Geometrization for Energy Levels of Isotropic Hyperfine Hamiltonian Block and Related Central Spin Problems for an Arbitrarily Complex Set of Spin-1/2 Nuclei

Description of interacting spin systems relies on understanding the spectral properties of the corresponding spin Hamiltonians. However, the eigenvalue problems arising here lead to algebraic problems too complex to be analytically tractable. This is already the case for the simplest nontrivial (Kmax−1) block for an isotropic hyperfine Hamiltonian for a radical with spin-12 nuclei, where n nuclei produce an n-th order algebraic equation with n independent parameters. Systems described by such blocks are now physically realizable, e.g., as radicals or radical pairs with polarized nuclear spins, appear as closed subensembles in more general radical settings, and have numerous counterparts in related central spin problems. We provide a simple geometrization of energy levels in this case: given n spin-12 nuclei with arbitrary positive couplings ai, take an n-dimensional hyper-ellipsoid with semiaxes ai, stretch it by a factor of n+1 along the spatial diagonal (1, 1, …, 1), read off the semiaxes of thus produced new hyper-ellipsoid qi, augment the set {qi} with q0=0, and obtain the sought n+1 energies as Ek=−12qk2+14∑iai. This procedure provides a way of seeing things that can only be solved numerically, giving a useful tool to gain insights that complement the numeric simulations usually inevitable here, and shows an intriguing connection to discrete Fourier transform and spectral properties of standard graphs.


