Rigorous Formalization of Orbital Functionals: Addressing the Noninteracting $v$-Representability Problem
Neil Qiang Su




[Corresponding author. ]nqsu@nankai.edu.cn
Frontiers Science Center for New Organic Matter, State Key Laboratory of Advanced Chemical Power Sources, Key Laboratory of Advanced Energy Materials Chemistry (Ministry of Education), Department of Chemistry, Nankai University, Tianjin 300071, China





Functionals that explicitly depend on occupied, unoccupied, or fractionally-occupied orbitals are rigorously formalized using Clifford algebras, and a variational principle is established that facilitates orbital (and occupation) optimization as a formal implementation method. Theoretically, these methodologies circumvent the limitations encountered in the original Kohn-Sham and related methods, particularly when the interacting system's electron density does not match that of any noninteracting reference system. This work redefines orbital (and occupation) functionals from a novel perspective, positioning them not merely as extensions of traditional density functionals, but as superior, rigorous alternatives.

Rigorous Formalization of Orbital Functionals: Addressing the Noninteracting v-Representability Problem
    Neil Qiang Su
    April 29, 2024
=======================================================================================================




The Hohenberg-Kohn (HK) theorems fundamentally established the basis for the original density functional theory (DFT) by demonstrating the one-to-one correspondence between the external potential and electron density for nondegenerate ground states, and by establishing a variational principle to compute the ground-state energy directly from densities <cit.>. However, the HK theorems and their extensions to ensemble degenerate states <cit.> do not resolve all challenges associated with the implementation and application of functional theories for many-electron systems. Notably, issues related to density v-representability are crucial—a density is deemed v-representable if it is a ground-state density for some external potential<cit.>.

Overall, v-representability issues stem from two critical perspectives: the Hohenberg-Kohn (HK) variational principle<cit.> and the Kohn-Sham (KS) method<cit.>. The HK variational principle restricts densities to those of the interacting system under some external potential, thereby introducing the interacting v-representable (IVR) condition<cit.>. However, the absence of explicit constraints for IVR densities hampers the variational process, necessitating the inclusion of non-IVR trial densities. Levy addressed this issue through the constrained search formulation <cit.>, proposing a universal functional over easily realizable N-representable (NR) densities—those derivable from any antisymmetric wavefunction. Further efforts on this v-representability issue have been made <cit.>, with notable contributions from Yang, Ayers, and Wu, who approached the issue through potential functional theory <cit.>.

Another v-representability issue emerges in the KS method, which presupposes the existence of a noninteracting reference system governed by an effective potential such that its density matches the ground-state density of the interacting system<cit.>. This prerequisite, known as the KS v-representable (KSVR) assumption, is essential for the original KS method and related methodologies <cit.>, facilitating the simplification of the variational determination of density and energy. However, the KSVR assumption reveals its limitations when the interacting system's ground-state density is non-KSVR, potentially corresponding to an excited state or failing to match any stationary state of a noninteracting system <cit.>. This KSVR problem thus compromises the reliability of computational results. Moving beyond traditional functionals like the generalized gradient approximations (GGAs) <cit.> and hybrid functionals <cit.>, emerging functionals that explicitly depend on occupied, unoccupied, or fractionally-occupied orbitals demonstrate superior performance <cit.>. Although it remains feasible to adhere to the KSVR assumption and use the optimized effective potential (OEP) method <cit.>, to obtain density and orbitals for these new functionals, the KSVR problem persists. Beyond the KS method, new implementations such as orbital optimization <cit.> (or equivalently, generalized OEP (GOEP) <cit.>), alongside combined orbital and occupation optimization <cit.>, have emerged. However, the theoretical foundations of these new functionals and implementations require further demonstration. This work aims to provide a solution to the KSVR problem and establish a rigorous theoretical basis for various orbital functionals and orbital (and occupation) optimization implementations.

For an N-electron interacting system subject to an external potential v (Ĥ_v=T̂ + V̂_ee+∑_i^N v(r_i)), the Levy constrained search formulation determines the ground-state energy through the following variational principle <cit.>:

    E_v = min_ρ E_v^ℒ[ρ],

