Calculating the Raman Signal Beyond Perturbation Theory for a Diatomic Molecular Crystal

We calculate the eigenstates of a diatomic molecule in a range of model mean-field potentials, and evaluate the evolution of their associated Raman spectra with field strength. We demonstrate that dramatic changes in the appearance of the Raman spectrum for a diatomic molecule occur without any associated change in the symmetry of the surrounding potential. The limiting cases of the quantum eigenstates correspond, in the classical sense, to free rotation, and libration of well-oriented molecules. However, there are also many mixed modes which are neither rotons nor librons. The consequence for Raman spectroscopy is a series of complications - the non-harmonic potential splits the Raman active modes, and breaks the selection rules on forbidden modes. The mass-dependence of the various states is different - rotors, oscillators and reorientations have $1/m$, $1/\sqrt{m}$ and weaker mass dependence respectively. This may allow one to identify the character of the mode with isotope spectroscopy. However it is complicated by mixed modes and transitions between two different eigenstates with different character. We conclude that significant changes in the Raman spectrum of molecular systems are insufficient to demonstrate a phase transition since such changes can also occur in a fixed symmetry potential upon increasing field strength.

Pure rotons: R=2 Hindered rotor re-orientations: R=1 Librons: R=√2 Oscillator re-orientations: R≈1 • Selection rules are valid only in high and low potential limits: for an intermediate strength potential many more modes are spectroscopically allowed.

Introduction
Optical spectroscopy is one of the most important methods with which one can study solid materials. The two main methods are Infrared (IR), in which light is absorbed by the sample, and Raman, in which light is reemitted at a different frequency. The energy of the absorbed IR photon, or the shift in energy of the Raman signal, correspond to quantised excitations in the material.
In a free molecule, these quantized excitations are vibrations or rotations. Selection rules can be derived using group theory. Each molecular mode contributes equally, so the signal is proportional to the number of molecules, N . In a harmonic, crystalline solid, the excitations are the phonons, and we can apply Bloch's theorem to calculate the phonon spectrum, assigning each normal mode to a reciprocal space q-vector. Conservation of momentum means that the only allowed modes -those which produce a spectroscopic signal -have q = 0. The number of modes is less than three times the number of atoms in the unit cell, but because of constructive interference throughout the crystal each mode is amplified by the total number of molecules, so again these modes produce a signal proportional to the number of molecules, N .
Molecular crystals lie between these two extremes. If the molecules are freely rotating, Bloch's theorem applies only to the molecular centre of mass. So there will be a combination of lattice phonons and free-molecule vibrations and rotations. If the intermolecular forces are strong enough to inhibit the rotation, but not strong enough to create long-ranged orientational order, then the assumption of free rotors breaks down. Coupling and reduced symmetry means that all 3N modes are Raman active.
For strong interactions the molecular orientation is fixed and the modes are harmonic oscillations about that direction, known as librons.
This mixed character regime between a free rotor and harmonic oscillator is often described as "anharmonicity", which is trivially correct, in the sense that the modes are not harmonic. Unfortunately, "anharmonic" is often taken to mean that modes can simply be described with perturbation theory on a harmonic oscillator, and their selection rules remain largely unchanged. For strong interactions third-order anharmonic corrections to the librons may give a reasonable approximation, but as the amplitude of oscillation increases this perturbative approach will fail.
Even if the perturbation approach breaks down, we can retain the approximation that the modes are localised on a single molecule. The wavefunctions remain localized even in the presence of strong interactions, provided the crystal has sufficient disorder [1,2]. This disorder can be in either molecular orientation, or in mass [3,4]. Here we consider the calculation of Raman spectra in the limit of sufficiently strong disorder where the nuclear wavefunctions are localized on a single molecule with the surrounding molecules described by a mean-field potential.
We will demonstrate that for a local potential with more than one inequivalent minimum, the eigenmodes include states which look like rotors and/or oscillators about different energy minima. The spectroscopic signal then contains rotor and oscillator modes, but also detrapping (oscillator → rotor) excitations, and reorientation modes (oscillator → oscillator) where the molecule moves between different minima pointing in different directions. Rotor and oscillator modes have precise mass-dependencies where the transitions between energy levels scale with 1/m and 1/ √ m, respectively. Reorientation modes are slightly more complicated because the transition energy comprises of two components: one related to the potential energy difference between the two minima, which is mass independent, and the other related to the zero point energy difference, which scales with 1/ √ m. Therefore, the exact mass relationship of the reorientation modes depends on both the relative depth and curvature of the two potential minima: when the depths are very different there is no mass scaling.
For homo-atomic molecules mode parity leads to strong selection rules, which must be considered for the correct interpretation of the spectra.
We consider two diatomic rotors of mass 1 and mass 2 which we treat as proxies for hydrogen and deuterium respectively. Hydrogen forms a variety of solid phases, several of which are diatomic. The lowest pressure phase, phase I, consists of free rotors on a hexagonal lattice [5]. In phases II and III at higher pressure, the molecules exhibit preferred orientations but their centres of mass retain the hexagonal symmetry seen in phase I [6] [7] [8]. The precise symmetry of these ordered phases is a topic of intense interest and many candidate structures have been proposed. [9] [10] [11] [12] [13].
The parity of the rotor wavefunction also plays a significant role in the solid phases of hydrogen [14]. The overall symmetry of the wavefunction which is given by a combination of the nuclear and orientational wavefunction of the molecule, must remain antisymmetric (ferimions) or symmetric (bosons). Therefore for para-hydrogen, which has nuclear spin groundstate |↑↓ − ↓↑ , the orientational wavefunction must have even parity. In the free rotor case this corresponds to states with J = 0, 2, 4.... Similarly for orthohydrogen, where the nuclear spins are aligned, the orientational wavefunction must have odd parity, corresponding to states with J = 1, 3, 5... for the free rotor.
Several of the results we present have ramifications for interpretation of the Raman spectra of phases of hydrogen. However the potential operators applied here are far simpler than those present in hydrogen and deuterium solids and we therefore refrain from direct comparison with experiment.

