Field-induced Bose-Einstein condensation and supersolid in the two-dimensional Kondo necklace

The application of an external magnetic field of sufficient strength to a spin system composed of a localized singlet can overcome the energy gap and trigger bosonic condensation and so provide an alternative method to realize exotic phases of matter in real materials. Previous research has indicated that a spin Hamiltonian with on-site Kondo coupling may be the effective many-body Hamiltonian for $\text{Ba}_2\text{NiO}_2\text{(AgSe)}_2$ (BNOAS) and here we study such a Hamiltonian using a tensor network ansatz in two dimensions. Our results unveil a phase diagram which indicates the underlying phases of BNOAS. We propose, in response to the possible doping-induced superconductivity of BNOAS, a fermionic model for further investigation. We hope that our discovery can bring up further interest in both theoretical and experimental researches for related nickelate compounds.


I. INTRODUCTION
Bose-Einstein condensation (BEC) is a phenomenon where a finite fraction of the quasiparticles in the system condenses into a single quantum mechanical entity on a macroscopic scale as a consequence of quantum statistical effects [1,2]. It leads to an exotic phase of matter, superfluid, which is a fluid with zero viscosity. The superfluid was originally discovered in the liquid helium-4 [3]. On the other hand, it was found that the fraction of BEC is strongly suppressed due to the strong interaction between 4 He atoms. Since then, great effort has been devoted to the search for weakly interacting or dilute Bose gases.
Theoretical studies [18][19][20][21][22] also unveiled the underlying mechanism for such effects. In a spin-singlet material, the local singlet state serves as the ground state with a finite gap, separating itself from triplet excitons. Upon applying a magnetic field, the three-fold excited states split into three triplon bands with S z =+1, 0, −1. The branch with spins aligning along the external field becomes soft and gradually reduces the gap to complete closure [7,23]. After the gap closing, a condensation of triplons that breaks U(1) symmetry takes place, leading to the effective spin superfluid (SSF) phase.
For such effects to take place, clearly we need a spin-singlet phase as the precursor. A recent study by Jin et al. [24] of * wepickett@ucdavis.edu † hyunyong@korea.ac.kr a magnetic material with layered nickelate Ba 2 NiO 2 (AgSe) 2 (BNOAS), recently synthesized under high pressure and high temperature [25], proposed the origin of its peculiar susceptibility χ sp . It is constant above 150 K and the same constant below 110 K in zero field, with a peak at T * =130 K. Thus there are no free moments at high or low temperature, yet a magnetic reconstruction occurs at 130 K giving the peak in χ sp . This behavior can be explained as arising from local spin singlets (contributing nothing to susceptibility) at high and low T, with some reconstruction occurring at T * [24]. Correlated first principles calculations predict a ground state consisting of a novel singlet within the Ni e g subshell, made of spins with the local "Kondo-like" spin texture: the d x 2 −y 2 electron (or hole, depending on viewpoint) is coupled with the d z 2 electron (or hole) to an unusual spin singlet with internal orbital structure and highly anisotropic exchange coupling [24], with the first signature of such an on-site, orbital entangled Lee-Pickett singlet having been seen in calculations on the infinite layer nickelate LaNiO 2 [26]. Due to the abovementioned facts, in the studies of BNOAS' quantum effect one might need to focus on the area where spin singlet state serves as the ground state.
Jin et al. [24] proposed an effective spin Hamiltonian, named the Kondo sieve, for describing the spin behavior of BNOAS: The local spin moment σ i (from the d x 2 −y 2 orbital of BNOAS) is Kondo-coupled with the τ i spin moment (from the same-site d z 2 orbital) with exchange coupling K. σ i and τ i fulfill the commutation relation that [ξ a , ξ b ] = i abc ξ c with ξ = σ or ξ = τ , and abc being the Levi-Civita symbol.
Within a layer, for the nearest neighbor, denoted by i, j , σ moments are coupled by the Heisenberg J term. Due to the multilayered nature of BNOAS, neighboring NiO 2 layers have J z coupling between τ -spin neighbors [i, j] alongẑ. The isite total spin operator S i = σ i +τ i is coupled with the external field through a Zeeman field h. Note that the model enjoys global U(1) symmetry with e iα i S z i , for total spin and with an arbitrary α. Eq. (1) as a whole stands for the Kondo sieve model H KS , while its first part, for each layer, represents the Kondo-necklace model H KN in two dimensions (2D) [27]. For the on-site K term, we have ψ singlet |Kσ i ·τ i |ψ singlet = − 3 4 K and ψ triplet |Kσ i · τ i |ψ triplet = 1 4 K. Therefore, an on-site singlet-triplet splitting with magnitude equal to K takes place. This system provides a good platform, by introducing an external magnetic field h, for field-induced BEC phases.
After gap closing due to the magnetic field, the mechanism can be thought of as a bosonic system, which makes it intriguing to consider the possibility of hosting a state with coexisting diagonal and off-diagonal orders, the so-called spin supersolid state [28]. For the usual hard-core Bose-Hubbardlike Hamiltonian on a square lattice, one needs more than the nearest-neighbor interaction to stabilize the spin supersolid phase [29,30]. Various spin supersolid phases can be detected by introducing frustration [31,32] or dipole-dipole interaction [33][34][35][36][37][38][39]. Related researches for SrCu 2 (BO 3 ) 2 which can be addressed by the Shastry-Sutherland model [12][13][14]16] or spin dimers [40,41] also indicated the formation of spin supersolid. However, despite the recently proposed spin supersolid phase induced by the spin-orbital coupling [42][43][44], due to the difficulty of experimental realization such spin supersolid states cannot be easily observed. On the other hand, field-induced condensation might provide a better platform for such exotic phases. It has been shown previously that with such a two-spin Hamiltonian, by breaking the SU(2) symmetry with an anisotropic strength ∆ that controls the coupling of σ z i σ z j , formation of a spin supersolid can be triggered by external field [45][46][47]. In this work, we study the Kondo necklace model H KN using a 2D tensor network ansatz called infinite projected entangled-pair states (iPEPS) [48]. As one will see, our results reveal that not only the BEC but also spin supersolid can be triggered by the magnetic field, suggesting a good potential of BNOAS for studying exotic phases of matter.

