text:
# GRECP/MRD-CI calculations of the spin-orbit splitting in the ground state of Tl and of the spectroscopic properties of TlH.

## I Introduction

During the last few years a large number of publications have dealt with calculations of the $`{}_{}{}^{2}P_{1/2}^{o}{}_{}{}^{2}P_{3/2}^{o}`$ splitting in the ground state of the Tl atom and spectroscopic constants for the $`0^+`$ ground state of TlH. Such interest to these systems arises because of their relatively simple electronic structure in the valence region. This makes them very convenient objects for testing methods for the description of relativistic and correlation effects. We can mention some recent papers  in which the electronic structure of thallium was studied and papers  in which the calculation of spectroscopic constants for TlH was carried out. With the exception of the atomic RCC calculation by Eliav et al.  and the atomic CI/MBPT2 calculation by Dzuba et al. , the published results cannot be considered to be very accurate and reliable, however, primarily because of the rather small basis sets and the small numbers of correlated electrons.

In calculations of Tl and TlH with the use of the relativistic effective core potential (RECP) approximation , in which only 13 thallium electrons are treated explicitly (13e-RECPs), one more problem appears. The correlation of the outer core (OC) and valence (V) electrons, occupying the $`5d`$ and $`ns,np,nd`$ ($`n=6,7,\mathrm{}`$) orbitals, respectively, cannot be satisfactorily described, mainly because the smoothed V-pseudoorbitals (pseudospinors) have the wrong behaviour in the OC region. One-electron functions $`\varphi _{x,k}^{corr}(r)`$, being some linear combinations of virtual orbitals, correlate to occupied orbitals $`\varphi _x^{occ}`$ (where $`x=c,v`$ stands for the OC and V orbital indices) and are usually localized in the same space region as $`\varphi _x^{occ}`$. Therefore, the original “direct” Coulomb two-electron integrals describing the OC-V correlation of $`\varphi _c^{occ}`$ and $`\varphi _v^{occ}`$ can be well reproduced by those with the pseudoorbitals, despite their localization in different space regions. However, a two-electron integral describing the “exchange” part of the OC-V correlation,

$$_𝐫𝑑𝐫\varphi _{c,k^{}}^{corr^{}}(𝐫)\varphi _v^{occ}(𝐫)_𝐫^{}𝑑𝐫^{}\varphi _{v,k}^{corr^{}}(𝐫^{})\varphi _c^{occ}(𝐫^{})\frac{1}{𝐫𝐫^{}},$$
(1)
cannot be well reproduced because the V-pseudoorbitals are smoothed in the OC region where the OC-pseudoorbitals are localized (for more theoretical details, see Ref. ).

The first RECPs for Tl with the $`5s,5p`$ shells treated explicitly (21e-RECPs) for which this disadvantage of the earlier “semicore” RECPs was overcome were generated and tested in single-configurational calculations by Mosyagin et al. . Some other inherent problems of the “nodeless” RECPs were also solved with the 21-electron Generalized RECP (21e-GRECP) version presented in Ref. . In Ref. , for the case of the 21e-GRECP it was also shown that the $`5s,5p`$ pseudospinors could be frozen while still providing significantly higher accuracy than 13e-RECPs because the valence and virtual $`ns`$ and $`np`$ ($`n=6,7,\mathrm{}`$) pseudoorbitals in the former case already have the proper nodal structure in the OC region.

## II The GRECP operator in the spin-orbit representation

In most existing quantum-chemical codes for molecular calculations with RECPs (as well as in the MRD-CI code used in the present work) spin-orbit basis sets are used. In these versions the number of the two-electron integrals is substantially smaller than in the case of spinor basis sets providing the same level of correlation treatment. Spin-orbit basis sets are preferable in the calculations in which correlation effects give a higher contribution to the properties of interest than those of a relativistic nature. This is usually the case for valence and outermost core electrons, which mainly determine chemical and spectroscopic properties of molecules.

Together with the spin-orbit basis set, the GRECP for Tl should also be employed in the spin-orbit representation. Following Ref. , the components of the spin-averaged part of the GRECP operator called the averaged relativistic effective potentials (AREP) are written in the form :

$$U_{n_vl}^{AREP}(r)=\frac{l+1}{2l+1}U_{n_vl+}(r)+\frac{l}{2l+1}U_{n_vl}(r),$$
(2)
$$𝐔_{n_cl}^{AREP}(r)=\frac{l+1}{2l+1}𝐕_{n_cn_vl+}(r)+\frac{l}{2l+1}𝐕_{n_cn_vl}(r),$$
(3)
$`𝐕_{n_cn_vl\pm }(r)`$ $`=`$ $`\left[U_{n_cl\pm }(r)U_{n_vl\pm }(r)\right]\stackrel{~}{𝐏}_{n_cl\pm }(r)+\stackrel{~}{𝐏}_{n_cl\pm }(r)\left[U_{n_cl\pm }(r)U_{n_vl\pm }(r)\right]`$ (4)
$``$ $`{\displaystyle \underset{n_c^{}}{}}\stackrel{~}{𝐏}_{n_cl\pm }(r)\left[{\displaystyle \frac{U_{n_cl\pm }(r)+U_{n_c^{}l\pm }(r)}{2}}U_{n_vl\pm }(r)\right]\stackrel{~}{𝐏}_{n_c^{}l\pm }(r),`$ (5)
where $`U_{nl\pm }(r)`$ are the potentials generated for the $`\stackrel{~}{\phi }_{nl\pm }(r)`$ pseudospinors by means of the Goddard scheme ; $`n_x`$ is the principal quantum number of an outercore ($`n_c`$), valence ($`n_v`$) or virtual ($`n_a`$) pseudospinor; $`l`$ and $`j`$ are angular and total electron momenta; $`\pm `$ stands for $`j=l\pm 1/2`$; $`\stackrel{~}{𝐏}_{n_cl\pm }(r)`$ is the radial projector on the OC pseudospinors:

$$\stackrel{~}{𝐏}_{n_cl\pm }(r)=\underset{m}{}\stackrel{~}{|n_c,l,\pm ,m}\stackrel{~}{n_c,l,\pm ,m|}.$$
(6)
Clearly, the AREP component of the GRECP may be used in calculations with nonrelativistic quantum-chemical codes in order to take account of spin-independent relativistic effects.

The operator of the effective spin-orbit interaction can be derived following the expression for the spin-angular projector $`𝐏_{l\pm }`$ from Ref. :

$$𝐏_{l\pm }(\mathrm{\Omega },\sigma )=\frac{1}{2l+1}\left[\left(l+\frac{1}{2}\pm \frac{1}{2}\right)𝐏_l(\mathrm{\Omega })\pm 2𝐏_l(\mathrm{\Omega })\stackrel{}{𝐥}\stackrel{}{𝐬}𝐏_l(\mathrm{\Omega })\right].$$
(7)
Its components, called the effective spin-orbit potentials (ESOP), can be written as

$$\mathrm{\Delta }U_{n_vl}(r)=U_{n_vl+}(r)U_{n_vl}(r),$$
(8)
$$\mathrm{\Delta }𝐔_{n_cl}(r)=𝐕_{n_cn_vl+}(r)𝐕_{n_cn_vl}(r),$$
(9)
$$𝐔_{nl}^{ESOP}=\frac{2\mathrm{\Delta }𝐔_{nl}(r)}{2l+1}𝐏_l\stackrel{}{𝐥}\stackrel{}{𝐬},$$
(10)
$$𝐏_l=\underset{m=l}{\overset{l}{}}|lmlm|,$$
(11)
where $`|lmlm|`$ is the projector on the spherical function $`Y_{lm}`$.

Neglecting the difference between $`U_{n_vL}^{AREP}`$ and $`U_{n_vLJ}`$ for virtual pseudospinors with $`l>L`$ (for theoretical details see Ref. ), one can write the GRECP operator $`𝐔`$ as

$`𝐔`$ $`=`$ $`U_{n_vL}^{AREP}(r)+{\displaystyle \underset{l=0}{\overset{L1}{}}}\left[U_{n_vl}^{AREP}(r)U_{n_vL}^{AREP}(r)\right]𝐏_l`$ (12)
$`+`$ $`{\displaystyle \underset{n_c}{}}{\displaystyle \underset{l=0}{\overset{L}{}}}𝐔_{n_cl}^{AREP}(r)𝐏_l+{\displaystyle \underset{l=1}{\overset{L}{}}}\left[𝐔_{n_vl}^{ESOP}+{\displaystyle \underset{n_c}{}}𝐔_{n_cl}^{ESOP}\right]𝐏_l.`$ (13)
Note that the nonlocal terms with the projectors on the most important correlation functions $`\stackrel{~}{\phi }_{n_xl\pm ;n_kl_k\pm }^{corr}(r)`$ (1) (where $`x=c,v`$) localized mainly in the OC and V regions and with the corresponding potentials $`U_{n_kl_k\pm }^{corr}(r)`$ can be taken into account in the considered expressions for the GRECP operator additionally to those with the OC projectors. Obviously, the non-local GRECP terms for the frozen OC pseudospinors can be omitted in the sum over ($`n_cl`$) in Eq. (13).

We should emphasize that in spite of the rather complicated form of the above GRECP operator, the main computational effort in calculating matrix elements with the GRECP is caused by the standard radially-local operator, which is also a part of conventional RECP operators, and not by the non-local GRECP terms. Thus, the additional complications in calculations with GRECPs are negligible in comparison with treatments employing conventional semi-local RECPs if comparable gaussian expansions are used for the partial potentials. The more critical point is that the effort in the calculation and transformation of two-electron integrals is always substantially higher than that in the computation of RECP integrals for all known RECP versions (including GRECPs) when appropriately large basis sets are employed in the precise calculations.

## III Frozen-core approximation for the outer-core shells

To perform precise calculations of chemical and spectroscopic properties, correlations should be taken into account not only within the valence regions of heavy atoms and heavy-atom molecules but in the core regions and between the valence and core electrons as well. In practice, the goal is to achieve a given level of accuracy by correlating as small a number of electrons as possible, thus reducing the computational effort. However, as discussed in the Introduction, the accuracy of the RECPs generated for a given number of explicitly treated electrons cannot always satisfy the accuracy requirements expected from correlating all these electrons in the corresponding all-electron calculation. This is true, in particular, for calculations of Tl, having a $`5d^{10}6s^26p^1`$ leading configuration in the ground state, and its compounds.

To attain an accuracy level of 400 cm-1 for the $`{}_{}{}^{2}P_{1/2}^{o}{}_{}{}^{2}P_{3/2}^{o}`$ splitting in the ground state and for excitation energies to low-lying states of Tl and to take account of the core polarization, one should correlate at least 13 electrons, i.e. include the $`5d`$ shell. This is achieved in the present MRD–CI calculations with $`f`$ and $`g`$ basis functions describing mainly polarization of the $`5d`$ shell (for other recent results see, e.g., ). Some data from our 13e-CI calculations of the SO-splitting in the ground state of Tl are collected in Table I in comparison with the 3e-CI results, which in our DF/CI (Dirac-Fock calculations followed by CI) and the GRECP/CI calculations have errors of about 600 cm-1.

We also should mention the recent relativistic coupled-cluster (RCC) results of Landau et al. , in which 35 electrons are correlated and a decrease of close to 90 cm-1 in the above mentioned SO splitting is due to the Breit interaction. Note that this interaction is not yet taken into account in the RECPs considered in the present work.

