High energy vibrational excitations of nitromethane in liquid water

The dynamics of vibrational energy flow is central to many areas of physical chemistry, ranging from time dependent spectroscopy¹ to chemical reactivity.² A full understanding of this dynamics requires a detailed picture of the pathways through which the excess energy is dissipated, a complex task given that the relaxation dynamics generally involves the participation of a mixture of many degrees of freedom,^3,4 particularly in condensed matter systems. For the broad spectrum of problems for which a fully classical approach for the nuclear degrees of freedom is a sensible choice, it is, in principle, feasible to obtain a highly detailed picture of the energy pathways. Several measures of energy flow are available, among which the computation of energy fluxes, i.e., the time integrated power after an excitation, turns out to be particularly powerful. Initially devised for vibrational relaxation problems,^5–10 it was later formalized and generalized,¹¹ allowing to explore in detail a wide variety of excitations (translational,¹² rotational,^12,13 vibrational,^9–11,14 and electronic^15–19).

In this paper, we apply the energy flux methodology to the analysis of high energy vibrational excitations. Such conditions arise, for instance, in combustion reactions,^20,21 and of more interest to the present work, following the rapid, radiationless internal conversion of an electronic excitation.²² Photosynthesis and photostability of DNA under UV radiation are two prominent examples. These processes are well suited to classical mechanical approaches. The high density of states^3,21 results in vibrational flow between states whose energy separation is smaller than the typical thermal energy, this being a hallmark of a sound classical approximation. As an example of the methodology’s capabilities, the case of bend relaxation in neat water illustrates the level of detail attainable.^9,11 For example, it is found that the energy transfer to the surrounding water molecules is markedly asymmetric: the two h-bond accepting neighbors receive three times the energy that goes into the two h-bond donor neighbors, a quantitative result that could not be obtained with alternative methods.^9,11

This work explores the creation of spatial anisotropies in the energy fluxes that follow high energy vibrational excitations. The presence of such anisotropies may be crucial for applications at the nanoscale. There is a substantial activity in the design of active particles (usually at the micrometric scale), which, by the effect of one of several possible mechanisms (e.g., mechanical deformations, chemical reactions at their surface, and gradients in their vicinity), are capable of self-propulsion and, as such, are usually referred to as swimmers.²³ A case in point is that of Janus swimmers: it has been experimentally demonstrated²⁴ that for half-metallic micrometric spheres, laser heating of the metallic side produces a temperature anisotropy around the particle that results in self-propulsion (a phenomenon known as self-thermophoresis, one of several possibilities for propulsion based on local gradients). Recently, we have demonstrated that self-thermophoresis is also feasible at the molecular scale with a proof-of-concept model,²⁵ inspired, in this case, by the energy fluxes observed during the solvation dynamics that follows electronic excitation of a model solute in water.¹⁵ The swimmer is here a fullerene (diameter of $≈1$ nm) decorated with a dipolar moiety, which, upon electronic excitation by external radiation, inverts its dipole moment (a standard model in computational studies of solvation relaxation^26,27). For a liquid water solvent, this process indeed results in strong heating around the dipole (i.e., anisotropic heating) and self-propulsion, with a velocity dependent on the inversion rate.²⁵ 

Having demonstrated the feasibility of molecular self-propulsion, we now consider whether it can also occur for more realistic models. To this end, we have chosen to investigate the vibrational relaxation of nitromethane (CH₃NO₂). The hydrophilic end of this amphiphilic molecule has a stronger coupling with water than the hydrophobic end and, as a result, could be significantly more efficient in the ability to transfer energy to the surrounding solvent. Based on our earlier studies,²⁵ we recognize that additional criteria for self-propulsion must be met beyond anisotropic energy flow. First, the excitation energy deposited in the molecule needs to be substantial (∼50 kcal/mol). This energy is characteristic of electronic transitions and, therefore, available in sub-picosecond times in the form of vibrational energy for systems undergoing internal conversion. Second, to produce a significant effect, such excess energy must dissipate into the solvent in a few picoseconds. Again, this timescale is typical for vibrational relaxation in water, as shown in Ref. 28 for low energy vibrational excitations (⁠$≤10$ kcal/mol).

Our study of nitromethane builds on several previous works. Its small size has facilitated the development of an accurate intramolecular force field, extensively tested across different phases (gas,^21,29 liquid,^30,31 and crystalline solid^32–36). The vibrational relaxation of highly excited vibrational states has been addressed as well: for the crystalline solid³⁶ and in a supercritical bath^21,29 (argon). This previous work (with some additional extensions to be reported here) will facilitate a contrasted picture of the results for nitromethane in liquid water, a system not yet addressed.

The outline of the remainder of this paper is as follows: Our methodology is described in some detail in Sec. II, while Sec. III is devoted to a brief outline of the phenomenological models employed to describe the different molecules involved as well as describe their interactions. The results obtained together with a detailed discussion are reported in Sec. IV, which is divided Subsections IV A–IV C, respectively, devoted to isolated nitromethane molecule, nitromethane in argon, and, central to this work, nitromethane in liquid water. Concluding remarks are offered in Sec. V.

