Simulating electric field interactions with polar molecules using spectroscopic databases

Ro-vibrational Stark-associated phenomena of small polyatomic molecules are modelled using extensive spectroscopic data generated as part of the ExoMol project. The external field Hamiltonian is built from the computed ro-vibrational line list of the molecule in question. The Hamiltonian we propose is general and suitable for any polar molecule in the presence of an electric field. By exploiting precomputed data, the often prohibitively expensive computations associated with high accuracy simulations of molecule-field interactions are avoided. Applications to strong terahertz field-induced ro-vibrational dynamics of PH3 and NH3, and spontaneous emission data for optoelectrical Sisyphus cooling of H2CO and CH3Cl are discussed.

line lists for PH 3 42 and NH 3 43 (available from the ExoMol database; the much larger hot line lists 43,44 are not necessary for the present applications). The ro-vibrational basis was truncated at the maximal rotational quantum number J = 20 and the maximal vibrational quanta v i corresponding to the ground and first excited vibrational states, i.e. v i = 0, 1. For NH 3 , due to strong anharmonicity of the large amplitude umbrella-type vibrational motion we extended the vibrational basis set to include v umb = 0, … , 7.
For both molecules only one quantum wavepacket simulation corresponding to the initial ground ro-vibrational state was performed. For PH 3 this was the v = 0; J = 0, k = 0 state, where the quantum number k is the projection of the total angular momentum onto the molecule-fixed z axis. For NH 3 , due to nuclear spin statistics, the initial state was the v = 0 + , J = 0, k = 0 state. In order to investigate the influence of temperature T on alignment and orientation, simulations were carried out at T = 0 K and T = 10 K for PH 3 . At T = 10 K a total of 81 initial states are populated. Final results were therefore averaged over all quantum trajectories from the different initial states.
The single-pulse THz electric field was represented by the analytical function 39 π τ πν τ τ = = − = E t E t t t E t ( ) cos ( / ) sin (2 ), [ /2, /2], ( ) 0 otherwise, where E 0 is the amplitude of the pulse, τ is the pulse duration, and ν is the central frequency. In all calculations we fixed τ to 4 ps. For simulations of the mixed Raman-THz field-free orientation of PH 3 , the lack of a polarizability model in the ExoMol database meant the simulations of the initial Raman excitation were carried out using a more general variational approach implemented in the TROVE 6 and RichMol program packages. The resulting wavepacket was then used as the initial wavefunction in subsequent THz-driven dynamics simulations, modelled using the ExoMol line list data. Note that all fields are polarized along the Z axis in calculations.
PH 3 . The results of our numerical simulations of PH 3 at T = 0 K and T = 10 K are shown in Figs 1 and 2, respectively. Four time-delayed THz pulses with central frequency ν = . 0 5 ps −1 and peak field strengths at E 0 = 250 kV/ cm and E 0 = 1 MV/cm have been used to induce transient alignment 〈 cos 2 θ〉 (t) and orientation 〈 cosθ〉 (t), where θ is the Euler angle. Such intense THz pulses in the few-cycle regime can be routinely generated in a convenient tabletop setup [45][46][47] . As expected, increasing the number of pulses significantly enhances the degree of orientation and alignment for both temperatures and peak field strengths. The periodic behaviour of the orientation dynamics follows the quantum rotational revival pattern; strong peaks are observed at times t ≈ 2.0, 4.2, 5.7, 7.9, 9.5 ps etc., which is approximately 1/2, 1, 3/2, 2, 5/2, … of a revival time T rev = 1/2 Bc ≈ 3.7 ps (B ≈ 4.5 cm −1 is the rotational constant of PH 3 and c is the speed of light). The enhancement of the alignment and orientation at delay times τ ≈ 4, 8, 12 ps etc., reveals recurrences separated by approximately Δ τ = T rev .
Increasing the peak field strength of the THz pulses results in a more complex pattern of alignment and orientation with additional local minima and maxima. These secondary peaks are always present but attenuated in the weak field regime. Alignment and orientation are enhanced in stronger electric fields, hence higher absolute peak values of 〈 cos 2 θ〉 (t) and 〈 cosθ〉 (t). However, raising the temperature from T = 0 K to T = 10 K activates several new quantum trajectories via thermally populated states and this leads to overall decoherence of the wavepacket. This effect is responsible for the lower absolute values of the orientation and alignment parameters at T = 10 K Scientific RepoRts | 7:45068 | DOI: 10.1038/srep45068 in Fig. 2. Therefore, losses in orientation or alignment at higher temperatures can be compensated for with a stronger peak intensity of the THz pulse.
As an alternative to multiple THz pulses, in ref. 38 it was demonstrated that the degree of alignment and orientation can be improved by applying a short intense Raman pulse to coherently prepare molecules in highly excited rotational states prior to THz exposure. This allows a larger molecular population to occupy states for which transition frequencies match the peak of the THz pulse envelope. Raman excitation follows Δ J = ± 2 selection rules and hence couples states of the same parity. Therefore, when starting from an initial J = 0 state, the Raman pulse predominantly excites J = 2 states. The subsequent THz pulse induces J → J′ : 0 → 1, 2 → 3 transitions and produces a mixed-parity wavepacket, thus enabling effective net orientation. The results of our simulations for PH 3 are shown in Fig. 3. The alignment pattern is insensitive to the time delay between the Raman and THz pulse and it appears the THz pulse contributes only marginally to alignment. For a given time delay the orientation peaks are separated by ≈ 3.9 ps, which corresponds to T rev . NH 3 . As a second example we explore strategies to establish control over rotational coherences and populations in NH 3 by applying optimally timed THz pulses. Recently it was shown both theoretically and experimentally 41 that a pair of THz pulses separated in time can induce larger population transfer (larger transient emission responses) than two interactions within a single short THz pulse.
To adjust the THz pulse parameters to the energy level structure of NH 3 , we first look at the Fourier intensity of a pair of THz pulses at a few characteristic frequencies as a function of the time delay τ between the pulses and the pulse central frequency ν (see Eq. (1)). The results are shown in Fig. 4 for the three rotational transitions, |0, 1〉 → |1, 0〉 , |1, 0〉 → |2, 1〉 , and |2, 1〉 → |3, 0〉 , where states of NH 3 are labelled as |J, v umb 〉 , and the quantum numbers k and m, which correspond to the projection of the total angular momentum onto the molecule-fixed and laboratory-fixed axes, respectively, are zero. It is evident that the optimal values of ν are 0.4, 1.1, and 1.6 ps −1 for the three different rotational excitations of NH 3 . Figure 5 depicts the populations of the ro-vibrational wave packets of NH 3 , initially in the |0, 1〉 , |1, 0〉 , and |2, 1〉 ro-vibrational states, following the excitation by a time-delayed THz pulse pair with E 0 = 500 kV/cm and respective pulse central frequencies ν = . 0 4, 1.1, and 1.6 ps −1 . Population transfer is significantly enhanced by two properly delayed THz pulses when compared to a single, but twice as intense, THz pulse (τ = 0). Slight adjustment of the central frequency of the THz pulses increases the selectivity of excitation, suppressing (or enhancing) two-photon transitions. As expected, population modulations for the three excited rotational states with J = 1, 2 Based on the results of Fig. 5 we can combine pairs of THz pulses to create coherences between different pairs of rotational states with even and odd parities. This is shown in Fig. 6. The first pulse at t ≈ 1.1 ps induces a J → J′ : 0 → 1 transition, creating a coherent superposition | 〉 ± | 〉 ( 0, 1 1, 0 ) 1 2 . The second pulse is tuned to affect population of the |1, 0〉 state and transforms the two-state coherence into a new one: | 〉 ± | 〉 ( 0, 1 2, 1 ) 1 2 . Parameters of the subsequent pulses (upper inset in Fig. 6) are optimized to maintain uniform division of molecular population between the two states. This illustrates the high level of control which can be exerted over the ro-vibrational wavepacket with the use of THz pulses.
To further investigate the validity of our approach we have computed the alignment and orientation of NH 3 using three different Hamiltonian models. The first, given in Eq. (2), is the electric field Hamiltonian in the dipole approximation. The second Hamiltonian we use incorporates the polarizability α into the model, while the third Hamiltonian also includes the hyperpolarizability β. Latter calculations were perfomed using our recently developed variational approach RichMol. The THz field was composed of three simultaneously applied pulses with respective central frequencies ν = .
0 3, 0.6, and 0.9 ps −1 each with duration τ = 4 ps. The results are illustrated in Fig. 7 for electric field strengths up to 10 MV/cm. Calculations were performed on an initial ground vibrational state with J = 1, k = 1, m = 0 and mixed parity + /− ; that is ψ ( 1, 1, 0, 0 1, 1, 0, 0 ) 1 2 . As shown in Fig. 7 it is only for the very intense 10 MV/cm THz pulses that we see contributions from the polarizability and hyperpolarizability, otherwise their effects are negligible. In general, alignment and orientation caused by the hyperpolarizability is an order of magnitude smaller than the polarizability and contributes with opposite sign. Field strengths as strong as 10 MV/cm, where the non-linear THz effects of α and β, although noticeable, are small enough to be neglected, can therefore be modelled using the Hamiltonian in the dipole approximation in conjunction with ExoMol line list data.
Intense THz pulses offer a promising approach for controlling large-amplitude motions in floppy molecules. In Fig. 8 we investigate the effect of an intense THz-pulse-train on the dynamics of the large-amplitude vibrational coordinate of NH 3 associated with the umbrella inversional motion. This was performed with 50 MV/cm intense quasi half-cycle THz pulses, obtained in Eq. (1) by setting the single-pulse central frequency ν = .
0 1 ps −1 and pulse duration τ = 4 ps. The left panel of Fig. 8 shows the expectation value of the inversion coordinate ρ inv (zero at planar molecular geometry) as a function of time and delay time τ between single pulses. In the region where the delay time is smaller than the single pulse duration, i.e. τ τ ≤ = 4 ps, the inversion tunnelling is signif- icantly extended or prohibited by the effects of the strong electric field. For many experiments on controlled molecules the presence of an intense electric field can lead to undesirable side effects. To avoid them, the delay time between THz pulses must exceed the single pulse duration, τ τ > , creating short windows in time with no strong electric field present. As seen in Fig. 8, for delay times just above the pulse duration it is possible to modulate, at least to some extent, the inversion tunnelling rate. For larger delay times the tunnelling dynamics quickly approach field-free behaviour.
To again validate the use of ExoMol line list data with the Hamiltonian in the dipole approximation for this type of observable, we have run calculations with an extended Hamiltonian containing non-linear contributions   from the polarizability and hyperpolarizability tensors. The results are shown in the right panel of Fig. 8 and confirm the validity of the ExoMol dipole-only approach to predict non-rigid dynamics with THz fields below 50 MV/cm.

