In this work, we study the magnetic phases of a spatially modulated chain of spin-1 Rydberg excitons. Using the Density Matrix Renormalization Group (DMRG) technique, we study various magnetic and topologically nontrivial phases using both single-particle properties, such as local magnetization and quantum entropy, and many-body ones, such as pair-wise Néel and long-range string correlations. In particular, we investigate the emergence and robustness of the Haldane phase, a topological phase of anti-ferromagnetic spin-1 chains. Furthermore, we devise a hybrid quantum algorithm employing restricted Boltzmann machine to simulate the ground state of such a system that shows very good agreement with the results of exact diagonalization and DMRG.

## I. INTRODUCTION

Quantum spin systems have been the subject of intense studies for a long time as they offer a playground for exploring fascinating physics.^{1} In 1983, Haldane showed that chains of half-integer and integer spins are profoundly different.^{2,3} In particular, while the Heisenberg antiferromagnetic chains have a ground state with an excitation gap and exponentially decaying correlations for integer spins, the ground state of half-integer spin chains has no excitation gap and the correlation decays algebraically. Since then, the rich physics of spin chains have been investigated extensively, leading to the discovery of intriguing features, such as half-integer spin edge modes^{4,5} and non-local string orders.^{6,7}

Within the last two decades, significant advances in the control of various quantum systems have led to the emergence of the quantum simulator plethora in well-controlled laboratory setups based on ultracold atoms in optical lattices, superconducting circuits, and quantum dots, to name a few.^{8,9} Such controllable quantum systems have been successfully used for the study of spin systems properties, such as magnetism,^{10–13} transport,^{14} and topology.^{15–18} Recently, there have been several attempts to simulate various condensed-matter models using an array of Rydberg atoms, i.e., highly excited atoms with macroscopic sizes and very susceptible to the external fields.^{19} The enhanced polarizability of Rydberg atoms results in large, long-range interactions and makes them a suitable platform for the investigation of various spin systems.^{20–23}

On the other hand, excitons, i.e., the quasi-particles of bounded electron–hole pairs in semiconductors, offer an alternative platform for the study of many-body quantum systems in a low-dimensional integrable and scalable platform at high densities not accessible with ultracold atoms. Like atoms, unique favorable scaling, such as strong long-range interaction, and Rydberg blockade effect are expected for Rydberg excitons, as well. The first observation of Rydberg excitons up to *n* = 25 in cuprous oxide (Cu_{2}O) in 2014^{24} revived the idea of excitonic quantum simulators and the field of semiconductor Rydberg physics.^{25}

Here, we study the emergent magnetic phases of a generalized spin-1 chain. The Hamiltonian parameters and the coupling coefficients are based on values attainable with Rydberg excitons in Cu_{2}O, where the spin-1 degrees of freedom can be mimicked with optically active *p*-excitons. Furthermore, we develop a hybrid quantum–classical algorithm to map this generalized spin-1 Hamiltonian to an spin-$12$ Ising Hamiltonian using Restricted Boltzmann Machine (RBM) architecture and simulate the ground state of some prototypical examples. It must be noted that recently such hybrid quantum simulations to explore excitonic properties of spin-$12$ systems have begun to gain attention, such as in Ref. 26, which creates an entangled state of photon–hole pairs representing an exciton condensate on a 53-qubit quantum computer using a Lipkin Hamiltonian.

The paper is organized as follows: In Sec. II, we introduce the model describing the most general dynamics of a spin-1 chain with nearest–neighbor interactions in a spatially modulated potential. In Sec. III, we employ the density matrix renormalization group (DMRG) to study various magnetic phases emerging in this system with physical parameters based on Rydberg excitons in Cu_{2}O. We classify the aforementioned many-body phases using pair-wise and long-range string correlations, highlighting the emergence of Néel order and topological Haldane phase at different modulated potentials. Section IV puts forward a hybrid algorithm for the efficient simulation of studying stationary states on noisy-intermediate-scale quantum (NISQ) hardware. Since we are studying spin-1 particles, the Hilbert space grows as 3^{N}, making the problem intractable rather quickly. We propose an algorithm that makes this scaling linear in the system size *N*, which allows us to simulate the exact dynamics of the model on a quantum computer. We summarize the results and conclude the work in Sec. V and discuss the follow-up directions.

The results of this work pave the way for the study of spin-1 chains, supporting diverse topological and magnetic phases, on an analog excitonic quantum simulator. Furthermore, it puts forward an efficient algorithm to study the stationary states of such systems on a gate-based digital quantum computer for the first time.

