Laurent H Coudert

Bérenger Gans

Gustavo Garcia

Jean-Christophe Loison

Renner-Teller effects in the photoelectron spectra of CNC, CCN, and HCCN

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Photoelectron spectroscopy allows us to study complicated rovibronic transitions, between neutral and cationic species, requiring in many cases dedicated models. This is well illustrated by the photoelectron spectra of ethylene 1 and methane 2 which required new models to compute the torsional energy levels of C 2 H + 4 and to treat the Jahn-Teller effect in CH + 4 . This is also well illustrated by the photoelectron spectrum of acetylene 3,4 which could be analyzed once the unusual tetra-atomic Renner-Teller coupling 5 in C 2 H + 2 was accounted for.

Recently, the photoelectron spectra of several radicals, including CNC, CCN, and HCCN, have been reported. 6 Just as the acetylene molecule, these species or their cation display a Renner-Teller (RT) coupling, but modeling their photoionization spectra is expected to be theoretically more challenging for several reasons. A larger RT interaction characterized by a Renner parameter ϵ of the order of 0.5 in CNC and CCN, and larger than 1 in HCCN + should be accounted for. In the quasilinear CNC + , CCN + , and HCCN, the bending mode should be dealt with as a large amplitude motion since it is characterized by a very low frequency on the order of 100 cm -1 . More importantly, although theoretical calculations of the line strength of photoelectron spectra are available for diatomic [7][8][9][10] and polyatomic 1,[11][12][13][14][15] molecules, there are no results for molecules displaying a RT coupling.

In this paper, a theoretical formalism aimed at accounting for the photoionization spectrum of non-rigid molecules displaying a RT coupling is developed and applied to the simulation of the photoelectron spectrum 6 of CNC, CCN, and HCCN. Starting from ab initio calculations, bending potentials were derived for these three neutral species and their cation.

The rovibrational energy levels of CNC + , CCN + , and HCCN, characterized by a nondegenerate sigma electronic state, were computed using an approach in which the overall rotation and the large amplitude bending mode are treated simultaneously [16][17][18][19][20][21] in order to account for the quasilinearity. Using previous theoretical results, [22][23][24][25][26][27] this approach is extended to include the RT effect and allows us to retrieve the rovibronic energy levels of CNC, CCN, and HCCN + . Finally, an expression for the line strength of photoelectron spectra when one of the electronic states displays the RT coupling is derived and written in terms of Franck-Condon factors involving the bending wavefunctions of the neutral and cationic species.

Since none of the photoionization spectra dealt with are rotationnaly resolved, 6 only the K structure is accounted for in the energy level calculation and the line strength is averaged over rotational levels.

These theoretical results are first applied to the ionizing transition of H 2 O which provides us with a test for the theoretical approach since the cation H 2 O + is strongly affected by the RT coupling. The threshold photoelectron spectrum (TPES), calculated from available potential energy surfaces 28,29 and using the newly derived line strength expression, compares fairly well with the experimental one. 30,31 The theoretical results are then applied to the modeling of the photoelectron spectra of CNC, CCN, and HCCN. These spectra turned out to be in good agreement with the experimental ones, 6 even for HCCN where the RT coupling in its cationic species is the largest.

This paper has five remaining sections. In Section II, the results of the ab initio calculations and their fitting are presented. Section III is concerned with the effective approaches used to account for the quasilinearity and the RT coupling. In Section IV, the expression of the photoionization cross section when one of the electronic states displays the RT coupling is introduced. These results are applied to H 2 O, CNC, CCN, and HCCN in Section V.

Section VI is the discussion.

II. AB INITIO CALCULATIONS

The ab initio calculations of the electronic states of the CNC, CCN, and HCCN rad- Numerical values were least squares fitted to a polynomial expansion F (γ) written in term of π -γ, the supplement of γ:

F (γ) = n ∑ i=0 f i (π -γ) i (1)
where f i , with 0 ≤ i ≤ n, are the fitted constants. 4 such constants were retrieved for each electronic state and for each parameter allowing us to reproduce accurately the ab initio values. For potential energies and bond lengths, the expansion in Eq. ( 1) was restricted to even values of i. For a pair of RT electronic substates, the f 0 parameters of both substates were constrained to be equal to ensure the required degeneracy for the linear configuration.

Due to a convergence problem, the ab initio calculations could not be carried out for γ < γ min , where the angle γ min is between 90 and 115 • depending on the species and the electronic state. In order to obtain physically meaningful results, Eq. ( 1) should not be used for values of γ smaller than γ min .

For the CNC radical, fitting of the ab initio potential energy points with the expansion in Eq. ( 1) yielded root-mean-square (RMS) deviations of 2.5 and 2.4 cm -1 for the neutral and cationic species, respectively. It was also possible to retrieve an accurate value for the Renner parameter ϵ. A value of 0.562 was obtained and turned out to be in good agreement with that reported by Merer and Travis, 34 0.549.

Similarly for the CCN radical, the fit of the ab initio potential energy points led to root-mean-square (RMS) deviations of 1.5 and 26 cm -1 for the neutral and cationic species, As can be seen in Table I, the HCCN radical and the HCCN + cation display a nearly linear equilibrium geometry. The ab initio calculations also confirm that the neutral and cationic species are planar 32,39,40 and that when γ = 180 • , the bending angle β = ∠CCN is 180 • . When fitting the ab initio values of this angle, the angle supplement π -β was fitted to the expansion of Eq. ( 1) which was restricted to odd values of i. This ensures the required linearity for γ = 180 • . The angle β is plotted in Fig. 4 as a function of γ. For γ ≥ 110 • , the results are consistent with the molecule assuming a Z-shaped trans geometry for the X 3 A ′′ electronic state of the neutral radical and the X + 2 A ′ electronic substate of the cation as the dihedral angle ∠HCCN is 180 • . For the A + 2 A ′′ electronic substate, the angle β displays almost no change and remains close to 180 • . The RMS deviations of the fits of the ab initio potential energy points are 12 and 4.1 cm -1 for the neutral and the cationic species, respectively.

The results of the fit of the ab initio bond lengths with the polynomial expansion in Eq. ( 1)

is summarized for all species in the form of figures available as supplementary material.

III. EFFECTIVE BENDING AND RT HAMILTONIANS

The systems dealt with in the present work display either a quasilinearity problem like the CNC + and CCN + cations, 36 and the HCCN radical 32,42 or a strong RT effect like the CNC and CCN radicals, 26,34,38 and the HCCN + cation. 39 In this section, effective Hamiltonians are introduced to deal with both effects and account simultaneously for the large amplitude bending mode and the overall rotation.

A. Effective bending Hamiltonian

The effective bending-rotation Hamiltonian H b-r introduced to treat the quasilinearity accounts for the large amplitude nature of the bending mode and relies on the Bending-Rotation approach [16][17][18][19][20][21] developed to deal with the anomalous centrifugal distortion of the water molecule. 43 The Bending-Rotation approach [16][17][18][19][20][21] accounts exactly for the fact that the A rotational constant goes to infinity for the linear configuration. H b-r is written in term of the large amplitude coordinate t = cos γ. Following Hougen et al. 44 and Mekhtiev et al., 45 a reference configuration represented by atom positions a i (t) is used and a four-dimensional effective Hamiltonian describing the large amplitude bending motion and the overall rotation is derived. A molecule fixed xyz axis system is chosen such that its origin is the molecular center of mass, its xz plane is the molecular plane, and the molecule is along the z axis for the t = -1 linear configuration. Obtaining atom positions a i (t) and their derivatives ∂a i (t)/∂t from Section II, the 4×4 symmetrical generalized inertia tensor 44 I(t) is computed using Eqs. [5] of Mekhtiev et al. 45 The effective Hamiltonian H b-r is then expressed with the generalized inverse inertia tensor 44,45 µ(t) = [I(t)] -1 as:

H b-r = 1 2 P t µ tt (t)P t + 1 2 ∑ δ=x,y,z µ δδ (t)N 2 δ + 1 2 µ xz (t){N x , N z } + 1 2 {µ yt (t), P t }N y + V (t), (2)
where µ tt (t), µ xx (t), µ yy (t), µ zz (t), µ xz (t), and µ yt (t) are components of the generalized inverse inertia tensor; 44,45 {, } is the anticommutator; P t is the momentum conjugate to t; N δ , with δ = x, y, z, are molecule fixed components of the rotational angular momentum N;

