Miguel De La Puente

Damien Laage
email: damien.laage@ens.psl.eu

How the acidity of water droplets and films is controlled by the air-water interface

Acidity is a key determinant of chemical reactivity in atmospheric aqueous aerosols and in water microdroplets used for catalysis. However, many fundamental questions about these systems have remained elusive, including how their acidity differs from that of bulk solutions, the degree of heterogeneity between their core and surface, and how the acid-base properties are affected by their size. Here, we perform hybrid DFT-quality neural network-based molecular simulations with explicit nuclear quantum effects and combine them with an analytic model to describe the pH and self-ion concentrations of droplets and films for sizes ranging from nm to µm. We determine how the acidity of water droplets and thin films is controlled by the properties of the air-water interface and by their surface-to-volume ratio. We show that while the pH is uniform in each system, hydronium and hydroxide ions exhibit concentration gradients that span the two outermost molecular layers, enriching the interface with hydronium cations and depleting it with hydroxide anions. Acidity depends strongly on the surface-to-volume ratio for system sizes below a few tens of nm, where the core becomes enriched in hydroxide ions and the pH increases as a result of hydronium stabilization at the interface. These results obtained for pure water systems have important implications for 1 our understanding of chemical reactivity in atmospheric aerosols and for catalysis in aqueous microdroplets.

Introduction

Many chemical reactions of major importance take place in aqueous systems with high surface-to-volume ratios. These include, for example, climate-critical atmospheric reactions at the surface of aerosols, 1,2 reactions presenting a dramatic rate-enhancement in microdroplets, 3,4 and the distinct reactivity of systems in nanoconfinement. 5 In these systems, pH is a key determinant of chemical reaction pathways and kinetics. [6][7][8][9] However, questions such as the extent to which their pH differs from that of bulk solutions and how it is impacted by the presence of the interface have remained elusive.

Previous extensive studies have concluded that the acidity at the air/water interface is different from that in the bulk. [10][11][12][13][14][15][16][17][18][19][20] However, these studies have mostly considered the air surface of a bulk solution. Characterizing the acidity of small droplets and thin films not only presents many experimental challenges, but also raises specific questions. 7 First, the relative influence of the interface is expected to be enhanced in these smaller systems, leading to a dependence of acidity on their size. Second, in the case of microdroplets, acidity has been suggested to be non-uniform 21,22 and to change between surface and core due to their different solvation conditions.

Since traditional pH probes for bulk solutions, such as reference electrodes, are obviously not suitable for very small droplets, experiments typically estimate the pH indirectly, using spectroscopic methods to measure the ratio between the acidic and basic forms of an indicator molecule. Recent pioneering experiments have used Raman microscopy to map the acidity within microdroplets and to investigate whether pH gradients exist within aerosols. [21][22][23] However, these studies have yielded conflicting results and a consensus on this key question is still lacking. Furthermore, it is important to immediately emphasize an important limitation of indicator-based spectroscopic approaches: they cannot provide an unambiguous measure of pH, because the acid-base equilibrium also depends on other factors, including the indicator pKa which can significantly change from its bulk value due to partial desolvation of the acid and conjugate base at the interface. 24,25 to identify its molecular-level determinants. However, in order to describe the water selfdissociation reaction equilibrium that controls the pH and to account for the quantum nature of the high-frequency proton motions, simulations must combine high-level electronic structure calculations with an explicit description of nuclear quantum effects. This results in a very high computational cost which severely limits the length of simulations and the precision of the calculated water self-dissociation equilibrium. However, this limitation can now be overcome with deep neural network potentials (NNPs), which are trained to provide DFT-quality reactive force fields at a fraction of the computational cost and which can be combined with path-integral molecular dynamics to account for nuclear quantum effects. 26 Here, we study the acidity of water droplets and films, focusing on pure water samples to determine how their size and surface-to-volume ratio affect their acidity. We first use NNP-based path-integral molecular dynamics simulations of water slabs to determine how the water self-dissociation equilibrium and the hydronium and hydroxide self-ion stabilities change as the air-water interface is approached. We then incorporate these results into a general analytic model to determine the pH and self-ion concentration profiles in droplets and films of decreasing size.

Our model reveals how the air-water interface affects the acidity of slabs and droplets.

Consistent with previous studies, 11,[14][15][16][27][28][29] our results show that the air-water interface is enriched in hydronium ions. However, in interfacial systems this hydronium concentration increase cannot be straightforwardly interpreted as a decrease in pH. A critical point is that while pH and hydronium concentration can be considered interchangeably to discuss the acidity of bulk solutions, these two quantities behave very differently in presence of an interface. The latter makes these systems heterogeneous and induces opposite gradients in the hydronium concentration and thermodynamic activity coefficient. Acidity can then be considered from two different perspectives. From a thermodynamic point of view relevant for example for equilibrium chemical composition, the pertinent quantity is the pH: it is governed by the chemical potential and is therefore uniform throughout the system at equilibrium. 30 In contrast, for reaction rate constants and for experiments whose signal is sensitive to the abundance of a self-ion, the hydronium concentration should be considered, and the latter exhibits a gradient between core and surface due to the different local solvation properties. We show that, unexpectedly, the pH increases in droplets and slabs of decreasing size, although the surface is always enriched in hydronium ions compared to the core. This pH increase in smaller systems is mostly due to the stabilization of hydronium ions at the interface which lowers their chemical potential, with only a minor contribution from the reduced dissociation equilibrium constant. Finally, the implications for chemical equilibria and reaction kinetics in atmospheric aerosols and microdroplets used in catalysis are discussed.

Methods

Our molecular dynamics simulations employ a NNP built using the DeePMD-kit code 31 on data referenced at the revPBE0-D3 [32][33][34] hybrid density functional level of theory (DFT), whose adequacy for describing the properties of bulk 35 liquid water and the air-water interface 36 has been established by previous studies. Our training set was generated through an iterative process and contains more than 11,000 configurations describing multiple system sizes and covering each stage of the water self-dissociation process. The final converged NNP reproduces the revPBE0-D3 energies and forces with a computational cost reduced by more than 20,000 times.

