Phase structure of the CP(1) model in the presence of a topological $\theta$-term

We numerically study the phase structure of the CP(1) model in the presence of a topological $\theta$-term, a regime afflicted by the sign problem for conventional lattice Monte Carlo simulations. Using a bond-weighted tensor renormalization group method, we compute the free energy for inverse couplings ranging from $0\leq \beta \leq 1.1$ and find a CP-violating, first-order phase transition at $\theta=\pi$. In contrast to previous findings, our numerical results provide no evidence for a critical coupling $\beta_c<1.1$ above which a second-order phase transition emerges at $\theta=\pi$ and/or the first-order transition line bifurcates at $\theta\neq\pi$. If such a critical coupling exists, as suggested by Haldane's conjecture, our study indicates that is larger than $\beta_c>1.1$.


I. INTRODUCTION
The CP(N − 1) models in 1+1 dimensions share many properties with QCD in 3+1 dimensions, among them confinement, asymptotic freedom, instantons, a 1/N expansion, a topological charge, and a θ-term. Thus, they serve as benchmark models for developing and testing both new numerical techniques and new proposed solutions to open questions of QCD.
Both the CP(N − 1) models and QCD contain many nonperturbative phenomena, which cannot be addressed with conventional lattice techniques. One particular example is the above-mentioned θ-term, which is the origin of the strong CP problem in QCD [1,2]. The θ-term gives rise to an imaginary contribution to the action in the Euclidean formulation of lattice gauge theories. This complex action problem, also called the sign problem [3], prevents the successful application of Markov chain Monte Carlo (MCMC) methods for large values of the topological vacuum angle θ. To evade this problem, new numerical techniques have to be developed.
The phase structure of different CP(N − 1) models has been intensively investigated with different methods in the past [34,35,[46][47][48][49][50]. In particular, the phase diagram of the CP(1) model with a θ-term has been studied with a strong coupling analysis [48], Monte Carlo simulations [50], and TRG studies [34,35]. Here, it was shown that a first-order phase transition occurs at small values of β and θ = π, where the CP symmetry spontaneously breaks (see Fig. 1). For larger values of β, corresponding to the weak-coupling regime, the numerical studies in Refs. [34,35,50] indicated that the phase transition at θ = π becomes a second-order transition. Since the CP(1) model is equivalent to the O(3) model [51], these results are in agreement with Haldane's conjecture [52,53] that the O(3) model at weak coupling becomes gapless at θ = π, corresponding to a second-order phase transition. At the same time, it has also been suggested that the first-order transition line bifurcates at θ = π for large values β, thus exhibiting a non-trivial dependence on the coupling [34,35,48,50]. However, the exact value of the critical coupling, β c , remains unknown for both scenarios. While a strong coupling analysis suggests a bifurcation at β c = 0.56 [48], Monte Carlo simulations yield a second-order transition beyond β c = 0.5 [50], and the TRG analysis finds that the second-order phase transition occurs around β c ≈ 0.3 − 0.4 [34,35].
In the present work, we reanalyze the phase diagram of the lattice CP(1) model in the presence of a θ-term using high-precision TRG simulations. We provide a detailed analysis of the different errors contributing to the current and previous numerical studies. Interestingly, we find no indication of a bifurcation point or the onset of a secondorder phase transition up to β = 1.1, which is much larger than suggested by the previous studies. If the bifurcation point really exists, as suggested by Haldane's conjecture, our results indicate it occurs at β c > 1.1.
The paper is organized as follows. In Sec. II, we will describe the CP(1) model and the numerical TRG methods. In Sec. III, we will present our numerical results. In Sec. IV, we will summarize and discuss our results.  (1) model, where β = 1/g 2 is the inverse coupling and θ is the free parameter of the topological θ-term. The weak-coupling limit β → ∞ corresponds to the continuum limit of the lattice model. For small values of β, the model undergoes a first-order phase transition at θ = π (solid gray line), up to a certain critical coupling βc (yellow shaded region). For β > βc (blue region), different predictions exist and the phase structure is not entirely clear. There are hints toward a secondorder phase transition at θ = π (dotted brown line), in agreement with Haldane's conjecture, and toward a bifurcation of the critical first-order line, such that the model undergoes two first-order transitions in θ = π (dashed dark purple lines).

II. MODEL AND METHODS
In order to study the phase structure of the CP(1) model with the TRG [54][55][56], we use a Euclidean time lattice formulation. The discretized action of the model is given by where β = 1/g 2 is the inverse coupling constant and θ ∈ [0, 2π] is the free parameter of the topological θ-term. The two-component complex scalar fields z x ∈ C 2 reside on the sites x = (x 1 , x 2 ) of a two-dimensional lattice and are fixed to unit length, z * x z x = 1 ∀x. The gauge variables U x,µ reside on the links connecting the lattice sites x and x +μ along the direction µ (see Fig. 2). They are related to the components A x,µ of an auxiliary vector field via U x,µ = exp(iA x,µ ). The topological charge q x is defined as In order to explore the phase diagram of the CP(1) model using TRG, we need to obtain the tensor network representation of the partition function where we have used a shorthand notation for the integration measure over the complex scalar fields, To this end, we first derive an expansion for exp(−S θ=0 ) which factors into two parts: a part dependent on the complex scalar fields z x and a part depending on the vector fields A x,µ . Inserting this expansion into Eq. (3) and taking the integral over the scalar fields, we arrive at a tensor representation of the matter part. Subsequently, we derive a corresponding representation for the θ-term. The final tensors for Z θ can then be obtained by multiplying both tensors for the individual parts together and integrating over the gauge fields. To obtain a tensor network representation for θ = 0, we use the characterlike expansion of the Boltzmann factor [48] where d l,m are the expansion dimensions, h l,m are the expansion coefficients, f l,m are the expansion characters, and Z 0 (β) is a normalization factor. The non-negative integers l and m will become the indices of the tensor representation after integrating out the degrees of freedom z x and A µ . Note that due to the product in front of the right-hand side of Eq.
As the explicit form of the dimensions d l,m and the characters f l,m is rather complicated, we refer to Refs. [48,56] for details and only quote their normalization conditions here, In order to construct the tensor representation, we introduce an additional decomposition of f l,m [see Fig. 3 where {a} = a 1 , ..., a l , a 1 , ..., a m , and the indices a i (a i ) range over 1, 2 (see Refs. [48,56] for details).
With this expansion, we can divide , such that we can subsequently integrate out the complex scalar fields z x site by site. Note that the θ-term does not depend on z x , which implies that only the first term in Eq. (1) is relevant for integrating out z x . Finally, we define the tensor [see Fig. 3 The tensor W x represents a rank-4 tensor, whose bonds are given by the four multi indices (l s , m s , {a}), , which allows for representing exp(−S θ=0 ) as a tensor network up to the factor e i(m−l)Ax,µ depending on the gauge field. Next, we derive a tensor representation of the topological θ-term starting again from its character expansion, which reads In the expression above, n p is an integer, and we have defined In order to take the integral over the gauge field A x,µ , we consider the following integrals of A x,1 and A x−2,2 , Similarly, the integral over A x,2 yields δ lt−mt up−vp . Using δ ab δ bc = δ ac , we finally get the following tensor representation, The partition function of the system can now be represented as the contraction of the rank-4 tensors T x stuv sitting at the different lattice sites. Compared to W x , the multi-index for each bond of the tensor T x comprises an additional index resulting from the character expansion of the topological term [see also Fig. 3(c)].
Since the indices l, m, and n p in the character expansions in Eq. (1) and Eq. (11) range over a countable infinite set, the bonds of the tensor T x are infinite dimensional. Hence, for numerical calculations we need to truncate the bonds to a finite dimension. To this end, we limit the number of terms that we keep in the character expansions. For the indices l, m we define a cutoff k max and restrict ourselves to values fulfilling l + m ≤ k max . This results in a bond size of χ β = 1 + k max 2 kmax+1 for the tensor representing the partition function without the θ-term. Similarly, we only keep a subset of χ θ terms in the expansion for the topological term in Eq. (11). Since our numerical analysis primarily focuses on the region θ = π, we choose the terms with the largest absolute value of G np (π). For example, for χ θ = 2, we take into account n p = 0, −1 because the absolute value of |G 0,−1 (π)| = 2/π is larger than all other values of |G np (π)|. The size of the truncated tensor T x is thus given by χ β × χ θ .
These cutoffs in the tensor T x introduce systematic errors. Since h l,m is given by a ratio of the modified Bessel functions (see Eq. (6)), the cutoff at k max truncates only an exponentially small contribution to the tensor T x stuv . The truncation in the expansion for the θ-term leads to an error of the order O(1/χ θ ), as shown in Eq. (12). We will estimate these systematic errors by using different k max and χ θ our numerical calculations.
In addition, the TRG algorithm introduces another systematic error. Starting from the initial, finite dimensional tensor T x , the coarse-gaining step of the algorithm would produce a tensor which has the same structure as the original one, but with a bond size corresponding to the square of the one of the original tensor. Hence, for the computation to be sustainable, one has to truncate the tensor during the coarse-graining step and limit the bond size to a maximum value D. We will estimate this additional systematic error by performing calculations for various values of the bond dimensions D for a fixed set of parameters. In particular, D should be chosen larger than the size of the initial tensor χ β × χ θ to capture the relevant physics of the model.

III. RESULTS
In order to explore the phase structure of the model, we compute the partition function using TRG methods, which in turn allows us to obtain the free energy density where V is the dimensionless volume of the lattice. Moreover, we can also compute the specific heat from the partition function using the relation In particular, studying the free energy will allow us to determine the phase structure of the model and to detect a possible phase transition as discontinuities in the derivatives of the free energy. As a first step, we focus on the case θ = 0 to explore systematic errors due to the truncation in k max over a large range β. This allows us to identify the range of couplings in which our TRG computations yield reliable results. As a second step, we turn to the model in the presence of a topological θ-term. Using again the TRG approach, we study the model in a region around θ = π, where we carefully estimate our systematic errors due to the finite values of (D, k max , χ θ , V ).

A. CP(1) model without θ-term
To begin with, we benchmark the TRG method for θ = 0 and explore the effect of systematic errors due to the truncation of the expansion in Eq. (5). To this end, we calculate the free energy and specific heat of the CP(1) model using the anisotropic TRG (ATRG) method with randomized singular value decomposition (SVD) [57,58]. In order to avoid errors due to the randomized SVD, we set the oversampling parameter to the value of D, meaning that we compute 2D singular values at each coarse-graining step and select the D dominant ones from those. This ATRG method is faster than the usual TRG algorithm, while the precision is comparable to the one of TRG for same D. Since we choose the oversampling parameter large enough, the systematic error from the randomized SVD is negligible. Figure 4 shows our results for the free energy and the specific heat for D = 80, k max = 2, 4, and a volume of V = 2 40 . Focusing on the free energy, we see that the systematic error due to the finite value of k max is negligible in the range of small values of the inverse coupling. In particular, for β ≤ 1.5 there is essentially no difference between results with k max = 2 and 4. Only as we enter the region of larger β, the curves start to deviate slightly while the qualitative behavior of the free energy is still the same for both values of k max .
Looking at the specific heat, in Fig. 4, we see that truncation effects due to finite k max are more pronounced. This is not too surprising, as the specific heat is obtained via Eq. (16) by approximating the second derivative with finite differences, where we choose ∆β = 0.5 for the data shown in the inset of Fig. 4. Thus, the systematic errors in log Z θ are enhanced by a factor of (∆β) 2 leading to significantly larger deviations in the regime of larger β. Focusing on smaller values of the inverse coupling, our data show that truncation effects are nevertheless small despite this effect, and we can reliably determine the specific heat with k max = 2 in regimes β ≤ 1.1. Hence, for all the following we will solely focus on that region and use k max ≥ 2, for which systematic errors due to the truncation of the expansion in Eq. (5) are very small. We note that previous work using TRG [35] observed large fluctuations in the specific heat in the regime of large β. In contrast, our ATRG results are stable even for large values of β. While Ref. [35] used a bond dimension D = 21, our calculations rely on a larger bond dimension D = 80. Thus, we suspect that these fluctuations were caused by truncation effects due to the small bond dimension.

B. CP(1) model with θ-term
In order to study the CP(1) model in the presence of a θ-term, we use the bond-weighted TRG method [59] in conjunction with regular SVD to truncate the tensor at each coarse-graining step. Compared to the original TRG method, bond-weighted TRG introduces bond weights on the edges of the tensor network, allowing for more precise results within the same calculation time. Figure 5 shows our results for the free energy as a function of θ for the volume V = 2 24 , bond dimension D = 80, k max = 2, a truncation of the character expansion of the topological term χ θ = 2, and three values of the inverse coupling β = 0.1, 0.6, 1.1. Focusing first on β = 0.1 in Fig. 5(a), we observe a clear cusp at θ = π. Hence, the first derivative of F (β = 0.1, θ) with respect to θ is discontinuous at θ = π, indicating that the model undergoes a first-order phase transition at this point. This observation is consistent with previous findings of a first-order transition for β = 0.1 [34,35,48,50].
For larger inverse couplings of β = 0.6 and β = 1.1, the absolute values of the free energy change, as Fig. 5(b) and Fig. 5(c) reveal, respectively. However, for both β = 0.6 and β = 1.1, we still find the same cusp structure as for the smaller inverse coupling of β = 0.1. Thus, even for these large values of β, the transition is still of first order. These results disagree with previous findings [34,35,48,50] that for β c 0.4 − 0.56, the transition at θ = π becomes of second order and the first-order transition line bifurcates at θ = π.
In order to ensure that the discrepancies for the phase structure obtained from our numerical results and previous studies are not caused by systematic errors due to the choice of k max , D, and χ θ , we carefully examine the effect of these parameters. Figure 6 shows the difference in the results for the free energy, ∆F (β, θ) = |F (β, θ; D, k max , χ θ ) − F (β, θ; 80, 2, 2)| , between (D, k max , χ θ ) = (80, 2, 2) and various larger choices of k max , D, and χ θ for a fixed volume of V = 2 24 . Note that changing χ θ and k max affects the size of the initial tensor as discussed above. While for (k max , χ θ ) = (2, 2) the bond size of the initial tensor is 34, for (2,4) and (3,2) it grows to 68 and 98 respectively. Since D should be chosen well above the initial tensor size, we do not only increase k max or χ θ , but simultaneously also the value of D. Figure 6(a) shows our results for the smallest value of the inverse coupling corresponding to β = 0.1. In particular, we observe that increasing the bond dimension while keeping (k max , χ θ ) = (2, 2) has an extremely small effect on the results, with an absolute change in the range below 10 −14 . Similarly, increasing k max does not lead to noticeable deviations and the absolute change is in the same range. Increasing the value of χ θ to 4 yields a slightly larger deviation, which is still very small and in  For the larger inverse couplings of β = 0.6 and β = 1.1 in Fig. 6(b) and Fig. 6(c), respectively, we observe a qualitatively similar picture. While the absolute deviation is larger than for the case of β = 0.1, these effects are well below the percentile range even for our largest value of β. Just as before, increasing D and k max has the least impact, whereas increasing χ θ to 4 now clearly has a slightly larger effect on the results. In relation to the values of the free energy shown in Fig. 5(b) and Fig. 5(c), which are on the order of 10 −1 and 3, respectively, the systematic effects due to finite (D, k max , χ θ ) are again negligible.
To further study the possible origin of the discrepancies between our and previous numerical results, we now investigate the cusp in the free energy as a function of the volume V . To this end, we plot a close-up around the region of θ = π for various volumes and various values of the inverse coupling in Fig. 7. Looking at our results for β = 0.1 [see Fig. 7(a)], we observe that for volumes below 2 14 , the free energy is smooth at θ = π. Thus, the first derivative does not show any discontinuity for small volumes. Only when increasing V to larger values, the cusp at θ = π eventually emerges, signaling the onset of the first-order transition directly in the free energy. For larger values of the inverse coupling, corresponding to β = 0.6 and 1.1, we see a qualitatively similar behavior, as Figs. 7(b) and 7(c) reveal. Interestingly, the finite-volume effects seem to be slightly less pronounced for larger values of β, as a comparison between Figs. 7(a) and 7(c) shows. For β = 1.1, we observe that the data for V = 2 14 already has a clear signature of a cusp, and increasing the volume only leads to a slightly sharper peak.
In summary, the analysis of the systematic errors from finite (V, D, k max , χ θ ) further supports the reliability of our findings throughout the entire range of β that we study. We conclude that the phase transition at θ = π is still of first order up to an inverse coupling of β = 1.1.

IV. DISCUSSION AND CONCLUSIONS
In this paper, we numerically studied the phase diagram of the lattice CP(1) model with a θ-term. We observe a CP-violating, first-order phase transition at θ = π and small β, in agreement with previous findings [34,35,48,50]. However, our results do not confirm previous indications for a bifurcation in the first-order transition line at some critical value of β, ranging from β c = 0.56 [48] to β c = 0.5 [50] to β c = 0.4 [34,35]. Instead, we find that the first-order transition at θ = π persists up to β = 1.1 without any bifurcation, and we observe no indication that the phase transition at θ = π becomes of second order for β > β c . The existence of a second-order transition at θ = π and β > β c was suggested by the previous studies based on a strong coupling analysis [48], Monte Carlo simulations [50], and first TRG studies [34,35]. Compared to the strong coupling analysis of the CP(1) model in Ref. [48], our method does not require an expansion in the inverse coupling β = 1/g 2 , but instead performs truncations and character expansions to construct the tensor representation. Thus, while for small β the results in [48] may be more accurate than ours, the strong coupling analysis becomes more challenging for larger values of β. Since our TRG results clearly show that the effects due to the truncations of the character expansions are very small throughout the entire range of β we study, our results suggest that the strong coupling expansion does not correctly capture the physics of the model up to β = 1.1.
Similarly, the Monte Carlo studies in Ref. [50] suffer from the sign problem, which makes it difficult to reliably examine the phase diagram of the CP(1) model for large values of θ. In particular, the Monte Carlo approach becomes more challenging when approaching the phase transition at θ = π, whereas our TRG calculation only shows small systematic errors, which are essentially independent of the value of θ for the parameter range we study.
Compared to previous TRG-based studies of the model, our approach to determine the phase structure is more direct. In particular, Refs. [34,35] did not study the free energy of the model at large volumes, but rather explored the free energy at small volumes and used finite-size scaling behavior of the topological susceptibility as an indicator to determine the order of the phase transition. In contrast, our computations are performed at a much larger volume, for which the results are well converged, which allows us to study the behavior of the free energy without having to rely on finite-size scaling. Indeed, as we have shown in Fig. 7, a clear signature of the first-order transition only occurs at volumes of V > 2 14 . We note that our bond dimension of D = 80 is only slightly larger than the bond dimension D = 68 used in Ref. [35], but increasing the bond dimension only has a minor effect on our results, see Fig. 6.
To summarize, our results demonstrate that the phase transition at θ = π is of first order up to β ≤ 1.1, and no second-order transition occurs. The discovery of such a second-order transition would be crucial for confirming Haldane's conjecture [52,53] that the O(3) model becomes gapless at θ = π and weak coupling. As our study reveals, the possibility of a second-order transition at θ = π and a bifurcation of the first-order transition line at θ = π is only given for β c > 1.1. Compared to previous studies, our method does not suffer from the sign problem as Monte Carlo simulations, and our results have small systematic errors of 10 −3 due to larger volumes than in earlier TRG studies. Thus, our study reveals that testing Haldane's conjecture in the weak-coupling regime requires both major numerical efforts and detailed systematic error analysis. In the future, we will extend our study to larger values of β > 1.1. While this task is more challenging due to enhanced truncation effects in the weak-coupling regime, our numerical make us positive that this parameter range can be reliably accessed with TRG methods. In the main text, we focused on the shape of the free energy and used the emergence of a cusp toward the thermodynamic limit as an indication for a first-order phase transition. Alternatively, following Refs. [60,61], we can also examine the finite-size scaling of the peak value χ peak of the susceptibility at θ = π. For a first-order phase transition, it is expected that χ peak is proportional to the volume, whereas for a second-order transition, one expects a scaling of χ peak ∝ V γ with γ < 1 [62]. In principle, the susceptibility could be directly determined from our TRG results for Z θ , by computing the derivative numerically. However, this would require a very fine resolution in order to avoid large errors. Thus, we follow Ref. [61] and approximate our data for the free energy close to θ = π with a polynomial. Since finite-size effects round off the peak in the free energy, we choose the functional form where we consider only even powers to ensure the symmetry around θ = π. An example of this interpolation is shown in Fig. 8. The peak value of the susceptibility can then be obtained from the polynomial interpolation using Eq. (A1), and is given by χ peak = −2βc 1 . In order to determine the volume dependence of χ peak , we repeat the analysis described above for various volumes. The values for the volume are chosen such that they are large enough to avoid excessive finite-size effects and small enough to avoid a cusp in the free energy, as we observed for large volumes in the main text. Our results for this analysis for various values of the coupling are shown in Fig. 9. Focusing on β = 0.6 first [see Fig. 9(a)], we clearly observe that χ peak is proportional to the volume and γ is in good agreement with 1. Increasing the inverse coupling to β = 0.8, we see a very similar picture, as the data in Fig. 9(b) reveals. Again, the peak value of the susceptibility scales linearly with the volume, and the value for γ obtained from our fit is compatible with 1 within error bars. For β = 1.1, we observe that it becomes increasingly harder to find a window where finite-size effects are not excessive and the free energy is still far enough away from showing a cusp. As a result, the data points for small values of the volume in Fig. 9(c) show slightly larger error bars than for the previous values of β. When determining the volume dependence of χ peak for this case, we again obtain a value of γ that agrees with 1 within error bars. In summary, the finite-size scaling of the peak value χ peak of the susceptibility indicates that the phase transition is of first order throughout the entire range of β under consideration. This is consistent with our results in the main text, as obtained from the direct examination of the free energy in the limit of large volumes.