where the trial density ρ can be any NR density, and E_v^ℒ[ρ] is the energy functional given by

    E_v^ℒ[ρ] = min_Ψ→ρ⟨Ψ | Ĥ  | Ψ⟩,

in which the minimization searches all normalized antisymmetric wavefunctions that yield ρ to identify the minimizing wavefunction Ψ_ρ^min. Upon successful minimization of Eqs. <ref> and <ref>, the resultant energy E_v and the corresponding ground-state density and wavefunction, denoted ρ_v and Ψ_v, respectively, are obtained.


Indeed, Eqs. <ref> and <ref> correspond respectively to the implementation method and the energy functional required. In practice, an effective strategy is essential to search reasonable NR densities for approaching ρ_v. The KS method restricts densities to those corresponding to a noninteracting reference system with an effective potential w_s, i.e., the KSVR densities ρ_w_s <cit.>. This method not only simplifies the search but also leads to a set of orbitals, known as KS orbitals, which are used to generate ρ_w_s and to construct both the kinetic energy functional and the exchange-correlation functional in Eq. <ref> <cit.>. However, since ρ_v is not necessarily a KSVR density, thus potentially rendering the KS implementation ineffective.


This work abandons the traditional use of a noninteracting reference system to address the limitations associated with KSVR densities. Subsequent discussions transition from a generalized Slater determinant to focusing on orbital functionals. Unlike the ground-state wavefunction of noninteracting systems, which is composed of occupied KS orbitals, the determinant in this study is constructed from an arbitrary set of orthonormal one-electron functions. To explore a more comprehensive scenario, hypercomplex orbitals are employed <cit.>:

    φ_n,i(𝐫) = ϕ_i^0(𝐫) + ∑_μ=1^nϕ_i^μ(𝐫) e_μ,

with {ϕ_i^μ} as real functions and {e_1, …, e_n} as a basis in a Clifford algebra of dimension n <cit.>, satisfying:

    e_μ^2 = -1;    e_μ e_ν = -e_ν e_μ.

The determinant constructed from {φ_n,i} is denoted Φ_n.  At n=0, {φ_n,i} and Φ_n revert to traditional real orbitals and determinant. With this construction, the density for Eq. <ref> can be obtained via:

    ⟨Φ_n | ρ̂ | Φ_n ⟩ = 2∑_i |φ_n,i(𝐫)|^2 =2 ∑_p λ_n,p |χ_n,p(𝐫)|^2,

where the last equality was derived in Ref. , and {χ_n,p, λ_n,p} represent hierarchically correlated orbitals (HCOs) and their occupations <cit.>. As indicated by Eq. <ref>, the density remains invariant under the unitary rotation of {φ_n,i} and of {χ_n,p} with the same occupation.



Theorem 1: The mappings among Φ_n, {φ_n,i}, and {χ_n,p, λ_n,p} are surjective for a given n, as shown in:

    Φ_n ⇔{φ_n,i}⇔{χ_n,p, λ_n,p}.

Specifically, there exists a one-to-one mapping between Φ_n and {φ_n,i}, and a many-to-one mapping between {φ_n,i} and {χ_n,p, λ_n,p}, within the unitary rotation of {φ_n,i} and of {χ_n,p} with the same occupation. Refer to Fig. <ref> for visualization.





Proof 1: The proof of this theorem is detailed in the supplementary material <cit.>. □

Theorem 1 establishes that the densities derived from Φ_n, {φ_n,i}, and {χ_n,p, λ_n,p}—collectively referred to as ρ_n—share the same domain. The subsequent theorem will delve deeper into the properties of ρ_n and delineate their advantages over ρ_w_s.


Theorem 2: The domain of ρ_n expands as n increases, exhibiting a matryoshka-like structure. The relationships among the domains of ρ_w_s, ρ_v, and ρ_n from n=0 to n=∞ are shown as:

    domρ_w_s⊆domρ_n=0⊆domρ_n=1⋯⊆domρ_n=∞⊇domρ_v