Obviously, the $`5d`$ shell should also be explicitly treated in calculations of molecules containing Tl to take into account core relaxation and polarization effects with satisfactory accuracy. For these calculations it would be optimal to use the RECPs with 13 electrons of Tl treated explicitly (13e-RECPs) such as the RECP of Ross et al.  or our valence RECP version . None of the known nodeless 13e-RECPs can provide the aforementioned accuracy, however. Although single-configurational tests  give errors of 100 cm-1 or somewhat more for excitation energies to low-lying states, they are dramatically increased for 13e-RECPs if all 13 electrons are correlated. The reasons are discussed in the Introduction (one can also see the results of the 13e-RECP/MRD-CI calculations in Ref.  and of the 13e-PP/MRCI calculations in Ref. ).

To overcome this disadvantage, one should use RECPs with at least 21 electrons, e.g. 21e-GRECP  and 21e-PP  for Tl. The $`5s`$ and $`5p`$ pseudospinors can be treated as frozen, however, while still providing the aforementioned accuracy. The $`5p`$ orbitals have energies about four times higher and their average radii are 1.4 times shorter than those for the $`5d`$ orbitals. Moreover, their angular correlation is supressed as compared with the $`5d`$ shell because the most important polarization functions ($`5d`$ for the $`5p`$ orbitals and $`5p`$ for the $`5s`$ orbitals) are completely occupied in the lowest-lying states. Therefore, the $`5s,5p`$ orbitals are substantially less active in chemical processes.

In order to freeze the $`5s`$ and $`5p`$ pseudospinors, one can apply the energy level shift technique . Following Huzinaga et al. , one should add the matrix elements of the SCF field operators (the Coulomb and spin-dependent exchange terms) over these OC pseudospinors to the one-electron part of the Hamiltonian together with the level shift terms

$$\underset{n_{c_f},l,\pm }{}B_{n_{c_f}l\pm }\stackrel{~}{𝐏}_{n_{c_f}l\pm }(r),$$
(14)
where $`B_{n_{c_f}l\pm }`$ is at least of order $`|2\epsilon _{n_{c_f}l\pm }|`$ and $`\epsilon _{n_{c_f}l\pm }`$ is the orbital energy of the OC pseudospinor $`\stackrel{~}{\varphi }_{n_{c_f}l\pm }(r)`$ to be frozen. Such nonlocal terms are needed in order prevent collapse of the molecular orbitals to the frozen states (the $`5s_{1/2},5p_{1/2,3/2}`$ pseudospinors for Tl). All terms with the frozen core pseudospinors described here (the Coulomb and exchange interactions, and the level shift operator) can easily be presented in spin-orbit form with the help of eq. (7), as was done above for the GRECP operator.

More importantly, these OC pseudospinors can be frozen in calculations with spin-orbit basis sets and they can already be frozen at the stage of calculation of the one-electron matrix elements of the Hamiltonian, as implemented in the MOLGEP code . Thus, any integrals with indices of the frozen spinors are completely excluded after the integral calculation step.

In single-configurational calculations with the numerical HFJ code  we have seen that the SO splitting of the $`5p`$ shell increases the resulting SO splitting of the $`{}_{}{}^{2}P_{}^{o}`$ ground state by about 400 cm-1, whereas the SO splitting of the $`5d`$ shell decreases the final SO splitting by almost the same value. Therefore, it is important to freeze the $`5p_{1/2}`$ and $`5p_{3/2}`$ (pseudo)spinors and not some averaged $`5p`$ (pseudo)orbitals if the SO interaction is to be taken into account in the $`5d`$ and valence shells.

In Ref. , the 21e-“energy-adjusted” Pseudopotential (PP) having the features which have been emphasized  as inherent for GRECPs (different potentials for the $`5p`$ and $`6p`$ pseudospinors in the case of Tl) is generated and applied to the calculation of the SO splitting in Tl, with the core correlations described by the core polarization potential (CPP). Some average OC pseudoorbitals are frozen and the SO splitting of 7810 cm-1 obtained in their 21e-PP/MRCI calculation is quite different than our result.

After applying the projection operator of eq. (7) to the level shift (14), Coulomb and exchange terms with the frozen core pseudospinors, the AREP and ESOP parts of the GRECP operator are to be modified to include these new contributions. This technique was successfully employed in our earlier calculations of the spin-rotational Hamiltonian parameters in the BaF and YbF molecules .

The freezing technique discussed above can be efficiently applied to those OC shells for which the spin-orbit interaction is clearly more important than the correlation and relaxation effects. If the latter effects are neglected entirely or taken into account within “correlated” GRECP versions , the corresponding OC pseudospinors can be frozen and the spin-orbit basis sets can be successfully used for other explicitly treated shells. This is true for the $`5p_{1/2,3/2}`$ subshells in Tl, contrary to the case of the $`5d_{3/2,5/2}`$ subshells. Freezing the OC pseudospinors allows one to optimize an atomic basis set only for the orbitals which are varied or correlated in subsequent calculations, thus avoiding the basis set optimization for the frozen states and reducing the number of the calculated and stored two-electron integrals. Otherwise, if the $`5p`$ shell should be correlated explicitly, a spinor basis set can be more appropriate than the spin-orbit one.

## IV The MRD-CI method

In the multireference single- and double-excitation CI approach , the $`\mathrm{\Lambda }S`$-basis sets of many-electron spin-adapted (and space symmetry-adapted) functions (SAFs) are employed. This method makes use of configuration selection and perturbative energy extrapolation techniques  and employs the Table CI algorithm  for efficient handling of the various open-shell cases which arise in the Hamiltonian matrix elements. Some new features of the selection scheme used in this work are considered below. The higher excitations in the CI treatment has been assessed by applying the generalized multireference analogue  of the Davidson correction  to the extrapolated $`T=0`$ energies of each root.

After selecting the $`\mathrm{\Lambda }S`$-sets of SAFs for a chosen threshold $`T_i(i=1,2)`$, they are collected together in accord with the relativistic double-group symmetry requirements and a spin-orbit CI (SO-CI) calculation is performed with these SAFs to obtain some SO-roots ($`\mathrm{\Psi }_I^{SO,T_i}`$) and their energies ($`_I^{SO,T_i}`$) which are of interest in a considered double group irreducible representation (irrep). Then the linear $`T=0`$ correction is evaluated in the basis of the calculations with the $`T_1`$ and $`T_2`$ thresholds. Finally, the generalized Davidson (or full CI) correction is applied to each root of interest.

The stage of the molecular spectroscopic constants calculation begins with the fitting of the relativistic CI potential curves to polynomials which are employed to construct appropriate Born-Oppenheimer nuclear motion Schrödinger equations solved by the Dunham method with the help of the DUNHAM-SPECTR code of Mitin .

### A Features of the spin-orbit selection procedure

Let us define a Hamiltonian H for a molecule as

$$𝐇=𝐇^{(0)}+𝐕^{corr}+𝐇^{SO},$$
(15)
where $`𝐇^{(0)}`$ is an unperturbed spin-independent Hamiltonian, $`𝐕^{corr}`$ is a two-electron operator describing correlations, and $`𝐇^{SO}`$ is a one-electron spin-orbit operator (ESOP in our case). Let us choose an orthonormal basis set of SAFs $`\{\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}\}`$ in the $`\mathrm{\Lambda }S`$-coupling scheme (or “spin-orbit” basis set). In particular, these SAFs can be solutions of Hartree-Fock equations with a spin-averaged RECP for the molecule considered. The $`𝐇^{(0)}`$ Hamiltonian is constructed to be diagonal in the given many-electron basis set:

$$𝐇^{(0)}\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}=E_I^{(n)\mathrm{\Lambda }S}\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S},$$
(16)
where $`n=0,1,\mathrm{}`$ (see below the description of the indices in more detail). Additionally define $`𝐇^{(0)}`$ so that

$$<\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}|𝐇^{(0)}|\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}><\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}|𝐇|\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}>$$
(17)
in order to exclude the first-order PT contributions to total energies of molecular states (this corresponds to the Epstein-Nesbet PT form).

We will ignore the two-electron spin-dependent (Breit) interactions which ordinarily can be neglected when studying chemical and spectroscopic properties. Breit and other quantum electrodynamic (QED) effects are relatively large for lanthanides and actinides, but for the V and OC shells they can be efficiently represented by the one-electron $`j`$-dependent RECP terms.

Let us distinguish the following types of many-electron functions which are considered in a double-group symmetry:

* $`\{\mathrm{\Phi }_I^{(0)\mathrm{\Lambda }S},E_I^{(0)\mathrm{\Lambda }S}\}_{I=0}^{N^{(0)\mathrm{\Lambda }S}}`$ are reference SAFs (“Mains”) and their energies

$$E_I^{(n)\mathrm{\Lambda }S}=<\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}|𝐇^{(0)}|\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}>$$
(18)
at $`n=0`$ for those $`\mathrm{\Lambda }S`$-irreps which are of interest for the final spin-orbit CI (SO-CI) calculation;
* $`\{\mathrm{\Psi }_I^{(0)\mathrm{\Lambda }S},_I^{(0)\mathrm{\Lambda }S}\}_{I=0}^{𝒩^{(0)\mathrm{\Lambda }S}}`$ are some of the CI solutions (“$`\mathrm{\Lambda }S`$-roots”) and their energies

$$_I^{(0)\mathrm{\Lambda }S}=<\mathrm{\Psi }_I^{(0)\mathrm{\Lambda }S}|𝐇^{(0)}+𝐕^{corr}|\mathrm{\Psi }_I^{(0)\mathrm{\Lambda }S}>$$
(19)
in the $`\mathrm{\Lambda }S`$-irrep which diagonalize the $`(𝐇^{(0)}+𝐕^{corr})`$ in the subspace of Mains only;
* $`\{\mathrm{\Psi }_I^{(0)SO},_I^{(0)SO}\}_{I=0}^{𝒩^{(0)SO}}`$ are some of the SO-CI solutions (“SO-roots” which are of interest) and their energies

$$_I^{(0)SO}=<\mathrm{\Psi }_I^{(0)SO}|𝐇^{(0)}+𝐕^{corr}+𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>$$
(20)
which diagonalize the complete $`𝐇`$ Hamiltonian in the subspace of all Mains collected from all the $`\mathrm{\Lambda }S`$-irreps considered;
* $`\{\mathrm{\Phi }_I^{(1)\mathrm{\Lambda }S},E_I^{(1)\mathrm{\Lambda }S}\}_{I=0}^{N^{(1)\mathrm{\Lambda }S}}`$ are the singly-excited SAFs (SE-SAFs) and their energies (18) at $`n=1`$, i.e.

$$\mathrm{\Phi }_I^{(1)\mathrm{\Lambda }S}\{𝐏^{\mathrm{\Lambda }S}𝐚_p^+𝐚_q\mathrm{\Phi }_J^{(0)\mathrm{\Lambda }^{}S^{}}\}\{\mathrm{\Phi }_K^{(0)\mathrm{\Lambda }S}\}(p,q;J,K),$$
(21)
where $`𝐏^{\mathrm{\Lambda }S}=|\mathrm{\Lambda }S><\mathrm{\Lambda }S|`$ is a projector on the subspace of the $`\mathrm{\Lambda }S`$-states, $`𝐚_p^+(𝐚_q)`$ are the creation (annihilation) operators of one-electron states (spin-orbitals) $`\varphi _p(\varphi _q)`$. The SE-SAFs can be automatically selected because of their relatively small number;
* $`\{\mathrm{\Phi }_I^{(2)\mathrm{\Lambda }S},E_I^{(2)\mathrm{\Lambda }S}\}_{I=0}^{N^{(2)\mathrm{\Lambda }S}}`$ are the doubly-excited SAFs (DE-SAFs)

