Ensemble Ground State of a Many-Electron System with Fractional Electron Number and Spin: Piecewise-Linearity and Flat-Plane Condition Generalized

Description of many-electron systems with a fractional electron number (Ntot) and fractional spin (Mtot) is of great importance in physical chemistry, solid-state physics, and materials science. In this Letter, we provide an exact description of the zero-temperature ensemble ground state of a general, finite, many-electron system and characterize the dependence of the energy and the spin-densities on both Ntot and Mtot, when the total spin is at its equilibrium value. We generalize the piecewise-linearity principle and the flat-plane condition and determine which pure states contribute to the ground-state ensemble. We find a new derivative discontinuity, which manifests for spin variation at a constant Ntot, as a jump in the Kohn–Sham potential. We identify a previously unknown degeneracy of the ground state, such that the total energy and density are unique, but the spin-densities are not. Our findings serve as a basis for development of advanced approximations in density functional theory and other many-electron methods.

Modelling of materials and chemical processes from first principles crucially depends on efficient, accurate, yet approximate methods (see, e.g., [1][2][3][4][5][6]).The quality of a given approximation can be assessed not only versus the experiment, but also by comparison to known exact properties of many-electron systems; in the context of density functional theory (DFT), see, e.g., [7][8][9] One such exact property is the piecewise-linear behavior of the energy versus electron number [10]: when varying the number of electrons in the system, N tot , in a continuous manner, allowing both integer and fractional values, the energy E(N tot ) is linear between any two integer values of N tot .As N tot surpasses an integer, the slope of E(N tot ) may abruptly change.Piecewise linearity exists also for the electron density n(r) and any property that is an expectation value of an operator.
Spin, which is a fundamental property of an electron, has to be taken into account in many-electron systems, both when magnetic fields are present, but also in their absence [11][12][13].Occurrence of fractional z-projection of the total spin, M tot , and the dependence of the energy on M tot , has been extensively studied in the literature [14][15][16][17][18][19][20][21][22].In Ref. 15 it was shown that the exact total energy must be constant for all M tot ∈ [−S min , S min ] (where S min is the equilibrium value of the total spin, S, i.e., that value of S for which the system has the lowest energy).In Ref. 16 the flat-plane condition has been developed, unifying piecewise-linearity [10] and the constancy condition [15], which are both completely general principles, and apply to any many-electron approach.Focusing on the H atom, Ref. 16 described its exact total energy, while continuously varying N tot from 0 to 2 and M tot accordingly, while M tot ∈ [− 1  2 , 1  2 ]. [23] The exact energy graph of H in the N ↑ − N ↓ plane (where N ↑ = 1  2 N tot + M tot and N ↓ = 1  2 N tot − M tot ) is comprised of two triangles whose vertices lie at the points of integer N ↑ and N ↓ .Approximate exchange-correlation (xc) functionals may strongly violate this condition, leading to delocalization error and static correlation error.These developments were followed by numerous studies (see, e.g., [24][25][26][27][28][29] and references therein) which focused on: (a) the constancy condition [15], examined by changing M tot between −S min and +S min while keeping N tot constant and integer; or (b) varying both N tot and M tot for systems with S min = 1  2 .In addition, Refs.14, 17, 20 and 22 analyzed the flat-plane behavior of the energy for general N ↑ and N ↓ , suggesting formulae describing the energy profile.They assumed that the total energy (and analogously other properties) for any fractional N tot and M tot is determined by the energies of three out of the four states with nearby integer numbers of ↑ and ↓ electrons.The overall energy graph is therefore a collection of triangles.Moreover, Refs.30 and 18 raised the possibility of a discontinuity in E versus M tot , which stems from the nonuniqueness of the external magnetic field in spin-DFT.
This work goes beyond previous results, and rigorously describes the exact ground state of a general, finite, many-electron system with an arbitrary fractional electron number, N tot , and a simultaneously fractional zprojection of the spin M tot , as long as the total spin, S, is at its equilibrium value, S = S min ; S min itself can be of any value.Such a ground state has to be described by an ensemble, attained by minimizing the total energy, while constraining N tot and M tot .
Our findings generalize the piecewise-linearity principle [10] and the flat-plane condition [16].We find which pure ground states contribute to the ensemble, and which do not.The contributing pure states are not necessarily among the four neighboring states to the point (N tot , M tot ), contrary to previous findings.Notably, the energy graph consists not only of triangles, but also of trapezoids.Furthermore, we identify a new type of a derivative discontinuity, which manifests in the case of spin variation, as a jump in the Kohn-Sham (KS) potential and compare it to the discontinuities defined in Ref. 12. Surprisingly, we reveal a degeneracy in the ensemble ground state, where the energy and the total den- sity are determined uniquely, but the spin-densities are not.Our findings provide new exact properties of manyelectron systems, useful for development of advanced approximations in DFT and in other many-electron methods.We start our derivation by considering an ensemble ground state Λ of a system with a possibly fractional N tot and fractional M tot , where no magnetic fields are present.The (non-relativistic) Hamiltonian is given by Ĥ = T + V + Ŵ , where T is the kinetic energy operator, V is a multiplicative external potential operator and Ŵ is the electron-electron repulsion operator.Subsequently, Λ minimizes the expectation value of the energy, Tr Λ Ĥ , under the constraints and Here N is the electron number operator and Ŝz is the operator for the z-projection of the total spin.Consider now |Ψ N,S,M ⟩, the lowest-energy eigenstate of the operators N , Ŝ2 , Ŝz , with eigenvalues N , S(S + 1) and M , respectively.N is a non-negative integer, M and S are integers (for even N ) or half-integers (for odd N ), and −S ⩽ M ⩽ S. For a given value of N , there exists one such value of S that minimizes the energy; we denote it S min .Subsequently, there exists a multiplet of (2S min + 1) states, all with the same energy (in absence of a magnetic field), and with M = −S min , . . ., S min .We assume that for all other values of S (each one with its own multiplet of states) the energy is strictly higher.We further assume, for simplicity, that for specified N , S and M , the pure ground state is not degenerate any further.
In the following we denote |Ψ N,M ⟩ as the lowestenergy eigenstate of N and Ŝz , with eigenvalues N and M , respectively.For M ∈ [−S min , S min ], we have, of course, |Ψ N,M ⟩ = Ψ N,Smin(N ),M ; varying M does not affect the energy E(N, M ) = ⟨Ψ N,M | Ĥ|Ψ N,M ⟩.For other values of M , the ground state energy is higher.The energy ⟨Ψ N,S,M | Ĥ|Ψ N,S,M ⟩ as a function of S and M is as illustrated in Fig. 1.Moreover, as long as M ∈ [−S min (N ), S min (N )], the density n N (r) = ⟨Ψ N,M |n(r)|Ψ N,M ⟩ does not depend on M (see Supporting Information (SI) for details).
For general, possibly fractional, values of N tot and M tot , the ground state Λ is expressed as Λ = [31] Here and below the sum over N includes all integers for which the N -electron system is bound, and the sum over M includes all possible (integer or half-integer) values of M , for a given value of N : In the following, we denote N tot = N 0 + α, where N 0 is the integer part and α ∈ [0, 1) is the fractional part of N tot .The equilibrium spin values for N 0 and N 0 + 1 electrons are, respectively, S 0 := S min (N 0 ) and S 1 := S min (N 0 + 1).
We now consider the important case of and prove that any ensemble state, which consists of only |Ψ N0,−S0 ⟩ , . . ., |Ψ N0,S0 ⟩ and |Ψ N0+1,−S1 ⟩ , . . ., |Ψ N0+1,S1 ⟩, and satisfies Constraints (1)-(3), is a ground state of the system.By explicitly writing the aforementioned ensemble state as and making use of Eqs. ( 1) and (3) (for Γ), we see that is the ground state energy for a given N .Notably, there is no dependence on the z-projection of the spin; Constraint (2) was not used.Now we show that Tr Γ Ĥ is indeed the ground state energy.
First, from the piecewise-linearity principle of Ref. 10, it follows that min Ξ →Ntot Tr Ξ Ĥ = (1 − α)E(N 0 ) + αE(N 0 + 1), namely that the ground state energy for any system with N tot electrons, represented by a state Ξ, is linear in N , between N 0 and N 0 + 1.Second, note that E ens (N tot , M tot ) := min Ξ →Ntot,Mtot Tr Ξ Ĥ ⩾ min Ξ →Ntot Tr Ξ Ĥ , merely stating the ensemble groundstate energy for given N tot and M tot is greater or equal to the ensemble ground-state energy given only N tot , because adding an additional constraint in a minimization yields a result larger or equal the original one.Third, since Γ satisfies Eqs. ( 1 This inequality holds only if all its terms equal each other, meaning that Tr Γ Ĥ = min Ξ →Ntot,Mtot Tr Ξ Ĥ , i.e., Γ is indeed a ground state. Next, to show that every ground-state is of the form of Γ (Eq.( 5)), still within Region (4), we assume, by way of contradiction, that this is false for a ground state Λ, which is explicitly written as Λ = to deduce properties of Λ.The sum on M in Eq. ( 7) runs over both positive and negative values, as before.Λ′ satisfies Eqs. ( 1) and (3), Tr Λ′ Ŝz = 0, and Tr Λ′ Ĥ = Tr Λ Ĥ = (1 − α)E(N 0 ) + αE(N 0 + 1).Therefore, Λ′ is itself a ground state, of N tot electrons and spin 0.
To find the coefficients λ N,M that bring the energy of Λ′ to minimum, while retaining its expectation values of 1l, N , Ŝz , we notice that, for each value of N (in the external sum of Eq. ( 7)), the energy is minimized by setting λ N,M = 0 for all M > S min (N ) and M < −S min (N ).Then, Λ′ = . This expression for the ensemble energy is familiar from Ref. 10: Using the convexity conjecture of the energy[3, 19,32], we state that this energy is minimized by still under the aforementioned constraints for Λ′ .This means that the only nonzero coefficients of Λ are λ N0,−S0 , . . ., λ N0,S0 and λ N0+1,−S1 , . . ., λ N0+1,S1 , in contradiction to the original assumption.This concludes our proof: the ground state of a system with N tot electrons with a total spin M tot is described only by the ensemble ground state Γ of Eq. ( 5).It consists The ground-state energy Eens(N ↑ , N ↓ ) for the C atom, within the region defined in Eq. ( 4), and for 0 ⩽ Ntot ⩽ 7. The slope is zero along the orange line segments, which correspond to constant and integer Ntot, and changes discontinuously in directions perpendicular to these lines. of 2(S 0 + S 1 + 1) pure states, at most.The coefficients γ N,M ∈ [0, 1] have to be determined from Constraints (1)-(3), which can be expressed as: Several important consequences follow from the fact that Γ of Eq. ( 5) is a ground state.1.Total energy.Within Region (4), the total energy (10) is piecewise-linear with respect to α and does not depend on the z-projection of the spin, M tot .The graph of E ens (N ↑ , N ↓ ) forms a series of triangles and trapezoids, as presented in Fig. 2. For the H atom, with the number of electrons in each spin channel varying from 0 to 1, this result is precisely the flat-plane condition as described in Ref. 16.However, for larger N tot and particularly for S min > 1 2 , e.g., the C atom of Fig. 2, the behavior of E ens is richer.Furthermore, since E ens does not depend on the spin, the slope in Fig. 2 along lines of constant N tot is (∂E ens /∂M tot ) Ntot = 0.This result is in agreement with Ref. 15.

Frontier eigenvalues and derivative discontinuity.
We denote M B = (1 − α)S 0 + αS 1 , the highest value of M tot (for a given N tot ), at the boundary of Region (4), and consider |M tot | < M B .In this case, (∂E ens /∂N ↑ ) N ↓ = (∂E ens /∂N ↓ ) N ↑ = (∂E ens /∂N tot ) Mtot , because (∂E ens /∂M tot ) Ntot = 0. Hence, by Janak's theorem [33], the highest occupied (ho) KS energy eigenvalues are the same for both spin channels, constant strictly within each trapezoid/triangle as in Fig. 2, and equal (the negative of) the corresponding ionization potential: . However, at and beyond the boundary of Region (4), the value of ε σ ho may differ from E(N 0 + 1) − E(N 0 ), as it is determined by the slope of the energy just outside and adjacently to the boundary.Whereas the full discussion of the energy profile outside Region (4) is beyond the scope of this work, we note that just outside the boundary the energy is a plane determined by three pure-state energy values, one of which is outside Region (4).Therefore, a change in the energy slope is expected as we cross the boundary M tot = M B of Region (4).Specifically, for spin migration (variation of M tot at constant N tot ), we predict a derivative discontinuity for the KS energy levels, which is manifested in a uniform jump in the KS potential v σ KS (r) at M tot = M B [a similar logic applies for the boundary We define the spin-migration derivative discontinuity around M tot = M B as where δ → 0 + , N tot is kept constant and i refers to any of the KS energy levels.In the following, quantities evaluated at M B + δ are denoted with a prime (e.g.ε ′ σ ho ) and those evaluated within Region (4) are without a prime.
As an illustration, choose N tot / ∈ N and consider a common example, where immediately outside Region (4), for M tot > M B , the ensemble ground state consists of |Ψ N0,S0 ⟩, |Ψ N0,S0+1 ⟩ and |Ψ N0+1,S1 ⟩.Then, the ensemble energy there is E ′ ens (N tot , M tot ) = KN tot +JM tot +C, where J = E(N 0 , S 0 + 1) − E(N 0 , S 0 ) -a spin-flip energy, and K = (E(N 0 + 1) − E(N 0 )) − (S 1 − S 0 )J.Therefore, the slope changes around M B equal Notably, the ↑-slope is continuous for , we obtain a discontinuity only in the ↓-channel: We note that as M tot surpasses M B , N ↓ crosses an integer, while N ↑ does not; the ↓-ho orbital within Region ( 4) is the ↓-lu orbital immediately outside it.The term in the parentheses in Eq. ( 14) is (the negative of) the KS gap immediately outside Region (4).It equals ε ↓ ho−1 − ε ↓ ho in terms of quantities inside Region (4).
In the uncommon case of S 1 > S 0 + 1 2 , we expect discontinuities in both spin channels, while neither N ↑ nor N ↓ cross an integer.
As to the value of the ho level at the boundary itself, it is determined by recalling that all the derivatives with respect to N σ are taken from below.Then, from geometric considerations in the In other words, in the above cases the ho level at the boundary has the same value as strictly within (4).Otherwise, at the boundary it has the value immediately outside (4).
∆ σ sm is related to the xc correction to the spin stiffness defined in Ref. 12, for the case of integer N tot : if immediately outside Region (4) the ensemble ground state consists of |Ψ N0,S0 ⟩ and |Ψ N0,S0+1 ⟩, the two quantities coincide.However, in the general case they differ, as Ref. 12 expressly makes use of excited-state spin ensembles, following Ref.34, and here we treat ground-state ensembles.3. Non-uniqueness and degeneracy.Remarkably, the coefficients γ N,M of Eq. ( 5) are confined by only three constraints, Eq. ( 9), whereas the number of coefficients is 2(S 0 + S 1 + 1).This means that the ground state Γ may not be uniquely defined.Strictly within Region (4) (i.e. ∈ N, all 2(S 0 + S 1 + 1) coefficients may be nonzero.Then, since there are 3 constraints (Eq.( 9)), the degeneracy of the ground-state is then we are left with 2S 0 + 1 coefficients and 2 constraints, hence the degeneracy is (2S 0 − 1)-fold.We discuss this non-uniqueness here in detail, starting with the cases where Γ is defined uniquely.
Case (a).If S 0 = 0 and S 1 = 1 2 (e.g., adding an electron to H + towards H, while varying M tot ), we are left with three uniquely determined γ N,M 's: 2 and S 1 = 0 (e.g., adding an electron to H towards H − , while varying M tot ), we also have three uniquely defined coefficients: At the boundary of Region (4), where the spin is at its maximal value, M tot = M B , we find that the expression for M tot , , reaches its maximum (under the aforementioned constraints on γ N,M 's) only when γ N0,S0 = 1 − α, γ N0+1,S1 = α and all the other γ N,M 's vanish; the ground state is determined uniquely.If, for example, S 1 = S 0 + 1 2 , this corresponds to the common scenario of ionization (via the up spin channel): fractionally varying N ↑ while N ↓ is a constant integer.A unique solution is obtained also when M tot reaches its minimal value, i.e., for Case (c).For the scenario of spin migration, when N tot = N 0 = const.and M tot varies between −S 0 and S 0 , the ground state is not defined uniquely, for S 0 ⩾ 1.From Eq. ( 9) we see that S1 M =−S1 γ N0+1,M = 0 and therefore γ N0+1,M = 0, for all M .We are then left with 2S 0 + 1 coefficients.
In particular, for S 0 = 1 (e.g., spin migration for the C atom), the three coefficients γ N0,1 , γ N0,0 and γ N0,−1 can be expressed as γ N0,1 , where x remains undetermined.To understand the meaning of the above ambiguity in the ground state, focus on the case , meaning that to get the average ⟨ Ŝz ⟩ = 0, the system must have an equal probability of x 2 to be in the M = 1 and M = −1 states, and a complementary probability of (1 − x) to be in the state with M = 0.Here we see a clear distinction of constraining the average spin of the system being 0 versus requiring each and every replica in a macroscopic statistical ensemble to have a spin of 0. The latter is obtained by demanding that the ground state of the system is an eigenstate of Ŝz , with eigenvalue 0 (which is equivalent here to setting x = 0).Surprisingly, despite the ambiguity in Γ, in this Case the density, as well as the spin-densities, are determined unambiguously and do not depend on x.Since the purestate densities n N,M (r) do not depend on M (see the SI for details), n ens (r) = n N0 (r).For the spin-densities, , where δ ↑ = 1 and δ ↓ = −1.The spin distribution, defined by Q(r) := n ↑ (r) − n ↓ (r), is expressed for our ensemble as being linear in M tot , and independent of x.Here we used the fact that for any pure state Case (d).For spin migration with S 0 = 3 2 (e.g., for the N atom), even the spin-densities are not determined uniquely anymore.Here we have four coefficients (dropping the index N 0 for brevity): , with two undetermined parameters, x and y.Whereas the total energy and the total density are determined unambiguously, the spin-densities linearly depend on y and are independent of x.The ensemble spin distribution is linear in y.
Case (e).For a fractional N tot = N 0 + α, S 0 = 1 2 and S 1 = 1 (e.g., adding an electron to C + towards C, while varying M tot ), we end up with five coefficients, which depend on two parameters, x and y: γ The total ensemble density n ens (r) = (1 − α)n N0 (r)+αn N0+1 (r) is piecewise-linear in α and remains fully determined, whereas the spin-densities are not.The ensemble spin distribution is being linear in y, but independent of α (in full contrast to n ens (r)).Cases (d) and (e) clearly show that, for a given electron number and spin, we have a set of degenerate ground states Γ(x, y), with the same total density, but with different spin-densities, n ↑ ens (r) and n ↓ ens (r).Notably, even in the spin-unpolarized case M tot = 0, the spin-densities are not determined, and are not necessarily equal to each other.4. Removing non-uniqueness.To further explore and interpret the surprising result of non-uniqueness at the ground state, we offer two ways for its full or partial removal.
First, to determine the ground state unambiguously, one can introduce additional constraints, to complement (1)-(3).For example, we can set ⟨ Ŝ2 z ⟩, ⟨ Ŝ3 z ⟩ and/or higher moments of Ŝz .In the SI we analyze how many such moments are sufficient to determine Γ and/or Q ens (r), in the general case.Specifically, for Case (c) above, ⟨ Ŝ2 z ⟩ = x, i.e., setting ⟨ Ŝ2 z ⟩ fully determines Γ.Alternatively, one can require the standard deviation ∆S z := ⟨ Ŝ2 z ⟩ − ⟨ Ŝz ⟩ 2 to be minimal.This requirement sets the ground state uniquely; a mathematical proof of this statement is provided in the SI.In particular, in Case (c) the deviation is In Case (d), ⟨ Ŝ2 z ⟩, as well as all even moments of Ŝz equal ⟨ Ŝ2k z ⟩ = (S 0 − 1) 2k (1 − x) + S 2k 0 x, being linear in x and independent of y (k ∈ N).Conversely, all the odd moments ⟨ Ŝ2k+1 z ⟩ = (S 0 −1) 2k y +S 2k 0 (M tot −y) are linear in y and independent of x.Therefore, it is sufficient to set ⟨ Ŝ2 z ⟩ and ⟨ Ŝ3 z ⟩, or to minimize ∆S z (see the SI for the technical details).In the latter case, we obtain a gradual transition from with two adjacent pure states involved at each point.The spin distribution is unambiguously determined and piecewiselinear in M tot : Second, to remove non-uniqueness of the ground state it is natural to introduce a weak magnetic field [35].Surprisingly, while lifting the degeneracy between pure states |Ψ N,M ⟩ with different M , a magnetic field does not fully remove the non-uniqueness.
Applying this magnetic field in Case (d), the ensemble energy becomes Ẽens (N 0 , M tot ) = E(N 0 ) + µ B B 0 M tot F N0,3/2 + µ B B 0 (F N0,1/2 − F N0,3/2 )y, being linear in y.One can then minimize Ẽens (N 0 , M tot ) with respect to y: If, say, for a given system and field profile F N0,1/2 > F N0,3/2 , then for B 0 > 0, we seek for the lowest possible y and for B 0 < 0 -for the highest.The range of y values is confined by the require- and Fortunately, minimizing y in this Case also uniquely sets x, and therefore the ground state is determined uniquely.As in the case of minimal ∆S z , also here we are left with only two pure states involved in the ensemble at each point.Now, however, we perform a gradual transition from In contrast, in Case (c), the ensemble energy in presence of the magnetic field becomes Ẽens (N 0 , M tot ) = E(N 0 ) + µ B B 0 M tot F N0,1 , being independent of x.Therefore, the non-uniqueness of the ground state cannot be removed.
In Case (e), the ensemble energy becomes Ẽens )y, and can be minimized with respect to y.However, as opposed to Case (d), here a minimal/maximal value of y can correspond to a range of possible x values; the ground state is not fully determined.This is not surprising, as Case (e) has Case (c) as a particular limit, for α = 1 (see the SI for details).
Finally, we note that applying a strong magnetic field could further remove the ambiguity in the ground state.However, this scenario may change which pure states contribute to the ensemble ground state, and therefore it is outside the scope of this Letter.
To conclude, in this Letter we rigorously described the exact ground state of a finite many-electron system, with N tot electrons and with spin M tot , fully employing the ensemble approach.Our results hold for any value of N tot , any values of the equilibrium spin, S min (N ), and for M tot confined by Eq. ( 4).Description outside of Region ( 4) is a subject for future work.
First, we found that the ground state is an ensemble, which is comprised of 2(S 0 + S 1 + 1) pure states, as indi-cated in Eq. ( 5).These states are not necessarily the four integer neighboring points to (N tot , M tot ) on the N ↑ − N ↓ grid.The graph of the total energy, E ens (N tot , M tot ) consists therefore not only of triangles, but also of trapezoids, when S min ⩾ 1.The total energy, the total density, and any quantity A, whose pure-state expectation values ⟨Ψ N,M | Â|Ψ N,M ⟩ do not depend on M , are piecewiselinear in N tot (i.e., in α) and independent of M tot .
Second, the highest-occupied KS ↑and ↓-orbital energies are found to be equal strictly within Region (4), and the KS potentials experience a spatially uniform "jump" at the boundary of this region (Eq.( 11)), directly related to the spin-migration derivative discontinuity, ∆ σ sm (Eq.( 11)), introduced here for the first time.
Third, we discovered that, generally speaking, the ground state Γ is not uniquely defined.Unlike the total energy and the total density, which are determined uniquely, the spin-densities are not (Eqs.( 16) and ( 17)).Two ways of removing this non-uniqueness are suggested: adding constraints, e.g., on higher moments of Ŝz , and minimizing ∆S z .Applying a weak magnetic field does not always fully remove the non-uniqueness.
Our findings generalize the piecewise-linearity and the flat-plane conditions and provide new exact properties for many-electron systems.These are expected to be useful in the development of advanced approximations in density functional theory and in other many-electron methods.This document provides further technical details on four subjects.In Sec.I we show that certain pure-state expectation values, including the density, ⟨Ψ N,S,M |n(r)|Ψ N,S,M ⟩, do not depend on M .In Sec.II we show how the ensemble ground state Γ and/or the spin distribution Q ens (r) can be uniquely determined by setting the higher moments of the operator Ŝz .In Sec.III we prove that enforcing the requirement that the deviation ∆S z is minimal always fully determines the ground state.In Sec.IV we provide a detailed analysis of Cases (c), (d) and (e), which were introduced in the main text, in absence and in presence of a weak, inhomogeneous magnetic field.
In the following, reference to an equation from the main text starts with the letters MT (e.g., Eq. (MT:4)), whereas references to equations within this document are by a single number.