This indicates that the domain of ρ_n fully encompasses KSVR densities and progressively includes all IVR densities, as depicted in Fig. <ref>.




Proof 2: Initially, consider the relationship between ρ_w_s and ρ_n=0. Since the ground state of a noninteracting reference system is typically represented by a Slater determinant, all ρ_w_s inherently reside within the domain of ρ_n=0. However, certain ρ_n=0 densities are identified as non-KSVR, either corresponding to excited states or not aligning with any stationary state of a non-interacting system, establishing the first relation in Eq. <ref>. This analysis assumes ρ_w_s is derived from conventional real KS orbitals. Should complex KS orbitals be used, ρ_n=0 transitions to ρ_n=1, with no additional modifications required.


The matryoshka-like structure of the domains of ρ_n for different n levels is elucidated by Eq. <ref>. As n increases, the degrees of freedom for {φ_n,i} increase, broadening the domain of ρ_n. Moreover, {φ_n,i} at lower n can be dimensionally reduced from those at higher n, ensuring that the domain of ρ_n at higher n fully encompasses that at lower n. Additionally, the constraints on {χ_n,p, λ_n,p} <cit.> illustrate another aspect: at n=0, {χ_n,p, λ_n,p} are subject to stringent constraints—{χ_n,p} are orthonormal and {λ_n,p} are integers, with N being 1 and the rest 0; as n increases, these constraints relax, thereby extending ρ_n's domain.


Regarding the domain of ρ_n=∞, as n approaches infinity, the constraints on {χ_n,p, λ_n,p} become minimal, described by <cit.>:

    ⟨χ_n,p | χ_n,q⟩ = δ_pq,


    0 ≤λ_n,p≤ 1,   ∑_pλ_n,p = N/2.

As detailed in Eq. <ref>, the density ρ_n in terms of {χ_n,p, λ_n,p} is:

    ρ_n(𝐫) = 2 ∑_p λ_n,p |χ_n,p(𝐫)|^2.

With these constraints, Eq. <ref> can reproduce any NR density, conclusively establishing that the domain of ρ_n=∞ encompasses that of ρ_v. □


The above two theorems confirm that Φ_n, {φ_n,i}, and {χ_n,p, λ_n,p} are equivalent for generating trial densities. The resulting densities, ρ_n, are notable for their scalability, extending their domain with increasing n to cover a broader range of plausible NR densities. Additionally, as derived from Ref. , under the ideal complete basis set limit, ρ_n=0 can produce any NR density, making the domain of ρ_n uniform across all n. Consequently, ρ_n are termed meta-NR densities, highlighting their ability to transcend traditional density search paradigms. Building on these insights, further discussions will elucidate how Φ_n, {φ_n,i}, and {χ_n,p, λ_n,p} serve as fundamental descriptors of the system.

Taking {χ_n,p, λ_n,p} as an example, the following variational principle is constructed:

    E_n,v = min_{χ_n,p, λ_n,p} E_n,v^pw[{χ_n,p, λ_n,p}],

with the energy functional defined as:

    E_n,v^pw[{χ_n,p, λ_n,p}] = min_Ψ→ρ_{χ_n,p, λ_n,p}⟨Ψ | Ĥ_v | Ψ⟩,

where the superscript pw stands for present work. Eq. <ref> seeks to minimize across the domain of {χ_n,p, λ_n,p}, with the density of each trial {χ_n,p, λ_n,p} given by ρ_{χ_n,p, λ_n,p} = 2 ∑_p λ_n,p |χ_n,p(𝐫)|^2.




Theorem 3: Assume ρ_v is a meta-NR density for a given n. Then, E_n,v^pw[{χ_n,p, λ_n,p}] is stationary with respect to variations in {χ_n,p, λ_n,p} when ρ_{χ_n,p, λ_n,p} = ρ_v, reaching its minimum and providing the ground state energy:

    E_v = E_n,v.


