A tensor network study of the complete ground state phase diagram of the spin-1 bilinear-biquadratic Heisenberg model on the square lattice

Using infinite projected entangled pair states, we study the ground state phase diagram of the spin-1 bilinear-biquadratic Heisenberg model on the square lattice directly in the thermodynamic limit. We find an unexpected partially nematic partially magnetic phase in between the antiferroquadrupolar and ferromagnetic regions. Furthermore, we describe all observed phases and discuss the nature of the phase transitions involved.


Introduction
Strongly correlated materials-even those described by relatively simple Hamiltonians-can display very rich behavior. In frustrated magnets, for example, classical ground state [1] degeneracies due to the inability of the system to minimize each term in the Hamiltonian individually give way for quantum fluctuations to play the predominant role in determining the actual ground state of the system. As such, frustrated magnets are a playground for all kinds of exotic states of matter, such as spin liquids [2,3] and spin-nematic states [4,5], the latter of which break spin-rotational symmetry while preserving time-reversal symmetry.
In this paper, we study the spin-1 bilinear-biquadratic Heisenberg (BBH) model on the square lattice: a model that consists of spin-1 particles, one per lattice site, that interact only with their nearest neighbors. The nearest-neighbor two-particle interaction is a combination of two types of competing interactions: the ordinary bilinear Heisenberg coupling, and the biquadratic coupling, which is the Heisenberg coupling squared. The corresponding coupling constants appearing in the Hamiltonian are commonly parametrized by an angle θ: where S i = (S x i , S y i , S z i ) is the vector of spin-matrices for the spin-1 particle on site i, The BBH model on the square lattice has been subject to a lot of interest in recent years, for a number of reasons. For one, at θ = π/4, the Hamiltonian is equivalent to the SU(3) Heisenberg model [6,7], which can be realized using cold-atom experiments [8][9][10][11]. Second, it was proposed that the nematic phases of the BBH model on the triangular lattice could be related to the unusual behavior of NiGa 2 S 4 [12][13][14] and Ba 3 NiSb 2 O 9 [15][16][17].
From a theoretical perspective, the BBH model is the most general lattice-translation, lattice-rotation and spin-rotation-symmetric Hamiltonian with nearest-neighbor interactions [18]. Moreover, the more interesting part of the phase diagram-the top half 0 < θ < π-is beyond the reach of quantum Monte Carlo [19] because of the sign problem, and has up to now only been studied either for small system sizes using exact diagonalization [20], or on larger systems using semi-classical [5,20] and perturbative [21] approaches.
In this article we continue upon our work on the emergence of the Haldane phase in the spin-1 BBH model on the square lattice [22] and provide a complete study of the ground state phase diagram using state-of-the-art tensor network algorithms. Contrary to previous studies, we find the occurrence of a half nematic half magnetic 'm = 1/2' phase that was previously predicted to appear only in the presence of an external magnetic field [20]. Additionally, we describe all phases and determine the nature of the corresponding phase transitions. The ground state phase diagram of the spin-1 BBH model according to our computations, which in itself is the main results of this paper, can be found in Fig. 8.
This paper is organized as follows. In Section 2 we provide a brief description of the numerical method used. Thereafter, Section 3 sets the stage by discussing some of the relevant properties of spin-1 particles-in particular in the context of spin-nematic order-as well as the symmetries of the Hamiltonian, followed by an overview of previous studies of the BBH model in Section 4. We then proceed to discuss our own findings and compare those with previous results in Section 5, and conclude with Section 6.
An iPEPS consists of five-legged tensors, one per lattice site, each with four auxiliary legs connecting to the four neighboring tensors in accordance with the square lattice structure, and one physical leg corresponding to the local Hilbert space of a single spin-1 particle (see Fig. 1). The auxiliary vector spaces all have the same dimension, called the bond dimension (D). The bond dimension controls the accuracy of the ansatz: the higher D, the more entanglement can be captured by the iPEPS, yielding better approximations to the ground state as D increases. In practice, we run simulations up to at most D = 10 (D = 16 if the state has U(1)-spin-rotation symmetry) and then extrapolate D → ∞, or, equivalently, truncation error to zero [50], when computing expectation values.  If the ground state is translationally invariant, the iPEPS can be represented by a single tensor repeated all over the lattice. However, if translational invariance is (partially) broken, then we have to use an iPEPS with a larger unit cell to reproduce the ground state; e.g., the antiferromagnetic ground states requires a 2x2 unit cell with two different tensors [51].
One of the main obstacles in dealing with infinite two-dimensional tensor networks is that, contrary to the one-dimensional case, it is impossible to contract the infinite network exactly. Instead, several approximation schemes have been developed over the last years, such as: transfer-matrix-based contraction schemes [23], coarse-graining-based contraction schemes [52,53] and corner-transfer matrix (CTM) methods [54][55][56] based on a formalism derived by Baxter [57,58]. For this paper, we have used a modified version [42,59] of the CTM algorithm developed by Orús and Vidal [60].
In the CTM contraction scheme, for each tensor in the unit cell additional environment tensors are introduced that represent the contraction of all other tensors surrounding the tensor in question. These environment tensors come with their own environment bond dimension χ, and are obtained as a fixed point of a row and column insertion and truncation procedure [60]. For this paper, we have chosen χ(D) > D 2 to be large enough to yield negligible variations in energy compared to those due to the use of finite D.
Given a randomly initialized iPEPS |ψ with a predefined unit cell, we evolve |ψ in imaginary time by applying exp(−τ H) for large enough τ to project it onto the ground state. At the cost of a controllable Trotter error, using a second order Trotter decomposition the imaginary time evolution operator exp(−τ H) can be decomposed into a sequence of gates that evolve a single bond of |ψ a small step in imaginary time. Application of each such gate increases the bond dimension of the updated bond. Truncating the bond dimension back to D after updating a single bond can be done using the simple or the full update algorithm.
The simple update algorithm [53] is inspired by its one-dimensional counterpart, known as infinite time-evolving block decimation (iTEBD) [61]. In one dimension, cutting a single bond corresponds to splitting the Hilbert space into two, and truncating the bond dimension back to D amounts to keeping the D largest Schmidt values corresponding to the cut in question, a procedure that yields the best rank-D approximation to the updated state in the 2-norm. Similarly, the simple update keeps the D largest values of a singular value decomposition involving only the tensors and accompanying weights connected to the to-be-truncated bond, which in this case does not correspond to an actual Schmidt decomposition, because cutting a single bond does not split the Hilbert space into two. Regardless, the simple update, albeit being somewhat ad hoc, is computationally very efficient and yields reasonably accurate results in many cases.
In the full update scheme [24] the truncated tensors are found by variationally minimizing the squared distance to the iPEPS with the updated bond. This requires computing the environment at every step, unless the fast full update [62] is used, which speeds up the algorithm by recycling previous environments. The full update algorithm is computationally more costly than its simple counterpart, but it does not suffer from uncontrollable approximations, and it is systematic, in the sense that in the D → ∞ limit the iPEPS should converge to the true ground state of the system in question.