## II. THE MODEL

In this section, we introduce the model of spin-1 chains characterized by the number of spins, *N*. The most general Hamiltonian with nearest–neighbor interaction is as follows:^{27}

where $Si\alpha $ are the components of the usual *S* = 1 spin operators with *α* ∈ {*x*, *y*, *z*}. For the periodic boundary conditions (PBC), *j*_{max} = (*N* − 1) with *j*_{max} + 1 = 0, whereas for the open boundary conditions (OBC), *j*_{max} = (*N* − 2). Several well-known spin models such as the transverse field Ising model,^{28,29} commonly studied Heisenberg spin models such as XXZ^{30,31} or XYZ,^{32,33} and the Affleck–Kennedy–Lieb–Tasaki (AKLT) model with bi-quadratic interaction terms^{34–38} are special cases of Eq. (1) with specific choices of coefficients ${ci}i=06$. Models beyond nearest neighbor interaction can also display a rich phase diagram as has been investigated recently^{39} in the context of fermion-pair condensation and exciton-pair condensation.

The realization of such a spin model on an array of Rydberg excitons in cuprous oxide has been proposed recently^{40,41} with the promise of attaining topologically non-trivial phases. Such excitons are created through optical excitations into Rydberg states of varying principal quantum numbers *n* with the azimuthal quantum number *l* = 1 (*p*-shell). The latter endows the array to behave like an effective pseudospin-1 particle with the van der Waals interaction between excitons leading to the terms in Eq. (1) under the assumption that only one exciton occupies each trapping site and the direct exchange interaction is negligible. The coefficients *c*_{0}–*c*_{6} for the said platform are listed in Table I in terms of a suitable energy scale *ϵ*.^{40}

. | Coefficients . | Value/ϵ
. | . | |
---|---|---|---|---|

c_{0} | −5.58 | |||

c_{1} | 9.53 | |||

c_{2} | −8.97 | |||

c_{3} | 1.27 | |||

c_{4} | 6.59 | |||

c_{5} | −3.18 | |||

c_{6} | 5.04 |

. | Coefficients . | Value/ϵ
. | . | |
---|---|---|---|---|

c_{0} | −5.58 | |||

c_{1} | 9.53 | |||

c_{2} | −8.97 | |||

c_{3} | 1.27 | |||

c_{4} | 6.59 | |||

c_{5} | −3.18 | |||

c_{6} | 5.04 |

Figure 1(a) shows that the energy scale *ϵ* escalates steeply with increasing principal quantum number *n* (∝*n*^{11}) for the Rydberg excitons at a fixed *R*_{0}, i.e., the inter-exciton distance. This quantity decays with increasing *R*_{0} ($\u221dR0\u22126$) due to the reduction in the van der Waals interaction between the excitons.

For Rydberg excitons in Cu_{2}O considered in this work, due to the same symmetry of the valence and conduction bands, the *p*-orbitals are optically active; hence, the three orbital angular momenta, *l* = 0, ±1, behave as a pseudo spin-1. A simple anisotropic potential behaving as $Sz2$ can be used to partially remove the degeneracy (between *m*_{l} = 0 and *m*_{l} = ±1) of these three orbitals. For a system of *N* excitons, if one assumes a parabolic harmonic trap for each site, lifting of degeneracy of the triplet *p*-excitons allows one to study richer and interesting physical systems with effects inaccessible without the use of such local perturbation. Indeed, as we shall see, a rich phase diagram emerges wherein the experimenter can switch from one-phase to another by just tuning this on-site anisotropy term. As a result, herein, we further investigate the same spin model with locally modulated on-site potential defined as **D** = {*D*_{j} = *D*_{1} + (−1)^{j}*D*_{2}|∀*j* ∈ **Z**_{N}, (*D*_{1}, *D*_{2}) ∈ **R**^{2}}.

The resultant spin system is similar to a Creutz ladder^{42–44} with two mutually interacting spin graphs with different local anisotropies characterized by (*D*_{1}, *D*_{2}) (in units of the scale *ϵ*) in an overall Hamiltonian of the form

