Cheikh T Bop

François Lique

Alexandre Faure

Ernesto Quintas-Sánchez

Richard Dawes

Non-LTE modelling of cyanoacetylene: evidence for isomer-specific excitation

Keywords: ISM, molecules -molecular data -molecular processes -ISM, abundances

Cyanoacetylene molecules are widespread in the interstellar medium (ISM) and several of its isomers have been detected in cold molecular clouds and circumstellar gas. Accurate estimates of the abundance ratio of cyanoacetylene isomers may provide deep insight into their environment. Such knowledge requires rigorous modeling of the emission spectra based on non-Local Thermodynamic Equilibrium (LTE) radiative transfer calculations. To this end, we computed excitation cross sections of HC 2 NC and HNC 3 induced by collision with para-and ortho-H 2 , using a quantum mechanical close-coupling method. Then, by thermally averaging these data, we derived rate coefficients for the first 31 low-lying rotational levels of each isomer for temperatures up to 80 K. For the para-H 2 collider, the propensity rules are in favor of rotational transitions involving Δ 1 = 2 for both isomers; while for the ortho-H 2 collider, Δ 1 = 2 and Δ 1 = 1 rotational transitions are favored for HC 2 NC and HNC 3 , respectively. A comparison of rate coefficients for the HC 3 N isomers shows differences up to an order of magnitude, especially at low temperatures. Finally, we performed non-LTE radiative transfer calculations to assess the impact of such variations in the analysis of observations. Our simulation suggests that the lack of collisional data specific to each isomer could lead to errors up to a factor of 2-3 in the excitation temperatures. We expect that these data could help in better understanding the cyanoacetylene chemistry and constraining the nitrogen chemistry in the ISM.

INTRODUCTION

The cyanopolyynes, namely HC 2 +1 N ( = 1, 2, 3, ...), are rod-like shaped nitriles which are made up of 2 carbon atoms bonded at the ends with a hydrogen atom (H) and a cyanide (CN) group. Since they usually have a low rotational constant and a relatively large dipole moment, which make them ideal barometers and thermometers of their astronomical environments (Cordiner et al. 2011(Cordiner et al. , 2013)), they are known as good gas-tracers. For instance, cyanoacetylene (HC 3 N), which has a rotational constant of ∼ 0.152 cm -1 (Creswell et al. 1977) and a dipole moment of ∼ 3.7 D (DeLeon & Muenter 1985), is found to be a good tracer for warm and dense regions (Costagliola & Aalto 2010). In addition, HC 3 N can also be used to understand the evolution of physical conditions of molecular clouds, ★ cheikhtidiane.bop@ucad.edu.sn † francois.lique@univ-lehavre.fr as recently done by Al-Edhari et al. (2017). Indeed, these authors used HC 3 N and HC 5 N emission spectra to unveil a crucial part of the IRAS 16293-2422 history, as they found that the source underwent a moderately rapid collapse and then a very rapid heating of the cold matter to 80 K. Cyanoacetylene is then of particular interest in astrochemistry. Like most cyanides, cyanoacetylene displays structural isomerism. So far, four stable isomers of HC 3 N-namely C 2 (H)CN, HNC 3 , HCNC 2 , and HC 2 NC-have been detected in laboratory experiments. Among them, only HNC 3 and HC 2 NC were observed in the interstellar and circumstellar gas (Kawaguchi et al. 1992a,b;Gensheimer 1997). The abundance ratio of isomers constitutes a useful tool to probe the physical and chemical conditions in the interstellar medium (ISM), as illustrated by Hacar et al. (2020). Therefore, accurately determining the abundance ratios between HC 3 N and its isomers may reveal important information about the media in which they are observed.

Very recently, determination of the abundance ratio of the HC 3 N isomers has been performed in several sources. In the L1544 prestellar core, Vastel et al. (2018) derived [HC 3 N] / [HC 2 NC] and [HC 3 N] / [HNC 3 ] abundance ratios of 28 -73 and 320 -800, respectively. In L483, Agúndez et al. (2019) derived the [HC 3 N] / [HC 2 NC] 74 and [HC 3 N] / [HNC 3 ] 840 abundance ratios. Cernicharo et al. (2020) derived, respectively, values of 77±8 and 442 ± 70 for these ratios in TMC-1, and values of 392 ± 22 and 1305 ± 45 in the IRC+10216 circumstellar envelope. All these studies concluded that HC 3 N is the most abundant isomer, followed by HC 2 NC and HNC 3 , respectively. However, the abundance ratio of these isomers varies greatly from one source to another. Vastel et al. (2018) used an astrochemical model to explain the abundance ratios derived from observations. They considered that the formation of HC 3 N isomers is likely to take place through dissociative recombination (DR) of HC 3 NH + : HC 3 NH + + e -→ HC 3 N + H

(1)

→ HC 2 NC + H (2) → HNC 3 + H (3)
Based on the branching ratios for the DR of DC 3 ND + measured by Vigren et al. (2012), they estimated ratios of 22 %, 4 %, and 22 % for the reactions (1), (2), and (3), respectively. HC 3 N is also formed through neutral-neutral reactions. This model reproduces fairly well the observed abundances of HC 3 N isomers in L1544, although it tends to overestimate the abundances of HC 2 NC, and especially HNC 3 , with respect to HC 3 N. The same conclusion was found by the astrochemical model of Osamura et al. (1999), performed for TMC-1. On one hand, such disagreements between astrochemical models and observations may indicate that different formation mechanisms could be involved. On the other hand, the abundance of HC 3 N isomers derived from the observations may also be uncertain because of the lack of available collisional data. Indeed, molecular abundances derivation should be performed using radiative transfer models beyond the Local Thermodynamic Equilibrium (LTE) approximation and therefore, requires the prior study of collisional excitation. Such studies currently exist only for the HC 3 N molecule.

The collisional excitation of HC 3 N by He and H 2 was first studied by Morris et al. (1976). Few time after, Green & Chapman (1978) determined HC 3 N-He rate coefficients using a potential energy surface (PES) obtained from an approximate 'electron gas' model combined with a Monte Carlo quasi-classical trajectory (QCT) method. The first quantum study was performed by Wernli et al. (2007); they investigated the scattering of HC 3 N induced by collisions with He and para-H 2 ( 2 = 0), 2 being the H 2 rotational momentum. These calculations were extended to the ortho-H 2 ( 2 = 1) collider by Faure et al. (2016).

As the HC 3 N isomers were not the object of any collisional excitation studies, rate coefficients for the HC 3 N molecule could only be used as a surrogate to model the emission spectra of its isomers (Cernicharo et al. 2020). It was recently shown (Hernández Vera et al. 2017) that isomerism effects can be important in collisional studies, and that rate coefficients can differ significantly for different isomers. Hence, we recently reported the first four dimensional (4D) PESs for the HC 2 NC-H 2 and HNC 3 -H 2 systems, as well as excitation rate coefficients due to para-H 2 ( 2 = 0) [hereafter p-H 2 ( 2 = 0)] collisions for temperature below 20 K (Bop et al. 2019). These data highlighted different excitation efficiencies for the HC 3 N isomers.

In this work, we extend our previous collisional excitation study of the HC 3 N isomers to higher temperatures and to the ortho-H 2 ( 2 = 1) [hereafter o-H 2 ( 2 = 1)] collider. Brief descriptions of the PESs and the quantum scattering methodology are presented in section 2. Rotational cross sections and rate coefficients of HC 2 NC and HNC 3 induced by collisions with both o-H 2 ( 2 = 1) and p-H 2 ( 2 = 0) are presented in section 3. In section 4 we study the excitation of the HC 3 N isomers in astrophysical media, and concluding remarks are given in section 5.

METHODS

HC 2 NC-H 2 and HNC 3 -H 2 potential energy surfaces