The model 3.1 Spin-1 nematic order
Individual spin-1 particles can display more types of order than their spin-1/2 counterparts. For example, every spin-1/2 single-particle state-assuming it is a pure state, which we take all single-particle states in this subsection to be-has a fixed magnetization m = S x 2 + S y 2 + S z 2 of exactly 1/2 (setting = 1), and is fully described by its magnetic dipole moment S = (S x , S y , S z ). This is a consequence of the fact that any spin-1/2 single-particle state can be obtained by applying a specific rotation to the up state in the z-basis (or any other reference state for that matter). By contrast, spin-1 particles can have a magnetization of anywhere between 0 and 1. Moreover, since spin-1 particles live in a higher-dimensional space, the dipole moment is not sufficient to fully describe a given spin-1 single-particle state, for which in addition also the quadrupole moment is required.
The quadrupole moment is measured by products of spin operators, which we can combine into the matrix T αβ = S α S β , with α, β ∈ {x, y, z}. However, T contains more than just the quadrupole moment, as the anti-symmetric part of T is just the vector of spin operators itself due to the spin commutation relations. The remaining symmetric part of T can, because Tr(T ) = S(S + 1) = 2 is constant, be captured by five independent operators that can conveniently be organized into the vector of quadrupolar operators. Typically, spin-1 single-particle states display both magnetic and quadrupolar order. However, when referring to a quadrupolar state, what we mean is a state that is solemnly described by Q. In particular, we require that its magnetic dipole moment m is zero. An example of a quadrupolar state is the |0 state in the S z -basis. As can be checked directly: m = 0. Interestingly, even in the absence of magnetization, the |0 state does break spin-rotation symmetry because (S z ) 2 = 0, whereas the fluctuations in x and y-direction are non-zero: (S x ) 2 = 1 = (S y ) 2 . Classically, this state can be thought of as a dipole moment that fluctuates in the x-y plane in such a way that it is zero on average. We will picture it by a disc portraying the plane of fluctuations (see Fig. 2). A normal vector to the plane of fluctuations is called a director (±e z in the case of |0 ). Any quadrupolar single-particle state is fully described by the orientation of its directors.
x y z Figure 2: The |0 quadrupolar state. The disc represents the magnetic dipole fluctuations in the x-y plane, and the dotted line represents a director perpendicular to the plane of fluctuations pointing along the z-axis.
A convenient on-site basis for investigating quadrupolar order is the time-reversal invariant basis: where {| ↑ , |0 , | ↓ } is the standard S z -eigenbasis. Invariance of the above vectors under the time-reversal operator T , which is anti-unitary, follows immediately from the fact that T interchanges | ↑ and | ↓ , and adds a sign to |0 : i.e T | ↑ = | ↓ , T | ↓ = | ↑ , and T |0 = −|0 . Any real linear combination of the above basis states α=x,y,z u α |α is also a quadrupolar state, with a director pointing along u. More generally, an arbitrary spin-1 single-particle state can be expanded in the basis of (3) as |d = α=x,y,z d α |α , where d has real and imaginary parts d = u + iv. In terms of u and v, the magnetic moment of |d is given by S = 2 u × v. We can normalize the state and use global phase invariance to have u and v satisfy u 2 + v 2 = 1 and u · v = 0. Assuming the last two equations hold, then so does the following: fully magnetized m = 1 states are precisely those for which u = v, whereas non-magnetic m = 0 states correspond to (u, v) = (0, 1) or (1, 0). The latter case (m = 0) describes purely quadrupolar states with a director d pointing along the direction of whichever of u or v is non-zero. Whenever u and v are both non-zero but not of equal magnitude, and consequently 0 < m < 1, the state is of mixed character: neither fully magnetic nor purely quadrupolar.
Because the basis of (3) is invariant under the anti-unitary time-reversal operator T , we can immediately conclude that the time-reversal invariant single-particle states are precisely those for which the coefficients in the above basis do not change under T , up to a global phase. Since T is anti-unitary, this means that the time-reversal invariant states are precisely those for which d is either completely real (v = 0), or purely imaginary (u = 0); i.e. |d is a quadrupolar state. (This is most obvious in the above gauge u · v = 0, where u and v both being non-zero results in S being flipped under T .) In other words, following Andreev's definition [4] of a spin-nematic state as a state that breaks spin-rotational symmetry but preserves time-reversal symmetry, we conclude that, for a single spin-1 particle, the notions of quadrupolar and spin-nematic are equivalent.
Comparing to nematic order in liquid crystals, which is the alignment of rod-like particles in a liquid, we observe that quadrupolar states have the exact same symmetry properties as rods. Indeed: two directors d and −d correspond to the same quantum state (they differ by an overall phase). Contrast this to magnetic states, where the magnetic moment is described by the magnetic dipole vector, which is certainly not equal to minus itself. More explicitly, the magnetic the | ↑ and the quadrupolar |0 states are invariant (the first up to a phase) under rotations about the z-axis, but, in addition, the |0 state is (up to a minus sign) also invariant under a π-rotation about any axis that lies in the x − y plane. Finally, we speak of nematic order when neighboring directors order in space, similar to the spatial ordering of rods in liquid crystals.

