Inverse Hamiltonian design of highly-entangled quantum systems

Solving inverse problems to identify Hamiltonians with desired properties holds promise for the discovery of fundamental principles. In quantum systems, quantum entanglement plays a pivotal role in not only characterizing the quantum nature but also developing quantum technology like quantum computing. Nonetheless, the design principles of the quantum entanglement are yet to be clarified. Here we apply an inverse design framework using automatic differentiation to quantum spin systems, aiming to construct Hamiltonians with large quantum entanglement. We show that the method automatically finds the Kitaev model with bond-dependent anisotropic interactions, whose ground state is a quantum spin liquid, on both honeycomb and square-octagon lattices. On triangular and maple-leaf lattices with geometrical frustration, it generates numerous solutions with spatially inhomogeneous interactions rather than converging to a specific model, but it still helps to construct unprecedented models. The comparative study reveals that bond-dependent anisotropic interactions, rather than isotropic Heisenberg interactions, amplify quantum entanglement, even in systems with geometrical frustration. The present study paves the way for the automatic design of new quantum systems with desired quantum nature and functionality.

Much of the existing research has focused on assessing quantum entanglement measures of specific quantum systems.Yet, for applications like the synthesis of quantum materials with unique properties and the development of devices suitable for quantum computation, the primary requirement is designing systems with the desired entanglement properties.Nevertheless, there remains a limited understanding of the design principles of quantum entanglement.
An efficient way to obtain the systems that exhibit the desired properties is to solve the inverse problem.Thus far, numerous methods have been proposed for such inverse design.In particular, for quantum systems, several techniques were developed, including the quantum state preparation [26] and the parent Hamiltonian identification with specific wave functions [27][28][29][30][31][32][33].In these methods, however, it is not straightforward to directly derive systems from tangible physical quantities or entanglement measures.Machine learning-driven approaches, such as the Bayesian optimization [34], the generative models [35], and the random forests [36], have also been explored.While these methods are potentially powerful, they require extensive data collections and face challenges in ensuring precision.An alternative method using automatic differentiation, which was recently developed by the authors [37], eliminates the need for data and addresses the challenges mentioned above.Thus, it can be a versatile tool for designing the quantum nature of the system like quantum entanglement, but it has been applied only to quantum systems that can be dealt with independent particle approximations.
In this study, we develop a method to design Hamiltonians of quantum many-body systems characterized by large quantum entanglement, leveraging the inverse design framework utilizing automatic differentiation [37].While the framework itself is versatile, the straightforward application by using, e.g., the EE, does not work efficiently, leading to instability of the optimization.To circumvent the difficulty, we develop two techniques: (i) We define and use the thermal ensemble of entanglement entropy (TEEE) as the objective function, and (ii) we incorporate a symmetrization in the optimization procedure such that the result is spatially homogeneous.It is noteworthy that our method can be extended to other entanglement measures such as logarithmic negativity (LN) [38], mutual information (MI) [39], and topological entanglement entropy [19,20].In addition, it can be adapted to design systems with predetermined entanglement properties instead of maximizing them.
We apply this method to quantum spin models with two-spin interactions on various lattice structures with and without geometrical frustration.For the models on the honeycomb and square-octagon lattices, both of which are free from geometrical frustration, we show that the optimized results are equivalent to the Kitaev model [40], a renowned exactly-solvable quantum spin liquid model and a potential candidate for topological quantum computation.In contrast, for the models on the triangular and maple leaf lattices with geometrical frustration, the optimization does not settle on a specific model, and instead, it yields numerous solutions with spatially inhomogeneous interactions, which show similar quantum entanglement.Even in these cases, however, we demonstrate that the present method is useful to find a homogeneous Hamiltonian with large entanglement; as an example, we obtain a model only with anisotropic interactions on the maple-leaf lattice.This underscores the proficiency of our method in autonomously generating models with large quantum entanglement.
Furthermore, we investigate the characteristics of interactions that are dominant in the obtained models with large entanglement.Our findings emphasize that systems with bond-dependent anisotropic interactions like the Kitaev-type interaction are more inclined to have enhanced entanglement, irrespective of the presence or absence of geometrical frustration.For instance, on the triangular lattice, large entanglement is achieved not by isotropic Heisenberg interactions expected for the resonating valence bond state [41], but by strongly anisotropic interactions depending on the bonds.This sheds light on the intricate relationship between quantum entanglement and the form of exchange interactions.Advancing and refining this method could pave the way for the discovery of systems that are pivotal for cutting-edge applications, such as quantum computing.
This paper is organized as follows.In Sec.II, we introduce the method.
After introducing the general framework, we propose the objective function in Sec.II A and describe the actual calculation flow in Sec.II B. In Sec.III, we present the results for quantum spin systems on honeycomb (Sec.III A), square-octagon (Sec.III B), triangular (Sec.III C), and maple-leaf lattices (Sec.III D).We discuss the results for different lattices from the view point of the ratio of anisotropic interactions in Sec.IV.Finally, we give the concluding remarks in Sec.V.