Proof 3: The variational principle in Eq. <ref> stipulates that when the trial density equals the ground state density (ρ = ρ_v), the Euler-Lagrange equation is satisfied <cit.>:

    δ E_v^ℒ[ρ]/δρ(r)|_ρ(r)=ρ_v(r) = μ,

where μ is the chemical potential ensuring the density is NR. For E_n,v^pw[{χ_n,p, λ_n,p}], consider small variations within {χ_n,p, λ_n,p}'s domain:

    δ E_n,v^pw[{χ_n,p, λ_n,p}]/δ{χ_n,p, λ_n,p} = ∫δ E_n,v^pw[{χ_n,p, λ_n,p}]/δρ_{χ_n,p, λ_n,p}(r)δρ_{χ_n,p, λ_n,p}(r)/δ{χ_n,p, λ_n,p} dr.

Upon comparing the definitions of E_v^ℒ[ρ] and E_n,v^pw[{χ_n,p, λ_n,p}] (namely Eqs. <ref> and <ref>), when ρ_{χ_n,p, λ_n,p} = ρ_v, the above equation simplifies to:

    δ E_n,v^pw[{χ_n,p, λ_n,p}]/δ{χ_n,p, λ_n,p} = μ∫δρ_{χ_n,p, λ_n,p}(r)/δ{χ_n,p, λ_n,p} dr = 0,

confirming that E_n,v^pw[{χ_n,p, λ_n,p}] is stationary. Furthermore, the following can be obtained:

    E   _n,v^pw[{χ_n,p, λ_n,p}]|_ρ_{χ_n,p, λ_n,p}≠ρ_v = min_Ψ→ρ_{χ_n,p, λ_n,p}⟨Ψ | Ĥ_v | Ψ⟩
       ≥min_Ψ→ρ_v⟨Ψ | Ĥ_v | Ψ⟩ = E_n,v^pw[{χ_n,p, λ_n,p}]|_ρ_{χ_n,p, λ_n,p} = ρ_v,

indicating that E_n,v^pw[{χ_n,p, λ_n,p}] reaches the minimum E_n,v when ρ_{χ_n,p, λ_n,p} = ρ_v and equals E_v. □


Theorem 3 thus demonstrates that {χ_n,p, λ_n,p} fundamentally defines the system, with the ground state energy obtainable by minimizing the functional E_n,v^pw[{χ_n,p, λ_n,p}]. For simplicity, E_n,v^pw[{χ_n,p, λ_n,p}] is decomposed as:

    E   _n,v^pw[{χ_n,p, λ_n,p}] = T_n[{χ_n,p, λ_n,p}] + J[ρ_{χ_n,p, λ_n,p}] 
       + E_n,xc[{χ_n,p, λ_n,p}] + ∫ v(r) ρ_{χ_n,p, λ_n,p}(r)   dr,

where T_n represents the kinetic energy <cit.>,

    T_n[{χ_n,p, λ_n,p}] = - ∑_pλ_n,p⟨χ_n,p | ∇^2 | χ_n,p⟩.

J denotes the Coulomb energy, E_n,xc is the exchange-correlation energy accounting for the remaining kinetic and electron interaction energies not covered by T_n and J, and the final term represents the external potential energy.



The above derivations and conclusions require only minor modifications to apply to Φ_n and {φ_n,i}. Specifically, in Eqs. <ref>-<ref>, replace {χ_n,p, λ_n,p} with Φ_n and {φ_n,i} respectively, substituting ρ_{χ_n,p, λ_n,p} with ρ_Φ_n = ⟨Φ_n | ρ̂ | Φ_n ⟩ and ρ_{φ_n,i} =2 ∑_i |φ_n,i(𝐫)|^2. The energy functionals  E_n,v^pw[Φ_n] and E_n,v^pw[{φ_n,i}] follow the decomposition in Eq. <ref>, with T_n of Eq. <ref>, expressed respectively in terms of Φ_n and {φ_n,i}, given by<cit.>:

    T_n[Φ_n] = ⟨Φ_n | T̂ | Φ_n ⟩,


    T_n[{φ_n,i}] = - ∑_i⟨φ_n,i | ∇^2 | φ_n,i⟩.




