Orientation correlation and local field in liquid nitrobenzene

Analysis¹ of a recent hyper-Rayleigh scattering (HRS, or second-harmonic light scattering) experiment for water² finds that the transverse and longitudinal orientation correlation functions for the water molecules at long range are almost completely determined from the observed HRS polarization and angle dependence. Such direct experimental information is important because intermolecular orientation correlation of dipolar molecules has a strong effect on the dielectric properties of simple polar liquids and more complex systems.^3,4 HRS, mediated by the molecular first hyperpolarizability tensor β,^5,6 is sensitive to these orientation correlations since the vector part of β transforms the same as the permanent molecular dipole moment that makes the dominant contribution to the dielectric response of polar liquids.

The dielectric response of a liquid composed of dipolar and polarizable molecules is a problem of long standing,^7–11 which is not entirely solved even now.^12,13 Difficulties arise due to the long range of the dipole interaction, the short-range and long-range orientation correlations of the molecules, and the non-locality introduced into the local field calculation by the molecular polarizability. The static dielectric constant is sensitive to all these factors. Onsager considered a dipolar molecule in a cavity of molecular size in the liquid continuum and introduced the distinction between the reaction field and the directing field acting on the molecule.⁷ Kirkwood then replaced Onsager’s cavity with a mesoscopic spherical cavity containing the molecule and all its neighbours, and treated the cavity contents by statistical mechanics.⁸ The Kirkwood-Onsager relation^3,9 for the static dielectric constant ε_s

combines continuum electrostatics and statistical mechanics. It is expressed in terms of the high frequency dielectric constant ε_∞, molecular number density ρ, permanent molecular dipole moment μ, the Kirkwood correlation factor g_K, and temperature T, where ε₀ is the vacuum permittivity and k_B is the Boltzmann constant. The Kirkwood factor

accounts for the short-range intermolecular orientation correlations (θ_j1 is the angle between the dipoles on molecules j and 1), and ε_∞ is the contribution at zero frequency due to the librational, vibrational, and electronic polarizabilities of the molecules, usually estimated from ε(ω) measured at a frequency too high for molecular orientation relaxation, or by extrapolating the optical refractive index. Eq. (1) is not exact for polar, polarizable molecules, but it is rigorous for non-polarizable molecules where ε_∞ = 1, and it reduces to ε_s = ε_∞ for non-polar molecules. The dipole moment of a molecule in the liquid differs from that of a molecule in the gas phase due to the molecular dipole induced by the electric field of the polarized surrounding liquid acting on the molecular polarizability, and this effect is accounted for by ε_∞ in Eq. (1). However, the permanent dipole μ appearing in Eq. (1) also differs from the gas phase dipole μ₀ as a result of distortion of the molecular geometry and electronic structure by short range forces, higher multipole fields, and electronic overlap with neighbour molecules, the effects of which are not accounted for in the derivation of Eq. (1). Both μ and g_K are functions of the molecular environment and are not separately determined by dielectric measurements.

Molecules in the liquid respond to the microscopic local electric field and not the macroscopic Maxwell field. At high frequencies or for non-polar liquids, the factor f relating the local and applied fields is the Lorentz local field factor

whereas the Onsager local field factor

is widely used for polar liquids at low frequencies.^14,15 Eq. (4) is obtained using continuum electrostatics for a molecule in a spherical cavity of molecular size in the liquid, neglecting short-range interactions and correlations of the molecules. Eq. (4) can be modified by including the factor g_K to account for short-range orientation correlations between molecules, or by using an ellipsoidal cavity to account for molecular shape. Local field factors enter the analysis of many experiments but usually cannot be separately measured.

An alternative to continuum model expressions such as Eq. (4) is to obtain the liquid dielectric properties from a molecular dynamics (MD) simulation including the explicit molecular details for the system of interest.^16–26 Such MD simulations have become feasible for systems ranging from simple polar liquids to protein molecules and larger biomolecular systems.^27–32 Electrostatic interactions often control the behaviour of complex, heterogeneous systems, but may be poorly represented by continuum model results.³³ Molecular dynamics simulations give easy access to the microscopic details, but the simulations usually employ semi-empirical and most often non-polarizable models, have difficulties with long-range correlations, and need validation. In particular, the dielectric constant is a liquid property determined by such MD simulations which often disagrees with experiment.^27,28 Recent MD simulation results find local fields^18,19 much lower than the Onsager local field from Eq. (4), indicating inadequacy of the continuum model results for dielectric properties, and the need for independent experimental determinations.

