Guangfeng Liu

Jie Liu

Andrew S Dunn

Peter Nadazdy

Peter Siffalovic

Roland Resel

Mamatimin Abbas

Guillaume Wantz

Yves Henri Geerts

Directional Crystallization from the Melt of an Organic p-Type and n-Type Semiconductor Blend

come

L'archive ouverte

INTRODUCTION

Pristine organic semiconductors possess a low electrical conductivity, and usually, their dielectric constant is in the order of 3-4. However, a new situation occurs when they are doped, as shown by Heeger, McDiaramid, and Shirakawa in 1977. [START_REF]benzothiophene derivatives[END_REF] Doping, which consists of the external introduction of charges into semiconductors with chemicals, electric field, or light, gives rise to an increase of electrical conductivity, up to 11 orders of magnitude in the case of polyacetylene. [START_REF] Yamamoto | Electrical Properties of Phthalocyanine-Halogen Complexes[END_REF] Doped organic semiconductors have been the topic of intense research activities for more than four decades. [START_REF] Jacobs | Controlling Molecular Doping in Organic Semiconductors[END_REF][START_REF] Lüssem | Doped Organic Transistors[END_REF][START_REF] Salzmann | Molecular Electrical Doping of Organic Semiconductors: Fundamental Mechanisms and Emerging Dopant Design Rules[END_REF] However, some fundamental questions remain unanswered, even today. Why, for example, the ratio of free charge carriers to the number of dopant molecules, ηdop, is so low, reaching only an order of magnitude of 10% in most systems? This might be caused by the low dielectric constant of organic semiconductors that opposes charge separation, the clustering of dopant molecules, the trapping of charge carriers, or a strong hybridization of dopant and matrix molecules. [START_REF] Lüssem | Doped Organic Transistors[END_REF][START_REF] Salzmann | Molecular Electrical Doping of Organic Semiconductors: Fundamental Mechanisms and Emerging Dopant Design Rules[END_REF] While relatively uniform dispersion and distribution of dopants can be achieved in conjugated polymers, the situation is more complicated for small molecular semiconductors, especially when they crystallize. The dispersion and distribution of doping agents depend on various factors, such as: i) the crystal lattice of the host semiconductor and the size of the guest atoms or molecules, ii) the stoichiometry, iii) the crystallization conditions, iv) the occurrence of structural defects, v) the presence of grain boundaries and vi) the miscibility and extent of charge transfer between donors and acceptors. [START_REF] Schweicher | Molecular Semiconductors for Logic Operations: Dead-End or Bright Future?[END_REF] So far, only a few doped single crystals of molecular semiconductors have been reported. For example, Hiramoto et al. have reported single crystals of rubrene doped with minute amounts of iron(III) chloride. 7 In a different field of research, the directional crystallization of thin films of organic molecules in the presence of a deliberately introduced impurities of various amounts has been extensively studied for the occurrence of Mullins-Sekerka instabilities and pattern formation at the growth front. [8][9][10][11][12] Typically, eutectic mixtures forming one miscible phase in melt and two separate crystal phases have been subjected to directional crystallization in a temperature gradient, ∇T, using an experimental set-up as described in Figure 1. 13 At the growth front, the interface between the liquid and solid phases, crystal growth and diffusion of the phases compete. The length-scale of phase separation between the two crystal phases is linked to the cooling rate at the growth front, Cfront (°C min -1 ), that is externally imposed by the lateral withdrawal velocity, Vp (µm s -1 ), and the magnitude of the thermal gradient, Gexp, (°C mm -1 ). As a function of experimental conditions, digitated, dendritic, or even more complex growth morphologies are essentially observed. [13][14][15] In this paper, we combine these two aforementioned seemingly separate research fields to control the morphology of thin films of 2,7-dioctyl [1]benzothieno [3,2-b][1]benzothiophene (C8-BTBT-C8) with various amounts of tetracyanoquinodimethane (TCNQ). These two organic semiconductors have been selected because the crystal structure of their charge transfer complexes has already been reported 16,[START_REF] Tsutsumi | Ntype field-effect transistors based on layered crystalline donor-acceptor semiconductors with dialkylated benzothienobenzothiophenes as electron donors[END_REF] (Figure S1 and Table S1), allowing a detailed structural elucidation of thin-film morphology. The latter is defined here as comprising of crystal shape, crystal size, crystal orientation, crystal alignment, less ordered regions, defects, eventual polymorphs, and grain boundaries. Our goal is to obtain parallel crystalline stripes of C8-BTBT-C8 separated by stripes of C8-BTBT-C8-TCNQ charge transfer complexes (C8-BTBT-C8-TCNQ CT).

EXPERIMENTAL SECTION

1. Materials. C8-BTBT-C8 and 2,7-bis(dodecyloxy) [1]benzothieno [3,2-b]-benzothiophene (C12-BTBT-C12) were synthesized according to literature. [START_REF] Ebata | Highly Soluble [1]Benzothieno[3,2-b]benzothiophene (BTBT) Derivatives for High-Performance, Solution-Processed Organic Field-Effect Transistors[END_REF] Purification of the materials was performed through a chromatographic column using hexane and the purity assessed by NMR. [START_REF] Ruzié | Design, synthesis, chemical stability, packing, cyclic voltammetry, ionisation potential, and charge transport of[END_REF] 7,7,8,8-tetracyanoquinodimethane (TCNQ, >98%, Sigma-Aldrich), polydimethylsiloxane (PDMS, trade name: SYLGARD®184, Sigma-Aldrich), and HPLC-grade (>99.8%, Sigma-Aldrich) toluene solvent were used without further purification. Glass substrates of size 15 mm × 15 mm × 0.16 mm (D 263 Borosilicate thin glass from Marienfeld, Germany) were sonicated for 10 minutes in deionized water and acetone, then dried with nitrogen and exposed to an ultraviolet-ozone (Procleaner plus; Bioforce Nanosciences) treatment for 30 minutes before each experiment.

Glass surface modification by PDMS (SYLGARD®184):

PDMS was mixed with a curing agent in a 10:1 w/w ratio. Then, this solution was diluted in toluene in a 1:4 (w/w) ratio. Next, the diluted PDMS solution was spin-coated onto the pre-cleaned 15 mm × 15 mm × 0.16 mm glass substrates at a speed of 6000 rpm. Finally, glass substrates coated with PDMS were heated for 2 hours at 90℃ for completely curing. The thickness of the PDMS layer, measured by an Alpha-Step D-600 Stylus Profiler, KLA-Tencor, was found to be ~ 350 nm. The surface roughness (Ra) of the layer, measured using a Veeco NT9080 optical profilometer, was around 7 nm. For the preparation of powder blended samples, C8-BTBT-C8 and TCNQ were mixed in molar ratios 1:1, 10:1 and 20:1.

Then, each blend was ground with a mortar and pestle for 10 min at room temperature in order to disperse TCNQ molecules in the rich phase of C8-BTBT-C8. For the preparation of thin films, 0.5 mg of C8-BTBT-C8:TCNQ blended powder (either 10:1 or 20:1) was deposited on a 15 mm × 15 mm × 0.16 mm pre-cleaned thin glass substrate (or a PDMS/glass substrate) and covered by another identical glass substrate (or a PDMS/glass substrate). Then, the whole sample was placed on the hot stage at 130℃ in order to obtain a uniform liquid crystal phase. Each thin film was then produced by processing the sample by directional crystallization. After the treatment, the cover substrate was removed for X-ray diffraction analysis. For the thin films prepared on pure glass substrates, removing the cover glass resulted in almost equal amounts of material on both glass slides (both still usable for X-ray diffraction measurements). However, for the thin films prepared on PDMS/glass substrates, removing the cover PDMS/glass led most of the material on the bottom PDMS/glass substrate.

Temperature gradient apparatus.

The set-up shown in Figure 1, consists of a Linkam GS350 system presenting two distinct heating stages separated by a gap. The distance (gap) between the two stages where the thermal gradient was generated was 2 mm. A 76 mm × 26 mm × 1 mm microscope glass slide (Marienfeld Cat. No. 1000000) was intercalated between the stages and the sample to ensure a constant displacement velocity of the sample. During the ∇T experiments, the hot stage was set at a temperature Th = 130℃; while the cold stage was maintained at a temperature Tc = 80℃. The system is covered by a hermetic lid so that the system is thermally independent of the laboratory environment.

Calibration of the magnitude of the temperature gradient set-up (Gcal).

Previous work [START_REF] Schweicher | Directional crystallization of C8-BTBT-C8 thin films in a temperature gradient[END_REF] has shown that the effective magnitude of the temperature gradient (Gexp) that takes place between the hot and cold zones is less than the magnitude calculated by the equation as G = (Th -Tc)/x, where x = 2.0 mm (the gap between the hot and cold stages). The calibration of the ∇T set-up has been carried out as previously reported and is detailed in supporting information (Table S2 and Figure S2). Gexp ≈ 10°C mm -1 for the Th -Tc couple as 130°C -80°C.

Thermogravimetric analysis (TGA). TGA measurements were performed with a Perkin

Elmer Pyris 6 TGA instrument. The powder samples of pure C8-BTBT-C8, C8-BTBT-C8:TCNQ blends and pure TCNQ underwent a linear ramp of increasing temperature up to 500℃ at a rate of 30℃ min -1 in an inert atmosphere of nitrogen. Moreover, to further test the thermal stability, the C8-BTBT-C8:TCNQ blended powder samples were also heated to 160℃ and 200℃ and maintained under flowing nitrogen gas for 1 hour.

Differential scanning calorimetry (DSC).

DSC measurements were performed for pure C8-BTBT-C8 powder, C8-BTBT-C8:TCNQ blended powders and pure TCNQ under a dry nitrogen gas atmosphere in a NETZSCH POLYMA DSC 214 instrument (temperature range 25℃ -160℃).

First, each sample was measured using 10℃ min -1 heating-cooling rates with three heating-cooling cycles. Then, to investigate the effect of heating and cooling rate on DSC peak shapes, a heating-cooling rate of 1℃ min -1 with three heating-cooling cycles was adopted for each sample. The data were analyzed using the software Netzsch Proteus version 7.1. 21

Polarized Optical Microscopy (POM).

The liquid crystals and crystallization behaviors of the powder samples of pure C8-BTBT-C8, C8-BTBT-C8:TCNQ blends were recorded at different temperatures by a hot-stage-POM (Linkam LTS350-Nikon Eclipse 80i). Another hot-stage-POM (Linkam GS350-Nikon Eclipse 80i) allowed in situ observations of the temperature gradient crystallization of the thin films. Images of the thin film samples were recorded during and after the thermal gradient treatments.

Crystal morphology calculation.

The crystal morphology of C8-BTBT-C8 and C8-BTBT-C8:TCNQ CT was modeled from the crystal structure obtained from Cambridge Structural Database (CSD) with codes YOKBIK [START_REF] Ebata | Highly Soluble [1]Benzothieno[3,2-b]benzothiophene (BTBT) Derivatives for High-Performance, Solution-Processed Organic Field-Effect Transistors[END_REF] and WOTWUZ, [START_REF] Tsutsumi | Ntype field-effect transistors based on layered crystalline donor-acceptor semiconductors with dialkylated benzothienobenzothiophenes as electron donors[END_REF] respectively, using Mercury 3.10.3 software based on the Bravais-Friedel-Donnay-Harker (BFDH) method. [START_REF] Bravais | Etudes Cristallographiques[END_REF][START_REF] Friede | Crystal Morphologies[END_REF][START_REF] Donnay | A new law of crystal morphology extending the law of Bravais[END_REF] 8. X-ray measurements. Specular X-ray diffraction (sXRD): sXRD patterns of the thin films were collected on a Bruker D8 Advance diffractometer using Cu Kα radiation (λ = 1.5418 Å). Measurements were conducted at room temperature in the 2θ range between 1.6° and 30° with an increment of 0.02° per step and a counting time of 3s per step. XRD curves are plotted as the scattering intensity versus the scattering wave vector, qz (Å -1 ), using the equation q = 4π × sinθ/λ, where λ = 1.5418 Å. Powder XRD patterns (PXRD) for the powder samples of pure C8-BTBT-C8 and C8-BTBT-C8: TCNQ blends were collected on a Rigaku Ultima IV diffractometer using Cu Kα radiation. Measurements were conducted at room temperature in the 2θ range between 1.6° and 50° with an increment of 0.02° per step and a counting time of 0.24 s per step. The variable temperature PXRD patterns of the samples were recorded using a DHS 1100 (Anton Paar company) heating stage on the Rigaku instrument, with the same scan parameters and 10 ℃ min -1 heating-cooling rates.