Our study is based on a detailed analysis of classical all-atom molecular dynamics (MD) simulations of a vibrationally excited nitromethane molecule (CH₃NO₂) in three different environments: vacuum (isolated nitromethane molecule), argon gas, and liquid water. In this section, we provide a detailed account of the various methodologies implemented (with in-house codes) for the analysis of the trajectories generated from all-atom MD simulations.

In a simulation of a molecule immersed in a liquid phase, the molecule undergoes translations, rotations, and internal distortions. To evaluate the normal mode displacements and velocities and calculate the corresponding vibrational energy of the molecule at a given time, its instantaneous configuration needs to be compared with an underlying undistorted configuration. Determining such configuration is a classic problem in the theory of molecular vibrations, solved by Eckart for non-linear molecules,³⁹ which requires the use of a non-inertial system of reference (Eckart frame).

Using this methodology, we can partition the energy following Eq. (3). We note that the projector method⁴⁶ has also been used in the literature for similar purposes, albeit it has been later established⁴⁵ that it is subject to a variety of drawbacks, particularly so for condensed phase simulations where there is no guarantee that the modes correspond to the Eckart frame.

A full understanding of molecular energy relaxation requires the identification of the pathways by which the excess energy flows. Identifying the pathways is crucial in the present context because we want to determine whether vibrational relaxation of nitromethane takes place anisotropically.

There are several possible approaches for following energy flow. As already emphasized, we will use the “energy flux” approach as it has several advantages over two indirect approaches in common use, which rely on energy partitioning. These approaches include the analysis of the correlation of energy fluctuations at equilibrium and/or direct monitoring of excess energy dynamics.^9,14 Unfortunately, neither of those approaches is fully satisfactory as they can miss key steps, provide ambiguous information on spatial extent, and be difficult to interpret. One finds, for instance, that in equilibrium simulations, energy time correlation functions are plagued with cross terms, making it difficult to disentangle individual contributions. As for nonequilibrium approaches (i.e., relaxation of excited states), monitoring of the excess energy as it spreads among degrees of freedom will miss intermediate fast relaxing states even though they may actually channel most of the excess energy. For both, the spatial extent of the process is ambiguous due to rapidly increasing noise with increasing distance from the initial excitation. In contrast, an energy flux approach has been shown to be free of these limitations, as demonstrated for a wide variety of excitations in neat liquid water.^{9,10,12,13,15–19} In each case, its use pinpointed the precise contribution of each motion involved (including ultrafast modes) and provided a precise map of the spatial extent of the process.

The complexity of the energy flux method depends sensitively on the desired number of pathways to be interrogated.¹¹ The analysis of the vibrational relaxation of water illustrates this point. To study relaxation of the lowest frequency mode (bending mode),^9,11 one only needs to consider energy fluxes into self-rotation and into neighboring solvent molecules. The solvent contribution can, in turn, be partitioned into first shell and the remainder. If instead one wants to study the whole manifold of fluxes after an excitation of the (higher frequency) stretching mode, the complexity increases substantially.¹⁰ A detailed picture for nitromethane, which would consider all intramolecular inter-mode channels plus transfer into self-rotation and into the solvent, while feasible, would be considerably more complex. However, since our main interest is to track the excitation energy transfer into the solvent, a simpler version will suffice.

For the case of nitromethane, an asymmetric energy transfer into the solvent would manifest itself as a difference between the work exerted by the molecule on the solvent surrounding the nitro and methyl groups. To capture this asymmetry, we partition the solvent into two regions. Using cylindrical coordinates r, ϕ, z in which the z-axis is parallel to nitromethane’s C–N bond and the origin is at the molecule’s center of mass, those solvent atoms with z_j > 0 are on the nitro group side and defined to be at the “front,” and those on the methyl group side are defined to be at the “back.” Correspondingly, work that nitromethane exerts on these two regions has contributions from W_front and W_back. For the liquid water case, a more detailed spatial analysis is carried out (see within).

The results of the flux calculations will depend sensitively on the forces. For this reason, it is essential to have confidence in the quality of the force fields that are employed. In this section, we provide a summary of the models used to describe the molecules involved and their mutual interactions for the different systems studied: isolated nitromethane, nitromethane in argon gas, and nitromethane in water.

For nitromethane’s intramolecular interactions, the Sorescu–Rice–Thompson (SRT)^30,32,33,49 force field was used. This potential energy surface (PES) is expressed in terms of a standard combination of stretching (Morse), bending (harmonic), and torsional (cosine) potential terms. Its parametrization is of hybrid type: a fit to (scaled) ab initio density functional frequencies is complemented with dissociation energies taken from previous phenomenologically based potentials and an increase in the torsional barrier (for methyl rotation) with respect to pure ab initio results. The frequencies that result from a normal mode analysis (see Sec. II B) can be found in Table I. One should be aware that the original SRT parameters are presented in Ref. 30, and a corrected set is reported in a later publication.³³ The parameters and functional forms used are detailed in Sec. I of the supplementary material.