HRS experiments and analysis, such as that previously presented for water,^1,2 provide direct experimental measurement of the orientation correlations in polar liquids. In the following work, similar measurements and analysis of HRS angle and polarization dependence will be presented for liquid nitrobenzene (NB). Analysis of the HRS experimental data is complicated by the presence of several overlapping contributions to the polarization and spectral functions, which are dissected and described using models with consistent interrelated parameters. The HRS results will be combined with measurements of ion-induced HRS, MD simulation, and the results of dielectric measurements, to determine short-range and long-range orientation correlations, and the static local field factor in liquid nitro benzene.

The HRS measurements in this work were made with linearly polarized light at scattering angles in the range from 0° to 180°. Scattering configurations with incident and scattered light polarized either perpendicular or parallel to the horizontal scattering plane are denoted as VV, HV, VH, and HH, where V denotes vertical polarization, H denotes horizontal polarization, and the first and second letters refer to the incident and scattered light, respectively. The HRS intensities at scattering angle θ_s in the horizontal plane are given by the expressions^1,2

The intensity coefficients A_T and A_L account for HRS from that part of the third rank tensor β which transforms as a vector, while the A₀ intensity coefficient accounts for the remaining HRS contributions, including those induced by short-range multipole intermolecular interactions. In the case that the molecules have no long-range orientation correlations, one has A_T = A_L, and the angular dependence of the combined A_T and A_L contribution reduces to exactly the same form as the A₀ contribution.

The third rank Cartesian molecular hyperpolarizability tensor β_ijk(−2ω; ω, ω) for second-harmonic scattering is symmetric in the last two indices j and k since the corresponding applied optical fields are indistinguishable, and far from resonance β also has approximate permutation symmetry in all indices (Kleinman symmetry, which becomes exact in the limit ω → 0).^5,6 An alternative expression for β is the direct sum of four irreducible spherical tensors,⁵ 

where $βm[ν, l]$ is a spherical tensor of rank l with 2l + 1 components m, index ν labels the symmetry under permutation of the Cartestian tensor indices (ss is totally symmetric, while ms is non-symmetric for first index permutations), and the mixed symmetry ν = ms tensors vanish when Kleinman symmetry holds. The vector (l = 1) β^[1] tensors transform the same as the dipole vector under rotations and are affected by the same orientation correlations as the dipole vector.

For a homogeneous, isotropic random vector field $U→(r→)$⁠, the most general correlation function has the form^1,34

where B_T(0) = B_L(0). The tensor B_ij is diagonal when one coordinate axis is aligned along $r→$⁠, with the transverse and longitudinal correlation functions B_T and B_L as the diagonal components. The corresponding spatial spectrum is

where the tensor S_ij is diagonal when one coordinate axis is aligned along the scattering vector $K→$⁠, with the transverse and longitudinal spectra S_T and S_L as the diagonal components. The scattering wavevector is $K→=2k→i−k→s$⁠, where $k→i$ and $k→s$ are the wavevectors of the incident and scattered light, respectively. Expressing the HRS intensities due to vector β in terms of correlation functions between molecular orientation vector components, one may show that the vector β^[1] HRS intensities have the form given in Eqs. (5)–(8), with intensity coefficients A_T and A_L proportional to S_T(K) and S_L(K), respectively.¹ 

The polarization factor R for vector β^[1] HRS (for molecules with C_2v, C_3v, C_4v, or C_∞v symmetry) is given by

where $β||=(1/5)∑α(βzαα+2βααz)$⁠, $β⊥=(1/5)∑α(2βzαα−βααz)$⁠, $β[ν, l]2=∑mβm[ν, l]2$⁠, and R = 3 when Kleinman symmetry applies.^5,6,35 Orientation correlations with l = 1 do not affect HRS due to the irreducible spherical components of β with l = 2 or 3, and for the A₀ HRS terms due to these β components,

with P² = 3/2 when Kleinman symmetry applies.⁵ The A₀ HRS terms also have contributions from transient molecular hyperpolarizabilities induced by the quadrupole and octopole fields of neighbouring molecules during collisions, for which P² > 3/2 is possible.^36–39 

In addition to the A_T and A_L HRS intensities, modified by dipolar correlations intrinsic to the neat polar liquid, there is a contribution to the A_L HRS intensity coefficient due to radial orientation correlation of molecular dipoles induced by the Coulomb electric field around each dissolved ion.^40,41 The ion-induced HRS intensities are given by

where

is the Debye-Huckel screening parameter, $γ⊥=(1/15)∑αβ(2γαββα−γααββ)$ where γ_αβγδ(−2ω; 0, ω, ω) is the second hyperpolarizability,³⁵ and ρ₊ (ρ₋) is the number density of ions with charge Z₊e (−Z₋e). The factor C includes all other geometric, optical field, and instrumental factors. The scattering wavenumber K is given by

with maximum value