Rotor, Oscillator and Reorientational Modes
The Raman active transitions of diatomic molecules are generally modelled as two distinct types of excitation; rotational excitations or 'rotons' and vibrational excitations, either 'librons' or 'vibrons'. Rotational excitations are well described for a diatomic rotor; the roton energy levels are given in 2D by: and in 3D by: where r is the bond length, m is the atomic mass and J is an integer quantum number. The selection rules for Raman transitions for the free rotor are ∆J = 0, ±2. Vibrons and librons are described as quantum harmonic oscillators for which energy levels are given by: with ω the frequency and k the effective spring constant. The different energy-mass relationships for the free rotor (1/m) and harmonic oscillator (1/ √ m) provide a way to characterise Raman modes using different isotopes of a diatomic molecule. [15] [16] Here, we will compare the frequencies of two specific molecules, with atomic mass 1 ("hydrogen") and 2 ("deuterium"): we expect a rotor to give an isotopic energy ratio R = 2 and an oscillator to give R = √ 2, where E i and E j are the calculated energy levels and υ represents the experimentally-measurable Raman shift. For a reorientional mode we expect a weak mass dependence. These modes give rise to a transition from an oscillator state in one minimum in e.g. the y-direction, to a new oscillator state in an inequivalent minimum in e.g. the z-direction. We therefore expect the Raman shift for this mode to be given by: and the value of R to be given by: where ∆E min is the difference in energy between the two minima. If ∆E min is large with respect to k y − √ k z we expect a value of R close to 1. If ∆E min is small with respect to k y − √ k z we expect R ≈ √ 2. Whilst the free rotor or oscillator models can describe the majority of excitations in diatomic systems, there are certain solids for which this approach is insufficient. In particular solid phases of hydrogen where nuclear quantum effects are significant. Here we describe a model for a system that sits between the idealised cases of rotors and oscillators: the hindered rotor.

Modelling a hindered rotor
We consider only the angular dependent modes: rotons and librons. Starting with the angular Schroedinger equation: where r = |r| is the molecular bond length and m is the mass of the nucleus.
We introduce a simple mean-field potential to describe the angular dependence due to local crystal structure. This is conveniently expanded in spherical harmonics.
giving the Hamiltonian in the basis of the free rotor as: where the primes distinguish spherical harmonics used for the potential and wavefunctions.
We then diagonalise the Hamiltonian to find the energy levels n: with hindered-rotor eigenstates now given by: These states are associated with the appropriate nuclear spin isomer to form overall antisymmetric (fermionic) or symmetric (bosonic) states. The nuclear wavefunction does not affect the eigenstate energy, or the polarizability tensor, but it does affect the degeneracy and therefore the intensity of the Raman signal. For H 2 there are singlet and triplet nuclear states, D 2 has singlet, triplet and quintuplet.

Calculating the Raman signal
To evaluate the transitions between states we express the polarizability operator in the hindered rotor basis, this is done through a unitary transformation on the polarizability tensor of the free molecule. We then calculate the Raman signal using the response of the quantum system to a sudden excitation. We compute the time-dependent response to excitation using the Liouville-von Neumann (LvN) equation, and then evaluate the expectation value of the resulting polarization [15,17,18,19]. The Raman spectrum is then obtained by Fourier transform (FT) of the response in the time domain. Full details are given in our previous work [15] and references within. This approach automatically satisfies selection rules based on transition polarizability and nuclear exchange parity.

