key: cord-0986206-xb83gwb5
authors: Clark, Victoria H.J.; Yurchenko, Sergei N.
title: Modelling the non-local thermodynamic equilibrium spectra of silylene (SiH2)
date: 2021-04-28
journal: Physical chemistry chemical physics : PCCP
DOI: 10.1039/d1cp00839k
sha: 91b9efe4ebd914526b817861d8c1a039025ab79b
doc_id: 986206
cord_uid: xb83gwb5

This paper sets out a robust methodology for modelling spectra of polyatomic molecules produced in reactive or dissociative environments, with vibrational populations outside local thermal equilibrium (LTE). The methodology is based on accurate, extensive ro-vibrational line lists containing transitions with high vibrational excitations and relies on the detailed ro-vibrational assignments. The developed methodology is applied to model non-LTE IR and visible spectra of silylene (SiH$_2$) produced in a decomposition of disilane (Si$_2$H$_6$), a reaction of technological importance. Two approaches for non-LTE vibrational populations of the product SiH$_2$ are introduced: a simplistic 1D approach based on the Harmonic approximation and a full 3D model incorporating accurate vibrational wavefunctions of SiH$_2$ computed variationally with the TROVE (Theoretical ROVibrational Energy) program. We show how their non-LTE spectral signatures can be used to trace different reaction channels of molecular dissociations.

Normally, molecules are assumed to be in local thermal equilibrium pertaining to a given temperature with the internal degrees of freedom (electronic-rotation-vibration) characterized by the Boltzmann distribution. However many different physical chemical, experimental and technological processes produce molecules that do not satisfy the Boltzmann law and as consequence have unusual, non-local thermal equilibrium (non-LTE) spectroscopic signatures. Molecules produced in reactions do not necessarily obey the Boltzmann thermal equilibrium, at least if the reaction time is shorter than the collision time. Instead, their internal degrees of freedom are populated based on the reaction paths rather than on the temperature of the surrounding environment. These out-of-LTE (i.e. non-LTE) populations encode information about the structural reaction dynamics and can manifest in the molecular spectra. The field of non-LTE spectroscopy has great potential to study these processes as the properties of the molecules producing the non-LTE spectroscopic signatures can shed the light on the dynamics of chemical reactions. [1] [2] [3] [4] [5] [6] The socalled transition state (TS) spectroscopy is a technique already widely used that employs the high-resolution non-LTE spectra of products to observe reaction processes that are hidden for the conventional spectroscopic methods. 7, 8 The novel high resolution non-LTE spectroscopic techniques allow decoupling of the vibrational and rotational degrees of freedom of molecules and thus control their vibrational and rotational populations, e.g. with rich vibrational and simplified rotational structures. 6, [9] [10] [11] [12] [13] The modern day study of non-LTE spectroscopy can be traced as far back as the 1930s to the original papers of Milne 14 , and the many key papers from the decades following. [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] The non-LTE spectroscopic effects play important role in high-resolution applications and there exist a number of accurate non-LTE spectroscopic and radiative transfer codes, see van der Tak et al. 25 , Funke et al. 26 , Pannier and Laux 27 and references therein. As such, non-LTE spectra are often vital for the modelling of astrophysical problems, including planetary atmospheric properties 28 , stellar atmospheres of solar system and exoplanets [29] [30] [31] [32] [33] [34] [35] [36] and the ISM. 25, 37, 38 The non-local thermodynamic effects of the spectra of molecules has been of interest to chemists and astronomers alike for many years. A notable example of this is the 2011 work by Ferus et al. 2 who studied the isomers of HCN within acetonitrile, formamide, and BrCN discharge. Using the features for both the HCN and HNC molecules from these spectra, Ferus et al. 2 were able to calculate the ratios of molecules within the reactions and also the reaction path taken by the HNC molecule during the isomerization. In this work we explore this idea to study non-LTE spectral signatures of silylene (SiH 2 ) produced from disilane (Si 2 H 6 ).

Reaction properties of silylene, silane and disilane such as the rate constants for the formation, destruction and chemical pathways are important for plasma physics aspects such as silicon deposition. [39] [40] [41] The ease of hydrogen transfer and high barriers in the saturated silicon system, leading to the ready formation of three-centre interactions and consequently the isomerisation reactions of Si 2 H 6 , are just as important to study as the elementary reactions. 42 Silane containing reactions are also of importance for astrophysics, with the presence of SiH 4 in IRC +10216 discussed by Goldhaber and Betz 43 , Kaiser and Osamura 44 .

The complexity of the silane containing systems has been discussed elsewhere 42, 45 , with quantitative calculations proving particularly difficult. There were a number of ab initio studies of the structural properties of Si 2 H 6 e.g., [46] [47] [48] [49] [50] [51] [52] [53] as well as of the formation and reactions involving this molecule. [54] [55] [56] [57] [58] Agrawal et al. 53 produced a global ab initio potential energy surface of disilane and used it to investigate the dissociation dynamics with classical trajectories. Márquez et al. 50 reported a force field for Si 2 H 6 . The main isomer of Si 2 H 6 has a staggered, ethane-like, structure (see Fig. 1 ) with a low barrier (1.2 kcal/mol) to the eclipsed conformation 59 . Si 2 H 6 has been shown to have a local minimum as an inverted stable structure with one of the Si-H 3 'umbrellas' pointing to the center as well as a transition state with a similar inverted configuration. 46, 60 These structures are nominally asymmetric (C s symmetry) but essentially acquiring the C 3v symmetry. Thermal decomposition of Si 2 H 6 has been extensively studied, both theoretically (mostly using RRKM, Rice-Ramsperger-Kassel-Marcus) and experimentally, with the reaction Si 2 H 6 → SiH 2 + SiH 4 as the most common 53 67 ). However, it is also possible for the disilane molecule to dissociate homolytically, as was originally thought to be the main pathway owing to disilanes similarities with ethane, and form 2SiH 3 , 53, 56, 68 or to undergo dehydrogenation to H 2 Si-SiH 2 , H 3 Si-SiH, or H 3 Si-H 2 . 53,54,57 Disilane can even undergo double dehydrogenation to form Si 2 H 2 , however to our knowledge this has only been reported as the main product when undergoing photolysis at 193 nm. 69 Yoshida et al. 57 notes that the transition state for the Si 2 H 6 → SiH 2 + SiH 4 is 8.48 kcal mol −1 lower than the transition state for Si 2 H 6 → H 3 SiSiH + H 2 , at 43.38 kcal mol −1 compared to 51.86 kcal mol −1 .

