Influence of the excitation energy of a probe solute molecule on its solvation dynamics and emission spectrum in 1-ethyl-3-methylimidazolium hexafluorophosphate (EMI+PF6) is studied via molecular dynamics simulations using a coarse-grained model description. By exciting the probe at different energies, each with an extremely narrow distribution, ensuing solvent relaxation and its dynamic variance are monitored using the isoconfigurational ensemble method. Resulting Stokes shift function, S(t), indicates that long-time solvent relaxation becomes slower with the decreasing excitation energy and approaches the equilibrium correlation function, C(t), of solvent fluctuations. This suggests that the system excited at the red-edge of the spectrum observes linear response better than that at the blue-edge. A detailed analysis of nonequilibrium trajectories shows that the effect of initial configurations on variance of relaxation dynamics is mainly confined to short times; it reaches a maximum around 0.1 ≲ t ≲ 1 ps and diminishes as time further increases. The influence of the initial velocity distribution, on the other hand, tends to grow with time and dominates the long-time variations of dynamics. The emission spectrum shows the red-edge effect in accord with previous studies.

Room-temperature ionic liquids (RTILs) consisting of bulky asymmetric organic cations and anions have attracted broad interest because of their unique properties, such as high chemical, electrical, and thermal stabilities, low vapor pressure, and large electrochemical window.1–3 It has also been suggested that RTILs are a new type of solvents or electrolytes that can replace conventional aqueous and organic solvents in a wide range of applications in chemistry and electrochemistry. Due to their strong Coulomb interaction and disordered molecular structure, RTILs can exist in a liquid state at room temperature, but have very high viscosity, which induces distinctive properties in dynamics.4–6 

Among others, solvation dynamics in RTILs have received a great deal of experimental7–17 and theoretical18–33 attention in the recent years. Their relaxation over a wide range of time scales spanning hundreds of femtoseconds to tens and hundreds of nanoseconds reveals RTILs’ unique dynamic character. In the case of imidazolium-based ILs, for example, solvation dynamics show biphasic relaxation—ultrafast subpicosecond dynamics followed by a very slow non-exponential decay.

Interestingly, solvation free energy profiles of RTILs were found to be approximately parabolic, indicating that the Marcus free energy relationship is generally valid for electron transfer reactions in RTIL environments.34–36 From this perspective of nearly harmonic potentials, nonequilibrium solvent relaxation is closely related to equilibrium solvent fluctuation dynamics.19 Recently, a considerable experimental effort has been focused on excitation-wavelength dependence of solvation dynamics in RTILs.37–40 Related dependence of steady-state emission spectra on the excitation energy, so-called red-edge effect (REE), has also been examined for RTILs both experimentally41–43 and theoretically.44–46 Time-resolved spectroscopy measurements indicate that solvation dynamics vary with the excitation wavelength.37–40 In some experiments, solvation dynamics were found to become faster with increasing excitation energy37,40 but in some others, the opposite trend was obtained.38,40 Thus while it is generally accepted that heterogeneity of RTILs is responsible for excitation-wavelength dependent solvation dynamics, there is a clear need for further analysis to better understand this interesting phenomenon.

In most of the previous simulation studies, initial configurations needed for solvent relaxation subsequent to the electronic excitation of the probe solute were chosen randomly from configurations equilibrated to the probe in the ground electronic state. The resulting initial configurations are characterized by a broad distribution in energies and thus do not allow a detailed study of the dependence of nonequilibrium solvation dynamics on the excitation energy. In this article, we circumvent this difficulty by considering multiple sets of initial configurations, each set with an extremely narrow distribution of excitation energy with a different mean. The comparison of solvent relaxation obtained with different sets of initial configurations enables a direct analysis of its dependence on the excitation energy. To gain insight into the effects of structural heterogeneities on solvation dynamics, we also investigate the time evolution of variance of emission energies using the isoconfigurational ensemble method. As a prototypical RTIL, we consider 1-ethyl-3-methylimidazolium hexafluorophosphate (EMI+PF6) using a coarse-grained model description and compared with aprotic acetonitrile.

This paper is organized as follows: Model descriptions and simulation methods are briefly explained in Section II. In Section III, we analyze the Stokes-shift function and the steady-state emission spectrum, focusing on their dependence on the excitation wavelength. Comparison of nonequilibrium and equilibrium solvation dynamics is also made there. Concluding remarks are offered in Section IV.

The simulation box comprises a single diatomic solute immersed either in 256 pairs of EMI+ and PF6 or in 512 molecules of CH3CN. Two constituent atoms of the rigid solute molecule are separated by 0.55 nm and interact with the solvent molecules through the Lennard-Jones and Coulomb potential. The atomic mass of each atom of the solute is 154.64 amu and LJ parameters are σ = 0.4 nm and ϵ = 0.83 kJ mol−1. We considered two types of solute charge distributions: a neutral pair (NP) with no charges and an ionic pair (IP) with charges q = ± 0.6364e, which yield a dipole moment of 16.8 D (e = elementary charge). For EMI+PF6, we used the coarse-grained description of Ref. 4, in which EMI+ are modeled as four sites connected through rigid bonds and PF6 is described as a united atom. For acetonitrile, the all-atom description of Ref. 47 was employed.