$$\mathrm{\Phi }_I^{(2)\mathrm{\Lambda }S}\{𝐏^{\mathrm{\Lambda }S}𝐚_p^+𝐚_q^+𝐚_r𝐚_s\mathrm{\Phi }_J^{(0)\mathrm{\Lambda }^{}S^{}}\}(\{\mathrm{\Phi }_K^{(1)\mathrm{\Lambda }S}\}\{\mathrm{\Phi }_L^{(0)\mathrm{\Lambda }S}\})(p,q,r,s;J,K,L)$$
(22)
and their energies (18) at $`n=2`$; a SAF $`\mathrm{\Phi }_I^{(2)\mathrm{\Lambda }S}`$ should be selected in accordance with some selection criteria to be used in the final SO-CI calculation. In principle, triple and higher excited sets of SAFs can be similarly defined.

The correlation operator, $`𝐕^{corr}`$, has the symmetry of the molecule and, therefore, can be rewritten as

$$𝐕^{corr}\underset{\mathrm{\Lambda }S}{}𝐏^{\mathrm{\Lambda }S}𝐕^{corr}𝐏^{\mathrm{\Lambda }S}.$$
(23)
It normally gives the most important contribution through the second-order Brillouin-Wigner PT energy correction in the basis set of $`\mathrm{\Phi }_J^{(n)\mathrm{\Lambda }S}`$ (after appropriate redefinition of $`𝐇^{(0)}`$ in the subspace of Mains, see Ref. ):

$$\underset{n=1,2}{}\underset{J}{}\frac{|<\mathrm{\Phi }_J^{(n)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Psi }_0^{(0)\mathrm{\Lambda }S}>|^2}{_0^{\mathrm{\Lambda }S}E_J^{(n)\mathrm{\Lambda }S}}$$
(24)
for the non-degenerate ground state $`\mathrm{\Psi }_0^{\mathrm{\Lambda }S}`$ with the exact energy $`_0^{\mathrm{\Lambda }S}`$ in the $`\mathrm{\Lambda }S`$-irrep (obviously, terms with $`n3`$ are automatically equal to zero because $`𝐕^{corr}`$ is a two-electron operator). A similar expression with the replacements $`\mathrm{\Psi }_0^{(0)\mathrm{\Lambda }S}\mathrm{\Psi }_I^{(0)\mathrm{\Lambda }S}`$ and $`_0^{\mathrm{\Lambda }S}_I^{\mathrm{\Lambda }S}`$ can be applied for excited states $`\mathrm{\Psi }_I^{\mathrm{\Lambda }S}`$ (some precautions should be taken concerning the degenerate states and the orthogonality constraints with respect to the lower-lying states with $`J<I`$). As a result, the first rows, columns and energies on the diagonal of the Hamiltonian matrix

$$<\mathrm{\Phi }_J^{(n)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Phi }_I^{(0)\mathrm{\Lambda }S}>,<\mathrm{\Phi }_J^{(n)\mathrm{\Lambda }S}|𝐇|\mathrm{\Phi }_J^{(n)\mathrm{\Lambda }S}>$$
(25)
for $`n=1,2`$ are usually employed in the selection procedures for SAFs $`\{\mathrm{\Phi }_I^{(1,2)\mathrm{\Lambda }S}\}`$ based on the nonrelativistic $`A_k`$ and $`B_k`$ approximations (when $`𝐇^{SO}`$ is not taken into account)  or on the multi-diagonalization scheme  for subsequent calculations of $`\mathrm{\Psi }_I^{\mathrm{\Lambda }S}`$. In spite of some differences between these selection schemes, they are not very essential for the final CI results if a high quality reference set (set of Mains) and a suitably small threshold are chosen.

For molecules with heavy and very heavy atoms, the $`𝐇^{SO}`$ operator can give large contributions to the energy both in second and in higher PT orders if a non-optimal set of Mains, $`\{\mathrm{\Phi }_I^{(0)\mathrm{\Lambda }S}\}`$, is chosen after an SCF calculation with the SO-averaged potentials (AREPs). The latter is the usual practice and the set of Mains generated in such a manner can be smaller than optimal for the case of large SO interaction. Therefore, not only second but third and maybe even higher PT order(s) can be important in the selection procedure for a “bad” set of the starting roots $`\mathrm{\Psi }_I^{(0)SO}`$. This means that the off-diagonal matrix elements of $`𝐇`$ between secondary many-electron basis functions (SE-, DE-SAFs) may be introduced into the selection procedure because $`𝐇^{SO}`$ is a substantially off-diagonal operator contrary to $`𝐕^{corr}`$:

$$<\mathrm{\Phi }_I^{(n)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Phi }_J^{(n^{})\mathrm{\Lambda }^{}S^{}}>(\mathrm{\Lambda }S)\text{ and }(\mathrm{\Lambda }^{}S^{})\text{ can be different,}n^{}\{n,|n\pm 1|\}.$$
(26)
In particular, $`𝐇^{SO}`$ gives zero matrix elements between SAFs belonging to the same $`\mathrm{\Lambda }S`$-irrep in the $`D_{2h}`$ or $`C_{2v}`$ symmetry groups.

For simplicity, let us consider the selection scheme based on the $`A_k`$ approximation (24).

In the nonrelativistic-type selection scheme, a SAF $`\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}`$ is selected in a $`\mathrm{\Lambda }S`$-irrep if

$$\frac{|<\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Psi }_I^{(0)\mathrm{\Lambda }S}>|^2}{E_J^{(1,2)\mathrm{\Lambda }S}_I^{(0)\mathrm{\Lambda }S}}\delta E_T^{\mathrm{\Lambda }S},$$
(27)
where $`IN^{SO}`$ and $`\delta E_T^{\mathrm{\Lambda }S}`$ is a threshold criterion for the energy selection scheme in the $`\mathrm{\Lambda }S`$-irrep. In (27) we have replaced the exact $`_I^{\mathrm{\Lambda }S}`$ energies by the approximate $`_I^{(0)\mathrm{\Lambda }S}`$ values (19) that corresponds to the Rayleigh-Schrödinger PT case. Such a simplification is justified for small $`\delta E_T^{\mathrm{\Lambda }S}`$ and good reference states.

In a SO-CI calculation within some relativistic double-group irrep, substitutions for the reference state $`(\mathrm{\Psi }_I^{(0)\mathrm{\Lambda }S}\mathrm{\Psi }_I^{(0)SO})`$ and the perturbation $`(𝐕^{corr}𝐕^{corr}+𝐇^{SO})`$ should be used in the previous expression, so that

$$\frac{|<\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐕^{corr}+𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>|^2}{E_J^{(1,2)\mathrm{\Lambda }S}_I^{(0)SO}}\delta E_T^{SO},$$
(28)
where $`\delta E_T^{SO}`$ is a selection threshold for $`\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}`$ to be used in the subsequent SO-CI calculation.

In more detail, the matrix element in the PT numerator of the above formula can be rewritten as

$`|<\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Psi }_I^{(0)SO}>|^2`$ (29)
$`+|<\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>|^2`$ (30)
$`+2\mathrm{}(<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>)`$ (31)
by taking into account eq. (23) in the calculation of the matrix elements for $`𝐕^{corr}`$, contrary to those for $`𝐇^{SO}`$. In spite of mixing different $`\mathrm{\Lambda }S`$-states due to $`𝐇^{SO}`$, the number of non-zero matrix elements with $`𝐇^{SO}`$ in eq. (31) is usually relatively small because the SO interaction is a one-electron operator (see eq. (26)) which is very localized compared with the long-range Coulomb interaction. Thus, one can see that the nonrelativistic-type selection due to $`𝐕^{corr}`$ with respect to $`\{𝐏^{\mathrm{\Lambda }S}\mathrm{\Psi }_I^{(0)SO}\}`$ in each considered $`\mathrm{\Lambda }S`$-irrep and automatic selection of all SE-SAFs $`\{\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}\}`$ (21) can be efficiently applied instead of eq. (28). It must be emphasized that contrary to the selection schemes in the nonrelativistic case, SE-SAFs should be generated with respect to the Mains from all the used $`\mathrm{\Lambda }^{}S^{}`$-irreps. In a more simplified treatment, the automatic selection of SE-SAFs can be done with respect to a subset of the most important Mains, e.g. having largest CI-coefficients in the $`\mathrm{\Psi }_I^{(0)SO}`$ roots.

Next let us consider the terms from the third-order PT energy (PT-3) for SAFs $`\{\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}\}`$ which can be essential for the SO selection procedure. Below we shall discuss only matrix elements in the PT numerators of the corresponding PT-3 terms because specific expressions for the energy denominators are not essential for our analysis and conclusions. For simplicity, we shall omit the terms conjugate to those considered.

The first two types of the PT-3 matrix elements are:

$`<\mathrm{\Psi }_I^{(0)SO}|𝐇^{SO}|\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Phi }_K^{(1)\mathrm{\Lambda }^{}S^{}}><\mathrm{\Phi }_K^{(1)\mathrm{\Lambda }^{}S^{}}|𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>,`$ (32)
$`<\mathrm{\Psi }_I^{(0)SO}|𝐇^{SO}|\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Phi }_L^{(1)\mathrm{\Lambda }S}><\mathrm{\Phi }_L^{(1)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>.`$ (33)
The first intermediate state, $`\mathrm{\Phi }_J^{(1)\mathrm{\Lambda }S}`$, is a test SE-SAF and the indices for other intermediate SAFs run over all the allowed ones. The PT-3 terms summed over the indices of the second intermediate state give contributions (together with the conjugate terms) for the selection of the test SE-SAF. However, the SE-SAFs can be selected automatically and these terms are out of our particular interest.

The following matrix element type

$$<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Phi }_K^{(1)\mathrm{\Lambda }^{}S^{}}><\mathrm{\Phi }_K^{(1)\mathrm{\Lambda }^{}S^{}}|𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>.$$
(34)
can be used for the selection of $`\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}`$ and $`\mathrm{\Phi }_K^{(1)\mathrm{\Lambda }^{}S^{}}`$ when summing over another set of intermediate states in the PT-3 expression. As one can see, this term can be used for the selection of both SE-SAFs and DE-SAFs. The above expression is quadratic in the (large) $`𝐇^{SO}`$ interaction contrary to the remaining terms considered below. The contribution of the terms with matrix elements (34) can be essential and their use for the selection of DE-SAFs $`\mathrm{\Phi }_J^{(2)\mathrm{\Lambda }S}`$ can be important for a subsequent SO-CI calculation.

The following matrix element types contain a second order perturbation in $`𝐕^{corr}`$ and, therefore, we can suggest that in general they are less important for our consideration than the above terms:

