Geometric and energetic considerations of surface fluctuations during ion transfer across the water-immiscible organic liquid interface

Most of our current understanding of ion transfer across the interface between two immiscible liquids (typically water and a weakly polar solvent) is based on electrochemical measurements (see Refs. 1–3 for recent reviews) and the corresponding theoretical and computational studies of the free energy change associated with this process (see Refs. 4 and 5 for recent reviews). Considerably more challenging is the experimental determination and interpretation of the dynamics of ion transfer. Specifically, while one can describe a voltage-current relationship for the transfer of ions across the interface in terms of first order rate constants, the published values of these rate constants over the last three decades seem to depend on the experimental method used. The various reasons provided for this situation include limitations in the instruments’ time resolution, uncertainty in the value of the electric potential drop across the interface, and unclear separation between diffusion through the bulk phases and the crossing of the interface, among others.^2,4

Putting aside the experimental uncertainty in the values of rate constants, development of a quantitative predictive theory for the rate of ion transfer requires a reasonable understanding of the mechanism by which an ion exchanges its solvation environment as it crosses the interface. Existing theoretical models such as those based on the Nernst-Planck equation (diffusion in an external field) have been successful in simulating experimentally observable cyclic voltammetry curves.⁴ However, ion transfer rate constants derived from this description give values that are up to two orders of magnitude larger than measured.^4,6 The slower rate constant suggests a barrier or some type of retardation effects. Better agreement with experiments can be obtained when the diffusion equation is numerically solved with a molecularly derived potential of mean force (PMF) as an input external force.⁷ If the PMF exhibits a barrier, one could employ transition state theory to estimate the rate constant.⁸ 

While the PMF calculated from a molecularly detailed model can be useful for some clues about the mechanism of ion transfer (and the nature of the barrier), it is an equilibrium mean field description in which the role of dynamic fluctuations is mostly ignored. Molecular dynamics simulations provide ample evidence that these fluctuations are critical for elucidating the mechanism of the transfer process. In particular, simulations suggest the importance of water molecules’ protrusions (water “fingers”) into the organic phase as a mechanism that, while necessary, results in slowing down the transfer rate into the organic phase and enhancing the reverse process by producing water protrusions in response to the approaching ion.^9–14 Theoretical models based on this picture provide additional support for this mechanism. For example, Marcus has estimated the rate of ion attachment and detachment from a water “finger” and diffusion along it and concluded that the rate limiting step for ion transfer is the motion of an ion along a “solvation coordinate.”¹⁵ Urbakh, Kornyshev, and co-workers developed a coupled, two-dimensional Langevin equation model of ion transfer in which one degree of freedom describes the ion’s motion and the second describes dynamic interfacial fluctuations modeled as capillary waves.^16,17 Numerical solutions and the analytical expression derived for the rate constant suggests that ion-surface coupling can result in a slowing down of the ion motion relative to simple diffusion.

In order to obtain molecular-level insight into the ion transfer process, it is useful to describe the collective ion-solvent interactions using a single (or a few) degree(s) of freedom “reaction coordinate(s).” In this paper, we utilize two independent collective coordinates and use them to examine both the thermodynamic and the dynamic fluctuations during the transfer of a hydrophilic Cl⁻ ion across the water/nitrobenzene interface. One of these coordinates, suggested by Schweighofer and Benjamin,¹⁰ is based on the solvation energy of the ion. The other, recently proposed by Morita and co-workers, is geometrically based and attempts to characterize the finger structure formed during the ion transfer.¹⁸ 

The rest of the paper is organized as follows. In Section II, we describe the simulation details, the umbrella sampling free energy calculations with biasing potential, and the two reaction coordinate calculations. In Section III, we discuss the results of these calculations and in Section IV, we summarize this work and present our conclusions.

The water/nitrobenzene liquid/liquid interfacial system consists of two adjacent slabs of 986 water molecules and 252 nitrobenzene molecules in a 31.28 × 31.28 × 300.0 Å rectangular box. The liquid/liquid interface is located in the X-Y plane at Z ≈ 0, with the water phase in the region of Z < 0 and the nitrobenzene phase in the Z > 0 region. Each liquid phase is in equilibrium with its respective vapor phase; only one liquid/liquid interface is present. Periodic boundary conditions are applied in all directions and a soft reflecting potential wall is located 5 Å from the simulation box boundaries in the Z-direction to prevent molecules from crossing into the adjacent vapor phase.

Intermolecular interaction potentials are represented as the pairwise sum of Lennard-Jones and Coulomb terms,