where sin²(θ₀/2) = (n_i − n_s)²/(4n_in_s), λ₀ is the vacuum wavelength for the incident light, and n_i and n_s are the refractive indices for the incident and scattered light, respectively. The ion-induced HRS contribution at low ion concentration comes mainly from molecules at distance r ≈ π/K from the ion⁴¹ and is due to the small orientation correlation induced in many molecules located at a large distance from the ion where the electric field of the ion is slowly varying in space and time. For nitrobenzene at T = 25 °C with λ₀ = 1064 nm, one has n_i = 1.5262, n_s = 1.5589, θ₀ = 1.2°, and K_m = 2π/172 nm.^42–44 For ε_s = 34.82 and Z₊ = Z₋ = Z = 1, the screening parameter $KD2=Km2$ will be obtained for ρ₊ + ρ₋ = 1.09 × 10⁻⁴ mol/L.⁴⁵ 

The maximum ion-induced HRS intensities, obtained from Eqs. (15) and (16) in the limit as θ_s → 0 and K_D → ∞, are

For comparison, the HRS intensity from the same liquid with unperturbed and uncorrelated randomly oriented molecules is

The HRS experimental apparatus and methods are similar to those previously described.^{2,40,41,46–50} Linearly polarized pulses from an injection-seeded single-longitudinal-mode Nd:YAG (yttrium aluminum garnet) laser (operating at λ₀ = 1064 nm, 4.3 kHz repetition rate, 100 ns pulse duration) are focused to a 4.5 μm beam waist radius in the liquid sample of deuterated nitrobenzene (C₆D₅NO₂, Aldrich, 99.5 atom % D, 0.2 μm filtered, deuterated to minimize thermal lens effects due to absorption of the laser beam). Scattered light is collected and collimated by an aspheric lens (f = 13.8 mm), analyzed by a linear polarizer, focused into an optical fiber, and fiber-coupled to a spectral filter or spectrometer followed by the photon counting detector. The sample temperature was T = 25.0 °C, and the laser beam average power in the sample was typically P_av = 1.5 W.

Ratios of HRS intensities were measured (I_{V V}/I_HV, I_HV/I_{V H}, and I_HH/I_{V H}). The effect of intensity drift was cancelled by using several hundred alternate 10 s measurements of the two polarization configurations for each ratio. Rapid switching between polarization configurations was enabled using a liquid crystal variable wave plate (LCVWP) to control the laser polarization and a fast rotator to control the analyzing polarizer for the scattered light. The collection numerical aperture (NA = n sin θ) was controlled by a circular aperture following the collection lens, and the HRS intensity ratio at NA = 0 was obtained by extrapolating measurements in the range of 0.04 < NA < 0.11, made using several apertures, to zero collection aperture.

HRS measurements were made at scattering angles θ_s in the range from 0° to 180°. The sample cell was either a special 8-window cell for 45°, 90°, or 135° measurements or a standard square 10 mm fluorimeter cuvette for θ_s measurements at or near 0°, 90°, or 180°. The laser beam was normally incident on the entrance window in all cases, but off-normal incidence at the exit window of the cuvette was required for some of the measurements. Measurements at θ_s = 0° and 180° were made with the collection lens at θ_s = 0°, with the incident laser beam reflected by a plane dichroic mirror following the collection lens and the HRS light passing through the mirror. The mirror is one focal length from the lens and tilted 1° so that the reflected beam, refocused by the lens, is parallel to the incident beam in the sample but displaced by 0.5 mm. Collection of HRS from either the forward or retro-reflected beam waist (θ_s = 0° or 180°) is selected by a 0.5 mm transverse displacement of the collection optics.

The HRS polarization ratio measurements were made with an interference filter selecting a 60 cm⁻¹ (2 nm full width at half maximum, FWHM) spectral band centered on the 532 nm second-harmonic wavelength. Since the ion-induced HRS signal is difficult to eliminate, it was separately measured and then subtracted. The ion-induced HRS contribution was measured by inserting a fiber-coupled confocal Fabry-Perot interferometer (10 cm mirror spacing, 750 MHz free spectral range) and scanning the spectrum with 13 MHz resolution. The spectral broadening of the ion-induced HRS signal due to the diffusive motion of the ions is typically <1 MHz, as compared to spectral broadening >10 GHz for the other components of the HRS spectrum, so the ion-induced signal appears as a sharp spike riding on the flat background produced by the multiple overlapping orders of the much broader HRS spectrum. The ratio of the integrated intensities of the spike and background (S/B) was determined from the Fabry-Perot spectrum and used to subtract the ion-induced contribution. Typically, 10⁴ alternate scans of the reference spectrum [laser second-harmonic produced by a potassium titanyl phosphate (KTP) crystal] and the HRS spectrum were recorded, where the reference scans were used to lock the center of the Fabry-Perot scans to the laser second-harmonic, and to provide the instrument response function used for de-convolution of the ion-induced HRS spectrum.