and V (t) is the potential energy function. In Eq. ( 2), only two non-diagonal components of the generalized inverse inertia tensor arise due to the symmetry plane. A term involving the determinant of the µ(t) tensor giving rise to a mass dependent potential analogous to those in Eq. ( 36) of Hougen et al. 44 and in Eq. [10] of Mekhtiev et al. 45 has been omitted in Eq. ( 2). The volume element to be used for the effective Hamiltonian H b-r is sin θdθdϕdχdt.

Although there are no analytical expression for the components of the generalized inverse inertia tensor, in the linear limit, when t → -1, the two non-diagonal components µ xz (t)

and µ ty (t) go to zero and the four diagonal components display the following behavior:

lim t→-1      µ zz (t) = A 1 + t , µ xx (t) = 2B, µ tt (t) = 2B e (1 -t 2 ), µ yy (t) = 2B, (3)
where A, B e , and B are three kinetic energy parameters. Although the meaning of A and B e is not obvious, B clearly is the rotational constant of the linear configuration. The results in Eqs. ( 3) are well illustrated by the CNC + cation. Taking into account the above requirements, the molecule fixed axis system is attached to the molecule so that its x axis bisects the ∠CNC bending angle. Figure 5 depicts the variations of

µ tt (t)/(1 -t 2 ), µ xx (t),
µ yy (t), and (1 + t)µ zz (t). All the curves in this figure display a smooth behavior and it can be deduced that A = B e = 4.9 cm -1 and B = 0.451 cm -1 .

Equations (3) mean that when t → -1, the limiting behavior of the terms in the effective Hamiltonian of Eq. ( 2) is the same as those in the Bending-Rotation Hamiltonian defined in Eqs. ( 1)-(3) of Coudert. 16 This suggests that the Schrödinger equation for the bending mode should be solved using the same θ αβ n (t) basis set functions:

θ αβ n (t) = (1 -t) α/2 (1 + t) β/2 P (α,β) n (t)/ √ h n , ( 4
)
where n is an integer with 0

≤ n ≤ n Max ; P (α,β) n (t), with α, β ≥ -1, is a Jacobi polynomial; 46
and h n is a normalizing factor given in Table 22.2 of Abramovitz and Stegun's book 46 and in Eq. ( 8) of Coudert. 16 The Schrödinger equation for the effective Hamiltonian in Eq. ( 2) is solved [16][17][18][19][20][21] using rovibrational wavefunctions of the form:

ψ N,k (t)|N, k⟩, ( 5
)
where ψ N,k (t) is a vibrational bending wavefunction and |N, k⟩ is a symmetric top rotational function defined as in Wigner [START_REF] Wigner | Group Theory[END_REF] and characterized by N the quantum number of the total angular momentum and by k the eigenvalue of its molecule fixed component N z . The symmetric top rotational function depends also on M , the eigenvalue of the laboratory fixed component N Z , but this quantum number is omitted as the energy does not depend on M .

The vibrational bending wavefunction ψ N,k (t) is a solution of the Schrödinger equation for the bending Hamiltonian:

H N,k b = 1 2 P t µ tt (t)P t + 1 2 k 2 µ zz (t) + V (t) + 1 4 [N (N + 1) -k 2 ][µ xx (t) + µ yy (t)]. (6)
The matrix of the bending Hamiltonian in Eq. ( 6) is set up and diagonalized taking the basis set functions in Eq. ( 4) with β such that: 16-21

β 2 = k 2 lim t→-1 4µ zz (t)(1 + t) 2 µ tt (t) . ( 7
)
This choice accounts for the singularity of the bending Hamiltonian H N,k b for the linear configuration and ensures that the basis set functions in Eq. ( 4) are eigenfunctions of H N,k b -V (t) when t is close to -1. In the case of the C 2v symmetry CNC + cation, Eq. ( 7) leads to β = |k|. The value of α can be estimated matching the n = β = 0 basis set function of R-T effect in photoelectron spectra Eq. ( 4) with the ground state eigenfunction of the harmonic oscillator bending Hamiltonian. This leads to:

α = 4 √ f 2 /B e , ( 8
)
where f 2 is the i = 1 expansion coefficient of the potential energy function with Eq. ( 1) and B e is defined in Eqs. (3). Using the results of Section II in the case of the CCN + cation, we obtain B e = 5.01 and f 2 = 475 cm -1 leading to α = 39. With such a value, the basis set functions of Eq. ( 4) with n ≤ 10 are vanishingly small outside the range -1 ≤ t ≤ 0.5, that is, 60

• ≤ γ ≤ 180 • .
The matrix elements of the operators in Eq. ( 6) between two basis set functions of Eq. ( 4)

were computed using Gaussian quadrature and accounting for the singularity of µ zz (t) and

P t µ zz (t)P t at t = -1.
For the former operator, the singularity is due to its behavior in the ). This means that values of γ, with γ < γ min , for which the results in Section II are not valid, will not be used. This provides us with an alternative way of choosing α.

Bending-rotation energies and eigenfunctions of the Hamiltonian in Eq. ( 2) will be written:

E(v, N, k) and Ψ v,N,k = ψ N,k v (t)|N, k⟩, ( 9
)
where v is the bending quantum number and ψ N,k v (t) is the corresponding eigenfunction. Symmetry adapted bending-rotation wavefunctions for the C s symmetry group can be built starting from Eq. ( 9), using the invariance of the bending Hamiltonian in Eq. ( 6

|N Kδ⟩ = (|N, K⟩ + δ|N, -K⟩)/ √ 2(1 + δ K,0 ), ( 10
)
where K ≥ 0 and δ = ±1, the symmetry adapted wavefunctions take the form:

Ψ v,N,K,δ = ψ N,K v (t)|N Kδ⟩. ( 11
)
These wavefunctions belong to the symmetry species A ′ and A ′′ of C s when δ(-1) N +K is +1 or -1, respectively. The symmetry species of the rovibronic wavefunction should be calculated replacing this term by δ(-1) N +K g(X Σ), where the g(X Σ) depends on the C s symmetry species of the nondegenerate X Σ electronic state and is +1 (-1) for an A ′ (A ′′ ) electronic state.

B. Effective RT Hamiltonian

The RT effect in systems with a nonzero projection Λ = 1 of the molecule fixed component L z of the electronic angular momentum L are considered. The two RT electronic substates that become degenerate in the linear limit belong to the A ′ and A ′′ symmetry species of C s and their potential energy functions, denoted V ′ (t) and V ′′ (t), respectively, can be obtained from Section II. In the present treatment, we only consider the large amplitude bending mode, the electron spin is ignored, and we use most of the results in Section III A. The molecule fixed axis system is attached to the molecule in the same way and a unique reference configuration is chosen for both electronic substates. The corresponding atom positions a i (t) are determined from the results of Section II and calculated from internal coordinates taken as the average of the A ′ and A ′′ substates values. This choice leads to a single generalized inertia tensor calculated as in the previous section. In agreement with previous investigations, [22][23][24][25][26][27] we use electronic wavefunctions of the form |±Λ⟩ = (|A ′ ⟩ ± i|A ′′ ⟩)/ √ 2 which, in the linear limit, are eigenfunctions of L z with eigenvalues ±Λ. The Schrödinger equation for the effective RT Hamiltonian used in this investigation is solved using rovibronic wavefunctions of the form:

ψ +Λ,N,k (t)|+Λ⟩|N, k⟩ + ψ -Λ,N,k (t)|-Λ⟩|N, k⟩, ( 12
)
where |N, k⟩ are symmetric top rotational functions defined as for Eq. ( 5) and ψ ±Λ,N,k (t) are two vibrational functions. For this pair of functions the Schrödinger equation is expressed as a 2 × 2 matrix:

  H +Λ,N,k b (t) V -(t) V -(t) H -Λ,N,k b (t)   , (13)
H ±Λ,N,k b = 1 2 P t µ tt (t)P t + 1 2 (k ± Λ) 2 µ zz (t) + 1 4 [N (N + 1) -k 2 ][µ xx (t) + µ yy (t)] + V + (t), (14)
where µ δδ (t), with δ = t, x, y, z, are diagonal components of the generalized inertia tensor and

V + (t) = [V ′ (t) + V ′′ (t)]/2.
The matrix of the Hamiltonian in Eq. ( 13) is set up and diagonalized expanding the vibrational functions ψ ±Λ,N,k (t) in terms of the basis set functions of Eq. ( 4). The corresponding values of α and β are denoted α ± and β ± . As Eqs. (3) are also valid, an equation similar to Eq. ( 7) holds:

β ± 2 = (k ± Λ) 2 lim t→-1 4µ zz (t)(1 + t) 2 µ tt (t) . ( 15
)
The value of α ± can be estimated from the potential energy functions of the A ′ and A ′′ substates. Using Eq. ( 14), an equation similar to Eq. ( 8) arises:

α ± = 4 √ f + 2 /B e , ( 16
)
where f + 2 is the average value of the f 2 expansion coefficients of Eq. ( 1) for the potential energy function of the A ′ and A ′′ substates; and where B e is defined in Eqs. (3). Evaluation of the matrix elements needed to set up the matrix of the Hamiltonian in Eq. ( 13) should be carried out as in Section III A. For H ±Λ,N,k b (t), Eqs. (A4)-(A6) should be used; for the potential energy function term V -(t), the Gauss-Jacobi quadrature for the weight function

(1 -t) (α + +α -)/2 (1 + t) (β + +β -)/2 should be utilized.
Symmetry adapted rovibronic functions for the C s point group, denoted Ψ v,Λ,N,K,δ , where v is the bending quantum number, K ≥ 0, and δ = ±1, are expressed in terms of the electronic wavefunctions of the A ′ and A ′′ substates; and of the Wang-type rotational functions |N Kδ⟩ of Eq. ( 10) as:

Ψ v,Λ,N,K,δ = 1 √ 2 {|A ′ ⟩[ψ +Λ,N,K v (t) + ψ -Λ,N,K v (t)]|N Kδ⟩ + i|A ′′ ⟩[ψ +Λ,N,K v (t) -ψ -Λ,N,K v (t)]|N K-δ⟩}, ( 17
)
where ψ +Λ,N,K v (t) and ψ -Λ,N,K v (t) are solutions of Eq. ( 13) characterized by the bending quantum number v and fulfilling:

∫ +1 -1 dt {|ψ +Λ,N,K v (t)| 2 + |ψ -Λ,N,K v (t)| 2 } = 1. ( 18
)
The wavefunctions in Eq. ( 17) belong to the symmetry species A ′ and A ′′ of C s when δ(-1) N +K is +1 or -1, respectively. This equation is physically more meaningful than Eq. ( 12) because it is expressed in terms of the actual |A ′ ⟩ and |A ′′ ⟩ electronic wavefunctions of the RT substates.

IV. RT COUPLING AND PHOTOIONIZATION CROSS SECTION

Theoretical calculations of the photoionization cross section were performed for diatomic [7][8][9][10] and polyatomic 1,11-15 molecules. Simplified results are available for molecules displaying the RT coupling. 22,[START_REF] Jungen | [END_REF] In this work the results derived by Willitsch et al. 15 for an asymmetric-top polyatomic molecule are used. The photoionization cross section derived by these authors is given in the first of their Eqs. (7). This equation, obtained ignoring the electron spin, involves electronic and rovibrational terms. The latter, given in the second of their Eqs. (7), consists of a Franck-Condon factor (q v ) 2 = |⟨v + |v ′′ ⟩| 2 and a rotational factor Q(l ′′ ). These factors are evaluated in this section when the neutral and cationic species are described by the large amplitude approaches introduced in Section III.

A. RT coupling in the X Π state of the neutral species

The electronic state of the cation, denoted X + Σ, is then either an A ′ or A ′′ state. The results of Willitsch et al. 15 are used to calculate the total photoionization cross section σ tot of the X + Σ ← X Π ionizing transition by adding the contribution of each RT component.

Below, these contributions are distinguished using a ± sign where the upper (lower) sign is for the X + Σ ← X A ′ (X + Σ ← X A ′′ ) transition. For either component, the lower level labeled v ′′ , N ′′ , K ′′ , δ ′′ is described by the wavefunction in Eq. ( 17) and connects to the upper level of the cation labeled v + , N + , K + , δ + described by the wavefunction in Eq. ( 11). It can be shown that the following Franck-Condon factor arises:

(q ± v ) 2 = ∫ +1 -1 ψ N + ,K + v + (t)[ψ +Λ,N ′′ ,K ′′ v ′′ (t) ± ψ -Λ,N ′′ ,K ′′ v ′′ (t)]dt 2 , ( 19
)
where the vibrational wavefunctions

ψ N + ,K + v + (t) and ψ ±Λ,N ′′ ,K ′′ v ′′
(t) should be taken from Eqs. ( 9) and ( 12), respectively. The Q(l ′′ ) factor can be written:

Q ± (l ′′ ) = (2N + + 1)     N + l ′′ N ′′ -K + λ ′′ K ′′   ± δ ′′   N + l ′′ N ′′ -K + λ ′′ -K ′′   + δ +   N + l ′′ N ′′ K + λ ′′ K ′′   ± δ + δ ′′   N + l ′′ N ′′ K + λ ′′ -K ′′     2 /h(K + , K ′′ ). ( 20
)
where

h(K + , K ′′ ) = 4(1 + δ K + ,0 )(1 + δ K ′′ ,0
). Taking into account the electronic terms and the contributions from both RT components, the total cross section can be written: 15

σ tot = ρ ′′ [(q + v ) 2 ∑ |λ ′′ |≤l ′′ Q + (l ′′ )B (A ′ ) l ′′ λ ′′ 2l ′′ + 1 + (q - v ) 2 ∑ |λ ′′ |≤l ′′ Q -(l ′′ )B (A ′′ ) l ′′ λ ′′ 2l ′′ + 1 ]/2, ( 21
)
where the electronic terms B

(A ′ ) l ′′ λ ′′ and B (A ′′ )
l ′′ λ ′′ are defined in Eq. ( 10) of Willitsch et al. 15 and correspond to either RT electronic substates; and ρ ′′ is a weighing factor depending on the lower level. The Franck-Condon factors in Eq. ( 21) should be calculated with a Gauss-Jacobi quadrature. Symmetry considerations based on the nature of the partial wave describing the photoelectron 13,15 put some restrictions on the quantum numbers of the upper and lower levels. Equation (3) of Willitsch et al. 15 and the results at the end of Sections III A and III B lead to:

δ + (-1) N + +K + g(X + Σ) δ ′′ (-1) N ′′ +K ′′ = -(-1) l , ( 22
)
where l is the orbital angular momentum of the photoelectron partial wave.

Since the rotational structure is not expected to be resolved in the photoionization spectrum of the heavy molecules studied in this work, the cross section should be averaged over rotational levels. It will be averaged over N + and N ′′ , but not over v + , K + and v ′′ , K ′′ as resolved vibrational and K structure might be observed. The averaged cross section will be calculated as:

σ vib tot (v + , K + ← v ′′ , K ′′ ) = ∑ N ′′ ,δ ′′ ∑ N + ,δ + σ tot . ( 23
)
Using the same ideas as in Section 5 of Buckingham et al. 7 and neglecting in Eq. ( 19) the dependence of the Franck-Condon factors on N + and N ′′ , evaluation of the summation over the rotational quantum numbers of the upper level N + and δ + leads to:

∑ N + ,δ + Q ± (l ′′ ) = (δ K + ,|K ′′ +λ ′′ | + δ K + ,|K ′′ -λ ′′ | )/4, ( 24
)
where N + and δ + are such that Eq. ( 22) is fulfilled and λ ′′ is nonzero. As the result in Eq. ( 24) does not depend on the rotational quantum numbers of the lower level, N ′′ and δ ′′ , the summation over these quantum numbers in Eq. ( 23) reduces to:

∑ N ′′ ,δ ′′ ρ ′′ = exp[-E rv (v ′′ , K ′′ )/kT ]Z r /Z = ρ ′′ (v ′′ , K ′′ )Z r /Z, ( 25
)
where Z is the partition function of the molecule; T the temperature; k the Boltzmann constant; E rv (v ′′ , K ′′ ) the rovibrational energy; and Z r the rotational partition function.

This partition function is assumed to be independent of v ′′ and K ′′ and to depend only on the temperature. The averaged cross section can now be written:

σ vib tot (v + , K + ← v ′′ , K ′′ ) = ρ ′′ (v ′′ , K ′′ )Z r /Z × ∑ |λ ′′ |≤l ′′ [(q + v ) 2 B (A ′ ) l ′′ λ ′′ + (q - v ) 2 B (A ′′ ) l ′′ λ ′′ ]/(2l ′′ + 1) × (δ K + ,|K ′′ +λ ′′ | + δ K + ,|K ′′ -λ ′′ | )/8. (26)
When using this equation, rovibrational energies for given v and K values should be computed setting N equal to K in Eqs. ( 6) and ( 14). The vibrational functions thus obtained should be used to obtain the Franck-Condon factors.

B. RT coupling in the X + Π state of the cationic species

The electronic state of the neutral, denoted X Σ, is then either an A ′ or A ′′ state. The total photoionization cross section σ tot of the X + Π ← X Σ ionizing transition is obtained a Levels are identified by their electronic state, the bent molecule vibrational quantum number v 2 , and the rotational quantum number K as defined in Sections III A and III B. b X 1 A 1 is the ground electronic state of H 2 O, X + 2 B 1 and A + 2 A 1 are respectively lower and upper RT substates of H 2 O + . c Energy in cm -1 calculated by Bunker and Stone 50 for H 2 O and Wu et al. 29 for H 2 O + . For K = 1, the average value of the K-type doublet was taken. d Energy in cm -1 calculated in Section V A taking the same energy origin as Bunker and Stone, 50 and Wu et al. 29 reached setting n Max , the maximum value of n for the basis set functions in Eq. ( 4), to 30, P the number of points of the Gauss-Jacobi quadrature to 33, and α, α ± to 10. Table II displays a comparison between the present results and those previously calculated. 29,50 For H 2 O, the discrepancies range from 0 to 10 cm -1 . For H 2 O + the discrepancies are larger, especially for the A + 2 A 1 state, but the theoretical approach reproduces well the large energy increase when K goes from 0 to 1.

The photoionization cross sections were calculated with Eq. ( 26) adding the contributions of lines with v + ≤ 25, v ′′ ≤ 10, and K + , K ′′ ≤ 12. The adiabatic ionization energy of the X + 2 B 1 state of H 2 O + was set to 12.621 eV, the value determined by Truong et al. 30 from the experimental TPES of water measured with synchrotron radiation between 12 and 40 eV photon energy. Assuming a temperature of 100 K and a Gaussian line shape with a half width at half maximum of 90 cm -1 , the TPES was computed and is plotted in Fig. 6. For emphasis purposes, three cases were treated. In the first and second ones, only X + 2 B 1 ← X 1 A 1 and A + 2 A 1 ← X 1 A 1 transitions were considered setting the (q - v ) 2 and (q + v ) 2 Franck-Condon factors to zero, respectively; in the third case all transitions were considered. The upper panel of Fig. 6 reveals that strong transitions spanning a small energy range arise near the ionization energy due to favorable Franck-Condon factors between the lower X + 2 B 1 RT substate and the ground X 1 A 1 state. The spectrum in this upper panel resembles that shown in Fig. 3 of Truong et al. 30 except that the one calculated in this work does not display the Franck-Condon progression due to the stretching ν 1 mode. As emphasized by the middle panel, unfavorable Franck-Condon factors between the upper A + 2 A 1 RT substate and the ground X 1 A 1 state lead to higher energy transitions. This calculated spectrum should be compared with Fig. 4 of Truong et al. 30 Even though, the ν 1 stretching mode is ignored, we can see that the Franck-Condon progression due to the ν 2 mode spans the same region than that in the observed spectrum. A comparison between the lower panel of Fig. 6 and the TPES in Fig. 2 of Truong et al. 30 reveals an agreement between the line intensities of the low and high energy portions of the spectrum.

The computed TPES could also be compared with the experimental photoelectron spectrum recorded by Ford et al. 31 using a higher resolution than Truong et al. 30 However such a comparison is outside the scope of the present paper as it requires a treatment of the photoionization cross section with values of l ′′ larger than one in Section IV.

B. CNC radical and CNC + ion

The rovibronic energy levels of the CNC radical in its ground X 2 Π g electronic state were computed using the ab initio results in Section II and the procedure described in Section III B. Converged results were obtained setting n Max , P , and α, α ± to 35, 38, and 20, respectively. For low lying rovibronic energy levels, Table III displays a comparison between Observed minus calculated energies are within a few cm -1 except for the highest lying levels which is not involved in any of the vibronic transitions observed by Merer and Travis. 34 The results of Sections II and III A were used to retrieve rovibrational energies for the CNC + cation. No experimental data are available for this species and energies could only be compared with those calculated by Jensen and Kraemer. 36 For the ν 0 2 , 2ν 0 2 , 3ν a Rovibronic levels of the CNC radical are identified with the bent molecule vibrational quantum number v 2 and their symmetry species label Γ. b Observed energy in cm -1 reported by Merer and Travis. 34 c Energy in cm -1 calculated in Section V B taking the same energy origin as Merer and Travis. 34 and are much larger than those calculated in this work: 96, 220, 350, and 494 cm -1 . The large differences stem from the fact that the bending potential retrieved in this work is much shallower and more anharmonic than that used by Jensen and Kraemer. 36 The photoionization cross section were calculated with Eq. ( 26) for lines with v + ≤ 15, v ′′ ≤ 30, and K + , K ′′ ≤ 21. A Gaussian line shape with a half width at half maximum of 90 cm -1 was used when computing the photoelectron spectrum (PES). The energy difference between the cation and the neutral for the linear geometry, given in Section II, was shifted by +0.085 eV in order to obtain the best agreement between observed and calculated PES. This leads to an adiabatic ionization potential of 9.792 eV which is 0.012 eV larger than the value of Garcia et al. 6 Figure 7 shows experimental and calculated spectra. Unlike in H 2 O, all transitions occur near the adiabatic ionization potential. This figure emphasizes that there is a small mismatch between the main peak maximum and the adiabatic ionization potential, which leads to the small discrepancy between the adiabatic ionization potential of this work and that of Garcia et al. 6 Almost no Franck-Condon progressions can be seen in the observed and calculated spectra because for all three electronic states the equilibrium configuration is the linear configuration. When using the bending potential of Jensen and Kraemer 36 for the cation, the agreement with the experimental spectrum is not as good because the main peak in the calculated PES is narrower. A figure showing this calculated PES and the experimental one is available in the supplementary material.

C. CCN radical and CCN + ion

The rovibronic energy levels of the CCN radical and the CCN + cation were calculated as in the previous section. For the radical, Table IV lists the energies up to 1000 cm -1 along with those computed by Hill et al. 26 It can be seen that both sets of values are within 10 cm -1 .

Merer and Travis 38 report for the energy difference of the (010) 2 ∆ and (010) 2 Σ -levels a value of 113.55 cm -1 which compares well with that calculated in this work, 109 cm -1 . For the (020) 2 Φ and (010) 2 Σ -levels, the energies determined by Kohguchi et al. 41 a Levels are identified with the bent molecule vibrational quantum number v 2 and their symmetry species label Γ. Only levels below 1000 cm -1 appear. b Calculated energy in cm -1 reported by Hill et al. 26 For doubly degenerate levels, the average value of the spin orbit components was taken. c Energy in cm -1 calculated in Section V C.

For the CCN + cation, Jensen and Kraemer 36 determined the vibrational energies of the ν 1 2 , 2ν 0 2 , 3ν 1 2 , and 4ν 0 2 states to be 153, 318, 482, and 650 cm -1 , respectively. These values are much larger than those obtained in this work: 104, 208, 318, and 430 cm -1 . Just as for CNC + , this is due to the fact that the bending potential retrieved in this work is much shallower and more anharmonic than that used by Jensen and Kraemer. 36 The photoionization cross sections were calculated as in the previous section and the TPES is shown in Fig. 8. The energy difference between the cation and the neutral for the linear geometry, given in Section II, was shifted by +0.080 eV in order to obtain the best agreement between observed and calculated TPES. This leads to an adiabatic ionization potential of 10.822 eV, which agrees well with the value of Garcia et al. 6 Most transitions are located near this energy. For the calculated spectrum, a Franck-Condon progression can clearly be seen on the high energy side of the main peak and a matching weak feature can be observed in the experimental spectrum. Lowering the temperature leads to a narrower main peak but does not really improve the agreement between observed and calculated spectra on the high energy side of the main peak.

D. HCCN radical and HCCN + ion

Setting n Max , P , and α, α ± to 40, 43, and 20, respectively, rovibronic energies were obtained for HCCN and HCCN + . For the former species, vibrational states are labeled using the linear molecule quantum numbers v 5 l 5 and their energies are given up to v 5 = 4 in Table V. These values should be compared with those in Tables 6 and9 of Koput. 32 Experimental energies of 128.9, 341.7, and 625 cm -1 are available for the 1 1 , 2 2 , and 3 3 states, respectively, and are in better agreement with the present calculated energies than with those calculated by Koput. 32 The vibronic energies calculated for the cation are consistent with a large RT interaction. To our knowledge no spectroscopic investigation of the cation is b Energies in cm -1 calculated in Section V D.

available in the literature.

The photoionization cross section was calculated with Eq. ( 26) for lines with v + ≤ 40, v ′′ ≤ 20, and K + , K ′′ ≤ 30. A Gaussian line shape with a half width at half maximum of 90 cm -1 was taken when computing the TPES, plotted in Fig. 9. About half of the transitions are located near 10.65 eV. The other half are spread over the region from 10.8 to 11.5 eV corresponding to a superposition of many Franck-Condon progressions. Unlike in CNC and CCN, the bending potential of the ground X 3 A ′′ state is quite different from those of the two upper electronic states, leading to a large number of such progressions.

The energy difference between the cation and the neutral for the linear geometry, given in Section II, was decreased by 0.010 eV so as to obtain the best agreement between observed and calculated TPES. The value calculated for the adiabatic ionization potential is then 10.621 eV. This value is 0.118 eV above that reported by Zhao et al. 39 and 0.021 eV below that obtained by Garcia et al. 6 This discrepancy, larger than in the case of CNC, is also due to the fact that the adiabatic ionization potential differs from the experimental energy of the main peak maximum, as emphasized by Figure 9. This latter value was incorrectly taken as the adiabatic ionization potential in Garcia et al. 6

VI. DISCUSSION

Calculated photoelectron spectra are computed for the CNC, CCN, and HCCN radicals based on a derived expression of the total photoionization cross section when one of the electronic states is affected by the RT effects. Two approaches aimed at spectroscopically modeling the quasilinearity and the RT coupling displayed by these species are also derived. The total photoionization cross section obtained in this work, appearing in Eq. ( 21), was derived making use of the results of Willitsch et al. 15 for the total photoionization cross section of an asymmetric top molecule. Their expression, obtained using the orbital approximation, is conveniently expanded in terms of products of rotational, vibrational, and electronic terms. In this work, the RT coupling taking place in one of the electronic states is dealt with by adding the contribution of both RT components of the ionizing electronic transition. The vibrational term is written in terms of Franck-Condon factors involving the bending wavefunctions.

In the approaches accounting for the quasilinearity and the RT coupling, introduced in Section III, the overall rotation and the bending mode are treated together, accounting for the large amplitude nature of this mode. These approaches allow us to describe well the bending-rotation and the rovibronic couplings, but they only provide us with a limited

  icals and their cation were carried out using the CCSD(T) (Coupled-Cluster with Single and Double and perturbative Triple excitations) method explicitly correlated (RCCSD(T)-F12) as well as internally contracted multireference configuration interaction method with Davidson correction (MRCI+Q) with complete active space self-consistent field (CASSCF) wavefunctions. All calculations were performed using the MOLPRO 2012 package and the Dunning augmented triple and quadruple zeta basis. The potential energy curves used for the photoelectron calculations are those calculated using RCCSD(T)-F12.Ground state potential energy curves, displayed in Figs.1, 2, and 3, were computed as a function of a bending angle denoted γ, optimizing all other structural parameters. For the triatomic CNC and CCN, γ was taken as the angle ∠CNC and ∠CCN, respectively; for the tetra-atomic HCCN, γ is the ∠HCC angle corresponding to the ν 5 mode.32 We have chosen to use the curves calculated at the RCCSD(T)-F12 level rather than those at MRCI level because the geometry optimization is performed only at CASSCF for MRCI calculations. Both methods, RCCSD(T)-F12 and MRCI, lead to very similar results. The CASSCF and MRCI calculations were performed at full valence active space for CNC and CCN, namely with 19 (18 for the cations) electrons distributed in 15 orbitals with the 1s orbitals of carbon and nitrogen atoms kept doubly occupied but fully optimized. The CASSCF and MRCI calculations were performed at smaller active space for HCCN due to convergence problems, namely with 20 (19 for HCCN + ) electrons distributed in 13 orbitals with the 1s and 2s orbitals of carbon and nitrogen atoms kept doubly occupied but fully optimized. The optimized value for the geometries at the RCCSD(T)-F12/AVTZ level for CNC, CCN, and HCCN and their cation are listed in

Figure 2 .

 2 Figure 2. Variations with the bending angle γ = ∠CCN of the potential energy function V in eV for the X 2 A ′′ and A 2 A ′ electronic substates of the CCN radical resulting from the RT coupling and for the nondegenerate X 1 A ′ electronic state of the CCN + cation. Ab initio potential values are indicated by dots. Solid lines are the fitted bending potentials calculated with Eq. (1).

Figure 3 .

 3 Figure 3. Variations with the bending angle γ = ∠HCC of the potential energy function V in eV for the X + 2 A ′ and A + 2 A ′′ electronic substates of the HCCN + cation resulting from the RT coupling, top panel, and for the nondegenerate X 3 A ′′ electronic state of the HCCN radical, bottom panel. Ab initio potential values are indicated by dots. Solid lines are the fitted bending potentials calculated with Eq. (1).

Figure 4 .

 4 Figure 4. Variations with the bending angle γ = ∠HCC of the angle β = ∠CCN in degrees for the X 2 A ′′ and A 2 A ′ electronic substates of the HCCN + cation and for the X 3 A ′′ electronic state of the HCCN radical. Ab initio values are indicated by dots. Solid lines are fitted values calculated with Eq. (1).

  linear limit, described by Eqs. (3). For the latter operator, as revealed by Eq. (A2), the singularity arises because of the form of the basis set functions. A P -point Gauss-Jacobi quadrature suited for a weight function of the form (1 -t) α ′ (1 + t) β ′ was used. As shown by Eqs. (A4)-(A6), taking α ′ = α -1 and β ′ = β -1 accounts for both singularities and ensures the best accuracy. Equation (8) puts an upper bound for the node values t i of the quadrature. When P = 21, α = 39, and β = 0, as obtained for the CCN + cation, the largest (smallest) value of t i (γ i ) is 0.40 (66.4 •

  ) under the transformation k → -k, and remembering that the effects of the inversion operation E* of C s on a symmetric top rotational function 48 are E*|N, k⟩ = (-1) N +k |N, -k⟩. Using Wang-type rotational wavefunctions:

  (t) are vibrational operators and V -(t) = [V ′ (t) -V ′′ (t)]/2. Neglecting the geometry variations of the expectation value of L z and L 2 z , the vibrational operators H ±Λ,N,k b (t)take the following form:[22][23][24][25][26][27]

Figure 6 .

 6 Figure 6. Simulated TPES of the X + 2 Π u ← X 1 A 1 ionizing transition of H 2 O. In the upper and middle panel only X + 2 B 1 ← X 1 A 1 and A + 2 A 1 ← X 1 A 1 transitions were considered, respectively.All transitions are included in the lower panel.

  the observed energies reported by Merer and Travis 34 and those calculated in this work.

Figure 7 .

 7 Figure 7. Observed 6 (Calculated) PES of the X + 1 A 1 ← X 2 Π g ionizing transition of CNC in the upper (lower) panel. In the lower panel, the vertical arrow is the adiabatic ionization potential.

  with respect to the (000) 2 Π level are 573.84 and 179.27 cm -1 and also agree well with the present values: 574 and 183 cm -1 .

Figure 8 .

 8 Figure 8. Observed 6 (Calculated) TPES of the X + 1 A ′ ← X 2 Π ionizing transition of CCN in the upper (lower) panel. In the lower panel, the vertical arrow is the adiabatic ionization potential.

a

  Levels are identified using the linear molecule vibrational quantum numbers v 5 l5 .

Figure 9 .

 9 Figure 9. Observed 6 (Calculated) TPES of the X + 2 Π ← X 3 A ′′ ionizing transition of HCCN in the upper (lower) panel. In the lower panel, the small vertical arrow is the adiabatic ionization potential.

  Table I where, when available, they are compared with previous experimental or theoretical values. The calculated energy differences

between the cation and the neutral for the linear geometry, not given in this table, are 9.734, 10.768, 10.574 eV, for CNC, CCN, and HCCN, respectively.

Table I .

 I Experimental and calculated structural parameters a In CNC, the RT coupling leads to lower 2 A 2 and upper 2 B 2 substates; for CCN, to lower 2 A ′′ and upper Figure 1. Variations with the bending angle γ = ∠CNC of the potential energy function V in eV for the X 2 A 2 and A 2 B 2 electronic substates of the CNC radical resulting from the RT coupling and for the X 1 A 1 electronic state of the CNC + cation. Ab initio potential values are indicated by dots. Solid lines are the fitted bending potentials calculated with Eq. (1).

	Species State b	Parameter c Experimental d Theory	Theory (this work) e
	CNC	X 2 Π g	∠CNC	180 34	180 35	180
			r(CN)	1.245 34	1.2534 35	1.2462
	CNC +	X + 1 Σ + g	∠CNC		180 36,37	180
			r(CN)		1.2427, 36 1.2534 37	1.2468
	CCN	X 2 Π	∠CCN	180 38	180 35	180
			r(CC)		1.4045, 35 1.3749 26	1.3821
			r(CN)		1.1889, 35 1.1847 26	1.1847
	CCN +	X + 1 Σ + ∠CNC		180 36,37	180
			r(CC)		1.3803, 36 1.381 37	1.3821
			r(CN)		1.1815, 36 1.2534 37	1.1906
	HCCN	X 3 A ′′	∠HCC	180 33	144.9, 32 137.5, 39 144.95 40 145.36
			∠CCN	180 33	175.4, 32 175.6, 39 175.45 40 175.26
			r(HC)	0.998 33	1.069, 32 1.063, 39 1.0689 40 1.0707
			r(CC)	1.323 33	1.328, 32 1.372, 39 1.3267 40 1.3285
			r(CN)	1.195 33	1.186, 32 1.184, 39 1.1850 40 1.1898
	HCCN + X + 2 A ′ ∠HCC		152.7 39	163.59
			∠CCN		174.1 39	175.92
			r(HC)		1.095 39	1.0842
	et al., 41 0.489.	r(CC) r(CN)		1.317 39 1.212 39	1.2981 1.2048
		A + 2 A ′′ ∠HCC		180 39	180
			∠CCN		180 39	180
			r(HC)		1.092 39	1.083
			r(CC)		1.300 39	1.291
			r(CN)		1.219 39	1.207
	a Bond lengths are in Å and bond angles in degrees.

b

Table II .

 II Calculated rovibronic energies a for H 2 O and H 2 O +

	State b	v 2	K	Cal. c	Cal. d
	X 1 A 1	0	1	39.7	39.2
		1	0	1594.3	1587.6
		1	1	1637.5	1630.3
		2	0	3151.7	3143.2
		2	1	3199.5	3190.7
		3	0	4666.1	4660.9
		3	1	4720.8	4715.5
	X + 2 B 1	0	1	37.3	38.8
		1	0	1417.8	1410.0
		1	1	1459.2	1452.9
		2	0	2784.4	2776.2
		2	1	2831.7	2825.0
		3	0	4097.3	4092.9
		3	1	4153.1	4150.3
	A + 2 A 1	0	0	8379.5	8381.9
		0	1	9033.9	9006.2
		1	0	9884.6	9901.0
		1	1	10647.1	10654.8
		2	0	11564.3	11605.8
		2	1	12498.4	12563.9

Table III .

 III Observed and calculated energies a for CNC

	v 2	Γ	Obs. b	Cal. c
	0	Π g	-24	-24
	1	Σ -u	120	117
	1	∆ u	248	245
	2	Π g	356	359
	1	Σ + u	473	479
	2	Φ g	496	500
	2	Π g	855	802
				1 2 , and 4ν 0

Table IV .

 IV Calculated rovibronic energies a for CCN

	(0v 2 0)	Γ	Cal. b	Cal. c
	(000)	Π	0	0
	(010)	Σ -	182	183
	(010)	∆	293	292
	(020)µ	Π	440	441
	(010)	Σ +	459	462
	(020)	Φ	575	574
	(030)	Σ -	669	672
	(030)µ	∆	713	715
	(020)κ	Π	824	832
	(040)µ	Π	924	929
	(040)µ	Φ	993	996

Table V .

 V Calculated vibrational energies a for HCCN

A ′ substates; and for HCCN + to lower 2 A ′ and upper 2 A ′′ substates. c For HCCN and HCCN + the dihedral ∠HCCN is 180 • , consistent with trans conformations. d For HCCN, Brown et al.

assumed a linear geometry. e Obtained at the RCCSD(T)-F12/AVTZ level of theory. respectively. The value obtained for the Renner parameter ϵ, 0.424, is in reasonable agreement with that reported by Merer and Travis, 38 0.44, but in better agreement with that calculated by Hill et al., 26 0.429. It is in poor agreement with that retrieved by Kohguchi

states, the energies reported by these authors are 165,

349, 530, and 720 cm -1 , respectively,

by adding the contribution of the RT components X + A ′ ← X Σ and X + A ′′ ← X Σ distinguished using + and -signs, respectively. The resulting total and averaged photoionization cross section can be obtained making a few changes in Section IV A. In Eqs. (21) and (26), the electronic terms B l ′′ λ ′′ . The Franck-Condon factor (q ± v ) 2 should be evaluated making the substitution v ′′ , N ′′ , K ′′ ↔ v + , N + , K + in Eq. (19). In the rotational factor Q(l ′′ ) of Eq. ( 20), the substitutions δ ′′ → ±δ ′′ and δ + → ±δ + should be made. At last in Eq. ( 22), g(X + Σ) becomes g(X Σ).

For the molecules dealt with in the present investigation, the photoelectron is ejected from a π orbital. This means that λ ′′ = 1. In Eqs. ( 21) and ( 26), the first term in the summation over l ′′ will be assumed to be the dominant term and this summation will be restricted to the l ′′ = 1 term. It will also be assumed that electronic terms B

of Section IV A are equal.

V. NUMERICAL RESULTS

Using the two previous sections, photoionization cross sections are calculated for the ionizing transitions of H 2 O, CNC, CCN, and HCCN. The first molecule provides us with a test of the theoretical approach and also allows us to understand the contribution of the various terms in Eq. ( 26) to the TPES. For the three remaining molecules, comparisons with the experimental photoionization spectra 6 are carried out. With the experimental setup used in this reference, the energy resolution is 17 meV and the rotational temperature, fairly well known, is around 200 K. The vibrational temperature is, however, not as well defined and can be quite high since the consecutive R-H + F --→ R + HF reactions are exothermic.

Below, we have chosen a temperature of 500 K for the simulations of the CNC, CCN, and HCCN spectra.

A. Neutral H 2 O and H 2 O + cation

Starting from the potential energy surfaces of Partridge and Schwenke 28 for the neutral species and of Wu et al. 29 for the cation, bending potentials and bond lengths were retrieved for the ground stretching states. The results of Sections III A and III B were then applied to the calculation of the rovibronic energies of H 2 O and H 2 O + . A sufficient accuracy was description as they ignore the electron spin and the stretching modes. This approximate description should nonetheless be sufficient for the qualitative results sought in this paper. A more satisfactory description would be achieved with treatments accounting for all vibrational modes. 27,51,52 In order to understand qualitatively the contribution of the various terms in the new expression of the total photoionization cross section, the theoretical treatment was first applied, in Section V A, to the simulation of the photoelectron spectrum of H 2 O for which both the neutral 28 and the cationic 29 species are spectroscopically well characterized; the latter displaying a strong RT coupling with a Renner parameter larger than 1. The calculated TPES, displayed in Fig. 6, is in qualitative agreement with the experimental one. 30 The theoretical approach was then applied to the computation of the photoelectron spectra of CNC, CCN, and HCCN. The bending potentials of these species were retrieved using the ab initio calculations presented in Section II and their rovibronic energies, computed in Sections V B, V C and V D, turned out to be in good agreement with spectroscopic data 34,38,41 and calculations. 26,32 Adjusting the adiabatic ionization potential of each molecule, a qualitative agreement between their experimental 6 and calculated PES or TPES was achieved and can be seen in Figs. 789. For CNC, both spectra display a single peak. Matching the energy of these peaks, the adiabatic ionization potential, deduced from Section V B, was increased by 0.085 eV and was found to be 9.792 eV. For CCN, both spectra span a larger energy range and the improved value of the vertical ionization potential, 0.080 eV higher than the value from Section V C, is 10.822 eV. For HCCN, both spectra display a narrow peak and a broad feature. Matching these spectra leads to an adiabatic ionization potential of 10.621 eV, which is 0.010 eV smaller than the value calculated in Section V D.

Appendix A: Matrix elements of inverse inertia tensor components

The matrix element of the diagonal component µ zz (t) of the inverse inertia tensor between two basis set functions of Eq. ( 4) is equal to:

The matrix element of the kinetic energy operator P t µ tt (t)P t can be obtained from:

where, in agreement with Partridge and Schwenke, (t) as:

Evaluating the matrix elements in Eqs. (A1) and (A2) with a P -point Gauss-Jacobi quadrature suited for the weight function (1 -t) α-1 (1 + t) β-1 , the quadrature leads to:

where ω i and t i are respectively the weights and nodes. When Op = µ zz (t), the function

When Op = P t µ tt (t)P t , the function f (t) is:

Equations ( 3) and (A3) show that in both cases f (t) is a finite smooth function of t and the quadrature should lead to accurate results.

SUPPLEMENTARY MATERIAL

See supplementary material for a PDF file containing four figures numbered S1 to S4.


Christophe Fond

Shape of the surface of a liquid under indentation by a flat punch -effect of surface tension for a semi-infinite axisymmetric media in the presence of gravity

Keywords: nanomechanics, indentation, surface tension, contact, liquid, incompressible, gravity

published or not. The documents may come

Shape of the surface of a liquid under indentation by a flat punch -effect of surface tension for a semi-infinite axisymmetric media in the presence of gravity

Preamble

This article follows two articles, in English and French, concerning the indentation of an elastic medium for which surface tension acts considerably. In the situation where the elasticity of the medium disappears in front of the surface tension, as is the case for microscopic indenters and compressible soft media, the problem tends towards the indentation of a stretched membrane (Fond, 2018c), (Fond, 2018d). In the situation where the medium is almost incompressible and very soft in shear, the problem tends towards that of a liquid. In finite medium the solutions are analyzed and it is shown in (Fond, 2018a) et (Fond, 2018b) that gravity can be neglected, always in the case of microscopic indenters. For semi-infinite environments, it will be shown further on that the effect of gravity can no longer be neglected. Considering a semi-infinite media and the attempt to superimpose the effects of the response of the bulk [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF] and that of the surface, it is therefore necessary to know the solution for these two situations considered independently. It has been shown that the insertion of the punch into a liquid tends towards infinity when gravity is neglected. The aim here is to provide a suitable solution with consideration of gravity, a case that does not admit an analytical solution. Many articles have been dedicated to this subject of indentation, for example [START_REF] Chakrabarti | Direct measurement of the surface tension of a soft elastic hydrogel: Exploration of elastocapillary instability in adhesion[END_REF], [START_REF] Cortat | How closely can a solid approach an air-water surface without becoming wet?[END_REF], [START_REF] Cortat | Using stable and unstable profiles to deduce deformation limits of the air-water interface[END_REF], [START_REF] Butt | Force measurements with the atomic force microscope: Technique, interpretation and applications[END_REF], [START_REF] Colchero | Observation of liquid neck formation with scanning force microscopy techniques[END_REF], [START_REF] Cappella | Force-distance curves by atomic force microscopy[END_REF], (Forcada et al., 1991), the article closest to the present analysis being that of [START_REF] Chan | Forces between a rigid probe particle and a liquid interface: I. the repulsive case[END_REF].

Introduction

The liquid is supposed to be incompressible. The static situation is considered, i.e. there is no inertial effect. An established and stable state is considered so that the viscosity of the fluid is neglected. The surface tension γ, assumed constant, balances a local pressure pρgδ exerted on the surface at a horizontal distance r from the center of the flat punch (Fond, 2018b). The static equation of equilibrium in Fig. 1 is therefore formulated by :

p -ρgδ = -γ( sin (atanδ ′ ) r + δ" (1 + δ ′2 ) 3/2 ) (1)
where p is the pressure under the indenter, ρ the density of the liquid, g the acceleration of gravity and δ the height of the surface at the distance r. For a classic free surface problem, the maximum height δ asympt. corresponds to an infinite radius from where p = ρgδ asympt. . This equation can be rewritten

p -ρgδ = -γ( δ ′ r √ 1+δ ′2 + δ ′ (1+δ ′2 ) 3/2
). For calculation convenience it was chosen δ = 0 at the base of the punch. If the zero height reference is the free surface then the punch base will be at the depth corresponding to -δ asympt. . This change of origin will lead to the depth δ(r) -δ asympt. . The equilibrium of the flat punch is given by (see Fig. 2) :

F = 2πaf + πa 2 p (2)
where f denotes the Laplace's tension at the corner of the punch given by f = γsin(β(a)).

Numerical integration of the equilibrium equation

Since the equilibrium equation 1 does not admit an analytical solution, a calculation by numerical integration has been developed. A distance must be chosen D >> a that ensures to be on the numerical asymptote1 . The line segment defined by the two points a and D at z = 0, i. e. (a, 0) -(D, 0), is discretized into sub-segments (r i , 0) -(r i+1 , 0), i ranging from 0 to n with r 0 = a and r n = D. Noticing b a = log 10 (a) and b D = log 10 (D), the spatial discretization is done so that r i = 10 ba+(i/n)(b D -ba) . The surface is composed of line segments (r i , δ i ) -(r i+1 , δ i+1 ). Given an angle β(a), the first segment is given by (r 0 , 0) -(r 1 , δ 1 ) where δ 1 denotes (r 1 -a)(tan(β(a))). In order to know the following points (r 2 , δ 2 ) then (r i , δ i ) a dichotomy procedure looks for the value of δ i which respects the equation. 1 with a precision based on the value abs(

δ m i -δ m-1 i δ i-1 -δ i-2 )
where m denotes the iteration by dichotomy to find δ i . For the calculations presented here the precision chosen is 10 -6 and n = 1000 for the number of segments. The curves δ ′ and δ" are calculated from the second degree polynomial which passes through the three points (r i , δ i ), (r i+1 , δ i+1 ) and (r i+2 , δ i+2 i + 2). The slope δ ′ is the slope of the polynomial function in the middle of the segment. 2Whatever the force exerted on the object, as long as β(a) < π/23 , the surface disturbance is typically exerted at a distance of 1cm as shown in Fig. 3 for γ = 0.07d/m 2 , ρ = 1000kg/m 3 , g = 9.81m/s 2 , a = 10nm. For values of r large enough that sinβ(r) ≈ sinβ(r) we find the values predicted by the approximate analytical solutions and the curves normalized by sinβ(a)4 are superimposed.

Fig. 4 shows that the reaction force to indentation varies almost linearly with the indentation depth. Indeed, less than 6% error is made on the prediction of force by considering that the original slope is kept up to sinβ(a) = 0.9999, i. e.β(a) = π/2.02. For the parameters chosen for the illustration, the contribution of hydrostatic pressure under the punch, term πto 2 p of the equation 2 where p = ρgδ asympt. , is negligible.

A dimensional analysis allows to adimension by introducing the capillary length λ = γ ρg (see Appendix A). The Fig. 5 provides the results of the numerical integration described above for sinβ(a) = 0.55 . Around r = λ the effect of gravity is felt and it goes from a trend δ(r) = asinβ(a)log(r/a) to the asymptote δ ≈ 0 for r > 6 5 λ. Whatever the force exerted on the object, the surface disturbance is typically less than 10λ.

For the intended use, essentially the maximum depression is sought. From the previous equations it comes:

δ asympt. ≈ a sinβ(a) log(6λ/5a) (3)
or even: so that:

δ asympt. ≈ F 2πγ log(6λ/5a) (4)
δ asympt. max < a log(6λ/5a) (5)
for a microscopic flat punch floating on a liquid of similar property to water on earth, the contribution to the pressure balance under the punch is always negligible compared to the contribution of surface tension.

Approximative function

The analysis of the solutions of the second order differential equation in Annex Appendix A provides information to find a suitable approximate function. In addition, a dimensional analysis highlights the term ρg γ commonly referred to as capillary length. For the intended use, essentially the maximum depression is sought and a slight error in the shape of the surface will be tolerated.

δ(r) ≈ asinβ(a) 2e -r/Λ 1 + e -r/Λ log(r/Λ) (6) où Λ ≈ 6 5 γ ρg .
Rappelons que F = 2πγa sinβ(a). The approximate function given by the eq. 6 and illustrated in Fig. 6 shows an underestimation of the

-7×10 -8 -6×10 -8 -5×10 -8 -4×10 -8 -3×10 -8 -2×10 -8 -1×10 -8 0 1×10 -8
1.0e-08 1.0e-07 1.0e-06 1.0e-05 1.0e-04 1.0e-03 1.0e-02 1.0e-01 sink in the vicinity of the asymptotic value and it should be considered that the depth -δ(r) < 0 ∀r. Note that Λ should be chosen slightly higher than λ for a better approximation.

δ(r) - δ asympt. (r) (m) r (m) legend λ=1.8mm a=10nm Fit Λ = λ Fit Λ = 6/5 λ

Case β = π/2 3.2.1. Non-submerged object

Fig. 7 illustrates the limit case for the angle at the corner of the punch. Let's first consider the case where no part is submerged, i.e. δ imm. = 0. When β = π/2, the contribution to the equilibrium of the surface tension is saturated and is worth f = 2πaγ. One of the two radii of curvature is exactly a, the other is at least a if p ≥ 0. If the pressure p under the flat punch is negligible then the sum of the curves must cancel each other out. Hence, the two radii of curvature are worth a. For the above numerical calculations for which the object is microscopic, i.e. a << λ, and sinβ = 0.9999, a simple geometrical calculation shows that the difference for δ asympt. to the sinβ = 1 situation is about 0.014a. In view of the results shown in Fig. 5, this value appears negligible.

Submerged object

When δ imm. > 0, if the object has been immersed monotonously I assume that β(a) = π/2, then the maximum force that can be balanced is worth:

F = 2πaγ + πa 2 ρg(δ asympt. + δ imm. ) (7)
At the limit of immersion, using the eq. 5 valid for typically a < λ 10 , we get:

F = 2πaγ + πa 3 ρglog(6λ/5a) (8)
Fig. 8 is a numerical computation for microscopic objects such as a < λ for water at room temperature, i. e. ρ = 1000kg/m 3 , g = 9.81m/s 2 and γ = 0.07 hence λ = 2.67mm. It is clear that the contribution of the pressure under the punch is always negligible compared to that of the surface tension. On the other hand, the pressure related to gravity is used to calculate δ asympt. and renders finite the depth considering an infinitely large liquid surface.

Objects with a size greater than the capillary length λ

The objective was achieved for objects of microscopic size for values of ρg typically in the order of 10 +4 P a and surface tensions typically in the order of 5 10 -2 N/m. Nevertheless, for a better understanding let's analyze the situation where the object has a dimension greater than λ. It is obvious that the approximate solution provided in eq. 6 can no longer be satisfactory for objects larger than λ. Let's illustrate this situation in cases where ρ = 10 +3 kg/m 3 , g = 9.81m/s 2 and γ = 0.07N/m 6 , i. e. λ = 2.67mm for a = 6 typically the case of water at room temperature on earth 8mm ≈ 3λ and a = 1dm ≈ 37λ in Fig. 9. To ensure proper integration of the differential equation, an angle of β(D) to the horizontal, very small, i. e. β(D) << 1, is given at a distance of D far enough from the edge of the punch, located at r = a, to be considered on the horizontal asymptote. Given this perturbation, integration is processed from D to a. When the depth δ(a) and the angle β(a) have been determined by this way, the integration from a to D is performed and it is checked that the same surface profile is restored in both directions of integration. For the examples in Fig. 9, the perturbation was chosen so that the slope at the edge of the flat punch is maximum, i. e. sinβ(a) = 1. Note that the force and depression increase almost linearly with this perturbation angle β(D).

We notice in Fig. 9 that the disturbance length of the surface 7 is about λ from the punch edge. The contribution of the hydrostatic pressure below the punch surface to the total reaction force is about 64% for a = 8mm and 96% for a = 1dm. It has been shown previously that this is negligible for microscopic 8 but in these cases it is considerable. The capillary length λ for the a radius therefore marks the transition between the sustentation by the 7 length from the punch to find a nearly horizontal surface 8 objects whose order of magnitude is much smaller than the order of magnitude of λ surface tension and the sustentation by the hydrostatic pressure generated under the tip by its depth. When β(a) has reached π/2, if the punch continues to sink the "triple point" moves along the punch wall and the object floats or sinks.

Discussion

The aim was to provide a simple solution for the case of indentation of a flat punch in a semi-infinite planar medium, in the case where the mechanical behaviour of this medium is similar to that of a liquid. The medium is considered to be liquid, of negligible viscosity and incompressible. Calculations confirm that the latter assumption is appropriate. Indeed, the extremely small pressure variations are not likely to considerably change the density of the liquid whose compressibility is typically in the order of 2GP a. The results provided here can be easily used and adapted for spherical or conical indenters for which the angle β(a) will be limited by the geometry of the indentor. Indeed, since it has been shown here that for microscopic objects the contribution of pressure under the indent to equilibrium is a second order term, it will not be necessary to redeploy numerical integrations in these cases. Only objects with a capillary length greater than λ would require a more precise quantification, but this does not concern the purpose of this article.

Conclusion

The calculations presented in this article for a liquid are consistent with those of (Fond, 2018c), (Fond, 2018a) for highly stressed membranes. They make it possible to locate the domain of validity in order to neglect certain parameters. They provide validation arguments for finite elements numerical indentation computations taking into consideration surface tension as well as elasticity for microscopic indentation. Given the relatively large volume of information required to properly present this numerical model, it will be presented in a separate article. This finite element numerical model calculates the effects of elasticity and surface tension in a coupled manner by covering the entire range of the ratio γ aµ where µ refers to the shear modulus of the solid. It was therefore essential to have solutions in cases where γ aµ << 1 and γ aµ >> 1.

Figure 1 :Figure 2 :

 12 Figure 1: Pressures and curvatures for an axisymmetric model (left) and gravity-related pressure variation (right). a F

  1e-07 1e-06 1e-05 1e-04 1e-03 1e-02 1e-

Figure 3 :

 3 Figure3: Normalized surface profiles for γ = 0.07d/m 2 , ρ = 1000kg/m 3 , g = 9.81m/s 2 , a = 10nm and sinβ(a) = 0.1, 0.5, 07 and 0.9. At the top left is a zoom on the fairly non-linear part for which typically sinβ(r) ≈ β(r).

Figure 4 :

 4 Figure 4: Force versus indentation depth for γ = 0.07d/m 2 , ρ = 1000kg/m 3 , g = 9.81m/s 2 , a = 10nm. The dotted line shows the slope at the origin. The relative error is calculated with the values of this line.

Figure 5 :

 5 Figure 5: Results of numerical integration and normalization by λ =

Figure 6 :

 6 Figure 6: Approximate function to describe the indentation depth -δ as a function of distance r for sinβ(a) = 0.5 and λ = 1.8mm.

Figure 7 :

 7 Figure 7: Flat punch equilibrium.

Figure 8 :

 8 Figure 8: Reaction force versus punch radius. The dotted lines indicate the contributions of surface tension and pressure under the punch.

Figure 9 :

 9 Figure 9: Surface profiles for a = 8mm and a = 1dm corresponding to λ = 2.67mm. The punches are materialized by the black dotted lines.

for the results presented here in double precision, i. e. to the nearest 10 -15

The procedure can also do the "reverse path" from D to a from a perturbation defined in D by (r n , δ n ) and (r n-1 , δ n-1 ) but in this direction it is not possible to predict a priori sinβ(a).

for β(a) → π/2, the contribution of surface tension can no longer increase and the object would sink into the liquid until the pressure under the punch can balance the force exerted. This situation is discussed below.

remind us that normalizing by sinβ(a) amounts to normalizing by the applied force since the contribution of surface tension to equilibrium is f = 2πγsinβ(a) and that F ≈ f if πa 2 p << f .

the maximum force will be reached for sinβ(a) → 1

Appendix A. Analyse des solutions de l'équation différentielle du second ordre

The equilibrium equation 1 does not admit an analytical solution if all its terms are non-zero. In order to find a suitable approximation of numerical solutions, it is convenient to analyze analytical solutions when certain terms of the equation of equilibrium vanish.

Appendix A.1. ρg = 0

For small angles, i. e. sinβ(a) ≈ β(a), when ρg = 0, there is an analytical solution to the equilibrium equation. For p = 0 it is formulated by:

with r ≥ a.

Appendix A.2. δ ′ = 0 When δ ′ = 0, there is an analytical solution to the equilibrium equation formulated by: δ

where c 1 and c 2 are constants.

Appendix A.3. δ" = 0 When δ" = 0, there is an analytical solution to the equilibrium equation formulated by: δ

where c 3 is a constant.

Fond, C., 2018b. Forme de la surface d'un liquide ou d'un solide se comportant quasiment comme un liquide sous l'action de l'indentation par un poinçon plat -effet de la tension de surface.

Fond, C., 2018c. Indentation of a highly prestressed membrane by a flat punch: analytical solution.

Fond, C., 2018d. The shape of the surface of a liquid or a solid that behaves almost like a liquid under the action of indentation by a flat punch -the effect of surface tension.

Forcada, M.L., Jakas, M.M., Gras-Marti, A., 1991. On liquid-film thickness measurements with the atomic-force microscope. The Journal of Chemical Physics 95, 706-708. https://doi.org/10.1063/1.461420.