Investigation of the N\'eel phase of the frustrated Heisenberg antiferromagnet by differentiable symmetric tensor networks

The recent progress in the optimization of two-dimensional tensor networks [H.-J. Liao, J.-G. Liu, L. Wang, and T. Xiang, Phys. Rev. X ${\bf 9}$, 031041 (2019)] based on automatic differentiation opened the way towards precise and fast optimization of such states and, in particular, infinite projected entangled-pair states (iPEPS) that constitute a generic-purpose ${\it Ansatz}$ for lattice problems governed by local Hamiltonians. In this work, we perform an extensive study of a paradigmatic model of frustrated magnetism, the $J_1-J_2$ Heisenberg antiferromagnet on the square lattice. By using advances in both optimization and subsequent data analysis, through finite correlation-length scaling, we report accurate estimations of the magnetization curve in the N\'eel phase for $J_2/J_1 \le 0.45$. The unrestricted iPEPS simulations reveal an $U(1)$ symmetric structure, which we identify and impose on tensors, resulting in a clean and consistent picture of antiferromagnetic order vanishing at the phase transition with a quantum paramagnet at $J_2/J_1 \approx 0.46(1)$. The present methodology can be extended beyond this model to study generic order-to-disorder transitions in magnetic systems.


I. INTRODUCTION
The spin-S antiferromagnet, with isotropic coupling J 1 between nearest-neighbor spins located on the sites of a square lattice, represents one of the most paradigmatic models of quantum magnetism. At zero temperature, the system develops long-range antiferromagnetic (Néel) order for any value of S: while for S ≥ 1 there are analytical arguments [1,2], for the extreme quantum case with S = 1/2, this has been numerically proven thanks to quantum Monte Carlo simulations on large systems [3][4][5]. Instead, any finite temperature will restore spin rotation symmetry, in agreement with the Mermin-Wagner theorem [6]. A magnetically disordered ground state may be also achieved by including further super-exchange couplings, most notably a next-nearest-neighbor interaction J 2 , which destabilizes the Néel order driving towards a quantum phase transition. In this respect, much effort has been spent to understand the ground-state properties of the J 1 − J 2 model defined by: Here, . . . and . . . stand for nearest-neighbor and next-nearest-neighbor sites on the square lattice, respectively; S i = (S x i , S y i , S z i ) represents the spin-1/2 operator on the site i. Both the spin-spin interactions are taken positive.
In the presence of finite J 2 a severe sign problem is present (especially in the local basis with z-component * Juraj.hasik@irsamc.ups-tlse.fr defined on each site), which prohibits quantum Monte Carlo algorithms from assessing large system sizes [7][8][9]. Over the last three decades several alternative methods have been introduced and kept improving such as exact diagonalizations, density-matrix renormalization group (DMRG), functional-renormalization group (fRG), and variational Monte Carlo (VMC) approaches. The ground-state properties of the J 1 − J 2 model have been intensively investigated, with contradicting results, supporting the existence of a valence-bond solid (with either columnar or plaquette order) [10][11][12] or a spin liquid (either gapped or gapless) [13][14][15][16], or even both [17][18][19][20][21]. One important aspect emerging in the latest calculations is the existence of a continuous quantum phase transition between the antiferromagnetic and the paramagnetic phases for J 2 /J 1 ≈ 0.5, where the staggered magnetization (hereafter named simply "magnetization") goes to zero.
Recently, borrowing concepts from quantum information, tensor-network methods have been introduced [22][23][24]. In one dimension, the so-called matrix-product states (MPS) offer a convenient and elegant rephrasing of previous DMRG ideas. MPS evolved into the method of choice and provide very accurate approximations of the exact ground-state properties. Generalizations in two dimensions are more problematic. The prominent example, projected entangled-pair states (PEPS), provide the correct entanglement structure of most quantum ground states of local spin Hamiltonians [25], however, they suffer from a steep scaling of computational effort when enlarging the system size. For this reason, their application has been limited to ladder geometries with small number of legs [26,27] and finite 2D clusters with open boundary and up to ≈ 200 − 300 sites [28]. In order to overcome this computational barrier and avoid boundary effects, algorithms that work directly in the thermodynamic limit (dubbed iPEPS) have been introduced and developed [29,30]: here, only a small number of tensors is explicitly considered and embedded into an environment that is self-consistently obtained (e.g., within the so-called corner-transfer matrix approaches [31] or channels [32]). The size of these tensors, and in turn the number of variational parameters of the wave function, is characterized by the so-called bond dimension D. The iPEPS are systematically improved by enlarging the bond dimension, accounting for increasingly entangled states.
In recent years, iPEPS have been applied to assess the nature of the ground state of the J 1 − J 2 model, mainly focusing on the highly-frustrated regime J 2 /J 1 ≈ 0.5 [33][34][35]. However, these attempts were not completely satisfactory, since they either used a simplified tensor structure, limited to the description of paramagnets, or suffer from optimization problems, arising in methods that are not fully satisfactory and consistent (e.g., the so-called simple and full update [29,36]). In this respect, a breakthrough in the field has been achieved by performing the tensor optimization using the ideas of algorithmic differentiation, or better the adjoint algorithmic differentiation (AAD) technique, which allow a very efficient optimization even in the presence of large number of parameters [37]. Here, Liao and collaborators limited their application to the unfrustrated Heisenberg model (with J 2 = 0), showing that extremely accurate and completely stable results may be obtained for both the ground-state energy and magnetization.
Even though PEPS (and iPEPS) Ansätze are designed to describe both gapped and gapless states (following the entanglement entropy's area law, up to additive corrections), it remains an open question whether generic optimization can reliably reproduce highly-entangled ground states, as the ones that are possibly emerging in the frustrated regime J 2 /J 1 ≈ 0.5 [15,18,20,21]. Therefore, in this work, we do not directly address the question of the nature of the magnetically disordered phase; instead, we focus our attention to the magnetically ordered phase with J 2 /J 1 ≤ 0.45 and perform an accurate determination of the magnetization curve as a function of the frustrating ratio. In addition to its conceptual importance, the problem of the disappearance of antiferromagnetic order under increasing frustration offers a stringent test to most numerical methods, in general, and to tensor network methods, in particular. To this end, we apply the same ideas of AAD to optimize the iPEPS Ansatz for the J 1 − J 2 model of Eq. (1). Importantly, unlike the previously proposed gradient-based optimizations [32,38], the AAD can be effortlessly extended beyond nearestneighbour Hamiltonians. The energy and magnetization are obtained for different values of the bond dimension D, from 2 up to 7. Then, the estimates for D → ∞ are obtained for each frustration ratio J 2 /J 1 . Note however that, this is not realized by a crude extrapolation in 1/D (for which the results for different values of D are considerably scattered) but, instead, by performing a correlation-length extrapolation, which is motivated by the finite-size scaling analysis that is well established in the Néel phase, as recently proposed in Refs. [39,40]. Despite the fact that this mode of extrapolation requires the calculation of the correlation length ξ, which may not be as accurate as other thermodynamic quantities (e.g., energy and magnetization), it has been shown to give remarkably good results for the unfrustrated Heisenberg model. In fact, as we have mentioned earlier, even though iPEPS can describe certain gapless phases, their generic optimization instead leads to states with finite correlation lengths. e.g., in the Néel phase, and the bond dimension D turn out not to be the correct object to quantify this aspect. As we will show, also in the presence of frustration, the analysis based on the correlation length gives reliable thermodynamic estimates, even though no exact results are available. Our calculations are compatible with a vanishing magnetization for J 2 /J 1 ≈ 0.45, which is in close agreement with recent calculations [15,[18][19][20][21] and give a reference for future investigations.
The paper is organized as follows: in section II, we will describe the iPEPS method; in section III, we present the results; in section IV, we finally draw our conclusions and discuss the perspectives.
The tensor a is chosen (and constructed such as) to be invariant under a number of symmetries. First of all, it belongs to the A 1 irreducible representation of the C 4v point group, thus enforcing all the spatial symmetries of the square lattice on the iPEPS. The antiferromagnetic correlations are incorporated in the ansatz by unitaries −iσ y , which rotate the physical S z basis at every site of one sublattice. We absorb these unitaries into observables leaving the definition of the wave function untouched [see Eq. (12)].
Secondly, the tensor a also possesses a further structure by requiring certain transformation properties under the action of U (1) group (see below). Such choice is motivated by the remaining U (1) symmetry in the ordered phase, which manifests itself as equivalence between different magnetizations connected by transverse (Goldstone) modes. As defined below, U (1) tensor classes are defined by assigning specific "charges" to the virtual and physical degrees of freedom.
When considering A 1 -and U (1)-symmetric states, the tensor a = a( λ) is taken to be a linear combination of (fixed) elementary tensors {t 0 , t 1 , . . .} (named a tensor "class") such that with coefficients λ being the variational parameters. The elementary tensors {t 0 , t 1 , . . .} are different representatives of the A 1 irreducible representation for some choice of the U (1) charges. Given an iPEPS defined by tensor a, the evaluation of any observable O amounts to a contraction of infinite double-layer network composed of tensors a together with the tensor representation of O. Such tensor network is the diagrammatic equivalent of usual expression O = ψ(a)|O|ψ(a) . A central aspect of iPEPS method is an approximate contraction of such networks. In this work, we realize them by finding the so-called environment tensors C and T of dimension χ, dubbed environment dimension, by the means of corner-transfer matrix (CTM) procedure [31]. These tensors compress the parts of the original infinite network in approximate but finitedimensional objects. Afterwards, the desired reduced density matrices can be constructed from C and T , together with the on-site tensor a. Ultimately, the exact value of any observable is recovered taking χ → ∞, which we extrapolate from the data for increasingly large χ.
The optimization of tensor a (or equivalently in the U (1)-symmetric approach, λ) is carried out using standard gradient-based method L-BFGS supplemented with backtracking linesearch. The gradients are evaluated by AAD, which back-propagates the gradient through the whole process of energy evaluation for fixed χ [37]: Starting with a given CTM, followed by assembling the reduced-density matrices from converged C, T tensors and finally evaluating the spin-spin interaction between nearest and next-nearest neighbors.

B. Extracting the relevant U (1) charges
For small enough frustration, in the Néel phase, the unconstrained optimization of tensor a leading to correct U (1)-symmetric iPEPS would be a desirable outcome. Key steps of the CTM algorithm for single-site iPEPS with C4v symmetry. (i) Initial tensors at iteration i: (ii) Construction of enlarged corner and its reshaping into matrix of dimensions D 2 χ×D 2 χ. (iii) Symmetric eigenvalue decomposition of enlarged corner and truncation down to leading χ eigenpairs by magnitude of the eigenvalues. Truncation is always done at the boundaries between degenerate eigenvalues (see text). (iv) Absorption and truncation with isometry P from step (ii) for half-row/-column tensor T .
Under circumstances, AAD optimization can arrive at an almost U (1)-symmetric tensorã. In such case, a direct and robust evidence can be seen in the nearly degenerate pairs of leading eigenvalues of the transfer matrix. Importantly, such iPEPS states provide an unbiased information about the energetically favourable U (1)-charge structure of tensor a. We are concerned with inferring these charges from the elements of tensorã. Obtaining the correct charge assignment of the smallest D tensors allows (i) to perform an efficient variational optimization over a greatly reduced number of parameters λ, (ii) to obtain truly U (1)-symmetric environments via CTM and, finally, (iii) to predict the correct charge content of higher-D a tensors and, hence, enable to perform (i) and (ii) for larger D.
Before describing how to achieve the goal of obtaining the charges from the almost symmetricã tensor, let us first briefly review the expected properties of the resulting U (1)-symmetric a tensor. In practice, one has to assign U (1) charges u = (u ↑ , u ↓ ) and v = (v 0 , . . . , v D−1 ) to the two physical spin-1/2 components and the D virtual degrees of freedom on each of the four auxiliary indices. Without loss of generality we take them to be integers. In this language, U (1) invariance is realized by simply enforcing a selection rule for the non-zero tensor elements a s uldr which should exhibit a local charge conservation where N is some fixed integer. Notice that in order to preserve C 4v symmetry the same v, associated to all the virtual indices, is taken on the four legs of the tensor. Note also that there is some freedom in the definition of the charges since shifts like u s → u s + α, v σ → v σ + β, and N → N + α + 4β, with α and β ∈ Z, leaves Eq. (5) invariant. It is easy to connect this charge conservation to the U (1) invariance of the a tensor. Indeed, the action of any element g ∈ U (1) on a is given by: where U and V are diagonal matrices depending on g, and all auxiliary indices are transformed by the same V : with the phase θ g ∈ R and γ = 0, . . . , D − 1. Therefore, the non-zero elements of the tensor a transform according to (ga) s uldr = a s uldr e iθg(u s +vu+v l +v d +vr) Hence, Eq. (5) implies that a is indeed invariant up to global phase under the action of U (1). Once the relevant U (1) charges u and v are known (see below), practically, Eq. (5) is used in the construction of the elementary tensors {t 0 , t 1 , . . . } by filtering out their non-zero elements.
Let us now describe how to infer the charges from an unrestricted tensor optimization that has produced an almost symmetric on-site tensorã, with bond dimension D. To identify the dominant (at least for small D) U (1)-symmetric component ofã, and then ultimately derive the hidden U (1) charges, we have to first perform a higher-order singular value decomposition ofã: with unitary matrices Z, Y , and the so-called core tensor c. The same unitary Y is associated to different auxiliary legs due to the enforced C 4v symmetry. The core tensor c plays an analogous role to singular values in standard singular value decomposition of a matrix. The untruncated core tensor c by itself defines a physically equivalent iPEPS to the one given byã. A good lower-rank approximation ofã can be obtained by truncation of the smallest elements of the core tensor c. The basic premise, supported by nearly degenerate transfer matrix spectrum for small D, is that the relative magnitude of symmetry breaking elements of tensor c is small. Therefore, we assume that the largest elements of tensor c respect the U (1)-symmetry constrain associated to an unknown set of charges u and v.
For the last step in identifying the charges, we reformulate the problem in terms of linear algebra. First, taking a set of n largest tensor elements (modulo C 4v symmetry), and writing down Eq. (5) for each of them will result in a set of n coupled linear equations (with integer coefficients) of the D + 2 unknown charges. Whenever n > D + 2, the linear system becomes over-complete and, increasing n still allows the same solution for the charges, unless n is taken too large so that (small) nonzero tensor elements breaking U (1)-symmetry are included. To solve this linear problem it is convenient to recast the constraints into a n × (D + 2) matrix. The matrix, containing integer matrix elements, is obtained by simply counting the total number of charges of each type γ and s on the virtual and physical legs, respectively. More precisely, we define vectors n(c s uldr ) of integer coordinates that count the number of times specific index value appears among the indices of a given tensor element. Expressing each individual element constraint (5) as n(c s uldr ) · ( u, v) = N and recasting them into matrix form, the linear system can be written in a compact fashion as M · ( u, v) = N .
To be explicit, let us consider the case D = 3 for which all charges can be obtained using only the n = D + 2 = 5  largest tensor elements of tensor c: If the tensor c possesses an (approximate) U (1) symmetry structure (as in the example above), then the linear system has a non-trivial solution in terms of charges u and v. To solve it, it is known that one needs to bring the matrix M into its Smith normal form (see Appendix A). Note here that the integer N can, in fact, be changed arbitrarily. Although, the explicit values of the charges will depend on N , the U (1) class of a tensors will not. In other word, there is some "gauge" freedom to determine each U (1) class. For the example with D = 3 considered here, we get integer charges, u ↑ = +1, u ↓ = −1, v 0 = 0, v 1 = 2 and v 2 = 0, as can be checked by direct substitution in Eq. (11) choosing N = 1. A complete list of the relevant charges are shown in Table I for bond dimension up to D = 9.
C. Reduced density matrices, CTM algorithm, and implementation details The evaluation of energy is realized through two distinct reduced-density matrices (RDM), ρ (N N ) and ρ (N N N ) , for nearest and next-nearest neighbour sites respectively. Their diagrammatic definition is shown in Fig. 1. The energy per site is then given by: withS α = −σ y S α (σ y ) T , as these are the only nonequivalent terms of Hamiltonian (1) acting on the singlesite iPEPS with C 4v symmetry.
The two RDMs are obtained by substituting the environment of a 2 × 2 cluster within the infinite tensor network with the CTM approximation and tracing out all but two (nearest-neighbor or next-nearest-neighbor) sites. The leading computational cost in contraction of these networks is O[(χD 2 ) 3 p 2 ] with p = 2 being the dimension of the physical index s. A more complete alternative is to consider a RDM of all four spins contained within the cluster. However, contracting such network with eight open physical indices is more expensive in terms of computational complexity and memory requirements, as both are amplified by a factor of p 2 .
The most demanding part of the calculations is the CTM algorithm. Given the highly constrained nature of our iPEPS, in particular the C 4v symmetry imposed on tensor a, we can utilize the efficient formulation of the algorithm of Ref. [31]. The C 4v symmetry of the on-site tensor a is reflected in the corner matrix C which is taken to be diagonal and half-row/-column tensor T which is symmetric with respect to the permutation of its environment indices. We show the diagrammatic description of the main steps within single CTM iteration in the Fig. 2.
There are few more remarks to be made regarding the implementation of the CTM algorithm. The initial C and T tensors are given by partially contracted double-layer tensor, e.g. C (dd )(rr ) = sul a s uldr a s uld r . In addition, after each step of the CTM the tensors C and T are symmetrized accordingly and normalized by their largest element. To establish the convergence of the CTM, we use the RDM of nearest neighbors ρ (N N ) 2×1 computed just from the 2 × 1 cluster at each CTM step. Once the difference (in sense of Frobenius norm) between ρ (N N ) 2×1 from two consecutive iterations becomes smaller than CT M , we consider the CTM converged. During optimization we set CT M = 10 −8 , which typically requires at most O(70) iterations to converge for largest (D, χ opt ) = (7, 147) simulations considered. For scaling of observables of optimized states we instead iterate CTM until CT M = 10 −12 . Remarkably, the U (1) symmetry is preserved along the CTM procedure, whenever we adjust the truncation as to never break the multiplet structure of the enlarged corner.
Finally, a peculiar complication is present in the process of computing gradients by AAD, with two distinct aspects. First, the standard definition of adjoint function of eigenvalue (or singular value) decomposition relies on computing the full decomposition [41]. Hence, in this context one cannot resort to significantly faster partial decompositions such as Lanczos (at least during gradient computation). This sets the leading complexity of CTM iteration to O[(χD 2 ) 3 ]. Recently, a developed differentiable dominant eigensolver tries to address this shortcoming by alternative adjoint formula [42]. The second, more fundamental aspect is the ill-defined adjoint in the case of degenerate eigenvalues stemming from the terms proportional to the inverse of spectral gaps. We use a smooth cutoff function [37] to tame this problematic terms. Although doing so, the accidental crossings of eigenvalues in course of CTM sometimes result in erroneous gradients. In general, we found this occurrence, manifested by the failure of linesearch, to be rare. The formulation of AAD applied to gauge-invariant scalars (such as energy), whose computation however involves eigendecomposition with degenerate spectrum, still remains an open problem.
The complete algorithm is available as a part of the open-source library peps-torch [43] focused on AAD optimization of iPEPS.

III. RESULTS
Our analysis is based upon an extensive set of calculations for various bond dimensions, ranging from D = 2 to 7, and different values of the frustrating ratio J 2 /J 1 up to 0.5. For the large bond dimensions considered, the optimizations have been performed with environment di- Due to C4v symmetry imposed on the ansatz, the transfer matrix is symmetric and can be diagonalized. Eigenvalues are ordered with descending magnitude with the leading eigenvalue λ0 normalized to unity. Bottom: RDM for two-point correlation functions, defined for r ≥ 1, and its connection to transfer matrix E. mensions up to χ opt = 4D 2 in the case of D = 5, 6 and up to χ opt = 3D 2 for D = 7. Here, we want to highlight a few important aspects of iPEPS that are crucial for the investigation of the magnetically ordered phase. First of all, within optimizations with no imposed symmetries, there is a generic tendency to break the physical U (1) symmetry of the Néel state (corresponding to global rotations around the axis of the spontaneous magnetization), leading to a slight (spin) nematic order, e.g., different values of the nearest-neighbour S x S x and S y S y correlations. This effect becomes more severe with increased frustration. For example, for most of the states with D > 3 and J 2 0.3, there is a sensible (e.g., 5 − 10% and even larger) difference in the correlation lengths corresponding to the transverse directions. Connected to this issue, we observe that it is possible to stabilize distinct "families" of local minima for various bond dimensions D, in particular D = 3 and 4, with substantial differences in their magnetization, correlation length, and the degree of nematic order. Every family corresponds to a specific way the quantum fluctuations are built on top of the classical Néel state, e.g., by converging towards one of the possible choices of U (1) charges or breaking the symmetry completely. Given the limited number of bond dimensions that are available within our AAD optimization, it is then of utmost importance to identify the family of minima that are connected and lead to a smooth and physically sound extrapolation in the D → ∞ limit. Therefore, using the scheme introduced in Sec. II B, we take the optimized and almost U (1)-symmetric states from unrestricted simulations (typically for J 2 ≈ 0) and infer their charge structure. The charges revealed by this analysis are listed in Table I and define the correct classes of C 4vsymmetric U (1) iPEPS for D ranging from 2 to 7, which best describe the Néel phase.
In order to obtain the thermodynamic estimates of the ground-state energy and magnetization (within the magnetically ordered phase), we compute these quantities for increasing values of the bond dimension D. A brute-force extrapolation in 1/D provides poor estimates, given the fact that the data are usually scattered, see for example the case of the magnetization reported in Appendix B. Instead, we follow the recent proposal that has been put forward in Refs. [39,40]. In this respect, for every value of D used, we compute the dominant correlation length ξ which is defined by the so-called transfer matrix E of iPEPS, see Fig. 3: where λ 1 is the second largest eigenvalue of the transfer matrix (without the loss of generality we assume that the largest one is normalized to 1). We remark that the value of ξ obtained in this way coincides with the correlation length of the usual spin-spin correlation function (or, more precisely, the transverse correlations): where ρ (2) (r), defined in Fig. 3, is the two-point RDM.
To obtain the χ → ∞ limit of the correlation length, we use the scaling formula [39,44]: which allows for more precise extrapolation of ξ than the usual 1/χ scaling across all ratios of J 2 /J 1 [45]. Finally, the thermodynamic estimates of the energy and magnetization (squared) are obtained by a suitable fit in powers of 1/ξ: where m = |Tr[ρ (1) S]| and ρ (1) is the single site RDM. Let us start discussing the ground-state energy, shown in Fig. 4. For the unfrustrated case J 2 = 0, our results are fully compatible with what has been previously obtained in Refs. [39,40]. The data points align perfectly according to the theoretical expectations and the extrapolated values are in very good agreement with quantum Monte Carlo results [4,5]. For example for D = 7 (after extrapolation in the environment dimension χ) we get e(D = 7) = −0.669432, which is identical to the linear extrapolation in 1/ξ 3 from D = 3 to 7. Including the subleading term 1/ξ 4 , the extrapolation gives e(∞) = −0.669437(2) (to be compared with the exact value e QMC = −0.669437(5) [4]). For future comparisons, the energies for D = 7 and different J 2 /J 1 ratios are reported in Table II. Under increasing the frustrating ratio, a remarkably smooth behavior persists up to J 2 /J 1 ≈ 0.3; then, for larger values, small fluctuations on the fourth digit of the energy, are visible, possibily indicating that the scaling regime moves to larger values of ξ (or D), not reachable within our current possibilities. Still, the quality of the results is sufficient to obtain reliable extrapolations for ξ → ∞. Our calculations show that the expected scaling is not limited to the unfrustrated case, but persists in the whole antiferromagnetic region, thus corroborating the ideas put forward in Refs. [39,40]. One remarkable feature is that, while for small values of D (i.e., for D = 2 and 3), the correlation length ξ clearly increases by increasing J 2 /J 1 , for larger values of D (i.e., for D = 4, 5, 6, and 7), it is essentially constant, or even slightly decreasing with J 2 /J 1 . This aspect will be discussed in connection to the magnetization curve that is presented below.
Then, we move to the central part of the present work, which deals with the magnetization, see Fig. 5. Here, we report m 2 (ξ) for different values of J 2 /J 1 (including 0.5) for D ranging from 2 to 7. Furthermore, the   [4]. For comparison, the variational Monte Carlo calculations of Ref. [46] are also included.
raw data for D = 7 are also shown in Table II. In the unfrustrated case, we get m 2 (D = 7) = 0.0994 and m 2 (∞) = 0.0948 (2), to be compared with the exact value m 2 QMC = 0.0942(2) [4]. In Fig. 5, we attempt both linear and quadratic fits. As in the case of energy extrapolations, we exclude the results with D = 2 from the fitting procedure, since they are clearly off, especially for intermediate and large values of J 2 /J 1 . According to our fits, the linear one looks more trustable than the quadratic one, which serves to give an upperbound to the value of the magnetization. Within the linear fit, we observe vanishing magnetization for J 2 /J 1 ≈ 0.46 (1), giving rise to a continuous transition to a magnetically disordered phase, whose nature is beyond the scope of the present work. We would like to emphasize that the results for J 2 /J 1 = 0.5 are clearly incompatible with a smooth be- havior in 1/ξ, strongly suggesting that at this point the ground state is already outside the magnetically-ordered phase. The final magnetization curve is shown in Fig. 6. For comparison, the variational Monte Carlo calculations, which have been obtained by using Gutzwillerprojected fermionic states, are also shown [46]. In the latter case, a quantum critical point for J 2 /J 1 ≈ 0.48, separating the antiferromagnetic phase and a gapless spin liquid, has been reported. The present results are expected to improve the accuracy of the magnetization (e.g., the accuracy of m 2 for the unfrustated case is smaller than 1%). Still, these two independent calculations give very similar behavior, with almost compatible values for the location of the quantum critical point. We would like to mention that, recent numerical calculations, including DMRG [19], neural-network approaches (based upon restricted Boltzmann machines on top of fermionic states) [21], and finite size PEPS calculations [47] also pointed out that the Néel phase survives up to J 2 /J 1 in the range 0.45 ÷ 0.47, a value that is considerably larger than the one predicted by linear spin-wave theory [48].
Finally, we would like to comment on the J 2dependence of the correlation length, which is clearly different at small (i.e. D = 2, 3) and larger (i.e. D = 4, · · · , 7) bond dimensions. A possible explanation of the rapid increase of ξ, for D = 2 and 3, when approaching the critical point, may be attributed to the fact that, for these very small bond dimensions, the antiferromagnetic state is poorly approximated as a "dressed" product state, having a finite magnetization but lacking the correct transverse (Goldstone) fluctuations. When approaching the phase transition, the magnetization decreases and the state starts to build up long-range entanglement (for D = 3 a short-range resonating-valence bond state can be constructed [49]). Therefore, a larger correlation length can be attained. Once the basic (low D) structure of tensor is established, optimizing at increasingly higher D further improves the description of the antiferromagnetic state and allows correlation length to grow, becoming large even in the presence of significant frustration. Then, no appreciable change of ξ is detected when approaching the quantum critical point. In this respect, we expect that ξ → ∞ in the whole Néel phase, including the critical point. Remarkably, despite optimized iPEPS being finitely correlated, the correct exponent of the power-law decay of transverse spin-spin correlations, i.e., S x 0 S x r 1/r (assuming magnetization along z-spin axis), can already be obtained, see Appendix C for the case with J 2 = 0.
As mentioned above, ξ corresponds to the correlation length of transverse spin-spin correlations. In addition to that, it is possible to evaluate, by a direct fitting procedure of the correlation function itself, the correlation length ξ L of the longitudinal correlations. We find also this quantity to be relatively large, i.e., ξ L ≈ ξ/2, see Fig. 7. Moreover, as for transverse spin-spin correlations, the short-range behavior of the longitudinal correlations reveals their power-law decay (see Appendix C), which then becomes rapidly cut off above the finite-D induced length scale ξ L . These findings show that our optimized iPEPS are even able to approximately capture the powerlaw behavior of transverse and longitudinal spin-spin correlations of the Néel phase.

IV. CONCLUSIONS
In this work, we have investigated the antiferromagnetic phase of the spin-1/2 J 1 − J 2 model on the square lattice, evaluating with unprecedent accuracy the energies and magnetizations for J 2 /J 1 ≤ 0.45. The results point towards the existence of a quantum critical point at J 2 /J 1 ≈ 0.46 (1), which separate the Néel antiferromagnet and a quantum paramagnet, whose nature is beyond the scope of the present study. The importance of our findings is twofold. From the methodological side, we combined state-of-the-art optimization techniques (based upon the AAD scheme [37]), clever parametrizations of the tensor network (based upon the underlying residual U (1) symmetry that exists in the Néel phase), and recently developed extrapolation analyses (based upon the correlation-length scaling [39,40]). In particular, the construction of U (1)-symmetric tensor is pivotal to a straight optimization procedure and correlation-length scaling to solid extrapolations to thermodynamic limit. With these tools in hand, it is possible to get reliable estimations for the ground-state energy but, most importantly, also for the magnetization within the frustrated regime, for which no exact methods can be applied. Therefore, the main outcome of the present work is to provide the magnetization curve for the spin-1/2 J 1 − J 2 model on the square lattice up to relatively large values of the frustrating ratios. In particular, the magne-tization curve shows a smooth behavior, which strongly suggest the existence of a continuous phase transition towards a quantum paramagnet.
Here, our calculations have been limited to the magnetically ordered phase, where relatively entangled states have been achieved. Indeed, rather long correlation lengths are obtained, indication that the tensor network may approximately describe the existence of gapless excitations (i.e., Goldstone modes). The magnetically disordered phase still remains elusive, presumably because of its high-entangled nature due to fractional excitations (spinons and visons). In this respect, the recentlydeveloped method to impose SU (2) symmetry [34,50] would be beneficial to the final understanding of the full phase diagram of the spin-1/2 J 1 − J 2 model.

ACKNOWLEDGMENTS
We thank Fabien Alet, Zhengcheng Gu, Andreas Läuchli, Wenyuan Liu, Pierre Pujol, Anders Sandvik, and Sandro Sorella for helpful discussions. This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2020-P1231. This project is supported by the TNSTRONG ANR-16-CE30-0025 and the TNTOP ANR-18-CE30-0026-01 grants awarded by the French Research Council. be described by a simple linear function in 1/D with appreciable accuracy. Considerable deviations are present, especially for larger bond dimensions across the studied range of J 2 /J 1 , preventing a smooth extrapolation in the D → ∞ limit.
Appendix C: Spin-spin correlations in the J2 = 0 limit In Fig. 9, assuming magnetization along z-spin axis, we show the decay of both transverse S x 0 S x r and longitudinal S z 0 S z r correlations for J 2 = 0 and D = 2, . . . , 7. Due to imposed U (1) symmetry the transverse correlations along x and y spin axes, S x 0 S x r and S y 0 S y r , are identi-cal. The extrapolated values are obtained by performing, for each distance r, an extrapolation in 1/ξ of the finite-D results using the three largest available bond dimensions D = 5, 6, and 7. Then, the extrapolated correlations are fitted in the short-distance region r ∈ [2,11] (excluding the nearest-neighbor case) with a power law f (r) ∝ r −β . The final result gives an exponent β ≈ 1.02(1) for transverse and β L ≈ 1.90(5) for longitudinal correlations. Linear extrapolations in 1/ξ, up to r = 20, are performed using the D = 5, 6, 7 data. The dashed lines are power-law fits to short-distance behavior, see text.