format: text

text:
# Geometric numerical integration of nonholonomic systems and optimal control problems

## 1 Introduction

Standard methods for simulating the motion of a dynamical system usually ignore many of the geometric features of this system (simplecticity, conservation laws, symmetries…). However, new methods have been recently developed, called geometric integrators, which are concerned with some of the extra features of geometric nature of the dynamical system (see \[HaLuWa:02\]).

In the first part of the paper, we propose a class of geometric integrators for nonholonomic systems \[Leomar:96D,NeiFuf:72\] based on a discretization of the Lagrangian function (in a more precise sense, we discretize the action function) and a coherent discretization of the constraint forces (see \[LeMaSa:02\]). These equations will be conceptually equivalent to the proposed for systems with external forces (see \[MarWes:01\]). Finally, second part corcerns with the construction of symplectic integrators for optimal control theory by using generating functions of the second kind.

## 2 Nonholonomic systems

### 2.1 Geometrical formulation of nonholonomic systems

Let $`Q`$ be a $`n`$-dimensional differentiable manifold, with local coordinates $`(q^i)`$ and tangent bundle $`TQ`$, with induced coordinates $`(q^i,\dot{q}^i)`$. Consider a Lagrangian system, with Lagrangian $`L:TQ\text{}`$, subject to nonholonomic constraints, defined by a submanifold $`D`$ of the velocity phase space $`TQ`$. We will assume that $`dimD=2nm`$ and that $`D`$ is locally described by the vanishing of $`m`$ independent functions $`\varphi ^a`$ (the “constraint functions”), satisfying the rank condition $`\text{rank }\left({\displaystyle \frac{\varphi ^a}{\dot{q}^i}}\right)=m`$. In the sequel, we will follow a Hamiltonian point of view. The canonical coordinates on $`T^{}Q`$ (the cotangent bundle of $`Q`$) are denoted by $`(q^i,p_i)`$. Assume, for simplicity, that the Lagrangian $`L`$ is hyperregular, that is, the Legendre transformation $`Leg:TQT^{}Q,(q^i,\dot{q}^i)(q^i,p_i=L/\dot{q}^i)`$, is a global diffeomorphism. The constraint functions on $`T^{}Q`$ become $`\mathrm{\Psi }^a=\varphi ^aLeg^1`$, i.e. $`\mathrm{\Psi }^a(q^i,p_i)=\varphi ^a(q^i,{\displaystyle \frac{H}{p_i}}),`$ where the Hamiltonian $`H:T^{}Q\text{}`$ is defined by $`H=E_LLeg^1`$. Here, $`E_L`$ denotes the energy of the system, locally defined by $`E_L=\dot{q}^i\frac{L}{\dot{q}^i}L`$. Since locally $`Leg^1(q^i,p_i)=(q^i,{\displaystyle \frac{H}{p_i}})`$, then $`H=p_i\dot{q}^iL(q^i,\dot{q}^i),`$ where $`\dot{q}^i`$ is expressed in terms of $`q^i`$ and $`p_i`$ by using $`Leg^1`$.

The equations of motion for the nonholonomic system on $`T^{}Q`$ can now be written as follows (see \[CaLeMa:99,Marl:95\] and references therein)

$$\{\begin{array}{ccc}\hfill \dot{q}^i& =& \frac{H}{p_i}\hfill \\ \hfill \dot{p}_i& =& \frac{H}{q^i}\lambda _a\frac{\mathrm{\Psi }^a}{p_j}_{ji},\hfill \end{array}$$
(1)
together with the constraint equations $`\mathrm{\Psi }^a(q,p)=0`$, where $`_{ij}`$ are the components of the inverse of the matrix $`(^{ij})=(^2H/p_ip_j)`$. Note that

$$(\frac{\mathrm{\Psi }^a}{p_j}_{ji})(q,p)=(\frac{\varphi ^a}{\dot{q}^i}Leg^1)(q,p).$$

Let $`M`$ denote the image of the constraint submanifold $`D`$ under the Legendre transformation, and let $`F`$ be the distribution on $`T^{}Q`$ along $`M`$, whose annihilator is given by $`F^o=Leg_{}(\stackrel{~}{F}^o))`$. Here, $`\stackrel{~}{F}^o`$ represents the constraint forces subbundle, locally defined by

$$\stackrel{~}{F}^o=\text{span}\{\mu ^a=\frac{\varphi ^a}{\dot{q}^i}dq^i\}$$

The Hamiltonian equations of motion of the nonholonomic system can be then rewritten in intrinsic form as

$$\begin{array}{ccc}\hfill (i_X\omega _QdH)_{|M}& & F^o\hfill \\ \hfill X_{|M}& & TM,\hfill \end{array}$$
(2)
where $`\omega _Q=d\theta _Q=dq^idp_i`$ (with $`\theta _Q=p_idq^i`$) is the canonical symplectic form on $`T^{}Q`$. Suppose in addition that the following compatibility condition $`F^{}TM=\{0\}`$ holds, where $`\mathrm{`}\mathrm{`}\mathrm{"}`$ denotes the symplectic orthogonal with respect to $`\omega _Q`$. Observe that, locally, this condition means that the matrix $`(C^{ab})=\left({\displaystyle \frac{\mathrm{\Psi }^a}{p_i}}_{ij}{\displaystyle \frac{\mathrm{\Psi }^b}{p_j}}\right)`$ is regular. The compatibility condition is not too restrictive, since it is trivially verified by the usual systems of mechanical type (i.e. with a Lagrangian of the form kinetic minus potential energy), where the $`_{ij}`$ represent the components of a positive definite Riemannian metric. The compatibility condition guarantees, in particular, the existence of a unique solution of the constrained equations of motion (2) which, henceforth, will be denoted by $`X_{H,M}`$ on the Hamiltonian side and $`Leg_{}^1(X_{H,M})=\xi _{L,D}`$ on the Lagrangian side.

Moreover, if we denote by $`X_H`$ the Hamiltonian function of $`H`$, i.e., $`i_{X_H}\omega _Q=dH`$ then, using the constraint functions, we may explicitely determine the Lagrange multipliers $`\lambda _a`$ as $`\lambda _a=𝒞_{ab}X_H(\mathrm{\Psi }^b).`$ Next, writing the 1-form $`\mathrm{\Lambda }=𝒞_{ab}X_H(\mathrm{\Psi }^b)\frac{\mathrm{\Psi }^a}{p_j}_{ji}dq^i`$ then, the nonholonomic equations are equivalently rewritten as

$$\{\begin{array}{ccc}\hfill \dot{q}^i& =& \frac{H}{p_i},\hfill \\ \hfill \dot{p}_i& =& \frac{H}{q^i}\mathrm{\Lambda }_i,\hfill \end{array}$$
(3)
for initial conditions $`(q_0,p_0)M`$ and $`\mathrm{\Lambda }=\mathrm{\Lambda }_idq^i`$. We also denote by $`\stackrel{~}{\mathrm{\Lambda }}=Leg^{}(\mathrm{\Lambda })`$ the 1-form on $`TQ`$ wich represents the constraint force once the Lagrange multipliers have been determined.

Now, consider the flow $`F_t:MM`$, $`tI\text{}`$ of the vector field $`X_{H,M}`$, solution of the nonholonomic problem. Since (3) is geometrically rewritten as

$$i_{X_{H,M}}\omega _Q=dH+\mathrm{\Lambda },$$

then

$$L_{X_{H,M}}\theta _Q=d(i_{X_{H,M}}\theta _QH)\mathrm{\Lambda },$$

or, equivalently, $`L_{X_{H,M}}\theta _Q=d(LLeg^1)\mathrm{\Lambda }.`$ Therefore, integrating

$$F_h^{}\theta _Q\theta _Q=d\left(_0^hL\stackrel{~}{F}_t𝑑t\right)_0^hF_t^{}\mathrm{\Lambda },$$
(4)
where $`\stackrel{~}{F}_t`$ is the flow of the vector field $`\xi _{L,D}`$.

### 2.2 “Generating functions” and nonholonomic mechanics

In what follows, we will follow similar arguments for the construction of generating functions for symplectic or canonical maps \[Arn:78\]. However, because of equation (4), we have that the nonholonomic flow is not a canonical transformation; i.e.,

$$F_h^{}\omega _Q\omega _Q=d\left(_0^hF_t^{}\mathrm{\Lambda }\right).$$
(5)
This description will allow us to construct a new family of nonholonomic integrators for equations (19). Denote by $`\pi _i:T^{}Q\times T^{}QT^{}Q`$, $`i=1,2`$, the canonical projections. Consider the following forms

$`\mathrm{\Theta }`$ $`=`$ $`\pi _2^{}\theta _Q\pi _1^{}\theta _Q,`$
$`\mathrm{\Omega }`$ $`=`$ $`\pi _2^{}\omega _Q\pi _1^{}\omega _Q=d\mathrm{\Theta }.`$
Denote by $`i_{F_h}:\text{Graph}(F_h)T^{}Q\times T^{}Q`$ the inclusion map and observe that $`\text{Graph}(F_h)M\times M`$. Then, from (4) $`i_{F_h}^{}\mathrm{\Theta }`$ is equal to

$$(\pi _{1}^{}{}_{|\text{Graph}(F_h)}{}^{})^{}\left[d\left(_0^hL\stackrel{~}{F}_t𝑑t\right)_0^hF_t^{}\mathrm{\Lambda }\right].$$

Let $`(q_0,p_0,q_1,p_1)`$ be coordinates in $`T^{}Q\times T^{}Q`$ in a neighborhood of some point in $`\text{Graph}(F_h)`$. If $`(q_0,p_0,q_1,p_1)\text{Graph}(F_h)`$ then $`\mathrm{\Psi }^a(q_0,p_0)=0`$ and $`\mathrm{\Psi }^a(q_1,p_1)=0`$. Moreover, along $`\text{Graph}(F_h)`$, $`q_1=q_1(q_0,p_0)`$, $`p_1=p_1(q_0,p_0)`$ and

$`p_1dq_1p_0dq_0`$ $`=`$ $`d\left({\displaystyle _0^h}L(q(t),\dot{q}(t))𝑑t\right)`$ (6)
$`{\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)),`$
where $`(q(t),\dot{q}(t))=\stackrel{~}{F}_t(q_0,\dot{q}_0)`$ with $`Leg(q_0,\dot{q}_0)=(q_0,p_0)`$. Here, $`\stackrel{~}{F}_t`$ denotes the flow of $`\xi _{L,D}`$. Equation (6) is satisfied along $`\text{Graph}(F_h)`$.

Assume that, in a neighborhood of some point $`x\text{Graph}(F_h)`$, we can change this system of coordinates to a new coordinates $`(q_0,q_1)`$. Denote by

$$S^h(q_0,q_1)=_0^hL(q(t),\dot{q}(t))𝑑t,$$

where $`q(t)`$ is a solution curve of the nonholonomic problem with $`q(0)=q`$ and $`q(h)=q_1`$ and an adequate extension of $`S^h`$. It is easy to show that this solution always exists for adequate values of $`q_0`$ and $`q_1`$.

Thus, we deduce that

$$\{\begin{array}{c}p_0=\frac{S^h}{q_0}+_0^h\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t))\frac{q}{q_0},\hfill \\ p_1=\frac{S^h}{q_1}_0^h\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t))\frac{q}{q_1},\hfill \end{array}$$
(7)
where $`(q_0,q_1)`$ verifies the constraint functions $`\phi ^a(q_0,q_1,h)=0`$, explicitely defined by

$`\phi ^a(q_0,q_1,h)=`$
$`\mathrm{\Psi }^a(q_0,{\displaystyle \frac{S^h}{q_0}}(q_0,q_1)+{\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)){\displaystyle \frac{q}{q_0}}),`$ (8)
where $`q(t)`$ is a solution of the nonholonomic problem with $`q(0)=q_0`$ and $`q(h)=q_h`$.

Next, we will show how the group composite law of the flow $`F_h`$, $`F_{Nh}=\underset{N}{\underset{}{F_h\mathrm{}F_h}}`$, is expressed in terms of the corresponding “generating functions” $`S^h`$. Moreover, the following Theorem will result in a new construction of numerical integrators for nonholonomic mechanics when we change the “generating function” and the constraint forces by appropriate approximations.

###### Theorem 2.1

The function $`S^{Nh}`$, the “generating function” for $`F_{Nh}`$, is given by

