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 optical cycling center (OCC) is a recent term describing two electronic states within a quantum object undergoing repeated optical laser excitation and spontaneous decay, while being isolated from its environment. The authors present ab-initio calculations of the ground and excited electronic and vibronic states of the polyatomic molecule SrOH providing detailed understanding of the complex molecular cooling processes with an OCC.

L aser cooling and trapping of atoms, enabled by the existence of closed optical cycling transitions, have revolutionized atomic physics and led to breakthroughs in several disciplines of science and technology 1 . These advances enabled the simulation of exotic phases in quantum-degenerate atomic gases, the creation of a novel generation of atomic clocks, matter-wave interferometry, and the development of other highly sensitive sensors. Temperatures below tens of microkelvin have also allowed the confinement of diatomic molecules, built from or associated with laser-cooled atoms, in electric, magnetic, and/or optical traps, where they are isolated from their environment and can be carefully studied.
Achieving similar temperatures for polyatomic molecules, however, remains challenging. Since polyatomic molecules are characterized by multiple degrees of freedom and have correspondingly more complex structures, it is far from obvious whether there exists polyatomic molecules that have the nearlyclosed optical cycling transitions required for successful laser cooling. Such transitions could then repeatedly scatter photons.
A diverse list of promising applications for ultracold polyatomic molecules exists. This includes creating novel types of sensors, advancing quantum information science, simulation of complex exotic materials, performing precision spectroscopy to test the Standard Model of particle physics, and, excitingly, the promise of control of quantum chemical reactions when each molecule is prepared in a unique rovibrational quantum state. Moreover, the de Broglie wavelength of colliding ultracold molecules is much larger than the range of intermolecular forces and, thus, the science of the breaking and making of chemical bonds has entered into an unexplored regime.
Over the decades many spectroscopic studies of molecules consisting of an alkaline-earth metal atom (M) and a ligand have been performed 2-10 , predominantly to determine their structure. The simplest polyatomic molecules of this type, the triatomic alkaline-earth monohydroxides M-OH, have attracted particular attention after the discovery of their peculiar ionic chemical bond. When a ground-state alkaline-earth atom is bound to OH one of its two outer-most ns 2 electrons is transferred to the ligand, leaving the second electron in an open shell molecular orbital localized around the metal atom. This electron can be optically excited without disturbing the atom-ligand bond leading to so-called highly diagonal FCFs and efficient optical cycling. High-resolution laser-spectroscopy experiments on CaOH (and CaOD) 5 and SrOH (and SrOD) 11 were first to deduce the strong ionic character of the metal-hydroxide bond. Reference 12 showed that the remaining valence electron of the metal atom can be easily promoted to any excited orbital.
Recently, Doppler laser cooling 13 of the simpler CaF, SrF, YbF, and YO dimers has been demonstrated [14][15][16][17][18] . Their success is associated with nearly-diagonal Franck-Condon factors (FCFs) on an allowed optical molecular transition. Diagonal FCFs ensure that spontaneous emission from vibrational state v of the excited electronic state populates with near unit probability vibrational state v 0 ¼ v of the ground electronic state. This closed two-level system can absorb and emit thousands of photons per molecule to achieve cooling.
Following these successful molecular cooling experiments, Isaev and Berger 19 suggested that similar near-diagonality of FCFs exists in other classes of polyatomic molecules. It is now understood that metal-monohydroxides and even larger metalmonoalkoxide molecules with metal atom as Sr, Ca, or Ba are promising candidates for laser cooling. In 2017 the first optical cycling transitions in the polyatomic monohydroxide molecule SrOH was demonstrated by Dr. Doyle's group 20 . They succeeded in reducing the translational motion of SrOH to the record low temperature of 750 microkelvin starting from 50 milikelvin and using the near unit values of the FCFs. To eliminate rotational branching during the photon cycling process experimentalists use of the rotationally closed J 00 ! J 0 À 1 angular momentum transitions (see for details refs. 20,21 ). A significant effort from the scientific community, however, is needed to identify and study the classes of polyatomic molecules with closed optical cycling transitions. Laser cooling of polyatomic molecules is relevant to those molecules that are able to scatter hundreds of photons (stimulated absorption followed by a spontaneous emission) between two vibrational states without loss to other vibrational states, i.e., have a cycling transition where for each cycle the kinetic energy of the molecule is reduced by ΔE=k ¼ 1μK and 10 μK depending on the cooling process and k is the Boltzmann constant. The requirement of 100 scattering cycles is to a certain degree arbitrary but reasonable and implies that FCF larger than 0.99 are needed. A larger FCF will allow better cooling.
In fact, a comprehensive analysis of level structures that are amenable to laser cooling has to be conducted. Special attention must be given to finding excited-state molecular potentials that have: i) a shape that is similar to that of the ground-state potential; ii) a strong dipole electronic transition with the ground state; and iii) weak transitions to dark states that are not optically coupled to the excited state. It is to be expected that only a small number of polyatomic molecules have these three properties.
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 higherexcited 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 rovibrational structure.
The relevant multi-dimensional PESs have been determined using a combination of coupled-cluster (CCSD(T)), equation-ofmotion coupled-cluster (EOM-CC), and multi-reference 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 1v for linear geometries and C s otherwise. In particular, we obtain four 2 A 0 and two 2 A 00 potentials using standard notation for the irreducible representations of C s , respectively. For linear geometries 2 A 0 states connect to either 2 Σ þ or 2 Π states of the C 1v symmetry group, while 2 A 00 states always become 2 Π states. Hence, we label potential surfaces by n 2 A 0;00 ðm 2 Λ ± Þ, where n ¼ 1; 2; 3; and m ¼ X; A; B; indicate the energy ordering of states for R and r close to their equilibrium separations and θ ¼ 180 .
CCSD(T) and EOM-CC calculations lead to the most-accurate PESs for each irreducible representation and, in principle, should return the corresponding adiabatic PES with avoided crossings from excited potentials in the same irreducible representations. Numerically, however, we find this not to be true due to the fact that the method use only excitations from a single reference configuration and that the electronic wavefunctions rapidly change from ionic to covalently bonded with small changes in the Jacobi coordinates, in particular regions of R. In essence, for a given reference configuration the CCSD(T) or EOM-CC simulations return diabatic potentials with electronic wavefunctions that have either an ionic bond, where one of the outer valence electrons of Sr is now tightly bound to oxygen, or a covalent bond, where both valence electrons mostly remain in orbit around the Sr nucleus and barely couple with the OH electron cloud. The diabatic PESs have lines in the plane ðR; θÞ along which two potentials have the same energy. Figure 2a shows four diabatic PESs of SrOH as functions of R and θ obtained by CCSD(T) and EOM-CC methods. Details regarding electronic basis sets, coupled-cluster method for the individual 2 A′ and 2 A″ PESs, and interpolation can be found in Methods. The calculations are performed at ten angular and about 45 radial geometries. This data is then interpolated using an analytical functional form. The interpolated PESs are essential in recognizing system characteristics, such as minima and transition states, i.e., saddle points, as well as crossings. Near extrema the curvature or Hessian matrix is evaluated.
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 2 2 A 0 and 3 2 A 0 states also have ionic character and minima at θ ¼ 180 . In fact, the minima of these three potentials have nearly the same equilibrium coordinates and curvatures providing excellent conditions for optical cycling. These conditions are further described in subsection "Wavefunction overlaps".
The fourth diabatic PES is the shallow excited 4 2 A 0 ð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 0 ð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 0 ð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. Figure 2b shows the interpolated adiabatic 2 A 0 PESs of the SrOH molecule as functions of R for six values of θ and r ¼ 1:832a 0 obtained by diagonalizing the 4 4 Hamiltonian with the diabatic PESs as diagonal matrix elements and coupling matrix elements described in subsection "Non-adiabatic couplings and conical intersections" at each ðR; θÞ. These cuts through the PESs exhibit several intriguing features. First, for θ ¼ 0 and 180 CIs, where two potentials still touch, are apparent. For other values of θ the crossings are avoided. The location of the CIs at linear geometries is specific to SrOH, other molecules will have CIs at other geometries. CIs lead to geometric or Berry phases 29 , i. e. sign changes associated with an electronic wavefunction when transported along a closed path encircling a CI. These phases lead to interference that allows efficient non-adiabatic transitions between surfaces 30 , modify product rotational state distributions in chemical reactions, and influence ro-vibrationally averaged transition matrix elements.
We have similarly obtained two diabatic 2 A 00 PESs. Both are 2 Π states at linear C 1v geometries. On the energy scale of Fig. 2a their shapes are nearly indistinguishable from those of the 2 2 A 0 ðA 2 ΠÞ and 4 2 A 0 ðF 2 ΠÞ potentials. In fact, 1 2 A 00 ðA 2 ΠÞ and 2 2 A 0 ðA 2 ΠÞ as well as 2 2 A 00 ðF 2 ΠÞ and 4 2 A 0 ðF 2 ΠÞ are degenerate for linear geometries. Unlike the corresponding 2 A 0 potentials, however, the two 2 A 00 potentials do not possess CIs.
Non-adiabatic couplings and conical intersections. The adiabatic 2 A 0 PESs in Fig. 2b are Born-Oppenheimer (BO) potentials 31 that are typically used as a starting point in describing scattering and chemical processes. The BO approximation is based on the realization that the motion of nuclei and electrons occur on different time or energy scales and usually only a single PES is required. This leads to a significant simplification of the description of scattering and chemical processes. For SrOH the BO approximation breaks down when potentials of 2 A 0 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 less-accurate MRCI calculations at the SD level with a basis set similar to that used for our CCSD(T) calculations (details are given in Methods.) Our MRCI calculations rely on 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 hψ n ðR 0 ; θ 0 Þjψ n ðR; θÞi 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 0 and 4 2 A 0 states near R ¼ 8a 0 as well as those between 1 2 A 00 and 2 2 A 00 near R ¼ 6a 0 . We assume the coupling between the diabatic 2 2 A 0 and 4 2 A 0 states near R ¼ 6a 0 is the same as between the A 00 ones. For each pair of states the adiabatic wavefunctions can then be approximated as the superposition 33 ψ n ðR; θÞ ψ n 0 ðR; θÞ In the limit V m ðR; θÞ ! V m 0 ðR; θÞ this expression needs to be treated carefully. The mixing angle will approach π=4 and only at CIs W mm 0 ð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 using where R c ðθÞ is the curve where V m ðR; θÞ ¼ V m 0 ð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 jR ref À R c ðθÞj. We use R ref ¼ 11a 0 and 7a 0 for the pair of 2 A 0 and 2 A 00 states, respectively. Values for θ j , R c , and R 0 are available upon request. Figure 3a and b show mixing angles ϑðR; θÞ for the two diabatic 2 A 0 and the two 2 A 00 states, respectively. In both panels the mixing angle is seen to change relatively rapidly over a small range of R, i.e., determined by R 0 ðθÞ, from nearly 90 to 0 . When ϑðR; θÞ ¼ 45 we have R ¼ R c ðθÞ. For the 2 A 0 states in panel a, R 0 ðθÞ is largest for θ % 90 and is zero for θ ! 0 and 180 , indicative of CIs where ϑðR; θÞ has an infinitely sharp jump from 90 to 0 . For the 2 A 00 states in panel b the width R 0 ðθÞ is nearly independent of θ and there is no CI. Finally, Fig. 3c shows the adiabatic 2 A 0 potentials near the CI at R % 8a 0 and θ ¼ 180 as determined by diagonalizing the 2 2 matrix containing the relevant diabatic Fig. 2 The ground and excited potential energy surfaces VðR; θÞ of SrOH at separation r ¼ 1:832a 0 . a A two-dimensional cut through the four relevant 2 A 0 diabatic potential surfaces as functions of separation R and angle θ. The blue, orange, green, and red curves correspond to states 1 2 A 0 ðX 2 Σ þ Þ, 2 2 A 0 ðA 2 ΠÞ, 3 2 A 0 ðB 2 Σ þ Þ, and 4 2 A 0 ðF 2 ΠÞ, respectively. Seams, containing conical intersections at the two linear geometries, are lines along which two diabatic potentials are degenerate. b Adiabatic 2 A 0 potentials as functions of separation R for six values of angle θ. The absolute ground state minimum occurs at angle θ ¼ 180 at separation R % 4a 0 . The zero of energy is set at the energy of a ground-state Sr atom infinitely far away from the ground-state OH(X 2 Π) molecule at its equilibrium bond length. Potential energies are expressed in units of cm À1 , using Planck's constant h and speed of light in vacuum c.
PESs coupled by the coupling matrix element at each ðR; θÞ. The figure corresponds to a surface plot of the data shown in Fig. 2b.
The Renner-Teller effect 35,36 occurs in linear or quasi-linear, open-shell triatomic molecules and relies on non-adiabatic coupling between rovibrational and electronic motion in polyatomic molecules. It is associated with an energy degeneracy of electronic states for high-symmetry geometries that is lifted by vibrational bending motion. The motion causes degenerate potential energy surfaces to split and breaks the Born-Oppenheimer (BO) approximation. Stretching modes, which do not break linear symmetry, are not affected.
We calculate the RT parameter by estimating the non-adiabatic coupling between 2 2 A 0 and 1 2 A 00 PESs. These PESs are degenerate at co-linear geometries but split otherwise. Figure 4 shows the potentials near the co-linear equilibrium separation R as functions of the bending angle θ. The potentials are degenerate at θ ¼ 180 . The slight anisotropy away from linear geometry causes the bending vibrational motion in the two PESs to differ. Our calculated ð0; 1; 0Þ bending energy relative to ð0; 0; 0Þ for the 2 2 A 0 and 1 2 A 00 states are ωðA 0 Þ ¼ 384 cm À1 and ωðA 00 Þ ¼ 412 cm À1 , respectively, corresponding to a 28 cm À1 difference. In the harmonic approximation this differential splitting is characterized by the Renner-Teller parameter 27 ϵ ¼ ωðA 0 Þ 2 À ωðA 00 Þ 2 ωðA 0 Þ 2 þ ωðA 00 Þ 2 ¼ À0:0693 : Presunka and Coxon 27 estimated a Renner-Teller parameter of À0:0791 for these states from fitting this and other non-adiabatic parameters to spectroscopic data. Here, we calculate ϵ for the first time from ab-initio calculations and find agreement with spectroscopic estimate to within 15%.
Wavefunction overlaps. Laser cooling of SrOH 20 relied on optical cycling transitions 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 . Table 1 shows atomic orbital composition for the ground and A and B electronically excited states of SrOH. The X electronic ground state is mostly described by Sr-s orbital as expected. One can notice that in the case of the A and B states there is some participation of d-orbitals of Sr atom. For the A state that participation is quite small and the excited state is mainly formed by Sr-p orbital participating in s À p transition. The B state has quite a substantial component originating from Sr-d-orbitals indicating increased s À d type hybridization with some participation of O-s and H-s atomic orbitals.
Our computation of vibrational wavefunctions and their overlap for electronic transitions is described in Methods. We focus on vibrational levels near the global minima of the diabatic PESs at θ ¼ 180 , as shown in Fig. 5. Effects from diabatic couplings, CIs, fine-structure interactions, and coriolis forces among the A 0 and A 00 surfaces can then be omitted. The lowest vibrational and rotational states are uniquely labeled by the three vibrational quantum numbers v s , v b , and v OH , total trimer orbital angular momentum J, and parity p. Here, the vibrational quantum numbers correspond to the Sr-O stretch, the Sr-O-H bend, and the O-H stretch, respectively. We use the abbreviated notation ðv s ; v b ; v OH Þ and always have v OH ¼ 0 for frozen separation r. Our vibrational energies are discussed and compared to results from existing literature in Methods. Figure 6a examines our diagonal FCFs of stretching modes ðv 0