where r is the distance between atom centers i and j. Standard Lorentz-Berthelot combining rules, σ_ij = (σ_i + σ_j)/2 and ε_ij = (ε_iε_j)^1/2, are used to generate Lennard-Jones parameters for all mixed interactions. For water, we use a version of the flexible simple point charge (SPC) model¹⁹ with intramolecular potentials as described by Kuchitsu and Morino.²⁰ The nitrobenzene molecules are modeled using a 14-site all-atom model described previously.²¹ In brief, this nitrobenzene model uses intramolecular parameters taken from the Amber force field²² and modified intermolecular parameters taken from ab initio calculations and reproduces the experimental dipole moment and enthalpy of vaporization. The chloride ion Lennard-Jones parameters (σ = 3.934 Å and ε = 0.832 kcal/mol) were used in previous simulations and give reasonable agreement with the experimental hydration free energy.²³ 

Molecular dynamics simulations are performed with our in-house MD software that uses the velocity Verlet algorithm to integrate the equations of motion. Unless otherwise specified, all reported simulation data are averaged over 10 independently generated starting configurations, each run for 1 ns, for a total ensemble average over 10 ns of simulation time. All simulations were performed at 298 K and with an integration time step of 0.5 fs.

The local free energy of the single chloride ion is given by

where δ is the Dirac delta function and z_I is the ion’s position along the axis normal to the interface. P(Z) is the probability to find the ion at Z and the ensemble average is calculated over all possible solvent configurations while the ion is located at z_I = Z. The reversible work required to move the ion from position Z₁ to Z₂ would be A(Z₁)-A(Z₂). While direct sampling of the ion position can theoretically be used to compute A(Z), in practice if at a given interval A(Z) varies significantly (more than 5k_BT), and/or if significant structural relaxation of the solvents is necessary to sample equilibrium configurations (as is the case here with the formation of “water fingers”), a different approach is required. In the umbrella sampling method,^24–26 the interval of interest is divided into a series of N overlapping lamellae parallel to the liquid-liquid interface such that within each lamella a statistically meaningful sampling of z_I (using the 10 ns trajectory) can be obtained and used to calculate the corresponding A_n(Z),

where P_n(Z) is the chloride ion Z coordinate probability distribution within lamella n.

In the simulations described below, the lamellae are 3 Å wide, overlap by 1 Å, and span the simulation box from the region of bulk water to bulk nitrobenzene. The ion is constrained within a specified lamella by a window potential (which is zero when the ion is inside that window but rises rapidly when the ion attempts to escape the window) through the course of the 10 ns trajectory. The series of A_n(Z) segments (n = 1, 2, …, N) is combined by using their overlapping regions^27,28 to arrive at the complete free energy profile for the entire interval of interest.

If A(Z) varies rapidly in a given region, most often due to the water finger “pulling” the ion toward the aqueous phase, poor statistical sampling may arise. To accelerate the exploration of phase space and improve sampling statistics in this region, we apply a biasing potential U_bias(Z) to the ion, modifying the total ion-solvent interaction energy,

where U_bias(Z) is a function of the ion Z position only and is typically of the form

where the constants α, γ, and ζ are selected (and given in supplementary material³⁹) so that the biasing potential approximates A_n(Z). The free energy profile in lamella with applied bias is therefore calculated as follows:²⁴ 

Quantitatively describing and tracking the transfer of a single ion from the aqueous phase to an adjacent immiscible phase has been challenging since early work.²⁹ The most obvious coordinate would be the ion’s distance from the liquid/liquid interface along an axis normal to the interface. Extracting the ion’s position in simulation Z-space is trivial, but assigning a reference Z-coordinate value to a liquid/liquid interface is not particularly clear due to fluctuations inherent to the liquid/liquid interface, including interfacial distortion caused by the water finger itself. More importantly, describing the state of the system using a single coordinate assumes that the two solvents provide a mean field for the motion of the ion and ignores solvent fluctuations. A relatively simple approach to take into account these fluctuations is to add an additional coordinate describing the interface. For example, as mentioned in the Introduction, Urbakh, Kornyshev and co-workers proposed a geometrical coordinate describing interface deformation and fluctuations. In this paper, we consider in detail two coordinates: an energy based one introduced by Schweighofer and Benjamin¹⁰ and a geometric one recently introduced by Morita and co-workers.¹⁸ Both of these coordinates can be computed for instantaneous system configurations and represent, in some sense to be explained below, a complimentary description of the environment experienced by the ion.

For a two-phase system, we define a solvation coordinate for each solvent,

where the subscripts W and NB represent the water or nitrobenzene solvents, I represents the ion, U is the total intermolecular interaction potential energy between the ion and the specified solvent, and U^bulk is the average total ion-solvent interaction potential when the ion is in the bulk solvent at the same temperature. Values of U^bulk (⁠$UI−Wbulk=−155.5$ kcal/mol and $UI−NBbulk=−89.39$ kcal/mol) are calculated by molecular dynamics simulations of a single ion in the uniform bulk liquids using the same force fields and parameters described above. Additional details regarding U^bulk calculations are provided in supplementary material.³⁹ The solvation coordinates consider both electrostatic and Lennard-Jones potentials between the ion and all solvent molecules present in the simulation. Note that s_W and s_NB are dimensionless and implicitly dependent on the ion’s location Z, in addition to the coordinates of all solvent molecules (although mainly sensitive to the solvent molecules in the nearest solvation shells). We expect the following: s_W(z → bulk water) = 1, s_W(z → bulk liquid) = f, where f > 1 represents the interaction of the ion with (mostly) the fraction of the hydration shell that was co-transferred to the organic phase (f ≈ 0 for a hydrophobic ion).