The spectroscopy of SiH 2 has been used to monitor the SiH 2 + SiH 4 → Si 2 H 6 and Si 2 H 6 → SiH 2 + SiH 4 reactions (see Fig. 1 ) and measure the corresponding rate constants and Arrhenius parameters by spectroscopically tracking electronic (Ã-X) transitions of SiH 2 39,46,65,70 . In these studies, the reconstructions of the amount of SiH 2 relied on the assumption of the Boltzmann thermal distribution when estimating the population of the lower state. No account of the possible non-LTE population of SiH 2 molecules after dissociation was made, which could potentially hamper the count of the SiH 2 molecules and effect the reaction rates estimated. A similar experimental technique was used in Hertl and Jolly 71 to monitor SiH 2 in SiH 4 plasma. It is the second, dissociation, part of the reaction shown in Fig. 1 (Si 2 H 6 → SiH 2 + SiH 4 ) we study in this work. More specifically, we show that (i) the (vibrational) populations of the molecules produced in reactions can be very different from the Boltzmann distribution and is important to take into account when interpreting spectroscopic measurements. That (ii) spectral shapes of the dissociated SiH 2 can bear strong non-LTE character, very different from the LTE spectrum of an LTE SiH 2 sample making it possible to distinguish between different reaction stages and even between different dissociation channels the silylene molecules it is produced from. In this work non-LTE spectra of SiH 2 under conditions similar to dissociation processes expected in these experiments are modelled. Recently we have computed an accurate ro-vibrational line list for SiH 2 , named CATS. 72 It covers a large range of rotational and vibrational excitation, capable of modelling very hot spectra of this molecule (up to T = 2000 K) as part of the ExoMol database. 73 The CATS line list was produced using the program TROVE, 74, 75 which solves the nuclear motion Schrödinger equation variationally. The ro-vibrational energies and corresponding wavefunctions were computed using an accurate, empirically refined potential energy surface (PES) of silane and a high-level ab initio dipole moment surface (DMS). The ro-vibrational probabilities (in the form of Einstein A coefficients) were computed using a high level ab initio dipole moment surface.

The study by Clark et al. 72 forms the basis for the present work, where we utilize the CATS line list, wavefunctions, purposebuilt numerical basis set, and the CATS computational TROVE setup to model non-LTE spectroscopic properties of SiH 2 produced from dissociation of Si 2 H 6 through different reaction channels. Using a simplified 1D Harmonic oscillator wavefunctions (see Pastorek et al. 6 ) and more sophisticated 3D vibrational CATS wavefunctions from accurate variational calculations, the non-LTE ro-vibrational populations of SiH 2 are generated and used to produce non-LTE spectroscopic spectra of different dissociation channels of disilane. To this end we investigate reaction topology connecting the global minimum of Si 2 H 6 with the closest saddle points and local minima as well as the corresponding structural properties using a high level ab initio theory cc-pVTZ-F12/CCSD(T)-f12b 76,77 employing the program MOL-PRO2015. 78 Theoretically, the non-LTE properties of dissociating molecules were studied by Band and Freed 79 . In the present work we use general approach of Berry 80 and Band and Freed 79 , which assume no significant structural changes between the reactant and product nuclear configuration, along with the slow vibrational relaxation of the product 13 to investigate non-LTE spectroscopic signatures of SiH 2 produced from dissociation of Si 2 H 6 . In this paper we specifically consider situations where the vibrational relaxations are not achieved during the time of the experiment, so that the molecules still hold the memory of the structure during the reaction or dissociation. The rotation relaxation time however is much shorter and the rotational degrees of freedom can be usually assumed to satisfy the Bolzmann equilibrium. 13 We also investigate possible non-LTE impact on the electronic A (0,2,0)-X (0,0,0) spectrum of SiH 2 . This is a favorite spectroscopic system for the detection of SiH 2 due to the large Franck-Condon factor and the availability of suitable laser. 39, 70, [81] [82] [83] [84] [85] [86] [87] [88] [89] [90] [91] The non-LTE absorption spectra of SiH 2 are simulated using the (non-LTE) ExoCross program, 92 where a new feature of non-Boltzmann populations was added. ExoCross has been previously used to model spectra of molecules in environments that can be characterized using two temperatures, vibrational and rotational. [93] [94] [95] [96] The paper is structured as follows. In Sec. 2 we describe the calculations of potential energy surfaces for the disilane and silylene structures. The theory used in this paper is described in Sec. 3. In Sec. 4.1 we calculate the 1D harmonic wavefunction population and use them to produce non-LTE spectra of SiH 2 corresponding to different dissociation routes. In Sec. 4.2 we calculate the populations and subsequent non-LTE using the full 3D wavefunctions and describe the new TROVE methodology. A non-LTE electronic A (0,2,0)-X (0,0,0) spectrum of SiH 2 is presented in Sec. 4.4. Conclusions are offered in Sec. 5.

