The Holohedrization Effect in Ligand Field Models

Abstract

The ligand field theory is an early and yet perennial class of quantum models accounting for the optical and magnetic properties of metal ions as a function of their environment in compounds. In the context of modern quantum chemistry, in order to predict properties from first principles, the ligand field paradigm can serve to illuminate the black box of heavy calculations, extracting heuristic meaning and causal roots. The genuine ligand field models are tacitly affected by an artificial feature, so-called holohedrization. It induces an inversion symmetry, even in cases where the local geometry does not show this element. This aspect received little attention over decades of using the ligand field Hamiltonians. In this work, we systematically investigate, assisted by state-of-the-art ab initio computer experiments, whether holohedrization is a hidden drawback of early models or if it also appears in realistic modeling. We found that the holohedrization trend also appears when using data from modern ab initio calculations.

1. Introduction

This work concerns a class of models pertaining to quantum mechanics applied to problems regarding optical and magnetic properties due to metal ions embedded in molecular compounds or solid-state lattices. The molecular entities forming the environment of the metal ions in molecules and crystals are called ligands [1,2]. Aiming to establish a causal relationship between properties and the geometry of local structure around a given metal ion, the ligand field (also called crystal field) theory [3,4] is a crude phenomenological model that appeared immediately after the birth of quantum mechanics. The “phenomenological” epithet denotes a rather simplified conceptual framework, retaining the physical essence of a given problem, based on parameters that are not computed from first principles but are accessible by adjustment to certain experimental data. To designate the same situation, the ligand field (LF) models are said to be based on effective Hamiltonians, suggesting that the underlying parameters are efficiently accounting for the main features of the system. The LF paradigm forms the conceptual ground of the structure–property qualitative correlation in the chemistry of d-type transition metal elements [5] and f-type compounds (lanthanides) [6]. The electronic structure regime usually experienced by metal ions in molecular compounds (called complex combinations) or lattices (oxides, halides, etc.) is called coordination bonding [1,2]. It is characterized by a weak or moderate perturbation of the free ion features.

It must also be added that the LF effective Hamiltonian accounts only for the one-electron components of the electronic structure (i.e., the interaction of electrons with the environment and their own kinetic energy). In the completion of LF methodology, the two-electron part (equating inter-electron Coulomb and exchange effects) is borrowed from the spectroscopy of atoms [7], with certain adjustments of the so-called Slater–Condon or Racah parameters [8,9,10]. However, in the following, we will not deal with two-electron parameterization.

To put things into perspective, one must point out that, nowadays, quantum computation methods can accurately reproduce the properties of small- and medium-sized molecules in the so-called ab initio mode, i.e., using only theory and pure principles. At the expense of various sorts of approximations, one may tackle larger systems. Computational methodologies developed around Schrödinger’s equation, and its Hamiltonian representation can be grouped into a branch labeled Wave Function Theory (WFT) [11]. Another class, Density Functional Theory (DFT) [12,13], is grounded on a theorem [14] stating that the electron density can be set as a primordial factor, provided that an operator called an exchange–correlation functional is known. The trouble is that no exact expression of such a conceptual panacea is available, having, in turn, tens of empirical versions. The quest for new density functional approximations and the increasingly good performances have led to an explosion of DFT usage and popularity in the last few decades as a practical route for applications. While molecular systems are approachable in the WFT frame, the applications in solid-state chemistry and physics are heavily based on DFT [15]. Aside from dependence on empirical formulas, an intrinsic limitation of DFT is that it cannot account for systems with degenerate ground states, e.g., molecules or atoms whose lowest levels consist of a set of equal energies. In the spherical symmetry of atoms or high-symmetry molecules, the degenerate ground state may occur in important situations, making this DFT drawback a non-negligible impediment. On the other side, as statistics show, many systems of interest lack symmetry and have a non-degenerate ground state. DFT cannot account for excited states (e.g., the objects determining the optical properties), although a technical trick known as Time-Dependent (TD) DFT may emulate information about a limited number of excited states (realized by single-electron orbital promotions) [16,17].

In turn, by multi-configurational procedures, WFT describes the degenerate and excited states, a goal coming with high demands on technical resources (time, memory, disk space), which may become prohibitive. Despite intrinsic limitations, DFT exerts fascination by virtue of its conceptual fundaments because electronic density is a physical measurable object for which physical intuition has clear clues. By contrary, since WFT descends from Schrödinger’s equation [18], it may seem to inherit the alleged lack of intuitiveness. Although puzzling paradoxes appeared at the birth of quantum mechanics [19] and are still vivid in modern [20,21,22,23] and actual research and debates [24,25], the WFT methodology is quite clear; nowadays, computational chemists can go for applications without hindering turmoil about meaning.