A. Ground States in Zero Field
For discussing the phase diagram, we first generalize to the 2D Kondo-necklace XXZ model in the following way: where When ∆=1 in Eq. (3), Eq. (2) reverts to the ordinary Kondonecklace model. We utilize iPEPS as the variational ansatz and optimize it for obtaining the ground state using automatic differentiation [49,50]. Properties in the thermodynamic limit can be attained by exploiting the corner transfer matrix renormalization group (CTMRG) [51]. In this work, we choose the virtual bond dimension D = 4 and the dimension for environment tensors χ = 20, which are found to be sufficient to obtain the qualitative phase diagram. For demonstration, we perform the bond dimension scaling analysis to check the stability of each phase (see Supplementary Note 1 and Supplementary Figure 6). Further information on the iPEPS method is provided in the Method section. Note that the interlayer coupling between τ field in Eq. (1) is not included in Eq. (2). In Fig. 1(a), we demonstrate the zero-field phase diagram of 2D Kondo-necklace model (∆=1) with a pie chart. Starting from θ=0, we obtain five different phases. In the antiferromagnetic (AFM) phase, the on-site σ and τ moments are antiparallel, while nearest-neighbor σ (τ ) moments also stay in antiparallel. With increasing θ, ordering of both spin moments disappear and the system enters the spin-singlet phase through a continuous transition, with transition point at θ = 0.31π (recall that this is for D = 4).
After leaving the first quadrant, the first term in Eq. (2) favors ferromagnetic (FM) due to the sign change from plus to minus. With θ larger than 0.725π magnetic order continuously appears again but this time nearest-neighbor moments align in parallel, while σ and τ remain antiparallel on each site. Thus we name this phase as σ-FM.
Entering the third quadrant of the phase diagram, both nearest-neighbor and local coupling terms become FM, leading to the FM phase where all spins align in the same direction. Finally, in fourth quadrant the local coupling term is FM while nearest-neighbor coupling term becomes AFM. As a result, nearest-neighbor spin moments align antiparallel, but on-site moments are in parallel, becoming the σ-AFM phase. The phase transitions from σ-FM to FM, from FM to σ-AFM, and from σ-AFM to AFM are all of the first order since they break different translational symmetries.
These phase transitions can be described by measuring the following order parameters: Previous studies using an effective analytical approach (EAA) [52], quantum Monte Carlo (QMC) [53], and stochastic series expansion (SSE) [54,55] obtained the transition point from AFM to spin-singlet at K c /J = 1.37 ≈ tan(0.2993π), K c /J = 1.3888(1) ≈ tan(0.3013π), and K c /J = 1.4 ≈ tan(0.3026π) (at low T =0.05J), and they are indicated in Fig. 1(b). Note that finite D iPEPS tends to overemphasize the order parameter (see for example Figure 8 of the work by Hasik et al. [56]) and thus, with higher D we expect that the transition point predicted by iPEPS would get closer to the rest three.