Introduction
Solving a Hamiltonian for a physical system with n free parameters naturally leads to the need to solve an algebraic problem with n free parameters, e.g., derive and solve a secular equation to find the energies of eigenstates as the roots of an at least n-th order polynomial. However, unless the situation is very symmetric, allows a perturbative treatment, or can be physically factored into lower-dimensional subtasks, anything beyond a cubic equation is an analytical dead-end, and one is faced with the choice of either resorting to simple models with no more than three parameters, or having to go numeric.
This work focuses on the particular problem of describing the energy levels of the spin system of an organic radical with electronic spin-1 /2 interacting with spin-1 /2 nuclei via isotropic hyperfine interaction described by Hamiltonian This is a representative case of the central spin problem in science, one selected spin interacting with multiple non-interacting spins of a different kind, or a localized spin coupled to a "spin bath," and parallels can be found in many problems of condensed matter and chemical physics and quantum information processing. Some examples apart from the already quoted organic radical or radical pairs include electrons in quantum dots [1] and which has dimensions (n + 1) × (n + 1), is a separate entity of key importance for describing recombining spin-correlated radical pairs that feature one or both radicals with an arbitrary set of n spin-1 /2 nuclei. This importance has two aspects. The first one relates to a very general situation of photo-or radiation-induced radical ion pairs in low-viscosity liquids that are formed in singlet initial spin state of the two unpaired electrons and are observed via recombination also from singlet collective spin state [55][56][57]. Such pairs are reasonably well described by considering only dynamic spin evolution driven by intra-radical isotropic hyperfine interactions, omitting inter-radical spin-spin interactions and spin relaxation, and therefore are described by a sum of two independent Hamiltonians of the type (1), which conserves total spin projection of the pair. All such systems include a subensemble of pairs starting from initial singlet electron spin state with all nuclear spins aligned in the same direction. Such a subensemble is present for a radical pair of arbitrary hyperfine complexity, and under stated assumptions its evolution is independent of other subensembles that must have a different total spin projection of the pair. Therefore, its contribution to recombination probability is additive. Furthermore, this evolution can be described using only the state subspaces of the individual radicals with either maximum total (electron plus nuclei) spin projection, or total projection one less than maximum (2). Thus, to describe this subensemble it suffices to solve the "penultimate" block of the hyperfine Hamiltonian for a radical, as the state with the maximum total spin projection is automatically an eigenstate and needs no efforts. Although the statistical weight of the "all nuclear spins up" subensemble of radical pairs rapidly decreases with increasing the number of nuclei, its additive contribution to any experimental observable is always present, and as discussed in [54] its behavior is in a sense typical for more complex subensembles with lower total spin projection. In particular, a perturbative treatment of applied weak magnetic field is possible in the most general case and gives such insights as an expected linear skew in magnetic field dependence of recombination probability upon passing through zero field if the weights of "all nuclear spins up" and "all nuclear spins down" subensembles are not equal, i.e., nuclear spins are polarized. Such a skew indeed was implicitly present in model calculations [58], although it has yet to be observed in experiment. This would be an interesting possibility of introducing anisotropy (via polarized nuclear spins, i.e., via the initial state) in an otherwise intrinsically isotropic spin system.
The second aspect is related to the abovementioned option of having radicals with initially polarized nuclear spins. One of the crucial paradigm shifts in NMR over the last two decades is a realization that non-Boltzmannian nuclear population (hyperpolarization) is possible and can be created by several means. The already implemented approaches include dynamic nuclear polarization DNP [59,60], chemically induced dynamic nuclear polarization CIDNP [61,62], optical nuclear polarization ONP [63], optical pumping of noble gases [64] or semiconductors [65], para-hydrogen induced polarization (PHIP) [66] and Signal Amplification By Reversible Exchange (SABRE) [67,68]. Achieving an orders of magnitude boost in NMR sensitivity, these methods are capable of producing diamagnetic systems with nuclear polarizations of percents to dozens of percent persisting for times of the order of nuclear relaxation times, seconds for hydrogens to minutes for nuclei with lower magnetogyric ratios, or maintained indefinitely in case of continuous generation. Much longer-lived collective symmetry-protected nuclear spin states can be created [69][70][71][72][73], and can be used as reservoirs for accumulating spin polarization that can then be released into the spin system. For the simplest possible systems of two spin-1 /2 nuclei the symmetryprotected state is usually the non-magnetic singlet state, with net magnetization carried by triplet states, and methods are now available for efficient conversion of triplet states to singlet and back (magnetization-to-singlet, M2S, and singlet-to-magnetization, S2M, conversion) [74][75][76][77][78]. Larger systems of up to 4 spins have been successfully polarized [79,80]. While hyperpolarization of organic molecules was initially mostly confined to hydrogens due to availability of a convenient source of hyperpolarized hydrogens in the form of molecular para-hydrogen, the more recent trend is polarization transfer to heteronuclei-like 4 of 29 13 C and 15 N [81][82][83][84][85][86]. Furthermore, the most efficient redistribution of hyperpolarization over the system of coupled nuclei occurs in the conditions of strong J-coupling, when scalar coupling of spins is larger than or comparable to the difference of their Zeeman energies. For homonuclear systems, like J-coupled hydrogens, this usually happens in milliTesla fields, while efficient coupling of heteronuclei like 13 C-1 H pairs requires mi-croTesla to nanoTesla fields due to vastly different magnetogyric ratios, which was one of the main drivers for developing ZULF NMR. Thus, the new reality in the NMR community is availability of hyperpolarized diamagnetic systems in close to zero fields with amounts, lifetimes, and degrees of polarization already sufficient for in-vivo use as MRI agents. But then the new experimental reality for the radical/radical pair community is availability of hyperpolarized starting systems to generate radicals from, with substantial polarization at 13 C and 15 N spin-1 /2 nuclei that often show higher hyperfine coupling constants than protons and will thus dominate the hyperfine structure of the radical. To the best of our knowledge, no such experiments on radicals or radical pairs with initially hyperpolarized nuclei have been reported so far, but this is certainly an emerging experimental possibility.
Thus, the penultimate block of the hyperfine Hamiltonian is indeed a "physically reasonable part" of the problem that may be practically realizable by itself, provides useful insights on more complex parts of a more general problem, and still allows fairly deep analysis without resorting to numeric simulations. The paper [54] had derived the secular equation for this block and certain spectral bounds for its eigenvalues. Still for n nuclei the problem has n hyperfine coupling constants a i as free parameters, and denies a general analytic solution, although the eigenvalues can certainly be found numerically by diagonalizing the initial matrix.
It often happens that things that are easy to imagine or draw free-hand are difficult to describe analytically. The prime example of this is computer graphics, which has long become a research field in applied mathematics in itself. Occasionally this complexity may work the other way round, and things that are difficult to treat analytically may happen to have an intuitively natural and appealing visualization. This work will provide such a visualization for the energies of the penultimate Hamiltonian block described above: given n spin-1 /2 nuclei with arbitrary positive couplings a i , take an n-dimensional hyper-ellipsoid centered at the origin of the reference frame and having semiaxes √ a i aligned with the n Cartesian axes, stretch it by a factor of √ n + 1 along the spatial diagonal (1, 1, . . . , 1), and the semiaxes of the obtained new hyper-ellipsoid q i , augmented with one special value q 0 = 0, produce the n + 1 energy levels as E k = − 1 2 q 2 k + 1 4 ∑ i a i . While such a procedure cannot yield analytic expressions for the eigenvalues and these would be computed numerically anyway, as the roots of the corresponding high-order polynomial cannot be found analytically for algebraic reasons, it is believed that such thinking of Hamiltonians in terms of geometric shapes can be useful in supporting the intuition and providing insights to complement the usually inevitable numeric calculations.
This work does not compete with the existing spin dynamics simulation software [87][88][89][90][91][92][93][94][95][96][97] and pursues a different goal. Its raison d'etre is not how to compute, but rather how to not compute and still obtain useful insights about spin dynamics. The purpose is to obtain a picture of the core, spectral, properties of the central spin Hamiltonian that will eventually show up in any simulation, rather than to generate a visual representation of spin dynamics. The notion of (hyper-)ellipsoids central to this work may also sound misleadingly familiar. Any purely dynamic spin system completely described by Hamiltonian can in one sense or another be represented by some ellipsoids, as there is a set of purely real eigenvalues that can be all made positive by a suitable shift of energy origin, and these are algebraically equivalent to ellipsoids. On the other hand, ordinary 3D-ellipsoids are very familiar for magnetic resonance community, as they are used to represent, e.g., chemical shielding and shift tensors, dipole-dipole and quadrupolar couplings, and atomic anisotropic displacement parameters (thermal ellipsoids) derived from NMR crystallography in NMR [98], or anisotropic hfi tensors and g-tensors in EPR [99]. However, the point of this work is a very special arbitrary-dimensional ellipsoid with direct connection to initial parameters a i and a startlingly simple, universal visualization that does not require any computation at all for an arbitrary complex system.
All the derivations in this work rely on basic complex/matrix algebra covered in university textbooks, e.g., [100][101][102] and probably do not require specific references, several more specific simple facts that are used are reproduced just inline instead of citing them. Connections to more remote areas that are not further explored here, such as Laplacian matrices for graphs and Discrete Fourier Transform, will only be indicated at suitable places. The usual conventions of diag{a 1 , a 2 , . . . , a n },1, |M|, (. . .) T and (. . .) + will be used to denote a diagonal matrix with the only nonzero elements being a i at diagonal, identity matrix, determinant of matrix, matrix transpose and Hermitian conjugate, ket notation |φ k will also be used for vectors.
The paper is organized as follows. First we briefly introduce the used notation and summarize the necessary results from [54] that will be used as the starting point for the derivations of this work, and formulate the problem of obtaining the eigenvalues (roots) of a certain matrix derived from Hamiltonian block. Then a zero-splitting transform is built in two flavors, an explicitly complex unitary one and a real orthogonal one, which both remove the present zero eigenvalue from the matrix/secular equation and produce a more convenient positive definite matrix, either complex Hermitian or real symmetric, respectively, inheriting all the non-zero roots. The real zero-split matrix is then used to introduce the circulant form and as a bridge to the reference case of equivalent nuclei. The scaling and diagonalizing transforms are built for the complex zero-split matrix M that together produce a diagonal matrix diag{a 1 , . . . , a n } from which M can be reconstructed as a simple factorization, the inverse matrix for M is constructed, and matrix M is used to treat the loss of equivalence as a perturbation. In the last sections the stretching transform outlined above is built using n-dimensional hyper-ellipsoid as a geometric image of the positive definite matrix, and the circle is completed by showing that the original characteristic polynomial can be obtained by variation of a functional produced by the stretching. The concluding section summarizes the essential result and briefly discusses its connection to solving Hamiltonians for systems of interacting spins.

