Renormalizability of the gradient flow in the 2D $O(N)$ non-linear sigma model

It is known that the gauge field and its composite operators evolved by the Yang--Mills gradient flow are ultraviolet (UV) finite without any multiplicative wave function renormalization. In this paper, we prove that the gradient flow in the 2D $O(N)$ non-linear sigma model possesses a similar property: The flowed $N$-vector field and its composite operators are UV finite without multiplicative wave function renormalization. Our proof in all orders of perturbation theory uses a $(2+1)$-dimensional field theoretical representation of the gradient flow, which possesses local gauge invariance without gauge field. As application of the UV finiteness of the gradient flow, we construct the energy--momentum tensor in the lattice formulation of the $O(N)$ non-linear sigma model that automatically restores the correct normalization and the conservation law in the continuum limit.


Introduction and the summary
The Yang-Mills gradient flow or the Wilson flow [1] has attracted much attention in recent years in the context of lattice gauge theory. Its known application includes, the scale setting [1,2], the definition of the topological charge [1,3], the definition of the non-perturbative gauge coupling [4,5], chiral condensation [6], and the improvement of the step scaling [7], etc. Even its application to supersymmetric theory [8] and to the operator product expansion [9] is considered. Ref. [10] is a review of this notion and further related works can be found in a review [11] and in a most recent paper on the non-perturbative beta function [12].
A crucial property of the Yang-Mills gradient flow, underlying the above application, is its "ultraviolet (UV) finiteness" [1,13]. The gradient flow is a one-parameter (called the flow time) evolution of the gauge field, according to a "heat diffusion equation" (called the flow equation). A remarkable fact that can rigorously be proven [13] in all orders of perturbation theory is that any correlation function of the evolved (or flowed) gauge field becomes UV finite without the wave function renormalization, as long as parameters of the theory are properly renormalized. Moreover, any local product of the flowed gauge field remains UV finite without further (multiplicative as well as subtractive) renormalization. This remarkable property of the gradient flow facilitates, in particular, the construction of renormalized composite operators of the gauge field. That is, any simple product of the flowed (bare) gauge field as it stands is a renormalized (i.e., UV finite) quantity.
In Ref. [14], as possible application of the gradient flow, one of us (H.S.) considered the construction of the energy-momentum tensor in lattice gauge theory. This application of the gradient flow to the energy-momentum tensor was further developed from a somewhat different perspective in Ref. [15]. The construction was then generalized to gauge theories including the fermion field [16]. The genuine energy-momentum tensor cannot be defined on the lattice because the lattice structure breaks the translational invariance explicitly. Even the construction of a lattice operator that reduces to the correctly-normalized conserved energy-momentum tensor in the continuum limit is quite non-trivial as investigated in Refs. [17,18]. Ref. [19] is a pioneering work on this issue.
The basic idea of Refs. [14,16], that uses the UV finiteness of the gradient flow in an essential way, is recapitulated in Sect. 6 of the present paper. The aim of Refs. [14,16] is to construct a lattice operator that automatically reduces to the correctly-normalized conserved energy-momentum tensor in the continuum limit. Theoretically, there is only little room for doubt on the reasoning in Refs. [14,16]. Practically, however, it is not a priori clear presently-available lattice parameters are sufficient to extract physical information by using the construction. On this issue, the promising result in Ref. [20] for thermodynamical quantities in the quenched QCD is quite encouraging. Still, it is indispensable to numerically demonstrate the conservation law of the energy-momentum tensor by using lattice Monte Carlo simulations.
Under these situations, it seems useful to study a simpler system that would allow a similar construction of the lattice energy-momentum tensor using the gradient flow. One of the basic assumptions in Refs. [14,16] is that the theory is asymptotically free. Not so many field theories exhibit the asymptotic freedom, however. This was our original motivation for the present study on the gradient flow in the two-dimensional O(N ) non-linear sigma model [21][22][23]. It is well-known [24] that the physics of this systems possesses many similarities with the four-dimensional non-abelian gauge theory. These include, the asymptotic freedom, the dynamical generation of the mass gap, and, for N = 3, the topological term and the associated θ-parameter. See also Refs. [25,26]. This system is advantageous also from computational perspective (and thus from our original motivation), because there exists a very efficient cluster simulation algorithm [27,28]. The state of art in non-perturbative lattice study of the two-dimensional O(N ) non-linear sigma model can be found in Ref. [29].
In the present paper, we will show that there exists another surprising similarity between the two-dimensional O(N ) non-linear sigma model and the four-dimensional gauge theory: Any correlation function of the flowed N -vector field in the former becomes UV finite without the wave function renormalization, as long as parameters of the theory are renormalized. This UV finiteness persists also for any local product of the flowed N -vector field. This similarity is surprising, because the UV finiteness of the flowed gauge field is a non-trivial consequence [13] of the gauge BRS symmetry that acts non-linearly on the gauge field. In fact, matter fields such as the fermion field transform linearly under the gauge BRS symmetry and they do require the wave function renormalisation even after the flow [6]. In the two-dimensional O(N ) non-linear sigma model, however, it is not clear at first glance what plays the same role as this gauge BRS symmetry in the four-dimensional gauge theory. Our proof clarifies this point. On the other hand, happily, because of the UV finiteness of the gradient flow in the two-dimensional O(N ) non-linear sigma model, we can repeat the construction of the lattice energy-momentum tensor in Refs. [14,16].
The following is the organization of the present paper and the summary of the contents of each section.
In Sect. 2, we introduce the flow equation in the two-dimensional O(N ) non-linear sigma model. If one considers the application in lattice numerical simulations, this is the equation that should be solved numerically in conjunction with the conventional Monte Carlo simulations. We then formulate the perturbative expansion for the system defined by the combination of the two-dimensional O(N ) non-linear sigma model and the flow equation (the flowed system).
In Sect. 3, on the basis of the perturbative expansion developed in Sect. 2, we explicitly compute the two-point function of the flowed bare N -vector field to the one-loop order. This explicit calculation shows that the two-point function is made UV finite by the conventional parameter renormalization in the non-linear sigma model [30], but without the wave function renormalization. We carry out the computation in dimensional regularization and in lattice regularization and arrived at the same conclusion. Although this computation is only in the one-loop level, it strongly indicates that the gradient flow in the non-linear sigma model has a similar UV property as the gauge theory.
As the proof for the four-dimensional gauge theory in Ref. [13] and the renormalizability proof in the stochastic quantization [31,32], our proof in all orders of perturbation theory uses a local field theory with one spacetime dimension higher: We use a (2 + 1)dimensional field theoretical representation of the flowed system. In Sect. 4, we define this (2 + 1)-dimensional local field theory. Then we show that the system defined through the flow equation in Sect. 2 and the (2 + 1)-dimensional field theory have equivalent perturbative expansions. It is easy to see the rough equivalence. However, a closer look reveals that there are discrepancies between two systems; the measure term in the former is missing in the latter, while the former does not have the flow-line loop diagrams in the latter. Presumably, the step to show that those two apparently different elements are actually equivalent (Sect. 4.4) is the hardest part in our argument. We will find, to address this very subtle problem in a convincing manner, it is necessary to first discretize the flow-time derivative and then take the continuum limit for this discretization; this necessity of discretization is also counter-intuitive.
Once having obtained a local field theory that is (perturbatively) equivalent to the flowed system, a possible way to proceed is to write down a Ward-Takahashi relation or a Zinn-Justin equation [33] (see Ref. [34], for example) for the 1PI generating functional, 1 which restricts the possible form of counterterms, on the basis of a certain symmetry in the (2 + 1)dimensional system. This is the content of Sect. 5. Here, we encounter another surprise: The (2 + 1)-dimensional field theory possesses local gauge symmetries, although it does not contain any gauge field. Note that the unique internal symmetry in the original two-dimensional O(N ) non-linear sigma model is the global O(N ) symmetry. Because of these gauge symmetries, we have to fix the gauge. Even under the gauge fixing, there still remains a residual symmetry that acts non-linearly on the N -vector field. We will find the Zinn-Justin equation associated with this non-linear symmetry does the job. Then by listing up possible counterterms (by borrowing the information obtained in Sect. 4.4) and examining the restriction implied by the Zinn-Justin equation, finally we show that only counterterms required are those of the original two-dimensional O(N ) non-linear sigma model. In particular, the flowed N -vector field (and its composite operators) is not renormalized. This completes our proof for the UV finiteness of the gradient flow.
In Sect. 6, on the basis of the UV finiteness established in Sect. 5, we construct the energymomentum tensor in lattice formulation of the non-linear sigma model, following the line of reasoning of Refs. [14,16].
In summary, we have found another example in which the gradient flow exhibits a remarkable UV finiteness; in the two-dimensional O(N ) non-linear sigma model, any correlation function of the flowed N -vector field and its composite operators is UV finite without the multiplicative (as well as the subtractive) renormalization. Our proof in the present paper also clarified subtle but very interesting technical issues arising in the theoretical analysis of the gradient flow, such as the necessity of the discretization of the flow-time derivative and the emergence of gauge and/or non-linear symmetries in the corresponding local field theory with one dimension higher. The knowledge obtained here will be useful in considering the application of the gradient flow to a wider range of systems.
Also, going back to our original motivation, we hope to numerically test the idea of Refs. [14,16] by using the energy-momentum tensor constructed in Sect. 6 in the near future.
2. Gradient flow in the two-dimensional O(N ) non-linear sigma model 2.1. Two-dimensional O(N) non-linear sigma model and the flow equation The two-dimensional O(N ) non-linear sigma model is a field theory of an N component vector with the unit length. Its partition function is given by and g 0 is the bare coupling constant. Although the spacetime dimension D is 2 for our target theory, expressions for generic D is useful because we will extensively use dimensional regularization in what follows.
In the present paper, as an analogue of the Yang-Mills gradient flow [1], we consider the following t-evolution of the N -vector field (the flow equation): where the initial value is given by the N -vector field in the O(N ) non-linear sigma model, which is the subject to the functional integral (2.1). The projection operator P ij (t, x) in the right-hand side of the flow equation (2.2) is defined by (in Eq. (2.2) and in what follows, the sum over the repeated index is understood). The projection operator is introduced so that the flow is consistent with the constraint n(t, x) 2 = 1, where n(t, x) 2 ≡ N i=1 n i (t, x)n i (t, x), that is, ∂ t n(t, x) 2 = 0. The latter would be a natural requirement for the flow equation for the O(N ) non-linear sigma model. In fact, a flow equation identical to Eq. (2.2) has been advocated also in Appendix B of Ref. [8] from a perspective of the symmetry of the present system. 2

