\section{Introduction}
In reference \cite{Caron-Huot:2016cwu}, the following formula was proposed (see also \cite{EliasMiro:2020tdv} for a very transparent derivation)
\begin{equation}\label{chw}
e^{-i\pi D}F^*_{\cal O}=S F^*_{\cal O}\,,
\end{equation}
relating the action of the dilation operator $D$ on some renormalized form factor of an operator $\cal O$, called $F_{\cal O}$, to that of the $S$ matrix on the same $F_{\cal O}$.

Despite its indisputable beauty, the formula, as it is, is hard to put into use. The literature that followed this path has mainly focussed on its one-loop implementation \cite{EliasMiro:2020tdv,Baratella:2020lzz}, and on some special two-loop examples which, for a reason or another, are free from infrared (IR) divergences \cite{EliasMiro:2020tdv,Bern:2020ikv,EliasMiro:2021jgu}.
The purpose of this work is to make the first steps into the implementation of \eq{chw} at two loops in cases where IR divergences are present. The key to this will be
\begin{itemize}
\item[(i)] the use of form factors where $\cal O$ is a conserved operator and
\item[(ii)] the use of general results about the structure of IR divergences, especially factorization.
\end{itemize}

Since we are mainly interested in a proof of principle, we will only touch on the simplest examples of a quartic scalar theory, maximally supersymmetric Yang-Mills and Yukawa theory. However, insight into these examples will instruct us on the more general scenario.

Following the notation of \cite{EliasMiro:2020tdv}, where $F_{\cal O}({n})\equiv\, _{\rm out}\langle {n}|{\cal O}(0) |0\rangle$ and $S_{n,m}\equiv\, _{\rm out}\langle {n}| {m}\rangle_{\rm in}$, we can rewrite \eq{chw} slightly more explicitly as
\begin{equation}
e^{-i\pi D}F^*_{\cal O}({n})=\sum_{{m}\geq 2} S_{n,m} F^*_{\cal O}({m})\,.
\end{equation}
With an abuse of notation, we denoted by $n,m$ some asymptotic state made of $n$ and $m$ particles respectively. The sum over $m$ implies an integration over the phase space of the corresponding multi-particle state.
The dilation operator $D\equiv p_i\partial_{p_i}$ can be traded in massless theories -- like the ones we consider -- with $-\mu\partial_{\mu}$, so that its action on form factors can be deduced from the renormalization group equations \cite{Caron-Huot:2016cwu}
\begin{equation}\label{Ddef}
D F_{\cal O}=-\mu\partial_{\mu}F_{\cal O}=(\gamma_{\rm UV}-\gamma_{\rm IR}+\beta_g\partial_g)F_{\cal O}\,.
\end{equation}
Notice that \eq{chw} should be interpreted as a matrix equation in the space of all local operators $\cal O$, accounting for possible RG mixing effects. However, since (i) we consider conserved operators for which $\gamma_{\rm UV}$ is zero and (ii) IR divergences do not mix form factors (color mixing is a very peculiar, ``almost diagonal'' mixing), we can safely interpret \eq{chw} to be an equation for the form factor of some fixed $\cal O$, say the energy-momentum tensor $T^{\mu\nu}$.

The version of \eq{chw} we will work with reads as follows
\begin{equation}\label{master}
e^{i\pi (\gamma_{\rm IR}-\beta_g\partial_g)}F^*_{\cal J}({2})=\sum_{{m}\geq 2} S_{2,m} F^*_{\cal J}({m})\,,
\end{equation}
where $\cal J$ is some conserved current,\footnote{In general, that $\cal J$ satisfies a Ward identity is not enough to conclude that $\gamma_{\rm UV}=0$ \cite{Collins:2005nj,Baratella:2021guc}. We will come back to this in due time.} and the choice $n=2$ corresponds to working with minimal form factors (i.e. those that do not vanish in the free theory limit), which is surely convenient to do. The state `2' will typically be made of a particle-antiparticle pair, so that \eq{master} should be read as an equation for the IR divergences associated to this \emph{state} (at variance with $\gamma_{\rm UV}$, which is a property of the operator $\cal O$). Once $\gamma_{\rm IR}$ is extracted with \eq{master}, it will enter universally in all processes involving these particles. For notational simplicity, we will omit the suffix ${\cal J}$ from now on.

\section{One loop}\label{1loop}

Expanding \eq{master} to one-loop order, with $S=1+i{\cal M}$, we get
\begin{equation}\label{oneloop}
\gamma_{\rm IR}^{(1)} \,F_0(2)=\frac{1}{\pi}\,\sum_{2'}\,{\cal M}_{2,2'}^{(0)} \,F_0(2')\,,
\end{equation}
where $F_0(2)\equiv\lim_{g\to 0}F(2)$ is the free theory form factor, which is real. By $x^{(\ell)}$ we denote the $\ell$-loop contribution to $x$. The sum over $2'$ denotes a sum over all relevant intermediate two-particle states, integrated over their phase space, while $F_0(2')$ denotes the $\cal J$ form factor over the state $2'$. States with $m\geq 3$ provide a higher-loop effect and drop out of the one-loop equation. More explicitly, \eq{oneloop} can be rewritten as the following two-particle phase space integral against a conventionally normalized amplitude ${\cal A}_{2,2'}$ \cite{Baratella:2020lzz}
\begin{equation}
\gamma_{\rm IR}^{(1)} \,F_0(2)=\frac{1}{32\pi^3}\,\sum_{2'}\int_{4\pi}d\Omega_{2'}\,{\cal A}_{2,2'}^{(0)} \,F_0(2')\,.
\end{equation}
This is an equation for $\gamma_{\rm IR}$ at one loop. To extract it, we assume that it has the following structure
\begin{equation}\label{gammaIRdimreg}
\gamma_{\rm IR}={\gamma_{\rm cusp}}\,\ln \frac{-s}{\mu^2}+\gamma_{\rm coll}\,,
\end{equation}
which is certainly true at one loop in any theory, and conjectured to be valid at all orders in QCD \cite{Becher:2009cu,Becher:2014oda} and maximally supersymmetric Yang-Mills \cite{Bern:2005iz} (from which the name `cusp' is borrowed in a perhaps fanciful way).