Grazing-incidence X-ray diffraction (GIXD):

GIXD experiments for thin films were carried out on the XRD1 beamline at the Elettra synchrotron (Trieste), using standard set-up. The X-ray wavelength in these experiments was 1.4 Å. The scattered intensities were collected by a twodimensional X-ray detector (Pilatus 2M, Dectris). Calibration was performed using LaB6 standard.

An incident angle of αi = 3º was chosen to reduce the width of the diffraction peaks. During each GIXD experiment, the samples were rotated by 360°, integrating the intensity during a rotation angle of 2°. In total, 180 diffraction patterns were acquired for each sample. Data processing was performed by the software GIDVis [START_REF] Schrode | GIDVis: a comprehensive software tool for geometry-independent grazing-incidence X-ray diffraction data analysis and pole-figure calculations[END_REF] for phase analysis and calculation of pole figures. STEREOPOLE [START_REF] Salzmann | software for the analysis of X-ray diffraction pole figures IDL[END_REF] software was used for the interpretation of the pole figures.

Confocal Raman microscopy (CRM):

Spatially-resolved Raman spectra were measured by a commercial CRM system (Alpha300 R+, WITec). The excitation laser wavelength was 532 nm.

The spectrophotometer (UHTS 300, WITec), equipped with a 600 gr/mm grating (blazed at 500 nm) and coupled to an EMCCD camera (Newton DU970N BV 353, Andor), was employed to collect Raman spectra through a 100 m diameter multimode optical fiber that also served as a confocal pinhole. The measured spectra were acquired and analyzed using WITec Project Plus software.

RESULTS AND DISCUSSION

Preparation and thermotropic behaviors of C8-BTBT-C8: TCNQ blends.

The thermal stability of pure C8-BTBT-C8, TCNQ, and C8-BTBT-C8:TCNQ blends (10:1 and 20:1 molar ratios), up to 200°C, were assessed by TGA; blends of 1:1 tend to be slightly less stable (see detailed description in the supporting information, Figure S3, S4 and S5). DSC measurements indicate that TCNQ has no transition in the temperature range from 25°C to 160°C (Figure S6), whereas C8-BTBT-C8 exhibits the expected Cr to LC (liquid crystal phase) and LC to isotropic melt phase transitions: at 109.8°C and 126.0°C, respectively (Figure 2a). Variable temperature PXRD and POM results confirm the presence of the known crystal structure of pure C8-BTBT-C8 at room temperature and the smectic-A (SmA) phase at high temperature (see Figure S13a, and S14, see also Table S4). These results agree with previous reports. [START_REF] Ebata | Highly Soluble [1]Benzothieno[3,2-b]benzothiophene (BTBT) Derivatives for High-Performance, Solution-Processed Organic Field-Effect Transistors[END_REF][START_REF] Schweicher | Directional crystallization of C8-BTBT-C8 thin films in a temperature gradient[END_REF] For investigating the thermotropic phase behaviors of the blends, the starting point is 1:1 ratio of C8-BTBT-C8 and TCNQ. Grinding, which is an incomplete mixing process, affords a coarse blend with grains of C8-BTBT-C8, TCNQ, and C8-BTBT-C8-TCNQ CT, as shown by PXRD (Figure S11). Raising the temperature enhances mixing and allows electron transfer from C8-BTBT-C8 to TCNQ crystals. This is corroborated by the variable temperature PXRD patterns of 1:1 blend (Figure S13b). However, at 160 °C, some Bragg peaks are still observed in the XRD pattern, indicating that the blend does not form an isotropic melt. This is confirmed by the POM image at 160 °C, which shows some texture (Figure S15). The second and the third heating-cooling DSC cycles (10 °C min -1 and 1 °C min -1 scanning rates) of 1:1 blend are comparable with one another (Figures S7b andS7c, Figure S8b andS8c). However, the first heating cycle is different because of incomplete mixing by grinding (Figure S7a and Figure S8a). The similarity of the second and third heating curves indicates that C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT have reached phases with equilibrium compositions. This conclusion is also corroborated by the similarity of the relative crystallization enthalpies (Figure S7, S8 and Table S3). These results highlight that a partial electron transfer occurs between C8-BTBT-C8 and TCNQ in 1:1 blend. An equilibrium takes place as described in Scheme 1, which depends upon temperature but also on the phase behavior. Since the phase behavior of 1:1 blend is particularly complex, we have decided to investigate the thermal behavior of 10:1 and 20:1 blends, expecting that they would deviate less from that of pure C8-BTBT-C8. Figure S12a shows the comparison of the PXRD pattern of the powder sample of pure C8-BTBT-C8, 10:1 and 20:1 blends, after grinding, at room temperature. Only Bragg peaks belonging to the crystal structure of C8-BTBT-C8 crystals can be observed in the PXRD patterns, indicating that TCNQ molecules were molecularly dispersed in the C8-BTBT-C8 rich phase.

Notably, C8-BTBT-C8-TCNQ CT co-crystals were not detected in these two blends at room temperature. Figure S12b and Figure 2c show the variable temperature PXRD patterns of 10:1 and 20:1 blends, which exhibit similar phase behavior, i.e. at 90°C, C8-BTBT-C8-TCNQ CT cocrystals can be detected. At 130℃, the blends transform into a SmA liquid crystal phase. At 160℃, the isotropic melt formed. Cooled down to room temperature, C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT co-crystals are detected but no TCNQ is observed (see detail description in the supporting information, also see Figure S12b, and Table S4). These diffraction results match the POM observations of the two blends at both 130 ℃ and 160 ℃ very well (Figure S16). Comparing the DSC curves of 1:1 blend, 10:1 and 20:1 blends, the DSC traces are much less complex (Figure 2b, Figure S9 andS10). In the first heating cycles of the blends, the LC to the isotropic melt phase transition temperatures (TLC-iso) for 10:1 and 20:1 blends are 135.0°C and 131.6°C, respectively, which are higher than for pure C8-BTBT-C8 (126.0°C). Impurities are known to depress melting points because the chemical potential of a mixture is always lower than that of its individual constituents. However, some electrostatic interactions arise from charge transfer complex formation, which increases TLC-iso. The second and third heating cycles, as well as the first, second and third cooling cycles, of the blends that are comparable with one another, demonstrate reproducibility and consistency of the system (see detailed description in the supporting information, also see Figure 2b, Figure S9, S10 and Table S3). In conclusion, the above results are consistent with the view that the blends have reached phases with equilibrium compositions after the first heating. Blends of 10:1 and 20:1 exhibit a more simple and comparable phase behavior to pure C8-BTBT-C8, meaning they qualify as good candidates for directional crystallization.

Directional crystallization of C8-BTBT-C8: TCNQ blends in ∇T.

Crystallization conditions.

The samples were prepared by adopting a sandwiched geometry, that is, depositing 0.5 mg of C8-BTBT-C8:TCNQ blended powder between two glass (or PDMS/glass) substrates (Figure 1).

No dewetting occurred for sandwiched films. Uncoated glass and PDMS coated glass were chosen to investigate the influence of the substrate on crystallization. The molar ratio of C8-BTBT-C8:TCNQ, and the cooling rate were varied to conduct the directional crystal growth of blended samples in ∇T, and to explore a broad set of crystallization conditions. For these, the cooling rates at the growth front Cexp were varied from 0.6 to 60 °C min -1 (films 1-14, Table S5).

Growth morphology and thin-film morphology of the blends on the glass substrate.

The crystal growth processes of samples 1-10 (on pure glass substrates) were monitored in situ by POM during directional crystallization (Figure 3). The digitated growth morphology [START_REF] Jackson | Kinetic Processes: Crystal Growth, Diffusion, and Phase Transitions in Materials[END_REF] was observed in films 3, 4, 8, and 9. On the other hand, films 1, 2, 6, and 7 show branched sea-weed type growth morphology. In addition, films 5 and 10 exhibit dendritic growth morphology. The results indicate that the cooling rate at the growth front Cexp (°C min -1 ) plays a vital role in determining the crystal growth morphology: moderate Cexp (6 and 12°C min -1 ) leads to the digitated growth, but too slow (≤ 3°C min -1 ) or too fast (≥ 60°C min -1 ) Cexp perturbs the growth (Figure 4).

As shown in Figure S17, after the ∇T crystallization, the crystalline stripes can be observed again in films 3, 4, 8, and 9. The width of the larger stripes ranges from 30 μm to 80 μm, and the width of the narrower stripes ranges from 5 μm to 30 μm (Table S6). Only faceted or mosaic crystals were found in other films (1, 2, 5, 6, 7 and 10, Figure S17). The thickness values of films, measured by a step profiler, range between 2~3 μm. However, many cracks through the thin films were found after cooling the films from Tc to room temperature. These originate from the differences in thermal expansion coefficients between organic crystals and the rigid glass substrate. Another problem is that after separating the two glass slides, almost equal amounts of material remained on the top and bottom glass substrates due to the relatively large adhesion of organic crystals to glass. Using PDMS/glass substrates improved these problems (vide infra).

To investigate the crystal structures, out-of-plane preferred orientation and in-plane alignment of the thin films, diffraction studies have been conducted on these films with in-house sXRD and synchrotron-based GIXD based facilities. Figure S21 shows the sXRD patterns of the thin films 1-10, indicating that the known crystal phase of C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT (Figure S1 and Table S1) that have been detected in the bulk crystals are also present in thin films.

No other polymorphs have been detected in the thin films at room temperature. Diffraction data demonstrate a high preferential orientation of the (00l) lattice planes parallel to the glass substrates for both C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT crystallites. Synchrotron-based 2D-GIXRD patterns displayed in Figure S24 (for films 3, 4, 8 and 9) exhibit a texture with a low mosaicity.

All spots in GIXD data can be indexed by the known crystal phase of C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT. [START_REF] Tsutsumi | Ntype field-effect transistors based on layered crystalline donor-acceptor semiconductors with dialkylated benzothienobenzothiophenes as electron donors[END_REF][START_REF] Ebata | Highly Soluble [1]Benzothieno[3,2-b]benzothiophene (BTBT) Derivatives for High-Performance, Solution-Processed Organic Field-Effect Transistors[END_REF] Synchrotron results also corroborates a lack of polymorphism in the thin films and a high preferential orientation of the (00l) planes for both crystallites, which are consistent with the observations made by sXRD.

Growth morphology and thin-film morphology of the blends on the PDMS/glass substrate.

The formation of cracks is a common problem when preparing organic thin films from melt crystallization. Therefore, a soft PDMS layer was used to coat glass substrates. Figure 5a and Figure S18 show that films 11-14 on the PDMS/glass substrate keep a digitated growth front during the ∇T crystallizations. As expected, after directional crystallization, fewer cracks are present (Figure 5b and Figure S19) in the samples compared to the samples prepared on pure glass substrates, resulting from the soft viscoelastic nature of PDMS. The comparison of the dark-field microscopy picture of film 9 (on the glass substrate) with film 14 (on the PDMS/glass substrate) also indicates that the soft PDMS layer can improve the cracking problem (Figure S20). Samples of 20:1 ratio (films 13 and 14) show better parallel striped crystals than 10:1 ratio samples (films 11 and 12). The widths of crystallite stripes in films 13-14 correspond to the values of the thin films obtained on glass substrates (Table S6).

The morphology of stripes for sample 14 was further analyzed using confocal Raman imaging.