2 Geometry optimisation and reaction topology of Si 2 H 6

In order to better understand the reaction process of breaking Si 2 H 6 , the topology of Si 2 H 6 has been investigated by performing a structural analysis of Si 2 H 6 using a high level ab initio theory. This includes finding the global minimum (GM), local minima (LM), transition states (TS), reaction barriers as well as reaction paths, as detailed below. A reaction slice through the global PES of Si 2 H 6 helps to indicate how likely local minima or transitions states were to be formed based on corresponding topology. These properties of Si 2 H 6 were obtained using the geometry optimization and reaction path finder implemented in MOLPRO2015 78 using the explicitly correlated coupled cluster method CCSD(T)-F12b 76,77 with the F12-optimized correlation consistent basis set, VTZ-F12 97 in the frozen core approximation. The calculations employed the diagonal fixed amplitude ansatz 3C(FIX) 98 and a Slater geminal exponent value of β = 1.0 a −1 0 . 99 The auxiliary basis sets were chosen to be the resolution of the identity OptRI 100 basis and the aug-cc-pV5Z/JKFIT 101 and cc-awCV5Z/MP2FIT 102 basis sets for density fitting. In the following this level of theory will be referenced to as VTZ/CCSD(T)-F12b.

We shall refer to different disilane isomers as dGM, dLM and dTS to distinguish them from the SiH 2 fragments GM, LM and TS as discussed below.

The global minimum of Si 2 H 6 (dGM) has a symmetrical, staggered D 3d structure. The closest local minimum (dLM) has an inverted, C 3v structure. The lowest transition state, which will be referred to as dTS (TS1 in Becerra et al. 46 and TS2 in Tonokura et al. 103 ), has also inverted structure, just a few kJ/mol above dLM. 46 These structures together with the corresponding opti-mized parameters for three geometries most relevant to our work are collected in Table 1 . Our structural parameters compare well with that from the literature. The structure of the dGM has also been determined spectroscopically, with the equilibrium bond lengths r Si−H = 1.492 Å, r Si−Si = 2.331 Å and bond angles β ∠HSiSi = 110.3 • and α ∠HSiSi = 108.6 • . 104 The structure of the deuterated isotopologue Si 2 H 5 D was reported as r Si−H = 1.4874(17) Å, r Si−Si = 2.3317(15) Å, β ∠HSiSi = 110.66 (16) • . 105 The reaction path connecting the disilane isomers dGM, dLM and dTS is shown in Fig. 2 . A zoom of the dLM side is shown as inset. The energies of the global and local minima are 178 kJ/mol and 4.1 kJ/mol below the transition state, respectively. The energy and geometry information for dLM, dGM, dTS are collected in Table 2 . The results calculated compare well with both the results of Becerra et al. 46 and Sakai and Nakamura 60 , albeit both vary for the dGM structure by 36 kJ mol −1 The dLM isomer of Si 2 H 6 has a shallow potential with a very low barrier to dTS of 344 cm -1 . It can be also recovered using lower levels of theory, for example, using MP2/6-311G(d, p) 46 and even with the UFF force fields implemented in Avogadro 1.2.0 106 via the steepest descent method and 4 steps per update.

With the aim to give more quantitative information on the structural and dynamical properties of five SiH 2 fragments from Si 2 H 6 , the fragments are described as follows. GM is an SiH 2 fragment from the global minimum structure (dGM); LM-L is an SiH 2 fragment from the left hand side (LHS) of the local minimum structure (dLM) with the Si-H 3 umbrella group pointing outside; LM-R is an SiH 2 fragment from the right hand side (RHS) of the dLM structure with Si-H 3 pointing inside; TS-L is an SiH 2 fragment from the LHS of the transition state structure (dTS), Si-H 3 umbrella group points outside; and TS-R is an SiH 2 fragment from the RHS of the dTS structure, Si-H 3 umbrella group pointing inside. The structural parameters and structures are shown in Table  1 and the energies are shown in Table 2 . The harmonic frequencies for the disilane molecules were computed with MOLPRO using the TVZ/CCSD(T)-F12b level of theory. The results from the frequency analysis are shown in Table 3 along with a comparison with literature (both experimental, if available, or theoretical data is shown). The symmetry of dGM is D 3d whereas the symmetry for the dLM and dTS are both C 3v . The columns titled "Theory" were calculated in this work, and the degenerate states have been removed, with an average calculated if there were any differences owing to computational errors associated with the lower symmetry used by MOLPRO. All frequencies of dLM are positive thus confirming that it is a minimum with a stable structure. The 'negative' (or imaginary) harmonic frequency of dTS is -551.6 cm -1 .

We now consider the decomposition reaction Si 2 H 6 → SiH 2 + SiH 4 shown in Fig. 1 and model the vibrational populations of the product SiH 2 assuming that the corresponding relaxation time to LTE is much longer than the time of the spectroscopic experiment. We aim at simulating non-LTE IR spectra of SiH 2 using these populations to demonstrate their impact on the spectral shape of dissociated species.

