A Simulation Study of the Fundamental Vibrational Shifts of HCl Diluted in Ar, Kr, and Xe: Anharmonic Corrections Effects

We have calculated the vibrational solvent shifts of the fundamental bands of HCl diluted in Ar, Kr, and Xe solutions at different thermodynamic conditions by means of the molecular dynamics technique and a model for the isotropic part of the interaction depending on the vibration. The theoretical vibrational shifts, which were compared with the available experimental data, have been determined by considering both, the usual linear Buckingham terms and the nonlinear anharmonic corrections, and the latter omitted in a previous work for the HCl in Ar and Kr. We have found that the Buckingham contributions dominate the solvent shifts of the fundamental bands of HCl in Ar, Kr, and Xe, although the anharmonic shifts’ present significant greater values than those obtained previously for N2 diluted in liquid Ar and pure liquid N2, both at normal conditions. We have analyzed the solvent shifts influence of the linear and quadratic (in the vibrational coordinate) oscillator-bath interaction terms and also the Dunham intramolecular potential effects on the anharmonic contributions.


Introduction
The study of the fundamental and overtones infrared and Raman vibrational solvent shifts gives relevant information about the molecular interaction depending of the intramolecular vibration [1][2][3][4][5][6][7]. In an earlier work, according to Buckingham, it was deduced that the diatomic solvent shift associated to the vibrational transition 0 → ν of a cubic anharmonic oscillator presents a linear dependence with the vibrational quantum number ν [1]. This result, which is valid in principle for the fundamental band or lower overtones of several diatomics, would break down for the fundamental bands of some molecules and high overtones of much ones. Recently, Alessi et al. [6,7] has deduced the non-linear (ν (ν + 1)) correction (anharmonic correction) to the linear Buckingham solvent shift for a diatomic molecule perturbed by an atomic or molecular solvent. Alessi et al. [6,7] analyzed the solvent shift of N 2 diluted in liquid Ar and pure N 2 liquid, both at normal conditions, founding that the non-linear solvent correction, also called anharmonic shift [6,7], represents only the 1% of the fundamental shift of the isotropic Raman band. This result is due in part to the small anharmonicity parameter χ e of N 2 (6.0732 · 10 −3 ) [6,7]; however, because the HCl anharmonicity constant χ e is clearly much larger than the N 2 one, we expect a major contribution of non-linear solvent shift for the fundamental band of HCl diluted in Ar, Kr, and Xe.
We have made a previous study [8] of the vibrational solvent shifts of HCl diluted in Ar and Kr considering the usual linear Buckingham approach. In the present work we will study the solvent shifts of HCl also in Ar and Kr and additionally in liquid Xe by considering the non-linear corrections to the solvent shift deduced by Alessi et al. [6,7]. Like in the previous work [8] both the diatomic vibration and rotation will be treated with a quantum approach (small system S), while the translational degree of freedom, both of the diatomic and the solvent molecules, will be treated classically (the bath B). Unlike the previous study [8] where the vibration was modeled by means of a cubic anharmonic oscillator, in this work the vibration will be characterized by using a fifth-order Dunham anharmonic oscillator [6,7]. Fourth-and fifth-order terms of the intramolecular potential must be considered because they have relevant effects on the anharmonic solvent shift. On the other hand, 2 International Journal of Spectroscopy classical evolution of the diatomic and solvent translational coordinates (the bath B) will be represented by using molecular dynamic simulation techniques for fixed values of temperature and density.
A key element in the solvent shift calculation is to know the functional form of the isotropic part of diatomic-solvent interaction potential depending on the vibration. Like in the previous work [8] we use for such part of the HCl-rare atom potential the model proposed by Marteau el al. [9], which depends only on a fitting parameter, that will be determined by comparing the theoretical and the experimental solvent shifts, at the different thermodynamic states in which the HCl absorption bands were measured.
For the solvent shifts of HCl diluted in dense Ar we have used the same site-site HCl-Ar potential of the previous work [8] (obtained by fitting the site-site potential to the functional form proposed by Holmgren et al. [10]), while for HCl diluted in Kr and Xe we have employed the sitesite HCl-rare atom potentials obtained by fitting the sitesite functional form to the accurate potentials proposed by Hutson and Howard [11].
In the present study we have analyzed the relative contributions of the linear Buckingham shift and the nonlinear anharmonic term to the vibrational solvent shift of HCl in Ar, Kr, and Xe, for the different thermodynamic states considered in the experimental study [12]. Additionally, we have analyzed the contributions to the solvent shift of the linear and quadratic (in the diatomic vibrational coordinate) parts of the oscillator-bath interaction. The effects on the anharmonic solvent shift of the third, fourth, and fifth Dunham potential terms also are discussed.
The organization of the paper is as follows. In Section 2 we develop the basic theoretical aspects related with the vibrational solvent shifts. In Section 3 we introduce the model for the isotropic part of the HCl-rare atom interaction depending of the vibration, the fitting of the site-site potentials, and the details of the molecular dynamics simulations. In Section 4 we compare the theoretical and the experimental solvent shift and discuss different aspects related with this magnitude. Finally, a brief summary and some conclusions are collected in Section 5.

