Peeking into the $\theta$ vacuum

We propose a subvolume method to study the $\theta$ dependence of the free energy density of the four-dimensional SU($N$) Yang-Mills theory on the lattice. As an attempt, the method is first applied to SU(2) Yang-Mills theory at $T=1.2\,T_c$ to understand the systematics of the method. We then proceed to the calculation of the vacuum energy density and obtain the $\theta$ dependence qualitatively different from the high temperature case. The numerical results combined with the theoretical requirements provide the evidence for the spontaneous CP violation at $\theta = \pi$, which is in accordance with the large $N$ prediction and indicates that the similarity between 4d SU($N$) and 2d CP$^{N-1}$ theories does not hold for $N$=2.


I. INTRODUCTION
The θ parameter of the 4d Yang-Mills theory controls relative weights of different topological sectors in the path integral. Despite long history, it still remains as a challenging problem to identify the effect of the θ parameter on the non-perturbative dynamics of the theory.
For the special value θ = π the Lagrangian has CP symmetry, and we can ask whether or not the CP symmetry is spontaneously broken. In the large N limit [1] spontaneous CP violation at θ = π was demonstrated in Refs. [2][3][4]. For finite N a mixed anomaly between the CP symmetry and the Z N center symmetry shows that the CP symmetry in the confining phase has to be broken [5]. A similar conclusion was derived by studying restoration of the equivalence of local observables in SU(N ) and SU(N )/Z N gauge theories in the infinite volume limit [6]. (See also Refs. [7][8][9].) While these theoretical developments have narrowed down possible scenarios, an explicit nonperturbative calculation is necessary to unambiguously settle the fate of the CP symmetry at θ = π. Any direct numerical simulation at θ = π, however, has been difficult due to the notorious sign problem [37].
In Ref. [12] three of the authors of the present paper studied the vacuum energy density of the 4d SU(2) Yang-Mills theory around θ = 0 by lattice numerical simulations. The case of SU(2) gauge group is of particular interest since N = 2 is farthest away from the large N limit: there is a well-known parallel between 2d CP 1 and 4d SU(2) model [38], and the known vacuum at θ = π in the former [14] alludes to the appearance of gapless theory in the latter. By observing that the first two numerical coefficients in the θ expansion obey the large N scaling [39], it was inferred in Ref. [12] that the 4d SU(2) Yang-Mills theory at θ = π has spontaneous CP breaking, contrary to the naive expectation from the 2d CP 1 model.
In this work, we develop a subvolume method to explore the θ dependence of the free energy without any series expansion in θ and apply to 4d SU(2) theory. We find that it indeed works and show an evidence of spontaneous CP violation at θ = π at zero temperature, in consistency with the results of Ref. [12]. The subvolume method here is inspired by Ref. [22] and similar to that introduced in Ref. [23] to study 2d CP N−1 model.

II. SUBVOLUME METHOD AND LATTICE SET UP
In what follows, the subvolume method is described. After generating a number of gauge configurations at θ = 0 and implementing the n APE steps of APE smearing [24], the topological charge density q(x) is calculated with the five-loop improved topological charge operator [18] on the lattice. Defining the subvolume topological charge Q sub , we calculate the free energy density for the subvolume as [22,23], where U denotes the link variable and S g the gauge action. Note that the expectation value · · · is estimated on the θ = 0 configurations. The free energy density is then obtained as the infinite volume limit of f sub (θ): The goal of this study is to answer whether the spontaneous CP violation does occur at θ = π. Thus, the order parameter should be useful and is calculated through  f and df /dθ are evaluated separately and later used to make a consistency check. Some cautions are in order. First, the size of the subvolume V sub cannot be arbitrary. Obviously it has to be large enough to cover the typical correlation length of the system, which is considered to be around 1/(aT c ) = 9.50 [25] in the lattice unit. We thus restrict the size of the subvolume to V sub ≥ 10 4 . Suppose that V sub is large enough but V sub ≪ V full , where V full denotes the full volume. Then, the resulting f (θ) is expected to be independent of V full , and f sub would shows a scaling behavior as a function of V sub . As V sub grows and approaches V full , Q sub becomes close to an integer global Q, and f sub starts to converge a full volume result, in which we are not interested as it becomes clear later. The critical goal in this method is to identify the scaling behavior from which the infinite volume limit is extracted.
Secondly, the observables (2) and (6) contain the expectation value cos(θQ sub ) . Either observable becomes incalculable when cos(θQ sub ) fluctuates across zero. This is the sign problem in this method and sets the θ dependent upper limit on the size of V sub . Thus, another crucial point is whether the scaling behavior is realized before V sub reaches the upper limit. We will return to this point in the discussion section.
The free energy density of 4d SU(2) Yang-Mills theory is calculated at high and zero temperatures below. We employ the configurations generated in our previous work at T = 0 [12], while the ensemble corresponding to T = 1.2 T c is newly generated at the lattice coupling β = 1.975, the same as the T = 0 case. The simulation parameters are summarized in Tab. I.
The topological charge density measured at each smearing step is uniformly shifted as q(x) → q(x) + ǫ at every configurations so that the global topological charge Q = x∈V full q(x) takes the integer closest to the original value. The calculation of Q sub is carried out every five smearing steps. For the APE smearing, we take α = 0.6 in the notation of Ref. [27].
Topological observables on the lattice can be distorted by topological lumps originating from lattice artifacts. One can take away those by the smearing procedure, but at the same time the smearing may deform physical topological excitations, too. We studied this point in detail in Ref. [12], and developed the procedure to restore the physical information. The procedure consists of the extrapolation of the observables to the zero smearing limit by fitting over a suitable interval of the smearing steps. The fit range is fixed in advance by examining the response of the global topological charge to the smearing as done in Ref. [12]. The resulting fit ranges and the topological susceptibilities χ = Q 2 /V full thus determined using the global topological charge Q are shown in Tab. I for later use.
The θ-dependence is explored in the range of θ = k π/10 with k ∈ [1,20]. Each configuration is separated by ten Hybrid Monte Carlo (HMC) trajectories. In the following analysis, statistical errors are estimated by the single-elimination jack-knife method with the bin size of 500 and 100 configurations for zero and high temperatures, respectively. We mainly show the analysis of f (θ), but df (θ)/dθ is analyzed in parallel with similar quality.