Notation and Previous Results
This section introduces the notation and adopts the necessary results from [54] for subsequent derivations. We consider a radical having n spin-1 /2 nuclei with arbitrary positive hyperfine coupling constants 0 < a 1 ≤ a 2 ≤ . . . ≤ a n described by isotropic hyperfine Hamiltonian (1) and its restriction to state subspace with total projection K = S z + ∑ i I iz onto an arbitrarily chosen axis z that is one step less than the maximum possible value (2), spanned by n + 1 product-basis functions |ξ 0 = |−; + + + · · · + ; |ξ 1 = |+; − + + · · · + ; . . . |ξ n = |+; + + · · · + − . (3) Here ± stand for projections m = ± 1 2 , the first position in the notation separated with semicolon is allotted to electron spin, and the following n positions correspond to n nuclear spins, with the subscript indicating the only flipped spin in the entire set. All matrices and vectors will be written with basis ordered by increasing index. The restriction to blocks with fixed K is a natural approach for this problem, as the Hamiltonian (1), either by itself or augmented with Zeeman term, conserves the total spin projection and is thus block-diagonal with respect to K. As discussed in detail in [54], the block with K = K max = (n + 1)/2 is trivial, the blocks with |K| < K max − 1 are too complex for an arbitrary system, and the penultimate block K = K max − 1 is still reasonably simple to analyze for a general system, while already rich enough to build a representative physical implementation.
The restricted Hamiltonian for this subspace is represented by an (n + 1) × (n + 1) real symmetric matrix It is convenient to shift and scale the matrix and change its overall sign as follows: All embellishments were removed from the operator symbol H for simplicity. We refer to eigenvalues λ k of (5) as "roots," and they are related to the corresponding energies E k of the original Hamiltonian via an obvious relation The eigenvectors are not affected by transform (5). It can be readily seen that the vector consisting of all ones is an eigenvector of (5) with λ = 0 and E = 1 4 ∑ i a i , as the sum of elements in every line is equal to 0. The adopted form differs from the form used in [54] by overall sign change, which makes all non-zero λ k positive as is more convenient for geometrization. It may be further noted that matrix (5) looks like Laplacian, or Kirchhoff matrix for a certain graph [103], with degrees of vertices at the diagonal and minus weights of edges as off-diagonal elements. This particular matrix corresponds to a star graph shown in Figure 1, and more complex subspaces with lower |K| values will be represented by more general bipartite graphs with two classes of vertices corresponding to electron spin up or down in the relevant ket vector of the type (3), with only inter-class connections by edges with weights given by the hyperfine coupling constants for a nuclear spin flipped when flipping the electron spin to pass from a function of one class to a function of another class. The simple transform (5) is thus useful not only because it provides a convenient picture of the Hamiltonian block, but rather because it reduces it to a mathematically well-defined and thoroughly studied object [103]. For example, a Laplacian matrix for any graph has a zero eigenvalue with multiplicity equal to the number of its connected components (which is one in our case), all other eigenvalues being positive. The eigenvectors are not affected by transform (5). It can be readily seen that the vector consisting of all ones is an eigenvector of (5) with 0   and 1 4 , i i E a   as the sum of elements in every line is equal to 0. The adopted form differs from the form used in [54] by overall sign change, which makes all non-zero k  positive as is more convenient for geometrization. It may be further noted that matrix (5) looks like Laplacian, or Kirchhoff matrix for a certain graph [103], with degrees of vertices at the diagonal and minus weights of edges as off-diagonal elements. This particular matrix corresponds to a star graph shown in Figure 1, and more complex subspaces with lower K values will be represented by more general bipartite graphs with two classes of vertices corresponding to electron spin up or down in the relevant ket vector of the type (3), with only inter-class connections by edges with weights given by the hyperfine coupling constants for a nuclear spin flipped when flipping the electron spin to pass from a function of one class to a function of another class. The simple transform (5) is thus useful not only because it provides a convenient picture of the Hamiltonian block, but rather because it reduces it to a mathematically well-defined and thoroughly studied object [103]. For example, a Laplacian matrix for any graph has a zero eigenvalue with multiplicity equal to the number of its connected components (which is one in our case), all other eigenvalues being positive. We also frequently meet the following symmetric combinations of the coupling constants: In the adopted notation the needed results from [54] that we build upon can be summarized as follows. The secular equation to find the roots can be written in three equivalent forms: We also frequently meet the following symmetric combinations of the coupling constants: In the adopted notation the needed results from [54] that we build upon can be summarized as follows. The secular equation to find the roots can be written in three equivalent forms:

•
The starting determinant: • The expanded and reduced determinant: The polynomial form: Equations (8) and (10) are explicitly valid for any selection of the coupling constants a i , while Equation (9) needs some care with equivalent nuclei (a i = a j for some i, j), but is also valid for any set of coupling constants after all arising uncertainties have been removed. In [54] this was demonstrated in a somewhat awkward way, but, in a hindsight, it can be readily seen from continuity arguments, as both the origin (8) and the result (10) of (9) are obviously continuous in a i , and a set with equivalent nuclei can be obtained as a limit of a set with distinct nuclei.
Equations (8)-(10) have one zero root mentioned above, corresponding to eigenvector of all ones obtained by lowering the projection of the function with the maximum possible total spin and projection all the remaining n roots are non-zero and positive. If all coupling constants are distinct, 0 < a 1 < a 2 < · · · < a n , all λ k are also distinct and are separated by consecutive a i as and all the corresponding eigenvectors are fully entangled, i.e., all their coordinates in the basis (3) are non-zero. This case corresponds to the last line in (9). If equivalent nuclei with multiple equal coupling constants are allowed, and a radical has n spin-1 /2 nuclei with constants 0 < a 1 < a 2 < · · · < a r , r ≤ n with multiplicities m i ≥ 1, i = 1, 2, . . . , r ≤ n, the n + 1 roots of the characteristic polynomial are organized as follows: the remaining r roots λ i satisfy : The roots λ i = a i appearing for each multiple constant give rise to the ambiguity in (9) mentioned above and related to formal division by zero compensated by multiplication with zero pre-factor, but analysis showed that all three forms are valid for an arbitrary set of hyperfine coupling constants.
In practice the most useful of the three is the polynomial form (10), where when writing A i the constants a i should be counted for each of n nuclei considering the multiplicities. Still this is an n-th order algebraic equation with up to n free parameters A i equivalent to n original free coupling constants a i . The apparent simplicity of the polynomial in (10) is deceptive, and although only regular numeric prefactors (k + 1) in each term (k + 1)A k (−λ) n−k distinguish it from the trivial equation having λ i = a i by Vieta's formulae, these "simple numbers" lead to catastrophic consequences in terms of complexity. The following sections will show how this problem can be at least partially overcome by geometric visualization.