Spontaneous emission data for Sisyphus cooling of H 2 CO and CH 3 Cl. Optoelectrical Sisyphus
cooling is a robust and general method for producing ultracold polyatomic molecules [48][49][50] . This technique works by moving molecules through a closed system of trapped ro-vibrational states in the presence of an electric field. Molecules lose kinetic energy as they travel up the electric field gradient of the trapping potential. The    cycle is repeated until sufficient cooling has taken place with sub-millikelvin temperatures possible. The speed and efficiency of Sisyphus cooling is dictated by the Einstein A coefficients and decay channels of the involved ro-vibrational energy levels; information which is readily available and computed to a high degree of accuracy in the ExoMol database. It is possible then to identify suitable transitions in different molecular systems which could be used for Sisyphus cooling. This type of analysis can aid future experiments if there is interest in a particular molecule, potentially leading to improved rates of cooling.
A general overview of the Sisyphus cooling level scheme is shown in Fig. 9. Starting from an initial, preferably highly populated ro-vibrational state, molecular population is transferred to an excited bridge state via an infrared (IR) transition. Spontaneous decay from the bridge state should be dominated by at most two channels to long-lived low-field seeking states. Large branching ratios are essential in this step with the Einstein A coefficient of the primary decay channel dictating the speed of the cooling cycle. As molecules penetrate the high-field region they climb up the Stark induced energy gradient and lose kinetic energy. For larger values of the quantum number M, the gradient is steeper and more energy is removed. It is important that the Stark split M sublevels are low-field seeking states (shifted towards higher energies) as this ensures the molecules remain trapped. Applied radio-frequency (RF) radiation stimulates emission to lower M quantum number states and the molecules gain a small amount of kinetic energy as they move down the electric field gradient. However, this gain is smaller than the previous loss in energy and the temperature of the molecules is reduced. A microwave (MW) field connects the two decay channels and closes the cycle. The process is repeated until sufficient cooling has occurred.
Sisyphus cooling has been realised experimentally for methyl fluoride (CH 3 F) 49 and formaldehyde (H 2 CO) 50 . However, a number of symmetric top molecules with strong parallel vibrational transitions have been proposed as suitable candidates for cooling using this method 48 . The closed system of ro-vibrational energy levels should satisfy the following criteria: • The initial state is highly populated and long-lived. Energy levels in the vibrational ground state with rotational quantum number J in the range 2 ≤ J ≤ 4 are generally good candidates. • There must be a sufficiently fast (intense) IR transition from the initial state to an excited bridge state with an Einstein A coefficient large enough that population transfer can occur in the time frame of effective trapping. • State-selective spontaneous decay is possible from the bridge state with one or two dominant decay channels.
The branching ratio of the two channels should sum to near unity. • Target states of spontaneous decay must be low-field seeking states with sufficiently long lifetimes to ensure the molecules remain trapped. • RF transitions from target states to lower M quantum number states are available.
• MW radiation can be used to move population from the secondary decay channel back to the initial state to close the cycle.
To illustrate how the ExoMol database, or in general any reasonably "complete" spectroscopic library, can be used to identify molecular transitions for Sisyphus cooling we have performed analysis on H 2 CO and CH 3 Cl. These systems have been treated within the ExoMol framework 19 and we refer the reader to the relevant publications for details of the ro-vibrational calculations [51][52][53][54] . Candidate ro-vibrational states were identified by sorting the .trans files according to upper state energy. There are a huge number of energy levels to process so to reduce the workload we only considered states involved in at least one intense transition (Einstein A coefficient of the order 10 1 or greater) and with J ≤ 4. Once identified, all transitions to the upper state were collected and branching ratios computed. Energy levels with two dominant decay channels or less were then analysed to check that the lower states (or target states) were low-field seeking. In Table 1 we list six collections of suitable states for cooling including those used by ref. 50 (third and forth row). Using the ExoMol database we found that the |v 1 , 3, 3〉 level has 41 decay channels, of which only transitions to the |0, 3, 3〉 and |0, 4, 3〉 levels were dominant. We predict a slightly smaller Einstein A coefficient for the |v 1 , 3, 3〉 → |0, 3, 3〉 transition and hence a slightly slower cooling rate than ref. 50. However, the difference is small (≈ 10 Hz) and our computed branching ratios are identical to those of ref. 50.
For other possible candidate states in the v 1 manifold, a slightly higher branching ratio and Einstein A coefficient are predicted for the |v 1 , 4, 4〉 → |0, 4, 4〉 transition. Alternatively, the v 2 C-O stretching mode appears suitable for Sisyphus cooling but using these states would half the speed of cooling compared to the v 1 levels. All lower states E″ in Table 1 belong to the vibrational ground state and thus remain accessible either thermally or by methods of selective state preparation. At T = 100 K for example, lower states are thermally populated to around 30%. Note that only transitions allowed by selection rules were considered, i.e. Δ J, Δ M = 0, ± 1 and Δ K = 0. CH 3 Cl. Methyl chloride has not yet been considered for Sisyphus cooling but we expect a similar experiment to CH 3 F 49 could be performed. As shown in Table 2, ro-vibrational levels of the v 1 symmetric CH 3 stretching mode and vibrational ground state appear the best candidates for cooling. Similar to methyl fluoride, which had a spontaneous decay rate of about 15 Hz for the |v 1 , 3, 3〉 → |0, 3, 3〉 transition, the respective rate in CH 3 Cl is only marginally larger and would hence provide a similar speed of cooling. Note that we have only considered CH 3 35 Cl for the present analysis but we would expect the same conclusions for CH 3 37 Cl.

