Relativistic aspects of orbital and magnetic anisotropies in the chemical bonding and structure of lanthanide molecules

The electronic structure of magnetic lanthanide atoms is fascinating from a fundamental perspective. They have electrons in a submerged open 4f shell lying beneath a filled 6s shell with strong relativistic correlations leading to a large magnetic moment and large electronic orbital angular momentum. This large angular momentum leads to strong anisotropies, i. e. orientation dependencies, in their mutual interactions. The long-ranged molecular anisotropies are crucial for proposals to use ultracold lanthanide atoms in spin-based quantum computers, the realization of exotic states in correlated matter, and the simulation of orbitronics found in magnetic technologies. Short-ranged interactions and bond formation among these atomic species have thus far not been well characterized. Efficient relativistic computations are required. Here, for the first time we theoretically determine the electronic and ro-vibrational states of heavy homonuclear lanthanide Er2 and Tm2 molecules by applying state-of-the-art relativistic methods. In spite of the complexity of their internal structure, we were able to obtain reliable spin-orbit and correlation-induced splittings between the 91 Er2 and 36 Tm2 electronic potentials dissociating to two ground-state atoms. A tensor analysis allows us to expand the potentials between the atoms in terms of a sum of seven spin-spin tensor operators simplifying future research. The strengths of the tensor operators as functions of atom separation are presented and relationships among the strengths, derived from the dispersive long-range interactions, are explained. Finally, low-lying spectroscopically relevant ro-vibrational energy levels are computed with coupled-channels calculations and analyzed.


I. INTRODUCTION
A challenging question of molecular chemistry is an accurate description of inter-atomic and inter-molecular bonding at the quantum-mechanical level. This problem has attracted much attention but is not always resolved. Over the last decades, novel perspectives on the problem have relied on ultracold atoms and molecules. For example, quantum degenerate gases of atoms offer a unique platform on which to build and form small molecules in single internal state as they avoid unwanted system complexity. Ultracold gasses of atoms and molecules typically also allow for a high level of control and tunability and are well isolated from their surroundings.
As part of these developments experimental breakthroughs in realizing quantum gases of atoms with large magnetic moments [1][2][3][4][5][6][7][8][9] have also contributed. These atomic species tend to have a far more complex electronic structure than that of alkali-metal or alkalineearth species most often studied in the field. The magnetic lanthanides from dysprosium to thulium with their exceptionally large magnetic moments and large orbital momenta are extreme examples of such species. This experimental research relied on controllable and tunable anisotropic dipolar interactions between the atoms. The highly anisotropic short-range interactions between lanthanide atoms, however, remain poorly understood as they require knowledge of their chemical bonds. These systems form an excellent environment for explorations at the interface between quantum chemistry and atomic and molecular physics.
In previous research, we developed a successful model Hamiltonian to study the anisotropic interactions of bosonic Dy and Er in an external magnetic field [10,11] and in collaboration with the experimental groups of Drs. Ferlaino and Pfau we found and analyzed hundreds of magnetic Feshbach resonances in their collisions [12][13][14]. These resonances can be used to convert an atomic gas into a gas of highly-magnetic molecules as well as to study the threshold properties or the ultracold collisionenergy dependence of three-body relaxation [15]. These atom-atom interactions have also been studied in thulium (Tm) [16,17].
In spite of advances in the simulation of ultracold collisional interactions between heavy lanthanide atoms, the fundamental nature of the relativistic bond and shortrange electronic states in lanthanide dimers as well as in the even-heavier actinide dimers remains mostly unexplored. Precise knowledge of these interactions is clearly desirable for predicting their quantum vibrations and rotations. There exists an exception though. Substantial progress has been made in understanding interaction in the homonuclear diuranium molecule U 2 [18][19][20][21]. The latest studies [20,21] paid particular attention to the chemical bond of U 2 with its multi-orbital character. Relativistic and correlation effects using the Dirac equation for the electrons were fully incorporated by the authors of Ref. [21] and enabled them to determine the energies of the lowest electronic states of U 2 in the vicinity of the equilibrium separation. In addition, accurate ground-state potentials for heternuclear dimer molecules that include one open 4f-shell lanthanide atom and one non-lanthanide atom have become available [22][23][24][25][26][27][28].
The bond between two ground-state Er and two ground-state Tm atoms is the focus of this paper. The interactions between the j = 6 Er atoms and between the j = 7/2 Tm atoms are anisotropic and orientation dependent. Here, j is the total electron angular momentum of an atom. The anisotropy is a consequence of potential energy differences for different relative orientations of the electron angular momenta in the open 4f 12 and 4f 13 shells of Er and Tm, respectively. These 4f electrons lie beneath a closed 6s 2 shell so that these molecules are chemically similar but have distinct physical properties. Electron motion in lanthanides is strongly correlated and relativistic and spin-orbit coupling is strong.
The ground-state manifold of Er 2 and Tm 2 has a large number of electronic states. They are labeled by projection quantum number Ω with values up to 2j of the total dimer electron angular momentum on the symmetry axis of the molecule and well as other selection quantum numbers. Because of this complexity, the intermolecular interactions until now have not been accurately characterized. To fulfill these objectives we have performed, for the first time, relativistic configuration-interaction calculations of all Ω states as a function of interatomic separation R for Er 2 and Tm 2 using the DIRAC code [29]. These configuration-interaction calculations determine the short-range energy splittings among the 91 and 36 distinct adiabatic potentials of the Er 2 and Tm 2 dimers, respectively. Furthermore, we have setup an analytical spincoupling or spin-tensor representation of the short-range electronic potential surfaces for their use in determining rotational-vibrational levels in this paper and future improved simulations of the scattering of ultra-cold Er and Tm atoms. This representation has seven spin tensor operators and follows from the analytic form of the long-range anisotropic dispersion or van-der-Waals interaction. We find that the splittings among the Er 2 and Tm 2 potentials are dominated by a single anisotropic dipolar coupling between one of the atomic angular momenta and the mechanical rotation of the atom pair. We have also computed the long-range coupling strengths for the seven tensor operators based on all known atomic transition energies and transition dipole moments of Er and Tm. In fact, we find simple relationships among the seven spin-tensor operators contributing to the longrange interaction Hamiltonian.
Finally, we predict the relativistic Hund's case (c) structure of the energetically-lowest rotationalvibrational levels of the homonuclear Er 2 and Tm 2 dimers using a discrete-variable representation for the vibrational motion. We hope that our predictions will pave the way to spectroscopic studies of these complex and interesting molecules in the near future. Er 2 (Ω σ =12 g ) FIG. 1. Potential energies of the energetically-lowest "spinstretched" states of Er2 (black curve, Ωσ = 12g) and Tm2 (red curve, Ωσ = 7u) as functions of internuclear separation R.
The calculations are based on non-relativistic configurationinteraction and non-relativistic coupled-cluster theory for Er2 and Tm2, respectively. The zero of energy is at the dissociation limit of two ground-state atoms.

