Like-charge anion pairing as a mechanism for ion co-localization in charged aqueous droplets

Abstract

Aqueous droplets serve as microreactors for chemical reactions in atmospheric and laboratory-generated aerosols. These droplets frequently carry excess charge, ranging from a small fraction of the Rayleigh limit to nearly the maximum charge sustainable before spontaneous fission. The effect of the droplet charge on chemical reactivity remains insufficiently understood. Here we identify like-charge contact ion-pairing of OH⁻ ions as a mechanism for transient co-localization of pre-reactive species. Using molecular dynamics simulations of charged water nanodroplets containing multiple OH⁻ ions, we quantify the equilibrium constant, lifetime, and formation rate of contact ion-pairs, with Cl⁻ and Na⁺ ions included for comparison. We find that OH⁻ and Cl⁻ exhibit a significantly higher propensity to form contact ion-pairs than Na⁺ ions. Hydroxide contact ion-pairs are stabilized by hydrogen-bonded water bridges forming transient [OH⁻(H₂O)_nOH⁻] structures (n = 1–3) with lifetimes of ≈18 ps. These configurations are enriched in the droplet subsurface region, where ion density is elevated and water self-diffusion is up to twofold faster than in the bulk-like interior. Although diffusion enhancement alone does not account for reported rate accelerations in droplets, it increases ion–ion encounter frequencies in regions where pairing is most probable. We suggest that rather than directly producing radicals, these contact ion-pairs act as transient, pre-reactive configurations that locally concentrate charge and restructure the hydrogen-bond network facilitating proton-transfer events within the hydrogen-bonded bridges and lowering barriers for electron detachment due to the presence of nearby charges.

---

1 Introduction

Aqueous droplets provide confined environments in which chemical processes can differ significantly from those in bulk solution. In particular, enhanced reaction rates^1,2 and altered product distributions^3–7 have been widely reported in microdroplet experiments across diverse systems. Despite extensive investigation, the molecular-level origins of these differences remain under active debate. Several mechanisms have been proposed, including interfacial effects, enhanced reactant confinement, strong electric fields, and modified solvation environments.^1,2,8–11

Here we show that contact ion-pairing of OH⁻ and Cl⁻ ions provides a general mechanism for pre-reactive co-localization in charged aqueous droplets. One experimentally important example where pre-reactive co-localization may be relevant is the reported formation of H₂O₂ in microdroplets. The mechanism of H₂O₂ formation has been the subject of considerable debate in recent years,^{3–7,10,12,13} but it remains unresolved. If H₂O₂ formation takes place via hydroxyl radicals, the formation requires not only the generation of hydroxyl radical intermediates, but also the effective detachment of the electron to avoid rapid recombination with the radical¹⁴ and finally the encounter of the hydroxyl radicals. Because of both the high reactivity of the radicals and the rarity of their formation, such encounters are unlikely unless precursor species are transiently co-localized. In the mechanism proposed by Colussi¹² co-localization has been implicitly considered for H₃O⁺ and OH⁻ species upon collision of droplets that carry excess H₃O⁺ ions with those with excess OH⁻ ions. Here we suggest that like-charge contact ion-pairing of OH⁻ provides an alternative and more general route to pre-reactive organization. Importantly, these OH⁻–OH⁻ contact ion-pairs are not assumed to undergo direct chemical transformation. Rather, they act as transient, correlated configurations that locally concentrate charge and reorganize the surrounding hydrogen-bond network. This co-localization may facilitate subsequent multi-body processes, such as interactions with charge carriers, interfacial electric-field fluctuations, lowering of ionization energy because of the presence of a neighboring negative charge or proton-transfer events, which may ultimately lead to the formation of reactive intermediates.

We study aqueous droplets primarily charged with OH⁻ ions, while Cl⁻ and Na⁺ ions are included for comparison. We examine the equilibrium constant between the contact ion-pair (CIP) and the “free” ions, the lifetime, and rate of formation of CIPs as well as the diffusion of the hydrated OH⁻ and the self-diffusion coefficient of water using “classical” molecular dynamics methods. The diffusion of OH⁻ in aqueous solution proceeds through both vehicular and structural transport mechanisms.^15–17 The structural diffusion occurs via proton transfer following the Grotthuss mechanism.¹⁸ Here, we perform simulations that do not allow for the breaking of chemical bonds, therefore, we only consider the vehicular diffusion. To our knowledge, this is the first study to quantify like-charge contact ion-pairing in multi-ion droplets and connect it to spatially heterogeneous diffusion and encounter kinetics.

While the solvation structure of individual alkali, halide, and hydroxide ions in bulk water has been extensively studied,^19–29 significantly fewer studies have examined the contact and solvent-separated ion-pairing of like-charged ions in confined or interfacial environments.^30,31 Existing work has largely focused on isolated ion-pairs, leaving open the question of how such interactions manifest in multi-ion, nanoscale droplets.

In addition to the thermodynamic stability of contact ion-pairs, the rate at which ions encounter one another is governed by molecular diffusion. While the self-diffusion coefficient of H₂O in bulk solution has been extensively characterized,^32–41 corresponding studies in droplets, where diffusion is expected to vary spatially, are still limited. Quantifying this position-dependent diffusion is essential for understanding structural diffusion for OH⁻ and H₃O⁺ ions and their effect in ion–ion encounter frequencies under nanoconfinement.

Although many experiments involve nearly neutral microdroplets,² droplets can also carry substantial net charge depending on their method of preparation. A droplet can sustain only a finite amount of net charge before spontaneous division, determined by the balance between surface tension, which stabilizes the spherical shape of the droplet, and electrostatic repulsion among like charges, which promotes droplet fragmentation. The maximum charge a droplet can hold is reached at the Rayleigh limit, where these forces exactly counterbalance one another. The proximity of a charged droplet to the Rayleigh limit is commonly expressed by the fissility parameter X given by

where Q is the total droplet charge, γ is the surface tension, ε₀ is the vacuum permittivity, and R₀ is the droplet radius. The Rayleigh limit is at X = 1. In this work, we examine droplets both near and well below this limit in order to (i) characterize ion spatial distributions as a function of charge density, a study that has not been performed thus far, and (ii) determine whether like-charge ion encounters remain frequent even at lower charge states.

2 Systems and simulation methods

Molecular dynamics (MD) simulations of aqueous charged droplets were conducted by utilizing the Nanoscale molecular dynamics (NAMD) ver. 2.14 software⁴² and visualized by visual molecular dynamics⁴³ (VMD) ver. 1.9.4. The composition of the systems and relevant physical parameters are shown in Table 1.

Table 1 Droplet composition and relevant physical parameters at T = 300 K. The first column presents the approximate number of H₂O molecules (N_H2O) in the droplet, the second column the droplet's equimolar radius (R_e), the third column the radius of the bulk-like region (R_b) computed as in our previous work^44,45 the fourth column the Rayleigh charge (Q_R) of a droplet estimated at X = 1 using eqn (1) in the main text where radius R = R_e and surface tension equals to 69.3 mN m⁻¹ for the TIP4P/2005 model⁴⁶ and 48.8 mN m⁻¹ for the TIP3P-CHARMM model.⁴⁷ The Q_R values with the asterisk refer to the TIP3P-CHARMM model. The fifth to seventh column present the number of ions (N_Na⁺, N_Cl⁻, N_OH⁻, respectively) and the corresponding uniform concentration in a spherical droplet estimated by dividing the number of ions with the droplet volume. The eighth column shows the X value for the maximum number of ions simulated in the droplet sizes reported in the first column by the TIP4P/2005 model

NH2O	Re [nm]	Rb [nm]	QR [e]	NNa^+; CNa^+ [M]	NCl^−; CCl^− [M]	NOH^−; COH^− [M]	X
1.0 × 10^3	1.93	0.95	10.4; 9.1*	2, 8; 0.11, 0.44	2, 8; 0.11, 0.44	2, 8; 0.11, 0.44	0.59
1.0 × 10^4	4.15	3.2	32.9; 27.6*	19, 33; 0.11, 0.18	33; 0.18	18, 33; 0.10, 0.18	1.0
1.5 × 10^4	4.75	3.7	40.2	38; 0.14	38; 0.14	—	0.89