$`<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Phi }_L^{(1)\mathrm{\Lambda }S}><\mathrm{\Phi }_L^{(1)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Psi }_I^{(0)SO}>,`$ (35)
$`<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Phi }_K^{(1,2)\mathrm{\Lambda }^{}S^{}}><\mathrm{\Phi }_K^{(1,2)\mathrm{\Lambda }^{}S^{}}|𝐕^{corr}|\mathrm{\Psi }_I^{(0)SO}>.`$ (36)
These terms, together with the conjugate ones, can be used for the selection of $`\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}`$ and $`\mathrm{\Phi }_L^{(1)\mathrm{\Lambda }S}`$. The term

$$<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(2)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(2)\mathrm{\Lambda }S}|𝐇^{SO}|\mathrm{\Phi }_K^{(2)\mathrm{\Lambda }^{}S^{}}><\mathrm{\Phi }_K^{(2)\mathrm{\Lambda }^{}S^{}}|𝐕^{corr}|\mathrm{\Psi }_I^{(0)SO}>$$
(37)
can be analyzed separately because it contains both the intermediate states as DE-SAFs. In general, it is more difficult to take such terms into account in the selection procedure, because of the large number of tested DE-SAFs. We should note, however, that when a tested $`\mathrm{\Phi }_J^{(2)\mathrm{\Lambda }S}`$ DE-SAF is fixed, the other intermediate states, $`\{\mathrm{\Phi }_K^{(2)\mathrm{\Lambda }^{}S^{}}\}`$, are those DE-SAFs which are only singly excited with respect to the tested one. Therefore, the number of them will not be very high.

For completeness, the matrix element type which is cubic in the $`𝐕^{corr}`$ perturbation should be listed:

$$<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(1,2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Phi }_K^{(1,2)\mathrm{\Lambda }S}><\mathrm{\Phi }_K^{(1,2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Psi }_I^{(0)SO}>.$$
(38)
This term is of nonrelativistic type and it is out of our particular interest because it does not contain the $`𝐇^{SO}`$ perturbation. Again, we can separate the term

$$<\mathrm{\Psi }_I^{(0)SO}|𝐕^{corr}|\mathrm{\Phi }_J^{(2)\mathrm{\Lambda }S}><\mathrm{\Phi }_J^{(2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Phi }_K^{(2)\mathrm{\Lambda }S}><\mathrm{\Phi }_K^{(2)\mathrm{\Lambda }S}|𝐕^{corr}|\mathrm{\Psi }_I^{(0)SO}>.$$
(39)
from the previous one only because the latter contains both the DE-SAF intermediate states.

We should emphasize that the terms containing SE- or DE-SAFs in the intermediate states of the PT-3 expressions are not taken into account in the $`B_k`$ and multi-diagonalization selection procedures, although these schemes include, in fact, contributions of higher than the second-order PT terms.

When analyzing the above PT-3 terms, it can be concluded that if one replaces the reference SO roots, $`\mathrm{\Psi }_I^{(0)SO}`$, by new reference states, $`\mathrm{\Psi }_I^{(0+1)SO}`$, which diagonalize the complete Hamiltonian $`𝐇`$ for the sets of both Mains and SE-SAFs taken together, and applies the selection criterion based on the second-order PT (28), then the main part of the above PT-3 terms will be taken into account in such a selection. An exception occurs for terms (37) and (39), but in general they are thought to be less important than the other third-order PT terms.

In a more sophisticated treatment, the reference $`\mathrm{\Psi }_I^{(0+1^{})SO}`$ SO states can be generated when diagonalizing $`𝐇`$ for the sets of Mains and those SE-SAFs ($`\{\mathrm{\Phi }_J^{(1^{})\mathrm{\Lambda }S}\}`$), which are automatically generated with respect to the most important subset of Mains ($`\{\mathrm{\Phi }_I^{(0^{})\mathrm{\Lambda }S}\}`$). The latter subset can be selected from a preliminary CI calculation for the set of Mains, e.g. in a basis of configurations with the highest CI coefficients in $`\mathrm{\Psi }_I^{(0)SO}`$, and so on. This is worthwhile in order to reduce the number of SAFs in the resulting reference states $`\mathrm{\Psi }_I^{(0+1^{})SO}`$ rather than in $`\mathrm{\Psi }_I^{(0+1)SO}`$, thus reducing the selection time which can otherwise be very large.

We should also note that the trial SE- and DE-SAFs, which are tested in the above selection procedure, are generated only for the set of Mains and not for the $`\{\mathrm{\Phi }_I^{(0+1^{})\mathrm{\Lambda }S}\}`$ set. Therefore, the number of the tested configurations and the selection time are reasonably limited. If the number of configurations used in $`\{\mathrm{\Phi }_I^{(0+1^{})\mathrm{\Lambda }S}\}`$ is not high, one can extend the set of Mains by including the above subset of SE-SAFs, thus obviously enlarging the set of the consequently generated and tested SE- and DE-SAFs.

Again we should emphasize that it is not necessary to use the third-order PT or the suggested automatic selection of SE-SAFs in a selection procedure if a fairly good set of the reference roots $`\mathrm{\Psi }_I^{(0)SO}`$ is used, i.e., if they provide good approximations to the required solutions $`\mathrm{\Psi }_I^{SO}`$. In particular, if the $`\{\mathrm{\Psi }_I^{(0)SO}\}`$ set is obtained from a preliminary series of SO-CI calculations of the studied states, this can be superfluous.

As an alternative to the above selection schemes with respect to the PT energy, the PT expressions for the CI coefficient of a trial SE- or DE-SAF can be also explored. Applying the above PT analysis to the case of the $`\mathrm{\Psi }_I^{(0+1^{})SO}>`$ reference state, a $`\mathrm{\Phi }_J^{(1^{\prime \prime },2)\mathrm{\Lambda }S}`$ SAF is selected if its CI coefficient $`C_J^{(1^{\prime \prime },2)\mathrm{\Lambda }S}`$ satisfies the inequality

$`|C_J^{(1^{\prime \prime },2)\mathrm{\Lambda }S}|C_{min},`$

where $`C_{min}`$ is the selection threshold for the CI coefficients and

$$C_J^{(1^{\prime \prime },2)\mathrm{\Lambda }S}=\frac{<\mathrm{\Phi }_J^{(1^{\prime \prime },2)\mathrm{\Lambda }S}|𝐇|\mathrm{\Psi }_I^{(0+1^{})SO}>}{_I^{(0+1^{})SO}E_J^{(1^{\prime \prime },2)\mathrm{\Lambda }S}},$$
(40)
is the first-order PT value for the CI coefficient of a tested SAF which is not included in the subset of the $`\mathrm{\Phi }_J^{(1^{})\mathrm{\Lambda }S}`$ reference SE-SAFs.

Such a means of selection can be preferable if those properties of primary interest cannot be calculated from potential energy curves or surfaces. Moreover, the PT selection with respect to both the energy and the CI coefficients can be applied simultaneously if the properties are of different nature.

## V Calculations

In the CI calculations of Tl and TlH we used the MRD-CI package  combined with the SO selection codes based on the scheme described above. Our test calculations have shown that spin-orbit selection is very helpful for preparation of appropriate sets of Mains and for reducing effort in the final CI calculations with an optimal set of selected SAFs.

### A Spin-orbit splitting in the ground state of Tl

Calculations for the Tl atom were performed to optimize the basis set and the level shift GRECP parameters for the 21e/8fs-GRECP, i.e. the 21 electron GRECP with 8 electrons occupying the frozen OC pseudospinors, $`5s_{1/2}`$ and $`5p_{1/2,3/2}`$. The quality of the generated basis set is analyzed by calculating the $`{}_{}{}^{2}P_{1/2}^{o}{}_{}{}^{2}P_{3/2}^{o}`$ splitting for the ground state.

Before discussing the present results, it is worthwhile to to make some brief comments concerning numerous values for the Tl($`{}_{}{}^{2}P_{}^{o}`$) spin-orbit splitting calculated and published in the last years. It is well known that this quantity calculated at the one-configuration Dirac-Fock level agrees very well, within 100 cm-1, with the experimental value of 7793 cm-1 . However, such a good agreement results from the fortuitous cancelation of a number of large errors caused by the DF approximation. The situation changes dramatically when even the three outermost $`6s^26p^1`$ electrons are correlated. Some Tl($`{}_{}{}^{2}P_{}^{o}`$) splitting values from our 3e-CI calculations employing different codes and basis sets are given in Table I and lie between 7130 and 7210 cm-1. The corresponding 3e-CI results obtained by other groups after 1996 range from 6800 to 7800 cm-1 and such a large divergency can not be considered as satisfactory, because the ground state of Tl has a very simple configuration structure as compared to other heavy elements. We consider our calculated values of about 7200 cm-1 for this splitting as reliable for an approach in which the $`5d`$ spinors are frozen after the DF calculation of the nonrelativistically averaged $`6s^26p^1`$ configuration. Taking into account that a contribution of aproximately -100 cm-1 arises from the Breit terms, the deviation from the experimental value for the splitting is around 600–700 cm-1 (this size of an error can be justified theoretically). As will be shown below, the computated value for the Tl($`{}_{}{}^{2}P_{}^{o}`$) splitting can be significantly improved if $`5d`$ electrons are explicitly included in the calculations and the corresponding basis set contains functions with sufficiently high angular momenta.

The optimal basis set was selected in a series of MRD-CI calculations for Tl (with different sets of primitives and numbers of contracted $`s,p,d,f`$ and $`g`$ functions) to minimize the sum of energies for the ground $`{}_{}{}^{2}P_{1/2}^{}`$ and $`{}_{}{}^{2}P_{3/2}^{}`$ states. In these calculations, the SAFs were selected in the $`{}_{}{}^{2}B_{1u}^{}`$, $`{}_{}{}^{2}B_{2u}^{}`$, and $`{}_{}{}^{2}B_{3u}^{}`$ irreps of the $`D_{2h}`$ group (nonrelativistic-type degenerate $`{}_{}{}^{2}P`$ ground states belong to these irreps) because these doublets are strongly mixed by the SO interaction, resulting in the splitting of the ground $`{}_{}{}^{2}P`$ state. We have found that two $`g`$ functions should be added to the basis set, giving a contribution of about 9000 cm-1 to the $`{}_{}{}^{2}P_{}^{o}`$ ground state total energies. The resulting  basis set and GRECP parameters for Tl can be found on http://www.qchem.pnpi.spb.ru.

For the  basis set we have also performed MRD-CI calculations including SAFs from the $`{}_{}{}^{2}A_{u}^{}`$ irrep and SAFs with quartet multiplicity ($`{}_{}{}^{4}B_{1u}^{}`$, $`{}_{}{}^{4}B_{2u}^{}`$, $`{}_{}{}^{4}B_{3u}^{}`$, and $`{}_{}{}^{4}A_{u}^{}`$). In our calculations with different basis sets, their contributions have decreased the SO splitting by about 170 cm-1 and the total energy by about 2000 cm-1. One can see from Table I that this decrease is mainly caused by the $`p`$\- and $`d`$-components which arise from reexpansion of the leading spinor configuration in terms of the spin-orbit configurations. For good accuracy we can recommend the inclusion of $`\mathrm{\Lambda }^{}|S\pm 1|`$-irreps for the calculation of states having leading configurations in $`\mathrm{\Lambda }S`$-irreps.

In Table I some of our final MRD-CI results are collected together with the atomic relativistic coupled-cluster (RCC) results  obtained with a very large basis set. In these MRD-CI calculations altogether 627 Mains in three basic irreps and about 100 Mains in five additional irreps were involved and SE-SAFs were automatically generated for three Mains to prepare the reference $`\{\mathrm{\Psi }_I^{(0+1^{})SO}\}_{I=1}^3`$ states. Relatively small thresholds, $`T_1`$=0.03 and $`T_2`$=0.01 $`\mu E_h`$, are used in the final runs with the  basis set (for the $`T`$=0 threshold and full-CI extrapolations), thus selecting respectively about 190000 and 450000 SAFs altogether.

One can see that the best SO splitting calculated in the present work underestimates the experimental result  by about 400 cm-1 (recall that additionally about 90 cm-1 is due to the Breit interaction ). Analyzing our previous GRECP/RCC calculations of Hg  it can be concluded that this occurs due to the neglect of the OC-V correlations with the OC $`5p`$ and $`4f`$ shells, and to a lesser extent with $`5s`$ rather than due to the atomic basis set incompleteness, the GRECP errors or the restricted CI approximation. The OC-V correlation (contribution to the total energy) in Tl and Hg will have the same order of magnitude for respective pairs of correlated electrons (spinors).

We also studied the reliability of the linear $`T0`$ extrapolation procedure currently used in the MRD-CI code. In the final results of our MRD-CI calculations the corresponding correction gives the highest contribution to the cumulative error. So this is a bottleneck of the present Tl and TlH calculations with the large number of Mains.

### B Spectroscopic constants of the ground state in TlH

The explicit treatment of $`5d`$ electrons in precise TlH (TlX) calculations is necessary not only due to the strong correlation between these and the valence electrons of Tl, but also because of the substantial influence of relaxation-polarization effects in this shell on the bond formation. This cannot be very accurately taken into account by employing a polarization potential  in combination with, e.g., 3e-RECPs . The influence of other atoms (X) in a TlX molecule on the $`5p`$, $`5s`$ and $`4f`$ shells of Tl is significantly smaller and can be neglected if an accuracy of a few hundreds of wavenumbers for the excitation energies of low-lying states is sufficient. We neglected their contributions in calculation of the TlH spectroscopic constants.

In calculating spectroscopic properties for the TlH ground state (Table II) we used the contracted  basis set for thallium discussed above and the  set for hydrogen (see http://www.qchem.pnpi.spb.ru) contracted from the primitive (6,3,1) gaussian basis set of Dunning . The SAFs were selected in the $`{}_{}{}^{1}A_{1}^{}`$, $`{}_{}{}^{3}B_{1}^{}`$, $`{}_{}{}^{3}B_{2}^{}`$ and $`{}_{}{}^{3}A_{2}^{}`$ $`\mathrm{\Lambda }S`$-irreps of the $`C_{2v}`$ group because the triplet states are most strongly admixed by the SO interaction to the nonrelativistic $`{}_{}{}^{1}A_{1}^{}`$ (or $`{}_{}{}^{1}\mathrm{\Sigma }_{}^{+}`$ in $`C_\mathrm{}v`$) ground state producing the relativistic $`0^+`$ ground state in the double $`C_\mathrm{}v^{}`$ group.

We have performed three series of TlH calculations for 16 interatomic distances. In these runs, the $`\mathrm{\Psi }_0^{(0+1^{})SO}`$ reference SO states are generated with the MRD-CI code by diagonalizing $`𝐇`$ for the set of Mains and the SE-SAFs which are automatically selected with respect to the single configuration SCF ground state (calculated with the SO-averaged GRECP), giving a contribution of more than 90 % to the final wave function.

The first run is used for preparing an optimal set of Mains for the second series of SO-CI calculations. Only one SCF configuration which has the lowest energy in each $`\mathrm{\Lambda }S`$-irrep is included into the subspace of $`\mathrm{\Lambda }S`$ Mains and, consequently, the SO reference state consists of these SCF configurations and the automatically selected SE-SAFs with respect to the SCF configuration from the $`{}_{}{}^{1}A_{1}^{}`$ irrep.

Those SAFs were selected as Mains for the second run which had the highest CI coefficients in the first run. As a result, 37 Mains in all irreps together are employed in the second run. Relatively small thresholds, $`T_1`$=1.0 and $`T_2`$=0.1 $`\mu E_h`$, are used in the second run (for the $`T0`$ extrapolation ), thus causing about 20000 and 85000 SAFs to be selected in the $`\mathrm{\Lambda }S`$-irreps altogether.

In the most computationally consuming third run (with the set of Mains consisting of the SAFs having the largest CI coefficients in the wave function from the second run), about 320 Mains are used altogether and the thresholds are set at $`T_1`$=0.1 and $`T_2`$=0.05 $`\mu E_h`$. About 70000 and 130000 SAFs, respectively, were used in the $`\mathrm{\Lambda }S`$-irreps altogether in the final SO-CI calculations.

One can see from Table II that the basis set superposition error (BSSE) (see  and references) must be taken into account for an accurate computation of spectroscopic constants. The BSSE was studied in the Tl+ ion calculations for the same interatomic distances as in TlH and estimated also in the Tl- calculations for three distances, i.e. with the ghost H atom. The same molecular basis set as in TlH was used for both the Tl and H atoms. The contribution from BSSE to the total energy is decisive for the $`5d^{10}`$ and $`6s^2`$ shells considered in the case of Tl+, while its changing due to addition of the $`6p`$ electrons (which are bonding in TlH) can be considered as relatively small, because the difference in BSSE for Tl+ and Tl- is not significant in comparison with other errors. In the calculations of the spectroscopic properties with the counterpoise corrections (CPC), the calculated TlH points on the potential curve were corrected with the calculated BSSE for Tl+, i.e. for the $`5d,6s`$ shells taken into account.

One can see that after applying the T=0, FCI, and counterpoise corrections, the calculated properties are in very good agreement with the experimental data both in the second and third runs. The accuracy obtained is notably better than for other existing results for TlH (and not only for those presented in Table II). We suggest, however, that the very good agreement of the calculated $`D_e`$ with the experimental value can be fortuitous and the “real” (full CI) value can be notably different from the listed one because of the approximations made.

## VI Resume

The SO splitting in the ground $`{}_{}{}^{2}P`$ state of Tl is calculated by the MRD-CI method with the 21e-GRECP when $`5d^{10},6s^2`$ and $`6p^1`$ electrons are correlated and the $`5s^2`$ and $`5p^6`$ pseudospinors are frozen in the framework of the level shift technique. A  basis set is optimized for Tl and an underestimation of about 400 cm-1 is found for the SO splitting as compared with the experimental data.

Further improvement of the accuracy can be attained when correlations with the outer core $`4f,5p`$ and $`5s`$ shells of Tl and Breit effects are taken into account. We expect that this can be efficiently done in the framework of the “correlated” 21e/8fs-GRECP version in which 13 electrons are treated explicitly as in the present calculation. The inclusion of $`h`$-type functions is also desirable, as has been demonstrated for Hg in Ref. .

Fourteen electrons are correlated in the calculation of spectroscopic constants for the $`0^+`$ ground state of TlH and very good agreement with the experimental data is found.

The developed spin-orbit selection scheme and code are demonstrated to be efficient when large sets of basis functions and reference configurations are required in high-precision electronic-structure calculations.

###### Acknowledgements.

This work was supported by the DFG/RFBR grant N 96–03–00069 and the RFBR grant N 99–03–33249. AVT is grateful to REHE program of the European Science Foundation for fellowship grants (NN 14–95 and 22–95) to visit the laboratory of one of us (RJB), where part of the work was done. We are thankful to K. Shulgina and T. Isaev (PNPI) for writing some codes used for automatic generation of Mains. We are grateful to Dipl.-Ing. H.-P. Liebermann for the help in combining the MOLGEP and MRD-CI codes. We are also grateful to Dr. G. Hirsch (deceased) for his kind hospitality and invaluable help during visits to Wuppertal by AVT and NSM. The main part of the present calculations was carried out at the computer center of the Bergische Universität GH Wuppertal. JECS codes developed by PNPI quantum chemistry group were used for remote control of the calculations.


format: text

text:
# 1 Introduction

## 1 Introduction

The last decade has seen some major progress towards a solution of the Cosmic X-ray Background (CXRB) problem. As a result of instrumentation with improved sensitivity and energy resolution and the realization that a combination of heavily absorbed and unabsorbed active galaxy populations are required to make the CXRB , , it is now possible to easily reproduce the shape of the CXRB, at least in the $`3`$-100 keV band (e.g. ,  and other papers in these proceedings). Previously, in the 1980s, spectral measurements of heavily absorbed AGN spectra were poor and rare since such spectra were associated with X-ray weak Seyfert 2 galaxies, so models of the CXRB predicted a spectrum which was too steep compared to that which was observed. X-ray absorption provides the flattening of the spectra of some of the sources necessary to account for the CXRB spectrum. More details of the state of the art of the latest models can be found within these proceedings and references therein.

However, a few nagging and important problems remain. The first problem is what to identify this X-ray absorbed population with. Are they just the regular Seyfert 2 galaxies that we see in our local Universe? Increasing evidence since the early 1990s has pointed to a new, as yet unobserved, population of flat X-ray spectra, Narrow Emission-Line Galaxies (NELG), dominating the CXRB contribution at the faintest fluxes (e.g. , , ). There has been a massive effort in the past decade by many different groups around the world to try and directly resolve the X-ray background with deep surveys using X-ray imaging telescopes like ROSAT and ASCA followed up by optical identification programs, complemented with fluctuations analyses and cross-correlation studies (e.g. , , , , , , , , , , , , and these proceedings). The deepest ROSAT HRI observations (1.4 Ms) combined with optical follow-ups indicate that the earlier associations of the faintest CXRB contributions with NELG may be questionable ; but see also . About 80% of the soft X-ray background has been directly resolved into discrete sources in the ROSAT band and $`>80\%`$ of these are claimed to be broad-line AGN from HRI deep field data  but only $`31\%`$ from the deepest PSPC data .

In the hard band all except one study resolves $`30\%`$ of the 2–10 keV CXRB into a discrete source population with flat X-ray spectra. The one exception is the ASCA fluctuations study by  which goes down to a limiting 2–6 keV flux of $`10^{14}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$ and claims to resolve $`60\%`$ of the CXRB in this energy band. The measured photon index of the fluctuations is $`\mathrm{\Gamma }=1.26\pm 0.12`$. As pointed out by  the deeper the flux limit of the hard surveys, the flatter the X-ray slope of the CXRB contributors. Optical follow-up studies of all these hard surveys are still in progress but preliminary findings all come up with the same result: that a large percentage of the sources, maybe as much as 50% turn out to be broad-line flat X-ray spectra sources. This may be surprising to some because broad-line AGN are traditionally associated with steep X-ray spectra with little absorption. We shall show in this paper that this view is antiquated and thus it is premature to conclude from this that there must be a new, undiscovered population of sources making the CXRB. We note that at this flux level, all that can be said about the X-ray spectra is that they are flat; it is possible to model the X-ray spectra with just a single power-law only, which may in fact be mimicking a much more complex spectrum.

The second problem concerns the source counts, $`\mathrm{log}N/\mathrm{log}S`$ (i.e. the differential number of sources contributing to the CXRB at a given flux), or equivalently, the X-ray luminosity function (XLF). If we use the soft X-ray band XLF (0.5–2 keV) measured with, for example, ROSAT and use it to predict the XLF in the hard X-ray band (2–10 keV) we get a significant discrepancy compared with observation if we use the same assumptions about the source spectra which are used to model the CXRB spectrum. Likewise if we use the hard band XLF to predict the soft band XLF the result does not agree with observation. This problem has been noted by many authors (e.g. , , , , other papers in these proceedings). It has even been demonstrated that the discrepancy can be resolved if we do not insist that the broad-line objects contributing to the CXRB in the ROSAT band have steep X-ray spectra  in the hard band. Again, due to a misconception about the X-ray spectra of broad-line AGN, a new source population is invoked.

In summary, the two major problems are that a significant fraction of the flat X-ray spectrum sources making the CXRB are turning out to be broad-line AGN rather than narrow-line AGN, and that the soft XLF does not correctly predict the hard XLF and vice-versa. In the remainder of this paper we examine the critical assumptions which go into modeling the CXRB spectrum and XLF and show that we need look no further than our own backyard to find the broad-line AGN with flat X-ray spectra which are required to explain the CXRB.

## 2 Critical Assumptions

Currently, the latest AGN-synthesis models of the CXRB spectrum and XLF make a number of assumptions about the source spectra which may not be altogether correct or justified. We examine the most important ones below.

### 2.1 Intrinsic AGN continuum: Hard X-ray Slope

The intrinsic hard X-ray spectral photon index of both the type 1 and type 2 AGN is universally assumed to be $`\mathrm{\Gamma }=2`$ in AGN synthesis models. While it is true that the mean photon index measured for samples of AGN peak at this value (e.g. ), the dispersion is large. This is most dramatically illustrated in Fig. 6 of  which shows a range of photon indices ($`\mathrm{\Delta }\mathrm{\Gamma }>1`$; see Fig. 5 in  ) found for a sample of PG quasars. Aside from these examples, we know there are a significant number of sources which do not have $`\mathrm{\Gamma }=2`$ and may have intrinsically flat X-ray spectra, perhaps too many to ignore. This is not new of course: we have known for a long time that the one of the brightest quasars, 3C 273, has $`\mathrm{\Gamma }1.5`$. Also, it appears that high-redshift quasars may generally have flatter X-ray spectra .

The point is that we do not need to be constrained to assume $`\mathrm{\Gamma }=2`$ for AGN synthesis models of the CXRB because we do not know what value (or range) is relevant for the bulk of the sources that make the CXRB.

### 2.2 Intrinsic AGN continuum: Soft Excess

The latest AGN synthesis models still assume that a soft excess above the basic $`\mathrm{\Gamma }=2`$ power law is a ubiquitous feature of the intrinsic continuum in both type 1 and type 2 AGN (e.g. , ). Specifically, the models assume that the intrinsic power law steepens to $`\mathrm{\Gamma }=2.3`$ below 1.5 keV. Yet the latest studies of quasars show that this is simply not borne out by observation. For example  find only 5 out of 14 PG quasars show a soft excess in the ASCA band. Also,  find at most 2 out of 24 radio-loud and radio-quiet quasars with a soft excess in the ASCA band. Sure, ROSAT / PSPC studies all show that quasar spectra are significantly steeper in the 0.1–2.4 keV band than in the 2–10 keV band. Aside from systematic problems with the PSPC response (e.g. , ) there may well be a soft-excess component that is common below 0.5 keV. However, in the 0.5–2 keV band, which is relevant here, there is no evidence whatsoever from ASCA or ROSAT that a soft excess in the 0.5–2 keV band is ubiquitous. The lack of evidence for an intrinsic soft excess is even more severe in type 2 AGN, in which the intrinsic soft X-ray spectrum cannot be seen in most cases. It is either cut-off or is dominated by complex line-emission.

The ubiquity of a soft X-ray excess in AGN is a myth left over from the days of low-energy resolution spectroscopy in the eighties, before it was realised that photoionized absorption edges are common in AGN and that these can mimic a steepening of the soft X-ray spectrum with low-energy resolution data. However, the mimicking is only approximate and the two types of spectra differ in important ways (e.g. see ).

### 2.3 Intrinsic AGN continuum: High-energy cutoff

AGN synthesis models generally assume an exponential cutoff of around 300 keV or so (e.g. , ). The observational situation is that this is a very difficult parameter to measure with current instrumentation. A good review can be found in . It can be seen that the error bars are large and there is quite a range in the allowed cutoff values. In fact, one of the best measurements (for the bright AGN NGC 4151) is low at, $`70\pm 15`$ keV (see Table 1 in ). For the sources which make the bulk of the CXRB, we do not know whether low values or high values are relevant because we have no spectral information in that energy band for the weak sources identified in deep surveys. The distribution of the high-energy turnover in the sources that make the CXRB is important for reproducing the spectral bump at $`30`$ keV.

### 2.4 Absorption in Type 2 AGN

Several attempts have been made recently to characterize the absorbing column distribution in type 2 AGN (e.g. , , ), to use as input to AGN synthesis models of the CXRB. It has been known that absorbing columns deduced using data only below 10 keV, not taking into account higher-energy data, can be wrong. However, type 2 AGN have generally been too weak to obtain good spectra above 10 keV. Only recently has it been demonstrated, using BeppoSAX data, that the absorbing system parameters deduced previously for some well-known sources have been completely wrong. For example, the X-ray spectrum of Circinus galaxy was thought to be due to pure Compton reflection. However, BeppoSAX discovered a transmission component at hard X-ray energies . Also, Mkn 3 was previously thought to be Compton-thin but BeppoSAX utilizing higher energy data actually shows it to be Compton-thick . Since the absorbing column distributions compiled so far (e.g. in , , ) are based on previous data usually without the benefit of hard X-ray data, it has to be said that the true column density distribution is actually unknown. In any case, the distributions that have been compiled are only for local sources and not the ones that contribute to the bulk of the CXRB.

Another uncertainty in AGN synthesis models is how to treat type 2 QSO. A few have been discovered (e.g. , ) but their existence has also been disputed (e.g. ).

### 2.5 The ratio of Type 2/Type 1 AGN

AGN synthesis models either assume a value for this ratio or predict it. There are three problems with this ratio. The first is that it is based on optical classification and it is explicitly assumed in the models that type 1 AGN are not subject to X-ray absorption. The latter is simply not true (see §2.6). Thus, the model ratio and the observed ratio are two completely different quantities. The model ratio is an X-ray ratio but the observed ratio is an optical ratio. The second problem is that reports of the observed ratio cover an order of magnitude in range, from $`110`$ (e.g. see discussion in ). The third problem is that the ratio may be completely different at high redshifts. Currently, the type 2/type 1 ratio as it stands is not a very useful parameter.

### 2.6 Absorption in Type 1 AGN

All AGN synthesis models ignore X-ray absorption in type 1 AGN. Here we wish to demonstrate that even a small amount of absorption can make the 0.5–10 keV X-ray spectrum flatter than the CXRB spectrum. Complexity in the form of photoionized absorbers in type 1 AGN is common (e.g. , ) yet this is ignored in modeling the CXRB. Consider NGC 3516 for example, which is a well-known broad-line AGN. Fig. 1 shows the best-fitting absorption plus scattering model (including photoionization edge features) and residuals to ASCA data for NGC 3516. Such a model universally describes the main features of the ASCA spectra of many type 1 and type 2 AGN (e.g. ). The absorbing column is only $`3\times 10^{22}\mathrm{cm}^2`$. While the intrinsic power law is steep ($`\mathrm{\Gamma }1.8`$), the overall spectrum is extremely flat. In fact if this source was much fainter, say with a 2–10 keV flux of $`10^{13}\mathrm{ergs}\mathrm{cm}^2\mathrm{s}^1`$, a 100 Ks simulation shows that we would measure a slope of $`\mathrm{\Gamma }=1.03_{0.17}^{+0.18}`$ (assuming no extra absorption, which is what is done with very faint sources detected in deep surveys).

NGC 3516 is not an isolated example. There are many other examples in the literature of broad-line AGN with overall flat X-ray spectra. NGC 3227, NGC 4151 are two other notable examples. If these sources were very faint, simulations similar to those for NGC 3516 show that we would measure $`\mathrm{\Gamma }=0.40_{0.25}^{+0.20}`$ and $`\mathrm{\Gamma }=0.07_{0.28}^{+0.14}`$ respectively. We also note that  found a heavily absorbed quasar, PG $`1411+442`$, in their PG example with $`N_H10^{23}\mathrm{cm}^2`$ which consequently has a very flat X-ray spectrum.

### 2.7 Scattering Fraction

The value of the fraction of the observed continuum in type 2 AGN which is due to scattering and not subject to the absorption suffered by the direct continuum, is typically assumed to be $`24\%`$. Although authors of AGN synthesis models admit there might be a small but real range in this parameter, we suggest that the range is very large. For example we measure $`40\%`$ for NGC 2992 and $`0.2\%`$ for NGC 7172; these numbers represent approximately the extreme values we have measured for a large number of objects.

## 3 Alternative Models of the CXRB

Other categories of objects have been proposed to explain some or all of the CXRB. Notably starbursts (e.g. , , ), LINERS (e.g. ), elliptical galaxies, and naked ADAFs (, ). It remains to be seen whether starbursts can achieve a sizable contribution to the CXRB. Only some LINER X-ray spectra are flat but most are steep . Likewise, although elliptical galaxies may have a hard X-ray tail , their overall spectrum compared to the CXRB is too steep. As for naked ADAFs, none have been recognizably observed so far and they require fine-tuning of the accretion rate.

## 4 Conclusions

We conclude that the sources having broad optical lines and flat X-ray spectra, which are showing up in the deep surveys at the faintest fluxes, need not be a new, mysterious population of undiscovered types of X-ray source. The source spectral assumptions which go into current AGN synthesis models have not caught up with improved spectral measurements of AGN over the last decade. For example, X-ray absorption in type 1 AGN is common and soft excesses are not ubiquitous. Intrinsic type 1 X-ray spectra can also be flat. The flat X-ray spectrum type 1 AGN we already know of may simply be more abundant at the faintest fluxes. It will be interesting to see what sources Chandra will show up in deep hard surveys .

Questions A couple of important questions were raised.

(1) Isn’t this what people are doing anyway? What have you done that is different? While it is true that AGN synthesis models already include both absorbed and unabsorbed AGN, and it may be simply a matter of which are called type 1 and type 2, what people are doing is modeling all type 1 AGN without X-ray absorption and then try to predict the hard counts from the soft counts (or vice-versa) using $`\mathrm{\Gamma }=2`$. They then find a discrepancy which cannot easily be explained if one adheres to incorrect assumptions about the source spectra (§2). It has been demonstrated that the discrepancy can be resolved if one empirically assumes a flatter slope . This flatter slope for type 1 AGN may be real, as has been discussed above. What people are also doing is to attempt to give interpretation to the model type 2/type 1 ratio. One cannot take seriously interpretations of models of the CXRB spectrum and source counts if the source spectra are not modeled correctly.

(2) Earlier we saw some spectral fitting results of faint sources in deep surveys, but these had steep X-ray slopes, not flat.
Those spectral slopes quoted were obtained by fitting a model with intrinsic absorption included. This is not what is relevant for comparison with the slope of the CXRB. One must fit the spectrum of those sources without extra absorption if one is to compare with the CXRB. The spectra are flat.


text:
# 1 Bulk-mediated neutrino flavor oscillations with 𝑚=0.01 and 𝑚_{1,2}=1∓𝛿𝑚/2, for various 𝛿𝑚. As discussed in the text, these flavor oscillations are essentially “maximal” even though all mixing angles on the brane are vanishing.

Over the past few years, experimental evidence that neutrinos have mass and undergo oscillations has been pointing towards the existence of new physics beyond the Standard Model. Recent SuperK data on solar neutrinos  and atmospheric neutrinos , combined with the LSND accelerator measurements , has provided a new challenge for our understanding of the lepton sector of the Standard Model.

In some sense, the observation of lepton mixing restores a parallel structure between the quark and lepton sectors of the Standard Model. However, although the quark sector of the Standard Model is characterized by relatively small flavor mixing angles, recent data suggest that the lepton sector experiences large flavor mixings. Explaining the origin of these large mixings thus becomes a theoretical challenge. One important observation in this regard is the fact that right-handed quarks carry both color and electromagnetic charge, while right-handed neutrinos are neutral under all Standard-Model gauge symmetries. Thus, it is plausible that right-handed neutrinos have completely different origins than their quark counterparts.

One interesting possibility  is that the right-handed neutrino may be realized as a closed- (rather than open-) string state, and hence may propagate in the higher-dimensional bulk corresponding to a large extra spacetime dimension. The possibility of having large extra spacetime dimensions has been proposed  as a method of lowering the fundamental high-energy scales of physics such as the GUT, Planck, and string scales. Allowing the right-handed neutrino to feel these extra spacetime dimensions would have profound consequences for neutrino physics; indeed, it might then be possible to use neutrinos as a method of probing the physics of large extra dimensions. Various phenomenological aspects of higher-dimensional neutrino physics have been explored in Refs. .

One key issue concerns the embedding of flavor into such higher-dimensional scenarios. All flavor models in the literature thus far introduce one bulk neutrino for each of the three brane neutrinos, thereby extending flavor into the bulk. However, in such cases, the brane/bulk couplings become arbitrary $`3\times 3`$ mixing matrices whose parameters are undetermined. Moreover, the three bulk neutrinos can in principle correspond to different extra spacetime dimensions with different radii. Thus, one obtains a scenario with many undetermined parameters governing neutrino masses and mixing angles.

By contrast, we shall here introduce a “compact” model in which only one bulk neutrino is required. In other words, we shall consider flavor to be a property internal to the Standard Model, one which is restricted to the brane and which therefore does not extend into the bulk. Furthermore we shall assume flavor-blind couplings between our brane and bulk fields. Thus, in this model, only the physics directly on the brane is flavor-sensitive. Moreover, for simplicity we shall also assume that our brane theory is completely flavor-diagonal. Despite these facts, we shall see that neutrino flavor oscillations can still arise.

We begin by briefly describing the model. On the brane, we introduce three left-handed neutrinos $`\nu _i`$ ($`i=1,2,3`$); these are our flavor eigenstates. We assume that these left-handed neutrinos have corresponding Majorana masses $`m_i`$ on the brane, and we shall take the $`m_i`$ to be unequal. It is in this way that we shall distinguish between the different flavors on the brane. Note that for our purposes, we shall simply consider these Majorana masses to be arbitrary input parameters, and we shall not speculate on their origins or sizes. Despite the presence of these masses, however, we shall not introduce any explicit flavor mixings between the left-handed neutrinos on the brane. Thus, our theory on the brane will be completely flavor-diagonal.

Our bulk theory will be identical to that considered in Ref. . Specifically, we shall consider a single five-dimensional Dirac fermion $`\mathrm{\Psi }`$, which in the Weyl basis can be decomposed into two two-component spinors: $`\mathrm{\Psi }=(\psi _1,\overline{\psi }_2)^T`$. This Dirac fermion does not carry any flavor indices, and is therefore completely flavor-neutral. We impose the orbifold relations $`\psi _{1,2}(y)=\pm \psi _{1,2}(y)`$ where $`y`$ is the coordinate of the fifth dimension. In order to have a diagonal mass matrix, we also introduce the linear combinations $`N^{(n)}(\psi _1^{(n)}+\psi _2^{(n)})/\sqrt{2}`$ and $`M^{(n)}(\psi _1^{(n)}\psi _2^{(n)})/\sqrt{2}`$ for all $`n>0`$. Assuming that the brane is located at the orbifold fixed point $`y=0`$, we see that $`\psi _2`$ vanishes on the brane. The most natural brane/bulk coupling is therefore simply between $`\nu _i`$ and $`\psi _1`$. Although in principle each left-handed neutrino $`\nu _i`$ on the brane can have a different coupling to $`\psi _1`$, we shall consider a simple model in which this brane/bulk coupling, labeled by $`\widehat{m}`$, is completely flavor-universal.

Given these assumptions, our Lagrangian has the form

$`_{\mathrm{brane}}`$ $`=`$ $`{\displaystyle }d^4x{\displaystyle \underset{i=1}{\overset{3}{}}}\{\overline{\nu }_ii\overline{\sigma }^\mu D_\mu \nu _i+m_i(\nu _i\nu _i+\mathrm{h}.\mathrm{c}.)\}`$
$`_{\mathrm{bulk}}`$ $`=`$ $`{\displaystyle d^4x𝑑yM_s\left\{\overline{\psi }_1i\overline{\sigma }^\mu _\mu \psi _1+\overline{\psi }_2i\overline{\sigma }^\mu _\mu \psi _2\right\}}`$
$`_{\mathrm{coupling}}`$ $`=`$ $`{\displaystyle }d^4x{\displaystyle \underset{i=1}{\overset{3}{}}}(\widehat{m}\nu _i\psi _1|_{y=0}+\mathrm{h}.\mathrm{c}.).`$ (1)
Here $`M_s`$ is the mass scale of the higher-dimensional fundamental theory. By compactifying the Lagrangian (1) down to four dimensions we obtain

$``$ $`=`$ $`{\displaystyle }d^4x\{{\displaystyle \underset{i=1}{\overset{3}{}}}\overline{\nu }_ii\overline{\sigma }^\mu D_\mu \nu _i+\overline{\psi }_1^{(0)}i\overline{\sigma }^\mu _\mu \psi _1^{(0)}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(\overline{N}^{(n)}i\overline{\sigma }^\mu _\mu N^{(n)}+\overline{M}^{(n)}i\overline{\sigma }^\mu _\mu M^{(n)})`$ (2)
$`+\{{\displaystyle \underset{i=1}{\overset{3}{}}}m_i\nu _i\nu _i+\frac{1}{2}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}[\left({\displaystyle \frac{n}{R}}\right)N^{(n)}N^{(n)}\left({\displaystyle \frac{n}{R}}\right)M^{(n)}M^{(n)}]`$
$`+m{\displaystyle \underset{i=1}{\overset{3}{}}}\nu _i(\psi _1^{(0)}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}N^{(n)}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}M^{(n)})+\mathrm{h}.\mathrm{c}.\}\}`$
where $`m\widehat{m}/\sqrt{2\pi M_sR}`$ is the volume-suppressed brane/bulk coupling resulting from the rescaling of the individual $`\psi _1^{(0)}`$, $`N^{(n)}`$, and $`M^{(n)}`$ Kaluza-Klein modes.

Given the Lagrangian (2), we see that the Standard-Model flavor-eigenstate neutrinos $`\nu _i`$ will mix with the entire tower of Kaluza-Klein states of the higher-dimensional $`\mathrm{\Psi }`$ field, even though they do not mix directly with each other. Defining

$$𝒩^T(\nu _1,\nu _2,\nu _3,\psi _1^{(0)},N^{(1)},M^{(1)},N^{(2)},M^{(2)},\mathrm{}),$$
(3)
we see that the mass terms in the Lagrangian (2) take the form $`\frac{1}{2}(𝒩^T𝒩+\mathrm{h}.\mathrm{c}.)`$ where $``$ takes the symmetric form

$$=\left(\begin{array}{ccccccccc}m_1& 0& 0& m& m& m& m& m& \mathrm{}\\ 0& m_2& 0& m& m& m& m& m& \mathrm{}\\ 0& 0& m_3& m& m& m& m& m& \mathrm{}\\ m& m& m& 0& 0& 0& 0& 0& \mathrm{}\\ m& m& m& 0& 1/R& 0& 0& 0& \mathrm{}\\ m& m& m& 0& 0& 1/R& 0& 0& \mathrm{}\\ m& m& m& 0& 0& 0& 2/R& 0& \mathrm{}\\ m& m& m& 0& 0& 0& 0& 2/R& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right).$$
(4)
In this model, the bulk theory is flavor-neutral and the brane/bulk couplings are flavor-blind. Moreover, the theory on the brane is flavor-diagonal. Nevertheless, it is immediately apparent the three brane neutrinos will undergo flavor oscillations as a result of their indirect mixings with the bulk Kaluza-Klein neutrinos. In order to demonstrate this explicitly, we can determine the eigenvalues and eigenvectors of this mass matrix. It turns out that the eigenvalues $`\lambda `$ of the matrix (4) are given exactly as the solutions to the transcendental equation

$$\mathrm{tan}\pi \lambda R=\pi m^2R\underset{i=1}{\overset{3}{}}\frac{1}{\lambda m_i}.$$
(5)
For each solution $`\lambda `$ to (5), the corresponding mass eigenstate $`|\stackrel{~}{\nu }_\lambda `$ is then exactly given by

$`|\stackrel{~}{\nu }_\lambda `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N_\lambda }}}[\left({\displaystyle \underset{j=1}{\overset{3}{}}}{\displaystyle \frac{m}{\lambda m_j}}\right)^1{\displaystyle \underset{i=1}{\overset{3}{}}}{\displaystyle \frac{\lambda }{\lambda m_i}}|\nu _i+|\psi _1^{(0)}`$ (6)
$`+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\lambda }{\lambda k/R}}|N^{(k)}+{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\lambda }{\lambda +k/R}}|M^{(k)}]`$
where $`N_\lambda `$ is an overall normalization constant.