All simulations were performed in the canonical (NVT) ensemble with a Nosé-Hoover thermostat using the GROMACS packages.48 The box size and the temperature employed were, respectively, L = 4.364 nm and T = 350 K for EMI+PF6 and 3.555 nm and 300 K for acetonitrile. Periodic, cubic boundary condition was applied and long-range electrostatic interactions were computed with the particle-mesh Ewald summation method. The trajectories were integrated via the Verlet algorithm using a time step of 2 fs with data saved every 20 fs. For equilibrium simulations, 5 different trajectories were considered for both EMI+PF6 and CH3CN. For each trajectory, the production run was 25 ns long for EMI+PF6 and 30 ns long for CH3CN.

The Franck-Condon (FC) energy between solute electronic states a and b with different charge distributions,

ΔEabEbEa,
(1)

is widely used as a collective molecular variable49 in simulation studies to make direct contact with time-resolved spectroscopy measurements. Nonequilibrium simulations are usually performed by instantaneously switching the solute electronic state from b to a (i.e., FC excitation) in the presence of a solvent configuration equilibrated to b and by averaging the ensuing time evolution of ΔEab over an ensemble of initial configurations. Thus the preparation of a reliable set of initial configurations is crucial for a nonequilibrium study. Since our main goal is to scrutinize the excitation-wavelength-dependence of solvation dynamics, we considered different excitation energies and their corresponding sets of initial configurations. Specifically, we selected 8 different excitation energies by shifting ΔEab from its equilibrium value to lower energies, i.e., ΔEab = 〈ΔEab〉 − ΔEshift, where ΔEshift measures the degree of red shift of the excitation wavelength. Here 〈ΔEab〉 (≡〈ΔEaba) is the equilibrium absorption energy and 〈⋯〉a denotes an equilibrium ensemble average in the presence of the solute in state a (=NP or IP). In the present study, we considered ΔEshift values of 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol.

For each of the ΔEshift values, we selected 400 configurations from equilibrium configurations in such a way that they are at least 50 ps apart from one another and satisfy |ΔEab − (〈ΔEab〉 − ΔEshift)| ≤ σE with σE = kBT/300. It is well known that if the ΔEshift value exceeds the energy fluctuations at equilibrium, sampling using the standard Boltzmann distribution becomes unreliable due to limited molecular dynamics (MD) statistics. We thus introduced solute charges intermediate between NP and IP à la free energy perturbation method50 to generate reliable ΔEab distributions centered at differing absorption energies. Three intermediate charge states were considered: QP with q = ± 0.1591e, HP with ±0.3182e, and TP with ±0.4773e. For illustration, the probability distributions of ΔEIP→NP obtained with these intermediate charges as well as with NP and IP are shown in Fig. 1. It is worthwhile to note that the distributions obtained with the NP and IP states hardly overlap, clearly exposing the difficulty with the standard Boltzmann sampling when ΔEshift is large. The distributions obtained with QP, HP, and TP ensure that a broad region of ΔEIP→NP from its equilibrium value, 〈ΔEIP→NP〉, to lower energies and thus a wide range of red shifts in the excitation energy are covered reasonably well statistically.

FIG. 1.

Distribution function of ΔEIP→NP of equilibrium simulation for five solute state, NP, QP, HP, TP, and IP, in EMI+PF6 (solid line) and acetonitrile (dashed line).

FIG. 1.

Distribution function of ΔEIP→NP of equilibrium simulation for five solute state, NP, QP, HP, TP, and IP, in EMI+PF6 (solid line) and acetonitrile (dashed line).

Close modal

Once a set of 400 configurations was generated for each of the 8 different values of ΔEshift using equilibrium trajectories obtained with NP, QP, HP, TP, and IP solutes, we constructed 50 different sets of velocities according to the Maxwell distribution and generated a velocity or isoconfigurational ensemble5,51 for each of the 400 configurations. Different members of an isoconfigurational ensemble have different sets of initial momenta assigned to the molecules/ions but they all share exactly the same configuration of initial positions for the molecules/ions. This enables to separate variances arising from initial structure and initial dynamics of the solute-solvent system and analyze their respective influence on subsequent relaxation. All nonequilibrium simulations were conducted for 1 ns in EMI+PF6 and for 200 ps in acetonitrile. A total of 20 000 different microstates in phase space were employed to compute averages of nonequilibrium relaxation for given initial excitation wavelength, i.e., a given ΔEIP→NP or ΔENP→IP value.

Simulation results for the time-correlation function, Ca/b(t), of solvent fluctuations at equilibrium and dynamic Stokes shift function, Sa/b(t), of the nonequilibrium FC energy relaxation defined as

Ca/b(t)δΔEab(t)δΔEab(0)(δΔEab(0))2
(2)

and

Sa/b(t)ΔEab(t)¯ΔEab()¯ΔEab(0)¯ΔEab()¯,
(3)

are presented in Fig. 2. In Eqs. (2) and (3), δΔEab is the deviation of ΔEab from the equilibrium ensemble average 〈ΔEab〉, and ΔEab(t)¯ denotes the nonequilibrium ensemble average of the ab transition at time t after an instantaneous change of solute charge state from b to a at t = 0. As mentioned above, an ensemble of 20 000 trajectories is used to calculate Sa/b(t) in our study.