Another key feature of our simulations is the explicit consideration of nuclear quantum effects, the importance of which for the water pH is evidenced by the pronounced pH increase between H 2 O and D 2 O heavy water (from pH=7.0 to pD '7.4). Nuclear quantum effects for all nuclei in our simulated system are described by the PIGLET [START_REF] Ceriotti | Efficient First-Principles Calculation of the Quantum Kinetic Energy and Momentum Distribution of Nuclei[END_REF] combination of pathintegral molecular dynamics with a generalized Langevin equation for improved convergence.

Our simulations are propagated at 300 K with the i-PI [START_REF] Kapil | i-PI 2.0: A universal force engine for advanced molecular simulations[END_REF] simulation package for air-water slab and bulk reference systems all containing 128 water molecules for a total cumulative simulation time of approximately 20 ns. Further details regarding the simulation methodology, the training set preparation, the NNP validation and the impact of nuclear quantum effects are provided in the SI.

Results and Discussion

We first determine how the air-water interface impacts the water self-dissociation reaction that produces the hydronium and hydroxide ions governing the acidity. We perform NNP molecular dynamics simulations of bulk water and of a slab of water with a liquid-vapor interface, and we calculate the dissociation reaction free energy profile at different depths relative to the interface and in the bulk for reference. We select a water molecule that we maintain at a certain depth, and we force its dissociation using a biasing potential (see details in SI). This yields a hydroxide ion kept at the same depth and a proton that is free to diffuse in the system (no assumption is made on the solvation structure of the proton and whether it adopts the Eigen or Zundel limiting forms). The water self-dissociation free energy is then determined as the free energy difference between the stable undissociated and dissociated states at a given depth.

Figure 1 shows that water self-dissociation is more unfavorable at the interface than in the bulk: the dissociation free energy increases by approximately 2 kcal/mol when moving from the bulk to the interface. This change is localized in the immediate vicinity of the interface and the bulk reaction free energy is recovered beyond two molecular layers from the interface. Further analysis of the trajectories used for the free energy calculations shows that our simulations are long enough to allow for the pair of self-ions to dissociate (see Figs. S16-S17) and confirm that the self-dissociation mechanism at the interface is unchanged with respect to the mechanism in the bulk. Figure 1: Water self-dissociation free energy profiles along the hydrogen coordination number around the reactive oxygen atom in the bulk (black) and at two different heights s relative to the air-water interface, respectively at the interface (green) and about one molecular layer below the interface (orange). Heights s are defined as the distance from the reacting oxygen atom to the plane of the Gibbs dividing surface, oriented such that negative distances indicate the liquid phase. The dissociating water molecule and the resulting hydroxide are fixed at height s while the proton-accepting molecule and the hydronium product are free. Thick colored lines represent 95% confidence intervals.

We therefore describe the change in the water dissociation free energy between its bulk standard state reference r G and its r G(s) value in the slab at a depth s relative to the interface as arising from the changes in solvation free energies for the reactants and products, using the thermodynamic cycle presented in Fig. 2a,

r G(s) r G = G H 2 O hyd (s) h G H 2 O hyd i + G HO hyd (s) + h G H 3 O + hyd i (1)
This scheme corresponds to the dissociation simulated in Fig. 1: the proton-donating water molecule is held at depth s, leading to a hydroxide ion at the same position, while the protonaccepting water molecule and the resulting hydronium ion are free to diffuse. Equation 1thus involves the changes in hydration free energy between the standard state and at depth s

in the slab for H 2 O and OH , G H 2 O,OH hyd (s) = G H 2 O,OH hyd (s) G H 2 O,OH hyd
, and the averages over the sample of the hydration free energy changes with respect to the standard state for

H 2 O and H 3 O + , h G H 2 O,H 3 O + hyd i (see eq. S4 in SI).
The interfacial stabilities of both water self-ions are therefore essential to understand changes in the self-dissociation equilibrium. Many studies have addressed this fundamental question, 10,12,[17][18][19][20]27,[START_REF] Wick | The behavior of NaOH at the air-water interface: A computational study[END_REF][START_REF] Mundy | Hydroxide anion at the air-water interface[END_REF][START_REF] Jagoda-Cwiklik | Behavior of the Eigen Form of Hydronium at the Air/Water Interface[END_REF] with contrasting results that are highly dependent on the method employed. Here, we compute the changes in hydration free energies for the reactants (H 2 O) and products (H 3 O + , OH ) as they move across the interface into the water phase from a series of independent biased NNP simulations on a slab system (see SI). In the following, we consider that the bulk standard-state hydration free energies can be approximated by the hydration free energies in the middle of the slab (see SI).

Our NNP simulations at the revPBE0-D3 level with nuclear quantum effects show that the hydronium cation is stabilized by -1.2±0.4 kcal/mol at the air-water interface with respect to the bulk (see Fig. 2b). This result is in excellent agreement with the -1.3±0.2 kcal/mol 29 and -0.9±0.2 kcal/mol 28 values determined by recent SFG measurements, and is consistent with previous calculations at different levels of theory. [15][16][17]25 The behavior of the hydroxide anion strongly contrasts with that of hydronium, since our simulations show that it is destabilized by +2.7 ± 0.7 kcal/mol at the interface. This repulsion was already observed in previous simulation studies performed with empirical valence bond 16 and polarizable 15 force fields and is consistent with the very low hydroxide surface activity measured by SFG. [START_REF] Das | The Surface Activity of the Hydrated Proton Is Substantially Higher than That of the Hydroxide Ion[END_REF] Further analysis of our trajectories show that the molecular factors causing the different surface affinities of hydronium and hydroxide in our NNP simulations are similar to these identified in previous computational studies. Hydronium ions behave as amphiphilic molecules and are stabilized at the interface due to the loss of a weak and destabilizing accepted hydrogen bond upon desolvation and the gain of almost a full hydrogen bond for the water molecule that is replaced at the interface 11,16 (see Figs. S6-S7). In contrast, the strong interfacial destabilization of hydroxide arises from the loss of stabilizing interactions with water (loss of a weak donated hydrogen bond and weakening of a strong accepted hydrogen bond, see Figs. S6-S8) 16 and the reduction of orientational entropy 15 (Fig. S9).

We then use eq. 1 to determine the r G(s) water self-dissociation reaction free energy continuously along the distance s to the interface, from the bulk standard-state reaction free energy r G (Fig. 1) and the hydration free energy profiles (Fig. 2b). Figure 2c shows that the model eq.1 provides an excellent description of the reaction free energy calculated explicitly at some selected depths by forcing the dissociation in Fig. 1. This strongly suggests that changes in water self-dissociation free energy at the interface are essentially determined by solvation effects.

This solvation picture of interfacial effects on water dissociation free energy (eq. 1) is thus supported by all-atom NNP molecular dynamics simulations of a paradigm air-water interface. We now use it to propose an analytic model capable of addressing systems not accessible via molecular simulations to describe the pH, hydronium and hydroxide ion concentrations in aqueous systems of various sizes and various shapes, including spherical droplets and slabs.

The spherical geometry is of crucial interest for atmospheric aerosols 1,2 and microdroplet catalysis, 3,4 while the slab geometry is relevant for, e.g., aqueous systems nanoconfined within hydrophobic slits. [START_REF] Muñoz-Santiburcio | Nanoconfinement in Slit Pores Enhances Water Self-Dissociation[END_REF] In the following, we show that in these interfacial systems pH is uniform, 30 but unlike in the bulk, local self-ion concentrations are not and present large differences between the core and the surface of the solution. We then determine the pH value and the self-ion concentration profiles as the system size changes, and discuss the impact of the interface.

By definition, the pH is determined by the thermodynamic activity of hydronium ions, which can be expressed as the product of the hydronium concentration and of the hydronium activity coefficient as

pH = log 10 (a H 3 O + ) = log 10 ✓ H 3 O + (s) [H 3 O + ] (s) c ◆ ( 2
)
where c =1 mol/L is the standard state concentration. The thermodynamic activity is determined by the chemical potential µ as a = exp(µ/k B T ), where k B is the Boltzmann constant and T the temperature. Since the chemical potential is the same everywhere in the system at equilibrium, the pH is unique and uniform in each system.

However, local self-ion concentrations are not uniform because the ions' hydration free energies depend on the distance s to the interface: ions tend to accumulate in regions where they are better solvated, and the local hydronium ion concentration follows the Maxwell-Boltzmann distribution proportional to exp

⇣ G H 3 O + hyd (s)/k B T ⌘
, where we recall that G H 3 O + hyd (s) is the hydronium hydration free energy change between the standard state and the distance s from the interface.

The thermodynamic activity is constant over the system because a local increase in concentration is compensated by a decrease in the activity coefficient. This can be shown by using the solvation thermodynamics perspective [START_REF] Ben-Naim | Solvation Thermodynamics[END_REF] to determine the hydronium chemical potential

µ H 3 O + = µ H 3 O + + k B T ln ⇥ H 3 O + ⇤ (s)/c + G H 3 O + hyd (s) (3)
Here, we assume that at the low hydronium (and hydroxide) concentrations being considered, the only deviation from ideality arises from the distinct solvation properties of the interfacial system, and not from ion-ion correlations. Eqs. 2-3 imply that (s) = exp

⇣ G H 3 O + hyd (s)/k B T
⌘ , which exactly compensates the spatial variations of the hydronium concentration (Fig. S19).

We now determine how the pH and the self-ion concentrations are affected by the system size. Our model relies on the following three key ingredients. First, the ions are distributed according to the Maxwell-Boltzmann probability based on their local hydration free energy.

Second, charge conservation in these closed systems requires the total numbers of hydronium and hydroxide ions are equal. Finally, the dissociation probability of a water molecule at a distance s from the interface is determined by the dissociation reaction free energy r G(s)

at that position (eq. 1) and the total number of hydronium ions in the system is obtained by averaging this dissociation probability over all distances, weighted by the probability to find a water molecule at this distance. We approximate the self-ion hydration free energy profiles by their numerical determinations in the above-described slab system and consider that hydration free energies do not change significantly beyond 7 Å from the interface (see SI for details and validation of these approximations).

As detailed in the SI, our model thus yields the pH of aqueous systems with a liquid-vapor interface,

pH = log 10 ✓ h[H 3 O + ]i c ◆ D G H 3 O + hyd E k B T ln 10 (4)
The first term on the right-hand side is the familiar bulk term involving the hydronium ion concentration, here averaged over the sample. The interface changes the solvation properties and introduces an additional term involving h G H 3 O + hyd i, the change in hydronium hydration free energy with respect to the standard state, averaged over the sample, the boundaries of which are defined by the Gibbs dividing surface.

The average hydronium concentration is (see SI):

⌦ [H 3 O + ] ↵ = c q h⇢ H 2 O i/⇢ H 2 O exp [ h r Gi/k B T ] (5)
where h r Gi is the water self-dissociation reaction free energy averaged over the positions of both reactants and both products and h⇢ H 2 O i/⇢ H 2 O is the ratio of the average water density in the system and of the bulk liquid water density in standard conditions. The change in water density between the standard conditions and the slab being very small, the structure of eq. 5 is very similar to that of the bulk concentration obtained from the bulk reaction free energy via the law of mass action in ideal conditions. It shows that in an interfacial system, the average hydronium concentration differs from its bulk water value essentially because of the change in the average water dissociation free energy.

Finally, the hydronium local concentration profile along the distance to the interface is obtained from the thermal distribution in presence of the hydration free energy profile as

[H 3 O + ](s) = ⌦ [H 3 O + ] ↵ exp G H 3 O + hyd (s) h G H 3 O + hyd i k B T ! (6)
We note here that while our numerical validation of the solvation model eq. 1 in Fig. 2 had broken the symmetry of the water self-dissociation reaction by fixing the water molecule that model eqs. 4-6. Straightforward extension of eqs. 4-6 to hydroxide ions provides the pOH and the hydroxide concentration profiles (see SI).

We now use eq. 4 to determine the pH of neat water spherical droplets and neat water slabs of decreasing sizes. Determining absolute pH values from simulations is challenging since in the bulk reference, limitations of the selected DFT level underestimate the dissociation free energy ( r G =13.1 kcal/mol) and yield pH'4.8 for bulk neat water instead of the wellknown experimental value of 7. Our present study thus instead focuses on pH shifts relative to the bulk value: the latter are determined by interfacial solvation of water self-ions, for which we have shown that our simulations provide a quantitative agreement with experiments, including SFG. 28,29 Strikingly, the results in Fig. 3a-b show that the pH increases, i.e., the solution becomes less acidic, when the system size decreases and the surface to volume ratio increases. This effect becomes significant for droplets smaller than 100 nm in radius and slabs thinner than 200 nm. For the smaller systems, e.g., spherical droplets with a radius of approximately 1 nm which are relevant for atmospheric aerosol nucleation, the pH increases by close to a full pH unit with respect to the bulk value.

The situation in interfacial systems is more complex than in the bulk, because the interface changes not only the degree of dissociation, which affects both pH and pOH similarly, but also the stabilization of each self-ion, which affects separately pH and pOH. It is therefore instructive to complement our pH results with pOH calculations. Figures 3c-d show that the pH increase in smaller systems is accompanied by a pOH decrease, indicating that these systems become increasingly basic.

While this may at first seem surprising that the presence of an interface where the concentration of hydronium ions is enhanced leads to a system which is less acidic and more basic, the decomposition of the pH shift provided by eq. 4 reveals the molecular origin of this effect. The interface can affect the pH in two ways: first by changing the average degree of water dissociation which determines the average hydronium concentration, and second via the ion hydration free energy by providing a favorable solvation environment for the hydronium. While the degree of dissociation is influenced by both hydronium and hydroxide solvation at the interface (eq. 1) and affects pH and pOH in the same way, the hydration free energy term is specific to each ion. Figure 3a-d show that the reduction in hydronium and hydroxide overall concentrations caused by the more unfavorable dissociation at the interface increases both pH and pOH, leading to an increase in pKw, but that it only accounts for a minor fraction of the pH and pOH shifts. Instead, our calculations show that most of the pH increase in these systems with high surface-to-volume ratio is due to the stabilization of the hydronium ions at the interface, which lowers their chemical potential, decreases their thermodynamic activity and increases the pH. Vice-versa, the destabilization of hydroxide ions at the interface increases their activity and decreases the pOH.

The pH increase in small samples is a general result that applies to all neat water systems with large surface to volume ratio. Some experiments [21][22][23] have recently used Raman microscopy techniques to determine the acidity of water droplets. However, they have employed a pH probe and measured the acid (AH) and conjugate base (A ) abundances. In these interfacial systems, the Henderson-Hasselbach equation

pH = pKa + log 10 A [A ] AH [AH] (7)
involves concentrations and activity coefficients which depend on the distance to the interface. The acid/base concentration ratio therefore does not directly report the pH since the pKa is modified with respect to its bulk value by the interface, 24,[START_REF] Wellen | Surface pK a of octanoic, nonanoic, and decanoic fatty acids at the air-water interface: applications to atmospheric aerosol chemistry[END_REF][START_REF] Eugene | Enhanced Acidity of Acetic and Pyruvic Acids on the Surface of Water[END_REF] and the spatiallydependent hydration free energy contribution to the activity coefficients needs to be considered. We further note that most of these experiments consider systems whose typical size exceeds that where our calculations predict the pH increase to be important.

We now examine the spatial distribution of hydronium and hydroxide ions in aqueous systems of various sizes and various shapes. While the pH is the relevant measure of acidity for the composition of a system at thermodynamic equilibrium, the local concentration in hydronium ions is a more useful quantity for chemical reaction kinetics which are of great importance, e.g., in the atmospheric chemistry context.

The concentration profiles calculated via eq. 6 in Fig. 4a show that although the total numbers of hydronium and hydroxide ions in the system are equal, the two self-ions have very different distributions. Both self-ions exhibit concentration gradients between the interface and the core of the aqueous system and these gradients have opposite signs for the two ions, as anticipated from the hydration free energies (Fig. 2b): the interface attracts hydronium ions but repels hydroxide ions. The interface is thus always enriched in hydronium ions but depleted in hydroxide ions compared to the core, which exhibits the opposite behavior due to charge conservation in the system. This leads to a strong excess (by more than two orders of magnitude) of hydronium ions relative to hydroxide ions at the interface (Fig. 4b). For both ions the gradient is localized and extends over no more than two molecular layers from the interface. A change in the system size or shape affects the average ion concentration and the average hydration free energy, which causes a uniform scaling of the concentration profile (leading to a shift on the logarithmic scale in Fig. 4a). Because our model assumes that hydration free energy profiles along the distance to the interface are not affected by the system size, the ratio between the interface and core concentrations is unchanged. Figure 4 shows a strong spatial heterogeneity between the core and interface concentrations. This suggests that adsorbed species on aqueous aerosols undergoing kinetically controlled reactions will display different reactivities depending on their surface activity and on the size of the aerosol. We now characterize this heterogeneity between interface and core, and how it depends on the surface-to-volume ratio of the aqueous phase. Figure 5 reports the changes in the hydronium and hydroxide concentrations relatively to the bulk both for the core and interfacial regions of systems of decreasing surface-to-volume ratios (or equivalently for a spherical droplet of increasing radius). Below a 0.1 nm 1 surface/volume ratio (i.e., above a 30 nm radius for a spherical droplet), the core and surface concentrations do not depend on the size or shape of the sample and have the values found for a macroscopic water phase with an air interface, i.e., the core is bulk-like and the interface is enriched in hydronium ions and depleted in hydroxide ions compared to the bulk. For smaller systems, as the surface-to-volume ratio grows, the hydronium and hydroxide concentrations become increasingly different in the core (Fig. 5a), while they become more similar at the interface (Fig. 5b): the core is progressively enriched in hydroxide ions and depleted in hydronium ions, thus becoming increasingly different from bulk water, while at the interface, the hydronium concentration decreases and the hydroxide concentration increases. This latter result may at first seem surprising: while the interface preferentially stabilizes hydronium ions, it is less enriched in hydronium ions when its relative weight in the system increases. This effect arises from the core. While in systems with vanishing surface-to-volume ratios the core acts as an infinite source of hydronium ions and an infinite sink of hydroxide ions, in small closed systems the weight of the core is reduced relative to that of the interface, so that it can no longer be considered as a reservoir of hydronium ions (see SI).

These results have major consequences for the understanding of chemical reactivity at interfaces under kinetic control. At the interface of macroscopic aqueous systems, hydronium ions are almost an order of magnitude more abundant than in the bulk, and almost three orders of magnitude more abundant than hydroxide anions at the interface, which clearly makes the local reactivity acidic. In small aerosols, hydrophilic species that preferentially reside in the core experience an environment enriched in hydroxide ions and follow basecatalyzed reaction pathways, while more hydrophobic and surface-active species are expected to follow acid-catalyzed pathways.

Conclusion

Our neural-network based molecular simulations of hybrid DFT quality with explicit nuclear quantum effects combined with an analytic model describing the pH, pOH and self-ion concentrations have determined how the acidity of water droplets and thin films is controlled by the properties of the air-water interface and by their surface-to-volume ratio. Our results show that the description of acidity in these heterogeneous systems is more complex than in the bulk. In closed aqueous systems such as droplets, the number of hydronium and hydroxide ions (determined by the mean dissociation free energy of water) is the same, but the ions are distributed differently due to differences in solvation. While the interface attracts hydronium cations, it repels hydroxide anions. A striking result is that as the size of the system decreases, the pH becomes less acidic while the pOH becomes more basic. This is due not so much to the small change in the self-dissociation mean free energy, but to the stabilization of hydronium ions by the interface, which lowers their chemical potential and raises the pH (and vice-versa for the interfacial destabilization of hydroxides, which lowers the pOH). The pH is uniform in each system, because it is determined by the hydronium chemical potential, which is the same at the interface and in the core. Similarly, the pOH is also uniform. In contrast, the hydronium and hydroxide concentrations show gradients that extend over several molecular layers at the interface. Both pH and concentration profiles depend on the size of the aqueous phase, reaching their macroscopic liquid-vapor interface limit for sizes above a few tens of nm. While we have focused here on the acidity of pure water systems, our model can be easily extended to situations in which multiple acidic species are present. These future developments to include other compounds affecting the acidity will have major implications for catalysis in microdroplets and for chemical reactivity in atmospheric aerosols which have a profound effect on global climate 1 and whose properties are radically affected by their acidity. [6][7][8][9][START_REF] Angle | Acidity across the interface from the ocean surface to sea spray aerosol[END_REF][START_REF] Angle | Direct Quantification of Changes in pH within Single Levitated Microdroplets and the Kinetics of Nitrate and Chloride Depletion[END_REF] Deep Neural Network Potentials

Neural Network Potential Training Parameters

Neural Network Potentials (NNP) were trained with the DeePMD-kit 4 code, version 2.0.3, and followed a typical two-neural-network architecture. First, atomic configurations are processed by an embedding network of three layers of 20, 40 and 80 neurons, which builds atomic descriptors of the environment. These descriptors are then given as input to a fitting network of four layers of 240 neurons that computes the atomic forces in the configuration.

The embedding follows the smooth-edition scheme 10 with conservation of angular and radial information within a cutoff of 6 Å. A cosine weighting function is applied for atoms beyond 5.5 Å to ensure a smooth cutoff. Although no explicit treatment of long range effects was included in the NNP workflow, atomic configurations at the air-water interface were added to the training set, which has been shown to be enough to correctly reproduce the liquid-vapor equilibrium with a short-range approach. 11 The parameters of the neural networks are set during the training procedure by minimizing a loss function based only on the mean squared error on the atomic forces predicted by the network (compared to the reference data). While it is customary to include the error on atomic energies as well within the loss function, we found this to be counter-productive when including systems with different numbers of atoms in bulk and at the interface in the training set. We show below that this approach yields correctly converged NNP. The minimization is performed iteratively, with a total number of batches ranging from 4•10

S5

Iterative Construction of the Training Set

The training set over which the NNP are trained is an ensemble of atomic configurations (the inputs) for which the atomic forces (labels) have been computed at a given level of theory, taken as the reference. In order to build a set of configurations that is representative of the phase space region of interest of the system, a homemade concurrent learning procedure inspired by the DP-GEN scheme 12 was used. This procedure proceeds by iterations composed of three phases: training, exploration and labelling.

In the training phase, three NNP are trained over the full available data set, with different random parameter initializations. These NNP are then used to propagate MD trajectories in the exploration phase, which allows to generate new atomic configurations. Since the self-dissociation reaction is a rare event in bulk water under standard conditions we performed biased simulations (steered molecular dynamics and metadynamics) to ensure that the training set contains enough reactive configurations to correctly reproduce the behavior of water far from the equilibrium states (reactants and products). The maximal deviation between the atomic forces computed by the three NNP over a given configuration is used as a measure of the error on said configuration. If this error is larger than a certain threshold (here, 0.1 eV/Å), but small enough that the configuration might be physically relevant (inferior to 0.5 eV/Å), it is considered as a candidate, as its description by the NNP can still improve. The maximum error bound is set because, since there are no physical guidelines to the NNP description, the error can grow indefinitely large, as any conceivable configuration can be generated when the MD trajectory is outside of the phase space region on which the training was performed. After the exploration phase, a maximum number of 500 random configurations are selected for labelling: the atomic forces are computed at the reference level of theory and the inputs (configurations) and labels (forces) are included in the training set for the next iteration. Iterations are performed until there are no candidate configurations after an exploration phase of at least 100 ps. It is important to note that, for every type of system that must be ultimately described by the NNP (chemical composition, bulk or S6 interface, etc.), this procedure must be repeated until convergence.

In this work, the reference level of theory was the hybrid DFT with dispersion corrections revPBE0-D3, [13][14][15] with GTH pseudo-potentials 16 for describing the core electrons and a triple-⇣ quality basis set (TZV2P 17 ) for the valence electrons. The ADMM approximation was used in conjunction with the FIT3 auxiliary basis set. 18 All the electronic structure calculations were performed with the CP2K 19 code.

In order to obtain NNP that are able to perform simulations with explicit accounting of NQE via path-integral MD, we found it essential to use the configurations of the ringpolymer "beads" as inputs, instead of using only "centroid" configurations (or configurations obtained with classical nuclei MD).

Performance and Validation of the Neural Network Potential

NNP-based MD simulations were found to be ⇡20,000 times more efficient than DFT-based simulations with the machines at our disposal. For a system of 128 water molecules and a timestep of 0.25 fs, production rates of ⇡ 2 ns/day and ⇡ 0.15 ns/day were obtained for classical nuclei and path-integral MD respectively on 1 Nvidia Tesla V100 SXM2 GPU with 16 GB of RAM. These efficient production rates made possible to simulate a total accrued time of ⇡ 20 ns and ⇡ 60 ns with path-integral and classical nuclei MD respectively.

In order to validate the description given by the NNP, we generated an independent set S7 S3, the NNP yielded results that were very similar for every type of system considered in this study (bulk vs interface, acidic vs basic vs neutral, classical nuclei vs path-integral, etc.), and in excellent agreement with those of the reference method.

Furthermore, the learning curve presented in Fig. S4 shows no overfitting (the difference in relative error over the training and validation sets does not increase with training set size), suggesting that every region of chemical space was identically (and correctly) described.

Finally, in order to ensure that all stages of the self-dissociation process are correctly reproduced, we extracted configurations from different windows of our self-dissociation umbrella sampling simulations (see below for further details about these calculations) and computed the relative error of the NNP predictions with respect to DFT. As can be seen in Fig. S5, the error between our NNP predictions and reference hybrid DFT calculations remain below 5% all along the reaction path. The NNP are therefore suitable to study the self-dissociation reaction in the different environments considered.

We can therefore conclude that our simulation method yielded statistically converged S8 results of hybrid DFT quality and with explicit accounting of nuclear quantum effects.

Free Energy Profile Calculations

All free energy profiles presented were obtained via Umbrella Sampling (US) calculations with a homemade weighted histogram analysis method (WHAM) 20 implementation. The errors correspond to 95% confidence intervals, which were obtained with a procedure described by Zhu and Hummer. 21 Two-dimensional free energy surfaces were also obtained with Umbrella Sampling and the WHAM implementation of Alan Grossfield. 22 The bias potentials required for the US procedure were included in the MD simulations with the PLUMED code 23 interfaced with either i-PI or LAMMPS.

Reaction Free Energy Profiles

The free energy of water self-dissociation was studied with a collective variable describing the coordination number n of all hydrogen atoms of the system around the reactive oxygen atom (O A ):

n = X j2{H} 1 ⇣ r Aj r 0 ⌘ 12 1 ⇣ r Aj r 0 ⌘ 24 (S1) S10
where the cut-off radius r 0 was set to 1.3 Å. With this collective variable, the reactant state was found to lie at n ⇡ 1.9 (undissociated water) and the product (deprotonated water) state corresponds to a flat basin below n ⇡ 1.3.

For systems with quantum nuclei, the reaction US simulations were performed with a total number of 53 window simulations of 40 ps where the first 10 ps were used as equilibration.

Since classical nuclei MD were significantly more efficient, the corresponding US calculations were performed with 53 window simulations of 70 ps, where the first 20 ps were used as equilibration. Since no significant changes in the errors of the WHAM procedure (or in the shape of the free energy profiles) were observed, we concluded that these simulations were long enough to converge the free energy in every given window. The target coordination number of every window was obtained by applying a harmonic potential with a force constant of 20 eV per coordination unit squared and, for simulations in a slab system, an additional harmonic restraint on the position of the center of mass of O A along the z direction was added, with a force constant of 9.72 eV/Å 2 . Thus, the simulation windows were separated by 0.025 units of coordination and spanned the [1.0, 2.3] interval.

We note that with such a coordination number interval we deliberately chose to deprotonate a water molecule (of fixed identity). While the opposite choice (protonating a molecule of fixed identity) would be equivalent in the bulk, it could give different results in an interfacial system (where the z position of one reactive water molecule is fixed but the other diffuses in an anisotropic environment). We studied these effects with classical nuclei MD by adding 28 additional window simulations to the US calculation. The window simulations were of same length and had target n values in the [2.3, 3.0] interval (see the results below).

Solvation Free Energy Profiles

Solvation free energy profiles were calculated by biasing along the distance d com between the center of mass of the water slab and the oxygen atom of the solute of interest (hydronium cation, hydroxide anion or water molecule). For the US with NQE we used a total number of S11 interval (uniformly separated by 0.125 Å). For classical nuclei MD, we extended the US with 8 additional windows in the [11.0, 12.0] Å interval. Target distances were enforced with a harmonic constraint of 0.5 eV/Å 2 . In order to prevent reactive events during these simulations, harmonic walls were added on the coordination number of the oxygen atom of interest, with a force constant of 20 eV per coordination unit squared. These acted when below n = 2.7 for the hydronium, above n = 1.3 for the hydroxide, and above and below n = 2.3 and n = 1.7 respectively for the water molecule.

The choice of the distance to the center of mass d com as a collective variable allowed efficient biasing during the US simulations, but it is not immediately indicative of the solvation state of the solute and its values are heavily dependent on the simulation slab thickness.

Therefore, we chose to represent the free energy as a function of the distance to the air-water interface (or Gibbs dividing surface), defined as the (xy) plane where water density reaches half of its bulk value. This choice of collective variable has been reported to yield very similar results to the more complex instantaneous interface 24 for the solvation free energy of the hydronium ion, 25 and it has the advantage of being independent of the chemical nature of the solute (whereas deformations of the instantaneous interface can strongly change between different solutes, especially when those are charged). Furthermore, the free energy profiles along this distance to the interface s can be straightforwardly obtained by translating the local frame of reference along the z axis. For the slab systems with 128 water molecules, the interfaces were found to be 7.55 Å above and below the center of mass of the slab, so that s = d com 7.55 Å.

Evolution of the Solvation Environment along the Distance to the

Air-Water Interface

Free energy calculations show that the hydronium cation is stabilized by ⇡ 1.2 kcal•mol 1 at the air-water interface, whereas the hydroxide anion is destabilized by ⇡ +2.7 kcal•mol 1 S12 (see Fig. 2b of the main text). The solvation free energy profile of a water molecule is almost flat when going from bulk to the interface and increases only when going to the vapor phase.

This does not mean that the water molecule is indifferent to the air-water interface, but rather that placing a water molecule at the interface in our biased simulation corresponds on average to replacing another identical molecule, which was originally at the interface.

Therefore, the solvation free energy profiles presented here can be understood as indicating if a solute is preferred at the interface when compared to the water molecules that it would replace, and effectively indicate the propensity of the species within the slab.

In order to understand the relative interfacial stabilities of both self-ions compared to water we studied the solvation structure of each species along the distance to the air-water interface (see Fig. S6 to S8). First, when the hydronium cation is moved from the middle of the slab to the air-water interface, there is no significant change in the positions of the peaks of the radial distribution functions around the ion oxygen atom (Fig. S6a andS7a).

Furthermore, only a small decrease in the coordination number around the ion oxygen atom at the position of its first solvation shell is observed (going from ⇡ 3.4 to 3.2). These results suggest that the solvation environment of H 3 O + changes only slightly at the interface, with the only remarkable difference being the loss of a weak, partially accepted hydrogen bond.

This fourth accepted hydrogen bond is known to destabilize the positively charged hydronium cation and thus is only transiently observed in bulk conditions. Therefore, the loss of this partial, destabilizing hydrogen bond, as well as the gain of almost one full (stabilizing) hydrogen bond for the water molecule when going from the interface to the bulk (Fig. S6c and S7c), explains the interfacial stabilization of the hydronium cation.

The changes in the first solvation shell of the hydroxide anion are more drastic. As can be seen in Fig. S6b and S7b, the coordination number of HO after the first peak of the radial distribution functions decreases by a full unit when going from bulk to interface. A close inspection of the radial distribution function of the hydroxide hydrogen atom with the water oxygen atoms (Fig. S8a) shows that upon desolvation, the hydroxide anion loses a S13 weak donated hydrogen bond. However, the loss of this weak hydrogen bond is accompanied by only a moderate decrease of the coordination number (Fig. S8b), meaning that there is as well a partial loss of an accepted hydrogen bond by the hydroxide oxygen atom. These findings are in good agreement with other ab initio MD results. 26 The loss of both a weak donated hydrogen bond and the weakening of a strong accepted hydrogen bond (or a moderate weakening of all accepted hydrogen bonds) at the interface could explain the unstable character of the hydroxide anion at the air-water interface. To complement this picture, we also studied the change in the distribution of the OH orientations of the hydroxide anion, as a decrease of orientational entropy has been proposed in the literature as the main cause for the hydroxide destabilization at the interface. 27,28 The results (Fig. S9) show that the orientation of the OH bond of the hydroxide anion at the air-water interface is indeed more constrained than in a bulk-like environment, with a preferential orientation of the OH bond towards the vapor phase. This induces an increased entropic cost of ⇡ 0.7 kcal•mol 1 at 300 K for the desolvation process. While important (⇡ k B T ), entropy cannot therefore account for the full interfacial destabilization of the hydroxide anion, which indicates that the loss of solvation plays a more important role than previously thought. We note that although water shows a similar desolvation at the interface (Fig. S6c and Fig. S7c), the associated free-energy cost should be lower than for the charged hydroxide anion.

S17

Comparison with Classical Nuclei MD Results

In order to estimate the impact of NQE, of symmetry breaking in the reaction coordinate, and of the simulation system size on the principal results of this study, we performed a set of complementary simulations with classical nuclei MD. We obtained similar reaction and solvation free energy profiles, and verified the validity of the local acidity model (eq. 1 of the main text). We also extended the dissociation free energy profiles and the acidity model to increasing values of the reaction coordinate (the target water molecule is protonated by a free molecule) and showed that this leads to identical conclusions about the role of solvation free energy for the change in water acidity at the interface. Finally, we performed the same type of calculations with a system of 256 molecules to show that our conclusions are not due to simulation cell or system size artifacts. We now describe these results.

Reaction and Hydration Free Energy Profiles

We computed free energy profiles for both the solvation of water and its self-ions, and for the self-dissociation of water in bulk and at different depths relative to the air-water interface with classical nuclei MD. These free energy profiles (Fig. S10) lead to the same qualitative conclusions as those obtained with path-integral MD (presented in the main text) concerning the relative stabilities of all species at the air-water interface. Furthermore, they can be used in the same manner (eq. 1 to 6 of the main text) to obtain continuous variations of the reaction free energy along the distance to the interface (Fig. S11), global pH values (Fig. S12), and spatially resolved acid/base concentration profiles. The agreement between model and reactive simulations remains excellent despite the absence of NQE, which induce a significant change in the absolute free energy values obtained for both solvation and reaction.

This shows that our model is robust with respect to the simulation method employed for determining the solvation free energies and it could therefore be used to determine acid-base reactivity of any given species from inexpensive classical MD (provided that an adequate S18 forcefield is available). The global pH values (Fig. S12) remain above bulk for all system sizes and the effect of hydronium solvation is still largely dominant compared to changes in total concentration, meaning that the conclusions of this study are also robust with respect to the simulation method employed. Figure S12: Total (black curves) and decomposed (green and purple curves, see eq. 4) pH values of systems of increasing typical size with respect to bulk pH. Plain (dashed) curves represent systems with planar (spherical) interfaces, ie. slabs (droplets). Typical sizes of slab and droplet systems correspond to the slab half-thickness and droplet radius respectively, obtained from the Gibbs dividing surface (GDS). Results are obtained with classical nuclei MD, errorbars are not presented for clarity, but are comparable to those shown in Fig. 3 of the main text.

nuclei and path-integral MD. Moreover, contrary to what could be expected, the explicit inclusion of NQE further reduces the self-dissociation free energy, which suggests that the reference DFT method that we use (revPBE0-D3) is not ideal for the study of the acidbase reactivity of water. This result is all the more surprising since revPBE0-D3 is known to be one of the best hybrid functionals for the structure and dynamics of (undissociated) liquid water. 29,30 Elucidating the origins of the shortcomings of the revPBE0-D3 potential for dissociating OH bonds when increased zero-point energies are included via path-integral MD is beyond the scope of this work. Nonetheless, we remain confident in the adequacy of the solvation-based model proposed here to study acid-base reactivity and on the conclusions that we have so far obtained from it, as none of these depend on the specificities of either the density functional or the type of MD simulation.

Symmetric Choice: Fixed Hydronium and Diffusing Hydroxide

As we stated above, by studying the deprotonation of a singularized water molecule (n < 2)

we broke the natural symmetry of the self-dissociation process. This choice is completely ir-S20 relevant in the bulk environment, since all water molecules are equivalent, but it becomes important at the air-water interface. Indeed, the simulation protocol employed at the interface needs to restrain the depth s of one water molecule, and thus of one of the self-dissociation products after reaction. Since neither reactants (water) nor products (H 3 O + and HO ) are indifferent to the air-water interface (see solvation structures above), the choice of protonating or deprotonating a water molecule at a fixed given depth will have an impact on the reaction free energy obtained in the simulation. The local acidity model presented in the main text (eq. 1) can be trivially adapted to the case where the water molecule that is kept fixed is protonated (see below for a more detailed discussion of the model):

r G 0 (s) = r G G H 2 O hyd (s) h G H 2 O hyd i + h G HO hyd i + G H 3 O + hyd (s) (S2)
where we employed the same notations as in eq. 1 in the main text. We note that formally all the free energies computed in our simulations correspond to canonical Helmholtz free energies, but the G notation is used for simplicity.

In order to verify the validity of eq. S2, and thus that acid-base reactivity at the interface can indeed be simply explained by the solvation environment of the reagents, we extended the dissociation free energy profiles to larger values of the protonation collective variable and compared the simulation results with the predictions of the model.

As can be seen in Fig. S13, the choice of protonating or deprotonating the fixed water molecule has indeed a drastic effect on the variations of the reaction free energy near the air-water interface. Indeed, while bulk reaction free energies obtained with either simulation choice are equal within the errorbars (0.2 kcal•mol 1 ), the variations of reaction free energy when the hydronium cation is left at a fixed interfacial depth (n ⇡ 2.8) are less pronounced than in the opposite situation (fixed hydroxide anion, n < 1.3). This is due to the overall destabilization of the hydroxide anion in the slab system being compensated by the stabilization of the hydronium cation near the interface, and to the latter being significantly less S21 a) b)