B. Field-induced BEC
The spin-singlet ground state near θ= π 2 provides a good platform for the field-induced BEC after turning on the magnetic field. To see this, note that the elementary excitations over the singlet ground state are the mobile triplons carrying the quantum numbers S z =0, ±1. With increasing field strength, the triplon band with S z =+1 becomes more favorable and finally crosses the energy of the spin-singlet ground state, leading to the condensation of the S z =+1 triplon. The local densities of triplons can be measured by the following operators: where ρ ± i and ρ 0 i are densities for the local triplon with S z i =±1 and S z i =0 quantum numbers, respectively. The local density of the singlet is determined by Utilizing the triplon density operator ρ + i , we can define an order parameter as follows: where |ρ + | > 0 reflects the translational symmetry breaking. For a complete characterization of the quantum phases, we define the condensate density order parameter that detects the U(1) symmetry breaking: FIG. 2. Field-induced Bose-Einstein condensation for HKN with θ=0.4π. h stands for the strength of field while SSF and FM represent the spin superfluid and ferromagnetic phases distinctly. |ρ + | and n0 come from the definition in Eqs. (7) and (8). The small inset demonstrates n0 for bond dimension D = 4 and D = 6 (dimension for environment tensors χ = 36) separately near the transition point from spin singlet to SSF. We can see that the transition point is located at around h = 0.58 from both bond dimensions.
One can also see that its finite value reflects the mixture of the spin singlet and the triplon with that occurs in the field-induced BEC. Recall that the magnetic field favors the S z =+1 triplon band, causing the gap closing and eventually crossing with the singlet band. This fulfills the general understanding of field-induced BEC that after band crossing, S z =+1 band and singlet band are hybridized, resulting in BEC. Fig. 2 displays two order parameters, i.e.,ρ + and n 0 , as a function of the external field strength at θ=0.4π (K ≈ 3J).
X X   TABLE I. The symmetry table for each phase. z-antiferromagnetic, spin supersolid, spin superfluid, and ferromagnetic phases are abbreviated into z-AFM, SS, SSF, and FM separately. indicates the symmetry is present and X indicates it is broken. Orders that appear after translation and U(1) symmetry breaking are also indicated in the brackets and their definitions come from Eqs. (7) and (8). Notice that the translational symmetry breaking denotes the appearance of spatial inhomogeneity in the triplon density, while the U(1) symmetry breaking indicates the condensation of the S z =+1 triplon. z-AFM and Solid both possess U(1) symmetry but have different quantum numbers S z : zero for z-AFM while 0.5 for Solid. It is clear to see from the symmetry that a first-order transition takes place at the phase boundary between Solid and SSF, where phases break distinct symmetries on different sides.
The spin-singlet phase remains stable before the ground state enters into the SSF phase around h ≈ 0.58 at which the S z =+1 triplon starts condensing. If we further enhance the magnetic field, finally spin moments are fully polarized and the asymptotic FM plateau appears, restoring U(1) symmetry. Note that the triplon density remains uniform (ρ = 0), and both transitions are expected to be continuous due to the U(1) symmetry breaking.

