Hydration structure and dynamics of phosphoric acid and its anions—Ultrafast 2D-IR spectroscopy and ab initio molecular dynamics simulations

Phosphate groups are key structural units and hydration sites of nucleotides, such as adenosine di-and tri-phosphate, the backbone of DNA and RNA, and/or the head groups in phospholipid cellular membranes.^1–3 In biochemistry, numerous metabolic and signaling processes rely on phosphorylation reactions for storing and releasing energy.^4,5 Beyond their fundamental relevance in chemistry and biology, phosphates have found widespread applications in pH buffering of solutions, materials treatment, fertilizers, or food processing. Monomeric phosphates exist as the charge neutral phosphoric acid H₃PO₄ or anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻, depending on the pH value of an aqueous environment. The different species display markedly different hydration structures, which are characterized by a pattern of hydrogen bonds with water molecules in the first solvation shell and electric interactions with solvent molecules beyond the first shell.

Phosphate hydration has been studied in detail by linear vibrational and dielectric spectroscopy.^6–10 In parallel, extensive theoretical work based on density functional theory, quantum mechanics/molecular mechanics (QM/MM) methods, and molecular dynamics simulations has addressed the basic molecular interactions and spatial arrangements of solute and solvent molecules.^7,11–17 Such work has mainly focused on time-averaged structures and interactions in thermal equilibrium. Hydration dynamics on the intrinsic molecular length and time scales has been addressed for a variety of ions by nonlinear femtosecond infrared spectroscopy.^18–23 Recently, such methods have been applied to elucidate structural fluctuations around phosphates and the related fluctuating electric fields on their intrinsic ultrashort time scale.^24–27 In a most direct approach, phosphate vibrations serve as specific probes for mapping local interactions and the structure and dynamics of the water shell.

The asymmetric phosphate stretching vibration ν_AS(PO₂) mainly involves displacement of the PO₂ unit of the phosphate group and has been used as a local probe in ultrafast infrared spectroscopy²⁸ to study model systems with hydrated PO₂ groups, e.g., (CH₃)₂PO₄ or KH₂PO₄, and larger DNA and RNA backbone structures. Two-dimensional infrared (2D-IR) spectroscopy has given detailed insight into vibrational line shapes, which exhibit a moderate inhomogeneous broadening due to site variations in hydration and line broadening by vibrational dephasing on a 300 fs time scale. Moreover, the formation of solvent-separated and contact ion pairs of PO₂ with Na⁺, Ca²⁺, and, particularly Mg²⁺ cations has been detected via characteristic frequency shifts and line shape changes.^29–31 The theoretical analysis of such results by DFT cluster calculations and QM/MM simulations of ultrafast molecular dynamics has demonstrated a moderate slowing-down of water dynamics in the first few hydration shells compared to bulk water and fluctuating electric fields with amplitudes of several ten megavolts/cm from the aqueous environment. The first hydration shell around the PO₂ unit consists of up to 6 water molecules forming hydrogen bonds with the two oxygen atoms. In contact ion pairs, repulsive short-range interactions of PO₂ with the cation lead to a spectral blue shift of the ν_AS(PO₂) absorption band.

There is limited and partly conflicting information on the hydration structure and dynamics of phosphoric acid H₃PO₄ and its anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻.^{7,8,10,12,13,16} So far, the ionic species have mainly been studied in aqueous solutions of sodium and potassium phosphate salts at high concentration, where, beyond separately solvated ions, cation/anion pairs are formed and enhance the structural heterogeneity. Such behavior becomes increasingly complex at high molar concentrations of phosphate ions, with a direct impact on vibrational⁸ and dielectric spectra.¹⁰ 

