Emulating optical cycling centers in polyatomic molecules

An optical cycling center (OCC) is a recently coined term to indicate two electronic states within a complex quantum object that can repeatedly experience optical laser excitation and spontaneous decay, while being well isolated from its environment. Here we present a quantitative understanding of electronic, vibrational, and rotational excitations of the polyatomic SrOH molecule, which possesses a localized OCC near its Sr atom. In particular, we describe the vibrationally-dependent trends in the Franck-Condon factors of the bending and stretching modes of the molecular electronic states coupled in the optical transition. These simulations required us to perform electronic structure calculations of the multi-dimensional potential energy surfaces of both ground and excited states, the determination of vibrational and bending modes, and corresponding Franck-Condon factors. We also discuss the extent to which the optical cycling center has diagonal Franck-Condon factors.

An engineering approach, however, can be used to add optical cycling centers (OCCs) to a variety of polyatomic molecules with increasing complexity [22]. For this approach to work we need to ensure that the electrons holding the optically active atom and the molecule together remain sufficiently localized so that the whole molecular system can be translationally cooled. This stringent requirement motivates our intent to conduct an assessment of the role of electron density and localization in the coupling to the ligand molecule.
Recent theoretical studies on the cooling properties of alkaline-earth monohydroxides other than SrOH can be found in Refs. [19,23]. They only determined molecular parameters and Frank-Condon factors for vibronic transitions between ground and excited states near equilibrium geometries.
Our objectives are to better quantify the extent to which polyatomic SrOH is an ideal platform for cooling with laser light as well as to determine the "global" shape of its ground and excited potential energy surfaces (PES) in anticipation of future research on dissociating SrOH into Sr and OH and on ultracold scattering between Sr and OH. In this paper, we describe a complex electronic structure determination of local and some global properties of four multi-dimensional intramolecular PESs of SrOH as well as calculate their corresponding vibrational structures. Our effort involves locating potential minima, potential avoided crossings as well as conical intersections (CIs), where two adiabatic PESs of the same electronic symmetry touch. We represent the PESs and thus the CIs both in terms of adiabatic and diabatic representations. We also determine Frank-Condon factors for not-only the lowest vibrational states but also higher-excited vibrational states of the potentials, where the transitions are no longer diagonal due to non-adiabatic and anharmonic corrections. At near linear geometries, we evaluate the Renner-Teller (RT) parameter [24] for some excited-state potentials.