Figure S13: a) Extension of the water self-dissociation free energy profiles from Fig. S10a to larger values of the coordination number. b) Predictions of the local acidity model for the continuous variations of the reaction free energy when the fixed water molecule is either deprotonated (leaving a fixed hydroxide anion, eq. 1) or protonated (fixed hydronium cation, eq. S2). Plain lines and colored dots represent the model predictions and reactive simulation results respectively, the thick colored lines represent 95% confidence intervals.

pronounced (⇡ 1.0 kcal•mol 1 with classical nuclei MD) than the destabilization of the hydroxide at the same position (⇡ +2.0 kcal•mol 1 with classical nuclei MD). The model predictions still agree with the US simulation results for the three interfacial depths considered (Fig. S13b) and, as will be shown below, the choice between using r G 0 (s) (eq. S2) or r G(s) (eq. 1) has no impact on the resulting pH and hydronium concentrations.

Simulation System Size Effects

In order to examine possible artifacts due to the finite size of our simulation system we repeated the same type of MD simulations with classical nuclei and a system that is twice larger (256 water molecules, see Fig S1c,d).

We first computed the self-dissociation free energy in a bulk-like cubic cell of side 19.76 Å, following the same US procedure as for the smaller system (described above). We then proceeded to compute the self-dissociation free energies at two different interfacial depths in a 256 water molecules slab system, where both the thickness and surface had been enlarged (21⇥21⇥63 Å 3 ). We also computed the solvation free energy profiles of water and its self-ions S22 in this system. This allows to detect effects coming from both incomplete bulk-like solvation due to a too thin slab and from difficult surface deformation induced by the increased surface tension associated with small interfacial area (only small capillary wavelengths are allowed due to the periodic boundary conditions along the x and y directions).