Further clarifications on the theorems are provided here. While current discussions primarily address pure ground states of interacting systems, extensions to ensemble states are feasible. This involves generalizing ρ_v to ensemble IVR densities and redefining energy functionals in Eqs. <ref> and <ref>, as detailed in Ref. . Specifically, for a g-fold ensemble state, Eq. <ref> is reformulated as:

    E_n,v^pw[{χ_n,p, λ_n,p}] = min_{Ψ^k,d^k}→ρ_{χ_n,p, λ_n,p}∑_k=1^g d^k ⟨Ψ^k | Ĥ_v | Ψ^k ⟩,

where the search for g anti-symmetric wavefunctions {Ψ^k} (each with its density ρ^k) and their weights {d^k} (satisfying 0 ≤ d^k ≤ 1 and ∑_k=1^g d^k = 1) ensures the ensemble density matches ρ_{χ_n,p, λ_n,p}, i.e., ∑_k=1^g d^k ρ^k(𝐫) = ρ_{χ_n,p, λ_n,p}. Concerning g, by assigning it a sufficiently large number, non-contributing weights will reduce to zero during the minimization process, thereby effectively capturing the system's true degeneracy. The energy functionals E_n,v^pw[Φ_n] and  E_n,v^pw[{φ_n,i}] are adapted similarly. Notably, even degenerate systems managed with a single Slater determinant Φ_n benefit from using hypercomplex orbitals {φ_n,i}, equating to the application of HCOs and their occupations {χ_n,p, λ_n,p}, thus effectively representing ensemble IVR densities and strong correlation effects <cit.>.



Moreover, the theoretical framework can be extended to handle excited states, including ensemble excited states. According to the Levy constrained search method and its extension <cit.>, when the trial density corresponds to an excited state, it identifies a stationary point rather than the global minimum. Thus, by modifying the minimization process in Eqs. <ref> and <ref> to target these stationary points, the theorems can be applied to compute the energies and densities of excited states.


While Φ_n, {φ_n,i}, and {χ_n,p, λ_n,p} all serve as fundamental variables of the system, only HCOs and their occupations, {χ_n,p, λ_n,p}, are real values <cit.>. Additionally, Theorem 1 indicates that {χ_n,p, λ_n,p} possess the smallest domain, enhancing their operability and optimizations. Furthermore, the fractional occupations of HCOs can effectively capture strong correlation effects, offering advantages in functional development.
Notably, at n=0, HCO occupations are strictly binary—either 1 or 0—clearly delineating occupied from unoccupied orbitals. This scenario allows the above three theorems to lay a rigorous theoretical foundation for various (occupied and unoccupied) orbital functionals <cit.> and corresponding orbital optimization <cit.> implementations. In fact, even at n=0, this approach can circumvent many issues associated with KSVR densities, thereby explaining the favorable outcomes observed in existing orbital-optimization calculations <cit.>.
When n=∞, the constraints on {χ_n,p, λ_n,p} imposed by Eqs. <ref> and <ref> align precisely with the ensemble N-representability constraints on natural orbitals and their occupations, thereby enabling these theorems to demonstrate the rigor of reduced density matrix functional theory <cit.> from another perspective.


The above addresses the spin-compensated case of N/2 electrons for each spin. As for spin-polarized systems containing N_α and N_β electrons with differing spins, Φ_n then comprises N_α orbitals {φ_n,i^α} and N_β orbitals {φ_n,i^β}, each linked to distinct sets of HCOs and their occupations, specifically {χ_n,p^α, λ_n,p^α} and {χ_n,p^β, λ_n,p^β}. In this scenario, the densities—including ρ_v, ρ_w_s, and ρ_n—are accordingly decomposed into components for α and β spins. Consequently, similar derivations can be carried out for both spins, details of which are not repeated here.