Evaluation of Frank-Condon factors.
For v s 4 the FCFs decrease linearly with increasing v s . In this regime, stretching modes with different bending quantum numbers are well separated in energy and FCFs based on the harmonic approximation of the potentials are in reasonable agreement with our more precise evaluation. The FCFs are largest for the 1 2 A 0 ðX 2 Σ þ Þ to 3 2 A 0 ðB 2 Σ þ Þ transition.
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 the expected angular extents of vibrational wavefunctions relative to the θ ¼ 180 equilibrium linear geometry as functions of v s . For v b ¼ 0 states the vibrational wavefunctions along θ are nodeless. We observe that the 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 0 ðX 2 Σ þ Þ and 2 2 A 0 ðA 2 ΠÞ states remains close but differ significantly from those of the 3 2 A 0 ð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 1 2 A 0 ðX 2 Σ þ Þ to 2 2 A 0 ðA 2 ΠÞ transitions remain diagonal. Figure 7 compares vibrational wavefunctions as functions of stretching-direction R for linear geometry θ ¼ 180 for the 1 2 A 0 ðX 2 Σ þ Þ to 2 2 A 0 ðA 2 ΠÞ and 1 2 A 0 ðX 2 Σ þ Þ to 3 2 A 0 ðB 2 Σ þ Þ 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 0 ðX 2 Σ þ Þ and 3 2 A 0 ð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. Table 2 presents our diagonal and off-diagonal FCF matrix elements. For simplicity FCFs smaller than 10 À3 have been omitted. Overall the 1 2 A 0 ðX 2 Σ þ Þ to 3 2 A 0 ðB 2 Σ þ Þ vibronic transitions have higher diagonal FCFs and smaller off-diagonal FCFs than those for the 1 2 A 0 ðX 2 Σ þ Þ to 2 2 A 0 ðA 2 ΠÞ transition. 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 0 ðX 2 Σ þ Þ to 2 2 A 0 ðA 2 ΠÞ and 1 2 A 0 ðX 2 Σ þ Þ to 3 2 A 0 ð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- Table 1 Orbital composition analysis for the ground X 2 Σ þ and lowest lying electronic A 2 Π and B 2 Σ þ states of the SrOH molecule. For each of the three molecular states columns correspond to amplitudes of the s, p, and d valence orbitals of the strontium cation (Sr-s, Sr-p, Sr-d), the s valence orbital of oxygen (O-s), and the s orbital of hydrogen (H-s).