High-symmetry points
It shall not come as a surprise to the reader that the biquadratic term (S i · S j ) 2 in the Hamiltonian (1), which involves on-site products of spin operators, is related to quadrupolar order. However, part of the biquadratic term also includes the ordinary spin matrices (the anti-symmetric part of T defined above), which are more naturally absorbed into the bilinear term. To better separate the two competing magnetic and quadrupolar interactions, let us rewrite the Hamiltonian in terms of S and Q. Doing so yields up to an irrelevant θ-dependent constant, with the spin and quadrupolar coupling constants given by J S (θ) = cos(θ) − sin(θ)/2 and J Q (θ) = sin(θ)/2. We should remark that quadrupolar order is formally captured by the symmetric part of T , and the set of operators in (2) is just one of many possible sets of operators that can be chosen to describe quadrupolar order. However, using the set of operators in (2) does unveil the at first sight hidden points of higher symmetry that the BBH Hamiltonian possesses. When expressing all three spin matrices and all five quadrupolar operators given by (2) in terms of the time-reversal invariant basis of (3): we observe that the spin matrices become equal to the three imaginary Gell-Mann matrices, whereas the above five quadrupolar operators equal the remaining five real Gell-Mann matrices. Recall that the Gell-Mann matrices are the generators of SU(3). Consequently, whenever the spin coupling constant equals the quadrupolar coupling constant, i.e. J S (θ) = J Q (θ) = J, the Hamiltonian reduces to an equal-weighted product over all Gell-Mann matrices: [63]. This occurs for θ = π/4 and θ = 5π/4. Additionally, when J S (θ) = −J Q (θ) = J, which happens for θ = π/2 and θ = 3π/2, all matrices again have equal weight, except that the imaginary matrices come with a relative extra minus sign. Since the square lattice is bipartite, we can compensate for this extra sign by taking the antifundamental representationλ = −λ * of SU(3) on one sublattice. The Hamiltonian then takes the form H = J i,j λ i ·λ j (where i is in the A-sublattice, and j in the B-sublattice), which is also SU(3) symmetric.
At the SU(3)-symmetric points, there is a larger set of operators commuting with the Hamiltonian, which means that the ground state degeneracy increases. Moreover, because the SU(3)-symmetric points are exactly the points at which the coupling constants are equal in magnitude (|J S | = |J Q |), the regions in between are precisely those for which one of the two coupling constants dominates the other. Hence, we can naively expect to have four different phases separated by the SU(3)-symmetric points, with magnetic or quadrupolar order depending on which of the two coupling constants is larger in magnitude.