Recently, Kikkawa et al. employed mathematical graph theory to develop a new geometrical coordinate useful for describing ion transfer between water and an adjacent immiscible phase.¹⁸ To calculate the “water finger coordinate,” w, the ion and water molecules are treated as vertices in an undirected graph. The edges of the graph are the geometrical distances between the vertices. A connected path between the ion and bulk water is defined by the requirement that all edges along the path are shorter than a threshold distance. At each system configuration, the coordinate w is equal to the minimum threshold distance that will give rise to a connected path. For example, a “naked” ion with Z = + 10 Å in a system where the interface is constrained at Z = 0 Å would have w ≈ 10 Å. If a water molecule was present with the same X and Y coordinates as the ion and Z = 4 Å, w ≈ 6 Å. Our implementation of this algorithm varies slightly; we assign the position of the water oxygen molecule as the vertex representing each water molecule instead of the water molecule’s center of mass. This assignment allows us to directly compare the water finger coordinate to the commonly used geometrical definition of water-water hydrogen bonds, where r_OO < 3.4 Å.^30–32 A hydrated (or partially hydrated) ion “connected” to the bulk water phase via an intact water finger would have a value of w less than or equal to this upper bound of 3.4 Å, which we shall refer to as $rOOmax$⁠. For simplicity, we define the condition of “connected to bulk water” when the ion is connected to a water molecule whose oxygen molecule is located at Z < 0. This is sufficient since finger breakup always occurs at Z > 0.

It is useful to note at this point several features of this coordinate: When the ion is in bulk water, 〈w〉 ≈ peak position of O–O radial distribution function (or ion-oxygen RDF peak position if longer) regardless of the position Z_I of the ion. When the ion is in the organic phase connected via water finger to the aqueous phase, w ≈ O–O distance corresponding to the longest hydrogen bond. As the water finger breaks, w corresponds to the distance between the two nearest water molecules—one that belongs to the ion partial hydration shell and one to the water phase, regardless of the size of the ion’s hydration shell. If the ion in the organic phase is “naked,” w is equal to the distance between the ion and the closest water molecule.

To gain further insight into the role of solvent fluctuations during the ion transfer process, we also consider a system in which the water/nitrobenzene interface is constrained to be a molecularly sharp (“flat”) interface. This is done by adding a soft reflecting potential in the X-Y plane at Z = 0 that acts on any water or nitrobenzene molecule whose center of mass attempts to cross the Z = 0 plane. This external potential essentially restricts the water molecules to the Z < 0 region and the nitrobenzene molecules to the Z > 0 region and gives rise to a molecularly sharp interface, where capillary fluctuations are strongly dampened. The reflecting potential is a function of the distance from the center of mass of each molecule so that perturbations to the molecular orientational distribution are minimized. The implemented potential is of the form

where $UnW$ is the reflecting potential on water molecule n, whose center of mass is located at z_COM along the interface normal and H(x) is the Heaviside step function (H(x) = 0 for x < 0; H(x) = 1 for x > 0). Similarly, $UnNB$ is the reflecting potential on nitrobenzene molecule n. γ is a force constant whose exact value is not important. In the simulations described below we use γ = 1200 kcal mol⁻¹ Å⁻². The reflecting potential prevents the ion from carrying any part of its hydration shell across the interface and prevents the formation of the water “finger.”

Figure 1 shows the density profiles of water and nitrobenzene with the ion in bulk water, confined to a window −17.5 Å < z_I < − 14.5 Å. We chose Z = 0 to be the location of the Gibbs Dividing Surface (GDS, the plane parallel to interface where the water density is approximately 50% of the bulk value. For the exact definition see Ref. 33). The interface region is quite narrow (the distance over which the density of water varies from 90% to 10% of the bulk value is 5 Å). The properties of the neat water/nitrobenzene interface have been described elsewhere.²¹ The ion transfer equilibrium calculations described below are carried out with the ion location extending from the bulk water (Z = − 15 Å) to the bulk nitrobenzene regions (Z = 20 Å).

Figure 2 depicts several snapshots from the simulations, which illustrate some of the important structural considerations involved in the ion transfer process. In an unconstrained system (Figure 2(a)) as the chloride ion moves from the aqueous phase toward the nitrobenzene phase, it carries with it most of its hydration shell, which creates a local disturbance at the interface. Hydration shell water molecules remain hydrogen bonded to neighboring water molecules, and the ion remains tethered to the aqueous phase by the water “finger” structure. As the ion moves further into the organic phase, surface distortions are increased, the water “finger” elongates and finally breaks, typically when it has extended about 1 nm into the organic phase, leaving the ion with a partial hydration shell (Figure 2(b)). In contrast, when the interface is constrained to be flat, the ion is transferred as a single species into the organic phase (Figure 2(c)). Panel d is a snapshot from a simulation in which the chloride ion and three water molecules Cl⁻(H₂O)₃ are transferred across the flat interface, which is to be discussed below.