The Zero-Splitting Transform
The matrix (5), although written in a very symmetric form, has a zero eigenvalue with the known eigenvector of all ones. The zero λ 0 is an inconvenience, as it makes the matrix degenerate and limits possible manipulations with it, while the known eigenvector suggests that the subspace that it spans can be split away. In this section we shall construct such a zero-splitting transform.
We need an unitary transform to another (n + 1)−dimensional basis, of which one vector is already known: and we want the remaining n vectors to reflect the equal roles of all n nuclei, so that the new basis retains this symmetry. A convenient and regular way of building such a basis is to generate it by consecutive powers of n distinct n-th roots of unity ε k = ε k , where ε = exp(2πi/n) is the n-th root of unity with the smallest positive argument, as where the first element and the normalizing prefactor in each vector have to be chosen to ensure orthonormality of the entire set. The corresponding unitary transform matrix C 1 has the form and taking the transform by multiplying out H 1 = C + 1 HC 1 , we obtain where the following complex combinations of parameters a i weighted with consecutive powers of the corresponding root of unity have been introduced: We have indeed split off the subspace spanned by the eigenvector of λ 0 = 0 by obtaining the zero "hook" in the transformed matrix and a smaller n × n but explicitly complex Hermitian matrix that is in its way symmetric with respect to all n nuclei. It can also be seen that the second function of the new basis |φ 2 = (−n, 1, 1, . . . , 1, 1) T / n(n + 1) stands out from the set (15), as coefficients for elements in the second row and second column are larger than for other nonzero rows and columns. This is a reflection of the entanglement of the eigenfunctions of Hamiltonian (5) mentioned in the previous Section: at least for distinct nuclei they all must contain a nonzero contribution from function |ξ 0 = |−; + + + · · · + of the set (3), that is, their first component must be non-zero. However, the only basis vector of the new set (15) with non-zero first component is |φ 2 , which endows it with a special role.
In splitting away the zero eigenvalue we have obtained a smaller complex matrix (17), but sacrificed the real symmetry of the original matrix (5). This can be restored for the smaller matrix by making another unitary transform and "shuffling back" the basis in the following way. Again, use complex vectors (15) to generate a new basis, but this time confine it only to the subspace with nonzero eigenvalues leaving alone the split-off one-dimensional subspace. Arrange it in the columns of an unitary matrix as and apply the transform backwards by multiplying out H 2 = C 2 H 1 C + 2 to obtain with The real orthogonal transform that produces H 2 directly from H given by the matrix product C 1 C + 2 has the form with If only the real split-zero transformed matrix is needed, the matrix (20) can be obtained directly from (5) using real-only matrices by taking the transform H 2 = C + 12 HC 12 and completely bypassing three complex matrices in (16)-(19). We also note that the suggested transforms have turned a sparse matrix (5) into dense but regular matrices using also dense transform matrices: while the starting matrix had non-zero elements only at diagonal and in the hook, the transformed matrices have not identically zero elements everywhere but in the hook.
The transforms (16) and (20) are unitary by construction and preserve the spectrum of the transformed matrix, therefore the matrices (17) and (20) have equation (10) as their characteristic polynomial. The zero root corresponds to the zero hook of the first row and first column, and all the positive roots are inherited by the smaller positive definite n × n matrices that have the characteristic polynomial (24) The traces of matrices (17) and (20) are equal to 2A 1 , the coefficient at the second highest degree of the characteristic polynomial, and determinants of the n × n matrices are equal to (n + 1)A n , the free term. In the following sections we shall omit the zero hook and work with the obtained n × n matrices, referring to them as the result of the zero-splitting transform of this section. The natural ordering of the basis vectors generated by n-th roots of unity given above will be used throughout.

The Real Symmetric Zero-Split Matrix
The n × n real symmetric matrix in (20) is useful for making a connection to the wellknown reference case of all equivalent nuclei and introducing circulants. Omitting the zero hook and setting all a i = a, we obtain the following n × n matrix: It can be seen that all its rows are generated by a single line (2, 1, 1, . . . , 1) circularly shifted one position to the right when going one row down. The eigenvalues of such a matrix, which will be our roots for the case of n equivalent nuclei, are found by guessing its eigenvectors based on this circulant symmetry: these are again given by consecutive powers of an n-th root of unity. It can be seen that the system of linear equations is satisfied identically, as multiplying the first equation by ε produces the second equation and so on, and the last equation multiplied by ε gives again the first one, closing the circulant pattern across the rows. Taking n different n-th roots of unity ε k , k = 0, . . . , n − 1, we thus obtain n different orthogonal eigenvectors, and the corresponding eigenvalue is directly read off from the last equation and is given by Since for n-th roots of unity ε k the eigenvalues of Equation (26) are λ = a, multiplicity n − 1, λ = (n + 1)a, multiplicity 1.
This dichotomy of obtaining either 1 or n + 1 depending on whether summing all consecutive powers of ε = 1 or ε = 1 will be met in numerous places. The special role of the basis vector corresponding to ε = 1 (and thus comprised of all ones) as opposed to vectors generated by other roots of unity ε = 1 is also encountered quite often.
The obtained roots are consistent with the results (13), giving λ = a with multiplicity n − 1 and a nondegenerate root λ > a for n equivalent nuclei. The obtained roots are also consistent with energy levels that can be directly obtained from consideration of spin Hamiltonian (1) in case of n equivalent spin-1 /2 nuclei with restriction to subspace with total spin projection K max − 1. The Hamiltonian and its energies in this case read Three types of states can enter the considered subspace with K = K max − 1. One of these is the state with I Σ = I max = n/2 and J = I Σ + 1 2 = (n + 1)/2, i.e., the function (11) from the multiplet with maximum possible total spin that gave us the zero root split away in the preceding section. The others are one state with I Σ = I max = n/2 and J = I Σ − 1 2 = (n − 1)/2, and n − 1 states with I Σ = I max − 1 = n/2 − 1 and J = I Σ + 1 2 = (n − 1)/2, the latter corresponding to n − 1 possible ways of combining n nuclear spins-1 /2 into total nuclear spin I Σ = I max − 1 = n/2 − 1 with maximum projection, i.e., to n − 1 possible nuclear multiplets with I Σ = n/2 − 1. Taking into account the connection (6) between the roots λ k and energies E k , we obtain Table 1 that reproduces (29). Table 1. Energies E and corresponding roots λ for subspace K = K max − 1.
The characteristic polynomial (24) factors into a much simpler form On the other hand, knowing from the general properties that for n equivalent nuclei we must have n − 1 roots equal to a leaves only one unknown λ to figure out, and this can be done by considering the simplest possible relation connecting the sum of all roots with the coefficient at the second highest degree of λ In fact, it suffices to know that we must have a degenerate root λ 1 with multiplicity n − 1 and a non-degenerate root λ 2 , and consider the relations for the trace and determinant to obtain the sought roots: it can be seen that the system of equations is satisfied with λ 1 = a and λ 2 = (n + 1)a. Not much can be done with the real symmetric matrix to move beyond the equivalent nuclei, as the shifting pattern symmetry of the type seen in (20) is lost when the couplings become not all equal, and the circulant property is thus lost as well. However, further advances are possible with the complexified matrix of (17) that upon inspection also shows a form of circulant symmetry, but without reference to specific coupling constants a i , which are now hidden inside the complex "coupling constants" A ε k .