HRS spectra were also measured at θ_s = 90° with lower resolution using a grating spectrometer (Jobin-Yvon U1000), with reference and HRS scans alternated in the sequence ref-VV-ref-HV (100 repeats). The reference spectrum was used for de-convolution of the HRS spectra, and the spectral variation of the polarization ratio was obtained from the alternate scans of the VV and HV spectra. The VV and HV HRS spectra are not contaminated by ion-induced HRS.

Table I and Figure 1 show measurements of the HRS polarization ratios I_{V V}/I_HV, I_HV/I_{V H}, and I_HH/I_{V H} for nitrobenzene as a function of scattering angle. The measurements have been corrected for the ion-induced HRS contribution and extrapolated to zero collection angle. Photon counting statistics account for the stated error bars on most of the data points. The ion-induced HRS contribution is <1% for all measurements with θ_s ≥ 90°, but it becomes large at small forward scattering angles, especially for HH polarization. The measured ratio of the ion-induced spike to the spectrally broader intrinsic solvent background HRS (spike/background, S/B) increases from 30% at 45° to 200% at 15.7° for HH polarization. The accuracy of I_HH/I_{V H} measured at small forward angles is limited by the uncertainty of the ion-induced spike correction. The intensity of the ion-induced spike is a strong function of scattering angle and sample purity. Figure 2 shows the angular and ion-concentration dependence of the ion-induced HRS intensity, and the relative size of S/B for HH and VH polarizations.

Fig. 1(a) shows that Eqs. (5)–(8) with A_L = 0 are a good fit to the HRS polarization ratio data. The K dependence of A_T can be more closely examined by setting A_L = 0 in Eqs. (6) and (7), and solving for A_T/A₀, giving

This expression gives an independent determination of A_T(K)/A₀ at the corresponding value of K ≈ K_msin(θ/2) for each I_HV/I_VH measurement (assuming A_L = 0). These results are shown in Fig. 1(b) and are consistent with A_T independent of K. The A_T term as a fraction of I_HV is given by

where F_HV,T = 0.5285 ± 0.0012 using A_T/A₀ from Fig. 1. The molecular dipoles in the liquid are distributed as a homogeneous, isotropic, random vector field. As shown in Ref. 1, a pure transverse spectrum which is finite and flat at K = 0 requires that the orientation correlation function varies as r⁻³ at long-range. Over the observed K range, A_T is flat within ±4%. This requires a correlation function varying as r^−3±δ at long-range with |δ| ≤ 0.02, which is consistent with the long-range orientation correlations produced by dipole-dipole interactions.

The polarization factor R = 3.12 from the fit in Fig. 1(a) is consistent with expected small deviations from Kleinman symmetry for the vector β^[1] HRS contribution, but the value P² = 1.81 from the fit is not possible with l = 2, 3 β HRS contributions alone, with or without deviations from Kleinman symmetry. However, additional contributions due to multipole-induced hyperpolarizabilities for colliding molecules can account for the observed P² value. The HRS spectrum contains more information about these contributions.

Figure 3(a) shows VV HRS spectra for nitrobenzene obtained from 100 cm⁻¹ scans at 3.8 cm⁻¹ resolution and 10 cm⁻¹ scans at 0.6 cm⁻¹ resolution (FWHM). Also shown for comparison in Fig. 3(a) is the VV HRS spectrum for benzene (C₆D₆, Aldrich, 99.6 at. % D). The VV HRS spectra and the corresponding HV spectra (not shown) for nitrobenzene were simultaneously fit with the sum of two Lorentzian components and one exponential component,

where ν is the frequency shift from the laser second-harmonic frequency 2ν_L. Eqs. (25) and (26) are multiplied by the factor [(ν + 2ν_L)/2ν_L]⁷exp(−hcν/2k_BT) to account for the spectrometer spectral slit width and HRS frequency dependence and the Boltzmann factor, and then convolved with the measured instrument spectral response function before comparison with the data. The benzene spectrum was similarly fit including only the exponential contribution. The NB fit parameters are summarized in Table II, and the decomposition into spectral components for the fit is shown in Figs. 3(b) and 3(c).