Results and Discussion
We present results for four different potential energy surfaces (Fig 1), gradually building up in complexity. In each case, the weak field limit is a perturbed rotor, while our model also explores the strong field regime.
The simplest, case 1, is a potential with a single minimum along zdirection with equal curvature in both x-and y-directions; in the limit of strong field we expect to find a single oscillator mode degenerate along x and y. In case 2 we introduce a potential with a minimum along y where the curvature is different in spatial dimensions x and z: here we expect to see an asymmetric 2D oscillator with modes with different frequencies in the two dimensions. Case 3 represents a combination of the previous two cases, there are now 2 inequivalent minima, parallel to the y and z axes. In case 4 we change the relative depth of these inequivalent minima.
All of these model potentials have higher symmetry than the hcp structure of solid hydrogen phase I, so they will have simpler spectroscopic signatures. : Potential isosurfaces for the four cases considered. Distance of the surface from the origin represents the strength of the potential, green is attractive, red repulsive. Case 1 has a potential minimum for the molecule pointing in the z-direction. Case 2 favours the y-orientation, case 3 has two minima with different depths in y-and z-directions, case 4 has two minima with similar depth.
3.1. Case 1: V = P Y 20 We start with a simple potential V = P Y 20 that only has variation in θ and is symmetric in φ. In this potential, the lowest energy state has the molecule pointing along the z-direction. To first order in perturbation theory, the potential is quadratic and symmetric. The scale of the potential, P, is arbitrary, but it is understood that it will vary monotonically with pressure.
The energy levels as a function of potential strength P are shown in Fig. 3 with the corresponding wavefunctions in Fig. 4. At zero potential, we recover the free rotor states, which can be well described by quantum numbers |J, m J and have degeneracies 1, 3, 5 ... (2J+1) with increasing energy. Figure 2 shows that in intermediate field, J is no longer a good quantum number and all states acquire mixed character. However in case 1, mixing in m J is restricted to basis functions with equal |m J | due to the rotational symmetry of the potential around the z-axis. All states can be assigned a clear parity and consist of entirely odd or entirely even J basis functions. States of opposite parity pair up with asymptotically degenerate energies towards the strong field (pure oscillator regime) as shown in Fig. 3. In the limit of high potential the molecule is oriented along z (θ = 0 in polar coordinates), and its excitations can be described using eigenstates of a 2D harmonic oscillator |n 1 , n 2 . The ground state is |0, 0 -a doublet between even and odd parity. The first excited states are degenerate |1, 0 , |0, 1 , a quadruplet when accounting for parity degeneracy. The second group of six excited states are |2, 0 , |0, 2 , and |1, 1 . For a symmetric harmonic oscillator these are still degenerate, but because the potential is anharmonic the |1, 1 state has slightly lower energy. The third group of eight states: |3, 0 , |0, 3 , |2, 1 and |1, 2 , comprises two quadruplets. States with the same absolute value of m J remain degenerate because the potential is rotationally symmetric around z. As a result, each original state with 2J + 1 degeneracy appears to split into J + 1 separate states. Top: isotope energy ratios R for case 1 for first ten excited states against the potential strength P (see eq. 4 for definition). As P is increased the ratio R transitions from 2 for a free rotor to √ 2 (dashed line) for a harmonic oscillator. Spectroscopy cannot measure the eigenstate energies, rather it measures energy differences between eigenstates: the Raman shifts relate to selected (E i − E j ) differences. As noted above, the energies depend on isotope mass, and in Fig. 3 we show the isotope energy ratio for transitions E 0 → E i for the first ten excited states i, in otherwise-identical molecules with atomic masses 1 and 2, respectively. Values of R move smoothly from 2 → √ 2, for all transitions shown, some of which are Raman active. Because the potential is softer than harmonic at large amplitudes, R overshoots and approaches √ 2 from below. Frequency (cm -1 ) Figure 5: Calculated Raman spectra vs P for deuterium in case 1. The right hand yaxis indicates the value of P on a logarithmic scale. We assume a temperature of 300 K.
At low values of P the spectra resemble the free rotor regime with a series of ∆J = 2 roton transitions visible. At high values of P a single peak is visible corresponding to an oscillator transition (|0, 0 → |0, 1 shown in Fig. 4). Figure 5 shows the calculated Raman spectrum for the simplest case 1. At low potential we see the familiar form of the quantum rotors, and at high potential the single peak of a harmonic oscillator (libron). The intermediate potential region is characterised by a very complicated spectrum and the appearance of many new peaks, despite the symmetry of the potential remaining unchanged. These arise from both splitting of the rotor peaks, and from transitions which acquire Raman activity due to the symmetry-breaking potential lifting the rotor selection rules.