C. Combination of Field and Anisotropy
In the previous sections we have demonstrated the phase diagram for the 2D Kondo-necklace model in zero field and the field-induced BEC out of the spin-singlet phase due to the magnetic field. We expect a richer phase structure may emerge by introducing the XXZ anisotropy in the Heisenberg interaction. In the boson language, σ z i σ z j is mapped to the repulsive interaction between neighboring bosons. While the interaction prefers the density wave or low density of triplons, the external field stabilizes dense populations of the S z =+1 triplon. Indeed, we find that these competing effects give rise to various quantum phases, and the phase diagram is presented in Fig. 3(a). To characterize the phases, we also show the order parameters in Fig. 3(b) as a function of h at ∆=3.
In the absence of the field, the strong anisotropy results in a trivial magnetic state, named after z-antiferromagnetic (z-AFM) state, out of the spin-singlet phase where σ spins form the Néel configuration inẑ direction, and τ spins align antiparallel to the on-site σ spin as well. The spin-singlet and z-AFM phases share the same U(1) quantum number while the uniformity of the triplon density breaks seemingly continuously, suggesting the continuous transition (see Supplementary Note 1 and Supplementary Figure 5(b)). On the other hand, as increasing ∆ or the repulsive interaction between triplons, the system evolves into a Solid phase out of the SSF phase in a wide region of the phase diagram. Note that the Solid phase is characterized by a checkerboard pattern of the S z =+1 triplon density with a fractional number per unit cell, i.e., ρ + =1/2 in Fig. 3(b), which resembles the Mott phase for a fermionic system with large Hubbard U interaction. Otherwise this is a factional magnetization plateau phase in the spin language. If we further enhance the strength of external field, the Solid state will melt down and SSF appears again. The first-order transition between Solid and SSF phases can be likely described by the XXZ model on the square lattice, which serves as the leading order mapping from the original Hamiltonian [57]. It is worth noting that the S z =−1 triplon is completely suppressed, i.e., ρ − =0 throughout the phase (See Supplementary Note 1 and Supplementary Figure 5(a)). We also find that a spin supersolid phase, which breaks both the uniformity of the triplon density and U(1) symmetry simultaneously, appears in a narrow region indicated in Fig. 3(a) and (b). The strong repulsive interaction introduced by XXZ anisotropy stabilizes the density wave of condensed triplons and thus gives rise to the translational symmetry breaking out of the SSF. The density wave is characterized by a (π, π) wave vector. In order to visualize each phase, we illustrate the order parameters schematically in Fig. 3(c), and summarize the symmetry properties of phases in Table I.