$$S^{Nh}(q_0,q_N)=\underset{k=0}{\overset{N1}{}}S^h(q_k,q_{k+1}),$$

where $`q_k`$, $`1kN1`$, are points verifying

$`D_2S^h(q_{k1},q_k)+D_1S^h(q_k,q_{k+1})=`$ (9)
$`{\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)){\displaystyle \frac{q}{q_1}}+{\displaystyle _h^{2h}}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)){\displaystyle \frac{q}{q_0}},`$
and $`q(t)`$ is a solution curve of the nonholonomic problem with $`q(0)=q_{k1}`$ and $`q(h)=q_k`$ (respectively, $`q(h)=q_k`$ and $`q(2h)=q_{k+1}`$) for the first integral (resp., second integral) of the right-hand side.

Proof: It is suffices to prove the result for $`N=2`$; that is,

$$S^{2h}(q_0,q_2)=S^h(q_0,q_1)+S^h(q_1,q_2),$$

where $`q_1`$ verifies condition (9).

Since

$`p_1dq_1p_0dq_0`$ $`=`$ $`dS^h(q_0,q_1){\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)),`$
$`p_2dq_2p_1dq_1`$ $`=`$ $`dS^h(q_1,q_2){\displaystyle _h^{2h}}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)),`$
then

$`p_2dq_2p_0dq_0=d\left(S^h(q_0,q_1)+S^h(q_1,q_2)\right)`$
$`{\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)){\displaystyle _h^{2h}}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)).`$
Since the variables $`q_1`$ do not appear on the left-hand side term, we obtain expression (9). Moreover, for this choice of $`q_1`$ then $`S^{2h}(q_0,q_2)=S^h(q_0,q_1)+S^h(q_1,q_2)`$ is a “generating function of the first kind” of $`F_{2h}`$.

Equations (9) determine an implicit system of difference equations which permit us to obtain $`q_2`$ from the initial data $`q_0`$ and $`q_1`$.

### 2.3  Nonholonomic integrators

In the sequel and, for simplicity, assume that $`Q`$ is a vector space. Since $`S^h(q_0,q_1)=_0^hL(q(t),\dot{q}(t))𝑑t`$, where $`q(t)`$ is a nonholonomic solution with $`q(0)=q_0`$ and $`q(h)=q_1`$, we can obtain nonholonomic integrators by taking adequate approximations of the “generating function” $`S^h`$ and the extra-term $`_0^h\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t))`$.

Consider, for instance, the approximation

$$S_\alpha ^h(q_0,q_1)=hL((1\alpha )q_0+\alpha q_1,\frac{q_1q_0}{h}),$$
(10)
for some parameter $`\alpha [0,1]`$. (In general, we will write $`S_\alpha ^h(q_0,q_1)S^h(q_0,q_1)`$.)

A natural approximation of the constraint forces adapted to our choice of approximation for $`S^h`$ are

$`{\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)){\displaystyle \frac{q}{q_0}}`$
$`(1\alpha )h\stackrel{~}{\mathrm{\Lambda }}((1\alpha )q_0+\alpha q_1,{\displaystyle \frac{q_1q_0}{h}}),`$
$`{\displaystyle _0^h}\stackrel{~}{\mathrm{\Lambda }}(q(t),\dot{q}(t)){\displaystyle \frac{q}{q_1}}\alpha h\stackrel{~}{\mathrm{\Lambda }}((1\alpha )q_0+\alpha q_1,{\displaystyle \frac{q_1q_0}{h}}).`$
Consequently, we obtain the following numerical method for nonholonomic systems

$`D_2S_\alpha ^h(q_{k1},q_k)+D_1S_\alpha ^h(q_k,q_{k+1})=`$
$`\alpha h\stackrel{~}{\mathrm{\Lambda }}((1\alpha )q_{k1}+\alpha q_k,{\displaystyle \frac{q_kq_{k1}}{h}})`$
$`+(1\alpha )h\stackrel{~}{\mathrm{\Lambda }}((1\alpha )q_k+\alpha q_{k+1},{\displaystyle \frac{q_{k+1}q_k}{h}}),`$
with $`1kN1`$ and initial condition satisfying

$`\stackrel{~}{\phi }^a`$ $`(`$ $`q_0,q_1,h)=\mathrm{\Psi }^a(q_0,{\displaystyle \frac{S_\alpha ^h}{q_0}}(q_0,q_1)`$
$`+(1\alpha )h\stackrel{~}{\mathrm{\Lambda }}((1\alpha )q_0+\alpha q_1,{\displaystyle \frac{q_1q_0}{h}}))=0.`$
###### Example 2.2

Nonholonomic particle.

Consider the Lagrangian $`L:T\text{}^3\text{}`$

$$L=\frac{1}{2}(\dot{x}^2+\dot{y}^2+\dot{z}^2)(x^2+y^2),$$

subject to the constraint $`\varphi =\dot{z}y\dot{x}=0`$. Taking $`\alpha =1/2`$ in (10) we obtain a geometric integrator for the continuous nonholonomic problem. The first figure compares the method introduced here to the traditional Runge-Kutta method of fourth order, showing an improvement in several orders of magnitude. Observe that, in this scale, the value of the energy in each step of our algorithm is practically undistinguishable from the initial value of the energy.

The second figure is a comparison between our method and the one proposed in \[CorMar:01\]. A similar behaviour is observed. Nevertheless, a slightly better behaviour can also be appreciated, where the proposed algorithm shows on average a better preservation of the original energy.

## 3 Optimal control theory

### 3.1 Geometric formulation of optimal control problems

A general optimal control problem consists of a set of differential equations

$$\dot{q}^i=\mathrm{\Gamma }^i(q(t),u(t)),\mathrm{\hspace{0.33em}1}in,$$
(11)
where $`q^i`$ denote the states and $`u`$ the control variables, and a cost function $`L(q,u)`$. Given some boundary conditions (usually $`q_0=q(t_0)`$ and $`q_F=q(t_f)`$) the aim is to find a $`C^2`$-piecewise smooth curve $`c(t)=(q(t),u(t))`$, satisfying the control equations (11) and minimizing the functional

$$𝒥(c)=_{t_0}^{t_f}L(q(t),u(t))𝑑t.$$
(12)
In a global description, one assumes a fiber bundle structure $`\pi :UQ`$, where $`Q`$ is the configuration manifold with local coordinates $`q^i`$ and $`U`$ is the bundle of controls, with local coordinates $`(q^i,u^a)`$, $`1in`$, $`1am`$.

The ordinary differential equations (11) on $`Q`$ depending on the parameters $`u`$ can be seen as a vector field $`\mathrm{\Gamma }`$ along the projection map $`\pi `$, that is, $`\mathrm{\Gamma }`$ is a smooth map $`\mathrm{\Gamma }:UTB`$ such that the diagram

is commutative. This vector field is locally written as $`\mathrm{\Gamma }=\mathrm{\Gamma }^i(q,u){\displaystyle \frac{}{q^i}}`$.

The solutions of such problem are provided by Pontryaguin’s maximum principle. If we construct the Hamiltonian function

$$H(q,p,u)=L(q,u)+p_i\mathrm{\Gamma }^i(q,u)$$
(13)
where $`p_i`$, $`1in`$, are now considered as Lagrange’s multipliers, then a curve $`\gamma :\text{}U`$, $`\gamma (t)=(q(t),u(t))`$ is an optimal trajectory if there exists functions $`p_i(t)`$, $`1in`$ such that they are solutions of the Hamilton equations:

$$\{\begin{array}{c}\dot{q}^i(t)=\frac{H}{p_i}(q(t),p(t),u(t))\hfill \\ \dot{p}_i(t)=\frac{H}{q^i}(q(t),p(t),u(t))\hfill \end{array}$$
(14)
and

$$H(q(t),p(t),u(t))=\underset{𝑣}{\text{min}}H(q(t),p(t),v),$$
(15)
with $`t[t_0,t_f]`$. This last condition is usually replaced by

$$\frac{H}{u^a}=0,1am$$
(16)
when we are looking for extremal trajectories.

It is well known that the Pontryaguin’s necessary conditions for extremality have a geometric interpretation in terms of presymplectic hamiltonian system. The total space of the system will be $`T^{}Q\times _QU`$. Let $`\omega _Q`$ be the canonical symplectic form on $`T^{}Q`$ and consider the canonical projection $`\text{pr}_1:T^{}Q\times _QUT^{}Q`$. Denote by $`\omega =\text{pr}_1^{}\omega _Q`$ the induced closed 2-form on $`T^{}Q\times _QU`$. The 2-form $`\omega `$ is degenerate and its characteristic distribution is locally spanned by $`/u^a`$, $`1am`$. Define the Pontryaguin’s hamiltonian function $`H:T^{}Q\times _QU\text{}`$ as follows $`H(\alpha _q,u_q)=L(u_q)+\alpha _q(\mathrm{\Gamma }(u_q))`$ where $`\alpha _qT_q^{}Q`$ and $`u_q\text{pr}^1(q)`$. Obviously, the coordinate expression of $`H`$ is (13).

Equations (14) (15) and (16) are intrinsically written as

$$i_X\omega =dH$$
(17)
Applying the Dirac-Bergmann-Gotay-Nester algorithm to the presymplectic system $`(T^{}Q\times _QU,\mathrm{\Omega },H)`$ we obtain that equations (16) correspond to the primary constraints for the presymplectic system: $`\varphi ^a=\frac{H}{u^a}=0`$. The equations have solution along the first constraint submanifold $`P_0`$ determined by the vanishing of the primary constraints. On the points of $`P_0`$ there is at least a pointwise solution of Equation (17) but such solutions are not, in general, tangent to $`P_0`$. These points must be removed leaving a subset $`P_1P_0`$ (it is assumed tan $`P_1`$ is also a submanifold). Then, we have to restrict $`P_1`$ to a submanifold where the solutions of (17) are tangent to $`P_1`$. Proceeding further, we obtain a sequence of submanifolds

$$\mathrm{}P_k\mathrm{}P_2P_1P_0T^{}Q\times _QU$$

If this algorithm stabilizes, i.e. there exists a positive integer $`k\text{}`$ such that $`P_k=P_{k+1}`$ and $`dimP_k0`$, then we will obtain an stable submanifold $`P_f=P_k`$, on which a vector field exists such that

$$\left(i_X\omega =dH\right)_{|P_f}$$
(18)
The constraints determining $`P_f`$ are known in the control literature as higher order conditions for optimality. Therefore, a necessary condition for optimality of the curve $`\gamma :\text{}U`$, $`\gamma (t)=(q(t),u(t))`$ will be the existence of a lift $`\stackrel{~}{\gamma }`$ of $`\gamma `$ to $`P_f`$ such that $`\stackrel{~}{\gamma }`$ will be an integral curve of a solution of Equations (18).

In the regular case, the final constraint algorithm is $`P_0`$ (that is, $`P_0=P_f`$) and all the constraints are second class following the classical classification of Dirac. In such case $`(P_0,\omega _0)`$ is a symplectic manifold, where $`\mathrm{\Omega }_0`$ denotes the restriction of the presymplectic 2-form to the constraint submanifold $`P_0`$. Locally, the symplecticity of $`(P_0,\omega _0)`$ is equivalent to the regularity of the matrix $`\left({\displaystyle \frac{^2H}{u^au^b}}\right)_{1a,bm}`$. The dynamical equations for the optimal control problem will be

$$i_{X_{P_0}}\omega _0=dH_{|P_0}$$
(19)
Taking coordinates $`(q^i,p_i)`$ on $`P_0`$, then the dynamical equations are:

$$\{\begin{array}{c}\dot{q}^i(t)=\frac{H_{|P_0}}{p_i}(q(t),p(t))\hfill \\ \dot{p}_i(t)=\frac{H_{|P_0}}{q^i}(q(t),p(t))\hfill \end{array}$$
(20)
where we have substituted in (14) the control variables $`u^a`$ by its value $`\overline{u}^a=f^a(q,p)`$ applying the implicit function theorem to the primary constraints $`\varphi ^a=0`$. In such case, there exists a unique solution $`X_{P_0}`$ of Equation (19) and its flow preserves the symplectic 2-form $`\omega _0`$, i.e. it is a canonical transformation.

### 3.2 Generating functions of the second kind

Let $`(,\mathrm{\Omega })`$ be an exact symplectic manifold ($`\mathrm{\Omega }`$ is symplectic and exact, $`\mathrm{\Omega }=d\mathrm{\Theta }`$) and suppose that $`F:`$ is a transformation from $``$ to itself and $`\text{Graph}(F)`$ the graph of $`F`$, $`\text{Graph}(F)\times `$. Denote by $`\pi _i:\times `$, $`i=1,2`$ the canonical projections and the forms:

$`\overline{\mathrm{\Theta }}`$ $`=`$ $`\pi _2^{}\mathrm{\Theta }\pi _1^{}\mathrm{\Theta }`$
$`\overline{\mathrm{\Omega }}`$ $`=`$ $`\pi _2^{}\mathrm{\Omega }\pi _1^{}\mathrm{\Omega }=d\overline{\mathrm{\Theta }}`$
Denote by $`i_F:\text{Graph}(F)\times `$ the inclusion map. Then, $`F`$ is a canonical transformation if and only if $`i_F^{}\overline{\mathrm{\Omega }}=0`$, that is, if $`\text{Graph}(F)`$ is a lagrangian submanifold of $`(\times ,\overline{\mathrm{\Omega }})`$. In such a case, $`i_F^{}\overline{\mathrm{\Omega }}=di_F^{}\overline{\mathrm{\Theta }}=0`$ and, at least locally, there exists a function $`S:\text{Graph}F\text{}`$ such that

$$i_F^{}\overline{\mathrm{\Theta }}=dS$$
(21)
Taking $`(q^i,p_i)`$ as natural coordinates in $`\text{Graph}(F)`$ and $`(q^i,p_i,𝐪^i,𝐩_i)`$ the coordinates in $`\times `$, then, along $`\text{Graph}(F)`$, $`𝐪^i=𝐪^i(q,p)`$ and $`𝐩_i=𝐩_i(q,p)`$ and $`𝐩_id𝐪^ip_idq^i=dS(q,p).`$ Suppose that $`(q^i,𝐩_i)`$ are independent local coordinates on $`\text{Graph}(F)`$ (see \[Arn:78\]); i.e. $`S=S(q,𝐩)`$ Since

$$𝐩_id𝐪^ip_idq^i=𝐪^id𝐩_i+d(𝐪^i𝐩_i)p_idq^i=dS,$$

if we define $`S_2(q,𝐩)=𝐪^i𝐩_iS(q,𝐩)`$, where $`𝐩`$ is expressed in terms of $`p`$ and $`𝐪`$, then $`𝐪^id𝐩_i+p_idq^i=dS_2(q,𝐩)`$

###### Definition 3.1

The function $`S_2(q,𝐩)`$ will be called a generating function of the second kind of the canonical transformation $`F`$.

Now, suppose that $`(,\mathrm{\Omega },H)`$ is a hamiltonian system and $`X_H`$ its hamiltonian vector field, say $`i_{X_H}\mathrm{\Omega }=dH`$. Denote by $`F_h:`$ its flow.

###### Theorem 3.2

Let a function $`S_2^{Nh}`$ be defined by

$$S_2^{Nh}(q_0,p_{Nh})=\underset{k=0}{\overset{N1}{}}(S_2^h(q_k,p_{k+1})q_{k+1}p_{k+1})$$

where $`q_k`$, $`1kN`$, and $`p_k`$, $`0kN1`$, are stationary points of the right-hand side, that is

$`q_{k+1}`$ $`=`$ $`{\displaystyle \frac{S_2^h}{p}}(q_k,p_{k+1}),0kN1`$
$`p_k`$ $`=`$ $`{\displaystyle \frac{S_2^h}{q}}(q_k,p_{k+1}),0kN1`$
then $`S_2^{Nh}`$ is a generating function of the second kind for $`F_{Nh}:`$.

Proof: It is similar to that of Theorem 2.1.

Finally, we have the following

###### Proposition 3.3

A generating function of the second kind for $`F_h`$ is given by

$$S_2^h(q_0,p_h)=p_hq_h_0^h\left(pdqHdt\right)$$

where $`t(q(t),p(t))`$ is an integral curve of the Hamilton equations such that $`q(0)=q_0`$ and $`p(h)=p_h`$.

### 3.3 Generating functions of the second kind and discrete optimal control problems

From Proposition 3.3 a generating function of the second kind for the Hamiltonian system $`(P_0,\mathrm{\Omega }_0,H_{|P_0})`$ which determines the dynamics of the optimal control problem given by (11) and (12) is

$`S_2^h(q_0,p_h)=p_hq_h`$
$`{\displaystyle _0^h}\left(p(t)\dot{q}(t)H_{|P_0}(q(t),p(t))\right)𝑑t`$ (22)
where $`t(q(t),p(t))`$ is an integral curve of the vector field $`X_{P_0}`$ with $`(q(0),p(0))=(q_0,p_0)`$ and $`(q(h),p(h))=(q_h,p_h)`$.

We now turn to the construction of a numerical integrator for the Hamiltonian system $`(P_0,\omega _0,H_{|P_0})`$ by using an approximation of the generating function. The proposed methods also realize the integration steps by canonical transformations; therefore, they are symplectic integrators.

###### Example 3.4

Consider, for instance, the following approximation to $`S_2^h`$:

$`\stackrel{~}{S}_2^h(q_k,p_{k+1})=p_{k+1}q_{k+1}hp_{k+1}\left({\displaystyle \frac{q_{k+1}q_k}{h}}\right)`$
$`+h\stackrel{~}{L}_d(q_k,p_{k+1})+hp_{k+1}\stackrel{~}{\mathrm{\Gamma }}_d(q_k,p_{k+1})`$
where $`\stackrel{~}{L}_d`$ and $`\stackrel{~}{\mathrm{\Gamma }}_d`$ are adequate approximations to $`L_{|P_0}`$ and $`\mathrm{\Gamma }_{|P_0}`$, respectively.

Denote by $`\stackrel{~}{f}(q_k,p_{k+1})`$ the function $`\stackrel{~}{f}(q_k,p_{k+1})=h\mathrm{\Gamma }_d(q_k,p_{k+1})+q_k`$. Since $`{\displaystyle \frac{q_{k+1}q_k}{h}}=\stackrel{~}{\mathrm{\Gamma }}_d(q_k,p_{k+1})`$ then,

$$\stackrel{~}{S}_2^h(q_k,p_{k+1})=\stackrel{~}{L}_d(q_k,p_{k+1})+p_{k+1}\stackrel{~}{f}(q_k,p_{k+1})$$

and hence the equations

$$\{\begin{array}{c}p_k=\frac{\stackrel{~}{S}_d^h}{q}(q_k,p_{k+1})\hfill \\ q_{k+1}=\frac{\stackrel{~}{S}_d^h}{p}(q_k,p_{k+1})\hfill \end{array}$$
(23)
are exactly the discrete equations corresponding to the classical discrete optimal control problem (see \[Lew:86\]), determined by the control equations: $`q_{k+1}^i=\stackrel{~}{f}^i(q_k,u_k),\text{(}(q_0)\text{ given)}`$ and with associate perfomance index: $`J=_{k=0}^{N1}\stackrel{~}{L}_d(q_k,u_k)`$ Observe that this discrete optimal control problem is symplectic in the sense explained in the subsection above.


text:
# On the evaluation of the improvement parameter in the lattice Hamiltonian approach to critical phenomena

## I Introduction

The quest for better analytical and numerical methods in the theoretical evaluation of measurable physical quantities, like critical exponents and amplitude ratios, is one of the lasting tasks of statistical field theory.

In recent years substantial progress in this field has been made by the introduction of a method based on the strong coupling lattice expansion of improved Hamiltonians CPRV-impr-O1  (for a review, cf. Ref. PV-rept ). The essential feature of this method is the possibility of removing all leading nonscaling contributions to physical quantities, in the neighborhood of criticality, by a specific choice of a parameter in the lattice Hamiltonian (critical coupling). The convergence of analytical and/or numerical evaluations is therefore impressively faster than in any other variant form of the models under inspection belonging to the same universality class.

The main limitation of this method lays in the absence of an efficient analytical technique for the determination of the critical parameter. Conceptual reasons for this limitation may be found in the impossibility of an $`ϵ`$ expansion for the value of the critical parameter, both in the $`4ϵ`$ and in the $`2+ϵ`$ expansion schemes. In practice one must resort to an extrapolation from numerical Monte Carlo finite-size evaluations of some physical quantity, typically the Binder cumulant CPRV-impr-O1 ; PV-rept .

Another well-known analytical approach to the study of critical lattice models is the $`1/N`$ expansion, which applies in particular to the physically very important class of $`O(N)`$ spin models in three dimensions. However it is known that, in exactly three dimensions and for nearest-neighbor interactions, the critical parameter in the large-$`N`$ limit is a negative number ZJ . Since the critical parameter controls the large-field behavior of the interaction potential, a negative value would naively imply an unbounded Hamiltonian. In practice this would prevent a Monte Carlo simulation of the system, thus seriously jeopardizing the conceptual meaning of the whole approach.

Nevertheless we still seriously believe the $`1/N`$ expansion to be, from a theoretical point of view, probably the most relevant expansion scheme that can be applied to any quantum and statistical field theory, in that there is no known obstruction to summability of the series expansion in powers of $`1/N`$ for the values of physical quantities.

Therefore we decided to explore the conceptual and numerical consequences of performing a systematic $`1/N`$ expansion of the critical parameter, for the class of three-dimensional $`O(N)`$ spin models, in order to check the actual relevance of the drawbacks that we mentioned above.

As a consequence of our analysis we found that the critical parameter can be formally computed within the expansion with no limitation related to the sign of its large-$`N`$ value and for space dimensionalities in the interval $`2<d<4`$. In particular we found that the sign of the first $`1/N`$ correction is positive, and one may then expect to find a value $`N_c`$ such that the parameter itself vanishes. It is however still unclear that the predictions of the $`1/N`$ expansion may be extended to the region $`N<N_c`$.

As a byproduct of our analysis we obtained new more efficient expressions for the asymptotic expansions of many important functions entering our calculations. We presented these results with some details because they might be relevant to other computations of critical and subcritical quantities.

In Sec. II we introduce the lattice $`O(N)`$ models and their $`1/N`$ expansion. In Secs. III, IV, and V we discuss on general grounds the relevant asymptotic epansions (gap equation, effective propagator and renomalized coupling). In Secs. VI and VII we specialize these expansions to the case of the standard nearest-neighbor interaction, with the help some useful integral representations. A few numerical results are presented in Sec. VIII. In Sec. IX we compute the $`1/N`$ correction to the improvement parameter and finally in Sec. X we discuss the meaning and relevance of our results.

## II The effective Hamiltonian and the graph expansion

Our starting point will be the usual $`N`$-component $`\varphi ^4`$ lattice Hamiltonian in $`d`$ dimensions:

$$H=\underset{x}{}\left[\frac{1}{2}_\mu _\mu \mathit{\varphi }(x)_\mu \mathit{\varphi }(x)+\frac{1}{2}\mu _0^2\mathit{\varphi }^2(x)+\frac{1}{4!}g_0\left(\mathit{\varphi }^2(x)\right)^2\right],$$
(1)
where $`_\mu \mathit{\varphi }(x)`$ is some (local) form of the lattice gradient; in the standard nearest-neighbor formulation $`_\mu \mathit{\varphi }(x)=\mathit{\varphi }(x+\mu )\mathit{\varphi }(x)`$.

Following Ref. CPRV-ON  we define the rescaled couplings

$$\beta =\frac{6\mu _0^2}{g_0N},\gamma =\frac{3}{g_0N}$$

and introduce an auxiliary field $`\alpha `$ in order to eliminate the quartic term in the Hamiltonian.

After a trivial Gaussian integration the resulting effective Hamiltonian is

$$H_{\mathrm{eff}}=\frac{N}{2}[TrLn\beta (_\mu _\mu +i\alpha )i\beta \alpha +\gamma \alpha ^2].$$
(2)
In the limit $`\gamma 0`$ this Hamiltonian reduces to the usual effective large-$`N`$ expression for the nonlinear $`\sigma `$ model. In the nearest-neighbor formulation, it is also known as the $`O(N)`$ Heisenberg model; its large-$`N`$ limit was investigated in Ref. MR .

The saddle-point condition on the effective Hamiltonian leads to the so-called gap equation

$$\beta +2\gamma m_0^2=_\pi ^\pi \frac{d^dp}{(2\pi )^d}\frac{1}{\overline{p}^2+m_0^2},$$
(3)
where $`\overline{p}^2`$ is the Fourier transform of the lattice Laplacian operator $`_\mu _\mu `$, which in the nearest-neighbor case takes the form $`\widehat{p}^2=2_\mu (1\mathrm{cos}p_\mu )`$.

The gap equation allows for the elimination of $`\beta `$ in favor of the new parameter $`m_0`$ (large-$`N`$ inverse correlation length) in the Feynman graph expansion.

In the large-$`N`$ limit criticality corresponds to the vanishing of $`m_0^2`$, and the criticality condition may then take the form

$$\beta _c=_\pi ^\pi \frac{d^dp}{(2\pi )^d}\frac{1}{\overline{p}^2},$$

corresponding to a finite value of $`\beta _c`$ for all $`d>2`$.

The graph expansion for this model, in the formulation based on the effective Hamiltonian, requires defining the (bare) propagator $`\mathrm{\Delta }`$ for the effective field $`\alpha `$; by standard manipulations one obtains

$$D(k,m_0,\gamma )\mathrm{\Delta }^1(k,m_0,\gamma )=\frac{1}{2}_\pi ^\pi \frac{d^dp}{(2\pi )^d}\frac{1}{\overline{p}^2+m_0^2}\frac{1}{(\overline{p+k})^2+m_0^2}+\gamma .$$
(4)
The above approach is quite general, and it leads to a systematic $`1/N`$ expansion of the correlation functions and of the physical quantities for arbitrary values of $`m_0`$ and $`\gamma `$.

However we want to focus our attention on the critical domain, and in particular we want to evaluate the coupling $`\gamma ^{}`$ such that the first nontrivial corrections to scaling turn out to vanish in the computation of physical quantities in the scaling region.

To this purpose it is convenient to parametrize the cutoff dependence of correlation functions and renormalized couplings in terms of the lattice spacing $`a`$, which can be made to appear explicitly in calculations by a rescaling of the coupling and momentum dependence.

Around criticality the dependence on $`a`$ is not analytic, and as a consequence we need asymptotic expansions in order to identify the scaling and leading nonscaling contributions to any computable quantity.

The basic technique for asymptotic expansions in powers of $`m_0a`$ is described in Ref. CR ; here we shall discuss its applications to the cases of interest for the present paper.

We only recall that, in order to regularize the generic lattice integral

$$I(k;m_0a)_\pi ^\pi \frac{d^dp}{(2\pi )^d}F(k;m_0a,p),$$

where $`k`$ is any collection of external momenta, we can make use of the formal identity

$$I(k;m_0a)=I_{\mathrm{lat}}(k;m_0a)+I_{\mathrm{con}}(k;m_0a),$$

where

$$I_{\mathrm{lat}}(k;m_0a)=_\pi ^\pi \frac{d^dp}{(2\pi )^d}F(k;m_0a,p)_{\mathrm{}}^+\mathrm{}a^d\frac{d^dp}{(2\pi )^d}T^{(\mathrm{IR})}F(k;m_0a,pa)$$

and

$$I_{\mathrm{con}}(k;m_0a)=a^d_{\mathrm{}}^+\mathrm{}\frac{d^dp}{(2\pi )^d}[F(k;m_0a,pa)T^{(\mathrm{UV})}F(k;m_0a,pa)].$$

The $`TT^{(\mathrm{IR})}+T^{(\mathrm{UV})}`$ operation amounts to a Taylor series expansion of the integrand in powers of $`m_0a`$, $`T^{(\mathrm{IR})}`$ and $`T^{(\mathrm{UV})}`$ corresponding respectively to the IR and UV singular terms in the expansion.

It is possible to prove that the expansions of $`I_{\mathrm{lat}}`$ (the “lattice contribution”) and $`I_{\mathrm{con}}`$ (the “continuum contribution”) are individually and fully regular; the nonanaliticity of the expansion is factored out in the $`a^d`$ term multiplying the continuum contribution.

## III Asymptotic expansion of the gap equation

For the purposes of the present paper and in order to show an explicit example of the asymptotic expansion procedure let us consider the gap equation in the large-$`N`$ limit.

Let us assume a generic lattice Laplacian such that $`\overline{p}^2p^2+cp^4+O(p^6)`$, where $`p^{2n}=_\mu p_\mu ^{2n}`$: in the nearest-neighbor version $`c=\frac{1}{12}`$.

From the previously derived results we obtain the following relationship:

$$\beta _c\beta =m_0^2a^2\left[_\pi ^\pi \frac{d^dp}{(2\pi )^d}\frac{1}{\overline{p}^2(\overline{p}^2+m_0^2a^2)}+2\gamma \right].$$

Let us now define

$$G(p,m_0a)\frac{m_0^2a^2}{\overline{p}^2(\overline{p}^2+m_0^2a^2)}$$

and perform the relevant expansions up to the first few nontrivial terms:

$`TG(p,m_0a)`$ $``$ $`{\displaystyle \frac{m_0^2a^2}{(\overline{p}^2)^2}}{\displaystyle \frac{m_0^4a^4}{(\overline{p}^2)^3}}+O\left(m_0^6a^6\right),`$
$`G(pa,m_0a)`$ $``$ $`{\displaystyle \frac{1}{a^2}}\left[{\displaystyle \frac{1}{p^2}}{\displaystyle \frac{1}{p^2+m_0^2}}\right]cp^4\left[{\displaystyle \frac{1}{(p^2)^2}}{\displaystyle \frac{1}{(p^2+m_0^2)^2}}\right]+O\left(a^2\right),`$
$`TG(pa,m_0a)`$ $``$ $`{\displaystyle \frac{m_0^2}{a^2}}\left[{\displaystyle \frac{1}{(p^2)^2}}{\displaystyle \frac{m_0^2}{(p^2)^3}}\right]m_0^2cp^4\left[{\displaystyle \frac{2}{(p^2)^3}}{\displaystyle \frac{3m_0^2}{(p^2)^4}}\right]+O(m_0^6,a^2).`$
By grouping together the IR singular terms we therefore obtain the lattice contributions:

$$2m_0^2a^2(\gamma \gamma _0)m_0^4a^4\delta _0+O\left(m_0^6a^6\right),$$

where we defined the following numerical constants:

$`\gamma _0`$ $``$ $`{\displaystyle \frac{1}{2}}\left[{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}+{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(p^2)^2}}\right],`$
$`\delta _0`$ $``$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(\overline{p}^2)^3}}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d^dp}{(2\pi )^d}}\left({\displaystyle \frac{1}{(p^2)^3}}3c{\displaystyle \frac{p^4}{(p^2)^4}}\right).`$ (5)
On the other side, by grouping together the UV singular terms, we obtain the following continuum contributions:

$$_{\mathrm{}}^+\mathrm{}\frac{d^dp}{(2\pi )^d}\left(a^{d2}\left[\frac{1}{p^2+m_0^2}+\frac{1}{p^2}\right]+a^dcp^4\left[\frac{1}{(p^2+m_0^2)^2}\frac{1}{(p^2)^2}2\frac{m_0^2}{(p^2)^3}\right]+O(a^{d+2})\right).$$

We can perform the continuum integrals by standard dimensional regularization techniques. We then finally find

$$\beta _c\beta b_0(m_0a)^{d2}+2(\gamma \gamma _0)m_0^2a^2\frac{3}{2}cb_0m_0^da^d\delta _0m_0^4a^4+O\left((m_0a)^{d+2}\right),$$
(6)
where

$$b_0\frac{\mathrm{\Gamma }(1\frac{d}{2})}{(4\pi )^{\frac{d}{2}}}.$$

From the asymptotic expansion of the gap equation we immediately learn the following lesson: it is possible to choose for the quartic coupling $`\gamma `$ the special value $`\gamma _0`$ such that the first nonscaling contribution to the large-$`N`$ saddle-point condition vanishes.

By this procedure we have identified the large-$`N`$ critical coupling for these versions of the model. In fact any change in the form of the local interaction, reflecting itself in the detailed $`p`$ dependence of the bare lattice massless propagator without changing its singular part, leads to a finite change in the numerical value of $`\gamma _0`$ without affecting its formal representation.

It is especially interesting to consider the class of models characterized by next-to-nearest neighbor interactions, and whose propagator of the fundamental excitations is obtained in the form

$$\overline{p}^2=\underset{\mu =1}{\overset{d}{}}\left(6c+\frac{5}{2}\right)8\left(c+\frac{1}{3}\right)\mathrm{cos}p_\mu +2\left(c+\frac{1}{12}\right)\mathrm{cos}2p_\mu =\widehat{p}^2+\left(c+\frac{1}{12}\right)\widehat{p}^4,$$

where

$$\widehat{p}^n\underset{\mu =1}{\overset{d}{}}\widehat{p}_\mu ^n,\widehat{p}_\mu 2\mathrm{sin}\frac{p_\mu }{2}.$$

Note that $`c=\frac{1}{12}`$ corresponds to the standard nearest-neighbor interaction, while $`c=0`$ corresponds to the $`O(a^2)`$ Symanzik tree-improved version of $`O(N)`$ models Symanzik .

In three dimensions we have numerically explored the range $`\frac{1}{12}c0`$: our results are presented in Figs. 1 and 2.

Let us notice in particular that the choice $`c=0.03332110\mathrm{}`$ corresponds to vanishing $`\gamma _0`$, and it is therefore, at least in the large-$`N`$ limit, an alternative version of a spin model where the leading corrections to scaling are automatically made to vanish. In turn when $`c=0`$ we obtain $`\gamma _0=0.003328210\mathrm{}`$, which implies a small but nonvanishing leading scaling violation, and $`\delta _0=0.002181406\mathrm{}`$ . Compared to the standard $`c=\frac{1}{12}`$ case this version has however the advantage of being numerically testable also by Monte Carlo methods, since $`\gamma _0>0`$.

## IV Asymptotic expansion of the effective propagator

In order to perform an asymptotic expansion of the $`O(1/N)`$ contributions to physical quantities we must compute one-loop graphs involving the effective propagator $`\mathrm{\Delta }(k,m_0,\gamma )`$.

We therefore need to evaluate the asymptotic expansion of $`\mathrm{\Delta }`$ or, more conveniently, of $`D(k,m_0,\gamma )`$. It is easy to recognize, from the definition of $`D`$ and from the general rule of the asymptotic expansion, that in general we may express the result in the form

$$D(k,m_0a,\gamma )=\underset{n=0}{\overset{\mathrm{}}{}}\left[A_n(k)(m_0a)^{2n}+B_n(k)(m_0a)^{2n+d2}\right],$$
(7)
where $`A_n(k)`$ have the form of lattice contributions (and only $`A_0`$ depends on $`\gamma `$), while $`B_n(k)`$ are continuum contributions that can be analytically computed, e.g., in dimensional regularization.

We shall not repeat the derivation (some details can be found in Ref. CR ), but only quote the relevant results:

$`A_0(k,\gamma )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{\overline{p}^2}}{\displaystyle \frac{1}{(\overline{p+k})^2}}+\gamma ,`$
$`A_1(k)`$ $`=`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}{\displaystyle \frac{1}{(\overline{p+k})^2}}+{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(p^2)^2}}{\displaystyle \frac{1}{\overline{k}^2}},`$
$`A_2(k)`$ $`=`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(\overline{p}^2)^3}}{\displaystyle \frac{1}{(\overline{p+k})^2}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}{\displaystyle \frac{1}{((\overline{p+k})^2)^2}}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(p^2)^3}}{\displaystyle \frac{1}{\overline{k}^2}}`$
$``$ $`3c{\displaystyle \frac{p^4}{(p^2)^4}}{\displaystyle \frac{1}{\overline{k}^2}}+{\displaystyle \frac{1}{(p^2)^2}}C_0(k)+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{(p^2)^2}}+{\displaystyle \frac{1}{((p+k)^2)^2}}\right){\displaystyle \frac{1}{(\overline{k}^2)^2}};`$
$`B_0(k)`$ $`=`$ $`{\displaystyle \frac{b_0}{\overline{k}^2}},B_1(k)=b_0\left({\displaystyle \frac{1}{(\overline{k}^2)^2}}+C_0(k){\displaystyle \frac{3}{2}}{\displaystyle \frac{c}{\overline{k}^2}}\right),`$ (8)
where we defined

$$C(k,m_0)\frac{1}{2d}_\mu \frac{^2}{k_\mu ^2}\frac{1}{\overline{k}^2+m_0^2},C_0(k)C(k,0)\frac{1}{2d}_\mu \frac{^2}{k_\mu ^2}\frac{1}{\overline{k}^2}.$$

Trivial manipulations allow to express $`A_n(k)`$ as pure lattice integrals:

$`A_1(k)`$ $`=`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}\left[{\displaystyle \frac{1}{\overline{k}^2}}{\displaystyle \frac{1}{(\overline{p+k})^2}}\right]+2{\displaystyle \frac{\gamma _0}{\overline{k}^2}},`$
$`A_2(k)`$ $`=`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}\left[{\displaystyle \frac{1}{(\overline{p}^2)^3}}{\displaystyle \frac{1}{(\overline{p+k})^2}}{\displaystyle \frac{1}{(\overline{p}^2)^3}}{\displaystyle \frac{1}{\overline{k}^2}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}C_0(k)\right]+{\displaystyle \frac{\delta _0}{\overline{k}^2}}2\gamma _0C_0(k)`$
$`+`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dp}{(2\pi )^d}}\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}{\displaystyle \frac{1}{((\overline{p+k})^2)^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{(\overline{p}^2)^2}}{\displaystyle \frac{1}{(\overline{k}^2)^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{((\overline{p+k})^2)^2}}{\displaystyle \frac{1}{(\overline{k}^2)^2}}\right]2{\displaystyle \frac{\gamma _0}{(\overline{k}^2)^2}}.`$
From our general considerations it should be by now clear that we shall also need a different expansion of $`D`$, homogeneous in powers of $`m`$ and $`k`$.