Case 2:
Next we consider a potential of the form V = P (Y 22 + Y 2−2 ) that has a single minimum pointing in the y-direction and with different curvature in θ and φ. This has a ground state with the molecule pointing along the potential minimum (Fig. 7). The reduced symmetry of this potential lifts all m J degeneracies in the rotor states.
As with case 1, in the strong potential limit there are degenerate odd and even parity states which can be labelled with 2D SHO quantum numbers |n 1 , n 2 , in energy order: |0, 0 , |1, 0 , |0, 1 , |2, 0 , |1, 1 , |0, 2 , |3, 0 and |2, 1 (Fig. 6, 7). In this limit, the energy levels all become even-odd parity pairs, and energies approximate to (n 1 + 1 2 ) ω 1 + (n 2 + 1 2 ) ω 2 .  Figure 6: Bottom: energy levels for the first twenty-five states of deuterium for case 2 (V = P (Y 22 + Y 2−2 )) relative to the ground state. The degeneracy in |m J | seen in case 1 is now broken by the reduced symmetry of the case 2 potential. This lower symmetry also lowers the degeneracy of the oscillator states in the strong field limit, from four-fold (seen in case 1) to two-fold. Top: isotope energy ratios (R) for first ten excited states in case 2 against the potential strength P. As with case 1, R transitions from 2 for a free rotor to √ 2 (dashed line) for a harmonic oscillator. Frequency (cm -1 ) Figure 8: Raman spectra for case 2 with deuterium. Right-hand y-axis indicates the offset of the spectra by the value of P on a log scale. As with case 1 (Fig. 5), the spectrum undergoes a dramatic transformation from rotor states to libration states upon increasing pressure. Unlike in case 1, the 2D SHO degeneracy is lifted by the different curvatures in θ and φ, resulting in two distinct libron modes at high potential strength.

Case
The third case is a potential with two distinct minima, along the z-and y-directions. In this case the classical energy minimum has the molecule pointing along the y-axis, with a metastable minimum for molecules pointing along z. Although this is still simpler than any real system, we can see in Fig. 9 that for intermediate strength potential the energy levels have no clear structure. The high energy limit still tends to be harmonic oscillator-like, and is sensitive to the relative stability of the two minima. In this case the global minimum is deep enough that the lowest three doublet levels are all related to it and can be labelled |n 1 , n 2 y : |0, 0 y , |1, 0 y and |0, 1 y . Fig. 10 shows that the metastable state, pointing in the z-direction, has an energy of 1800 cm −1 higher than the ground state. This state can be described with two further quantum numbers |n 1 , n 2 z shown as |0, 0 z in Fig. 10. In the energy plot (Fig. 9) we observe that the energy of this state relative to the global ground state rises approximately linearly with the potential, starting from J = 1 and J = 2 states in the free rotor regime. Fig. 12 shows that the |00 y → |00 z transition is Raman active and corresponds to the re-orientation of the molecule from one minimum to another. For this model potential, the re-orientation can be observed at low temperature and intermediate field strength. The isotopic mass ratio for this transition is close to 1 because the transition energy is dominated by the potential energy difference between the two minima.   Figure 9: Bottom: Energy levels for first twenty-five states of deuterium for case 3 (V = P (Y 22 + Y 2−2 − Y 40 )) relative to the ground state. A metastable state exists corresponding to orientation along the z-axis (|00 z in Fig. 10). The system now consists of two sets of independent oscillator states, one for each minimum. Top: isotope energy ratios (R) for first ten excited states for case 3 against the potential strength P. The value of R now transitions from 2 for a free rotor to either 1 or √ 2 for a re-orientation or harmonic oscillator, respectively. The dashed line indicates the re-orientational transition corresponding to the |00 y → |00 z transition shown in Fig. 12. Frequency (cm -1 ) Figure 11: Raman spectra for case 3 with deuterium at 300 K. Right-hand y-axis indicates the offset of the spectra by the value of P on a log scale. Despite the system being composed of two independent oscillators, the Raman spectra broadly resemble case 2, with a single oscillator with two libron modes at high field. This is because the metastable |00 z state is not populated at 300 K. Frequency (cm -1 ) Figure 12: Raman Spectrum at P= 800 cm −1 with reorientaional mode labelled. The temperature in this case is set to 10 K to allow clear identification of the re-orientational mode.