The SRT force field provides a good match of the experimental spectrum, although one can note that there is a slight misordering for three modes. Specifically, M11 (basically a nitro group symmetric stretch) presents a higher frequency than the twin methyl bending modes (M9, M10), at variance with the observed ordering⁵⁰ and with ab initio calculations.^32,51 In an attempt to check the origin of such misordering, a new parametrization of the same functional form has been developed from purely density functional ab initio computations using a larger basis set than in the original SRT parametrization.⁴⁹ The resulting harmonic spectrum is displayed in the last column (“NEW” header) of Table I, and a detailed account of the procedure followed can be found in the supplementary material. The resulting parameters are similar to those in the SRT force field, which translates into rather similar frequencies. As expected, the slight discrepancies between the force fields are found for the methyl stretchings (related to not having used scaled values in this new parametrization), dissociation energies (not fixed), and in the torsional CH₃–NO₂ barrier (again not fixed). However, the misordering of frequencies is also present so that its origin should probably be attributed to a shortcoming of the functional form chosen for the PES.

In a separate set of simulations, the flexible SPC/fw water model⁵⁷ has been employed. Intramolecular interaction coefficients, charges, and LJ intermolecular interaction coefficients can be found in Ref. 57.

In Subsections IV A–IV C, we describe vibrational relaxation of nitromethane as an isolated molecule and solvated in supercritical Ar and liquid water. The simulation protocols are specific for each environment and are provided in each of the respective subsections. All numerical MD simulations have been performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)⁵⁸ MD package running on a local cluster.

Since the onset of anisotropic energy fluxes depends on a balance of timescales between intramolecular and intermolecular energy transfer, we begin our analysis with an investigation of energy transfer in the isolated nitromethane molecule. We start by describing previous works and then compare our results to them.

Recently, ab initio molecular dynamics (AIMD) calculations⁵¹ have been reported for low energy excitations of an isolated nitromethane molecule. Internal energy redistribution was found to occur on timescales ranging from ∼0.2 ps to times slightly greater than 1 ps, depending on the mode initially excited. These calculations provide an excellent benchmark for testing the SRT potential and for assessing the dependence of our relaxation results on force field parameterization. We will then extend our study to consider the extent to which internal dynamics might depend on the excess energy initially deposited.

The ab initio calculations in Ref. 51 consist of a series of time dependent density functional calculations, where the relaxation of each mode is analyzed in turn. Starting from an equilibrium configuration, an excitation of one (for high frequency modes) or two quanta (for low frequencies) of vibrational energy is injected into a parent mode as purely kinetic energy, while the rest of the modes are assigned thermal energies. The dynamics is averaged in each case over a set of 50 independent trajectories.

In our calculations, most of the parameters employed (excitation protocol, length of simulations) have been chosen to closely mimic those of the ab initio calculations. For initial conditions, a NVT trajectory is run with a Nosé–Hoover thermostat set at 400 K, with a coupling constant of 10 fs. Starting from configurations taken from this equilibration run, 200 independent NVE trajectories have been launched following the initial excitation of the selected mode as described in Sec. II C. The equilibration and production runs timesteps are 0.1 and 0.05 fs, respectively. Finally, we note that in the ab initio calculations, a projector method was employed for the partitioning of energy among the internal modes; here, the Eckart based partition, described in Sec. II B, has been used.

The quantities displayed in Fig. 1 can be compared to the results of Fig. 1 in Ref. 51. It should be noted that in Ref. 51, modes are ordered in terms of decreasing frequency, which is the inverse of the one used here (following Ref. 49). To minimize confusion, the panels in Fig. 1 are ordered to mirror the corresponding ones of Ref. 51.

A mode by mode comparison of the 14 panels shows that the same qualitative (and quantitative in most cases) behavior is observed. A noticeable difference is that those modes undergoing fast relaxation have shorter relaxation times with the SRT force field. While the shortest times in ab initio calculations are of the order of 0.2 ps, for SRT results, we find sub-0.1 ps decay for several modes. These differences at the very short timescales are unlikely to be relevant to the liquid phase simulations where the overall vibrational cooling (VC) takes place on a picosecond timescale.

We consider particularly relevant, in this comparison, that in both cases, it is possible to roughly divide the modes into two sets. A first set (basically those modes related with the methyl group) is characterized by a dominant sub-picosecond relaxation (⁠$∼0.2ps$⁠, a time after which they reach a plateau with about one sixth of the initial energy). In contrast, there is a second set (mostly related to the nitro moiety, lower two rows and first panel of second row in Fig. 1), which contains markedly slower modes. Therefore, irrespective of where the excess energy is initially injected, the nitro group, with its lower frequency modes, will necessarily be a major recipient of the excess energy on a picosecond timescale.

The results for the newly developed force field can be found alongside the SRT results in Fig. 1. We find very little differences between both force fields. Given this good overall accord, in the rest of this work, we will report the results for the SRT force field.

Finally, we note that the ab initio calculations (and our classical simulations) of the isolated were performed for relatively low energy excitations (⁠$∼200meV$ for the highest energy excitations, see upper leftmost panel in Fig. 1). However, in the work to be reported for argon gas and liquid water, we will be interested in excitations of $∼2000meV$⁠, i.e., one order of magnitude larger. It is therefore of interest to explore whether the intramolecular energy fluxes may be significantly altered for the isolated molecule. Simulations similar to the ones previously described have been performed, although now the same energy (200 kJ/mol–2000 meV) has been deposited on each mode of interest. The results are presented in Fig. S7 of the supplementary material. No substantial changes are found, modes centered on the nitro group (generally slow) are basically unaltered, while those located on the methyl group maintain a similar fast initial decay. This substantial similarity for high/low initial excitations energies will still be present when the molecule is placed in liquid water (see Sec. IV C).