This overview of the large realms of computation and modeling of quantum systems serves to shed clear light on phenomenological models and effective Hamiltonians, like ligand field theory. According to the initial description above, these may look obsolete. After the emergence of realistic and powerful calculations via WFT or DFT, the LF approaches apparently receded from the first line of the theoretical description of metal-based compounds. However, LF came to a new life, rebranded as a tool to illuminate the black box of complicated WFT calculations [26,27,28,29]. Another practical utility of LF models can be observed with respect to the DFT frame. As previously briefed, DFT is not suited for certain problems with degenerate or quasi-degenerate ground states and excited levels [30], as needed in the magnetism and spectroscopy of metal ion systems. However, DFT can be wisely conducted to set numeric experiments yielding LF parameters. Then, subsequently, the spectrum of states needed in the description of properties can be emulated in a phenomenological LF frame [31,32,33].

From such a perspective, one may foresee the perennial value of LF phenomenological models. Divagating a bit, one may say that the LF models, already at the venerable age of almost one century, are here to stay for a long time. In several decades from now, DFT will not look as we know it today. In the future, either breakthroughs in the quest for the “holy Grail” of the exact exchange–correlation functional will create a new generation of DFT methods, or WFT will take over the whole realm of calculations on the wave crest of extremely powerful new computers. But, a phenomenological method, like LF, containing the quintessence of the physical problem at hand, will survive overall changes in the landscape of computational techniques, safeguarding the heuristic insight aimed by scientific inquiry.

2. Methods

The quantum calculations were performed with the GAMESS [34] suite (General Atomic and Molecular Electronic Structure System). The calculations are intentionally kept at the simplest level, choosing a 6–31 G basis set for Fe, F, and Cl atoms. The Gd was treated on a SARC-ZORA basis [35]. For computer experiments, we chose metal ions with half-filled d and f shells, i.e., d⁵ and f⁷ configurations, taking Fe³⁺ and Gd³⁺ in idealized molecular compounds with fluoride and chloride anions. The quantum chemical calculations were performed at the restricted open-shell Hartree–Fock level; although, for improved convergence, we worked in Complete Active Space Self Consistent Field (CASSCF) mode, within the trivial situation of a single high-spin state. We did not consider DFT, aiming to maintain compatibility with tentative further analyses of systems with other dⁿ or fⁿ configurations, implying degenerate or quasi-degenerate ground states and a multi-configuration WFT treatment. The mathematical and graphical handling was performed in Matlab-Octave [36,37] and Mathematica^TM environments [38,39]. The Mathematica^TM computer code is presented in the Supplementary Materials.

3. Results and Discussion

3.1. The Holohedrization in Standard Ligand Field Models

There are several varieties of ligand field models, some of which directly parameterize the content of the Hamiltonian matrix [40,41,42,43,44]. Here, we will be interested in the ligand field models proposing a Hamiltonian operator, expanded in spherical harmonics [45]. This principle meets in a fortunate way the fact that the basic ingredients of electronic structure methods are the atomic orbital wavefunctions, which are factorized into a radial part (with respect to the distance from the atomic nucleus) and an angular part, consisting of spherical harmonics. Another basic LF principle is that it addresses a limited basis, which is usually made of d-type or f-type atomic orbitals. Then, the basis on which the LF operator obtains its matrix representation is also made of spherical harmonics, Y_l,m(θ, φ), with the l = 2 or l = 3 quantum numbers. The radial part is absorbed in the parameters of the model so that the LF handles only the polar coordinates:

$${\widehat{V}}_{LF}\left(\theta ,\phi \right)={\sum }_{{k=k}_{min} }^{k_{max}} {\sum }_{q=-k }^k{B}_q^k Z_{k,q}\left(\theta ,\phi \right) ,$$

(1)

where $B_q^k$ values are adjustable to a specific problem and Z_k,q are combinations of spherical harmonics, conventionally defined as follows [46]:

$$Z_{k,0}\left(\theta ,\phi \right)=\sqrt{\frac{4\pi }{2k+1}}{Y}_{k,0}\left(\theta ,\phi \right),$$

(2a)

$$Z_{k,q}\left(\theta ,\phi \right)=\sqrt{\frac{4\pi }{2k+1}}\left(Y_{k,q}\left(\theta ,\phi \right)+{\left(-1\right)}^q{Y}_{k,-q}\left(\theta ,\phi \right)\right),$$

(2b)