Theoretical Background
In order to deduce the basic expressions of the solvent shift we make reference to a spectral theory previously developed [13], although the results obtained are more general than the proper theory. This theory [13] considers the evolution of a small system S, formed by the diatomic rotation and vibration, which interacts with a thermal bath B, constituted by the translational degree of freedom, both of the diatomic and solvent molecules. In this theory it is assumed that the bath average of the S-B interaction H is zero. This condition always can be achieved by redefining the S-B interaction H and the system S Hamiltonian H S as H − H , and H S + H respectively. It is precisely the new Hamiltonian H S + H who defines the frequencies of the vibrorotational spectral lines, introducing the term H the frequency shifts which are the object of our study.
The S-B interaction Hamiltonian H can be divided into the isotropic V I and the anisotropic V A parts: Because the bath is macroscopically isotropic, in lowest order approximation, when the anisotropic intermolecular forces are neglected in the bath average, the mean value of the anisotropic part of the interaction is zero V A = 0 [14]; therefore, in this approach, the frequency shift produced by the solvent is basically due to the isotropic part of the interaction H S + V I and hence to the diatomic vibration. Although Alessi et al. [6,7] included the anisotropic force effects on bath averages in their study of pure N 2 and N 2 diluted in liquid Ar, we have neglected these effects for HCl in rare gases in the present study; they will be considered in a future work. The isotropic part of the S-B interaction can be written as the sum of the binary isotropic potentials depending of the vibration v (Q, r p ): where Q is the vibrational coordinate and r p is the norm of the vector joining the center of mass of the diatomic and the p-solvent molecule. The binary potential will be approximated by means of the first two terms of its power expansion in the vibrational coordinate: so the isotropic S-B interaction can be written as where with n = 1, 2. Then, the bath average of the isotropic part of the S-B potential (4) is given by where F n is the bath average of the coupling coefficient (5). In [7] also the cubic and quartic terms of V I were considered, but it was found that their contributions to the solvent shift were small. The vibrorotational Hamiltonian H S is modeled as a rigid rotor H R coupled H RV to an anharmonic oscillator H V : The anharmonic oscillator is represented initially by the usual Morse model: International Journal of Spectroscopy 3 where P Q is the vibrational momentum, μ is the reduced mass of the diatomic, and D and a are the parameters of the Morse potential. The energy spectrum of the Morse Hamiltonian is [7] where the harmonic frequency is and the anharmonicity parameter is The expression (9) is usually used to fit spectroscopy data, by means of which obtain the values of the vibrational harmonic frequency ω h and the anharmonicity ω h χ e . When the Morse potential (8) is expanded in powers of the vibrational coordinate Q up to the fifth order we obtain the Dunham model of the anharmonic oscillator [7]: where In order to relating more easily the perturbation coefficient f j with the vibrational spectroscopy constants ω h and χ e , we rearrange the expressions ((10), (11)) as [7] By making second-order quantum mechanical perturbation theory calculations with the Hamiltonian H V + V I ( (6), (12)) we obtain the following expression for the solvent shift associated to the vibrational transition 0 → ν [7]: where Δω B 0 → ν is the linear solvent shift predicted by Buckingham's theory [1], while Δω A 0 → ν is the non-linear solvent shift, also called anharmonic shift, deduced by Alessi et al. [7]. Using the relations ( (13), (14)) we can write the linear and non-linear solvent shifts ( (16), (17)) in terms of the spectroscopic constants ω h and χ e , that is, defining the coefficients: where R e is the equilibrium intermolecular distance of the diatomic; the solvent shifts can be expressed as where B e is the diatomic rotational constant. The first summand in (20) corresponds to the two f 3 f 4 terms of F 1 in (17), the second summand corresponds to the f 5 term in (17), and the last summand corresponds to the two f 3 f 4 terms of F 2 in (17). The first term of (20) is the result of the large cancelation effect of the two f 3 f 4 terms of F 1 (17) [7], being the unique summand that presents opposite sign to F n .