In our description of the non-LTE vibrational population of the dissociated molecule we follow the polyatomic Franck-Condon type approximation by Band and Freed 79 and Berry 80 based on the structural differences between reactant and product assuming no significant change in nuclear configuration of the molecule. In order to connect the product (gas phase SiH 2 ) to an initial structure of Si 2 H 6 through the dissociation process, we assume that the dissociation happens instantaneously, i.e. the initial configuration of the product SiH 2 corresponds to the structural parameters (bond lengths Si-H and inter-bond angles ∠HSiH) of SiH 2 as a fragment of Si 2 H 6 , for which the parameters collected in Table 1 are used. For example, for the dissociation from the dGM structure, the initial configuration of the gas phase SiH 2 is assumed to be r Si−H = 1.482 Å, α ∠HSiH = 108.6 • . Naturally, this is a very deformed geometry comparing to the equilibrium structure of the gas phase SiH 2 , r e = 1.5137 ±0.003 Å and α e = 92.04 ±0.05 • . 107 After being dissociated, in relaxing to be a free molecule, the fragment SiH 2 has added vibrational energy and is hence in a situation when it's vibrational populations do not match the LTE distribution for corresponding temperature of the surroundings, at least for the vibrational degrees of freedom.

Experience shows that the rotational equillibration time is usually very short and we can thus safely assume the LTE conditions for the rotational degrees of freedom with the rotational temper- ature the same as the temperature of the surroundings (see also e.g. Dudás et al. 13 ). The vibrational population however is not in the thermal equilibrium and therefore no sensible vibrational temperature could be associated with the corresponding population.

The non-LTE population N J,k, (T ) of a ro-vibrational state |J, k, is then given by: 6

where N vib is a non-LTE vibrational population, T rot is the rotational temperature, J is the total angular momentum quantum number; k is a generic rotational quantum number, e.g. the projection of the total angular momentum on the molecular z axis; is a generic vibrational quantum number/label, e.g. a combination ( 1 , 2 , 3 ) to describe vibrational states of a triatomic molecule; g ns i is the nuclear spin degeneracy; T is the temperature, c 2 is the second radiation constant. In Eq. (1),Ẽ rot J,k is the rotational part of the ro-vibrational energy term value approximated asẼ

where hc ·Ẽ vib is the vibrational (J = 0) energy ('band center') and hc ·Ẽ J,k, is the total ro-vibrational energy. The non-LTE partition function Q nLTE in Eq. (1) is given by

Our aim is to calculate the vibrational populations N v ≡ N vib v of SiH 2 as produced by instantaneous (vertical) dissociation from three structures, dGM, dLM and dTS. In case of dLM and dTS, the SiH 2 fragment can originate from any of the two different sides of their inverted structures, which should be taken into account. We therefore have to consider five different fragments, as shown in Table 4 .

Let us assume that the Si 2 H 6 is LTE and hence is in its ground vibrational state at the moment of dissociation, while the fragment SiH 2 can end up in any vibrationally excited state ≡ ( 1 , 2 , 3 ) with some transition probability giving rise to the vibrational population N vib . On top of that we also assume a full separation of the stretching Si-H and bending ∠HSiH modes inside Si 2 H 6 in its ground vibrational state. Possible consequences of deviation from these approximations are discussed below.

Under the assumptions made we define the vibrational population of the SiH 2 fragment as a Franck-Condon factor for a vertical transition from the ground vibrational state of disilane Si 2 H 6 to SiH 2 + SiH 4 with the gas phase (g.ph.) SiH 2 transferred to some vibrational state | (g.ph.) = | 1 , 2 , 3 . In the approximation of the full separation of the fragment SiH 2 from the rest of Si 2 H 6 , the population of SiH 2 can be represented as an overlap between the ground state wavefunction | = 0(fragment) of a fragment SiH 2 and that of the corresponding vibrational state | = (g.ph.) of the gas phase SiH 2 as given by

Here the spacial wavefunction of the SiH 2 fragment (i.e. a combination of two adjacent Si-H bonds in disilane with an angle between them forming the dissociating SiH 2 ) is projected on vibrational eigenfunctions of the gas phase SiH 2 to give the corresponding populations of SiH 2 . The highest populated energy level will have the largest overlap between these wavefunctions.

The calculated temperature dependent populations N J,k, (T ) in Eq. (1) can be then combined with a molecular line list for SiH 2 to simulate absorption or emission spectra of this molecule under the non-LTE conditions as defined by N J,k, (T ) in Eq. (1). Here we use the ExoMol line list CATS by Clark et al. 72 as provided by Exo-Mol (www.exomol.com). Technically this is done by incorporating the non-LTE vibrational densities N into the ExoMol States file as described in Section 4 (the ExoMol file formats are discussed extensively elsewhere 73 ). A non-LTE spectrum of SiH 2 for given T and P is then calculated using CATS' Einstein-A coefficients with the ExoCross program, 92 where a new non-LTE option has been implemented as part of this work. The rotational populations are assumed to be in LTE according with Eq. (1).

Two approaches were used for the calculation of the population densities of the fragment SiH 2 . One approach -named the decoupled 1D approach -is where the 3D wavefunctions of the fragment as well as of free SiH 2 are represented by products of 1D parts with the harmonic oscillators as wavefunctions. This simplified model is mainly used to illustrate the idea of our non-LTE treatment. The second, more accurate approach -named the 3D approach -is based on the full 3D vibrational wavefunctions computed using the variational program TROVE. 74 Both approaches are presented in the following in order to assess and compare the accuracy achieved. 16 . The dTS 16 can be found at 960 cm -1 . The other two 925 cm -1 line were assigned to 9 and 10 (E). 

A vibrational state | = | 1 , 2 , 3 of SiH 2 is characterized by the three (normal mode) quantum numbers 1 , 2 and 3 corresponding to the two stretching modes ( 1 and 3 ) and one bending mode ( 2 ) of the SiH 2 molecule. The 1D approach considers the stretching Si-H 1 , Si-H 2 and bending ∠HSiH modes, both of the molecular fragment and gas phase SiH 2 molecules, as fully independent and described by one dimensional (1D) wavefunctions under the harmonic approximation, as given by:

Here x is a dimensionless coordinate describing either the stretching r = r Si−H or bending α = α ∠H−Si−H coordinate as follows:

with

and ω s = ω Si−H , ω b = ω ∠HSiH . In Eq. (5) H (x) is a Hermite polynomial and C is the corresponding normalization constant.

The constants a s and a b correspond to inverse masses of the vibrational part of a free three-atomic molecule expressed in terms of the internal coordinates r 1 , r 2 , α (see, e.g. Sutcliffe and Tennyson 108 , Yurchenko et al. 109 ).

A 1D population for the corresponding mode of the gas phase SiH 2 molecule is given by Eq. (4) with | = Ψ (x). The different disilane fragments have different structural parameters r e , α e , ω s and ω b , see Tables 1 and 3 , and thus lead to different ground state vibrational 1D wavefunctions |0(fragment) (stretching or bending) and hence result in different vibrational populations N of the gas phase SiH 2 according with Eq. (4).

For each of the three modes (two stretching and one bending), 1D wavefunctions of the gas phase SiH 2 for 30 vibrational states from 0 up to 29 were calculated. These 1D wavefunctions were then numerically integrated with the corresponding ground state 1D wavefunctions | = 0 of the fragment in question.

The total vibrational population N 1 , 2 , 3 in this approximation is then given by a product

where N str 1 and N str 3 are obtained using the stretching harmonic oscillators wavefunctions | 1 and | 3 , while N bnd 2 is obtained using the corresponding bending harmonic oscillator wavefunction | 2 . The independent treatment of the two stretching populations is partly justified by the local mode character of the vibrational degrees of freedom of SiH 2 due to the 90 • bond angle (see, e.g. Jensen 110 and Clark et al. 72 ). The asymmetric vibrational modes of SiH 2 (B 2 in C 2v ) are non populated in this 1D approximation. This is because for the parallel nature of the Franck-Condon transitions from the ground vibrational state of disilane, which is fully symmetric (A 1 ) and the excitation, only symmetric states of SiH 2 give rise to non-zero integrals in Eq. (4). For example, the vibrational population of the (2,1,0) state (2 stretching and 1 bend-ing quanta of A 1 ) of the GM fragment is obtained as a product of N str 1 = 0.274, N bnd 2 = 0.333 and N str 3 = 0.967, respectively, resulting in N (2,1,0) = 0.088, while the population of the B 2 -type (0,0,1) vibratitonal state is assumed to be zero. The populations N 1 , 2 , 3 are pre-calculated for each vibrational state ( 1 , 2 , 3 ) of SiH 2 and added to the CATS State file to be used in non-LTE simulations (see below, Section 4.3). All the vibrational populations computed and used as part of this work are provided in the supplementary material.

Examples of overlapping bending mode Harmonic wavefunctions used in calculations of populations for = 0, 1 and 2 of the five fragments are shown in Figures 3 and 4 for the bending and stretching modes respectively. In all cases the black curve represents the ground state wavefunction for the non-LTE fragment, and the blue, green and red line show the = 0, = 1 and = 2 LTE wavefunctions of gas phase SiH 2 . 72 The corresponding 1D populations as an integral of the overlaps between the LTE and non-LTE wavefunctions are plotted in Figure 5 for the bending and stretching modes of the five fragments. The stretching populations exhibit a typical Boltzmann-like distribution with the ground vibrational state = 0 as the mostly populated in all five cases. This is expected because their equilibrium bond lengths are rather similar to that of the gas phase SiH 2 . In case of the bending populations, only LM-R and TS-R have = 0 to be with the highest populations, while for LM-L, TS-L and GM the N distributions exhibit strong non-LTE character with = 2 to be almost as populated as the ground vibrational state = 0.

The shapes and positions of the curves in Figures 3 and 4 match with the parameters from Table 1 . The larger ω used for the bending modes manifests itself as wider curves, while the black curves are all centred around the equilibrium bond angles and length listed in Table 1 .

The similarity between the populations of the LM-L, TS-L and GM fragments is expected owing to their similar structural parameters. It is interesting to see the most populated vibrational levels of LM-R and TS-R are always lower than the vibrational levels of LM-L, TS-L and GM. Figures 6 and 7 show the 1000 cm -1 and 2000 cm -1 (10 µm and 5 µm) bands for the SiH 2 absorption spectrum, respectively, simulated using the non-LTE densities from Fig. 5 for all five cases considered and compared to the LTE scenario assuming the (rotational) temperature of T = 296 K and using the CATS line list. The strongest bands are indicated using different colours. Figure 6 focuses on the 1000 cm -1 band. Most of the non-LTE spectra contain bending hot bands (020)-(010) and (020)-(010), which are stronger than the fundamental band (000)-(000). It can be seen that the P and R branches of the non-LTE spectra are shifted to lower wavenumbers in the GM, LM-L and TS-L spectra. In the TS-R and LM-R spectra the bands are not shifted, with the LM-R spectrum having only the fundamental (010)-(000) band visible.

The plots in Fig. 7 show the 2000 cm -1 band in the region of the polyad (100)/(020)/(001) for the five fragments, with the strongest fundamental band (001)-(000). The non-LTE intensities of the hot bands (011)-(010), (200)-(020) are found to be comparable to the intensities of the (001)-(000) band. The Q branch is clearly shifted for the GM, TS-L and LM-L molecules. The band is less shifted for the TS-R and LM-R fragments, but owing to the increased similarity between the fragment and molecular structures with TS-R and LM-R this is to be expected. Only the main polyad system (100)/(020)/(001) is visible for the TS-R spectrum (indicated as (001)-(000) in Fig. 7) .

