Energetic analysis of succinic acid in water droplets: insight into the size-dependent solubility of atmospheric nanoparticles

Size-dependent solubility is prevalent in atmospheric nanoparticles, but a molecular level understanding is still insufficient, especially for organic compounds. Here, we performed molecular dynamics simulations to investigate the size dependence of succinic acid solvation on the scale of ~1-4 nm with the potential of mean forces method. Our analyses reveal 15 that the surface preference of succinic acid is stronger for a droplet than the slab of the same size, and the surface propensity is enhanced due to the curvature effect as the droplet becomes smaller. Energetic analyses show that such surface preference is primarily an enthalpic effect in both systems, while the entropic effect further enhances the surface propensity in droplets. On the other hand, with decreasing droplet size, the solubility of succinic acid in the internal bulk volume may decrease, imposing an opposite effect on the size dependence of solubility as compared with the enhanced surface propensity. 20 Meanwhile, structural analyses, however, show that the surface to internal bulk volume ratio increases drastically, especially when considering the surface in respect to succinic acid, e.g., for droplet with radius of 1 nm, the internal bulk volume would be already close to zero for the succinic acid molecule.


Introduction
Organic compositions in aerosol particles play a vital role in new particle formation and growth, as well as aerosol activation 25 and chemical reactions in the atmosphere (Matsumoto et al., 1997;Riipinen et al., 2012;Metzger et al., 2010;Kirkby et al., 2011;Zhang et al., 2004;Kroll et al., 2011). Organic compositions participate in chemical and physical processes in the internal bulk volume (referred to as the bulk of confined volumes except surface volumes) as well as at the surface area of aerosol particles, and the reactions at the surface are evidenced to be unique compared to the reactions in the internal bulk volume (Finlayson-Pitts, 2009;Donaldson and Valsaraj, 2010;Hayase et al., 2011;Martins-Costa et al., 2012;Monge et al., 2012). probabilities, and phase separation, is still not well understood experimentally and theoretically (Rudich, 2003;Donaldson and Vaida, 2006;Fuzzi et al., 2006;Rudich et al., 2007;Bluhm and Siegmann, 2009;Werner et al., 2016).
It has been found that particle size plays an important role in regulating the phase transition of atmospheric nanoparticles. Bulk thermodynamics have shown that chemical composition, relative humidity, temperature, and pressure are the main factors 35 controlling aerosol phase state, but they cannot fully explain the observed phenomenon for nanometer-sized particles (Biskos et al., 2006a;Biskos et al., 2006b;Virtanen et al., 2010;Krieger et al., 2012;You et al., 2012;Veghte et al., 2013;Cheng et al., 2015;Werner et al., 2016;Freedman, 2017). Several studies have demonstrated that the phase transitions of nanoparticles have significant size dependence. For instance, Biskos et al. (2006a) found that deliquescence and efflorescence occur at higher relative humidity for sodium chloride nanoparticles using Humidified Tandem Differential Mobility Analyzer (HTDMA) 40 experiments. And Cheng et al. (2015) developed Differential Köhler Analysis to explain the size dependence of deliquescence and efflorescence for the sodium chloride nanoparticles. The morphology of particles was also observed through experiments, where phase separation of aerosol particles is dependent on both particle size and components (You et al., 2012;Veghte et al., 2013;Cheng et al., 2015).
Size dependence of organic molecule solvation in the nanoparticles as one of the phase transitions has been suspected to be 45 prevalent. Favorable surface solvation of organic molecules was known to enhance the solubility of solutes at the surface of the nanodroplet with molecular dynamics (MD) simulation by Hub et al. (2012). Similar results were shown in Sayou et al. (2017) and attractive and repulsive interactions of the organic solute with water molecules were analyzed to explain the surface propensity of organic molecules in the droplet. Besides, Werner et al. (2016) probed the surface activity of succinic acid in sodium chloride and ammonium sulphate solutions by X-ray Photoelectron Spectroscopy (XPS) experiments. The result 50 suggests that succinic acid has a higher solubility at the surface of bulk solutions compared with the total solubility of bulk solutions because of fewer water molecules for hydration at the surface. Furthermore, Werner et al. (2016) developed a partitioned particle model, where the surface region used the surface concentrations estimated from XPS experiments and the internal bulk volume had the bulk-like particle concentration, indicating that the size-dependent solubility of surface-active organic molecules can further enhance cloud activation and chemical reactions in the atmosphere. Despite of the importance 55 of the size dependence of organic molecule solvation, how the solubility in the internal bulk volume as well as at the surface part of nanodroplets changes, particularly with changing droplet size and the underlying driving mechanism is still unclear.
In this study, size dependence of organic composition solvation in the nanoparticles is investigated based on the molecular dynamics (MD) simulation and energetic analysis. Succinic acid is chosen as a surrogate in the simulation because it is present in the atmospheric aerosol particles (Hyvärinen et al., 2006;Werner et al., 2014). The simulations are performed for one single 60 succinic acid molecule in different sizes of water nanodroplets (i.e., droplet radius from ~1 nm to 4 nm) and different thickness of water planar slabs (i.e., half slab thickness from ~1 nm to 4 nm). We compare the potential of mean forces (PMFs) derived from the simulations, as well as the enthalpic and entropic components decomposed from PMFs, to identify the potential driving forces for the size dependence of succinic acid solvation and discuss the size dependence of succinic acid solubility. https://doi.org/10.5194/acp-2020-1329 Preprint. Discussion started: 4 February 2021 c Author(s) 2021. CC BY 4.0 License.