Discussion
Comprehensive spectroscopic line lists generated as part of the ExoMol project can be used to model ro-vibrational phenomena for small polyatomic molecules in the presence of external electric fields. By doing so we avoid repeating the same, often expensive calculations. Currently the ExoMol database contains complete sets of high accuracy ro-vibrational energies, transition frequencies, Einstein A coefficients, and complex dipole moment phases for many important diatomic and small polyatomic molecules valid for temperatures up to T = 1500 K. These data can be straightforwardly utilized to model molecule-field interactions in the dipole approximation. In the future we plan to extend the ExoMol database with Raman and electric quadrupole transition moments for selected molecules, such as NH 3 and H 2 O, thereby extending the range of possible applications.
Two illustrative examples of molecule-field interaction applications were presented. Strong THz field induced ro-vibrational dynamics of PH 3 and NH 3 were simulated and these represent the first high accuracy calculations on polyatomic molecules which have been reported in the literature. For optoelectrical Sisyphus cooling, suitable collections of states were identified in H 2 CO and CH 3 Cl. In both molecules the v 1 symmetric stretching mode provided the fastest cooling rates and most suitable ro-vibrational energy levels for Sisyphus cooling.
Although we have chosen relatively straightforward processes to look at, our approach can simulate any case with no constraint on the availability of data (provided the molecule is in the ExoMol database). At present, we are unaware of any other line list database which stores the dipole moment phase factors. This information is always produced with the dipole line strengths or Einstein coefficients but is discarded when the line lists are compiled, thus making them unusable for modeling field-driven effects. We therefore encourage other theory-based data compilations to retain this information in future.