Results
Electronic potentials of SrOH. The computation of PESs is crucial for defining the landscape in which the nuclei transverse upon interacting with one another. These surfaces can have many features in the context of their topology that are important for the internal transition mechanisms. We begin by describing our calculation of the ground and excited PESs of strontium monohydroxide SrOH. Past theoretical studies on M-OH were devoted to either spectroscopic characteristics near equilibrium geometries [25] or the PESs for the lighter BeOH and MgOH with fixed O-H separation [26]. Experimental studies of CaOH [10] and SrOH [2,7,8] predominantly focussed on their near-equilibrium ro-vibrational structure.
The relevant multi-dimensional PESs have been determined using a combination of coupled-cluster (CCSD(T)), equation-of-motion coupled-cluster (EOM-CC), and multireference configuration interaction (MRCI) methods. This combination allows us to overcome the limitations of the CCSD(T) and EOM-CC methods associated with their single reference nature and the MRCI method with its limits on the size of active space and characterize multiple avoided crossings between relevant PESs. Since the molecule contains one heavy atom, relativistic effects are embedded into the description of the core electrons.
Potentials are presented in the Jacobi coordinates R and r defined in Fig. 1. For our purposes it is sufficient to determine the PESs, which are only functions of (R, θ, r) in the two-dimensional plane with the OH separation fixed to its diatomic equilibrium separation as the energy required to vibrationally excite OH is nearly seven times larger that than those of the Sr-OH stretch and bend. Coupling to the OH stretch can then be ignored for our objectives. (We always use r = 1.832a 0 and also suppress the r dependence in our notation.) The PESs and corresponding electronic wavefunctions are labeled by their total electronic spin angular momentum, here always a doublet or spin 1/2, as well as the irreducible representations of point groups C ∞v for linear geometries and C s otherwise. In particular, we obtain four 2 A and two 2 A potentials using standard notation for the irreducible representations of C s , respectively. For linear geometries 2 A states connect to either 2 Σ + or 2 Π states of the C ∞v symmetry group, while 2 A states always become 2 Π states. Hence, we label   The optimized geometry for the ground-state SrOH molecule occurs at θ = 180 • . Its electronic wavefunction has X 2 Σ + symmetry and is ionically bonded. These observations are consistent with spectroscopic data [27] and the semi-empirical analyses of Ref. [28]. The The fourth diabatic PES is the shallow excited 4 2 A (F 2 Π) potential with a covalently bonded electronic state that dissociates to Sr( 1 S) and OH(X 2 Π). For large R the ionic PESs correlate to electronically-excited states of the Sr atom. Crucially, we observe that the three ionic diabatic PESs each have a curved line in the (R, θ) plane along which the 4 2 A (F 2 Π) potential equals the energy of the ionic potential. It is worth noting that to good approximation this curve is independent of θ. The 4 2 A (F 2 Π) potential is expected to play a crucial role in the zero-energy dissociation of SrOH to create ultracold OH fragments, an important radical for various scientific applications. For SrOH the BO approximation breaks down when potentials of 2 A symmetry come close and even become degenerate for one or more geometries. Couplings between these potentials can not then be neglected. The corresponding non-adiabatic transition probability for "hopping" between potentials increases dramatically especially for conical intersections and can greatly affect chemical properties [32].
Our CCSD(T) and EOM-CC calculations only resulted in diabatic PESs V m (R, θ) and electronic wavefunctions |φ m , due to their reliance on a single reference configuration. Here, index m labels states and corresponding PESs. As the diabatic electronic wavefunctions are predominantly described by this single reference configuration, their dependence on R and θ is weak and negligible, and suppressed in our notation.
Coupling matrix elements between the diabatic potentials are constructed by performing  Mixing angles similarly determined from our calculations for the 1 2 A and 2 2 A diabatic states.
No conical intersection is present. c) Adiabatic 1 2 A (blue surface) and 4 2 A (red surface) potential energy surfaces V as functions of separation R and and angle θ near the conical intersection at two or more reference configurations and thus do lead to adiabatic BO potentials U n (R, θ) and adiabatic electronic wavefunctions |ψ n (R, θ) , where index n labels states. Adiabatic wavefunctions strongly depend on R and θ near crossings or avoid crossings between potentials. We also compute diagonal overlap matrix elements ψ n (R , θ )|ψ n (R, θ) at different geometries within the MRCI method.
Numerically, we find that the general shapes of the CCSD(T) and MRCI potentials are the same. Specifically, the geometries where diabatic potentials in Fig. 2a cross and adiabatic MRCI potentials avoid are in good agreement. We can then assume that away from the crossings |ψ n (R, θ) ≈ |φ m for some diabatic index m. Near avoided crossings the adiabatic wavefunction is a superposition of diabatic wavefunctions.
We focus on the coupling matrix element between the diabatic 1 2 A and 4 2 A states near R = 8a 0 as well as those between 1 2 A and 2 2 A near R = 6a 0 . We assume the coupling between the diabatic 2 2 A and 4 2 A states near R = 6a 0 is the same as between the A ones. For each pair of states the adiabatic wavefunctions can then be approximated as the with indices n, n and m, m , respectively, where mixing angle ϑ(R, θ) is to be determined.
To determine the mixing angle we, first, identify a reference geometry R ref for each θ away The mixing angle at other radial geometries is then given by and the coupling matrix element between the two diabatic states m and m is In the limit V m (R, θ) → V m (R, θ) this expression needs to be treated carefully. The mixing angle will approach π/4 and only at CIs W mm (R, θ) = 0.
In practice, we parametrize ϑ(R, θ) as we perform MRCI calculations only on a restricted set of geometries (R i , θ j ). We ensure a smooth functional dependence on R at fixed θ by where R c (θ) is the curve where V m (R, θ) = V m (R, θ) and R 0 (θ) ≥ 0 is a coupling width that is fitted to reproduce the MRCI overlap matrix at θ j (at other θ it is found with the Akima interpolation method [34].) By construction R 0 (θ) is much smaller than |R ref − R c (θ)|. We use R ref = 11a 0 and 7a 0 for the pair of 2 A and 2 A states, respectively. Values for θ j , R c , and R 0 are available upon request. We calculate the RT parameter by estimating the non-adiabatic coupling between 2 2 A and 1 2 A PESs. These PESs are degenerate at co-linear geometries but split otherwise. difference. In the harmonic approximation this differential splitting is characterized by the Renner-Teller parameter [27] = between either X 2 Σ + ↔ A 2 Π or X 2 Σ + ↔ B 2 Σ + vibronic states. The success of the experiment is associated with near-diagonal FCFs between the sets of levels supported by these electronic potentials, which ensures that spontaneous emission from a rovibrational state of the excited electronic state populates with near unit probability the corresponding rovibrational state of the ground electronic state. This two-level system absorbs and emits many photons to achieve cooling [13].
We observe a dramatic decrease of the FCFs for v s = 6 and 7, emphasized by a blue shaded area in Fig. 6a. We explain this behavior in terms of the complex overlap of vibrational wavefunctions in R and θ. Figure 6b shows bending range for v s < 4 vibrational levels in all three electronic states is ≈ 15 • and that for v s > 5 the wavefunction is more extended. In fact, for v s > 5 the angular extent for the 1 2 A (X 2 Σ + ) and 2 2 A (A 2 Π) states remains close but differ significantly from those of the 3 2 A (B 2 Σ + ) state. This irregularity in the v b = 0 series coincides with the corresponding irregularity in FCFs, but clearly is not the whole story as we could conclude that the transitions. We observe that for the (3,0,0) vibrational level the wavefunction of the three electronic potentials have similar shape, although the overlap for the 1 2 A (X 2 Σ + ) and 3 2 A (B 2 Σ + ) states is significantly better. For the (7,0,0) level neither overlap is good.
Combining our observations on the angular extend in Fig. 6b with those on the stretching direction then explains the dramatic decrease of the FCFs for v s = 6 and 7.
The (0, 0, 0) to (0, 0, 0) FCFs for these electronic states was measured in Ref. [28]. They found 0.957(3) and 0.977(2) for the 1 2 A (X 2 Σ + ) to 2 2 A (A 2 Π) and 1 2 A (X 2 Σ + ) to 3 2 A (B 2 Σ + ) transition, respectively. These value are smaller than those found in our calculations. Better agreement might be found by allowing the OH bond to change. Extended tables of the FCF matrices are presented in Supplementary Note 1.