Perturbative expansion
As usual, for the perturbative treatment of the O(N ) non-linear sigma model, we parametrize the constraint n(x) 2 = 1 in Eq. (2.1) in terms of N − 1 independent components (the π-field) as and then expand expressions regrading π(x) as a small fluctuation. In this perturbative treatment, the partition function becomes (2.7) The above arbitrary choice of the perturbative branch, Eq. (2.6) with small π(x), however induces the infrared (IR) divergences in the perturbative expansion of O(N ) non-invariant quantities [35]. To regularize the IR divergences, we introduce the mass term and take the massless limit m 0 → 0 in the very end of calculation. With this mass term, the particular perturbative branch (2.6) is favored for a weak coupling. Also for the flowed field n i (t, x), since n(t, x) 2 = 1 holds along the flow evolution, we set Then the perturbative expansion of the flow equation (2.2) is obtained from the integral representation, and (2.14) Noting that (∂ t − ∂ µ ∂ µ )K t (x) = 0 and K t=0 (x) = δ D (x), we see that Eq. (2.11) solves Eq. (2.2) with the initial condition (2.3). By iteratively solving Eq. (2.11) in terms of the initial value π k (y), therefore, we have a perturbative solution of the flow equation. This expansion can be represented diagrammatically (the flow Feynman diagram [13]) and, throughout this paper, we represent the heat kernel (2.13) by a doubled wavy line in Fig. 1. This line is also called the "flow-line propagator" or simply the "flow line". The initial value of the flow, π k (y) in Eq. (2.11) is a quantum field subject to the functional integral (2.7). From Eq. (2.11) and Eq. (2.7) (with the mass term (2.8)), one then sees that the quantum free propagator of the flowed field is given by π k (t, x)π l (s, y) Note that in this propagator, the flow times at the end points appear in the sum (not the difference). Throughout this paper, this free propagator will be denoted by a single wavy line (Fig. 2). Finally, the functional integral (2.7) generates interaction vertices among π k (x)'s. The interaction vertices in the action integral will be denoted by a filled circle (see Fig. 3 for the example). On the other hand, the interaction vertices arise from the functional measure in Eq. (2.7), the "measure term", will be represented by a cross as Fig. 5.