The local anisotropy profile for a particular realization of (*D*_{1}, *D*_{2}) is schematically represented in Fig. 1(b). This periodic modulation can be thought as two interlaced 1D exciton traps with different depths. Such 1D arrays can be easily realized experimentally using different approaches, e.g., via a periodic arrangement of static stressors^{45,46} or surface acoustic waves^{47} to locally deform the crystal lattice and trap the excitons, or using spatial light modulators^{48} to create excitons at desired locations,^{49} or through the optical confinement.^{50} If a spatial light modulator (SLM) is used, the connectivity graph of the Hamiltonian can be changed *in situ,* which is an advantage over the strain traps that are static and hence not programmable. To realize the strong-limit of the Creutz ladder, one can envision the even (odd) sites in Fig. 1(b) as the nodes of the upper (lower) rungs.^{42} The flexibility of changing the interaction strength between adjacent sites through changing the excitonic Rydberg states, as depicted in Fig. 1(a), allows one to tune the interaction strength between the nodes *in situ* and explore various phases of the ladder model. Similar anisotropic potentials with modulation after two sites can be used to simulate pair-wise interacting models, such as the Su-Schreiffer Heeger (SSH) model,^{51} as well. Properties of such models can be directly revealed using single-site as well as non-local optical correlation measurements native to the platform.

Since the *c*_{1} > 0, the ground state of Eq. (1) has an anti-ferromagnetic character,^{40} the predominance of which can be conveniently tuned using a local anisotropy profile as in Eq. (2). The modulated potential considered in this case serves to probe deeper into the physics of the model by demarcating the regime wherein such spin phases can be stabilized at will, including the topologically non-trivial Haldane phase, thereby highlighting its robustness. Furthermore, it provides a pathway to transition from trivial phases to topological ones by controlling experimentally accessible parameters by traversing a richer phase diagram hitherto unexplored.

Various phases would be characterized by one-body properties, such as local magnetization and von Neumann entropy, followed by many-body properties, such as Néel correlation^{52} and long-range string correlations.^{34} In Sec. III, we first explicate the classical resources used for the simulation and then present a linear-scaling quantum algorithm designed for the study of such systems.

## III. CLASSICAL SIMULATIONS

The classical simulations of the spin model have been performed using the DMRG algorithm introduced by White^{53} and later applied to many physical^{54–57} and chemically relevant problems.^{58–60} The algorithm is extremely popular to circumvent the exponential scaling of many-body Hilbert spaces,^{61} especially in 1D, which is also the scope of this work. A detailed overview of such algorithms can be found elsewhere.^{61,62}

We used the Julia version of ITensor as well as the DMRGPY wrapper^{63} from the ITensor library for computation with maximum bond dimensions of the required Matrix Product State (MPS) set to be 250 for each sweep. 100 sweeps are used over the linear array with a cutoff energy of convergence at 10^{−9} per sweep. For excited state computations (not discussed in this paper), a weight parameter *ω* to penalize the overlap with the ground state may be used with steadily decreasing noise profile $\u226410\u22126$ for a crisp convergence.

Figures 2(a) and 2(b) shows the ground state energy profile in the (*D*_{1}, *D*_{2}) plane, subject to PBC for two chain sizes *N* = 4 and *N* = 16, respectively. We see that the ground state energy is lower for smaller values of *D*_{1} irrespective of the modulation depth *D*_{2} (in the range studied) and rises when *D*_{1} is escalated. We analyze the properties of the spin chain at the designated points (A–I) marked in Fig. 2(b) for *N* = 16. (Similar trends for different chain lengths can be observed.) Using spin correlations, in Sec. III A, we discuss the implications of this result and show that for (|*D*_{1}|, |*D*_{2}|) ≈ (0, 0), e.g., point E, there is a distinct Haldane phase, whereas for (*D*_{1}, *D*_{2}) ≪ 0, e.g., point A, we have an anti-ferromagnet, and for (*D*_{1}, *D*_{2}) ≫ 0, e.g., point I, we have the large |*D*|-limit disordered phase.^{64–67}

### A. Single and many-body properties

To determine different phases obtainable in this periodically modulated spin-1 chain, we study certain order parameters such as local magnetization $(\u27e8Sz(i)\u27e9)$, its associated variance $\sigma Sz,i$, and the von Neumann entropy of the reduced single-spin density matrix $\rho i1$ along the chain as well as the pair-wise Néel [*C*_{N}(*i*)] and the long-range string [*C*_{S}(*i*)] correlation defined as

and

The local average of the magnetization $\u27e8Szi\u27e9$ is plotted in Fig. 3(a), and the von Neumann entropy is shown in Fig. 3(b), respectively, for points (*A*, *B*, *D*, *E*, *G*, *H*) defined in Fig. 2(b). Since a similar local anisotropy profile is established in (*C*, *F*, *I*) as in (*A*, *D*, *G*), respectively, with the odd and the even sites swapped [see Eq. (2)], the former points are omitted from further discussion.