II. RESULTS AND DISCUSSION
A. Ground electronic states of Er2 and Tm2.
In this section we provide the relevant information on molecular electronic properties for two homonuclear lanthanide molecules, Er 2 and Tm 2 , using a two-step approach to determine short-range electronic potential surfaces for all molecular states that dissociate to Er or Tm atoms in the electronic ground states [Xe]4f 12 6s 2 ( 3 H 6 ) and [Xe]4f 13 6s 2 ( 2 F 7/2 ), respectively. These electronic configurations contain partially-filled or open submerged 4f and chemically-active 6s atomic orbitals. Computational details and justification of the two step process are presented in Sec. III as well as the Appendices.
In the first step of our study, we focus on the spinstretched states of Er 2 with Ω σ = 12 g and Tm 2 with Ω σ = 7 u , where subscripts σ = g and u for gerade or ungerade indicate the inversion symmetry of the electron wavefunction with respect to the center of charge. These states have the maximum allowed total electron spin quantum number S and the maximum projection quantum number Λ of the total electron orbital angular momentum along the internuclear axis, corresponding to the S = 2, Λ = 10 and S = 1, Λ = 6 states for Er 2 and Tm 2 , respectively. We have used non-relativistic configuration-interaction or coupled-cluster calculations to determine an accurate depth for the potential energy of these spin-stretched states. Figure 1 shows the spin-stretched potential energy curves for Er 2 and Tm 2 as functions of interatomic separation R. The relatively shallow potential depths of just under hc × 800 cm −1 reflects the covalent bond of the two closed 6s 2 orbitals. The equilibrium separations at the potential minima are R e = 8.7a 0 and 8.6a 0 for Er 2 and Tm 2 , respectively. Here, h is Planck's constant, c is the speed of light in vacuum, and a 0 = 0.0529177 nm is Relativistic Ω ± g potential energy curves with gerade symmetry for Er2 (panel a) and Tm2 (panel c) as functions of internuclear separation R near the equilibrium separation as obtained from electronic structure calculations. All potentials approach zero energy for R → ∞. Panels (b) and (d) show potential energies from panels (a) and (c) using the same line colors at the equilibrium separation as functions of projection quantum number Ω for Er2 and Tm2, respectively. Gerade Ω = 0 states are 0 + states. Potentials of states with Ω and −Ω are degenerate.
the Bohr radius. The depth and shape of these potentials is similar to that of the X 1 Σ + g state of the non-magnetic Yb 2 [30].
In a second step, we determine energy splittings among the Er 2 and Tm 2 potentials, using fully-relativistic configuration-interaction method implemented into DIRAC 2019 [29] that includes spin-orbit and anisotropic short-range interactions between atoms with two open 4f shells. The states are described by the projection quantum number of the total electronic angular momentum  el =  1 +  2 on the internuclear axis Ω, the gerade and ungerade symmetry, and a parity symmetry for Ω = 0 states. Here, Ω ranges from 0 to 12 for Er 2 and 0 to 7 for Tm 2 and labeling Ω ± = 0 + or 0 − indicates a symmetry with respect to the reflection of the electron wavefunction through a plane containing the internuclear axis. Figures 2(a) and (c) show the gerade relativistic potential energy surfaces (PES) for Er 2 and Tm 2 as functions of internuclear separation near their equilibrium separation. There are 49 gerade potentials for Er 2 and 16 for Tm 2 . For the separations shown in the figure and, in fact, for larger separations the splittings are less than 10 % of the depth of the potentials relative to the dissociation energy. The figures for the ungerade states is qualitatively similar and reproduced in the Appendices. There are 42 and 20 ungerade potentials for Er 2 and Tm 2 , respectively.
The splittings among the gerade relativistic potentials seem at first glance nontrivial. A pattern, however, emerges when we plot the potential energies at a single R near the equilibrium separation as functions of projection quantum number Ω, see Figs. 2(b) and (d). For both dimers the energetically lowest potential is a 0 + g state. We also observe that the splittings among states with the same Ω gradually decrease with increasing Ω. In fact, the potential energies at the equilibrium separation are arranged in parabolic shapes. Finally, for Er 2 the 0 + g state with the smallest well depth is nearly degenerate with the spin stretched 12 g state. A discussion of the origin of this pattern is given in the following subsection. Please note that the Ω = 7 state for Tm 2 has ungerade symmetry.