The self-dissociation and solvation free-energies obtained are presented in Fig. S14, and it is easy to see the great resemblance with those obtained in the 128 molecules system (Fig. S10). Solvation free energy differences between the middle of the two slabs and the corresponding interfaces are in perfect agreement with each other, proving that there are no significant size artifacts. Furthermore, we used the solvation free energy profiles in Fig. S14b and the reference bulk self-dissociation free energy in Fig. S14a in conjunction with our local acidity model (eq. 1), and compared the results to those obtained with reactive US simulations at the interface. As can be seen in Fig. S15, the agreement remains excellent. This is another indication of the robustness of the model proposed here, which can be used to study a broad range of system sizes.

S24

Water Self-Dissociation Mechanism

The self-dissociation mechanism of water has recently attracted great attention both in the bulk 31 and at solid-liquid interfaces, 32,33 due to its fundamental importance for the understanding of acid-base chemistry and to its implications for the energy sector. A recent study has shown the importance of the role played by the solvent in the separation of water selfions in the bulk, 31 making the explicit consideration of ion-ion distances essential for the correct convergence of self-dissociation free energies. Here our NNP-based simulations are long enough to converge one-dimensional dissociation free energy profiles along the coordination number collective variable (eq. S1) without explicitly biasing this separation coordinate.