$$Z_{k,-q}\left(\theta ,\phi \right)=\sqrt{\frac{4\pi }{2k+1}}i\left(Y_{k,-q}\left(\theta ,\phi \right)-{\left(-1\right)}^q{Y}_{k,q}\left(\theta ,\phi \right)\right),$$

(2c)

where q ≥ 0. The above expressions are meant to lead to real spherical harmonics, converting the exp(−iqφ) factor in sine and cosine, respectively. In standard ligand field models, there is a limitation on the k indices in the above summation. Namely, confined to the goal of using the above operators on a 2l + 1-dimensional basis consisting of the Y_l,m set (with the integer m running from −l to l), the ligand field operator is confined to k components, from 2 to 2l. More specifically, LF is applied either on d orbitals (l = 2, with k running only on k = 2 and k = 4) or the f shell (l = 3, with k = 2, 4, and 6).

The filtering of k indices occurs because of selection rules at the integration of the Y^*_l,m·Y_k,q·Y_l,m′ triads over angular coordinates m and m′ forming the line and column entries in the Hamiltonian matrix of the ${\widehat{V}}_k$ operator. Such integrals are resolved with the help of Clebsch–Gordan coefficients or the so-called Wigner’s 3j symbols [47]. Thee k = 0 term is discarded in the LF treatments because it gives a constant shifting of all the eigenvalues. Eliminating the k = 0 component leads to eigenvalues having their average fixed at zero.

In more detail, the construction of the ligand field Hamiltonian by spherical harmonics comes from the initial idea that it describes Coulomb interactions between one electron at the metal ion and the charge distribution in the environment, equated via the multipolar expansion of the electronic cloud on the surrounding atoms [3,4].

Observing that only spherical functions with even k are entering the Hamiltonian, we arrive at the title problem of holohedrization [40]. The inversion symmetry refers to the change in a function when the Cartesian coordinates (x, y, z) are switched to (−x, −y, −z). On a sphere of unitary radius, or at the transformation to solid spherical harmonics, when performing the conversion to Cartesian coordinates, r^lY_l,m(θ, φ) → C_l,m(x, y, z) ≡ Y_l,m(x, y, z), the sets can be labeled according to their behavior at inversion. Thus, the odd l indices show odd parity Y_l,m(−x, −y, −z) = −Y_l,m(x, y, z), while the even l functions are invariant at inversion, Y_l,m(−x, −y, −z) = Y_l,m(x, y, z). Then, the Hamiltonian in Equation (2) is invariant at inversion, irrespective of the set of $B_q^k$ quantities. Since the $B_q^k$ coefficients are enciphering the nature of the local geometry around the described metal ion (the so-called coordination polyhedron) [1,2], there is no way to account for the lack of an inversion center. The Hamiltonian cannot discriminate whether a perturbation exerted from a vertex comes from its place, (x, y, z) or its antipode, (−x, −y, −z). This means that the LF model artificially doubles a polyhedron by its inversion image. In other words, LF may describe the systems at a higher symmetry than their actual geometry.

The ${\widehat{V}}_{LF}\left(\theta ,\phi \right)$ operator can be drawn as a colormap on the surface of the sphere, as a function of angular coordinates. Although quite obvious, this sort of illustration is not used in textbooks or research approaches of LF problems. Figure 1 offers a picturesque demonstration of holohedrization. Namely, the symmetry of the LF potential for a half of octahedron, {MX₃}, provide that the XMX angles are kept at 90° (ortho-axiality), is the same as a rigorous octahedron {MX₆}. On the left side in Figure 1, one observes six equivalent maxima (represented in red) under the ligands. The blue areas correspond to minima between ligands. For the ortho-axial {MX₃}, LF yields six equivalent areas with maxima. Because holohedrization is unable to distinguish between opposed directions on a line, the perturbation exerted by a ligand on an M-X axis is smeared as if produced by half of a ligand on both ends, (X/2)-M-(X/2). Then, the above-described {MX₃} pyramid can be formalized as an {M(X/2)₆} octahedron. Dividing the perturbation power of three ligands on six positions, the overall potential becomes shallower, as suggested by the paler coloration of the map on the right-side panel compared with the left one.