We see that for points (*D*, *E*, *G*, *H*), the magnetization $\u27e8Szi\u27e9$ at each site, depicted with red stars, is predominantly zero with a non-vanishing variance of $Siz$, locally. This essentially is a consequence of attaining a many-body quantum state as a superposition of many contributing configurations. However, that is not the case for (*A*, *B*) due to a preference of one of the symmetry-broken configurational states as is indicated by the low on-site variance of *S*^{z}.

This is further bolstered in the corresponding von Neumann entropy profile presented in Fig. 3(b) calculated from the reduced density matrix at each site $(\rho i1)$. As can be seen, the entropy attains a low value for (*A*, *B*) due to the preference of one of the uncorrelated configurations, leading to a pure one-site reduced density matrix $(\rho i1)$ at all sites *i*. The single-site entropy for other points (*D* and *E*, *G* and *H*), however, are near maximal, indicating an almost maximally mixed reduced density matrix $(\rho i1)$ at each site arising from a many-body state as the superposition of computational basis states.

To determine the phase of the spin chain at different points, we investigate the correlations as shown in Fig. 4, where the blue stars show the Néel order and the red circles indicate the string correlation. For points (*A*, *B*) wherein *D*_{1} = −8, the Néel correlation is sustained and non-decaying irrespective of *D*_{2}. The modulating Hamiltonian, i.e., the second term in Eq. (2), leads to a stabilized ground state where sites with $Spz=\xb11$ are favored over $Spz=0$.

In the disordered phase at *D*_{1} = 8, however, e.g., (*G*, *H*), the local anisotropy at each site is positive, and hence, unlike the previous case, configurations with $Spz=0$ are favored due to the non-negative penalty of $Spz=\xb11$ sites enforced by the local modulation. As a result, many different configurations contribute significantly to the ground state, and the overall anti-ferromagnetic order is lost with the lack of exclusive dominance from Néel-ordered configurations.

When *D*_{1} = 0, e.g., at (*D*, *E*), the local anisotropy is ±8 for alternating spins, and hence, partial remnants of anti-ferromagnetic order as for *D*_{1} = −8 damped quickly due to the sites with local $Spz=0$ being favored, owing to the non-negative penalty as in (*G*, *H*). We also see a non-decaying string correlation in points (*D*, *E*) of Fig. 4. It must be emphasized that due to the nature of the exponential term in Eq. (4), which can contribute −1 for sites with $Szp=\xb11$ and 1 for sites with $Szp=0$, the string correlation remains negative for all points (*A*–*H*), except at *j* = 0. The said order between site indexed at “0” and at “*j*” can only acquire non-zeros values for configurations wherein the terminal spin $(Sz0,Szj)$ are in any one of the four states (±1, ±1) since the *S*_{z} = 0 on either spin would destroy the correlation order defined in Eq. (4). The non-negativity of string order can be qualitatively understood for certain phases using the number of contributing configurations at least for the case of spin indices at the terminus of the excitonic array i.e., *j* = *N* − 1. For all such configurations in which the terminal spins are of opposite spin-type, the intermediary spins that contribute to the exponential term in the string order must have $\u2211p=1j\u22121Szp=0$ to ensure an overall $\u2211jSzj=0$, which is the symmetry of the ground as well as the first excited state (cf. Sec. S1 in the supplementary material). Therefore, the total number of such configuration is proportional to $2(\u2211i=1(j\u22121)/2j\u22121ij\u22121\u2212ii+1)$, all having the string correlation of −1. On the other hand, when the terminal spins are in the same spin state (non-zero) for which the string order is +1, the intermediary spins in such configurations must have $\u2211p=1j\u22121Szp=\xb12$, which is proportional to a total count of $2(j\u221212\u2211i=1(j\u22123)/2j\u22123ij\u22123\u2212ii+1)$. For (*A*, *B*), the dominant configuration(s) are of the first kind with terminal spins at (±1, ∓1), which explains their larger negative string correlations compared to other points (cf. Sec. S2 in the supplementary material).

The above discussion together with the decaying Néel order and the sustenance of string correlation at points (*D*, *E*, *F*) (especially for *E*) supports the emergence of the topological Haldane phase for this trap configurations. At point *E*, there is no perturbation; hence, we recover the properties of the Haldane phase discussed in Ref. 40. In Sec. IV, we discuss how such spin-1 system can be mapped to an RBM based off spin-1/2 models with a hybrid-algorithm linearly scaling in qubit numbers with the chain length *N*, thereby allowing it to be studied on NISQ era quantum computers.