The Complexified Coupling Constants
The complexified "coupling constants" A ε k defined by Equation (18) are the third form of parametrizing the problem of finding the roots with the input parameters of the original physical problem, the hyperfine coupling constants a i . The set {a i } itself is the first and obvious parametrization used in the Hamiltonian (1) and its matrix representation (5). The second one is the set {A i } defined by (7), which is natural for the polynomial representation, but is not linearly related to {a i }. The set A(ε k ) can be viewed as a result of transforming the "vector" consisting of a i with an unitary matrix built from vectors again generated by consecutive powers of n-th roots of unity: where C 3 is the unitary matrix and the √ n multiplier has been added since the vectors corresponding to the original definition of (18) were not normalized. As opposed to the set {A i }, the mapping {a i } → A ε k is linear bijective and can be easily reversed if needed using the Hermitian conjugate of matrix from (34) and dividing the result by n. The inverse relation has the form A useful view of the mapping {a i } → A ε k is interpreting the original set {a i } as a function of an evenly spaced discrete argument {i}. Then the transform (34) is recognized as a discrete Fourier transform (DFT) for the "function" {a i }, and transform (35) is inverse DFT, which opens a connection to the vast mathematical field of fast Fourier transforms [104] and more algebraically abstract subject of Fourier transform in finite groups [105] and Gaussian sums over finite fields [106]. FFT is widely used in signal processing for efficient manipulation of polynomials with coefficients {a i } [107], and is not directly employed here, still it is a huge resource to draw from.
A useful bit of intuition from the connection with Fourier transform is that the slower the variation of the function, the faster its Fourier coefficients decay. Thus, a useful property of the set A(ε k ) is that for equivalent nuclei all A ε k vanish identically except for A(1) = na = A 1 as a direct consequence of the sum rule (28), since for equivalent nuclei A ε k turns simply into the sum of all powers of ε k times the coupling constant. The matrix (17) (with omitted zero hook) then turns into a diagonal one H eq = a × diag{n + 1, 1, 1, . . . , 1, 1}, reproducing eigenvalues (29) and factored characteristic polynomial (31). For small deviations from equivalence the matrix will be almost diagonal, with first superdiagonal and subdiagonal given by small A(ε) and its complex conjugate A ε n−1 , respectively, and progressively farther super-/subdiagonals given by further consecutive components of the Fourier transform A(ε k ) that will progressively decay if variations in the original set {a i } are modest. Thus, the complex split-zero matrix is natural for describing the loss of equivalence. However, before doing this, one additional step will be useful.

The Scaling Transform
Let us take the complex split-zero matrix (17), omit the zero hook, and denote the resulting n × n matrix as M : As mentioned earlier, M has inherited all the nontrivial positive roots for the original problem and has (24) as its characteristic polynomial. The form of matrix in (37) invites to re-scale the first basis vector as |φ 1 → |φ 1 / √ n + 1, which will remove the factors (n + 1) and √ n + 1 in the corresponding matrix elements of the first line and first column. Such a transformation is equivalent to matrix transform to produce This scaling transform is not unitary, and it does not preserve the spectrum. However, it changes the problem in a manageable and, in this place, in a convenient way. We also note that the scaling transform has modified the trace and determinant of the complex zero-split matrix M in the following way: The first relation is seen directly from the definition of M 1 in (39), and the second one comes from taking determinants at both sides of the scaling transform (38) and taking |M| = (n + 1)A n as the last term in characteristic polynomial (24).