H₃PO₄, H₂PO₄⁻, and HPO₄²⁻ can act as hydrogen bond acceptors and hydrogen bond donors, while PO₄³⁻ acts as hydrogen bond acceptor only. The impact of such dual hydrogen bond functionality on hydration dynamics is largely unexplored and insight from femtosecond nonlinear vibrational spectroscopy has remained very limited. In this article, we present a systematic 2D-IR study of hydrated phosphoric acid and its anions in the frequency range of the symmetric and asymmetric phosphate stretching vibrations. The relative abundance of the different species is tuned by changing the pH value of diluted aqueous solutions, thus avoiding the issue of ion pair formation. In the samples, there are ∼280 water molecules per phosphate, in order to suppress interactions between solvated ions and provide a sufficient fraction of bulk-like water. The fluctuating structure of the water environment around the phosphates gives rise to spectral diffusion of vibrational excitations, which is mapped onto the 2D-IR line shapes of symmetric and asymmetric phosphate stretching vibrations. The 2D-IR spectra demonstrate a pronounced slowing-down of water shell dynamics when going from H₃PO₄ to PO₄³⁻. The experimental results are complemented by ab initio molecular dynamics simulations of hydration geometries and dynamics for H₃PO₄ and its anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ in water.

Solutions of phosphoric acid H₃PO₄ (Thermo Fisher) in ultrapure water (Rotipuran Ultra, Roth) were prepared with a H₃PO₄ concentration of c = 0.2M. At this concentration, there are some 280 water molecules per PO₄ unit. The large water excess allows for safely suppressing interactions between PO₄ moieties at different solvation sites.

The concentration of Na⁺ ions at different pH values is presented at the bottom line of Table I. For pH = 1.0, 9.7, and 13.2, the species H₃PO₄, HPO₄²⁻, and PO₄³⁻ are predominant, whereas appreciable amounts of both H₃PO₄ and H₂PO₄⁻ exist at pH = 2.3. At pH = 13.2, there are similar concentrations of PO₄³⁻ and Na⁺ ions, which could result in the formation of different ion pair species, contributing to the spectroscopic observables. The results on PO₂²⁻/Na⁺ contact ion pairs reported in Ref. 30 suggest that the fraction of such ion pairs and their contributions to the linear infrared absorption and 2D-IR spectra is minor for equimolar solutions with c ≈ 0.2M.

Linear infrared absorption spectra were recorded with a Fourier transform infrared spectrometer (Bruker Vertex 80, spectral resolution 2 cm⁻¹). The liquid samples were held in a demountable cell (Harrick) with two 1 mm thick CaF₂ windows separated by a Teflon spacer of 12 or 25 µm thickness. The infrared absorption of neat water was measured under the same conditions and subtracted from the spectra of the solutions (cf. Fig. 1).

2D-IR spectra were measured in 3-pulse photon-echo experiments with heterodyne detection of the nonlinear signal. Infrared pulses with a center frequency of ∼1100 cm⁻¹, a spectral width of ∼190 cm⁻¹ (intensity FWHM), a duration of 120 fs, and a pulse energy of 5 µJ are generated by optical parametric amplification with a 1 kHz repetition rate.³² Two passively phase-stabilized pulse pairs with wavevectors (k₁, k₂) (mutual delay τ, coherence time) and (k₃, k_LO) are generated by reflection from a diffractive optic element.

In a box-CARS beam geometry, three pulses are focused into the sample by an off-axis parabolic mirror to generate a photon echo signal. The third-order signal emitted in the phase-matching direction −k₁ + k₂ + k₃ = k_Sig is overlapped with a local oscillator pulse (k_LO) for spectrally resolved detection by a 64 pixel mercury cadmium telluride detector array (frequency ν₃, spectral resolution 2 cm⁻¹). Signals measured for different coherence times τ are Fourier-transformed along τ to generate the excitation frequency coordinate ν₁. In Fig. 3, the absorptive 2D-IR signal, which is given by the real part of the sum of the rephasing and nonrephasing signals, is plotted as a function of ν₁ and ν₃. Femtosecond pump–probe experiments were performed in the 2D-IR setup with pulse 1 as the pump and the attenuated pulse 3 as the probe pulse.