Using the CCSD(T)-F12b1 approach (Adler et al. 2007), along with the VTZ-F122 basis set (Hill et al. 2010), we derived 4D PESs for the HC 2 NC-H 2 and HNC 3 -H 2 van der Waals (vdW) systems within the rigid rotor approximation. As we have done in the past for other vdW linear dimers (Wang et al. 2016;Donoghue et al. 2016;Barclay et al. 2018;Desrousseaux et al. 2019;Castro-Juárez et al. 2019;[START_REF] Quintas-Sánchez | Physical Chemistry Chemical Physics Quintas[END_REF], the 4D PES analytical representation was constructed using an automated interpolating moving least squares methodology (Dawes et al. 2010;Majumder et al. 2016), implemented in the recently released software-package AUTOSURF (Quintas-Sánchez & Dawes 2019). The global estimated root-mean-squared fitting error tolerance was 0.33 cm -1 for both PESs, and the total number of symmetry unique points needed to reach that target was 6095 for HC 2 NC-H 2 and 4172 for HNC 3 -H 2 . As usual (Dawes & Quintas-Sánchez 2018), a local fit was expanded about each data point, and the final potential is obtained as the normalized weighted sum of the local fits. The MOLPRO electronic structure package (Werner et al. 2012) was used to perform all the ab initio calculations used to fit the PESs. These vdW complexes present each two minima; with well-depths of 174.1 cm -1 and 178.0 cm -1 for HC 2 NC-H 2 , and 270.8 cm -1 and 467.9 cm -1 for HNC 3 -H 2 . More details concerning the computational details, as well as the analytical representation, well-depths and its precise locations for each PESs, can be obtained in Bop et al. (2019).

In order to perform quantum scattering calculations, the angular dependence of the PES was refitted over spherical harmonics. Due to the high-anisotropy, the regularization procedure of Wernli et al. (2007) was employed. The final expansion included 380 angular basis functions and each expansion coefficient was interpolated radially in the range = 5-30 a 0 with a smooth connection to standard short-and long-range extrapolations. Full details can be found in Bop et al. (2019).

Scattering calculations

In this work, we focus on the rotational transitions of HC 2 NC and HNC 3 induced by collision with H 2 molecules:

HC 2 NC( 1 ) + H 2 ( 2 ) → HC 2 NC( 1 ) + H 2 ( 2 ) (4) HNC 3 ( 1 ) + H 2 ( 2 ) → HNC 3 ( 1 ) + H 2 ( 2 ) (5)
where 1 and 2 stand for the rotational angular momenta of the HC 3 N isomers and H 2 , respectively. Separate calculations were done for p-H 2 (even 2 ) and o-H 2 (odd 2 ) since ortho-para-H 2 conversion is not allowed during inelastic collisions. The splitting of the HC 2 NC and HNC 3 rotational lines due to the nuclear spin of the N atom is not taken into account here as hyperfine transitions of HC 2 NC and HNC 3 are currently not resolved in astrophysical media. Hyperfine resolved rate coefficients can be approximately derived from our data, in case they are needed, using the procedure of Faure & Lique (2012) and Lanza & Lique (2014).

Using the 4D HC 2 NC-H 2 and HNC 3 -H 2 PESs, we derived inelastic cross sections for the 31 first low-lying rotational levels ( 1 = 0-30) of HC 2 NC and HNC 3 . The energy levels of the HC 2 NC, HNC 3 and H 2 molecules were determined using the experimental rotational and centrifugal distortion constants, 0 and 0 . For HC 2 NC and HNC 3 , the two employed constants ( 0 , 0 ) are (0.166 cm -1 , 2.091 × 10 -8 cm -1 ) and (0.156 cm -1 , 2.063 × 10 -8 cm -1 ), respectively; while for H 2 , 0 = 59.322 cm -1 and 0 = 0.047 cm -1 (Huber & Herzberg 1979). The maximum total energy used in the scattering calculations was set to 500 cm -1 and 620 cm -1 for p-and o-H 2 calculations (The ground o-H 2 ( 2 = 1) level is at about 118 cm -1 above the p-H 2 ( 2 = 0) level). We used the exact close-coupling (CC) quantum mechanical method (Green 1975) as implemented in the MOLSCAT computer code (Hutson & Green 1994). The modified log derivative-airy integrator (Alexander & Manolopoulos 1987) was used to solve the coupled equations.

For both of the HC 2 NC-H 2 and HNC 3 -H 2 collision systems, for both p-and o-H 2 calculations, the close-coupling equations were propagated from from 2 Å, to at least 25 Å. In order to ensure convergence of the state-to-state inelastic cross sections ( 1 2 → 1 2

), the rotational bases of the HC 2 NC, HNC 3 and H 2 molecules were adjusted as a function of the total energy. For instance, at the highest total energy, 1 of HC 2 NC was fixed at 33 and 34 for p-H 2 and o-H 2 collisions, respectively; while 1 of HNC 3 remained at 34 for collisions with both spin-isomers (p-H 2 and o-H 2 ). The cross sections obtained with such rotational basis of the cyanopolyynes differ by less than 1% compared to cross sections obtained with rotational basis of the cyanopolyynes containing 41 levels (e.g. increasing 1 up to 40).

Only the lowest rotation levels of p-H 2 ( 2 = 0) and o-H 2 ( 2 = 1) were included in the calculations. Table 1 shows the impact of the p-H 2 rotational basis on the magnitude of the cross sections as a function of the energy. Two series of test were performed, the first one was for one single total angular momentum ( = 0) and the second test consists of summing the partial cross sections over the first 11 total angular momenta ( = 0 -10). We found that the convergence increases when summing over the total angular momentum. For example, in the cases of HC 2 NC-p-H 2 , summing partial cross-sections up to = 10 leads to average biases less than ∼10% for all considered energies whereas the difference could be greater than a factor 2 for = 0 only. For HNC 3 -p-H 2 and HC 2 NC-p-H 2 , we found that the convergence with respect to the H 2 basis increases with the increase of the energy and is relatively independent to the energy, respectively.

For HC 2 NC-p-H 2 , the global convergence of the results considering the minimal H 2 basis is good enough, i.e. better than 10% for all the energy grid considered here. The difference can be higher for state-to-state data but remains below 20-30% for all transitions. For HNC 3 -p-H 2 , because of the larger well depth in the PES, the deviation is greater but less than 20% beyond 100 cm -1 . The large bias (> 50%) obtained at 50 cm -1 is due the resonances which are significant because of the deep potential well of this collision system.

From this table, we can see that our cross sections are globally Rate coefficients (cm 3 s Rate coefficients (cm 3 s converged to better than 15-20%. In the case of the o-H 2 collider, setting j 2 = 1 -3 instead of 2 = 1 did not affect much the cross sections (less than 10% of deviation) so that the convergence of the o-H 2 basis is estimated to be better than that of the p-H 2 ones. We smoothly scanned the total energy range up to 250 cm -1 (370 cm -1 ) for collisions with p-H 2 (o-H 2 ), with a fine grid step of 0.2 cm -1 . Thereafter, we gradually increased this step up to 20 cm -1 , which was used in the total energy interval 400-500 cm -1 (520-620 cm -1 ) for the HC 2 NC-p-H 2 and HNC 3 -p-H 2 (HC 2 NCo-H 2 and HNC 3 -o-H 2 ) collisional systems. Note that the rotational energy at 1 = 30 for HC 2 NC and HNC 3 is 144.14 cm -1 and 144.80 cm -1 , respectively. Hence, the kinetic energy range used in this work is essentially the same for both isomers.

-1 ) 1, 0 → 0, 0 2, 1, → 1, 1 2, 0 → 0, 0 2, 0, → 1, 0 2, 1 → 0, 1 1, 1 → 0, 1 (a)
-1 ) 1, 1 → 0, 1 2, 1, → 1, 1 1, 0 → 0, 0 2, 0, → 1, 0 2, 1 → 0, 1 2, 0 → 0, 0 (b)
From the computed cross sections, we derive the corresponding rotational rate coefficients ( 1 2 → 1 2

) using the Maxwell-Boltzmann thermal average over the kinetic energy:

1 2 → 1 2 ( ) = ( 8 ) 1/2 2 ∫ ∞ 0 1 2 → 1 2 ( ) - (6)
In equation 6, = ( × ) -1 , with being the Boltzmann con-Table 1. The effect of the H 2 rotational manifold in the cross-section ( ˚ 2 ) calculations for selected energies and total angular momenta ( ). The values in brackets represent the relative errors (in percent) when restricting 2 to zero and < Δ > stands for the average deviation (in percent) over all inelastic transitions. stant and being the reduced mass of the collisional system. Since we focus on the 31 low-lying rotational levels of the complexes of interest, the kinetic energy range used in the calculations permitted computation of rate coefficients for temperatures up to 80 K.

HC 2 NC-para-H 2 = 50 cm -1 = 100 cm -1 = 500 cm -1 = 0 = 0 -10 = 0 = 0 -10 = 0 = 0 -10 1 1 2 = 0 2 = 0 -2 2 = 0 2 = 0 -2 2 = 0 2 = 0 -2 2 = 0 2 = 0 -2 2 = 0 2 = 0 -2 2 = 0 2 = 0 -2 1 0 8.03E-3 1.04E-2 6.36E-1 9.10E

ROTATIONAL CROSS SECTIONS AND RATE COEFFICIENTS

We display in Fig. 1 the temperature variation of the HC 2 NC-H 2 and HNC 3 -H 2 rate coefficients. Comparing the plots in the two panels of Fig. 1 reveals that rate coefficients of HNC 3 -H 2 are globally larger than those of HC 2 NC-H 2 . Such behavior can be related to the fact that the minimum of the HNC 3 -H 2 interaction potential is ∼ 2.5 times deeper than that of HC 2 NC-H 2 . It is interesting to note that, for both isomers, the temperature variation of the rate coefficients is weak, the rate coefficients being almost constant for temperatures above ∼ 20 K. We also observe that the rate coefficients for o-H 2 collisions are larger than those for p-H 2 collisions. The difference is however moderate (within a factor of three), as previously found for HC 3 N-H 2 collisions (Faure et al. 2016).

To gain further insight into propensity rules, we plot in Fig. 2 the de-excitation rate coefficients of the HC 3 N isomers out of the 1 = 10 level at 50 K. The HC 3 N rate coefficients were taken from Faure et al. (2016). For collisions with p-H 2 , for both isomers, one can see that the propensity in favor of even Δ 1 transitions observed in our previous work (Bop et al. 2019) persists at higher temperatures. In contrast, for collisions with o-H 2 , the propensity rules are different for the two isomers. A propensity rule in favor of even Δ 1 transitions is found for HC 2 NC (however, less marked than for p-H 2 collisions) while for HNC 3 , the magnitude of rate coefficients typically follows an energy-gap law, (i.e., Δ 1 = 1 > Δ 1 = 2 > Δ 1 = 3 ...). The change of propensity rules with respect to the rotational states of H 2 in the case of HNC 3 -H 2 collisions, was also found for the HC 3 N-H 2 and HCN-H 2 collisional systems (Faure et al. 2016;Hernández Vera et al. 2017). Such behavior originates from the quadrupole moment of H 2 ( 2 > 0) which modifies the interaction with the target. Indeed, the dipole-quadrupole interaction governs the long-range interactions but vanishes for H 2 ( 2 = 0). The difference between o-and p-H 2 is marked in the case of HNC 3 since HNC 3 exhibits a strong permanent dipole moment of 5.665 D (Botschwina et al. 1992). Such inversion of the propensity rules is not seen for HC 2 NC because of the weaker dipole moment (2.93 ), (Krüger et al. 1991).

In order to check the validity of using the rate coefficients of one isomer to derive the abundance of the other isomers, we depict in Fig. 3 the rate coefficients of HC 2 NC and HNC 3 as a function of those of HC 3 N and the rate coefficients of HNC 3 as a function of those of HC 2 NC. For each isomer, the relative deviation to others is less than a factor of 3 when considering the o-H 2 collider. The deviation is even less than a factor of 1.5-2 when considering only the dominant rate coefficients ( ( ) ≥ 10 -11 cm 3 s -1 ).

On the other hand, larger differences exist between the excitation of the different isomers induced by collisions with p-H 2 . Indeed, the difference can be up to an order of magnitude for the weaker rate coefficients and up to a factor 3 for the dominant ones, regardless of temperature. It is interesting to note that the rate coefficients of HC 2 NC and HNC 3 globally outweigh those of HC 3 N whereas a slight dominance of HNC 3 with respect to HC 2 NC is found. The stronger magnitude of the HNC 3 -H 2 rate coefficients can be related to the larger well-depth in the interaction potential compared to the two other isomers since all isomers exhibit a similar rotational structure.

9,1 8,1 7,1 6,1 5,1 4,1 3,1 2,1 1,1 0,1 j ′ 1 ,j ′ 2 0.0 0.5 1.0 1.5 2.0 2.5
Rate coefficients (cm 3 s -1 )

×10 -10 j 1 , j 2 = 10, 1 (a)
HC 2 NC HNC 3 HC 3 N 9,0 8,0 7,0 6,0 5,0 4,0 3,0 2,0 1,0 0,0 j ′ 1 ,j ′

THE EXCITATION OF HC 3 N, HC 2 NC, AND HNC 3 IN ASTROPHYSICAL MEDIA

The comparison of rate coefficients made in the above paragraph showed that HC 3 N and its isomers have different excitation behaviours when they collide with molecular hydrogen. It is then of significant interest to simulate the impact of using the appropriate collisional data when deriving the abundance of the HC 3 N isomers in astrophysical media. Non-LTE radiative transfer calculations were performed with the RADEX code ( Van der Tak et al. 2007), using the escape probability approximation. Our collisional data were combined with spectroscopic data from the Jet Propulsion Laboratory (JPL) catalog, as downloaded from the Cologne Database for Molecular Spectroscopy (CDMS) portal (Müller et al. 2005) 3 . In practice, we calculated the excitation and brightness temperatures ( and , respectively) for the 1 = 5 → 4, 8 → 7 and 11 → 10 emission lines, that are frequently observed in the interstellar and circumstellar gas (Vastel et al. 2018;Agúndez et al. 2019;Cernicharo et al. 2020). In our simulation, the cosmic microwave background used as 3 https://cdms.astro.uni-koeln.de/cdms/portal/ background radiation field was set to 2.73 K and the line width Δ to 1 km s -1 . The kinetic temperature was fixed at 10 K to simulate the physical conditions of cold molecular clouds such as TMC-1 or L1544, and at 50 K to reach out the temperature of the outer layers of circumstellar envelopes such as IRC+10216. Furthermore, we chose, for the ortho-to-para ratio of H 2 , values of 0 (only p-H 2 ; model 1); and 3, corresponding to the high temperature limit (model 2). However, any o-and p-H 2 proportion can be used in models since we provided state-to-state rate coefficients for both o-and p-H 2 colliders. In practice, the total density of molecular hydrogen (H 2 ) was spanned from 10 2 to 10 8 cm -3 and the column densities of HC 3 N, HC 2 NC, and HNC 3 was fixed at = 10 12 cm -2 . This value is close to the HC 2 NC and HNC 3 column density derived in cold molecular clouds (Vastel et al. 2018). 4 Figs. 4, 5, and6 show the excitation temperature obtained in our simulations as a function of (H 2 ).

First, we note a supra-thermal excitation for the 1 = 5 → 4 line at 50 K for densities in the 10 4 -10 5 cm -3 range. Indeed, the excitation temperature ( ) is larger than the kinetic temperature ( ). Then, as expected, for all isomers the excitation temperature varies from ∼ 2.73 K to the kinetic temperature (10 K or 50 K). At low H 2 density, the medium is dilute and the excitation temperature is reduced to the cosmic background temperature ( ≈ CMB = 2.73 K) while at high H 2 density, tends asymptotically to the adopted kinetic temperature [ = 10 K for panels (a), (b), and (c) and = 50 K for panels (d), (e), and (f)]. An agreement between the and the kinetic temperature indicates that LTE conditions are fulfilled. According to our calculations, this is true when the H 2 density exceeds 5×10 4 to 10 6 cm -3 depending on the isomers, the oto-p-H 2 ratios, the emission lines, and the kinetic temperatures. The density at which = is called the LTE density below. Because of the weak temperature dependence of the rate coefficients, the LTE densities are not significantly influenced by the kinetic temperature. The o-to-p-H 2 ratio also slightly impact the LTE densities because of the moderate difference between o-and p-H 2 rate coefficients. At the opposite, because of the significant variation of the Einstein coefficients with the increasing rotational state, the LTE density changes by up to an order of magnitude when changing from 1 = 5 → 4 to 1 = 11 → 10. For all isomers, the LTE density is found to be generally larger than the typical densities in both cold molecular clouds and circumstellar envelopes indicating clearly that non-LTE modeling of the emission spectra is mandatory.

We then focus in the intermediate region where radiative and collisional processes are in competition. The corresponding density is called the critical density. The comparison of the of the different isomers reveals that the excitation conditions of HC 3 N, HC 2 NC and HNC 3 are different for all lines. In particular, we found that (HC 2 NC) > (HC 3 N) > (HNC 3 )] and the differences for the can be up to a factor 2-3. This clearly demonstrates that the excitation conditions of the different isomers are quite different and that it is very risky to assume similar for all isomers. In Table 2, we show the 1 = 8 → 7 brightness temperature ratios of the three HC 3 N isomers ( [HC 2 NC]/ [HC 3 N],

[HNC 3 ]/ [HC 3 N], and [HNC 3 ]/ [HC 2 NC]) for two H 2 density and two temperatures.