First, Raman spectra of pure TCNQ and C8-BTBT-C8 materials were measured (Figure S21). For mapping the spatial distribution of TCNQ molecules, the following vibrational modes at 1212 cm -1 , 1460 cm -1 and 1607 cm -1 can be used, which belong to C=CH bending, C-CN stretching and C=C ring stretching, respectively [START_REF] Gucciardi | Nano-Raman imaging of Cu-TCNQ clusters in TCNQ thin films by scanning near-field optical microscopy[END_REF] . However, the first two overlap with the C8-BTBT-C8 vibrational modes at 1208 cm -1 and 1472 cm -1 , while only the vibrational mode at 1607 cm -1 has no overlap with any of the C8-BTBT-C8 vibrational modes [START_REF] Mai | Effects of Ambient Gases on the Electrical Performance of Solution-Processed C8-BTBT Thin-Film Transistors[END_REF][START_REF] Li | In situ purification to eliminate the influence of impurities in solution-processed organic crystals for transistor arrays[END_REF] . For this reason, the peak intensity at 1607 cm -1 was used to map the spatial distribution of TCNQ shown in Figure 6a. The well-defined horizontal stripes indicate regions of increased TCNQ concentration in agreement with the POM images (Figure 5). In addition, the numerically integrated vertical line scans across the entire Raman map are shown in Figure 6b. The graph reveals that the concentration contrast between the TCNQ-rich and TCNQ-depleted regions is at a level of some 50%. For assessing the in-plane alignment of the crystallites in the thin films, pole figures were calculated based on the GIXD data. The particular {113} reflection for C8-BTBT-C8 and {01-1} reflection for C8-BTBT-C8-TCNQ CT were chosen because of their strong intensity. The {113} reflection of C8-BTBT-C8 crystallites in film 14 (Figure 9a) shows that four symmetrical discrete poles are located at the tilt angle of about 60° and the rotation angles of 60°, 165°, 240° and 345°.

Only four narrow discrete spots indicate that the C8-BTBT-C8 crystallites show a uniaxial alignment in film 14. The rotation angles of the discrete spots also reveal that the C8-BTBT-C8 molecules in film 14 were well in-plane aligned along the [010] or [0-10] direction, which is parallel to the gradient direction of the ∇T approach (Figure 10a). This result agrees well with the POM images shown in Figure 5a and 5b. It should be noted that processing pure C8-BTBT-C8 films by directional crystallization leads to no uniaxial alignment. [START_REF] Schweicher | Directional crystallization of C8-BTBT-C8 thin films in a temperature gradient[END_REF] The absence of uniaxial alignment has also been observed for C8-BTBT-C8-TCNQ CT in Figure 9b. Pole figures of films 3, 4, 8 and 9 (glass substrates) and films 11-14 (PDMS/glass substrates) are shown in Figure S26a and Figure S26b, respectively. These results corroborate those of film 14, except that some discrete {113} reflections are missing, most likely due to weak signal intensity. In summary, both C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT crystallites show a high preferential orientation of the (00l) planes, which means that C8-BTBT-C8 molecules and C8-BTBT-C8-TCNQ CT molecules preferentially stand upright on the glass substrate or on the PDMS/glass substrate (Figure 10b and

CONCLUSIONS

The directional crystal growth of semiconductors has been demonstrated so that separated areas of p-and n-type semiconductors are formed. Variable temperature powder XRD patterns (PXRD), differential scanning calorimetry (DSC) and polarized optical microscopy (POM) results show that 10:1 and 20:1 molar ratios of C8-BTBT-C8:TCNQ blends exhibit a phase behavior rather comparable to pure C8-BTBT-C8, qualifying them for directional crystallization. The results indicate that the crystal growing modes (digitated, dendritic, or branched sea-weed type) of the blends are linked to the cooling speed, which is externally imposed by the lateral withdrawal velocity and the magnitude of the thermal gradient. On both glass and PDMS substrates, the parallel crystalline stripes of C8-BTBT-C8 and C8-BTBT-C8-TCNQ CT have been obtained at a moderate cooling rate (6 or 12°C min -1 ), and a clear assignment of the different stripes to the different crystalline phases are obtained by Raman microscopy. The comprehensive characterization of the structure and the morphology of blend thin films was completed by POM, in-house XRD and synchrotron-based GIXD measurements. The results confirm a preferential orientation of the crystallites with the (001) plane parallel to the substrate surface for both types of phases. Moreover, the in-plane alignment of both phase types along the thermal gradient direction has been confirmed. While C8-BTBT-C8 shows high crystalline alignment with low in-plane mosaicity, the order is less pronounced for the charge transfer complex C8-BTBT-C8-TCNQ CT.

Figure 1 .

 1 Figure 1. Schematic diagram of the ∇T set-up. Tiso and TCr are the isotropic melt and crystallization

Scheme 1 .

 1 Scheme 1. Partial electron transfer occurring between the donor C8-BTBT-C8 and the acceptor

Figure 2 .

 2 Figure 2. DSC traces of (a) pure C8-BTBT-C8 and (b) C8-BTBT-C8:TCNQ blend with molar

Figure 3 .

 3 Figure 3. The crystal growth front of the blended films 1-10 on glass substrates under different

Figure 4 .

 4 Figure 4. A diagram shows the sample zones with branched sea-weed growth, dendritic growth,

Figure 5 .

 5 Figure 5. POM images of film 14 prepared on PDMS coated glass with a molar ration 20:1, a

Figure 6 .

 6 Figure 6. Confocal Raman mapping. a) Distribution map of the TCNQ Raman signal at 1607 cm -1 .

Figure 7 .

 7 Figure 7. Specular X-ray diffraction patterns of film 14 together with calculated diffraction

Figure 8 .

 8 Figure 8. (a) Reciprocal space maps calculated from GIXD investigations for film 14 plotted in

Figure

  Figure S27).