B. Spin tensor decomposition of Er2 and Tm2
PESs.
A tensor decomposition of the potential energy surfaces enables us to write PESs as weighted sums of spinspin coupling terms. It removes the need for a complicated evaluation of non-adiabatic couplings among potentials. We believe that the tensor format is essential for our molecular systems with their tens to hundred adiabatic channel potentials in the ground configuration.
Here, we apply the tensor decomposition technique developed in our previous study of scattering dynamics between ultracold Dy atoms [10] assuming that the molecular electronic wavefunction is well represented by superpositions of (anti-)symmetrized, parity-conserving products of atomic electronic ground states |j i m i or |j i Ω i for atom i = 1, 2. Atomic states are labeled by eigenvalues of the total electronic angular momentum operator  i , where projection quantum numbers along a space-fixed quantization axis are denoted by m i and those along the body-fixed internuclear axis by Ω i . For homonuclear systems j 1 = j 2 ≡ j. Nevertheless, subscripts 1 and 2 on operators  i and atomic states are kept to indicate the appropriate atom. (As always we omit the reduced Planck constanth in describing the eigenvalues of angular momentum operators.) The atom-atom interactions are then expressed as a sum of isotropic and anisotropic spin-tensor interactions. In principle, an infinite number of such interactions of ever increasing complexity exist. We, however, only include the seven low-rank tensors that describe the vander-Waals interaction at large interatomic separations [10]. These are with rank-k spherical tensor operators T (i) kq with components q, spherical harmonic functions C kq (R) with C kq (0) = δ q0 , andR is the orientation of the interatomic axis. Here, δ ij is the Kronecker delta. The seven T (i) kq correspond to three isotropic rank-0 tensorŝ three anisotropic rank-2 tensorŝ and a single anisotropic rank-4 tensor Thus N 0 = 3, N 2 = 3, and N 4 = 1 in Eq. (1). We have followed the ⊗ notation of Ref. [31] for combining spherical tensor operators, which is equivalent to the notation used in Ref. [32]. Then I is the identity operator, and [j 1 ⊗ j 2 ] kq denotes a tensor product of angular momentum operators  1 and  2 coupled to an operator of rank k and component q.
We obtain the strengths V (i) k (R) by a least-squares procedure minimizing the differences of the splittings with respect to the spin-stretched potential. We only do so for separations R ≤ 12a 0 for which DIRAC 2019 converged. All V k=2 (R) are consistent with zero. Thus, the dominant spin-tensor operators are either spin independent or a rank-2 tensor operator, corresponding to an R-dependent effective atomic quadrupole moment coupled to the rotation of the dimer. The spin-independent strength V      k=2 (R) can be found in the Appendices. The dominance of the isotropic spin-independent tensor operator is a consequence of the fact that the spatial extent of the 4f orbital is much smaller than that of the 6s orbital. In fact, the 4f orbital of one atom does not significantly overlap with that of a nearby Er or Tm atom even for R near the equilibrium separation.
We can now explain the origin of the patterns seen in Fig. 2. Including only the two dominant spin tensor operators in V ( R), eigenstates of V ( R) correspond to single |j 1 Ω 1 |j 2 Ω 2 + σ|j 2 Ω 2 |j 1 Ω 1 states. They have eigenenergies (9) with quadratic or parabolic dependences on Ω 1 and Ω 2 (and thus Ω). For Er 2 with positive V (1) k=2 (R) and integer j, Eq. (9) predicts that the Ω 1 = Ω 2 = 0 state and thus an Ω = 0 state has the lowest potential energy. In fact, for Er 2 this is a 0 + g state. In addition, Eq. (9) implies that the energetically-highest Ω = 0 state is degenerate with the sole spin-stretched Ω = 12 state. In fact, multiple degenerate adiabatic states with the same value for Ω 1 but opposite-signed values for Ω 2 exist.
For Tm 2 also with positive V (1) k=2 (R) but now halfinteger j, the model is quite satisfactory as well. Equation (9) predicts that states with |Ω 1 | = |Ω 2 | = 1/2 have the lowest energy. In this case, Ω σ is either 0 + g or 0 − u (both with Ω 1 = −Ω 2 = 1/2) or 1 u (with Ω 1 = Ω 2 = 1/2) and the ground state should be three-fold degenerate. In fact, the energetically-lowest level from the DIRAC calculations is a 0 + g state. Any removal of degeneracies is due to one or more of the five weaker spin-tensor operators not accounted for in Eq. (9).
Although the five weaker spin-tensor operators could not be reliably extracted from the least-squares adjustment to all splittings among the relativistic potentials of Tm 2 , additional analyses show that the spin-tensors in Eqs. (3) and (6) are the most important of the five. The first-order correction to the energy due to these two spin-tensor operators is with a positive value within the parenthesis. Thus the 0 + g state has a lower potential energy than the 1 u state.
C. Spin-tensors for long-range interactions.
The results shown in Figs. 2 and 3 focused on the deepest parts of the potentials. For scattering of ultracold atoms the long-range or large R part of the potential is equally important. The long-range form involves the van-der-Waals as well as magnetic and quadrupolar interactions.
In our model for V ( R) all seven strengths V 2 /R 3 describes the magnetic dipole-dipole interaction between the magnetic moments of the lanthanides atoms. Its strength D where g j is the electronic g-factor of the atomic ground state, α is the fine-structure constant, and E h = 4.359 74 × 10 −18 J is the Hartree energy. The magnetic dipole-dipole interaction is not captured by electronic structure calculations, but is relevant for scattering calculations.
The rank-4 strength V (1) k=4 (R) has a second long-range contribution as well. It approaches Q for homonuclear dimers that is solely determined by the atomic quadrupole moment Q = j i j i |Q 20 |j i j i of the m i = j i spin-stretched state of Er or Tm [33] and e is the positive elementary charge. The quadrupole moment for erbium was calculated in our previous paper [12] and equals 0.029ea 2 0 . For thulium Q is not available, but expected to be equally small compared to ea 2 0 . For rovibrational simulations with thulium we use Q = 0.
We have determined the C (i) k coefficients for Er 2 and Tm 2 from second-order perturbation theory in the electric dipole-dipole interaction using experimentallydetermined atomic transition frequencies and oscillator strengths or Einstein A coefficients as well as their reported uncertainties. We closely follow the calculations in Refs. [10,34] for the dysprosium dimer. The evaluation of the seven van-der-Waals coefficients is described in the Appendices. Values are given in Table I, while correlation coefficients are given in the Appendices. The relative for Er + Er and Tm + Tm sorted by value and then by rank k. The first three columns label the seven tensor operators. The last two columns give their value and its one-standard-deviation statistical uncertainties. The strength of tensor operator [j1 ⊗ j2]2 is √ 2 times larger than that of [j1 ⊗ j2]0. Similarly, the strengths of the last three tensor operators of each dimer are related by simple algebraic factors discussed in the text. For numerical convenience the strengths are given with additional digits.
Homonuclear Erbium dimer k i Operator T −0.009 080 527 0.0053 For Er 2 we observe that the absolute value of the magnetic dipole-dipole interaction |D 35a 0 the magnetic dipole-dipole interaction is the strongest anisotropic interaction, while for smaller R the effective quadrupole interaction of Eq. (5) is the strongest anisotropic interaction.
Less obvious from Table I is that we have been able to derive non-trivial algebraic relationships among the C (i) k thereby reducing the number of independent dispersion coefficients from 7 to 4. We find that showing that the spin-exchange strength and the effective dipole-dipole strength multiplying Eqs. (3) and (6), respectively, are related. Similarly, we find that relating the strengths of the three spin-tensors constructed from the two effective atomic quadrupole op- The derivation of these relations can be found in the Appendices.