The Loss of Equivalence as a Perturbation
Let us take the system of equations for the eigensystem of matrix (37) M → c = λ → c and scale it using the transform (38) by writing the following sequence of determinant identities: is the modified identity matrix with first element scaled down, to obtain the equation Taking all nuclei to be equivalent, we have A(1) = na, A ε k = 1 = 0, and (43) becomes This equation has unit vectors → e k having only one 1 in the k-th position, k = 1, . . . , n, as its solution ("eigenvectors"), and has one non-degenerate root λ = (n + 1)a with eigenvector → e 1 and (n − 1)−degenerate root λ = a with eigenvectors → e k , k = 2, . . . , n. This is of course again the solution for equivalent nuclei (29), but in a much more convenient way to apply the equivalence-lifting perturbation, because algebraically the problem now is to apply a perturbation to identity-like matrix, a well-developed topic covered in matrix theory books [108]. We shall only discuss first-order perturbation here, underlining its similarity to conventional quantum mechanical perturbation theory.
Equation (43) up to first order in perturbation can be written as where M 0 1 is the matrix from (44), δM 1 is the matrix from (43) with all A ε k , k = 0, . . . , n − 1 changed to the corresponding δA ε k , and unperturbed eigenvectors and eigenvalues are the solution of (44). As in quantum mechanics, we have two subproblems, for single nondegenerate state with with first component equal to zero to keep the perturbed vector → c 0 + δ → c normalized to 1 to first order. Equation (46) is satisfied identically in zero order, and in the first order we obtain Expanding, this gives As expected, the terms with δλ and δc i separate, and we obtain the first-order corrections to the non-degenerate root and its eigenvector: It can be seen that the correction to the root, i.e., to the energy of the eigenstate, depends only on δA(1) = δA 1 = ∑ i δa i , i.e., on the shift of the center of gravity of the set of coupling constants, while the correction to the function depends on δA ε k , i.e., on lifting the degeneracy of the set of coupling constants.
Since the state remains non-degenerate, its function must be explicitly real, and this is indeed so. The correction δ → c to the state vector written in basis {|η i } means that the function is corrected by |δψ = ∑ i δc i |η i . We recall that the complex split-zero matrix M is written in complexified shuffled basis of type (15), of which the needed part {|η k }, k = 2, . . . , n can be expressed via original basis (3) with an (n − 1) × n matrix as Taking δc i from (50) and |η i from (51), after rearrangement we obtain the following result for the correction and the original unperturbed function: Note that the function of the original basis with flipped electron spin |ξ 0 = |−; + + · · · + is present only in the unperturbed function ψ 0 1 and does not enter the perturbative correction |δψ 1 that only redistributes the functions {|ξ k }, k = 2, . . . , n.
Turning to the problem of the (n − 1)−degenerate subspace, we first note that this subspace was not modified by the scaling transform (38), and therefore the problem reduces to the ordinary quantum mechanical degenerate perturbation theory, i.e., to finding the eigenvalues and eigenvectors of the perturbation matrix restricted to the degenerate subspace. Since the basis of this subspace is just the set of unit vectors → e k , k = 2, . . . , n, we have to find the eigensystem of matrix from (43) with omitted first row and first column and all A ε k , k = 0, . . . , n − 1 again changed to the corresponding δA ε k , i.e., the (n − 1) × (n − 1) matrix written in basis (51) as The secular equation for the matrix (53) will produce an algebraic equation of degree n − 1, only one step lower than the original Equation (24). Since the matrix has δA(1) at diagonal, we can say that all eigenvalues will be shifted by δA(1)/n and then split around it as the roots of the polynomial with coefficients derived from the off-diagonal elements δA ε k , i.e., the splitting depends on lifting the degeneracy of the set of coupling constants. A finer point to note is that, as discussed at the end of Section 2.4, the matrix (53) depends not directly on δa i , but rather on components of the Fourier transform of their sequence, and thus the splitting depends not only on lifting the degeneracy, but also on how abrupt it occurs in the ordered sequence of coupling constants.
The splitting may be non-zero in linear in δa i order, for example, for the simplest possible case of n = 3 and a pair of degenerate states we obtain which may produce linear splitting: e.g., for δa 1 = δa, δa 2 = δa 3 = 0 we obtain δλ = 0, 2 3 δa. The corresponding eigenfunctions in this case are (|η 2 ± |η 3 )/ √ 2 that produce the familiar for three-fold symmetry (which in fact is not present here, as the nuclei are no longer assumed equivalent) combinations The case of four nuclei and three degenerate states yields a cubic equation with three real roots having their center of gravity at δA(1)/4, and even in this case the obtained explicit expressions are already too awkward to be useful and give no new insight, and increasing n pushes the problem beyond reasonable effort, while the pattern of lifting the degeneracy is given by the general result (13).
We also note that the degenerate perturbation, similar to the non-degenerate case discussed above, does not mix functions with different projections of electron spin, and all eigenfunctions of matrix (53) remain confined to subspace spanned by basis (51) and thus include only functions {|ξ k }, k = 1, . . . , n with electron spin projection m e = + 1 2 . This has far-reaching consequences for the original physical problem of describing recombining spin-correlated radical pairs as set forth in [54] in the sense that for first-order perturbations the set of degenerate states remains unreachable from the initial state of the electron-singlet radical pair with all nuclear spins up. The same is true for the case of a separate radical with polarized nuclei. The more difficult to describe degenerate states simply do not participate in spin evolution. Thus, small perturbations to equivalence only change the contribution of the non-degenerate state, which is a relatively mild correction with first-order changes to energies and wavefunctions, in comparison to possible lifting of degeneracy and not small, zero order, rotation of basis. As we see, the latter does not happen here.

The Diagonalizing Transform
The scaling transform (38) has produced the matrix (39) that also has a circulant pattern to it. It can be seen that all its rows are generated by the single line A(1), A(ε), A ε 2 , . . . , A ε n−1 circularly shifted one position to the right when going one row down. The eigenvectors of such a matrix are again given by consecutive powers of an n-th root of unity. Writing the eigenproblem for M 1 as it can be seen that vectors of the form satisfy the system (56), and the first line produces the expression for eigenvalue Substituting expressions for A(ε m ) (18) and rearranging terms, we obtain the sum The expression in parentheses for each term is the sum of all n consecutive powers of one of the n-th roots of unity and vanishes identically unless the root is unity, so the only surviving term in the sum is for k + j − 1 = 0(modn), or j = n + 1 − k(modn). We thus obtain the following eigenvalues for M 1 : k = 1 : µ 1 = a n ; k = 2 : µ 2 = a n−1 ; k = 3 : µ 3 = a n−2 ; k = j : µ j = a n−j+1 ; k = n : µ n = a 1 .
The backwards-going sequence of coupling constants a i in (60) can be reversed by taking in (57) ε n−k = ε k as the root of unity that generates the vector → c k , which finally produces eigenvalues and eigenvectors of M 1 as and the matrix serving as the obtained diagonalizing transform for M 1 is recognized as matrix C + 3 from Equation (35): . . , a n }.
Recalling that M 1 was itself produced from the complex zero-split matrix M via the scaling transform (38), we have where we again stress that the scaling transform with matrix C −1 4 is not unitary, and so the entire transform M → D is not unitary. This matrix relation can be reversed to reconstruct the complex zero-split matrix M from diagonal matrix D as Matrix M from (64) has equation (24) as its characteristic polynomial and gives all the nonzero roots of the original problem. It can be seen that it is constructed from a simple diagonal matrix of the original parameters a i by two simple transforms, (complex) rotation C 3 followed by stretching C 4 . The transform (64) can also be thought of as changing the characteristic polynomial in the following way: generating the Equation (24) from the trivial equation with λ i = a i .