Figure 9 .

 9 Figure 9. (a) Pole figure of the {113} reflection of C8-BTBT-C8 crystallites, and (b) of the {01-

Figure 10 .

 10 Figure 10. (a) Schematic representation of C8-BTBT-C8 molecular orientation versus the

  The Supporting Information is available free of charge on the ACS Publications website at DOI: XXX, including the crystal structure information, calibration of the magnitude of the temperature gradient set-up, TGA and DSC curves, variable temperature PXRD patterns, optical microscopy images, Raman spectra, in-house sXRD patterns, synchrotron-based GIXD and Pole-figures. Corresponding Author * guangfeng.liu@ulb.be (7) Ohashi, C.; Izawa, S.; Shinmura, Y.; Kikuchi, M.; Watase, S.; Izaki, M.; Naito, H.; Hiramoto, M., Hall Effect in Bulk-Doped Organic Single Crystals. Adv Mater. 2017, 29 (23), 1605619. (8) Mullins, W. W.; Sekerka, R. F., Morphological Stability of a Particle Growing by Diffusion or Heat Flow. J. Appl. Phys. 1963, 34 (2), 323-329. (9) Mullins, W. W.; Sekerka, R. F., Stability of a Planar Interface During Solidification of a Dilute Binary Alloy. J. Appl. Phys. 1964, 35 (2), 444-451. (10) Jackson, K. A.; Hunt, J. D., Transparent compounds that freeze like metals. Acta Met. 1965, 13 (11), 1212-1215. (11) Langer, J. S., Instabilities and pattern formation in crystal growth. Rev. Mod. Phys. 1980, 52 (1), 1-28.(12) Caroli, B.; Caroli, C.; Roulet, B., The Mullins-Sekerka instability in directional solidification of thin samples. J. Crystal Growth. 1986, 76 (1), 31-49. (13) Godreche, C., Solids far from equilibrium. Cambridge University Press, 2011 (14) Akamatsu, S.; Faivre, G.; Ihle, T., Symmetry-broken double fingers and seaweed patterns in thin-film directional solidification of a nonfaceted cubic crystal. Phys. Rev. E 1995, 51 (5), 4751-4773. (15) Ginibre, M.; Akamatsu, S.; Faivre, G., Experimental determination of the stability diagram of a lamellar eutectic growth front. Phys. Rev. E 1997, 56 (1), 780-796. (16) Shibata, Y.; Tsutsumi, J. y.; Matsuoka, S.; Matsubara, K.; Yoshida, Y.; Chikamatsu, M.; Hasegawa, T., Uniaxially oriented polycrystalline thin films and air-stable n-type transistors based on donor-acceptor semiconductor (diC8BTBT)(FnTCNQ) [n = 0, 2, 4]. Appl. Phys. Lett. 2015, 106 (14), 143303.

BTBT n° 2.4565.11, Phasetrans n° T.0058.14, Pi-Fast n° T.0072.18, and 2D to 3D No. 30489208.

Financial supports from the French Community of Belgian (ARC n° 20061) and by the Walloon

Region (WCS No. 1117306, SOLIDYE n° 1510602) are also acknowledged. We thank the synchrotron Elettra, Trieste for provision of synchrotron radiation and Luisa Barba for her support to use the beamline XRD1.


Andrei V Nikitin

Alexander E Protasevich

Alena A Rodina

Michael Rey

Attila Tajti

Vladimir G Tyuterev

V E Zuev

Vibrational levels of formaldehyde: Calculations from new high precision potential energy surfaces and comparison with experimental band origins

Keywords: formaldehyde, potential energy surface, vibrational energy levels, kinetic energy operator

Introduction

Accurate knowledge of spectroscopy of formaldehyde is important for the monitoring of this species in the atmosphere and also for detection of the presence of this molecule as an environmental pollutant indoors [START_REF] Barbe | Measurements of tropospheric and stratospheric H2CO by an infrared high resolution technique[END_REF] [2] [3] [START_REF] Zhu | Validation of satellite formaldehyde (HCHO) retrievals using observations from 12 aircraft campaigns Atmos[END_REF]. Improving spectral parameters and line-by-line analyses are essential for reliable interpretation of near infrared spectra returned by the groundbased and orbiting observatories [START_REF] Steck | Retrieval of global upper tropospheric and stratospheric formaldehyde (H2CO) distributions from high-resolution MIPAS-Envisat spectra[END_REF] [START_REF] Dufour | Global uppertropospheric formaldehyde: Seasonal cycles observed by the ACE-FTS satellite instrument[END_REF]. An assignment of infrared spectra corresponding to high-energy overtone and combination bands is known to be a tedious and difficult task. Recent ab initio calculations of potential energy surfaces (PES) and of dipole moment surfaces (DMS) together with a progress in variational nuclear motion calculation have permitted a significant progress in predictions of vibration-rotation bands. Previous experience with some 4-5 atomic molecules [START_REF] Rey | Accurate Theoretical Methane Line Lists in the Infrared up to 3000 K and Quasi-continuum Absorption/Emission Modeling for Astrophysical Applications[END_REF] [8] [START_REF] Wong | Atlas of Experimental and Theoretical High-temperature Methane Cross Sections from T=295 to 1000K in the Near-infrared The[END_REF] [10] [START_REF] Rey | Understanding global infrared opacity and hot bands of greenhouse molecules with low vibrational modes from first-principles calculations: The case of CF4[END_REF] have shown that ab initio based approaches could be useful for better modelling and assignments of high-resolution spectra. In the present work, we report theoretical vibrational levels of formaldehyde computed from new accurate ab initio PES using the algorithms similar to one recently applied for the methane molecule [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF] and compare our results with other recent calculations [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] [START_REF] Yachmenev | A new "spectroscopic" potential energy surface for formaldehyde in its ground electronic state[END_REF].

From the theoretical standpoint, the formaldehyde molecule is an interesting case because of the rapid basis set convergence of variational calculations and a possibility to compare different methods conducted in various coordinate systems using high-level quantum-chemical calculations. Burleigh et al. [START_REF] Burleigh | An accurate quartic force field for formaldehyde[END_REF] have been used the quadratic force field for formaldehyde obtained by refining the ab initio quartic force field CCST(T)/cc-pVTZ of Martin, Lee, and Taylor [START_REF] Martin | An accurate ab initio Ouartic Force Field for Formaldehyde and its Isotopomers[END_REF] by means of iteratively fitting a subset of the coefficients of a Taylor-series expansion of the potential-energy surface to the set of the observed transition frequencies.

Yachmenev et al. [START_REF] Yachmenev | A new "spectroscopic" potential energy surface for formaldehyde in its ground electronic state[END_REF] have calculated ab initio PES using CCSD (T) / aug-cc-pVQZ method in a wide range of geometries. A detailed review of earlier works refs [START_REF] Burleigh | A random matrix approach to rotation-vibration mixing in H2CO and D2CO[END_REF] [15] [START_REF] Martin | An accurate ab initio Ouartic Force Field for Formaldehyde and its Isotopomers[END_REF] [START_REF] Mladenovic | Discrete variable approaches to tetratomic molecules Part II: application to H2O2 and H2CO[END_REF] [START_REF] Yagi | Ab initio potential energy surface for vibrational state calculations of H2CO[END_REF] [START_REF] Al-Refaie | ExoMol line lists -VIII. A variationally computed line list for hot formaldehyde[END_REF] that reported calculations of the formaldehyde PES and energy levels can be found refs. [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] [14].

The accuracy of ab initio calculations is in general lower than that of high-resolution experimental methods, except for studies of diatomic molecules with a small number of electrons. In the case of triatomics, very accurate ab initio calculations of potential energy surfaces (PES) and of dipole moment surfaces (DMS) led to a breakthrough in the extension of spectra analyses towards higher energy ranges. This was, for example, the case with the spectroscopy of water [START_REF] Partridge | The determination of an accurate isotope dependent potential energy surface for water from extensive ab initio calculations and experimental data[END_REF] [START_REF] Polyansky | Water on the sun:line assignments based on variational calculations[END_REF] [23] [START_REF] Polyansky | Calculation of Rotation-Vibration Energy Levels of the Water Molecule with Near-Experimental Accuracy Based on an ab Initio Potential Energy Surface[END_REF], carbon dioxide [START_REF] Huang | An isotopic-independent highly accurate potential energy surface for CO2 isotopologues and an initial (C16O2)-C-12 infrared line list[END_REF] [26] [START_REF] Yurchenko | A Variationally Computed T = 300 K Line List for NH3[END_REF], hydrogen sulfide [START_REF] Cours | Ab initio dipole moment functions of H232S and intensity anomalies in rovibrational spectra[END_REF] [29] and sulfur dioxide [START_REF] Huang | Highly accurate potential energy surface, dipole moment surface, rovibrational energy levels, and infrared line list for (SO2)-S-32-O-16 up to 8000 cm-1[END_REF], as well as with the ozone spectra analyses [START_REF] Vl | New analytical model for the ozone electronic ground state potential surface and accurate ab initio vibrational predictions at high energy range[END_REF] [32] [START_REF] Vl | Does the "reef structure" at the ozone transition state towards the dissociation exist? New insight from calculations and ultrasensitive spectroscopy experiments[END_REF]. In the latter case, ab initio predictions were mandatory for precise intensity quantification [START_REF] Tyuterev | Ab initio predictions and laboratory validation for consistent ozone intensities in the MW, 10 and 5 µ m ranges[END_REF], for understanding the PES properties in the transition state range towards the dissociation threshold [START_REF] Dawes | Highly accurate ozone formation potential and implications for kinetics[END_REF] [START_REF] Vl | Does the "reef structure" at the ozone transition state towards the dissociation exist? New insight from calculations and ultrasensitive spectroscopy experiments[END_REF], the interactions between potential wells [START_REF] Kokoouline | Localized and delocalized bound states of the main isotopologue 48O3 and of 18O-enriched 50O3 isotopomers of the ozone molecule near the dissociation threshold[END_REF], and modeling of the isotopic exchange reactions [START_REF] Yuen | The Role of Ozone Vibrational Resonances in the Isotope Exchange Reaction 16O16O+18O→18O16O+16O: The Time-Dependent Picture[END_REF]. Many ab initio PESs have been subsequently refined by a fit to experimental data to achieve better accuracy in line positions. Theoretical line lists for ammonia [START_REF] Huang | Rovibrational spectra of ammonia. I. Unprecedented accuracy of a potential energy surface used with nonadiabatic corrections[END_REF] [START_REF] Yurchenko | A Variationally Computed T = 300 K Line List for NH3[END_REF] and phosphine [START_REF] Yurchenko | Ab initio dipole moment and theoretical rovibrational intensities in the electronic ground state of PH3[END_REF] [40] [START_REF] Nikitin | High order dipole moment surfaces of PH3 and ab initio intensity predictions in the Octad range[END_REF] are successful examples of this trend for four-atom species. New PESs and DMSs have been recently reported for five-atom [START_REF] Owens | Accurate ab initio vibrational energies of methyl chloride[END_REF] [START_REF] Owens | A global potential energy surface and dipole moment surface for silane[END_REF], six-atom molecules [START_REF] Delahaye | A new accurate ground-state potential energy surface of ethylene and predictions for rotational and vibrational energy levels[END_REF] [START_REF] Delahaye | Accurate 12D dipole moment surfaces of ethylene[END_REF], [START_REF] Li | High-Level, First-Principles, Full-Dimensional Quantum Calculation of the Ro-vibrational Spectrum of the Simplest Criegee Intermediate (CH2OO)[END_REF] and seven-atom [START_REF] Nikitin | First Full-Dimensional Potential Energy and Dipole Moment Surfaces of SF6[END_REF] molecules.

Among the large number of ab initio calculations for many-electron systems, only two works reported calculations of the vibrational levels with the accuracy of about 0.1 cm -1 , at least in the low and medium energy ranges: ref [START_REF] Polyansky | Calculation of Rotation-Vibration Energy Levels of the Water Molecule with Near-Experimental Accuracy Based on an ab Initio Potential Energy Surface[END_REF] [START_REF] Polyansky | ExoMol molecular line lists XXX: a complete high-accuracy line list for water[END_REF] for the water molecule, and ref [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF] for methane.

While both of these molecules contain 10 electrons, the formaldehyde molecule has 16 electrons.

In the recent study by Morgan et al. [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] the high accuracy vibrational fundamentals of formaldehyde have been reported, but not all calculated vibrational energy levels were published.

In contrast to methane, formaldehyde is an experimentally much less studied molecule. One of the possible reasons for this is the unavailability of the cold spectra caused by the high temperature of congelation and a complex technique for creating the formaldehyde gas. While for methane about 200 assigned vibrational energy levels have been measured with high accuracy, in case of formaldehyde, no more than 20 vibrational levels could be accurately obtained from the high resolution rotationally resolved spectra. The calculations reported in this paper can be considered as an additional source of information that could be useful for extended assignment and modeling of high-resolution experimental spectra of formaldehyde and its isotopologues.

This paper is structured as follows. The choice of the coordinates, the kinetic energy operator (KEO) and the analytical representations for the PES including their benchmark comparisons are discussed in Section 2. The dependence of the theoretical vibrational levels on the electronic basis set, the impact of relativistic corrections, diagonal Born-Oppenheimer corrections (DBOC), and of high-order electron correlation including connected triple and quadruple excitations in the coupled-cluster method is considered in Section 3. A detailed comparison with available experimental data in Sections 3 and 4 shows that the resulting fulldimensional ab initio PES provides currently the most accurate theoretical band origins for formaldehyde. We also discuss uncertainties and discrepancies in published experimental levels obtained from low-and high-resolution spectra. An empirically optimized PES obtained with a fine tuning of six PES parameters and the corresponding predicted vibrational levels are also reported.

Coordinates and analytical PES representation

For calculations of the vibrational-rotational energy levels of the formaldehyde molecule we used the kinetic energy operator (KEO) and PES written in terms of three sets of internal coordinates based on the angle-bisector molecule fixed axes [START_REF] Lukka | A simple method for the derivation of exact quantum-mechanical vibration-rotation Hamiltonians in terms of internal coordinates[END_REF], [START_REF] Colwell | The derivation of vibration-rotation kinetic energy operators in internal coordinates[END_REF] . The sets of coordinates differ only in the definition of the internal vectors: (I) bond angle-bond length coordinates in which KEO was derived in ref [START_REF] Colwell | The derivation of vibration-rotation kinetic energy operators in internal coordinates[END_REF], (II) Jacobi-like coordinates (see later), (III) Radau-Jacobi coordinates from refs [START_REF] Bramley | A general discrete variable method to calculate vibrational energy levels of three-and four-atom molecules[END_REF] and [START_REF] Poulin | Calculation of vibrational (J=0) excitation energies and band intensities of formaldehyde using the recursive residue generation method[END_REF]. The internal coordinates (I) and (II) are non-orthogonal sets while (III) is an orthogonal coordinate set. The kinetic energy operator of the nuclei (using the notation from ref [START_REF] Makarewicz | Rovibrational Molecular Hamiltonian in Mixed Bond-Angle and Umbrella-Like Coordinates[END_REF]) in the laboratory frame r starts from the CH 1 H 2 center of mass. Figure 1 shows these three sets of the internal vectors. Elements of the matrix of transformation to internal vectors and the effective masses for three sets of coordinates.

Set of coordinates

ik A , 1 ij m (I) 4 ik k   4 11 ij i mm   (II)   3 1 2 4 1 k ik k m m m m      33 1 2 4 11 1 ij i j i m m m m        (III)         3 3 1 2 4 3 3 4 4 12 1 1 1 1 k ik i k k i k k k m m m m m mm                          a) 3 1 2 4 11 ij i ij i m m m m      a) the parameter   4 1 2 4
m m m m    

. Here we applied the same tensor techniques for an optimal sampling of the grid of nuclear configurations as described in the previous works for non-abelian symmetry groups [START_REF] Nikitin | Vibration energy levels of the PH3, PH2D, and PHD2 molecules calculated from high order potential energy surface[END_REF] [55]. This permitted accounting for the full symmetry of the molecule in order to reduce the number of geometrical nuclear configurations for ab initio calculations of electronic energies. To build the corresponding grid of geometries, we used at the first step six one-dimensional sections(Figure 2 show three angular sections). The extent of the grid points was chosen in a way which ensures that a maximum number of parameters of our analytical PES representation would be well-defined in the least-squares fits to ab initio electronic energies. PES parameters responsible for the coupling of various vibrational modes were systematically included. The corresponding initial reference grid (G (R) ) contained 22509 points. This relatively dense grid provides a possibility for testing various high-order PES expansions in elementary functions of the symmetrized coordinates. The grid points were arranged in the order of increasing energy with the cut-off corresponding to the extent of the studied infrared bands. For this grid, the initial calculations were carried out with the coupled-cluster CCSD(T) method using the frozen core approximation and the quintuple-zeta cc-pV5Z basis set (that will be referred to as V5Z sets in abbreviated notation) and the MOLPRO program [START_REF] Werner | Molpro: a general purpose quantum chemistry program package[END_REF] [57].

The analytical PES was constructed using the irreducible tensor formulation described in refs [54] [55]. The PES can be expressed as a power series in elementary functions of the symmetrized coordinates involving Morse-type functions for "radial" and a sine function of angular variables: ( , , , , , )