In summary, this study integrates various functionals—those dependent on occupied, unoccupied, or fractionally-occupied orbitals—into a unified framework. Each type of functional corresponds to a Clifford algebra with a specific dimensional basis, with their theoretical rigor thoroughly demonstrated. Specifically, Theorem 3 establishes a variational principle that lays a robust theoretical foundation for both orbital optimization and simultaneous orbital and occupation optimization. Furthermore, Theorem 2 investigates the domains of densities determined by the basis dimensions within Clifford algebras, confirming that these methodologies can effectively overcome the limitations posed by the KSVR assumption. Moreover, when properly implemented, these methodologies can also resolve common issues encountered in GGA and hybrid functional calculations, particularly when the desired solutions incorrectly target excited states of noninteracting systems—this work ensures the validity of such implementation.







Support from the National Natural Science Foundation of China (Grants No. 22122303 and No. 22073049) and Fundamental Research Funds for the Central Universities (Nankai University, Grant No. 63206008) is appreciated. 




§ DETAILS OF THE PROOF OF THEOREM 1



Theorem 1: The mappings among Φ_n, {φ_n,i}, and {χ_n,p, λ_n,p} are surjective for a given n, as shown in:

    Φ_n ⇔{φ_n,i}⇔{χ_n,p, λ_n,p}.

Specifically, there exists a one-to-one mapping between Φ_n and {φ_n,i}, and a many-to-one mapping between {φ_n,i} and {χ_n,p, λ_n,p}, within the unitary rotation of {φ_n,i} and of {χ_n,p} with the same occupation. 

Proof 1: The one-to-one surjective mapping between Φ_n and {φ_n,i}, within the unitary rotation of {φ_n,i}, is evident. The demonstration below illustrates the many-to-one surjective mapping between {φ_n,i} and {χ_n,p, λ_n,p}, within the unitary rotation of {φ_n,i} and of {χ_n,p} with the same occupation. 


For completeness, a brief review of relevant concepts is provided. Hypercomplex orbitals are defined as <cit.>

    φ_n,p(𝐫) = ϕ_p^0(𝐫) + ∑_μ=1^nϕ_p^μ(𝐫) e_μ,

where {e_1, e_2, …, e_n} form a basis in a Clifford algebra of dimension n <cit.>. For constructing Φ_n, N/2 hypercomplex orbitals are utilized; without loss of generality, the first N/2 orbitals, denoted {φ_n,i}, are selected. The set {ϕ_p^μ} consists of real functions expanded on a set of orthonormal functions {ξ_p}, which does not compromise generality:

    ϕ_p^μ(r) = ∑_q=1^K ξ_q(r) V_pq^μ,

where K is the dimension of the basis set, and {V^μ} are a set of n+1 K × K matrices, each V^μ representing the expansion coefficients for the μ-th component of the hypercomplex orbitals. These matrices preserve the orthonormality of the hypercomplex orbitals, leading to the following conditions on {V^μ}:

    ∑_μ=0^nV^μV^μ T = I_K,  
    ∑_μ=0^nV^μ TV^μ = I_K,  
    V^μV^ν T = V^νV^μ T,  
    V^μ TV^ν = V^ν TV^μ,

where I_K is the K × K identity matrix, and the superscript T denotes the matrix transpose.
Through a series of derivations <cit.>, the density can be represented on the set {χ_n,p, λ_n,p}, as given by:

    ρ_n(𝐫) = 2 ∑_p λ_n,p |χ_n,p(𝐫)|^2.

Here, {χ_n,p, λ_n,p} correspond to the eigenvectors and eigenvalues, respectively, of the symmetric matrix D_n, defined by:

    D_n = ∑_μ=0^nV^μ TI^N/2_K V^μ,

where I^N/2_K is a K × K diagonal matrix, with the first N/2 diagonal entries set to 1, and the rest set to 0. The spectral decomposition of D_n is given by:

    D_n = U Λ U^T.