One-loop calculation of correlation functions of the flowed field
An explicit one-loop calculation of correlation functions of the flowed field is quite instructive, because it shows a remarkable UV property of the gradient flow. As the UV regularization, we first adopt dimensional regularization, setting where γ E is Euler's constant.
On the other hand, the contribution of the another one-loop diagram, diagram 02 in Fig. 4, that contains the flow vertex is We note that the measure term in Eq. (2.16) vanishes identically in dimensional regularization with which δ D (0) ≡ 0. Thus, in total, we have π k (t, x)π l (s, y) Now, the parameter renormalization in the original O(N ) non-linear sigma model (2.7) with the mass term (2.8) is known to be (in the minimal subtraction (MS) scheme) and (3.6) where Z 3 is the wave function renormalization factor for the un-flowed π-field, π k (x) = Z 1/2 3 π k R (x). 4 If we make these substitutions in Eqs. (3.4), we obtain the following completely UV finite expression, π k (t, x)π l (s, y) Remarkably, when expressed in terms of renormalized parameters, the two-point function of the flowed π-field is UV finite without the multiplicative wave function renormalization. 5 This UV finiteness of the flowed field is similar to that of the four-dimensional gauge field flowed by the Yang-Mills gradient flow, a property first observed in Ref. [1] in lower order perturbative computations and then proven in all orders of perturbation theory in Ref. [13]. The above result indicates that by a similar mechanism as the four-dimensional gauge theory, the Nvector field flowed to positive flow times is UV finite only with the parameter renormalization. It is also instructive to repeat the above calculation by using lattice regularization instead of dimensional regularization. We adopt the prescription that in Eq.
x , where a denotes the lattice spacing, and the derivative ∂ µ is replaced by the forward difference operator. The Laplacian in the flow equation where ∂ µ and ∂ * µ are the forward and the backward difference operators, respectively. Then the contribution of Fig. 3 is which is quadratically divergent. The quadratic divergence in the last term is cancelled by the measure term (2.16) with δ D (0) → 1/a 2 for lattice regularization. In fact, the contribution of the measure term to the two-point function (Fig. 5) is (3.10) On the other hand, the contribution of Fig. 4 is Thus, we have in total π k (t, x)π l (s, y) It is obvious that all UV divergences are removed by the parameter renormalization (3.5) and (3.6) with the replacement 1/ǫ → − ln a; again, remarkably, no wave function renormalization is required. Although the two-point function (3.7) is UV finite, it contains IR divergences (i.e., it diverges for m → 0) because it is not an O(N ) invariant "physical" quantity [35]. As a simple example of an IR finite O(N ) invariant observable, we can consider the "energy density", defined by which is analogous to the energy density introduced in Ref. [1] for the gauge theory. For the vacuum expectation value, there are four flow Feynman diagrams to the next-to-leading order, as depicted in Figs. 6-9 (the cross denotes the operator E(t, x)). A straightforward calculation using dimensional regularization yields This is IR finite as expected and UV finite in terms of the renormalized coupling constant, again indicating the UV finiteness of the flowed field. If this UV finiteness persists to all orders (we will prove this in a later section), the result (3.15) shows that the combination t E(t, x) provides a possible non-perturbative definition of a renormalized coupling as the gradient flow scheme in the four-dimensional gauge theory (see Refs. [4,5], for example). That is, we can set Then it must be interesting to investigate the running of this non-perturbative coupling in numerical lattice simulations, in view of the expected conformal and walking behaviors of the O(3) non-linear sigma model with non-zero θ-parameters [26].