At high H 2 density (∼ 10 8 cm -3 ), when LTE conditions are fulfilled, the brightness temperature ratios are not equal to unity due

p-H 2 X = HC 2 NC Y = HC 3 N p-H 2 X = HNC 3 Y = HC 3 N p-H 2 X = HNC 3 Y = HC 2 NC o-H 2 X = HC 2 NC Y = HC 3 N o-H 2 X = HNC 3 Y = HC 3 N o-H 2 X = HNC 3 Y = HC 2 NC (a) (b) (c) (d) (e) (f)
Rate coefficients of X Rate coefficients of Y Table 2. Comparison of the brightness temperature ratios and the ratios of the square of the electric dipole moment. The brightness temperatures are calculated at 10 and 50 K. For each line, the first entries and the second entries (values in brackets) were obtained at = 10 8 cm -3 and = 3 × 10 4 cm -3 , respectively. ratios of ratios of 2 (Krüger et al. 1991) dipole moment of HNC 3 (Botschwina et al. 1992) dipole moment of HC 3 N (DeLeon & Muenter 1985) to the difference between the dipole moments of the HC 3 N isomers. Indeed, the Einstein coefficients associated to the radiative transitions are proportional to the square of the electric dipole moment ( 2 ). We compare in Table 2 the ratios obtained at = 10 8 cm -3 to the squared ratios of the dipole moments. For both kinetic temperatures, the agreement is very good, the minor deviations coming from the different frequencies of the cyanopolyynes isomers. In Table 2, we also provide the ratios at = 3 × 10 4 cm -3 , the typical H 2 density in molecular clouds. As one can see, the deviation compared to the squared ratios of the dipole moments is large (up to a factor 2) demonstrating again the importance of non-LTE effects.