III. RELATION TO BNOAS
Our work was stimulated by the report of the BNOAS insulator built on a d 8 Ni ion. It was noted in the Introduction that the Kondo sieve model is the minimal spin Hamiltonian for BNOAS, which is three-dimensional (3D) but contains a layered structure. While one needs interlayer interaction for the 3D coupling, our 2D Kondo necklace model is expected to establish the underlying intralayer behavior. In addition, this is a good platform for study of the field-induced magnetic orders, besides the well-known examples TlCuCl 3 or SrCu 2 (BO 3 ) 2 . More interestingly, by introducing the symmetry required XXZ anisotropy the field-induced spin supersolid can be realized. Note that realization of supersolids has been an active research topic and recently reported to arise from a BEC made of dipolar atoms [58][59][60][61][62][63]. However, such realization with the cold-atom equipment is not an easy task and therefore we do not have many examples so far, considering the date when this concept of supersolid was first proposed [28].
As suggested by Ng and Lee [64], instead of searching for a supersolid phase at very low temperature, magnetic materials with spin singlets in their ground states provide a more promising scenario for its formation. For BNOAS, the "local" spin and "Kondo" spin moments both arise from electrons in e g subshell of Ni. Earlier work had revealed that the Ni ion in 2D materials is suitable for the study of XXZ-type antiferromagnetism [65][66][67]. These discoveries indicate that BNOAS has the potential to serve as a good platform for the realization of a spin supersolid phase, making it worthwhile to study BNOAS further, theoretically and experimentally, and exploit more of its underlying physics.
As mentioned by Jin et al. [24], doping electrons into BNOAS leads towards the region of possible high-T c superconductivity, considering the similarities to hole-doped  2) and (3). h stands for the field strength and ∆ represents the strength of XXZ anisotropy. z-antiferromagnetic, spin supersolid, spin superfluid, and ferromagnetic phases are abbreviated into z-AFM, SS, SSF, and FM separately. The SS phase lies in the range 0.57 h 0.7 with ∆ 2. Filled and empty circles indicate the first-order and second-order phase boundaries, separately. A more careful investigation for those continuous boundaries is non-trivial and thus we will leave it for the future works. (b) Order parameters along with h for ∆=3. Their definitions can be found in Eqs. (5), (7), and (8). The inset shows the results obtained using bond dimension D=5 (dimension for environment tensors χ=25), focusing on the SS area. (c) The triplon configurations for each phase. Within the 2 × 2 unit cell, the size of each lattice site stands for the magnitude of ρ + i . Blue (white) color reflects nonzero (zero) condensation (n0). Here, with blue color it indicates U(1) symmetry breaking. Note that for z-AFM, ρ + i has a checkerboard pattern with a small but nonzero on-site magnitude, and thus the homogeneity is broken and we have small but nonzeroρ + order.
cuprates and Sr-doped NdNiO 2 superconductors [68]. In addition, the spin-singlet state serves as another scenarioan unusual one -for a self-doped Mott insulator [69]. For Nd 1−x Sr x NiO 2 , the sparse Nd 5d conduction carriers may couple with Ni 3d x 2 −y 2 electrons to form Kondo singlets dynamically. Such singlets will suppress the AFM order, leading to a paramagnetic ground state which can be metallic. For BNOAS, on the other hand, without doping we have a magnetically inert ground state with a Kondo singlet occupying every site. Because of the hard-core nature of the on-site singlets, it should be an insulator. Upon doping BNOAS, however, the singlet density decreases while the long-range AFM order is hindered unless a high doping level. Therefore, with an intermediate doping level we suggest that a insulator-metal transition could also take place. Extending these similarities, this self-doped superconducting transition can be modeled by a t−J-like Hamiltonian [69]. Here we propose an effective Hamiltonian for describing the microscopic mechanism for electron-doped BNOAS. We start from the 3d x 2 −y 2 spins (σ) and adopt an anisotropic t − J model where c iα (c † iα ) is the annihilation (creation) operator for 3d x 2 −y 2 electrons with α =↑, ↓ and H.C. denotes the Hermitian conjugate. The spin operator is connected to fermionic operator by σ β i = 1 2 α,α c † iα ρ β α,α c iα where ρ β is the Pauli matrix with β=x, y, z. P G =Π i (1 − n i↑ n i↓ ), with n iα =c † iα c iα , is the Gutzwiller projection operator to prevent the double occupancy on each site [70]. Longer-range hopping can be also included to better explain experimental observations but here we only demonstrate the nearest-neighbor hopping. While the t − J model is often applied for hole-doped cuprate superconductors, its electron-doped counterpart only requires a sign change for the hopping constant due to the particle-hole transformation and thus we can directly borrow the same form here [71]. Additionally, to preserve the degree of freedom for the XXZ anisotropy, we retain the form of ( · ) ∆ for the superexchange term.
The 3d z 2 spin τ is coupled to the d x 2 −y 2 spin σ and can be described by H K : The final effective Hamiltonian is the sum H eff = H J + H K . This Hamiltonian, which we call the t − J − K model, can be numerically solved with various methodologies [72][73][74][75][76][77][78] by Monte Carlo [79] or renormalized mean-field theory [80] besides tensor network, for its properties in real and momentum space. Thus, we consider this to be our future challenge.