I. DEPENDENCE OF PURE-STATE QUANTITIES ON M
In this section we clarify which pure-state quantities are independent of the value of M , as long as M ∈ [−S min , S min ], and which are not.Consider a quantity described by an operator P and compare P (M ) = ⟨Ψ N,M | P |Ψ N,M ⟩ with, say, P (M − 1).Recall that Ŝ± |Ψ N,M ∓1 ⟩ = const.|Ψ N,M ⟩, where Ŝ± = Ŝ † ∓ = Ŝx ± i Ŝy are the ladder operators, which raise/lower M by 1. Furthermore, note that |Ψ N,M ⟩ is an eigenstate of the operator Ŝ− Ŝ+ .Then, It then directly follows that for any operator P that commutes with Ŝ+ , P (M − 1) = P (M ), i.e., the value of P is M -independent.Such is the case for the Hamiltonian Ĥ, in absence of a magnetic field, for the electron density operator n(r), for the reduced density matrix, ρ1 (r, r ′ ), and for any operator that depends only on spatial, and not on spin-coordinates.However, this is not the case for the Hamiltonian in presence of a magnetic field, for the spin-density operator nσ (r) or for the spin-dependent reduced density matrix γ1 (rσ, r ′ σ ′ ).

II. HIGHER MOMENTS OF Ŝz
In this section we discuss the possibility to determine the ground state Γ uniquely by choosing the values of moments of Ŝz (i.e.⟨ Ŝz ⟩, ⟨ Ŝ2 z ⟩, ⟨ Ŝ3 z ⟩, etc.), in addition to Eqs. (MT:1)-(MT:3), and the possibility to determine the spin distribution Q ens (r) = n ↑ ens (r) − n ↓ ens (r), by choosing the values of odd moments of Ŝz .Recall that the first moment is M tot = ⟨ Ŝz ⟩.
First, consider N tot ∈ N (as relevant in Cases (c) and (d) in the main text).There are 2S 0 + 1 ensemble coefficients, and the moments of Ŝz are given by  The (2S 0 + 1) × (2S 0 + 1) matrix of Eq. ( 1) is a Vandermonde matrix [2], where no two columns are identical, hence it is invertible.Therefore, the moments of Ŝz up to Ŝ2S0 z determine the ensemble coefficients, and thus Γ. * Author to whom correspondence should be addressed: eli.kraisler@mail.huji.ac.il arXiv:2308.03465v3[cond-mat.mtrl-sci]20 Feb 2024 Next we show that Q ens (r) is determined by the odd moments, ⟨ Ŝk z ⟩, where 1 ⩽ k ⩽ 2S 0 .Importantly, the density n ens (r) = n ↑ ens (r) + n ↓ ens (r) is unique, depending only on N tot .Thus, if Q ens (r) is determined uniquely (while Γ may remain undetermined), the spin-densities n σ ens (r) are determined from Q ens (r) uniquely, as well.If for instance N tot is an odd integer, then 2S 0 is an odd integer, and the odd moments of Ŝz up to ⟨ Ŝ2S0 z ⟩ can be expressed as The matrix in Eq. ( 2) is similar to that of Eq. ( 1), however the powers in each row are not consecutive integers, hence it is not a Vandermonde matrix but a generalized Vandermonde matrix [2].This matrix is also invertible, meaning that the odd moments ⟨ Ŝz ⟩, . . ., ⟨ Ŝ2S0 z ⟩ determine the differences (γ N0,1/2 − γ N0,−1/2 ), . . ., (γ N0,S0 − γ N0,−S0 ).Since Q N,M (r) = −Q N,−M (r), the ensemble spin distribution can be expressed as r), meaning it is also determined by the odd moments ⟨ Ŝz ⟩, . . ., ⟨ Ŝ2S0 z ⟩.One can show in a similar manner that if N tot is an even integer then Q ens (r) is determined by ⟨ Ŝz ⟩, ⟨ Ŝ3 z ⟩, . . ., ⟨ Ŝ2S0−1 z ⟩.Next we consider the case of N tot / ∈ N, where there are now 2(S 0 + S 1 + 1) ensemble coefficients.We express the moments of Ŝz , as well as α, by If this matrix is invertible, then indeed the ground state is determined by the values of the moments of Ŝz , up to ⟨ Ŝ2S0+2S1 z ⟩, together with α.We have explicitly verified that this is the case for S 0 , S 1 ⩽ 8, using the NumPy [1] library in Python.This covers essentially all the chemically relevant scenarios.However, if this matrix is not invertible, then it is sufficient to choose the value of ⟨ Ŝ2S0+2S1+1 z ⟩ as an additional constraint (consistently with all the other constraints).Then, The matrix of Eq. ( 4) is another Vandermonde matrix, hence it is invertible.
To determine Q ens (r), now with N tot / ∈ N, assume that N 0 is odd.Then, analogously to Eq. ( 2), The matrix in Eq. ( 5) is another generalized Vandermonde matrix, meaning it is invertible.Since Q ens (r) = r), we see that Q ens (r) is fully determined by the odd moments of Ŝz , up to ⟨ Ŝ2S0+2S1 z ⟩.This can be shown for the case of N 0 being even, in the same way.