Previous studies
In 1988 Papanicolaou [5] conducted a product-state analysis of the spin-1 BBH model, a concise summary of which can be found in [20]. Because the square lattice is bipartite, finding the product ground state reduces to a two-body problem (one product state per sublattice). Papanicolaou's results agree with our analysis above; he found four different phases separated by the SU(3) points: the familiar ferromagnetic (FM) and antiferromagnetic (AFM) phases, and in between a ferroquadrupolar (FQ) phase wherein neighboring sites have aligned directors, and a phase he called semi-ordered (see Fig. 3). The latter has a degenerate product ground state: minimizing the two-particle energy leads to the condition that one sublattice state has to be quadrupolar, whereas the other sublattice state can be either quadrupolar with a director perpendicular to that of the first sublattice, or magnetic with a magnetic moment aligned with the director of the first sublattice, or a combination of the two.
As an interesting side note, let us mention that exactly at each SU(3) point both product ground states left and right of the SU(3) point in question are simultaneous ground states of the system, and can be rotated into one another through some element of SU(3).
The product ground state degeneracy in the top right octant of the phase diagram means that we can expect quantum fluctuations to play an important role there. Moreover, the types of ground states found open up a lot of possibilities for interesting types of order [64]; for example: the fact that neighboring directors can be perpendicular in three different ways (e.g., pointing in x, y and z direction, respectively), allows for the possibility of a threesublattice ground state (Fig. 4), which is surprising for a model with nearest-neighbor inter-  : Different types of product ground states in semi-ordered phase (π/4 ≤ θ ≤ π/2). Left: three-sublattice order with neighboring directors being perpendicular (antiferroquadrupolar 'AFQ3' state-discs with differing colors are quadrupolar states viewed at different angles); right: two-sublattice order with directors aligning parallel to neighboring magnetic moments (half nematic half magnetic 'm = 1/2' state). Note that the color coding is only meant to highlight which particles are in the same quantum state; see also [51] .
The energy per site for the product ground state is plotted in Fig. 5: it shows clear kinks at the SU(3)-symmetric points, suggesting that the phase transitions are all of first order.
The lower half −π ≤ θ ≤ 0 of the phase diagram is free of the sign problem and has been studied using quantum Monte Carlo in 2002 by Harada and Kawashima [19], who found that the actual ground state phase diagram agrees with the product ground state phase diagram. Moreover, they observed jumps in the quadrupolar order parameter at the SU(3) points θ = 5π/4 and 3π/2, which agrees with the product state analysis and confirms that the phase transitions from the antiferromagnetic to the ferroquadrupolar and from the ferroquadrupolar to the ferromagnetic phases are both of first order.
The Monte Carlo study was followed up by an exact-diagonalization and linear-flavor-wave study at the SU(3)-symmetric point θ = π/4 by Tóth et al. [6] in 2010, who concluded that the ground state is actually a three-sublattice state that is favored over the two-sublattice state because it has a lower zero-point energy (the energy of the wave-excitations on top of the classical ground state averaged over the Brillouin-zone). The three-sublattice order of the ground state at the SU(3) point was confirmed by DMRG and iPEPS simulations by Bauer et al. [7]. A further study by Tóth et al. [20] in 2012 based on exact diagonalization and linearflavor-wave theory examined the question whether the above-mentioned three-sublattice phase extends beyond the SU(3) point θ = π/4. They found that, in the semi-ordered phase, the ground state is the three-sublattice antiferroquadrupolar (AFQ3) state with perpendicular directors on neighboring sites (the analogous product state is depicted in Fig 4: left). Interestingly, Tóth et al. also found that the three-sublattice pattern extends into the antiferromagnetic regime, where, for 0.2π θ < π/4, instead of ordinary antiferromagnetism, a 120 • magnetic order emerges (Fig. 6). This result was thereafter confirmed by series expansion [21].