Holohedrization appears when the LF Hamiltonian is tailored for a basis made of a single type of pure atomic orbitals, either d or f. Thus, in the matrix representation, the standard LF operator implies either $⟨d_i\left|{\widehat{V}}_{LF}\right|d_j⟩$ or $⟨f_i\left|{\widehat{V}}_{LF}\right|f_j⟩$ elements. To have non-vanishing matrix elements, the symmetry of the operator must match those of d × d or f × f product representations. The d orbitals show even parity at inversion (labeled by subscript g), while the f functions undergo a sign change (labeled by u). Despite the different parity of d versus f shells, their products behave with g-type parity (because g = g × g and g = u × u), demanding this symmetry for the operator itself. Thus, holohedrization appears from the empirical premise of having as a basis pure d or f orbitals. At the same time, in the spirit of phenomenological approaches, it is tacitly hoped that the Hamiltonian will include ligand contributions; the occurrence of holohedrization is largely overlooked despite early recognition of this drawback [40]. Working within a non-empirical electronic structure theory, where the requirement of pure d or f basis is not imposed, does not imply holohedrization as a necessary outcome. To check this issue, we proposed the case studies debated in the following.

The octahedral LF Hamiltonian, tailored for a d-type orbital basis, i.e., limiting the Equation (1) to k_max = 4, is:

$${\widehat{V}}_{LF}^d[\mathrm{M}{\mathrm{X}}_6{]}_{O_h}^{}=14\sqrt{\pi }Dq\left(Y_{\mathrm{4,0}}\left(\theta ,\phi \right)+\sqrt{\frac{5}{14}}\left(Y_{\mathrm{4,4}}\left(\theta ,\phi \right)+Y_{4,-4}\left(\theta ,\phi \right)\right)\right).$$

(3)

Using symmetry, many parameters from Equation (1) are vanishing or are simplified to a common factor, traditionally denoted as Dq. Expanding the spherical harmonics in their trigonometric functions, one obtains:

$${\widehat{V}}_{LF}^d[\mathrm{M}{\mathrm{X}}_6{]}_{O_h}^{}=\frac{21}{8}Dq\left(3-30Cos(\theta {)}^2+35Cos(\theta {)}^4+5Cos\left(4\varphi \right)Sin(\theta {)}^4\right),$$

(4)

observing that this function is invariant at changing the coordinates to their antipodes, namely θ → 2π − θ, φ → φ + π. Because of holohedrization, the LF Hamiltonian of the ortho-axial {MX₃} is half of the octahedral one:

$${\widehat{V}}_{LF}^d[\mathrm{M}{\mathrm{X}}_3{]}_{\mathrm{X}\mathrm{M}\mathrm{X}=90°}^{}={\frac{1}{2}\widehat{V}}_{LF}^d[\mathrm{M}{\mathrm{X}}_6{]}_{O_h}^{}$$

(5)

Figure 1 was obtained by Equations (4) and (5) for the left- and right-side panels, respectively, taking Dq = 1000 cm⁻¹. In the following, we are going to explore, however, whether more realistic representations of the effective field from ligands are still undergoing a holohedrization effect.

The ligand field operator has the meaning of an effective field perturbing the electrons located in the d or f orbitals of the metal ion. The true electric field exerted by the charged ligand is not expected to show any holohedrization. Figure 2 shows the color maps of six fluoride anions forming an octahedron compared with the tri-fluoride set, which was obtained by keeping the ions placed on the positive sides of the Cartesian axes. The octahedron shows the same qualitative pattern as those shown in Figure 1 for the effective operator in Equations (3) or (4). The ortho-axial pyramid shows a polarized map, convening to set the positive zone in the quadrant marked by the three ligands.

The eigenvalues of the holohedrized {MX₃} Hamiltonian model for d electrons, three times −2Dq and two times 3Dq, are also half of the octahedron, −4Dq and 6Dq, respectively, assuming that the perturbation power of the X ligands, proportional with the Dq quantity, is transferable between coordination units. Holohedrization is reflected in higher symmetry of the eigenvalues pattern in the LF models which do not rely on the Hamiltonians in Equation (1). For instance, in the Angular Overlap Model (AOM) [41,43,44,45], the eigenvalues for the {MX₆} octahedron are 3e_σ^X for the doubly degenerate e_g orbital set and 4e_π^X for the triply degenerate t_2g functions. The AOM parameters are classified according to the bonding capability (σ,π) and ligand type (X). In the AOM, an ortho-axial {MX₃} obtains the (3/2)e_σ^X and 2e_π^X eigenvalues with double and triple orbital degeneracies, i.e., experiencing artificial higher symmetry, by holohedrization. Obtaining an octahedral pattern (O_h) for a case that should obey trigonal symmetry (C_3v point group, for the pyramid) is clearly a drawback of the model.

3.2. The Partial Holohedrization Effect in Orbital Energies from Ab Initio Quantum Calculations

