For conical intersections of two states (*I,J* = *I +* 1) the vectors defining the branching or g-h plane, the energy difference gradient vector **g**^{I,J}, and the interstate coupling vector **h**^{I,J}, can be made orthogonal by a one parameter rotation of the degenerate electronic eigenstates. The representation obtained from this rotation is used to construct the parameters that describe the vicinity of the conical intersection seam, the conical parameters, s^{I,J}_{x} (**R**), s^{I,J}_{y} (**R**), g^{I,J}(**R**), and *h*^{I,J}(**R**). As a result of the orthogonalization these parameters can be made continuous functions of **R**, the internuclear coordinates. In this work we generalize this notion to construct continuous parametrizations of conical intersection seams of three or more states. The generalization derives from a recently introduced procedure for using non-degenerate electronic states to construct coupled diabatic states that represent adiabatic states coupled by conical intersections. The procedure is illustrated using the seam of conical intersections of three states in parazolyl as an example.

## I. INTRODUCTION

Recently there has been renewed interest in the role of accidental three state conical intersections in nonadiabatic processes.^{1–13} Studies dealing with: (i) the location of points of conical intersection;^{4,5,14} (ii) the seam and branching spaces which are of dimension *N*^{int}-5 and 5, respectively;^{15} (iii) the degeneracy space, a subspace of the branching space in which the degeneracy of two states is preserved;^{16} (iv) the geometric phase;^{1,9,17,18} and (v) dynamical studies of actual^{11,12} and model^{13} systems, have been reported.

The description of the branching, or *g-h* plane for intersections of two states of the same symmetry is key to interpreting and determining the dynamics near two state conical intersections. As described in Sec. II, see also Ref. 19, for two state intersections the branching plane is defined by two vectors **g**^{I,J}, the energy difference gradient vector and **h**^{I,J}, the interstate coupling vector. Since there are two degenerate states at a point of conical intersection these states are defined up to a one parameter rotation. When this flexibility is used to make the **g**^{I,J} and **h**^{I,J} vectors orthogonal,^{20} the parameters that define the topography of the intersection seam, the conical parameters^{21} |$s_x^{I,J}$| $ s x I , J $ (**R**), |$s_y^{I,J}$| $ s y I , J $ (**R**), *g*^{I,J}(**R**), and *h*^{I,J}(**R**) can be made continuous functions of **R**, the atom centered cartesian coordinates. This continuity is essential for the polynomial representations of diabatic states we use to describe adiabatic states coupled by conical intersection seams.^{22} The orthogonalization procedure, which defines the arbitrary rotation of the two degenerate states, can be exploited in algorithms for locating conical intersections.^{23} From a pedagogic perspective the rotated **g**^{I,J} and **h**^{I,J} vectors reflect the relevant molecular symmetry an attribute which is usually absent in the original or nascent vectors. In this work we generalize these ideas to conical intersections of three or more states. Section II outlines the theory. Section III demonstrates its application to pyrazolyl, an azolyl for which previous work has demonstrated the existence of low-lying three state intersections.^{5,24,25} Section IV summarizes and concludes.

## II. THEORETICAL APPROACH

### A. Background: Two state degeneracies

Let the **c**^{J}(**R**) satisfy the electronic Schrödinger equation in a configuration state function (CSF) basis,

where **R = R**^{0} + δ**R**. In the vicinity of an *N* state conical intersection at **R**^{0}, taken to be the zero of energy, through first order in nuclear coordinates, Eq. (1) becomes^{6}

where

As described in Ref. 26 (see also Ref. 27) the **c**^{I}(**R**^{0}) are used as a geometry independent transformation of the generally applicable **H**^{CSF}, to **H**^{d} an *N*x*N* matrix whose spectrum agrees with that of **H**^{CSF} through first order in δ**R** near **R**^{0}. For both **H**^{d} and **H**^{CSF} the *N* degenerate states are determined up to an [*N*] ≡ *N*(*N* − 1)/2 parameter orthogonal geometry dependent transformation, **U**. This indeterminacy makes the **h**^{I,J}(**R**) discontinuous functions of **R**. We have previously discussed^{20} the two state conical intersection case where *N* = 2 and [*N*] = 1. In this case we introduced the orthogonality requirement,