Discussion
We have quantitatively analyzed the unique localized chemical bond between valence electrons of Sr and OH radical. We have shown that one of the two Sr valence electrons remains localized on the Sr + cation and can be optically excited without disturbing the atom-ligand bond leading to highly-diagonal FCFs and efficient optical-cycling conditions. The goal of our theoretical work has been to analyze the cooling process based on a rigorous calculation of the electronic structure of this important molecule previously totally unknown, as a means to look into the nature of optical cycling centers in polyatomic molecules and help extend its generality. Furthermore, our work indicates that there are several electronic and vibrational states that can be used for the laser cooling with FCFs close to unity. We also explain why the X-B transitions have more diagonal FCFs than X-A transitions.
Our analysis of the electronic and vibrational energy landscape of SrOH also points to intriguing molecular features that require further explorations. For instance, we located conical intersections in SrOH that can have important implications for non-adiabatic transitions between PESs, collisional dynamics within this molecule, as well as its fragmentation.
Our characterization of the geometric Renner-Teller effects provides better understanding of the coupling between electronic and vibrational motion and accounts for the possible nominally-forbidden loss channels.
Finally, Ref. [21] proposed that replacing the hydrogen atom in the monohydroxide with larger ligand should not significantly disturb the valence electron of metal-cation so that optical cycling and cooling even larger polyatomic molecule might remain possible. Detailed confirmation of this idea also awaits future research.