Ab initio (first principles) molecular dynamics (AIMD) simulations of H₃PO₄, H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ in water were performed with the CP2K program package (version 8.2).³³ AIMD simulations start from a pre-equilibrated box of 64 water molecules (TIP5P force field, cell dimensions 12.4138 × 12.4138 × 12.4138 ų; the box is taken from Ref. 34) by replacing one water molecule at the center of the simulation box by a H₃PO₄ molecule, followed by relaxation of the simulation box. Trajectories of the ions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ in water were generated by successively removing a proton after every ∼5 ps of simulation time. After an initial 5 ps equilibration period of every trajectory that was discarded, AIMD simulations cover a simulation time of 80 ps for H₃PO₄ and HPO₄²⁻ and 77 ps for H₂PO₄⁻ and PO₄³⁻ (total simulation time: 314 ps).

AIMD simulations were performed on the hybrid density functional level of theory using the PW6B95 meta-GGA functional³⁵ that uses 28% of exact exchange in the exchange correlation functional. The PW6B95 functional has shown an excellent performance for the thermodynamics and structure of main group elements³⁶ and yields the radial distribution function of liquid water with excellent accuracy [cf. inlay Fig. 5(e)].

AIMD simulations used a high energy cutoff of 280 a.u. and a relative cutoff of 30 a.u. The screening threshold EPS_SCHWARZ for exact exchange was set to 1.0 × 10⁻⁶, and the interaction potential was truncated at 4.5 Å. The TZV2P-GTH basis set was used for hydrogen, oxygen, and phosphorous atoms together with the auxiliary cFit3 basis and the PBE pseudopotential was used for core electrons. Simulations were performed in the NVT ensemble at 300 K using a simulation time step of 0.5 fs and a generalized Langevin equation thermostat (ndim = 5, relaxation time scale 1.0 ps⁻¹). All the AIMD simulations were performed on the CooLMUC-2 Linux Cluster of the Leibniz Supercomputing Centre.

For analysis of the simulation results, radial distribution functions were calculated with the VMD program³⁷ that was also used to produce graphical representations of the solvation shell. Electric fields imposed by the water solvent were evaluated for H₃PO₄, H₂PO₄, HPO₄²⁻, and PO₄³⁻ after mapping the real space trajectory to the initial simulation box, centered around the phosphorus atom. The simulation of electric field employs fixed charges of the TIP3P-FB water model,³⁸ thereby neglecting polarization effects in the solvation shell due to interactions with the solute. Simulation with a polarizable model¹⁵ has shown that induction contributions to the electric field are a minor contribution (10%–15%) in the solvation forces of phosphates. Moreover, it was demonstrated that the pronounced librational fluctuation dynamics are dominated by the electrostatic contribution to the correlation function. Electric fields have been evaluated at the midpoint of P=O (H₃PO₄) and P–O_H bonds (H₃PO₄, H₂PO₄⁻, HPO₄²⁻), followed by projection along the bond axis. For PO₂ groups of H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻, electric fields have been evaluated at the bisector midpoint of the PO₂ units, followed by projection along the PO₂ bisector axis. This procedure has been shown to quantitatively explain the spectral properties of the asymmetric phosphate stretching vibration ν_AS(PO₂) via the vibrational Stark effect.¹⁵ 

Infrared absorption spectra of the different phosphate species are shown in Fig. 1. The absorbance A = −log(T) (T: sample transmission) measured with a sample thickness of 25 µm is plotted as a function of wavenumber between 950 and 1250 cm⁻¹ (solid lines). The dashed line in panel (b) gives the absorption spectrum of KH₂PO₄ dissolved in water for a concentration of c = 1.0M, at which KH₂PO₄ fully dissociates in hydrated K⁺ and H₂PO₄⁻ ions.²⁶ This spectrum is normalized to the peak absorbance at 1077 cm⁻¹ of the H₃PO₄ solution at pH = 2.3.