PMF calculations
Potential of Mean Forces (PMFs), the free-energy profiles along the chosen coordinate (along the center of mass of water droplets or slabs to the air), were computed using the method of umbrella sampling with a weighted histogram analysis method (Kumar et al., 1995;Souaille and Roux, 2001;Hub et al., 2010). This method is based on the protocol described in Caleman et al. (2011) andSayou et al. (2017). 70 The MD simulations were computed with GROMACS 2016.5 (Van Der Spoel et al., 2005;Abraham et al., 2015). The solute (succinic acid molecule) was treated with the Generalized Amber force field (GAFF) with the electrostatic potential (ESP) calculated charges (B3LYP/aug-cc-pVTZ). Although aerosol pH in the atmosphere may vary in a large range (Zheng et al., 2020), the majority of fine particles (diameter < 1 m) is with pH below ~3-4 ( Fig. 2 of Pye et al. (2020)). In this pH range, the form of succinic acid is more relevant as hydrogen-succinate and succinate dominate only at higher pH. Besides, recent 75 study shows that the surface of pure water droplets and slabs is slightly enriched in hydronium, which might lead to some portion of succinic acid having a different protonation state in the internal bulk volume as on the surface as well (Hub et al., 2014), but the influence is not expected to be large to the aerosol pH. The solvent (water) was described with the TIP4P/2005 force field, as compared to four other water models the combination of TIP4P/2005 with above force field of succinic acid produce the closest hydration free energy (∆ hyd ) of single succinic acid molecule to the experimental value ( Fig. S1 and SI). 80 Newton's equation of motion was integrated with the leapfrog integration scheme at a time step of 2 fs (Van Gunsteren and Berendsen, 1988). The temperature was controlled at 298.15 K using the stochastic dynamics integrator at a coupling constant of 0.1 ps (Van Gunsteren and Berendsen, 1988). The covalent bonds to hydrogen atoms were fixed by Linear Constraint Solver (LINCS) algorithm (Miyamoto and Kollman, 1992;Hess et al., 1997).
The droplets without periodic boundary conditions were simulated as isolated systems (Fig. 1a). The simulations consisted of 85 a single succinic acid molecule with 140, 1117, 3768, and 8932 water molecules, corresponding to aqueous droplets with approximate radius 1 nm, 2 nm, 3 nm and 4 nm, respectively. The center of mass of the water molecules was fixed at the origin, and no cutoffs were applied to any intermolecular interactions (Fig. 1a).
The umbrella potential was implemented as where was set to 1000 kJ mol -1 nm -2 , is the radial distance of the solute center of mass from the center of mass of the water molecules. Along with the reaction coordinate , 20 windows with c of 0, 0.2, …, 4 nm, 25 windows with c of 0, 0.2, …, 5 nm, 35 windows with c of 0, 0.2, …, 7 nm, and 35 windows with c of 0, 0.2, …, 7 nm were prepared for the droplet systems (i.e., droplet radius from ~1 nm to 4 nm). Since the simulations were performed at the room temperature, to prevent individual https://doi.org/10.5194/acp-2020-1329 Preprint. Discussion started: 4 February 2021 c Author(s) 2021. CC BY 4.0 License.
water molecules from frequent evaporation from the droplet surface that may lead to an ill-defined droplet center of mass, a 95 spherical flat-bottom restraining potential was applied to each water molecule. The potential was given by where is the radial distance of the oxygen atom of the water molecule from the center of mass of the water droplet, fb is the radius of the sphere around the center of mass without any additional force, which was taken to be 0.5 nm larger than the radius of the water droplet. The force constant fb was set to 1000 kJ mol -1 nm -2 , and is the Heaviside step function. The energy of 100 each structure was minimized, and the umbrella sampling was then conducted for 51 ns for each window of the solute location.
After removing the first 1 ns of each trajectory for equilibration, the PMFs for droplet systems were computed by the integration of the mean force. An analytical entropy correction for the PMFs of the droplet systems (Neumann, 1980) was computed as where B is the Boltzmann factor, is the radial distance of the solute center of mass from the center of mass of the water 105 molecules. This correction is due to the fact that for the case when the succinic acid molecule is at a distance , it can be anywhere in the volume of a spherical shell of size 4 2 ∆ (where ∆ = 0.2 is the thickness of the shell) and hence the associated entropy, which scales as the logarithm of the volume of this shell times Boltzmann's constant, yields a term as written in Eq. (3).
The slab systems were simulated under the periodic boundary condition in x, y, z directions with the canonical (NVT) ensemble. 110 The single succinic acid molecule was solvated in the center of rectangular water box with edge lengths of 3 nm to the x and y directions and of 2, 4, 6 and 8 nm to the z direction. The slab systems ( Fig. 1b) consisted of 560, 1200, 1800, and 2400 water molecules, which correspond to 1, 2, 3, and 4 nm of half slab thickness, respectively. The systems were then elongated to 30 nm in the z direction to create liquid-vapor interfaces. The smooth particle-mesh Ewald (PME) method with a real-space cutoff of 1.35 nm, a spline order of 6, and relative tolerance of 10 -5 was used to the Coulombic interactions. The PME with a real-115 space cutoff of 1.35 nm was also used for the van der Waals interactions. The umbrella potential was implemented as where was set to 1000 kJ mol -1 nm -2 , z is the distance along the z direction of the solute center of mass from the center of mass of the water molecules. Along with the reaction coordinate z, 20 windows with c of 0, 0.2, …, 4 nm, 25 windows with where z is the distance along the z direction of the oxygen atom of the water molecule from the cell center. fb was taken to be 0.5 nm larger than the half of thickness of the slab in the direction. The force constant fb was set to 1000 kJ mol -1 nm -2 . 125 Figure 1c depicts the schematic free energy profiles of a succinic acid along the radius in the droplet system and along the zaxis in the slab system. As the molecule moves from the center of mass of water toward the gas phase, its free energy goes through a minimum at the water-vapor interface, then increases again, and finally reaches a plateau in the gas phase. The free energy difference of the transition of the solute from the gas to the internal bulk volume of liquid, ∆ ba , represents the free energy of hydration of succinic acid (Sayou et al., 2017). ∆ sb is the free energy of the transition from the internal bulk volume 130 to the liquid-vapor interface, and ∆ sa is the free energy of the transition from the gas phase into the liquid-vapor interface.
To separate the energetic contributions that govern surface propensity, the PMFs of droplet systems and slab systems were further decomposed into enthalpic and entropic components. Under the canonical (NVT) ensemble, the enthalpic part ∆ was computed as the average potential energy of the system in each simulation. The entropic component was computed as ∆ was further decomposed into water-water and solvent-water interactions in each system.