FIG. 2.

Left: Comparison of the normalized time-correlation function, Ca/b(t), and the Stokes-shift function, Sa/b(t), for NP/IP and IP/NP cases in EMI+PF6 (top) and acetonitrile (bottom). Middle: Sa/b(t) for NP/IP case in EMI+PF6 (top) and acetonitrile (bottom) with various ΔEshift, 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol. Right: Sa/b(t) for IP/NP case with the same red-shift as the middle panels.

FIG. 2.

Left: Comparison of the normalized time-correlation function, Ca/b(t), and the Stokes-shift function, Sa/b(t), for NP/IP and IP/NP cases in EMI+PF6 (top) and acetonitrile (bottom). Middle: Sa/b(t) for NP/IP case in EMI+PF6 (top) and acetonitrile (bottom) with various ΔEshift, 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol. Right: Sa/b(t) for IP/NP case with the same red-shift as the middle panels.

Close modal

We first consider the results in the left panels, where Ca/b(t) and Sa/b(t) (with ΔEshift = 0) are compared. Overall, our results compare reasonably well with those of an earlier study19 though the all-atom and united-atom descriptions employed there for EMI+PF6 are quite different from our coarse-grained model description. Both Ca/b(t) and Sa/b(t) are characterized by biphasic relaxation, i.e., ultrafast subpicosecond initial relaxation and a very slow nonexponential decay. The former and the latter are attributed, respectively, to hindered translational and diffusive motions of RTIL ions.18–20,24–26,30,31 We notice that Ca/b(t) and Sa/b(t) in the presence of the solute in the same ground state exhibit similar long-time behaviors in EMI+PF6. This confirms the conjecture made in Ref. 19 that equilibrium and nonequilibrium solvation dynamics, if they occur under similar RTIL environments, will share similar temporal behaviors. By contrast, the results for Ca/b(t) and Sa/b(t) in acetonitrile are almost superimposable regardless of their solute electronic states, revealing that linear response holds very well in this solvent.

Properties of equilibrium solvation dynamics we obtained are summarized in Table I. Our result for 〈ΔEab〉 for EMI+PF6 is 58.7 kcal/mol, which is smaller than that in Ref. 19 obtained with the IP solute charges of ±1.0e. This is as expected because 〈ΔEab〉, which gauges the solute-solvent electrostatic interaction, is nearly proportional to the magnitude of the diatomic solute charges. This trend was also observed in Ref. 22. The solvent frequency, ωS((δΔĖab)2(δΔEab)2)1/2, which measures the inertial relaxation of Ca/b(t), increases as the solute ground state is changed from NP to IP. This behavior arises primarily from 〈(δΔEab)2〉 that becomes reduced in the presence of a charged solute, compared to a neutral solute. The enhancement in solvent structure near the solute induced by strong solute-solvent Coulomb interactions (“electrostriction”) and accompanying increase in structural rigidity make solvent fluctuations difficult, thereby lowering 〈(δΔEab)2〉.52 Because (δΔĖab)2 is similar for IP and NP, a decrease in 〈(δΔEab)2〉 yields an increase in ωS. Therefore, the short time relaxation of Ca/b(t) becomes accelerated with increasing solute charge distribution.

TABLE I.

Equilibrium property.a

Solventa/b〈ΔEab〈(δΔEab)2(δΔĖab)2ωS
EMI+PF6 NP/IP −0.141 55.1 2431 6.6 
 IP/NP 58.7 38.2 2443 8.0 
CH3CN NP/IP −0.165 34.0 12 119 18.9 
 IP/NP 55.6 27.8 11 971 20.7 
Solventa/b〈ΔEab〈(δΔEab)2(δΔĖab)2ωS
EMI+PF6 NP/IP −0.141 55.1 2431 6.6 
 IP/NP 58.7 38.2 2443 8.0 
CH3CN NP/IP −0.165 34.0 12 119 18.9 
 IP/NP 55.6 27.8 11 971 20.7 
a

NP and IP states are assumed to be degenerate in vacuo. Energy unit: kcal/mol, time unit: ps.

The middle and right panels of Fig. 2 show Sa/b(t) with different ΔEshift. Though there are discernible differences, the overall relaxation behavior of Sa/b(t) does not vary too much with ΔEshift. This appears to be in agreement with experimental results. Similar to Ca/b(t), Sa/b(t) show biphasic relaxation irrespective of initial excitation energies.

Though not huge, we quantitate the influence of ΔEshift on solvent relaxation by fitting Ca/b(t) and Sa/b(t) as a combination of a Gaussian and a stretched exponential,30 

f(t)=fGexp(ωG2t2/2)+(1fG)exp{(t/τ)β}.
(4)