Finally, weak maser emissions have been observed for the 1 = 1 → 0 (Hunt et al. 1999) and 1 = 4 → 3 lines (Ellingsen et al. 2017) of HC 3 N. It is then of interest to see if such maser action is predicted for HC 3 N isomers. We calculated the excitation temperatures for the 1 = 1 → 0 and 4 → 3 line at different temperatures. Below 40 K, no population inversion is predicted for the 4 → 3 line, whereas population inversion occurs even at 10 K for the 1 → 0 line. In Fig. 7, we computed the logarithm of the excitation temperature (see caption) for the 1 = 4 → 3 assuming a thermal o-to p-H 2 ratio of ∼ 0.3. As can be seen in the figure, population inversion is present for all isomers, depending on the H 2 density. The inversion occurs at larger H 2 densities for HNC 3 than for HC 3 N and HC 2 NC because of stronger Einstein coefficients. Hence (weak) maser emissions may help to estimate the H 2 density. It is also interesting to note that recent abundance ratio determinations (Cernicharo et al. 2020) were based on the 1 = 4 → 3 emission line so that the analysis of such emission spectra in the absence of collisional isomer-specific rate coefficients might be quite uncertain.

CONCLUSION

State-to-state excitation cross-sections of HC 2 NC and HNC 3 induced by collision with p-H 2 and o-H 2 have been calculated using a close-coupling method. These results were then thermally averaged using a Maxwell-Boltzmann velocity distribution to determine (de)excitation rate coefficients involving the first 31 low-lying rotational levels of both HC 3 N isomers for temperatures up to 80 K.

The computed data exhibit some interesting propensity rules. [panels (d), (e), and (f)] and a fixed column density of 10 12 cm -2 for all isomers. For model 1, only p-H 2 is considered, while in model 2 oand p-H 2 are in the 3:1 proportion. and HNC 3 colliding with p-H 2 ; while for collisions with o-H 2 , the propensity is found in favor of Δ 1 = 2 and Δ 1 = 1 transitions for HC 2 NC and HNC 3 targets, respectively. We found that the HC 2 NC and HNC 3 data globally agree within a factor of 3 with those computed for HC 3 N by Faure et al. (2016), despite larger differences between some state-to-state transitions. Globally, the HNC 3 -H 2 are found to be larger than those of the other isomers.

(a) (b) (c) (d) (e) (f) n (cm -3 ) T ex (K)
We also performed radiative transfer calculations for the three cyanoppolyynes isomers using the RADEX code. The comparison of excitation temperatures showed that the excitation of the three isomers is different and put in evidence the need for specific analysis for the three isomers. In particular, we found that the excitation temperatures of HNC 3 are stronger than that of the other two and clearly demonstrates that special care should be taken when deriving its abundance from the observational spectra. We anticipate that the abundance of HNC 3 derived from the observation may be reduced by the use of an accurate non-LTE model, leading to an even stronger disagreement with astrochemical models. Numerous studies on excitation of interstellar isomers such as HCN/HNC (Sarrasin et al. 2010;Hernández Vera et al. 2017), AlCN/AlNC (Hernández Vera et al. 2013), MgCN/MgNC (Hernández Vera et al. 2013), SiCN/SiNC (Hernández Vera et al. 2015) and HCO + /HOC + (Santander et al. 2019) have been reported over the last decades. In all these studies, differences ranging from few tenths of a percent up to one order of magnitude were found. It is then clear that ideally, excitation conditions of all isomers should be studied individually, and we strongly recommend to compute specific rate coefficients of the astrophysical isomers in order to accurately determine their abundance ratios. This will lead to better insights about their mechanism of formation in the ISM and will put more robust constraints on astrochemical models. Among the interstellar isomers, the cyanopolyynes are of great interest, and a possible target for future excitation studies could be the HC 4 NC molecule, an isomer of the ubiquitous HC 5 N that was recently detected in TMC-1 (Xue et al. 2020).

Figure 1 .

 1 Figure 1. Rotational rate coefficients of HC 2 NC [panel (a)] and HNC 3 [panel (b)] induced by collision with p-H 2 ( 2 = 0) and o-H 2 ( 2 = 1). In both panels, dashed lines stand for the results obtained in the p-H 2 ladder, while solid lines refer to the o-H 2 symmetry. The transitions are labeled as 1 , 2 → 1 , 2 .

Figure 2 .

 2 Figure 2. Downward rate coefficients of HC 2 NC, HNC 3 and HC 3 N (Faure et al. 2016) induced by collision with o-H 2 [panel (a)] and p-H 2 [panel (b)] at T = 50 K. The transitions are labeled as 10, 1 → 1 , 1 and 10, 0 → 1 , 0 in panels (a) and (b), respectively.

Figure 3 .

 3 Figure 3. Comparison of rate coefficients of HC 2 NC [panels (a) and (d)] and HNC 3 [panels (b) and (e)] to those of HC 3 N as well as of HC 2 NC to those of HNC 3 [panels (c) and (f)], at T = 10 K and 50 K. The collisions treated in the p-H 2 [o-H 2 ] ladder are represented in the panels (a), (b) and (c) [(d), (e), and (f)]. In all panels, the diagonal lines are distant from each other by a factor of three and the solid one stands for y = x.

Figure 4 .

 4 Figure 4. Excitation temperature for the 1 = 5 → 4 line of HC 3 N, HC 2 NC and HNC 3 calculated for temperatures = 10 K [panels (a), (b), and (c)] and = 50 K [panels (d), (e), and (f)] and a fixed column density of 10 12 cm -2 for all isomers. For model 1, only p-H 2 is considered, while in model 2 oand p-H 2 are in the 3:1 proportion.

Figure 5 .

 5 Figure 5. Same as Fig. 4 but for the 1 = 8 → 7 line.

Figure 6 .

 6 Figure 6. Same as Fig. 4 but for the 1 = 11 → 10 line.

Figure 7 .

 7 Figure 7. Excitation temperature for the 1 = 4 → 3 line of HC 3 N, HC 2 NC and HNC 3 calculated for a temperature of 50 K and a fixed column density of 10 12 cm -2 , for all isomers.

  Transitions involving even Δ 1 values are favored for both HC 2 NC

					model 1: 100% of p-H 2
		10 1.0
		8	(a)				(b)			(c)
		0.8
		6
	T ex (K)	4 150 0.4 0.6	(d)	HC 2 NC HC 3 N			(e)	HNC 3 HC 3 N	(f)	HC 2 NC HNC 3
		100
		0.2
		50
		10 2 0 0.0 0.0	10 4	10 6 0.2	10 8	10 2	10 4 n (cm -3 ) 10 6 10 8 10 2 0.4 0.6	10 4 0.8	10 6	10 8 1.0
				model 2: 75% of o-H 2 and 25% of p-H 2
		10 1.0
		8
		0.8
		6
		4 0.6		HC 2 NC HC 3 N					HNC 3 HC 3 N			HC 2 NC HNC 3
		150 0.4
		100
		0.2
		50
		10 2 0 0.0 0.0	10 4	10 6 0.2	10 8	10 2	0.4	10 4	10 6	0.6	10 8 10 2	10 4 0.8	10 6	10 8 1.0

© 2020 The Authors

explicitly correlated coupled cluster with single, double and non-iterative triple excitation

polarized valence triple zeta Gaussian basis set which explicitly treats the correlation of electrons MNRAS 000, 1-?? (2020)

MNRAS 000, 1-?? (2020)

An increase of the column density by one or two orders of magnitude does not change the excitation temperature , since the opacity is still low ( << 1).MNRAS 000, 1-?? (2020)