The relaxation of high energy vibrational excitations of nitromethane in supercritical argon has received considerable attention.^21,29,47 While not the main focus of the present work, we believe that it is worth revisiting as the previous considerations for the isolated molecule allow for a different perspective, which, in addition, will be helpful for the analysis of relaxation in liquid water. The aforementioned studies relied on a microcanonical excitation protocol of nitromethane, equivalent to the reasonable assumption that relaxation is a two-step process, where fast internal equilibration is achieved prior to any substantial transfer into the solvent. An important outcome of the analysis of the isolated molecule is that both ab initio and classical calculations (for low and high excitation energies) coincide in that there are different timescales for the intramolecular dynamics. In particular, it seems possible that the slower internal energy relaxation of nitro modes (of a few ps) could induce an anisotropic energy transfer to the colliding argon atoms. To test this possibility, we have simulated relaxation in argon focusing on single mode excitations of nitromethane (see Sec. II C), an unexplored issue for the argon bath. Our goal is to explore the energy pathway in the process of vibrational cooling (VC, in which the molecule returns to the vibrational ground state⁵⁹).

Except for the initial excitation, the simulation protocol largely follows that of Refs. 21, 29, and 47. A single nitromethane molecule is placed in a cubic box with 1000 argon atoms with periodic boundary conditions. All atoms in the simulation box are given random uniformly distributed velocities re-scaled to be consistent with a temperature of 300 K. A 10 ps equilibration simulation is then performed with both a thermostat and barostat, where the Nosé–Hoover thermostat is set to 300 K with a 100 fs coupling constant and the barostat (also Nosé–Hoover) is set to 100 atm with a 10 fs coupling constant. An atom to atom cutoff of 15 AÅ was enforced for the Buckingham pair interactions. Finally, a timestep of 0.1 fs has been chosen. After equilibration, a long trajectory is generated during which the thermostat and barostat are maintained. Every 1 ps a configuration is saved as a seed for an independent simulation of vibrational relaxation. For these side simulations, the thermostat and barostat are disabled, and a (high energy) 200 kJ/mol instantaneous excitation of a specific normal mode of nitromethane is introduced (see Sec. II C). Each of these (NVE) simulations is 1000 ps long, with a timestep of 0.05 fs, with a total of 200 trajectories for each excitation.

We start by comparing the relaxation dynamics of two different excitations: M14 (a methyl stretch with fast relaxation for the isolated molecule, see the second panel of the top row in Fig. 1) and M2 (a nitro rocking mode with slow relaxation, see the second panel of the bottom row in Fig. 1). The results are displayed in Fig. 2, which contains the time dependent kinetic energy for modes M2 and M14 and the sum of all other modes up to a time of 5 ps for each of the two initial excitations. A non-negligible difference in relaxation can indeed be observed for short times, as the characteristic isolated molecule dynamics of each excited mode is still present: M14 decays rapidly to the point that after 5 ps, the excess energy is barely noticeable, while for M2, one fourth of the initial energy still remains. For both excitations, the rest of the modes increase their energy above their initial thermal energy, an increase not indicative of VC but of internal equilibration, as very little energy has actually been transferred into the gas on a timescale of a few picoseconds. This is clearly illustrated in Fig. 3, where total vibrational energies are displayed for each initial excitation. Here, the initial energy is given, on average, by the thermal equilibrium energy $(∼45kJ/mol)$ plus the excitation energy of 200 kJ/mol. On the 1000 ps time scale displayed, the curves are almost indistinguishable, including the one for the microcanonical excitation protocol. The slight dispersion between curves can be attributed to statistical indeterminacy, originated from the limited number of trajectories (200) used for each mode excitation.

As expected, the small differences at short times (few picoseconds) disappear altogether in the timescale of VC (hundreds of picoseconds), i.e., all modes relax with the same timescale. The cause lies in the inefficiency of energy transfer in the gas phase,²¹ estimated to be on average of $∼0.75kJ/mol$ per collision (to be compared with a total excess energy of 200 kJ/mol), and the low collision rate (a rough kinetic theory estimate predicts $∼1$ collision every 10 ps for a pressure of 100 atm). This is illustrated by a sample of one of the trajectories, displayed in Fig. 4, which in addition shows how the energy flux method can provide a vivid picture of the collisions. The top panel contains the total power exerted by nitromethane: collisions are easily located, being characterized by oscillations of the power (with a mixture of low frequencies, characteristic of center of mass atom/molecule vibrations, and higher frequencies resulting from nitromethane’s internal dynamics). The bottom panel displays the net work that the (excited) nitromethane molecule exerts on the surrounding gas, simply obtained from integration of the curve in the top panel. Over a time of a few picoseconds (characteristic of internal equilibration), energy transfer is small (a few kJ/mol for the most energetic collision at $∼3$ ps), to the point that it can be negative in some cases (for times up to $∼3$ ps in Fig. 4 bottom panel), even though the nitromethane molecule is highly excited.