where Λ is a diagonal matrix containing the eigenvalues {λ_n,p}, and U, a unitary matrix, defines {χ_n,p} as χ_n,p = ∑_q=1^K ξ_q U_qp. Thus, {χ_n,p} are orthonormal:

    ⟨χ_n,p | χ_n,q⟩ = δ_pq.

Interestingly, the constraints on {λ_n,p} depend on the dimension n of the Clifford algebra. Specifically, for any n, all {V^μ} that satisfy Eq. <ref> define the domain Ω_n for {χ_n,p, λ_n,p}. Notably, at n=0, Ω_n=0 encompasses {χ_n,p} that meet Eq. <ref>, with N/2 of them having an occupation of 1 and the rest 0, thereby clearly delineating occupied from unoccupied orbitals. For further discussion, see Refs.  and .

The surjectivity of the mapping between {φ_n,i} and {χ_n,p, λ_n,p} is first established. For any orthonormal set {φ_n,i}, a corresponding set {V^μ} exists that satisfies Eq. <ref>, thus ensuring the presence of {χ_n,p, λ_n,p} in Ω_n. Conversely, for each {χ_n,p, λ_n,p} in Ω_n, there must be a set {V^μ} compliant with Eq. <ref>, defining a compatible set {φ_n,i}. Without such a set {V^μ}, {χ_n,p, λ_n,p} would not fall within Ω_n. This confirms the surjective nature of the mapping between {φ_n,i} and {χ_n,p, λ_n,p}.


For any unitary rotation of the first N/2 hypercomplex orbitals {φ_n,i}, it corresponds to the following transformation on all matrices in {V^μ}:

    V'^μ = W V^μ,

where W is defined as:

    W=[            Q O_N/2, K-N/2; O_K-N/2, N/2      I_K-N/2 ],

with Q being a N/2 × N/2 unitary matrix, and O_I,J denoting an I × J zero matrix. It can be derived that for any {V^μ} satisfying Eq. <ref>, and its corresponding D_n defined by Eq. <ref>, applying the transformation in Eq. <ref> to all matrices in {V^μ} yields a new set {V'^μ} that also satisfies Eq. <ref>, producing the same D_n as defined by Eq. <ref>. Therefore, the transformed and original matrices correspond to the same {χ_n,p, λ_n,p}. Moreover, if D_n has identical eigenvalues, i.e., multiple orbitals in {χ_n,p} have the same occupations, any unitary rotation among these orbitals, by the spectral decomposition defined in Eq. <ref>, leads to the same D_n before and after the transformation, thus corresponding to the same {V^μ} and their defined {φ_n,i}. In fact, the invariance under these transformations also applies to the energy functionals defined in this work. The mappings between {φ_n,i} and {χ_n,p, λ_n,p} discussed below do not account for the changes these transformations impart to {φ_n,i} and {χ_n,p, λ_n,p}.


Lastly, the many-to-one mapping between {φ_n,i} and {χ_n,p, λ_n,p} is confirmed. For a given set of {φ_n,i}, there is a specific set of {V^μ} that satisfies Eq. <ref> and a distinct matrix D_n, thus associating a unique set of {χ_n,p, λ_n,p} with it. Additionally, if the roles of two bases in all hypercomplex orbitals (Eq. <ref>), such as e_μ and e_ν, and the functions ϕ_p^μ and ϕ_p^ν are interchanged, this is equivalent to swapping the coefficient matrices V^μ and V^ν. Clearly, such a transformation cannot be achieved by a unitary rotation of {φ_n,i}; however, it is evident that the new matrices {V'^μ} also satisfy Eq. <ref>, and the D_n defined by Eq. <ref> remains unchanged after the transformation. Therefore, the mapping between {φ_n,i} and {χ_n,p, λ_n,p} is many-to-one. This discussion only mentions a simple transformation; there are evidently other more complex transformations involving hypercomplex orbitals that leave {χ_n,p, λ_n,p} unchanged, hence the domain of {χ_n,p, λ_n,p} is considerably smaller.  □

aip

 