II. METHOD
In this section, we introduce a method to obtain the model with large quantum entanglement, based on the framework of inverse Hamiltonian design using automatic differentiation [37].In this framework, the Hamiltonian is parameterized by a number of parameters θ as H(θ).We also construct an objective function, L(θ), specifically designed such that its minimum corresponds to the realization of a desired state.Utilizing the gradient ∂L(θ)/∂θ, efficiently computed through automatic differentiation, we iteratively update θ using the gradient descent method to minimize L(θ).This facilitates the design of a Hamiltonian having the desired proper-ties.In this study, we apply the method to the design of quantum many-body systems with large quantum entanglement.In the following subsections, we detail the objective function and provide the comprehensive calculation flow.

A. Objective function
A prevalent measure of quantum entanglement is the EE, derived from the density matrix, ρ = |ψ⟩⟨ψ|, where |ψ⟩ represents the wave function of the systems [42].The EE between bi-partitioned subsystems A and B is given by with ρ A = Tr B ρ. Here, Tr A(B) indicates the partial trace over the subsystem A (B).We find, however, that straightforward application of the EE to the objective function in the framework of the inverse Hamiltonian design leads to instability in the optimization process.This is because the energy levels of the ground state and the first excited state are swapped during the optimization process, resulting in abrupt changes in the value of the objective function.Another problem is that the value of EE depends on how the system is bi-partitioned.
To address these issues, we define the TEEE, which aggregates the EE with the Boltzmann distribution, as where β is the inverse temperature, and E n (θ) and S n,g,ξ (θ) denote the energy and EE for each eigenstate n, respectively; g represents the groups of bipartitioning patterns associated with translational and rotational symmetries of the system, and ξ represents the bi-partition index within a group g.For instance, for a cluster with four sites {s 1 , s 2 , s 3 , s 4 } in a periodic onedimensional chain, there are two groups of bi-partitioning patterns, one with two elements and the other with one element, represented as (g = 1, ξ = 1) = (s 1 , s 2 |s 3 , s 4 ) and (g = 1, ξ = 2) = (s 2 , s 3 |s 1 , s 4 ), and (g = 2, ξ = 1) = (s 1 , s 3 |s 2 , s 4 ).Employing the TEEE mitigates the computational instability mentioned above, as it incorporates not merely the ground state but also the excited states.In addition, to resolve the problem that the EE depends on the bi-partitioning pattern, we construct the objective function L including a penalty to suppress the deviation of TEEE between different ξ within each group g as where with Here, N g denotes the number of bi-partitioning patterns within each group g, and λ ≥ 0 is a hyperparameter for the penalty term.The first term in Eq. ( 3) aims to enhance the overall TEEE, while the second term ensures uniformity of TEEE within the group.We also tried other objective functions to increase the MI or the LN instead of the TEEE, and found that the TEEE is most suitable for the present purpose; this is because it allows solutions with degeneracy between states connected by symmetry, whereas the MI and LN only allow nondegenerate solutions (see Appendix A).The overall calculation procedure is illustrated in Fig. 1.First, we define a Hamiltonian depending on parameters θ.Next, the Hamiltonian is diagonalized, yielding the energy E n (θ) and the wave function |ψ n (θ)⟩ for each eigenstate n.Then, we compute the density matrix as ρ n (θ) = |ψ n (θ)⟩⟨ψ n (θ)| for every state n.Subsequently, the EE S n,g,ξ (θ) is computed for each state with respect to each bi-partitioning pattern g, ξ.Through the calculation of S T g,ξ in Eq. ( 2), we compute the objective function L(θ) as given in Eq. (3).Once L(θ) is obtained, its gradient with respect to θ, ∂L(θ)/∂θ, is computed using automatic differentiation.Then, θ is updated to minimize L(θ) using ∂L(θ)/∂θ by the gradient descent method.Iteratively applying this procedure, we are able to construct a Hamiltonian, H(θ opt ), with large and spatially uniform quantum entanglement, where θ opt is parameters after the optimization.
While our approach is applicable to arbitrary Hamiltonians, we focus on a quantum spin Hamiltonian with two-spin exchange interactions between spin- 1  2 moments, which is given by where µ spans {x, y, z}, i, j are indices designating individual sites, and σµ i represents the Pauli operator at site i.We parametrize the exchange constants by θ = {θ µ ij } as which satisfies the condition µ (J µ ij ) 2 = 1.We note that the parametrization is not unique, and the optimized results may depend on the way of parametrization as well as the objective function; see Appendix B. While the present study deals with only the diagonal interactions in spin space, it can be straightforwardly extended to include off-diagonal components such as J xy ij , as well as other forms of interactions rather than the two-spin exchange interactions for arbitrary spin length S.
We apply the following numerical techniques to make the optimization better and more efficient.First, a tiny random hermitian matrix, whose element is in the order of 10 −8 , is added to the Hamiltonian to avoid the instability caused by degeneracy in the energy levels.This step is crucial because the current version of the JAX library lacks support for the automatic differentiation of diagonalization when involving degeneracy in eigenvalues.The minor modification slightly lifts the degeneracy, thereby stabilizing the optimization process.See Appendix C for the details.Next, in the calculation of S T g,ξ , only N states from the lowest energy are used out of 2 N S total states for the number of spins N S .It is because high energy states do not contribute to the value of TEEE if β is large enough.We set N = 40 throughout the paper; increasing N has a negligible effect on the results.To compute the EE, we use the relation Trρ A log ρ A = i γ i log γ i , where {γ i } represents the eigenvalues of ρ A .For the gradient decent method, we employ AdaBelief optimizer [43] with hyperparameters β 1 = 0.9 and β 2 = 0.999.The hyperparameter λ and the learning rate η are adjusted in accordance with optimization steps.
We conduct calculations with 100 different initial conditions.In case the optimized Hamiltonians have spatially inhomogeneous interactions, we may further optimize to purify the Hamiltonian with another objective function, as shown in Appendix D. Automatic differentiation is performed using the JAX library in Python [44], and optimization is carried out with the optax library.All computations are executed on an NVIDIA A100 GPU.