Structural analysis
For structural analysis, the water density profiles were calculated along the radius for the droplet systems and along the zdirection for the slab systems. For droplet, the simulation cell was divided into multiple shells with a thickness of 0.02 nm, and the density of each shell was statistically calculated based on the cumulative radial distribution functions of water 140 molecules. For slab, the simulation cell was divided into multiple slices with a thickness of 0.02 nm, and the density of each slice was statistically averaged over the simulation time. Both water density profiles from droplets and slabs were further fitted to a hyperbolic tangent function (Vega and de Miguel, 2007) as where L and V are densities of water and vapor, 0 is the position of the Gibbs-dividing surface, and is a parameter about 145 surface thickness.

Size dependence of surface propensity
In Fig. 2, the calculated potential of mean forces (PMFs) profiles show a significant surface propensity of succinic acid in the liquid phase with a minimum at the surface region for both water droplet and slab systems at all sizes. It indicates that the 150 succinic acid is more preferably located at the surface region than inside the internal bulk volume and the gas phase.
We further calculated the differences of free energy between succinic acid in the internal bulk volume, surface and air to identify the size dependence of surface propensity of succinic acid. As shown in Fig. 3a, ∆ sb , the free energy difference of the transition of succinic acid from the internal bulk volume to the surface, is negative in every droplet and slab system, reflecting the surface active of succinic acid. ∆ sb in the droplet system is smaller than that in the slab system of the same size, 155 suggesting that the surface propensity in the droplet system is more enhanced than that in the slab system due to the curvature effect. Similar results were observed in Sayou et al. (2017), where for ethane, benzene, methanol, and H2S, curvature dependence is appreciable and the magnitude is mainly related to the structure of organic molecules. Hub et al. (2012) also supported this view in their study. Besides, based on XPS experiments, Werner et al. (2016) found the enrichment of surface partitioning of succinic acid at the surface region of the bulk solutions. 160 Meanwhile, ∆ sa , the free energy difference of the transition of succinic acid from gas to the water surface, which is negative as well, is larger in the droplet system than that in the slab system of the same size (Fig. 3b), showing that the propensity of succinic acid at the surface layer (relative to the gas phase) in the droplet system is less than that in the slab system. Furthermore, the magnitude of the propensity is affected by the nature of the organic molecules as well as the temperature (Sayou et al., 2017). As a result, the decrease of ∆ sb and increase of ∆ sa from slab to droplet with the uniform size indicates 165 that succinic acid is "pushed" from the internal bulk volume to the surface region more strongly but "pushed" from the air to the surface region more slightly for the droplet systems. These two opposite forces further affect succinic acid to prefer to be exposed more to the air in the droplet rather than in the slab, which further demonstrates that the droplet curvature enhances the surface propensity of succinic acid.
When we compare the slab systems of different thicknesses, both ∆ sb (Fig. 3a, blue circles) and ∆ sa (Fig. 3b, blue triangles) 170 are almost constant with the increase of the thickness, suggesting a near constant surface propensity. In other words, the change of slab thickness only has a very slight effect on the surface propensity of succinic acid. This is different from the results of nanodroplet with a curvature/Kelvin effect. Without a curvature, the pressure inside the slab does not depend on the slab thickness leading to almost unchanged surface propensity. And at the considered temperature, this pressure to very good approximation is zero. However, for the droplet systems, when the radius of the droplet increases but the curvature of the 175 droplet decreases at the same time, ∆ sb decreases first from radius 1 nm to 3 nm but then increases (Fig. 3a, orange circles), showing that the surface propensity relative to the internal bulk volume is enhanced first and then weaken. ∆ sa shows a similar trend (Fig. 3b, orange triangles), meaning that the surface propensity relative to the air is also enhanced first and then https://doi.org/10.5194/acp-2020-1329 Preprint. Discussion started: 4 February 2021 c Author(s) 2021. CC BY 4.0 License.
weakened. Overall, the surface propensity is weakened as the droplet becomes larger and the difference between the droplet and slab systems becomes smaller when the droplet size increases. This result shows that different from thin slabs: for a nano-180 sized droplet, due to the Kelvin equation, the pressure inside the droplet is enhanced when the droplet size decreases.