However, the projection of these free energies along a second collective variable provides additional information concerning the reaction mechanism and its potential modifications at the air-water interface. We therefore employed a very similar collective variable to that used in the literature 31 to monitor the separation of water self-ions in solution:

R ion = 1 (n / n c ) pcut 1 (n / n c ) qcut r ion (n; n c ) (S3)
where r ion (n) gives the distance between both self-ions if the coordination number n is below the cut-off value n c = 1.56, and the average water-water distance (2.8 Å) otherwise. The parameters p cut and q cut were set to 12 and 24 respectively. Eq. S3 corresponds to a simplified version of the R ion collective variable used in the literature. 31 Since the reactive simulations are long enough that no biasing is needed to converge the R ion variable, this expression based on the explicit identification of both self-ions' positions is less computationally cumbersome for post-processing purposes than the computation of all coordination numbers required by the traditional version. The two-dimensional free energy surfaces resulting from the projection of the one-dimensional free energy profiles along this variable are presented in Fig. S16 and S17 for path-integral and classical nuclei MD respectively.

Dissociation in bulk proceeds almost sequentially, first the coordination number of the S25 reactive water molecule is lowered (a proton is transferred to a molecule in the first solvation shell), and then the self-ion pair can diffuse and separate in a rather flat free energy landscape.