Overview
To get a first rough picture of the ground state phase diagram, and also as a consistency check, we initialize D = 1, 2, 3, . . . , 6 simple update simulations for θ at 40 evenly spaced points covering all of [0, 2π] using 2x2 and 3x3 unit cells ( Fig. 7; only D = 1, 2 and 6 shown). The top plot of Fig. 7 shows the energy per site, which is seen to decrease as the bond dimension increases, except in the FM region π/2 < θ < 5π/4, where the ground state can be represented by a product state (i.e. D = 1). Moreover, we observe a transition from a twosublattice to a three-sublattice ground state around θ ≈ 0.2π, with the 0.2π three-sublattice simulations showing 120 • magnetic order observed in the AFM3 phase.
The bottom plot of Fig. 7 shows the magnetization per site (for D = 6 simulations), which is of order one in the FM and AFM phases, and zero in the quadrupolar phases, as expected. Contrary to magnetic order, quadrupolar order cannot be captured by a single order parameter, because it is described by a matrix: the traceless symmetric part of S α S β , denoted by Q αβ [65]. The convention used for the quadrupolar operators in the vector Q given by Eq. (2) has a preferred choice of direction in spin-space, as can be seen from the first two components of Q which represent the diagonal part of Q αβ . Using Q works well for nematic order along the z-direction; however, we do not control the direction of spontaneous symmetry breaking in our simulations. Therefore, when measuring quadrupolar order, we prefer to work with invariants of Q αβ such as its eigenvalues, or equivalently, the matrix invariants I Q = tr(Q) = 0, II Q = 1 2 tr(Q) 2 − tr(Q 2 ) = − 1 2 tr(Q 2 ) and III Q = det(Q) rather than the individual components of Q. Fig. 7 shows the two non-zero matrix invariants II Q and III Q , which are clearly larger in magnitude in the quadrupolar phases than they are in the magnetized phases [66]. Moreover, the determinant III Q changes sign when going from a nematic to a magnetic state.
The D = 1, 2, . . . , 6 simple update simulations reproduce the AFM, AFQ3, FM and FQ phases separated by the SU(3) points, and hint at the existence of the AFM3 phase [67] found by Tóth et al. [20] in between the AFM and AFQ3 phases. Because the lower half −π < θ < 0 of the phase diagram has already been established by quantum Monte Carlo [19], and in the top left quadrant π/2 < θ < π we only found ferromagnetic translationally invariant product states that have the same energy independent of D or unit cell size, the only section of the phase diagram that remains to be thoroughly investigated is the top right quadrant 0 < θ < π/2. In order to obtain more accurate data on this interesting region, we will use finer-spaced higher-D full update simulations rather than the less accurate simple update simulations shown in Fig. 7.
Part of the top right quadrant we have already tackled in our recent paper [22], where we attempted to determine the extent of the AFM3 phase, and in the process stumbled upon the paramagnetic (extended) Haldane phase that appears in between the AFM and AFM3 phases. Next, we turn our attention to the quadrupolar part of the top right quadrant, and find yet another unexpected phase: the partially nematic partially magnetic 'm = 1/2' phase that we shall discuss in detail below.
Combining all of our iPEPS results, we arrive at the phase diagram shown in Fig. 8. Note that the Haldane, AFM3 and m = 1/2 phases are not visible in the simple update plots in Fig. 7 because the θ-grid used is too coarse, and, in case of the Haldane phase, the bond dimension D is too low. In the following section, we will go through all phases in Fig. 8 in anti-clockwise order and discuss their properties, and, where necessary, elaborate on our findings. In the section thereafter, we will discuss the nature of all occurring phase transitions-also shown in Fig. 8. The analysis presented below involves simulations with a bond dimension as high as D = 10 (D = 16 for simulations in the Haldane phase), followed by a D → ∞ extrapolation where necessary. As a consequence, the expectation values may differ slightly from those shown in Fig. 7 which are only meant to give a rough overview of the type of order that can be expected.

AFM
The familiar antiferromagnetic phase can be described by a 2x2 two-sublattice unit cell, with spins pointing in opposite directions on each sublattice. The magnetization per site varies  throughout the phase from m ≈ 0.5 − 0.6 close to the FQ phase (Fig. 14) to m ≈ 0.8 around θ = 0, after which it monotonically decreases to zero or almost zero when approaching the Haldane phase [22]. States in the AFM phase are U(1)-spin-rotation symmetric around the axis of magnetization.

Haldane
As described in our recent paper [22], we set about to nail down the precise value of θ c for which the proposed antiferromagnetic to 120 • magnetically ordered phase transition was supposed to occur. While doing so, we discovered that, as we increased θ within the AFM phase, the magnetization vanished before the transition to the AFM3 phase, hinting at the existence of an intermediate phase in between the AFM and AFM3 phases. Further investigation showed that the lowest energy per site was obtained by a 1x1 unit cell iPEPS which converged to a paramagnetic state that preserves spin-rotation (hence also time-reversal) and lattice translation symmetry, but breaks lattice rotation symmetry-in the sense that there is an energy difference between the x and y-bonds-reminiscent of a system of coupled one-dimensional chains. This led us to investigate the anisotropic BBH model, in which we introduced an additional coupling parameter (0 ≤ J y ≤ 1) for the y-bonds only that allowed us to interpolate between the isotropic two-dimensional system and the one-dimensional system of decoupled chains. We then identified the intermediate paramagnetic phase as the Haldane phase, by showing that it can be adiabatically connected to the one-dimensional Haldane phase of decoupled spin-1 BBH chains by sending J y to zero. As a side note, we should remark that the Haldane phase we find is not related to the two-dimensional AKLT state. The latter can be obtained by taking four spin-1/2 particles per site-each forming a singlet bond with a spin-1/2 particle on a neighboring site-and projecting onto the on-site spin-2 subspace, which results a state that has a spin-2 particle on each lattice site. Rather, the analogous picture for the state we find is that of (decoupled) one-dimensional AKLT chains, each having two spin-1/2 particles per site that form singlet bonds only in one direction before projecting onto the on-site spin-1 subspace.
The location of the AFM to Haldane phase transition (at θ = 0.189(2)π) we have determined by investigating for what value of θ the magnetization per site in the AFM phase vanishes. The phase transition from the Haldane to the AFM3 phase (at θ = 0.217(4)π) follows from an energy per site comparison of simulations in both phases, in the same spirit as the AFQ3 to m = 1/2 phase transition discussed below.