III. RESULTS
We apply the method to the model in Eq. ( 7) on four different lattices: the honeycomb lattice, the squareoctagon lattice, the triangular lattice, and the maple leaf lattice.The former two are bipartite, while the latter two include triangles, i.e., geometrically frustrated.We present the results on each lattice one by one in the following sections.

A. Honeycomb lattice
First, we apply the method to the honeycomb lattice consisting of six-site cluster with periodic boundary conditions, illustrated in Fig. 2(c).This lattice comprises nine distinct bonds and 27 parameters.We restrict our consideration to bi-partitions where the system is divided into two interconnected subsystems, with each subsystem having three sites, as detailed in Table I.In this case, we have only one group which contains nine different elements.The inverse temperature is maintained at β = 60 throughout the calculation.3) and the learning rate η are shown in the lower panel.We note that L is reduced faster by taking λ nonzero on the way of the optimization than by taking it nonzero from the beginning.The inset shows the changes of ST and ∆S T ; ST reaches around 1.7, while ∆S T reaches a sufficiently small value.This suggests that the quantum entanglement is optimized with being almost identical for all bipartition patterns.Note that ∆S T remains larger when optimizing with λ = 0, leading to the spatially inhomogeneous interactions.
Figure 2(b) displays the evolution of J µ ij during the optimization process.One can observe that they eventually converge to one of the values 0, 1, or −1.Since each bond satisfies the condition µ (J µ ij ) 2 = 1, this means that only one of the µ ∈ {x, y, z} becomes nonzero for each J ij , resulting in an Ising-type anisotropic interaction.The optimized interactions of each bond are plotted by color on the lattice in Fig. 2(c), with the color code in the inset.Note that we plot the absolute values of the optimized J ij since their signs do not matter to the value of the objective function; we confirm that changing all these signs to negative does not change the value of the objective function.Interestingly, the interactions are all Ising-type and bond dependent, realizing the Kitaev model, which is known to exhibit spin liquid behavior with fractional Majorana excitations enabling topological quantum computations [40].
In Fig. 3(a), we present the optimization process of L for 100 different initial conditions.The histogram of the values of L at the end of the optimization is shown in the panel on the right.For comparison, we also show the value of L for the antiferromagnetic (AFM) Heisenberg model with J µ ij = 1/ √ 3 for all ⟨ij⟩ and µ, which is expected to show an AFM long-range order on this bipartite lattice.We find that all the 100 results have much smaller values of L than this value.We also find that 57 out of 100 converge to the Kitaev model with the smallest value of L. In addition, there is a gap between the smallest value of L and other values, indicating that the solutions of the Kitaev model are uniquely selected out.
Thus, the inverse design by maximizing the TEEE automatically finds the celebrated Kitaev model on the honeycomb lattice.The Kitaev model has the so-called bond frustration, where the local energy-optimal configuration cannot be satisfied on the whole due to the bonddependent Ising-type anisotropic interactions [45].This makes the wave function spatially nonlocal, leading to large quantum entanglement in the spin-liquid ground state.

B. Square-octagon lattice
Next, we apply the method to the square-octagon lattice, which is 1/5-depleted square lattice, as illustrated in Fig. 4.This lattice structure is realized in real materials [46], and has also been studied as a surface code in quantum computing [47,48].We consider an eight-site  4) and ∆S T in Eq. ( 5).The lower panel shows the schedule of λ in Eq. ( 3) and the learning rate η.(b) Changes in each component of spin interactions.(c) Optimized interaction parameters.The shades denote six-site cluster unit cells.The color of the bonds represents each spin interaction according to the color code in the inset which shows the ratio of the elements of J, e.g., red denotes that J x = 1 and J y = J z = 0.While the optimized results in Fig. 2(b) have a mixture of both positive and negative values, we plot the absolute values in (c) since the value of L remains consistent when all these signs of J µ ij are entirely negative.
cluster with periodic boundary conditions, comprising 12 distinct bonds and 36 parameters.We restrict our consideration to bi-partitions where the system is divided into two interconnected subsystems, with each subsystem having four sites, as detailed in Table II.In contrast to the honeycomb case in Table I, there are five groups, each with one to four elements.The inverse temperature is maintained at β = 160, which is larger than that for the honeycomb lattice case, because the larger cluster size increases the number of eigenstates.
In Fig. 3(b), the optimization process of L is depicted for 100 different initial conditions.The value of L of the AFM Heisenberg model is also plotted as a comparison, which has a lower value than the honeycomb lattice case and even lower than some optimized solutions.However, akin to the honeycomb lattice case, out of the 100 solutions, 10 converge to the minimal L value, distinctly separated from the other solutions by a discernible finite gap.
Figures 4(a) and 4(b) illustrate two of these optimal solutions.The color scheme is the same as the inset of Fig. 2(c).Remarkably, all the bonds shown in both cases exhibit bond-dependent Ising-type anisotropy.The solution in Fig. 4(a) is nondegenerate, while that in Fig. 4(b)  is doubly degenerate.Note that we modify the signs of J µ ij in Figs. 4, following the rules that we found: In Fig. 4(a), the value of L remains consistent whether all signs of bonds are uniformly negative or uniformly positive.Conversely, in Fig. 4(b), L remains the same value if an odd number of the four bonds that form a square are positive, with the remaining bonds being negative; the signs of any other bonds do not affect the value of L.
We note that some bonds are not completely anisotropic, e.g., J 36 ∼ 0.96 in Fig. 4(a).However, using the purification method outlined in Appendix D, they are transformed into fully anisotropic Hamiltonians with keeping the value of L almost intact.The remaining eight optimal solutions are almost identical to either of the two solutions.
The result in Fig. 4(a) is known as the Kitaev model, which manifests a spin liquid in the ground state [49,50].It is understood that this has large entanglement due to the bond frustration, like the honeycomb lattice case.Note that the solution in Fig. 4(b) is different from the Kitaev model, but it also has bond frustration due to the anisotropic interactions with different signs in the fourbond squares as mentioned above.

