Finite-Field Calculations of Transition Properties by the Fock Space Relativistic Coupled Cluster Method: Transitions between Different Fock Space Sectors

: Reliable information on transition matrix elements of various property operators between molecular electronic states is of crucial importance for predicting spectroscopic, electric, magnetic and radiative properties of molecules. The ﬁnite-ﬁeld technique is a simple and rather accurate tool for evaluating transition matrix elements of ﬁrst-order properties in the frames of the Fock space relativistic coupled cluster approach. We formulate and discuss the extension of this technique to the case of transitions between the electronic states associated with different sectors of the Fock space. Pilot applications to the evaluation of transition dipole moments between the closed-shell-like states (vacuum sector) and those dominated by single excitations of the Fermi vacuum (the 1 h 1 p sector) in heavy atoms (Xe and Hg) and simple molecules of heavy element compounds (I 2 and TlF) are reported.


Introduction
In recent years, there has been a marked increase of interest in highly accurate theoretical data on excited electronic states and electronic transitions in molecules of heavy element compounds. Such data are of key importance for ultra-low temperature physics [1][2][3][4], searches for violation of time-reversal and spatial-parity symmetries of fundamental interactions in low-energy spectroscopic experiments [5] and high-resolution spectroscopy of short-lived radioactive atoms and molecules [6]. The preparation of experimental studies and interpretation of their results require detailed knowledge of energetic, electric, magnetic, radiative and other molecular properties. Besides the information on pure electronic state properties, transition property values, i.e., matrix elements of property operators between the electronic state wavefunctions, is necessary to predict the excited-state dynamics.
The Fock space version of the relativistic coupled cluster method (FS RCC) is one of the most powerful tools of electronic structure modelling for heavy element compounds, widely used for obtaining accurate and reliable information on excited state potential energy surfaces and transition energies in molecules [7][8][9][10]. The prediction of expectation and transition values of various property operators within the FS RCC framework is usually a more complicated task than energy spectra calculations. Due to the exponential form of the wave operators, the expressions for property matrix elements between the FS CC wavefunctions in terms of cluster amplitudes and model eigenvectors have the form of (quasi-) nonterminating series. Truncating these series, one arrives at computational schemes which had been successfully applied to ab initio calculations of transition probabilities in small systems (see, e.g., [11][12][13] and references therein). The simplest scheme of this kind includes abandoning of all the amplitude-dependent terms, and estimating the required property matrix elements between the FS (R)CC wavefunctions by the transition amplitudes between the corresponding model space functions (e.g., left and right eigenvectors of the effective Hamiltonian) [14,15]. This scheme describes properly indirect influence of the couplings between the model space and its orthogonal complement (usually referred to as outer space manifold) on the transition property values. The latter circumstance is of a critical importance for the qualitatively correct description of transition properties in presence of weakly avoided crossings. At the same time, the accuracy of this approach remains limited due to the complete neglect of outer-space contributions to effective property operators.
Apparently the most universal and rigorous approach to evaluating transition matrix elements of first-order properties [16][17][18][19] is based on the use of the constrained-variational technique of Jørgensen et al. [20]. This approach implies solving the so-called Λ equations for each pair of states generated in a single FS CC calculation and leads to rather involved computations; to our knowledge, until now it has only limited implementation within the non-relativistic FS CC framework (for a recent review, see [21]).
A practical solution of the problem might consist in combining the potential surfaces, constructed using multireference coupled cluster methods, with the transition property surfaces obtained within the multireference configuration interaction (MRCI) schemes (see, e.g., [22]). Although the accuracy of MRCI estimates of some transition properties (in particular, transition dipoles) is generally sufficient for practical purposes, mismatches of avoided crossings and/or conical intersection's positions in MRCI and coupled cluster calculations can give rise to serious distortion of corresponding matrix elements between vibronic states. Another serious problem can arise from typically large number of electrons and shells which should be explicitly included into the correlation treatment for heavy element compounds, since the MRCI method is not size-consistent.
The expectation values of numerous important properties can be straightforwardly determined via using the finite-field technique, i.e., by numerical differentiation of the calculated energies with respect to the amplitude of appropriate perturbation (see [23][24][25][26][27] for recent applications). The extension of this approach to evaluating the transition matrix elements of first-order property operators proposed in [28,29] is based on the analysis of variations of the effective Hamiltonian eigenvectors induced by appropriate external field. Although no explicit information on the composition of the wavefunctions outside of the model space is used, the resulting transition matrix element estimates implicitly include the leading contributions from couplings between the model space and the outer space manifolds to the effective property operator [28] (see also Appendix A), thus generally having a significantly better accuracy than those obtained as property matrix elements between effective Hamiltonian eigenvectors. Note, however, that, in contrast to the finite-field coupled-cluster method for expectation values, its counterpart for off-diagonal matrix elements does not ensure the convergency of the results to exact values with appropriate extension of FS (R)CC computational scheme (inclusion of higher excitations in the cluster operator, approaching the complete basis set limit, etc.). The transition matrix elements are determined simultaneously for all pairs of states obtained in FS (R)CC effective Hamiltonian calculations. Applications of the finite-field technique to atomic and molecular relativistic calculations of E1 transition probabilities and off-diagonal magnetic hyperfine interactions were described by the authors of [4,[29][30][31] and the authors of [32,33], respectively.
Up to now, the applications of the finite-field technique were restricted to the cases of single Fock space sector and complete model spaces. If the model spaces are incomplete and the model states belong to the same Hilbert space but to different Fock space sectors, the straightforward use of the finite-field procedure in its present form [29] would always lead to unphysical zero values of the corresponding off-diagonal matrix elements. In the present paper, we describe the extension of the technique to the complicated cases mentioned above, focusing on the calculations of transition matrix elements between the closed-shell-like ground state and excited states dominated by single excitations (0h0p-1h1p-type transitions). Pilot applications to transition dipole moment calculations in heavy atoms and diatomic molecules of heavy element compounds are reported.

Theory
Consider a first-order property operator D defining the dependence of the total many-electron Hamiltonian H on some external field strength parameter F, H(F) = H(0) + DF. Let us choose a field-independent model space L P as the linear span of several Slater determinants and denote the projectors onto L P and its orthogonal complement by P and Q, respectively. Within the FS RCC approach, the field-dependent many-electron problem is solved through constructing the wave operator Ω and the effective Hamiltonian H. The diagonalizaition of H operating within L P , yields low-lying (target) eigenvalues E i , i = 1, . . . , dim L P of the total Hamiltonian H, whereas the wave operator Ω should reconstruct the target eigenfunctions ψ i of the total Hamiltonian from the corresponding eigenvectors of H (model eigenstates), ψ i = Ω ψ i . The wave operator is assumed to be a single normal-ordered exponential of a cluster operator T (linear combination of excitation operators, usually classified according to their total ranks k and the numbers of annihilated holes n h and particles n p ): Naturally, both H and Ω depend on the parameter F. We suppose that { ψ i } are normalized, so that the norms N i : N 2 i = ψ i |ψ i of the full Hamiltonian eigenvectors exceed 1.
It is convenient to express the transition property matrix elements in terms of the eigenvectors ofH and the effective property operator D. We adopt the following definition of D [34]: where Ω is the "left" counterpart of the wave operator Ω, (P denotes the orthogonal projector onto the subspace spanned by the target eigenfunctions of H and the inversion is restricted to the model space). The matrix elements (4) are related to those of the effective operator (5) by [34]: Here, { ψ ⊥⊥ i } are the left eigenfunctions of H, biorthonormalized to its right eigenfunctions, The derivative of the effective Hamiltonian with respect to the parameter F does not coincide with the effective property operator [28]: Rewriting Equation (8) in terms of matrix elements as one can notice that the additional term vanishes for the diagonal matrix elements but survives for the off-diagonal ones. However, a simple perturbation theory analysis [28] demonstrates that the approximation of D by the effective Hamiltonian derivative incorporates the leading contributions from the couplings to the orthogonal complement of L P . Combined with the off-diagonal Hellmann-Feynman theorem [35] for the effective Hamiltonian, this approximation yields The finite-field technique is based on estimating the derivative matrix element in the rhs of (10) within the central finite-difference approximation, evaluating and diagonalizing the effective Hamiltonians H(−∆F) and H(∆F) for two small finite values of the field strength parameter: Absolute values of transition matrix elements can be immediately obtained from Equations (7) and (11): for all pairs of target states simultaneously. It is worth underlining that the choice of the wave operator normalization is essential for the reliability of the basic approximation (10). If the condition (3) is fulfilled for any value of field strength parameter, (∂T/∂F) has no closed component. Taking into account the exponential form of Ω and Equation (6), one can represent the operator determining the error of our basic approximation (see Equation (9)) by the series and notice that in the case of intermediate normalization the first term (which can be large even in the case when all amplitudes for F = 0 are small) vanishes.
The situation with incomplete model spaces and transitions between the states associated with different Fock space sectors seems to be more complicated. For the sake of simplicity, let us focus on the case of transitions from the closed-shell like state to those dominated by the creation of single hole-particle pairs and described within the valence-universal FS (R)CC approach for the 1h1p sector; the generalization to other cases of "intersector" transitions is straightforward. The model space should include the vacuum closed-shell state and the subspace spanned by all possible combinations of one active hole and one active particle; the corresponding projectors are denoted by P (0h0p) and P (1h1p) , The conventional valence-universal FS (R)CC formulation supposes that the 1h1p sector [7,36] is fully decoupled from the vacuum sector, P HP = P (0h0p) HP (0h0p) + P (1h1p) HP (1h1p) . (15) so that the model wavefunctions ψ i (F) and ψ j (F ) for the ground state i and an excited state j are strictly orthogonal for any field strengths F and F and the finite-field estimate of transition property matrix elements within L P seems to be senseless. Note that the intermediate normalization of Ω is incompatible with the connectivity of the effective Hamiltonian for incomplete model spaces [37,38]. Two types of cluster components contributing to PΩP are usually introduced: single excitations creating an active hole and an active particle, T (see Figure 1) explicitly contains the terms linearly and quadratically dependent on cluster amplitudes. Let us perform a posteriori transformation of the effective Hamiltonian (15), corresponding to the restoration of the intermediate normalization of the wave operator, (cf. [39]). Taking into account Equation (16) and the normalized exponential form of the original wave operator, one readily realizes that where T Op is the purely open (PT Op P = 0) part of the cluster operator and O(T 2 ) stands for the sum of terms of second and higher orders in cluster amplitudes. In a strict analogy to Equation (13), the first term in the similar expansion of P Ω (∂Ω /∂F)P should disappear. One can thus expect that the model-space finite-field transition property estimates obtained according to Equations (11) and (12) with the eigenfunctions of the transformed Hamiltonian H are generally as accurate as in the case of transitions within a single sector. The computational efforts for transforming the effective Hamiltonian (or its eigenvectors) are negligible since the number of amplitudes required for constructing PΩP is small. Although in the case of the scalar property operator the CC equations are to be solved at least twice, for F = −∆F and F = ∆F, the choice of a common (field-independent) set of one-electron spinors enables one to avoid the doubling of the amount of computations at least for one-electron properties, using the same transformed and sorted two-body integrals and other intermediate data in both calculations. Furthermore, due to the smallness of the step size ∆F, the amplitudes obtained with one field strength value provide an excellent initial guess for computing the amplitudes for another field strength, reducing dramatically the number of iterations. Unfortunately, the substantial increase of computational work is unavoidable when the applied field lowers the symmetry of the system under study.
Let us notice finally that the simplest finite-field scheme does not make use of the information concerning the dependence of the cluster amplitudes on F which can be readily extracted from the results of calculations. Obviously, one could employ this information to estimate the derivative ∂(ΩP)/∂F in the rejected term in the rhs of Equation (9) and to evaluate approximately this term, truncating in some way the expansions in powers of the cluster operator.

Pilot Applications to Transition Dipole Moment Calculations
In this section, we apply the finite-field FS RCC scheme to evaluate transition dipole moments in heavy atoms (Xe, Hg) and molecules containing heavy atoms (I 2 , TlF). The calculations described below employed the accurate relativistic shape-consistent semilocal pseudopotential model [40][41][42] for heavy atoms; the inner core shells with principal quantum numbers 1-3 for fifth-row atoms and 1-4 for sixth-row atoms were excluded from the explicit treatment. The light fluorine atom was treated in the nonrelativistic all-electron manner. The cluster operator expansions were always restricted to single and double excitations (FS RCCSD). One-electron spinors were generated by solving the twoor four-component spin-orbit-coupled SCF equations. The component of spinors were expanded in contracted Gaussian basis sets. All explicitly treated electrons were correlated. The valence-universal approach to constructing the effective Hamiltonian for the 1h1p sector for the neutral systems with zero or small electron affinities gives rise to certain difficulties in choosing the set of active "particle" spinors. Normally one faces a rather objectionable choice between problematically large amplitudes of T (0h1p) 1 if the active particle spinors are compact (e.g., taken from the solutions of the SCF problem for a positive ion) and large T (1h1p) 2 amplitudes when these spinors are diffuse (obtained as virtual SCF spinors for the neutral system). We bypassed this difficulty via choosing quite large sets of active particle spinors and using the adjustable denominator shift technique [43] to suppress the possible divergencies in the presence of intruder states (more specifically, we employed the simulated imaginary shifts [32]). To conserve the core separability of the original FS RCC scheme, shifts were never applied to the energy denominators in the equations for T (0h0p) amplitudes; a balanced treatment of excitation and deexcitation contributions to PΩP implied the use of non-shifted equations for T (1h1p) 1 amplitudes as well. Shift amplitudes (s K in Equation (8) of [32]) were always defined by a single parameter s 2 , universal denominator shift amplitude for all double excitations in the 1h0p, 0h1p and 1h1p sectors; for single excitations in the 0h1p and 1h0p sectors, we used s K = s 2 /2. The shift attenuation parameter (m in Equation (8) of [32]) was assumed equal to 3.
The construction of spinors and evaluation of one-and two-electron integrals was performed using the DIRAC19 program package [15,44], whereas all FS RCC calculations were carried out with the help of the EXP-T code [45][46][47]. The VIBROT code [48] was used to solve the vibrational problem and compute the excited-state lifetime of TlF.

Transition Dipoles for Excitations of Closed-Shell Atoms
Experimental-based estimates of transition dipole moments in heavy-element-containing species accurate to about 1% and better are generally available only for atoms and atomic ions. Reliable estimates of transition dipole moments for several transitions in atomic Xe and Hg (Tables 1 and 2) were obtained from Einstein coefficients from Sansonetti and Martin [49]. Our FS RCC calculations for Xe employed the adaptation of the aug-cc-pVQZ-PP basis of Peterson et al. [50,51] to the pseudopotential [42], augmented by the p 3/2 − p 1/2 -like functions to improve the description of spin-orbit splittings, core-valence correlation functions from aug-cc-pwCVTZ-PP set [52] and additional diffuse functions in order to better reproduce the Rydberg nature of the excited states. The total basis size was [12s13p11d6 f 3g] (see Supplementary Materials for details). The one-electron spinors were obtained by solving the SCF equation for the ground-state configuration of the neutral atom (5p 6 ) or positive ion (5p 5 ); the corresponding results are marked in Table 1 by the symbols (Xe 0 ) or (Xe + ). Three pairs of highest-energy spinors occupied in the neutral Xe and 43 pairs of lowest-energy unoccupied spinors were considered as active. Rather small shift amplitudes (s 2 = −0.2 a.u.) were found to be sufficient to ensure the numerical stability of solving the FS RCC equations.
The resulting finite-field transition dipole moment estimates D FF ij along with the experimental-based values and matrix elements of the dipole operator between the FS RCC model , are listed in Table 1. The deviation of D FF ij from the empirical values never exceeds 0.04 a.u.; for intensive transitions, the error is within 4 %. It is worth noting that D FF ij only slightly depend on the choice of one-electron spinors (which, in turn, defines the model space). The dipole matrix elements between the H eigenvectors, despite their strong dependence on this choice, can be considered as less accurate but rather reasonable preliminary estimates of the transition dipole moments. Similar calculations were carried out for the Hg atom. The employed primitive Gaussian basis set [11s11p10d5 f 4g3h] was essentially an uncontracted version of the augmented quadruple zeta basis [53] adapted to the use with the shape-consistent pseudopotential [42]. We performed two series of calculations with one-electron spinors taken from the neutral atom (configuration 5d 10 6s 2 ) or positive ion (5d 10 6s 1 ) SCF problem. Six Kramers pairs of highest occupied spinors and 31 pairs of lowest unoccupied spinors spanned active-hole and active-particle subspaces, respectively. The same shift amplitudes as in the case of the xenon atom were used. The results compared to the experimental values are listed in Table 2. The finite-field scheme underestimates the 1 S 0 ↔ 3 P o 1 transition moment by 16% and the 1 S 0 ↔ 1 P o 1 transition moment by 11%. This level of accuracy is actually close to that reported in [29] for the finite-field calculations within the (0h2p) sector. It might be worth noting that the difference between the finite-field estimates of ψ ⊥⊥ i | D| ψ j and ψ ⊥⊥ j | D| ψ i according to Equation (11) were quite large (up to 20 %). This can be presumably attributed to significant differences in the norms of the ground-state and excited-state wavefunctions outside of the model space (see Equation (7)), in other words, to a non-balanced description of these states within the chosen model space. The results are stable with respect to further extensions of active spinor subspaces. A possible way of improving the D FF ij estimates in such cases might consist in accounting for the dependence of the wave operator on the external field strength Equation (9). Similarly to the case of Xe, the D FF ij values for Hg exhibit only a very weak dependence on the choice of the SCF problem for generating one-electron spinors.

Transition Dipole Moment Functions in I 2 and TlF
Optical transitions in the iodine molecule are among the most studied molecular electronic transitions in heavy element compounds (see [54] for a recent review). The rather intensive X0 + g −B0 + u transition is of particular interest for testing relativistic transition dipole calculation techniques since it corresponds to the forbidden singlet-triplet ( 1 Σ + g − 3 Π u ) excitation at the non-relativistic (or scalar relativistic) limit. The empirical dependence of the X0 + g − B0 + u transition dipole moment on the internuclear separation derived from a large amount of experimental data by Tellinghuisen [55] is expected to be accurate within a few percent at least in the vicinity of the ground-state equilibrium point. Our calculations were performed using the [9s 10p 8d 5 f 3g] basis set for iodine constructed in the same way as for Xe, except for the absence of Rydberg-like functions which are not needed to describe the valence-like states under study. The active space comprised five Kramers pairs of occupied spinors and 17 pairs of virtual spinors. The shift amplitude values corresponded to s 2 = −0.30 a.u.
The finite-field FS RCC transition dipole moment curve ( Figure 2) exhibits a good agreement with its empirical counterpart for moderate internuclear separations where the single-reference treatment of the ground state remains reasonable, i.e., approximately for R ≤ 3 Å. At larger R, an unusually rapid increase of amplitudes of double excitations in the vacuum sector (first of all, σ 2 g → σ * 2 u -like ones), see the dotted curve in Figure 2, gradually lowers the efficiency of the single-reference description of the ground state. The deviation of D FF X−B values from their empirical counterparts in the vicinity of the ground-state equilibrium point is about 4%. Except for the domain of large R, our results also agree with those of quasirelativistic many-body multipartitioning perturbation theory calculations [56] which ensured equal-footing and inherently multireference treatment of both states involved. In contrast with the atomic cases, the quality of transition dipole estimates by the dipole matrix elements between the model space vectors seems fully unacceptable for benchmark calculations.
The transition between the ground X(1)0 + and second excited B(1)1 (sometimes also referred to as a 3 Π 1 ) states of the thallium monofluoride molecule attracts great attention due to its possible use for laser cooling of TlF, which is considered as a perspective system for searches for "new physics" in tabletop-scale experiments (see [2] and references therein). The lifetime of the lowest vibrational level (v = 0) of the B(1)1 state, 99 ± 9 ns, was measured by Hunter et al. [2]. It should be noted that the accuracy of the lifetime measurements in [2] is of an unprecedented level for molecules of heavy element compounds. Due to nearly coinciding equilibrium separations R e and vibrational frequencies of this state with those of the ground one and the negligible probability of its radiative decay into the first excited A(2)0 + state (energy factor of the branching rate is about 10 −4 ), one can immediately obtain an estimate of the transition dipole at R e (B) ≈ R e (X): D X−B = 0.315 ± 0.014 a.u. The computational study of radiative transitions in TlF was performed recently by Liu et al. [57] using the scalar relativistic multireference configuration interaction method with subsequent incorporation of the direct spin-orbit interactions between low-lying scalar eigenstates (MRCI + SO). The resulting B(1)1, v = 0 lifetime (85 ns) and the transition dipole moment value at R e (X) ≈ R e (B) (0.321 a.u.) agree perfectly with their experimental counterparts. The present calculation used the [9s10p8d6 f 4g1h] Tl basis set with core-valence-correlation and diffuse functions from Zaitsevskii and Eliav [43] and the standard aug-cc-pVQZ fluorine basis [58,59], the solutions of the ground-state SCF problem as one-electron spinors, the model space defined by 4 pairs of active hole spinors and 16 pairs of active particle spinors, and the same shift parameters as for the I 2 molecule. Our estimate of vertical transition moment (0.297 a.u.) is quite close to both experimental-based and MRCI + SO values. In accordance with the MRCI + SO data [57], D FF X−B as function of the internuclear separation is rather flat near the equilibrium points and becomes decreasing at larger R; the difference between D FF X−B and its MRCI + SO counterpart becomes more significant at short R (see Figure 3). Combining the present transition dipole moment and energy difference function along with the accurate RCCSD(T) ground state potential from Zaitsevskii and Eliav [43], we evaluated the radiative lifetime of B1, v = 0 according to the Tellinghuisen's formula [60]. The resulting value, 112 ns, is only slightly longer than the measured one.
Similar to the case of transition in I 2 , the dipole matrix elements between the model states provide much worse estimates of the transition moment (however, the correct trend is reproduced).

Concluding Remarks
The simple finite-field technique of molecular transition property matrix element calculations within the framework of the Fock space relativistic coupled cluster method, deriving the matrix element estimates from the analysis of variations of the model-space projections of electronic eigenstates under the influence of an external field, is extended to the case of transitions between the states associated with different Fock space sectors. In contrast to the case of transitions between the states of the same sector, these projections are not immediately obtained as eigenstates of FS RCC effective Hamiltonian; an additional transformation of these eigenstates is required. The transformation matrix is readily constructed from the cluster amplitudes corresponding to the excitations which are closed with respect to the direct sum of the involved model subspaces of the Fock space. Pilot applications to off-diagonal electric dipole matrix elements for excitations of closed-shell-like ground states in systems containing heavy-element atoms yielded the estimates approximately at the same level of accuracy as for "single-sector" transitions, normally within 10 % for intensive transitions. This level of accuracy is comparable with the experimental one for the molecular measurements. The reliability of results is clearly superior to that achieved in calculation of dipole matrix elements between the eigenvectors of FS RCC effective Hamiltonians.
Due to the use of field-independent model spaces and one-electron spinors for solving the many-electron problem for a molecule in the finite-strength external field, the additional work required to evaluate transition property values remains moderate. However, in applications to highly symmetric systems and in the cases when one focuses on a small number of transitions but has to estimate off-diagonal matrix elements of numerous or multicomponent (tensor) properties the technique can become computationally extensive, so that the development of relativistic counterparts of analytical methods for transition density matrix construction remains desirable.
An important advantage of the finite-field technique consists in the possibility to simultaneously evaluate the transition matrix elements for all pairs of states corresponding to eigenvalues of the FS RCC effective Hamiltonian. This advantage can be decisive for studying intensity distributions in dense spectra, especially in those of actinide and lanthanide compounds. The work on assessment of reliability of the technique as a tool for describing the excitations in molecules and cluster models of crystals containing f -elements [61] is in progress in our group. This expression does not coincide with that for the second-order one for the field-free effective property operator: D (2) = PDΩ (2) P + P Ω (2) DP + P Ω (1) DΩ (1) P = PDΩ (2) P − PVG 2 VPDP + PVGVGDP − PVPVG 2 DP + PVGDGVP, but the difference is restricted to only one class of renormalization-type terms (underlined).