Energetic analysis
To separate the energetic contributions that govern surface propensity of succinic acid, we decompose the PMF into the enthalpic and entropic components for both droplet systems (radius from ~1 nm to 3 nm) and slab systems (half thickness from ~1 nm to 3 nm) (Fig. 4). Energetic decomposition for larger droplets and slabs is excluded due to the large uncertainties in the 185 enthalpic component. Enthalpic curves are further divided into succinic acid-water interactions and water-water interactions ( Fig. 5). For both droplet and slab systems, enthalpies show a significant minimum at the surface, demonstrating that the surface preference of succinic acid in both systems is primarily an enthalpic effect. However, the enthalpic minimum at the surface is partially compensated by the entropy, indicating that the entropy affects the magnitude of surface preference of succinic acid. 190 On the contrary, the enthalpies show an increase in both droplets and slabs when being ~0.5 nm away from the surface, which may be caused by a reduced number of water-water interactions (Fig. 5, lines with open triangles). However, this increase of enthalpy is offset partially by the decreased entropy due to the increased rotational freedom of water molecules (Hub et al., 2012), resulting in a slight jump in the internal PMFs close to the surface. Furthermore, this jump becomes more significant when the size of the droplet decreases, showing that enthalpy plays a more important role in that region when the droplet size 195 decreases. In addition, enthalpy plays a leading role in the performance of the succinic acid molecule in the internal bulk volume of droplets as well as in the internal bulk volume of slabs. However, the internal enthalpy and entropy inside droplets fluctuate more intensely than those in the slab systems due to the more active molecule interactions (Fig. 5).