The predominant phosphate species at pH = 1.0, 9.7, and 13.2 is H₃PO₄, HPO₄²⁻, and PO₄³⁻, respectively (cf. Table I). Their absorption spectra are shown in panels (a), (c), and (d) of Fig. 1. At pH = 2.3 [Fig. 1(b)], there is a mixture of 72.5% H₂PO₄⁻ and 27.5% H₃PO₄. The corresponding infrared spectrum displays strong bands of H₂PO₄⁻ at 1077 and 1160 cm⁻¹ and weak bands of H₃PO₄ at 1010 cm⁻¹ and between 1130 and 1210 cm⁻¹. The spectrum of the KH₂PO₄ sample [dashed line in panel (b)] exhibits the H₂PO₄⁻ bands only. The asymmetric and symmetric P–O_H vibrations of H₂PO₂⁻ are centered at 940 and 880 cm⁻¹, i.e., outside the frequency range shown in Fig. 1(b) and have been analyzed in detail in Ref. 27. The absorption bands of the symmetric and asymmetric phosphate stretching vibrations ν_S,AS are marked in Fig. 1 and the respective frequencies are presented in Table II. The symmetric stretching mode ν_S of H₃PO₄ shows a very small transition dipole but is clearly identified in the Raman spectrum.⁸ For the phosphate species H₃PO₄, H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻, we observe a successive redshift from 1170 cm⁻¹ [ν(P = O)] to 1008 cm⁻¹ [ν_AS(PO₄)], reflecting the weakening of the bonding potential upon the increase in negative charge density. In the present study, we focus on the asymmetric phosphate stretching vibrations ν_AS as probes of hydration dynamics.

The symmetric and asymmetric stretching vibrations of PO₂⁻ groups in (CH₃)₂PO₄, KH₂PO₄, DNA, and RNA embedded in an aqueous environment display lifetimes between 250 and 400 fs, as has been discussed in Refs. 26 and 28. The fast population decays of phosphate stretching vibrations and their independence from the particular molecular species and environment point to a relaxation via a population transfer to PO₄ vibrations at lower frequencies, which are anharmonically coupled to the stretching vibrations.

The 2D-IR spectra of the different phosphate species are shown in Fig. 3. The spectra were measured with phosphoric acid samples at different pH values, except for panel (j), which displays the 2D-IR spectrum of a KH₂PO₄ sample.²⁶ In panels (e)–(l), the absorptive 2D-IR signal is plotted as a function of excitation frequency ν₁ and detection frequency ν₃ for different waiting times T. At the early waiting time T = 300 fs [panels (e)–(h), (j)], there is a negligible temporal overlap of the three femtosecond pulses in the photon-echo experiment. The range toward longer waiting times is limited to some 600 fs [panels (i), (k), (l)] because of the femtosecond population decay times of the vibrational excitations. The 2D-IR spectra display diagonal peaks at detection frequencies around the frequency positions of all the vibrational absorption bands in the linear infrared absorption spectra [panels (a)–(d); please note the different ν₃ ranges of the spectra of different phosphate species (columns)].

The 2D-IR line shapes of the v = 0 → 1 signals (yellow-red peaks) display elliptical contours with a different tilting angle relative to the frequency diagonal ν₁ = ν₃. To quantify the orientation of the elliptic peaks, we derive center lines (CL), which connect the frequency positions of the maximum signals in (horizontal) cuts along the detection frequency axis ν₃ at a fixed excitation frequency ν₁.⁴⁰ The thick blue lines shown in Fig. 3 give the CL for the asymmetric stretching modes ν_AS. The inverse slope of such lines, the centerline slope (CLS), is a measure for the correlation of excitation (ν₁) and detection (ν₃) frequencies. A CLS = 1 corresponds to CL parallel to the ν₁ = ν₃ frequency diagonal and, thus a perfect correlation of ν₁ and ν₃, such as in the case of static inhomogeneous broadening. In contrast, a CLS = 0 corresponds to vertical CL (parallel to the ν₁ axis) and total loss of correlation, such as induced by spectral diffusion.