D. Ro-vibrational eigenstates.
We finish our analyses of Er 2 and Tm 2 by computing their energetically-lowest ro-vibrational eigenstates. That is, we compute eigenstates of −h 2 ∇ 2 /(2µ) + V ( R), where µ = m/2 and m is the mass of the Er or Tm atom [35]. We discretize the radial component of the kinetic energy operator −h 2 ∇ 2 /(2µ) using the discrete-variable representation of Ref. [36]. Details regarding the selection of the spin and angular momentum basis and, in particular, the orbital or partial-wave angular momentum and total molecular angular momentum J of the two rotating atoms are given in Sec.. We present results from two calculations. One is based on potential V ( R) including only the two dominant spin tensorsT 2 , constructed by joining the electronic structure data with the long-range potentials, and one where all seven tensor operators are included. For latter case, the R-dependence of the remaining five weaker spin tensor operators is given by their long-range form for all R. Figure 4 shows three views of ro-vibrational eigenenergies of bosonic 168 Er 2 states near the minimum of the adiabatic potentials. As the nuclear spin of 168 Er is zero gerade basis states have even values for partial wave . Ungerade basis states require odd . The pattern of the energy levels in the hc×150 cm −1 energy range in Fig. 4a) can be understood from the seven 0 + g adiabatic potentials shown in Figs. 2a) and b) and Eq. (9). The vibrational energy spacing based on the harmonic approximation around the minima of the nearly-parallel adiabatic potentials is hc × 27.0 cm −1 so that vibrational levels v = 0, 1, · · · , 5 are visible in the figure. For each v the spacings among the seven 0 + g states are to good approximation found from √ 6V k=2 (R e ) = hc × 2.3 cm −1 at the equilibrium separation. The Ω 1 fine-structure thus overlaps with the vibrational structure. Figure 4b) shows v = 0 Ω ± g/u 168 Er 2 eigenstates over an energy region of only hc×15 cm −1 versus total molecular angular momentum J. based on calculations that include all spin tensors. Panel c) shows equivalent data including only the two strongest spin tensors. The level density in both two panels is large, although the level patterns are distinct with differences around hc × 1 cm −1 . In other words, the weaker spin-tensors can not be ignored in an accurate analysis of the lowest energy states of Er 2 . Surprisingly, we predict that the J = 10 rotational state of the v = 0 0 + g ground state has the absolute lowest energy when all interactions are included. The reason for this and other unexpected rotational progressions is the degeneracies of different Ω states inherent in the model of Eq. (9) as well as more accidental degeneracies due to the high level density. The coriolis force in a rotat- ing molecule breaks these (near-)degeneracies for states with Ω values that differ by one unit. For large J coupling matrix elements are on the order of B e Jj 1 , which can easily reach values of order hc × 1 cm −1 comparable to or larger than V (1) k=2 (R e ), even when B e is not. Degenerate perturbation theory then predicts "rotational" progressions of the form ±A e J + B e J(J + 1), where energy A e with |A e | B e can be computed on a case by case basis. Figure 5 shows the lowest eigenenergies of 169 Tm 2 near the minimum of its adiabatic potentials versus total molecular angular momentum J. Both gerade and ungerade states are shown and assigned Ω ± g/u labels. The nuclear spin ı of 169 Tm is 1/2, making the atoms bosons, and hyperfine interactions of the form a hf  1 · ı 1 /h 2 etc. mix gerade and ungerade states, although even and odd remain uncoupled. We do not include such interactions and gerade states with either even and odd must be shown. Similarly, ungerade states with even and odd exist. It is worth noting that the 169 Tm hyperfine constant a hf = hc × 0.0062 cm −1 [16] and hyperfine couplings can mostly be neglected in the analysis of the bound states on the scale shown in the figure.
The Tm 2 level structure is significantly simpler than that for Er 2 , even though the vibrational spacings and rotational constants of the two dimers are nearly the same. The dominant anisotropic spin-tensor interaction in Tm 2 is close to three times larger than in Er 2 , as seen in Fig. 3, thereby reducing the number of Ω ± σ states between the v = 0 and 1 vibrational states. The vibrational and rotational spectroscopic constants for Tm 2 are hc × 27.2 cm −1 and hc × 0.0096 cm −1 , respectively. Unlike for Er 2 the remaining five weaker spin-tensor operator play no significant role. Figure 5b) shows that the v = 0, J = 0 0 + g state is the absolute ground state. This 0 + g state has a "simple" B e J(J + 1) rotational progression. The nearby ungerade states have a less-conventional rotational progression. Here, this follows from the four-fold degeneracy for the lowest-energy states implied by Eq. (9). The four states have |Ω 1 | = |Ω 2 | = 1/2 and labels Ω ± σ = 0 + g , 0 − u and twice 1 u . In fact, a careful analytical analysis of the cen-trifugal and coriolis interactions within the degenerate manifold shows that in the limit J → 0 the ungerade levels have energies that lie tens of B e above that of the 0 + g state. For large J two of the three ungerade states become equal mixtures of 1 u and 0 − u . The ungerade state with the second-lowest energy remains of pure 1 u symmetry. Finally, we observe that the energetically-lowest ungerade state has J = 3.