Interaction Potentials and Simulations
The HCl-solvent binary interaction potential depending on the vibration is a key element in the calculation of the vibrational solvent shift of the diatomic molecule. Not many models for this part of the isotropic interaction depending of the vibration are available in the literature. For this kind of potential it is usual to consider site-site models which take into account the presence of the vibrational coordinate in the site-site distances, also it is assumed that the Lennard-Jones parameters and σ present a linear dependence with the vibrational coordinate [5][6][7]. However, in this work, like in the previous one [8], we consider the model introduced by Marteau et al. [9] for the binary isotropic interaction potential depending on the interaction. By writing the vibrational components v (n) (r) (3) as the binary isotropic components of Marteau et al. can be expressed: where and σ are the usual 6-12 Lennard-Jones parameters of the diatomic-solvent interaction. S 2 and T 2 are given by International Journal of Spectroscopy  where α is the diatomic polarizability and Q = Q/R e , while S 1 and T 1 are given by where f r is a phenomenological parameter, called the repulsive factor, which will be obtained by comparing the theoretical solvent shift with the experimental data [13]. Unless f r , all the parameters involved in the isotropic components (22) are collected in Table 1.
For the Molecular dynamics simulations we consider a single HCl molecule and 250 solvent atoms collected in a cubic box with periodic boundary conditions. In each case the size of the box was chosen to give the experimental solvent density ρ. The motion equations were solved by using a Leap-Frog-Verlet algorithm which includes a SHAKE procedure to fix the diatomic bond length to R e = 1.275Å and a coupling to a thermal bath to keep the system temperature near to the experimental solvent value T. The time step was chosen as Δt = 2.5 fs, including each simulation run with a total of 10 6 time steps. Previous to each simulation run, a long equilibration period was considered.
As it is usual, the Molecular dynamics simulations were carried out considering that the interaction between the solvent atoms can be described by means of the 6-12 Lennard-Jones potentials, in our case with the values of the parameters σ = 3.405Å, = 119.8 K for Ar-Ar, σ = 3.575Å, = 191.4 K for Kr-Kr, and σ = 3.89Å, = 282.4 K for Xe-Xe [16].
Also, as it is usual in the Molecular dynamics simulations, the binary interaction between the diatomic molecule and the solvent atom will be modeled by means of a sitesite potential. Because actually there are more accurate potential models available for the HCl-Ar (Kr,Xe) interaction [10,11]; the site-site potentials were fitted to these more precise potential models. For the HCl-Ar we adopt the result obtained by Medina et al. [16] which fitted the site-site model to the Holmgren et al. potential [10]; the values obtained for the site-site parameters are shown in Table 2. Logically, the results of the fitting depend of the interatomic interval considered for the potential, so as to have similar fitting conditions to the HCl-Ar. In the HCl-Kr case we have fitted the Hutson and Howard potential [11] considering the r interval [3.5,10.5]. For the HCl-Xe case we have fitted the Hutson potential [11] by considering the [3.7,11] interval. In Table 2 are shown the values of the site-site parameters obtained for HCl-Kr and HCl-Xe with this fitting procedure.
In the Molecular dynamics simulations it is necessary to specify the values of the solvent density and temperature. However, our experimental information is referred to the pressure and temperature [12], so the density must be obtained by an indirect method. For the liquid-vapor coexistence states the densities were determined by means of a cubic splines interpolation of the experimental data (ρ and T) reported by Bertsev and Zelikina [18]. While for the supercritical states the densities were calculated by using the vdW-711 equation of state [19]. The experimental values of the pressure and temperature and the values obtained for the density are shown in Tables 3, 4, and 5 for Ar, Kr, and Xe, respectively. In all the cases the superscripts in the temperatures denote supercritical states, while the remaining ones are in the liquid-vapor coexistence curve.