With this expression for $\gamma_{\rm IR}$, how do we match the left- and right-hand-side of \eq{oneloop}? Surprisingly enough, we learned from \cite{Caron-Huot:2016cwu} that, in order to make sense of its appearance in \eq{oneloop}, the `cusp' contribution to the infrared anomalous dimension {has to} be identified with a certain phase space integral, having all in all
\begin{equation}\label{gammaIR}
\gamma_{\rm IR}=\frac{\gamma_{\rm cusp}}{2}\int_0^{\pi}\frac{d\cos\theta}{\sin^2(\theta/2)}+\gamma_{\rm coll}\,,
\end{equation}
where the polar angle $\theta$ parametrizes, together with the azimuthal angle $\phi$, the phase space of particles $2'$ with respect to some reference frame defined by the 2 system, so that $\theta\to 0$ corresponds to a forward scattering. Notice that the integral in \eq{gammaIR} is divergent.\footnote{\eq{gammaIR} is the appropriate one for distinguishable particles in the cut. For identical particles, we need to substitute $\sin^{-2}(\theta/2)\to 2\sin^{-2}\theta$. This is thoroughly discussed in \cite{Baratella:2020dvw}.} By enforcing the validity of \eq{oneloop}, we can now uniquely fix the finite coefficients $\gamma^{(1)}_{\rm cusp}$ and $\gamma^{(1)}_{\rm coll}$. Of course, the non-vanishing of $\gamma^{(1)}_{\rm cusp}$ is possible only when the phase space integral on the right-hand-side is also divergent. This point of view on \eq{oneloop} was not particularly emphasized in previous related works.

A rigorous way to go from \eq{gammaIR} to \eq{gammaIRdimreg} is via dimensional regularization. By promoting the divergent integral multiplying $\gamma_{\rm cusp}$ in \eq{gammaIR} to its finite version in $4-2\epsilon$ dimensions, we find
\begin{align}
\int_0^\pi \frac{d\cos\theta}{\sin^2(\theta/2)}\to n(\epsilon)\left(\frac{-s}{\mu^2}\right)^{-\epsilon}\int_0^\pi \frac{d\cos\theta \,(\sin\theta)^{-2\epsilon}}{\sin^2(\theta/2)} =-2\left[\frac{1}{\epsilon}-\ln\frac{-s}{\mu^2}+\ldots \right],
\end{align}
where $n(\epsilon)$ is a scheme-dependent quantity that tends to 1 for $\epsilon\to 0$.\footnote{It is important to keep track of the scheme dependence when studying the RG at two loops. However, the detailed study of scheme dependence is not the main focus here.} After subtraction of divergences, we precisely match \eq{gammaIR} to \eq{gammaIRdimreg}. In practice, two approaches are possible which of course give equivalent results. Either one adopts \eq{gammaIR} and computes the right-hand-side of \eq{oneloop} in four dimensions (in general with divergences in the phase space integral for $\theta\to 0$) or, alternatively, one can calculate the phase space integral on the right-hand-side of \eq{oneloop} in dimensional regularization, and match the result to $\gamma_{\rm IR}F_0(2)$ after subtraction, with $\gamma_{\rm IR}$ as in \eq{gammaIRdimreg}.

In a sense, the second approach is more satisfactory, as it treats coherently the divergences appearing in phase space integrals and those appearing in loop integrals (which will become relevant when expanding \eq{master} to higher loop orders). However, the important thing in practice is to consistently compare apples to apples and oranges to oranges.

\vspace{.5cm}

Collinear and cusp divergences have a markedly different structure, that allows one to separate them. The first ones are related to a rational contribution to \eq{oneloop}, while the second ones are purely logarithmic. (The possibility of having logarithmic pieces is special to $\gamma_{\rm IR}$.)

There is indeed an elegant way to project out the `cusp' contribution from \eq{oneloop}, which consists in keeping only its rational part. We can express this as
\begin{equation}\label{gammacoll}
\gamma_{\rm coll}^{(1)} \,F_0(2)=\frac{1}{\pi}\,{\cal R}\,\sum_{2'}\,{\cal M}_{2,2'}^{(0)} \,F_0(2')\,.
\end{equation}
There are several efficient ways to implement the action of $\cal R$, which is defined abstractly as the projector onto the rational part of the function it acts on \cite{Mastrolia:2009dr,Arkani-Hamed:2008owk}. The possibility of separating rational and logarithmic yields can be very useful in practice, since one is often interested in UV divergences of form factors, and not in IR ones. In this case, logarithmic contributions to $\gamma$ can be simply projected out, leaving out in general UV and collinear divergences together. Since the two are structurally identical (for comparison with more standard approaches, both of them appear as simple poles in dimensional regularization), it is important to be able to separate them. This is where the study of conserved operators, with $\gamma_{\rm UV}=0$, becomes relevant, together with the fact that IR divergences are universal and depend only on long distance physics, i.e. on the particles partaking in the process.

In summary, if one is interested in capturing UV divergences in a massless on-shell scheme, it is crucial to have knowledge of the collinear divergences in order not to mistake them for UV ones. Instead, cusp divergences can be projected out since they are structurally peculiar and cannot be mistaken.

\section{Two loops}

The aim of this section is to present the general procedure for extracting IR anomalous dimensions up to two loops. The following simple assumption plays a crucial role in simplifying the problem: that the (one-loop) renormalized form factor $F(2)$ can be factorized, i.e.
\begin{equation}\label{factorizable}
F(2)=C(\lambda)F_0(2)\,,
\end{equation}
where $F_0$ is the free-theory limit, and $\lambda$ is some squared coupling (e.g. $g^2$ in gauge theories). See \cite{Yang:2019vag} for an example in super-Yang-Mills.
Using this, and throwing away contributions from the three-loop order on, we find
\begin{align}
e^{i\pi (\gamma_{\rm IR}-\beta_\lambda \partial_\lambda)}F^*(2)&=\sum_{2'}S_{2,2'} F^*(2')+\sum_{3}S_{2,3} F^*(3)+\ldots \nonumber \\
e^{i\pi (\gamma_{\rm IR}-\beta_\lambda \partial_\lambda)}C^* F_0(2)&=\sum_{2'}(C')^* S_{2,2'} F_0(2')+\sum_{3}S_{2,3} F^*(3)+\ldots \nonumber \\
\left[ (C^*)^{-1} e^{i\pi (\gamma_{\rm IR}-\beta_{\lambda} \partial_{\lambda})} C^* \right] F_0(2) &=\sum_{2'}\left(\frac{C'}{C}\right)^* S_{2,2'} F_0(2')+({C^*})^{-1}\sum_{3}S_{2,3} F^*(3)+ \ldots \nonumber
\end{align}
Since the three particle cut is already by itself a two-loop effect, and considering that
\begin{equation}
C=1+a+{\cal O}(\lambda^2)\,,
\end{equation}
(and similarly for $C'$) we get that, up to three-loop effects, the master formula \eq{master} simplifies to
\begin{equation}\label{master2}
e^{i\pi \left(\gamma_{\rm IR}-\frac{1}{\lambda}\,a^* \beta_\lambda^{(1)}\right)}F_0(2)=\sum_{2'}e^{\left(a'-a\right)^*} S_{2,2'}F_0(2')+\sum_3 S_{2,3}\,F(3)+{\cal O}(\lambda^3)\,,
\end{equation}
where $\beta_\lambda\equiv \beta_\lambda^{(1)}+{\cal O}(\lambda^3)$, and we used that $\beta^{(1)}_\lambda\propto\lambda^2$ and $a\propto \lambda$. Finally, we dropped the * from $F^*(3)$ since it is a real (tree-level) object in the two-loop approximation. \eq{master2} is the central formula of this work.
Before moving to an example of its implementation, let us collect the main observations which are general in nature.
\begin{enumerate}
\item The separation in loop orders is pretty clear in \eq{master2}. The zero-loop order is just the identity $F_0=F_0$ (on the right-hand-side, take the 1 from $S_{2,2}=1+i{\cal M}_{2,2}$), while the one-loop order is easily seen to reduce to \eq{oneloop}, as the terms involving $a^*$ and $S_{2,3}$ are two-loop by themselves.
\item To go to two loops, \eq{master2} tells us that $F(2)$ should be known up to one loop, including finite terms that cannot be computed with \eq{oneloop}. The relevant information is captured by $a$, which is of the form
\begin{equation}\label{aeq}
a\equiv a_{\rm fin}^{(1)}-\frac{1}{2}\,\gamma_{\rm coll}^{(1)}\ln\frac{-s}{\mu^2}-\frac{1}{4}\,\gamma_{\rm cusp}^{(1)}\left(\ln\frac{-s}{\mu^2}\right)^2\,.
\end{equation}
\item Besides the full one-loop computation of $F(2)$, the other new ingredients at two loops are (i) the three-particle phase space integral, here schematically written as ${S}_{2,3}F(3)$, with all quantities at tree-level, and (ii) the two-particle phase space integral of $F_0(2)$ with the one-loop order $2\to 2$ amplitude ${S}^{(1)}_{2,2}$. This is shown in Fig.~\ref{fig:2loop}.
\item While the zero- and one-loop orders of \eq{master} produce respectively a real and a purely imaginary equation, starting from two loops we can split the equations into a real and an imaginary part, leading to further simplifications on which we will come back in Sec.~\ref{SYM}.
\item Depending on the properties specific to the theory under consideration, \eq{master2} can be simplified even further, for example in the maximally supersymmetric Yang-Mills theory, where $\beta_g=0$. In this case, the correction proportional to $\beta$ in the left-hand-side of \eq{master2} is absent, and it is not necessary to compute the one-loop finite terms of $F(2)$.
\end{enumerate}
The two-loop master formula \eq{master2} is powerful, as it allows to compute two-loop RG effects with only one-loop information as an input. This is analogous to the computation of one-loop RG effects from tree-level amplitudes with \eq{oneloop}.

\subsection{Infrared renormalization of $\varphi^4$ theory up to two loops}

\begin{figure}[t]
\centering
\begin{tikzpicture}[line width=1.1 pt, scale=.7, baseline=(current bounding box.center)]

	\begin{scope}[shift={(7.0,0)}]
	\draw (-1.5,0) to (-.2,1) ;
	\draw (-1.5,0) to (-.2,-1) ;
	\draw (1.5,.5) to  (.2,1) ;
	\draw (1.5,-.5) to  (.2,-1) ;
	\draw (1.5,.5) to  (2.8,1) ;
	\draw (1.5,-.5) to  (2.8,-1) ;
	\node at (-1.5,0) {$\times$} ;
	\node at (-2.2,0) {${\cal J}$} ;
	\draw[dotted] (0,1.4) -- (0,-1.4) ;
	\draw[fill=gray!15] (1.5,0) ellipse (.7 cm and .9 cm) ;
	\draw[fill=white] (1.5,0) ellipse (.21 cm and .27 cm) ;
	\node at (0,-2) {${S}^{(1)}_{2,2}F_{0}(2)$} ;
	\end{scope}

	\begin{scope}[shift={(15.0,0)}]
	\draw (-1.5,0) to (-.2,0) ;
	\draw (-1.5,.5) to  (-.2,1) ;
	\draw (-1.5,-.5) to  (-.2,-1) ;
	\draw (1.5,0) to (.2,0) ;
	\draw (1.5,.5) to  (.2,1) ;
	\draw (1.5,-.5) to  (.2,-1) ;
	\draw (1.5,.5) to  (2.8,1) ;
	\draw (1.5,-.5) to  (2.8,-1) ;
	\draw[fill=gray!15] (-1.5,0) ellipse (.7 cm and .9 cm) ;
	\node at (-2.2,0) {$\times$} ;
	\node at (-2.9,0) {${\cal J}$} ;
	\draw[dotted] (0,1.4) -- (0,-1.4) ;
	\draw[fill=gray!15] (1.5,0) ellipse (.7 cm and .9 cm) ;
	\node at (0,-2) {$S_{2,3}F(3)$} ;
	\end{scope}

 \end{tikzpicture}
 \caption{\emph{Two-loop diagrams entering the right-hand-side of \eq{master2}.}}
 \label{fig:2loop}
 \end{figure}

As an example of the method, we now consider the simplest theory of a neutral scalar interacting through a quartic interaction, such that, in standard amplitude normalization, we have up to one loop
\begin{equation}
{\cal A}_{2,2} = -\left(\lambda+\delta^{(1)}_\lambda\right)-\frac{\lambda^2}{32\pi^2}\left( \ln \frac{-s}{\mu^2}+ \ln \frac{-t}{\mu^2}+ \ln \frac{-u}{\mu^2}\right)\,,
\end{equation}
with $\delta^{(1)}_\lambda$ accounting for a possible one-loop rational contribution which, however, will not play a role in the following.

The relevant form factor is instead defined by
\begin{equation}
F(2)\equiv \,  _{\rm out}\langle 1_\varphi,2_\varphi |T^{\mu\nu}(0) |0\rangle\,,
\end{equation}
with the further implicit assumption that we project $F(2)$ onto its $J=2$ component, for example by taking the outgoing particles to be in some $d$-wave state. This extra care is required in order to deal with the fact that $T^{\mu\nu}$ can and in general \emph{does} mix in the UV. However, to our understanding, all mixing effects happen in the $J=0$ sector. This is because there is no $J=2$ identically conserved operator with which $T^{\mu\nu}$ can mix, unless one turns on gravity. In this case, the energy-momentum tensor can mix with $\partial_\rho \partial_\sigma C^{\mu\rho\nu\sigma}$ (where $C$ is the Weyl tensor), which is purely $J=2$. However this mixing has to vanish in the limit $M_P \to \infty$, i.e. when gravity is neglected. Notice that, if the renormalization of the form factor were not multiplicative, \eq{factorizable} could not have been invoked.

We want now to compute the right-hand-side of \eq{master2} in the $\varphi^4$ theory. We first observe that $S_{2,3}=0$ in this case, so that the last term of \eq{master2} is zero. Second, to greatly simplify the two-particle phase space integral, we note that the free-theory form factor can be written, in angular coordinates, as
\begin{equation}
F_0(2)=F_0(\theta,\phi)=\sum_{m=-2}^2 \tau_{m}^{\mu\nu}d^2_{m,0}(\theta)e^{im\phi}\,,
\end{equation}
where the explicit form of $\tau_m$ is not relevant. The angles $\theta$ and $\phi$ parametrize the direction of the outgoing $\varphi$ particles with respect to some reference frame. The appearance of just the $d^2$ Wigner functions is a consequence of having projected out the $J=0$ component. Notice also that $F_0(0,0)=\tau_0^{\mu\nu}$. When we apply \eq{master2} to the problem at hand, we get (the factor 1/2 in the phase space integral is for identical internal particles)
\begin{align}\label{scalar}
e^{i\pi \left(\gamma_{\rm IR}-\frac{1}{\lambda}a^* \beta_\lambda^{(1)}\right)}\tau_0^{\mu\nu}&=\tau_0^{\mu\nu}+\frac{i}{64\pi^2}\int_{4\pi}s_\theta d\theta \,d\phi \,{\cal A}_{2,2}(\theta) \,F_0(\theta,\phi)\nonumber \\
&=\tau_0^{\mu\nu}\left(1+\frac{i}{32\pi}\int_0^\pi s_\theta d\theta\,{\cal A}_{2,2}(\theta)d^2_{0,0}(\theta) \right) \nonumber \\
&=\tau_0^{\mu\nu}\left(1+\frac{i}{32\pi}\int_0^\pi s_\theta d\theta\,\left[K-\frac{\lambda^2}{16\pi^2}\ln s_\theta\right]d^2_{0,0}(\theta) \right) \nonumber \\
&=\tau_0^{\mu\nu}\left(1+\frac{2i\pi\lambda^2}{3(32\pi^2)^2}\right).
\end{align}
In the third equality we introduced
\begin{equation}
K\equiv -\lambda-\delta^{(1)}_\lambda-\frac{3\lambda^2}{32\pi^2}\ln\frac{-s}{\mu^2}\,.
\end{equation}
Since $K$ has no angular dependence, it is annihilated by the $d^2$ phase space integral. Notice that this includes the one-loop contribution to \eq{master2} (coming from the tree-level ${\cal A}_{2,2}$), the two-loop contribution coming from the one-loop rational term of $\cal A$, i.e. $-\delta^{(1)}_\lambda$, and finally the two-loop contribution coming from the $\ln\mu$ term of $\cal A$.
Before extracting the result, let us make some remarks about \eq{scalar}.
\begin{enumerate}
\item Perhaps the most important point is that there is a certain degree of ``bootstrapping'' in \eq{scalar}. First of all, in the fact that both sides of the equation are proportional to the same (free-theory) tensor $\tau_0^{\mu\nu}$.
\item More interestingly, we know that the physics that feeds $a^*$ includes contributions to $\gamma_{\rm IR}$ at one loop, which are divergent (see \eq{aeq}). If these were non-zero (and since $\beta^{(1)}_\lambda$ does not vanish in this theory), we would need similar divergent two-loop contributions on the right-hand-side. But we find just a number. Consistently, we see that $a$ is indeed finite, since the one-loop contribution in the right-hand-side is zero (it is all in $K$). In this respect, $\varphi^4$ is special.
\item Since $a^*$ contains a finite, scheme dependent term, the two-loop anomalous dimension will in general be scheme dependent as well. The scheme dependence is precisely dictated by \eq{scalar}, and in general by \eq{master2}. However, as we said, by expanding \eq{scalar} we get, at one loop, $\gamma_{\rm IR}^{(1)}=0$. Since there is no need for a renormalization at one loop (there are no UV nor IR divergences) $a\equiv a_{\rm fin}^{(1)}$ is physical. In order not to spoil the universal properties of the energy-momentum tensor, it must vanish.
\end{enumerate}
We therefore get at two loops
\begin{equation}
\gamma_{\rm IR}^{(2)}=\frac{2\lambda^2}{3(32\pi^2)^2}\,.
\end{equation}
Since we obtain a pure number, we rediscover that only collinear divergences are present at two loops in $\lambda \varphi^4$.

\subsection{Real and imaginary part of the two-loop equation}\label{SYM}

\begin{figure}[t]
\centering
\begin{tikzpicture}[line width=1.1 pt, scale=.7, baseline=(current bounding box.center)]

	\begin{scope}[shift={(5.4,-4)}]
	\draw (-1.5,0) to (-.2,1) ;
	\draw (-1.5,0) to (-.2,-1) ;
	\draw (1.5,.5) to  (.2,1) ;
	\draw (1.5,-.5) to  (.2,-1) ;
	\draw (1.5,.5) to  (2.8,1) ;
	\draw (1.5,-.5) to  (2.8,-1) ;
	\node at (-1.5,0) {$\times$} ;
	\node at (-2.2,0) {${\cal J}$} ;
	\draw[dotted] (0,1.4) -- (0,-1.4) ;
	\draw[fill=gray!15] (1.5,0) ellipse (.7 cm and .9 cm) ;
	\draw (4.5,.5) to  (3.2,1) ;
	\draw (4.5,-.5) to  (3.2,-1) ;
	\draw (4.5,.5) to  (5.8,1) ;
	\draw (4.5,-.5) to  (5.8,-1) ;
	\draw[dotted] (3,1.4) -- (3,-1.4) ;
	\draw[fill=gray!15] (4.5,0) ellipse (.7 cm and .9 cm) ;
	\node at (1.6,-1.9) {${\cal M}_{2,2}^{(0)}\, {\cal M}_{2,2}^{(0)} \,F_0(2)$} ;
	\node at (-0.4,1.35) {\footnotesize{$\hat{n}'$}} ;
	\node at (2.6,1.35) {\footnotesize{$\hat{n}$}} ;
	\node at (5.5,1.32) {\footnotesize{$\hat{z}$}} ;
	\end{scope}

 \end{tikzpicture}
 \caption{\emph{Class of diagrams entering into the real part of \eq{master2} at two loops. On general grounds, they do not provide genuinely two-loop information, as they are $({\rm one-loop})^2$. These diagrams can have `cusp' divergences in each of the two cuts, requiring a double subtraction. See \eq{doubleint}.}}
 \label{fig:1loop^2}
 \end{figure}

Let us now consider another example which, from the point of view of \eq{master2}, is structurally simple. The maximally supersymmetric Yang-Mills theory in the limit $N_c\to\infty$ enjoys a number of properties, including (i) the fact that $\beta_g=0$ to all orders in perturbation theory, (ii) that there is only one species (or better multiplet) and finally (iii) that $\gamma_{\rm coll}^{(1)}=0$. This means that the two-loop order in \eq{master2} simplifies to
\begin{equation}
\left[i\pi \gamma_{\rm IR}^{(2)}-\frac{\pi^2}{2}\left(\gamma_{\rm IR}^{(1)}\right)^2\right]F_0(2)=i{\cal M}_{2,2}^{(1)}F_0(2)+i{\cal M}_{2,3}F(3)\,.
\end{equation}
As explained in section 2.3 of \cite{Bern:2020ikv}, the two-loop equation is in general complex and can be further split into two. We get in this case
\begin{align}
\pi \gamma_{\rm IR}^{(2)}\,\,F_0(2)&={\rm Re} {\cal M}_{2,2}^{(1)} \,F_0(2)+{\cal M}_{2,3}F(3)\,,\label{real}\\
\frac{\pi^2}{2}\left(\gamma_{\rm IR}^{(1)}\right)^2 F_0(2)&={\rm Im} {\cal M}_{2,2}^{(1)} \,F_0(2)\,.\label{imaginary}
\end{align}
We see that the informative equation is the first one, as the second one contains only one-loop quantities. However, the consistency of the formalism requires it to be identically true. Since \eq{imaginary} provides a non-obvious identity, we study in detail how it happens to be true. The key observation is that, by virtue of the optical theorem, imaginary parts of loop amplitudes are computed by cuts. In formulas, we have
\begin{equation}
{\rm Im} {\cal M}_{2,2}^{(1)}=\frac{1}{2}{\cal M}_{2,2}^{(0)}{\cal M}_{2,2}^{(0)}\,,
\end{equation}
where, as usual, in the right-hand-side of the equation a phase space integration over the intermediate state is understood, and the (crucial) factor 1/2 is dictated by the optical theorem \cite{Bern:2020ikv}. All in all, we find
\begin{equation}\label{doublesub}
\left(\pi\gamma_{\rm IR}^{(1)}\right)^2 F_0(2)={\cal M}_{2,2}^{(0)}\, {\cal M}_{2,2}^{(0)} \,F_0(2)\,,
\end{equation}
with a double phase space integration in the right-hand-side. How do we make sense of this equation? Remember that, on general grounds, the interpretation of the `cusp' contribution appearing in the dilation operator $D$ is in terms of a divergent integral, like in \eq{gammaIR}. What is the meaning of its square?
In light of the consistency of the master formula, we argue that it must be an integral being capable of absorbing the divergences of the right-hand-side of \eq{doublesub}, the divergences coming from the forward region of ${\cal M}_{2,2}^{(0)}$. If we parametrize with $\hat{n}$ ($\hat{n}'$) the direction of the two-particle state in the most external (internal) cut of Fig.~\ref{fig:1loop^2}, and we call $\hat{z}$ the direction of the external two-particle state (we are in the center of momentum frame of the system, where each pair of particles is back-to-back), we find that the divergences of the cut in Fig.~\ref{fig:1loop^2} can only be absorbed by an integral of the following form
\begin{equation}\label{doubleint}
(\gamma_{\rm IR}^{(1)})^2=(\gamma_{\rm cusp}^{(1)})^2\int d{\hat{n}} \,\frac{1}{1-\hat{z}\cdot\hat{n}}\,\int d{\hat{n}'}\, \frac{1}{1-\hat{n}\cdot\hat{n}'}\,,
\end{equation}
where, in the same notation, the divergent integral in \eq{gammaIR} would read (we set $\int d\hat{n}=2$)
\begin{equation}
\gamma_{\rm IR}^{(1)}=\gamma_{\rm cusp}^{(1)}\int d{\hat{n}} \,\frac{1}{1-\hat{z}\cdot\hat{n}}\,.
\end{equation}
In a general theory, the structure of the (one-loop)$^2$ cancellation can be more intricate, with contributions coming from $(\gamma_{\rm coll}^{(1)})^2$, that in super-Yang-Mills is zero, and also mixed ones.

\vspace{.5cm}

From a practical point of view, it is surely more interesting to say something about the structure of the genuinely two-loop contributions, i.e. \eq{real}, before concluding the section. \eq{gammaIR} is valid at all orders. Therefore we know that the right-hand-side of \eq{real} must necessarily take the form of a rational contribution, plus possibly a term with the same divergences as the phase space integral in the left-hand-side of \eq{gammaIR}. Alternatively, if one prefers to dimensionally regulate the cuts, the result of the phase space integral must have the form of \eq{gammaIRdimreg}. This is quite a restrictive statement on the right-hand-side of \eq{real}. Moreover, as we explained at the end of section \ref{1loop}, the two contributions to $\gamma_{\rm IR}$ (i.e. $\gamma_{\rm cusp}$ and $\gamma_{\rm coll}$) can be structurally distinguished, allowing one to project out the cusp contribution (as we said, this is especially useful if the focus is on the computation of UV divergences in the on-shell scheme). The concrete implementation of these ideas is deferred to future work.

\subsection{Cancellation of non-local terms (Yukawa theory)}

In this section we would like to consider an example which, from the point of view of \eq{master2}, is slightly more complicated, in that both $\beta^{(1)}$ and $a$ are different from zero, with $a$ being infinite due to the presence of collinear divergences at one loop: the Yukawa theory of a scalar $\phi$ and a fermion $\psi$. In particular, we want to study some of the two-loop aspects of the form factor $F_\phi(2)\equiv\, _{\rm out}\langle 1_\phi2_\phi|T^{\mu\nu}(0) |0\rangle$.

Without going into the full computation, we want to test the consistency of the two-loop equation. In particular, we want to show how certain non-local terms proportional to $\ln(-s/\mu^2)$ cancel among themselves in \eq{master2}. These are terms that can be recognized as (one-loop)$^2$ contributions, but cannot be `projected out' by simply taking the imaginary part of \eq{master2}, like we did in Sec.~\ref{SYM}. These terms can appear {explicitly} in \eq{master2} at two loops (e.g. through $a$ or $S^{(1)}$) or in phase space integrals like the one in Fig.~\ref{fig_coll}, and should not be confused with the genuinely two-loop logarithms which have, instead, to be interpreted as contributions to $\gamma_{\rm cusp}$ (cf. Sec.~\ref{1loop}).

By inspection of \eq{master2}, and knowing that
\begin{equation}
F_\phi(2)=\left(1-\frac{1}{2}\gamma^{(1)}_{\phi}\ln \frac{-s}{\mu^2}+\ldots\right)F_{0,\phi}(2)\,,
\end{equation}
and similarly for $\psi$ (the ellipses stand for higher-order and finite contributions, whilst $\gamma_{\rm cusp}^{(1)}=0$), we find several places where terms proportional to $\ln \mu$ appear. All such terms are of second order in the loop expansion. On the left-hand-side of \eq{master2}, we find
\begin{equation}
i\pi \frac{1}{2}\beta_y^{(1)}\partial_y \gamma_\phi^{(1)}\ln\frac{-s}{\mu^2} F_{0,\phi}=i\pi \frac{1}{y}\beta_y^{(1)}\gamma_\phi^{(1)}\ln\frac{-s}{\mu^2} F_{0,\phi}\,,
\end{equation}
where we used that $\gamma_\phi^{(1)}\propto y^2$. On the right-hand-side, upon expanding the exponential, we find the term $(a'-a)^*i{\cal M}_{2,2'}^{(0)}F_0(2')$ which, rewritten in standard amplitude normalization, yields
\begin{equation}
\frac{\gamma_\phi^{(1)}-\gamma_\psi^{(1)}}{2}\ln\frac{-s}{\mu^2}\times\frac{i}{32\pi^2}\int_{4\pi} d\Omega \,{\cal A}^{(0)}_{\phi\phi\to\psi\psi}\,F_{0,\psi}=\frac{\gamma_\phi^{(1)}-\gamma_\psi^{(1)}}{2}\ln\frac{-s}{\mu^2}\times i\pi\gamma_\phi^{(1)}F_{0,\phi}\,,
\end{equation}
where we invoked \eq{oneloop} (and the specific properties of the model) to equate the phase space integral to the collinear anomalous dimension of $\phi$. Another possible source for terms containing $\ln\mu$ in the right-hand-side of \eq{master2} comes from the convolution of the relevant one-loop $S$ matrix elements with $F_{0,\psi}$ and $F_{0,\phi}$. In the model at hand, the relevant $S$ matrix elements are respectively $S^{(1)}_{\phi\phi\to\psi\psi}$ and $S^{(1)}_{\phi\phi\to\phi\phi}$. Regarding the second one, its $\ln\mu$ pieces are in the $s$-wave and do not survive the integration against $T^{\mu\nu}$. The first one, instead, has $\ln\mu$ contributions that enter non-trivially. Thanks to the RG equation that the $S$ matrix elements need to satisfy, we have\footnote{For amplitudes that are defined by minimal subtraction in the UV and in the IR, like the ones we consider here, the appropriate RG equation reads
\begin{equation}
\dot{S}=\left(\beta_g\partial_g+\mu\partial_\mu-\gamma_{\rm IR}\right)S=0\,.
\end{equation}}
\begin{equation}
{\cal A}^{(\ln\mu)}_{\phi\phi\to\psi\psi}=\left(\frac{\beta^{(1)}_y}{y}-\frac{\gamma_\phi^{(1)}+\gamma_\psi^{(1)}}{2}\right)\ln\frac{-s}{\mu^2}{\cal A}^{(0)}_{\phi\phi\to\psi\psi}\,,
\end{equation}
so that all in all we find a contribution to the right-hand-side of \eq{master2} that reads
\begin{equation}
\frac{i}{32\pi^2}\int_{4\pi}d\Omega\,{\cal A}^{(\ln\mu)}_{\phi\phi\to\psi\psi}\,F_{0,\psi}=\left(\frac{\beta^{(1)}_y}{y}-\frac{\gamma_\phi^{(1)}+\gamma_\psi^{(1)}}{2}\right)\ln\frac{-s}{\mu^2}\times i\pi\gamma_\phi^{(1)}F_{0,\phi}\,.
\end{equation}
By comparing the left- and right-hand-sides, we find remarkably that there are many cancellations, and one is left only with
\begin{equation}\label{countert}
-i\pi\gamma_\phi^{(1)}\gamma_\psi^{(1)}\ln\frac{-s}{\mu^2}F_{0,\phi}
\end{equation}
on the right-hand-side of \eq{master2}. Consistency requires that this non-local piece should cancel eventually. By a more careful study of the model, we indeed find a diagram for which \eq{countert} provides a suitable ``counterterm''. It is given by the last term of \eq{master2} which, in this model, schematically reads $S_{\phi\phi\to\phi\psi\psi}F_{\phi\psi\psi}$. As a matter of fact, such three-particle phase space integral features, as we explicitly checked, collinear $\psi$ divergences coming from diagrams like the one in Fig.~\ref{fig_coll}.

\begin{figure}[t]
\centering
\begin{tikzpicture}[line width=1.1 pt, scale=0.9, baseline=(current bounding box.center)]

\begin{scope}[shift={(7.0,0)}]
	\draw (-1.5,0) to (-.1,1) ;
	\draw (-1.5,0) to (-.1,-1) ;
	\draw (1.5,0) to (.1,1) ;
	\draw (1.5,0) to (.1,-1) ;
	\draw[dashed]  (1.75,0) to  (3.15,1) ;
	\draw[dashed]  (1.75,0) to  (3.15,-1) ;
	\draw[dashed]  (-.8,.5) to  (-.1,.05) ;
	\draw[dashed]  (.8,.5) to  (.1,.05) ;
	\node at (-1.5,0) {$\times$} ;
	\node at (-2.2,0) {${\cal J}$} ;
	\draw[dotted] (0,1.4) -- (0,-1.4) ;
	\draw[fill=gray!15] (1.625,0) ellipse (.8 cm and .7 cm) ;

	\node at (0,-2) {${S}_{2,3}F(3)$} ;
	\end{scope}

 \end{tikzpicture}
 \caption{\emph{Diagram with infrared divergences in Yukawa theory, coming from the collinear emission of a scalar $\phi$ (dashed line) from $\psi$ (solid line).}}
 \label{fig_coll}
 \end{figure}

By working in dimensional regularization we find that, in accordance with the general theorems on IR singularities, the relevant phase space integral behaves as
\begin{equation}
\frac{i}{32\pi^2}\int_{4\pi}d\Omega \,{\cal A}^{(0)}_{\phi\phi\to\psi\psi}\,F_{0,\psi} \times \gamma_{\psi}^{(1)}\ln\frac{-s}{\mu^2}=i\pi\gamma_\phi^{(1)}\gamma_\psi^{(1)}\ln\frac{-s}{\mu^2}F_{0,\phi}\,,
\end{equation}
which is precisely canceled by the counterterm \eq{countert}.

\section{Conclusions}

\eq{chw} is an all order statement that, when expanded in loops, reveals to be extremely powerful, in that it allows to reconstruct the RG at $n$-loop order in terms of ($n-$1)-loop objects (in a sharply defined way). Like any computational technique, it requires practice to become operative. The main purpose of this work is to begin a discussion on the two-loop zoology of \eq{chw} when IR divergences are present, to show its consistence and disclose its simplicity.

Since UV and IR divergences are on the same ground in dimensional regularization (we consider here scaleless theories), one needs a way to disentangle them. For this purpose, we followed the path that was proposed in the literature of taking $\cal O$ to be conserved, to ensure that $\gamma_{\cal O}\equiv\gamma_{\rm UV}=0$ (strictly speaking, one has to also make sure that there are no identically conserved operators with the same quantum numbers of $\cal O$ \cite{Collins:2005nj}).

After some practice, one comes to agree with the experts on the statement that IR divergences are simple. In the language of \eq{chw}, this is essentially because they do not mix different form factors (up to color mixing in Yang-Mills). This property was the key for a simplification of the master formula at two-loops: in practice, for what concerns the two-particle cuts, we are instructed to compute only phase space integrals with the free theory (or tree-level) form factor $F_0$ (while, in general, one-loop $S$ matrix elements produce genuinely two-loop contributions).

We tried to let \eq{chw} talk. One (perhaps marginal, technical) suggestion that we captured is that we should better perform the phase space integrals in dimensional regularization, consistently with what we do for $S$ and $F$. This allows to make sense of some apparent ambiguities which are encountered in a straightforward interpretation of the master equation. See for example the discussion about \eq{gammaIRdimreg} and \eq{gammaIR}.

The discussion was intentionally kept at the level of the structural analysis of \eq{master2}, limiting explicit computations to the minimum necessary. This was according to the belief that the main difficulties in computing two-loop effects with on-shell methods comes from the presence of spurious (one-loop$)^2$ effects. The main purpose here was to show how these have a clear structure and can be often easily subtracted following general principles. Getting to the end of the computation is of course also useful to fully taste the flavor of the method, but this is postponed to future work.

The method that we explored is completely general. This is a common feature in the on-shell approach, and the key to its versatility. Of course, some model-dependent property can make these techniques even more powerful, allowing them to be efficiently extended to higher loops (clearly, there is no intrinsic limitation in generalizing the formulas in this direction).

Even though what was reported here is conceptually independent from the spinor-helicity formalism \cite{Cheung:2017pzi}, the practical implementation of the method greatly benefits from it. The connection between form factors and massive amplitudes is also illuminating and computationally useful \cite{Arkani-Hamed:2017jhn}.

Knowing the structure of IR divergences is necessary if one needs to subtract them, like in the study of the UV running of low energy EFTs in the on-shell massless scheme \cite{EliasMiro:2020tdv,Baratella:2020lzz,Baratella:2021guc,Jiang:2020mhe}. For what concerns the study of EFTs (and especially of the Standard Model EFT), a lot of recent activity was dedicated to the identification of amplitude bases for higher-dimensional operators \cite{Shadmi:2018xan,Ma:2019gtx,Aoude:2019tzn,Durieux:2019siw,AccettulliHuber:2021uoa}, and the study of various kinds of selection rules for mixing among operators, like helicity and angular momentum selection rules \cite{Cheung:2015aba,Craig:2019wmo,Jiang:2020rwz}. These sets of ideas, put together with the methods that were reported here, can be of great practical use for the systematic study of the RG at two loops \cite{Bern:2020ikv}.

The zoology of \eq{master2} is richer than what was reported in this work, but it surely is ``finite''. We leave for the future a more complete investigation of the formula, with its implications and simplifying strategies.

\subsection*{Acknowledgements}
This research has been partially supported by the DFG Cluster of Excellence 2094 ORIGINS and the Collaborative Research Center SFB1258. The author warmly thanks, for the hospita-lity, the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP), which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy -- EXC-2094 -- 390783311. The author is also greatly indebted to the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA$^+$, for providing support and, especially, an extremely stimulating atmosphere during the workshop ``Amplitudes Meet BSM 2022''.