Without delving into the details, we find that

$$D(ka,m_0a,\gamma )\frac{1}{2}a^{d4}_{\mathrm{}}^+\mathrm{}\frac{d^dp}{(2\pi )^d}\frac{1}{p^2+m_0^2}\frac{1}{(p+k)^2+m_0^2}+(\gamma \gamma _0)+O(a^{d2}).$$

Notice that the integral can be analytically computed in all dimensions, and the final result is

$$D(ka,m_0a,\gamma )d_0\left(\frac{1}{4}k^2a^2+m_0^2a^2\right)^{\frac{d}{2}2}{}_{2}{}^{}F_{1}^{}(2\frac{d}{2},\frac{1}{2},\frac{3}{2},\frac{1}{\xi ^2})+\gamma \gamma _0D_c,$$
(10)
where $`d_0=\frac{1}{2}(1\frac{d}{2})b_0`$, $`{}_{2}{}^{}F_{1}^{}`$ is the Gauss hypergeometric function and $`\xi =\sqrt{1+4m_0^2/k^2}`$.

Let us in general denote by the label “$`c`$” the quantities occurring in the leading order in the homogeneous (continuum) expansion of $`D`$, and in particular $`\gamma _c\gamma \gamma _0`$.

Following Ref. DRG , we can exploit identities between hypergeometric functions to recast the above result into the form