## IV. QUANTUM SIMULATION USING RESTRICTED BOLTZMANN MACHINE (RBM)

Quantum machine learning algorithms^{68} based on parameterized quantum circuits have begun to gain enormous attention for understanding the physical properties of atomic and molecular systems.^{69–72} One such network, commonly known as the Restricted Boltzmann Machine,^{73} has been employed successfully as a neural-network ansatz for a quantum state^{74} in many different applications, such as in fermionic assemblies,^{75–77} in anyonic symmetries,^{78} in supervised phase-classification tasks,^{79,80} in dynamical evolution,^{81} and in chemistry^{82–85} or even for classification of classical data.^{86,87} The network is a universal approximator for any probability density^{73,88} and can efficiently simulate a volume-law entangled state even when sparsely parameterized.^{89} It has been provably demonstrated that evaluating the full distribution classically would require exponential resources as long as polynomial hierarchy is not collapsed;^{90} however, recently some of the authors have reported a quantum algorithm training the said ansatz to efficiently simulate the ground and excited state properties of materials and molecules with quadratic resources.^{82}

### A. The generative model/ansatz for the task

The specific description of the network can be found in Fig. 5(a). The core idea involves encoding the target state of the desired Hamiltonian (to be called the driver system) into a neural-network architecture that resembles a bipartite spin graph *G* with tunable connectivity. Formally, the graph can be envisioned as *G* = (*V*, *E*). The set of vertices *V* (henceforth to be called neurons) is further divided as $V={\nu}i=1m+n={\sigma}i=1n\u22c3{h}j=1m\u22c3{p}k=12$ with $(m,n)\u2208Z++$, where $Z++$ means all strictly positive integers, i.e., $Z++={x>0|\u2200x\u2208Z}$. The neurons of the type $\sigma \u20d7\u2208{1,\u22121}n$ [shown in blue in Fig. 5(a)] collectively form the visible node register. For a given Hamiltonian matrix $H\u2208Cd\xd7d$ of the driver system, the number *n* is chosen such that *n* = ⌈ log_{2}(*d*)⌉. The neurons of the type $h\u20d7\u2208{1,\u22121}m$ [shown in green in Fig. 5(a)] constitute the hidden node register. This number *m* can be chosen arbitrarily to accomplish the desired accuracy threshold (Here, *m* ∼ *n*). The role of the hidden node register is to escalate the expressibility of the generative model by increasing the number of controllable parameters and inducing higher-order hidden correlation among the spins of the visible register.^{73,91,92} Two neurons of the type *p*_{k} ∈ [−1, 1] shown in orange and red in Fig. 5(a) constitute the phase nodes whose role will be discussed shortly.^{93}

It is also assumed that all possible pairs of neurons *σ*_{i} and *h*_{j} are interconnected. Furthermore, each *σ*_{i} shares edges with the two neurons ${p}k=12$, also. Collectively, all such edges are labeled as ${eij}i=1,j=1n,m+2$ and constitute the set *E* where |*E*| = *mn* + 2*n*. Associated with the vertices *ν*_{i} ∈ *V*, we define a bias vector $(a\u20d7,b\u20d7,c\u20d7)\u2208Rm+n+2$, where $a\u20d7={a}i=1n$ is for neurons in the visible-node register $\sigma \u20d7$, $b\u20d7={b}j=1m$ for neurons in the hidden-node register $h\u20d7$, and (*c*_{1}, *c*_{2}) for the two neurons ${p}k=12$ in the phase-node register. Similarly, associated with *e*_{ij} ∈ *E*, where *i* ∈ ⌈*n*⌉ and *j* ∈ ⌈*m*⌉, we define a weight matrix $W\u20d7\u2208Rn\xd7m$. Together with the tunable parameters $(a\u20d7,b\u20d7,W\u20d7)$, we define the energy function $E(a\u20d7,b\u20d7,W\u20d7,\sigma \u20d7,h\u20d7)$ as

This subset of the network encodes a probability distribution, which is essentially the classical thermal state of the corresponding Ising spin Hamiltonian in Eq. (5) defined as follows:^{96–99}

This subset of *G* thus constitutes the usual RBM network.

Any realization of spin configuration $(\sigma \u20d7,h\u20d7)$ of the combined registers of (*m* + *n*) spins is sampled from the distribution in Eq. (6). Since the parameter vector $(a\u20d7,b\u20d7,W\u20d7)$ is real-valued, the functional form in Eq. (6) mimics a probability distribution, only. To account for the phase of the target wavefunction, herein, we encode the following complex-valued function within the neurons ${pk}k=12$ through the following activation function:^{93}