The CL of the H₃PO₄ spectra shown in Figs. 3(e) and 3(i) are vertical within the experimental accuracy, corresponding to CLS≈0. This observation points to a very fast spectral diffusion within the initial time window of up to T = 300 fs. The 2D-IR peak of H₂PO₄⁻ at (ν₁, ν₃)=(1156, 1159) cm⁻¹ [T = 300 fs, Fig. 3(f)] is superimposed by a broad component originating from H₃PO₄, which affects the CLS. Thus, we analyze the 2D-IR spectrum of KH₂PO₄ [Fig. 3(j)], displaying a small CLS ≈ 0.25 ± 0.05/0.1 of the ν_AS peak and again pointing to pronounced spectral diffusion within the first 300 fs. The results presented in Ref. 26 suggest a predominant 50 fs decay of the frequency-fluctuation correlation function, which governs the experimental 2D-IR line shape (Fig. 6(b) of Ref. 26).

The CLS for HPO₄²⁻ [Figs. 3(g) and 3(k)] has a value of 0.36 ± 0.05, whereas the PO₄³⁻ peaks at T = 300 and 500 fs display a CLS = 0.62 ± 0.05. The frequency cuts of the PO₄³⁻ spectra shown in Fig. 4(d) exhibit a similar width of the positive 2D-IR signal at T = 300 and 500 fs, pointing to a limited spectral diffusion up to 500 fs. A similar behavior is borne out by cuts along the diagonal ν₁ = ν₃ for the two waiting times (not shown). The CLS for the different species suggest that the loss of frequency correlation becomes slower with a decreasing number of hydrogen atoms attached to the PO₄ unit, which corresponds to a decreasing number of OH donor groups from the phosphate in hydrogen bonds with the first water solvation shell. Moreover, the Coulomb interaction between the phosphate and the water shell becomes stronger with increasing ion charge. Such two mechanisms appear to make the water environment more rigid and to slow down spectral diffusion via the fluctuating electric force from the solvent.

The time averaged hydration structure of phosphoric acid H₃PO₄ and its anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ obtained from high-level AIMD simulations is shown in Fig. 5. The solvation shell is quantified via the phosphorous–water–oxygen (P…O_w) radial distribution function [rdf, Fig. 5(e)]. Representative snapshots of first hydration shell geometries are shown in Figs. 5(a)–5(d). For H₃PO₄, the average simulated phosphorous-oxygen bond lengths are 1.491 Å (P=O) and 1.572 Å (P–O_H), in excellent agreement with bond lengths derived from EXAFS experiments⁴¹ (1.49 and 1.54 Å, respectively). Similar good agreement is found for all ions (e.g., PO₄³⁻: mean P–O distance 1.558 Å; exp: 1.54 Å), suggesting a high quality of AIMD simulations of phosphates in liquid water with the PW6B95 functional.

For H₃PO₄, the P…O_w rdf shows a moderately structured first solvation shell (P…O_w distance <4.0 Å) and negligible long-range structure beyond the first hydration shell (P…O_w distances >4.0 Å). The first solvation shell is made up of 4.7 water molecules (integrated first maximum of P…O_w rdf). For the ions H₂PO₄⁻ and HPO₄²⁻, the solvation shell is strikingly similar and substantially more structured than for H₃PO₄. The first hydration layer is composed of 8 and 9 water molecules, respectively (integrated first maximum of rdf P…O_w until 4.5 Å), and a pronounced second maximum appears in the rdf. For PO₄³⁻, the first hydration layer is composed of 13 water molecules (integrated first maximum of rdf P…O_w until 4.5 Å). As reflected in the different height of the first peak of the rdf of HPO₄²⁻ and PO₄³⁻, the first solvation shell of PO₄³⁻ contains a substantially larger number of water molecules. On average, each of the P=O units is surrounded by about 3 water molecules in a tetrahedral arrangement. Apparent is a long-range structure in the hydration shell around PO₄³⁻ with a pronounced second minimum and recognizable third maximum in the P…O_w rdf (P…O_w distance ∼8.2 Å, not shown). Comparison of the first maxima of the P…O_w rdf yields that maxima are shifted to slightly larger P…O_w distances in the order H₃PO₄, H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ (3.68, 3.71, 3.76, and 3.78 Å, respectively), presumably due to increasing steric hindrance with increasing water content of the first hydration shell.