$$D_c=a_0\xi ^{d3}(ka)^{d4}+\gamma _c+\frac{b_0}{k^2a^2\xi ^2}{}_{2}{}^{}F_{1}^{}(\frac{d1}{2},1,\frac{d}{2},1\frac{1}{\xi ^2})(m_0a)^{d2},$$
(11)
where

$$a_0\frac{1}{2}\frac{\mathrm{\Gamma }(\frac{d}{2}1)^2\mathrm{\Gamma }(2\frac{d}{2})}{(4\pi )^{\frac{d}{2}}\mathrm{\Gamma }(d2)}.$$

For $`d=3`$, Eq. (11) reduces to

$$D_c=\frac{1}{16ak}+\gamma _c\frac{1}{8\pi ak}arccot\frac{k}{2m_0}.$$

Eq. (11) is especially appropriate for the asymptotic expansion of $`D_c`$, because all the nonanalytic dependence on $`m_0a`$ is explicitly factored out in the last term. In particular (after a rescaling $`kak`$) we obtain the following behaviors:

$$A_{0c}=a_0(k^2)^{\frac{d}{2}2}+\gamma _c,A_{nc}=\frac{1}{n!}\frac{4^n\mathrm{\Gamma }(\frac{d1}{2})}{\mathrm{\Gamma }(\frac{d1}{2}n)}a_0(k^2)^{\frac{d}{2}2n}.$$

Finally notice also that the zero-momentum value of $`D`$ is related to the derivative of the gap equation with respect to the mass, and we can obtain the relationship

$$D_0D(0,m_0a,\gamma )d_0(m_0a)^{d4}+\gamma _c\frac{3d}{8}cb_0(m_0a)^{d2}\delta _0m_0^2a^2,$$
(12)
and as a consequence

$$D_{0c}d_0(m_0a)^{d4}+\gamma _c.$$

## V Asymptotic expansion of the renormalized coupling

We now recall from the literature the expression of the $`O(1/N)`$ contribution to the (unrenormalized) self-energy of the fundamental quanta:

$$\mathrm{\Sigma }_1(p,m_0)=\mathrm{\Sigma }_{1a}(p,m_0)+\mathrm{\Sigma }_{1b}(p,m_0),$$

where

$`\mathrm{\Sigma }_{1a}(p,m_0)`$ $`=`$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dk}{(2\pi )^d}}{\displaystyle \frac{\mathrm{\Delta }(k,m_0,\gamma )}{(\overline{p+k})^2+m_0^2}},`$
$`\mathrm{\Sigma }_{1b}(m_0)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }(0,m_0,\gamma ){\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dk}{(2\pi )^d}}\mathrm{\Delta }(k,m_0,\gamma ){\displaystyle \frac{}{m_0^2}}\mathrm{\Delta }^1(k,m_0,\gamma ).`$
Again we might perform an asymptotic expansion of this expression, on the lines traced in Ref. CR , recovering in the scaling limit the (unrenormalized) continuum contribution, and in principle we might evaluate the first nonleading contribution.

However for our purposes it is much more convenient to work directly with quantities chosen in such a way that all renormalization effects are automatically removed, i.e., quantities whose scaling limit is a finite, renormalization-group invariant, amplitude.

The simplest such object is the so-called “renormalized coupling” $`g_r`$, whose continuum (scaling) value $`g_r^{}`$ has been computed to $`O(1/N)`$ in Ref. CPRV-ON .

Actually the formal expression derived in Ref. CPRV-ON  is correct also for the lattice versions of the model, when continuum propagators are replaced by their lattice counterparts. Parametrizing the result in terms of the renormalized mass $`m`$ we therefore obtain

$$g_r(m,\gamma )=m^{d4}\mathrm{\Delta }(0,m,\gamma )\left[1+\frac{1}{N}g_r^{(1)}(m,\gamma )+O\left(\frac{1}{N^2}\right)\right],$$
(13)
where $`g_r^{(1)}`$ is given by

$`g_r^{(1)}(m,\gamma )`$ $`=`$ $`\mathrm{\Delta }(0.m,\gamma ){\displaystyle \frac{\mathrm{\Delta }^1(0,m,\gamma )}{m^2}}(\mathrm{\Sigma }_{1a}(0,m)+\mathrm{\Sigma }_{1b}(m)m^2{\displaystyle \frac{\mathrm{\Sigma }_{1a}}{p^2}}|_0)`$ (14)
$`+`$ $`\left(2{\displaystyle \frac{\mathrm{\Sigma }_{1a}(0,m)}{m^2}}+{\displaystyle \frac{\mathrm{\Sigma }_{1b}(m)}{m^2}}2{\displaystyle \frac{\mathrm{\Sigma }_{1a}}{p^2}}|_02\mathrm{\Delta }^1(0,m,\gamma )T(m)\right),`$
and we defined

$$T(m_0)=_\pi ^\pi \frac{d^dk}{(2\pi )^d}\left(\frac{\mathrm{\Delta }(k,m_0,\gamma )}{\overline{k}^2+m_0^2}\right)^2.$$

It is now a matter of trivial algebraic manipulations to show that the expression for $`g_r^{(1)}`$ can be cast into the form

$`g_r^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{D_0}}{\displaystyle \underset{\pi }{\overset{\pi }{}}}{\displaystyle \frac{d^dk}{(2\pi )^d}}\{{\displaystyle \frac{1}{D}}({\displaystyle \frac{1}{2}}{\displaystyle \frac{}{m^2}}[{\displaystyle \frac{D}{m^2}}+{\displaystyle \frac{2D_0}{\overline{k}^2+m^2}}]{\displaystyle \frac{D_0}{(\overline{k}^2+m^2)^2}}[2D_0+m^2{\displaystyle \frac{D_0}{m^2}}]C)`$ (15)
$`{\displaystyle \frac{1}{2D^2}}[{\displaystyle \frac{D}{m^2}}+{\displaystyle \frac{2D_0}{\overline{k}^2+m^2}}]^2\}.`$
In the asymptotic expansion of $`g_r^{(1)}`$ we may again identify a continuum and a lattice contribution.

Renormalization-group theory insures us about the expected properties of these contributions. In particular the leading continuum term should be finite and should therefore require no UV counterterms. Its evaluation will however allow us to identify the IR counterterms needed in order to regularize the lattice term.

The lattice contribution in turn should not affect the leading (scaling) order, since $`g_r^{}`$ is an invariant amplitude which should not depend on the detailed form of the Hamiltonian. In turn, trivial power counting in $`a`$ shows us that the first correction to scaling is generated by the leading lattice contribution, which is $`O\left((ma)^{4d}\right)`$.

We can prove the following identities, corresponding to similar results of Ref. CPRV-ON :

$`{\displaystyle \frac{D_c}{m^2}}`$ $`=`$ $`{\displaystyle \frac{2}{k^2+4m^2}}\left[(d3)D_cD_{0c}+(4d)\gamma _c\right],`$
$`{\displaystyle \frac{D_{0c}}{m^2}}`$ $`=`$ $`\left({\displaystyle \frac{d}{2}}2\right){\displaystyle \frac{D_{0c}\gamma _c}{m^2}},`$
$`C_c`$ $`=`$ $`{\displaystyle \frac{1}{(k^2+m^2)^2}}\left({\displaystyle \frac{4}{d}}{\displaystyle \frac{k^2}{k^2+m^2}}1\right).`$
As a consequence we obtain in leading order the following continuum contribution, depending on the single (dimensionless) variable $`x\gamma _c(ma)^{4d}`$:

$`g_{r,\mathrm{con}}^{(1)}(x)`$ $``$ $`{\displaystyle \frac{1}{D_{0c}}}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{d^dk}{(2\pi )^d}}`$
$`\times (`$ $``$ $`2\left[{\displaystyle \frac{d3}{k^2+4m^2}}+{\displaystyle \frac{\gamma _c}{D_c}}{\displaystyle \frac{4d}{k^2+4m^2}}+{\displaystyle \frac{D_{0c}}{D_c}}{\displaystyle \frac{3m^2}{(k^2+m^2)(k^2+4m^2)}}\right]^2+2{\displaystyle \frac{(d3)(d5)}{(k^2+4m^2)^2}}`$ (16)
$`+`$ $`{\displaystyle \frac{\gamma _c(4d)}{D_c}}\left[{\displaystyle \frac{2(d5)}{(k^2+4m^2)^2}}+{\displaystyle \frac{3}{2}}{\displaystyle \frac{1}{(k^2+4m^2)(k^2+m^2)}}+{\displaystyle \frac{1}{2}}\left(1{\displaystyle \frac{4}{d}}\right){\displaystyle \frac{1}{(k^2+m^2)^2}}\right]`$
$`+`$ $`{\displaystyle \frac{2\gamma _c}{D_c}}{\displaystyle \frac{(4d)m^2}{d(k^2+m^2)^3}}`$
$`+`$ $`{\displaystyle \frac{D_{0c}}{D_c}}{\displaystyle \frac{m^2}{k^2+m^2}}[{\displaystyle \frac{6(d5)}{(k^2+4m^2)^2}}+{\displaystyle \frac{3}{2}}{\displaystyle \frac{d8}{(k^2+4m^2)(k^2+m^2)}}+{\displaystyle \frac{2}{(k^2+m^2)^2}}]),`$
which can be shown to correspond exactly to the result presented in Ref. CPRV-ON .

In particular the fixed-point value is obtained by setting $`x=0`$, corresponding to $`\gamma _c=0`$, that is the condition for the removal of the leading nonscaling behavior in the large-$`N`$ limit.

From the above result we may immediately read off the structure of counterterms, just by taking the power series expansion in powers of $`m^2`$, and in particular by replacing $`D_c`$ with $`A_{0c}`$ and $`D_{0c}`$ with $`d_0(ma)^{d4}`$. The resulting (singular) expression is

$$\frac{(ma)^{4d}}{d_0}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{d^dk}{(2\pi )^d}\left(2\left[\frac{d3}{k^2}+\frac{\gamma _c}{A_{0c}}\frac{4d}{k^2}\right]^2+2\frac{(d3)(d5)}{(k^2)^2}+2\frac{\gamma _c}{A_{0c}}\left(d4\frac{1}{d}\right)\frac{4d}{(k^2)^2}\right).$$

Let us now compute the lattice contribution. By applying the already defined asymptotic expansions we obtain:

$$\frac{D}{m^2}+\frac{2D_0}{\overline{k}^2+m^2}A_1+2\frac{\gamma _c}{\overline{k}^2}b_0\left(\frac{d}{2}C_0+\frac{1}{(\overline{k}^2)^2}\right)m^{d2}+2\left(A_2\frac{\gamma _c}{(\overline{k}^2)^2}\frac{\delta _0}{\overline{k}^2}\right)m^2+O(m^d)$$

and

$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{m^2}}\left[{\displaystyle \frac{D}{m^2}}+{\displaystyle \frac{2D_0}{\overline{k}^2+m^2}}\right]{\displaystyle \frac{D_0}{(\overline{k}^2+m^2)^2}}\left[2D_0+m^2{\displaystyle \frac{D_0}{m^2}}\right]C`$
$``$ $`A_2{\displaystyle \frac{\delta _0}{\overline{k}^2}}2\gamma _c\left(C_0+{\displaystyle \frac{1}{(\overline{k}^2)^2}}\right)+O(m^{d2}).`$
Notice that the most singular contributions, $`O(m^{d4})`$, have been removed in the above combinations.