The Inverse Matrix
The multiplicative decomposition (64) for the complex zero-split matrix makes it easy to obtain the inverse matrix to M by taking the inverse of both sides of (64): The inverse matrix is thus obtained from diagonal matrix D −1 built from reciprocal coupling constants b i = 1/a i via the same rotation C 3 and compression C −1 4 instead of stretching C 4 . Introducing reciprocal complexified coupling constants B ε k similar to (18) as we obtain for the inverse matrix W an expression similar to expression (37) for M : that differs from (37) only by changing all A ε k to their reciprocal counterparts B ε k and all multiplicative prefactors n + 1, √ n + 1 also to their reciprocals. This reciprocity-type connection between the matrix and its inverse continues to their spectral properties. If a non-degenerate matrix A has the set {λ i } as its spectrum, its inverse A −1 has spectrum µ i = λ −1 i with the same ordering of eigenvalues, and is diagonalized by the same transform. Indeed [101], if take Apparently, So, we can write and from uniqueness of the inverse matrix, we get with the same transforming matrix C, i.e., the same set of eigenvectors taken in the same order, and with reciprocal eigenvalues taken in the same order. Finally, the reciprocal connection continues to the characteristic polynomial (24). Taking advantage of the fact that λ = 0 is not a root, dividing Equation (24) by (−λ) n A n (n + 1), and introducing µ = λ −1 , we obtain the following pair of equations: where we naturally obtain the reciprocal counterpart of definitions (7) The roots of the two equations in (74) are related as µ i = λ −1 i . The essence of last sections can be more usefully summarized if we treat the coupling constants a i simply as numbers rather than dimensional physical parameters, so that a i and a −1 i can coexist in the same space and be plotted in the same axes. Then we can take two diagonal matrices, written in the same basis and having either a i or a −1 i at diagonal, stretch one matrix by a factor of √ n + 1 and compress the other one by the same factor along the same direction, the spatial diagonal → e d = (1, 1, . . . , 1) T / √ n, and obtain the matrix that produces the sought nonzero roots for the original physical problem and its inverse, with their eigenvalues related reciprocally. This duality and the possibility to obtain simultaneously both the sought roots and their inverses will prove very useful for geometrization that we introduce in the following section.