With the equilibrium structures of the TS-R and LM-R fragments being similar to the equilibrium structure of SiH 2 , their non-LTE spectra are expected to be a similar spectrum to LTE. Indeed, for the 1D harmonic approach their P, Q and R branches maintain the expected LTE intensities for both the 1000 cm -1 and 2000 cm -1 bands.

In a full 3D approach, the ground state wavefunction |0, 0, 0(fragment) represents an SiH 2 fragment of an 18D ground state vibrational wavefunction of Si 2 H 6 :

φ Si 2 H 6 = |0, 0, 0(SiH 2 ) |0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0(SiH 4 ) .

Here we assume the approximation that the corresponding three modes (Si-H 1 , Si-H 2 and ∠H 1 SiH 2 ) are independent from the rest of the molecule so that all other modes, not relevant for the gas phase SiH 2 , can be eliminated (integrated out), including the reaction coordinate and vibrational modes of SiH 4 . This is in line with the assumptions used previously by Band and Freed 79 and Berry 80 . Apart from this approximation we will treat the SiH 2 fragment as accurate as possible. The corresponding wavefunction |0, 0, 0(fragment) is obtained by solving a 3D vibrational Schrödinger equation for these three degrees of freedom with a realistic PES obtained using a high level of ab initio theory (the same as above, VTZ/CCSD(T)-F12b with MOLPRO). The vibrational populations are then modelled using the Franck-Condon integrals as follows:

with the sum of all populations over all states totalling 1. In this equation, | 1 , 2 , 3 (g.ph.) is an accurate vibrational wavefunction of a gas phase SiH 2 molecule, obtained by solving the vibrational Schrödinger equation with an accurate PES. We use the TROVE variational program and the refined PES of SiH 2 by Clark et al. 72 to generate | 1 , 2 , 3 (g.ph.) for all vibrational excitations required. For the details on the TROVE calculations see below and also Clark et al. 72 In order to simplify the 3D integration in Eq. (11), the variational wavefunction |0, 0, 0(fragment) is obtained using the same vibrational basis set as the variational solution of the gas phase SiH 2 . By taking advantage of the compatibility of the orthogonality of the basis sets, the Franck-Condon factors are then given by as a sum of products 0, 0, 0(fragment)| 1 , 2 , 3 (g.ph. ) = The 1D bending mode harmonic wavefunctions for the five SiH 2 fragments described in Table 1 as a function of bond angle in degrees. Ground state fragment (black) compared with the LTE | = 0 (blue), | = 1 (green), | = 2 (red)

of the corresponding eigen-coefficients C i 1 ,i 2 ,i 3 (fragment) and C i 1 ,i 2 ,i 3 (g.ph.), obtained variationally in independent calculations using a new implementation in TROVE.

TROVE uses optimized non-standard vibrational basis sets, generated numerically by solving 1D Schrödinger equations for realistic 1D potentials. 74 This procedure allows producing compact basis functions optimized for a specific problem. In our case, the PESs of the corresponding five fragments and of the gas phase SiH 2 are different and therefore the generated basis sets would be different and even not orthogonal. We therefore implemented a feature in TROVE allowing to read and use externally generated basis functions. Of course all relevant calculation setups must be compatible, including the numerical grids used for the stretching and bending modes and their sizes. Using foreign basis sets certainly degrades their quality. However, since we are only interested in fragments' ground state wavefunctions, this degradation can be mitigated by including enough basis functions. Our typical 1D basis sets contain 12 -24 functions (see details below), which should be more than enough to obtain a converged ground state solution even with non-optimized basis sets.

For our new 3D populations corresponding to dissociations from the fragments, five PESs were generated as follows. We assume that PESs of a dissociating Si 2 H 6 molecule can be approximated as a sum of two independent fragments:

where the individual stretching and bending modes of SiH 2 fragments are fully separable:

The stretching part of the potential is given by a Morse-like expansion

while the bending part is a Taylor-type expansion in terms of the displacement from the corresponding equilibrium value:

The expansion constants f i and g i representing fragments' potential energies V SiH 2 (r 1 , r 2 , α) were obtained by fitting Eqs. (14, 15) to the ab initio data computed as 1D slices on the global surfaces for the five fragments from Si 2 H 6 (dGM, dLM and dTS) using VTZ/CCSD(T)-F12b consisting of 24 bending and 34 stretching geometries, distributed around the corresponding equilibria. Figure 8 illustrates the ab initio PESs of different fragments as 1D cuts for the stretching and bending modes compared to the corresponding cuts of the gas phase SiH 2 molecule. The bending cuts have especially different shapes with shifts to larger equilibrium angles and much steeper PESs. The differences in the stretching cuts are less pronounced. These features are important for the non-LTE behaviour of the corresponding excited states populations, with the bending degree of freedom to have stronger non-LTE character than stretching.

The vibrational wavefunctions of SiH 2 were computed using the variational nuclear motion program TROVE with the same setup as in Clark et al. 72 . Details of the TROVE methodology are discussed extensively elsewhere. 74, [111] [112] [113] Here, we give a brief outline of the main calculation steps. The TROVE kinetic energy operator is Taylor expanded up to sixth order around the SiH 2 equi- Fig. 4 The 1D stretching mode harmonic wavefunctions for the five SiH 2 fragments described in Table 1 as a function of bond length in Angstrom. Ground state fragment (black) compared with the LTE | = 0 (blue), | = 1 (green), | = 2 (red) librium geometry in terms of linearized coordinates. 114 The primitive basis set is constructed from 1D mode numerical basis functions using the Numerov-Cooley approach 115,116 by solving three 1D Schrödinger equations, for each vibrational degree of freedom. The stretching basis functions are then improved by solving a 2D Schrödinger equation for a reduced stretching Hamiltonian. The resulting stretching eigenfunctions are contracted, classified according with the C 2v (M) symmetry group 114 using an optimized symmetrization procedure 112 and combined with the bending primitive basis functions to form our final, symmetryadapted 3D vibrational basis set. The basis set coverage is defined by the polyad number cut-off