Note that the different flavor eigenstates $`|\nu _i`$ consist of different linear combinations of the different mass eigenstates $`|\stackrel{~}{\nu }_\lambda `$ for each $`\lambda `$. Thus, the different flavor eigenstates will experience a relative oscillation with each other. This oscillation is entirely “bulk-mediated” in the sense that there are no explicit flavor mixings on the brane; it is the presence of the higher-dimensional bulk which is completely responsible for inducing the flavor oscillations on the brane. Note that these oscillations are similar in spirit to the so-called “indirect” four-dimensional neutrino oscillations discussed in Ref. , except that here we have an infinite tower of higher-dimensional Kaluza-Klein states mediating our flavor oscillations and we have no explicit mixing angles on the brane.

We shall now discuss several phenomenological features of the “bulk-mediated” flavor oscillations inherent in this model. In general, we are interested in the probabilities $`P_{ij}(t)`$ that $`|\nu _i`$ oscillates into $`|\nu _j`$ as a function of time $`t`$. For $`i=j`$, this refers to flavor preservation on the brane, while for $`ij`$ this refers to flavor conversion on the brane. Note that in general, $`P_{ij}=P_{ji}`$ for all $`(i,j)`$. For phenomenological purposes, we are mostly interested in the probabilities $`P_{11}(t)`$, $`P_{22}(t)`$, and $`P_{21}(t)`$ where the subscripts $`(1,2,3)`$ signify $`(\nu _e,\nu _\mu ,\nu _\tau )`$ respectively. These are ultimately the three probabilities which are relevant for addressing constraints from solar, atmospheric, and long-baseline neutrino experiments.