Discussion on size dependence of solubility
With XPS detected surface enrichment, Werner et al. (2016) predicted a higher solubility of succinic acid in the small 200 nanodroplet, based on the speculative assumptions that (1) with decreasing droplet size, the thickness of surface is assumed to be constant of ~0.4 nm, a value determined as the full width at half maximum of carbon profile peak with MD simulation in a slab system (3.93.93.9 nm 3 succinic acid aqueous solution in a 3.93.913 nm 3 slab, simulated for 1ns at 278.15 K) (Werner et al., 2014); (2) the surface region of the nanodroplet was supposed to have the same enhanced surface concentration of succinic acid in the bulk solutions; and (3) the internal bulk volume of the nanodroplet was recognized to have the same 205 solubility of succinic acid as that in the internal bulk volume of bulk solutions. With our results, we may revisit these assumptions and gain more insight into the size dependence of solvation of succinic acid in nanodroplets compared to slab systems. The free energy of hydration (∆ ba ) gives implications on the size dependence of solubility of succinic acid (Fig. 3c). As shown in Fig. S1, the free energy of hydration of succinic acid in bulk solutions ∆ hyd is −44.94 ± 1.4 kJ mol -1 by TI method 210 calculation based on the GAFF/ESP force field for succinic acid and the TIP4P/2005 force field for water molecules (SI). For slab systems, ∆ ba increases very slightly along with the increase of the slab thickness but not beyond the range of ∆ hyd (Fig.   3c, pink block). It means that the potential of succinic acid hydration into the slab systems is basically unchanged no matter how the slab thickness changes. This further indicates that the solubility of succinic acid in the internal bulk volume as well as at the surface region of the slab systems is almost constant. 215 Interestingly, for droplet systems, the hydration free energy ∆ ba , under the interactions of ∆ sb and ∆ sa by the curvature effect, are significantly larger than ∆ hyd at all sizes. The unfavorable succinic acid solvation in the internal bulk volume of the nanodroplets decreases the solubility of succinic acid in the internal bulk volume of the nanodroplet compared to the solubility in the internal bulk volume of the slabs. Furthermore, the hydration free energy is larger in the internal bulk volume of small droplets than that in large ones, suggesting a decrease of succinic acid solubility in the internal bulk volume of 220 nanodroplet when the droplet size becomes smaller.
Although a larger hydration free energy usually suggests a smaller solubility in water, i.e., a decrease of succinic acid solubility in the internal bulk volume of nanodroplet when the droplet size becomes smaller, the favorable surface solvation of succinic acid enhances the surface solubility of succinic acid in the nanodroplets. Therefore, the overall size dependence of succinic acid solubility with decreasing sizes of the nanodroplets is controlled by the competition between these two effects (i.e., 225 decreased solubility in the internal bulk volume and increased surface enrichment) as well as the drastic change of surface to internal bulk volume ratio.
To examine how the surface thickness and surface to internal bulk volume would change when the size of the systems decreases, we calculated the water density profiles (Fig. 2, dash lines) fitted by the hyperbolic tangent function (Vega and de Miguel, 2007;Wang et al., 2019) based on the distributions of water density from MD simulations (Fig. 2, open circles). The 230 surface thicknesses of different sizes of nanodroplets and slabs can then be obtained (Fig. 6a). It should be noted that the surface thickness of the slab profiles systematically depends on the lateral linear dimension of the slab, which is 3 nm here.
Due to the capillary wave broadening it diverges logarithmically with this linear dimension, while in contrast, no such ambiguity occurs for the surface thickness of droplets (Wang et al., 2019). The structural analyses show that the surface to internal bulk volume ratio in water droplet increases significantly with decreasing droplet size (Fig. 6c), although the surface 235 thickness of water decreases at the same time. When droplet radius changes from 4 nm to 1 nm The water surface thickness slightly decreases from ~0.34 nm to ~0.25 nm (Fig. 6a). This result is consistent with the previous researches involving nanodroplets and slab water (Thompson et al., 1984;Ismail et al., 2006;Vega and de Miguel, 2007;Wang et al., 2019).
However, how to define the surface is another important issue here. Regarding succinic acid molecule, thermodynamic analysis shows that when the droplet is smaller, enthalpy and entropy ranges in the internal bulk volume become narrower. As shown in Fig. 5, for droplets of all sizes, the interaction of succinic acid molecule and water molecules shows a significant increase in both droplets and slabs when being ~1 nm away from the surface, indicating the surface thickness is effectively ~1 nm in respect to succinic acid molecule and is independent to the change of droplet size. Moreover, relatively stable interactions exist in the internal bulk volume of droplets with radius of 2 nm (Fig. 5b) and 3 nm (Fig. 5c), but not in the droplet with radius of 1 nm (Fig 5a). In this regard, for droplet with radius of 1 nm, the internal bulk volume is close to zero for succinic acid molecule, 245 but in terms of water molecule, the radius of internal bulk volume is still ~0.75 nm (Fig. 6b). Thus, in respect to succinic acid molecule, the surface to internal bulk volume ratio would increase much stronger than that for water molecule. These results suggest that the width of the PMF profiles is not simply controlled by the width of the interface, what also matter is the effective range of intermolecular forces between the succinic acid molecule and the water molecules. Since these forces contain longrange contributions, it is understandable that a saturation of the PMF profiles in the center of the droplet or slab in Fig. 2 (at 250 the same horizontal plateau) is not observed, unlike the behavior of the density profiles.