where 1 and 3 are the stretching and 2 is the bending quantum numbers, with the maximal excitations 12, 12 and 24 respectively.

For the gas phase SiH 2 calculations we employed the empirically refined PES by Clark et al. 72 For the five SiH 2 fragments only the vibrational ground state wavefunctions |0, 0, 0 were computed using the same setup and utilizing the basis functions from the g.ph. calculations as described above, but for the fragments' ab initio PESs from Eq. (13).

The 3D vibrational populations N 1 , 2 , 3 in Eq. (11) were computed using the corresponding TROVE eigenfunctions via Eq. (12) .

The 3D results generally agree with the 1D. As an example, Fig. 9 shows the TS-L case: The bending primitive wavefunctions of the gas phase SiH 2 are compared to the = 0 bending wavefunction generated for the PES of TS-L. Figure 10 shows the corresponding non-LTE vibrational populations as a function of the corresponding energies for the same TS-L scenario of the 3D vibrational populations. These populations are very different from the Boltzmann distribution, also shown on this figure for the SiH 2 vibrational states at T = 1000 K, which exhibits an exponential decay with = 0 at its maximum. It nicely demonstrates that it would not be possible to associate a single vibrational tempera- 

A non-LTE absorption line intensity I fi (cm/molecule) can be calculated as

where A fi is the Einstein-A coefficient (s −1 ),ν fi is the transition wavenumber (cm −1 ), Q nLTE (T ) is the non-LTE partition function defined in Eq. (3). In order to simulate absorption spectra of the gas phase SiH 2 assuming a non-LTE vibrational populations N 1 , 2 , 3 , we use the line positions and Einstein A coefficients from the ExoMol CATS line list 72 employing the ExoCross program. 92 The ExoMol line lists are formatted as two files, a States file and a Transition file. It is described extensively elsewhere 73 and in this paper we shall only discuss how the ExoMol format has been adapted for use in non-LTE situations.

To adapt the CATS States file 72 for non-LTE applications an additional 'density' column N 1 , 2 , 3 was added as a final column, see an extract from the States file in Table 6 . This column contains the weightings to the transition probabilities as populations of the vibrational levels occupied by the gas phase SiH 2 . The 'density' column is specific for the calculation of line lists for non-LTE molecules is not routinely included into the ExoMol States files. This column is read by ExoCross and used to give the population weighting to each line intensity. Evaluation of the non-LTE population N J,k, (T ) given for each ro-vibrational state as in Eq. (1) is based on the knowledge of the corresponding vibrational state ( 1 , 2 , 3 ) as well as the rotational energy contributioñ E rot J,k . Therefore for this approach to work it is mandatory for all ro-vibrational states to be vibrationally assigned in order to be able to subtract the vibrational contributionẼ vib from the total energy according with Eq. (2). All vibrational quantum numbers ( 1 , 2 , 3 ) are not required, only a vibrational index indicating the vibrational state in question. In our model, all asymmetric vibrational states (B 2 ) are not populated due to the zero overlap with the ground state of the A 1 symmetry in Eq. (11) , as part of the completely vertical Franck-Condon approximation. Figure 11 shows non-LTE spectra of SiH 2 in the two main spectroscopic regions, 1000 and 2000 cm -1 for all five fragments considered (GM, TS-L, TS-R, LM-L and LM-R). The corresponding vibrational non-LTE populations were generated with the 3D TROVE approach. The rotational populations assume the Boltzmann distribution with the rotational temperature of T rot = 296 K. These non-LTE spectra are compared to the corresponding LTE spectra of SiH 2 at T = 296 K, which comprise mainly two fundamental bands, (0, 1, 0) ← (0, 0, 0) (1000 cm -1 region) and (0, 0, 1) ← (0, 0, 0) (2000 cm -1 region), with the hot bands suppressed due to the relatively low temperature. The non-LTE spectra are dominated by the hot bands (0, 2, 0) ← (0, 1, 0), (0, 3, 0) ← (0, 2, 0) (1000 cm -1 band), (1, 1, 0) ← (0, 1, 0) and (1, 1, 0) ← (0, 2, 0) (2000 cm -1 band). The centers of the hot bands are systematically red shifted compared to the fundamental band centres and serve as distinct signatures of the non-LTE effects. The Q-branch of (011) ← (020) is especially distinct compared to the Q-branch of the LTE fundamental band in the 2000 cm -1 region, with the difference of about 15 cm -1 . Only the LM-R spectra of non-LTE are very similar to the LTE spectra. This is not surprising considering that the equilibrium values of r SiH and α ∠H−Si−H of the LM-R structure are very similar to the corresponding equilibrium parameters of the gas phase SiH 2 .

The stark differences of the non-LTE spectra offer an ability for experiment to distinguish between SiH 2 molecules produced from fragmenting Si 2 H 6 , and even to indicate dissociation channels involved.

The visible electronic bandÃ-X of SiH 2 has often been used to study different reactions involving leading to silylene. 91 Here we used the program RENNER 117 to simulate a non-LTE electronic spectrum of SiH 2 with the spectroscopic model by Yurchenko et al. 118 for theà 1 B 1 -X 1 A 1 system. The model includes two empirically adjusted PESs, for the A and X states, respectively, and an ab initio (MRCI)Ã-X transition dipole moment surface (TDMS).