The fitting parameters of Ca/b(t) and Sa/b(t) are compiled in Table II. There are several noteworthy aspects. First, as ΔEshift increases, the short time relaxation of Sa/b(t) becomes faster (i.e., its ωG increases) in the IP/NP case but becomes slower (and its ωG decreases) in the NP/IP case. This is closely related to the observation made above that the initial relaxation of Ca/b(t) becomes faster as the solute charge increases. Specifically, in the case of SIP/NP(t), an increase in ΔEshift makes the initial configurations more “IP-like” environments and thus accelerates the solvent relaxation, while it is the opposite in the case of SNP/IP(t). Since increase in ΔEshift weakens the perturbation to the system caused by the initial excitation, linear response is better followed with increasing ΔEshift. One consequence of this is that ωG of Sa/b(t) becomes closer to that of Ca/b(t) with increasing ΔEshift for both the IP/NP and NP/IP cases.

TABLE II.

Fitting parameters of Sa/b(t) with various ΔEshift and Ca/b(t) in EMI+PF6 and acetonitrile.a

Solventa/bΔEshiftfGωGτβa/bΔEshiftfGωGτβ
EMI+PF6 NP/IP 0.65 8.26 22.9 0.49 IP/NP 0.64 7.66 35.7 0.33 
  2.5 0.64 8.19 23.5 0.49  2.5 0.65 7.73 40.3 0.34 
  0.63 8.15 25.9 0.50  0.66 7.82 46.6 0.35 
  7.5 0.64 8.05 30.6 0.52  7.5 0.66 7.84 50.7 0.36 
  10 0.65 7.92 32.1 0.53  10 0.66 8.05 61.6 0.36 
  15 0.64 7.82 32.9 0.52  15 0.68 8.17 88.0 0.39 
  30 0.64 7.57 46.0 0.53  30 0.74 8.41 116.1 0.47 
  45 0.68 7.88 45.6 0.55  45 0.76 9.32 137.9 0.53 
  Ca/b(t0.70 7.65 49.2 0.65  Ca/b(t0.71 9.80 75.9 0.38 
CH3CN NP/IP 0.51 13.9 1.20 0.44 IP/NP 0.51 14.4 1.13 0.43 
  2.5 0.51 13.9 1.28 0.45  2.5 0.50 14.5 1.14 0.43 
  0.52 14.2 1.46 0.46  0.51 14.8 1.22 0.43 
  7.5 0.52 14.3 1.47 0.45  7.5 0.52 14.9 1.34 0.44 
  10 0.52 14.6 1.54 0.46  10 0.52 15.5 1.35 0.44 
  15 0.52 14.2 1.63 0.46  15 0.53 15.7 1.44 0.44 
  30 0.52 14.5 2.06 0.50  30 0.56 17.0 2.32 0.48 
  45 0.56 15.8 2.74 0.54  45 0.57 17.4 2.34 0.46 
  Ca/b(t0.57 14.7 1.95 0.54  Ca/b(t0.58 17.8 2.65 0.52 
Solventa/bΔEshiftfGωGτβa/bΔEshiftfGωGτβ
EMI+PF6 NP/IP 0.65 8.26 22.9 0.49 IP/NP 0.64 7.66 35.7 0.33 
  2.5 0.64 8.19 23.5 0.49  2.5 0.65 7.73 40.3 0.34 
  0.63 8.15 25.9 0.50  0.66 7.82 46.6 0.35 
  7.5 0.64 8.05 30.6 0.52  7.5 0.66 7.84 50.7 0.36 
  10 0.65 7.92 32.1 0.53  10 0.66 8.05 61.6 0.36 
  15 0.64 7.82 32.9 0.52  15 0.68 8.17 88.0 0.39 
  30 0.64 7.57 46.0 0.53  30 0.74 8.41 116.1 0.47 
  45 0.68 7.88 45.6 0.55  45 0.76 9.32 137.9 0.53 
  Ca/b(t0.70 7.65 49.2 0.65  Ca/b(t0.71 9.80 75.9 0.38 
CH3CN NP/IP 0.51 13.9 1.20 0.44 IP/NP 0.51 14.4 1.13 0.43 
  2.5 0.51 13.9 1.28 0.45  2.5 0.50 14.5 1.14 0.43 
  0.52 14.2 1.46 0.46  0.51 14.8 1.22 0.43 
  7.5 0.52 14.3 1.47 0.45  7.5 0.52 14.9 1.34 0.44 
  10 0.52 14.6 1.54 0.46  10 0.52 15.5 1.35 0.44 
  15 0.52 14.2 1.63 0.46  15 0.53 15.7 1.44 0.44 
  30 0.52 14.5 2.06 0.50  30 0.56 17.0 2.32 0.48 
  45 0.56 15.8 2.74 0.54  45 0.57 17.4 2.34 0.46 
  Ca/b(t0.57 14.7 1.95 0.54  Ca/b(t0.58 17.8 2.65 0.52 
a

The numbers in the third and the ninth column indicate ΔEshift for Sa/b(t) and the last rows of each solvent indicate fitting parameters of Ca/b(t). The energy unit is kcal/mol and the time unit is ps.

The long time relaxation of Sa/b(t), gauged by τ and β, also generally converges to that of Ca/b(t) as ΔEshift increases. As above, this is ascribed to the decreasing perturbation with increasing ΔEshift. Interestingly, Sa/b(t) with ΔEshift = 0 tends to relax faster than Ca/b(t) at long times regardless of the solute charge states. As such, the long-time relaxation of Sa/b(t) becomes slower as ΔEshift and thus the excitation wavelength increases for both the NP/IP and IP/NP cases. This is different from short-time dynamics considered above, which show differing trends of ωG for the IP/NP and NP/IP cases. We note that the long-time relaxation of both Sa/b(t) and Ca/b(t) in acetonitrile shows a non-exponential behavior (β ≈ 0.5) according to our simulations, whereas an exponential decay was obtained with a coarse-grained model in Ref. 30. We ascribe this difference to two factors: slower diffusion in the all-atom model than in the coarse-grained model at the same temperature and the partial charges of hydrogen atoms of the methyl groups that tend to enhance short-range specific interactions.

We turn to the effects of the initial configuration and velocity distributions on nonequilibrium solvent relaxation. To estimate their respective influence, we consider three different variances53 of ΔEab(t), Δtotal(t), Δvel(t), and Δcon(t) defined as

Δtotal(t)=ΔEab(t)2x,pΔEab(t)x,p2,
Δvel(t)=ΔEab(t)2pΔEab(t)p2x,
Δcon(t)=ΔEab(t)p2xΔEab(t)x,p2=Δtotal(t)Δvel(t),
(5)

where 〈⋯〉p and 〈⋯〉x indicate a nonequilibrium average over the isoconfigurational ensemble of 50 trajectories (i.e., average over different initial velocities for a given configuration) and over 400 different initial configurations, respectively. Δtotal is the total variance of ΔEab(t) determined with the full ensemble (i.e., 20 000 trajectories), while Δvel arising from differing initial velocity assignments measures the influence of the initial velocity distribution on later dynamics. The contribution from the distribution of initial configurations, Δcon, is determined as the difference of Δtotal and Δvel. The results for the three variances for different ΔEshift are presented in Fig. 3 for EMI+PF6.

FIG. 3.

Three kinds of variance, Δtotal (left), Δvel (middle), and Δcon (right), of ΔEab(t) in nonequilibrium simulations for NP/IP (tom) and IP/NP (bottom) cases with various ΔEshift, 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol, in EMI+PF6.

FIG. 3.

Three kinds of variance, Δtotal (left), Δvel (middle), and Δcon (right), of ΔEab(t) in nonequilibrium simulations for NP/IP (tom) and IP/NP (bottom) cases with various ΔEshift, 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol, in EMI+PF6.

Close modal

In all the cases we studied, Δtotal (left panels in Fig. 3) start with 0 because the ΔEab(t) value at t = 0 is identical for all 20 000 trajectories for a given excitation energy. Their increase in time reveals that the solvent relaxation progresses in different pathways. In many cases, Δtotal reaches its maximum in less than ∼1 ps, indicating that relaxation pathways become most heterogeneous in the same time scale as ultrafast subpicosecond decay in EMI+PF6. By contrast, several cases of NP/IP relaxation are characterized by a maximum in Δtotal located near 20 ps (see below). However, even in these cases, Δtotal has either a local minimum or a shoulder structure near 1 ps. As expected, Δtotal converges to the equilibrium value of 〈(δΔEab)2〉 at long times (see Table I) regardless of the initial conditions. This demonstrates explicitly that nonequilibrium fluctuations in the long time limit become equivalent to equilibrium fluctuations.

A further insight can be gained by decomposing Δtotal into Δvel and Δcon. The results in the right panels in Fig. 3 show that Δcon reaches its maximum value at t ≲ 1 ps. This maximum, responsible for the main peak or shoulder structure observed in Δtotal, shows that Δcon is the major contributor to Δtotal at short times. While Δcon shows a rather steep decrease after maximum, Δvel gradually increases. As a result, Δvel becomes totally dominant in Δtotal at long times. These differing behaviors of Δcon and Δvel suggest that initial configurations can influence subsequent dynamics only for a short period of time. This is consistent with a previous study of a binary Lennard-Jones liquid system;54 there it was also found that initial configurations cannot predict long time dynamics. One salient aspect of Δvel is that for ΔEshift ≲ 30 kcal/mol, it is a non-monotonic function of time with a maximum occurring near t = 100 ps in the NP/IP case. This is ascribed to the absence of the solute-solvent Coulomb interactions with the NP solute charge distribution. Specifically, abrupt elimination of the IP-solvent Coulomb interactions with the initial excitation in the NP/IP case subjects solvent ions to a strong nonequilibrium force, which helps the solvent to explore a large “volume” of relaxation pathways, i.e., a large degree of variations in relaxation pathways. As ΔEshift increases, the initial solvent environment becomes less “IP-like” and the nonequilibrium force the solvent feels after FC excitation becomes reduced. As a result, the volume of reaction pathways that can be explored and hence Δvel decrease. The variances have generally lower values with greater ΔEshift, suggesting that the degree of variations in relaxation pathways becomes smaller as the initial nonequilibrium state is generated closer to the fully relaxed equilibrium state. This means that at given time, the linewidth of the fluorescence emission spectrum decreases as the excitation wavelength increases. It would be interesting to see if this behavior would be borne out in experiments.

In this section, we study the influence of solute dipole moment on solvation dynamics by considering the QP and HP charge distributions (see Sec. II above). The diatomic solute model employed in the present work has been extensively utilized over the past decade18–23,30,35,52 because solvation dynamics in room temperature ionic liquids are well captured through this model. However, since the dipole difference of the IP and NP states used here is somewhat larger than that of typical dye molecules, it is worthwhile to investigate how the dipole difference of the ground and excited states of the probe solute affects its solvation dynamics. To do so, we analyzed Ca/b(t) and Sa/b(t) for the QP (4.2 D) and HP (8.4 D) solutes that have a smaller dipole moment than the IP (16.8 D) solute. Ca/b(t) and Sa/b(t) for the QP/NP and HP/NP cases were obtained by Eqs. (2) and (3); 10 000 initial configurations selected from those used for the IP/NP and NP/IP cases were employed for Sa/b(t) calculations. We note that Ca/b(t) for the NP/QP and NP/HP cases are exactly the same as that for the NP/IP case because equilibrium dynamics occur in the presence of the same NP state in all three cases. Sa/b(t) for the NP/QP and NP/HQ cases are similar to Sa/b(t) for NP/IP with ΔEshift = 45 and 30 kcal/mol, respectively.

The results for Ca/b(t) and Sa/b(t) for QP/NP and HP/NP are exhibited in Fig. 4. For comparison, the results for the IP/NP case are also shown there. Overall solvation dynamics do not vary strongly with the dipole moment of the active state a. Despite some minor differences, especially in Ca/b(t), all three dipolar states, IP, HP, and QP, show very similar relaxation behaviors, including the biphasic character and time scale of Ca/b(t) and Sa/b(t) dynamics. Furthermore, according to a recent study,55 a diatomic solute similar to ours captures main features of solvation dynamics of a dye molecule, such as coumarin 153, correctly. These findings indicate that the IP/NP solute system provides a very reasonable framework to study solvation dynamics in RTILs and our analysis based on it is robust.

FIG. 4.

Comparison of the normalized time-correlation function, Ca/b(t), and the Stokes-shift function, Sa/b(t), for QP/NP (left), HP/NP (middle), and IP/NP cases (right) in EMI+PF6.

FIG. 4.

Comparison of the normalized time-correlation function, Ca/b(t), and the Stokes-shift function, Sa/b(t), for QP/NP (left), HP/NP (middle), and IP/NP cases (right) in EMI+PF6.

Close modal

Finally, we examine the steady-state emission spectrum of the diatomic solute given by45 

Ia/b(ΔEem)=0δ(Eb(t)Ea(t)ΔEem)×et/τfdt,
(6)

where τf is the fluorescence lifetime of the solute and a and b denote its ground and excited states. Since solvent relaxation occurs in the presence of the solute in state b in Eq. (6), solvent dynamics relevant to Ia/b are given by Sa/b(t).

The results determined with various ΔEshift are displayed in Fig. 5. As mentioned above, the larger the ΔEshift value, the lower the excitation energy and the closer the initial nonequilibrium FC state to the free energy minimum of the excited state. Each spectrum is obtained by averaging over 20 000 trajectories which have equal excitation energy and area normalized. The lifetime of the solute is arbitrarily chosen as 100 ps, which roughly corresponds to the lifetime of 2-amino-7-nitrofluorene (ANF) probe molecule.43,45

FIG. 5.

The steady-state emission spectra for NP/IP (left) and IP/NP (right) cases with various ΔEshift, 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol, in EMI+PF6 (top) and acetonitrile (bottom).

FIG. 5.

The steady-state emission spectra for NP/IP (left) and IP/NP (right) cases with various ΔEshift, 0, 2.5, 5, 7.5, 10, 15, 30, and 45 kcal/mol, in EMI+PF6 (top) and acetonitrile (bottom).

Close modal

The emission spectra show significant dependence on the absorption energy in EMI+PF6 for both NP/IP and IP/NP, whereas they are nearly identical in acetonitrile. Our extensive nonequilibrium simulation study based on a large number of trajectories thus provides convincing evidence for the presence of REE in RTIL systems, in good agreement with Ref. 45.

In the present article, we studied solvation dynamics of a probe solute and related fluorescence spectra in EMI+PF6 with a special focus on the effect of the excitation energy by performing MD simulations with a coarse-grained model description. Despite differences in potential models, our results for overall solvation dynamics share common features with previous studies.19,30 Specifically, the time correlation function Ca/b(t) of solvation energy fluctuations and the time-dependent Stokes shift function Sa/b(t) are characterized by two regimes, i.e., ultrafast relaxation followed by a slow decay. These functions show deviations at sort times but become nearly congruent at long times in EMI+PF6. By contrast, Ca/b(t) and Sa/b(t) in acetonitrile obtained with an all-atom potential model show a good agreement in the entire time range we studied.

One of our main findings was that the more the excitation energy is red-shifted, the slower the long time relaxation of Sa/b(t) becomes in EMI+PF6. Furthermore, Sa/b(t) gradually converges to Ca/b(t) as ΔEshift increases. The solvent frequency ωG and the time-scale τ, gauging the short time and the long time dynamics, respectively, become closer to those of Ca/b(t), complying with the linear response theory. We have also examined the effect of the initial configurations and velocity distributions on the nonequilibrium solvent relaxation by analyzing the variances of 20  000 different trajectories but with same initial excitation energy. The contribution of the initial configuration to the variance diminishes while that of the initial velocity grows as time evolves. We found that the steady-state emission spectrum varies with the excitation energy in EMI+PF6 due to the slow relaxation of the solvent compared with the fluorescence lifetime of the probe.

One novel aspect of our analysis is that we employed iso-excitation energy ensemble, i.e., we constructed initial configurations that yield the same excitation energy. This enabled us to examine the dependence of nonequilibrium dynamics and emission spectrum on the excitation energy, in contrast to previous works where the excitation energies are sampled randomly.19,30,45 In addition, we simulated a large number of trajectories, 20 000 for each ΔEshift, to study nonequilibrium relaxation, adopting the isoconfigurational ensemble approach. This not only yields excellent statistics but also permits analysis of variances of relaxation dynamics arising from initial structure and dynamics. The use of a computationally-efficient coarse-grained RTIL model made this analysis possible (the total cumulative simulation time for RTIL was 320 μs). To our knowledge, the present study is the first simulation to carry out excitation-energy-selective nonequilibrium solvation dynamics using the isoconfigurational ensemble method. Moreover, our study demonstrates that using the concept of dynamics propensity, the isoconfigurational ensemble method provides useful insight into solvation dynamics in a way similar to prior studies of dynamic heterogeneity of supercooled liquid51 and RTIL systems.5 

In the present study, we used the fluorescence lifetime of 100 ps, pertinent to ANF probe molecules. It would be desirable in the future to extend this by considering a range of different solute lifetimes to examine the solute dependence of the emission spectrum and REE. It would be also worthwhile to study the time dependent emission spectra to investigate the time scale at which the REE appears distinctly.

This work was supported in part by the National Science Foundation through NSF Grant No. CHE-1223988. This work was also supported in part by National Research Foundation (NRF) grants funded by the Korean Government (MEST) (Grant Nos. NRF-2009-0070517 and NRF-2014R1A1A2056119), and by grants funded by the KISTI Supercomputing Center (Grant Nos. KSC-2011-C1-15 and KSC-2012-C2-47). Y.S. acknowledges financial support from the BK 21 Program of Korea.

1.
N. V.
Plechkova
and
K. R.
Seddon
,
Chem. Soc. Rev.
37
,
123
(
2008
).
2.
J. P.
Hallett
and
T.
Welton
,
Chem. Rev.
111
,
3508
(
2011
).
3.
H.
Weingärtner
,
Angew. Chem., Int. Ed.
47
,
654
(
2008
).
4.
D.
Jeong
,
M. Y.
Choi
,
H. J.
Kim
, and
Y.
Jung
,
Phys. Chem. Chem. Phys.
12
,
2001
(
2010
).
5.
D.
Kim
,
D.
Jeong
, and
Y.
Jung
,
Phys. Chem. Chem. Phys.
16
,
19712
(
2014
).
6.
S.-W.
Park
,
S.
Kim
, and
Y.
Jung
,
Phys. Chem. Chem. Phys.
17
,
29281
(
2015
).
7.
R.
Karmakar
and
A.
Samanta
,
J. Phys. Chem. A
106
,
4447
(
2002
).
8.
R.
Karmakar
and
A.
Samanta
,
J. Phys. Chem. A
106
,
6670
(
2002
).
9.
R.
Karmakar
and
A.
Samanta
,
J. Phys. Chem. A
107
,
7340
(
2003
).
10.
J. A.
Ingram
,
R. S.
Moog
,
N.
Ito
,
R.
Biswas
, and
M.
Maroncelli
,
J. Phys. Chem. B
107
,
5926
(
2003
).
11.
S.
Arzhantsev
,
H.
Jin
,
N.
Ito
, and
M.
Maroncelli
,
Chem. Phys. Lett.
417
,
524
(
2006
).
12.
B.
Lang
,
G.
Angulo
, and
E.
Vauthey
,
J. Phys. Chem. A
110
,
7028
(
2006
).
13.
H.
Jin
,
G. A.
Baker
,
S.
Arzhantsev
,
J.
Dong
, and
M.
Maroncelli
,
J. Phys. Chem. B
111
,
7291
(
2007
).
14.
A.
Paul
and
A.
Samanta
,
J. Phys. Chem. B
111
,
4724
(
2007
).
15.
A.
Paul
and
A.
Samanta
,
Indian J. Chem., Sect. A: Inorg., Bio-inorg., Phys., Theor. Anal. Chem.
49
,
649
(
2010
).
16.
A.
Samanta
,
J. Phys. Chem. Lett.
1
,
1557
(
2010
).
17.
X.-X.
Zhang
,
M.
Liang
,
N. P.
Ernsting
, and
M.
Maroncelli
,
J. Phys. Chem. B
117
,
4291
(
2013
).
18.
Y.
Shim
,
J.
Duan
,
M. Y.
Choi
, and
H. J.
Kim
,
J. Chem. Phys.
119
,
6411
(
2003
).
19.
Y.
Shim
,
M. Y.
Choi
, and
H. J.
Kim
,
J. Chem. Phys.
122
,
044511
(
2005
).
20.
Y.
Shim
and
H. J.
Kim
,
J. Phys. Chem. B
112
,
11028
(
2008
).
21.
D.
Jeong
,
M. Y.
Choi
,
Y.
Jung
, and
H. J.
Kim
,
J. Chem. Phys.
128
,
174504
(
2008
).
22.
Y.
Shim
and
H. J.
Kim
,
J. Phys. Chem. B
114
,
10160
(
2010
).
23.
Y.
Shim
and
H. J.
Kim
,
J. Phys. Chem. B
117
,
11743
(
2013
).
24.
M. N.
Kobrak
and
V.
Znamenskiy
,
Chem. Phys. Lett.
395
,
127
(
2004
).
25.
M. N.
Kobrak
,
J. Chem. Phys.
125
,
064502
(
2006
).
26.
M. N.
Kobrak
,
J. Chem. Phys.
127
,
184507
(
2007
).
27.
C. J.
Margulis
,
Mol. Phys.
102
,
829
(
2004
).
28.
B. L.
Bhargava
and
S.
Balasubramanian
,
J. Chem. Phys.
123
,
144505
(
2005
).
29.
Y.
Zhao
and
Z.
Hu
,
Chem. Commun.
48
,
2231
(
2012
).
30.
D.
Roy
and
M.
Maroncelli
,
J. Phys. Chem. B
116
,
5951
(
2012
).
31.
Z. L.
Terranova
and
S. A.
Corcelli
,
J. Phys. Chem. B
117
,
15659
(
2013
).
32.
M.
Schmollngruber
,
C.
Schröder
, and
O.
Steinhauser
,
J. Chem. Phys.
138
,
204504
(
2013
).
33.
M.
Schmollngruber
,
C.
Schröder
, and
O.
Steinhauser
,
Phys. Chem. Chem. Phys.
16
,
10999
(
2014
).
34.
Y.
Shim
and
H. J.
Kim
,
J. Phys. Chem. B
111
,
4510
(
2007
).
35.
Y.
Shim
and
H. J.
Kim
,
J. Phys. Chem. B
113
,
12964
(
2009
).
36.
R. M.
Lynden-Bell
,
J. Phys. Chem. B
111
,
10800
(
2007
).
37.
H.
Jin
,
X.
Li
, and
M.
Maroncelli
,
J. Phys. Chem. B
111
,
13473
(
2007
).
38.
A.
Adhikari
,
K.
Sahu
,
S.
Dey
,
S.
Ghosh
,
U.
Mandal
, and
K.
Bhattacharyya
,
J. Phys. Chem. B
111
,
12809
(
2007
).
39.
A.
Adhikari
,
S.
Dey
,
D. K.
Das
,
U.
Mandal
,
S.
Ghosh
, and
K.
Bhattacharyya
,
J. Phys. Chem. B
112
,
6350
(
2008
).
40.
Y.
Kimura
,
K.
Suda
,
M.
Shibuya
,
Y.
Yasaka
, and
M.
Ueno
,
Bull. Chem. Soc. Jpn.
88
,
939
(
2015
).
41.
P. K.
Mandal
,
M.
Sarkar
, and
A.
Samanta
,
J. Phys. Chem. A
108
,
9048
(
2004
).
42.
A.
Paul
,
P. K.
Mandal
, and
A.
Samanta
,
J. Phys. Chem. B
109
,
9148
(
2005
).
43.
A.
Samanta
,
J. Phys. Chem. B
110
,
13704
(
2006
).
44.
Z.
Hu
and
C. J.
Margulis
,
J. Phys. Chem. B
110
,
11025
(
2006
).
45.
Z.
Hu
and
C. J.
Margulis
,
Proc. Natl. Acad. Sci. U. S. A.
103
,
831
(
2006
).
46.
Z.
Hu
and
C. J.
Margulis
,
Acc. Chem. Res.
40
,
1097
(
2007
).
47.
A. M.
Nikitin
and
A. P.
Lyubartsev
,
J. Comput. Chem.
28
,
2020
(
2007
).
48.
S.
Pronk
,
S.
Páll
,
R.
Schulz
,
P.
Larsson
,
P.
Bjelkmar
,
R.
Apostolov
,
M. R.
Shirts
,
J. C.
Smith
,
P. M.
Kasson
,
D.
van der Spoel
,
B.
Hess
, and
E.
Lindahl
,
Bioinformatics
29
,
845
(
2013
).
49.
A.
Warshel
,
J. Phys. Chem.
86
,
2218
(
1982
).
50.
G.
King
and
A.
Warshel
,
J. Chem. Phys.
93
,
8682
(
1990
).
51.
A.
Widmer-Cooper
,
P.
Harrowell
, and
H.
Fynewever
,
Phys. Rev. Lett.
93
,
135701
(
2004
).
52.
Y.
Shim
,
M. Y.
Choi
, and
H. J.
Kim
,
J. Chem. Phys.
122
,
044510
(
2005
).
53.
L.
Berthier
and
R. L.
Jack
,
Phys. Rev. E
76
,
041509
(
2007
).
54.
J. A.
Rodriguez Fris
,
L. M.
Alarcón
, and
G. A.
Appignanesi
,
Phys. Rev. E
76
,
011502
(
2007
).
55.
E. C.
Wu
and
H. J.
Kim
,
J. Phys. Chem. B
120
,
4644
(
2016
).