The free energy profiles corresponding to the different ion transfer processes examined here are shown in Figure 3(a). These free energy profiles were calculated as explained in Section II B. Additional details, including the values of the biasing potentials used, are given in supplementary material.³⁹ 

The free energy profile for the transfer of the chloride ion across the normal water/nitrobenzene interface has a shallow minimum at about Z = − 2 Å. From that point, it monotonically increases as the ion is transferred across the interface to bulk nitrobenzene. The net free energy of transfer is 12 ± 1 kcal/mol. We estimate the error in this calculation by considering the standard deviation of ΔA values from 10 free energy profiles generated from each of the 10 independent initial configurations. Statistical fluctuations are additive as the adjacent windows are “stitched together” to form a continuous free energy profile. This approach to error estimation is therefore quite conservative since the curve reported in Figure 3(a) averages data from all 10 configurations within each lamella before consolidating them into one continuous curve describing the ion transfer event, dampening statistical fluctuations between the independent transfer events. This result may be compared with the estimated experimental value³⁴ of 9.5 ± 0.5 kcal/mol. The calculated value of ΔA would presumably be lower and closer to experimental results with the implementation of a polarizable nitrobenzene force field. Fig. 3(b) depicts the corresponding average number of hydration shell water molecules (defined as the number of water molecules whose oxygen-chloride ion distance is less than the first minimum of the bulk O–Cl⁻ radial distribution function r_min = 3.85 Å). The change in the average hydration number as a function of Z tracks the change in the free energy profile. The decline (relative to the bulk value of 7.70 ± 0.05) begins at around Z = − 3 Å and reaches the value of 4.7 ± 0.6 when the ion is located within the range 15 Å < Z < 18 Å. Within this region the water finger fluctuation has broken and the ion and its remaining hydration shell move toward the bulk organic phase. The increasing fluctuation in hydration shell population is due to the breakup and formation of the water finger, which is most likely to happen when the ion is near Z ≈ 12 Å when the finger is fully extended, as will be discussed below.

When the interface is constrained to be molecularly sharp and the ion is transferred without any water molecules, the free energy profile is significantly different (red curve in Figure 3(a)): The shallow local minimum is gone, the rise in the free energy occurs at smaller values of Z compared with the normal interface and the net free energy of transfer is significantly larger (23 ± 1 kcal/mol). Clearly, transferring the ion into the organic phase is much more favorable when some number of water molecules accompany the ion. Since the initial increase in the free energy correlates quite closely with the reduced number of hydration shell water molecules, we can attribute the steeper rise in the free energy to the more rapid decrease in the hydration number (green line in Figure 3(b)), made necessary by the constraint imposed on the interface. A continuum viewpoint of this process can be described by invoking the model of a sphere moving toward a flat surface in a viscous fluid.¹⁵ In this model, the resistance deviates from Stokes’ formula when the sphere approaches the surface, increasing with 1/d, where d is the distance between the sphere and the surface. In the normal, unconstrained system this additional resistance results in the distortion of the interface, caused by the approaching ion. In the unconstrained system this resistance of the surface to the ion’s approach is considerably smaller than the approach toward a flat, solid surface since the unconstrained liquid/liquid interface may distort in response to the ion’s motion. Since the water molecules are unable to cross the Z = 0 plane, the ion must approach the interface as it would a solid surface (and then pass through it). The dramatically increased free energy requirement of the flat interface system in the region of −3 > Z > + 3 illustrates the energetic value of the surface distortion. These results demonstrate that, in the unconstrained (“normal”) simulation, interfacial distortions caused by the approach of the ion to the interface significantly reduce the free energy required for ion transfer.