III. TESTING THE METHOD AT HIGH TEMPERATURE
We first apply the subvolume method to the calculation of the free energy density above T c , where the instanton prediction, f (θ) ∼ χ(1 − cos θ) [28][29][30], is believed to be valid and numerically supported for SU(N ) with N ≥ 3 [31,32] as well. Using the translational invariance, a single-sized subvolume V sub /a 4 = l 3 × N T with l = 10, 12, · · · , 24 is taken from 64 places per a configuration, and the results are averaged.
The l dependence of f sub is shown in Fig. 1. It is found that in general the measured l dependence is not constant and the leading correction linear in 1/l is inevitable. The linear dependence is seen for θ = π/2 for l ≤ 22, whereas it ends around l = 20 for θ = π and it is hard to determine the linear region for θ = 3π/2. Thus, it turns out to be difficult to identify the scaling region unambiguously especially when θ is large, and hence we give up the precise determination. Instead, we choose three fit ranges in each extrapolation and try to estimate the potential size of the systematic uncertainty due to the ambiguity of the scaling region. Three fit ranges, l ∈ [12,16], [14,18] and [16,20], are examined when fitting to the expected scaling behavior where s(θ) denotes the surface tension of the nonzero θ domain and a the lattice spacing. All fits performed in this analysis yield χ 2 /dof < 3. It is interesting to see that the relative relation f sub (3π/2) > f sub (π) at small l flips toward the large l limit and f (θ) ends up with non-monotonic function. Since the data of f sub smoothly (but sometimes rapidly) connect the full volume (l 3 ×6 = 24 3 × 6) results from the above, the extrapolation fitting the data near the full volume tends to be smaller than that fitting the data far from the full volume. As a result, the discrepancy, i.e. the potential size of the systematic uncertainty, turns out to be larger at larger θ.
The results thus obtained are then extrapolated to n APE = 0 at each value of θ with the fit range shown in Tab. I. In the extrapolation, the linear fit goes well with χ 2 /dof < 3. The stability against small shifts of the fit range is seen in Fig. 2.
Finally, the free energy density obtained with three fit ranges are shown in Fig. 3 together with the full volume result (filled squares), where f (θ) is normalized by the topological susceptibility in Tab. I. The prediction from the dilute instanton gas approximation, 1 − cos(θ), is shown by the dashed curve. The function, θ 2 /2, is also shown as the solid curve for comparison. Taking into account the uncertainty arising from the ambiguity of the scaling region, the numerical results are consistent with the instanton prediction. Note that non-monotonic behavior of f (θ) seems robust at high temperature but is far from obvious before the extrapolations, as the surface tension term in Eq. (7) is monotonic.  Fig. 1, l ∈ [12, 16], [14,18] and [16,20], are shown in triangle-up, circle and triangledown, respectively. f (θ) can also be obtained from the numerical integration of df (θ)/dθ as shown by the dotted curves. The agreement with those curves supports that the two nontrivial extrapolations included in the whole analysis do not pick up accidental fluctuations and are stable.
The result with full volume is found to well agree with the instanton prediction. One may think that this is the simplest way to obtain f (θ). However, we will see that it does not work at T = 0. From the test, assuming that the instanton prediction is valid at high temperature, we learn that the scaling behavior of f sub would be linear and the region showing such a behavior starts around the dynamical length scale (∼ 1/(aT c )).