III. METHODS
A. Spin-stretched electronic potentials.
The spin-stretched potentials V ss (R) for Er 2 and Tm 2 have been obtained from (partially-) spin-restricted single-reference non-relativisitic coupled-cluster calculations that included single, double and perturbative triple (RCCSD(T)) excitations [37]. For the calculations the total electron spin S is 2 and 1 for Er 2 and Tm 2 , respectively. In addition, the projection quantum number Λ of the total electron orbital angular momentum along the internuclear axis is 10 and 6 for Er 2 and Tm 2 , respectively.
We use the Stuttgart/Cologne "small-core" quasirelativistic effective core potentials (ECPs) developed for rare-earth elements (ECP28MWBSO) [38] to describe the twenty eight 1s, 2sp, and 3spd electrons of the Er and Tm atoms. The remaining electrons are described by the Relativistic Small Core Segmented (RSCSEG) [39] atomic basis set of quadruple-zeta quality developed for the ECPs. We extend the basis set with three diffuse Gaussian s functions with exponents 0.1495, 0.01, and 0.00412; two p functions with exponents 0.04895 and 0.0211; one d function with exponent 0.02799; and one f and g function both with exponent 0.1068, respectively. All exponents are in units of a −2 0 . In order to converge the preliminary self-consistent-field (SCF) calculations of neutral Er 2 and Tm 2 , we start from the SCF orbitals for molecular ions Er 4+ 2 and Tm 4+ 2 . Based on a comparison of results found with different basis set size, the one-standard-deviation uncertainty of the spin-stretched potential is about hc × 50 cm −1 at the equilibrium separation and drops to less than hc×1 cm −1 at R = 20a 0 . See tables in the Appendices for precise data on the spin-stretched potentials. A description of the extrapolation to R > 20a 0 can also be found there.

B. Relativistic calculation of potential splittings.
We use the direct relativistic configuration interaction (DIRRCI) method [40] in DIRAC [29] to determine the energy splittings between the relativistic adiabatic potential curves of Er 2 and Tm 2 dissociating to two groundstate atoms.
As in the non-relativistic calculations converging SCF calculations for the neutral dimers start from SCF or-bitals for Er 4+ 2 and Tm 4+ 2 . A reordering of the occupied 6s and open-shell 4f orbitals is then required in the input data for DIRAC. In practice, we only determine the SCF orbitals at R = 12a 0 in this manner. SCF orbitals for R < 12a 0 are found starting from orbitals for the neutral dimer obtained for a slightly larger R. We repeat the scheme to small R until the potentials are repulsive.
The active space in the DIRRCI calculations is solely composed of molecular orbitals arising from the 4f atomic shells. The 6s orbitals are kept doubly occupied and are not part of the active space. In addition, 5d orbitals remain unoccupied. These constraints balance the need for reasonable estimates of splittings among the relativistic potentials and a reasonable run-time and memory usage for the calculations. For state assignment it proved useful to determine the expectation values of the z-components of the total electronic orbital angular momentum and spin operators along the internuclear axis. DIRRCI calculations have been performed for 25 and 28 internuclear separations between 7a 0 < R ≤ 12a 0 for Er 2 and Tm 2 , respectively. At each of these R values, we determine the 91 relativistic potential energies and eigenstates for Er 2 . For Tm 2 we find 36 eigenpairs. An eigenstate is uniquely labeled by n, Ω ± g/u with n = 1, 2, 3, . . . for Ω ± g/u states ordered by increasing potential energy. We denote relativistic eigenenergies by U rel (R; n, Ω ± g/u ). For the 12 g and 7 u spin-stretched state for Er 2 and Tm 2 , respectively, there exists just one potential dissociating to two ground state atoms. For our identical ground-state atoms gerade states are 0 + states, while ungerade states are 0 − states. A discussion of di-atomic symmetries and molecular state labels can be found in Refs. [41] and [42].
With the constraints on the active space in the DIRRCI calculations, we sacrificed on the accuracy of the depth of the potentials. Those are mainly determined by excitations of electrons out of the 6s orbitals. In order to obtain accurate adiabatic potential energy curves V rel (R; n, Ω ± g/u ), we assume that the non-relativistic spinstretched potential V ss (R) is a good representation of the potential for the relativistic spin-stretched state. The adiabatic potentials for the other adiabatic states are then − U rel (R; n, Ω ± g/u ) − U rel (R; 1, Ω ss ) when R ≤ 12a 0 and Ω ss = 12 g and 7 u for Er 2 and Tm 2 , respectively. The uncertainties in the R ≤ 12a 0 calculations and extrapolation to R > 12a 0 using the long-range dispersion potentials are discussed in the Appendices. The potentials in Eq. (13) are used in the least-squares fitting to spin-tensor operators as described in the main text.
C. Basis sets in ro-vibrational state calculations.
We use the unit-normalized coupled spin and angular momentum basis for the calculation of ro-vibrational states of Er 2 and Tm 2 with |(j 1 j 2 )j el m el = m1m2 |j 1 m 1 |j 2 m 2 j 1 j 2 m 1 m 2 |j el m el and spherical harmonic functions Y m (R).
Here, j 1 j 2 m 1 m 2 |jm are Clebsch-Gordan coefficients. The total molecular angular momentum J = +  el is conserved,  el =  1 +  2 , and projection quantum numbers m x and M are with respect to a space-fixed coordinate system. Basis states with even and odd j el contribute to gerade and ungerade molecular states, respectively, and are not mixed by the molecular Hamiltonian. Similarly, the Hamiltonian does not couple basis states with even partial wave with those with odd . Atomic masses have been taken from Refs. [43,44] and atomic g-factors from Ref. [45].