Comparing H₃PO₄ and its anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻, we find increasingly long-range structured hydration shells around the phosphates with increasing Coulomb interactions, i.e., phosphate ion charge (Fig. 5), which also increases the structuring in O_w…O_w rdf in comparison with neat water. Such increased water structuring correlates with the structure making properties of the highly charged ions HPO₄²⁻ and PO₄³⁻ in the Hofmeister series.

The rdf of hydroxy oxygen–water oxygen atoms [PO_H…O_w rdf, Fig. 5(f)] characterizes the water population at the hydroxy group and the hydrogen bond donor properties of H₃PO₄, H₂PO₄⁻, and HPO₄²⁻. For H₃PO₄, we find particularly strong hydrogen bonds with the first hydration shell [first maximum of PO_H…O_w rdf = 2.56 Å]. For H₂PO₄⁻ and HPO₄²⁻, the first maxima of the PO_H…O_w rdf are shifted to larger distances (2.68 and 2.76 Å) that become comparable in strength with the water–water hydrogen bond [2.78 Å, inlay Fig. 5(e)]. The population of the first hydration shell is lower for H₂PO₄⁻ than for HPO₄²⁻.

The rdf of phosphate oxygen–water oxygen atoms (P=O…O_w rdf) were calculated to quantify the hydrogen bond strength of P=O groups and water molecules [Fig. 5(g)]. We find an increasing hydrogen bond strength in the order H₃PO₄, H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ with respective P=O⋯OH₂ distances of 2.76, 2.74, 2.72, and 2.66 Å (first maximum of P=O…O_w rdf). On average, the P=O unit of H₃PO₄ is populated by about 2.1 water molecules (integrated first maximum of rdf P=O…O_w until 3.0 Å), while the P=O units of the ions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ are populated by 2.5–2.75 water molecules.

Compared to previous ab initio MD simulations,¹⁶ a very similar hydration structure and comparable number of water molecules is found in the first solvation shell of the ions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻, but a substantially reduced water population of the P=O unit of H₃PO₄ is observed in the current study. A possible origin of such discrepancies is given by the current higher level, hybrid density functional treatment, compared to a local GGA functional in Ref. 16. The hybrid density functional level of theory is found to provide a substantially improved water structure [inlay, Fig. 5(e)] and P=O bond length in H₃PO₄ [Figs. 6(a) and 6(b)]. Such balanced description of interactions is expected to be particularly relevant for describing the hydration shell of H₃PO₄, thus, the presented simulations are considered more reliable.

While PO₄³⁻ has tetrahedral symmetry in the gas phase, instantaneous geometries in a fluctuating liquid are subject to symmetry distortions due to stochastic interactions with the solvation shell. Figure 6 shows the fluctuations of bond length of the POx groups (x = 1, ..., 4) of H₃PO₄ and the ions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻. Figures 6(b), 6(d), 6(f), and 6(h) show a shift of the bond length distributions to longer phosphorous oxygen distances with increasing negative charge density, reflecting the gradual weakening of bond order from a P=O to a P–O⁻ unit. For a particular solvent configuration, the oxygen atoms of the POx groups are not equivalent but differ due to the local environment. Nevertheless, the time evolution of the two PO₂ bond length of H₂PO₄⁻ [Fig. 6(c)] is strongly (anti-)correlated, reflecting the superimposed motion of symmetric and asymmetric PO₂ stretching vibrations. For HPO₄²⁻ and PO₄³⁻ [Figs. 6(e) and 6(g)], the correlated time evolution of two oxygens is preserved on short time scales (∼200 fs), while the bond length of the additional oxygen atoms of the POx unit is decoupled from the motion of the PO₂ unit, as reflected in the larger amplitude bond length excursions of a single alternating bond phosphorous oxygen bond length. In other words, the symmetry of degenerate PO₂ units is broken in HPO₄²⁻ and PO₄³⁻, leading to decoupling of P–O and PO₂ units. The latter are randomized by the fluctuating forces arising from librational motions and hydrogen bond rearrangements in the solvation shell, eventually making the four oxygens indistinguishable on a time scale of tens of picoseconds.