(D + 1)-dimensional field theoretical representation of the gradient flow
In the next section, we reveal the renormalization structure of the flowed system defined in Sect. 2. We prove in particular that the flowed N -vector field does not require the wave function renormalization. Our strategy is identical to the case of the four-dimensional gauge theory [13]; we seek a (D + 1)-dimensional local field theory that reproduces the flow Feynman rules in the preceding sections and use this to show the renormalizability. We neglect the IR regulating mass term (2.8) in this section, because it complicates the argument destroying the O(N ) symmetry. We will consider the effect of the mass term at the very end of the next section.

Partition function
As Refs. [6,13], we consider a (D + 1)-dimensional (D = 2 for our target theory) field theory defined in the half space, (t, x) ∈ [0, ∞) × R D that (at least perturbatively) is equivalent to the gradient flow in the two-dimensional O(N ) non-linear sigma model. We will find that, to resolve subtleties associated with the measure term and the flow-line loop (see Sect. 4.4), it is necessary to specify a prescription for the flow-time derivative. We will use the forward difference prescription (with the discretization length ǫ) for this. 6 The regularization for the D-dimensional "spacetime" direction is on the other hand arbitrary and we may assume for instance dimensional regularization or lattice regularization. The partition function of the (D + 1)-dimensional field theory we consider is defined by where t = 0, ǫ, 2ǫ, . . . , and 2) can equivalently be written as P ij (t, x)∂ t n j (t, x) with the projection operator P ij (t, x) in Eq. (2.4). The integration over the another Lagrange multiplier ξ i (x) in Eq. (4.1), on the other hand, imposes the initial condition (2.3).
It can be shown that, with the factor 1 − n ⊥ (t + ǫ, x) 2 in the integration measure, the partition function Z (4.1) can be obtained from the original partition function Z O(N ) (2.1) by inserting unity (up to infinite gauge volume; see below). However, since 1 − n ⊥ (t + ǫ, x) 2 = 1 + O(ǫ 2 ) → 1 for ǫ → 0, this factor can be neglected in the ǫ → 0 limit and we do not explicitly include this factor in what follows.

Symmetries and the gauge fixing
The above (D + 1)-dimensional system possesses following symmetries. One is the global O(N ) symmetry that is inherited from the original O(N ) non-linear sigma model: where ǫ ji = −ǫ ij are infinitesimal constant parameters. Other, somewhat unexpected ones are local gauge symmetries: and g(x) and h(t, x) are local parameters which can depend their arguments. These local symmetries, which exist even with the discretized flow-time and D-dimensional regularization, follow from the constraints n(x) 2 = n(t, x) 2 = 1 in the functional integral and the property n i (t, x)P ij (t, x) = 0. Because of these gauge symmetries, the partition function (4.1) itself is infinite. This is not a problem in our present context, because what we need at this moment is a generating functional of the perturbative expansion of the flowed system. To formulate perturbation theory in the above (D + 1)-dimensional field theory, we thus have to first fix the gauge symmetries (4.5). For this, we adopt the following gauge fixing conditions, and follow the Faddeev-Popov procedure. Thus we insert unity into the functional integral (4.1). Then, using the invariance of the action and the functional measure under the transformations (4.5), we can mod out the gauge volume from the partition function (4.1). We further solve the constraints n(x) 2 = n(t, x) 2 = 1 in terms of N − 1 independent components, as Eqs. (2.5) and (2.6) and Eqs. (2.9) and (2.10). Then, after the gauge volume is factored out, the partition function is given by where and where the combination R k (t, x) is defined by Eq. (2.14) and