C. Triangular lattice
Here, we apply the method to the triangular lattice known to have geometrical frustration.The quantum spin models with geometrical frustration have been studied as a possible route to spin liquid states, such as the resonating valence bond state [51].We consider eight-site cluster with periodic boundary conditions as illustrated in Fig. 5.This comprises 24 distinct bonds and 72 parameters.We restrict our consideration to bi-partitions where the system is divided into two interconnected sub- systems, with each subsystem having four sites, as detailed in Table III.The inverse temperature is maintained at β = 160.
Figure 3(c) shows the optimization processes for 100 different initial conditions.There is only one solution that gives the minimum value of L, and several different solutions appear without showing an apparent gap, forming rather continuous spectrum.This is in stark contrast to the cases of honeycomb and square-octagon lattices where multiple results starting from different initial conditions converge to the specific optimal solutions well separated from the others.In addition, the value of L for the AFM Heisenberg model is close to the lowest L solution, which is also in contrast to the previous two cases.The ground state of the AFM Heisenberg model on the triangular lattice is expected to exhibit a noncollinear 120 • order [52,53], where the magnetic moment ∼ 0.20 [53,54] is reduced compared to that for the honeycomb lattice case ∼ 0.27 [55].This suggests larger quantum entanglement in the triangular lattice case than the honeycomb lattice case.Notably, our method can generate models with larger entanglement than the AFM Heisenberg model.In Fig. 5, we plot the optimized interactions for the solution with the lowest L value in Fig. 3(c).The interactions are spatially inhomogeneous without showing any clear pattern, but the anisotropy is relatively prevalent.
The lack of convergence to a single solution starting from different initial conditions in the triangular lattice case can be attributed to possible coexistence of two different types of frustration, namely, geometrical frustration and bond frustration.In contrast to the bipartite honeycomb and square-octagon with no geometrical frustration, the triangular lattice is non-bipartite and can exploit both two frustrations, for which multiple metastable solutions may arise rather than converge to a single solution.Indeed, the representative solution in Fig. 5 appears to reconcile two frustrations by adopting the spatiallyinhomogeneous anisotropic interactions.

D. Maple-leaf lattice
Finally, we apply the method to the maple-leaf lattice, which is a 1/7-depleted triangular lattice, and hence, has a slightly weaker geometrical frustration than the triangular lattice.Quantum spin systems with this lattice structure are realized in real materials [56,57].We consider a six-site cluster with periodic boundary conditions, as illustrated in Fig. 6.This lattice comprises 15 distinct bonds and 45 parameters.We restrict our consideration to bi-partitions where the system is divided into two interconnected subsystems, with each subsystem having three sites, as detailed in Table IV.The inverse temperature is maintained at β = 60.As depicted in Fig. 3(d), there is only one solution with the smallest L, without showing a clear gap to the FIG. 6.
Representative optimized interaction parameters for the maple-leaf lattice model in the solutions with (a) L ≃ −1.738 and (b) L ≃ −1.708.The shades denote sixsite cluster unit cells.The notations are common to those in Fig. 4. In (a), the signs appear to be random, while in (b) the signs are all taken to be negative though the optimized solutions have a mixture of both positive and negative, since the value of L remains consistent in both cases.
other solutions, among the calculations for 100 initial conditions.This is similar to the triangular lattice case, but also different from the results of the honeycomb and square-octagon lattices, where multiple calculations converge to a single optimal solution.We plot the result for the smallest value of L ≃ −1.738 in our calculations in Fig. 6(a).It appears that the interactions are asymmetric and spatially inhomogeneous.Furthermore, the signs of the interactions do not show a specific spatial pattern.This may be due to the fact that, similarly to the triangular lattice case, both geometrical frustration and bond frustration are allowed to coexist.Meanwhile, L of the AFM Heisenberg model is not that small, as shown in Fig. 3(d), presumably because the maple-leaf lattice has weaker geometrical frustration than the triangular lattice.
We also plot the solution with a larger L ≃ −1.708 in Fig. 6(b).Contrary to the other solutions, it has a unique structure constructed solely by Ising-type anisotropic interactions.The structure includes the quadrangle by four J x (J y , J z ) enclosing J y (J z , J x ), and is symmetric for not only C 2 rotation but also operations that simultaneously perform a C 3 rotation in real space and a 2π/3 rotation about the [111] axis in spin space.Although the results in Fig. 6(b) have a mixture of both positive and negative values for the value of J µ ij , we confirm that changing all these signs to negative does not change the value of L. Thus, although this solution is not the best one in terms of quantum entanglement, it gives a concise model with bond-dependent Ising-type interactions compatible with the structural unit cell.To the best of our knowledge, such a model has not been reported thus far.This demonstrates that our framework is useful to develop an unprecedented, interesting model.