This mechanism is consistent with what has been observed in the literature and suggests that the only product state of the reaction is a solvent-separated ion pair. The contact ion pair state, frequently observed as a local free energy minimum in salt dissociation, is not observed here, probably due to the very favorable character of the ion-ion recombination reaction. When the reaction is undertaken in a system presenting an air-water interface, the diffusion profile of the ions (surface along R ion at fixed coordination number) becomes slightly rougher, with the apparition of a true barrier at R ion ⇡ 3 Å. A more complex diffusion landscape for both self-ions is to be expected in an interfacial system, due to the rapidly changing solvation environment (Fig. 2 of main text and S10b). Despite this change in the shape of the barrier and in the relative free energies of reactant and products, the selfdissociation mechanism itself does not change (deprotonation followed by diffusion leading to a solvent-separated ion pair), suggesting that self-dissociation thermodynamics are governed by relative solvation free energies in the reactant and product states. Nonetheless, despite the moderate changes observed for the free energy barrier height, the changes in barrier shape might be accompanied by a change in friction, which could greatly impact kinetics.

Although interesting, a detailed study of the reaction kinetics is beyond the scope of this work.

S26

Figure S16: Two-dimensional reaction free energy surfaces obtained with a second collective variable R ion (eq. S3) measuring the ion-ion distance in solution for the bulk and interface dissociations studied with path-integral MD.

S27

Figure S17: Two-dimensional reaction free energy surfaces obtained with a second collective variable R ion (eq. S3) measuring the ion-ion distance in solution for the bulk and interface dissociations studied with classical nuclei MD.

S28

Size Dependence of pH and Self-Ion Distributions: eqs 1,

4-6

In this work we use a previously developed model 34 relating the acidity of a given species across the air-water interface to the solvation free energies of all the reagents. An important result was the explicit dependence of the acid dissociation fraction on the total interfacial system size. Here, we extend this model to the case of neat water, where the dissociating species and the solvent are the same. We use it to obtain a new analytical expression of the pH value of an interfacial system, which explicitly depends on the total system size. Since reactivity in such a heterogeneous medium cannot be described solely by this average value, we derive the spatial distributions of water self-ions near the interface and characterize how they depend on system size. Finally, the counter-intuitive increase of hydronium concentration at the interface with system size is explained by the interplay between surface saturation and size of the hydronium reservoir. Below we present the derivations for all these results.

Local Acidity Model for Water Self-Dissociation

The local acidity model (eq. 1 of the main text) provides a representation of the reaction free energy of an acid-base reaction as a function of the solvation free energies (both local and averaged) of all reagents and of the bulk reaction free energy. It was originally derived for reactions of acidic species with water 34 and here we modify it slightly for the special case of water self-dissociation. The main assumption of the model is the non-interaction of the different reagents in solution (which is exact for infinitely diluted solutions) so that the reaction can be decomposed into the simple thermodynamic cycle of Fig. 2a, graphically summarized in Fig. S18. In our US simulations at the interface, both one of the reactant water molecules and the hydroxide anion have a constrained depth s within the slab, whereas the other water molecule and the hydronium cation are free to diffuse. The model postulates that the only difference with the bulk reaction comes from the change in the solvation S29 environment of all reagents, so that the reaction free energy at any given position within the slab can be obtained as given in eq. 1.

Figure S18: Schematic depiction of the thermodynamic cycle underlying the local acidity model (eq. 1, reproduced in the figure for the case of a system in slab geometry). The self-dissociation reaction of a water molecule at fixed depth s in a slab of thickness L can be decomposed as follows: first the reactants (water molecules) are brought from the slab system to a bulk environment (one is brought from depth s, the other from its average state within the slab), the reaction proceeds in bulk under standard conditions, and finally the products are again brought to the slab system (the hydroxide is brought to position s, the hydronium is left free within the slab). As discussed in the text, this can easily be adapted to different system geometries.

The average solvation free energy changes (with respect to bulk) involved in the model are obtained from the ensemble of local free energies. For any given species i in an arbitrary geometry, this yields:

⌦ G i hyd ↵ = k B T ln  1 ⇤ Z ds (s)e G i hyd (s)/k B T + k B T ln  1 ⇤ Z ds (s) (S4)
with ⇤ the thermal De Broglie wavelength, s the distance to the interface and (s)ds the volume element of the system geometry (for a slab parallel to the xy plane s = z and (s) = L x Ly is the surface area of the slab, and for a spherical droplet s = r and (s) = 4⇡s 2 ). In the limit of an infinite system, ⌦ G i hyd ↵ = 0 because the interfacial region where G i hyd (s) 6 = 0 becomes negligible compared to the bulk-like region in the core of the system.

The local solvation free energy difference between any given distance to the interface s and the bulk environment G i hyd (s) is approximated by the difference between the free S30 energy at the center of the finite-size system under consideration and that at distance s from the interface in the same system. When systems larger than the ones used in our simulations were considered (Fig. 3-5 in the main text), we extended the solvation free energies from our MD simulations by zero-padding in the middle of the simulation slab. Both of these approximations rely on the assumption that the solvation free energies are already converged with respect to bulk in the middle of our simulation slabs (this was thoroughly studied in our previous work and is also confirmed by the simulations performed with the slab of 256 water molecules).