Conclusion
In conclusion, the size dependence of succinic acid solvation was probed based on the MD simulations and energetic analysis.
The simulations were conducted for a single succinic acid molecule in the different sizes of water nanodroplets (i.e., droplet radius from ~1 nm to 4 nm) and different thickness of water planar slabs (i.e., half slab thickness from ~1 nm to 4 nm). On the 255 one hand, energetic analyses reveal a stronger surface preference for a droplet than the slab of the same size, and the surface propensity is enhanced due to the curvature effect as the droplet becomes smaller. The energetic driving force underlying such surface preference is primarily an enthalpic effect in both systems, while the entropic effect further enhances the surface propensity in droplets. On the other hand, the energetic analyses gives implications on the size dependence of solubility of succinic acid in nanodroplets. The solubility of succinic acid in the internal bulk volume may decrease when the droplet 260 becomes smaller, imposing an opposite effect on the size dependence of solubility as compared with the enhanced surface propensity. Meanwhile, the overall size dependence of succinic acid solubility is also controlled by the drastic change of surface to internal bulk volume ratio. However, structural analyses show that the surface to internal bulk volume ratio increases drastically, especially when considering the surface in respect to succinic acid, e.g., for droplet with radius of 1 nm, the internal bulk volume would be already close to zero for the succinic acid molecule. 265

Data availability
Readers who are interested in the data should contact the authors: Yafang Cheng (yafang.cheng@mpic.de), or Chuchu Chen (c.chen@mpic.de).

Supporting Information
The supporting information is available free of charge on the website. 270

Author contributions
Y.C. and H.S. conceived and led the study. C.C. performed the MD simulation and analyzed the data. D.S. and M.G provided the force field of succinic acid. C.C., Y.C., K.B., H.S., and X.W. interpreted the results. K.B., U.P., D.S., and M.G. discussed the results and commented on the manuscript. C.C., Y.C., and H.S. wrote the manuscript with inputs from all coauthors.

Competing interests 275
The authors declare no competing financial interest.