IV. SUMMARY
In this work, we make use of iPEPS, a 2D tensor network ansatz, to solve the Kondo necklace model in two dimensions. Without the XXZ anisotropy, we obtain the zero-field phase diagram and locate the region of spin-singlet formation. Upon turning on an external magnetic field, the spinsinglet phase goes through a phase transition into spin superfluid, a well-known phenomenon called field-induced BEC. By adding XXZ anisotropy, we argue that now the triplet state with S z =0 is more favorable and closes the gap with large enough anisotropy ∆. With external field, an exotic spin supersolid phase appears between two magnetization plateaus. We provide a ∆-h phase diagram and relocate the region where we believe spin supersolid can be realized with 2D Kondo-necklace model.
Since BNOAS has (weakly) coupled infinite-layer nickelate planes, and adopting the 3D Kondo sieve model to be its effective Hamiltonian, we expect such field-induced BEC/spin supersolid phases can be realized within its nickel-oxide layer. For a more careful investigation in the future we will consider probing the continuous phase boundaries in more detail for their universality class or critical exponents. Moreover, QMC provides another useful tool, especially for finite temperature. Such thermalization can also be approached by iPEPS through a purification process and it shows good consistency with QMC [81]. It will be of great interest to study the thermal properties of Kondo-necklace model with both numerical techniques and thus one of our future works.
Looking toward extending theoretical work, we propose an effective t − J − K Hamiltonian that can be used to describe the potential superconductivity arising in BNOAS after doping. Further studies for this effect from both experimental and theoretical sides are again of great interest and regarded as our future goal, too.

A. Infinite Projected Entangled-pair State
In this work we adopt the 2D tensor network ansatz, the infinite projected entangled-pair states, for our calculation. The iPEPS tensor network ansatz is an effective numerical method for dealing with various quantum problems in two dimensions [48,74,75,82]. It has many merits that some other numerical techniques do not have, such as its attainability to the thermodynamic limit and freedom from the vicious sign FIG. 4. Some basics for the infinite projected entangled-pair state ansatz. (a) The four bulk tensors a1, a2, a3, and a4, with bond dimension D for each bond. After tracing out the physical bonds with their complex-conjugate tensors (a † 1 to a † 4 ), we obtain the doublelayered tensors A1, A2, A3, and A4 shown in (b). Along with eight edge tensors (T1 to T8) and four corner tensors (K1 to K4), the norm of ansatz, Ψ|Ψ , can be calculated as the tensor graph in (b). The black thick bond is of dimension D 2 while dotted bonds have dimension χ. Area enclosed by red dashed line is the unit cell. In (c) we demonstrate the computation graph for the GS energy EGS from the initial input to be the bulk tensors. CTMRG stands for the corner transfer matrix renormalization group and RDM refers to the reduced density matrix. H is our target Hamiltonian. problem in quantum Monte Carlo. This ansatz contains two parts, the bulk tensors and environmental tensors for achieving the infinite size. Technically, the number of bulk tensors can be freely chosen, but a better choice will allow it to be able to reflect the real space modulation for the ground state (GS). Thus, as shown in Fig. 4(a), in this work we apply a 2 × 2 unit cell for the bulk (we have also tried other bulk sizes, see Supplementary Note 2), meaning that there are four different rank-5 bulk tensors from a 1 to a 4 . For each tensor it has four virtual bonds (bonds that connect nearby tensors), capturing the entanglement between site to site, with bond dimension D. The remaining one leg of each tensor, connecting a n and a † n , is the physical bond with dimension d, equal to the dimension of local Hilbert space. The accuracy of iPEPS method is determined by the bond dimension since larger D captures better the entanglement among sites. For critical systems we usually need a larger D since the correlation length diverges. However, for quantum states that are less entangled, such as the scenario in this work, we have finite correlation length and thus the results remain nearly the same after an adequate D.
With bulk tensors, the iPEPS ansatz is complete with the environmental tensors in Fig. 4(b), where rank-2 K n tensors stand for the corners and rank-3 T n tensors for the edges (n = 1 to 4). A 1 to A 4 tensors correspond to double-layered tensors by tracing out the physical bonds, as indicated in Fig. 4(a). The dimension for K tensors is χ · χ while for T tensors it is D 2 · χ · χ. The environmental tensors can be obtained from double-layered tensors through a process called CTMRG [51,75,83]. With corner and edge tensors, not only can we extrapolate the system size to the infinity, but the physical observables can also be calculated by constructing the corresponding reduced density matrix.