Feynman rules in the (D + 1)-dimensional system
Next we derive the Feynman rules in the above system (4.10)-(4.13). To write down the free propagator, we introduce the heat kernel with the discretized flow time, by which fulfills Clearly, K ǫ t (x) reduces to the heat kernel (2.13) in the continuum flow-time limit, K ǫ t (x) ǫ→0 → K t (x). By using this object, we change the integration variables from π k (t, x) to p k (t, x) as [13], Then the action becomes where abbreviated terms are cubic or higher in fields. It is then straightforward to find free propagators and the result is where ϑ(t) is a "regularized" step function, Note that ϑ(0) = 0 (not 1/2, for example). Since other free propagators among π k (x), p k (t, x), λ k (t, x), and ξ k (x) vanish, Eqs.
In passing, we note x) in the limit ǫ → 0, combined with the above πλ-propagator, precisely reproduces the last term of the integral equation (2.11) (i.e., the flow vertex).
Thus, we have observed that our present (D + 1)-dimensional system basically reproduces the perturbative expansion of the flowed system defined in Sect. 2; they seem to be basically equivalent. Nevertheless, we should note that the equivalence appears not quite complete. The measure term (2.16) is missing in Eqs. (4.10)-(4.13) (the factor x 1 − π(x) 2 / 1 − π(x) 2 becomes unity under the integration over ξ and this is not the measure term). Although the measure term (2.16) identically vanishes when one uses dimensional regularization, the measure term plays an important role for other regularization, such as lattice regularization. If the equivalence including the measure term does  not hold, then the renormalizability proof in the next section, that is based on the present (D + 1)-dimensional field theory, does not apply to the gradient flow with, for example, lattice regularization. Then, the UV finiteness of the gradient flow with lattice regularization, that we observed through an explicit calculation in Sect. 3, is not explained by the proof.
We will find that, rather surprisingly, the measure term is generated from naively-O(ǫ) terms in the action (4.12). The aim of the next subsection is to clarify this point and to establish the perturbative equivalence between the above (D + 1)-dimensional system and the flowed system in Sec. 2.
Next we note that the perturbative expansion of Eq. (4.10) generates loop diagrams consisting solely of the flow-line propagator (4.23). Such "flow-line loop diagrams" are depicted in Figs. 10 and 11. 7 It is now very important to recognize that there is no counterpart of the above flowline loop diagrams in the perturbative expansion of the flowed system in Sec. 2. This is a consequence of the retarded nature of the flow equation and one can confirm this by drawing flow-line diagrams starting from Eq. (2.11). Thus, there appears some (apparent; see below) discrepancy between perturbative expansions of the above two systems.
Let us begin our investigation from the flow-line loop diagram in Fig. 10 which starts and ends at the same flow vertex. First note that the πλ-propagator (4.23) vanishes when the flow time of λ is greater or equal to the flow time of π. Therefore, in Eq. (4.12), the genuine flow vertex containing the non-linear term R k (t, x) does not contribute to the flow-line loop diagram in Fig. 10. What contribute is the self-contraction in the combination E (4.13). The 7 The flow-line loops cannot become higher than one-loop, because the flow vertex is linear in λ.
If we Taylor expand π k (t + ǫ, x) in this expression with respect to ǫ, we find Since this is a total derivative, only the boundary field π k (t = 0, x) = π k (x) is contained. Then, remarkably, Eq. (4.28) coincides with the measure term (2.16) for ǫ → 0. The above result (4.28) can be obtained in a somewhat different manner. We first Taylor expand −E which yields and the delta function at the equal flow-time is interpreted as δ(0) = 1/ǫ. The self-contraction in Eq. (4.29) thus cancels the factor ǫ and leaves the O(1) result, Eq. (4.28). Next, we see that a flow-line loop diagram which contains a plurality of flow vertices, such as the diagram in Fig. 11, vanishes as ǫ → 0. A little thought shows that all vertices in such a flow-line loop diagram must be the vertex arising from E (4.13). This is again because πλ-propagator (4.23) vanishes when the flow time of λ is greater than or equal to that of π. The vertex is O(ǫ) as Eq. (4.29).
The integration of the flow time of each vertex eliminates one delta function and finally one is left with an overall integration and δ(0) = 1/ǫ. In the present case of a plurality of flow vertices, however, the power of ǫ coming from vertices is always greater than two; thus the flow-line loop diagram vanishes for ǫ → 0. The conclusion is that Eq. (4.28) is the unique contribution of the flow-line loop diagrams for ǫ → 0.
By similar reasoning, it can be confirmed that Eq. (4.28) is the unique place in which an apparent O(ǫ) term in the action contribute in the ǫ → 0 limit. The integration of the flow time of each vertex eliminates one delta function and the singularity δ(0) = 1/ǫ can arise only from the flow-line loop diagrams, the case already considered above. This observation justifies the Taylor expansion with respect to ǫ and neglecting O(ǫ) terms besides that particular term in Eq. (4.29).
Thus, we have observed that the perturbative expansions in the above two systems are equivalent by a remarkable mechanism: A flow-line loop diagram in the (D + 1)-dimensional system, which is not generated in the perturbative expansion of the original flow equation, reproduces the measure term (2.16) which is absent in the original partition function of the (D + 1)-dimensional system, Eq. (4.10). The mechanism is remarkable, because an apparent O(ǫ) term in the action, i.e., E (4.13), plays the crucial role through the flow-line loop. Now, having established the equivalence between the (D + 1)-dimensional field theory (4.10)-(4.13) and the flowed system in Sec. 2, we are ready to prove the UV finiteness of the gradient flow in Sect. 2.