III. MINIMIZING ∆Sz
In this section we prove that Γ is determined uniquely by adding the constraint of minimizing ∆S z , or equiv-alently, minimizing ⟨ Ŝ2 z ⟩, as ⟨ Ŝz ⟩ = M tot is given by Eq. (MT:2).Consider an ensemble Γ (Eq.(MT:5)), which minimizes ⟨ Ŝ2 z ⟩.Since the ensemble coefficients are all between 0 and 1, we may write them as γ N,M = sin 2 (θ N,M ).Using Lagrange multipliers, with Eqs.(MT:1)-(MT:3) as constraints, we find where κ 1 , κ 2 , κ 3 are some constants.Taking the derivatives, we have 2 sin Hence, for each N, M such that γ N,M ̸ = 0, 1, we have Since Eq. ( 7) is quadratic in M , there are at most two solutions for each value of N .This means there are at most four nonzero coefficientsγ N0,a , γ N0,b , γ N0+1,c , γ N0+1,d .Here M = a, b are the solutions to Eq. ( 7) for N = N 0 , and M = c, d are the solutions for N = N 0 + 1.If Eq. ( 7) has only one solution for N = N 0 or N = N 0 + 1, then a = b or c = d, respectively.In addition, from Eq. ( 7) we have Without loss of generality, we assume that b < a.Furthermore, since M tot = γ N0,a a + γ N0,b b, and all γ N,M 's are within [0, 1], we realize that b < M tot < a.By rewriting Eq. ( 9) as Therefore a is the smallest integer or half-integer (in accordance with N 0 ) greater than M tot , and b is the largest integer/half-integer smaller than M tot .In particular, a = b + 1.Therefore, the ground state minimizing ⟨ Ŝ2 z ⟩ is unique in this case.
(ii) Consider the case N tot / ∈ N.
We define Γ0 =  8).Therefore, by minimizing ∆S z , we are left with at most three nonzero γ N,M 's, which can then be derived from Eqs. (MT:1)-(MT:3), uniquely.Particularly, in case γ N0+1,d = 0 and is ruled out, Eqs.(MT:1)-(MT:3) read: The determinant of the matrix above equals a − b = 1.Therefore, the matrix is invertible.If instead both γ N0+1,c and γ N0+1,d are present, but γ N0,b vanishes, the corresponding matrix is invertible, as well.Finally, we show that a, b and c are uniquely determined.Assume by contradiction that there exists another ensemble state Γ′ , which also minimizes ⟨ Ŝ2 z ⟩ under the constraints (MT:1)-(MT:3).Thus, there is a different triad a ′ , b ′ and c ′ , which corresponds to a minimum.But then we may refer to Γ′′ = 1 2 ( Γ + Γ′ ): it also minimizes ⟨ Ŝ2 z ⟩ under the same constraints, but consists of more than three pure states, which is a contradiction.In case the nonzero ensemble coefficients are, e.g., γ N0,a , γ N0+1,c and γ N0+1,d , the same logic applies.