We have thus seen that differences between modes are limited to very short times. Necessarily, any bias on the spatial distribution of the average energy transferred into the surrounding medium will only occur at short times as well. The comparison of the work exerted on the argon atoms located in the “front” and the “back” of nitromethane, W_front and W_back (see Sec. II), is displayed in Fig. 5, with no apparent imbalance between both quantities up to a time of 10 ps. Additionally, Fig. 5 further reinforces the picture of inefficient transfer of the energy transfer from the excited nitromethane to the argon gas (now averaged over all the trajectories), with only $∼3kJ/mol$ being transferred up to a time of 10 ps.

To summarize, the intramolecular dynamics for isolated nitromethane is characterized by a sub-picosecond transfer into low energy modes, which are mostly located on the nitro moiety (see Sec. IV A), an effect that should, in principle, be conducive to anisotropic relaxation. Nevertheless, this accumulation of energy only lasts for a few picoseconds before internal equilibration sets in (helped by the collisions with the gas). Once this equilibrium is reached, no prevalence of transfer from the methyl or nitro groups can be expected. As a result, at longer times, any initial excitation (microcanonic, mode-specific, kinetic, etc.) results in exactly the same relaxation behavior (see Fig. 3). To have anisotropic dissipation, the energy transfer from the nitro group into the surrounding gas should be extremely fast (≲1 ps) in order that a substantial amount of energy is transferred before internal equilibration is reached.

We finally turn to the central issue of vibrational relaxation of nitromethane in liquid water following a high energy vibrational excitation. It was argued in the Introduction that nitromethane should be a promising candidate if one seeks a molecule capable of sustaining non-isotropic thermal gradients in its immediate environment, but this entails substantial demands. The work in Ref. 25 showed that in order to create temperature gradients around the excited molecule, the lifetime for vibrational relaxation in water should be of the order of a few picoseconds. For nitromethane, we argued in Sec. IV B that to create an asymmetric dissipation, the energy transfer from the nitro moiety into the solvent should occur in a timescale of ≲1 ps. This is a short timescale, considering that for neat liquid nitromethane, the vibrational relaxation time is still one order of magnitude larger (⁠$≈13$ ps).³¹ 

To address this issue, in this section, we calculate the relaxation rate of nitromethane in liquid water from all-atom MD simulations and determine the role played in this process by rotating/translating modes of solvation water molecules. Finally, we show that such energy relaxation occurs anisotropically, with a significantly larger transfer of energy into the solvent surrounding the nitro moiety than into water molecules surrounding the methyl moiety.

Similarly to the argon gas simulations, a single nitromethane molecule is placed in a cubic box of side 21.7 AÅ solvated with 343 water molecules, which corresponds to a density of 1.0034 g/cm³, with periodic boundary conditions. Water molecules are kept rigid using SHAKE with a tolerance of 10⁻⁷. Long range interactions are computed using the particle–particle particle-mesh (PPPM) solver implemented in LAMMPS⁵⁸ with a relative force error of 10⁻⁵, while a cutoff of 10 AÅ is enforced for Lennard–Jones pair interactions. All atoms in the simulation box are given random uniformly distributed velocities, re-scaled to be consistent with a temperature of 300 K. Simulations start with a long (NVT) equilibration run, with a time step of 0.1 fs. We use a Nosé–Hoover thermostat set to 300 K with a 10 fs coupling constant. After 10 ps of equilibration, the thermostat is maintained in a long simulation run from which seed configurations, 1 ps apart, are stored. During this long run, the dihedral angles H–C–N–O were constrained as a measure to avoid mixing of modes related to the methyl group for the initial configuration; more details can be found in Sec. VI of the supplementary material. These seed configurations are used as initial configurations for runs for which a 200 kJ/mol excitation is injected into a specific normal mode of nitromethane (see Sec. II C), in exactly the same way as was done in argon gas. For each one, 5 ps NVE (no thermostating) simulation is carried out with a timestep of 0.05 fs. Results in this work correspond to an average over 500 trajectories for each mode. A representative set of excitations has been repeated using a flexible model for water (see Sec. III), without any changes of the procedure.

Results for the vibrational energy relaxation of nitromethane in liquid (SPC/E) water are displayed in Fig. 6. Each (solid) curve corresponds to the total molecular vibrational energy after excitation of a given mode (see Sec. II C), while the dashed curve corresponds to the case obtained with a microcanonical excitation protocol. Similar to the argon simulations, the initial energy on average is $∼45+200kJ/mol$⁠. Two aspects stand out. First, full relaxation is now attained for all cases in a picosecond timescale, a three orders of magnitude speed-up with respect to supercritical argon and one magnitude faster than for neat liquid nitromethane.³¹ The second column in Table II provides a quantitative estimation of the decay time for each individual mode, with the average being 1.5 ps (first row, corresponding to the microcanonical average). A second feature that attracts attention is the existence of a substantial dispersion in relaxation times, which encompasses one order of magnitude (with a minimum of 0.23 ps for M5 and a maximum of 2.56 ps for M15), in stark contrast again with the corresponding results for argon gas (see Fig. 3). Two main groups can be distinguished: fast modes (M2–M5) with sub-ps relaxation times and slow modes (M1 and M6–M15). Fast modes share two main features: they have nitro-associated deformations and low frequencies. These are the modes with slow internal relaxation to other modes, as found in the study of the isolated molecule (see Sec. IV A). Correspondingly, slow relaxing modes in liquid water are associated with methyl-centered distortions, with relatively higher frequencies and fast internal relaxation to other modes for the isolated molecule.