Let us consider some calculation examples that are idealized yet realistic. Although more advanced options are possible, we will take as the ligand field sequence the frontier orbitals with predominant metal ion character, which are occupied by unpaired electrons, in the chosen computation samples. The [FeF₆]³⁻ unit has the following relative orbital energies in the LF set: {0, 0, 0, 9788.57, 9788.57}, measured in cm⁻¹. One observes the well-known pattern that is specific for the octahedron, with a triple degenerate t_2g as the lowest level and a higher e_g couple. The gap, 10Dq = 9788.57, falls in the expected range for d-type compounds [5]. The neutral [FeF₃] (enforced at the same Fe-F bond lengths, 1.984Å, like in octahedron and right F-M-F angles), yields the following LF sequence: {0, 899.85, 899.85, 7747.45, 7747.45} cm⁻¹. This system has C_3v symmetry, with d-type orbitals spanning the a₁ + 2e representations, recognizing in the above result the two doubly degenerate e representations. Since the first three energy values (corresponding to a₁ + e representation in C_3v) from the above list are mutually close, they can be roughly assimilated to a triple degenerate t_2g set in a quasi-octahedral pattern. Then, the energy scheme can be taken as indicating a holohedrization trend. The Supplementary Materials show the molecular orbitals of the [FeF₃] unit, illustrating the preponderance of d-type orbitals, aside from small ligand tails. This is a prerequisite of manifesting holohedrization, which would be rigorously true if the active orbitals were pure d atomic functions.

Similar conclusions are drawn by examining the case of chloride systems. The [FeCl₆]⁶⁻ (with Fe-Cl = 2.556 Å bond length) shows the {0, 0, 0, 6759.82, 6759.82} cm⁻¹ LF energies, while the [FeCl₃] one obtains {0, 855.96, 855.96, 6386.72, 6386.72} cm⁻¹. Here, one may see that by taking the average energy of a₁ + e as a pseudo-octahedral t_2g, one obtains a gap (5530.8 cm⁻¹) that is quite far from the half of 10Dq in the rigorous octahedron (6759.82 cm⁻¹). This does not literally obey the holohedrization constraint. On the other hand, the chlorine ligand contributions in the molecular orbitals of the ligand field sequence are relatively small, pointing to effective holohedrized behavior.

