The strong CP problem solved by itself due to long-distance vacuum effects

The vacuum of quantum chromodynamics has an incredibly rich structure at the nonperturbative level, which is intimately connected with the topology of gauge fields, and put to a test by the strong CP problem. We investigate the long-distance properties of the theory in the presence of the topological $\theta$ term. This is done on the lattice, using the gradient flow to isolate the long-distance modes in the functional integral measure and tracing it over successive length scales. The key point is that the vacuum splits into disconnected topological sectors with markedly different physical characteristics, which gives rise to a nontrivial behavior depending on $\theta$. We find that the color fields produced by quarks and gluons are screened, and confinement is lost, for bare vacuum angles $|\theta|>0$, thus providing a natural solution of the strong CP problem. The renormalized vacuum angle $\theta$ is found to flow to zero in the infrared limit, leading to a self-consistent solution within QCD.


Introduction
Quantum chromodynamics (QCD) describes the strong interactions remarkably well, from the smallest distances probed so far to hadronic length scales where quarks and gluons confine to hadrons. Yet it faces a problem. The theory allows for a CP-violating term S θ in the action, S " S 0`S θ . In Euclidean space-time it reads where Q is the topological charge, and θ is the vacuum angle. It is an arbitrary phase with valueś π ă θ ď π. A nonvanishing value of θ would result in an electric dipole moment d n of the neutron. The current experimental upper limit is |d n | ă 1.8ˆ10´1 3 e fm [1], which suggests that θ is anomalously small. This feature is referred to as the strong CP problem, which is considered as one of the major unsolved problems in the elementary particles field.
It is known from the case of the massive Schwinger model [2] that a θ term may change the phase of the system. Callan, Dashen and Gross [3] have claimed that a similar phenomenon will occur in QCD. The statement is that the color fields produced by quarks and gluons will be screened by instantons for |θ| ą 0. 't Hooft [4] has argued that the relevant degrees of freedom responsible for confinement are color-magnetic monopoles, realized by partial gauge fixing [5], which leaves the maximal abelian subgroup U(1)ˆU(1) Ă SU(3) unbroken. Quarks and gluons have color-electric charges with respect to the U(1) subgroups. Confinement occurs when the monopoles condense in the vacuum, by analogy to superconductivity. This has first been verified on the lattice by Kronfeld et al. [6], and tested in [7]. In the presence of a θ term the monopoles acquire a fractional color-electric charge proportional to θ, as has been noticed by Witten [8]. It is then expected that the color fields of quarks and gluons will be screened by forming bound states with the monopoles, driving the theory into a Higgs or Coulomb phase [4,9] for |θ| ą 0. A similar mechanism can be expected for instanton and center vortex [10] based models of the vacuum, which both can be connected to abelian monopoles [11,12].
This strongly suggests that θ is restricted to zero in the confining phase of the theory, which would mean that the strong CP problem is solved by itself. A direct derivation of the phase structure of QCD for nonvanishing values of θ from first principles has remained elusive. A crucial step in solving the strong CP problem is to isolate the relevant degrees of freedom. This is a multi-scale problem, which involves the passage from the short-distance weakly coupled regime, the lattice, to the long-distance strongly coupled confinement regime. The framework for dealing with physical problems involving different energy scales is the multi-scale renormalization group (RG) flow. Exact RG transformations are very difficult to implement numerically. The gradient flow [13,14] provides a powerful alternative for scale setting, with no need for costly ensemble matching. It can be considered as a particular, infinitesimal realization of the coarse-graining step of momentum space RG transformations [15,16,17,18]à la Wilson [19], Polchinski [20] and Wetterich [21], leaving the long-distance physics unchanged, and as such can be used to study RG transformations directly.
In this work we investigate the long-distance properties of the theory in the presence of the θ term (1) using the gradient flow, and show that CP is naturally conserved in the confining phase.