Let us now find out the leading IR singularities of the terms appearing in the integral representing $`g_{r,\mathrm{lat}}^{(1)}`$. We must only notice that, for all $`A_n(k)`$, the singular behavior when $`k0`$ is determined by the corresponding behavior of $`A_{nc}(k)`$, with corrections whose singularity is depressed by a factor $`k^2`$. As a consequence we obtain in the IR limit

$`{\displaystyle \frac{1}{A_0}}\left(A_1+2{\displaystyle \frac{\gamma _c}{k^2}}\right)`$ $``$ $`{\displaystyle \frac{2}{k^2}}\left[d3+(4d){\displaystyle \frac{\gamma _c}{A_{0c}}}\right],`$ (17)
$`{\displaystyle \frac{1}{A_0}}\left(A_2{\displaystyle \frac{\delta _0}{\overline{k}^2}}2\gamma _c(C_0+{\displaystyle \frac{1}{(\overline{k}^2)^2}})\right)`$ $``$ $`{\displaystyle \frac{2}{(k^2)^2}}\left[(d3)(d5)+\left(d4{\displaystyle \frac{1}{d}}\right)(4d){\displaystyle \frac{\gamma _c}{A_{0c}}}\right].`$
These singularities are perfectly matched by the terms coming from the $`T^{(\mathrm{IR})}`$ expansion of the continuum contribution, as expected. As a consequence we are able to write down an exact, finite representation of the leading lattice contribution to $`g_r^{(1)}`$, taking the form

$$g_{r,\mathrm{lat}}^{(1)}(m,\gamma _c)\frac{(ma)^{4d}}{d_0}\delta g^{(1)}(\gamma _c),$$
(18)
where

$`\delta g^{(1)}(\gamma _c)`$ $``$ $`{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dk}{(2\pi )^d}}\left({\displaystyle \frac{1}{A_0}}\left(A_2{\displaystyle \frac{\delta _0}{\overline{k}^2}}2\gamma _c\left(C_0+{\displaystyle \frac{1}{(\overline{k}^2)^2}}\right)\right){\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{A_0}}\left(A_1+2{\displaystyle \frac{\gamma _c}{\overline{k}^2}}\right)\right]^2\right)`$ (19)
$``$ $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{d^dk}{(2\pi )^d}}{\displaystyle \frac{2}{(k^2)^2}}([(d3)(d5)+(d4{\displaystyle \frac{1}{d}})(4d){\displaystyle \frac{\gamma _c}{A_{0c}}}]`$
$`[d3+(4d){\displaystyle \frac{\gamma _c}{A_{0c}}}]^2).`$
It is worth observing that $`g_{r,\mathrm{con}}^{(1)}(x)`$ showed a nonanalyticity whose leading dependence was proportional to $`x\mathrm{ln}x`$. Correspondingly $`\delta g^{(1)}(\gamma _c)`$ has a nonanalytic $`\gamma _c\mathrm{ln}\gamma _c`$ dependence. The coefficients of these singularities match properly in order to reproduce an overall dependence proportional to $`\gamma _c(ma)^{4d}\mathrm{ln}ma`$, as expected from general renormalization group arguments because of the anomalous dimension of the leading irrelevant operator.

## VI Integral representation for the gap equation

For the purpose of actual numerical calculations one must find an efficient way of performing lattice momentum integrals. In practice this may be obtained by resorting to parametric (Feynman and Schwinger) representations of the lattice propagators. These representations are especially useful in the case corresponding to the standard nearest-neighbor Hamiltonian.

Let us first consider the integral appearing in the gap equation

$$\beta +2\gamma m_0^2=_\pi ^\pi \frac{d^dp}{(2\pi )^d}\frac{1}{\widehat{p}^2+m_0^2}\mathrm{\Xi }(m_0^2).$$

Introducing Schwinger’s proper time representation we obtain

$$\mathrm{\Xi }(m_0^2)=_0^{\mathrm{}}𝑑\alpha e^{\alpha m_0^2}_\pi ^\pi \frac{d^dp}{(2\pi )^d}e^{2\alpha _\mu (1\mathrm{cos}p_\mu )}=_0^{\mathrm{}}𝑑\alpha e^{\alpha m_0^2}\left[e^{2\alpha }I_0(2\alpha )\right]^d,$$
(20)
where $`I_0`$ is the standard modified Bessel function, admitting for large values of its argument the following asymptotic expansion:

$$e^{2\alpha }I_0(2\alpha )\frac{1}{(4\pi \alpha )^{\frac{1}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(1)^n\mathrm{\Gamma }(n+\frac{1}{2})}{n!\mathrm{\Gamma }(\frac{1}{2}n)}\frac{1}{(4\alpha )^n}.$$

By a proper change of variables the numerical evaluation of the above integral is now possible even in the small $`m_0^2`$ regime.

When $`d=3`$ it is also possible to obtain an analytical expression for $`\mathrm{\Xi }(0)`$ MR , which was first derived in Ref. Watson :

$$\beta _c=\mathrm{\Xi }(0)=\frac{\varkappa ^2+1}{\pi ^2}K^2(\varkappa ),$$
(21)
where $`K`$ is the complete elliptic integral of the first kind and $`\varkappa =(2\sqrt{3})(\sqrt{3}\sqrt{2})`$, and the resulting numerical value is $`\beta _c=0.252731009858663\mathrm{}`$

We can compute $`\mathrm{\Xi }(m_0^2)`$ as an asymptotic expansion around $`m_0^2=0`$:

$$\mathrm{\Xi }(m_0^2)=\underset{n=0}{\overset{\mathrm{}}{}}a_nm_0^{2n}+b_nm_0^{2n+d2}.$$
(22)
It is conceivable that when $`d=3`$ the coefficients $`a_n`$ of the expansion can be computed analytically in terms of elliptic integrals, but we contented ourselves with a numerical calculation of the “lattice terms” appearing in the analytic part of the expansion, while it is always possible to obtain closed form expressions for the coefficients $`b_n`$ of the nonanalytic part.

We can obtain explicit expressions for the coefficients in the expansion (22) by subtracting a proper number of terms of the asymptotic expansion of the Bessel function raised to the power $`d`$. Let us label the coefficients of this expansion according to the equation

$$\left[e^{2\alpha }I_0(2\alpha )\right]^d\frac{1}{(4\pi \alpha )^{\frac{d}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}c_n(d)\alpha ^n;$$
(23)
the coefficients can be computed recursively from the equations

$`c_0(d)`$ $`=`$ $`1,`$
$`c_n(d)`$ $`=`$ $`{\displaystyle \frac{1}{n}}{\displaystyle \underset{k=1}{\overset{n}{}}}(kdn+k){\displaystyle \frac{(1)^k\mathrm{\Gamma }(\frac{1}{2}+k)}{4^kk!\mathrm{\Gamma }(\frac{1}{2}k)}}c_{nk}.`$ (24)
Let us add to the integrand in Eq. (20) the formally vanishing term

$$\frac{1}{(4\pi \alpha )^{\frac{d}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}\frac{(\alpha m_0^2)^n}{n!}\underset{m=0}{\overset{n1}{}}c_m(d)\alpha ^m+\frac{1}{(4\pi \alpha )^{\frac{d}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}c_n(d)\alpha ^n\underset{m=n+1}{\overset{\mathrm{}}{}}\frac{(\alpha m_0^2)^m}{m!},$$

where we have only interchanged the order of the summations in the two contributions. It is now possible to group the first contribution with the original integrand and recognize that the resulting combination defines an analytic function of $`m_0^2`$, since each coefficient

$$a_n_0^{\mathrm{}}𝑑\alpha \frac{(\alpha )^n}{n!}\left\{\left[e^{2\alpha }I_0(2\alpha )\right]^d\frac{1}{(4\pi \alpha )^{\frac{d}{2}}}\underset{m=0}{\overset{n1}{}}c_m(d)\alpha ^m\right\}$$
(25)
in the corresponding power series expansion is a finite, UV and IR regulated integral for any $`d`$ in the range $`2<d<4`$. In turn, the integration of the second contribution may be represented, after trivial resummations and rescalings, in the form

$$\frac{1}{(4\pi )^{\frac{d}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}c_n(d)k_n(d)(m_0^2)^{n+\frac{d}{2}1},$$

where

$$k_n(d)_0^{\mathrm{}}𝑑xx^{n\frac{d}{2}}\left[e^x\underset{m=0}{\overset{n}{}}\frac{(x)^m}{m!}\right]$$

are finite UV and IR regulated integrals in the range $`2<d<4`$; integrating by parts we obtain a recursive equation for $`k_n`$, which can be solved immediately, obtaining $`k_n=\mathrm{\Gamma }(1\frac{d}{2}n)`$, and consequently

$$b_n=\frac{1}{(4\pi )^{\frac{d}{2}}}c_n(d)\mathrm{\Gamma }\left(1\frac{d}{2}n\right).$$
(26)
It is now trivial to set $`d=3`$ in Eq. (24) to obtain an explicit recursive expression for $`b_n`$.

To compute numerical values of $`a_n`$ in $`d=3`$, it is not practical to use Eq. (25), given the slow convergence of the integration for large $`\alpha `$. We found it more convenient to split the $`\alpha `$ integration in Eq. (20) at $`\alpha =1`$. For $`\alpha 1`$ we can expand in powers of $`m_0^2`$ under the integration sign and integrate term by term. For $`\alpha 1`$, we subtract $`n+l`$ terms of the expansion (23), expand the integrand in powers of $`m_0^2`$ up to $`O(m_0^{2n})`$, and integrate term by term; while the resulting integrals converge for any $`l0`$, it is useful to set $`l3`$ to ensure fast convergence of the integration for large $`\alpha `$. The integral of the regulator can be computed analytically in terms of the incomplete Gamma function, and the singular part of its asymptotic expansion reproduces the singular part of the expansion (22).

Collecting all numerical and analytical results, we obtain for $`d=3`$ the following asymptotic expansion:

$`\mathrm{\Xi }(m_0^2)`$ $`=`$ $`\beta _c{\displaystyle \frac{1}{4\pi }}m_00.012164158583022m_0^2+{\displaystyle \frac{1}{32\pi }}m_0^3+0.00083776240606293m_0^4`$ (27)
$``$ $`{\displaystyle \frac{11}{2560\pi }}m_0^50.000066743211781194m_0^6+{\displaystyle \frac{281}{430080\pi }}m_0^7+5.7884488124445\times 10^6m_0^8`$
$``$ $`{\displaystyle \frac{71}{655360\pi }}m_0^95.320777180475\times 10^7m_0^{10}+{\displaystyle \frac{7783}{403701760\pi }}m_0^{11}`$
$`+`$ $`5.107544566044\times 10^8m_0^{12}{\displaystyle \frac{70289}{19377684480\pi }}m_0^{13}+O(m_0^{14}),`$
which gives an error smaller than $`10^{14}`$ in the range $`0m_0^20.1`$.

It is immediate to extract from this expansion the 3-$`d`$ large-$`N`$ nearest-neighbor model values

$$\gamma _00.0060820792915113,\delta _00.00083776240606293.$$

The asymptotic expansion of $`\mathrm{\Xi }^{(n)}(m_0^2)`$, the $`n`$th derivative of $`\mathrm{\Xi }(m_0^2)`$ with respect to $`m_0^2`$, is easily obtained from the above expression. There is an obvious precision loss, but in the above-mentioned range the error in the second derivative is still smaller than $`10^{10}`$.

In the following, we will need to compute $`\mathrm{\Xi }^{(n)}`$ in a fast and accurate way, for generic values of $`m_0^2`$; to this purpose, we tabulated $`\mathrm{\Xi }`$ for values of $`m_0^2`$ on a uniform grid with step $`h=10^3`$ and compute $`\mathrm{\Xi }^{(n)}(m_0^2)`$ by $`n+4`$-point Lagrange interpolation.

## VII Integral representations for the effective propagator

Let us now recall from Ref. CR  the following basic result:

$$D(k,m_0,0)=\frac{1}{2}_0^1𝑑x_\pi ^\pi \frac{d^dq}{(2\pi )^d}\frac{1}{[m_0^2+2_\mu (1z_\mu \mathrm{cos}q_\mu )]^2},$$

where $`z_\mu =\sqrt{1x(1x)\widehat{k}_\mu ^2}`$.

Using Schwinger representation we then obtain:

$`D(k,m_0,0)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}𝑑x{\displaystyle _0^{\mathrm{}}}𝑑\alpha \alpha e^{\alpha m_0^2}{\displaystyle _\pi ^\pi }{\displaystyle \frac{d^dq}{(2\pi )^d}}e^{2\alpha _\mu (1z_\mu \mathrm{cos}q_\mu )}`$ (28)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}𝑑x{\displaystyle _0^{\mathrm{}}}𝑑\alpha \alpha e^{\alpha m_0^2}{\displaystyle \underset{\mu }{}}\left[e^{2\alpha }I_0\left(2\alpha z_\mu \right)\right].`$
Alternatively we might use the Schwinger representation directly and obtain