where $d\u20d7,f\u20d7\u2208Rn$ are the weights for the connections between ${p}k=12$ and ${\sigma}i=1n$ [see Fig. 5(a)]. The full ansatz thus consists of neurons of the RBM as well as the two neurons of the phase node with a complete set of parameter vector $X=(a\u20d7,b\u20d7,W\u20d7,d\u20d7,f\u20d7,c1,c2)$ that can be adjusted variationally during the training process.

### B. The outline of the algorithm and results

The purpose of using the aforesaid generative network is to construct an ansatz for the amplitude and the phase of the target state wavefunction for the Hamiltonian defined in Eq. (2) as follows:

wherein |*σ*_{1}, *σ*_{2}, …, *σ*_{n}⟩ represents the spin configurations of the visible node register ${\sigma}i=1n$, which forms the computational basis in which the target wavefunction is resolved. The parameter vector $X\u20d7=(a\u20d7,b\u20d7,c\u20d7,W\u20d7,d\u20d7,f\u20d7)$ is tuned until the energy of the ansatz state is variationally minimized. The wavefunction in Eq. (8) obtained at the end of this training process has a large overlap with the target state.

Recently, the authors have developed an algorithm to retrieve the distribution directly using a quantum circuit^{82} with a circuit width, a number of gates, and parameter count of *O*(*m* × *n*) for a visible (hidden) node register of *n*(*m*) qubits. Since for *S* = 1, conversion of the problem to qubit degrees of freedom already leads to wasteful circuit widths, using the aforesaid quadratic scaling algorithm for this system would be inefficient. To this end here, we develop a *linear* scaling algorithm *O*(*m* + *n*), which also retrieves the required distribution function based on an arbitrary Gibbs-state preparation scheme as described in the following. The overarching theme of our algorithm involves a double optimization protocol. The outer-optimization involves tuning the parameter vector $X\u20d7$, whereas the inner optimization involves constructing the Gibbs state, given an $(a\u20d7,b\u20d7,W\u20d7)$ defining a connected Ising spin system. In other words, since our protocol involves neural-network encoding of the target state [see Fig. 5(a)], the inner optimization loop is entrusted with preparation of the thermal state representing the ansatz in its desired structural form for a fixed incumbent instance of the parameter vector $X\u20d7$, whereas the outer loop trains the ansatz and changes the network parameters $X\u20d7$ following the usual gradient based updates. Preparation of such thermal states in the inner loop is, in general, a non-trivial task as it involves non-unitary transformation from a pure state with entropy production. We shall see that this is aided by mediating the interaction of the register of (*m* + *n*) qubits with an ancillary qubit, leading to the formation of a mixed state (thermal state) in the reduced space of the data register of (*m* + *n*) qubits. In the absence of inner optimization loop, convergence would be poor as the subsequent analytical gradients are reliant on the form of the structural ansatz, whereas in the absence of the outer-optimization loop, no training would happen at all.

The algorithm starts by randomly assigning values to the parameter vector $X\u20d7=(a\u20d7,b\u20d7,W\u20d7,c\u20d7,d\u20d7,f\u20d7)$ by sampling from a uniform distribution within [−0.02, −0.02] to avoid the problem of vanishing gradient.