It is important to note that the only difference in the hydration free energies obtained in systems with different geometries comes from the difference in volume elements (s)ds, i.e.

it is an entropic term dependent on the difference in available space for a given distance to the interface. This implicitly neglects more subtle effects of the interface geometry on the solvation of acid-base species, such as e.g. curvature. For the effects of interface curvature on the solvation of ions the validity of this hypothesis was verified independently and will be the object of a later publication. The impact of the Laplace pressure on the water selfdissociation equilibrium is neglected due to the extremely low pressure-dependence of the water self-ionization product. 35 pH and local concentration profiles of interfacial systems of varying sizes

The pH value of a system, defined in terms of the activity of the hydronium cation, is the quantity that determines the acid-base composition at thermodynamic equilibrium. As discussed in the main text, it can be written in terms of the hydronium ion concentration and activity coefficient:

pH = log 10 e (µ H 3 O + µ H 3 O + ) = log 10 ✓ H 3 O + (s) [H 3 O + ](s) c ◆ (S5) S31
where is the inverse temperature (k B T ) 1 and c = 1 mol/L is the standard state concentration.

In our simulations, the chemical potential is not fixed, but an average equilibrium value can nevertheless be obtained, as will be shown below. We will assume that ion concentrations are low enough that all deviations from ideality arise from the distinct solvation properties of the interfacial system, and not from ion-ion correlations, which is legitimate for a pure airwater interface and makes the derivation simpler. As stated in the main text, the hydronium chemical potential at equilibrium can then be obtained under the above assumptions from the solvation thermodynamics perspective:

36 µ H 3 O + = µ H 3 O + + G H 3 O + hyd (s) + 1 ln [H 3 O + ] (s) c (3)
which implies that

H 3 O + (s) = exp h + G H 3 O + hyd (s)
i . Here G H 3 O + hyd (s) is the hydration free energy change of the hydronium cation at distance s with respect to bulk standard conditions, which determines the relative ion concentration profile at thermal equilibrium (Maxwell-Boltzmann distribution):

[H 3 O + ] (s) h[H 3 O + ]i = e G H 3 O + hyd (s) 1 R ds (s) R ds (s)e G H 3 O + hyd (s) (S6)
where (s)ds corresponds to the volume element of the system geometry and h[H 3 O + ]i is the average hydronium concentration in the system. This quantity will be further defined below, but simply using the definition of the average hydration free energy eq. S4 already provides the simpler expression presented in the main text (eq. 6) for the local concentration profile of hydronium cations:

[H 3 O + ](s) = ⌦ [H 3 O + ] ↵ exp G H 3 O + hyd (s) h G H 3 O + hyd i k B T ! (6)
Combining this expression with that of the solvation thermodynamics chemical potential S32 eq. 3 makes apparent that this quantity is indeed uniform at chemical equilibrium: We note that a similar expression could be easily written for other acid-base equilibria, where the dissociation probability should be equally weighted by the acid density distribution within the system. However, while this formula is easily computable from the local acidity model eq. 1, it breaks the natural symmetry of the water-self dissociation process, since in the expression of r G obtained from eq. 1 the hydroxide anion is fixed while the hydronium cation can freely diffuse. In order to obtain the more symmetric eq. 5 presented in the main text we introduce the link between water density and free energy at the interface

µ H 3 O + = µ H 3 O + + G H 3 O + hyd (s
⇢ H 2 O (s)/⇢ H 2 O = exp[ G H 2 O
hyd (s)] into eq. S15 and do a few simple manipulations with the definitions of the local reaction free energy eq. 1 and the average hydration free energies S34 eq. S4 : Finally, we can introduce the fully averaged reaction free energy (average of eq. 1)

⌦⇥ H 3 O + ⇤↵ = c s h⇢ H2O i ⇢ H 2 O R ds (
h r Gi = r G ⌦ G H 2 O hyd ↵ ⌦ G H 2 O hyd ↵ + D G H 3 O + hyd E + D G HO hyd E (S20)
to obtain eq. 5 from the main text:

⌦⇥ H 3 O + ⇤↵ = c s h⇢ H2O i ⇢ H 2 O e h r Gi (5)
This quantity was used to obtain Fig. 3-5 in the main text, and also allows a direct comparison between the hydronium concentration and activity coefficient profiles, which is presented in Fig. S19.

Accumulation of Hydronium Cations at the Interface

Although one could have expected the concentration of hydronium cations at the interface to increase in systems with large surface/volume ratios, as the relative impact of the interface is larger, our results (Fig. 4 and 5 of the main text) show the opposite size-dependence.

S35

Figure S19: Hydronium cation concentration and activity coefficient as a function of the distance s to the air-water interface in our simulation slab system of thickness L = 15.1 Å. It is apparent that local increases in the ion concentration are exactly compensated by a decrease in the activity coefficient.

We now explain how this increase in the accumulation of hydronium cations at the interface with system size results from the interplay between the surfactant character of the hydronium cation and the finite size of the system.

We consider a model two-state interfacial system of typical size L, consisting of a bulk-like core region and an interface layer of fixed thickness and surface area A. The total number of hydronium cations in such a system is simply the product of system volume V = AL and total ion concentration:

N H 3 O + (L) = ⌦⇥ H 3 O + ⇤↵ AL (S21)
This total number corresponds to hydronium cations both in the core and interface regions of the system, so that we can write the total number of hydroniums at the interface as (assuming thermal equilibrium):

N int H 3 O + (L) = N H 3 O + (L) A e G int A e G int + e Gcore (AL A ) (S22) = ⌦⇥ H 3 O + ⇤↵ AL 1 1 + e + G hyd (L/ 1) (S23)
where G core and G int correspond to the hydronium solvation free energies in the core and S36 interface regions respectively and G hyd = G int G core is the hydronium hydration free energy change from bulk to interface conditions.

We can then obtain the interface concentration by dividing by the interface volume

V int = A : ⌦⇥ H 3 O + ⇤↵ int = ⌦⇥ H 3 O + ⇤↵ L 1 + e + G hyd L 1 (S24)
which can be straightforwardly adapted to spherical geometry to describe droplet systems:

⌦⇥ H 3 O + ⇤↵ int = ⌦⇥ H 3 O + ⇤↵ R 3 1 + e + G hyd R 3 1 (S25)
As can be seen in Fig. S20, eqs. S24 and S25 do capture the qualitative variations of the average interface hydronium concentration with total system size. It is therefore the combination of the interfacial stabilization ( G hyd < 0) of the hydronium cation with the closed nature of the small systems under consideration (fixed total number of cations) that produces the observed increase in interface concentration with system size.

Figure 2 :

 2 Figure 2: (a) Thermodynamic cycle used to model the water dissociation reaction free energy within the slab eq. 1. The upper row describes the dissociation in the slab, where a water molecule at a fixed depth s dissociates and yields a hydroxide ion at the same depth and a hydronium ion which is free to diffuse. The lower row shows the water dissociation reaction in the standard state. The hi and s subscript notations indicate whether a molecule is free to diffuse or at a fixed depth relative to the interface. (b) Hydration free energy profiles along the distance s to the Gibbs dividing interface for water and hydronium and hydroxide ions for a 128 H 2 O slab system of average thickness L = 15.1 Å; free energies are referenced with respect to their values in the middle of the slab, at s = L/2. (c) Dissociation free energy as a continuous function of the distance to the interface given by the local acidity model eq. 1 (blue line) and obtained from Fig. 1 (black dots). The bulk reference value is shown by the horizontal line and the insets show typical configurations at the corresponding dissociation depths. For both figures the 95% confidence intervals are represented by colored regions.

Figure 3 :

 3 Figure 3: pH shift relative to the bulk value for (a) slabs and (b) spherical droplets of decreasing sizes. The pH shift (black) obtained from eq. 4 and its contributions due to hydronium concentration (green) and hydronium stabilization (purple) are shown, together with the shift in pKw=pH+pOH (dashes). The slab thickness and droplet radius are determined from the Gibbs dividing surface. (c) and (d) idem for pOH.

Figure 4 :

 4 Figure4: (a) Logarithmic measure of self-ion concentrations relative to their bulk value as a function of the distance to the air-water Gibbs dividing interface for a range of systems with surface/volume (S/V) ratios going from 0 (macroscopic air-water interface system, plain line) to 3/7.65 Å 1 (spherical droplet of 7.65 Å radius equal to our simulation slab half-thickness, dashed lines). (b) Logarithmic measure of the acid/base concentration ratio as a function of the distance to the air-water interface for the same range of systems.

Figure 5 :

 5 Figure 5: Change in the local concentration of water self-ions relative to bulk in (a) the core and (b) at the Gibbs dividing surface of systems of increasing size, represented by the surface-to-volume ratio (lower axis) and the radius of a spherical droplet (upper axis). Errors in panel (b) are not represented for readability due to the logarithmic scale, but are comparable to those in panel (a).

5 to 2 •

 2 10 6 depending on the size of the training set (see Iterative Construction of the Training Set below) and with a learning rate that exponentially decays from 1•10 3 to 1•10 8 during the iterative procedure.

  Our initial training set consisted in 2188 configurations of a single hydronium ion in a bulk cubic box of 64 water molecules and 12.42 Å side length, generated with this same procedure for a previous work (in publishing). The final training set contains 11,148 configurations, separated into 9 different systems. The training set composition and the different systems are presented in Fig. S2. We note that, due to the iterative nature of the training set construction, there are more configurations with a hydronium ion than for the other chemical compositions but the NNP yield results of similar accuracy for every system under consideration (see next section).

Figure S2 :

 S2 Figure S2: Number of input configurations in the training set per system considered. Systems can either contain a hydronium cation (H 3 O + ), a hydroxide anion (HO ) or pure water (H 2 O). Configurations can be either bulk or interface, and contain either 64 or 128 water molecules in cubic boxes of either 12.42 Å or 16 Å side-lengths respectively for bulk simulations and in a 16⇥16⇥48 Å 3 box for the interface.

Figure

  Figure S3: a) Root mean square (RMSE) and mean absolute (MAE) errors on 10 different validation data sets made of 128 water molecules in either bulk or interface conditions, potentially containing a hydronium cation or hydroxide anion (indicated above the bars). Hatched / Plain colored bars indicate configurations generated with path-integral / classical nuclei MD. The dashed horizontal black line indicates the 0.1 eV/Å threshold used to select candidate configurations. b) Correlation density plot between the atomic forces predicted by the NNP and those obtained with the reference DFT method over the entire 1,000 configurations of the validation data set.