Proof of the renormalizability of the gradient flow
In this section, on the basis of the (D + 1)-dimensional field theory in the preceding section, we show that any correlation function of the flowed N -vector field in terms of the renormalized coupling is UV finite, without the wave function renormalization.

Residual non-linear symmetry
We first note that even with the gauge fixing (4.7), there remains a residual symmetry that is a particular combination of the global O(N ) symmetry (4.4) and the local symmetries (4.5). It is given by the requirement that it does not affect the gauge fixing conditions. That is, From these, we have Under this residual symmetry, other field components transform as where indices k and l run over only from 1 to N − 1. The interesting part in the above residual symmetry is the O(N )/O(N − 1) part corresponding to the choice of parameters ǫ kl = 0. Writing ǫ k ≡ ǫ kN , it induces the following non-linear transformations It can directly be confirmed that the integration measure and the action in Eqs.

Ward-Takahashi relation or the Zinn-Justin equation
We can express the invariance of the system under the non-linear transformation (5.5) as an identity for the generating functional of 1PI correlation functions. First, we introduce the source terms for elementary fields, except for the Lagrange multiplier field ξ k (x). It turns out that this omission of the ξ-source greatly simplifies the discussion of the renormalization. This implies that we omit correlation functions including ξ k (x) from our consideration. However, since the ξ-field appears only in the quadratic (i.e., free) part of the action S only linearly, if 1PI correlation functions of other elementary fields turn to be UV finite after renormalization, any correlation functions including the elementary ξ-field are also UV finite. Hence nothing is lost by the omission of the ξ-source for our present purpose.
To write down the Ward-Takahashi relation associated with the symmetry (5.5), we supplement also additional terms to the action, as where where the source K k,l1...ln (t, x) is symmetric in indices (l 1 , . . . , l n ) by definition. We now consider the variation of integration variables of the form of Eq. (5.5) in the partition function, We note Then by the standard argument, invariance of the integration measure and of S imply that the generating functional of 1PI functions, defined by the Legendre transformation, where π k (x), π k (t, x), and λ k (t, x) denote expectation values of elementary fields, follows an identity In writing down this identity, we have taken the continuum flow time limit ǫ → 0. This is justified because we have observed that the symmetry (5.5) is preserved by the flow-time discretization. Also, we have observed that the (D + 1)-dimensional system Eq. (4.10)-(4.13) with ǫ → 0 reproduces the perturbative expansion of the flow equation. Thus we can study the renormalizability of the flowed system in Sect. 2 by using the identity (5.15).

Structure of the renormalization
Our statement of the renormalizability is that the 1PI generating functional Γ can be made UV finite in terms of renormalized quantities, by appropriately choosing the constants Z and Z 3 in (5.16) order by order in perturbation theory. In particular, we claim that the flowed or "bulk" fields, π k (t, x) and λ k (t, x), do not require the multiplicative renormalization. Our argument proceeds by mathematical induction based on the loop expansion. We set where Γ (ℓ) is the generating functional in the ℓ-th loop order. The above assertion is certainly true for ℓ = 0 (tree level approximation) for which Z = Z 3 = 1 is sufficient. Then suppose that in perturbation theory with renormalized quantities fixed, the constants Z and Z 3 in Eq. (5.16) can be chosen so that Γ (0) , . . . , Γ (ℓ) , are UV finite in terms of renormalized quantities. Then consider the (ℓ + 1)-th loop order calculation on the basis of the above chosen Z and Z 3 . Since Z and Z 3 have already been chosen so that Γ (0) , . . . , Γ (ℓ) , are finite, by considering the UV divergent part of the identity (5.15) in the (ℓ + 1)-th loop order, we have where Γ (ℓ+1)div denotes the UV divergent part of Γ (ℓ+1) and and We next study the most general form of the divergent part Γ (ℓ+1)div . First of all, by a general theorem, the divergent part must be an integral of a local polynomial of fields and their derivatives. We then note that there is no divergence corresponding to a local term in the "bulk" t > 0, a term that is written as an ∞ 0 dt d D x integral of a local polynomial of fields and their derivatives: As we explained in detail in Sect 4.4, there is no loop diagram consisting solely of the "flow line" πλ-propagator (4.23), other than the diagram in Fig. 10 that reduces to the measure term at the boundary t = 0, Eq. (2.16). Then, since the ππpropagator (4.22) possesses the Gaussian damping factor e −(t+s)p 2 (for ǫ → 0), any loop diagram in which the flow times of the vertices (they must be the same for the divergent part) are positive is UV finite. Therefore, there is no divergence that is written as the bulk integral.
Any divergent part is thus written as the integral on the boundary t = 0. Noting that for D = 2 the fields π k R (x) and π k (t, x) possess the the mass dimension 0, H R (x), λ k (t, x), and K k,l1...ln (t, x) possess 2, and H(t, x) possesses 4, the most general possible form of the divergent part is (5.21) where B contains at most two derivatives and E k,l1...ln is symmetric in indices (l 1 , . . . , l n ). Note that we did not include the flow field at zero flow time, π k (0, x), in the possible form of diverges (5.21). The redundancy to use this field variable in addition to π k R (x) follows from 21 the relation π k (0, x) = π k (x), (5.22) that is the expectation value of the variation of the action with respect to the ξ-field. Note that here field variables denote the expectation values with the presence of source fields and not the integration variables in the functional integral. This identity shows that as the arguments of the 1PI generating functional, the variables π k (0, x) and π k (x) cannot be independent, because they cannot take different values for any configuration of the source fields. We note also that the combination and , n ≥ 2. (5.29) The above conditions for B and C, Eqs. (5.25) and (5.26), are completely identical to conditions on the divergent part in the original two-dimensional O(N ) non-linear sigma model [30]. The general solution to these is given by [30] where δZ and δZ 3 are constants. Next, from the linearly-realized O(N − 1) symmetry (that is preserved in our all steps), one can set D k = π k R (x)d(π R (x) 2 ). Then Eq. The 1PI generating functional in the (ℓ + 1)-th loop order, Γ (ℓ+1) , is thus made UV finite. This completes mathematical induction for the renormalizability. In particular, we showed that there is no need of the wave function renormalization for the flowed fields, π k (t, x) and λ k (t, x). We have shown that any correlation function of the flowed fields is UV finite under the conventional parameter renormalization, without the multiplicative wave function renormalization. Then, it is easy to see that, because of Gaussian damping factors in propagators, this UV finiteness holds even when some spacetime coordinates of the correlation function coincide, i.e., even in the equal-point limit, as long as all flow-time coordinates of flowed fields are strictly positive. 8 The local product of any number of flowed fields does not contain UV divergences. This robust UV finiteness, that makes the construction of renormalized composite operators straightforward, is the key property in application of the gradient flow in lattice field theory.
In renormalized perturbation theory, one uses the propagators and the vertices in terms of renormalized parameters and renormalized fields. This renormalized Feynman rule is obtained by making the substitution (5.16) in the action (4.12). The part including the ξ-field becomes (5.37) and the second term is regarded as the perturbation. In this renormalized perturbation theory, from the first term, the free propagator is given by while the second term is regarded as a counterterm. In this way, we can use Eq. (5.38) also for π k R (x) by identifying π k R (x) = π k (0, x). As the πξ-propagator (4.24) shows, the second term in Eq. (5.37) acts as a two-point vertex at the boundary t = 0 that connects between π k (t, x) and π l R (y). This counterterm thus plays the same role as the boundary counterterm ∆S bc in the gauge theory (Sec. 7.1 of Ref. [13]).
Finally, the IR regulating mass term (2.8) can readily be incorporated in the above argument by the substitution (5.39) In particular, from Eq. (5.16), we see that the generating functional becomes UV finite in terms of 3 )m 2 , as we already noted in Eq. (3.6).