B. Tensor optimization
With the structure of iPEPS explained above, next we discuss how to optimize this ansatz. Since the environmental tensors come from bulk tensors, the number of variational parameters is decided by the bulk. Thus, optimization of iPEPS concerns the optimization of the bulk tensor. Traditionally, people make use of the technique called simple or full update based on the imaginary-time evolution for the optimization [84]. On the other hand, as a variational ansatz, a direct minimization of GS energy by changing variational parameters systematically based on the gradients might be a more direct approach. Nonetheless, it is not an easy job to evaluate the energy gradient of each variational parameter [85].
Recently, a new way of calculating such energy gradients has been proposed by using the technique called automatic differentiation [49]. It is a numerical way of evaluating function derivatives to machine precision and is often applied for updating the neural network for machine learning [86,87]. The idea of automatic differentiation is based on the chain rule: it assumes that a numerical function is composed of elementary operations (addition, subtraction, multiplication, and division) and functions (sin, log, etc.), irrespective of the complexity. For a given function, we can visualize its computation process by constructing a computation graph composed of individual nodes, where each node represents a intermediate result. The inputs will go through the computation graph and produce the final output. Later, by applying the forward or backward propagation, we are able to obtain the derivatives of the outputs with respect to each input.

C. Working Flow
Since our problem here has multiple inputs (tensor elements as the variational parameters) and only one output (GS energy, E GS ), it is more natural to apply the backward mode. First we need to record down the computation graph from initial bulk tensors to the final energy, as demonstrated in Fig. 4(c). Starting from bulk tensors, we can construct double-layered tensors and initial environmental tensors (K 0 and T 0 ). By applying CTMRG until convergence, we obtain the effective final environment with K N and T N . Then we construct the reduced density matrix with bulk tensors and environment. Finally, along with the Hamiltonian, we can compute E GS . After completing one computation flow (one epoch), the energy gradients are evaluated by backward propagation, and we make use of these gradients to update the tensor elements with a desired degree (learning rate). After a large enough number of epochs, the energy converges and we obtain a good GS ansatz for further calculation of physical observables, or for other investigations of the mother Hamiltonian. Please refer to the open repository by Hasik et al. [50] for a practical package of this iPEPS method with automatic differentiation.

VII. DATA AVAILABILITY
The authors declare that the main data supporting the findings of this study are available within the article and its Supplementary Discussion. All relevant data in this paper are available from the authors upon reasonable request. An open access repository for the basic codes of iPEPS and datasets is available at https://doi.org/10.5281/zenodo.6420017. Here, we present the triplon densities as a function of the field strength h at ∆=3.0. The triplon densities can be measured as follows: Note that the densities do not change in the z-AFM and Solid phases while increasing the strength of the magnetic field. That is because the U(1) symmetry is preserved, and thus the averaged density of each triplon remains constant. It is worth noting that the density of S z =−1 tripon, ρ − , is completely suppressed in the Solid phase while the one of S z =+1 triplon is 1/2, indicating that the Solid phase is a magnetization plateaux phase with a fractional magnetization m z =ρ + − ρ − =1/2. The results are demonstrated in Supplementary Fig. 5(a). For the transition of spin singlet to z-AFM, we demonstrate the |ρ + | to ∆ plot in Supplementary Fig. 5(b). A continuous phase transition is indicated with our numerical results. However, due to the difficulty of examining the phase boundary with finite D, we leave the further investigation for the nature of this phase transition for the future works.
In Supplementary Fig. 6(a), we extrapolate the two order parameters, |ρ + | and n 0 at h = 0.6, to infinite D. We obtain the results that |ρ + | D→∞ = 0.2202 and n 0D→∞ = 0.1236, suggesting that these two orders will remain coexisting for spin supersolid phase. We also make the similar extrapolation of ρ + for Solid and SSF phases in Supplementary Fig. 6(b). The result shows consistency after increasing D and for SSF (Solid) phase ρ + D→∞ =0.1212 (0.5).   Table II. We can easily find out that for the 2×2 and 4×2 unit cells they share very close energies and are always favored. It is because they are both commensurate to the minimum unit cell required for our antiferromagnetic Hamiltonian. On the other hand, for 3×2 or 3×3 unit cell not only the energies are higher but also the accuracy is poor for those translation-symmetry-breaking phases. This effect is because of the incommensurability, which causes a very slow or even unconvergent CTMRG process. As a result, we argue that for our 2D Kondo-necklace model in the square lattice, due to the lack of frustration, a 2×2 unit cell is the most suitable and efficient for studying its behaviour.