IV. APPLYING TO ZERO TEMPERATURE
Next we apply the subvolume method to calculate the vacuum energy density. This time the subvolume is defined by V sub /a 4 = l 4 with l = 10, 12, · · · , 24 and taken from 512 places per configuration. The l dependence of f sub (θ) is shown in Fig. 4 as before. Due to the sign problem in this method, some results at large θ and large l could not be calculated. But the available data show linear behavior. Following the previous analysis, three fit ranges of l ∈ [10,14], [12,16] and [14,18] are taken in the fit to (7) to estimate the systematic uncertainty. Contrary to the high temperature case, f (θ) turns out to be stable against the variation of the fit range, and does not show any sign of the flip, indicating monotonic behaviors of f (θ) as a function of θ.
The linear extrapolation to n APE = 0 is carried out with the fit range shown in Tab. I, and the fit is found to work well with χ 2 /dof < 3 as shown in Fig. 5. The stability against shift of the fit range is also confirmed.
Finally, the resulting f (θ) and df (θ)/dθ are shown in    (1 − cos θ). The stability of the two extrapolations during the analysis is confirmed as f (θ) and df (θ)/dθ well agree with the dotted curves. While the full volume calculation works only in the vicinity of θ = 0, the subvolume method succeeds to calculate, at least, to θ ∼ π. There are crucial differences from the high temperature case. First, the different choices of the fit range in l yield consistent results, and hence the potential systematic error from the ambiguity of the scaling region seems to be under control. Second, f (θ) is a monotonically increasing function, at least, to θ ∼ π, and the direct calculation of df (θ)/dθ clearly shows d f (θ)/dθ| θ=π = 0. Since d f (θ)/dθ = −i q(x) is CP odd, we conclude that CP is spontaneously broken at θ = π in the vacuum of the 4d SU(2) Yang-Mills theory [40] and that there is a phase transition to recover the CP symmetry at some finite temperature. In other  Fig. 4, l ∈ [10,14], [12,16] and [14,18], are shown in triangle-up, circle and triangledown, respectively. words, it is found that the 4d SU(2) Yang-Mills theory is in the large-N class unlike the 2d CP 1 model [41].
The subvolume method is equivalent to modifying the value of θ inside the subvolume. If the difference of θ is a multiple of 2π and the calculation respects the 2πperiodicity, the free energy would scale as the surface area of the subvolume when the subvolume is large enough. The lack of 2π periodicity in the free energy density should thus be interpreted as the presence of a metastable vacuum for a fixed value of θ (except for θ = π where two vacua interchanged by CP are degenerate and stable). Thus, we expect that the meta-stable vacuum should eventually decay into the stable one by the creation of a dynamical domain wall that attaches to the interface.
The absence of the decay of the domain into the domain-wall in the lattice calculation has an analog in the calculation of the static potential [42]. The static potential is calculated by inserting a Wilson loop, and should show the string breaking for configurations with light dynamical quarks when the two test charges are distant enough. But it does not occur, at least, within naive methods, and the resulting potential sticks to the original branch even after passing the transition point. The probable reason is that the overlap between the original state with two static charges and another lower energy state with two mesons are extremely small. We infer that the same happens to the calculation of f (θ) for θ > π and the first order phase transition is missing [43]. It is clearly interesting to directly see the formation of the domain wall on the lattice though it would not be straightforward.
We have mentioned the θ dependent upper limit on the size of V sub in sec. II. We examine the relation between the limit and θ |Q sub | at T = 0. Figure 7 shows θ |Q sub | as a function of 1/l, where we have used the approximate relation |Q sub | = (l 4 /V full ) 1/2 |Q| and the measured value of |Q| at T = 0 by ignoring the corrections to the relation of O(1/l) and O(1/(χV sub )). In the figure, the filled symbols represent the points where the calculations are failed due to the sign problem. It is seen that the upper limit indeed decreases with θ. Numerical investigation suggests that the upper limit on the subvolume scales as 1/θ.

VI. SUMMARY
We developed the subvolume method, which enables us to extract the θ dependence of the free energy den-sity in 4d Yang-Mills theory not restricted to θ ∼ 0. At high temperature, the method yields θ dependence consistent with the instanton prediction, as expected, within the large uncertainty due to the ambiguity of the scaling region. To fix this ambiguity, it is necessary to go to larger lattices. On the other hand, at zero temperature the sign problem arises instead, but still f (θ) could be calculated to θ ∼ π with the systematic uncertainty under control. Combining the numerical result with the theoretical requirement leads to the conclusion that the vacuum of 4d SU(2) Yang-Mills theory undergoes spontaneous CP violation at θ = π as large N theory does. Although the overlap problem prohibits the domain-wall from being formed, it is interesting to learn that such a object actually exists in the Yang-Mills theory [22].
We have tested the stability of the results by exchanging the order of extrapolations and obtained the consistent results with enlarged uncertainties.
In order to promote this study to the quantitative level, it is necessary to perform lattice simulations with larger volumes and finer lattice spacings. Further studies will be presented in the forthcoming paper [35]. It is a fascinating question to ask if our method is applicable for other questions with sign problems, such as gauge theories with finite values of the chemical potential.
While numerical results are not accurate past θ ∼ 3π/2, there are indications that the derivative df (θ)/dθ decreases past θ = π, and becomes smaller near θ ∼ 2π. This is consistent with the expectation [13] that there are two meta-stable branches of the SU(2) theory, each of which has 4π periodicity.