i j k l m n CO ijklmn ijklmn V r r r f S S S S S S     
, where (2a) To validate various representations of the PES using the least-squares technique, an extended grid of geometries was used. This grid was augmented by adding supplementary sets of points corresponding to 180- values with respect to all values of  (for 180  

    1 1 2 2 5 1 2 ( ) ( ) / 2, ( ), ( ) ( ) / 2 CO S y r y r S y r S y r y r      , ( 2b
)     3 1 2 6 1 2 ( ) ( ) / 2, ( )
( ) / 2, ( ) sin( ) e S f f S f f f              . ( 2c
) in the initial grid.

Figure 2. One-dimensional PES sections for angular symmetrized coordinates.

To compare the efficiency of the analytical form (2) in various coordinate sets, the benchmark PES calculations were carried out using the least-squares fit to V5Z grid values, without the further corrections discussed in the following sections. Even though the V5Z basis set did not provide the best quantitative accuracy, we considered that it was preferable to use all PESs at the same theoretical level without corrections for a reliable comparison of the general qualitative shape of the surfaces. Ab initio potential energies were fitted using the analytical symmetry adapted representation (3) and the weight function employed by Schwenke and

Partridge in ref [58] . 2 2.00200200 2) 1.00200200 ) ( 0005 . 0 tanh( ) ( 1     E E E w (3)
This weight function decreases with the energy E (expressed in cm -1 ) in order to de-emphasize the contribution to the PES fit of large grid displacements for vibrational energies beyond

E 1 =15000cm -1 .
The least-squares fit of the PES in various sets of coordinates leads to quite similar results, but for the coordinate sets II and III, the RMS deviations slightly decrease at E 1 > 12000 cm -1 . Table 2 summarizes the RMS deviations of the fit in various coordinates. The PES expansion was extended up to sixth order in all S i -coordinates including the torsion function cos(τ)+1 but some eighth-order terms were added for torsion functions cos(τ /2) and τ-π.

The number of parameters, which were statistically well determined in the fit of the PES to ab initio electronic energies, versus the total number of parameters of the corresponding expansion orders, is given in the second row of Table 2. The results show that the coordinate sets (II) and

(III) provide a more accurate description of the PES when using the weight function (3). When fitting a six-dimensional PES, all three forms of the torsion function result in close lying standard deviation (STD) values but the use of (II) and (III) provides slightly lower RMS values.

Ref [START_REF] Yachmenev | A new "spectroscopic" potential energy surface for formaldehyde in its ground electronic state[END_REF] reported a more complicated PES representation, but the use of that representation was not found useful in our case.

Note: Abbreviation KEO stands for kinetic energy operator

To control the equilibrium geometry, we calculated the rotational energy levels combining the V5Z PES with the equilibrium geometries optimized using larger basis sets. This corresponds to a series of calculations using the aug-cc-pCVXZ core-valence basis sets augmented with diffuse functions for one-particle basis cardinal numbers X = 4 and 5 (denoted simply as ACVQZ, and ACV5Z). We have also accounted for the Douglas-Kroll (DK) type relativistic corrections. The parameters of the equilibrium geometry affect the vibrational energy levels, especially those associated with the angular motions [START_REF] Nikitin | Rotational and vibrational energy levels of methane calculated from a new potential energy surface[END_REF]. Table 3 gives the equilibrium geometry parameters, while Table 4 gives the energy levels at J=1 for the equilibrium geometries mentioned in Table 3. Note that changes up to the fourth decimal place do not affect the vibrational energy levels.

That is why, the calculations for all geometries except for the first one, led to almost identical results for the low lying rotational levels. The accuracy of the known experimental rotational levels was estimated to be about 10 -4 cm -1 . In this work, however, we did not have an intention to obtain highly accurate values of the rotational levels, but to find the appropriate parameters of the equilibrium geometry consistent with reliable vibrational predictions. The diagram shows that the difference between the fourth and the first calculation is considerable in the range over 8000 cm -1 . This is mainly caused by the different torsion parts.

Below 7000 cm -1 , all calculations lead to close results. Below 4500 cm -1 , the influence of the PES type on the energy levels is not significant: the differences in the level positions do not exceed 0.05 cm -1 .

Improvement of the basis set and corrections to the PES

Since all versions of calculations for levels up to 7000 cm -1 using the V5Z basis in the previous sections lead to close results, we shall use the internal coordinate set (I) and torsion function cos(τ/2) while improving the PES at higher levels of the theory. The choice of the shape is due to the fact that the internal coordinates are easier for calculating corrections.

First of all, we conducted calculations according to the procedure described for the methane molecule in ref [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF] using the augmented core-valence basis set ACV5Z at a sparser grid of geometries. As these calculations are more demanding, we computed electronic energies with the CCSD(T)/ACV5Z level of theory at the grid ( G (AC) ) containing only 3820 points in the range up to about 6000 cm -1 from the equilibrium geometry. Similarly to ref [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF], the energy differences Δ AC = ACV5Z-V5Z involving two basis sets were computed on this grid. We have checked that these differences behave quite smoothly with respect to the symmetrized coordinates S i and were possible to be modeled analytically using the fourth-order expansions Δ AC (S i ) of the type (2a-2c). This allowed us to extend the electronic energies according to the relation ACV5Z V5Z + Δ AC (S i ) for all remaining geometries and refitting the PES(ACV5Z) on the full reference grid G (R) of 22509 points using the analytical form (2a) as discussed in the previous section.

As the next step, we accounted for the following three corrections to the surface: "massvelocity-Darwin" (MVD) relativistic corrections, the diagonal Born-Oppenheimer correction (DBOC), and higher-order dynamic electron correlations (HODEC). [75], whereas Morgan et al. [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] have reported a very thorough study of the impact of triple and quadruple excitations to the equilibrium geometry and the quartic force field of formaldehyde.

In this work, we calculated MVD corrections with MOLPRO on the full reference grid G (R) of 22509 points. Much more expensive calculations for DBOC and HODEC corrections were carried out at 2000 points of a smaller grid G (Corr) of nuclear geometries and approximated by a 4 th order expansion of the type (2). This gave a quite smooth analytical function of all six symmetrized coordinates for the sum of these three corrections:

Δ (3corr) (S i ) = Δ Rel (S i ) + Δ DBOC (S i ) + Δ T(Q) (S i ). (4)
The DBOC and HODEC corrections were calculated using the CFOUR program suite [START_REF] Stanton | CFour program[END_REF] with the cc-pVTZ one-particle basis set. The HODEC corrections on the G (Corr) grid were computed using the noniterative quadruple CCSDT(Q) method [START_REF] Kucharski | The coupled-cluster single, double, triple, and quadruple excitation method[END_REF] [72] [73] [START_REF] Matthews | Accelerating the convergence of higher-order coupled cluster methods[END_REF]. At this step we have constructed the ab initio corrected surface

PES(ACV5Z+3Corr) = PES(ACV5Z) + Δ (3corr) (S i ) (5)
Table 5 gives the energy levels using the PES(ACV5Z) and PES(ACV5Z+3Corr) that account for the above considered corrections, analytically extrapolated from the sparser grid G (Corr) .

It is seen that a values of the HODEC dynamic correction for formaldehyde computed with the CCSDT(Q) method is several times larger than those for the methane molecule [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF] . The differences between the calculated and measured levels are considerably larger for the v 2 overtones. Although the PES (ACV5Z + 3Corr) for formaldehyde is constructed in a similar way as the PES for methane, its accuracy is significantly lower than that for methane. It was necessary to improve the potential energy surface, particularly along the CO bond.

Since the relativistic corrections and DBOC are relatively small, we have tried to explore further contributions related to HODEC correlations including quadruple excitations, as well as to the larger, aug-cc-pCV6Z basis set. To this end we used at the next step the iterative CCSDTQ [START_REF] Matthews | Non-orthogonal spin-adaptation of coupled cluster methods: A new implementation of methods including quadruple excitations[END_REF] [74] [START_REF] Stanton | CFour program[END_REF] method which is expected to provide more rigorous account for dynamic electronic correlations corresponding to quadruple excitations. However, this method converged very slowly for low-symmetry geometries or for large r CO bond distances. of radial coordinates (2b). The choice of analytical representations had a marginal effect on vibrational energy levels.

Our final ab intio PES including these four corrections was constructed as

PES(ACV5Z+4Corr) = PES(ACV5Z) + Δ (3corr) (S i ) + Δ Q+6Z (S 1 ,S 2 ,S 5 ) (6)
The fourth correction Δ Q+6Z improves the accuracy of the calculated fundamental frequencies and brings them closer to the measured values, especially for v 2 . Also, the correction diminished the calculated value of the equilibrium geometry r COе . Contrary to the r CO correction, the r CH correction does have a significant effect on the energy levels. The partial contributions of the considered corrections to vibrational energies and the values of the deviations are shown in Figure 6. The contributions of various corrections to vibrational levels and our final ab initio results are collected in Table 5. The experimental energy levels are taken from refs [START_REF] Perez | Observation of hot bands in the infrared spectrum of H2CO[END_REF] [78] [START_REF] Flaud | High resolution spectroscopy of H2 12C16O in the 1.9 to 2.56 μm spectral range[END_REF], and [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF].

Note that the (Exp.-Calc.) deviations exceed 1 cm -1 only for six levels among a total of 35 experimentally known levels. The deviation of -1 cm -1 for 2ν 2 is quite reasonable considering that the deviation for ν 2 is equal to -0.44 cm -1 . Four levels from the other five have been taken from ref [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF]. This work reported the level values with only one significant digit after the decimal point. Ref [START_REF] Flaud | High resolution spectroscopy of H2 12C16O in the 1.9 to 2.56 μm spectral range[END_REF] reported another vibrational band center ν 2 + ν 4 + ν 6 derived from 40 rotational-vibrational energy levels including J values up to 17 with an uncertainty of 0.33 cm -1 .

Unfortunately, this work [START_REF] Flaud | High resolution spectroscopy of H2 12C16O in the 1.9 to 2.56 μm spectral range[END_REF] did not mention a minimum value of J for experimental transitions in the ν 2 + ν 4 + ν 6 band. If the fitting procedure does not use a sufficient number of the assigned transitions with J < 6, the accuracy of the extrapolated empirical J=0 level might be poor. A large difference can be explained by the fact that the two calculations deviate in different directions. Four other band centers above 4400 cm -1 have been reported in ref [START_REF] Flaud | High resolution spectroscopy of H2 12C16O in the 1.9 to 2.56 μm spectral range[END_REF]: ν 1 + ν 2 at 4529.5025(50) cm -1 , ν 2 + ν 5 at 4571.6947(50) cm -1 , 2ν 2 + ν 6 at 4734.2078(50) cm -1 , and 3ν 2 at 5177.75952(70) cm -1 . The corresponding calculated levels at 4530. [START_REF] Yagi | Ab initio potential energy surface for vibrational state calculations of H2CO[END_REF], 4571.72, 4734.95, and 5179.37 cm -1 are slightly higher than the empirical values, but the most significant part of the overestimation is due to the deviation of ν 2 levels, which can be roughly estimated as -0.44n

where n is number of ν 2 quanta. ), a -from ref [START_REF] Perrin | New analysis of the v2, v3, v4, and v6 bands of formaldehyde H2 12C16O line positions and intensities in the 5-10 lm spectral region[END_REF], b -from ref [START_REF] Perrin | New analysis of the 2v4, v4+v6, 2v6, v3+v4, v3+v6, v1, v5, v2+v4, 2n3, v2+v6 and v2+v3 bands of formaldehyde H2 12C16O: Line positions and intensities in the 3.5 mm spectral region[END_REF], c -from ref [START_REF] Perez | Observation of hot bands in the infrared spectrum of H2CO[END_REF], dfrom ref [START_REF] Flaud | High resolution spectroscopy of H2 12C16O in the 1.9 to 2.56 μm spectral range[END_REF], with one decimal place from ref [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF]. #) Rel.= relativistic; T(Q) = CCSDT(Q) -CCSD(T) correlation; Q+6Z = quadruple Δ Q+6Z correction (see text for explanation).

##) Experimental minus calculated levels corresponding to the full ab initio PES(ACVZ6 + 4Corr).

$

) RMS without five suspicious levels given in Table 6.