where

For

Eq. (3a) leads to the requirement

where g^{I,J} = *∥***g**^{I,J}∥ and h^{I,J} = *∥***h**^{I,J}∥. The requirement (3a) is a natural choice since in that case **H**^{d} has the appealing form

where **s**^{I,J} = (**h**^{I,I} + **h**^{J,J})/2, δ*R*^{g} = **g**^{I, J}/*g*^{I, J} · δ**R**, and δ*R*^{h} = **h**^{I, J}/*h*^{I, J} · δ**R**.

### B. N state degeneracies

The generalization to three or more states is not obvious since there are not enough parameters in **U**^{(N)} (only [*N*]) to make [*N* + 1] − 1 mutually orthogonal linear combinations of the **h**^{I,J}. Here [*N* + 1] − 1 appears since the trace of **H**^{d} is excluded as it does not contribute to the state splitting. In this work we show that for an *N* state degeneracy of the states *J, J*^{(min}^{)} ≤ *J ≤ J*^{(max)} = *J*^{(min)} + (*N* − 1), the [*N*] conditions,

determine **U**^{(N)} and lead to continuous **h**^{I,J}(**R)**. A condition similar to Eq. (6) was introduced and carefully analyzed in our recent work describing the representation of adiabatic potential energy surfaces coupled by conical intersections.^{28} Here we provide a concise recapitulation of that analysis (Eqs. (7a)–(7c), (8), (9a)–(9c), (10a), (10b), (11a)–(11d), (12a)–(12c), and (13)–(16)) sufficient to explain the present procedure.

The desired orthogonal transformation **U**^{(N)} is constructed as a product of 2 × 2 rotations **U**^{N,i,j}(θ^{i, j}) (see Sec. III C),

where **U**^{N, i, j}(θ^{i, j}) is the Givens rotation matrix between basis *i* and *j*:

except for

Below we will suppress the *i,j* superscript on θ^{i, j} when no confusion will result. Instead of dealing directly with Eq. (6) we showed that those equations are equivalent to minimization of the function,

which is comprised of [*N*] (*K*,*L*) pairs.

To establish that equivalence, we proceed as follows. Denote the wave function of electronic state *I* in any basis, such as the CSF basis, as **d**^{I}, and the corresponding wave function after the Givens rotation as **d**^{I}(*θ*). Then the effect of one individual Givens rotation on the wave functions can be written in vector form as

so that the rotated vectors satisfy

Therefore,

where *K* ≠ *I*, *J*. Other **h** vectors remain unchanged.

For the current **d**^{I} to satisfy Eq. (6)**g**^{I,J} · **h**^{I,J} must vanish for all pairs. From Eq. (10a), Eq. (6) is equivalent to the condition that *S* is stationary with respect to any rotation, that is for any pairs of states *I* < *J*,

Since *S*, from Eq. (8), is a sum of norms, it is non-negative. Being a non-negative continuous function, *S* is guaranteed to have a minimum. The solution to Eq. (6) is therefore also guaranteed.

This solution can be obtained by a sequence of 2 × 2 rotations.^{28} In order to show this, *S*(*θ*) is first divided into 3 parts:

and

Since the last two parts are both constant with respect to *θ*, we only need to find the minimum of the first term ‖**h**^{I, J}(θ)‖^{2} with respect to angle *θ*.

Defining *C* and *α* as

Eq. (11b) becomes

Note that using Eqs. (12a)–(12c), Eq. (4b) is equivalent to

which has solutions

where we have chosen the values of θ for which cos(α + 4θ) = −1.

On the other hand, using Eq. (13), this choice of θ yields the relation

where min() denotes the smaller of the two arguments. The ≤ in Eq. (16) becomes an equality only when (**h**^{I, J} · **g**^{I, J})^{2} = 0. Thus (*h*^{I, J})^{2} and hence *S* are strictly decreasing unless **h**^{I, J} · **g**^{I, J} = 0. By choosing the pairs of states *I* and *J* with the largest rotation angle every iteration, the rotation will always lower the value of *S* until **h**^{I, J} · **g**^{I, J} = 0 holds for all pairs of states *I* and *J*, which is the desired result. Using Eq. (9a) and the definition of **h**^{J,K}, **d**^{I} is updated each iteration, until the solution to Eq. (10b) is obtained. As discussed below the resulting **h**^{i,j} are used to define the rotated branching space.

Equation (6) may have more than one solution corresponding to local solutions of Eq. (10b). Further changing the phase of a single eigenstate, a reflection not included in Eqs. (7a)–(7c) which changes **h**^{I,J} but not **g**^{I,J}, can change the final value of *S*. This perhaps counterintuitive result was verified by explicit numerical calculations. These two issues are attributable to the range of starting values of the **g**^{I,J} and **h**^{I,J}. Methods to deal with these situations are described in the Sec. III.

### C. The branching or g-h space

At **R**^{0} a point of *N* state conical intersection **H**^{d} in Eq. (2a) is, to first order in δ**R**, a rigorous diabatic representation. See Eq. (2b). However, it is inappropriate for defining the branching space of this conical intersection since it does not explicitly account for the trace. For points of conical intersection of 3 states, to be considered in Sec. III, the g-h or branching space is defined by the five vectors:^{8}

Also required is the trace vector

which does not contribute to the state splitting and therefore is not a member of the branching space. In this representation, **H**^{d} is given by^{8}

where

Note that **H**^{d} does not have the simple form of Eq. (5) since here only **h**^{[}^{1]} and **h**^{[4]} are orthogonal. The continuous representation of the branching or g-h space vectors using Eq. (10b) for a practical example is discussed below. For the *N* state generalization of Eq. (18) see Ref. 8 and also Ref. 29. These works use Lie group techniques to process the **B** matrices.

## III. APPLICATION TO PYRAZOLYL

The algorithm described in Sec. II B is used here to develop a continuous representation of the branching space for a portion of the three state conical intersection seam in pyrazolyl (CH)_{3}N_{2}, Ref. 5, pictured in Fig. 1. This seam has been implicated in the photoelectron spectrum of the pyrazolide anion, which has been simulated previously.^{24,25} However to date only isolated points on the seam of three state conical intersections have been discussed.

### A. Electronic structure description

The low-lying 1,2,3^{2}A states of pyrazolyl exhibit a seam of three-state conical intersections whose description is inherently multi-reference in character and so multi-reference configuration interaction singles and doubles (MR-CISD) wave functions will be employed to describe this seam. The active space for these demonstration calculations is comprised of three orbitals including the two *π* type orbitals with one additional nodal plane perpendicular to molecular plane, as well as the anti-symmetric linear combination of the two lone pair orbitals one on each of the nitrogen atoms, yielding a 5 electron/3 orbital complete active space (CAS) reference wavefunction, with 3 CSFs that correspond to the 3 electronic states. The molecular orbitals were obtained from a state-averaged multiconfigurational self consistent field (MCSCF) procedure that included all three electronic states with equal weights and employed polarized double-zeta^{30} atomic basis sets. This basis set, in conjunction with the second-order CI treatment of electron correlation and employing a generalized interacting space restriction approximation, yielded MRCI wavefunctions comprised of 1 127 824 CSFs. All calculations reported in this work were carried out using the COLUMBUS suite of electronic structure codes.^{31–33} All energies are reported relative to the ground state minimum of pyrazolyl for which E_{1} = −224.7651475041 a.u.

### B. Three state seam of conical intersection and the g-h space

The minimum energy three state conical intersection structure for pyrazolyl, **Q**_{0} with energy E_{1}(**Q**_{0}) = 4905.19 cm^{−1}, is pictured in Fig. 1. Starting from this point a number of additional points on the seam of three state intersections were determined along each of three paths: (1) the stretching of the C2-H7 bond, (2) the out-of-plane angle of H6 with respect to the molecular plane, and (3) the out-of-plane angle of H7. We desire conical intersections systematically displaced along the path parameter but otherwise randomly placed. Such paths were determined by incrementing the parameter describing the path and then running the conical intersection search without energy optimization. The produced three state conical intersections are degenerate to <0.5 cm^{‑1}. The degenerate energies along these three paths are plotted in Figs. 2(a)–2(c), respectively, to demonstrate their continuity. Nascent **g**, **h,** and **s** vectors and the corresponding vectors after the solution of Eq. (10b), indicated by superscript ∼, are juxtaposed in Figs. 3 and 4, respectively. Figs. 3(a)–3(f) (Figs. 4(a)–4(f)) report the nascent vectors **h**^{[j]}, j = 1–6 (the rotated vectors |${\bf \tilde h}^{[j]}$| $ h \u0303 [ j ] $ j = 1–6) at the first point along path (3). This point is chosen because while it is nearly identical to the minimum energy 3 state crossing point, **Q**_{0}, which has C_{2v} symmetry, its geometry carries only C_{1} point group symmetry, making it a more suitable example of a general point without symmetry. The differences between the corresponding vectors in these figures are seen to be dramatic, except for **h**^{[6]} which is invariant to any rotation. The orthogonality of |${\bf \tilde h}$| $ h \u0303 $ ^{1,2} and |${\bf \tilde g}$| $ g \u0303 $ ^{1,2}, which is expected, is clearly evident. Also evident is the near C_{2v} symmetry of the rotated vectors, |${\bf \tilde h}^{[j]}$| $ h \u0303 [ j ] $ j = 1–6, in Fig. 4 which again is expected. No evidence of this symmetry is present in Fig. 3 except of course for **h**^{[6]}.

### C. Continuous representation of the seam

Figs. 5–7 consider the continuity of the vectors describing the branching space along the seam paths (1), (2), (3), respectively. At each seam point **Q**_{i} these figures report (as individual markers) ||**h**^{[j]}(**Q**_{i}) − **h**^{[j]}(**Q**_{0})|| the norm of the difference between the nascent vectors **h**^{[j]},j = 1–6 and the corresponding reference vectors, **h**^{[j]}**(Q**_{0}**)**, as well as (solid lines) |$||{\bf \tilde h}^{[j]} ({\bf Q}_i) - {\bf \tilde h}^{[j]} ({\bf Q}_0)||$| $ \Vert h \u0303 [ j ] ( Q i ) \u2212 h \u0303 [ j ] ( Q 0 ) \Vert $ the difference norms for the corresponding rotated quantities. The continuity (lack of continuity) is clearly evident in the rotated (nascent) quantities. For each of the solid lines the difference norm is small and a continuous function of the displacement. For the individual markers the difference from the origin vector is a discontinuous function of the displacement, except for **h**^{[6]}. **h**^{[6]} = |${\bf \tilde h}^{[6]}$| $ h \u0303 [ 6 ] $ and are continuous as is expected of the trace.

A comment on the irregular behavior of the nascent states and the corresponding **g** and **h** vectors is in order. The “discontinuities” pictured in Figs. 5–7 are a necessary consequence of the search algorithm and the conical topography. Since exact degeneracies are never achieved with a numerical procedure the computed, approximately degenerate, well defined eigenstates, the nascent representation, fall on a continuum of locations near the apex of the three state cone. Because of the conical topography these points give rise to distinct values of **g**^{I,J}. From Eq. (9c) we see that the distinct values of **g**^{I,J} correspond to distinct mixings of the N states at the conical intersection. Since at neighboring points on the seam completely different positions on the cones are possible, indeed likely, the nascent states representations behave irregularly along the seam. Our orthogonalization returns all mixtures to a common representation.

The observed continuity is not achieved by simply solving Eq. (10b) as described in Sec. II. This is a consequence of the multiple local solutions to that equation. Instead, a preliminary transformation is first carried out to minimize the quantity

where |${\bf \tilde h}^{I,J} ({\bf Q}_{ref})$| $ h \u0303 I , J ( Q r e f ) $ are the rotated vectors at a reference geometry where an orthogonalized representation is already achieved through the aforementioned rotation procedure and **h**^{U, I, J} is **h**^{I,J} constructed from **U** transformed eigenvectors as in Eq. (9b). Every iteration in the determination of **U**, the phase of the **U** transformed wave functions may be adjusted to accommodate the need to alter the wave function phase as noted at the end of Sec. II B. Δ has a unique minimum which is determined in a very straightforward manner with numerical gradients and line search techniques. After this preliminary transformation convergence of Eqs. (10a) and (10b) is readily achieved. This preconditioning by ensuring that the same minimum of *S* is obtained results in a representation that is continuous along the whole seam. In this work only a single reference geometry was used. More sophisticated selection schemes are possible. In future work it will be interesting to study the topography, particularly the local minima, of *S.*

### D. Implications

The nascent and orthogonalized representations of the g-h space of a point of conical intersection of three states contain the same information about the local topography. However the orthogonalized form has several uses for which the nascent form is not suitable. Here we have focused on the geometry dependence of the representation. Only the orthogonalized representation offers the possibility of extrapolating or interpolating that information to infer changes in the local topography associated with motion in the seam space. In a qualitative sense the orthogonalized representation accomplishes this by aligning the coordinate axes at adjacent geometries evincing only the true geometric variations in the branching space vectors. The rotated vectors were also shown to evince the relevant point group symmetry.

The orthogonalized representation finds application in our construction of quasi diabatic representations of coupled adiabatic potential energy surfaces.^{28} In this case since we fit both energy gradients and derivative couplings we need a way to compare the conical parameters of the fit and *ab initio* representations of a single point of three state conical intersection. The orthogonalized representation enables that comparison. For more details on the use of conical intersections in coupled surface fitting, see Ref. 28. Further, as in the case of two state conical intersections, the orthogonalized g-h space vectors, by avoiding spurious changes due to the incomplete definition which afflicts the nascent representation, can be used to extrapolate the iterative procedure for locating conical intersections. See Ref. 23 for details in the two state case. This issue will be considered in detail in a future study.

## IV. SUMMARY AND CONCLUSIONS

At a point of conical intersection of *N* states, the electronic states are defined up to the *N(N* − 1*)*/2 parameters defining an orthogonal mixing of the degenerate states. This flexibility renders the branching space vectors discontinuous functions of nuclear coordinates along the seam. In this work the orthogonalization approach used remove this flexibility in the two state case is generalized to *N* state conical intersections. Using the three state seam of conical intersection in pyrazolyl as an example, the orthogonalization condition, Eq. (6), or its equivalent Eq. (10b), was shown to efficiently eliminate the discontinuous nuclear coordinate dependence of the branching space vectors in the seam space, that occurs for the nascent eigenstates. The possible issue of local solutions to the defining equations, which also exists in the two state case, where the **g** and **h** vectors can switch, was noted. A method to avoid the issue was demonstrated. Other uses of the rotated representation, including improving the algorithm for locating points of conical intersection, were briefly noted.

## ACKNOWLEDGMENTS

This work was supported by National Science Foundation Grant No. CHE-1361121 to D.R.Y.

## REFERENCES

*ab initio*Electronic Structure Program,