An interesting situation results if we consider lanthanide-based octahedrons and pyramids. The [GdF₆]³⁻ octahedral unit (with 2.556 Å bond length) has a {0, 263.37, 263.37, 263.37, 724.27, 724.27, 724.27} set of energies (in cm⁻¹) corresponding to the split of f orbitals in the {a_2u, t_2u, t_1u} octahedral representations. Note that the order of magnitude in the total gap is much smaller than the transition metal example, a situation that is characteristic for lanthanides [6] because the f shell radial profile has a maximum at a point shrunk below the ionic radius and the f-type electrons undergo only a small perturbation from the environment. The [GdF₃] hypothetical ortho-axial unit produces the following LF type energies: {0.//109.74, 109.74, 263.37//351.16, 351.16, 504.79} (in cm⁻¹). The “//” delimiters in the above list suggest the quasi-degenerate sequences assignable to the {a_2u, t_2u, t_1u} quasi-octahedral pattern. The active orbitals of the [GdF₃] unit, depicted in the Supplementary Materials, show negligible ligand contributions, a fact determining holohedrization.

3.3. Holohedrization in the Ab Initio Results

In the following, we will design an extended LF-like Hamiltonian using the orbital eigenvalues and eigenvectors from ab initio calculations. Knowing the ε_i eigenvalues and the ψ_i ortho-normal eigenvectors of a quantum object, the most general formalism [7] to design an effective Hamiltonian is the following expansion:

$${\widehat{H}}_{eff}={\sum }_{\mathrm{i}}\left|{\psi }_i\right.⟩{\epsilon }_i⟨{\psi }_i|.$$

(6)

We used the $⟨bra|$ and $\left|ket\right.⟩$ notations for the functions designated to form the left and right side components of the matrix elements of the operator, $⟨bra|\widehat{H}|ket⟩$, which are integrals over the 3D space coordinates. Because of ortho-normalization, we have $⟨{\psi }_i|\widehat{H}|{\psi }_i⟩={\epsilon }_i$, and the effective operator is obtained by the reverse engineering of the available solutions. Here, we will consider as eigenvalues the orbital energies discussed in the previous section. Their related eigenvectors, the ψ_i molecular orbitals, are obtained as a linear expansion over atomic components, which will be written separately for metal centers and ligands by μ and λ, respectively. The first ones are placed in origin and the others are at {X_L, Y_L, Z_L,} Cartesian coordinates of the atoms L forming the ligands:

$${\psi }_i(x,y,z)={\sum }_{M\in metal ion}{c}_{iM}\cdot {\mu }_M(x,y,z)+{\sum }_{L\in ligands}{c}_{iL}\cdot {\lambda }_L(x-X_L,y-Y_L,z-Z_L).$$

(7)

In the performed calculations, as in the vast majority of computational chemistry approaches, the atomic components consist of the so-called Gaussian-type orbitals (GTOs) [48,49]. The GTOs are handled in the format of the local Cartesian coordinates (x, y, z) of each atom, but they can also be formally expressed in the equivalent polar coordinates (r, θ, φ). Then, one may detail the dichotomy of the atomic components (with χ notation standing for both above μ and λ) in the radial and angular part:

$$\chi \left(x,y,z\right)\equiv \chi \left(r,\theta ,\phi \right)=R\left(l,\zeta ,r\right)\cdot Z_{lm}\left(\theta ,\phi \right),$$

(8)

the notation Z_lm represents, like in Equation (2), the conversion of Y_lm complex spherical harmonics to real forms. The GTO radial part is:

$$R\left(l,\zeta ,r\right)=\mathcal{N}\left(l,\zeta \right)\cdot r^l\mathit{exp}\left(-\zeta \cdot r^2\right),$$

(9a)

$$\mathcal{N}(l,\zeta )=2^{\frac{2l+5}{4}}{\zeta }^{\frac{2l+3}{4}}/\sqrt{\Gamma (l+3/2)} ,$$

(9b)

where Γ is the well-known generalization of the factorial (gamma function) [50]. The set of ζ parameters forms the definition of the selected basis sets, while $c_{iM}$ and $c_{iL}$ are the result of the calculation (although some pre-determined factors are also part of basis set formats).

To make an LF-like operator, one must integrate the ψ_i molecular orbitals in Equation (6) on radial coordinate, keeping the origin on the metal ion. As a first approximation, one may consider only the metal-based atomic orbitals since these are usually the main part of LF-type molecular functions. In this case, the integration can be performed analytically, obtaining a generalized ligand field operator:

$$V_{LF}^M(\theta ,\phi )={\sum }_{M_1}{\sum }_{M_2}{\epsilon }_i c_{i{M}_1}{c}_{i{M}_2}\cdot S_{M_1,M_2}{\cdot Z}_{l_{M_1}{m}_{M_1}}^{}(\theta ,\phi )Z_{l_{M_2}{m}_{M_2}}^{}(\theta ,\phi ),$$

(10)

where:

$$S_{M_1,M_2}\equiv S\left(l_{M_1},{\zeta }_{M_1},l_{M_2},{\zeta }_{M_2}\right)={\int }_{r=0}^{\infty }R\left(l_{M_1},{\zeta }_{M_1},r\right)\cdot R\left(l_{M_2},{\zeta }_{M_2},r\right)r^{2}dr,$$

(11)

and:

$$\left(l_1,{\zeta }_1,l_2,{\zeta }_2\right)=\frac{1}{2}\mathcal{N}\left(l_1,{\zeta }_1\right)\mathcal{N}\left(l_2,{\zeta }_2\right)\cdot \Gamma \left(\frac{l_1+l_2+3}{2}\right)/{({\zeta }_1+{\zeta }_2)}^{\frac{l_1+l_2+3}{2}}.$$

(12)

The products of two spherical harmonics with l₁ and l₂ quantum numbers can be expanded as a sum of spherical functions running between |l₁ − l₂| and l₁ + l₂. For complex standard spherical functions, the expansion is performed with well-defined coefficients, either by Wigner’s 3j symbols or Clebsch–Gordan factors [47]. For real functions, as formulated in Equation (10), this conversion can be adapted, leading to an expansion factored in Z_lm spherical harmonics.

In principle, one may also realize the analytic expansion for the full molecular orbitals, with ligand parts included. However, for practical reasons, we turned to numerical integration. In this view, we borrowed a technique employed in the making of pseudo-potentials used in solid-state calculations [51]. The point is to express a radial integral as a sum of the integrand over a fixed set of points, factored by their specific weights.

$${\int }_{r=0}^{\mathrm{\infty }}f\left(r\right)dr\approx {\sum }_{m=1}^{n_{max}}{w}_m\cdot f\left(r_m\right).$$

(13)

We took a radial grid [52,53] that is denser at the origin and sparser at long range:

$$r_m=\delta r_0\frac{\mathit{exp}\left(mh\right)-1}{\mathit{exp}\left(h\right)-1},$$

(14a)

$$w_m=h\cdot \delta r_0\frac{\mathit{exp}\left(kh\right)}{\mathit{exp}\left(h\right)-1}.$$

(14b)

The maximal extension, r_max, and the number of the points, n_max, are tuned by δr₀ and h parameters, as follows:

$$n_{max}=int\left[\frac{1}{h}\mathit{ln}\left(1+r_{max}\frac{\mathit{exp}\left(h\right)-1}{\delta r_0}\right)\right],$$

(15)

where int means taking the integer part. Here, we used δr₀ = 0.001 Bohr and h = 0.0211, leading to n_max = 254 and r_max = 10 Bohr. The functions subjected to numerical integration in Equation (13) are the extensions of Formulas (10) and (11) that are used for metal–ligand terms:

$$f\left(r\right)=R\left(l_M,{\zeta }_M,r\right)\cdot R\left(l_L,{\zeta }_L,|\mathit{r}-{\mathit{R}}_{\mathit{L}}|\right)\cdot Y_{LM}^{real}({\theta }_L,\phi ),$$

(16)

where R_L marks the position of the ligand atom L and the bold writing suggests Cartesian vector triads. Because of representation with respect to the L centers, the spherical harmonics of ligands acquire a dependence on the radial part. Actually, only the θ variable has to be shifted to the L center by the r·cos(θ) − |r − r_L|·cos(θ_L) = R_ML relationship (where R_ML is the metal–ligand atom distance) because φ is the same for the M and L sites. The ligand–ligand analogs of the terms in Equations (1) and (6) can be neglected because of their small coefficients and practically null overlap of their radial functions. Although the metal-only part can be performed analytically in Equations (10)–(12), we also worked it numerically for comparability with the extended approach, as sketched in Equations (13)–(16).

The radial integrations from 0 to r_max = 10 Bohr were performed on a mesh of polar coordinates with N_θ = 25 ticks for θ and N_φ = 48 for φ. Then, the set of {θ_i = (i − 1)π/(N_θ − 1), φ_j = 2(j − 1)π/(N_φ − 1)} grid points, containing the v_ij results of the radial integration in the different directions, is fitted by a generalization of Equation (1), allowing the spherical harmonics to run on a large range, from k_min = 1 to k_max = 10. We discarded the k = 0 function, giving only the average of the potential on the sphere, i.e., not leading to LF splitting.

In Figure 3, one may see the results of the analysis for the [FeF₃] and [GdF₃] case studies. We see, in the left-side panels, that the fingerprint of holohedrization is retained, with local maxima in the direction of metal–ligand axes and at their antipodes. As a fine detail, one may observe that in the [FeF₃] case, the potential below the ligand is slightly different than in the opposed part. For [GdF₃], one may note that the circular profiles corresponding to the perturbation are a bit off-centered with respect to the metal–ligand axes. Taking the difference between the complete integration of the metal + ligand effects and those limited to the metal ones, a non-holohedral local perturbation oriented toward the ligands is revealed. It turns out then that the polarized behavior may appear entirely due to ligands, but this contribution is, at least in the considered cases, relatively minor. For the [GdF₃] system, the difference map is slightly asymmetric and hearth-like shaped because of the observed deviation of the effective potential from the Gd-F axes. This may be because the f-type orbitals are contributing insignificantly to the chemical bonding [54]. Then, the ligands interact, in fact, with other atomic orbitals on lanthanide, the 6s, 6p, and 5d functions. An s-p or p-d hybridization on the metal may create a slight anisotropy, modulating the ligand contributions in the octant containing the {F₃} octahedral face versus the neighboring areas. Conversely, the d-type orbitals are directly involved in the chemical bonding of the system, maintaining the quasi-axial symmetry of the metal–ligand effects.

3.4. Constructing a Generalized Ligand Field Hamiltonian

The above results come from the direct representation of the v_ij data on the {θ_i, φ_j} grid of points. In the following, we will complete the treatment, fitting the numeric results into an extended series of spherical harmonics in Equation (1): the k_min = 1 and k_max = 10 limits. Taking a large spectrum of spherical functions, we arrive at an overdetermined linear system of equations:

$$\mathbf{𝒜}\cdot \mathit{b}=\mathit{u},$$

(17)

where the coefficients to be fitted are the generalized LF parameters, $B_k^q$, which are concatenated in one column for successive k sets, skipping the trivial k = 0:

$$b_{k^2+k+q }\equiv b_{k^2+k+q,1 }=B_k^q.$$

(18)

The data feeding the problem are the v_ij values of the potential estimated at the {θ_i, φ_j} mesh of polar coordinates, which are counted by merging the successive lines of the v_ij matrix in a single column:

$$u_{N_{\phi }\cdot \left(i-1\right)+j}\equiv u_{N_{\phi }\cdot \left(i-1\right)+j,1}=v_{ij}.$$

(19)

The structure factors are obtained by estimating the whole set of real spherical harmonics at the grid of polar points, which are arranged in the following matrix:

$${𝒜}_{{N_{\phi }\cdot \left(i-1\right)+j, k}^2+k+q }=Z_{kq}({\theta }_i,{\phi }_j).$$

(20)

In the discussed setting, there are 1200 v_ij points and 120 generalized LF coefficients $B_k^q$ of spherical harmonics, running from k_min = 1 to k_max = 10. The solution, in the least square sense, is obtained by pseudo-inverse:

$$\mathit{b}={\mathbf{𝒜}}^{\#}\cdot \mathit{u},$$

(21a)

$${\mathbf{𝒜}}^{\#}={\left({\mathbf{𝒜}}^{\mathit{T}}\cdot \mathbf{𝒜}\right)}^{-1}\cdot {\mathbf{𝒜}}^{\mathit{T}}.$$

(21b)

This analysis can be applied to the full computed potential (the v_ij points) or the model confined to metal-only contributions (say, v_ij^M points). Additionally, one may attempt the fit on the limited sets of spherical harmonics used in the traditional ligand field, namely k = 2 and 4 for d-type elements (the [FeF₃] system) or k = 2, 4, and 6 for f-type ones (the GdF₃ case). Figure 4 shows the concatenated lists of the $B_k^q$ coefficients with k from 1 to 10, and, inside each k domain, q ranges from −k to k. As equated in the above linear system, the $B_k^q$ parameter occupies the k² + k + q position.

Because the coefficients at large k are small (or null, by design), it is convenient to represent the results on a logarithmic scale, by the {log₁₀(k² + k + q), $B_k^q$} couples. In this way, the domains of large k indices are compressed at the right-side margin. The overview in Figure 4 shows a relative similarity between the different fitting choices. The poorer setting is for the graphics shown at the bottom line, with spherical harmonics limited to the traditional LF pattern. Thus, one may note the null parameters at k = 1, 3, and for k > 4 in the d-type case (left side of Figure 4). The f-type system yields null $B_k^q$ for k = 1, 3, 5, and k > 6. Going to the approximation of metal-only based ab initio fit, one notes the quenching of higher terms, with k > 4 in the Fe system and k > 6 in the Gd one. This modeling option enables the intermediate contributions from k = 1 and 3 (left side) and k = 1, 3, and 5 (right side) terms. However, the related $B_k^q$ magnitudes are small. The top graphs, with the fully extended model, show the main contributions in the same pattern as the previously limited fit attempts. This suggests that, despite intrinsic limitations, the traditional LF paradigm may work well, at least in the case of systems with a pronounced ionic character, as the chosen systems are. We collaterally reached a similar perspective using a different methodology and other case studies [55]. At the same time, we tested that the results on chlorine analogs, which were expected to have a lower ionic nature, behave in a similar way to the discussed fluorine cases. The situation may change if we take neutral polyatomic ligands, but this is a matter for further developments.

4. Conclusions

This article is a detective search that questions the issue of holohedrization, which affects the standard ligand field phenomenological models. It means that the artificial symmetry increases when treating a metal ion in a surrounding environment without an inversion center; the model behaves as if the system has one. In other words, a low-symmetry coordination polyhedron effectively acts in the ligand field approach as if it is doubled by inversion. This issue was early recognized but received very little attention for decades and is practically forgotten nowadays. Although ligand field models are semiempirical, in a practical sense, can help, at a conceptual level, decrypt physical meaning from the black box of modern quantum calculations on molecules and solid-state crystals. Aiming to see if the effect persists at the extension of the ligand field effective operator with the help of state-of-the-art ab initio calculations, we embarked on a quest for an interdisciplinary insight. For the sake of systematic progress, we selected simple case studies. The quantum chemical calculations on idealized (yet realistic) systems were analyzed in a non-routine manner, handling the output data with advanced mathematical tools, analytically and numerically. We found that holohedrization also appears in extended modeling, at least in ionic systems, when starting from orbital energies and computed molecular orbitals as objects assignable to the ligand field phenomenology. The inquiry remains open to debate, with us aiming to test, in the future, systems with less ionic character and a multi-configurational approach to the debated problem.