$$D(k,m_0,0)=\frac{1}{2}_0^{\mathrm{}}𝑑s_0^{\mathrm{}}𝑑te^{(s+t)m_0^2}\underset{\mu }{}\left[e^{2(s+t)}I_0\left(2\sqrt{s^2+t^2+2st\mathrm{cos}k_\mu }\right)\right],$$

that can be reduced to the previous one by the variable change $`s=x\alpha `$, $`t=(1x)\alpha `$.

The direct numerical evaluation of Eq. (28), and especially of its derivatives with respect to $`m_0^2`$, in $`d=3`$ is difficult, particularly for small values of $`k`$ or $`m_0`$; the convergence can be improved dramatically by adding and subtracting a symmetric combination of Bessel functions:

$`D(k,m_0,0)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}𝑑x{\displaystyle _0^{\mathrm{}}}𝑑\alpha \alpha e^{\alpha (6+m_0^2)}\left[{\displaystyle \underset{\mu }{}}I_0\left(2\alpha z_\mu \right)I_0^3\left(2\alpha \overline{z}\right)\right]`$ (29)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}𝑑x{\displaystyle \frac{1}{\overline{z}^2}}\mathrm{\Xi }^{}\left({\displaystyle \frac{m_0^2+6}{\overline{z}}}1\right),`$
where $`\overline{z}=\frac{1}{3}(z_1+z_2+z_3)`$; the subtracted integral and its first few derivatives with respect to $`m_0^2`$ can now be computed accurately by Gauss-Legendre integration on a small grid; the integration in $`x`$ of $`\mathrm{\Xi }^{}`$ and of its derivatives is also easy, once a few singular terms of Eq. (22) have been subtracted.

In order to identify the lattice contributions to the asymptotic expansion of $`D`$ we must expand the integrand in powers of $`m_0^2`$. The IR singularities turn into unsuppressed positive powers of $`\alpha `$, present in the large-$`\alpha `$ regime when $`z_\mu 1`$. These singularities become worse and worse with increasing powers of $`m_0^2`$. In order to classify them according to their degree we need to consider the homogeneous expansion of the integrand in powers of $`\alpha ^1`$ and of $`\sqrt{1z_\mu ^2}x(1x)\widehat{k}_\mu ^2`$.

Recalling once more the asymptotic expansion of the Bessel function, we can write down the homogeneous expansion in the form:

$$\underset{\mu }{}e^{2\alpha }I_0(2\alpha z_\mu )\frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}C_n[d,\alpha x(1x)\widehat{k}_\mu ^2]\alpha ^n,$$
(30)
where in turn

$$C_n[d,\alpha x(1x)\widehat{k}_\mu ^2]\underset{m=0}{\overset{2n}{}}\gamma _{nm}(d)[\alpha x(1x)\widehat{k}^2]^m$$
(31)
and $`\gamma _{nm}(d)`$ may show a dependence on the ratios $`\widehat{k}^{2p}/(\widehat{k}^2)^p`$.

Repeating the procedure developed in the previous section we may now obtain the following decomposition:

$$D(k,m_0,0)=\underset{n=0}{\overset{\mathrm{}}{}}R_n(k)m_0^{2n}+D_s(k,m_0),$$
(32)
where

$$R_n(k)\frac{(1)^n}{2n!}_0^1𝑑x_0^{\mathrm{}}𝑑\alpha \alpha ^{n+1}\left[\underset{\mu }{}e^{2\alpha }I_0(2\alpha z_\mu )\frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}\underset{m=0}{\overset{n}{}}C_m\alpha ^m\right]$$
(33)
is the part of $`A_n(k)R_n(k)+S_n(k)`$ which shows a regular dependence on $`k`$ in the $`k0`$ limit, and

$$D_s(k,m_0)\frac{1}{2}_0^1𝑑x_0^{\mathrm{}}𝑑\alpha \frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}\underset{n=0}{\overset{\mathrm{}}{}}\left[e^{\alpha m_0^2}\underset{m=0}{\overset{n1}{}}\frac{(\alpha m_0^2)^m}{m!}\right]C_n\alpha ^{1n}$$
(34)
is a calculable expression including all the singular dependence on $`k`$ and admitting an asymptotic expansion in the form

$$D_s(k,m_0)\underset{n=0}{\overset{\mathrm{}}{}}\left[S_n(k)m_0^{2n}+B_n(k)m_0^{2n+d2}\right].$$

In order to compute $`D_s(k,m_0)`$ let us notice that it may be also expressed in the form

$$D_s(k,m_0)=\underset{n=0}{\overset{\mathrm{}}{}}\underset{m=0}{\overset{2n}{}}\gamma _{nm}(d)I_{nm}(m_0^2)$$
(35)
where

$$I_{nm}(m_0^2)\frac{1}{2}_0^1𝑑x[x(1x)\widehat{k}^2]^m_0^{\mathrm{}}𝑑\alpha \frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}\left[e^{\alpha m_0^2}\underset{m=0}{\overset{n1}{}}\frac{(\alpha m_0^2)^m}{m!}\right]\alpha ^{1+mn}$$
(36)
enjoys the property

$$\frac{dI_{nm}}{dm_0^2}=I_{n1,m}.$$

Therefore we only need to evaluate

$`I_{0m}(m_0^2)`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}𝑑x[x(1x)\widehat{k}^2]^m{\displaystyle _0^{\mathrm{}}}𝑑\alpha e^{\alpha x(1x)\widehat{k}^2\alpha m_0^2}{\displaystyle \frac{\alpha ^{1\frac{d}{2}+m}}{(4\pi )^{\frac{d}{2}}}}`$ (37)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Gamma }(2\frac{d}{2}+m)}{(4\pi )^{\frac{d}{2}}}}{\displaystyle _0^1}𝑑x{\displaystyle \frac{[x(1x)\widehat{k}^2]^m}{[m_0^2+x(1x)\widehat{k}^2]^{2\frac{d}{2}+m}}}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Gamma }(2\frac{d}{2})}{(4\pi )^{\frac{d}{2}}}}(\widehat{k}^2)^m\left({\displaystyle \frac{}{\widehat{k}^2}}\right)^m{\displaystyle _0^1}{\displaystyle \frac{dx}{[m_0^2+x(1x)\widehat{k}^2]^{2\frac{d}{2}}}}.`$
A trivial comparison with previous results shows that the last integral is directly related to the continuum propagator by the replacement $`k^2\widehat{k}^2`$, and therefore, setting $`\widehat{\xi }^2=1+4m_0^2/\widehat{k}^2`$ we obtain

$$I_{0m}(m_0^2)=(\widehat{k}^2)^m\left(\frac{}{\widehat{k}^2}\right)^m\left[a_0\widehat{\xi }^{d3}(\widehat{k}^2)^{\frac{d}{2}2}+\frac{b_0}{\widehat{k}^2\widehat{\xi }^2}{}_{2}{}^{}F_{1}^{}(\frac{d1}{2},1,\frac{d}{2},1\frac{1}{\widehat{\xi }^2})m_0^{d2}\right],$$
(38)
where we must appreciate that the above expression is naturally decomposed into an analytic and a nonanalytic term, which implies that we can immediately relate all coefficients $`S_n(k)`$ to the derivatives with respect to $`\widehat{k}^2`$ of $`\widehat{\xi }^{d3}(\widehat{k}^2)^{\frac{d}{2}2}`$.

Straightforward manipulations lead to the general relationship

$$S_n(k)=\frac{a_0}{n!}(\widehat{k}^2)^{\frac{d}{2}2}\left(\frac{4}{\widehat{k}^2}\right)^n\underset{p=0}{\overset{n}{}}\frac{(1)^p\mathrm{\Gamma }(\frac{d1}{2})}{\mathrm{\Gamma }(\frac{d1}{2}n+p)}\left(\frac{\widehat{k}^2}{4}\right)^p\underset{q=0}{\overset{2p}{}}\frac{(1)^q\mathrm{\Gamma }(\frac{d}{2}1n+p)}{\mathrm{\Gamma }(\frac{d}{2}1n+pq)}\gamma _{pq}(d).$$
(39)
The above expression brings into evidence a peculiar feature exhibited by the functions $`S_n(k)`$ when $`d=3`$. In this case the arguments of the function $`\mathrm{\Gamma }(1n+p)`$ appearing in the denominator are integer nonpositive numbers whenever $`pn`$, and therefore the corresponding contributions to the sum vanish, giving

$$S_n(k)=(1)^n\frac{a_0}{n!}(\widehat{k}^2)^{\frac{1}{2}}\underset{q=0}{\overset{2n}{}}\frac{(1)^q\mathrm{\Gamma }(\frac{1}{2})}{\mathrm{\Gamma }(\frac{1}{2}q)}\gamma _{nq}(3).$$
(40)
When $`d=3`$ the singular contributions to $`A_n(k)`$ are proportional to $`(\widehat{k}^2)^{\frac{1}{2}}`$ for all $`n`$, and this is basically a consequence of the vanishing dependence on $`\xi `$ in the analytic part of the continuum propagator.

The simplest example of this procedure is offered by the integral representation of $`A_0(k,0)`$:

$$A_0(k,0)=a_0(\widehat{k}^2)^{\frac{d}{2}2}+R_0(k),$$
(41)
where

$$R_0(k)\frac{1}{2}_0^1𝑑x_0^{\mathrm{}}𝑑\alpha \alpha \left[\underset{\mu }{}e^{2\alpha }I_0(2\alpha z_\mu )\frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}\right].$$
(42)
This representation is worth a few observations. Notice that the explicit term might be simply obtained from the expression of $`A_{0c}`$ by the replacement $`k^2\widehat{k}^2`$, and therefore it reproduces exactly the singular behavior of $`A_0`$ when $`k0`$. Moreover $`R_0(k)`$ has, by construction, a finite limit when $`k0`$. It is easy to check that $`R_0(0)=\gamma _0,`$ thus verifying explicitly that the choice $`\gamma =\gamma _0`$ leads to the $`O(a^{4d})`$ improvement of the lattice massless propagator.

By the same technique we obtain

$$A_1(k)=a_0\left[\frac{2(d3)}{\widehat{k}^2}+\frac{(6d)(4d)}{16}\frac{\widehat{k}^4}{(\widehat{k}^2)^2}\frac{(8d)}{16}\right](\widehat{k}^2)^{\frac{d}{2}2}+R_1(k),$$
(43)
where

$`R_1(k)`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}dx{\displaystyle _0^{\mathrm{}}}d\alpha \alpha ^2[{\displaystyle \underset{\mu }{}}e^{2\alpha }I_0(2\alpha z_\mu )`$ (44)
$``$ $`{\displaystyle \frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}}(1+{\displaystyle \frac{d}{16\alpha }}+{\displaystyle \frac{1}{4}}x(1x)\widehat{k}^2{\displaystyle \frac{1}{4}}\alpha x^2(1x)^2\widehat{k}^4)].`$
Again we notice that the explicit term reproduces the leading singularities, and moreover $`R_1(0)=\delta _0`$.