Figure S4 :Figure S5 :

 S4S5 Figure S4: Relative prediction error (root mean square error relative to standard deviation) obtained with NNP trained on different numbers of configurations randomly drawn from the total training set. Plain / dashed orange lines indicate the errors on the path-integral / classical nuclei configurations of the validation set. Constant error indicates that the NNP description is converged with respect to training set size, even for small amounts of training samples (up to 5% of total training set obtained with the iterative procedure), showing the importance of having dispersed data over the chemical space.

Figure S6 :Figure S7 :Figure

 S6S7 Figure S6: Oxygen-oxygen radial distribution functions (g OO ) and average coordination numbers (n OO ) for a series of distances to the air-water interface for the oxygen atom of a) a hydronium cation, b) a hydroxide anion, and c) a water molecule.

Figure

  Figure S9: a) Orientation probability density distributions of the hydroxide OH vector (see inset) for a series of distances to the air-water interface. b) Contribution of the orientational entropy of HO (S ori. , obtained from the probability density ⇢(cos ✓) as shown in the inset) to the total solvation free energy for a series of distances to the air-water interface.

Figure

  Figure S10: a) Self-dissociation free energy profiles in bulk and at different distances (d) to the air-water interface of the slab system containing 128 water molecules. b) Solvation free energy profiles of water and its self-ions along the distance to the air-water interface. In both cases negative/positive distances indicate the liquid/vapor phase and thick colored lines represent the 95% confidence intervals. Results were obtained with classical nuclei MD.

Figure S11 :

 S11 Figure S11: Water self-dissociation free energy as a function of the distance to the interface given by the local acidity model (blue line) and obtained from the reactive US simulations of Fig. S10a (black dots), both with classical nuclei MD. The bulk reference is represented by a horizontal black line and the insets show typical configurations at the corresponding dissociation depths

FigureS23Figure S15 :

 S15 FigureS14: a) Self-dissociation free energy profiles in bulk and at different distances (d) to the air-water interface of the slab system containing 256 water molecules. b) Solvation free energy profiles of water and its self-ions along the distance to the air-water interface. In both cases negative/positive distances indicate the liquid/vapor phase and thick colored lines represent the 95% confidence intervals.

Figure S20 :

 S20 Figure S20: Hydronium cation concentration at the surface of systems of increasing size in either a) slab or b) droplet geometries as obtained from the pH model eq. S6 (blue dashed line) and from eqs. S24 and S25 (orange plain lines). The interface thickness used was set to 1 Å.

  All that is left to compute the global pH (eq. 4) and the ion concentration profile (eq. 6) is obtaining the average hydronium concentration h[H 3 O + ]i. We start from a local expression of the reaction free energy measured in our simulations, akin to a local version of the law of mass action where we further assimilate the activity of water in our system to that of a⇢ H 2 O (s) h[H 3 O + ]i ⇢ H 2 O e H 2 O (s)is the density of water at depth s (for large values of s we get the standard state density ⇢ H 2 O so that we recover an activity of 1 for the solvent).An instructive reformulation of this expression of the hydronium concentration can be obtained by introducing an average of the reaction free energy over the water density

	from which we can obtain the global hydronium concentration:
		) + E + H 3 O + ⇤↵ 2 = c 2 = µ H 3 O + + D G H 3 O + hyd ⌦⇥ ⇢ H 2 O 1 ln e 1 ln h[H 3 O + ]i G H 3 O + hyd c Z 1 R ds (s) = c 2 h⇢ H2O i ⇢ H 2 O R ds (s)⇢ H 2 O (s)e (s) + 1 ln ds (s)⇢ H 2 O (s)e h[H 3 O + ]i e c + r G(s) ⌧ r G(s) R ds (s)⇢ H 2 O (s)	G	H 3 O + hyd	(S7) (S8) (S12) (S13)
	and immediately gives eq. 4 from the main text when combined with the definition eq. 2:
			pH = log 10 ⌦ e r G ↵ H 2 O = ✓	h[H 3 O + ]i c R ds (s)⇢ H 2 O (s)e ◆ D G H 3 O + hyd k B T ln 10 E r G(s) R ds (s)⇢ H 2 O (s)	(S14)
	yielding a global hydronium concentration that is thus obtained by a careful average of the
	dissociation probability (linked to exp[
	solvent:		⌦⇥	H 3 O + ⇤↵	= c	s	h⇢ H2O i ⇢ H 2 O	he	r G i H 2 O	(S15)
				c e	r G(s) =		[HO ] (s) h[H 3 O + ]i c ⇢ H 2 O (s)/⇢ H 2 O	(S9)
		[HO ](s) = c 2 From this last expression combined with electroneutrality we can write the total number
	of ions as:
	Z	ds (s)[HO ](s) =	Z	ds (s)[H 3 O + ](s) =	c 2 h[H 3 O + ]i ⇢ H 2 O	Z	ds (s)⇢ H 2 O (s)e	r G(s) (S11)
									S33

r G(s) (S10)

where ⇢ r G(s)]) weighted by the density of the dissociating species (water in this case):

  s)⇢ H 2 O (s)e

	= c	v u u u t h⇢ H2O i ⇢ H 2 O	R	r G(s) hyd (s) e ds (s)⇢ H 2 O (s) ds (s)e R G H 2 O ✓ R	r G ds (s)e G H 2 O hyd (s) h G H 2 O hyd i+ G HO hyd (s)+ G H 2 O hyd (s)	⌧	(S16) G H 3 O + hyd	◆
										(S17)
	= c = c	v u u u t h⇢ H2O i ⇢ H 2 O s h⇢ H2O i ⇢ H 2 O	e e	✓ ✓	r G r G	h G H 2 O hyd i+ e h G H 2 O ⌧ G H 3 O + hyd hyd i R 2h G H 2 O hyd i+ ⌧ G H 3 O + hyd	◆ ds (s) e h G HO hyd i R +h G HO hyd i	ds (s)	(S18)

◆ (S19)

Acknowledgement

This work was supported by an HPC allocation from GENCI at IDRIS (Grant 2023-A0130707156).

Supporting Information Available

The following files are available free of charge. Simulation details, NNP training and validation, free energy profile calculations, characterization of dissociation mechanism, derivation of eqs. 4-6. The NNP file is available from the authors upon reasonable request.

Supplementary Information

How the acidity of water droplets and films is controlled by the air-water interface

Simulated Systems in the Bulk and at the Interface

Four simulation systems were considered in this work, which are presented in Fig. S1. Except when explicitly stated otherwise, all results were obtained with the bulk and slab systems containing 128 water molecules (Fig. S1a,b). The effect of the simulation system size on the principal results of this study was considered (see below) by performing similar calculations in bulk and slab systems containing 256 water molecules. In order to model the air-water interface, 3D-periodic slab systems were used (see Fig. S1b,d), in which approximately three quarters of the simulation cell were empty of particles along the z direction, as is customary.

MD with Classical Nuclei

MD simulations with classical nuclei and neural-network potentials (NNP) were performed with the LAMMPS software 2,3 interfaced with DeepMD-kit. 4 Simulations were carried out in the canonical (NVT) ensemble at a target temperature of 300 K. Constant temperature was S3 enforced through a Nosé-Hoover thermostat 5-7 with a time constant of 0.1 ps. The timestep was set to 0.25 fs for all simulations.

MD with Nuclear Quantum Effects

Nuclear Quantum Effects (NQE) were included into the MD simulations through the pathintegral formalism. We used the Path-Integral Generalized Langevin Equation Thermostat 8 implementation of the i-PI software 9 with a target temperature of 300 K. As in the simulations with classical nuclei, the timestep was 0.25 fs.