AFM3
The three-sublattice 120 • magnetically ordered phase partially breaks translational symmetry in a 3x3 unit cell pattern, and has the spins on the three sublattices aligned in a plane such that they are at relative angles of 120 • (see Fig. 6). The magnetization per site varies from m ≈ 0.4 at θ = 0.22π [22] to m ≈ 0.1 − 0.3 near θ = π/4 (Fig. 10). AFM3 states have no residual spin-rotation symmetry.

AFQ3
Similar to the AFM3 phase, the three-sublattice antiferroquadrupolar phase is also described by a 3x3 unit cell, but now the magnetization per site is near zero at π/4 (Fig. 10, [68]), and exactly zero from about θ ≈ 0.27π up to the m = 1/2 phase (see Fig. 11 for θ = 0.487π). At the product state level, in addition to being time-reversal symmetric, the AFQ3 state has a residual spin-rotation symmetry given by π-rotations around the (mutually perpendicular) axes of nematic polarization.

m = 1/2
Next, let us focus on the extent of the three-sublattice AFQ3 phase in the semi-ordered region π/4 < θ < π/2. Recall that, at the product state level, there is a multitude of distinct types of ground states, two of which are depicted in Fig. 4. Tóth et al. [20] predicted that the three-sublattice (AFQ3) state is the ground state for π/4 < θ < π/2. However, their product state analysis in this same region shows that the addition of an infinitesimal external magnetic field favors the two-sublattice half nematic half magnetic m = 1/2 state as the ground state. Augmented by linear-flavor-wave theory, they continue to show that the AFQ3 phase extends into the small but non-zero magnetic field h > 0 region for π/4 < θ < π/2, but that, as θ approaches π/2 from below, the phase boundary between the m = 1/2 and AFQ3 phases moves down to h → 0 as θ → π/2. Since their exact diagonalization results for the AFQ3 phase with non-zero magnetic field do not extrapolate well to infinite system size, and linear-flavor-wave theory is semi-classical, it is not clear what happens close to θ = π/2. We have thoroughly investigated the parameter range of 0.48π < θ < π/2 using iPEPS. To get an idea of what types of states to expect, we initialized simulations using eleven different types of unit cells (to also allow for possibilities such as stripe order, as well as the more conventional types of order that can be represented by unit cells up to size 4x4 [69]) and evolved them in imaginary time using the simple update algorithm up to bond dimension D = 6. We then picked the minimal unit cell for each type of state obtained (e.g. the 4x2 state displayed the same pattern as the 2x2 state, so we discarded the former), and evolved those in imaginary time using the full update algorithm up to D = 8. This left two competitive states: the three-sublattice AFQ3 state described by a 3x3 unit cell iPEPS, and the two-sublattice m = 1/2 state described by a 2x2 unit cell iPEPS. Interestingly, the m = 1/2 state turned out to have a lower energy than the AFQ3 state for 0.49π θ < π/2. Finally, we ramped up the bond dimension to D = 10 to nail down the precise location of the phase transition between the AFQ3 and m = 1/2 phases.
As can be seen in Fig. 9, for θ = 0.487π, the AFQ3 simulation is lower in energy than the m = 1/2 simulation, whereas for θ = 0.490π the opposite is true. Taking the error bars from the D → ∞ extrapolation into account, we can conclude that-in contrast to previous predictions-the AFQ3 phase does not persist all the way up to the SU(3)-symmetric point π/2, but remains stable only up to θ = 0.4886(7)π, after which the system transitions into the m = 1/2 phase.  In case the reader is wondering why the m = 1/2 phase occupies such a small portion of the phase diagram, it is insightful to remark that the size in θ-space is not a physical quantity, as [cos(θ), sin(θ)] is but one of the many possible parameterizations of the coupling constants of the Hamiltonian in Eq. (1). Additionally, there is a reason for the m = 1/2 phase to only appear this close to π/2. The actual ground state is not a product state, and on the quadrupolar sublattice a small magnetic moment parallel to that of the magnetized sublattice can be detected. This is only energetically favorable when the Heisenberg coupling parameter J S (θ) in the Hamiltonian expressed in terms of S and Q (see Eq. (4)) becomes of similar magnitude as the quadrupolar coupling parameter J Q (θ), which only happens around 0.49π where J S (θ) decreases rapidly from 0 towards −1/2.
As its name suggests, the magnetization per site of m = 1/2 states is exactly 1/2 throughout the entire phase (e.g., see Figs. 11 and 12 for θ close to the transition points). It has a residual U(1)-spin-rotation symmetry around the axis of magnetization. As θ approaches π/4, the m = 1/2 state gradually turns into a product state (see Section 5.3.5).

FM
The familiar ferromagnetic phase is described by a translationally invariant product state represented by a 1x1 unit cell. It has a residual U(1)-spin-rotation symmetry around the axis of magnetization, and a magnetization per site of exactly m = 1 (see Figs. 7, 12 and 13).