Note that most of the "outlier" levels in our (exp.-calc.) with deviations larger than 1 cm -1 also showed significant deviations from empirical values in calculations of Burleigh et al. [START_REF] Burleigh | An accurate quartic force field for formaldehyde[END_REF] and of Yachmenev et al. [START_REF] Yachmenev | A new "spectroscopic" potential energy surface for formaldehyde in its ground electronic state[END_REF]. Table 6 compares three calculations for these five levels. [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF], [START_REF] Flaud | High resolution spectroscopy of H2 12C16O in the 1.9 to 2.56 μm spectral range[END_REF] Note that empirical levels reported by Bouwens et al. [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF] in 1996 have been obtained from low resolution dispersed fluorescence spectra with an average uncertainty of 1 cm -1 . However, some of these empirical levels could have larger errors as they differ considerably from values obtained from high-resolution spectra analyses in more recent papers. For example, the study of Bouwens et al [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF] reported the ν 3 + ν 6 level at 2496.11 cm -1 and the ν 1 + ν 4 level at 3940.21 cm - 1 . In a more recent paper of A. Perin et al. [START_REF] Perrin | New analysis of the 2v4, v4+v6, 2v6, v3+v4, v3+v6, v1, v5, v2+v4, 2n3, v2+v6 and v2+v3 bands of formaldehyde H2 12C16O: Line positions and intensities in the 3.5 mm spectral region[END_REF] [81] the same levels were obtained from analyses of vibration-rotation spectra at 2494.3543 cm -1 and at 3941.5295 cm -1 , respectively. Figure 7 shows the diagram illustrating the differences between empirical levels obtained from dispersed fluorescence spectra [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF] and those obtained from infrared rotationally resolved spectra [START_REF] Perrin | New analysis of the 2v4, v4+v6, 2v6, v3+v4, v3+v6, v1, v5, v2+v4, 2n3, v2+v6 and v2+v3 bands of formaldehyde H2 12C16O: Line positions and intensities in the 3.5 mm spectral region[END_REF] [81]. A comparison of fundamental band centers computed from our final ab initio PES including all corrections with the results of the recent, very thorough study by Morgan et al [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] is given in Table 7. According to notations of ref [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF], here MAE ( 11

n ii i XY n    ) is the mean (Exp.-Calc.) deviation, while MAPE ( 1 100% n ii i i XY nX   
) is the mean absolute percent error defined by Eq.( 17) of ref [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF]. In order to improve the accuracy of vibrational predictions, we have also produced an empirically optimized PES, by a fitting of 6 harmonic second order terms of the potential expansion to 20 levels experimentally determined with three decimal digits (reported in Table 5 excepting the level at 4163.289 cm -1 ). With this fine tuning of the six PES parameters, the RMS (exp.-calc.) deviation for 20 accurate levels derived from high-resolution spectra fall down from 0.21 cm -1 to 0.008 cm -1 . Figure 8 shows the diagram of differences between the experimental levels and those calculated from our best ab initio PES and from the empirically optimized PES. The experimental values derived from low resolution dispersed fluorescence spectra of Bouwens et al [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF], for which the accuracy is considerably lower than for levels obtained from analyses of rotationally resolved infrared spectra using an effective Hamiltonian technique, are shown in the rectangular frame. It is seen that without the levels in the frames, the errors of our ab initio levels are smaller than the discrepancies between the low resolution data of Bouwens et al [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF] and the experimental levels obtained in more recent works (see Figure 7). It can thus be concluded that the accuracy of our final ab initio calculations is not worse than that of low resolution experimental levels reported in [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF], at least for the bands with doubly and triply excited vibrational quanta. Energy levels up to 7000 cm-1, calculated from our empirically fitted PES in orthogonal coordinates are shown in the first column of Table 8. These coordinates permit using the KEO, which is considered as "exact" (not requiring series expansion), however the latter term is only meaningful within the framework of the Born-Oppenheimer approximation. In order to calculate energy levels in normal coordinates the energy values generated from the analytical form of the PES in orthogonal coordinates were computed on our G (R) grid of geometries (see paragraph 2) and then re-fitted in internal coordinates with the torsion function of (τ -π). The TENSOR code [START_REF] Rey | Group-theoretical formulation of an Eckart-frame kinetic energy operator in curvilinear coordinates for polyatomic molecules[END_REF] [84] [85] [START_REF] Rey | Convergence of normal mode variational calculations of methane spectra: Theoretical line list in the icosad range computed from potential energy and dipole moment surfaces[END_REF] was then used to calculate the energy levels in normal coordinates using the Watson-Eckart KEO. Most likely, the accuracy of the levels computed from empirically optimized PES for doubly and triply excited quanta is higher than the accuracy of the levels obtained using the pure ab initio PES. However, in case of the bands corresponding to four or more vibrational quanta, the error caused by an empirical correction might be larger than that for levels obtained using the initial ab initio PES, because only 20 low energy experimental levels were included in the fit. Note that four band centers ν 1 + ν 2 , ν 2 + ν 5 , 2ν 2 +ν 6 , and 3ν 2 were not used in fit of the PES, but the corresponding energy levels calculated from empirically optimized PES became noticeably (RMS four times lower) closer to the empirical values. The difference between the two calculations in Table 8 is several times bigger than the initial difference between the corresponding Calc 3 and Calc 4 results for the ab initio V5Z PES (figure 3), although the calculations are converged in both cases. This could occur due to the limited order of the PES expansion and to numerical issues in the procedure of the coordinate transformations on a final grid of points. For this reason, the predicted values of high-energy levels in Table 8 have to be taken with caution. They could be useful in the initial steps of the assignment, if no other, more precise values are available.

Deuterated D2CO isotopologue

The vibrational levels of D 2 CO given in Table 9 were calculated, using nuclear masses, from two potential surfaces. The calculation (CalcD 2 CO 1) corresponds to our best ab initio PES(ACV5Z+4Corr) and (CalcD 2 CO 2) corresponds to the empirically optimized PES fitted to twenty vibrational levels of the main formaldehyde isotopologue as described in the previous section. Experimental levels were taken from refs [START_REF] Perrin | New High-Resolution Analysis of the v3 , v4 , and v6 Bands of D2CO Measured by Fourier Transform Spectroscopy[END_REF] (the first 3 levels), [START_REF] Lohilahti | The v2 Band of Formaldehyde-d2[END_REF] (ν2), and [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] (the remaining 8 levels). Note that ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] reports both the vibrational levels obtained from experimental transitions from the ground state to the upper level J=0 (see Table 3 in ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF]), as well as the spectroscopic parameter of the band center derived with the least squares technique.

The spectroscopic parameter of the band center and vibrational level in the same units (cm-1) should have closely lying values in the absence of anharmonic Fermi-type resonances, but could considerably differ however, in the case of a strong Fermi coupling. Table 10 of ref [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] erroneously gives the spectroscopic parameter 2054.694 cm -1 instead of the experimental vibrational level 2060.91963 cm -1 for ν 1 . The experimental vibrational level v 1 is considerably closer to the calculated value of ref [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF]. In Table 9, we compare the experimental and calculated levels. Note that all calculated levels, except for 2ν 4 , are close to experimental ones. The choice of various representations of potentials, coordinates, and the size of the vibrational basis do not have any considerable effect on the calculated energy levels given in Table 9. The D 2 CO isotopologue is heavier than H 2 CO, and the DBOC correction for D 2 CO is expected to be less important than this correction for H 2 CO. In case of D2CO (CalcD 2 CO 1), the DBOC was not recalculated. To estimate the contribution of DBOC, we presented it in the form of an expansion in a second-order Taylor series. Comparison of parameters of the DBOC of the two isotopologues D2CO/H2CO shows that the contribution of DBOC_D2CO is ~ 0.6 *DBOC_H2CO. Since the contribution is weak, the -0.4 contribution for the main isotopologue (see column DBOC of Table 5) could be added to the calculated D2CO levels (CalcD 2 CO 1). In this case, levels ν 1 (2060.887cm -1 ), ν 2 (1701.258 cm -1 ), ν 5 (2162.657 cm -1 ) would become a little closer to the empirical values.

Since a considerable number of transitions in the 2ν 4 band have been assigned and accurate combination differences have been derived in ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF], an outlier in the 2ν 4 band center printed in Table 3 of ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] could likely be interpreted as a misprint. On the other hand, the authors of ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] have stated (the second sentence of the fourth paragraph): "One can clearly see the three strong fundamental bands ν 2 , ν 1 , and ν 5 with the centers located near 1701.6 cm -1 , 2060.9 cm -1 , and 2162.9 cm-1, respectively. The weak bands 2ν 4 , 1867.8 cm -1 , ν 4 + ν 6 , 1930.0 cm -1 , and 2ν 6 , 1974.3 cm -1 , can also be recognized very clearly in spite of their weakness". These band centers are close (within 0.1 cm -1 ) to reported experimental vibrational levels of ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF], except for the experimental vibrational level 2ν 4 , which was shifted by 2 cm -1 with respect to the band center cited above. Note, that the center of the 2ν 4 band reported in ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] (1867.8 cm -1 ) is very close to our calculated vibrational level of 2ν 4 (1867.6 cm -1 ). Figure 7 of ref [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] indicates the presence of a Fermi resonance between the 2ν 4 and ν 1 bands. However, Table 6 of the same paper [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF] missed the parameter of this resonance.

The last row of Table 9 gives the RMS values of differences between experimental and calculated levels. The empirically optimized PES provides the best description of the energy levels, except for the doubted 2ν 4 value. [START_REF] Perrin | New High-Resolution Analysis of the v3 , v4 , and v6 Bands of D2CO Measured by Fourier Transform Spectroscopy[END_REF] [87] [START_REF] Lohilahti | High resolution infrared study of D2CO in the region of 1780-2400 cm-1: assignment and preliminary analysis[END_REF]

Conclusion

The recent thorough study of the quartic force field of the formaldehyde molecule by Morgan et al [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF] obtained at 183 nuclear geometries near the equilibrium configuration has shown the importance of high-order dynamic electron correlation by including triple and quadruple excitations in the coupled-cluster hierarchy via the CCSDT(Q) method. This has permitted improving fundamental band origins with respect to the standard CCSD(T) approach with an RMS deviation of 0.52 cm -1 [START_REF] Morgan | Geometric Energy Derivatives at the Complete Basis Set Limit: Application to the Equilibrium Structure and Molecular Force Field of Formaldehyde[END_REF]. In this work, we confirm these conclusions using the algorithm of the construction of the full six-dimensional PES at an extended set of geometries. To this end, we used successive grids adapted to the level of the theory that was previously applied for the methane molecule [START_REF] Nikitin | First fully ab initio potential energy surface of methane with a spectroscopic accuracy[END_REF]. These methods, including also relativistic and adiabatic DBOC corrections, allowed us to approach the spectroscopic accuracy in ab initio calculations for a larger range of vibrational energies. For the first timein case of a four-atom molecule with 16 electrons like H2CO -it was possible to obtain the (ab initioexperimental) RMS deviation of 0.4 cm -1 including 33 experimental vibrational band origins up to three vibrational quanta (Table 6).

First-principle calculations provide an independent insight into uncertainties of available experimental data, particularly for those which had been deduced from low-resolution dispersed fluorescence spectra [START_REF] Bouwens | Pure vibrational spectroscopy of S 0 formaldehyde by dispersed fluorescence[END_REF]. A detailed comparison of various sources in Table 7 and Figures 678suggests that many experimental levels could be less accurate than the present ab initio calculations. By excluding five suspicious band origins corresponding to large experimental uncertainties, the RMS deviation of our ab initio PES for vibrational levels up to 3700 cm-1 (Table 6) drops down to 0.2 cm -1 . A similar accuracy is obtained for vibrational levels of the deuterated isotopologue D2CO (Table 9), confirming the reliability of the ab initio PES.

We have also produced an empirically optimized PES by adjusting six second order parameters to twenty experimental vibrational levels derived from high-resolution spectra with the (obs-calc) RMS deviation of 0.008 cm -1 . Computed from this PES 160 band origins of H2CO up to 7000 cm -1 are given in Table 8. Both the ab initio and the empirically optimized PES are provided in the Supplementary Materials as a C++ code.

For a further investigation of the corresponding accuracy issues, the ab initio PES can be used to derive effective spectroscopic models by the contact transformation method as it was done in the case of methane [89] [90] [91]. This will permit to accurately compute physically meaningful values of the resonance coupling parameters from ab initio surface for advanced analyses of high-resolution spectra.