IV. DISCUSSION
In Sec.III, we showed that our method successfully finds a model with optimized TEEE for each lattice.Interestingly, the solutions look qualitatively different between the bipartite lattices without geometrical frustration (honeycomb and square-octagon) and the nonbipartite lattices with geometrical frustration (triangular and maple-leaf).For the former, the optimization processes starting from different initial conditions converge onto a single solution with strongly anisotropic Kitaevtype interactions, while for the latter, they do not appear to single out a unique solution.Nevertheless, the interesting observation is that, even for the latter geometrically frustrated cases, the obtained solutions are not isotropic but rather anisotropic in spin space; the optimized TEEE has a lower value than that for the isotropic AFM Heisenberg model.This suggests that the systems try to maximize the quantum entanglement by exploiting the bond frustration even in the presence of geometrical frustration.
For more comparison, we estimate the degree of spin anisotropy in each lattice by computing the ratio of the anisotropic interactions by where f (x) = 1 for |x| > x thres and otherwise f (x) = 0 (we take the threshold x thres = 0.95 in the following calculations); N b is the number of distinct bonds in each cluster used in the calculations, and the average ⟨• • • ⟩ is taken over solutions with L in the range satisfying L < 0.98L min among the 100 calculations, where L min is the lowest value of L. The results are plotted in Fig. 7.
For the honeycomb and square-octagon lattices that are FIG.7. Changes in the ratio of anisotropic spin interactions, rani in Eq. ( 9), during the optimization process on each lattice.
free from geometrical frustration, the values of r ani are rapidly increased; r ani reaches almost 1 in the honeycomb lattice case, and it exceeds 0.95 in the square-octagon case.In contrast, the increase of r ani is less pronounced in the triangular and maple-leaf lattices with geometrical frustration; r ani is about 0.4 in the triangular lattice case, and it once increases but decreases down to ∼ 0.3 in the maple-leaf lattice.Nonetheless, the final values of r ani are nonzero and larger than those for the initial states with random J µ ij , indicating that the optimization of TEEE favors anisotropic interactions even for the geometrically frustrated cases.This is an important lesson from the automatic design: To optimize the quantum entanglement, it would be, in general, helpful to utilize anisotropic interactions.This is rather counterintuitive, since the anisotropic interactions suppress spin-singlet formation that is widely recognized to be crucial to stabilize quantum spin liquids like the resonating valence bond state.