IV. CONCLUSIONS
We have studied the electronic properties of two heavy homonuclear lanthanide molecules, Er 2 and Tm 2 . A hybrid non-relativistic/relativistic electronic structure approach was needed to overcome the computational challenges arising from the complexity of their open submerged 4f electronic shell structure partially hidden by a closed 6s 2 shell. This allowed us for the first time to determine a complete set of ground-state potentials for a wide range of interatomic separations.
A non-relativistic coupled-cluster calculation was used to determine the spin-stretched potential energy surfaces with the maximum allowed total electron spin S and projection quantum number Λ of the total electron orbital angular momentum along the internuclear axis. Then we used a relativistic multi-configuration-interaction calculation to determine the splittings among the potentials dissociating to two ground state atoms. There are 91 gerade/ungerade potentials for Er 2 (with Ωs from 0 to 12) and 36 potentials for Tm 2 (with Ωs from 0 to 7). We identified the splittings as due to different relative orientations of the angular momenta of 4f shell electrons.
To facilitate the application of our electronic structure predictions in spectroscopic and scattering dynamics studies we analytically expressed the potential energy operator for Er 2 and Tm 2 as a sum of a small number of spherical-tensor operators and elucidated the relationships between their electrostatic, relativistic, and magnetic dipole-dipole interactions. The most remark-able aspect of this analysis is that to good approximation the potential energy operator can be described with only two spin-tensor interactions, one isotropic and one anisotropic.
Finally, we computed the spectroscopically relevant lowest ro-vibrational eigenstates of Er 2 and Tm 2 . This data can be used as preliminary information for setting up spectroscopic studies of these exotic and technologically important systems.

Data availability statement
The data in this paper is self contained.

Appendix A: Table of content for Appendices
These appendices contain the input data and a derivation needed to reproduce the potentials presented in our article on the interactions between two ground-state erbium atoms and two ground-state thulium atoms. We give tables for non-relativistic spin-stretched potentials and relativistic anisotropic spin-tensor strengths. The main text has a graph of the potentials of gerade states of the two dimers. Here, we present the equivalent figure of potentials for ungerade states.
In addition, we derive the long-range van-der-Waals dispersion interactions for our high-spin Er and Tm atoms. We find that these interactions can be described in terms of seven spin-tensor operators, whose strengths or van-der-Waals coefficients are linearly dependent. In fact, only four independent dispersion coefficients exist. We also give tables of Er and Tm atomic transition frequencies and Einstein A coefficients or oscillator strengths on which the values of the dispersion coefficients in the table in the main text are based.
We use Planck's constant h and the speed of light in vacuum c in converting energies into wavenumbers.