^{84}If training with random initialization is not successful within the desired accuracy threshold (say*η*), an initial parameter set of a previously converged point in a similar optimization problem is used as the initial guess. This process is called warm starting^{82,83}and is routinely used in machine learning for non-convex functions as is the case here.The next step is using the $(a\u20d7,b\u20d7,W\u20d7)$ set to construct a quantum circuit that would return the thermal state $\rho th(a\u20d7,b\u20d7,W\u20d7)$ of (

*m*+*n*) qubits corresponding to the Ising energy function [defined in Eq. (5)]. The probability distribution that encodes the amplitude field of the ansatz as defined in Eq. (6) is then retrieved from the diagonal elements of this thermal state. To construct the thermal state for the specific RBM energy function we follow the protocol designed in Ref. 100 using an ansatz $U(\theta \u20d7)$ [see Fig. 5(b)] with a circuit width of*O*(*m*+*n*) and a depth of*D*. The depth*D*has to be tuned by the user for attaining the accuracy threshold of interest. For this work,*D*= 3 suffices for all examples. The number of gates required in the circuit is*O*(*D*(*m*+*n*)). The protocol relies on minimizing the Helmholtz free-energy of (*m*+*n*) qubits as the cost-function $F(\rho (\theta \u20d7,a\u20d7,b\u20d7,W\u20d7))$ with respect to $\theta \u20d7$, keeping $(a\u20d7,b\u20d7,W\u20d7)$ fixed at the input value. The extra ancillary qubit acts as a digital bath to facilitate the formation of a mixed state for the register with*m*+*n*qubits. (cf. Sec. S4 in the supplementary material for more information about this step).- With $P(a\u20d7,b\u20d7,W\u20d7,\sigma \u20d7,h\u20d7)$ from the previous step of the algorithm, the phase function defined in Eq. (7) was computed using the remaining part of the parameter set $(c\u20d7,d\u20d7,f\u20d7)$. Combining the two, the wavefunction $\psi (X\u20d7)$ can now be constructed and the cost-function for the outer loop of the optimization be evaluated. Since the objective is learning the eigenstates of the Hamiltonian $H\u0303$ in Eq. (2), we use the following cost function for tuning the parameter vector $X\u20d7$:where $\lambda \u2208R++$ is the penalty set as a hyper-parameter. The symbol $R++$ means all strictly positive real numbers, i.e., $R++={x>0|\u2200x\u2208R}$. $O\u0302$ is a user-defined symmetry operator $([H\u0303,O\u0302]=0\u0302)$ of the system whose eigenstate is desired with a specific eigenvalue(9)$C(|\psi (X\u20d7)\u3009,H\u0303,O\u0302,\lambda )=\u27e8\psi |H\u0303|\psi \u27e9+\lambda \u27e8\psi |(O\u0302\u2212\omega )2|\psi \u27e9,$
*ω*. The same prescription allows us to not only target specific states based on symmetry but also to explore excited states by sampling the orthogonal complement of the ground state projection operator*P*_{gr}, i.e., by choosing $O\u0302=Pgr$ and*ω*= 0. Quantum algorithms much like this one have now been formulated, which can simulate excited states efficiently.^{101,102}One must note that for the ground state*λ*= 0 is substituted in Eq. (9) as the usual minimization of energy suffices. Once the cost-function $C(|\psi (X\u20d7)\u3009,H\u0303,O\u0302,\lambda )$ is evaluated in the outer loop, one can check for convergence (if $C(|\psi (X\u20d7)\u3009,H\u0303,O\u0302,\lambda )\u2264\eta $ where

*η*is the convergence threshold) or if the maximum number of iterations has been completed. In either case, the resultant state $|\psi (X*\u20d7)\u3009$ and energy $\u27e8\psi (X*\u20d7)|H\u0303|\psi (X*\u20d7)\u27e9$ are printed wherein $X*\u20d7=argminX\u20d7C(|\psi (X\u20d7)\u3009,H\u0303,O\u0302,\lambda )$.

However, if either condition is not satisfied, then the parameter vector $X\u20d7$ is updated using a gradient-based algorithm (ADAM optimizer^{103} is used in this case) with learning rate $\alpha 2\u2208R++$ as

With the updated parameter set, one returns to the second step of the algorithm for the next iteration. The process is continued until the desired threshold convergence is reached or the maximum number of iterations set is exhausted.

As a proof of concept, we now use the algorithm described to compute the ground states for a linear array of *N* = 4 excitons in Eq. (2) for different values of (*D*_{1}, *D*_{2}) as indicated in Fig. 2(b) under the PBC. The corresponding ground state energies and their associated errors are displayed in Figs. 6(a) and 6(b), respectively. The energy errors are defined with respect to the results of the exact diagonalization (ED) with a convergence threshold of 2 × 10^{−3}*ϵ*. Note that the choice of this threshold is motivated by the fact that it corresponds to a relative error percentage of $<0.1%$ as can be established using the absolute values of the energy in Fig. 2(a) and the obtained errors from Fig. 6(b). The corresponding infidelities of the many-body wavefunctions obtained from the RBM network are plotted in Fig. 6(c). These infidelities are defined as $1\u2212F(\psi RBM,\psi ED)$, where $F(\psi RBM,\psi ED)=|\u27e8\psi RBM|\psi ED\u27e9|2$, and can only attain a value of zero if the states from *ED* and *RBM* are the same. Since the infidelities of the *RBM* wavefunctions are extremely low, this indicates that the state are nearly alike with the ones from *ED* at all points *A*–*I* and hence any observable property when computed from these states (including but not limited to one-body properties or many-body properties, such as correlation functions) would be in good agreement. We also present the Néel and string order correlation functions obtained from these states in Sec. S5 of the supplementary material. The circuit described has been simulated using the *Statevector Simulator* backend in Aer provider in Qiskit, which corresponds to IBM’s Quantum Information Software Kit (Qiskit).^{104} The details about mapping the system from psuedospin-1 degree of freedom to spin-$12$ for training the network *G* = (*V*, *E*) are presented in Sec. S3 of the supplementary material. The mapping techniques can also be used for quantum simulation of dynamical systems with several coupled potential energy surfaces for which a generalized mapping to spin systems has been derived and discussed recently.^{105} For *N* = 4, the total number of configurations accessible to the state-space would be *d* = 3^{4} = 81, which would require *n* = 7, *m* ∼ 7, and an additional ancillary qubit for the simulation of the circuit as described in the algorithm above. However, as the ground state is a singlet, symmetry restrictions allow us to disregard non-conforming configurations tapering the number of qubits to (*n* = *m* = 5). The total circuit width is 11, which is already significantly higher than the usual simulation requirements with the same number of spins in *S* = $12$ case. For all the (*D*_{1}, *D*_{2}) points studied, we have seen remarkable agreement with the exact results and DMRG calculations as presented in the above figure.