Lattice energy-momentum tensor
In the preceding section, we have shown that any local product (the composite operator) of the bare flowed N -vector field becomes UV finite under the conventional parameter renormalization in the two-dimensional O(N ) non-linear sigma model. As application of this fact, in the present section, we consider the construction of the energy-momentum tensor, the Noether current associated with the translational invariance, in lattice formulation of the non-linear sigma model. The idea is the same as Refs. [14] and [16]: Since lattice regularization explicitly breaks the translational invariance, the construction of the energy-momentum tensor is awkward. Instead of considering this construction directly, we construct a composite operator of the flowed field which, under dimensional regularization, becomes the energymomentum tensor. Since dimensional regularization preserves the translational invariance, the description of the energy-momentum tensor that fulfills the correct Ward-Takahashi relation is straightforward. On the other hand, since the composite operator of the flowed field is UV finite under the parameter renormalization, it must become independent of the regularization in the limit that the regulator is removed (after the renormalization; as long as the same renormalization conditions are adopted). In this way, low energy correlation functions of the energy-momentum tensor may be computed by using lattice regularization. This construction in Ref. [14] has been applied to the thermodynamics of the quenched QCD in Ref. [20] and promising results have been obtained.

Energy-momentum tensor with dimensional regularization
The energy-momentum tensor T µν (x) can be obtained from the variation of the action (2.1) under the infinitesimal translation with a localized parameter, The explicit form is given by Assuming that we are using dimensional regularization, which preserves the translational invariance, the above classical expression as it stands fulfills the correct Ward-Takahashi relation associated with the translational invariance: In this expression, D is a bounded integration region and O ext is an operator outside the region D and O int is an operator inside the region. We defined the renormalized energy-momentum tensor by subtracting the vacuum expectation value, {T µν } R (x) ≡ T µν (x) − T µν (x) . The Ward-Takahashi relation ensures that the bare quantity T µν (x) is not multiplicatively renormalized. Although naively the energy-momentum tensor (6.3) is traceless for D = 2, UV divergences in the composite operator (1/g 2 0 )∂ ρ n i (x)∂ ρ n i (x) being proportional to 1/ǫ makes this expectation invalid even for ǫ → 0. Instead, we have the the trace anomaly, where the MS scheme is assumed in the renormalized composite operator in the right-hand side and the coefficient is given by the β function, (here, the derivative with respect to the renormalization scale µ is taken while bare quantities are kept fixed) and [36][37][38],