FQ
Finally, the ferroquadrupolar phase is also translationally invariant, and can be represented by a 1x1 unit cell. Its magnetization is zero or very close to zero throughout the phase (see Figs. 13 and 14 for the transition points [68]). At the product state level, it is symmetric under time-reversal, residual U(1)-spin-rotations around the axis of nematic order, and π-rotations around any axis perpendicular to the axis of nematic order. Close to the FM phase, states in the FQ phase are product states, but, as we approach the AFM phase, the FQ states become gradually more and more entangled (see Sections 5.3.6 and 5.3.7).

Nature of the phase transitions
In the following section, we will discuss the nature of the phase transitions in the phase diagram shown in Fig. 8, including the transitions at 0.189(2)π and 0.217(4)π that have already been discussed in [22] (a brief summary of which will be presented below).
From the previous section, we gather that all neighboring phases break different symmetries, which suggests that all phase transitions are either first order, or unconventional second order transitions (see e.g. [70,71]). To support this hypothesis-and where possible distinguish between the two options-we will next look for kinks in the energy per site, or jumps in typical order parameters such as the magnetization or the Q-matrix invariants per site due to varying θ. At θ = 0.189(2)π, we have a transition between the AFM and Haldane phases. In the AFM phase, we observed that the magnetization per site goes to zero when approaching the Haldane phase suggesting a second order phase transition, which is unconventional considering the fact that both states break different lattice translation and rotation, as well as different spinrotation symmetries. Moreover, we did not find clear hysteresis behavior, which supports the claim that this transition is second order. However, due to the error bars in the magnetization close to the transition, we were not able to say with certainty whether the phase transition is a second or a weak first order transition; all we can say for sure is that this is not a clear first order transition as no jump in magnetization or kink in the energy were observed. 5.3.2 Haldane to AFM3: θ = 0.217(4)π The phase transition between the Haldane and the three-sublattice AFM3 phase displays a clear jump in magnetization, which goes from zero (Haldane) straight to 0.4 (AFM3) at the transition point. This implies that the transition is first order.

AFM3 to AFQ3: θ = π/4
At θ = π/4, we have a transition from the three-sublattice 120 • magnetically ordered phase to the three-sublattice antiferroquadrupolar phase. Precisely at the symmetric point θ = π/4, both the 120 • magnetically ordered state and the antiferroquadrupolar state are ground states of the system. Hence, we can simulate both at the phase transition by initializing the simulations from deep within the 120 • magnetically ordered and quadrupolar phases, respectively, and then moving towards the phase transition by initializing each simulation from the previous one. The resulting data at the critical point is shown in Fig. 10.
The left plot of Fig. 10 shows a subtle kink in the energy per site, which we observe for each fixed bond dimension simulation for D = 2, 3, . . . , 8 (only D = 4, 6, 8 shown), supporting the occurrence of a first order transition. The right plot shows the magnetization per site exactly at the phase transition (where we increased D up to 10). Despite the fact that the magnetization data points do not fit on a perfect straight line, a rough extrapolation of D → ∞ shows that the magnetic AFM3 state at the transition has a magnetization of at least 0.1, whereas the quadrupolar AFQ3 state has a magnetization of around zero [68], which agrees with the above observation that this is a first order transition. Because the magnetization does not extrapolate very nicely to D → ∞, we have also investigated the behavior of the eigenvalues of the quadrupole matrix (Appendix: Fig. 15 right), where a jump in spectrum can be observed. Additionally, we have also observed the occurrence of hysteresis (Appendix: Fig. 15 left). All of the above taken together, combined with the fact that, at the product state level, the two phases break different symmetries, lead us to conclude that this must be a first order phase transition.
5.3.4 AFQ3 to m = 1/2: θ = 0.4886(7)π Using the same simulations as the ones shown in Fig. 9 (but now only for θ = 0.487π and θ = 0.490π to prevent the figure from getting too cluttered), we see that the magnetization on both sides of the phase transition for the m = 1/2 states is exactly one-half, whereas it is exactly zero for the antiferroquadrupolar states (Fig. 11). Thus, at the transition from AFQ3 to m = 1/2 at θ = 0.4886(7)π, the magnetization jumps from zero to one-half, which means that the transition is of first order. This notion is confirmed by the right plot of Fig. 9, which shows that, at the intersection, the energy per site graphs for the AFQ3 and m = 1/2 phases have different slopes, indicating a kink in the energy per site at the transition. In this case, the transition is between two magnetized phases: the half-magnetized m = 1/2 phase and the ferromagnetic phase. States in the ferromagnetic phase are product states, which in our simulations can be seen by the fact that the energy does not improve as the bond dimension increases. Hence, any fixed D simulation will reproduce the exact ground state (D = 2 shown in Fig. 12).
As θ approaches π/2 from below, the m = 1/2 state also turns into a product state, which can be seen from the fact that the the energy per site graphs for different values of D all converge to the same value at θ = π/2 (Fig. 12 left). Hence, at the phase transition, no D → ∞ extrapolation is necessary, as the D = 2 results already give the exact ground state.
The left plot of Fig. 12 displays a clear kink in the energy, demonstrating that this is a first order transition. Moreover, Fig. 12 right shows that the m = 1/2 phase has a magnetization per site of exactly 1/2, whereas the ferromagnetic phase is fully magnetized, confirming that the transition is indeed of first order.  The remaining two critical points at 5π/4 and 3π/2 have previously been investigated using quantum Monte Carlo simulations [19]. Harada and Kawashima demonstrated that (S z ) 2 −2/3, which they used as the quadrupolar order parameter, exhibits a jump at θ = 5π/4 and at θ = 3π/2. This indicates that both transitions are of first order; a conclusion that also follows from our simulations, which we present below for completeness. Because the jump in the quadrupole moment has already been established, we will focus on the energy and magnetization in the following, but remark that we also observe clear jumps in the spectrum of the Q-matrix.
5.3.6 FM to FQ: θ = 5π/4 As before, we are dealing with product states on both sides of the phase transition. Thus, we can use any fixed D simulations to investigate the nature of the ferromagnetic to ferroquadrupolar phase transition (D = 2 shown). Both the kink in energy and the jump in magnetization displayed in Fig. 13 confirm that the phase transition at θ = 5π/4 is indeed of first order.