With the flat interface constraint, we are able to gain some clues about the energetic benefit of interfacial fluctuations, surface distortion, and water finger formation during the ion transfer event. However, this constraint also completely strips the ion of its water solvation shell. Co-transfer of part of this hydration shell is a well-established phenomenon, witnessed in both experiment and simulation (for a recent review see Ref. 5) that should contribute toward the reduction of free energy in liquid/liquid ion transfer.³⁵ We can provide a partial account of this effect by considering the transfer of the ion with a fixed number n of water molecules across the constrained flat interface. In these simulations, we set n = 3, which is close to the experimentally accepted number.³⁶ Oxygen atoms of the 3 water molecules assigned to the solvation shell are attached to the chloride ion using the interaction potential given by $γHxx2;x=rCl−−rO−rmin$⁠, where the constant γ is 1200 kcal mol⁻¹ Å⁻², H is the Heaviside step function, and r_min = 3.85 Å is the location of the first minimum of the chloride-oxygen water radial distribution function. The three water molecules tethered to the ion are not influenced by the external potential that keeps the interface flat and therefore are allowed to cross into the organic phase with the ion. Thus, we are in fact considering the transfer of the stable species Cl⁻(H₂O)₃ across the flat interface. The addition of this small solvation shell reduces the free energy of ion transfer across the flat, constrained interface by 17% compared to the transfer of the fully “naked” ion. The shape of the free energy curves is different as well (Figure 3(a)): When 3 water molecules accompany the ion, the free energy plateau begins slightly further into the organic phase due to the presence of hydration shell water. The hydration shell allows the ion to remain within hydrogen bonding distance to the aqueous phase by the length of an additional water-water hydrogen bond. This additional water-water hydrogen bond extends beyond the potential barrier imposed at Z = 0. The second feature considers the small peak prior to the plateau, which we assign to the energetics of rearranging interfacial water molecules to interact with approaching or leaving ion. This peak is not readily present in the Cl⁻(H₂O)₃ simulation, which is a reasonable result since the small solvation shell water molecules would require less energy to reintegrate with the aqueous phase than the naked chloride ion.

Given the relatively long lifetime of water protrusions and other surface fluctuations and given the uncertainty and variability in the number of water molecules that are co-extracted with the ion, it is reasonable to question whether the free energy curves depicted in Fig. 3 describe true equilibrium behavior. To partially explore this issue and gain additional insight into the ion transfer process, we have also attempted to calculate the free energy profile for the back transfer of the ion from the nitrobenzene to the aqueous phase. The challenge here is to correctly select the initial state of the ion in the organic phase. One approach is to start from a configuration reached during the forward (water → organic) simulations.^12,13 In reality, an ensemble average over properly weighted different hydration clusters, Cl⁻(H₂O)_n, would be required to truly describe the reverse process. In Fig. 3, we show the free energy profile calculated for the transfer of Cl⁻(H₂O)₃ from the organic phase to the aqueous phase across the normal interface as the process with the most likely largest contribution. Indeed the free energy profile (dashed line Fig. 3(a)) is quite close to the curve for the normal transfer of Cl⁻ from water to nitrobenzene. Note that the average hydration number for this reverse transfer (dashed line Fig. 3(b)) increases as the Cl⁻(H₂O)₃ approaches the interface due to the water surface fluctuations. The increase begins near Z = 12 Å, which is the point identified earlier as the average maximum length of water fingers before breakup in the transfer from water to nitrobenzene.

These results suggest that the energetics of interfacial ion transfer (for small ions) is largely dictated by surface distortion and the formation of the water finger structure. As the finger extends to near maximum length, just prior to breaking, the structure is typically a single chain of hydrogen-bonded water molecules surrounded by molecules of the adjacent immiscible phase. Breaking the elongated water finger requires significantly less free energy than stripping the ion from the hydration shell as nitrobenzene molecules that surround the water finger interact and interfere with the increasingly fragile chain of hydrogen-bonded water. Simulations of the constrained surface do not allow the distortion of the surface of the random motion of nitrobenzene molecules to reduce the energetic cost of ion transfer. However when we allow a fixed number of water molecules to co-transfer with the water, we partially account for the energetic benefit of this co-transfer. Crossing the interface is quite similar in both flat-interface simulations, but the hydrated ion demonstrates an “easier” transition into the organic phase.

Although these PMF curves are calculated using an equilibrium umbrella sampling technique with a single ion, we may use these results to comment on the observable dynamics of a macroscopic system. However, it should be again noted that the physical constraints applied in the flat and Cl⁻(H₂O)₃ systems to deconstruct the energetics of this ion transfer do not have realistic experimental analogues. If equilibrated systems of these designs could exist, we should expect that they would have Boltzmann-weighted ion populations related to the free energy differences calculated above. Further, the flux of the chloride ions across the interface should be related to the free energy curves in Figure 3(a) by

where j_I is the flux across the interface, β is 1/(k_BT), and ΔA is the difference between the lowest and highest values in the respective PMF curve. With this approach, the flat interface constraint reduces ion flux by about a factor of 10⁷. When the 3-water solvation shell is allowed to co-transfer in the Cl⁻(H₂O)₃ simulation, the flux increases by about 2 orders of magnitude, about 10⁵ times slower than in the unconstrained system.

The description of the ion transfer process using only the ion’s position along the interface normal averages out the detailed mechanistic picture qualitatively described above. As the mechanism of ion transfer is intimately connected with the deformation of the interface, formation of the water “finger,” and the co-transfer of some hydration shell water molecules, additional coordinates to describe this information are required to develop a quantitative account of the transfer process. In general, for any generalized coordinate Γ(r) (such as s_W, s_NB, or w defined above) where r represents all nuclear positions (including the ion location x_I, y_I, z_I), the local equilibrium average 〈Γ〉(Z) is given by