---

In droplets charged with Na⁺ or Cl⁻ ions the water molecules were modeled by the TIP4P/2005 (transferable intermolecular potential with four points/2005) model⁴⁸ combined with the ion parameters from OPLS (Optimized Potentials for Liquid Simulations).^49–51 Ion-oxygen parameters were calculated by the Lorentz–Berthelot rules.

In droplets charged with OH⁻ the H₂O molecules were modeled by the TIP3P-CHARMM model⁵² and the OH⁻ using the parameters reported by Jorgensen et al.⁵³ Force field parameters for OH⁻ are summarized in Table S1 in supplementary information (SI). Despite the simplicity of the model⁵³ used here, it reproduces hypercoordinated solvation structure^{20,22,23,26,54,55} of OH⁻ with 4.2 closest oxygen sites of H₂O at a distance of 2.55 Å (maximum of the O (OH⁻)–O (H₂O) radial distribution function) and two oxygen sites of H₂O molecules further away up to a distance of 3.05 Å as indicated by the radial distribution function of the systems that we calculated.

The Rayleigh charge (Q_R) was estimated by setting X = 1 in eqn (1). The number of ions in the studied droplets varied between Q_R to a smaller amount that corresponds to X ≈ 0.3, and X ≈ 0.5. Pristine aqueous droplets were used as reference systems.

The initial configurations for the charged droplets were prepared using well-equilibrated pristine droplet configurations after 300 ns of equilibration time at T = 300 K, by randomly replacing H₂O molecules with the ions. Extensive simulations with production time near to microsecond for sodiated and chlorinated droplets were also performed at T = 200 K in order to compare the room temperature H₂O orientation with that in the heterogeneous environment of a supercooled droplet.^44,45

The equilibrium MD simulations of droplets were performed in the canonical ensemble by placing the droplet at the center of a spherical cavity where the spherical boundary condition was applied. The center of mass (COM) of the droplet was restrained at the center of a spherical cavity. The restraint of the COM was applied by using the collective variables module, COLVARS.⁵⁶ The cavity radius was 10 Å larger than the droplet's R_e when X is approximately less than 0.4 and three times the droplet radius for droplets with X ≈ 1. This cavity size ensures that the development of droplet's shape fluctuations are not suppressed by the spherical walls. Direct electrostatics was implemented with a cutoff equal to the diameter of the cavity.

Newton's equation of motion for each atomic site was integrated using the velocity Verlet⁵⁷ algorithm with a time step of 2 fs. Temperature was maintained at 300 K using the Langevin thermostat applied to all oxygen sites and ions with a damping coefficient of 1 ps. The duration of the production runs were in the range of 150–200 ns at T = 300 K, while at T = 200 K the production run for N_H2O ≈ 10³ was 0.86 µs, and for N_H2O ≈ 10⁴, 0.4 µs.

To examine the H₂O dynamics by computing the diffusion coefficient, microcanonical ensemble (NVE, where N is the number of molecules, V the volume, and E the total energy) simulations of droplets containing ≈10⁴ H₂O molecules and 33 Na⁺ or 33 Cl⁻ ions were performed in the temperature range of 270–300 K. The vehicular diffusion coefficient of OH⁻ was computed in droplets composed of ≈10⁴ H₂O molecules and 18 OH⁻ ions. In the NVE simulations, the setup was the same as that for the NVT ensemble, but the Langevin damping coefficient was reduced to 0.01 ps⁻¹.

2.1 Umbrella sampling simulations

The potential of mean force (PMF) along the interionic distance (ξ) between two ions of the same sign was computed for 2Na⁺, 2Cl⁻, 2OH⁻ ions in a droplet of N_H2O ≈ 10³ at 300 K using umbrella sampling⁵⁸ (US). For OH⁻ the PMF was also computed in a bulk solution which was simulated with periodic boundary conditions (PBC).

The interionic distance, ξ, in the range of 3–10 Å was separated in windows of width of 0.2 Å. The last configurations of equilibrium simulations of a pair of Na⁺ or Cl⁻ or OH⁻ ions in N_H2O ≈ 10³ constituted the initial configuration for steered MD.⁵⁹ Steered MD generated the “seeds” for each window of US simulations. The length of the sampling in each window was 20 ns. The steering force constant of 5 kcal mol⁻¹ Å⁻² was applied to pull two ions from their initial position to their target location at constant velocity. Sampling windows were spaced uniformly at 0.2 Å with a biasing potential of 50 kcal mol⁻¹. These parameters were selected to ensure sufficient overlap between adjacent windows.

The US simulations for Na⁺ or Cl⁻ were performed in the TIP4P/2005 H₂O model with ion parameters taken from OPLS^49,50 with a timestep of 2 fs. The PMF of a pair of OH⁻ was computed with two force fields: (a) TIP3P-CHARMM with OH⁻ parameters from ref. 53 and (b) scaled charge⁶⁰ OH⁻ in TIP4P/2005 model (parameters are presented in Table S1 in SI). The use of different water and OH⁻ models reflects the availability of validated parameterizations for the respective ions and allows consistency with previous studies of hydroxide solvation. Both models for OH⁻ were applied in a droplet and bulk system. In the bulk simulations a simulation cubic cell with dimension 32 Å was prepared containing N_H2O ≈ 10³ molecules and two pairs of NaOH. Periodic boundary conditions were applied in the three dimensions and the electrostatic forces were treated with particle mesh Ewald. The sampled distributions of different windows were converted into the PMF profile using the weighted histogram analysis method⁶¹ (WHAM). The implementation of the WHAM code developed by Grossfield⁶² iteratively reconstructed the unbiased PMF from the sampling results of individual windows.

3 Results and discussion

3.1 Droplet structure and ion distributions

Molecular dynamics simulations of droplets composed of N_H2O ≈ 10³–1.5 × 10⁴ charged with OH⁻, Cl⁻ or Na⁺ ions were performed. The simulated systems spanned fissility parameters from X ≈ 0.3 to X ≈ 1. In this section we present the radial probability densities of OH⁻, Cl⁻ and Na⁺ ions.