Several constraints are placed on the fit. The ν₁, ν₂, and ν₃ components are assumed to be due to the β^[1], β^[2] ⊕ β^[3], and collision-induced contributions, respectively, so R² = 9.73 from Fig. 1(a) and $P22=1.5$ are assumed for the ν₁ and ν₂ spectral fit components. The ratio of integrated intensities I_VV/I_HV for the fit (in a 60 cm⁻¹ spectral band with center at ν = 0 cm⁻¹) was constrained to agree with the results from Fig. 1(a), I_VV/I_HV = 5.996 including all three components, and I_VV/I_HV = 1.81 excluding the ν₁ component, and the integrated intensity fraction F_HV,1 for the ν₁ component of I_HV was constrained to agree with F_HV,T = 0.5285 ± 0.0012. The value for I_VV(ν)/I_HV(ν) from the fit was also constrained to agree with the value 2.75 ± 0.05 measured at ν = − 30 cm⁻¹. The HRS spectrum for the centrosymmetric benzene molecule is entirely collision-induced. It is expected that the $P32$ and ν₃ fit parameters for nitrobenzene will be close to the values I_VV/I_HV = 3.29 ± 0.01 (from Ref. 48) and ν₃ = 28 ± 2 cm⁻¹ measured for benzene (although intensity of the ν₃ collision-induced contribution to I_HV is 49 times larger for nitrobenzene). The constrained fit to the nitrobenzene HRS spectral data is adequate and supports the interpretation of the polarization data. The fit parameters differ somewhat in this work from those obtained from the previous measurement for nitrobenzene,⁵¹ but the fraction F_HV,1 of the HV intensity for the ν₁ contribution is almost identical to the previous result (0.528 vs. 0.526).

Ion-induced HRS can provide additional information about orientation correlations. Table III and Fig. 4 show previous ion-induced HRS measurements for solutions of triflic acid (CF₃SO₃H, TfA) in nitrobenzene (C₆D₅NO₂)⁵⁰ re-analyzed with the solvent impurity correction removed and mole fraction of added TfA converted to molar concentration x. The anion concentration z, due to the triflate ions produced by dissociation of the added TfA, added to z_i for the solvent impurity ions initially present, is given by⁵⁰ 

and where ρ₊ + ρ₋ = 2z for a symmetric electrolyte. Fig. 4 shows fits to the HRS data using^40,52

with the fit parameters (S/B)_∞ = 0.384 ± 0.009, z₀ = 28.8 ± 2.0 μM, z_i = 0.4 ± 0.1 μM, K_a = 0.4 ± 0.1 mM, and ν_S,0 = 45 ± 15 kHz. Evaluating Eq. (17) for K_D at θ_s = 90° for nitrobenzene gives (ρ₊ + ρ₋)/2 = z₀ = 27.2 μM as compared to 28.8 ± 2.0 μM from the fit. The fit in Fig. 4 indicates that ionic contaminants at the 0.5 μM level can account for the ion-induced HRS in the present measurements, and that dissociation of TfA is >80% complete for x < 0.1 mM. The ionic dissociation constant K_a for TfA in nitrobenzene has not been measured, but the value from the fit is comparable to the measured values for other strong acids and their salts in nitrobenzene [20 mM for (CH₃)₄NSO₃F, 0.5 mM for (C₆H₅)₃CBF₄, and 6 μM for HSO₃F].^53,54 The parameter ν_S,0 = K²D/2π = 45 kHz corresponds to the ion diffusion coefficient D = 4.3 × 10⁻¹⁰ m² s⁻¹, as compared to 2.7 × 10⁻¹⁰ m² s⁻¹ for (C₆H₅)₃C⁺ in NB⁵³ and 5.7 × 10⁻¹⁰ m² s⁻¹ for SO₃F⁻ in NB,⁵⁴ estimated from conductivity measurements using the Nernst–Einstein equation. The overall agreement between ion-induced HRS theory and experiment is good enough that the ion-induced HRS measurements can be combined with Eqs. (15) and (22) to make a useful determination of the local field factor (see below).

The molecular dipole orientations in a liquid can be described using correlation functions B_L(r) and B_T(r) for a homogeneous, isotropic, random vector field, and the HRS intensities due to vector β^[1] can be expressed in terms of these correlation functions. The observation that the vector β^[1] HRS contribution is purely transverse restricts the allowed form for the correlation function, requiring that the correlation function has vanishing longitudinal spatial spectrum S_L(K) = 0, and a transverse spatial spectrum S_T(K) which is finite and non-zero at K = 0. Such a correlation function, with the required B ∝ r⁻³ asymptotic form due to the long-range dipole interaction, is¹ 

where B_L(0) = B_T(0) = 1 and B_L(∞)/B_T(∞) = − 2, and where the scale length a is the only undetermined parameter. At short range, the position and orientation correlations imposed by specific molecular size and shape become important, so at short range, the orientation correlation function for molecules in the liquid can be expected to differ from the simple generic form of Eqs. (30) and (31). The molecular distribution at short range in the liquid is most directly and completely addressed by molecular dynamics (MD) simulation, and the molecular orientation correlation function can be completely determined by combining the HRS and MD results.

