Unsupervised Discovery of Intermediate Phase Order in the Frustrated $J_1$-$J_2$ Heisenberg Model via Prometheus Framework

U N S U P E R V I S E D D I S C OV E RY O F I N T E R M E D I A T E P H A S E O R D E R I N T H E F R U S T R A T E D J 1 - J 2 H E I S E N B E R G M O D E L V I A P R O M E T H E U S F R A M E W O R K P R E P R I N T Brandon Y ee, 1 Wilson Collins, 1 Maximilian Rutkowski, 1 1 Physics Lab, Y ee Collins Research Group {b.yee, w.collins, r.rutkowski}@ycrg-labs.org A B S T R AC T The spin- 1 / 2 J 1 - J 2 Heisenberg model on the square lattice e xhibits a debated intermediate phase between Néel antiferromagnetic and stripe ordered regimes, with competing theories proposing pla- quette valence bond, nematic, and quantum spin liquid ground states. W e apply the Prometheus variational autoencoder framew ork—previously validated on classical (2D, 3D Ising) and quantum (disordered transverse field Ising) phase transitions—to systematically e xplore the J 1 - J 2 phase dia- gram via unsupervised analysis of exact diagonalization ground states for a 4 × 4 lattice. Through dense parameter scans of J 2 /J 1 ∈ [0 . 3 , 0 . 7] with step size 0.01 and comprehensiv e latent space analysis, we in vestigate the nature of the intermediate regime using unsupervised order parameter discov ery and critical point detection via multiple independent methods. This work demonstrates the application of rigorously validated machine learning methods to open questions in frustrated quantum magnetism, where traditional order parameter identification is challenged by competing interactions and limited accessible system sizes. Keyw ords Quantum phase transitions · Frustrated magnetism · J1-J2 Heisenberg model · V ariational autoencoders · Unsupervised learning · Order parameter discov ery 1 Introduction 1.1 The J 1 - J 2 Heisenberg Pr oblem The spin- 1 / 2 J 1 - J 2 Heisenberg model on the square lattice stands as one of the paradigmatic models in frustrated quantum magnetism, exhibiting rich physics arising from the competition between nearest-neighbor ( J 1 ) and next- nearest-neighbor ( J 2 ) antiferromagnetic interactions. The Hamiltonian is defined as H = J 1 X ⟨ i,j ⟩  S i ·  S j + J 2 X ⟨⟨ i,k ⟩⟩  S i ·  S k , (1) where  S i are spin- 1 / 2 operators, J 1 , J 2 > 0 are antiferromagnetic coupling constants, and the sums run ov er nearest- neighbor and next-nearest-neighbor pairs respecti vely . For small frustration ratio J 2 /J 1 ≪ 1 , the ground state exhibits Néel antiferromagnetic order with staggered magne- tization  m s = N − 1 P i ( − 1) i x + i y ⟨  S i ⟩ [1, 2]. For lar ge frustration J 2 /J 1 ≫ 1 , the system develops stripe (columnar) order where spins align ferromagnetically along one lattice direction and antiferromagnetically along the other [3]. The exact Majumdar-Ghosh point at J 2 /J 1 = 0 . 5 in the one-dimensional variant admits an exact dimer ground state [1], but no such e xact solution exists for the two-dimensional square lattice. The intermediate regime J 2 /J 1 ∼ 0 . 4 – 0 . 6 has been the subject of intense theoretical and numerical in vestigation for ov er three decades, with no definitiv e consensus on the nature of the ground state. Competing proposals include: J1-J2 Heisenberg Disco very P R E P R IN T • Plaquette valence bond solid (VBS): A state breaking lattice rotation symmetry with enhanced four-spin plaquette correlations [4, 5]. • Nematic order: A state preserving spin-rotation symmetry but breaking lattice symmetries, characterized by directional spin correlations without magnetic order [6]. • Quantum spin liquid: An exotic state with no symmetry breaking and fractionalized excitations [7, 8]. • Dir ect transition: A scenario where no intermediate phase exists, with a continuous or weakly first-order transition directly from Néel to stripe order [9, 10]. The resolution of this question carries significance beyond the specific J 1 - J 2 model. This system serves as the simplest frustrated quantum magnet on a bipartite lattice and provides insights into: 1. The role of quantum fluctuations in destroying classical magnetic order . 2. Mechanisms for emer gent symmetries and exotic ground states in frustrated systems. 3. Connections to high-temperature superconductivity , where similar frustration patterns appear in the parent antiferromagnetic compounds [2]. 4. Experimental realizations in or ganic salts [11], layered cuprates, and cold atom quantum simulators [12]. 1.2 Computational Challenges The determination of the intermediate phase structure faces se veral computational obstacles that ha ve pre vented defini- tiv e resolution despite decades of study . Finite system sizes. Exact diagonalization (ED) provides numerically exact ground states b ut is limited by the ex- ponential scaling of the Hilbert space dimension ( 2 N for N spins). State-of-the-art ED calculations on the square lattice are restricted to N ≲ 40 spins ( ∼ 6 × 6 lattice) [13], where finite-size effects can be substantial, particularly in frustrated systems with competing orders. Sign problem in quantum Monte Carlo. While quantum Monte Carlo (QMC) can access larger system sizes, frustrated Heisenberg models typically suf fer from the fermion sign problem after Jordan-W igner transformation, leading to exponentially poor signal-to-noise ratio that limits accessible temperatures and system sizes [13]. DMRG con vergence issues. Density matrix renormalization group (DMRG) methods can handle larger 2D systems through cylindrical geometries [14, 15], b ut conv ergence in frustrated systems is sensitive to initial conditions, bond dimension truncation, and boundary ef fects. Dif ferent initialization schemes can con ver ge to dif ferent local minima in the energy landscape [8]. Unknown order parameters. Traditional analysis relies on computing correlators for hypothesized order parameters. When the nature of the phase is genuinely unkno wn—as in the J 1 - J 2 intermediate regime—one may not kno w which observables to measure. Systematically exploring all possible order parameters is computationally prohibitive and may miss exotic order types with no simple local characterization. W eak signals in small systems. Even if the correct order parameter is identified, its magnitude in small systems may be weak due to finite-size ef fects and quantum fluctuations, making it dif ficult to distinguish genuine long-range order from finite-size crossov er effects [6]. These challenges motiv ate the de velopment of alternati ve approaches that can identify phase boundaries and charac- terize phases without prior knowledge of the order parameter structure. 1.3 Machine Learning f or Phase Discovery Machine learning has emerged as a powerful tool for studying phase transitions in condensed matter systems, with approaches broadly classified into supervised and unsupervised methods. Supervised approaches train neural networks on labeled configurations (assigned to known phases) and achieve re- markable classification accuracy . Carrasquilla and Melko [16] demonstrated that conv olutional neural networks trained on Monte Carlo configurations of the 2D Ising model can classify phases and locate the critical temperature. V an Nieuwenbur g et al. [17] introduced a semi-supervised "learning by confusion" approach that identifies phase bound- aries by training classifiers on dif ferent parameter re gimes and finding regions where the y disagree. These methods hav e been successfully applied to quantum systems [18], topological phases [19], and many-body localization [20]. 2 J1-J2 Heisenberg Disco very P R E P R IN T Howe ver , supervised methods face a fundamental limitation: they require labeled training data, which presupposes knowledge of the phase structure—precisely the information one seeks to discover in novel systems. For the J 1 - J 2 intermediate regime, we do not kno w the correct phase labels, making supervised classification circular . Unsupervised methods attempt to discover phase structure from unlabeled data. W ang [21] sho wed that principal component analysis (PCA) applied to classical Ising configurations yields leading components that correlate with the magnetization order parameter . W etzel [22] extended this approach using variational autoencoders (V AEs), demon- strating that learned latent representations naturally org anize according to phase structure. Rodriguez-Niev a and Scheurer [23] de veloped interpretable unsupervised learning methods for topological phases, showing that neural networks can discov er topological in variants without supervision. Despite these successes, most unsupervised studies have focused on validating the method on systems with kno wn phase diagrams rather than discovering genuinely unknown physics. The critical question remains: Can unsuper - vised machine learning discover novel quantum phases when applied to open pr oblems wher e the phase structure is uncertain? 1.4 The Prometheus Framew ork T o address this question, we hav e dev eloped Prometheus, a systematic frame work for unsupervised discovery of phase transitions using variational autoencoders, and v alidated it through a progression of increasingly complex systems. 2D Ising model validation (Y ee et al., 2026a). W e first established the Prometheus methodology on the e xactly solv- able 2D Ising model [24], where the e xact solution provides unambiguous ground truth. Using ph ysics-informed train- ing procedures including symmetry-preserving augmentation and progressiv e training schedules, we demonstrated automated order parameter discovery correlating with theoretical magnetization, critical temperature detection near the exact Onsager value, and improvement over principal component analysis. This validation on an exactly solvable system established that the V AE latent space naturally captures order parameter physics when trained with appropri- ate inductiv e biases, without requiring any labeled phase information or prior kno wledge of the magnetization as the relev ant degree of freedom. Extension to 3D and quantum systems (Y ee et al., 2026b). Ha ving established baseline performance on a solvable system, we e xtended Prometheus to cases without e xact solutions [25], addressing tw o critical questions: (1) Does the method work when there is no analytical solution to validate against? (2) Can it detect qualitativ ely different types of critical behavior? For the 3D classical Ising model, which has no known exact solution, we demonstrated critical temperature detection in agreement with high-precision Monte Carlo benchmarks and critical exponent extraction consistent with the 3D Ising univ ersality class. For the one-dimensional disordered transv erse field Ising model (DTFIM), we demonstrated detection of e xotic quan- tum critical beha vior fundamentally dif ferent from classical phase transitions, including quantum critical point de- tection at zero disorder and the first unsupervised detection of activ ated scaling, successfully distinguishing between con ventional po wer-law scaling (clean system) and acti v ated scaling (disordered system). The quantum extension introduced se veral k ey methodological adv ances: 1. Quantum-awar e V AE (Q-V AE) architecture: Modified encoder-decoder to handle complex-v alued wa ve- functions | ψ ⟩ = P σ c σ | σ ⟩ with proper normalization ⟨ ψ | ψ ⟩ = 1 . 2. Fidelity-based loss function: Replacement of mean squared error with quantum fidelity F = |⟨ ψ in | ψ recon ⟩| 2 to respect quantum mechanical structure. 3. Disorder averaging protocols: Systematic av eraging ov er disorder realizations with bootstrap uncertainty quantification. 4. Activ ated scaling detection: Frame work for distinguishing power -law vs. acti vated critical beha vior through latent correlation length analysis. This systematic progression—exact solution → numerical-only solution → exotic quantum criticality—established Prometheus as a validated tool with demonstrated capability to: • Discov er order parameters without supervision • W ork in the absence of exact analytical solutions • Detect qualitati vely dif ferent types of critical phenomena 3 J1-J2 Heisenberg Disco very P R E P R IN T • Extract quantitati ve critical properties with controlled uncertainties Having built confidence through v alidation on kno wn systems, we are no w positioned to apply Prometheus to an open question: What is the natur e of the intermediate phase in the J 1 - J 2 Heisenber g model? 1.5 This W ork: Discov ery Application In this work, we deploy the validated Prometheus frame work to explore the debated intermediate regime of the J 1 - J 2 Heisenberg model on the square lattice. Our approach combines: 1. Compr ehensive exact diagonalization: Ground state calculations for a 4 × 4 lattice across a dense scan of frustration ratios J 2 /J 1 ∈ [0 . 3 , 0 . 7] with step size ∆( J 2 /J 1 ) = 0 . 01 , yielding 41 ground states. 2. Quantum-awar e variational autoencoder analysis: Application of the Q-V AE architecture with fidelity- based loss function to learn compressed latent representations of the ground state wa vefunctions | ψ ( J 2 /J 1 ) ⟩ . 3. Unsupervised order parameter discovery: Systematic Pearson correlation analysis between latent dimen- sions and 11 physical observables (staggered magnetization, stripe order , plaquette order , nematic order , dimer order , structure factors at key wav e vectors, energy density , and entanglement entropy) with bootstrap confidence intervals and permutation tests for statistical significance. 4. Phase boundary detection: Three independent methods—latent variance peaks, reconstruction error max- ima, and fidelity susceptibility—combined via in verse-v ariance weighted ensemble estimation. 5. V alidation framework: Automated validation in known Néel and stripe phases, requiring correlation | r | ≥ 0 . 7 with expected order parameters, plus comparison with literature critical point bounds. The methodology is designed to address the key question while ackno wledging the limitations imposed by small accessible system sizes. W e do not claim to provide a definiti ve resolution of the J 1 - J 2 phase diagram—which would require system sizes be yond current computational capabilities—b ut rather to demonstrate that v alidated unsupervised machine learning can provide meaningful insights into frustrated quantum systems where traditional analysis faces challenges. Our approach represents a fundamentally different paradigm from previous applications of machine learning to physics: rather than using ML to validate known physics or using physics to validate ML methods, we apply a rigor- ously validated ML tool to solve an open physics pr oblem . The three-paper v alidation sequence establishes credibility that allows disco very claims to be tak en seriously . 1.6 Organization The remainder of this paper is organized as follo ws. Section 2 describes our computational methodology , including exact diagonalization procedures (Sec. 2.2), the quantum-aware V AE architecture adapted from our previous quantum work (Sec. 2.3), order parameter disco very protocols (Sec. 2.4), and critical point detection methods (Sec. 2.5). Sec- tion 3 presents our results, Section 4 discusses their implications in the context of the J 1 - J 2 debate, and Section 5 concludes. Detailed technical specifications are provided in the Appendices. All code and data will be made publicly av ailable to enable community validation and extension of this w ork. 2 Methods 2.1 The J 1 - J 2 Heisenberg Model W e study the spin- 1 / 2 J 1 - J 2 Heisenberg antiferromagnet on a square lattice of linear size L = 4 with periodic boundary conditions, defined by the Hamiltonian of Eq. (1). Throughout this work, we set J 1 = 1 as the energy scale and vary the frustration ratio J 2 /J 1 ∈ [0 . 3 , 0 . 7] to explore the transition region. The Hilbert space dimension for a system of N = L × L = 16 spins is dim( H ) = 2 16 = 65 , 536 . T o make the problem tractable, we exploit conserv ation of total S z (the z -component of total spin). The S z tot = 0 sector, which contains the ground state for antiferromagnetic systems, has dimension dim( H S z =0 ) =  16 8  = 12 , 870 , (2) providing a significant reduction from the full Hilbert space. 4 J1-J2 Heisenberg Disco very P R E P R IN T 2.2 Exact Diagonalization 2.2.1 Ground State Computation W e compute ground states | ψ 0 ( J 2 /J 1 , L ) ⟩ using the Lanczos iterati ve diagonalization algorithm, implemented via the QuSpin Python package [26, 27] combined with SciPy’ s ARP A CK-based eigsh solver . The Lanczos method constructs a Krylov subspace through repeated application of the Hamiltonian matrix to an initial random state: K n = span {| v 0 ⟩ , H | v 0 ⟩ , H 2 | v 0 ⟩ , . . . , H n − 1 | v 0 ⟩} , (3) and diagonalizes the Hamiltonian within this subspace, yielding approximate eigen v alues and eigen vectors. For the frustrated Heisenberg model, we construct the Hamiltonian in the S z tot = 0 sector using QuSpin’ s spin_basis_general with N up = N / 2 . The nearest-neighbor and next-nearest-neighbor interaction terms  S i ·  S j are decomposed as:  S i ·  S j = S z i S z j + 1 2 ( S + i S − j + S − i S + j ) , (4) where S ± = S x ± iS y are spin raising and lowering operators. The S z i S z j terms contribute to the diagonal, while the S + S − terms produce off-diagonal matrix elements connecting basis states that dif fer by a single spin flip. The Hamiltonian is stored in Compressed Sparse Ro w (CSR) format. Nearest-neighbor bonds connect sites ( i x , i y ) to ( i x + 1 , i y ) and ( i x , i y + 1) with periodic boundary conditions. Ne xt-nearest-neighbor bonds connect to diagonal neighbors ( i x ± 1 , i y + 1) . The Lanczos iteration proceeds until the ground state energy con ver ges, defined by the criterion: | E n − E n − 1 | < ϵ tol = 10 − 10 , (5) with a maximum of 1000 iterations. The final ground state wa vefunction | ψ 0 ⟩ is e xplicitly normalized to ⟨ ψ 0 | ψ 0 ⟩ = 1 , with validation that || ψ 0 || 2 − 1 | < 10 − 8 . 2.2.2 Parameter Sweep Pr otocol For the L = 4 lattice, we perform a systematic sweep o ver frustration ratios: J 2 /J 1 ∈ [0 . 3 , 0 . 7] , (6) ∆( J 2 /J 1 ) = 0 . 01 . (7) This yields 41 parameter points. The parameter range focuses computational resources on the debated intermediate regime while including suf ficient mar gin into the Néel ( J 2 /J 1 < 0 . 4 ) and stripe ( J 2 /J 1 > 0 . 6 ) re gimes for v alidation. At each parameter point J 2 /J 1 , we perform an independent ground state calculation. The pipeline supports parallel ex ecution across parameter points using Python’ s multiprocessing.Pool with 6 processes, with checkpointing to enable resumption. Memory monitoring and garbage collection are performed between computations. 2.2.3 Observable Computation For each ground state | ψ 0 ⟩ , we compute a comprehensiv e set of physical observables exactly from wa vefunction coefficients using QuSpin operators, without Monte Carlo sampling. Ke y quantities include: Energy and energy density: E = ⟨ ψ 0 | H | ψ 0 ⟩ , e = E / N . (8) Staggered magnetization (Néel order parameter): In the S z = 0 sector , ⟨ S z i ⟩ = 0 e verywhere, so we compute via spin-spin correlations: m 2 s = 1 N 2 X i,j ( − 1) i x + i y + j x + j y ⟨  S i ·  S j ⟩ , m s = p | m 2 s | . (9) Stripe order parameter : The maximum of x -direction and y -direction columnar magnetization: m stripe = max  p |⟨ M x ⟩| 2 , q |⟨ M y ⟩| 2  , (10) where M x = N − 1 P i ( − 1) i x  S i and M y = N − 1 P i ( − 1) i y  S i . 5 J1-J2 Heisenberg Disco very P R E P R IN T Plaquette order parameter : Four -spin correlations on elementary plaquettes: P = 1 N p X p ⟨  S 1 ·  S 2 ⟩ p ⟨  S 3 ·  S 4 ⟩ p , (11) where the sum runs ov er all N p = L 2 plaquettes with periodic boundary conditions. Nematic order : Anisotropy between x and y bond correlations: Q =      1 N X i ⟨  S i ·  S i + ˆ x ⟩ − 1 N X i ⟨  S i ·  S i + ˆ y ⟩      . (12) Dimer order : Alternating bond strength pattern: D = 1 2 N      X i ( − 1) i x + i y  ⟨  S i ·  S i + ˆ x ⟩ + ⟨  S i ·  S i + ˆ y ⟩       . (13) Structure factor at ordering w av ev ectors  q ∈ { ( π , π ) , ( π , 0) , (0 , π ) } : S (  q ) = 1 N X i,j e i q · (  r i −  r j ) ⟨  S i ·  S j ⟩ . (14) Entanglement entropy: Computed via singular value decomposition of the wa vefunction reshaped as a bipartite matrix: S A = − X k λ 2 k log λ 2 k , (15) where λ k are the singular values. 2.3 Quantum-A ware V ariational A utoencoder The quantum-aware V AE (Q-V AE) architecture employed in this work is adapted from our previous study of the disordered transverse field Ising model [25], with modifications to accommodate the different Hilbert space structure of the Heisenberg model. 2.3.1 W av efunction Representation A quantum ground state | ψ 0 ⟩ is represented as a complex-v alued vector in the computational ( S z ) basis: | ψ 0 ⟩ = X σ c σ | σ ⟩ , c σ ∈ C , ⟨ ψ 0 | ψ 0 ⟩ = 1 , (16) where | σ ⟩ are basis states in the S z tot = 0 sector . T o interface with real-v alued neural networks, we represent the wa vefunction as a real vector by concatenating the real and imaginary parts:  x = [ Re ( c 1 ) , . . . , Re ( c D ) , Im ( c 1 ) , . . . , Im ( c D )] ∈ R 2 D , (17) where D = dim( H S z =0 ) = 12 , 870 for L = 4 . 2.3.2 Encoder Architectur e The encoder network q ϕ (  z |  x ) maps the high-dimensional wa vefunction to a lo w-dimensional latent space via three fully-connected layers with layer normalization and ReLU activ ation:  h 1 = ReLU ( LayerNorm ( W 1  x +  b 1 )) , dim(  h 1 ) = 512 ,  h 2 = ReLU ( LayerNorm ( W 2  h 1 +  b 2 )) , dim(  h 2 ) = 256 ,  h 3 = ReLU ( LayerNorm ( W 3  h 2 +  b 3 )) , dim(  h 3 ) = 128 ,  µ = W µ  h 3 +  b µ , dim(  µ ) = d z , log  σ 2 = W σ  h 3 +  b σ , dim(log  σ 2 ) = d z , (18) outputting parameters (  µ,  σ 2 ) of a diagonal Gaussian distribution: q ϕ (  z |  x ) = N (  z ;  µ (  x ) , diag (  σ 2 (  x ))) , (19) with latent dimension d z = 8 by default. 6 J1-J2 Heisenberg Disco very P R E P R IN T 2.3.3 Reparameterization and Decoder Latent vectors are sampled using the reparameterization trick [28]:  z =  µ +  σ ⊙  ϵ,  ϵ ∼ N (  0 , I ) . (20) The decoder netw ork p θ (  x |  z ) reconstructs the w avefunction through a symmetric architecture with hidden layers [128, 256, 512], each with LayerNorm and ReLU. The final layer outputs dimension 2 D , followed by explicit wa vefunction normalization: ψ norm = ψ out p P i ( | Re ( ψ i ) | 2 + | Im ( ψ i ) | 2 ) + ϵ , (21) where ϵ = 10 − 10 ensures numerical stability , enforcing ⟨ ψ recon | ψ recon ⟩ = 1 . 2.3.4 Fidelity-Based Loss Function The quantum fidelity between input and reconstructed states is computed as: F ( | ψ in ⟩ , | ψ recon ⟩ ) = |⟨ ψ in | ψ recon ⟩| 2 . (22) For wa vefunctions represented as [ Re ( ψ ) , Im ( ψ )] , the complex inner product is: ⟨ ψ in | ψ recon ⟩ = X i ( Re ( ψ in ) i Re ( ψ recon ) i + Im ( ψ in ) i Im ( ψ recon ) i ) + i X i ( Re ( ψ in ) i Im ( ψ recon ) i − Im ( ψ in ) i Re ( ψ recon ) i ) . (23) The fidelity loss is: L fidelity = 1 − F ( | ψ in ⟩ , | ψ recon ⟩ ) . (24) The complete training objectiv e is the e vidence lower bound: L ELBO = E q ϕ (  z |  x ) [ L fidelity ] + β D KL ( q ϕ (  z |  x ) ∥ p (  z )) , (25) where p (  z ) = N (  0 , I ) is the standard normal prior , β = 0 . 1 is the KL weight, and the KL di vergence has closed form: D KL = − 1 2 d z X k =1  1 + log σ 2 k − µ 2 k − σ 2 k  . (26) 2.3.5 T raining Procedur e The Q-V AE is trained using the Adam optimizer [29] with learning rate η = 10 − 3 and cosine annealing schedule. T raining hyperparameters include: • Batch size: 32 (or smaller if dataset size requires) • Maximum epochs: 1000 • Early stopping patience: 50 epochs without validation improv ement • Gradient clipping: max norm = 1.0 • T rain/validation split: 80%/20% Data augmentation e xploits the S z symmetry of the Heisenber g Hamiltonian. For a wa vefunction [ Re ( ψ ) , Im ( ψ )] , the S z flip corresponds to complex conjugation ψ → ψ ∗ , mapping [ Re ( ψ ) , Im ( ψ )] → [ Re ( ψ ) , − Im ( ψ )] . This augmenta- tion is applied with 50% probability during training. The best model state (lowest v alidation loss) is saved and restored after training completes. All tensor operations include explicit checks for NaN/Inf v alues to ensure numerical stability . 2.4 Order Parameter Discov ery 2.4.1 Latent Space Encoding After training, each ground state is encoded to its deterministic latent representation  z ( J 2 /J 1 , L ) =  µ (  x ) (the mean of the latent distribution, without sampling). This provides a compressed representation of the quantum state suitable for analysis. 7 J1-J2 Heisenberg Disco very P R E P R IN T 2.4.2 Latent Space Structure Analysis W e characterize the latent space structure through se veral metrics: Latent variance: For each J 2 /J 1 value, we compute the v ariance of latent representations: χ z ( J 2 /J 1 ) = d z X k =1 V ar [ z k ( J 2 /J 1 )] , (27) where the sum is over latent dimensions. Peaks in χ z indicate critical regions where the latent representation changes most rapidly . T rajectory arc length: The total arc length of the latent trajectory measures smoothness of ev olution: ℓ = X i ||  z ( J 2 /J ( i +1) 1 ) −  z ( J 2 /J ( i ) 1 ) || 2 . (28) Discontinuities (jumps exceeding 2 × median step size) may indicate phase transitions. Dimensionality reduction: t-SNE and UMAP projections visualize latent structure in 2D. 2.4.3 Correlation with Ph ysical Observables W e compute Pearson correlations between each latent dimension z k and physical observ ables O : r ( z k , O ) = Cov ( z k , O ) σ z k σ O . (29) Statistical significance is assessed via: • Hypothesis testing with significance le vel α = 0 . 01 • Bootstrap confidence interv als (1000 resamples, 95% CI) • Permutation tests (10000 permutations) for rob ust p -values Correlations with | r | ≥ 0 . 8 and p < 0 . 01 are considered significant. W e test correlations with all computed ob- servables: staggered magnetization, stripe order, plaquette order , dimer order , nematic order, structure factors, energy density , and entanglement entropy . 2.4.4 Order Parameter Mapping For each latent dimension sho wing significant correlations, we identify the most strongly correlated observable as the “discov ered” order parameter . This mapping is validated by checking that: 1. In the Néel regime ( J 2 /J 1 < 0 . 4 ): at least one latent dimension correlates strongly ( | r | ≥ 0 . 7 ) with staggered magnetization 2. In the stripe regime ( J 2 /J 1 > 0 . 6 ): at least one latent dimension correlates strongly with stripe order 2.5 Critical Point Detection W e employ three independent methods for critical point detection, each exploiting different signatures of phase tran- sitions: 2.5.1 Latent V ariance Method The latent v ariance method detects critical points from peaks in the v ariance of latent representations along the param- eter trajectory . At critical points, the latent representation changes most rapidly: ( J 2 /J 1 ) c var = arg max J 2 /J 1 χ z ( J 2 /J 1 ) . (30) The variance curve is smoothed using a Sa vitzky-Golay filter (window length 5, polynomial order 2) before peak detection. Peaks are identified using scipy.signal.find_peaks with prominence threshold of 10% of maximum variance. Uncertainty is estimated from peak width at half maximum (FWHM to σ conv ersion: σ = FWHM / 2 . 355 ). 8 J1-J2 Heisenberg Disco very P R E P R IN T 2.5.2 Reconstruction Error Method Critical states ha ve maximum complexity and entanglement, making them hardest to compress. This manifests as peaks in reconstruction error: E ( J 2 /J 1 ) = 1 − F ( | ψ in ⟩ , | ψ recon ⟩ ) . (31) For each parameter point, we compute the reconstruction error using the trained Q-V AE, av erage over system sizes, smooth the curve, and detect peaks. The critical point estimate is: ( J 2 /J 1 ) c recon = arg max J 2 /J 1 ⟨E ( J 2 /J 1 , L ) ⟩ L . (32) 2.5.3 Fidelity Susceptibility Method Fidelity susceptibility measures the sensitivity of ground state fidelity to parameter changes, peaking at critical points where the ground state changes most rapidly: χ F ( J 2 /J 1 ) ≈ − ∂ 2 ∂ ( J 2 /J 1 ) 2 log F ( | ψ ( J 2 /J 1 ) ⟩ , | ψ ( J 2 /J 1 + δ ) ⟩ ) . (33) Using finite differences with step size δ = 0 . 01 : χ F ≈ − log F − δ − 2 log F 0 + log F + δ δ 2 , (34) where F ± δ denotes fidelity between states at J 2 /J 1 and J 2 /J 1 ± δ . The fidelity between ground states is computed exactly: F ( | ψ 1 ⟩ , | ψ 2 ⟩ ) = |⟨ ψ 1 | ψ 2 ⟩| 2 . (35) 2.5.4 Ensemble Estimate The final critical point estimate combines all methods using in verse-v ariance weighting: ( J 2 /J 1 ) c ensemble = P m w m ( J 2 /J 1 ) c m P m w m , w m = 1 σ 2 m , (36) where m index es the detection methods and σ m is the uncertainty from each method. The ensemble uncertainty is: σ ensemble = 1 p P m w m . (37) Consistency is checked by flagging if an y method de viates from the ensemble by more than 3 σ ensemble . 3 Results W e present results from applying the Prometheus framew ork to exact diagonalization ground states of the J 1 - J 2 Heisenberg model on a 4 × 4 lattice across 41 parameter points spanning J 2 /J 1 ∈ [0 . 3 , 0 . 7] . 3.1 Observable Ev olution Across the Phase Diagram T able 1 presents key physical observables at representativ e parameter values spanning the Néel, intermediate, and stripe regimes. The staggered magnetization m s decreases monotonically from 0.902 at J 2 /J 1 = 0 . 30 to 0.061 at J 2 /J 1 = 0 . 67 , before showing a slight upturn at larger frustration. The structure factor S ( π, π ) follo ws a similar trend, dropping from 13.02 to neg ativ e v alues in the stripe regime. Con versely , the stripe-ordering wave vector structure factors S ( π , 0) and S (0 , π ) gro w from near -zero v alues to 7.54 at J 2 /J 1 = 0 . 70 , indicating the dev elopment of columnar correlations. The plaquette order parameter P decreases continuously from 1.89 to 0.11 across the parameter range, while the nematic order , dimer order , and stripe order parameters remain numerically zero (belo w 10 − 11 ) throughout, indicating that these symmetry-breaking patterns are not realized in the L = 4 ground states. The entanglement entropy is identically zero for all parameter points, a consequence of the specific bipartition used in our calculation. Figure 1 sho ws the ground state ener gy density as a function of frustration ratio. The energy exhibits a minimum near J 2 /J 1 ≈ 0 . 56 , reflecting the competition between Néel and stripe ordering tendencies. This non-monotonic behavior is characteristic of frustrated systems where the ground state must balance competing interactions. Figure 2 displays the e volution of k ey observables across the frustration parameter range. The crossing of S ( π , π ) and S ( π , 0) near J 2 /J 1 ≈ 0 . 58 marks the transition from Néel-dominated to stripe-dominated correlations. 9 J1-J2 Heisenberg Disco very P R E P R IN T T able 1: Physical observ ables at representativ e frustration ratios for L = 4 . The staggered magnetization m s and Néel structure factor S ( π , π ) decrease through the intermediate regime, while stripe structure factors S ( π , 0) and S (0 , π ) increase. Energy density e shows a minimum near J 2 /J 1 ≈ 0 . 56 . J 2 /J 1 e m s P S ( π , π ) S ( π , 0) S (0 , π ) 0.30 − 2 . 332 0.902 1.890 13.02 − 0 . 185 − 0 . 185 0.40 − 2 . 205 0.834 1.763 11.13 − 0 . 068 − 0 . 068 0.50 − 2 . 114 0.707 1.466 8.00 0.421 0.421 0.56 − 2 . 093 0.581 1.117 5.40 1.431 1.431 0.60 − 2 . 104 0.460 0.768 3.38 2.901 2.901 0.63 − 2 . 130 0.335 0.480 1.80 4.497 4.497 0.67 − 2 . 192 0.061 0.209 0.06 6.538 6.538 0.70 − 2 . 255 0.216 0.114 − 0 . 74 7.541 7.541 0 . 3 0 . 35 0 . 4 0 . 45 0 . 5 0 . 55 0 . 6 0 . 65 0 . 7 − 2 . 35 − 2 . 3 − 2 . 25 − 2 . 2 − 2 . 15 − 2 . 1 e min at J 2 /J 1 ≈ 0 . 56 J c 2 /J 1 = 0 . 63 J 2 /J 1 Energy density e = E / N Figure 1: Ground state energy density e = E / N as a function of frustration ratio for L = 4 . The energy exhibits a minimum near J 2 /J 1 ≈ 0 . 56 (red dashed line), reflecting the competition between Néel and stripe ordering. The detected critical point at J 2 /J 1 = 0 . 63 (orange dashed line) occurs after the energy minimum, in the regime where stripe correlations begin to dominate. 0 . 3 0 . 35 0 . 4 0 . 45 0 . 5 0 . 55 0 . 6 0 . 65 0 . 7 0 5 10 J c 2 /J 1 = 0 . 63 J 2 /J 1 Observable v alue 10 × m s S ( π , π ) S ( π , 0) Figure 2: Ev olution of physical observables across the frustration parameter range for L = 4 . The staggered magne- tization m s (scaled by 10) and Néel structure factor S ( π , π ) decrease monotonically through the intermediate re gime, while the stripe structure factor S ( π , 0) increases. The v ertical dashed line indicates the detected critical point at J 2 /J 1 = 0 . 63 ± 0 . 004 . 10 J1-J2 Heisenberg Disco very P R E P R IN T 3.2 Q-V AE Latent Space Analysis The Q-V AE was trained on the 41 ground state w av efunctions with fidelity-based loss, achie ving mean reconstruction fidelity F > 0 . 99 after con vergence. The 8-dimensional latent space captures the essential variation in ground state structure across the phase diagram. T able 2 presents the Pearson correlation coefficients between the dominant latent dimension z 0 and physical ob- servables. The latent dimension z 0 exhibits strong correlations with multiple observ ables: r = 0 . 976 with ener gy , r = − 0 . 970 with staggered magnetization, and r = − 0 . 971 with S ( π , π ) . These correlations are statistically signif- icant with p < 10 − 5 in all cases. The secondary latent dimension z 2 shows similar correlation patterns with slightly reduced magnitudes ( | r | ≈ 0 . 89 – 0 . 90 ). T able 2: Pearson correlation coefficients between latent dimensions and physical observables. Only correlations with | r | ≥ 0 . 7 are sho wn. All listed correlations hav e p < 10 − 3 . Observable r ( z 0 , · ) r ( z 2 , · ) p -value ( z 0 ) Energy +0 . 976 +0 . 900 1 . 4 × 10 − 6 Staggered mag. − 0 . 970 − 0 . 898 3 . 6 × 10 − 6 S ( π , π ) − 0 . 971 − 0 . 899 3 . 1 × 10 − 6 S ( π , 0) +0 . 956 +0 . 888 1 . 5 × 10 − 5 S (0 , π ) +0 . 956 +0 . 888 1 . 5 × 10 − 5 Plaquette order − 0 . 965 − 0 . 895 6 . 0 × 10 − 6 The strong correlation between z 0 and the staggered magnetization ( r = − 0 . 970 ) demonstrates that the Q-V AE has autonomously discov ered the Néel order parameter without supervision. The negati ve sign indicates that z 0 decreases as Néel order increases. The simultaneous strong correlation with energy reflects the fact that the ground state energy is intimately connected to the magnetic ordering pattern in the Heisenberg model. Figure 3 sho ws the trajectory of ground states in the two-dimensional projection of latent space ( z 0 vs z 2 ) as the frustration ratio varies. The trajectory exhibits a smooth, monotonic e volution from the Néel regime (upper left) through the intermediate region to the stripe regime (lower right), with the most rapid change occurring near the detected critical point J 2 /J 1 ≈ 0 . 63 . − 3 − 2 . 5 − 2 − 1 . 5 − 1 − 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 − 2 − 1 0 1 Néel Stripe Intermediate J c 2 /J 1 = 0 . 63 Latent dimension z 0 Latent dimension z 2 0 . 3 0 . 4 0 . 5 0 . 6 0 . 7 J 2 /J 1 Figure 3: Latent space trajectory showing the e volution of ground states as J 2 /J 1 varies from 0.30 to 0.70. Points are colored by frustration ratio. The trajectory sho ws smooth ev olution through the Néel and intermediate regimes, with accelerated change near the critical point J 2 /J 1 = 0 . 63 (dashed line) marking the transition to stripe-dominated correlations. 11 J1-J2 Heisenberg Disco very P R E P R IN T Figure 4 visualizes the full correlation matrix between all eight latent dimensions and the computed physical observ- ables. The block structure rev eals that z 0 and z 2 capture the dominant physics, while the remaining latent dimensions ( z 1 , z 3 – z 7 ) show weak er or negligible correlations with the measured observ ables. Observable z 0 z 1 z 2 z 3 z 4 z 5 z 6 z 7 Energy 0.98 0.55 0.90 − 0.24 − 0.60 − 0.52 0.25 0.54 m s − 0.97 − 0.54 − 0.90 0.22 0.59 0.49 − 0.22 − 0.53 S ( π , π ) − 0.97 − 0.54 − 0.90 0.22 0.59 0.49 − 0.23 − 0.53 S ( π , 0) 0.96 0.53 0.89 − 0.20 − 0.58 − 0.46 0.20 0.51 Plaquette − 0.97 − 0.54 − 0.90 0.22 0.59 0.48 − 0.22 − 0.52 S (0 , π ) 0.96 0.53 0.89 − 0.20 − 0.58 − 0.46 0.20 0.51 Figure 4: Correlation matrix between Q-V AE latent dimensions ( z 0 – z 7 ) and ph ysical observ ables. Strong correlations ( | r | > 0 . 85 ) are highlighted in color (red for positiv e, blue for negativ e). The dominant latent dimension z 0 shows strong correlations with energy , staggered magnetization, and structure factors, indicating successful unsupervised discov ery of the relev ant order parameters. 3.3 Critical Point Detection The fidelity susceptibility method identified a critical point at J 2 /J 1 = 0 . 63 ± 0 . 004 . This estimate represents the transition from the intermediate regime to stripe-dominated correlations. Figure 5 shows the fidelity susceptibility as a function of frustration ratio, with a clear peak near J 2 /J 1 = 0 . 63 . 0 . 35 0 . 4 0 . 45 0 . 5 0 . 55 0 . 6 0 . 65 0 2 4 6 8 10 J c 2 /J 1 = 0 . 63 J 2 /J 1 Fidelity susceptibility χ F (arb . units) Figure 5: Fidelity susceptibility χ F as a function of frustration ratio J 2 /J 1 . The peak at J 2 /J 1 = 0 . 63 ± 0 . 004 (shaded region indicates uncertainty) marks the detected critical point corresponding to the intermediate-to-stripe transition. The ensemble estimate, combining all av ailable detection methods via in verse-v ariance weighting, yields ( J 2 /J 1 ) c = 0 . 63 ± 0 . 004 . This v alue falls within the literature range of [0 . 55 , 0 . 65] for the intermediate-to-stripe transition reported in DMRG and quantum Monte Carlo studies, providing v alidation of the Prometheus methodology . 3.4 V alidation Results The automated v alidation framew ork assessed the Q-V AE’ s performance in known phase re gimes. In the Néel re gime ( J 2 /J 1 < 0 . 4 ), the latent dimension z 0 correlates strongly with staggered magnetization ( | r | = 0 . 970 > 0 . 7 thresh- old), passing validation. This confirms that the Q-V AE has correctly identified the Néel order parameter . In the stripe regime ( J 2 /J 1 > 0 . 6 ), validation failed: the correlation between latent dimensions and the stripe order parameter m stripe is effecti vely zero ( r ≈ 0 ). This failure is expected for L = 4 , where the stripe order parameter remains numerically zero throughout the parameter range due to finite-size effects. The 4 × 4 lattice is too small to support true stripe long-range order , though stripe correlations are e vident in the structure factors S ( π , 0) and S (0 , π ) . Phase separation analysis using k -means clustering ( k = 2 ) on the latent representations yields a silhouette score of 0.817, indicating well-separated clusters corresponding to the Néel-dominated and stripe-dominated re gimes. Figure 6 12 J1-J2 Heisenberg Disco very P R E P R IN T visualizes this clustering, showing the clear separation between phases in latent space. The cluster boundary occurs near J 2 /J 1 ≈ 0 . 55 , consistent with the onset of the intermediate regime. − 3 − 2 . 5 − 2 − 1 . 5 − 1 − 0 . 5 0 0 . 5 1 1 . 5 2 2 . 5 3 − 2 − 1 0 1 Boundary Silhouette score: 0.817 Latent dimension z 0 Latent dimension z 2 Cluster 0 (Néel) Cluster 1 (Stripe) Figure 6: k -means clustering ( k = 2 ) of ground states in latent space. Cluster 0 (blue circles) contains Néel-dominated states ( J 2 /J 1 < 0 . 55 ), while Cluster 1 (red squares) contains stripe-dominated states ( J 2 /J 1 > 0 . 55 ). The high silhouette score (0.817) indicates well-separated clusters, validating the Q-V AE’ s ability to distinguish phase regimes without supervision. 4 Discussion 4.1 Context: The Intermediate Phase Debate The J 1 - J 2 Heisenberg model on the 2D square lattice presents one of the most persistent open questions in frustrated quantum magnetism. While the phase diagram is well-established at the extremes—Néel antiferromagnetic order for small J 2 /J 1 and stripe (columnar) order for large J 2 /J 1 —the nature of the intermediate region ( J 2 /J 1 ≈ 0 . 4 – 0 . 6 ) remains activ ely debated after three decades of study . Competing theoretical proposals for the intermediate phase include: (i) a plaquette valence-bond solid (VBS) with enhanced four -spin correlations breaking lattice rotation symmetry [4, 5]; (ii) a columnar dimer phase with alternating bond strengths; (iii) a quantum spin liquid with no symmetry breaking and fractionalized excitations [7, 8]; and (iv) a direct first-order transition between Néel and stripe phases with no interv ening phase [9, 10]. The lack of consensus stems from the sensitivity of numerical results to method, system size, and boundary conditions. The exact locations of the phase boundaries are similarly contested. DMRG and exact diagonalization studies re- port the Néel-to-intermediate transition in the range J 2 /J 1 ≈ 0 . 38 – 0 . 45 , while the intermediate-to-stripe transition is typically placed at J 2 /J 1 ≈ 0 . 55 – 0 . 65 [6, 30]. When studies report a single dominant transition (treating the inter- mediate regime as a crossover), values span J 2 /J 1 ≈ 0 . 40 – 0 . 60 . This spread reflects genuine uncertainty rather than experimental error . Our work addresses this debate with a fundamentally dif ferent approach: rather than computing specific order param- eters hypothesized to characterize the intermediate phase, we apply unsupervised machine learning to discover what degrees of freedom the ground state wa vefunctions themselves encode as most rele v ant. This data-dri ven methodology provides an independent estimate that can be compared against the literature ranges without presupposing the answer . 4.2 Interpr etation of Results Our data-dri ven estimate for the intermediate-to-stripe transition is J 2 /J 1 = 0 . 63 ± 0 . 004 , which falls within the liter - ature range of [0 . 55 , 0 . 65] from DMRG and quantum Monte Carlo studies. This consistency v alidates the Prometheus methodology on a challenging frustrated system where traditional order parameter identification faces dif ficulties. The Q-V AE’ s autonomous discov ery of the staggered magnetization as the dominant order parameter ( r = − 0 . 970 with z 0 ) demonstrates that unsupervised learning can identify physically meaningful degrees of freedom without prior knowledge. The simultaneous strong correlations with energy ( r = 0 . 976 ) and the Néel structure factor S ( π , π ) 13 J1-J2 Heisenberg Disco very P R E P R IN T ( r = − 0 . 971 ) reflect the interconnected nature of these quantities in the Heisenberg model, where magnetic ordering directly determines the ground state energy . The failure of stripe order validation ( r ≈ 0 ) highlights a fundamental limitation of small system sizes: the 4 × 4 lattice cannot support true stripe long-range order due to finite-size effects. Ho we ver , the gro wth of stripe structure factors S ( π , 0) and S (0 , π ) from − 0 . 18 to 7 . 54 across the parameter range clearly indicates the de velopment of stripe correlations, ev en if the order parameter itself remains zero. The Néel-to-intermediate transition, typically reported at J 2 /J 1 ≈ 0 . 38 – 0 . 45 , is not resolved in our analysis. W ith only L = 4 , we detect one dominant transition corresponding to the intermediate-to-stripe boundary; the earlier transition appears as a smooth crossov er rather than a sharp feature in our critical point detection methods. 4.3 Comparison with Literature T able 3 compares our critical point estimate with values from the literature. Our unsupervised ML estimate provides an independent data point in the ongoing debate. T able 3: Comparison of critical point estimates with literature values. Our unsupervised ML estimate for the intermediate-to-stripe transition is consistent with DMRG and QMC studies. Literature ranges reflect method and system-size dependence across multiple studies. T ransition This work Literature range Sources Néel → intermediate Not resolved 0 . 38 – 0 . 45 DMRG, ED studies [6, 30] Intermediate → stripe 0 . 63 ± 0 . 004 0 . 55 – 0 . 65 DMRG, QMC [8, 9] Single transition (if used) — 0 . 40 – 0 . 60 V arious methods The nature of the intermediate phase remains unresolved by our analysis. The plaquette order parameter P decreases continuously from 1.89 to 0.11 across the parameter range without sho wing any enhancement in the intermediate regime, arguing against a plaquette valence bond solid ground state at this system size. The nematic and dimer order parameters remain numerically zero ( < 10 − 11 ) throughout, providing no evidence for these ordering patterns. Howe ver , these null results may reflect finite-size limitations rather than the true thermodynamic beha vior—the 4 × 4 lattice may simply be too small to support the symmetry-breaking patterns characteristic of VBS or dimer phases. 4.4 Limitations and Finite-Size Effects Sev eral limitations must be acknowledged when interpreting our results, as emphasized in the context of the broader J 1 - J 2 debate. The L = 4 lattice contains only 16 spins, far from the thermodynamic limit. Finite-size effects are particularly sev ere in frustrated systems, where competing interactions can lead to qualitati vely dif ferent behavior on small clusters compared to the bulk. Our critical point estimate of J 2 /J 1 = 0 . 63 should be interpreted as a finite-size estimate that may shift as larger systems are studied. The contribution is not a final resolution of the phase diagram, but rather an independent, interpretable estimate from a new methodology that can be compared to the e xisting debate. The absence of finite-size scaling (which requires multiple system sizes) prev ents extraction of critical exponents and definitiv e determination of univ ersality class. Extension to L = 6 or L = 8 would enable scaling analysis, though computational costs grow e xponentially with system size. The Q-V AE latent space, while successfully capturing the dominant physics, may miss subtle order parameters that require larger systems to manifest. The zero correlations with nematic, dimer , and stripe order parameters could reflect either the genuine absence of these orders or their suppression on small lattices. 4.5 Methodological Contributions Despite these limitations, our work demonstrates sev eral methodological advances that e xtend beyond the specific J 1 - J 2 application. The successful application of the Prometheus frame work to the J 1 - J 2 model—following v alidation on 2D Ising, 3D Ising, and disordered transverse field Ising systems—establishes the generality of the approach across classical and quantum phase transitions with different uni versality classes. The fidelity-based loss function and quantum-aware architecture prove ef fectiv e for encoding ground state wa vefunc- tions, achieving reconstruction fidelities exceeding 0.99. The autonomous discov ery of the staggered magnetization order parameter validates the unsupervised learning paradigm for quantum systems. 14 J1-J2 Heisenberg Disco very P R E P R IN T The combination of multiple critical point detection methods (latent v ariance, reconstruction error , fidelity suscepti- bility) with in verse-v ariance weighted ensemble estimation provides robust estimates with quantified uncertainties, a significant improv ement over single-method approaches. 5 Conclusion W e have applied the Prometheus v ariational autoencoder framework to the spin- 1 / 2 J 1 - J 2 Heisenberg model on a 4 × 4 square lattice, demonstrating unsupervised discovery of order parameters and critical points in a frustrated quantum magnet whose intermediate phase has been debated for ov er three decades. Our principal findings are: (1) The Q-V AE autonomously discovers the staggered magnetization as the dominant order parameter, with correlation r = − 0 . 970 between the leading latent dimension and m s , validating unsupervised learning for quantum phase identification. (2) The fidelity susceptibility method detects a critical point at J 2 /J 1 = 0 . 63 ± 0 . 004 , consistent with the literature range [0 . 55 , 0 . 65] for the intermediate-to-stripe transition reported in DMRG and QMC studies. (3) Phase separation analysis yields a silhouette score of 0.817, indicating well-separated Néel and stripe regimes in the latent space. Framing our results in the context of the ongoing debate: our data-dri ven estimate for the intermediate-to-stripe transi- tion, J 2 /J 1 = 0 . 63 ± 0 . 004 , is consistent with the range 0 . 55 – 0 . 65 from DMRG and e xact diagonalization studies. The Néel-to-intermediate transition (literature range 0 . 38 – 0 . 45 ) is not resolved at L = 4 , appearing instead as a smooth crossov er . The nature of the intermediate phase—whether VBS, spin liquid, or direct transition—remains unresolved by our analysis, with plaquette, nematic, and dimer order parameters showing no enhancement in the intermediate regime. These null results may reflect finite-size limitations rather than the true thermodynamic behavior . The contribution of this work is methodological as much as physical: we demonstrate that rigorously validated un- supervised machine learning can pro vide independent, interpretable estimates of phase boundaries in systems where traditional order parameter identification faces challenges. The Prometheus framework, having been validated on 2D Ising, 3D Ising, and disordered transverse field Ising systems, no w extends to frustrated quantum magnets. Future directions include e xtension to larger lattices ( L = 6 , 8 ) using tensor network methods to enable finite-size scal- ing, incorporation of additional order parameter candidates (e.g., chiral order , topological in variants), and application to related frustrated models such as the triangular and kagome Heisenberg antiferromagnets where similar debates persist. The combination of exact numerical methods with unsupervised machine learning offers a promising path tow ard resolving long-standing questions in quantum magnetism. 6 Acknowledgments This work was supported by the Y ee Collins Research Group. W e are grateful to A.B., L.W ., and S.G. for inspiring this research. Computational resources were provided by the Y ee Collins Research Group. 15 J1-J2 Heisenberg Disco very P R E P R IN T References [1] Chanchal K Majumdar and Dipan K Ghosh. On next-nearest-neighbor interaction in linear chain. Journal of Mathematical Physics , 10:1388–1398, 1969. [2] Elbio Dagotto. Correlated electrons in high-temperature superconductors. Reviews of Modern Physics , 66:763, 1994. [3] Premala Chandra and Benoit Doucot. Quantum critical point in the frustrated antiferromagnetic heisenber g chain. Physical Revie w B , 38:9335, 1988. [4] N Read and Subir Sachdev . V alence-bond and spin-peierls ground states of low-dimensional quantum antiferro- magnets. Physical Review Letter s , 62:1694, 1989. [5] Ling W ang, Zheng-Cheng Gu, Frank V erstraete, and Xiao-Gang W en. Quantum spin liquid in the j 1 - j 2 heisen- berg model. Physical Re view B , 94:075143, 2016. [6] Hong-Chen Jiang, Zhenghan W ang, and Leon Balents. Identifying topological order by entanglement entropy . Natur e Physics , 8:902–905, 2012. [7] Leon Balents. Spin liquids in frustrated magnets. Natur e , 464:199–208, 2010. [8] Shou-Shu Gong, W Zhu, DN Sheng, Ole xei I Motrunich, and Matthe w P A Fisher . Phase diagram of the spin-1/2 j 1 - j 2 heisenberg model on a square lattice. Physical Re view Letters , 113:027201, 2014. [9] Y in Hu, Shou-Shu Gong, W Zhu, and DN Sheng. Competing spin-liquid phases. Physical Revie w B , 92:140403, 2015. [10] Francesco Ferrari and Federico Becca. Dynamical structure factor of the j 1 - j 2 heisenberg model. Physical Revie w B , 98:100405, 2018. [11] Minoru Y amashita et al. Thermodynamic properties of a spin-1/2 spin-liquid state. Science , 328:1246–1248, 2010. [12] Christian Gross and Immanuel Bloch. Quantum simulations with ultracold atoms. Science , 357:995–1001, 2017. [13] Anders W Sandvik. Computational studies of quantum spin systems. AIP Confer ence Pr oceedings , 1297:135– 338, 2010. [14] Ste ven R White. Density matrix formulation for quantum renormalization groups. Physical Revie w Letters , 69:2863, 1992. [15] Ulrich Schollwöck. The density-matrix renormalization group. Annals of Physics , 326:96–192, 2011. [16] Juan Carrasquilla and Roger G Melko. Machine learning phases of matter . Natur e Physics , 13:431–434, 2017. [17] Evert PL v an Nieuwenburg, Y e-Hua Liu, and Sebastian D Huber . Learning phase transitions by confusion. Natur e Physics , 13:435–439, 2017. [18] K enny Ch’ng, Juan Carrasquilla, Roger G Melko, and Ehsan Khatami. Machine learning phases of strongly correlated fermions. Physical Review X , 7:031038, 2017. [19] Y i Zhang et al. Machine learning in electronic-quantum-matter imaging experiments. Nature , 570:484–490, 2019. [20] F Schindler , N Regnault, and T Neupert. Probing man y-body localization with neural networks. Physical Revie w B , 95:245134, 2017. [21] Lei W ang. Disco vering phase transitions with unsupervised learning. Physical Revie w B , 94:195105, 2016. [22] Sebastian J W etzel. Unsupervised learning of phase transitions. Physical Revie w E , 96:022140, 2017. [23] Joaquin F Rodriguez-Niev a and Mathias S Scheurer . Identifying topological order through unsupervised machine learning. Natur e Physics , 15:790–795, 2019. [24] Brandon Y ee, W ilson Collins, Caden W ang, and Mihir T ekal. Prometheus: Unsupervised discov ery of phase transitions and order parameters in the 2d ising model using variational autoencoders. In Pr oceedings of the AAAI Confer ence on Artificial Intelligence , 2026. Accepted. [25] Brandon Y ee, W ilson Collins, and Maximilian Rutkowski. From classical to quantum: Extending prometheus for unsupervised discov ery of phase transitions in three dimensions and quantum systems. 2026. [26] Phillip W einberg and Marin Bukov . Quspin: a python package for dynamics and exact diagonalisation. SciP ost Physics , 2:003, 2017. [27] Marin Bukov , Ale xandre GR Day , Phillip W einberg, Anatoli Polkovnik ov , P ankaj Mehta, and Dries Sels. Quspin: Part ii. SciP ost Physics , 7:020, 2019. 16 J1-J2 Heisenberg Disco very P R E P R IN T [28] Diederik P Kingma and Max W elling. Auto-encoding v ariational bayes. In ICLR , 2014. [29] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR , 2015. [30] Ling W ang, KSD Beach, and Anders W Sandvik. Plaquette valence-bond crystal in the frustrated heisenberg quantum antiferromagnet on the square lattice. Physical Review B , 73:014431, 2006. 17 J1-J2 Heisenberg Disco very P R E P R IN T A Q-V AE Ar chitecture Specifications The Q-V AE architecture is implemented with the following specifications: A.1 Network Ar chitecture Encoder: Three fully-connected hidden layers with dimensions [512, 256, 128], each follo wed by LayerNorm and ReLU activ ation. The final hidden layer projects to two output heads:  µ (latent mean) and log  σ 2 (log variance), each of dimension d z = 8 (configurable). Decoder: Symmetric architecture with hidden layers [128, 256, 512], each with LayerNorm and ReLU. The final layer outputs dimension 2 × dim( H S z =0 ) (real and imaginary parts), followed by e xplicit wa vefunction normalization: ψ norm = ψ out p P i ( | Re ( ψ i ) | 2 + | Im ( ψ i ) | 2 ) + ϵ , (38) where ϵ = 10 − 10 ensures numerical stability . A.2 T raining Hyperparameters • Learning rate: η = 10 − 3 with cosine annealing schedule • Batch size: 32 (or smaller if dataset size requires) • Maximum epochs: 1000 • Early stopping patience: 50 epochs without validation improv ement • KL di ver gence weight: β = 0 . 1 • Gradient clipping: max norm = 1.0 • T rain/validation split: 80%/20% • Optimizer: Adam with default momentum parameters A.3 Data A ugmentation W e exploit the S z symmetry of the Heisenber g Hamiltonian. For a w avefunction represented as [ Re ( ψ ) , Im ( ψ )] , the S z flip corresponds to complex conjugation ψ → ψ ∗ , mapping [ Re ( ψ ) , Im ( ψ )] → [ Re ( ψ ) , − Im ( ψ )] . This augmentation is applied with 50% probability during training. B Physical Obser vable Definitions All observables are computed exactly from wav efunction coef ficients using QuSpin operators. The spin correlation function ⟨  S i ·  S j ⟩ is decomposed as: ⟨  S i ·  S j ⟩ = ⟨ S z i S z j ⟩ + 1 2 ⟨ S + i S − j + S − i S + j ⟩ . (39) B.1 Staggered Magnetization In the S z = 0 sector , ⟨ S z i ⟩ = 0 ev erywhere. W e compute the staggered magnetization via spin-spin correlations: m 2 s = 1 N 2 X i,j ( − 1) i x + i y + j x + j y ⟨  S i ·  S j ⟩ , (40) with m s = p | m 2 s | . B.2 Stripe Order The stripe order parameter measures columnar magnetization in both lattice directions: m stripe = max  p |⟨ M x ⟩| 2 , q |⟨ M y ⟩| 2  , (41) where M x = N − 1 P i ( − 1) i x  S i and M y = N − 1 P i ( − 1) i y  S i . 18 J1-J2 Heisenberg Disco very P R E P R IN T B.3 Plaquette Order The plaquette order parameter measures four-spin correlations on elementary plaquettes: P = 1 N p X p ⟨  S 1 ·  S 2 ⟩ p ⟨  S 3 ·  S 4 ⟩ p , (42) where the sum runs over all N p = L 2 plaquettes with periodic boundary conditions, and sites { 1 , 2 , 3 , 4 } are the four corners of each plaquette. B.4 Nematic Order The nematic order parameter measures the anisotropy between x and y bond correlations: Q =      1 N X i ⟨  S i ·  S i + ˆ x ⟩ − 1 N X i ⟨  S i ·  S i + ˆ y ⟩      . (43) B.5 Dimer Order The dimer order parameter measures alternating bond strength: D = 1 2 N      X i ( − 1) i x + i y  ⟨  S i ·  S i + ˆ x ⟩ + ⟨  S i ·  S i + ˆ y ⟩       . (44) B.6 Structure F actor The spin structure factor at wa ve vector  q is: S (  q ) = 1 N X i,j e i q · (  r i −  r j ) ⟨  S i ·  S j ⟩ . (45) W e compute S (  q ) at key wa ve vectors: ( π , π ) (Néel ordering), ( π, 0) and (0 , π ) (stripe ordering). B.7 Entanglement Entropy For a bipartition into subsystems A and B , the von Neumann entropy is computed via singular v alue decomposition of the wa vefunction reshaped as a matrix ψ i A ,i B : S A = − X k λ 2 k log λ 2 k , (46) where λ k are the singular values. C Computational Details and Reproducibility C.1 Software Stack • Python 3.8+ • QuSpin 0.3.6+ for exact diagonalization • PyT orch 1.9+ for Q-V AE implementation • NumPy , SciPy for numerical operations • scikit-learn for clustering and dimensionality reduction • Pydantic for configuration v alidation • HDF5 (h5py) for data storage 19 J1-J2 Heisenberg Disco very P R E P R IN T C.2 Exact Diagonalization Details The Hamiltonian is constructed in the S z tot = 0 sector using QuSpin’ s spin_basis_general with N up = N / 2 . Ground states are computed using SciPy’ s eigsh (ARP A CK Lanczos) with con vergence tolerance 10 − 10 and maxi- mum 1000 iterations. The Hamiltonian is stored in CSR sparse format. For L = 4 : dim( H S z =0 ) =  16 8  = 12 , 870 C.3 Checkpointing and Resumption The pipeline supports checkpointing at multiple stages: • ED checkpoints: Ground states saved per lattice size after completion • Q-V AE checkpoints: Model weights, training history , and best validation state • Analysis checkpoints: Latent representations, observables, and intermediate results All data is stored in HDF5 format with metadata for reproducibility . C.4 Parallel Execution ED computations support parallel ex ecution across parameter points using Python’ s multiprocessing.Pool . Mem- ory monitoring and garbage collection are performed between computations to manage resources for large Hilbert spaces. C.5 V alidation Framework The validation module performs automated checks: • Néel phase v alidation: Requires | r ( latent , m s ) | ≥ 0 . 7 for J 2 /J 1 ∈ [0 , 0 . 4] • Stripe phase v alidation: Requires | r ( latent , m stripe ) | ≥ 0 . 7 for J 2 /J 1 ∈ [0 . 6 , 1 . 0] • Literature comparison: Critical points compared against DMRG/QMC bounds All code and data will be released under open licenses at https://github.com/YCRG- Labs/ prometheus- heisenberg- discovery . D Complete Observable Data T able 4 presents the complete set of computed observables for all 41 parameter points. All quantities are computed exactly from e xact diagonalization ground states without Monte Carlo sampling. E Latent-Observable Corr elation Matrix T able 5 presents the complete Pearson correlation matrix between all eight latent dimensions and the ph ysical observ- ables with non-zero variance. Correlations with | r | ≥ 0 . 7 are highlighted. The correlation structure re veals that the Q-V AE has organized its latent space such that z 0 and z 2 capture the dominant variation associated with the magnetic ordering transition. The remaining dimensions ( z 1 , z 3 – z 7 ) show moderate to weak correlations, suggesting the y encode secondary features of the ground state wa vefunctions not directly related to the measured order parameters. 20 J1-J2 Heisenberg Disco very P R E P R IN T T able 4: Complete observable data for L = 4 across all parameter points. Energy density e , staggered magnetization m s , plaquette order P , and structure factors S ( π, π ) , S ( π, 0) are shown. Stripe order , nematic order, and dimer order are omitted as they are numerically zero ( < 10 − 11 ) for all points. J 2 /J 1 e m s P S ( π, π ) S ( π, 0) 0.30 − 2 . 332 0.902 1.890 13.02 − 0 . 185 0.31 − 2 . 318 0.897 1.882 12.87 − 0 . 179 0.32 − 2 . 304 0.892 1.873 12.72 − 0 . 172 0.33 − 2 . 291 0.886 1.863 12.56 − 0 . 164 0.34 − 2 . 278 0.880 1.852 12.38 − 0 . 155 0.35 − 2 . 265 0.873 1.840 12.20 − 0 . 145 0.36 − 2 . 252 0.866 1.827 12.01 − 0 . 133 0.37 − 2 . 240 0.859 1.813 11.81 − 0 . 120 0.38 − 2 . 228 0.851 1.798 11.59 − 0 . 105 0.39 − 2 . 216 0.843 1.781 11.37 − 0 . 088 0.40 − 2 . 205 0.834 1.763 11.13 − 0 . 068 0.41 − 2 . 193 0.825 1.743 10.89 − 0 . 045 0.42 − 2 . 183 0.815 1.722 10.62 − 0 . 019 0.43 − 2 . 172 0.804 1.699 10.35 +0 . 011 0.44 − 2 . 163 0.793 1.673 10.06 +0 . 046 0.45 − 2 . 153 0.781 1.645 9.75 +0 . 086 0.46 − 2 . 144 0.768 1.615 9.43 +0 . 134 0.47 − 2 . 136 0.754 1.582 9.10 +0 . 189 0.48 − 2 . 128 0.739 1.547 8.75 +0 . 254 0.49 − 2 . 121 0.724 1.508 8.38 +0 . 331 0.50 − 2 . 114 0.707 1.466 8.00 +0 . 421 0.51 − 2 . 109 0.689 1.420 7.60 +0 . 527 0.52 − 2 . 104 0.670 1.369 7.19 +0 . 652 0.53 − 2 . 100 0.650 1.314 6.77 +0 . 800 0.54 − 2 . 096 0.629 1.254 6.32 +0 . 976 0.55 − 2 . 094 0.606 1.188 5.87 +1 . 185 0.56 − 2 . 093 0.581 1.117 5.40 +1 . 431 0.57 − 2 . 094 0.554 1.039 4.91 +1 . 721 0.58 − 2 . 095 0.525 0.954 4.42 +2 . 061 0.59 − 2 . 099 0.494 0.864 3.90 +2 . 454 0.60 − 2 . 104 0.460 0.768 3.38 +2 . 901 0.61 − 2 . 110 0.422 0.670 2.85 +3 . 399 0.62 − 2 . 119 0.381 0.573 2.32 +3 . 937 0.63 − 2 . 130 0.335 0.480 1.80 +4 . 497 0.64 − 2 . 143 0.285 0.395 1.30 +5 . 058 0.65 − 2 . 158 0.229 0.321 0.84 +5 . 596 0.66 − 2 . 174 0.163 0.259 0.42 +6 . 093 0.67 − 2 . 192 0.061 0.209 0.06 +6 . 538 0.68 − 2 . 212 0.126 0.169 − 0 . 25 +6 . 926 0.69 − 2 . 233 0.180 0.138 − 0 . 52 +7 . 258 0.70 − 2 . 255 0.216 0.114 − 0 . 74 +7 . 541 T able 5: Complete correlation matrix between latent dimensions and physical observables. V alues with | r | ≥ 0 . 7 indicate strong correlations. Observ ables with zero v ariance (stripe order, nematic order , dimer order , entanglement entropy) are omitted. Observable z 0 z 1 z 2 z 3 z 4 z 5 z 6 z 7 Energy 0.98 0.55 0.90 − 0 . 24 − 0 . 60 − 0 . 52 0.25 0.54 m s − 0 . 97 − 0 . 54 − 0 . 90 0.22 0.59 0.49 − 0 . 22 − 0 . 53 S ( π , π ) − 0 . 97 − 0 . 54 − 0 . 90 0.22 0.59 0.49 − 0 . 23 − 0 . 53 S ( π , 0) 0.96 0.53 0.89 − 0 . 20 − 0 . 58 − 0 . 46 0.20 0.51 S (0 , π ) 0.96 0.53 0.89 − 0 . 20 − 0 . 58 − 0 . 46 0.20 0.51 Plaquette − 0 . 97 − 0 . 54 − 0 . 90 0.22 0.59 0.48 − 0 . 22 − 0 . 52 21