Appendix B: Non-relativistic spin-stretched states of Er2 and Tm2
The spin-stretched potentials for Er 2 and Tm 2 dissociating to two ground-state atoms have been obtained from non-relativistic single-reference coupled-cluster calculations that included single, double and perturbative triple excitations (RCCSD(T)). The orbital basis sets for these calculations have been described in Sec. III. The total electron spin S and orbital-angular-momentum projection quantum number Λ are conserved quantities and for this state have their maximum allowed value. We have (S, Λ) = (2, 10) and (1,6) for Er 2 and Tm 2 , respectively.
The spin-stretched potential U ss (R) for Er 2 has been determined with coupled-cluster theory as implemented in CFOUR [46] at 59 separations between R min = 7.3a 0 and R max = 20a 0 as well as at R ∞ = 200a 0 in order to determine the dissociation energy of the potential. Here, a 0 = 0.0529177 nm is the Bohr radius. The spinstretched potential for Tm 2 has been determined using coupled-cluster theory as implemented in Molpro [47] at 72 separations between R min = 6.25a 0 and R max = 20a 0 as well as at R ∞ = 60a 0 . Potentials of Er 2 and Tm 2 up to R = R max are given in Tables II  and III, respectively. The spin-stretched potential for R < R min is found by linear extrapolation using the first two separations larger than or equal to R min . For R > R disp with R disp > R max we use the dispersive form V disp (R) = C 6,ss /R 6 + C 8,ss /R 8 + C 10,ss /R 10 , (B2) where the van-der-Waals coefficient C 6,ss is with j = 6 and 7/2 for Er 2 and Tm 2 , respectively based on two relevant values for C (i) k are given in the table in the main text and the derivation in this Appendix B. Coefficients C 8,ss and C 10,ss are fixed such that the dispersive form agrees with the potential energy from coupled-cluster theory at the two largest radial points R ≤ R max . We use R disp = R max + 0.5a 0 for both Er 2 and Tm 2 , add (R disp , V disp (R disp )) to the coupled-cluster data, and for R ∈ (R min , R disp ) interpolate this extended coupled-cluster data set times R 6 using the Akima spline [48]. The function R 6 V ss (R) varies significantly less than V ss (R). The fitted C 8,ss are consistent with typical values based on the induced quadrupole-quadrupole interaction for other di-atomic molecules [49] and the contributions from the five omitted dispersion terms is small compared to the uncertainties in the potentials.
The uncertainty budget of V ss (R) as function of R has two components. The first is the complete basis set error of the RCCSD(T) calculations. The second is that four-electron excitations might need to be included in our open-shell molecules. This corresponds to accounting for non-perturbative triple as well as quadruple excitations. Basis-set superpositions errors increase the depth of the potentials V ss (R), while non-perturbative triple and quadruple corrections often are of opposite sign, nearly cancel, but lead to shallower potentials for dimers [50,51]. Here, based on the differences of calculations with triple-and quadruple-zeta accuracy basis sets, we assume that the one-standard-deviation uncertainties of V ss (R) is 2 × u(C   We have used the direct relativistic configuration interaction (DIRRCI) method in DIRAC 2019 [29] to determine the energy splittings among the relativistic adiabatic potential curves of Er 2 and Tm 2 for R ≤ 12a 0 . Basis sets have been described in Sec. ??. In Sec. ??, we also described how the spin-stretched potential V ss (R) and the energy splittings are used to construct relativistic adiabatic potential curves V rel (R; n, Ω ± g/u ). We find a common uncorrelated one-standard-deviation uncertainty u(R) = hc × 10 cm −1 independent of R for all potential energy splittings. This follows from a compar-ison of DIRRCI calculations with different basis set size. In addition, the uncertainty in the splittings and that of V ss (R) are uncorrelated.
The V rel (R; n, Ω ± g/u ) were fit to an expansion in terms of seven spin-spin tensor operators with strengths V (i) k (R) defined in the main text. Only V (1) 0 (R) and V (1) 2 (R) were found to be statistically relevant and we finally decided to only present results with those two strengths as fitting parameters with the remaining five strengths set to zero. The reduced chi-square χ 2 ν of this adjustment is less than one for all R < 12a 0 so that the fit is consistent.
Tables IV and V contain values of the spin tensor strengths V where j = 6 and 7/2 for Er 2 and Tm 2 , respectively. The spin-tensor strength V For R > 12a 0 the relativistic configuration-interaction calculations do not converge. As shown in Fig. 4 in the main text, however, the adjusted anisotropic strengths V (1) 2 (R) at R = 12a 0 for the two dimers are already consistent, i.e. within our uncertainties, with its asymptotic van-der-Waals C (1) 2 /R 6 behavior. We then use the vander-Waals behavior for R > R rel with R rel = 12a 0 +0.5a 0 . A smooth connection of the strength between 12a 0 and R rel is ensured by adding point (R rel , C      Long-range van-der-Waals dispersion interactions have been derived and studied in many settings [34,52,53]. We repeat part of these derivations in order to explain the relationships among the strengths of the seven spintensor operators contributing the atom-atom interaction V ( R) for large separations R. In this section we rely on Ref. [32] for notation regarding angular momentum operators as well as the manipulation of these operators with one clearly-stated exception regarding reduced matrix elements.
We consider interacting atoms with electronic eigenstates |n bβ with total electronic angular momentum quantum numbers b, projection quantum numbers β on a space-or laboratory-fixed axis, and energies E nb that are independent of β. Label n further uniquely specifies states. The ground state of the atoms is |g jm .
The van-der-Waals potential operator between two ground-state atoms, |g 1 j 1 m 1 , g 2 j 2 m 2 = |g 1 j 1 m 1 |g 2 j 2 m 2 , is derived from (degenerate) secondorder perturbation theory in the anisotropic electric dipole-dipole interaction V dd ( R) and is given by d i is the rank-1 electric dipole moment operator of atom i = 1 or 2, C kq (R) is a spherical harmonic, and rank-k spherical-tensor operator T kq (R, S) ≡ [R ⊗ S] kq with the ⊗ notation of [31] is constructed from spherical-tensor operators R and S with rank r and s, respectively. The prime on the sums in Eq. (D1) indicates that the sums exclude the term where both atoms are in the ground state and b i = |j i −1|, . . . , j i +1. The energy denominator is negative and does not depend on projection quantum numbers. Finally, 0 is the vacuum electric permittivity. The sums in Eq. (D1) can be rearranged in several steps using Appendix VI of Ref. [32]. As operators d 1 and d 2 commute, we first note that with T kq (C 2 , C 2 ) = k0|2200 C kq and Clebsch-Gordan coefficient j 3 m 3 |j 1 j 2 m 1 m 2 . We omitted the dependence on orientationR of the spherical harmonics for clarity. The right hand side of Eq. (D4) is only nonzero for even k = 0, 2, or 4. Secondly, we note that where the d i have been grouped by atom and l i = 0, 1, or 2 and l 1 + l 2 is even for a nonzero value of the nine-j symbol · · · · · · · · · . Next, we isolate the sums over projection quantum numbers and labels n in Eq. (D1). For atom i, we define spherical tensor operators with rank l = 0, 1, or 2 and b = |j − 1|, . . . , j + 1. Here, gj||d i ||nb and nb||d i ||gj are reduced matrix elements of the electric dipole moment operator between the atomic ground state and excited state |n bβ . Crucially, for a ground-state atom we derive that gjm|B lq (b, j; i)|gjm = (D6) jm|jlm q (2b + 1)(2l + 1)W (j1j1; bl) based on the Wigner-Eckart theorem, Eq. (3.12) of Ref. [32], and symmetries of the Racah symbol W (abcd; ef ). These matrix elements are independent of label n of the excited state, but still depend on its total angular momentum b. Moreover, using the Wigner-Eckart theorem again, we realize that the m, m , and q dependences of gjm|B lq (b, j; i)|gjm are identical to those for the identity operator, atomic angular momentum operator j q , and dipole operator T 2q (j, j) for l = 0, 1, and 2, respectively. In fact, we find with rank-l operator O lq (i) = I, j iq /h, and T 2q (j i , j i )/h 2 for atom i and l = 0, 1, and 2, respectively. Here,h is the reduced Planck constant and the function M (b, j; l) is given by and M (b, j; l = 2) = 2b + 1 2j + 1 .

(D10)
We put everything together to find for two groundstate atoms where the matrix Here, e is the elementary charge, E h is the Hartree energy, and a 0 is the Bohr radius. The prime in the sums over labels n 1 and n 2 excludes the case where both atoms are in their ground state and we have introduced the more symmetric reduced matrix elements (j||d||j ) = (−1) j−j (j ||d||j) * = √ 2j + 1 j||d||j used by, for example, Edmonds in Ref. [54].
The allowed values for k, l 1 , and l 2 and the operators on the first line of Eq. (D11) lead to the seven spin-tensor operators defined in the main text. The last three lines of Eq. (D11) correspond to the van-der-Waals coefficients C (i) k . For example, the choice O l1 (1) = j 1 /h and O l2 (2) = j 2 /h leads to spin-tensor operators [j 1 ⊗ j 2 ] k0 /h 2 with k = 0 or 2 in the main text using the ⊗ notation of Ref. [31] for combining spherical tensor operators. This notation is equivalent to the notation used in Ref. [32]. Further analysis of Eq. (D11) shows that the operators with the same l 1 and l 2 but different k have related vander-Waals coefficients as the k dependence is isolated in the nine-j symbol and the Clebsch-Gordan coefficient in the second line. This leads to the relationships between the two spin-tensors with l 1 = l 2 = 1 and the three spintensors with l 1 = l 2 = 2 given in the main text.
Lists of currently available atomic transition energies E nb − E gj and observed Einstein A coefficients or oscillator strengths f for erbium and thulium atoms are given in Tables VI, VII, and VIII below. The relevant relationships between the reduced matrix elements (gj||d||nb) in Eq. (D12) on the one hand and A and f on the other are respectively. Here, α is the fine-structure constant. For Er 2 atomic transition data or lines from 48 excited states to the ground state are available. The majority of the data is taken from Refs. [55,56]. References [57] and [58] each supply one line, while for two other lines we rely on private communications [59,60]. For Tm 2 atomic transition data from 65 excited states to the ground state are available. Data have been taken from Refs. [45,61,62]. When line strength information of a transition is available from more than one source, the most accurate datum is chosen.
The uncertainties of and correlations among the vander-Waals coefficients follow from error propagation on Eqs. (D11) and (D12) and are dominated by the uncorrelated uncertainties of the Einstein A coefficients and oscillator strengths. Uncertainties in the transition energies give negligible contributions. Our values for the van-der-Waals coefficients are listed in the Table in the main text. Their covariances can be found in Table IX below.
TABLE VI. Atomic transition energies and Einstein A coefficients from exited states of erbium to its j = 6 ground state. Columns labeled ∆E, A, u(A), and j give transition energies, the value and one-standard deviation uncertainty of the Einstein A coefficients, and the total electronic angular momenta j of the excited states, respectively. A reference to the original data is given in columns labeled by "Ref. TABLE VII. Some thulium excited atomic eigen energies with respect to its j = 7/2 ground state and oscillator strengths f from the ground state to these excited states. The first column gives the transition energy. The second and third column are the value and one-standard deviation uncertainty of the oscillator strength, respectively. The fourth column gives the total electronic angular momentum j of the excited state. A reference to the original data is given in the last column. Relevant thulium lines for which Einstein A coefficients are available can be found in Table VIII.
38342.570 0.00169 0.00338 7/2 [62] 39019.090 0.00098 0.00196 9/2 [62] 39259.920 0.00262 0.00525 5/2 [62] 39580.720 0.00822 0.01644 7/2 [62] 39847.040 0.00192 0.00384 7/2 [62] 40101.720 0.00121 0.00242 9/2 [62] TABLE VIII. Atomic transition energies and Einstein A coefficients from exited states of thulium to its j = 7/2 ground state. Columns labeled ∆E, A, u(A), and j give transition energies, the value and one-standard deviation uncertainty of the Einstein A coefficients, and the total electronic angular momenta j of the excited states, respectively. A reference to the original data is given in columns labeled by "Ref.". Relevant thulium lines for which oscillator strengths are available can be found in Table  VII Table I of the main text. Coefficients are labeled by rank k and index i following the notation in the main text. Correlation coefficients with the remaining three dispersion coefficients follow from the exact algebraic relationships among the dispersion coefficients.