V. CONCLUSION AND OUTLOOK
In summary, we have developed a method to design the Hamiltonian with large quantum entanglement using automatic differentiation.In the method, the parameters in Hamiltonians are optimized so as to reduce the objective function corresponding to the desired entanglement property of the system.Among several quantities representing the quantum entanglement, we introduced an extension of the entanglement entropy, which we call the thermal ensemble of entanglement entropy, as the objective function, to achieve stable and efficient optimization.Applying this method to the quantum spin systems on four different lattice structures, we showed that it successfully generates the Hamiltonian with large quantum entanglement.The optimization on the honey-comb and square-octagon lattices automatically finds the Kitaev model, which is known to provide exactly-solvable spin liquid ground states.On the triangular and mapleleaf lattices with geometrical frustration, the results do not converge to a specific model, suggesting a degenerate manifold of the models with spatially inhomogeneous interactions.Through the optimization, however, we found an unprecedented homogenous model with large entanglement for the maple-leaf lattice case.Furthermore, by comparing optimized interactions on these different lattices, we found that the introduction of anisotropic interactions plays an important role in the enhancement of quantum entanglement even in the presence of geometrical frustration.This finding is rather incompatible with the conventional wisdom that the spin-singlet formation by isotropic Heisenberg interaction is widely recognized as a key ingredient for the realization of quantum entangled states, such as the resonating valence bond state.
The optimized solutions in Sec.III A correspond to the gapless phase of the Kitaev model.This model is exactly-solvable and the ground state offers a quantum spin liquid with two types of fractional excitations, itinerant Majorana fermions and localized Z 2 fluxes [40].The entanglement entropy of this model was shown to have two separate contributions from the two fractional excitations [58].It is noteworthy that the framework automatically designs such an intriguing model with exact solvability.Meanwhile, the inverse design discovered new models on the square-octagon and maple-leaf lattices.These models may also have interesting properties of quantum entanglement.This is a subject for future research.
Our approach has a high flexibility, being effectively applicable to a diverse range of systems, including quantum spin systems, strongly correlated systems for both fermionic and bosonic cases, and even non-Hermitian systems.Nonetheless, an increase in the number of parameters tends to complicate the process of attaining optimal solutions.To address this issue, it might be necessary to develop new strategies for the efficient sampling of initial conditions, potentially employing techniques like Bayesian optimization.However, it is important to note that finding a globally optimal solution is not a prerequisite for the successful identification of unprecedented models, as exemplified by the finding of new models on the square-octagon and maple-leaf lattice cases.Therefore, despite the potential challenges in resolving global optimal issues, the current approach retains substantial utility for the discovery of new models in various complex systems.
In this research, the exact diagonalization was employed to calculate the entanglement entropy.However, the automatic differentiation of exact diagonalization is notably memory-intensive.When utilizing a single NVIDIA A100 GPU machine with 40 GB of memory, the maximum manageable matrix size is 2 12 × 2 12 .Given the versatility of automatic differentiation, it can be integrated with the Lanczos method or tensor networks to accommodate larger system sizes.It is noteworthy that the implementation of automatic differentiation in tensor networks has been documented recently [59].
Exploring applications to diverse entanglement measures promises to be intriguing.Maximizing of topological entanglement entropy could pave the way for novel topological quantum systems.For instance, applying this to the out-of-time-ordered correlation [60][61][62] and tripartite mutual information [63] might unveil distinctive quantum scrambling behaviors.Historically, the conception of models and systems has predominantly relied on the experience and intuition of researchers.The automation of system construction, possessing the target quantum functionality through an inverse design approach, would facilitate the uncovering of innovative principles and systems.The Hamiltonian and its bi-partitions are the same as in Sec.III A. The objective function L is given by Eq. ( 3) with the replacement of S T g,ξ with I g,ξ in Eqs. ( 3)-( 6), where I g,ξ is the MI for bi-partition with group g and index ξ.Ī and ∆I are defined accordingly.We set β = 60.
Figures 8(a) shows the optimization process of L for the calculations with 100 different initial conditions for the model in Eq. ( 7) on the honeycomb lattice.While the optimized values of L exhibit similar magnitudes, they do not converge to a specific optimal state, which stands in stark contrast to the results for the TEEE in Fig. 3(a), where more than half of the solutions converge to the same optimal state.Furthermore, we note that the magnitudes of ∆I in this scenario, being on the order of 10 −2 as shown in the inset of Fig. 8(a), are considerably larger than those observed in the case of the TEEE, indicating that the optimized results are not spatially homogenous.Figure 8(b) depicts the optimized interaction parameters for the solution with the smallest value of L. The result breaks both translational and rotational symmetry.We confirmed that most of the other solutions also have low symmetry compared with the result obtained in the optimization using the TEEE in Sec.III A.
Figure 8(c) displays a histogram of the populations of the ground state and the first excited state in the 100 optimized results.The population of a state |ψ⟩ is defined by ⟨ψ|ρ B |ψ⟩.The ground state clearly dominates the population distribution.This is attributed to the term S B tot in Eq. (A1), which tries to make the thermal density matrix as close to a pure state as possible and avoid degeneracy in the ground state during the optimization process.Such a trend is in stark contrast to the optimization process for the TEEE that allows degeneracy.This is the reason why the TEEE has an advantage over the MI in the current problem; the interaction parameters with high spatial symmetry are unlikely to emerge when attempting to increase the MI because it avoids degeneracy.Next, we discuss the results for another measure of quantum entanglement for mixed states, the LN, which is defined by where || • || denotes the trace norm of the matrix, and ρ B TB is the partial transpose over the subsystem B of ρ B .In the numerical calculation, ||ρ B TB (β)|| is computed by summing up the absolute values of the eigenvalues obtained by diagonalizing ρ B TB .Since the LN is an entanglement monotone [64], maximizing it can create a system with large entanglement.The optimization is performed similarly to the MI case above.
Figure 9 summarizes the results by using the LN in a similar manner to Fig. 8.The results are qualitatively similar to those obtained by the MI; the optimization does not single out a specific optimal state from the others, the interaction parameters break the lattice symmetry, and the population is dominated by the ground state.Thus, we conclude that the LN is inferior to the TEEE as well as the MI.8) and (b) Eq. (B1).The notations are common to those in Fig. 4. The thickness of the bond corresponds to the absolute value of its interaction.See the text for details.
in Sec.II B, the parametrization in Eq. ( 8) is not unique.Instead, one can take, for instance, where N b is the number of distinct bonds in each cluster used for the calculations.Note that Eq. ( 8) satisfies the local constraint on each bond, µ (J µ ij ) 2 = 1, while Eq.(B1) satisfies the global constraint, ⟨i,j⟩µ (J µ ij ) 2 = N b .Hence, Eq. (B1) allows different bond interactions to take different magnitudes; namely, it allows the optimization in wider parameter range than in Eq. ( 8).
The optimized results depend on not only the parametrization but also the objective function L. Indeed, we find that the results are almost intact for the parameterizations in Eq. ( 8) and Eq.(B1) when taking the TEEE in Eq. (3) for L, but depend on the parameterization when taking the LN in Eq. (A6) for L. We demonstrate this below.Firstly, we perform calculations based on the parametrization in Eq. ( 8), by using the objective function L(θ) defined as where E N g,ξ is the LN for the bi-partition labeled by g and ξ (see Table I).Note that here we set λ = 0 in Eq. ( 3), meaning different values of the LN are allowed for different bi-partition ξ of the same group g.