## V. CONCLUSION

In this work, we studied different topological and magnetic phases of a spatially modulated 1D chain of spin-1 particles with the most general nearest–neighbor interactions realizable with optically active Rydberg *p*-excitons in Cu_{2}O. Different phases have been categorized based on local spin magnetization and the von Neumann entropy as well as their pair-wise Néel and long-range string correlations. By studying the phases as a function of the local modulation, we examined both the phase transitions as well as their robustness to local perturbations. Later, we proposed a quantum algorithm that allowed us to simulate the many-body ground state of such a system on a quantum computer using a hybrid approach based on RBM. As mentioned, strongly interacting Rydberg excitons in Cu_{2}O are promising candidates for realizing analog quantum simulators in a solid-state platform. The ability to control the interaction strength is critical for investigating many-body systems, and in this case, the interaction mediated via Rydberg excitons can be controlled with the inter-particle spacing *R*_{0} as well as the principal quantum number of the state, *n*. Therefore, Rydberg excitons provide a highly controllable and tunable tabletop solid-state platform with direct optical access to the individual sites for coherent control as well as readout. Since the excitons are inherently out-of-equilibrium particles, an important follow-up study is to investigate the effect of drive and dissipation on the results presented in this work. Furthermore, it would be interesting to explore the implications of such phases on the photoluminescence spectra of the exciton chain to examine whether the attainable magnetic and topological phases can be directly studied by laser spectroscopy techniques. From the theoretical point of view, geometric analysis of the phase transitions associated with some of the phases in this work may also be undertaken.^{106} Another interesting direction is the study of long-range interactions readily achievable with Rydberg excitons in such spin systems. It would be interesting to explore the generalized models with interactions beyond the nearest neighbor and study their effects on the magnetic and topological phases obtained in this work.

## SUPPLEMENTARY MATERIAL

Supplementary material consisting of five sub-sections is available along with this paper. In Sec. S1, we describe the symmetries of the model like permutation, which is a *Z*_{2} symmetry of the model. Appropriate angular momentum symmetry, which has been used for qubit tapering, is also discussed. In Sec. S2, we present the explicit amplitudes and phases of the various components of the wavefunction in the configuration basis obtained for a linear chain of *N* = 4 excitons in various phases—anti-ferromagnet, Haldane phase, and large-|*D*| disordered phase. In Sec. S3, we describe two methods for mapping a spin-1 system to spin-$12$ for training the network *G* = (*V*, *E*). In Sec. S4, we present an explicit description of the quantum circuit for the hardware implementation of the quantum algorithm. In Sec. S5, we present the Néel and string order correlation function obtained by training the network *G* = (*V*, *E*) using the quantum algorithm and compare it against exact diagonalization.

## ACKNOWLEDGMENTS

The authors acknowledge the fruitful discussions with Stefan Scheel and Valentin Walther. S.K. and M.S. would like to acknowledge the financial support of the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, and Quantum Science Center. H.A. acknowledges the Purdue University Startup fund. We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors and do not reflect the official policy or position of IBM or the IBM Quantum team.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Manas Sajjan**: Conceptualization (lead); Data curation (lead); Formal analysis (equal); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (lead). **Hadiseh Alaeian**: Formal analysis (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Sabre Kais**: Formal analysis (equal); Funding acquisition (equal); Supervision (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

_{2}× Z

_{2}symmetry breaking in Haldane-gap antiferromagnets

_{2}O

_{2}O

_{2}O

_{2}O from strain traps

_{2}O