We also mention the result for $`A_2(k)`$, since it will be needed in the evaluation of $`\delta g^{(1)}`$:

$$A_2(k)=\frac{1}{2}a_0\left[\frac{a_{2,0}}{(\widehat{k}^2)^2}+\frac{a_{2,1}}{\widehat{k}^2}+a_{2,2}\right](\widehat{k}^2)^{\frac{d}{2}2}+R_2(k),$$
(45)
where

$`a_{2,0}`$ $`=`$ $`4(d3)(d5),`$
$`a_{2,1}`$ $`=`$ $`{\displaystyle \frac{(d3)(d12)}{8}}+{\displaystyle \frac{(d3)(d6)(d8)}{8}}{\displaystyle \frac{\widehat{k}^4}{(\widehat{k}^2)^2}},`$
$`a_{2,2}`$ $`=`$ $`{\displaystyle \frac{(d16)(d8)}{512}}+{\displaystyle \frac{(d4)(d6)(d8)}{256}}{\displaystyle \frac{\widehat{k}^4}{(\widehat{k}^2)^2}}`$
$`+`$ $`{\displaystyle \frac{(d4)(d6)(d8)}{64}}{\displaystyle \frac{\widehat{k}^6}{(\widehat{k}^2)^3}}+{\displaystyle \frac{(d4)(d6)(d8)(d10)}{512}}{\displaystyle \frac{(\widehat{k}^4)^2}{(\widehat{k}^2)^4}},`$
and we obtain

$$R_2(k)\frac{1}{4}_0^1𝑑x_0^{\mathrm{}}𝑑\alpha \alpha ^3\left[\underset{\mu }{}e^{2\alpha }I_0(2\alpha z_\mu )\frac{e^{\alpha x(1x)\widehat{k}^2}}{(4\pi \alpha )^{\frac{d}{2}}}\left(1+r_{2,1}(k)+r_{2,2}(k)\right)\right],$$
(46)
where in turn

$`r_{2,1}(k)`$ $`=`$ $`{\displaystyle \frac{d}{16\alpha }}+{\displaystyle \frac{1}{4}}x(1x)\widehat{k}^2{\displaystyle \frac{1}{4}}\alpha x^2(1x)^2\widehat{k}^4,`$
$`r_{2,2}(k)`$ $`=`$ $`{\displaystyle \frac{d(d+8)}{512\alpha ^2}}+{\displaystyle \frac{1}{64\alpha }}(d+2)x(1x)\widehat{k}^2+{\displaystyle \frac{1}{64}}\left(2(\widehat{k}^2)^2+(8d)\widehat{k}^4\right)x^2(1x)^2`$
$``$ $`{\displaystyle \frac{1}{16}}\alpha x^3(1x)^3(2\widehat{k}^6+\widehat{k}^4\widehat{k}^2)+{\displaystyle \frac{1}{32}}\alpha ^2x^4(1x)^4(\widehat{k}^4)^2.`$
An obvious relationship exists between the values $`R_n(0)`$ and the coefficients of the analytic part of the asymptotic expansion of the function $`\mathrm{\Xi }^{}(m_0^2)`$. One easily finds that $`R_n(0)=\frac{1}{2}(n+1)a_{n+1}`$.

It should be by now clear that a trivial generalization of the same technique allows for an explicit, albeit more cumbersome, evaluation of all the functions $`B_n(k)`$. We could verify that the correct expressions for $`B_0(k)`$ and $`B_1(k)`$ are reproduced, and we computed $`B_2(k)`$ for a better accuracy check of our numerical estimates. The result is too cumbersome to be reported here.

In practice, $`A_0`$, $`A_1`$, and $`A_2`$ can be computed numerically in $`d=3`$ exploiting the subtraction (29); since

$$e^{6t}\left[\underset{\mu }{}I_0(2\alpha z_\mu )I_0^3(2\alpha \overline{z})\right]e^{\alpha x(1x)\widehat{k}^2}x^2(1x)^2\frac{3\widehat{k}^4\left(\widehat{k}^2\right)^2}{48(4\pi \alpha )^{3/2}},$$
(47)
the subtraction is enough to regularize $`A_1`$ and $`A_2`$. Let us consider, e.g., $`A_1`$; by straightforward manipulations we obtain

$`A_1(k)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _0^1}dx\{{\displaystyle _0^1}d\alpha \alpha ^2e^{6\alpha }[{\displaystyle \underset{\mu }{}}I_0(2\alpha z_\mu )I_0^3(2\alpha \overline{z})]`$ (48)
$`+{\displaystyle \frac{1}{\overline{z}^3}}\mathrm{\Xi }^{\prime \prime }({\displaystyle \frac{6}{\overline{z}}}1){\displaystyle \frac{1}{16\pi \left(x(1x)\widehat{k}^2\right)^{3/2}}}\};`$
it is useful to overregulate the $`x`$ integration by subtracting the r.h.s. of Eq. (47), integrated over $`\alpha `$. For $`A_2`$ we follow the same procedure; the subtractions are more complicated and not worth writing here. The computation of $`A_0`$ is similar but of course easier.

We verified explicitly that $`A_n`$, computed in this way, and $`B_n`$ are consistent with Eq. (7) for $`n2`$.

## VIII Numerical results

It is worthwhile to present a selection of the numerical results that we obtained in $`d=3`$ for the nearest-neighbor formulation ($`c=\frac{1}{12}`$).

In Figs. 3, 4, and 5 we compare $`g_r^{(1)}(m,\gamma _c)`$, obtained by direct evaluation of Eq. (15), with $`g_{r,\mathrm{con}}^{(1)}(\gamma _c(ma)^{4d})+g_{r,\mathrm{lat}}^{(1)}(m,\gamma _c)`$, obtained from Eqs. (16), (18), and (19), for three values of $`\gamma _c`$; of special interest is the value $`\gamma _c=\gamma _0`$, i.e., $`\gamma =0`$, corresponding to the nonlinear $`\sigma `$ model.

$`\delta g^{(1)}(\gamma _c)`$ is plotted in Fig. 6; of special interest is the value

$$\delta g^{(1)}(0)=0.0049699\mathrm{}$$

## IX Construction of the improved Hamiltonian

As we mentioned in the Introduction, the improvement procedure aims at a systematic cancellation of the next-to-leading effects in the invariant amplitudes. The renormalization-group theory insures us that a single choice of $`\gamma `$ exists such that this cancellation occurs in all amplitudes. It is therefore sufficient to find the value $`\gamma =\gamma ^{}`$ for which the cancellation occurs in the renormalized coupling.

In the context of the $`1/N`$ expansion we may assume $`\gamma ^{}`$ to admit an expansion in powers of $`1/N`$:

$$\gamma ^{}=\gamma _0^{}+\frac{1}{N}\gamma _1^{}+O\left(\frac{1}{N^2}\right).$$

We have already recognized in Sec. III that $`\gamma _0^{}=\gamma _0`$. We may therefore define

$$\gamma _c^{}\gamma ^{}\gamma _0\frac{1}{N}\gamma _1^{}.$$

Substituting this result in the expression of $`g_r`$ and expanding in powers of $`m`$ we then obtain:

$$g_r(ma,\gamma ^{})\frac{1}{d_0}\left[1+\frac{1}{N}g_{r,\mathrm{con}}^{(1)}(0)+\frac{1}{N}\frac{1}{d_0}(\delta g_1(0)\gamma _1^{})(ma)^{4d}+O(m^2a^2)+O\left(\frac{1}{N^2}\right)\right].$$

We then recognize that the condition for the cancellation of the first nonleading contribution is

$$\gamma _1^{}=\delta g_1(0)_\pi ^\pi \frac{d^dk}{(2\pi )^d}\left(\frac{1}{\stackrel{~}{A}_0}\left(A_2\frac{\delta _0}{\overline{k}^2}\right)\frac{1}{2}\left(\frac{A_1}{\stackrel{~}{A}_0}\right)^2\right)+_{\mathrm{}}^+\mathrm{}\frac{d^dk}{(2\pi )^d}\frac{4(d3)}{(k^2)^2},$$

where $`\stackrel{~}{A}_0A_0(k,\gamma _0)`$.

## X Conclusions

In the case of the three-dimensional $`O(N)`$ models with standard nearest-neighbor interactions, our analytical results led us to the prediction

$`\gamma ^{}0.00608207+{\displaystyle \frac{1}{N}}\mathrm{\hspace{0.17em}0.0049699}+O\left({\displaystyle \frac{1}{N^2}}\right).`$ (49)
We might now obtain an estimate for the value $`N_c`$ for which $`\gamma ^{}`$ vanishes. The numerical value of the present estimate is admittedly not very promising, but, as we mentioned in the introduction, the essential feature of $`\gamma _1^{}`$ is its positive sign, suggesting that $`N_c`$ may exist and possibly be within the range of convergence of the $`1/N`$ expansion.

The physical interpretation of $`N_c`$ amounts to the statement that, for this special value of $`N`$, the lattice version of the nonlinear $`\sigma `$ model ($`\gamma =0`$) in three dimensions shows the absence of leading irrelevant operators in the expansion of the Hamiltonian into scaling fields.

For all values $`NN_c`$, one would get $`\gamma ^{}>0`$, and as a consequence one may proceed to analyze the models in their improved version, both by numerical Monte Carlo methods and by perturbative expansion techniques. Numerical evidence shows that this is actually the case for the physically interesting cases $`N3`$. It might be interesting to perform a numerical study with the purpose of estimating $`N_c`$. A naïve extrapolation from the known numerical values of $`\gamma ^{}(N)`$,

$$\gamma ^{}(1)=0.0159(3),\gamma ^{}(2)=0.0078(2),\gamma ^{}(3)=0.0043(4),\gamma ^{}(4)=0.0021(7),$$

obtained from Refs. CPRV-impr-O1 , CHPRV-impr-O2 , CHPRV-impr-O3 , and Hasenbusch-O4  respectively using the formula $`\gamma =\beta _c^2/(8\lambda N)`$, suggests $`N_c5÷6`$.

Concerning the extension of the $`1/N`$ expansion itself to the region $`N<N_c`$, we must cautiously mention that some dramatic change in the analytical behavior of the function $`\gamma ^{}(N)`$ may certainly occur at $`N=N_c`$. It is not possible to perform a strong-coupling expansion of the models when $`\gamma <0`$, as one may immediately realize from an analysis of the gap equation.

We would like to mention that the numerical evaluation of higher orders of the $`1/N`$ expansion is technically not beyond reach, along the lines traced by Ref. Flyvbjerg  and exploiting the more accurate results for the effective propagator obtained in the present paper.

Finally, as a consequence of the discussion of the previous sections, it should be clear that, in the class of models with next-to-nearest neighbor interactions, it is always possible to find a choice of Hamiltonian parameters such that improvement becomes possible for arbitrary values of $`N`$.