Appendix C: Computational instability
In this Appendix, we investigate the computational instability regarding the automatic differentiation of diagonalization.The JAX library does not currently support the automatic differentiation of diagonalization when involving degeneracy in the eigenvalues.To circumvent this issue, we add a tiny random Hermitian matrix, denoted as δH rand , to the Hamiltonian.This minor modification slightly lift the degeneracy.Each element of δH rand is generated from a Gaussian distribution with a mean of 0 and a variance of δ rand .In the optimization process, we use a different random matrix at each step, thereby facilitating more robust calculations.
Let us demonstrate the instability using a two-spin model as an example.The Hamiltonian is defined by Here, we use the LN discussed in Appendix A as the objective function, since the TEEE for the two-spin model takes almost always the maximum value for any interaction.Note that, however, the following arguments are applicable to any objective functions that requires the numerical diagonalization.
Figure 11(a) shows the optimization process of the objective function L = −E N , with the interaction parameters {J x , J y , J z } in the inset.We set β = 20 and δ rand = 10 −10 .Note that ĒN = E N since there is only one way for bi-partition in the two-spin problem.In this case, the LN converges to log 2, and the interaction parameters converge J x = J y = J z = 1 √ 3 ; namely, the optimization leads to the AFM Heisenberg model with isotropic interactions.In Fig. 11(b), we plot the deviation of the ground state energy, E 0 , from that of AFM Heisenberg model, E exact , while changing δ rand .It clearly shows that the deviation is proportional to δ rand , while the calculations for δ rand < 10 −15 becomes unstable as in- dicated by red shade.Figure 11(c) shows the changes in the energy differences between the adjacent states when δ rand = 0; here E i is the energy of the ith excited state.It shows that computational instability arises when the energy differences, for instance, E 3 − E 2 and E 2 − E 1 , are reduced below ∼ 10 −16 .This indicates that the optimization breaks down without δH rand due to the degeneracy.For the stable computation, we take δ rand in the order of 10 −8 in the calculations in the main text as well as other Appendices.where (D2) The smaller the parameter entropy, the fewer nonzero components in the interaction parameters.We optimize the parameter entropy to be smaller after the optimization of L using TEEE; namely, we perform an additional optimization using the following objective function, where λ is the hyperparameter, ReLU denotes the rectified linear unit function, L opt is the final value L after the optimization of the TEEE, and L(θ) in the input for ReLU is the same as Eq.(3).By optimizing θ to minimize the objective function L(θ), we can obtain a simpler model without increasing the value of L.
First, we apply this purification to the results of maximizing the TEEE on the square-octagon lattice in Sec.III B. The interactions in Fig. 4(a) are not completely homogeneous, with some bonds exhibiting slight deviations from the ideal Kitaev model.Figure 12(a) shows the optimization process of the parameter entropy.It is successfully reduced to zero, indicating that each bond has only one nonzero component among {J x , J y , J z } as shown in Fig. 12(b).Figure 12(c) presents the interactions after the purification.Compared to Fig. 4(a), all bonds become completely anisotropic as in the ideal Kitaev model.This comparison confirms that minimizing the parameter entropy is useful to purify the models, yielding more interpretable and tractable models.
Next, we apply the method to the results of the triangular lattice in Sec.III C. As indicated in Fig. 13(a), the optimization process starting from the result in Fig. 5 leads to a reduction in the value of S param , although it does not achieve a value of zero, which is different from the square-octagon lattice case.The optimization process for J µ ij in Fig. 13(b) shows that the interactions are not homogeneous, while there is a change in the values of J µ ij .The interactions after the purification lack any clear pattern as illustrated in Fig. 13(c).These results indicate that this purification method using the parameter entropy does not always successfully reduce the inhomogeneity of the model.The development of alternative objective functions for the purification may yield more robust results; it will be a focus for future research.

B. Calculation flow and computational details FIG. 1 .
FIG. 1.Flowchart of the inverse design of Hamiltonians with desired quantum entanglement properties using automatic differentiation.See the main text for details.
partition patterns of the six-site cluster of the honeycomb lattice into two three-site subsystems A and B for calculating the EE.The indices in A and B correspond to the site numbers in Fig.2(c).