where f(r) is the normalized phase space density (⁠$fr=exp−βUr/∫exp−βUrdr$⁠, U(r) is the total system potential energy, and β = 1/k_BT) and z_I is the ion’s location along the interface normal. The equilibrium simulations in each of the overlapping windows provide statistically reliable data to compute these ensemble averages, which we examine next.

The three panels in Figure 4 show the average values of the water solvation coordinate 〈s_W〉 and the nitrobenzene solvation coordinate 〈s_NB〉 as a function of the ion’s Z position for the three systems examined. Since most (about 70%) of the ion-water interaction energy is determined by the molecules in the first hydration shell, it is not surprising that the three 〈s_W〉(Z) curves in Figure 4(a) resemble the corresponding curves in Figure 3(b). Under normal conditions, 〈s_W〉 gradually decreases as the ion moves into the organic phase. The decrease is not dramatic since most ion-water interaction occurs between the ion and its solvation shell, which is mostly intact when the ion crosses the interface. The solvation coordinate continues to decrease as the ion moves into the organic phase and the average number of water molecules in the 1st hydration shell and contributions from the 2nd hydration shell and beyond decreases. In the flat-interface system, 〈s_W〉 decreases rapidly at the interface since the ion is stripped of its solvation shell water. The tail of the flat-interface curve represents the remaining ion-water long-range columbic attraction. The flat-interface, fixed hydration shell simulation (Cl⁻(H₂O)₃) represents an intermediate behavior. Since three water molecules are always present in the vicinity of the chloride ion, the water solvation coordinate plateaus once the ion leaves the short-range influence of the aqueous phase.

The nitrobenzene solvation coordinate is shown in Figure 4(b). In the normal, unconstrained system 〈s_NB〉 increases gradually as the ion moves to the organic phase. With the flat-interface constraint enabled, the ion is more exposed to the nitrobenzene phase as it approaches the interface since the constrained interface is unable to distort toward +Z. When the naked ion crosses into the organic phase, 〈s_NB〉 increases rapidly and plateaus. With three solvation-shell water molecules tethered to the ion in the flat-surface system (Cl⁻(H₂O)₃), the system plateaus at a lower value due to the presence of these solvation water molecules.

It is interesting to examine the sum 〈s_W〉 + 〈s_NB〉 for each of these systems. A fixed value of 1 for this sum, independent of Z, would suggest that the interaction energy of the ion with the water is replaced on a relative basis with the interaction energy of the ion with the less polar organic liquid. This was found to be the case for the transfer of tetra-alkyl ammonium ions and ion pairs across the water/chloroform interface.³⁷ Here, the small hydrophilic ion gives rise to a slightly different picture. While the sum does not vary much from 1 (only about 5%-15% total change depending on the system), the variations provide additional insight into the structure of the hydrated ion that is consistent with the discussion above. Thus, for the normal interface, 〈s_W〉 + 〈s_NB〉 varies from 1 in bulk water to 1.15 in bulk nitrobenzene. The slight increase suggests that the decrease in the total ion-water interactions is slightly made up for by the fact that the water molecules in the first hydration shell are held more tightly together when the ion is in bulk nitrobenzene than when the ion is in bulk water. This is also confirmed by noting the increased water-ion interaction energy per water molecule in the hydration shell and the longer water-ion residence time when the hydrated ion is in the bulk organic phase.³⁸ The slight increase in 〈s_W〉 + 〈s_NB〉 as the ion is transferred from water across the flat interface is initially the same as that of the normal interface until the ion reaches the Gibbs surface. At this point, the rapid drop in the ion-water interaction is seen as a sudden decrease in the sum of the solvation coordinates as the ion is stripped of its hydration shell. When the ion with a fixed hydration shell (Cl⁻(H₂O)₃) crosses the interface (Figure 4(c), green curve), the interaction energy between the ion and the three tethered water molecules increases, similar to the behavior of the unconstrained system’s hydration shell. The sum of the solvation coordinates approaches 1.1 as the ion and its fixed n = 3 shell move into the bulk nitrobenzene region.

