Mapping spins to fermions via the Jordan–Wigner (JW) transformation can render mean-field (Hartree–Fock, HF) descriptions effective for strongly correlated spin systems. As established in recent work, the application of such approaches is not limited by the nonlocal structure of JW strings or by site ordering because string operators can be absorbed into Thouless rotations of a Slater determinant, and the variational optimization of a unitary Lie-algebraic similarity transformation removes any ordering dependence. Leveraging these ideas, we develop a self-consistent field (SCF) scheme that expresses the mean-field energy as a functional of the single-particle density matrix, providing an alternative to gradient-based optimization of Thouless parameters. We derive the analytical orbital Hessian to diagnose HF stability and compute the ground-state correlation energy through the random-phase approximation (RPA). Benchmark results for the XXZ and J 1 –J 2 model on one- and two-dimensional lattices demonstrate that RPA significantly improves mean-field accuracy. 24 4 2026 22 9 4367 4367–4378 16 5 2026 © 2026 The Authors. Published by American Chemical Society This article is licensed under CC-BY 4.0 Abstract Mapping spins to fermions via the Jordan–Wigner (JW) transformation can render mean-field (Hartree–Fock, HF) descriptions effective for strongly correlated spin systems. As established in recent work, the application of such approaches is not limited by the nonlocal structure of JW strings or by site ordering because string operators can be absorbed into Thouless rotations of a Slater determinant, and the variational optimization of a unitary Lie-algebraic similarity transformation removes any ordering dependence. Although the use of spin and point-group symmetries to factor the Hamiltonian into invariant subspaces − extends the reach of exact diagonalization (ED), the exponential We express the mean-field energy as a functional of the single-particle density matrix ρ , interpret JW strings as Thouless rotations, and remove site-ordering dependence via unitary LAST (combined with HF, this constitutes the orbital-optimized uLAST method, oo-uLAST, see Section ). 7 ∑ m < n J m n s m z s n z = E const + ∑ k l t k l c k † c l + 1 2 ∑ k l m n [ k n | l m ] c k † c l † c m c n Thus, a Fock matrix F z for z -coupling can be constructed from the single-particle density matrix, ρ kl ≡ ⟨Φ| c l † c k |Φ⟩, in the usual way (see eq ), 8 F z = t + Γ with Γ defined in eq : 9 Γ k l = ∑ m n [ k l | m n ] ρ n m The respective contribution to the mean-field energy is given in eq : 10 E z [ ρ ] = Tr ( t ρ ) + 1 2 Tr ( Γ ρ ) + E const In contrast, xy -coupling involves string operators (strings vanish in open chains with nearest-neighbor interactions; in rings, a remaining string for the coupling between the first and last sites reduces to a sign factor in a definite fermion-number sector). The unitary orbital rotation R ∈ C N orb × N orb effected by a general string e i α q n q is given in eq , where C ∈ C N orb × N f and C R collect the orthonormal occupied orbitals defining the Slater determinants |Φ⟩ and |Φ R ⟩ ≡ e i ∑ q α q n q , respectively; N orb is the size of the single-particle basis, and N f is the fermion number. 12 C R = RC = ( e i α 1 0 0 0 ⋱ 0 0 0 e i α N orb ) C With the aim of deriving the xy -coupling contribution to the Fock operator, F xy , we now formulate the respective mean-field energy as a functional of the single-particle density matrix, ρ = CC † . OOO is connected to our MATLAB code through a small C++/MEX interface layer that implements a callback: given the current orbital coefficients and occupations proposed by OOO , the interface constructs the (generally complex) density matrix ρ = CC † and calls a MATLAB routine to evaluate the HF energy, E [ ρ ], and the Fock matrix, F [ ρ ]. These quantities are returned to OOO , which performs the diagonalization and updates the orbitals using its built-in acceleration strategy. This cycle is repeated until self-consistency is reached. For robustness, we first perform a short preconditioning stage with our original MATLAB SCF loop using simple density damping and then start OOO from the corresponding Fock matrix. 2.3. Gauge Freedom in oo-uLAST We note a redundancy in the parametrization of the real symmetric matrix Θ that defines the uLAST correlator γ (cf. eq ). Consider the separable shift of eq , 30 Θ p q → Θ p q ′ = Θ p q + χ p + χ q , χ p ∈ R which changes the generator Δγ = γ­( Θ ′ ) – γ­( Θ ), 31 Δ γ = i 2 ∑ p < q ( χ p + χ q ) n p n q = i 2 ∑ p χ p n p ∑ q ≠ p n q = i 2 ∑ p χ p n p ( N − 1 ) where N ≡ ∑ q n q . Overall, PBC and more compact clusters generally show larger residual errors, consistent with the increased correlation demands for a larger number of closed loops. The spikes and kinks in the RPA error curves have a well-defined physical origin: They arise at parameter values where competing Hartree–Fock solutions become nearly degenerate, and the lowest-energy solution undergoes a qualitative change. This is directly analogous to level crossings in electronic-structure theory. Such features correspond loosely to the existence of multiple phases in the thermodynamic limit and are a well-known consequence of finite-size mean-field calculations near phase boundaries. We note that the orbital Hessian derived here also provides the building blocks for the coupled-perturbed Hartree–Fock equations, enabling the computation of response properties in future work. Furthermore, while RPA already provides a substantial and computationally inexpensive improvement over the mean-field description, it represents a first rung on the correlation ladder. Exploring more sophisticated correlation methods in the fermionic JW/EJW framebearing in mind that the nonlocal string structure of the Hamiltonian may require adaptations beyond standard coupled-cluster formulationsis a natural direction for future work. Acknowledgments S.G.T.