In the following, we shall rescale all of our mass variables $`\{m_i,m\}`$ by $`R`$ so that these variables henceforth correspond to dimensionless quantities. Thus, our model has only four free dimensionless input parameters, $`\{m,m_1,m_2,m_3\}`$, while the radius $`R`$ serves as an overall length scale. These parameters ultimately determine not only the masses of the physical neutrinos, but also their relative effective mixing angles. We shall also define a dimensionless time variable, $`\stackrel{~}{t}t/(2pR^2)(1/2R^2)(L/E)`$, where $`Lct`$ is the spatial distance between the locations of neutrino production and neutrino detection, and $`Epc`$ is the neutrino energy. We can then plot our probabilities $`P_{ij}(t)`$ as functions of $`\stackrel{~}{t}`$.

Let us first consider the case when the brane/bulk coupling is extremely small, i.e., $`m1`$. In such cases, the mass matrix (4) is nearly diagonal, and there is relatively little mixing between the brane neutrinos and the bulk neutrinos. However, we find that large flavor oscillations can nevertheless be achieved for certain values of the input parameters $`\{m,m_i\}`$. As an explicit example, let us consider the case $`m=0.01`$, $`m_{1,2}=1\delta m/2`$, and $`m_3=5`$, with $`\delta m=m_2m_1`$ free to vary. The resulting neutrino oscillations are shown in Fig. 1. Note that as $`\delta m0`$, the resulting flavor oscillation becomes effectively maximal in the sense that we achieve full flavor conversion. In achieving this large flavor oscillation, we are exploiting a resonance between two approximately degenerate brane neutrinos $`|\nu _{1,2}`$ and the first excited Kaluza-Klein bulk neutrino $`|N^{(1)}`$, so that large effective flavor oscillations between the brane neutrinos are mediated indirectly through their oscillations with the bulk neutrino. Nevertheless, no mixing angles on the brane are required.