IV. ANALYSIS OF CASES (C), (D) AND (E)
In this section we present in detail Cases (c), (d) and (e), which were introduced in the main text.We start by considering spin migration at integer electron number, namely, the situation when N tot = N 0 = const., i.e. α = 0 and M tot ∈ [−S 0 , S 0 ].As mentioned in the main text, from Eq. (MT:9) we realize that S1 M =−S1 γ N0+1,M = 0, and since γ N,M ∈ [0, 1], it follows that γ N0+1,M = 0 for all M .As we are left with γ N0,M only, for brevity, we may omit in what follows the index N 0 .Yet, when it serves clarity we add it back.From Eq. (MT:9) we then get S0 M =−S0 γ M = 1 and S0 M =−S0 M γ M = M tot -two equations to determine (2S 0 + 1) coefficients.We therefore expect that for S 0 ⩾ 1 the coefficients γ N,M , and therefore the ensemble state Γ, would not be determined unambiguously.
We now refer to specific scenarios of spin migration, for selected values of S 0 .Trivially, for S 0 = 0, there is no spin migration at all, γ 0 = 1 and M tot must equal 0.
For S 0 = 1 2 (e.g., spin migration for the H or for the Li atom), which is a particular scenario within Case (a), we find that The coefficients are determined unambiguously, and so are all the quantities that describe the system.In particular, -independent of M tot .More generally, all the even moments of Ŝz , ⟨ Ŝ2k z ⟩ = ( 1 2 ) 2k are M tot -independent and all the odd moments ⟨ Ŝ2k+1 To find the ensemble spin-densities, n σ ens (r), we recall that: (i) for any pure-state, |Ψ N,M ⟩, the spin-densities relate as n ↓ N,M (r) = n ↑ N,−M (r); (ii) the total density is M -independent (see Sec.I above).Then, where and Q(r) := n ↑ (r) − n ↓ (r) is the spin distribution.From Eq. ( 12) it directly follows that the total ensemble density n ens (r) = n N0 (r) is M tot -independent and the spin distribution is linear in M tot .For all cases of integer N tot , the spin polarization, , has an analogous form to that of Q ens (r).
For S 0 = 1 (e.g., spin migration for the C atom), which is Case (c) of the main text, we obtain three coefficients that are not identically zero.They can be expressed as: where x remains undetermined.The value of x is not arbitrary, though: it is confined by the requirement that all γ N,M ∈ [0, 1].Applying this requirement to all three coefficients above, we find that Whereas the total energy, E ens (N tot , M tot ) = E(N 0 ) and all the odd moments of Ŝz , ⟨ Ŝ2k+1 z ⟩ = M tot , are determined uniquely, the even moments of Ŝz equal ⟨ Ŝ2k z ⟩ = x.Therefore, to fully determine the ensemble ground state for this Case, it is possible to add a requirement as to the value of, e.g., ⟨ Ŝ2 z ⟩ (in agreement with the conclusions of Sec.II above).
Alternatively, it is possible to restrict the deviation, ∆S z , to be minimal.Accounting for the aforementioned range of the variable x, we find that (∆S z ) min = |M tot | − M 2 tot , as depicted in Fig. 1.The ensemble coefficients in this case equal (see Fig. 2).Note that for every value of M tot we obtain two nonzero ensemble coefficients at most.Therefore, as we increase M tot from -1 to 0, we see a linear transition from the pure state |Ψ N0,−1 ⟩⟨Ψ N0,−1 | to |Ψ N0,0 ⟩⟨Ψ N0,0 |; further increase of M tot from 0 to 1 gives a linear transition from |Ψ N0,0 ⟩⟨Ψ N0,0 | to |Ψ N0,1 ⟩⟨Ψ N0,1 |.At the points M tot = −1, 0, 1, we have pure states (only one coefficient is nonzero) and therefore the deviation ∆S z vanishes (see Fig. 1).For completeness, consider the scenario where ∆S z is taken to be maximal.Then, x = 1 and (∆S z ) max = 1 − M 2 tot (depicted in Fig. 1, as well).For such a case, The spin-densities for this Case can be written at first as n σ ens (r) = 1 2 (1 − x)n N0,0 (r) + 1 2 x n N0,1 (r) + 1 2 δ σ M tot Q N0,1 (r), but since the pure-state densities are M -independent, it directly follows that, similarly to the S 0 = 1 2 case, n σ ens (r) satisfies Eq. ( 12) and does not depend on x, and Q ens (r) satisfies Eq. (14).
In presence of a magnetic field, the total energy in Case (c) changes linearly with M tot , as described in the main text.It does not depend, however, on x, and therefore cannot be minimized with respect to x. Surprisingly, introduction of a magnetic field does not remove the ambiguity in the ground state, for this Case.Therefore, in this respect, the requirement on ∆S z is stronger.
For S 0 = 3 2 (e.g., spin migration for the N atom), which is Case (d) of the main text, we obtain four coefficients, which are not identically zero.They can be expressed using the variables x and y: Whereas E ens (N tot , M tot ) = E(N 0 ) and ⟨ Ŝz ⟩ = M tot are determined uniquely, other quantities that describe the system are not.The even moments of Ŝz equal ⟨ Ŝ2k z ⟩ = S 2k 0 x + (S 0 − 1) 2k (1 − x), and are linear in x and independent of y, whereas the odd moments of Ŝz are ⟨ Ŝ2k+1 z ⟩ = S 2k 0 (M tot −y)+(S 0 −1) 2k y and are linear in y and independent of x.To fully determine the ensemble ground state for this Case, we can add two constraints, e.g., for ⟨ Ŝ2 z ⟩ and ⟨ Ŝ3 z ⟩.As opposed to previous cases, here even the spindensities are not determined unambiguously.They can be expressed as being independent of x and linear in y.The spin distribution Q ens (r) can then be written as it is linear in y, as well.For y = 0 it reduces to Eq. ( 14).
As we already mentioned above, one of the ways to fully determine the ground state is to require that the deviation ∆S z is minimal.In our Case, ∆S z = 1 4 − M 2 tot + 2x, which means that the minimal value of ∆S z is obtained for the minimal value of x.As before, the domain of x (and y) is restricted by the requirement that all the ensemble coefficients are within [0, 1].Superimposing all eight conditions that emerge from this requirement, we are left with two inequalities: |y| ⩽ 1 2 (1 − x) and |y − M tot | ⩽ 3 2 x.Graphically, the domain of x and y (which parametrically depends on M tot ) is illustrated in Fig. 3. From Fig. 3 we conclude that the lowest value of x is Observing the geometry of the x−y domain, we find that for each value of x min there is one corresponding value of y: Therefore, by minimization of ∆S z , the ground state is defined completely, in agreement with Sec.III.The dependence of x min and y 0 on M tot is depicted on Figs. 4  and 5, respectively.For completeness, let us note that the maximal value of x is x max = 1, for any M tot , and y 1 = y(x max ) = 0. Therefore, in this case also maximizing ∆S z sets both x and y and fully determines the ensemble ground state.The minimal value of ∆S z is and the maximal value is (∆S z ) max = 9 4 − M 2 tot ; see Fig. 6.
Finally, we discuss Case (e) of the main text: N tot = N 0 + α, α ∈ [0, 1), S 0 = 1 2 and S 1 = 1.This Case corresponds to adding an electron, e.g., to C + towards C, while varying M tot .In this Case we have five ensemble coefficients, with three constraints (Eq.(MT:9)).The coefficients can be expressed using two variables, x and y: (cf.Eq. ( 18)).The spin-densities are independent of x and linear in y.It directly follows that the total density n ens (N tot , M tot ) = (1 − α)n N0 (r) + αn N0+1 (r), (37) is piecewise-linear in α and unambiguously defined, as expected, while the spin distribution is given by It depends only on y, and can be determined by setting ⟨ Ŝ3 z ⟩, according to Sec.II.As in previous Cases, we address the scenario of minimal ∆S z .Here, ∆S z = 1 4 (1 + 3α) − x − M 2 tot ; a minimum is obtained for the highest possible x.The domain of x and y, which emerges from the requirement that all the ensemble coefficients are within [0, 1], can be described by the inequalities: 0 ⩽ x ⩽ 1, |y| ⩽ 1 2 (1 − α) and |y − M tot | ⩽ α − x.The x − y domain is illustrated in Fig. 13.
In Case (e), the highest value of x is .
(39)  , and for M tot ∈ [ 3α−1 2 , 1+α 2 , ], x = 0. Figs. 17 and 18 show y min and the range of x graphically.Therefore, while the magnetic field in this case can set the spin-densities, the odd spin moments, and other quantities, the ground state remains ambiguous (due to x).

1
Fritz Haber Research Center for Molecular Dynamics and Institute of Chemistry, The Hebrew University of Jerusalem, 9091401 Jerusalem, Israel (Dated: February 21, 2024)