In general, an important increase of the energy relaxation rate in the liquid phase can be attributed to the caging of the solute, which (in a gas-phase inspired picture) results in a much higher rate of collisions per unit time. Another important factor for the energy relaxation rate is the efficiency of the energy transfer per “collision” (where frequency matching between solute and solvent modes is an important factor). Neat liquid water stands out in this connection as it is characterized by extremely high relaxation rates: for low energy (one quantum) excitations, both the vibrational stretch⁶⁰ and bend⁶¹ decay in $∼0.2$ ps. The energy relaxation of excited polar solutes in liquid water can also be ultrafast (picosecond timescale), as shown in the pioneering theoretical studies for a diatomic model of (polar) methyl chloride.^5,6 The explanation of fast relaxation in water was found to be two-pronged,⁶ with a first factor being the strong Coulomb forces present, which speed-up relaxation by orders of magnitude compared to a non-polar liquid (compromising the utility of the gas-phase concept of short-lived collision). The second factor that arguably makes water special as solvent is the existence of an extremely broad librational band that can channel the excess energy, which results from water’s distinctly small moment of inertia.^5,6,62

In the case of nitromethane, the unique properties of liquid water permit an ultrafast vibrational cooling (of $≈1.5$ ps) one order of magnitude faster than for neat crystalline³⁶ and liquid³¹ nitromethane (⁠$≈13$ ps) and three orders of magnitude faster than for nitromethane in supercritical argon²¹(⁠$≈1750$ ps).

Figure 7 shows the results for the overall flux of energy into translations/rotations of the surrounding water solvent (again for a microcanonical excitation protocol). As expected, direct transfer into water rotations is dominant, accounting for $∼72$% of the initial excess energy. To gain some perspective, in the relaxation of the bending mode of a water molecule in neat water, $∼95$% of the transferred energy goes into rotations (⁠$∼2/5$ to rotation of the molecule’s neighbors and $∼3/5$ into self-rotation by centrifugal coupling).^9,11 While rotation participation is thus smaller for nitromethane, it is still comparable to other ultrafast processes for solutes in water. For the case of solvation relaxation after ionization of an initially neutral atomic solute in water, an entirely different process, it is found that a very similar $∼78$% of the initial excess energy is transferred into neighbor’s rotations.¹⁵ The present case is thus aligned with our previous findings, confirming the crucial role of librations in efficiently funneling the initial excess energy.

The energy flux method allows one, with no increase in computational complexity, to further extend this analysis to each individual mode excitation. Figure 8 compares the results for the fastest and slowest modes (M5 with a 0.23 ps relaxation time and M15 with 2.55 ps). While this substantial difference in rate is perfectly apparent at short times, at $∼5$ ps, the share of energy transferred to translations and rotations has almost converged to the (average) values obtained for microcanonical excitation (although saturation is not yet achieved for the slower mode). The results for all modes can be found in Fig. S8 of the supplementary material. The basic feature to extract from these plots is the prevalence of energy transfer into water librations, which applies to each individual mode. This does not preclude that a non-negligible dispersion exists, particularly clear when comparing M2 (nitro group rocking) and M3 (nitro group wagging): while the latter displays the largest weight of transfer into water rotations, M2 shows a minimal difference between transfer into rotations and translations (with transfer into rotations still prevailing).

The amphiphilic nature of nitromethane should result in a stronger coupling of the nitro moiety with the solvent as compared to the methyl end, contributing to an efficient and faster energy transfer on the nitro side of the molecule. The differential coupling of the nitro and methyl moieties with liquid water is determined by nitromethane’s polarity and chemical composition: it is characterized by a dipole moment (3.46 D, larger than that of the water solvent), directed from the nitro group (with two oxygens interacting with the solvent) to the methyl group (with three hydrogens in contact with the solvent). Consequently, it is the difference in the oxygen–hydrogen interactions at nitromethane’s ends, which determine the asymmetric coupling of the molecule with water. The strength of the different oxygen–hydrogen interactions can be illustrated by comparing the two distinct oxygen–hydrogen radial distribution functions (RDFs), displayed in the two panels of Fig. 9: (nitro)oxygen–(water)hydrogen (for the nitro end of the molecule, O–H of the top panel) vs (methyl)hydrogen–(water)oxygen (for the water-methyl pair, H–O of the bottom panel). Focusing first on the top dashed and bottom solid curves (O–H RDFs), it is obvious that the first peak, while located in both cases at $∼2$ Å, is considerably stronger on the nitro side. This feature contrasts with a barely noticeable peak on the methyl side and corroborates the picture of a substantially stronger hydration on the nitro end. The additional RDFs included in Fig. 9 facilitate a more nuanced picture. The top panel, which corresponds to the nitro moiety, also displays the O–O RDF (solid), with a strong first peak located at $∼3$ Å (i.e., 1 Å further away than the first peak for the O–H RDF, dashed). Considering that the internal O–H bond of a water molecule is also $∼1$ Å, both RDFs (O–H and O–O) taken together signal the typical co-linear O–H⋯O structure characteristic of strong hydrogen bonds. A different arrangement is found on the methyl side (bottom panel), where the O–H and H–H RDFs are barely distinguishable, a result consistent with the known tendency of water O–H bonds to be parallel to hydrophobic surfaces, with some pointing away.⁶³ 