Methods
ExoMol line list data structure. A complete description of the ExoMol data structure along with examples was recently reported in ref. 11. Therefore we provide only a brief description that is relevant for the present work. The two main files available for download from the ExoMol website (www.exomol.com) are the states . and .trans files, unique to each molecule. The states . file contains all computed ro-vibrational energy levels (in cm −1 ). For polyatomics, which are the focus of the present study, each energy level has a unique state ID with symmetry and quantum number labelling. The trans . files, which are split into frequency windows as to be more manageable, contain all calculated molecular transitions. Upper and lower state labels, transition frequencies and Einstein A coefficients are provided.
Modelling Stark associated phenomena requires information on the complex phase of the dipole moment matrix elements (discussed below). This information, simply a + or − , is not currently available for all molecules in the ExoMol database. However, work is underway to rectify this and the complex phase for each molecular transition can be extracted from the dipole . file. Several other molecule specific files related to atmospheric applications (e.g. pressure broadening of spectral lines) are available but since we do not use them in this work, we do not discuss them here.
Electric field interaction Hamiltonian. The complete Hamiltonian is written as a sum of the field-free ro-vibrational Hamiltonian and the molecule-field interaction Hamiltonian, which is approximated in this work as the permanent dipole interaction Here ε l,J and |l, J, m〉 denote the field-free ro-vibrational energy and wavefunction, respectively, of a state identified by the running number l (ID number) in the line list energy file (.states file), J is the rotational quantum number of the total angular momentum, and the quantum number m = − J, … , J is the projection of the total angular momentum on the laboratory-fixed Z axis. The coupling strength is given by the product of the ro-vibrational matrix elements of electronic ground state electric dipole moment µ ′ ′ ′ A l J m l J m ( , , , , , ) in the ro-vibrational basis and the time-dependent classical electric field E A (t), defined in the laboratory fixed axes system. To describe the radiative fields we use the form , where A is the polarization axis, E t t T ( ; , ) is the pulse time profile with a maximum at t = t 0 and duration (FWHM) T, and ω A is the carrier frequency. Calculation of the µ A ro-vibrational matrix elements from Einstein A coefficients, extracted from the ExoMol line list, is described in detail below. The selection rules for µ A allow coupling between ro-vibrational states with Δ J = J′ − J = 0, ± 1 and Δ m = m′ − m = 0, ± 1. Hence, the total Hamiltonian may be constructed in the form of a block-tridiagonal matrix, having blocks corresponding to Δ m = − 1, 0, or + 1 in the lower, main, and upper diagonals, respectively. Using this representation, various linear algebra operations can be performed efficiently and the memory requirements reduced by use of the band matrix storage scheme. For A = Z, µ Z can couple only states with Δ m = 0, thus the total Hamiltonian becomes factorized into independent blocks for each chosen m, which are processed independently.
The time-evolution of the wavefunction is given by the time-dependent Schrödinger equation with the Hamiltonian in Eq. (2), which may be time-independent or time-dependent. For static fields the problem reduces to solving an eigenvalue equation for the Hamiltonian. For radiative fields the time-evolution of the wavefunction is described by the time-evolution operator U(t, t′ ) as Ψ Ψ = ′ ′ t U t t t ( ) ( , ) ( ). Providing Δ t = t − t′ is sufficiently small compared to the characteristic oscillation period of the field perturbation, U(t, t′ ) can be evaluated, for example, using the split-operator method as  3) is the evaluation of the matrix exponential at each time step t = t′ + Δ t/2, which in our case is computed using the iterative approximation based on Krylov subspace methods, as implemented in the Expokit computational package 55 . Alternative and more sophisticated representations for the split-operator technique can be found in refs 56 and 57.
Dipole matrix elements. We aim to derive the relationship between the Einstein A coefficients, A fi , and the matrix elements of the electric dipole moment operator, ψ µ ψ 〈 | | 〉 f A i ( = A X Y Z , , ), defined relative to the laboratory fixed system. It is convenient to first express A fi in terms of the dipole line strength S fi as π πε  where v and J, k, m refer to the vibrational and rotational quanta, respectively, μ α (α = x, y, z) is the dipole moment in the molecule-fixed axes system, and T σα is the Cartesian-to-spherical tensor transformation matrix It should be noted that the actual wavefuncitons |ψ〉 are expansions in terms of |J, k, m〉 |v〉 as basis functions. We omit the expansion coefficients in Eq. (5) for the sake of simplicity, which does not alter the methodology presented below. Using the same variables, the matrix element of the dipole moment operator ψ µ ψ 〈 | | 〉