Case 4.
The final case also involves a potential with minima along y and z, but here the two orientations are close in energy. There is an unusual situation where, as the potential increases, the ground state switches its character from pointing along y to pointing along z. To illustrate the crossover we plot the eigenstate energies in Fig. 13, relative to the reference state |00 y rather than the difference from the ground state as in the previous figures. The crossover can be understood as a trade-off between the zero point energies and the potential energy of the two minima pointing in the y-and z-directions. The state pointing along y has a higher potential energy minimum, but it is stabilised by a lower zero point energy. As the potential energy difference between the two minima increases with field strength, while the zero point energy remains unchanged, the state pointing in the z-direction eventually becomes more stable and hence a crossover takes place at around P = 1000 cm −1 . Frequency (cm -1 ) Figure 15: Raman spectra for case 4 with D 2 at 300 K. Right-hand y-axis indicates the value of P on a log scale. The two independent oscillators are closer in energy in this case and consequently the spectra at high field is composed of contributions from each system. Indicated with triangles are two kinds of re-orientational modes: rotational axis re-orientation (orange) and directional re-orientation (green). Frequency (cm -1 ) Figure 16: Raman Spectrum at P= 800 cm −1 for D 2 in case 4 at 300 K. The reorientaional mode between the ground state and a meta-stable state pointing along the z axis ((|00 y → |00 z ) is indicated by the arrow. The mode appears at low Raman shift in this case as the two orientations are close in energy.
As in the previous cases we see a breaking of the 2J + 1 rotor degeneracy as P increases. This also means that as the m J degeneracy is lifted, transitions which in the zero potential limit are associated with the elastic Rayleigh scattering (∆J = 0), acquire a small non-zero frequency and are strongly Raman active. These transitions gradually emerge from the Rayleigh line as the strength of the potential is increased from P = 0 → 100 (see orange triangles in Fig. 15). They represent a more subtle re-orientation associated with the axis of rotation of a weakly hindered rotor with fixed angular momentum. At high field strengths a distinct type of re-orientational modes appear where the molecule transitions from one oscillator minimum to another as shown Fig. 16. In Fig. 15, indicated with green triangles, one of these modes can be seen emerging as a shoulder from the S 0 (0) (J = 0 → 2) peak and gradually moving to lower frequency with increasing P. These modes are Raman active at intermediate field strengths but the predicted Raman intensity gradually decreases as P is increased.

Conclusion
We have presented a detailed analysis of Raman spectra for a diatomic molecule in four different simple potentials composed of spherical harmonics. For weak potentials the eigenstates describe a nearly-free rotor, and the good quantum numbers are J and m J . With strong potential they describe molecules in well defined orientations, with harmonic oscillations (librons) where two SHO quantum numbers are associated with each inequivalent direction. The eigenstates in each direction are described by precisely two quantum numbers since each minimum has two degrees of freedom. For intermediate potential strength, the eigenstates are mixed, with no good quantum numbers. In addition, for homoatomic molecules, all states are symmetric to nuclear exchange, which creates a parity degeneracy at high field.
We can identify three types of Raman active modes: excitations of rotational states, excitations of oscillator states, and re-orientational excitations. Oscillator and rotor excitations have well defined selection rules, but all transitions between mixed states are Raman active. Consequently, as the potential increases, the simple roton spectra at low potential changes to a complex manifold containing many overlapping states. The roton, libron and reorientation modes all have different mass dependencies. Experimentally this can serve as a tool to identify the character of each mode by taking ratios of Raman peak frequencies between different isotopes.
Nuclear exchange symmetry is not explicit in the model, but appears implicitly through the well defined parity of the wavefunctions. It is interesting to note that in the harmonic oscillator limit the parity does not appear as a Raman selection rule because each SHO level is a doublet with odd and even parity. This suggests that at large field strength, changes in the Raman spectra no longer reflect the spin symmetry of the hydrogen molecule, therefore these spectra are agnostic to the breakdown of ortho-and para-states observed in NMR studies at intermediate pressures [20].
We note that even our most complex potential, case 4, has higher symmetry than known diatomic crystalline materials. For example, molecular hydrogen adopts an hcp phase which would be modeled by a more complicated mean-field potential that has one minimum in the c-direction, and six minima with molecules tilted slightly out of the basal plane. The Raman spectra shown here are simpler than could ever be obtained experimentally, however this allows us to easily understand and interpret the character of the different transitions involved 1 .

Data availability
Codes can be obtained from the authors on request. No experimental data was collected in this study.