FQ to AFM: θ = 3π/2
When increasing θ from 5π/4 to 3π/2, the ground state becomes more and more entangled while remaining in the ferroquadrupolar phase. Upon reaching 3π/2 the ground state is no longer a product state, and therefore we require higher D simulations to investigate the phase transition at 3π/2.  As shown by Fig. 14, there is a subtle but observable kink in the energy per site as a function of θ. Additionally, the magnetization per site at the phase transition, which does not extrapolate very well to D → ∞, does appear to saturate in between 0.5 and 0.6 in the AFM phase, whereas it extrapolates to zero in the FQ phase, confirming that this transition is also of first order.

Conclusion
We have studied the ground state phase diagram of the spin-1 BBH model on the square lattice by means of infinite projected entangled pair states (iPEPS). Using low bond dimension simple update simulations, we were able to reproduce all phases that were previously predicted to occur [5,20] in this model. We then turned our attention to the challenging top right quadrant of the phase diagram: the part that has the largest product ground state degeneracy, is beyond the reach of Monte Carlo, and had up to now only been studied on small systems with exact diagonalization [20] or by means of semi-classical approaches [20,21].
In this interesting region, we found two new phases: the paramagnetic extended Haldane phase (discussed in our previous work [22]) that preserves O(3)-spin and lattice translational symmetry, but breaks lattice rotational symmetry and can be adiabatically connected to the one-dimensional Haldane phase; and the m = 1/2 phase that had been known to exist in the presence of an external magnetic field [20], but, in contrast to previous predictions, our iPEPS data showed it also exists in a small θ-window without external field. We concluded our analysis by investigating the nature of all phase transitions of the BBH model, mainly by looking at kinks in the energy and jumps in the magnetization per site as a function of θ, while also taking symmetry considerations into account. We found clear or subtle kinks in the energy for all transitions but the AFM to Haldane transition. Similarly, the magnetization displayed θ/π  slight to clear jumps except at the above-mentioned transition, leading us to conclude that all phase transitions are of first order, except possibly the AFM to Haldane transition, which we predict to be either of second or weak first order.
In addition to demonstrating that the spin-1 BBH model on the square lattice exhibits various exotic phases, such as several (partially) nematic phases (m = 1/2, FQ and AFQ3), threesublattice ordering (AFM3 and AFQ3), and even a highly symmetric paramagnetic phase that only breaks lattice rotational symmetry (Haldane), which are interesting in their own right from a theoretical point of view, we hope to have paved the way for a future investigation of the BBH model on the numerically more challenging triangular lattice, possibly offering insight into the not-yet understood behavior of NiGa 2 S 4 [12][13][14] and Ba 3 NiSb 2 O 9 [15][16][17]. Finally, our results show that iPEPS is a competitive method for analyzing strongly correlated spin systems, especially where quantum Monte Carlo suffers from the sign problem.

A Additional data
The appendix contains additional data to support the main text. We have included plots that provide extra evidence for the occurrence of a first order transition at θ = π/4. For the other transitions, the jump in magnetization or kink in the energy is clear enough to draw conclusions from. Additional plots can be provided upon request.
At θ = π/4, hysteresis can be observed around the phase transition (Fig. 15 left): the quadrupolar phase specifically can be simulated at θ = 0.249π < π/4 (where the ground state is in the magnetized phase) before jumping to a magnetized state at θ = 0.245π. Also, the eigenvalues of the Q-matrix (Fig. 15 right) differ in both phases: a rough extrapolation of D → ∞ shows that in the antiferroquadrupolar phase we have eigenvalues of approximately 0.15 (twice) and −0.3 once, whereas in the 120 • magnetically ordered phase we have 0.15, 0 and −0.15, indicating a jump in the spectrum, which leads us to conclude that the transition is of first order.
Additionally, from the spectrum of the Q-matrix we obtain some extra information about the quadrupolar phase: because Q has two identical eigenvalues, the nematic order is uniaxial (as opposed to biaxial).