State
Sr-s Sr-p 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 For all excited diabatic potentials, 1 2 A 00 , 2 2 A 0 , 3 2 A 0 , 2 2 A 00 , and 4 2 A 0 , 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 0 symmetry and two roots of the A 00 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 0 and A 00 PESs and between the 2 2 A 0 and 3 2 A 0 PESs.
The CC based methods used in construction of the 2 A 0 and 2 A 00 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 0 and two 2 A 00 twodimensional 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 V mr ðR; θÞ ¼ 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 V lr ðR; θÞ ¼ 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 VðR; θÞ ¼ ½1 À f ðRÞ V mr ðR; θÞ þ f ðRÞ V lr ðR; θÞ ð 8Þ and f ðRÞ ¼ ð1 þ tanh α R À R sw ð Þ ½ Þ =2: The parameter values for describing U l ðRÞ, C nl , α, and R sw are available upon request.
Bound state energies. For each diabatic PESs the corresponding vibrational Hamiltonian includes the kinetic-energy operators for Jacobi coordinates R and r with masses μ ¼ m Sr m OH =ðm Sr þ m OH Þ and m OH , respectively. Here, m Sr is the mass of the Sr atom and m OH that of the OH molecule. The vibrational Hamiltonian commutes with the trimer orbital-angular-momentum operator J = j þ l as well as the parity operator p, the symmetry under spatial inversion of all electrons and nuclei around the center-of-mass position of the trimer. Here, j and l are the orbital-angular-momentum operators of ground-state OH and of Sr relative to the center of mass of OH, respectively.
We diagonalize the vibrational Hamiltonian by first creating a one-dimensional basis of bound states of the radial Hamiltonian À_ 2 =ð2μÞd 2 =dR 2 þ VðRÞ using the distributed Gaussian basis approach and where VðRÞ is the interpolated PES computed at θ ¼ 180 and r ¼ 1:832a 0 (_ is the reduced Planck constant.) Then the vibrational Hamiltonian is diagonalized in the product basis of radial bound states and superpositions of products of spherical harmonic functions for the orientation of R and r, such that the basis states are eigenstates of J 2 , J z , and parity p.
In Table 3 we compare some of our calculated energies with total angular momentum J ¼ 0 to other semi-empirical results as well as experimental data. Our results for the Sr-O stretch and Sr-O-H bend energies agree well with experimental data, even though the OH bond is not allowed to stretch and diabatic couplings are omitted. The energy difference between the ðv s ; v b ; 0Þ ¼ ð1; 0; 0Þ and ð0; 0; 0Þ states, the Sr-O stretching mode, are reproduced to better than 1% for the 1 2 A 0 ðX 2 Σ þ Þ and 1 2 A 00 ðA 2 ΠÞ states. The discrepancy increases to 2.5% for the 3 2 A 0 ð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 Table 2 Franck-Condon factors between pairs of ðv s ; v b ; 0Þ rovibrational SrOH states of J ¼ 02 2 A 0 ðA 2 ΠÞ, J ¼ 03 2 A 0 ðB 2 Σ þ Þ, and J ¼ 11 2 A 0 ðX 2 Σ þ Þ. In the table the three states are denoted by A, B and X, respectively. Here parity p = +1 or −1 corresponds to whether vb is even or odd Table 3 Comparison of energies for three ðv s ; v b ; 0Þ stretching and bending levels with total angular momentum J ¼ 0 and parity þ1 for the 1 2 A 0 ðX 2 Σ þ Þ, 1 2 A 00 ðA 2 ΠÞ, and 3 2 A 0 ðB 2 Σ þ Þ states of the 88 Sr- 16 46 536 − − Exp. 45 − 401 771 Energies are relative to that of the (0, 0, 0) mode and in units of cm −1 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 Figs. 1-3.

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.