Methods
Ab initio methods for the SrOH potential energy surfaces. We have used the coupled-cluster CCSD(T) method from CFOUR [37] to calculate the ground and first excited diabatic PESs with total electron spin 1/2 of A symmetry (ground state), A symmetry (first excited state), and the highly excited diabatic PES with spin 1/2 of A symmetry but with completely different reference configuration. They are labeled 1 2 A , 1 2 A and 2 2 A , respectively. The CCSD(T) calculations are based on the atomic bases sets def2-QZVPP (8s8p5d3f)/[7s5p4d3f] [38] paired with the Stuttgart ECP28MDF effective core potential [39] for Sr, the all electron aug-cc-pVQZ (13s7p4d3f2g)/[6s5p4d3f2g] for O, and aug-cc-pVQZ (7s4p3d2f)/[5s4p3d2f] for H [40].
For all excited diabatic potentials, 1 2 A , 2 2 A , 3 2 A , 2 2 A , and 4 2 A , we employed the equation-of-motion coupled-cluster (EOM-CCSD(dT)) approach, implemented in Q-CHEM, version 5.0 [41]. The excitation-energy root space for the EOM-CCSD(dT) method spans four roots of A symmetry and two roots of the A symmetry. We use the Stuttgart ECP28MDF relativistic effective core potential [39] along with Peterson's pseudopotential based, polarized valence, correlation consistent triple-zeta (aug-cc-pVTZ-PP) basis for the Sr atom [42] and corresponding Dunning's aug-cc-pVTZ basis for the O and H atoms [40].
These calculations allowed us to obtain the relative splittings between A and A PESs and between the 2 2 A and 3 2 A PESs.
The CC based methods used in construction of the 2 A and 2 A potentials are capable of providing highly accurate quantitative data about PESs. Unfortunately, these methods are unsuitable for the calculation of PESs near their crossing or avoided crossing. The way out is a MRCI calculation, which realistically, can only be done with an active space of limited size.
We characterize non-adiabatic couplings between states by using the MRCI calculations with single and double excitations, implemented in MOLPRO [43]. The ECP28MDF Köln pseudo-potential [39] with the augmented correlation-consistent triple-zeta basis aug-cc-pVTZ-PP [42] for Sr and the all-electron aug-cc-pVQZ basis for O and H [40] are applied.
The calculations are performed with reference orbitals constructed from state-averaged multi-reference self-consistent-field calculations. Finally, all electronic structure calculations were performed with assumption that the OH separation is fixed at r = 1.832a 0 .
Interpolation of potential energy surfaces. The four 2 A and two 2 A two-dimensional diabatic SrOH potentials have been computed on a grid in Jacobi coordinates R and θ. In order to find the potential energies at other coordinates we expand the angular dependence of each diabatic PESs in terms of Legendre polynomials P l (cos θ) via for R < 16.5a 0 with L = 9 and radial expansion coefficients U l (R). For separations R > 16.5a 0 the PESs are fit to with dispersion coefficients C nl . The radial coefficients U l (R) and C nl are determined by the Reproducing Kernel Hilbert Space (RKHS) interpolation method [44]. The two expansions are smoothly joined using Within these approximations we have computed several of the energetically-lowest vibrational levels with total angular momentum J p = 0 + and 1 ± of the diabatic 1 2 A (X 2 Σ + ), , and 3 2 A (B 2 Σ + ). We focus on these rotational states as J = 1 → J = 0 transitions from the 1 2 A (X 2 Σ + ) ground state are used in laser cooling of SrOH [20].
In Table III  reproduced to better than 1% for the 1 2 A (X 2 Σ + ) and 1 2 A (A 2 Π) states. The discrepancy increases to 2.5% for the 3 2 A (B 2 Σ + ) state.
For the bending mode with quanta v b our theoretical energy differences are larger than those found experimentally, although the difference is no more than 10% for the three electronic states. Adding a second quanta in the bending mode shows that anharmonic contributions are non-negligible. We also corroborate the DFT and 2D-DVR results of Nguyen et al. [28]. Nguyen's DFT results are based on the harmonic approximation.
The vibrational wavefunctions are used the compute FCFs among various pairs of ionic states. Their values and trends were presented in subsection "Wavefunction overlaps". Additional bound-state energies and wavefunctions are shown in Supplementary Figures 1-3.
Our determination of the permanent dipole moments of the 1 2 A (X 2 Σ + ), 2 2 A (A 2 Π), and fine-structure resolved measurement, i.e. they found permanent dipole moments of 0.590 (45) Debye and 0.424(5) Debye for Ω = 1/2 and 3/2, respectively, suggesting that our value is too small. (The numbers in parenthesis are one-standard-deviation uncertainties in the last digits.)

Data Availability
The data sets generated during the current study are partially included in this published article (and its Supplementary Note 1) and also available from the corresponding author on reasonable request.   , v b , 0) stretching and bending levels with total angular momentum J = 0 and parity +1 for the 1 2 A (X 2 Σ + ), 1 2 A (A 2 Π), and 3 2 A (B 2 Σ + ) states of the 88 Sr-16 O 1 H isotopologue using our diabatic and adiabatic potential energy surfaces, and data from the literature. Energies are relative to that of the (0, 0, 0) mode and in units of cm −1 .