Fig. 1 shows a typical snapshot of a droplet where the bulk-like interior (0 < r < R_b where r is the distance from the droplet's COM), sub-surface (R_b < r < R_v), and surface (R_v < r < R_e) regions have been marked. The method for determining the R_b, R_v which is the outer radius of the subsurface determined by using the Voronoi volume of the molecules, and R_e, which is the equimolar radius were described in our previous work.⁴⁴ The interior region, 0 < r < R_b, is characterized by bulk density and structure,⁶³ the sub-surface region by high ion density,^64,65 lower mass density, lower number of H-bonds per molecule,⁵¹ and enhanced diffusion that we will discuss in detail later. The surface region is characterized by high level of inhomogeneity due to shape fluctuations, more disrupted H-bonded network⁵¹ than the subsurface, more enhanced diffusion than the subsurface, and strong orientation of the H₂O molecules due to the presence of the vapor–liquid interface and the proximity of the ions to the interface.^64,66

Fig. 1 Typical snapshot of a droplet composed of N_H2O ≈ 10³ and N_Na⁺ = 8 ions, where for clarity only the oxygen sites of the H₂O molecules are shown. In the snapshot, differently colored spherical regions⁴⁴ are shown at certain distance from the droplet's center of mass (COM). The interior orange-colored sphere from the droplet's COM to a distance r = R_b displays the bulk-like region. The green region at R_b < r < R_v is the subsurface region, and the grey region at R_v < r < R_e is the surface region.

Fig. 2(a) shows the H₂O radial density profiles of droplets containing Na⁺ or Cl⁻ or OH⁻ ions as a function of the distance from the droplet's center of mass (COM) at 300 K. Fig. 2(b) and (c) show the Na⁺ and Cl⁻ ion radial distributions (the profiles have been divided by the volume element 4πr²dr), respectively, using the TIP4P/2005 H₂O model with OPLS ion parameters. Fig. 2(d) shows the ion radial distribution for OH⁻ ions using TIP3P-CHARMM H₂O model⁵² and OH⁻ parameters reported by Jorgensen et al.⁵³ Here we compare the OH⁻ distributions with those of Na⁺ and Cl⁻ using different force fields. We can make this comparison because the droplet is a conductor, therefore it is expected that the ions will transfer to the surface. In our previous work we showed that for multiple ions the usage of a polarizable force field for Na⁺ and Cl⁻ led to small changes in the ion distribution⁶⁶ relative to a non-polarizable force field.

Fig. 2 (a) Water density (ρ) as a function of distance (r) from the droplet's center of mass (COM) at 300 K when the droplet has charge near the Rayleigh limit. Dashed, dash-dotted, and solid lines represent droplet composed of N_H2O ≈ 10³, ≈10⁴, and ≈1.5 × 10⁴, respectively. (b) Na⁺ ion radial distribution profiles for N_Na⁺ = 2,8 in N_H2O ≈ 10³, N_Na⁺ = 19, 33 in N_H2O ≈ 10⁴, and N_Na⁺ = 38 in N_H2O ≈ 1.5 × 10⁴. (c) Same as (b) but for Cl⁻. (d) Same as (b) but for 2, 8, 18, 33 OH⁻ ions. The shaded area indicates the error bars. Error bars were computed following the standard error of the mean, with M blocks. Each block was ∼20 ns long. Configurations were sampled every 5 ps.

When only two ions are present in a droplet composed of N_H2O ≈ 10³, regardless of the ion type, the ions are uniformly distributed throughout the droplet's volume. For Na⁺ and Cl⁻ ions at X ≈ 1 and ≈0.33 the maximum of the ion radial distribution lies at the subsurface region and the distribution decays exponentially toward the droplet's center as we have discussed in previous works.^64,65,67 A noticeable difference is that in N_OH⁻ = 18 the OH⁻ ions are deeper below the droplet boundary defined by R_e than the Na⁺ and Cl⁻ ions and they are more uniformly distributed within the droplet at r < R_b than the Na⁺ and Cl⁻ ions. In order to examine the OH⁻ location further, the ion distribution was also computed for N_OH⁻ = 33 in N_H2O ≈ 10⁴, where the system is above the Rayleigh limit (X > 1). The charged droplets at X > 1 readily fragment by emitting clusters of solvated OH⁻ ions. We prevented fragmentation by enclosing the droplet in a spherical cavity with a radius 5 Å larger than R_e. We find that similarly to the N_OH⁻ = 18 case, the ions are more uniformly distributed throughout the droplet relative to Na⁺ and Cl⁻ ions, without a significant maximum in the subsurface region. The droplets with N_OH⁻ = 33 appear to have more ions near R_e because of the larger shape fluctuations due to the higher charge. In reality the ions have the same depth from the surface.

The observation of the deeper location of the OH⁻ from the interface is supported by the results of Jungwirth et al.⁶⁸ who demonstrated that the OH⁻ ions experience a weak interfacial repulsion by using data from both synchrotron photoelectron spectroscopy experiments and MD simulations with a polarizable force field. Even though the present simulations used a non-polarizable hydroxide model in TIP3P-CHARMM, the repulsion of the OH⁻ from the surface is captured in the ion distribution when compared with the Cl⁻ and Na⁺ ion distribution. While on the average the OH⁻ ions are found two to three H₂O layers below the surface, we also identified fluctuations where the OH⁻ are part of the outer surface layer. As we showed in previous research⁶⁶ for Na⁺ and Cl⁻ ions using polarizable model and OPLS parameters in droplets, when the ions are near the surface they still maintain on the average their first solvation shell with only a minute reduction relative to the fully solvated ion. Reduction in the solvation appears in the second solvation shell. Differently, when a OH⁻ ion is closest to the surface due to position and shape fluctuations, the oxygen site lies toward the droplet interior where it is well solvated while the H site points toward the vacuo without other H₂O molecules H-bonded with it.

Additional analysis of the droplet's shape fluctuations that provide insight into the ion location and of the orientation of the H₂O molecules^51,64,65 are presented in Section S2 and S3 in SI. Although these fluctuations do not directly determine contact ion-pairing that we discuss in the next section, they may contribute to the heterogeneity of the interfacial region where ion correlations are most pronounced.

3.2 Like-charge contact ion-pairing

Here we examine the likelihood of contact ion-pairing^30,69–71 between two ions of the same sign in droplets. First we examine the potential of mean force (PMF) of a like-charge single pair and then the radial distribution function (RDF) of multiple ions of the same sign. The focus of the discussion is on the OH⁻ contact ion-pairing, while Na⁺ and Cl⁻ pairing are used for comparison with OH⁻.

3.2.1. Single ion pair in a droplet. The PMFs for two OH⁻ in droplet and bulk solution using the TIP3P-CHARMM force field are shown in Fig. 3. The PMFs show a shallow minimum at ξ_{OH⁻–OH⁻} = 4.0 Å that corresponds to the CIP structure [OH⁻(H₂O)_nOH⁻] where n = 1,2 and rarely 3. The nH₂O molecules form H-bonded bridges with the two oxygens of the OH⁻ pair. A H₂O bridge between two molecules is defined if this H₂O molecule has simultaneously H-bonds with the two OH⁻ ion. In the H-bond the O–O distance is taken to be <3.5 Å and the OHO angle greater than 150°. A typical snapshot of the bridged structure is shown in Fig. 3. In other works, density functional theory geometric optimization identifies the bridging to have the lowest energy structure despite the total charge −2 of the system.⁵⁵ This metastable [OH⁻(H₂O)₂OH⁻] structure has predominantly been found in low hydration environment.⁵⁵ In models that allow for electronic polarization the CIP H₂O bridged structures have the potential to reduce direct electrostatic repulsion by distributing the charge through the hydrogen-bond network. We explore this direction with ab initio molecular dynamics (AIMD) calculations in future work.

Fig. 3 Potential of mean force (PMF) as a function of the oxygen–oxygen interionic distance, ξ_{OH⁻–OH⁻}, computed by umbrella sampling in a droplet composed of N_H2O ≈ 10³ and a single pair of OH⁻ ions and in bulk solution at 300 K. The TIP3P-CHARMM water model with OH⁻ parameters from ref. 50 is used. The shaded areas indicate error bars, calculated using the standard error of the mean. In the typical snapshots the oxygen site of the OH⁻ is colored blue. The snapshot at ξ_{OH⁻–OH⁻} ≈ 4.0 Å shows a H₂O molecule bridging OH⁻–OH⁻ CIP, and at ξ_{OH⁻–OH⁻} ≈ 6.4 Å shows a OH⁻–OH⁻ configuration with a few water molecules along the O–O interionic distance.

In Fig. 3 a second broad shallow minimum at ξ_{OH⁻–OH⁻} ≈ 6.0–6.4 Å is also found. The two minima are separated by a free energy barrier of ≈2.0 kcal mol⁻¹. A typical snapshot of the pair structure that corresponds to the second minimum is shown in Fig. 3. The second minimum is followed by an insignificant barrier that leads to the free ions that undergo relative diffusive motion.⁷² Since all the free barriers along ξ_{OH⁻–OH⁻} are low with the highest one being the first barrier, we define a two-state process, where the two states are the CIP (right hand-side) and the “free” (left hand-side) OH⁻ configurations in the following reaction

where n = 1–3, k_f denotes the forward and k_b the back rate constant. Because of the presence of negligible free energy barriers beyond the first barrier, we consider that the rate of formation of the CIP depends on the relative diffusion of the OH⁻ ions. The encounter rate of the OH⁻ to form the CIP is discussed in the next section.

Fig. 3 shows that the bulk and droplet PMFs along ξ_{OH⁻–OH⁻} coincide. The similarity between the droplet and bulk PMFs is explained by the OH⁻ ion location in the droplet. As shown in Fig. 2(d) in droplets composed of N_H2O ≈ 10³ and a single pair of OH⁻, the two ions are uniformly distributed in the subsurface and bulk-like interior. Probability density profiles of the distance of each OH⁻ from the droplet's COM when they are in the CIP form (Fig. S5 in SI) show that often the CIP forms near R_b, therefore the environment of the OH⁻ on the average is similar to that of the bulk solution.

It has been established in many studies that the form of a PMF strongly depends on the force field.⁴⁵ In order to test the reliability of the higher propensity of CIP formation between 2OH⁻, the PMF was also computed with the scaled charge force field (see methodology section). The PMFs with this force field are shown in Fig. S6 in SI. The scaled charge force field yields a more repulsive PMF than the TIP3P-CHARMM force field. While absolute PMF depths vary, the presence of a metastable CIP minimum persists across models. Moreover, both models show that the CIP is stabilized via H-bonds formed with H₂O molecules bridges in similar structures.

Now we compare the PMF of OH⁻ with that of Cl⁻ and Na⁺. Fig. 4 shows the potential of mean force (PMF) as a function of the Na⁺–Na⁺ and Cl⁻–Cl⁻ interionic distance, ξ, within a droplet composed of N_H2O ≈ 10³ and a single pair of Na⁺ or Cl⁻ ions.

Fig. 4 Potential of mean force (PMF) as a function of the interionic distance, ξ_Na–Na, computed by umbrella sampling in a droplet composed of N_H2O ≈ 10³ and a single pair of Na⁺ ions at 300 K. The snapshots show the H₂O molecules surrounding the Na⁺ ions (red-coloured spheres) within 3.5 Å from the COM of the pair. The inset shows the PMF as a function of the Cl⁻–Cl⁻ interionic distance in the same system as in the main graph but with a single pair of Cl⁻ ions, computed in the same way as for Na⁺–Na⁺. The axes have the same meaning and units as in the main graph.

It is found that the PMFs for these systems are overall repulsive, with a shallow minimum when the ions are found at ξ_Cl–Cl ≈ 5.1 Å for the Cl⁻ pair and ξ_Na–Na ≈ 4.0–4.5 Å for the Na⁺ pair. The PMF of Cl⁻ shows a more distinct short distance minimum than the Na⁺ ions. This shallow minimum is separated by a low barrier from a broad low-lying free energy basin. During the simulations, the two Cl⁻ are in the majority of the time in the droplet bulk-like interior, therefore, the environment that they experience is on the average bulk-like. These PMFs are in agreement with the PMFs of the same ions computed by other authors in bulk solution.^69–71

We cannot directly compare the PMFs for Cl⁻ and OH⁻ ions because of the different force fields used. However, these different force fields consistently show a higher propensity for Cl⁻ and OH⁻ contact ion-pairing than for the Na⁺ ions. The stronger contact ion-pairing propensity of anions reflects their lower charge density and the directional hydrogen-bonding network that facilitates solvent-bridged configurations that stabilize the CIP, particularly for OH⁻.

It is important to note that hydroxide ions exhibit complex dynamics involving proton transfer¹⁷ and charge delocalization, which are not captured by molecular mechanics force fields as those we use here. While AIMD can in principle describe these effects, their predictions may depend sensitively on the treatment of electronic structure and the relative timescales of ion approach, CIP lifetime, and solvent reorganization. In this context, the present simulations provide a complementary perspective by sampling long-time statistical behavior and ion encounter frequencies. Therefore, the ion-pair configurations identified here should be interpreted as transient, hydrogen-bond-mediated correlations whose thermodynamic stability may differ under more sophisticated models, while their role in enhancing ion co-localization is expected to remain qualitatively valid.

Based on experimental data of polarization-resolved femtosecond vibrational spectroscopy and terahertz time-domain dielectric relaxation measurements, the vibrational and orientational dynamics of water molecules outside the first hydration shell of OH⁻ have been found to take place at 2.5 ± 0.2 ps as in pure H₂O in a bulk solution.⁷³ This suggests that a single OH⁻ does not have a long-range effect on the hydrogen-bonded network of liquid water.⁷³ However, this may not be the case around the OH⁻–OH⁻ pair because of the higher amount of charge accumulated in a short region. This is a topic we explore in future work using AIMD.

3.2.2. Pairing in droplets with multiple ions. Fig. 5(a) shows the radial distribution function (RDF) between Na⁺ ions, and the inset between Cl⁻ ions in droplets with multiple ion. The Na⁺ ions show negligible preference for CIP formation⁵¹ in contrast to Cl⁻ and OH⁻ (Fig. 5(b)) where the CIP pairing is pronounced. The inset of Fig. 5(b) shows a typical snapshot of OH⁻–OH⁻ CIP surrounded by strongly oriented H₂O molecules. Similarly to the study of the single pair, due to the difference in the force fields we cannot directly compare the propensity of ion pairing between the Cl⁻ and OH⁻ ions. However, the higher propensity of contact ion-pairing for two different types of anions, Cl⁻ and OH⁻, using different force fields relative to Na⁺ ions indicates that their preference for pairing arises from fundamental effects that do not have a strong force-field dependence.

Fig. 5 RDF of (a) Na⁺–Na⁺ with the inset of Cl⁻–Cl⁻, in droplets composed of N_H2O ≈ 10³-8 ions and N_H2O ≈ 10⁴-33 ions. (b) O–O in OH⁻–OH⁻ in droplets of the same size as for (a) with N_OH = 8 and 18 ions, respectively. (c) Probability (P(r)) of the distance of the COM of the pairs with interionic distance ξ_{OH⁻–OH⁻} < 4.9 Å from the droplet's COM. The histogram is normalized to 1. Bin width is 0.4 Å.

Fig. 5(c) shows a histogram of the raw data (not divided by 4πr²dr) of the distance of the COM of the OH⁻ pairs with interionic distance <4.9 Å from the droplet's COM. The histogram shows higher probability of the contact ion-pairing in the subsurface region and near R_b, which is consistent with the probability distribution of OH⁻ in Fig. 2(c).

In the analysis that follows we focus on the OH⁻. In the formation of the CIP the hypercoordinated OH⁻ plays a key role.^{20,22,23,26,54,55} Within the approximations of the molecular mechanics force field that we use here, each OH⁻ is surrounded by four strongly oriented H₂O molecules where each acts as a donor of a single hydrogen, while the other hydrogen site points outward towards other H₂O molecules.

The distinct environment of the subsurface region leads to a variation of the dielectric constant relative to the bulk interior. It is expected that the static dielectric will be lower⁷⁴ than that of the bulk solution due to the reduction in the number of H-bonds.⁵¹ Even though the lower dielectric constant may imply stronger repulsion between the ions, the reduced number of H-bonds may facilitate the stabilization of the CIP.

When the two OH⁻ ions approach interionic distance ξ_{OH⁻–OH⁻} = 4.0 Å metastable structures [OH⁻(H₂O)_nOH⁻] where n = 1, 2 and rarely 3 are formed. Rarely, longer bridges of 2 H-bonded H₂O molecules may be also formed. We found that the H site of the OH⁻ does not form H-bonds with the H₂O molecules and that when OH⁻ lies near the surface, the H site points toward the vacuo. We monitored the bridged configurations to compute the equilibrium constant of the process shown in eqn (2), the lifetime of the CIP, and the kinetics of the process.

3.2.3. Equilibrium constant. For every OH⁻ pair we recorded the number of H₂O bridges, zero to three over all the configurations in the MD trajectory. By assuming a Markov process, we found the eigenvalues of the 4 × 4 transition probability matrix of H₂O bridge formation. The unit eigenvalue corresponds to the equilibrium distribution of the number of bridges per pair. We found that this distribution follows closely a Poisson distribution. To validate the Markov process assumption we also examined the distribution of the bridges per pair directly from the MD trajectory that also showed a Poisson distribution. The expectation value of the first moment of the Poisson distribution yields λ = 0.0189 ± 0.0002 per pair and per configuration for a droplet composed of N_H2O ≈ 10³–N_OH⁻ = 8, and λ = (12.2 ± 0.33) × 10⁻⁴ for N_H2O ≈ 10⁴–N_OH⁻ = 18. These values multiplied by the number of pairs yield probability 0.5 and 0.2, respectively, to encounter a pair among the saved configurations of the MD trajectory.

The equilibrium constant of the reaction in eqn (2) was estimated by

where N_A is Avogadro's number and V the droplet's volume. The estimates of K_eq are found to be (0.35 ± 0.1) L mol⁻¹ for N_H2O ≈ 10³–N_OH⁻ = 8 and (0.22 ± 0.1) L mol⁻¹ for N_H2O ≈ 10⁴–N_OH⁻ = 18. The difference in the values is attributed to the inhomogeneity of the droplet environment.

3.2.4. Lifetime of the CIP. We computed the lifetime of the H₂O bridged structures for a droplet composed of N_H2O ≈ 10³–N_OH⁻ = 8 and N_H2O ≈ 10⁴–N_OH⁻ = 18 in the following way. The eigenvalues, λ_k (k = 0–3) of the transition probability matrix for N_H2O ≈ 10³–N_OH⁻ = 8 are 1, 0.78409, 0.0499, and 0.00682 and for for N_H2O ≈ 10⁴–N_OH⁻ = 18 are 1, 0.754, 0.076, and −0.0086. The order of magnitude difference between the second and third eigenvalue shows that there is a time-scale separation between the chemical process of interest and the other molecular motions in the system. The lifetime, τ, of the H₂O bridged structures was estimated from τ = −5 ps/ln(λ₁), where the numerator equals 5 ps because the configurations in the MD trajectory were saved every 5 ps. For N_H2O ≈ 10³–N_OH⁻ = 8 was found τ = 20.5 ps and for N_H2O ≈ 10⁴–N_OH⁻ = 18 τ = 17.7 ps. The agreement between the two different system sizes is good. The 20 ps lifetime significantly exceeds typical hydrogen-bond lifetimes.^75–77 Indicating that these structures persist over multiple solvent reorganization events.

3.3 Diffusion dynamics

Because the formation of CIPs is governed by ion–ion encounter rates, which depend on molecular diffusion, we next quantify the position-dependent diffusion within the droplet. Since the contact ion-pairing depends on encounter frequency, first we examine the radially dependent self-diffusion coefficient of H₂O in a N_H2O ≈ 10⁴ and N_Na⁺ = 33 (or N_Cl⁻ = 33) droplet. We implemented a technique that we presented in a previous article⁷⁸ to estimate the position dependent self-diffusion coefficient. The steps we performed are described below. Firstly, at a specific temperature we collected an ensemble of molecule displacements squared over a finite interval of time (5 ps, 10 ps, 20 ps, and 40 ps) at which the position of the molecule does not change significantly. Examples of distributions of displacements over 40 ps at (a) T = 270 K, and (b) T = 300 K are shown in Fig. S7 in SI. A time interval of 40 ps was chosen to be sufficiently brief to maintain the molecule within its neighborhood, for example, within this time interval a molecule in the interior does not reach the surface boundary. In the temperature range of 230 K–300 K the convergence of the values of the self-diffusion coefficients in the droplet interior for 40 ps is shown in Fig. S8 and S9 in SI.

Secondly, we assumed a functional model of the self-diffusion coefficient dependence on the distance from the droplet's COM, r expresses as

where D₀ can be associated with the self-diffusion coefficient in the bulk-like droplet's interior and R_L is the distance from the COM of the droplet where the apparent self-diffusion coefficient is double that of the one in the center. This quadratic function is chosen because it is the lowest power polynomial that smoothly passes from the droplet's center. Thirdly, the ensemble of displacements was analyzed using the maximum likelihood estimate and the optimal parameters for the r dependence of the model (eqn (4)) of the self-diffusion coefficient were selected. The calculations were performed using the maximum likelihood estimate functions of the R statistical analysis package.⁷⁹ The optimal parameters are presented in Table 2.

Table 2 Parameters used in eqn (4) (second and third column) as a function of temperature (T) (first column) for TIP4P/2005 model at T = 300–230 K and TIP3P-CHARMM model for 300 K. The data for the TIP3P-CHARMM carry an (*)

T [K]	D0 [nm^2 ns^−1]	RL [Å]
300	1.97	49.9
300	3.18 (*)	55.7 (*)
270	0.83	45.5
250	0.37	40.9
230	0.084	29.1

---

The model (eqn (4)) was verified by comparing the computed probability distribution of the displacements squared with the theoretical prediction. For a diffusive process the probability, denoted as P_r (theoretical prediction), of Δ² for a fixed value of position r, is given by the gamma distribution (Γ) with the shape parameter α = 3/2:

where τ is the time interval of 40 ps, Δ is the magnitude of the displacement, and D (r) is the position dependent self-diffusion coefficient. The comparison is shown in Fig. S9 in SI.

Fig. 6 shows the radially dependent self-diffusion coefficients at various T computed using eqn (4) with the parameters presented in Table 2. First we discuss the values of D₀ in Table 2 that provide the self-diffusion coefficient at the bulk-like interior of the droplet. In the N_H2O ≈ 10⁴ (R_e = 4.15 nm) droplet we examine here, the interior pressure was estimated to range from ≈332 bar at 300 K to ≈371 bar at 230 K using the Young-Laplace equation with value of the surface tension for the planar interface. Ref. 80 that shows the convergence of the surface tension in droplets with size led us to consider that an aqueous droplet with R_e = 4.15 nm has reached the bulk value of surface tension. The data in ref. 39 for TIP4P/2005 bulk at 298.15 K and pressure 500 bar give the value of 2.30 [nm² ns⁻¹], which is near the interior values found in this study.

Fig. 6 Radially dependent self-diffusion coefficient of H₂O in a droplet composed N_H2O ≈ 10⁴–N_Na⁺ = 33 ions at different temperatures.

The temperature and pressure dependence of the self-diffusion coefficient of bulk H₂O has been discussed in many articles.³³ It has been found that at pressure 1 bar and above a cross-over temperature the self-diffusion of H₂O follows an Arrhenius equation while below shows significant deviation from the Arrhenius equation that increases with the degree of supercooling. The cross-over temperature has been placed at different values in the literature, and it may be as low as 225 K. We use the Arrhenius equation to fit the temperature dependence of the self-diffusion coefficient values reported in the second column of Table 2, that provides an estimate of the activation barrier of 6.14 ± 0.67 kcal mol⁻¹. The activation energy is ≈1.0 kcal mol⁻¹ higher than the values reported³³ for various H₂O models at 1 bar. The higher value may arise from the higher pressure at the droplet interior that prevents the escape of the H₂O molecule from a cage. Another possible explanation for the activation energy difference is that the prefactor in the Arrhenius equation has a temperature dependence.

The self-diffusion coefficient increases gradually from the droplet's bulk-like interior to the surface. In the temperature range T = 300–230 K near R_e the self-diffusion coefficient increases up to twice that in the bulk-like interior. We expect a similar increase to be present in microdroplets. The moderately increased surface self-diffusion is expected to contribute to the acceleration of reaction rates when the reaction is diffusion-limited such as that of the formation of the CIPs.

3.3.1. Diffusion of OH⁻. In a droplet composed of N_H2O ≈ 10⁴–N_OH⁻ = 18 the OH⁻ ions are found on an average distance of R_avg = 2.98 nm from the droplet's COM. This distance is at the interface between the bulk-like region and subsurface (Fig. 2). We estimate the OH⁻ diffusion coefficient (D_OH⁻) for their motion on a spherical surface using D_OH⁻ = R_avg²〈θ²〉/4t, where 〈θ²〉 = 0.483 and t = 0.5 ns. D_OH⁻ is estimated to be 2.14 nm² ns⁻¹. Using eqn (4) and the parameters in Table 2, the self-diffusion coefficient of TIP3P H₂O at R_avg = 2.98 nm was estimated to be 4.08 ± 0.2 nm² ns⁻¹. These estimates show that the vehicular diffusion of OH⁻ is almost half that of TIP3P-CHARMM H₂O. The slower diffusion is expected because of the tightly bound solvation layer surrounding a OH⁻ ion. The experimental value^15,16 of D_OH⁻ = 5.2–5.4 nm² ns⁻¹ is ≈2.5 times higher than the value estimated here because of the structural diffusion¹⁷ and also a possible different value of the vehicular diffusion. Regardless of the model-dependent values of the diffusion coefficient, the present study suggests that in the depth from the surface where the OH⁻ are located, the self-diffusion of H₂O is moderately higher than its bulk H₂O value and that the vehicular diffusion rate is half that of the self-diffusion of H₂O. In addition, the faster H₂O self-diffusion rate is expected to accelerate the diffusion of OH⁻.

Using D_OH⁻ we estimate k_f for the formation of the CIP using the Smoluchowski equation^81,82

k_f[OH⁻] = 4π(2D_OH⁻)(2R)(1/V) (6)

where 2R is the collision radius taken to be 0.4 nm, which is the first minimum in the PMF shown in Fig. 3, and V is the volume of the droplet taken to be 300 nm³ for a droplet composed of N_H2O ≈ 10⁴. Using the detailed balance condition, we find that the life-time (τ) of a OH⁻ CIP is given bywhere λ = (12.2 ± 0.33) × 10⁻⁴ as defined earlier. It was estimated that τ = 17 ps. This value is in good agreement with the lifetime we estimated in the previous section using the eigenvalues of the transition probability matrix of H₂O bridges. The agreement between independent kinetic estimates (see discussion in previous section) provides support for the robustness of the contact ion-pairing dynamics. Moreover, the agreement of the values suggests that OH⁻ contact ion-pairing occurs in the entire droplet volume. The use of the entire droplet volume is justified by the OH⁻ concentration profiles (Fig. 2) that decay slowly from the maximum of the distribution in the subsurface toward the bulk-like interior. The OH⁻ diffusion is strongly determined by structural diffusion, which is in turn is affected by the different self-diffusion coefficient of H₂O and H-bonded network in the droplet interior vs. in the surface layers. These effects are to be examined with AIMD.

3.4 Implications for microdroplet chemistry

For a droplet of R_e = 1 µm (N_H2O ≈ 1.4 × 10¹¹–N_ion = 1.23 × 10⁵) that has a more relevant size to that used in the experiments, the concentration of the ions if they were uniformly distributed in the entire volume would be 4.87 × 10⁻⁵ M. However, in the charged droplet the ions are mainly distributed in the surface and subsurface region^65,83 (thickness of approximately 1–1.5 nm) that increases the concentration by 333 times (0.016 M) relative to the uniform distribution. Therefore, if the ions participate in the reaction their concentration itself is a factor that is expected to contribute by at least two-three orders of magnitude in the reaction acceleration. This contribution in the rate of reaction arises directly from the mass action law.

Regarding the rate of diffusion of the H₂O molecules, similar to nanoscopic droplets, microdroplets are expected to exhibit enhanced diffusion at the subsurface relative to the interior. In a droplet with R_e = 1 µm the average interionic separation of the ions is ≈10 nm if they were all distributed uniformly on the surface. In practice this distance for OH⁻ ions may be shorter because these ions are located below the surface, which reduces their effective separation. If we consider the experimental value of OH⁻ diffusion,^15,16 that includes both vehicular and structural diffusion, 5.2–5.4 nm² ns⁻¹, the faster diffusion of H₂O in the subsurface relative to the interior and the fact that the ions are confined to diffuse in a narrow shell at the surface the encounter time is estimated to be <1.4 ns. These estimates indicate that ion–ion encounters occur within a nanosecond even at relatively modest subsurface concentrations in micron-sized droplets. Consequently, transient OH⁻–OH⁻ CIPs are expected to form frequently during the droplet lifetime. Although each CIP is short-lived, their continual formation provides repeated opportunities for pre-reactive configurations, effectively increasing the probability of interfacial chemical events. The faster diffusion of H₂O in the droplet surface combined with reduced H-bonding⁵¹ is expected to promote like-charge ion pairing because it can reduce the structural diffusion of OH⁻. We propose that this effect is analogous to the reduced diffusion of hydronium ions on the surface of ice due to local trapping by under-coordinated water molecules.⁸⁴ Moreover, the faster diffusion allows for faster solvent re-organization in capturing an electron potentially detached from a OH⁻ preventing rapid recombination.¹⁴

Fluctuations in the charge distribution, including both the positions of the free ions and the orientation of H₂O molecules, can produce local enhancements of the electric field. In previous work,⁶⁵ we estimated that for a droplet composed of ≈2 × 10⁴ H₂O molecules and 36 Na⁺ ions, the average value of the normal component of the electric field is 0.187 V Å⁻¹, with fluctuations reaching magnitudes approximately twice the average value. These large fluctuations were observed over roughly 2% of the droplet surface area. In nanodroplets, such transient electric fields can approach the magnitude required to induce H₂O dissociation. For droplets with charge near the Rayleigh limit, the average squared electric field scales inversely with the droplet radius. Consequently, as droplet size increases, the average electric field decreases substantially. The charge relaxation to shape fluctuations has been thoroughly described in ref. 85.

A distinct droplet region with high electric field is the Taylor cone appearing on the surface of a droplet near the Rayleigh limit^{66,83,86–88} regardless of the droplet size. These conical shapes have high local values of the electric field at their tips, which may potentially serve as catalytic centers for accelerating chemical reactions.

Although strong electric fields are frequently invoked to explain the enhanced reactivity observed in charged microdroplets,¹¹ the scaling of the average electric field with droplet radius suggests that this mechanism becomes significantly weaker in larger droplets. Consequently, electric-field-driven activation alone may not fully account for the substantial reaction-rate enhancements reported experimentally. This observation motivates the exploration of alternative molecular-level mechanisms that may operate independently of the macroscopic electric field strength. One such mechanism is the formation of correlated ion configurations such as the OH⁻ (or Cl⁻) like-charge contact ion-pairing observed in the present simulations that are present even when the average electric field is relatively small. Regardless of droplet size, the overlap between enhanced diffusion and elevated ion density in the subsurface region maximizes ion–ion encounter rates precisely where non-specific features of the surface region such as lower dielectric constant and reduced H-bonding favor contact ion-pairing.

4 Conclusion

In this work, we examined pre-reactive co-localization via like-charge ion contact ion-pairing. Our simulations revealed that OH⁻ and Cl⁻ ions exhibit a significantly stronger tendency to form CIPs than Na⁺ ions. Hydroxide CIPs, and similarly Cl⁻ pairs, are stabilized by hydrogen-bond-mediated bridging structures involving surrounding water molecules, forming transient configurations [OH⁻(H₂O)_nOH⁻] where n = 1,2 and rarely 3. These structures have the potential to partially delocalize charge through the hydrogen-bond network mitigating electrostatic repulsion and enabling the formation of metastable CIPs.

The molecular mechanics MD simulations performed here estimated that the lifetime of these hydroxide CIPs is ≈18 ps. This time is sufficiently long to allow repeated solvent reorganization and local proton-transfer events within their hydrogen-bonded environment. Although short-lived on macroscopic timescales, their continual formation and decay provide recurring opportunities to sample pre-reactive configurations. In this sense, we speculate that such transient ion pairs can act as precursors to chemically relevant processes, including electron detachment into the surrounding solvent or vacuum interface, field-assisted water dissociation, and radical formation.

An open question in the literature is the H₂O₂ formation^3–13 in aqueous microdroplets. We suggest that the formation of hydroxyl radicals without the assistance of any catalytic species, and their encounter to form H₂O₂ in a microdroplet is very rare. We propose that for H₂O₂ to be produced in an experimentally detectable quantity every binary or unitary step of the reaction mechanism should have the precursors to radicals in close proximity. The OH⁻ contact-ion-paring is a possible structural precursor that it is a center of initiating reactivity such as proton transfer from the surrounding H₂O molecules. The manner in which the presence of a negative ion in the vicinity of OH⁻ facilitates the electron detachment is to be examined by AIMD.

We also characterized the position-dependent self-diffusion of H₂O within droplets with equimolar radius of ≈4.15 nm, over the temperature range 230–300 K. This droplet size is large enough to be used as a representative example for droplets even with orders of magnitude larger diameter because it has a clearly defined bulk-like interior, subsurface and surface regions. We found that the diffusion coefficient increases gradually from the bulk-like interior toward the interface, reaching values up to approximately twice as large near the droplet boundary. While this enhancement can increase encounter frequencies between reactants, its magnitude alone is insufficient to explain the orders-of-magnitude rate accelerations reported in microdroplet experiments.

Our findings suggest that enhanced reactivity in charged droplets arises from the interplay of multiple factors rather than a single dominant mechanism. In particular, transient, hydrogen-bond-bridged contact ion-pairing provides a means of locally concentrating ionic species, namely OH⁻, and reorganizing the solvent environment, thereby creating reactive configurations in the form of dynamic ion correlations even in the absence of strong average electric fields. This mechanism complements other proposed effects, such as interfacial electric fields, pH gradients, and evaporation-driven concentration changes.

A limitation of the present study is the use of molecular mechanics, non-reactive force fields, which do not capture proton transfer, charge delocalization beyond fixed charges, or electronic polarization effects. These processes are particularly relevant for hydroxide ions, where structural diffusion and proton transfer play a central role in dynamics. As a result, the lifetimes and stability of the hydroxide CIPs reported here may differ quantitatively from those obtained using AIMD. Nevertheless, the classical simulations provide statistically robust insight into ion spatial distributions and encounter frequencies, which are expected to remain qualitatively valid. The present studies provide reference data to be compared with future AIMD data. The use of a polarizable molecular mechanics force field is expected to strengthen the stability of the OH⁻–OH⁻ CIP relative to the non-polarizable model by strengthening the polarization in the H₂O bridges. Application of AIMD should carefully account for the response time of the electron density to the dynamics of the H₂O molecules and to the approaching rate of the OH⁻ ions. When AIMD is used for this problem, we suggest to confine the distance between the OH⁻–OH⁻ ions in values near the maximum of the CIP and near the second incipient minimum to compute the ionization energies of the OH⁻ as well as the reactions that the pair initiates in the surrounding molecules. It is expected that the ionization energies of each of the OH⁻ ions will be lower than that of the free OH⁻ ions because of the presence of the neighboring negative charge. A neighboring negative charge may also prevent electron recombination and may initiate proton transfer reactions within the surrounding molecules. The neighbouring negative charge may not be limited to an anion, but it may also be a negatively charged surface. In future research we examine how these precursor configurations evolve into chemically reactive events using reactive and AIMD modeling.

Finally, the distinct radial distribution of hydroxide ions to alkali ions is the starting point to address the differences in ionization efficiency between negative and positive modes in electrospray ionization, where the depth of ions below the interface determine their reactivity with analytes, their likelihood of emission,^83,89 and the size of the emitted cluster.^83,89 We also suggest extension of both the experimental and computational studies to supercooled droplets and higher-temperature regimes and to different H-bonded solvents such as fluorinated alcohols to further clarify the role of hydrogen-bond dynamics and diffusion in governing interfacial chemistry.

Author contributions

H. N. formal analysis, visualization, validation, investigation. S. C. conceptualization, formal analysis, methodology, validation, investigation, visualization, writing – original draft, writing – review & editing, funding acquisition, supervision, resources.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data supporting this article have been included as part of the supplementary information (SI). Supplementary information: (S1) Force field parameters for OH⁻. (S2) Characterization of the droplet's shape fluctuations using the moment of inertia. (S3) Orientation of the H₂O molecules in droplets. (S4) Like-charges ion pairing. (S5) Diffusion dynamics. See DOI: https://doi.org/10.1039/d6cp01638c.

Acknowledgements

SC acknowledges an NSERC-Discovery grant (Canada) for funding this research. HN acknowledges the province of Ontario and The University of Western Ontario for an Ontario Graduate Scholarship. Digital Research Alliance of Canada is acknowledged for providing the computing facilities to perform this research.

References

1. Z. Wei, Y. Li, R. G. Cooks and X. Yan, Annu. Rev. Phys. Chem., 2020, 71, 31–51.
2. K. R. Wilson and A. M. Prophet, Annu. Rev. Phys. Chem., 2024, 75, 185–208.
3. J. K. Lee, K. L. Walker, H. S. Han, J. Kang, F. B. Prinz, R. M. Waymouth, H. G. Nam and R. N. Zare, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 19294–19298.
4. J. K. Lee, H. S. Han, S. Chaikasetsin, D. P. Marron, R. M. Waymouth, F. B. Prinz and R. N. Zare, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 30934–30941.
5. C. Zhu and J. S. Francisco, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 19222–19224.
6. M. T. Martins-Costa and M. F. Ruiz-López, J. Am. Chem. Soc., 2023, 145, 1400–1406.
7. S. F. Nami-Ana, M. A. Mehrgardi, M. Mofidfar and R. N. Zare, J. Am. Chem. Soc., 2024, 146, 31945–31949.
8. M. F. Ruiz-Lopez, J. S. Francisco, M. T. Martins-Costa and J. M. Anglada, Nat. Rev. Chem., 2020, 4, 459–475.
9. R. A. LaCour, J. P. Heindel, R. Zhao and T. Head-Gordon, J. Am. Chem. Soc., 2025, 147, 6299–6317.
10. A. J. Colussi, J. Phys. Chem. Lett., 2025, 16, 3366–3370.
11. H. Hao, I. Leven and T. Head-Gordon, Nat. Commun., 2022, 13, 280.
12. A. J. Colussi, J. Am. Chem. Soc., 2023, 145, 16315–16317.
13. M. A. Eatoo and H. Mishra, Chem. Sci., 2024, 15, 3093–3103.
14. C. Petersen, J. Thøgersen, S. K. Jensen and S. R. Keiding, J. Chem. Phys. A, 2007, 111, 11410–11420.
15. J. Sluyters and M. Sluyters-Rehbach, J. Phys. Chem. B, 2010, 114, 15582–15589.
16. B. Halle and G. Karlström, J. Chem. Soc., Faraday Trans. 2, 1983, 79, 1047–1073.
17. M. Chen, L. Zheng, B. Santra, H.-Y. Ko, R. A. DiStasio Jr, M. L. Klein, R. Car and X. Wu, Nat. Chem., 2018, 10, 413–419.
18. C. De Grotthuss, Ann. Chim., 1806, 58, 54.
19. M. Tuckerman, K. Laasonen, M. Sprik and M. Parrinello, J. Phys. Chem., 1995, 99, 5749–5752.
20. P. B. Balbuena, K. P. Johnston and P. J. Rossky, J. Phys. Chem., 1996, 100, 2706–2715.
21. A. Botti, F. Bruni, S. Imberti, M. Ricci and A. Soper, J. Chem. Phys., 2003, 119, 5001–5004.
22. M. Zapaowski and W. M. Bartczak, Comput. Chem., 2000, 24, 459–468.
23. E. Brodskaya, A. P. Lyubartsev and A. Laaksonen, J. Phys. Chem. B, 2002, 106, 6479–6487.
24. Y.-L. S. Tse, C. Chen, G. E. Lindberg, R. Kumar and G. A. Voth, J. Am. Chem. Soc., 2015, 137, 12610–12616.
25. R. Buchner, G. Hefter, P. M. May and P. Sipos, J. Phys. Chem. B, 1999, 103, 11186–11190.
26. D. Bucher, A. Gray-Weale and S. Kuyucak, J. Chem. Theory Comput., 2010, 6, 2888–2895.
27. N. Agmon, H. J. Bakker, R. K. Campen, R. H. Henchman, P. Pohl, S. Roke, M. Thämer and A. Hassanali, Chem. Rev., 2016, 116, 7642–7672.
28. W. Cao, H. Wen, S. S. Xantheas and X.-B. Wang, Sci. Adv., 2023, 9, eadf4309.
29. N. F. Van Der Vegt, K. Haldrup, S. Roke, J. Zheng, M. Lund and H. J. Bakker, Chem. Rev., 2016, 116, 7626–7641.
30. D. Laria and R. Fernández-Prini, J. Chem. Phys., 1995, 102, 7664–7673.
31. I. Palaia, A. Goyal, E. Del Gado, L. Samaj and E. Trizac, J. Phys. Chem. B, 2022, 126, 3143–3149.
32. K. Krynicki, C. D. Green and D. W. Sawyer, Faraday Discuss. Chem. Soc., 1978, 66, 199–208.
33. I. N. Tsimpanogiannis, O. A. Moultos, L. F. Franco, M. B. D. M. Spera, M. Erdös and I. G. Economou, Mol. Simul., 2019, 45, 425–453.
34. D. Rozmanov and P. G. Kusalik, J. Chem. Phys., 2012, 136, 044507.
35. T. I. Morozova, N. A. Garca and J.-L. Barrat, J. Chem. Phys., 2022, 156, 126101.
36. S. H. Lee and J. Kim, Mol. Phys., 2019, 117, 1926–1933.
37. P. Montero de Hijes, E. Sanz, L. Joly, C. Valeriani and F. Caupin, J. Chem. Phys., 2018, 149, 094503.
38. S. Tazi, A. Botan, M. Salanne, V. Marry, P. Turq and B. Rotenberg, J. Phys.: Condens. Matter, 2012, 24, 284117.
39. G. Guevara-Carrion, J. Vrabec and H. Hasse, J. Chem. Phys., 2011, 134, 074508.
40. X. Teng, B. Liu and T. Ichiye, J. Chem. Phys., 2020, 153, 104510.
41. P. Mark and L. Nilsson, J. Phys. Chem. A, 2001, 105, 9954–9960.
42. J. C. Phillips, D. J. Hardy, J. D. C. Maia, J. E. Stone, J. V. Ribeiro, R. C. Bernardi, R. Buch, G. Fiorin, J. Henin, W. Jiang, R. McGreevy, M. C. R. Melo, B. K. Radak, R. D. Skeel, A. Singharoy, Y. Wang, B. Roux, A. Aksimentiev, Z. Luthey-Schulten, L. V. Kale, K. Schulten, C. Chipot and E. Tajkhorshid, J. Chem. Phys., 2020, 153, 044130.
43. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graph., 1996, 14, 33–38.
44. S. M. A. Malek, V. Kwan, I. Saika-Voivod and S. Consta, J. Am. Chem. Soc., 2021, 143, 13113–13123.
45. V. Kwan, S. R. Maiti, I. Saika-Voivod and S. Consta, J. Am. Chem. Soc., 2022, 144, 11148–11158.
46. C. Vega and E. de Miguel, J. Chem. Phys., 2007, 126, 154707.
47. A. E. Ismail, G. S. Grest and M. J. Stevens, J. Chem. Phys., 2006, 125, 014702.
48. J. L. F. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 234505.
49. J. Aqvist, J. Phys. Chem., 1990, 94, 8021–8024.
50. J. Chandrasekhar, D. C. Spellmeyer and W. L. Jorgensen, J. Am. Chem. Soc., 1984, 106, 903–910.
51. V. Kwan and S. Consta, J. Am. Soc. Mass Spectrom., 2020, 32, 33–45.
52. B. R. Brooks, C. L. Brooks III, A. D. Mackerell Jr., L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York and M. Karplus, J. Comput. Chem., 2009, 30, 1545–1614.
53. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
54. A. Dutta and T. Lazaridis, J. Phys. Chem. B, 2024, 128, 12161–12170.
55. I. Zadok, H. Long, B. Pivovar, A. Roznowska, A. Michalak, D. R. Dekel and S. Srebnik, J. Mol. Liq., 2020, 313, 113485.
56. G. Fiorin, M. L. Klein and J. Hénin, Mol. Phys., 2013, 111, 3345–3362.
57. M. Allen and D. Tildesley, Computer Simulation of Liquids, Oxford University Press, 2017.
58. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
59. S. Izrailev, S. Stepaniants, B. Isralewitz, D. Kosztin, H. Lu, F. Molnar, W. Wriggers and K. Schulten, Computational Molecular Dynamics: Challenges, Methods, Ideas: Proceedings of the 2nd International Symposium on Algorithms for Macromolecular Modelling, Berlin, May 21-24, 1997, 1999, pp. 39-65.
60. P. Habibi, A. Rahbari, S. Blazquez, C. Vega, P. Dey, T. J. H. Vlugt and O. A. Moultos, J. Phys. Chem. B, 2022, 126, 9376–9387.
61. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
62. A. Grossfield, An implementation of WHAM: the Weighted Histogram Analysis Method Version 2.1.0, https://membrane.urmc.rochester.edu/sites/default/files/wham/doc.pdf, 2013.
63. V. Kwan, S. Consta and S. M. Malek, J. Phys. Chem. B, 2023, 128, 193–207.
64. V. Kwan, A. Malevanets and S. Consta, J. Phys. Chem. A, 2019, 123, 9298–9310.
65. V. Kwan and S. Consta, Chem. Phys. Lett., 2020, 746, 137238.
66. V. Kwan and S. Consta, J. Phys. Chem. A, 2022, 126, 3229–3238.
67. A. Malevanets and S. Consta, J. Chem. Phys., 2013, 138, 184312.
68. B. Winter, M. Faubel, R. Vácha and P. Jungwirth, Chem. Phys. Lett., 2009, 474, 241–247.
69. C. H. Choi, S. Re, M. H. Rashid, H. Li, M. Feig and Y. Sugita, J. Phys. Chem. B, 2013, 117, 9273–9279.
70. L. X. Dang, B. M. Pettitt and P. J. Rossky, J. Chem. Phys., 1992, 96, 4046–4047.
71. E. Guardia, R. Rey and J. Padró, J. Chem. Phys., 1991, 95, 2823–2831.
72. G. Ciccotti, M. Ferrario, J. T. Hynes and R. Kapral, Chem. Phys., 1989, 129, 241–251.
73. J. Hunger, L. Liu, K.-J. Tielrooij, M. Bonn and H. Bakker, J. Chem. Phys., 2011, 135, 124517.
74. C. Zhang, S. Yue, A. Z. Panagiotopoulos, M. L. Klein and X. Wu, Phys. Rev. Lett., 2023, 131, 076801.
75. A. Luzar and D. Chandler, Nature, 1996, 379, 55–57.
76. A. Luzar, J. Chem. Phys., 2000, 113, 10663–10675.
77. J. Liu, X. He, J. Z. Zhang and L.-W. Qi, Chem. Sci., 2018, 9, 2065–2073.
78. S. Consta, L. M. Wingen, Y. Qin, V. Perraud and B. J. Finlayson-Pitts, Phys. Chem. Chem. Phys., 2024, 26, 28220–28233.
79. R Core Team, R: A language and environment for statistical computing, R Foundation for Statistical Computing, Vienna, Austria, 2021.
80. S. Malek, P. H. Poole and I. Saika-Voivod, J. Chem. Phys., 2019, 150, 234507.
81. M. Smoluchowski, Z. Physik. Chem, 1917, 92, 129–155.
82. M. Doi, J. Phys. A-Math. Gen., 1976, 9, 1479–1495.
83. V. Kwan, R. O'Dwyer, D. Laur, J. Tan and S. Consta, J. Phys. Chem. A, 2021, 125, 2954–2966.
84. K. Park, W. Lin and F. Paesani, J. Phys. Chem. B, 2014, 118, 8081–8089.
85. A. I. Zhakin, Phys. Usp., 2013, 56, 141.
86. L. Rayleigh, Philos. Mag., 1882, 14, 184–186.
87. M. Sharawy and S. Consta, J. Phys. Chem. A, 2016, 120, 8871–8880.
88. S. Consta, J. Phys. Chem. B, 2022, 126, 8350–8357.
89. J. V. Iribarne and B. A. Thomson, J. Chem. Phys., 1976, 64, 2287–2294.

---