m 1 R , 2 R , 3 R , and 4 R 4 m 3 r 3 r 2 r

 12344332 can be written as a sum of the KEO of the center of mass and the T R relative motions of the nuclei, as follows: is the mass of the k th nuclei, and the linear transformation of the 4 us use the following numbering of indices: 1-for the nucleus of the first hydrogen atom H 1 , 2 for the nucleus of the second hydrogen atom H 2 , 3for oxygen O, 4for carbon С. The vectors start from the same point and define the positions of the H 1 , H 2 , O, and С nuclei, respectively. are masses of the H 1 , H 2 , O, and С nuclei, define the relative positions of the H 1 , H 2 , and O nuclei with respect to different centers of mass. The differences in the coordinate systems are as follows: for (I), start from the C nucleus position; for (II), start from the Radau canonical point for the CH 1 H 2 group, and 3

Figure 1 .

 1 Figure 1. Internal vectors (depicted in black) for three sets of the internal coordinates based on the angle bisector coordinate system for the the formaldehyde molecule : a) bond angle-bond length coordinates (I), b) Jacobi-like coordinates (II), c) Radau-Jacobi coordinates (III). The 1, 2, 3, and 4 indices denote the H 1 , H 2 , O, and C atom nuclei, respectively. Internal vectors with heads in the positions of nuclei 1, 2, and 3 have the following starting point: a) in the position of the nucleus 4, b) in the center of mass of the nuclei 1, 2,and 4(depicted in green) c) in the canonical point (depicted in yellow) of the nuclei 1, 2, and 4 (for vectors with heads in 1 and 2), and in the center of mass (depicted in green) of the nuclei 1, 2, and 4 (for the vector with head in 3).

Figure 3 .

 3 Figure 3. Comparison of the energy levels (the symmetry block A 1 ) for four calculations. The variational methods and the ab initio PES versions corresponding to Calc1, Calc2, Calc3 and Calc4 are defined in Table2.

Figure 3

 3 Figure 3 displays the diagram of differences in the calculated levels for the symmetry block A 1 .

  Full PES calculations at this level of the theory are extremely demanding, but it was possible to study the corresponding contributions on the most relevant one-dimensional cuts. Previous experiences in the PES corrections for various molecules showed that corrections performed along chemical bonds provide a maximum contribution. Note that for formaldehyde, all geometries keep the C 2v symmetry while changing the r CO bond length, whereas they have the C S symmetry while changing the r CH bond length. The correction Δ Q = CCSDTQ -CCSDT(Q) was evaluated for about twenty points along both bond stretching coordinates. The one-dimensional Δ AC6Z = ACV6Z-ACV5Z corrections were also computed at the same sets of geometries. The onedimensional corrections of the CO cut of the potential are shown in Figure4. These corrections result in an increase of the potential energy at distances longer than r e (CO). At the distances shorter than r e (CO), the potential energy slightly decreases. For the r CO bond, the Δ Q = CCSDTQ-CCSDT(Q) correction in the VTZ basis considerably exceeds the Δ AC6Z = ACV6Z-ACV5Z correction calculated with the CCSD(T) technique. For the r CH bond, the contributions of these corrections are of a similar magnitude. We summarized two corrections Δ Q+6Z = Δ Q + Δ AC6Z and approximated it analytically as a power series expansion in Morse-type functions Δ Q+6Z (S 1 ,S 2 ,S 5 )

Figure 4

 4 Figure 4 One-dimensional radial corrections along the CO bond(shown as solid line). The contributions Δ Q = CCSDTQ -CCSDT(Q) -the rigorous account for connected quadruple excitations in dynamic electron correlation -are shown as black crosses, while the Δ AC6Z =(ACV6Z-ACV5Z) and Δ V6Z = (V6Z-V5Z) corrections calculated at the CCSD(T) level are shown as circles and squares, respectively. The sum of two corrections Δ Q+6Z = Δ Q + Δ AC6Z shown as a solid line.

Figure 5

 5 Figure 5 One-dimensional radial corrections along the CH bond (shown as solid line). The contributions Δ Q = CCSDTQ -CCSDT(Q) -the rigorous account for connected quadruple excitations in dynamic electron correlation -are shown as black crosses, while the Δ AC6Z = (ACV6Z-ACV5Z) correction calculated at the CCSD(T) level is shown as circles. The sum of two corrections Δ Q+6Z = Δ Q + Δ AC6Z shown as a solid line.

Fig 6 .

 6 Fig 6. Upper panel(a): contribution of relativistic effects, DBOC, and high-order electron correlations T(Q) to vibrational energies up to 4500 cm -1 . Lower panel(b): Experimentalcalculated deviations for vibration levels corresponding to ab initio calculations using PES (ACV5Z+3Corrections) in blue and to PES(ACV5Z+4Corrections) in red. Both vertical scales are given in cm -1 .

Fig 8 .

 8 Fig 8. Experimental minus calculated deviations for vibrational levels computed from pure ab initio PES(ACV5Z+4Corr) and from the empirically optimized PES. Vertical scales are given in cm -1 . Levels from Bouwens et al. [80] and the suspicious level ν 2 + ν 4 + ν 6 from [79] are

  in cm -1 . **) RMS deviation without the suspicious level 2ν 4 . ***) See text in the beginning of the paragraph

Table 1

 1

gives the matrix elements ik A and values of , 1 ij m for these coordinate sets.

Table 1

 1

Table 2 .

 2 RMS and standard deviations of the PES fits in various coordinate sets.

Table 3 .

 3 Equilibrium geometry of the formaldehyde molecule, as optimized at various levels of ab initio theory

	Coordinates	CCSD(T)/	CCSD(T)/	CCSD(T)-DK/	Morgan [13] Empirical
		ACVQZ	ACV5Z	ACV5Z-DK**		values [59] [60]
	r e (CO) / Å	1.205240	1.20437	1.20426	1.20457	1.20461(19)
	r e (CH) / Å	1.100744	1.10041	1.10031	1.10052	1.10046(16)
	θ(HCH) / Deg	116.6186	116.6495	116.6546	116.694	116.722(93)
	θ(OCH)* / Deg 121.6906	121.6752	121.6726	121.653	121.639
	* θ(OCH)=π-θ(HCH)/2

** accounting for Douglas-Kroll-Hess relativistic corrections

Table 4 .

 4 The energy levels (in cm -1 ) obtained for the equilibrium geometries from Table3with the (V5Z) PES.

	J=1 Levels	Exp.*	ACV5Z**	ACV5Z**	Morgan	Empirical
			Watson KEO	"Exact" KEO	[13]	values [59] [60]
	A 2	2.4295	2.422452	2.426531	2.421981	2.422041
	B 1	10.5395	10.539508	10.537128	10.533067	10.531232

* Experimental values [61] [60]. ** Basis sets used for the equilibrium geometry. KEO stands for the kinetic energy operator used in the calculations

Table 5 . Vibrational levels with contributions of relativistic, DBOC corrections, and high-order electron correlations up to 4400 cm -1 .

 5

								This work				ref. [14]
				Contributions of		Theoretical levels with corrections
	Vib. State	С	Exp.	corrections	#						Calc	Exp-
			levels *	Rel.	DBO C	T(Q)	Q+6Z	ACV5Z	ACV5Z +3corr	ACV5Z +4corr	Exp-Calc ##	Calc
	ν4	B1	1167.25628(2) a	-0.14 -0.07 -3.82 0.55	1170.469	1166.43	1166.99	0.27	1166.1	1.16
	ν6	B2	1249.09568(2) a	-0.06 0.08	-2.24 0.65	1250.471	1248.25	1248.91	0.18	1245.6	3.49
	ν3	A1	1500.17574(12) a	0.06	-0.13 -2.28 0.76	1501.967	1499.62	1500.38	-0.21 1499.1	1.07
	ν2	A1	1746.00886(13) a	-0.91 -0.23 -8.86 3.75	1752.703	1742.69	1746.45	-0.44 1744.61 1.4
	2ν4	A1	2327.5239(5) b	-0.28 -0.14 -7.59 1.10	2333.903	2325.89	2326.99	0.53	2325.18 2.35
	ν4 + ν6	A2	2422.9701(50) b	-0.21 0.00	-6.10 1.22	2427.594	2421.29	2422.52	0.45	2418.43 4.54
	2ν6	A1	2494.3543(5) b	-0.13 0.15	-4.49 1.32	2497.191	2492.71	2494.03	0.32	2487.73 6.63
	ν3 + ν4	B1	2667.0481(20) b	-0.09 -0.19 -6.18 1.33	2672.132	2665.67	2667.00	0.05	2664.64 2.4
	ν3 + ν6	B2	2719.1550(5) b	-0.03 -0.15 -4.43 1.04	2722.697	2718.08	2719.13	0.03	2715.24 3.92
	ν1	A1	2782.4575(5) b	-0.25 -0.61 -4.78 0.34	2788.101	2782.46	2782.81	-0.35 2781.74 0.71
	ν5	B2	2843.3256(5) b	-0.37 -0.37 -6.14 0.90	2848.961	2842.08	2842.98	0.34	2842.37 0.95
	ν2 + ν4	B1	2905.9685(20) b	-1.06 -0.30 -12.8 4.34	2916.043	2901.80	2906.14	-0.18 2903.35 2.61
	2ν3	A1	2998.9873(5) b	0.11	-0.25 -4.69 1.55	3002.723	2997.89	2999.44	-0.46 2997.24 1.74
	ν2 + ν6	B2	3000.0659(5) b	-0.88 -0.19 -10.4 3.98	3007.897	2996.34	3000.32	-0.26 2995.91 4.16
	ν2 + ν3	A1	3238.4548(20) b	-0.89 -0.36 -11.3 4.62	3247.114	3234.49	3239.10	-0.65 3236.11 2.35
	2ν2	A1	3471.6	-1.81 -0.46 -17.8 7.58	3485.266	3465.10	3472.69	-1.09 3469.39 2.21
	3ν4	B1	3480.7	-0.41 -0.21 -11.3 1.65	3491.059	3479.11	3480.76	-0.06 3478.04 2.66
	2ν4 + ν6	B2	3586.6	-0.35 -0.06 -9.82 1.77	3594.654	3584.42	3586.19	0.41	3581.09 5.51
	ν4 +2ν6	B1	3673.5 $	-0.28 0.07	-8.42 1.86	3682.802	3674.16	3676.03	-2.53 3668.88 4.62
	3ν6	B2	-	-0.21 0.20	-6.79 1.89	3740.452	3733.65	3735.54	-	-	-
	ν3 + 2ν4	A1	3825.3	-0.24 -0.26 -9.98 1.89	3834.334	3823.86	3825.74	-0.45 3822.43 2.87
	ν3 + ν4+ ν6	A2	3886.5 $	-0.21 -0.25 -8.39 1.50	3895.439	3886.59	3888.10	-1.60 3883.37 3.13
	ν3+2ν6	A1	3937.4 $	-0.14 -0.14 -6.77 1.49	3940.654	3933.60	3935.09	2.31	3929.86 7.54
	ν1 + ν4	B1	3941.5295 c	-0.40 -0.66 -8.83 0.89	3950.568	3940.68	3941.57	-0.04 3939.85 1.68
	ν4 + ν5	A2	3996.518 c	-0.49 -0.40 -9.94 1.48	4005.302	3994.47	3995.94	0.57	3994.75 1.77
	ν1 + ν6	B2	4021.0806(6) d	-0.31 -0.50 -7.17 1.00	4028.157	4020.18	4021.18	-0.10 4017.7	3.38
	ν2 + 2ν4	A1	4058.3	-1.18 -0.36 -16.5 4.84	4071.474	4053.43	4058.28	0.02	4054.6	3.7
	ν5 + ν6	A1	4083.1	-0.52 -0.21 -9.30 2.07	4090.512	4080.47	4082.54	0.56	4078.88 4.22
	ν2 + ν4+ ν6	A2 4163.289(330) d,$	-1.06 -0.26 -14.6 4.66	4176.049	4160.08	4164.74	-1.45 4159.35 3.94
	2ν3 + ν4	B1	-	-0.05 -0.31 -8.65 2.13	4172.264	4163.25	4165.37	-	-	-
	2ν3+ ν6	B2	4192.3816(5) d	0.00	-0.33 -6.74 1.54	4198.047	4190.96	4192.51	-0.13 4189.62 2.77
	ν2 + 2ν6	A1	4248.7	-0.87 -0.15 -12.2 4.27	4256.837	4243.52	4247.79	0.90	4241.02 7.68
	ν1 + ν3	A1	4253.8 $	-0.18 -0.72 -7.00 1.07	4262.684	4254.79	4255.86	-2.06 4253.86 -0.06