With this understanding of the nature of the solvation coordinates s_W and s_NB, it is useful to consider the dependence of the system’s free energy on them. In general, the two-dimensional free energy A(Z, s) in which Z and s (where s stands for s_W or s_NB) vary independently, would contain a relatively complete description of the ion transfer. In practice, computing A(Z, s) would be challenging because values of s significantly different from the equilibrium values are extremely rare and irrelevant to the process of electrochemical ion transfer. When the ion is located in the aqueous phase, the probability of observing a value of s_W significantly different from 1 (or s_NB different from 0) is essentially zero, corresponding to a very high value of the free energy. When the ion is at or near the GDS, relatively large fluctuations in the values of s_W and s_NB are more likely, corresponding to lower free energy values. As the ion is moved to the organic phase, at each value of Z the two-dimensional free energy A(Z, s) has an approximate parabolic dependence on s_W or s_NB with the minimum (equilibrium value) of A attained at 〈s_W〉 or 〈s_NB〉. This equilibrium value is plotted in Figure 5 as a function of 〈s_W〉 (Figure 5(a)) and 〈s_NB〉 (Figure 5(b)). The fact that the top panel (where s_W varies from 1 to 0) and the bottom panel (where s_NB varies from 0 to 1) are almost identical is a direct result of the relationship 〈s_W〉 + 〈s_NB〉 ≈ 1 mentioned above. They describe the free energy change when the ion is transferred from water to nitrobenzene when the solvation coordinates are allowed to reach equilibrium at each value of Z. Unlike the significantly different behavior depicted in Figure 4 for the three systems, the similarity in the shape of the three curves in Figure 5 is quite striking, with the main difference being the final value of the free energy. This suggests that the local equilibrium free energy of the ion is more strongly correlated to its solvation state than to its location. The relatively small differences between the three curves in the region 0.5 < s_W < 0.7 (more clearly observed in the top panel) suggest that significant contribution to the free energy of transfer is the deformation of the interface and the creation of the water finger, but that this is less than the “cost” of transferring a “naked” ion.

While Figure 5 suggests that solvation coordinate is highly correlated with the position of the ion, it is interesting to examine in some detail the parabolic dependence of the free energy A(Z, s) on the solvation coordinate s. Specifically, for small deviations from equilibrium,

Figure 6 shows the values of the “solvent force constants” k_W and k_NB determined from a parabolic fit of the fluctuations in the solvent coordinates about their equilibrium values. A large value of the force constant indicates a tighter solvation environment, which requires larger free energy to deform. An example of this parabolic relation is shown in the inset of Figure 6. Here A(Z, s_W) of an ion confined to a window defined by −1.5 Å < Z < 1.5 Å and the parabolic fit from which k_W is obtained (dashed curve) are shown. As the ion is transferred from the aqueous phase to the organic phase the decrease in the k_W force constant corresponds to the large fluctuations in the hydration energy caused by the depletion of the hydration shell and the formation of the water finger. In contrast, k_NB is mostly constant when the ion is in the organic phase, with a rapid increase as the ion sheds the solvating nitrobenzene molecules upon entering the interface region and moving toward the aqueous phase. We should point out that when the ion is in the organic phase and the water finger is broken, several hydration shell structures as well as different finger structures may contribute to the fluctuations in s_W and these fluctuations are no longer well described by a single parabola.

We next consider the variation in the water finger coordinate w as the ion is transferred through the interface. Figure 7 shows the average 〈w〉 as a function of the chloride ion’s Z position, calculated by averaging the values of w in 0.1 Å bins along Z. As long as the ion is attached to the bulk water phase via an unbroken network of water molecules, 〈w〉 is near the O–O distance corresponding to a water-water hydrogen bond (about r_OO = 2.75 Å, which corresponds to the peak value of the bulk water O–O radial distribution function). The value of Z at which 〈w〉 begins to significantly increase beyond $rOOmax$ (3.4 Å, see Section II C 2) is the point where the ion breaks away from the aqueous phase. In the constrained surface systems, this separation from the bulk water begins at Z ≈ 3 Å for the transfer of the naked ion and at Z ≈ 5 Å for the transfer of Cl⁻(H₂O)₃ across the flat interface. The difference of about 2 Å between the two constrained systems is attributed to an additional hydration water between the ion and the aqueous phase during the transfer of Cl⁻(H₂O)₃. When the ion breaks free of the interface, the last interacting water remaining at the interface is no longer pulled toward the +Z direction, and it recedes into the aqueous phase, resulting in the rapid increase in 〈w〉 seen in both of the constrained simulations (Z ≈ 4 Å and Z ≈ 6 Å). As the ion moves away from the flat constrained surface, 〈w〉 increases linearly with the ion’s Z position and in this region, to a good approximation, 〈w〉 = Z for the naked ion and 〈w〉 = Z − Z_OCl for the transfer of Cl⁻(H₂O)₃, where Z_OCl is the projection along the interface normal of the position of the peak Cl–O radial distribution function, about 3.2 Å. The average shortening of 〈w〉, when compared to the naked ion at the same Z position, is approximately 2.5 Å due to the intervening of slightly less than a single water molecule (on average) between the ion and the flat interface.

In the unconstrained, normal system, the water finger appears to break when the ion is approximately at Z = 12 Å. Examination of individual trajectories (see below) suggests that the breakup follows elongation of the water finger and significant water surface fluctuations, which make the initial rise not as sharp as in the constrained systems. As the water finger elongates, the water finger structure changes from a more conical morphology to a linear chain of water molecules, finally resulting in the breakup of a hydrogen bond between a hydration shell water molecules and the other water molecules in the finger.