The experimentally observed differences in the spectral diffusion properties of H₃PO₄ and the anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻ are related to fluctuations and inhomogeneities of the hydration shells, sensed via electrostatic coupling of the phosphate vibrational oscillator with the solvent environment. Thus, the electric field acting on the phosphates represents a collective coordinate of solvation, allowing a characterization of the hydration dynamics. We have calculated the fluctuations of the electric field at the bond midpoint, projected on the bond axis of P=O and P–O_H bonds (H₃PO₄, H₂PO₄⁻, and HPO₄²⁻) and at the bisector midpoint of PO₂ groups (H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻) [see Fig. 7(a)]. Electric field fluctuation amplitudes and electric field correlation functions are shown in Fig. 7. As demonstrated in Fig. 6, the gas phase tetrahedral symmetry of PO₄³⁻ is broken in a fluctuating liquid, leading to a localization of modes in PO₂ units, a behavior similar to other tetrahedral ions.⁴³ We have thus chosen to evaluate the electric field at the bisector midpoint of successive PO₂ groups of H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻.

The trajectory of the electric field imposed on the P=O bond of H₃PO₄ [Fig. 7(b)] reveals fast and large amplitude, sub-100 fs fluctuations due to librational motions of water molecules together with slower modulations on the ∼1–10 ps time scale. On average, the imposed electric field is substantial [84 MV/cm, Fig. 7(c)] with large amplitude modulations reaching up to 200 MV/cm, providing a rationale for the modification of the gas phase molecular structure of H₃PO₄ in the aqueous solution, i.e., the contraction of P=O and P–O_H bond distances from 1.54 and 1.69 Ź⁶ to 1.491 and 1.572 Å in the liquid. The temporal behavior of field fluctuations is characterized via the electric field fluctuation correlation function [Fig. 7(d)]. For the P=O bond of H₃PO₄, the fast (∼100 fs) decay is followed by a slower decay extending into the picosecond time scale. In comparison, the amplitude of the initial decay of the field correlation function for the P–O_H groups of H₃PO₄ is more pronounced, and slower dynamics on the 400–1000 fs time scale almost negligible. The weak residual inhomogeneity among the three P–OH sites stems from the strong hydrogen bonds of hydroxy groups, with the first hydration shell that lock the conformation of the hydroxy groups on the low picosecond time scale.

The electric field fluctuation correlation functions of H₂PO₄⁻ and HPO₄²⁻ [Figs. 7(e) and 7(f)] are very similar; for both ions, the correlation function recorded at PO₂ groups shows a longer lived plateau extending into the picosecond time scale due to the inhomogeneity of hydrogen bond geometries, while the correlation functions recorded at P–O_H sites display a fast (∼100 fs) decay.

The averaged electric field fluctuation correlation function recorded at the bisector midpoints of the PO₂ units of PO₄³⁻ [solid line in Fig. 7(g)] shows an initial fast (50–100 fs) decay, followed by slower dynamics on the 400 fs to 1.2 ps time scale. Notably, the different PO₂ oscillators exhibit a similar temporal behavior, but the amplitudes of the fast and slower time scale in the individual correlation function of different PO₂ sites vary [dashed lines in Fig. 7(g)]. The varying amplitudes of the slow temporal component in the correlation functions reflect inhomogeneities of the water population and water hydrogen bond geometries at the different sites of the PO₄³⁻ tetrahedron leading to an inhomogeneity persisting up to picoseconds in the PO₄³⁻ case.