ACKNOWLEDGMENTS

The authors acknowledge the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of Centre National de la Recherche Scientifique (CNRS)/Institut National des Sciences de l'Univers (INSU) with Institut de Chimie (INC)/Institut de Physique (INP) co-funded by Commissariat a l'Energie Atomique (CEA) and Centre National d'Etudes Spatiales (CNES). FL acknowledges financial support from the Institut Universitaire de France and the French National Research Agency (ANR) through a grant to the Anion Cos Chem Project (ANR-14-CE33-0013). R.D. and E.Q-S. are supported by the US Department of Energy Office of Science, Office of Basic Energy Sciences (Award DE-SC0019740). This work was performed using HPC resources from GENCI-CINES (Grant No. A0070411036), as well as the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material.


H D Cornean
email: cornean@math.aau.dk

D Krejčiřík

T G Pedersen

AND N Raymond

E Stockmeyer

ON THE TWO-DIMENSIONAL QUANTUM CONFINED STARK EFFECT IN STRONG ELECTRIC FIELDS

We consider a Stark Hamiltonian on a two-dimensional bounded domain with Dirichlet boundary conditions. In the strong electric field limit we derive, under certain local convexity conditions, a three-term asymptotic expansion of the low-lying eigenvalues. This shows that the excitation frequencies are proportional to the square root of the boundary curvature at a certain point determined by the direction of the electric field.

1. Introduction 1.1. Physical motivation. There has been a lot of theoretical and experimental work on the properties of quantum confined semiconductor devices. These systems exhibit interesting features that may find applications in nano-technology (see e.g. [START_REF] Harrison | Quantum wells, wires and dots: theoretical and computational physics of semiconductor nanostructures[END_REF][START_REF] Harutyunyan | Cylindrical nanolayer in the strong uniform electrical field: The field localization of carriers and electrooptical transitions[END_REF]). Physically, the confinement may be achieved through an interface with the vacuum or an hetero-junction, i.e., by interfacing the semiconductor with an isolator or with another semiconductor of larger gap [START_REF] Harrison | Quantum wells, wires and dots: theoretical and computational physics of semiconductor nanostructures[END_REF]. The most prominent semiconductor devices are the socalled quantum wells, quantum wires and quantum dots, where the material is confined in one, two and three orthogonal directions, respectively. In the presence of symmetry along the non-confined directions the description of the system is reduced to the study of a Hamiltonian in dimension one, for quantum wells, and in dimension two for quantum wires (see e.g., [START_REF] Harrison | Quantum wells, wires and dots: theoretical and computational physics of semiconductor nanostructures[END_REF]Chapter 8]).

An interesting phenomenon is the behavior of the energy levels of the system when a uniform electric field is applied; this is known as the quantum confined Stark effect [START_REF] Miller | Bandedge electroabsorption in quantum well structures: The quantum-confined Stark effect[END_REF] (see also [START_REF] Pedersen | Stark effect in spherical quantum dots[END_REF] and references therein). In a simplified model one may consider independently electrons or holes. Then, in the so-called effective mass envelope function approximation [START_REF] Miller | Electroabsorption of highly confined systems: Theory of the quantum-confined Franz-Keldysh effect in semiconductor quantum wires and dots[END_REF], the effective Hamiltonian describing the quantum confined Stark effect is formally given by

´ 2 2m ∆ `q F ¨x `Vconf , (1.1)
where is Plank's constant divided by 2π, m ą 0 is the effective mass of the electron (or hole), ´q is the charge of the particle (electron or hole), F is the electric field, and V conf is some confining potential. In the simplest case one models the confining potential as infinite potential walls, i.e., one considers the first two terms of the Hamiltonian in (1.1) restricted to the domain of confinement with Dirichlet boundary conditions. The quantum confined Stark effect modeled as in (1.1) with Dirichlet boundary conditions have been considered, for instance, in [START_REF] Bastard | Variational calculations on a quantum well in an electric field[END_REF][START_REF] Miller | Bandedge electroabsorption in quantum well structures: The quantum-confined Stark effect[END_REF] for quantum wells, in [START_REF] Sukumar | Polarizability of a carrier in an isolated well of a quantumwell wire[END_REF][START_REF] Wang | Transverse Stark effect in the optical absorption in a square semiconducting quantum wire[END_REF][START_REF] Yi | Transverse Stark effect in a rectangular semiconducting quantum wire[END_REF][START_REF] Vázquez | Transverse Stark effect of electrons in a semiconducting quantum wire[END_REF][START_REF] Lee | Stark effect in the optical absorption in quantum wires[END_REF][START_REF] Harutyunyan | Cylindrical nanolayer in the strong uniform electrical field: The field localization of carriers and electrooptical transitions[END_REF] for quantum wires, and in [START_REF] Wei | Stark effect of electrons in semiconducting rectangular quantum boxes[END_REF][START_REF] Pedersen | Stark effect in spherical quantum dots[END_REF] for quantum dots.

In this work we are interested in the study of the low-lying eigenvalues of the following two-dimensional Hamiltonian restricted to an open set Ω Ă R 2 :

H " ´ 2 2m pB 2 x `B2 y q `qF x , (1.2)
acting on a dense subspace of the square integral functions L 2 pΩq with Dirichlet boundary conditions. This is a model Hamiltonian to describe the energy levels of a quantum wire with cross section Ω in the presence of an electric field perpendicular to the non-confined direction. Notice that we have chosen coordinates such that F is parallel to the x-axis.

The low-lying eigenvalues of this operator have been studied, partially numerically, for different geometries in several papers such as in [START_REF] Sukumar | Polarizability of a carrier in an isolated well of a quantumwell wire[END_REF][START_REF] Wang | Transverse Stark effect in the optical absorption in a square semiconducting quantum wire[END_REF] for squares, in [START_REF] Yi | Transverse Stark effect in a rectangular semiconducting quantum wire[END_REF] for rectangles, in [START_REF] Vázquez | Transverse Stark effect of electrons in a semiconducting quantum wire[END_REF][START_REF] Lee | Stark effect in the optical absorption in quantum wires[END_REF][START_REF] Pedersen | Stark effect in finite-barrier quantum wells, wires, and dots[END_REF] for disks, and in [START_REF] Harutyunyan | Cylindrical nanolayer in the strong uniform electrical field: The field localization of carriers and electrooptical transitions[END_REF] for annuli.

From now on we will work with domains satisfying the following conditions, see Figure 1.

Assumption 1.1. The set Ω is open, bounded, and connected. We assume that there exists a unique point A 0 P BΩ such that the first component of A 0 is given by

x min " inf px,yqPΩ
x " min px,yqPΩ

x .

We also assume that BΩ is smooth near A 0 , and that the curvature at A 0 , denoted by κ 0 , is positive.

Let us describe the notion of curvature we use in this paper. Consider a smooth counterclockwise parametrization s Þ Ñ γpsq by arc-length of the boundary near A 0 " γp0q. If npsq is the outward pointing normal to BΩ at γpsq, the curvature κpsq at γpsq is defined through the relation γ 2 psq " ´κpsqnpsq. We have κp0q " κ 0 ą 0. Then we can provide a three terms asymptotic expansion for the individual eigenvalues of H in the limit of strong electric fields. Our main result, Theorem 1.2, implies that for qF ą 0 the n-th eigenvalue of H behaves in the limit of strong electric field as

E n " qF x min `ˆpqF q 2 2m ˙1 3 z 1 `pn ´1 2 q ˆqF κ 0 m ˙1 2 `OpF 1{3 q , (1.3)
where z 1 « 2.338 is the absolute value of the smallest zero of the Airy function.