Vibrational Solvent Shifts
The experimental spectral profile of the fundamental absorption bands of HCl diluted in Ar, Kr, and Xe along the liquidvapor coexistence curve and in supercritical states has been reported by Pérez et al. [12]. From these absorption line shapes the band origin frequencies have been determined by fitting the vibrorotational spectral line frequencies of the rotational resonances [12]. From the band origin frequencies we have obtained the vibrational solvent shifts Δω 0 → 1 by using the isotope abundance-weighted HCl molecular constants [20] ω h = 2990.4 cm −1 and χ e = 0.01766; also we have considered the weighted value of the rotational constant B e = 10.589 cm −1 [20] which has been used in the theoretical calculations of the solvent shift ((19), (20)). The values of the experimental solvent shifts Δω 0 → 1 obtained in this way are collected in Tables 3, 4, and 5 for HCl-Ar, HCl-Kr, and HCl-Xe, respectively.
International Journal of Spectroscopy 5   We have calculated the theoretical solvent shifts Δω 0 → 1 of HCl by using ((15), (19), and (20)) calculating the bath averages F n (n = 1, 2) of the coupling coefficients (5) by means of the molecular dynamics simulations and employing the values of the repulsive factor f r (24) which give the best fit of the theoretical to the experimental solvent shifts Δω 0 → 1 . The repulsive factors obtained in the fitting procedure are f r = 2.1 for the Ar, f r = 1.9 for the Kr, and f r = 1.7 for the Xe. These values are greater than those obtained with the Buckingham cubic model of the previous paper [8], a 18% major for the HCl-Ar, and a 19% for the HCl-Kr. The non-linear Dunham model, ((15), (19), and (20)), produces red solvent shifts Δω 0 → 1 greater in magnitude than those of the linear Buckingham model of the previous paper [8], so it is necessary to increase the repulsive factor f r , in order to produce a decrease of the magnitude of the vibrational shifts Δω 0 → 1 and obtain theoretical solvent shifts near to the experimental values.
The theoretical solvent shifts Δω 0 → 1 of the HCl in Ar, Kr, and Xe solutions are collected in Tables 3, 4, and 5, respectively. Also, in Figure 1 we have plotted the theoretical and the experimental solvent shifts Δω 0 → 1 of HCl in Ar, Kr, and Xe for different values of the temperature. As it can be appreciated in general terms there is a good agreement between the theoretical and the experimental solvent shifts; the better accord is found for HCl-Ar where the theoretical experimental differences are less than 10% for all the temperatures. For the HCl-Kr these differences are less than 10% unless for its two maxima temperatures, while for the HCl-Xe the differences are less than 10% unless for its three maxima temperatures. In all the cases the vibrational shift Δω 0 → 1 is a red shift which increases its magnitude following the sequence Ar, Kr, and Xe, that is, the sequence in which increasing the diatomic-atomic interaction (see Table 1). On the other hand, we can observe that for the three solvents the vibrational shifts Δω 0 → 1 decrease their magnitude with increasing temperature and decreasing density, or in other words, with increasing mean interatomic distance. In all the cases studied in this work the non-linear anharmonic shift Δω A 0 → 1 (20) contributes an 8% to the solvent shift Δω 0 → 1 (15), so for the HCl in rare gases the linear Buckingham approximation Δω B 0 → 1 (19) dominates with a 92% of the solvent shift of the fundamental absorption bands. In Figure 2 we have plotted the theoretical solvent shifts Δω 0 → 1 together with the linear (19) and non-linear (20) contributions, and as can be appreciated, the Buckingham term is the principal contribution to the fundamental solvent shifts Δω 0 → 1 . However, the anharmonic non-linear terms for the HCl in rare gases are greater than the corresponding values Δω A 0 → 1 for pure N 2 liquid and N 2 diluted in liquid Ar at normal conditions, where the nonlinear term only represents the 1% of the vibrational shifts [7]. These differences are due in part to the anharmonicity parameter χ e [7] which is smaller for the N 2 molecule. So for the first ten overtones of N 2 the solvent shifts Δω 0 → ν present a relatively linear behavior [7], such as predicting the Buckingham theory [1]. However, for the first ten overtones of the HCl in rare gases the solvent shifts present a non-linear behavior such as is what shown in Figure 3, where we have plotted the vibrational shifts Δω 0 → ν of the HCl in Ar together with the Buckingham approach Δω B 0 → ν . In this figure we can observe as the solvent shift diverge of the Buckingham linear theory conforms the quantum number ν is increased, so we can verify that for high overtones the anharmonic terms Δω A 0 → ν make very relevant contributions to the solvent shifts. In this way, we must take into account the non-linear solvent shifts in the study of the HCl vibrational spectra, particularly for the vibrational overtones, as is the case of the first overtone absorption bands of HCl in fluid SF 6 , for which there are experimental data at different thermodynamics conditions [21].
In contrast to pure N 2 liquid and N 2 in liquid Ar where the F 1 term is negative and the F 2 term is positive [7], for the HCl in rare gases both F 1 and F 2 terms ( (19), (20)) are negative, so in the Buckingham approach the F 1 terms represent the 71% of the linear contribution Δω B 0 → 1 , while the F 2 terms represent the 29% of the of the Buckingham solvent shift. For the non-linear anharmonic shift Δω A 0 → 1 the F 1 terms represent the 33% and the F 2 terms represent the 67%. In this last case the f 3 f 4 positive contributions of F 1 represent the 2% of Δω A 0 → 1 , the f 5 contribution of F 1 represents the 36% of Δω A 0 → 1 , and finally the f 3 f 4 terms of F 2 represent 67% of Δω A 0 → 1 .

Summary and Conclusions
We have calculated the vibrational solvent shifts of the fundamental absorption bands of the HCl diluted in Ar, Kr, and Xe by means of the molecular dynamics simulation and the model proposed by Marteau et al. [9] for the isotropic part of the interaction depending of the vibration. The Marteau potential depends of a phenomenological parameter, the repulsive factor f r , which has been obtained by comparing the theoretical and the experimental solvent shifts. In this work, in contrast with the previous paper [8], we have taken into account the non-linear anharmonic contribution to the solvent shift, obtaining a reasonable agreement between the theoretical and the experimental solvent shifts. We have found that for HCl in rare gases the non-linear anharmonic term Δω A 0 → 1 makes a modest contribution to the fundamental solvent shift; however, this contribution is more important than that observed in pure N 2 liquid and N 2 diluted in Ar liquid [7]. On the other hand, for medium and high overtones the anharmonic shifts make an important contribution to the vibrational solvent shifts of HCl in rare gases.
Finally, we have quantified the effects on the vibrational solvent shifts of the linear and quadratic terms of the isotropic part of interaction depending on the vibration, and of the Dunham terms on the anharmonic solvent shifts.