Small flow time expansion and the energy-momentum tensor
We construct a composite operator of the flowed field which reduces to the two-dimensional composite operator (6.3) by using the small flow time expansion introduced in Ref. [13]. For this, we take an O(N ) invariant dimension 2 second-rank composite operator of the flowed field, According to the argument in Ref. [13], for t → 0, this composite operator of the flowed field can be expressed as a series of two-dimensional local operators with increasing mass dimensions, as Similarly, we have Inverting these relations with respect to the two-dimensional operators and substituting them into Eq. (6.3), we have where Hence, if we know the t → 0 behavior of the coefficients ζ IJ (t) in Eqs. (6.10) and (6.11), then the energy-momentum tensor (6.3) can be obtained by the t → 0 limit of the right-hand side of Eq. (6.12). Thus, we are interested in the t → 0 behavior of the coefficients ζ IJ (t). Since all the composite operators in the above expansions are bare ones, by the standard renormalization group argument, the expansion coefficients are independent of the renormalization scale q, if they are expressed in terms of the running parameterḡ(q). In particular, we may take q = 1/ √ 8t. Then because of the asymptotic freedom, the running coupling behavesḡ(1/ √ 8t) → 0 for t → 0 and ζ IJ (t) for t → 0 can be evaluated by perturbation theory. To find the coefficients ζ IJ (t) in Eq. (6.10), we consider correlation functions, and and compute the coefficients of combinations, in tensors M µν (p, q; t) and M µν (p, q). Then, we determine ζ IJ (t) so that the relation (6.10) holds in correlation functions in view of the combinations (6.17). Set ζ IJ (t) = ζ (0) IJ (t) + · · · , where superscripts denote the number of loops. In the tree level, ζ   We use dimensional regularization to regularize UV divergences and the mass term (2.8) to regularize IR divergences.
The contribution of each diagram to M µν (p, q; t) in Eq. (6.15) and to M µν (p, q) in Eq. (6.16) is tabulated in Tables 1 and 2. From these results, we have      Note that IR divergences are cancelled out in the coefficients. Using these results in Eqs. (6.13) and (6.14), we have c 1 (t) = 1 g 2 − 1 4π (N − 2) ln(8πµ 2 t) + O(g 2 ), (6.21) c 2 (t) = 1 4π (N − 2) + O(g 2 ) = b 0 + O(g 2 ), (6.22) where g is the renormalized coupling in the MS scheme (3.5). Note that these coefficients are UV finite in terms of the renormalized parameter. This must be so, because all composite operators appearing in Eq. (6.12) are renormalized ones. For c 2 (t) (6.22), one may proceed one step further [14] by requiring Eq. (6.12) reproduces the trace anomaly (6.5) to the two-loop order. By taking the trace of Eq. (6.3) and comparing it with Eq. (6.5), we find ∂ ρ n i ∂ ρ n i R (x) = 1 + O(g 4 ) ∂ ρ n i (x)∂ ρ n i (x) − ∂ ρ n i (x)∂ ρ n i (x) . (6.23) Using Eq. (6.11) with Eq. (6.20), we find that for Eq. (6.12) to reproduce the trace anomaly (6.5) to the two-loop order, (6.24) The expression for the energy-momentum tensor that is usable with lattice regularization is thus given by the t → 0 limit of Eq. (6.12) with the coefficients in Eqs. (6.21) and (6.24). As noted above, we can replace the renormalization constant g and the renormalization scale µ in Eqs. (6.21) and (6.24) by the running couplingḡ(q) with the renormalization scale q and set q = 1/ √ 8t. We may use for example the four-loop running coupling [39], in actual numerical simulations.

A facile computational method for ζ IJ (t)
In the above calculation of the matching coefficients ζ IJ (t), we have regularized IR divergences by introducing the bare mass m 0 for the N -vector field. The required computation is as the result somewhat troublesome. In this final subsection, we point out that at least in the one-loop level there exists a "facile method" that avoids the introduction of the IR regularizing mass [40]. This method has been particularly useful for gauge theories [14,16] because one can regularize IR divergences without intruding a gauge breaking mass parameter; IR divergences are regularized by "dimensional regularization".
We first note that for a Feynman diagram that contribute to Eq. (6.16), for example the diagram in Fig. 15, there always exists a corresponding flow Feynman diagram that contributes to Eq. (6.15). The topology of both diagrams is identical (Fig. 15 for the present example) but in the latter the propagators carry the Gaussian damping factor e −tℓ 2 , where ℓ is the loop momentum, as Eq. (2.15). As we have observed above, what is relevant to ζ IJ (t) is the difference of the values of these two diagrams which, by dimensional counting, has the structure, as the analytic continuation from Re(D) > 2. Now, interestingly, the right-hand sides of Eq. (6.27) and that of Eq. (6.29) are identical as a function of D. Thus, we may use the latter method instead of the former. The latter is computationally much simpler because only flow Feynman diagrams have to be computed and the IR regulator m 0 is not necessary. We tabulated the result of this facile method in Table 3. It can be confirmed that each entry coincides with the difference between corresponding entries of Table 1 and of Table 2 as it must be the case; the resulting matching coefficients ζ IJ (t) directly obtained from Table 3 are of course identical to the previous ones, Eqs. (6.18)-(6.20).