Additional insight into the utility of the water finger coordinate for describing the mechanism of ion transfer is provided in Figure 8, which shows the local equilibrium free energy of the ion plotted against the equilibrium value of the water finger coordinate. The nearly step-like behavior of the three systems suggests that, in all systems examined, by the time the last hydrogen bond is broken (when 〈w〉 ≈ 3.4 Å, the commonly used maximum oxygen-oxygen distance when geometrically defining the existence of water-water hydrogen bonds), most of the free energy change has been accounted for, especially when surface fluctuations are suppressed. This is more clearly observed in the normalized free energy change in the bottom panel, where ΔA_t is the net free energy of transfer from water to nitrobenzene. This behavior is consistent with our earlier discussion, which assigns most of the free energy change to the creation and elongation of the water finger in the normal system or the stripping of the water hydration shell in the constrained systems. The fact that 〈w〉 does not appreciably change during this critical phase of the finger formation and elongation reduces its utility for understanding the mechanism of the transfer (except for the obvious observation that transfer via formation of a water finger is much more energetically favorable).

Because the water finger coordinate w is (loosely speaking) a measure of the longest oxygen-oxygen distance, it may be used to monitor the integrity of the water finger. This is demonstrated in Figure 9, which depicts Z, w, s_W and s_NB vs. time for a selected 0.8 ns trajectory of the unconstrained system. In this particular simulation, the ion is constrained within the window 8.5 Å > Z > 11.5 Å (in the organic phase). The relatively flat line near w ≈ 3 Å indicates the region with an intact water finger, where the ion is tethered to the aqueous phase by at least a chain of hydrogen bonded water molecules with a maximum r_OO of approximately 3.4 Å. At 0.27 ns the water finger breaks, w increases to over 10 Å, and remains broken with rapid fluctuations for the next 0.3 ns until it reforms at t = 0.62 ns and remains intact for the rest of the trajectory. This spontaneous formation and breakup of the water finger while the ion is held above the interface quantitatively demonstrate the harpoon or “protrusion”-based transfer mechanisms suggested in earlier work.^9,15

Figure 9 also reveals the statistical fluctuations of the solvation and water finger coordinates during the most erratic part of the transfer process: separation of the ion from the water phase. The solvation coordinates s_W and s_NB have a fairly consistent standard deviation of 0.1 in the MD trajectory observed in Figure 9. The water finger demonstrates a significantly different behavior in regard to fluctuation, with a standard deviation of 0.07 Å when the water finger structure is intact (0.0 ns < t < 0.2 ns) and 1 Å when the water finger has broken (0.3 ns < t < 0.6 ns). This difference in σ_w values reflects the nature of the coordinate. When the water finger is intact, fluctuations in w are essentially fluctuations in hydrogen bond distances. When broken, fluctuations in w largely describe geometric fluctuations of the liquid-liquid interface.

Several interesting points can be gleaned from Figure 9. When the finger is intact (t < 0.27 or t > 0.62 ns), the ion is able to explore the whole range of Z values, suggesting that the water finger is able to keep up and fluctuate accordingly. When the finger is broken, w experiences large fluctuations while the ion is mostly near the “upper” (Z ≈ 11 Å) edge of the window. Without the anchor provided by the water finger, the ion is more likely to drift towards the bulk of the organic phase due to the applied biasing potential. This is also reflected by the variations observed in the values of the solvation coordinates (bottom panel): s_NB increases and s_W decreases during this time frame (while the sum fluctuates around a fixed average of 1.1). The large fluctuations in w correspond to large surface water deformation. For example at t ≈ 0.4 ns we see that w ≈ Z ≈ 11 Å corresponding to a situation where the ion’s partial hydration shell points away from the interface while the interface is locally flat.

Characterizing liquid surface fluctuations during ion transfer across the interface between water and an immiscible organic solvent and taking into account the partial co-transfer of the ion hydration shell are crucial for correctly describing the mechanism of this important process. Our free energy calculations for the transfer of Cl⁻ across the water/nitrobenzene interface show that a major contribution to the ion’s free energy of transfer from the aqueous to the organic phase is in the deformation of the interface and the creation of water molecules protrusions (“finger”) into the organic phase. However, a significant “savings” in free energy is gained by the co-transfer of several water molecules with the ion. The ion solvation energy with each of the two liquids (normalized by the corresponding bulk values) gives crucial information about the local thermodynamic state of the ion and is strongly correlated with the ion local free energy and position. In contrast, the geometric coordinate w does not appreciably change during the critical phase of the water finger formation and elongation. This coordinate is useful for monitoring the integrity of the water “finger” and for identifying the onset of the “harpoon”-like mechanism for the ion transfer.

We are grateful for several illuminating discussions with professor Morita. J.K. also thanks L. Wang, N. Kikkawa, and the Morita group for their hospitality. Financial support from the National Science Foundation through Grant No. CHE-1363076 is acknowledged.