MD simulations for N = 8000 nitrobenzene (C₆H₅NO₂) molecules in the NPT ensemble at T = 298 K and P = 101 kPa were performed using the gromacs software package (version 5.0.6).⁵⁵ The molecules were in an L = 11 nm cubic box with periodic boundary conditions, long-range Coulomb interactions were treated by the particle mesh Ewald (PME) method, while Lennard-Jones interactions were treated using a 1.0 nm cutoff radius. The MD simulations used molecular topology and structure from virtualchemistry.org (CAS 98-95-3 files) for models with OPLS/AA (Optimized Potential for Liquid Simulation for All Atoms) force fields,^27,28 CGenFF (CHARMM General Force Field) and GAFF (Generalized Amber Force Field) force fields. The OPLS/AA model used atom type opls-767 Ar-NO2 for N12, and the CGenFF model used parameters updated for the nitrobenzene NO₂ atoms (charmm36-jun2015.ff). The OPLS/AA partial charges are based on the results of quantum chemical calculation for the molecule. The CGenFF model partial charges are adjusted to make the molecular dipole moment better agree with the experimental gas phase dipole, whereas the GAFF model partial charges overestimate the permanent dipole moment to account for the effect of molecular polarizability. The CGenFF model was also used to simulate nitrobenzene-d5 (C₆D₅NO₂). The time step was 2 fs, and the system was equilibrated for 5 ns before the 10 ns production run. Similar results were obtained for 20 ns MD simulations of 1000 molecules.

The principal results calculated from the MD trajectories are the mean square fluctuation of the total dipole M, $GK=(M2−M2)/(Nμ2)$⁠, the time autocorrelation function $P1(τ)=μˆi(t+τ)⋅μˆi(t)$ and correlation time τ₁ for the molecular dipoles, the radial distribution function g(r), and the dipole-dipole correlation functions $F(r)=μˆi⋅μˆj=cosθij$⁠, $L(r)=(μˆi⋅rˆij)(μˆj⋅rˆij)=BL(r)$ and T(r) = [F(r) − L(r)]/2 = B_T(r), for molecular pairs ij with center-of-mass intermolecular separation $r→ij$⁠. Figure 5 shows the MD results for g(r), L(r), and T(r). The MD results for g(r) and F(r) are similar to the results from previous investigations using several model force fields for NB.^20,25 The pair distribution functions g(r) plotted in Fig. 5(a) show that other molecules are excluded for r < 0.32 nm and that positional correlations vanish for r > 2.0 nm. The orientation correlation functions r³L(r) and r³T(r) plotted in Figs. 5(b)-5(d) are expected to tend to constants at long range where the correlations tend to the form given by Eqs. (30) and (31). The large oscillations seen for r³L(r) flatten out at about r = 3.0 nm, and at larger distances the simulation results become noisy. The smaller oscillations in r³T(r) flatten out at about r = 1.5 nm.

The electrostatic lattice sum for the infinite periodic array of cubic simulation cells in the MD simulation is only conditionally convergent. It is evaluated as the absolutely convergent Ewald sum, plus the remaining conditionally convergent contributions which have the form of a boundary surface term.⁵⁶ For an infinite spherical sample with dielectric constant ε, surrounded by a medium with dielectric constant ε′, the surface contribution to the energy for the simulation cell is 2π M²L⁻³/(2ε′ + 1).^56,57 The present simulations use the (conducting) boundary condition ε′ = ∞ for which this contribution vanishes, but it does not vanish for the boundary condition ε′ = ε describing a sample embedded in an identical dielectric continuum. Corrections are required to obtain the results for ε′ = ε, for the infinite homogeneous system. The required corrections to the correlation functions are determined from the asymptotic (r → ∞) expressions,^57,58

where

No correction is required for $hε′112(r)$ since the asymptotic value is independent of ε′, but $h∞110(∞)=(3yN)−1(ε−1)2/ε$ must be subtracted from $hε′=∞110(r)$ to obtain the correct asymptotic result $hε110(∞)=0$⁠. Corrections $3δF=9δL=9δT=−h∞110(∞)=−1.2±0.2×10−4$ were applied, and it is the corrected results from the MD simulations that are shown in Figure 5.