In particular, we find an interesting behaviour of the energy splitting in terms of the geometry of Ω which is proportional to the F 1{2 . Such an equidistant splitting for the eigenvalues in the strong electric field regime has been observed in [START_REF] Lee | Stark effect in the optical absorption in quantum wires[END_REF] for disk shaped Ω using partially numerical methods (see Figure 7 from [START_REF] Lee | Stark effect in the optical absorption in quantum wires[END_REF]). However, its dependence on the curvature does not seem to have been reported before. Heuristically the third term in the expansion may be explained as follows: under a strong electric field the particle is pushed towards x min and behaves as an harmonic oscillator in the direction perpendicular to the field with elastic constant proportional to the curvature at x min . Let us mention here that a physics-oriented paper is in preparation. Our first investigations in this direction suggest that the spectral splitting (1.3) might be experimentally accessible.

Main result.

Let Ω be as in Assumption 1.1. Notice that by factorizing qF ą 0 in the expression of H in (1.2) we have L h :" pqF q ´1H " ´h2 pB 2 x `B2 y q `x , h " ? 2mqF .

(1.4)

We define L h as the unique self-adjoint operator defined through the quadratic form

Q h pϕq " h 2 ż Ω |∇ϕpxq| 2 dx `żΩ x |ϕpxq| 2 dx , ϕ P H 1 0 pΩq .
The operator L h has domain contained in H 1 0 pΩq and acts as in (1.4). This is the Dirichlet realization of the Stark Hamiltonian. We want to describe the first eigenvalues of this operator in the limit h Ñ 0. We obtain (1.3) from the following result. Denote by pλ n phqq ně1 the eigenvalues of L h in increasing order, where each eigenvalue is repeated according to its multiplicity. Theorem 1.2. Let n P t1, 2, . . . u. Then, we have as h Ñ 0

λ n phq " x min `z1 h 2 3 `hp2n ´1q c κ 0 2 `Oph 4 3 q .
Remark 1.3. Various extensions of our main theorem can be considered.

i. It is possible to prove asymptotic expansions in powers of h 1{6 of the low-lying eigenvalues by using a formal series analysis. ii. In our generic geometric situation, we could even prove that the eigenfunctions admit WKB expansions. iii. Our strategy may be adapted to deal with a finite number of non-degenerate minima, and it would even be possible to investigate tunnel effects when these minima have symmetries.

The proofs of such extensions can be adapted from a recent literature developed in the context of the Born-Oppenheimer approximation (see for instance [START_REF] Raymond | Bound states of the magnetic Schrödinger operator[END_REF]Chapter 11], or the generalizations [START_REF]Formules de Weyl par réduction de dimension[END_REF][START_REF] Martinez | A general effective Hamiltonian method[END_REF], and also [START_REF] Helffer | Tunneling for the Robin Laplacian in smooth planar domains[END_REF] where tunnelling estimates are also considered).

Remark 1.4. If we replace the potential x by ix in the Hamiltonian and one looks at the real part of the eigenvalues the second term in the asymptotic expansion has a factor 1{2 (see e.g., [START_REF] Grebenkov | The complex Airy operator on the line with a semipermeable barrier[END_REF][START_REF] Henry | On the semi-classical analysis of Schrödinger operators with purely imaginary electric potentials in a bounded domain[END_REF][START_REF] Almog | Spectral analysis of a complex Schrödinger operator in the semiclassical limit[END_REF][START_REF] Almog | The stability of the normal state of superconductors in the presence of electric currents[END_REF]). Concerning the analysis in the case of imaginary electric fields, the reader might also want to consider [START_REF] Grebenkov | On spectral properties of the Bloch-Torrey operator in two dimensions[END_REF][START_REF] Almog | Spectral semi-classical analysis of a complex Schrödinger operator in exterior domains[END_REF][START_REF] Almog | On a Schrödinger operator with a purely imaginary potential in the semiclassical limit[END_REF].

The rest of this article is organized as follows: In Section 2 we show that low energy eigenfunctions and its derivatives are exponentially well localized around x min . We use this to find an effective Hamiltonian Ă M h , whose low energy eigenvalues are those of H modulo an exponentially small error, this is done in Section 3. This operator is expressed in tubular coordinates and acts on a tiny domain around x min with Dirichlet boundary conditions. In the last section we provide asymptotic upper and lower bounds for the eigenvalues of Ă M h .

Localization near the potential minimum

The following proposition states that the eigenfunctions associated with the low-lying eigenvalues are localized in x near x min . Proposition 2.1. Let M ą 0. There exist ε, C, h 0 ą 0 such that, for all h P p0, h 0 q, for all eigenvalues λ such that λ ď x min `M h (2.2)

Proof. We write the Agmon formula, for all ψ P DompL h q and all bounded Lipschitz functions Φ:

xL h ψ, e 2Φ{h ψy " Q h pe Φ{h ψq ´}e Φ{h ∇Φψ} 2 .
Let ψ be an eigenfunction corresponding to an eigenvalue λ ď x min `M h

ż Ω ˆˆ1 ´9 4 ε 2 ˙px ´xmin q ´M h 2{3 ˙|e Φ{h ψ| 2 dxdy ď 0 .
We take ε such that 1 ´9 4 ε 2 ą 0 and fix R ą 0 such that

ˆ1 ´9 4 ε 2 ˙R ´M " 1 . (2.5)
We write

ż |x´x min |ěRh 2{3 ˆˆ1 ´9 4 ε 2 ˙px ´xmin q ´M h 2{3 ˙|e Φ{h ψ| 2 dxdy ď ´ż|x´x min |ăRh 2{3 ˆˆ1 ´9 4 ε 2 ˙px ´xmin q ´M h 2{3 ˙|e Φ{h ψ| 2 dxdy ,
and get (using (2.5)):

h 2{3 ż |x´x min |ěRh 2{3 |e Φ{h ψ| 2 dxdy ď ż |x´x min |ěRh 2{3 ˆˆ1 ´9 4 ε 2 ˙px ´xmin q ´M h 2{3 ˙|e Φ{h ψ| 2 dxdy ď M h 2{3 ż |x´x min |ăRh 2{3 e 2Φ{h |ψ| 2 dxdy ď Ch 2{3 }ψ} 2 .
This proves (2.1). Next we prove (2.2). First observe that repeating the calculation above without dropping the first term in (2.4) we get that ż

Ω h 2 |∇e Φ{h ψ| 2 dxdy ď Ch 2{3 }ψ} 2 .
Next, fix 0 ă ε ă ε. We have that

e ε|x´x min | 3 2 {h ph∇ψq " ´3ε 2
?
x ´xmin e ε|x´x min | we complete the proof of (2.2).

We denote by ABpA 0 , ηq the complement of the open disc BpA 0 , ηq.

Corollary 2.2. Let M, η ą 0. There exist c η , C η , h 0 ą 0 such that if h P p0, h 0 q, then any eigenfunction ψ corresponding to an eigenvalue λ ď x min `M h Proof. The set Ω X ABpA 0 , ηq is compact. The map Ω X ABpA 0 , ηq Q px, yq Þ Ñ x ´xmin is continuous and positive, thus it has a positive lower bound. The conclusion follows from Theorem 2.1.

Tubular coordinates and localized operator

We can now reduce our investigation to a neighborhood of A 0 . We use tubular coordinates in a neighborhood of A 0 (see for instance [START_REF] Fournais | Spectral methods in surface superconductivity[END_REF]). Due to Assumption 1.1 there exist δ 1 , δ 2 ą 0 such that (see Figure 2)

W δ " p´δ 1 , δ 1 q ˆp0, δ 2 q Q ps, tq Þ Ñ τ ps, tq " γpsq ´tnpsq P τ pW δ q induces a diffeomorphism. Here, n is the outward pointing normal and γ is the natural length-parametrization of the boundary BΩ; we set γp0q " A 0 . Denoting by θpsq the turning angle at the point γpsq P BΩ we may write npsq " pcos θpsq, sin θpsqq, the tangent vector γ 1 psq " p´sin θpsq, cos θpsqq, and the curvature κpsq " θ 1 psq. The Jacobian of τ " pτ 1 , τ 2 q " τ 1 e 1 `τ2 e 2 is given by mps, tq " 1 ´κpsqt .

(3.1)

We fix δ 2 ą 0 so small that the map τ induces a local diffeomorphism between a rectangle and the tubular neighborhood of the boundary. In view of the definition of A 0 we have

n 1 p0q ¨e1 " 0 , γ 1 p0q ¨e1 " 0 , γ 2 p0q " ´κ0 np0q " κ 0 e 1 , (3.2)
where κ 0 " κp0q.

3.2. Spectral reduction to a localized operator. For 0 ă δ ď δ 1 , δ 2 we define L h,δ to be the Dirichlet realization of ´h2 B 2

x ´h2 B 2 y `x on L 2 pτ pW δ qq with W δ " p´δ, δq ˆp0, δq. We denote by λ n ph, δq the corresponding eigenvalues of L h,δ . By using the decay estimates of Theorem 2.1, we get the following: Proposition 3.1. For every fixed n ě 1 there exist k, K ą 0 such that λ n ph, δq ´Ke ´k{h ď λ n phq ď λ n ph, δq .

Proof. The second inequality is a direct consequence of the min-max variational principle since the form domain of L h,δ is included in the one of L h .

We now prove the first inequality but only for n " 1. Let ψ 1 be a unit eigenvector of L h corresponding to its groundstate λ 1 phq. Let 0 ď χ δ ď 1 be a smooth cut-off function such that χ δ pτ ps, tqq " 1 if |s|, t ď δ{2, and χ δ pτ ps, tqq " 0 if |s| ě 3 4 δ or t ě 3 4 δ. Then χ δ ψ 1 belongs to the domain of L h,δ and since the support of the derivatives of χ δ is away from A 0 , from Corollary 2.2 we have }χ δ ψ 1 } " 1 `Ope ´k{h q, }p∇χ δ q∇ψ 1 } " Ope ´k{h q, }pL h,δ ´λ1 phqqχ δ ψ 1 } " Ope ´k{h q.

Then the min-max principle implies:

λ 1 ph, δq ď xχ δ ψ 1 , L h,δ χ δ ψ 1 y }χ δ ψ 1 } 2 ď λ 1 phq `Ope ´k{h q.
If n ą 1, one can construct n quasi-modes for L h,δ out of the first n eigenmodes of L h by using a similar argument.

Therefore, we can focus on the spectral analysis of L h,δ . The operator L h,δ is unitarily equivalent to the Dirichlet realization of Note that for δ small enough, there exists α ą 0 such that for all ps, tq P p´δ, δqˆp0, δq: τ 1 ps, tq ě x min `αps 2 `tq .

M h " ´h2 m ´1B s m ´1B s ´h2 m ´1B t mB t `
Proposition 3.2. Let M ą 0. There exist ε, C, h 0 ą 0 such that, for all h P p0, h 0 q, and for all eigenfunctions ψ corresponding to eigenvalues λ of M h with λ ď x min `M h Proof. If δ is small enough then mps, tq " 1 and x ´xmin " t. Therefore we can directly use the strategy of Proposition 2.1 applied for the eigenfunctions of L h,δ , but in both (3.4) and (3.5) we need to choose an ε which is smaller than the one in Proposition 2.1 in order to control the linear growth in t of mps, tq ´1 and B s mps, tq.

Therefore, the operator M h can be replaced by . By using this exponential decay, there exist C, h 0 ą 0 such that, for all h P p0, h 0 q,

Ă M h " ´h2 m ´1B s m ´1B s ´h2 m ´1B t mB t `
µ n phq ´1 C e ´Ch ´3η{2 ď λ n ph, δq ď µ n phq . (3.7)
Thus, modulo an exponentially small error, the asymptotic analysis of λ n phq is reduced to that of µ n phq.

By shrinking the spectral window, we can even get a localization with respect to the s variable, as stated in the next proposition. Proposition 3.3. Let M ą 0 and η P p0, 1 3 q. There exist ε, C, h 0 ą 0 such that, for all h P p0, h 0 q, and for all eigenfunctions ψ of Ă M h corresponding to eigenvalues λ ď x min `h 2 3 z 1 `M h we have:

ż W δ,h e εs 2 {h |ψ| 2 dsdt ď C}ψ} 2 .
Proof. The proof follows the same lines as that of Proposition 2.1. We let Φpsq " εs 2 {2 and write the Agmon formula:

ż W δ,h
´m´2 |hB s e Φ{h ψ| 2 `|hB t e Φ{h ψ| 2 `τ1 ps, tq|e Φ{h ψ| 2 ´pλ `|∇Φ| 2 q|e Φ{h ψ| 2 ¯mdsdt " 0 .

First, we drop the tangential derivative:

ż W δ,h
´|hB t e Φ{h ψ| 2 `τ1 ps, tq|e Φ{h ψ| 2 ´pλ `|∇Φ| 2 q|e Φ{h ψ| 2 ¯mdsdt ď 0 .

We observe from (3.3) that there exist k ą 0 such that τ 1 ps, tq ě x min `t `ks 2 . Introducing the last inequality in the above integral we get:

ż W δ,h
´|hB t e Φ{h ψ| 2 `t|e Φ{h ψ| 2 ¯mdsdt `żW δ,h p´λ `xmin `ks 2 ´|∇Φ| 2 q|e Φ{h ψ| 2 mdsdt ď 0 .

On W δ,h , for sufficiently small h , there exists C ą 0 such that mps, tq ě 1 ´Ch 2{3´η . By using this in the above integrals we have that

p1 ´Ch 2{3´η q ż W δ,h
´|hB t e Φ ψ| 2 `t|e Φ{h ψ| 2 ¯dsdt `żW δ,h p´λ `xmin `ks 2 ´|∇Φ| 2 q|e Φ{h ψ| 2 mdsdt ď 0 .

Then, with a possibly larger constant C we have that

p1 ´Ch 2{3´η q ż W δ,h ´|hB t e Φ ψ| 2 `t|e Φ{h ψ| 2 ´h 2 3 z 1 |e Φ{h ψ| 2 ¯dsdt `żW δ,h p´λ `xmin `h 2 3 z 1 `ks 2 ´|∇Φ| 2 ´Ch 4 3 ´ηq|e
Φ{h ψ| 2 mdsdt ď 0 .

By using the min-max principle and the Dirichlet bracketing (only with respect to t, s being fixed), we have

ż W δ,h ´|hB t e Φ ψ| 2 `t|e Φ{h ψ| 2 ´h 2 3 z 1 |e Φ{h ψ| 2 ¯dsdt ě 0 . Therefore, ż W δ,h ´p´λ `xmin `h 2 3 z 1 `ks 2 ´|∇Φ| 2 ´Ch 4 3 ´ηq|e
Φ{h ψ| 2 ¯mdsdt ď 0 , so that, using the assumption on the location of λ, we obtain

ż W δ,h p´M h `ks 2 ´|∇Φ| 2 ´Ch 4 3 ´ηq|e
Φ{h ψ| 2 ¯mdsdt ď 0 . Now if η P `0, 1 3 ˘is kept fixed and h is small enough, then: ż

W δ,h
p ks 2 ´|∇Φ| 2 ´2M hq|e Φ{h ψ| 2 mdsdt ď 0 , and ż

W δ,h ´pk ´ε2 qs 2 ´2M h ¯|e Φ{h ψ| 2 dsdt ď 0 .
For ε small enough, the conclusion follows as in the proof of Proposition 2.1.

This shows that the eigenfunctions of "low energy" are localized near A 0 at a scale h 1{2 in the s direction, and at a scale h 2{3 in the t direction.

Proof of the main theorem

In view of Equation (3.7) and Proposition 3.1 our main result Theorem 1.2 is a direct consequence of the next proposition that provides the asymptotic behavior of the lowlying eigenvalues, µ n phq, of the operator Ă M h defined in the previous section. In the remainder of this section we show the above proposition by obtaining suitable upper and lower bounds. M h and apply the min-max principle. Consider a smooth function χ with compact support equal to 1 near p0, 0q on a scale of order h 1{2´θ in the s direction, and on a scale of order h 2{3´θ in the t direction, for some θ P p0, ηq. Let us introduce the following family of test functions ϕ n,h ps, tq " χps, tqf n,h psqa h ptq , where i. f n,h psq " h ´1 4 f n ph ´1 2 sq and pf n q nPNzt0u is an L 2 pRq-normalized family of eigenfunctions of the harmonic oscillator ´h2 B 2 s `κ0 2 s 2 (rescaled Hermite functions), ii. a h ptq " h ´1 3 Aiph ´2 3 t ´z1 q where Ai is the L 2 pR `q normalized Airy function and ´z1 ă 0 is its first zero. In particular, we have that a h ptq satisfies p´h 2 B 2 t `tqa h ptq " z 1 h 2{3 a h ptq for t ą 0 with Dirichlet boundary conditions. An explicit computation shows that for any n P N there exist constants ε, c ą 0 such that ż Using this we immediately see that we find c 1 , c 2 ą 0 such that

1 ě }ϕ n,h } 2 ě 1 ´c1 ż t |ϕ n,h ps, tq| 2 dsdt ě 1 ´c2 h 2{3 . (4.2)
We are interested in estimating the matrix elements xϕ n 1 ,h , p Ă M h ´xmin qϕ n,h y. First, we notice that the support of χ ´1 and the supports of the derivatives of χ are located in a region where either h ´1{2 |s| ě Ch ´θ or h ´2{3 t ě Ch ´θ, thus all the integrals entering the matrix element which contain such derivatives will be of order h ´8 due to the exponential localization of the f n 's and of Ai. Second, the operators hB s and hB t acting on f n,h 's and a h respectively will generate a factor of at most h 1{3 , and each integral contains two such factors; thus each term in the scalar product has an order of magnitude of at most h 2{3 . Third, if we replace the function m in (3.1) by 1, the error contains an extra factor t which due to the decay of a h may be replaced by h 2{3 . Together with the a-priori decay of h 2{3 coming from the derivatives, this error term will grow at most like h 4{3 .

Moreover, replacing τ 1 ps, tq ´xmin with κ 0 2 s 2 `t will produce an error like h 3{2 . Thus we may write:

xϕ n 1 ,h , p Ă M h ´xmin qϕ n,h y "h 2 ż R ´Bs f n,h B s f n 1 ,h `κ0 2 s 2 f n,h f n 1 ,h ¯ds `δn,n 1 ż R ``|B t a h | 2 `t|a h | 2 ˘dt `Oph 4{3 q "δ n,n 1 ˆp2n ´1qh c κ 0 2 `z1 h 2{3 ˙`Oph 4{3 q. (4.3)
Using similar arguments, the Gram-Schmidt matrix elements xϕ n 1 ,h , ϕ n,h y will equal δ n,n 1 `Oph 2{3 q. Thus if N is fixed and h is small enough, the subspace span 1ďjďN ϕ j,h

will have dimension N and we may find an orthonormal basis tψ n,h u N n"1 such that

ψ n,h " N ÿ j"1 c j ϕ j,h , c j " δ n,j `Oph 2{3 q .
Thus the matrix elements xψ n 1 ,h , p Ă M h ´xmin qψ n,h y will obey the same estimate as in (4.3), which via the min-max principle imply

µ n phq ď x min `z1 h 2 3 `hp2n ´1q c κ 0 2 `Oph 4 3 q, 1 ď n ď N.
4.2. Lower bound. Let N ě 1 and consider a family of eigenfunctions pψ j,h q 1ďjďN associated with the eigenvalues pµ j phqq 1ďjďN . We let

E N phq " span 1ďjďN ψ j,h Ă L 2 pW δ,h , mdsdtq .
Note that the decay estimates of Propositions 3.2 and 3.3 can be extended to ψ P E N phq.

Let us choose any ψ P E N phq with norm one. Because mps, tq ď 1, we have the important inequality

1 " ż W δ,h |ψ| 2 mdsdt ď }ψ} 2 L 2 pW δ,h ;dsdtq . (4.4)
We also have

µ N phq ě x Ă M h ψ, ψy " ż W δ,h `m|hB t ψ| 2 `m´1 |hB s ψ| 2 `mτ 1 ps, tq|ψ| 2 ˘dsdt . (4.5)
We recall that mps, tq " 1 ´κpsqt so that, using (3.5) ) .

Thus, for h small enough, the N -th eigenvalue of A h b 1 `1 b H h , denoted by ν N phq, is given by

ν N phq " z 1 h 2{3 `p2N ´1qh c κ 0 2 . (4.7)
Notice that the set E N phq is contained in the form domain of A h b 1 `1 b H h and, seen as a subset of L 2 pW δ,h ; dsdtq, still has the dimension N for h small enough. This is because the ψ j,h 's are almost orthogonal in the "flat" space, up to an error of order h 2{3 . Thus, by the min-max principle, we have

ν N phq ď sup ψPE N phq A´A h b 1 `1 b H h ¯ψ, ψ E L 2 pdsdtq }ψ} 2 L 2 pdsdtq
ď pµ N phq ´xmin `Ch 4{3 q , (4.8)

where in the last inequality we used (4.6). This implies the desired lower bound and concludes the proof of Proposition 4.1.

Figure 1 .

 1 Figure 1. The domain Ω.

2 3 , 3 2Ω e ε|x´x min | 3 2

 333 and all corresponding eigenfunctions ψ, ż Ω e ε|x´x min | {h |ψ| 2 dx ď C}ψ} 2 , (2.1) and ż {h |h∇ψ| 2 dx ď Ch 2{3 }ψ} 2 .

3 2? t e ´pε´εqt 3 2

 33 {h ψ `h∇ ˆeε|x´x min | 3 2 {h ψ ˙. Using the triangle inequality, (2.3), (2.1) and the estimate sup tě0 {h ď Ch 1 3 ,

2 3

 2 satisfies the estimates ż ΩXABpA 0 ,ηq |ψ| 2 dx ď C η e ´cη{h }ψ} 2 and ż ΩXABpA 0 ,ηq |∇ψ| 2 dx ď C η e ´cη{h }ψ} 2 .

Figure 2 . 1 .

 21 Figure 2. The local tubular coordinates

3 2δ e εt 3 2

 33 {h |ψ| 2 dsdt ď C}ψ} 2 , (3.4) ż W {h |h∇ s,t ψ| 2 dsdt ď Ch 2{3 }ψ} 2 .(3.5)

Proposition 4 . 1 . 2 3

 412 Let n P t1, 2, . . . u. Then, as h Ñ 0 we have µ n phq " x min `z1 h

4. 1 .

 1 Upper bound. To get the upper bound in Theorem 1.2, it is sufficient to use convenient test functions in the domain of Ă

R e εs 2 3 2

 23 {h |f n,h psq| 2 ds ď c , {h |a h ptq| 2 dt ď c . (4.1)

  Φ{h ψ| 2 `px ´λ ´|∇Φ| 2 q|e Φ{h ψ| 2 dxdy " 0 ,

	2 3 . We get
	ż
	(2.3)
	and thus
	ż

Ω h 2 |∇e Ω h 2 |∇e Φ{h ψ| 2 `px ´xmin ´M h 2{3 ´|∇Φ| 2 q|e Φ{h ψ| 2 dxdy ď 0 .

(2.4) Now, we choose Φ " ε|x ´xmin | 3{2 and drop the first term above to get that

  τ1 ps, tq , with Dirichlet boundary conditions, acting on L 2 pW δ,h , mps, tqdsdtq, with W δ,h " p´δ, δq ˆp0, h ´m´2 |hB s ψ| 2 `|hB t ψ| 2 `τ1 ps, tq|ψ| 2 ¯mdsdt . (3.6) Let µ n phq be the asociated (ordered) eigenvalues of Ă M h . The decay estimates of Proposition 3.2 are still satisfied by the eigenfunctions of Ă M h with eigenvalues λ ď x min `M h

	2 3 ´ηq ,
	for some η P `0, 1 3 ˘with corresponding quadratic form
	ż
	xψ, Ă M h ψy "
	W δ,h
	2
	3

  in the second inequality below, `|hB t ψ| 2 `|hB s ψ| 2 `mτ 1 ps, tq|ψ| 2 ˘dsdt ´C ż `|hB t ψ| 2 `|hB s ψ| 2 `mτ 1 ps, tq|ψ| 2 ˘dsdt ´Ch 4{3 , " x min `żW δ,h `|hB t ψ| 2 `|hB s ψ| 2 `mpτ 1 ps, tq ´xmin q|ψ| 2 ˘dsdt ´Ch 4{3 ě x min `żW δ,h ˆ|hB t ψ| 2 `|hB s ψ| 2 `ˆt `κ0 s 2 2 ˙|ψ| 2 ˙dsdt ´Ch 4{3 p|s 2 t| `|t 2 |q|ψ| 2 dsdt,the last integral can be estimated to be of order h 4{3 as well using the exponential decay in the s and t variables (Propositions 3.2 and 3.3). Then, there is a C ą 0 such thatx Ă M h ψ, ψy ě x min `żW δ,h ˆ|hB t ψ| 2 `|hB s ψ| 2 `ˆt `κ0 s 2 2 ˙|ψ| 2 ˙dsdt ´Ch 4{3 . ˆ|hB t ψ| 2 `|hB s ψ| 2 `ˆt `κ0 s 2 2 ˙|ψ| 2 ˙dsdt ď µ N phq ´xmin `Ch 4{3 (4.6)ď pµ N phq ´xmin `Ch 4{3 q}ψ} 2 L 2 pdsdtq . On the left hand side of the above identity we recognize the quadratic form associated to the operator A

	ż
	x Ă M h ψ, ψy ě					t|h∇ s,t ψ| 2 dsdt
	W δ,h				W δ,h
	ż
	ě
	W δ,h
	ż
	´C
	W δ,h
	Now using the inequalities in (4.5) and (4.4) we have:
	ż
	W δ,h
	)	,	sppH h q "	! p2n ´1qh c	κ 0 2	, n ě 1

h b 1 `1 b H h , where A h " ´h2 B 2 t `t on L 2 pp0, `8q, dtq with Dirichlet boundary conditions and H h " ´h2 B 2 s `κ0 s 2 2 on L 2 pp´8, 8q, dsq. The spectrum of A h and H h is given by sppA h q " ! z n h 2{3 , n ě 1

Acknowledgment. D.K. has been supported by the EXPRO grant No. 20-17749X of the Czech Science Foundation (GACR). D.K. is also grateful to the Aalborg University for supporting his stay in 2019. E.S. has been partially funded by Fondecyt (Chile) project # 118-0355 and thanks Esteban Ramos Moore for stimulating discussions. H.C., E.S, and N.R. are deeply grateful to the Mittag-Leffler Institute where this collaboration was stimulated during the thematic semester "Spectral Methods in Mathematical Physics" in 2019. E.S. and N.R. are also grateful to the CIRM where the ideas of this paper were discussed ("Research in Pairs" session in 2019).