The gradient flow
The gradient flow describes the evolution of fields and physical quantities as a function of flow time t. The flow of SU(3) gauge fields is defined by [14] B t B µ pt, xq " D ν G µν pt, xq , G µν " B µ B ν´Bν B µ`r B µ , B ν s , D µ¨" B µ¨`r B µ ,¨s , (2) where B µ pt, xq " B a µ pt, xq T a , and B µ pt " 0, xq " A µ pxq is the original gauge field of QCD. It thus defines a sequence of gauge fields parameterized by t. The renormalization scale µ is set by the flow time, µ " 1{ ? 8t for t " 0, where ? 8t is the 'smoothing range' over which the gauge field is averaged. Correlators of the flowed fields are automatically renormalized [22]. The expectation value of the energy density defines a renormalized coupling [14] g 2 GF pµq " at flow time t in the gradient flow (GF) scheme. Varying µ, the coupling satisfies standard RG equations. The dependence on µ encodes the dynamical properties of the theory, from confinement in the infrared (IR) to asymptotic freedom at short distances. For a start we may restrict our investigations to the SU(3) Yang-Mills theory. If the strong CP problem is resolved in the Yang-Mills theory, then it is expected that it is also resolved in QCD. Quarks play no role in the IR if they are massive. We use the plaquette action [23] Re Tr U µν pxq¯ (5) to generate representative ensembles of fundamental gauge fields. For any such gauge field the flow equation (2) is integrated to the requested flow time t. We use a continuum-like version of the energy density Ept, xq obtained from a symmetric (clover-like) definition of the field strength tensor G µν pt, xq [14]. The simulations are done for β " 6{g 2 " 6.0 on 16 4 , 24 4 and 32 4 lattices.
The lattice spacing at this value of β is a " 0.082p2q fm, where we have taken ? t 0 " 0.147p4q fm to set the scale [24,25], with t 0 defined by t 2 0 xEpt 0 qy " 0.3. Our current ensembles include 4000 configurations on the 16 4 lattice and 5000 configurations on the 24 4 and 32 4 lattices each. First results on the two smaller lattices have been reported in [26]. In this work we will consider bulk quantities only, like the energy density Eptq, for example. Limits on the flow time are set by the lattice volume. For bulk quantities we expect finite size effects to become noticeable at ? 8t Á L, where (dimensionful) L is the linear extent of the lattice. To check this we plot the dimensionless quantity t 2 xEptqy on the 16 4 , 24 4 and 32 4 lattice as a function of t{a 2 in Fig. 1. The data on the 32 4 lattice fall on a straight line up to t{a 2 « 150, corresponding to ? 8t « 32. On the 24 4 and 16 4 lattices finite size effects show up at t{a 2 « 80 and 50, respectively, also in rough agreement with ? 8t « L.
Physical observables should be independent of the RG scale. In the Yang-Mills theory at zero temperature the choice of bulk quantities is limited. Two such observables are the topological susceptibility and the normalized Polyakov loop susceptibility. The topological susceptibility is defined by where V " L 4 . We define the normalized Polyakov loop susceptibility [27] by where with U 0 px 0 , xq being the (complex valued) link matrix in time direction and V 3 the spatial volume, V 3 " L 3 . Note that x|P| 2 y denotes the Polyakov loop correlator The Polyakov loop requires normalization to be interpreted as the free energy of static quarks. We use the field-theoretic definition (1) of Q. On the 24 4 lattice at flow time t{a 2 " 50, for example, the mean deviation from the nearest integer, ∆Q, turns out to be |∆Q{Q| « 0.03 for |Q| ą 0. In Fig. 2 we show χ t and χ P as a function of flow time for t{a 2 " 10 to 100. Both quantities are independent of t within the errorbars, as expected. A linear fit to the topological susceptibility gives ? t 0 χ 1{4 t " 0.162p3q. This number agrees precisely with the result of a high statistics calculation reported in [28], ? t 0 χ gives χ P " 0.289p7q. This number compares well with the result of a 2D Gaussian distribution of (real and imaginary) P, which gives χ P " 4{π´1 " 0.273. The flow-time independence of the topological susceptibility χ t is supported by two facts. Namely, that the vacuum splits into disconnected topological sectors, and that the field-theoretic definition (1) of Q is independent of the flow time. The former will be discussed in detail in Sec. 4. The latter follows from the fact that the charge density q " p1{32π 2 q G a µνG a µν can be written as a total derivative, q " B µ ω µ , where ω µ is the Chern-Simons density or 0-cochain [29]. The latter is gauge variant. However, its flow-time derivative, B t ω µ " p1{8π 2 q B t B a νG a µν , is gauge invariant, which leads to B t Q " 0 on the periodic lattice after partial integration.

Running coupling and linear confinement
The gradient flow running coupling α GF pµq " g 2 GF pµq{4π, introduced in (4), plays a key role in our investigations. In Fig. 3 we show α GF {π on the 32 4 lattice as a function of t{a 2 " 1{8 a 2 µ 2 . The striking result is that α GF is a linear function of t " 1{8µ 2 in the nonperturbative regime, starting well below α GF " 1 (t 2 xEptqy " 0.3). We have already seen that the linear behavior extends to ? 8t « L, Fig. 1, corresponding to µ À 100 MeV, and it may be assumed that it will extend to even larger values of t as the volume is increased. We will come back to this point later. The result is the gradient flow beta function In terms of µ, the RG equation (10) has the implicit solution  which leads to α GF pµq " The understanding is that at shorter distances the gradient flow beta function, as well as the running coupling, connect analytically with the perturbative expressions. A fit of (12) to the lattice data shown in Fig. 3 gives ? t 0 Λ GF " 0.475p16q. In any other renormalization scheme the beta function β S is given implicitly by the relation Solving (13) for β S pα S q and α S pµq in the nonperturbative regime, using (11) and (12), gives This states that the linear rise of the running coupling, Fig. 3, is universal, up to a scheme dependent constant factor. To make contact with phenomenology, we need to transform the gradient flow coupling α GF to a common scheme. A preferred scheme in the Yang-Mills theory is the V scheme [30]. In this scheme From the literature we know Λ V {Λ MS " 1.600 [30] and Λ MS {Λ GF " 0.534 [14]. This leads to ? t 0 Λ V " 0.406p14q and ? t 0 Λ MS " 0.217p7q. The latter number is in excellent agreement with the outcome of a recent dedicated lattice calculation [31], ? t 0 Λ MS " 0.220p3q. The linear rise of α V pµq with 1{µ 2 , which is commonly dubbed infrared slavery, effectively describes many low-energy phenomena of the theory. So, for example, the static quarkantiquark potential, which can be described by the exchange of a single dressed gluon, Vpqq " 4 3 α V pqq{q 2 . A popular example is the Richardson potential [32], which reproduces the spectroscopy of heavy quark systems, like charmonium and bottomonium, very well. Upon performing the Fourier transformation of Vpqq to configuration space, we obtain where σ, the string tension, is given by which is exactly what one expects from Regge phenomenology.
The consistently good agreement of our results so far with phenomenology provides an invaluable test of the use and potential of the gradient flow.
It is interesting to compare the beta function (10), matched with the perturbative expression at shorter distances, with its perturbative counterparts. In Fig. 4 we show the gradient flow beta function, together with the four-loop [33] and twenty-loop [34] perturbative results in the qq scheme (Λ qq {Λ V " 0.655). As was to be expected, the perturbative beta function gradually approaches the nonperturbative expression with increasing order.

Yang-Mills theory with θ term
We will show now that with increasing flow time the ensemble of gauge fields splits into quantum mechanically disconnected sectors of topological charge Q. A similar behavior has been found long time ago by 'cooling' lattice gauge field configurations [35,36]. This is expected to happen at ever smaller flow times as the lattice spacing is reduced, until the original lattice gauge field ensemble itself splits into isolated topological sectors [14]. Furthermore, we shall show that key observables depend sensitively on the charge Q of the respective sector. When Fourier transformed to the θ vacuum, this suggests a nontrivial phase structure of the theory.
We distinguish the topological sectors by the affix Q. A prominent example is the energy density xEpQ, tqy. In Fig. 5 we plot the average 'action' VxEpQ, tqy{8π 2 , normalized to one for a single classical instanton, on the 32 4 lattice as a function of t{a 2 . While VxEpQ, tqy{8π 2 develops a plateau for borderline charges at large flow time, the ensemble average, that is the statistical  average across all topological sectors, vanishes like 1{t. Unlike 'cooling', the topological sectors are remarkably stable, even for borderline charges. The probability distribution function PpQq for topological charge Q is given by PpQq " ş Q DU µ expt´S 0 u{ ş DU µ expt´S 0 u. From B t Q " 0 follows that PpQq is independent of the flow time t once the ensemble has settled into disconnected topological sectors. Thus, it is not to be feared that the gradient flow will end in a semi-classical ensemble of noninteracting instantons.

Running coupling
If the general expectation is correct and the color fields are screened for |θ| ą 0, we should find, in the first place, that the color charge is screened. This will show in the θ-dependence of the running coupling. From xEpQ, tqy we derive α V pQ, µq depending on Q. In Fig. 6 we plot α V pQ, µq divided by the ensemble average α V pµq as a function of t{a 2 and |Q|. Already at relatively small flow times α V pQ, µq begins to fan out according to Q. If multiplied by the lattice volume, the data in all three figures fall on top of each other, apart from finite volume effects. The transformation of α V pQ, µq to the θ vacuum is achieved by the discrete Fourier transform weighted by the topological charge distribution PpQq. On the finite lattice Zpθq is well defined. It is independent of t, like PpQq. We will return to Zpθq and the resulting free energy in Sec. 5. In (17) the parameter θ is the bare vacuum angle, that labels superselection sectors of the theory [37]. It is the parameter that appears in the lattice action and, hence, determines the topological properties of the vacuum. Limits are set by the precision of PpQq and α V pQ, µq. The statistics needed for consistent precision will increase relatively moderately with the root of the volume. The charge distribution PpQq is determined by the real part of the action, S 0 , which increases with |Q| and, thus, suppresses configurations with a large number of (anti-)instantons.
t/a 2  That makes it increasingly difficult to determine PpQq precisely for large values of |Q|. This circumstance is completely independent of whether we simulate at θ " 0 or any other value |θ| ą 0, which is to say that the situation would not improve if we could simulate the complex action. Fortunately, for our main conclusions we will need to know the Fourier sum for small values of |θ| only, which is rather insensitive to fluctuations at large values of |Q|. By comparing results on different volumes we can control statistical and systematic effects.
We plot α V pθ, µq for our three volumes in Fig. 7. The two entries on the left show finite size effects for t{a 2 Á 50 and t{a 2 Á 80, respectively, as expected (cf. Fig. 1). The shape of the curves does not change though. The calculations of α V pθ, µq are rather robust. The determining factor is that α V pQ, µq is a monotonically increasing function of |Q|, which makes the numerator of (17) fall off much faster than the denominator, Zpθq, while the charge distribution PpQq largely cancels out. An interesting feature of α V pθ, µq is that for fixed values of |θ| ą 0 it rises to a certain point before it drops towards zero as t " 1{8µ 2 is increased. A similar behavior is found for the running coupling α s pT, µq at finite temperature T ą T c [38]. With our current statistics we are not able to compute α V pθ, µq with confidence for t{a 2 À 10 and 20 on the 24 4 and 32 4 lattice, respectively. But there is no doubt that it will continue to flatten, as seen on the 16 4 lattice. The main hindrance is that fluctuations of α V pQ, µq with respect to Q can be severe on large volumes before the ensemble has settled into truly disconnected topological sectors, as can be seen in Fig. 6. The situation is expected to improve for smaller lattice spacings a. Figure 7 shows that for any fixed value |θ| ą 0 the running coupling α V pθ, µq is bounded by a finite value, αV " α V pθ, µ˚q, µ˚ą 0, however small θ is, which thwarts linear confinement, better called 'separation-of-charge confinement' [39]. When viewed from a fixed value of flow time, we find that α V pθ, µq drops to zero as |θ| is increased. This happens the faster, the larger t is. The scale parameter µ corresponds to the distance r " expt´γ E u{ ? 2 µ « 0.40{µ at which the color charge is probed, which derives from converting the potential Vpqq to Vprq in coordinate space [40]. The curves in Fig. 7 thus acquire a concrete physical meaning. For example, at |θ| " 0.8 the color charge is totally screened at distances r Á 0.6 [fm], while at |θ| " 0.4 it is screened for r Á 0. 8 [fm]. This means that for any fixed value |θ| ą 0 quarks and gluons can be separated, perhaps with increasing cost of energy, by being neutralized by collective gluonic excitations.
Analytically, α V pθ, µq in Fig. 7 can be expressed by [26] within the errors. On the 24 4 lattice D « 0.10 and λ « 0.75. The charge is totally screened when the right-hand side vanishes. We define the screening radius, λ S , at which the running coupling has dropped by a factor 1{e. From (18) we find λ S « 0.31{|θ| [fm]. The situation here is very similar to the finite temperature phase transition. While confinement is lost for T ą T c , the screening length [41] is consistently described by λ S 9 1{pT´T c q.
It is illustrative to compare our results to the predictions of the dual superconductor model of confinement and its extension to nonvanishing values of θ [4], which offers a dynamical explanation of the screening mechanism, as described in the Introduction. In the θ vacuum the monopoles acquire a color-electric charge q " θ{2π. We may assume that the monopole density ρ [42,43,44] does not change significantly for small values of θ. The Debye screening length of a charged particle in the (super)conducting vacuum is given by λ D " a pE F {ρ q 2 q, where E F is the Fermi energy. This leads to λ D " a p4π 2 E F {ρq{|θ|, which matches our result. From Fig. 7 it is evident that within QCD, θ can only assume values within the envelope of the curves, to maintain confinement, for example. This indicates that θ needs to be renormalized in the IR, along with α V . From (18) we derive RG equations, which govern the flow of both α V and θ as a function of the scale parameter µ. The result is a set of partial differential equations. For α V θ 4 ! 1 the equations simplify to which applies to the major part of Fig. 8. Here θpµq denotes the renormalized θ parameter, which can be thought of as the coupling that enters the effective Lagrangian averaged over virtualities larger than µ. In Fig. 8 we show the RG flow of α V pθ, µq and θpµq in the pθ, π{α V q plane for different initial values of θ, with t increasing from top to bottom. Assuming that there is no phase transition down to the IR limit, this leads us to conclude that both θpµq and π{α V pθ, µq flow to zero in the limit µ Ñ 0, very likely to an IR fixed point. Bar of any loop corrections, the properties of the theory can be directly read off from the fixed point couplings, to the end that confinement implies strong CP invariance. At the upper end of the curves, that is in the perturbative regime, CP is trivially conserved. Our result, that the θ parameter renormalizes to zero in the IR, does not come unexpected [45,46,47,48]. Largely identical RG equations for the running coupling and θ parameter follow from models of the vacuum based on instantons [37,46], considering their rather limited size [49,50], with basically the same conclusion for θpµq. It appears that the renormalization of θ is a generic property of instanton fluctuations. Furthermore, it has been shown analytically [47], based on an exact RG evolution equation [51], that θp0q " 0 when α S p0q " 8. This is exactly what we have found. The best known example though, where a RG flow as shown in Fig. 8 becomes effective, is the quantum Hall effect. In this case, described by the CP N´1 model at large N [48], the θ parameter has a precise parallel in the Hall conductivity, θ{2π " σ xy , which flows to quantization, θ " 0 rmod 2πs, at IR scales, viz. low temperature. For long, the quantum Hall effect has served as a model for the solution of the strong CP problem [45].

Hadron observables
Let us consider hadron observables now. The Polyakov loop susceptibility, which we considered already in Sec. 2, is a sensitive probe of the long-distance properties of the theory. The Polyakov loop describes a single static quark in the fundamental representation traveling around the periodic lattice. As such, it should be screened for nonvanishing values of θ. We show the integrated bare Polyakov loop correlator x|P| 2 y Q depending on Q as a function of t{a 2 in Fig. 9 (left panel).
Also shown is a scatter plot of P at t{a 2 " 60 (right panel). Both plots are from the 16 4 lattice, where we have the largest number of entries per charge. We find xPy " 0 in each sector. That implies center symmetry for all values of θ. The entry on the right shows that for small values of |Q| the Polyakov loop P populates the entire theoretically allowed region, while |P| stays small for larger values of |Q|. Similar results are found on the larger lattices. In the θ vacuum we have From (20) we derive the normalized Polyakov loop susceptibility which describes the connected part of the Polyakov loop correlator x|P| 2 y θ . We plot the Polyakov loop susceptibility in Fig. 10 for t{a 2 " 10 to 100. It shows that the Polyakov loop gets screened, as expected, within a narrow region around θ " 0. Furthermore, the susceptibility is independent of the flow time t, not only for θ " 0, Fig. 2, but for all values of θ. The Polyakov loop is somewhat special, as it can readily be screened by static dyon loops wrapped around the lattice [11,52]. A further key quantity is the correlation length. Above the vacuum, the energy density Ept, xq, Eq. (3), projects onto J PC " 0``glueball states. The lowest energy state, which we denote by m 0``, is called the mass gap. The inverse of the mass gap defines the correlation length, ξ " 1{m 0``, which describes the length scale over that fluctuations are correlated. The correlation length can be read off from the variance of the energy density, which is identical to the integrated glueball correlator. In the θ vacuum it reads (on the hypercubic lattice) where N " L 6 {16. In (22) the dependence on the flow time has been suppressed to avoid confusion with the Euclidean time. Explicitly, We have assumed that the correlator is dominated by the lowest glueball state, the mass gap, which should not make a difference though. In Fig. 11 we show xE 2 y θ´x Ey 2 θ as a function of θ on the 24 4 lattice for three different values of t{a 2 . Currently, we have the best control over the errors on our midsize lattice. Again, the three curves turn out to be independent of the flow time. Like the Polyakov loop susceptibility, the glueball correlator quickly drops to zero away from θ " 0. Likewise, the correlation length ξ vanishes, most probably in combination with the amplitude |xθ|E|0``y| 2 , and the theory ceases to have a (finite) mass gap.
A common feature of both the Polyakov loop and the energy density is that they are totally screened at |θ| Á 0.4 on the larger lattices, Figs. 10 and 11. It is to be noted that the Polyakov loops and energy densities in (20) and (23) overlap at small separations x. Thus, to be completely screened, the screening length must be smaller than the (color) charge radius of the static quark, for example. Our estimate of the screening length λ S at |θ| « 0.4 is λ S « 0.8 fm, assuming that the screening length is universal. This is to be compared to the typical size of a hadron of 1´2 fm, which makes our result seem credible.

Errors
For the sake of legibility, in particular to highlight the dependence on the flow time, we have not drawn error envelopes in Figs. 7, 10 and 11. The running coupling α V pθ, µq, as well as xEptqy θ , which is the Fourier transform of a convex function, is strictly positive [53]. The Polyakov loop susceptibility χ P pθq and the variance xE 2 y θ´x Ey 2 θ are positive by definition. However, after  17), is compared with a smooth fit.
the three quantities have dropped to zero, they start to oscillate around zero with a frequency of Op|Q| max q, where |Q| max is about the largest charge that has been detected, due to the finite volume and statistics. This is demonstrated in Fig. 12 for the partition function Zpθq on the 32 4 lattice. The oscillation amplitude is À 0.02 over the entire range of θ. Various techniques to filter unphysical high-frequency modes from discrete Fourier transforms have been proposed in the literature [54]. We fit the tail of the distributions to a smooth function. Alternatively, one can employ a low-pass filter, like the Savitzky-Golay smoothing filter, which practically gives the same result. With this in mind, we estimate the vertical error in Fig. 7 at 8%, but no smaller than 0.5, and the horizontal error at 10% for |θ| À 0.5 and up to Op30%q for |θ| Á 1, due to uncertainties at small flow times. In Fig. 10 the vertical error is 6%, but no smaller than 0.02, and the horizontal error is 15%. In Fig. 11 the vertical error is less than 0.2, and the horizontal error is estimated at 18%.

Free energy
The defining features of the θ vacuum are reflected in the free energy Fpθq, and the probability distribution PpQq of the topological charge Q it derives from. The free energy is given by For example, a first (second) order phase transition will arise when the first (second) derivative of Fpθq becomes discontinuous. While it will be difficult to compute Fpθq directly from the partition function Zpθq at larger values of θ, it can be reconstructed from the lower connected moments of the topological charge 1 V xQ n y θ,c "´i n B n Fpθq Bθ n . The higher the moments, the larger the volume needs to be. On the 16 4 lattice PpQq cannot be distinguished from a Gaussian distribution, PpQq " p1{ a 4πxQ 2 yq expt´Q 2 {2xQ 2 yu. The situation changes on the larger lattices. In Fig. 13 we plot the logarithm of PpQq on the 32 4 lattice. On this lattice L « 2.6 fm. We see a clear difference from a Gaussian distribution, let alone from a dilute instanton gas. The difference becomes noticeable for charges |Q| Á 16. The steeper (than Gaussian) slope of PpQq is an indication of instanton repulsion [49]. Deviations from a Gaussian distribution have been observed before in [28]. Let us now consider the ratio Note that xQ 2 y θ , being the Fourier transform of a convex function, is positive and that xQy 2 θ is inherently negative. Thus R ě 0. We plot the ratio R in Fig. 14 for our 16 4 and 32 4 lattices. A Gaussian distribution gives xQy θ "´i xQ 2 y and xQ 2 y θ,c " xQ 2 y, which leads to R " θ 2 xQ 2 y. This is approximately what we find on the 16 4 lattice. On the 32 4 lattice R starts to oscillate about R « 1{2 outside a narrow region around θ " 0, with a frequency of Op|Q| max q. The exact mean does not matter here. However, 1{2 is the most plausible value, in which case the topological charge Q is totally disconnected, considering reflection positivity. This results in the differential equation which has the solution BFpθq where ℘ is Weierstrass' elliptic function [55], and g 3 is a parameter that determines its period (2π in our case). If this continues to hold on larger volumes, the topological charge density, p1{32π 2 q xG a µνG a µν y θ " p1{Vq xQy θ , will drop to zero at |θ| Á 0 in the infinite volume limit. This would mean that the θ term will not have observable effects on hadron masses and matrix elements [56], as far as they exist. However, to substantiate this claim, further investigations on larger lattices will have to be done. That will be hard though, because Zpθq becomes an increasingly narrow function of θ and one quickly runs into the 'zero over zero' problem.
The ratio R vanishes identically at θ " 0. It needs to be seen if any firm conclusion on the order of the transition can be drawn from the behavior of the free energy in the vicinity of θ " 0. So far everything speaks for a second order phase transition, without spontaneous breaking of CP invariance [57].

Summary and outlook
The gradient flow proved a powerful tool for tracing the gauge field over successive length scales µ. It passed several tests and showed its potential for extracting low-energy quantities of the theory, highlighted by the topological susceptibility χ t , the lambda parameter Λ MS and the string tension σ. Under the gradient flow the path integral splits into quantum mechanically disconnected topological sectors, with the probability distribution PpQq for topological charge Q being independent of the flow time t. A key observation is that the long-distance properties of the theory depend sensitively on Q. It can be readily checked that already a moderate dependence of an observable on Q leads to a highly nontrivial dependence on θ.
So far we have concentrated on bulk quantities, most importantly the running coupling α V . For θ " 0 we found that α V increases quadratically with distance, α V » Λ 2 V {µ 2 9 r 2 , in accord with infrared slavery. This is an essential condition for our final conclusion. For |θ| ą 0 the color charge turned out to be screened with screening length λ S 9 1{|θ|. Thus, to maintain confinement, θ is restricted to an increasingly narrow region around zero, depending on µ, which means that the physical θ is renormalized. Analytically, the change of α V and θ with scale µ could be described by RG equations. The solution, shown in Fig. 8, connects the IR with the perturbative regime. In the IR limit α V and θ flow to 1{α V " θ " 0, possibly an infrared fixed point. This suggests that CP is conserved at all scales within QCD. To be true, the vacuum expectation value xG a µνG a µν y θ must be zero everywhere for |θ| ą 0, as has been objected in [56]. We have argued that this will be the case on large volumes, which needs to be substantiated though.
Not much changes in three-flavor QCD for heavy quarks. In Fig. 15 we plot the running coupling α GF pµq as a function of t{a 2 at the SU(3) flavor symmetric point [58], m π " m K « 410 MeV, for four different values of β. The individual curves fall onto a single straight line when transformed to a common scale. As before, the gauge fields split into disconnected topological sectors at larger flow times. In this case Λ MS is obtained from α GF pµq extrapolated to the chiral limit, followed by a continuum extrapolation. Preliminary results suggest agreement with the literature. A flavor-diagonal CP-violating phase of the quark mass matrix, which often enters phenomenological estimates of CP violating processes, can always be rotated into a θ term.
The search for an electric dipole moment d n of the neutron directly from QCD constitutes a crucial test for our results. The RG flow in Fig. 8 indicates that CP is conserved from perturbative scales down to the IR limit. The experimental upper bound on d n [1], which limits the separation of positive and negative charges within the neutron to less than 10´1 3 fm, agrees with that. It would be rather naïve to believe that any effect on that scale can be attributed to QCD, and its nonperturbative properties in particular. Recent lattice calculations do not provide a uniform picture. While the results of [59,60,61] and [62], corrected for spurious mixing effects of the Pauli and dipole form factors in [63], are consistent with d n " 0, the authors of [64] claim a nonvanishing value. One might find a nonvanishing condensate xG a µνG a µν y θ in the finite volume, as we have seen in Sec. 5. But before any lattice calculations can be trusted, the results need to be extrapolated to the infinite volume. This is particularly the case for small values of θ. In case of a vanishing dipole moment no limit on θ can be drawn from the experimental upper bound.
The bound on θ circulating in the literature is based on estimates of d n from chiral effective field theory [65,66]. However, there is the possibility that chiral symmetry is restored in the θ vacuum, given the fact that chiral symmetry breaking appears to be inseparably connected with confinement, which would deprive the calculations of their fundamental basis. Indeed, it has been argued that the chiral condensate ΣpQ, µq depending on Q grows with a power of |Q| [67,68,69], which would cause Σpθ, µq to vanish for |θ| Á 0, similar to α V pθ, µq and the Polyakov loop susceptibility. Calculations of Σpθ, µq on the lattice are in progress [70]. The nontrivial phase structure of QCD has far-reaching consequences for anomalous chiral transformations. Whatever the new phases are, our results are incompatible with the axion extension of the Standard Model. The vacuum will be unstable under the Peccei-Quinn [71] chiral transformation U PQ p1q " expti δQ 5 u, resulting in the shift symmetry θ Ñ θ`δ, which thwarts the axion conjecture.