The stronger coupling of the nitro group, together with vibrational frequencies that overlap with the librational band of water, explain the extremely fast relaxation of nitro-centered modes. We note that relaxation of modes M2 (nitro rocking), M3 (nitro wagging), and M5 (nitro symmetric stretching) are characterized by relaxation rates almost identical to that of the vibrational modes of neat water itself (⁠$∼0.2$ ps as noted before), with M4 (nitro symmetric bending) slightly slower but still in the sub-ps.

We now address the central issue of interest, namely, the onset of anisotropy during vibrational relaxation. This is an issue intimately linked to the relaxation time dispersion shown in Fig. 6 and summarized in Table II. It is worthwhile to emphasize that there is no a priori reason why the existence of such dispersion should imply per se the onset of anisotropy. However, it turns out that for the particular case of nitromethane (and it is to be expected for molecules that have similar traits), several factors concur to establish this link. First of all and as discussed in Sec. IV A, the internal dynamics is one that favors fast internal relaxation of the higher frequency modes, partly into lower frequency modes. Second, the modes in nitromethane have a marked local character (see Table I), with lower frequency modes preferentially involving nitro group motions and higher frequency modes related to motions of the methyl group. Finally, the amphiphilic nature of nitromethane induces a more efficient and faster energy transfer onto the solvent surrounding the nitro side of the molecule. It is the combined effect of these factors that should produce a distinct anisotropy for the energy flux into the surrounding water solvent.

To characterize the asymmetry of the energy dissipation into liquid water, we first use the simple measure employed for the argon bath, i.e., the comparison of the work exerted on the solvent at the nitro end [W_front(t)] and the methyl end [W_back(t)]. Figure 10 shows the time dependency of both quantities, still for a microcanonical excitation protocol. It unambiguously shows the existence of a clear asymmetry, with W_front doubling the value of W_back, a feature that is largely independent of time during the vibrational relaxation.

The monotonous increase of both curves allows one to further simplify the analysis by focusing on the difference between the net (final) work performed on the solvent at the front of nitromethane minus the one performed at the back region, W_front − W_back. A positive sign indicates that more energy is transferred to the nitro side, and a negative sign indicates a dominance of transfer to the methyl side. Finally, we define the degree of asymmetry as the index that results from normalizing this difference (i.e., dividing by W_front + W_back) and expressing it as a percentage. The results for each mode are listed on the second column of Table II. The first row corresponds to microcanonical excitation, with an asymmetry index of 30%, which indicates that W_front doubles W_back, as previously noted. It is particularly remarkable that all excitations show a positive asymmetry, even for excitations involving methyl centered modes. In addition, of the 15 vibrational modes, 12 of them show a degree of asymmetry close to $∼30$%. The only exceptions correspond to three low frequency modes: M2 (rocking of the nitro group) with a 47% asymmetry, M5 (stretching of the nitro group) with the maximum value of 67%, and finally M1 (relative rotation between methyl and nitro groups) with just a 4%. As we see, the first two cases, with maximum asymmetry, correspond to the nitro group as expected. The low asymmetry of M1 is consistent with its nature, a relative rotation between nitro and methyl moieties, i.e., a nonlocal concerted motion.

The asymmetry of energy transfer into the solvent can be obtained with higher spatial resolution. As described in Sec. II D, the energy flux method allows one to partition the work, exerted by nitromethane, in many possible ways. Instead of the coarse-grained partition (front/back) used so far, we will now focus on the average work exerted on each point in space around the nitromethane molecule. This can be obtained from adding (and averaging) the work on the water molecules at a small region (“pixel”) around each point of a conveniently defined mesh during the nonequilibrium simulations. Each point in space can be identified in terms of cylindrical coordinates: the distance ρ to the molecular axis (defined by the C–N bond), the height z along that axis (taking the center of mass as the origin), and the azimuth (angle around the axis). In order to simplify the visualization, we exploit the (almost) axial symmetry of nitromethane, condensing the results into a two-dimensional representation. This can be accomplished by summing over azimuths so that each pixel will contain the total work at a given radius ρ and height z.

The resulting two-dimensional map of the work exerted by nitromethane onto the solvent is displayed in Fig. 11 for a microcanonical excitation protocol. It further clarifies the nature of the asymmetry previously discussed in terms of a simple quantitative index. First and foremost, asymmetry is now clearly visualized, with most of the colored pixels corresponding to positive work located at z > 0, around the nitro moiety of nitromethane (“front” in our previous discussion). In addition, the work exerted in those regions is larger than the work done around the methyl side of the molecule (“back”), shown in Fig. 11 by pixels with higher color intensity. An additional feature to note, missed by the asymmetry index, is that energy transfer is highly local. Only a small subset of pixels (those in the immediate vicinity of nitromethane) show a net work, while the rest are characterized by a negligible net work (within statistical indeterminacy). It should be noted in addition that only a portion of the total simulation box is actually displayed in Fig. 11 (up to a radius of 7 Å, to be compared to a simulation box of side $∼22$ Å), further highlighting the small region with a non-null work. Again, this markedly local transfer of energy is in line with the behavior found in a variety of non-equilibrium relaxation processes in liquid water, for which it has systematically been found a transfer to librations of the immediate water solvent neighbors.