The dielectric constant and the Kirkwood factor can be calculated from the MD simulation results using the expressions ε₁ = 1 + 3y G_K and g_K,1 = G_K(2ε₁ + 1)/(3ε₁),^57,59,60 but the accuracy is limited by slow fluctuations of $M2$ which do not decrease with increasing system size. Relative uncertainty 5%–10% for G_K, ε₁ and g_K,1 is estimated using $σε−1/ε−1=[2τM/(3t)]1/2$ [Eq. (13) in Ref. 57]. Smaller statistical error may be obtained using the correlation function expression $gK,2=1+(4π/3)ρ ∫hε110(r) r2dr=1+4πρ ∫F(r) g(r) r2dr$ and solving (ε₂ − 1)(2ε₂ + 1)/(3ε₂) = 3yg_K,2  for ε₂.⁵⁷ A second alternative is to solve $(ε3−1)2/ε3=4πyρlimr→∞r3h112(r)=18π yρ a3$ for ε₃ and use g_K,3 = (ε₃ − 1)(2ε₃ + 1)/(9yε₃),⁵⁷ where a³ is the parameter in Eqs. (30) and (31) (determined below). The static dielectric constant ε_s including the effect of molecular polarizability is obtained using Eq. (1) with ρ, μ, and, g_K from the MD simulation combined with ε_∞ = 3.50 from dielectric measurements.⁶¹ The results of these calculations are compared with experimental measurements in Table IV. Good agreement with experiment for density ρ, dielectric constant ε_s, and correlation time τ₁ is obtained for the CGenFF model, especially using correlation function expressions for g_K and ε. This agreement indicates that orientation correlation is determined by short-range forces and the permanent dipole moment, and is insensitive to the molecular polarizability.

The parameter a is determined from the MD simulation by requiring that it produces a pure transverse HRS spectrum. The β^[1] HRS intensities are directly related to the spatial spectra of the correlation functions in the K = 0 limit, given by¹ 

where the integrals are evaluated piecewise over three regions. The delta function self-correlation in region 1 (r < r₁ = 0.32 nm) gives¹ 

where $4π r03/3=ρ−1$ is the volume per molecule, r₀ = 0.344 nm, and S_T,1 = S_L,1 = 56.8 × 10⁻³ nm³ for NB. The MD correlation functions integrated over region 2 (r₁ < r < r₂ = 3.0 nm) give S_T,2 and S_L,2, while the long-range contributions from the correlation functions given by Eqs. (30) and (31), integrated over region 3 (r > r₂ = 3.0 nm), give¹ 

A pure transverse spectrum with S_L(0) = 0, in agreement with the observed HRS spectrum, is obtained by adjusting the parameter a in Eqs. (38) and (39). The intensity contributions to Eqs. (35) and (36) are given in Table V for each MD simulation. The parameter a from Table V is used to evaluate Eqs. (30) and (31) for the long-range asymptotic correlation functions, which are plotted in Figure 5 for comparison with the MD correlation functions. The HRS intensity S_T(0) is larger by the factor C_T = S_T(0)/S_T,1 than the intensity S_T,1 that would be obtained for β^[1] HRS in the absence of orientation correlation. Using C_T from Table V, I_HV and I_VH can be expressed in terms of β_⊥ as

where all orientation correlations, short-range interactions, and β components are taken into account.

The static local field factor f(0) appearing in Eq. (15) for (S/B)_∞ can be determined by combining the results of HRS and dielectric measurements and MD simulations. The ion-induced HRS result (S/B)_∞ = 0.384 ± 0.009 for VH polarization at θ_s = 90° from Fig. 4 can be expressed using Eqs. (15), (40), and (41) as

Orientation correlations enter through g_K for ion-induced HRS³⁸ and through C_T for background HRS. The molecular hyperpolarizability β_⊥ would exactly cancel in this expression if γ_⊥ could be ignored. Using the gas phase measurements μ = 4.22 D,^45,β_|| = 6.3 ± 0.3 × 10⁻⁵¹ C³ m³ J⁻², and γ_|| = 1.2 ± 0.2 × 10⁻⁶⁰ C⁴ m⁴ J⁻³,⁶² assuming γ_⊥/β_⊥ = γ_||/β_||, and ignoring the difference between gas and liquid phase molecular properties, gives Q = 3k_BT γ_⊥/μ β_⊥ = 0.170 ± 0.034. Then using Eq. (42) with ε_s = 34.46 at T = 300 K, (S/B)_∞ = 0.384 ± 0.009, F_HV,T = 0.5285 ± 0.0012, and I_HV/I_VH = 1.359 ± 0.002, one obtains $(1+Q/gK)μ gKf(0)CT−1/2=2.631±0.031×10−29$ C m = 7.888 ± 0.09 D from the HRS results.