In the RENNER calculations, the size of the basis set originally used in Yurchenko et al. 118 was reduced in order to be able to increase the rotational excitations. The main purpose of this exercise is to show a qualitative impact of the non-LTE populations on the spectral shape of the electronic band and not so much the quality of the line positions, and therefore a smaller basis set is justified. We used 18 and 12 bending basis functions for theX andà electronic states, respectively, for every |k| block (where |k| ≤ J). TheX 1 A 1 electronic state basis set included N A =12 stretching functions of the A 1 symmetry and N B =10 stretching functions figure) and compared to the same spectra simulated using the LTE population at T = 296 K. The non-LTE populations were obtained using the 1D approach (see text). Only the strongest bands are shown. The CATS line list was used of the B 2 symmetry. For theà 1 B 1 state N A =10 and N B =8 stretching functions were used. These stretching functions were constructed from the Morse oscillator functions |n 1 |n 3 with n 1 + n 3 ≤ N stretch =12.

A rovibronic line list for theà 1 B 1 -X 1 A 1 of SiH 2 was generated covering the rotational excitations up to J max = 15 with the lower state energies (X) truncated at hc · 25000 cm −1 and the upper state energies truncated at hc · 28000 cm −1 .

For the non-LTE simulations we used the 1D vibrational population model with the structural parameters corresponding to TS-L from Table 1. Figure 12 shows a non-LTE electronic spectrum of SiH 2 in the region of the bandà (0,2,0) ←X (0,0,0), assuming the rotational temperature T rot = 500 K, compared to an LTE spectrum of T = 500 K. The non-LTE spectrum contains the hot band A (0,3,0) ←X (0,1,0) which can be used to identify the non-LTE character of the system. The rovibronic line 1 01 ← 1 10 belonging to this band used in a number of experimental studies involving SiH 2 as a reaction product 71, 89, 90 to estimate reaction rates. It is common for such studies to assume the Boltzmann equilibrium at different stages of the analysis of the measurements. For example, the partition function of SiH 2 is required to estimate the number density of SiH 2 in its lower, ground electronic state 70, 71, 119 , which is directly affected by the LTE assumption. As our calculations show, the number densities of SiH 2 as a reaction product can vary significantly under non-LTE depending on the reaction pathway and impact the experimental rates.

The focus of this paper is on the new features which have been added to TROVE to allow modelling the non-LTE populations of polyatomic molecules. We have demonstrated this capability by modelling the non-LTE line list of SiH 2 calculated with 3D wavefunctions and TROVE, and compared them to non-LTE spectra of SiH 2 modelled using a 1D harmonic approach, and the LTE line list calculated previously by ExoMol.

There are two stable isomers of disilane, a local minimum structure and a global minimum structure, with a third transition state structure also known. Non-LTE spectra of SiH 2 corresponding to dissociation of disilane from different sides of the three disilane isomer were computed. We have shown that the non-LTE spectra of SiH 2 are different in most cases. This is important as the spectrum of SiH 2 is used to monitor the quantity present in a reaction as a means to track the progress of SiH 2 + SiH 4 → Si 2 H 6 and Si 2 H 6 → SiH 2 + SiH 4 when calculating the corresponding rate constant. If Si 2 H 6 is decomposing at a rate slower than it is being formed, then tracking the quantity of SiH 2 can give a rate constant that is not reflective of the speed of reaction, and merely an indication of the equilibrium balance of the two species SiH 2 and Si 2 H 6 .

In two approaches considered, 1D and the 3D, we assume that the rotational degrees of freedom are equillibrated quickly once the dissociation from disilane occurs, hence we use the Boltzmann distribution for the rotational degrees of freedom. We also assume that the SiH 2 fragment during the instantaneous dissociation and is fully decoupled from the rest of the Si 2 H 6 molecule, i.e. can be described by a 3D wavefunction in its lowest, relaxed vibrational configuration and has the same structural parameters as the Si 2 H 6 molecule.

We have shown that the non-LTE spectra of SiH 2 can be calculated by the new TROVE methodology and existing ExoMol line list, and it compares well to the simpler 1D harmonic approximation published previously. The method could be applied to the non-LTE spectroscopy of other small molecules including SiH 4 , which has not been explored here.

We have also shown that despite the many approximations used in the 1D approximation (separability of the modes, Harmonic approximate etc.), the results compare well to the results obtained using the full 3D approach. This lends confidence in using the simplified but robust 1D approach in similar non-LTE studies, as Vibrational (normal mode) quantum numbers; Γ vib : Vibrational symmetry in C 2v (M); K a : Asymmetric top quantum number; K c : Asymmetric top quantum number; Γ rot :

Rotational symmetry in C 2v (M); C i : Largest coefficient used in the TROVE assignment; n 1 -n 3 : Vibrational (TROVE) quantum numbers; i vib : Vibrational state counting number; N vib Population density for each i based on 1 -3 . This column is not produced by the ExoMol format. It must be calculated separately.

e.g. we have used to model the CO non-LTE spectra, 6 recently, which are planning to explore in the future. The methods described here can be used model the intensity distribution of the reaction products and to ascertain from what molecule the SiH 2 dissociated from. The equilibrium structured parameters (bond lengths and angles) can be treated as effective parameters to be adjusted to reproduce the experimental spectra. 

Proc. Nat. Acad. Sci

Non-LTE radiative transfer J o u r n a l N a me

Radiation transport in spectral lines

Molecular Symmetry and Spectroscopy