Comparing the correlation decays of H₃PO₄ and the anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻, pronounced initial fast (50–100 fs) decay is followed by varying degrees of slower dynamics. The fast dynamics arise from librational motion predominantly of first solvation shell water molecules¹⁵ for a fixed hydrogen bond geometry. The slower (∼1 ps) time scale leads to a decay of residual correlation of correlation functions due to the randomization of hydrogen bond geometries that interchange on the 1–20 ps time scale and impose particularly pronounced inhomogeneities in the PO₄³⁻ case. The breaking and reformation of hydrogen bonds involves large angular jumps of water molecules stochastically occurring in the hydration shell.⁴⁴ The dual hydrogen bond functionality of H₃PO₄, H₂PO₄⁻, and HPO₄² is an additional source of fluctuations of the vibrational Hamiltonian due to the strong coupling of P–OH and PO₂ mode displacements. The detailed mechanism of vibrational couplings among these modes has been investigated in Ref. 26. In essence, the fast (quasi-homogeneous) fluctuations at a hydrogen bond donor site [cf. red lines shown in Figs. 7(e) and 7(f)] fully propagate into fluctuations of ν_AS(PO₂) vibrations, in addition to the field fluctuations experienced at PO₂ units. This behavior explains the absence of major inhomogeneities in the 2D-IR spectra of H₂PO₄⁻ and HPO₄²⁻. In contrast, due to the absence of H-bond donor functionality, the PO₂ units of PO₄³⁻ most clearly reveal inhomogeneities in the hydrogen bond arrangement around the ions with picosecond lifetime.

The experiments and simulations presented here are in the limit of independent phosphate hydration sites, i.e., interactions between neighboring sites play a negligible role. When increasing phosphate concentrations to the molar range, hydration shells of neighboring sites begin to interfere in space and influence each other. This scenario is particularly relevant for solvated PO₄³⁻ ions, where the first two hydration layers contain up to 40 water molecules in a structured arrangement.¹⁰ At high PO₄³⁻ concentrations, one expects a replacement of bulk-like water by an arrangement of solvent molecules structured by the ions.

It is interesting to compare the hydration dynamics of the tetrahedral ions PO₄³⁻ and SO₄²⁻ in the femtosecond time domain. While the current study reveals inherent inhomogeneities of hydrogen bond arrangement and hydration structure, inhomogeneities of the hydration shell of SO₄²⁻ on the ∼500 fs time scale are only observed in the presence of Mg²⁺ ions (but not for singly charged Na⁺ and NH₄⁺ ions).⁴³ The tightly bound solvation shell or PO₄³⁻ thus contrasts with a more labile and fluxional solvation shell of SO₄²⁻, where inhomogeneity is only induced via the formation of solvent shared ions pairs with Mg²⁺ ions.

In conclusion, the combination of femtosecond 2D-IR spectroscopy and ab initio molecular dynamics simulations reveals distinctly different hydration structures and dynamics of phosphoric acid H₃PO₄ and its anions H₂PO₄⁻, HPO₄²⁻, and PO₄³⁻. The first solvent shells around the solutes contain between 6 and 13 water molecules and a molecular arrangement with phosphate OH groups acting as hydrogen bond donors to water oxygens. The next two water layers show a different degree of structural ordering and rigidity. The latter is most pronounced around PO₄³⁻ ions. As a result of such structural differences, the fluctuating electric forces the water environment exerts on the vibrational phosphate oscillators display different time correlation functions. The predominant fast (50 fs) correlation decays around H₃PO₄ and H₂PO₄^- result in a quasi-homogeneous 2D-IR line shape of the asymmetric phosphate stretching vibration. The slower picosecond correlation decay components found for HPO₄²⁻ and PO₄³⁻ point to an incomplete randomization of electric forces. For PO₄³⁻ ions, such components manifest in an incomplete spectral diffusion on a 600 fs time scale and partly inhomogeneous 2D-IR line shapes.

This research received funding from the European Research Council (ERC) under European Union’s Horizon 2020 research and innovation program (Grant Agreements Nos. 833365 and 802817). The authors acknowledge the computational and data resources provided by the Leibniz Supercomputing Centre (www.lrz.de).

The authors have no conflicts to disclose.

T.E., A.K., and B.P.F. conceived the study. A.K. performed the experiments. A.K. and TE analyzed the experimental results. B.P.F. performed the theoretical calculations and simulations. The manuscript was written by all authors.

Achintya Kundu: Conceptualization (lead); Formal analysis (equal); Investigation (equal); Project administration (lead); Supervision (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Benjamin P. Fingerhut: Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Investigation (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Thomas Elsaesser: Conceptualization (lead); Formal analysis (equal); Investigation (equal); Project administration (lead); Supervision (lead); Validation (equal); Writing – original draft (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.