Analysis of dielectric response measurements using the Kirkwood-Onsager relation given by Eq. (1) determines the factor μ ²g_K. The high frequency dielectric constant ε_∞ entering the analysis can be estimated by extrapolating the optical refractive index (ε_∞ = n² = 2.3 at 20 °C),⁴³ by extrapolating dielectric response measurements for ν < 14 GHz (ε_∞ = 3.3 at 21.4 °C),⁶³ or by fitting dielectric response measurements for ν < 3 THz (ε_∞a = 3.50 at 20 °C).⁶¹ This last experiment reaches the high frequency limit for all orientation relaxation modes, and ε_∞a includes the directly measured 29 cm⁻¹ librational contribution in addition to the vibrational and electronic polarizability contributions. In that work,⁶¹ the value μ = 4.00 D for NB was obtained from measurements of the dielectric increment for dilute solutions of NB in the non-polar solvents CCl₄ and benzene, and then the value g_K = 0.65 for neat NB was obtained using Eq. (1) assuming that μ for NB is the same in the pure liquid and in dilute solution. The small value for g_K was interpreted as due to formation of anti-parallel NB dimers, consistent with negative B_L and B_T for r < 0.6 nm in the MD simulation. That experiment determined μ ²g_K = 10.4 D² for nitrobenzene at 20 °C.⁶¹ 

One can determine the local field factor f(0) by combining μ ²g_K = 10.4 D² obtained from the dielectric response measurements, with $(1+Q/gK)μ gKf(0)CT−1/2=7.888±0.09$ D and Q = 0.170 ± 0.034 obtained from the HRS analysis. The Onsager local field factor f_O(0) = 2.617 is given by Eq. (4) with ε_s = 34.46 (T = 300 K) and ε_∞ = 3.50. Inserting the values for C_T and g_K obtained from the MD simulation into these equations, and solving for μ and f(0)/f_O(0), gives with the results shown in Table VI. The stated error bars in Table VI account for statistical uncertainties for the experimental and MD simulation results. Solutions for all the MD simulations give f(0)/f_O(0) < 1, but only the analysis using input from the CGenFF MD simulation gives a nearly self-consistent value for μ. Taking the CGenFF results as the best estimates, one sees that the NB molecular dipole decreases from μ = 4.22 D in gas phase⁴⁵ to 4.00 D for dilute solution in non-polar solvent,⁶¹ to 3.90 D in pure liquid NB. These results indicate that μ is reduced by short range intermolecular interactions in the liquid, although previous second-harmonic generation experiments that have probed the effects of intermolecular interactions for NB do not give evidence for μ decreasing in a more polar environment.^64,65

Recent MD simulations for nonlinear optical chromophores dissolved in chloroform have found local field smaller than the Onsager value from Eq. (4), by a factor of 0.74 for disperse red (DR),¹⁹ 0.68 for p-nitroaniline (pNA), and 0.64 for p-nitro-N,N-dimethylaniline (pNDMA).¹⁸ These simulations, for 2 to 16 chromophore molecules in a box with 512 chloroform molecules, determine local field from the chromophore orientation order induced by an applied electric field. The results can be fit by a generalized Onsager local field function f = ε(L_aε_∞ + 1 − L_a)/[ε(1 − L_a) + L_aε_∞], where L_a is the depolarization factor for the molecular cavity, and where L_a is 1/3 for a sphere and goes to zero for a very narrow prolate ellipsoid.^14,18,19 The L_a and f that fit the MD results have a weak dependence on the strength of the applied field. The local field values from these MD simulations are much smaller than 0.85 × f_O(0) for NB from the present work. This could be the combined effect of f < f_O(0) and orientation correlation g_K < 1 for the probe molecules in solution for the MD simulations. These results indicate that local fields are smaller than the Onsager value, but more work is needed to disentangle the possible effects.

HRS can be partitioned into several contributions with different polarization, angle, and spectral dependence, which are sensitive to different molecular interactions, correlations, and dynamics. The vector β^[1] contribution is sensitive to correlations due to both long-range dipole-dipole coupling and short-range intermolecular interactions, and HRS is able to distinguish between transverse and longitudinal vector correlations. Orientation correlation functions over the entire range of r are determined by combining the results of HRS measurements and MD simulations for NB. The HRS intensity contributions due to short-range and long-range correlations nearly cancel in liquid NB, and the intensity enhancement for vector β^[1] HRS due to orientation correlations is only 5%. The assignment of the different HRS contributions made in the polarization analysis is validated when the HRS spectrum is decomposed into corresponding components. The narrow central component has the polarization of vector β^[1] HRS and has width consistent with the P₁ dipole autocorrelation time from the MD simulation. The intermediate-width spectral component has the polarization of octupolar β^[3] HRS, and the broad spectral wing has polarization and width consistent with collision-induced HRS. The long-range radial orientation correlations induced by the Coulomb field around dissolved ions produces another HRS contribution which is identified by its characteristic sharp spectrum and longitudinal polarization. The signal due to this long-range ion-induced orientation correlation has been employed to determine the molecular dipole μ, Kirkwood factor g_K, and local field factor in nitrobenzene.

This work was supported by the National Science Foundation (NSF) through Grant No. CHE-1212114.