Given the high energies involved, it seems fair to ask whether the relaxation dynamics depends on the initial excitation energy or, in other words, whether linear response holds in this regime. To this end, we have computed the relaxation behavior for a wide range of initial energies. The normalized results, shown in Fig. 12, correspond to initial excitation of mode M5. This mode was selected because it has the fastest initial decay (see Fig. 6) and, therefore, is expected to be more sensitive than other modes to increasing amounts of initial excess energy. However, the results Fig. 12 do not show any noticeable signs of an excitation energy dependence of the relaxation dynamics, i.e., substantially increasing the initial excess energy does not hinder the ability of the water solvent to channel it. To our knowledge, this is the first reported case supporting, for the relevant case of liquid water, a behavior similar to that reported in Ref. 64 for vibrational relaxation in liquid argon and postulated to have wide applicability. It constitutes an additional illustration of the remarkable range of validity of linear response in liquid water, an issue that has been recently addressed¹² in detail for rotational and translational excitations in neat liquids, including water. The absence of any dependence is actually very similar to the findings for translational excitations (one can see, for instance, Fig. S3.2 in Ref. 12, where the largest excitations are comparable to the present ones). This extends to rotational excitations as well, although in this case a certain loss of coherence is observed for increasing excitation energy¹² but still far from the observations of total loss of coherence observed for some systems in Ref. 65.

The results reported are obtained from all-atom MD simulations, which employ the SPC/E rigid model for water. Considering that nitromethane has several modes with frequencies close and above the frequency of the water bend mode, one might ask whether the use of a flexible model might have an impact on the results. We recall though that a more relevant factor is the high density of states and the concomitant small energy spacing, which should dispel any such effect. To conduct our test, we note that mode M12 is nearly resonant with the water bend mode, especially for the SRT model employed for nitromethane. Separate simulations using the flexible SPC/fw model of water (see Sec. III) were performed in which the M12 mode of nitromethane was excited. No noticeable difference with the results reported employing the rigid SPC/E water model is found for the energy relaxation dynamics (see Fig. S9 of the supplementary material). We note that a similarly negligible effect was found as well for the case of solvation relaxation,¹⁵ where comparatively high energies were dissipated without any noticeable differences between rigid and flexible models for the water solvent.

The core finding of this work is the existence of a marked energy flux anisotropy during relaxation of high energy vibrational excitations of a nitromethane solute in liquid water. This anisotropy has been shown to take place irrespective of the excited vibrational mode, and the main factors have been identified, particularly the crucial role of the nitro group. Its strong coupling with water, combined with mode frequencies that overlap with water’s spectrum, results in an extremely efficient transfer of energy into the water solvent librations (within a timescale of $∼0.2$ ps). While these features are similar to those found for low energy excitations in neat liquid water,⁹ the present case shows that ultrafast relaxation is still present for highly energetic excitations and, more importantly, opens the way to external control as the solute can be addressed (and actuated) with tailored electromagnetic radiation. It seems reasonable to assume that a similar behavior might be present, and exploited, not only for nitromethane but for molecules with similar traits.

Among possible applications, it is of particular interest to test whether periodic (high energy) excitations might facilitate the onset of self-thermophoresis.^23,25 It should be considered, though, that several other important factors are at play, most notably Brownian rotation of the solute. It is well known²³ that as the solute radius decreases, motion becomes more random and less directed so that the impact of the self-thermophoretic effect is muted. Nitromethane certainly constitutes a worst case scenario as far as size is concerned. However, there are several ways to still exploit self-propulsion. We showed in Ref. 25 that it is possible to attain full control of the particle’s dynamics in reduced dimensionalities. In this case, the direction of (one-dimensional) motion inside a nanotube could be selected at will.²⁵ It has recently been shown that nitromethane, when dissolved in liquid water, is preferentially absorbed to the interior and surface of carbon nanotubes,⁶⁶ suggesting that it might be possible to exploit nitromethane’s excited state dynamics in one-dimensional and two-dimensional environments, respectively.

The supplementary material includes force field parameters, fitting curves, details on the trend lines displayed on Figs. 1 and 2, results for high energy excitations on isolated nitromethane, additional results of relaxation and work of nitromethane in liquid water, a figure comparing the relaxation of a M12 excitation in rigid SPC/E and flexible SPC/fw water, and details on the constraint of the H–C–N–O dihedrals mentioned in Sec. IV C.

A.J.R. acknowledges the financial support from the Departament de Recerca i Universitats de la Generalitat de Catalunya. C.C. acknowledges support from Spanish Grant No. PID2021-124297NB-C31, funded by MCIN/AEI/10.13039/501100011033 and “ERDF A way of making Europe.” E.L.S. gratefully acknowledges support from the NSF via Grant No. CHE-1900095. A.J.R. and R.R. acknowledge support from Grant No. PID2021-124297NB-C32, funded by MCIN/AEI/10.13039/501100011033 and FEDER, and from the Generalitat de Catalunya (Grant No. 2021 SGR 01411).

The authors have no conflicts to disclose.

Arnau Jurado Romero: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Carles Calero: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Edwin L. Sibert III: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Rossend Rey: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

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