Figure 2 (
Figure 2(a) shows the optimization process of the objective function L. The changes in λ in Eq. (3) and the learning rate η are shown in the lower panel.We note

FIG. 2 .
FIG. 2. A representative optimization result for the honeycomb lattice model.(a) Changes in the objective function L through the optimization process.The inset shows changes in ST in Eq. (4) and ∆S T in Eq. (5).The lower panel shows the schedule of λ in Eq. (3) and the learning rate η.(b) Changes in each component of spin interactions.(c) Optimized interaction parameters.The shades denote six-site cluster unit cells.The color of the bonds represents each spin interaction according to the color code in the inset which shows the ratio of the elements of J, e.g., red denotes that J x = 1 and J y = J z = 0.While the optimized results in Fig.2(b) have a mixture of both positive and negative values, we plot the absolute values in (c) since the value of L remains consistent when all these signs of J µ ij are entirely negative.

FIG. 3 .
FIG. 3. Changes in L for 100 different initial conditions for (a) honeycomb, (b) square-octagon, (c) triangular, and (d) maple-leaf lattices.The right panel in each figure is the histogram of L after convergence.The dotted line indicates L for the AFM Heisenberg model on each lattice.Note that the abrupt changes at steps 1000 and 1300 are due to the changes in the learning rate η, shown below in Fig. 2(a).

FIG. 4 .
FIG. 4.Representative optimized interaction parameters for the square-octagon lattice model in the solution with (a) L ≃ −2.0070 and (b) L ≃ −2.0086, which are nondegenerate and doubly degenerate states, respectively.The shades denote eight-site cluster unit cells.The color of the bonds represents the value of each spin interaction according to the inset in Fig.2(c) and the solid (dotted) line indicates that the sign of the largest component of J µ ij is negative (positive).Note that though the optimized values of J µ ij have a mixture of both positive and negative values, we illustrate modified signs since the value of L remains consistent in both cases: (a) when all these signs are either entirely positive or entirely negative and (b) when an odd number of bonds forming the square are positive.

FIG. 5 .
FIG. 5. Representative optimized interaction parameters for the triangular lattice model in the solution with the lowest L ≃ −2.0451.The shades denote eight-site cluster unit cells.The notations are common to those in Fig. 4.
Bi-partition patterns of the six-site cluster of the maple-leaf lattice into two three-site subsystems A and B for calculating the EE.The indices in A and B correspond to the site numbers in Fig.6.

FIG. 8 .
FIG. 8. (a) Changes in L using the MI for 100 different initial conditions for the model on a honeycomb lattice.Note that the abrupt increase in L at the 300 step is due to the change in λ.The inset shows the changes in ∆I.The lower panel shows the schedule of λ and learning rate η.(b) Optimized interaction parameters.The notations are common to those in Fig. 4. (c) Histogram of the populations of the ground state and the first excited state for 100 solutions.

FIG. 9 .
FIG. 9. (a) Changes in L using the LN for 100 different initial conditions for the model on a honeycomb lattice.Note that the abrupt increase in L at the 300 step is due to the change in λ.The lower panel shows the schedule of λ and learning rate η.(b) Optimized interaction parameters.The notations are common to those in Fig. 4. (c) Histogram of the population of the ground state and the first excited state for 100 solutions.

FIG. 10 .
FIG. 10. (a),(b) The optimized interaction parameters for the honeycomb lattice model with the use of the objective function of Eq. (B2) under the constraint of (a) Eq. (8) and (b) Eq. (B1).The notations are common to those in Fig.4.The thickness of the bond corresponds to the absolute value of its interaction.See the text for details.
shows one of the optimized interactions with the smallest L(= − ĒN ) ∼ −1.78 among 30 calculations for the honeycomb lattice model.Six out of nine interactions are completely anisotropic, similar to the Kitaev model, while the remaining three take more isotropic values.Thus, this optimization procedure finds a model with a mixture of the Kitaev-type anisotropic interactions and the Heisenberg-like isotropic interactions.Next, we perform calculations based on the parametrization in Eq. (B1), keeping the other setup the same as above.Figure 10(b) shows the optimized interactions.Remarkably, in this case all the interactions are strongly anisotropic, showing different amplitudes of the interactions depending on the bonds.The amplitudes are represented by the thickness of each bond proportional to µ (J µ ij ) 2 .The optimized values are |J x 14 |, |J x 23 |, |J x 56 | ≃ 0.517 and |J y 16 |, |J y 25 |, |J y 34 |, |J z 12 |, |J z 36 |, |J z 45 | ≃ 1.168, and the other components of J µ ij are almost zero.This configuration corresponds to a spatially-anisotropic Kitaev model in which the J x bonds are weaker than J y and J z .Interestingly, the value of L for this solution is L(= − ĒN ) ∼ −1.80, which is smaller than the result with the parametrization using Eq.(8).This comparison demonstrates that better parametrization depends on the problem.Choices of the objective function and the model parametrization generate different Hamiltonians.Therefore, one should choose how to parameterize the interaction carefully according to the purpose.

FIG. 11 . 3 ,
FIG. 11.(a) Changes in the logarithmic negativity E N through the optimization process for the two-spin model.The inset shows changes in each component of spin interactions, J µ .We set η = 0.03 during the optimization.(b) Deviation of the ground state energy of the optimized result, E0, from the exact energy for J x = J y = J z = 1/ √ 3, Eexact = − √ 3/4, as a function of the variance of the random matrix elements, δ rand .(c) Changes in the energy difference between the lowenergy eigenstates.Ei is the energy of the ith excited state.In (b) and (c), the red area represents the region where the calculations become unstable.

FIG. 13 .
FIG. 13.(a) Changes in Sparam through the optimization process using Eq.(D3) starting from the result in Fig. 5 for the model on the triangular lattice.We set λ = 1000 and η = 0.03 during the optimization.(b) Changes in the interaction parameters.(c) Optimized interaction parameters.The notations are common to those in Fig. 5.

TABLE III .
Bi-partition patterns of eight-site cluster of the triangular lattice into two four-site subsystems A and B for calculating the EE.The indices in A and B correspond to the site numbers in Fig.5.