Table 6 . Five vibrational levels from Table 5 with our (Exp.-Calc.) deviations exceeding 1 cm -1 .

 6

	Vib. State С	Exp. *	This work	ref. [15]		ref. [14]
	ν 4 +2ν 6	B 1	3673.5 ?	3676.032 -2.53	3675.9	-2.4	3675.07 -1.57
	ν 3 + ν 4 + ν 6 A 2	3886.5 ?	3888.101 -1.60	3886.7	-0.2	3886.47 0.03
	ν 3 +2ν 6	A 1	3937.4 ?	3935.095	2.31	3935.4	2.0	3936.02 1.38
	ν 2 + ν 4 + ν 6 A 2 4163.289 ? 4164.739 -1.45	4163.6	0.3	4164.12 -0.83
	ν 1 + ν 3	A 1	4253.8 ?	4255.859 -2.06	4255.3	-1.5	4254.58 -0.78

*) experimental values reported in refs

Table 7 . Experimental minus calculated deviations for fundamentals of H 2 CO.

 7

	Vib. State	C	Exp.	Exp-Calc This Work	Morgan [13] CCSDT(Q) /CBS	Morgan [13] CCSDT(Q) /CBS+Rel
	ν 4	B 1	1167.256	0.27	-0.8	-0.95
	ν 6	B 2	1249.095	0.18	0.03	-0.08
	ν 3	A 1	1500.175	-0.21	-0.03	0.01
	ν 2	A 1	1746.009	-0.44	-0.32	-1.48
	ν 1	A 1	2782.457	-0.35	1.18	0.49
	ν 5	B 2	2843.326	0.34	1.06	0.25
	MAE			0.29	0.57	0.54
	MAPE(%)			0.016	0.029	0.034
	RMS			0.31	0.73	0.75
	All values are given in cm -1 .

Table . Vibrational levels ( in cm -1 ) of the main formaldehyde isotopologue H 2 CO 478 computed from our empirically optimized PES using two sets of coordinates with two 479 versions of variational calculations.

 .

		A1 B2 A1 000600	4974.735 6196.629 6909.239	4974.536 6196.618 6909.337	0.199 0.011 -0.098
		B1 B2 B1 020300	4978.413 6211.767 6909.809	4978.371 6211.549 6909.812	0.042 0.218 -0.003
	480 481	B2 A1 A1 201000 A1 B1 A2 110101 Sym, vib labels* B1 A1 B1 011102	5044.315 6216.906 6914.801 5092.702 6236.880 6922.631 Calc A # 5111.129 6254.983 6966.384	5044.276 6216.889 6915.485 5092.730 6236.997 6922.694 Calc B ## 5111.015 6255.058 6966.388	0.039 0.017 -0.684 -0.027 -0.117 -0.063 Difference (A-B) 0.114 -0.075 -0.003
		B1 B2 A1 A2 003101	1167.263 5140.810 6265.909 6978.675	1167.247 5140.722 6265.782 6978.675	0.016 0.088 0.127 0.000
		B2 B2 B2 A1 110002	1249.103 5151.476 6275.837 6990.034	1249.044 5151.302 6275.861 6990.107	0.059 0.175 -0.024 -0.074
		A1 A1 A1 A2 A2 B1 *) Vibrational state labels include C 2v symmetry types and normal mode vibrational assignment. 1500.140 1500.138 0.002 1746.002 1746.005 5178.050 5178.078 -0.028 5186.406 5186.361 0.045 6283.377 6283.364 0.013 6322.849 6322.896 -0.047 #) Calc A. corresponds to direct calculations using our empirically optimized PES in orthogonal coordinates (I) and -0.003 A1 2327.519 2327.492 0.027 A2 2422.989 2422.912 0.076 A1 2494.406 2494.293 0.113 B1 5205.475 5205.445 0.029 B1 5246.933 5246.868 0.066 A1 5253.556 5253.499 0.058 A2 6335.338 6335.170 the torsion functions cos(τ/2) with the exact kinetic energy operator in polyspherical coordinates. 0.168 B2 6340.507 6340.531 ##) Calc B. corresponds to calculations using our empirically optimized PES refitted in normal coordinates on a -0.024 A1 6346.545 6346.522 0.022 geometrical grid with the Watson-Eckart kinetic energy operator.
		B1 B2 A1	2667.030 5312.412 6360.567	2667.017 5312.323 6360.303	0.014 0.089 0.264
		B2 A1 A1	2719.222 5322.146 6401.440	2719.197 5322.106 6401.405	0.026 0.041 0.036
		A1 B2 B1	2782.492 5325.537 6426.649	2782.499 5325.463 6426.645	-0.008 0.074 0.004
		B2 A2 A1	2843.343 5359.048 6437.667	2843.311 5359.007 6437.681	0.032 0.041 -0.014
		B1 A1 B2	2905.969 5383.971 6451.851	2905.958 5383.855 6451.859	0.011 0.116 -0.008
		A1 B1 A2	2998.953 5414.636 6466.652	2998.941 5414.625 6466.599	0.012 0.011 0.053
		B2 B1 B1	3000.159 5419.717 6472.643	3000.110 5419.635 6472.596	0.049 0.082 0.047
		A1 B2 B2	3238.402 5433.111 6476.996	3238.407 5433.318 6476.944	-0.005 -0.207 0.052
		A1 A1 A2	3471.801 5463.113 6495.501	3471.806 5463.393 6495.445	-0.005 -0.281 0.056
		B1 A2 B2	3481.518 5490.220 6510.481	3481.485 5490.169 6510.441	0.033 0.051 0.039
		B2 B2 A1	3586.906 5490.308 6546.991	3586.813 5490.176 6546.894	0.093 0.131 0.097
		B1 B2 B1	3676.685 5530.358 6553.062	3676.568 5530.434 6553.068	0.117 -0.075 -0.006
		B2 A1 A1	3736.082 5544.288 6564.842	3735.927 5544.261 6564.913	0.155 0.026 -0.071
		A1 A1 B2	3826.023 5552.705 6578.550	3825.994 5552.692 6578.313	0.030 0.013 0.237
		A2 A2 A1	3888.514 5627.353 6580.280	3888.480 5627.327 6580.214	0.034 0.026 0.066
		A1 A1 A2	3935.451 5651.311 6584.443	3935.337 5651.317 6584.991	0.114 -0.007 -0.548
		B1 B1 B1	3941.539 5660.831 6612.308	3941.542 5660.816 6612.932	-0.004 0.015 -0.624
		A2 B2 A1	3996.548 5665.587 6629.339	3996.503 5665.536 6629.655	0.045 0.051 -0.317
		B2 B1 B2	4021.081 5681.295 6636.922	4021.039 5681.319 6636.909	0.041 -0.025 0.013
		A1 A1 A1	4058.360 5687.375 6653.251	4058.336 5687.256 6653.314	0.024 0.119 -0.062
		A1 A2 A2	4082.935 5718.950 6670.621	4082.874 5718.932 6670.512	0.061 0.018 0.109
		A2 A1 A2	4164.835 5729.450 6683.934	4164.773 5729.448 6684.163	0.062 0.002 -0.229
		B1 B2 B2	4165.153 5765.234 6684.586	4165.135 5765.242 6685.095	0.018 -0.008 -0.509
		B2 A1 B1	4192.437 5769.551 6693.205	4192.407 5769.541 6693.178	0.030 0.010 0.027
		A1 B1 B1	4247.860 5772.142 6708.592	4247.772 5772.133 6708.624	0.087 0.009 -0.032
		A1 A1 A1	4255.327 5809.358 6727.367	4255.329 5809.322 6727.179	-0.002 0.036 0.188
		B2 B2 B2	4335.077 5822.554 6758.622	4335.028 5822.520 6758.625	0.049 0.033 -0.003
		B1 B1 A1	4397.710 5887.355 6766.827	4397.703 5887.344 6767.112	0.007 0.011 -0.285
		B2 A2 B2	4466.891 5889.405 6778.065	4466.869 5889.373 6778.064	0.022 0.032 0.001
		A1 B2 B1	4495.549 5892.767 6795.388	4495.533 5892.686 6795.461	0.016 0.081 -0.074
		A1 B2 A1	4529.420 5935.857 6813.391	4529.457 5935.799 6813.333	-0.037 0.057 0.057
		B2 A1 A1	4571.525 5985.373 6825.366	4571.519 5985.326 6825.462	0.006 0.046 -0.096
		B1 A1 A2	4624.282 5988.654 6829.984	4624.281 5988.651 6829.955	0.001 0.003 0.030
		A1 A1 A1	4629.634 5997.360 6838.098	4629.602 5997.369 6838.072	0.032 -0.009 0.026
		A1 B1 B1	4729.001 6002.841 6851.401	4728.995 6002.727 6851.322	0.006 0.113 0.079
		B2 B2 B2	4734.453 6052.566 6856.265	4734.419 6052.564 6856.337	0.035 0.002 -0.071
		A2 B2 B2	4743.299 6097.490 6864.634	4743.198 6097.335 6865.029	0.101 0.155 -0.395
		A1 B1 A1	4843.901 6107.429 6865.517	4843.775 6107.441 6865.643	0.126 -0.012 -0.126
		A2 A1 B1	4928.842 6124.578 6888.866	4928.695 6124.543 6888.915	0.147 0.034 -0.049
		A1 B1 B2	4955.909 6180.301 6894.266	4955.924 6180.124 6894.350	-0.015 0.177 -0.083
		A2 B2	6190.937 6904.180	6190.917 6903.974	0.020 0.206

Table 9 . Vibrational levels for the D 2 CO isotopologue up to 3000 cm -1 .

 9

	Vib. State	C	Exp. levels *)	CalcD 2 CO 1 *** Exp.-CalcD 2 CO 1	CalcD 2 CO 2***	Exp.-CalcD 2 CO 2
	ν 4	B 1	938.03549	937.936	0.099	937.711	0.324
	ν 6	B 2	989.25028	989.174	0.076	989.053	0.197
	ν 3	A 1	1100.44254	1100.403	0.039	1100.536	-0.093
	ν 2	A 1	1701.619103	1701.166	0.453	1701.631	-0.012
	2ν 4	A 1	1865.84444	1867.602	-1.757	1867.21	-1.365
	ν 4 + ν 6	A 2	1930.02676	1929.852	0.174	1929.524	0.502
	2ν 6	A 1	1974.32541	1974.232	0.092	1973.993	0.332
	ν 3 + ν 4	B 1	2038.90981	2038.801	0.108	2038.708	0.202
	ν 1	A 1	2060.91963	2060.643	0.275	2060.852	0.067
	ν 3 + ν 6	B 2	2072.66611	2072.585	0.080	2072.509	0.157
	ν 5	B 2	2162.92330	2162.509	0.414	2162.118	0.805
	2ν 3	A 1	2201.73284	2201.727	0.005	2201.98	-0.247
	ν 2 + ν 4	B 1		2633.045		2633.258
	ν 2 + ν 6	B 2		2685.004		2685.288
	ν 2 + ν 3	A 1		2790.027		2789.501
	RMS deviation				0.55		0.51
	RMS ( -level 2ν 4 ) **				0.21		0.34
	*) Experimental levels from refs

ACKNOWLEDGMENTS

The support of RFBR grant (Grant No. 19-03-00581) is acknowledged. V.T. acknowledges the support from Academic D.Mendeleev program of Tomsk State University and M.R. form LEFE CHAT CNRS program.

SUPPLEMENTARY MATERIAL

See supplementary material for the ab initio and the empirically optimized PES in orthogonal coordinates.