The Geometrization
A geometric image of a real symmetric positive definite n × n matrix is an n-dimensional hyper-ellipsoid in real n-space obtained by constructing a quadratic form from the matrix. For simplicity we shall omit "hyper" and refer to them simply as ellipsoids. For a diagonal matrix D = diag{a 1 , a ,2 , . . . , a n } we construct the form and obtain an ellipsoid with semiaxes given by inverse square roots of the diagonal elements, aligned with the Cartesian axes of the n-space. For a general real positive definite matrix, we obtain a similar ellipsoid, but rotated and with semiaxes built from eigenvalues as p i = 1/ √ λ i . The same is true for a complex Hermitian positive definite matrix, although it may be not straightforward to visualize it. Still some unitary transform, i.e., pure rotation, converts it into a real diagonal matrix without modifying its spectrum, and the ellipsoid visualization remains valid.
Therefore, the transforms that simultaneously reconstruct the matrices M containing all the nonzero roots of the original problem (64) and its inverse W (66) from diagonal matrices D = diag{a 1 , a 2 , . . . , a n } and D −1 = diag{1/a 1 , 1/a 2 , . . . , 1/a n } correspond to scaling the ellipsoids with semiaxes 1/ √ a i and √ a i , respectively, to produce different ellipsoids with semiaxes 1/ √ λ i and √ λ i , respectively. While formally we have to stretch/compress a complex Hermitian matrix rather than real symmetric one, the transform actually proceeds along the only explicitly real basis vector → e d = (1, 1, . . . , 1) T / √ n and does not change anything in its complementary (n − 1)-dimensional space spanned by the complex basis vectors. If needed, an additional unitary transform can be performed in this complementary space to obtain a real basis and a real symmetric matrix with the same spectrum, and thus explicitly transform an ellipsoid to an ellipsoid. This is the "shuffling back" transform (19) that gave us the real symmetric matrix out of the complex zero-split one, but in doing so ruined the symmetry that was later used to diagonalize the matrix. We now see that this additional transform, confined to the complementary space, is immaterial, and the net result is geometrically stretching/compressing an n-dimensional ellipsoid by a factor of √ n + 1 along the direction → e d = (1, 1, . . . , 1) T / √ n. This is almost what we strived to achieve, the only nuisance being the recurring reciprocity that interferes with intuition. It is natural to work with matrix M corresponding to the original problem, and think in terms of its representing ellipsoid. However, the semiaxes of the starting ellipsoid for M are the inverse square roots of the corresponding parameters a i , the sought roots λ i are the inverse squared semiaxes of the resulting ellipsoid, and the stretching of matrix is geometrically the compressing of its representing ellipsoid. However, all this is taken care of if we consider the dual transform of the inverse matrix D −1 : the semiaxes of the starting ellipsoid are now √ a i , the squared semiaxes of the resulting ellipsoid are now λ i , and the compressing of the matrix implied by transform (66) is the spatial stretching of the ellipsoid.
We have thus constructed a simple geometric visualization procedure for the roots of the zero-split matrix (64): given n spin-1 /2 nuclei with arbitrary positive coupling constants a i , take an n-dimensional ellipsoid centered at the origin of the reference frame and having semiaxes √ a i aligned with the n Cartesian axes, stretch it by a factor of √ n + 1 along the spatial diagonal → e d = (1, 1, . . . , 1) T / √ n, and the semiaxes of the obtained new ellipsoid q i produce the n sought positive roots as λ i = q 2 i , i = 1, . . . , n. Augment this set with additional zero root λ 0 = q 2 0 = 0, and obtain the n + 1 energy levels of the penultimate Hamiltonian block as Figure 2 shows graphically how stretching a 2D ellipse by a factor of √ 3 along the (1, 1) T / √ 2 direction produces a transformed and rotated ellipse with larger semiaxes corresponding to the sought roots for a system with two nuclei, illustrating the solutions λ 1,2 = (a 1 + a 2 ) ± a 2 1 + a 2 2 − a 1 a 2 to equation λ 2 − 2(a 1 + a 2 ) + 3a 1 a 2 = 0.  Figure 2 shows graphically how stretching a 2D ellipse by a factor of 3 along the      Stretching the representing ellipse for a two-nuclei system with a 1 = 25 and a 2 = 9 relative units, corresponding to semiaxes of 5 and 3 r.u., along the direction (1, 1) T / √ 2 by a factor of √ 3 to obtain the roots as the squared semiaxes of the transformed ellipse. The figure shows the initial and resulting ellipses with their semiaxes as bold sticks, the stretch lines along the (1, 1) T / √ 2 direction and the reference line normal to them, the old and new Cartesian axes for the ellipses, and the correspondence of the points at the cross-section of the initial ellipse with the stretch lines and their images at the final ellipse that are a factor of √ 3 farther away from the reference line along the stretch lines. Figure 3 shows the operation of the developed transform for a three-nuclei system, when the representing ellipsoid is a conventional 3D ellipsoid that can still be relatively easily drawn and visualized. The figure shows four panels for the representative cases of the starting ellipsoid having the shape of a sphere (all three semiaxes are equal, a = b = c), oblong ellipsoid of revolution (one longer semiaxis and two equal ones, a > b = c), oblate ellipsoid of revolution (one shorter semiaxis and two equal ones, a = b < c), and a general ellipsoid (all three semiaxes are different). The starting ellipsoid with its semiaxes aligned with Cartesian axes is stretched by a factor of 2 along the spatial diagonal (1, 1, 1) T / √ 3. The images were generated in Mathematica by parametrizing the initial ellipsoid as applying the transform to the radius-vector as (79) and using the ParametricPlot3D function to draw the initial and the transformed ellipsoids shown in the figure.
and using the ParametricPlot3D function to draw the initial and the transformed ellipsoids shown in the figure.
When building the transform (79) to generate the figures we took advantage of the freedom to choose the basis vectors in the orthogonal complement to the stretching direction mentioned above, and used a pair of real vectors instead of complex ones generated by cubic roots of 1 to embed the figure in ordinary 3D space. A careful study of Figures 2 and 3 illustrates the general spectral bounds of (12) and (13): the stretching either leaves some semiaxes unchanged for the degenerate case of the ellipsoids of revolution, of increases them for differing semiaxes, but so as to not step over the next larger one, and the longest semiaxis of the initial ellipsoid always becomes longer. The first one would be n equivalent nuclei with the same coupling constant . a The Figure 3. The transform for several representative cases of a three-nuclei system: stretching a 3D ellipsoid a factor of 2 along the spatial diagonal (1, 1, 1) T / √ 3. Legends in the panels indicate semiaxes of the starting ellipsoid in r.u. that are equal to square roots of the original hyperfine coupling constants a i , while squared semiaxes q i of the stretched ellipsoids give the sought energies When building the transform (79) to generate the figures we took advantage of the freedom to choose the basis vectors in the orthogonal complement to the stretching direction mentioned above, and used a pair of real vectors instead of complex ones generated by cubic roots of 1 to embed the figure in ordinary 3D space. A careful study of Figures 2 and 3 illustrates the general spectral bounds of (12) and (13): the stretching either leaves some semiaxes unchanged for the degenerate case of the ellipsoids of revolution, of increases them for differing semiaxes, but so as to not step over the next larger one, and the longest semiaxis of the initial ellipsoid always becomes longer.
A couple of simple examples can also be provided for systems of higher dimensions. The first one would be n equivalent nuclei with the same coupling constant a. The ellipsoid then turns into an n-dimensional sphere of radius √ a, and stretching it in any direction by a factor of √ n + 1 produces an ellipsoid of revolution with one semiaxis equal to (n + 1)a and remaining n − 1 semiaxes unchanged at √ a, which gives the solution (29). The second example would be n − 1 equivalent nuclei with the coupling constants a 1 = · · · = a n−1 = a and one nucleus with a larger coupling a n > a. The starting ellipsoid for this configuration is an ellipsoid of revolution with one long semiaxis √ a n and the remaining n − 1 equal semiaxes √ a, which we have to stretch along direction → e d = (1, 1, . . . , 1) T / √ n. The stretching will be confined to plane σ dn spanning vectors → e d and → e n = (0, 0, . . . , 0, 1) T , with the latter corresponding to the long semiaxis. The stretched ellipsoid will have n − 2 semiaxes normal to σ dn staying at √ a that produce the (n − 2)−degenerate root λ 1 = a. The restriction of the ellipsoid to plane σ dn is an ellipse with semiaxes √ a and √ a n that we now have to stretch by a factor of √ n + 1 along the direction with angle γ to axis → e n = (0, 0, . . . , 0, 1) T given by cos γ = → e n · → e n = 1 √ n , sin γ = n − 1 n .
In these simples cases the solutions provided by the geometric approach are certainly more complex than the straightforward algebraic solutions to quadratic equations. For even slightly more complex systems both approaches fail to produce closed-form solutions. This is that heaping up complexity behind the simple to imagine or draw things referred to in the Introduction. However, while for the algebraic approach this is a plain dead end, the geometric visualization still remains valid, and it is earnestly believed, useful. simply by the coupling constants a i . Finding eigenvalues of a matrix can be viewed as going from the initial basis to the basis of its eigenvectors, and, although this can be, and numerically almost always is, done directly in one step, the ability to decompose it into a sequence of more comprehendible and visualizable steps by going through several intermediate bases can provide new insights about the general properties of the system. In our case the initial matrix was first shifted, uniformly scaled and changed in sign to obtain the Laplacian-like matrix (5) that is guaranteed to have one zero root and all remaining roots positive, which can be seen as a hyper-paraboloid, a higher-dimensional generalization of the shape of a wine glass. After such conditioning of the initial matrix to generate the starting representation for geometrization, all subsequent steps can be viewed as manipulations with geometric shapes. First the zero splitting transform was applied to obtain a smaller positive definite matrix (37), which can be seen as taking a cross-section of the hyperboloid normal to its axis, direction (1, 1, . . . , 1) T n+1 in (n + 1)−dimensional space, and obtaining a hyper-ellipsoid one step lower in dimension, similar to looking at the glass from the top to see a circular pattern. Then a stretching transform along direction (1, 1, . . . , 1) T n in the new space (38) produced a simpler hyper-ellipsoid, which could be rotated (62) to align it with Cartesian axes of its space, i.e., to diagonalize the corresponding matrix. This sequence of transforms is not unitary all way through and it does modify the spectrum, as it involves a non-unitary stretching transform in the middle, but it still can be seen as two legs of unitary rotations connected by an easily visualizable non-unitary bridge of the stretch, in which all the algebraic complexity of the problem is condensed. Such "click algebra" of going in simple steps between consecutive bases is believed to be useful for analyzing matrix-related problems common in science.