Even though this flavor mixing appears to resemble a simple two-state oscillation in the limit $`\delta m1`$, the time-averaged probabilities for each flavor are not 1/2, but 3/8. This reflects the fact that we are dealing with a three-state oscillation in which 1/4 of the total neutrino probability has been lost to the bulk neutrino $`|N^{(1)}`$, and is hence “sterile” from the point of view of flavor on the brane.

When all three of the brane neutrino flavors are approximately degenerate, so that $`m_1m_2m_3k\text{ZZ}`$, we find that the bulk neutrino mediates an effective four-state flavor oscillation which closely resembles a three-state oscillation on the brane. As an example, let us consider the case $`m=0.01`$, $`m_2=1`$, and $`m_{1,3}=1\delta m/2`$. The resulting probabilities for flavor oscillation and preservation are then shown in Fig. 2. Although the probabilities for flavor preservation and conversion are flavor-independent for short times (top plot of Fig. 2), a strong flavor dependence develops for later times (remaining plots). In fact, the corresponding time-averaged probabilities are given exactly by

$`\overline{P_{11}}=4/9,`$ $`\overline{P_{12}}=5/18,`$ $`\overline{P_{13}}=1/9,`$
$`\overline{P_{22}}=5/18,`$ $`\overline{P_{23}}=5/18,`$ $`\overline{P_{33}}=4/9.`$ (7)
Thus, if we start with $`|\nu _2`$, we see that we are equally likely to have conversion into each flavor, whereas if we start with $`|\nu _1`$ or $`|\nu _3`$, we are more likely to preserve the initial flavor, with equally spaced declining probabilities for conversion into nearby flavors. However, no matter which flavor of brane neutrino we begin with, we see that exactly 1/6 of the initial probability is ultimately lost into bulk neutrinos.

The above discussion pertains to only several specific choices of parameters $`\{m,m_1,m_2,m_3\}`$ in this model; a more complete discussion of these and other scenarios will be presented in Ref. . In all cases, however, we find that we can generate sizable neutrino oscillations even if the brane/bulk coupling parameter $`m`$ is extremely small and the brane theory is flavor-diagonal. In general, large brane/bulk oscillations are triggered when any of the $`m_i`$ take values that are near integers (in units of $`R^1`$). Note, however, that this does not require a large fine-tuning. When the extra dimensions are large, the corresponding Kaluza-Klein states are relatively closely spaced. Thus, regardless of the specific mass of a flavor neutrino on the brane, there is always a Kaluza-Klein state in the bulk which is effectively degenerate with it. Thus, such brane/bulk resonances are natural and can therefore mediate large brane/bulk neutrino oscillations. On the other hand, when $`m`$ is small, sizable flavor oscillations leading to flavor conversion still require approximate flavor degeneracies on the brane.

This situation changes dramatically as the brane/bulk coupling parameter $`m`$ becomes larger. In such cases, a larger population of bulk Kaluza-Klein neutrino states participates in the neutrino mixings, and it is no longer necessary to have flavor degeneracies on the brane in order to generate sizable bulk-mediated flavor oscillations.

As an example of this, let us consider the case when none of the $`m_i`$ are close to an integer or to each other. For example, we may take $`m_1=1/2`$, $`m_2=3/2`$, and $`m_31`$. For small values of $`m`$, this situation gives neither neutrino oscillations nor flavor conversion amongst the two lightest flavors. However, as $`m`$ increases, both neutrino oscillations and flavor conversions are generated. This is shown in Fig. 3, where we have plotted the time-averaged probabilities for neutrino flavor preservation and conversion as functions of the brane/bulk coupling $`m`$.

These effects can be further enhanced if any of the $`m_i`$ approaches an integer (in units of $`R^1`$), for in such cases there is an additional effect due to a strong resonance oscillation between the corresponding brane neutrino and its degenerate Kaluza-Klein partner in the bulk. This situation mirrors the phenomenon already discussed for the small-$`m`$ case. As the most extreme case of this, let us consider the situation in which all three of the brane neutrino flavors are approximately degenerate with a bulk neutrino. For concreteness we take $`m_2=k\text{ZZ}`$ and $`m_{1,3}=k\delta m/2`$. Choosing $`k=1`$ and $`\delta m=10^3`$, we obtain the time-averaged probabilities shown in Fig. 4.

This figure may be understood as follows. As $`m0`$, we expect that the preservation probabilities should go to 1 while the conversion probabilities should go to zero. This is indeed what happens, except that as $`m0`$, the time-averaged probability $`\overline{P_{22}}`$ goes not to $`1`$ but to $`1/2`$. This reflects the fact that for any non-zero value of $`m`$, there continues to exist a maximal two-state resonance between $`|\nu _2`$ and $`|N^{(1)}`$.

As $`m`$ increases into the range $`5\times 10^4\stackrel{<}{}m\stackrel{<}{}0.1`$, we see from Fig. 4 that the neutrino probabilities then settle into an effective four-state resonance oscillation pattern for which the time-averaged probabilities are given in (7). As we see from Fig. 4, these time-averaged resonance values persist over several orders of magnitude in $`m`$. Finally, as $`m`$ increases beyond $`0.1`$, larger numbers of Kaluza-Klein states begin to participate in the neutrino oscillations. This increases the probabilities for oscillations into bulk neutrinos, which in turn causes our time-averaged neutrino probabilities on the brane to decline relative to their values in (7). Remarkably, however, we see from Fig. 4 that even though these probabilities decline from the values given in (7), they still preserve their relative spacings as functions of $`m`$. Moreover, we see that the time-averaged probabilities on the brane do not fall to zero as $`m\mathrm{}`$, but instead take the non-zero asymptotic values

$`\overline{P_{11}}=7/18,`$ $`\overline{P_{12}}=2/9,`$ $`\overline{P_{13}}=1/18,`$
$`\overline{P_{22}}=2/9,`$ $`\overline{P_{23}}=2/9,`$ $`\overline{P_{33}}=7/18.`$ (8)
The fact that these asymptotic probabilities are non-zero illustrates that even when an infinite number of bulk Kaluza-Klein states contribute to the neutrino oscillations, and even when the brane/bulk coupling is infinite, we still do not lose all of our initial neutrino probability into the bulk states! Instead, as $`m\mathrm{}`$, we lose only 1/3 of the initial probability into bulk neutrinos. In other words, we see that 2/3 of the original probability remains “on the brane” even when $`m\mathrm{}`$. This is in fact a completely general feature of our model, and holds regardless of the specific values of the $`m_i`$. Moreover, this property generalizes: with $`n_f`$ flavors on the brane, we find that only $`1/n_f`$ of the initial neutrino probability is lost into bulk neutinos.*** Note that this property relies crucially on the fact that we have taken only one bulk neutrino in our model; in other words, we have not extended flavor into the bulk. By contrast, if we had taken a separate bulk neutrino for each brane neutrino, we would have instead found that all neutrino probability is ultimately lost into sterile bulk neutrinos as $`m\mathrm{}`$.

This surprising observation suggests that in our model, the large-$`m`$ case may be able to evade various four-dimensional bounds (such as those from supernova cooling rates) which ordinarily restrict the sizes of mixings with sterile neutrinos. Indeed, our results such as those in Eq. (8) apply for large values of $`m`$ for which it is not possible to employ constraints derived using a “perturbative” approach in which only the lightest few Kaluza-Klein bulk neutrinos mix significantly. Rather, an exact treatment is required in which the full Kaluza-Klein tower is incorporated simultaneously, for in our model it is the entire Kaluza-Klein tower that plays a role in reconverting bulk neutrinos back into neutrino flavors on the brane as the brane/bulk coupling becomes large. This issue, especially in conjunction with MSW matter effects, will be discussed further in Ref. .

Thus, to summarize, we see that it is possible to have large effective “bulk-mediated” flavor mixing on the brane even when all flavor mixing angles on the brane are zero. Moreover, this occurs even when the bulk theory is flavor-neutral and the brane/bulk couplings are flavor-blind. This illustrates that in higher dimensions, it may not be necessary to have large flavor mixing angles on the brane in order to accommodate experimental data; small (or even vanishing) flavor mixing angles on the brane may suffice. Furthermore, even though this model contains only five free parameters, it yields a surprisingly rich neutrino oscillation phenomenology, and exploits the higher-dimensional nature of the bulk in an essential way. Moreover, its essential structure involving only one bulk neutrino leads to a radically different neutrino phenomenology for large brane/bulk coupling than is possible in models with multiple bulk neutrinos, and suggests that this limit might actually evade supernova cooling constraints. Thus, this model might be profitably used as the basis of a detailed investigation of the experimental viability of various higher-dimensional neutrino oscillation mechanisms. This will be discussed further in Ref. .

Acknowledgments

We wish to thank C.P. Burgess, B. Campbell, L. Dixon, Y. Grossman, S. Lola, R. Mohapatra, K. Orginos, and M.H. Reno for discussions. This work was supported in part by the National Science Foundation under Grant PHY-0071054 and by the Department of Energy under Grants DE-FG02-95ER40906 and DE-FG03-93ER40792. We also wish to acknowledge the hospitality of the Aspen Center for Physics and the CERN Theory Division where portions of this work were completed.