aainstitutetext: Graduate School of China Academy of Engineering Physics, No.
10 Xibeiwang East Road, Haidian District, Beijing, 100193, China

# Five-loop anomalous dimensions of $\phi^{Q}$ operators in a scalar theory
with $O(N)$ symmetry

Qingjun Jin a Yi Li<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We compute the complete $Q$-dependence of anomalous dimensions of traceless
symmetric tensor operator $\phi^{Q}$ in $O(N)$ scalar theory to five-loop. The
renormalization factors are extracted from $\phi^{Q}\rightarrow Q\phi$ form
factors, and the integrand of form factors are constructed with the help of
unitarity cut method. The anomalous dimensions match the known results in
Badel:2019oxl ; Antipin:2020abu , where the leading and subleading terms in
the large $Q$ limit were obtained using a semiclassical method.

##  1 Introduction

Different physical systems around second order phase transitions exhibit
universal critical behaviors, and can be studied through field-theoretic
renormalization group approach Wilson:1971bg ; Wilson:1971dh . The critical
exponents can be obtained from the scaling dimensions of operators at the
Wilson-Fisher fixed point Wilson:1971dc  combined with resummation methods
LeGuillou:1979ixc . One of the most extensively studied quantum field theories
in this context is the scalar $\phi^{4}$ theory with $O(N)$ symmetry, which
describes critical phenomena of different statistical systems including
superfluid $\sideset{{}^{4}}{}{\mathop{He}}$, ferromagnets, etc.

Recently, a semiclassical method (see e.g. Son:1995wz ; Hellerman:2015nra ;
Alvarez-Gaume:2016vff ; Alvarez-Gaume:2019biu ; Arias-Tamargo:2019xld ;
Gaume:2020bmp ) was developed in Badel:2019oxl , which provides an convenient
way to study the anomalous dimensions of $\phi^{Q}$-type charged operators in
$\phi^{4}$ theory with $U(1)$ symmetry. The leading and subleading terms of
the anomalous dimensions in the large $Q$ limit can be evaluated to
arbitrarily high loops. The method was applied to the $O(N)$ scalar theory in
Antipin:2020abu , and the result was consistent with the traditional
perturbative calculation to 3-loop. The comparison was extended to 4-loop in
Jack:2021ypd , where the subset of Feynman diagrams which contribute to the
anomalous dimensions in the large charge limit were evaluated. Similar
computations and comparisons have also been performed for a $U(N)\times U(N)$
model in Antipin:2020rdw ; Jack:2021lja , and for a $U(N)\times U(M)$ model in
Antipin:2021akb .

In order to further test the anomalous dimensions from semiclassical method,
the perturbative computations should be extended to five and more loops. In
the perturbative approach, the anomalous dimensions are extracted from the
ultraviolet (UV) divergences of correlation functions or form factors.
However, there are two major challenges in this computation. First, the number
of contributing Feynman diagrams increase drastically as the loop number and
$Q$ increase. Second, it can be very difficult to evaluate the UV divergences
of 5-loop integrals in the form factors.

The first challenge can be overcome by using the unitarity cut method, which
has been playing a key role in the computation of gauge and gravity scattering
amplitudes and form factors (see e.g. Bern:1994zx ; Bern:2008qj ; Boels:2008ef
; Bern:2012uc ; Cachazo:2013hca ; Cachazo:2013iea ; Bern:2015ooa ;
Yang:2016ear ). The form factor with arbitrary $Q$ can be constructed by a few
simple "irreducible" cut integrands. The UV divergences are computed using the
UV decomposition method Jin:2022zcj , which direct extracts the UV
counterterms from the integrand.

In this work, we perform a perturbative computation of the anomalous
dimensions of $\phi^{Q}$-type operators in the $O(N)$ model to 5-loop. In
Section 2, we construct the form factor integrands using unitarity cut method.
In Section 3, we present the full $Q$-dependence of anomalous dimensions to
5-loop, and compare with known results. In Section 4, we discuss some possible
future directions.

##  2 Construction of $\phi^{Q}\rightarrow Q\phi$ form factors using unitarity
cut

The $\phi^{4}$ theory with $O(N)$ symmetry is given by the Lagrangian,

$\mathcal{L}=\frac{1}{2}(\partial_{\mu}\phi)^{2}+\frac{g}{4!}(\phi^{2})^{2}\
,$ (1)

in which $\phi=(\phi_{1},\cdots,\phi_{N})$ is a real scalar field with $N$
components. We will work in a $D=4-\epsilon$ dimensional Euclidean space111As
far as $\phi^{4}$ theory is concerned, there is no essential difference
between Euclidean space and Minkowski space, since the theory is free of IR
divergence., and use the minimal subtraction scheme. Following Antipin:2020abu
, we consider the symmetric traceless component of $\phi^{Q}$ operator:

$\displaystyle(\phi^{2})_{i_{1}i_{2}}$
$\displaystyle\equiv\phi_{i_{1}}\phi_{i_{2}}-\frac{1}{N}\delta_{i_{1}i_{2}}\phi^{2},\
$ (2) $\displaystyle(\phi^{3})_{i_{1}i_{2}i_{3}}$
$\displaystyle\equiv\phi_{i_{1}}\phi_{i_{2}}\phi_{i_{3}}-\frac{1}{N+2}\phi^{2}(\phi_{i_{1}}\delta_{i_{2}i_{3}}+\phi_{i_{2}}\delta_{i_{1}i_{3}}+\phi_{i_{3}}\delta_{i_{1}i_{2}}),\
$ $\displaystyle\cdots$

The anomalous dimensions of $\phi^{Q}$ operators with $Q=2,4$ were obtained in
Braun:2013tva ; Kompaniets:2019zes  and Calabrese:2002bm  respectively, and in
this paper we consider operators with general $Q$. The anomalous dimensions of
these operators can be extracted from the UV divergences of the following form
factor,

$\Bigl{(}\mathbf{F}_{Q}\Bigr{)}_{i_{1}\cdots i_{Q}}^{j_{1}\cdots
j_{Q}}\equiv\int d^{D}xe^{-ik\cdot
x}\Bigl{\langle}\phi_{j_{1}}(p_{1})\cdots\phi_{j_{Q}}(p_{Q})\Bigr{|}(\phi^{Q})_{i_{1}\cdots
i_{Q}}(x)\Bigr{|}0\Bigr{\rangle}\ .$ (3)

We will call it the $\phi^{Q}\rightarrow Q\phi$ form factor.

In Section 2.1, we study the structure of $O(N)$ indices, and show that these
indices can be removed by contracting the operator and external states with
the same zero-norm vector. In Section 2.2, we construct form factor integrands
with the help of unitarity cut method. A few "irreducible" cut integrands are
taken as building blocks of form factor integrands with general $Q$. In
Section 2.3, we discuss the scattering amplitudes which are required in the
unitarity cut computations.

###  2.1 $O(N)$ tensor structures

The first difficulty we encountered in the computation of these form factors
comes from the cumbersome $O(N)$ tensor structure in the operator. In the
simplest case $Q=2$, the Feynman rule already contains 3 terms,

$V_{i_{1}i_{2}}^{j_{1}j_{2}}=\delta_{i_{1}}^{j_{1}}\delta_{i_{2}}^{j_{2}}+\delta_{i_{1}}^{j_{2}}\delta_{i_{2}}^{j_{1}}-\frac{2}{N}\delta_{i_{1}i_{2}}\delta^{j_{1}j_{2}}\
.$ (4)

The number of terms in the Feynman rule increases exponentially as $Q$
increases.

In order to overcome this difficulty, we observe all Feynman diagrams
contributing to the $\phi^{Q}\rightarrow Q\phi$ form factor actually have the
same tensor structure, which equals the corresponding Feynman rule:

$\Bigl{(}\mathbf{F}_{Q}\Bigr{)}_{i_{1}\cdots i_{Q}}^{j_{1}\cdots
j_{Q}}=V_{i_{1}\cdots i_{Q}}^{j_{1}\cdots
j_{Q}}\mathcal{F}_{Q}(p_{1},\cdots,p_{n})\ ,$ (5)

in which $\mathcal{F}_{Q}(p_{1},\cdots,p_{n})$ is a loop integrand without
$O(N)$ indices. The factorization of the tensor structure is basically a
consequence of charge conservation, but it can also be understood on the level
of Feynman diagrams. A proof of this factorization based on Feynman diagrams
can be found in Appendix A.

This factorization allows us to contract the operator with a zero-norm complex
vector $q$ satisfying $q^{2}=0$ and $q\cdot\bar{q}=1$, without altering the
anomalous dimension222The most simple choice of $q$ is
$\frac{1}{\sqrt{2}}(1,i,0,\cdots,0)$. In the $U(1)$ case ($N=2$),
$\mathcal{O}_{Q}\propto(\phi_{1}+i\phi_{2})^{Q}$ is exactly the operator
$\phi^{n}$ ($n$ corresponds to $Q$ here) in Badel:2019oxl .:

$\mathcal{O}_{Q}\equiv\frac{1}{Q!}q_{i_{1}}\cdots
q_{i_{Q}}(\phi^{Q})_{i_{1}\cdots i_{Q}}=\frac{(q\cdot\phi)^{Q}}{Q!}\ .$ (6)

The Feynman rule of $\mathcal{O}_{Q}$ only contains a single term. The form
factor of $\mathcal{O}_{Q}$ is

$\int d^{D}xe^{-ik\cdot
x}\Bigl{\langle}\phi_{j_{1}}(p_{1})\cdots\phi_{j_{Q}}(p_{Q})\Bigr{|}\mathcal{O}_{Q}(x)\Bigr{|}0\Bigr{\rangle}=\bar{q}_{j_{1}}\cdots\bar{q}_{j_{Q}}\mathcal{F}_{Q}(p_{1},\cdots,p_{Q})\
,$ (7)

in which $\bar{q}$ appears instead of $q$ because the operator is "in-coming"
instead of "out-going".

The expression can be further simplified if we also contract the external
fields $\phi_{j}$ with $q$, and define $\varphi\equiv q\cdot\phi$. The form
factor is now free of $O(N)$ indices,

$\int d^{D}xe^{-ik\cdot
x}\Bigl{\langle}\varphi(p_{1})\cdots\varphi(p_{Q})\Bigr{|}\mathcal{O}_{Q}(x)\Bigr{|}0\Bigr{\rangle}=\mathcal{F}_{Q}(p_{1},\cdots,p_{Q})\
.$ (8)

We will call it the $\varphi^{Q}\rightarrow Q\varphi$ form factor. The
$L$-loop contribution to the $\varphi^{Q}\rightarrow Q\varphi$ form factor
will be denoted by $\mathcal{F}_{Q}^{(L)}$. The tree form factor is
simply333Strictly speaking, all form factors contain a
$\delta^{D}(k-\sum_{i=1}^{Q}p_{i})$ from the $x$-integration, which guarantees
the momentum conservation. So $\mathcal{F}_{Q}^{(0)}=1$ actually means
$\mathcal{F}_{Q}^{(0)}=\delta^{D}(k-\sum_{i=1}^{Q}p_{i})$. For compactness we
will neglect this overall factor in all tree and loop form factors.
$\mathcal{F}_{Q}^{(0)}=1$.

###  2.2 Unitarity cut

Unitarity based methods simplify the computation of scattering amplitudes(form
factors) in two ways. First, the intermediate expressions during the
computation can be substantially simplified by making use of on-shell
amplitudes, instead of off-shell interaction vertices, as the building block
of loop amplitude integrand. Second, different cut channels can be related by
permutations of the internal and external legs, so a lot repeated computations
can be avoided.

Since in this work we are interested in a pure scalar theory, the on-shell
constraint does not really simplify the expressions, but unitarity cut method
can still help us to exploit the permutation symmetry of the internal and
external legs.

In each Feynman diagram contributing to the $\varphi^{Q}\rightarrow Q\varphi$
form factor, there is a single vertex which contains $\varphi^{Q}$, and we
will only cut all the internal legs which are directly incident to this
vertex444One may consider other types of cuts which affect legs not
originating from $\varphi^{Q}$, but the cuts we consider are sufficient to fix
the complete integrand. (in other words, all internal legs originating form
$\varphi^{Q}$). After the cut, the form factor is split into a
$\varphi^{Q}\rightarrow Q\varphi$ tree form factor, and one or more scalar
amplitudes. For example, Figure 1 shows a unitarity cut of the
$\varphi^{5}\rightarrow 5\varphi$ form factor. The building blocks of the cut
are a $\varphi^{5}\rightarrow 5\varphi$ tree form factor (the blob in the
middle), and two 4-point scalar tree amplitudes (the other two blobs). Figure
2 shows a unitarity cut which contains two building blocks: a
$\varphi^{Q}\rightarrow Q\varphi$ tree form factor, and a $2n$-point scalar
loop amplitude.

Figure 1: A reducible cut. The blue line is the $\varphi^{Q}$ operator. Each
solid blob represents a tree amplitude (or tree form factor), while a blob
with several holes represents a multiloop amplitude (or multiloop form
factor). Figure 2: An irreducible $n$-line cut of $\varphi^{Q}\rightarrow
Q\varphi$ $L$-loop form factor. The l.h.s. of the cut is a
$\varphi^{Q}\rightarrow Q\varphi$ tree form factor, the r.h.s. represents a
$(L-n+1)$-loop $2n$-point scattering amplitude.

If a unitarity cut contains two or more scalar amplitudes (besides the tree
form factor), we will call it a reducible cut. If the cut contains a single
scalar amplitudes, we will call it an irreducible cut. The cut in Figure 1 is
a reducible cut, while the cut in Figure 2 is an irreducible cut. The
integrand of a reducible cut can be split into multiple pieces. For example,
Figure 1 corresponds to the following cut integrand555As discussed in Section
2.1, the tree form factor equals 1. Therefore the cut integrand does not
depend on the number of external legs which are originating from
$\varphi^{Q}$.:

$I=I_{1}I_{2},\ I_{1}=\frac{g}{l_{1}^{2}(l+p_{12})^{2}},\
I_{2}=\frac{g}{l_{2}^{2}(l_{2}+p_{34})^{2}}\ .$ (9)

Both $I_{1}$ and $I_{2}$ are identical to a cut integrand of a one-loop
$\varphi^{2}\rightarrow 2\varphi$ form factor. In general, each reducible cut
can be regarded as the product of several lower loop irreducible cuts.

It is not hard to see that:

  1. 1.

Each Feynman diagram contains at least two internal legs which are originating
from $\varphi^{Q}$, so it can always be captured666A Feynman diagram is
captured by a cut, if all cut propagators are also propagators of this Feynman
diagram. On the contrary, if a cut propagator is not the propagator of a
Feynman diagram, then this Feynman diagram gives zero under the cut, which
means it is not captured by the cut. by a reducible or an irreducible cut.
Therefore these two types of cuts will suffice to fix the complete integrand.

  2. 2.

Different cuts have no intersections. If a Feynman diagram is captured by one
cut, then it cannot be captured by any other cuts.

This means the loop integrand $\mathcal{F}_{Q}^{(L)}$ can be obtained by
summing over contributions of all reducible and irreducible cuts. More
precisely, let $\\{\mathcal{C}_{1},\cdots\mathcal{C}_{N_{c}}\\}$ be the
complete set of reducible and irreducible cuts, and
$\\{l_{1},\cdots,l_{n_{i}}\\}$ be the momenta of cut legs in
$\mathcal{C}_{i}$. The cut integrand will denoted by
$\mathcal{F}_{Q}^{(L)}\Bigr{|}_{\mathcal{C}_{i}}$, and we define the restored
cut integrand as

$\mathcal{F}_{Q}^{(L)}\Bigr{|}_{\mathcal{C}_{i}}^{\text{restored}}=\frac{\mathcal{F}_{Q}^{(L)}\Bigr{|}_{\mathcal{C}_{i}}}{l_{1}^{2}\cdots
l_{n_{i}}^{2}}\ .$ (10)

The restored cut integrand is actually the contribution of all Feynman
diagrams captured by the cut. Therefore

$\mathcal{F}_{Q}^{(L)}=\sum_{i=1}^{N_{c}}\mathcal{F}_{Q}^{(L)}\Bigr{|}_{\mathcal{C}_{i}}^{\text{restored}}\
.$ (11)

Let us examine the irreducible cuts in Figure 2. The l.h.s. of the cut is
$\mathcal{F}_{Q}^{(0)}=1$, so the cut integrand does not depend on $Q$ or
$p_{n+1},\cdots,p_{Q}$. Therefore the corresponding restored cut integrand can
be denoted by $\mathtt{F}_{n}^{(L)}(p_{1},\cdots,p_{n})$. Similar cut channels
can be obtained by permuting all $Q$ external legs, and in all $C_{Q}^{n}$
different $\mathtt{F}_{n}^{(L)}(p_{i_{1}},\cdots,p_{i_{n}})$ are generated.
However, the UV counterterm is independent of the external momenta, because
the integral has divergence degree 0 (see Appendix B for details). So from now
on we will not discriminate $\mathtt{F}_{n}^{(L)}$ with different external
momenta, and their total contribution is $C_{Q}^{n}\mathtt{F}_{n}^{(L)}$.

The value of $n$ in Figure 2 satisfies $2\leq n\leq L+1$, so the contribution
of all irreducible cuts is given by

$\mathcal{F}_{Q}^{(L),\text{irreducible}}=\sum_{n=2}^{L+1}C_{Q}^{n}\mathtt{F}_{n}^{(L)}\
.$ (12)

The contribution of reducible cuts can be constructed by the product of
several $\mathtt{F}_{n_{i}}^{(L_{i})}$ divided by the corresponding symmetry
factor. The complete form factor can be written as

$\mathcal{F}_{Q}^{(L)}=\sum_{k=1}^{L}\sum_{\text{allowed
}(n_{i},L_{i})}\frac{Q!}{(Q-\sum_{i=1}^{k}n_{i})!}\frac{1}{S^{L_{1}\cdots
L_{k}}_{n_{1}\cdots
n_{k}}}\prod_{i=1}^{k}\frac{1}{n_{i}!}\mathtt{F}_{n_{i}}^{(L_{i})}\ ,$ (13)

in which the allowed $(n_{i},L_{i})$ should satisfy the following condition:

$n_{i}\geq 2,\ \sum_{i=1}^{k}n_{i}\leq L+k,\ \sum_{i=1}^{k}L_{i}=L,\ $ (14)

and $\\{(n_{i},L_{i})|i=1,\cdots,k\\}$ should be an ordered set, which means
$n_{1}\leq n_{2}\leq\cdots\leq n_{k}$, and if $n_{i}=n_{i+1}$, then $L_{i}\leq
L_{i+1}$. $S^{L_{1}\cdots L_{k}}_{n_{1}\cdots n_{k}}$ is a symmetry factor
when there are duplicated $(n_{i},L_{i})$ pairs. For example,
$S^{31112}_{22335}=2!$ because there are two $(n_{i},L_{i})=(3,1)$ pairs.

As an example of (13), the $L=3$ expression is

$\displaystyle\mathcal{F}_{Q}^{(3)}=C_{Q}^{2}\mathtt{F}_{2}^{(3)}+C_{Q}^{3}\mathtt{F}_{3}^{(3)}+C_{Q}^{4}\mathtt{F}_{4}^{(3)}+6C_{Q}^{4}\mathtt{F}_{2}^{(1)}\mathtt{F}_{2}^{(2)}+10C_{Q}^{5}\mathtt{F}_{2}^{(1)}\mathtt{F}_{3}^{(2)}+15C_{Q}^{6}\Bigl{[}\mathtt{F}_{2}^{(1)}\Bigr{]}^{3}\
.$ (15)

Figure 3: The cuts corresponding to the last term and the first term in (15).

The last term in (15) corresponds to the reducible cut in Figure 3 (a). It has
$(n_{i},L_{i})=(2,1)$ for $i=1,2,3$, so the factor $S_{222}^{111}=3!$. The
overall factor for this cut is

$\frac{Q!}{(Q-6)!}\frac{1}{S^{111}_{222}(2!)^{3}}=\frac{Q!}{(Q-6)!}\frac{1}{3!(2!)^{3}}=15C_{Q}^{6}\
.$ (16)

The first term in (15) corresponds to the irreducible cut in Figure 3 (b), and
the corresponding factor is $\frac{Q!}{(Q-2)!}\frac{1}{2!}=C_{Q}^{2}$.

From (13) and (15) we observe that for a given $L$, $\mathcal{F}_{Q}^{(L)}$
can be constructed by the following building blocks:

$\mathcal{F}_{Q}^{(L)}:\
\\{\mathtt{F}_{n}^{(L_{1})}|n=2,\cdots,L+1;L_{1}=n-1,\cdots,L\\}\ .$ (17)

Therefore the complexity of $\mathcal{F}_{Q}^{(L)}$ only depends on $L$, while
$Q$ only appears in the overall factors.

###  2.3 Integrand of $2n$-point scattering amplitudes

The r.h.s. of the unitarity cut in Figure 2 is a $2n$-point $(L-n+1)$-loop
scattering amplitude $\mathcal{A}_{2n}^{(L-n+1)}$. The list of cut integrands
and scattering amplitudes needed to computed $\mathcal{F}_{Q}^{(L)}$ to 5-loop
are given in Table 1.

$\mathtt{F}_{n}^{(L)}$ | $\mathtt{F}_{2}^{(1)}$ |  $\mathtt{F}_{2}^{(2)}$, $\mathtt{F}_{3}^{(2)}$ |  $\mathtt{F}_{2}^{(3)}$, $\mathtt{F}_{3}^{(3)}$, $\mathtt{F}_{4}^{(3)}$ |  $\mathtt{F}_{2}^{(4)}$, $\mathtt{F}_{3}^{(4)}$, $\mathtt{F}_{4}^{(4)}$, $\mathtt{F}_{5}^{(4)}$ |  $\mathtt{F}_{2}^{(5)}$, $\mathtt{F}_{3}^{(5)}$, $\mathtt{F}_{4}^{(5)}$, $\mathtt{F}_{5}^{(5)}$, $\mathtt{F}_{6}^{(5)}$
---|---|---|---|---|---
$\mathcal{A}_{n}^{(L)}$ | $\mathcal{A}_{4}^{(0)}$ |  $\mathcal{A}_{4}^{(1)}$, $\mathcal{A}_{6}^{(0)}$ |  $\mathcal{A}_{4}^{(2)}$, $\mathcal{A}_{6}^{(1)}$, $\mathcal{A}_{8}^{(0)}$ |  $\mathcal{A}_{4}^{(3)}$, $\mathcal{A}_{6}^{(2)}$, $\mathcal{A}_{8}^{(1)}$, $\mathcal{A}_{10}^{(0)}$ |  $\mathcal{A}_{4}^{(4)}$, $\mathcal{A}_{6}^{(3)}$, $\mathcal{A}_{8}^{(2)}$, $\mathcal{A}_{10}^{(1)}$, $\mathcal{A}_{12}^{(0)}$
Table 1: Cut integrands and scattering amplitudes required to compute
$\mathcal{F}_{Q}^{(L)}$ to 5-loop.

When $n$ or $L$ are large, it can be very difficult to generate the complete
set of Feynman diagrams contributing to $\mathcal{A}_{2n}^{(L-n+1)}$.
Fortunately, as discussed earlier, all Feynman diagrams with the same topology
have the same contribution to the UV counterterm. Therefore we only need to
generate all different topologies, and determine the necessary combinatorial
factors from the permutation of external legs.

Some $\mathcal{A}_{2n}^{(L)}$ topologies are shown in Figure 4, in which we
have added outgoing(incoming) arrows to $\varphi$($\bar{\varphi}$) external
legs. These topologies can be constructed by finding all inequivalent ways of
adding arrows to the external legs in $\phi^{4}$ topologies.

Figure 4: Some $\mathcal{A}_{2n}^{(L)}$ topologies.

The $\phi^{4}$ topologies can be constructed recursively by adding vertices to
or gluing external legs in lower point or lower loop topologies. The
$(2n+2)$-point tree topologies can be generated by finding all inequivalent
ways of adding a $\phi^{4}$ vertex to $2n$-point topologies. An example is
given in Figure 5.

Figure 5: (a): the 6-point $\phi^{4}$ topology. (b): the 8-point $\phi^{4}$
topology. (c) and (d): the topologies contributing to $\mathcal{A}_{8}^{(0)}$.

The $L$-loop $2n$-point $\phi^{4}$ topologies can be generated by gluing two
external legs together in the $(L-1)$-loop $(2n+2)$-point $\phi^{4}$
topologies, and removing graphs containing snail subgraphs. An example is
given in Figure 6.

Figure 6: From 8-point tree topology to 6-point 1-loop topologies. (d) should
be removed because it contains a snail subgraph, therefore it has no
contribution to the UV counterterm.

The integrand of $\mathcal{A}_{2n}^{(L-n+1)}$ can be constructed from the
topologies, together with some combinatorial factors. Let $T$ be a topology of
$\mathcal{A}_{2n}^{(L-n+1)}$, and let $\\{G_{i}(T)|i=1,\cdots,N_{T}\\}$ be the
Feynman diagrams corresponds to $T$. $N_{T}$ depends on the symmetry factor
$S_{\varphi}$($S_{\bar{\varphi}}$) of the diagram under permutation of $n$
$\varphi$($\bar{\varphi}$) external legs:

$N_{T}=\frac{(n!)^{2}}{S_{\varphi}S_{\bar{\varphi}}}\ .$ (18)

The integrand corresponding to $G_{i}(T)$ will be denoted by the same $I(T)$,
then the total integrand can be written as

$\mathcal{A}_{2n}^{(L-n+1)}=\sum_{T}\frac{(n!)^{2}}{S^{int}_{T}S_{\varphi}S_{\bar{\varphi}}}I(T)\
,$ (19)

in which $S^{int}_{T}$ is the symmetry factor of the diagram under permutation
of internal legs.

We can also define $S_{T}\equiv S^{int}_{T}S_{\varphi}S_{\bar{\varphi}}$,
which is the symmetry factor of $T$ under the permutations of internal and
external legs all together. Then (19) becomes

$\mathcal{A}_{2n}^{(L-n+1)}=\sum_{T}\frac{(n!)^{2}}{S_{T}}I(T)\ .$ (20)

The restored cut integrand of the form factor is

$\mathtt{F}_{n}^{(L)}=\frac{1}{l_{1}^{2}\cdots
l_{n}^{2}}\sum_{T}\frac{n!}{S_{T}}I(T)\ ,$ (21)

in which we have added the cut propagators and removed a $n!$ factor to
compensate the symmetry factor of $n$ cut legs.

We also generated some loop form factors directly using QGRAF Nogueira:1991ex
, including all form factors with $Q\leq 6,L\leq 3$, and
$\varphi^{2}\rightarrow 2\varphi$, $\varphi^{3}\rightarrow 3\varphi$ to 5
loop. All the results are consistent with the unitarity cut method.

##  3 The anomalous dimensions

The bare operator is defined as $\varphi^{Q}_{bare}\equiv
Z_{\varphi^{Q}}\varphi^{Q}$, and the renormalization Z-factor can be
determined by the one-particle-irreducible (1PI) part of the form factor,
$\mathcal{F}_{Q}^{(L),1PI}$. Unitarity cut method produces the complete
scattering amplitudes (form factors) including the contributions of both 1PI
and one-particle-reducible(1PR) diagrams, but the contribution of 1PR diagrams
can be easily removed in a pure scalar theory777For $\varphi^{Q}\rightarrow
Q\varphi$ form factor, 1PR diagrams with
$\frac{1}{(p_{1}+p_{2}+\cdots+p_{n})^{2}}$ type poles are forbidden by charge
conservation unless $n=1$. So we only need to remove the diagrams with
$\frac{1}{p_{i}^{2}}$ type propagators.

We evaluate the UV divergences of the form factors using the UV decomposition
method. The UV divergences of multiloop integrals can be decomposed into the
local divergence and lower loop sub divergences, and the anomalous dimensions
are determined solely by the local divergences. Here we only demonstrate the
method by considering a two loop integral which corresponds to the Feynman
diagram in Figure 7. More details can be found in Jin:2022zcj .

Figure 7: An integral which contributes to $\mathtt{F}_{2}^{(2)}$. The arrows
in the diagram represent the directions of the momenta.

The local divergence of a generic integral can be reduced to that of vacuum
integrals by an asymptotic expansion around large loop momenta $l_{i}$. Since
this integral has divergence degree 0, we can simply set $p_{i}\rightarrow 0$
without changing the local divergence:

$\displaystyle\mathbf{L}\frac{1}{(l_{1}+p_{1})^{2}(l_{1}-p_{2})^{2}(l_{1}-l_{2})^{2}l_{2}^{2}}=\mathbf{L}\frac{1}{(l_{1}^{2})^{2}(l_{1}-l_{2})^{2}l_{2}^{2}}\
,$ (22)

in which $\mathbf{L}$ is the local divergence operator, and we have ignored
the integration symbol. Then the infrared divergence of the vacuum integral is
regulated by adding mass to a propagator:

$\displaystyle\mathbf{L}\frac{1}{(l_{1}^{2})^{2}(l_{1}-l_{2})^{2}l_{2}^{2}}=\mathbf{L}\frac{1}{(l_{1}^{2}+m^{2})^{2}(l_{1}-l_{2})^{2}l_{2}^{2}}\
.$ (23)

The local divergence of the massive vacuum integral can be obtained by
subtracting sub UV divergences from the total UV divergence of the
integral888The local divergence is different from Jin:2022zcj  because the
different conventions of $\epsilon$. In this work $D=4-\epsilon$, but in
Jin:2022zcj  $D=4-2\epsilon$. The $e^{\epsilon\gamma}$ factor in (24) is from
the convention of $\overline{\text{MS}}$ scheme.:

$\displaystyle\mathbf{L}\frac{1}{(l_{1}^{2}+m^{2})^{2}(l_{1}-l_{2})^{2}l_{2}^{2}}=\frac{1}{(l_{1}^{2}+m^{2})^{2}(l_{1}-l_{2})^{2}l_{2}^{2}}-\frac{1}{(l_{1}^{2}+m^{2})^{2}}\mathbf{L}\frac{1}{(l_{1}-l_{2})^{2}l_{2}^{2}}$
(24) $\displaystyle=$
$\displaystyle\frac{m^{2d-8}e^{\epsilon\gamma}\Gamma(4-d)\Gamma(2-\frac{d}{2})\Gamma(-1+\frac{d}{2})^{2}}{\Gamma(\frac{d}{2})}-\frac{e^{\frac{1}{2}\epsilon\gamma}}{\epsilon}m^{d-4}\Gamma(2-\frac{d}{2})=-\frac{2}{\epsilon^{2}}+\frac{1}{\epsilon}\
.$

The Z-factor is given by the local divergence of the form factor:

$Z_{\varphi^{Q}}=1-\sum_{L=1}^{\infty}\mathbf{L}\mathcal{F}_{Q}^{(L),1PI}\ .$
(25)

The UV decomposition method turns out to be quite efficient even at 5-loop,
and we were able to evaluate the local divergences of hundreds of 5-loop
integrals in a few hours on a laptop computer.

The anomalous dimensions $\gamma_{\varphi^{Q}}$ can be derived from the
combination $\mathcal{Z}_{Q}\equiv Z_{\phi}^{-\frac{Q}{2}}Z_{\varphi^{Q}}$:

$\displaystyle\gamma_{\varphi^{Q}}=$
$\displaystyle-\frac{\partial\ln\mathcal{Z}_{Q}}{\partial\ln\mu}=\Bigl{[}\epsilon
g-\beta(g)\Bigr{]}\mathcal{Z}_{Q}^{-1}\frac{\partial\mathcal{Z}_{Q}}{\partial
g}=Q\sum_{L=1}^{\infty}g^{L}\gamma_{L}\ .$ (26)

The wavefunction renormalization factor $Z_{\phi}$ and the beta function
$\beta(g)$ can be found in Kleinert:1991rg , and we have also reproduced them
using the UV decomposition method.

The 1PI form factors vanishes when $Q=1$, but (26) still holds:

$\gamma_{\varphi^{1}}=-\frac{\partial\ln
Z_{\phi}^{-\frac{1}{2}}}{\partial\ln\mu}=\frac{1}{2}\frac{\partial\ln
Z_{\phi}}{\partial\ln\mu}\ .$ (27)

The $\gamma_{L}$ in (26) are given by

$\displaystyle\gamma_{1}=$ $\displaystyle\frac{Q-1}{3},\ $ (28)
$\displaystyle\gamma_{2}=$
$\displaystyle\frac{2Q^{2}}{9}+\frac{(N-6)Q}{18}+\frac{2-3N}{36},\ $
$\displaystyle\gamma_{3}=$
$\displaystyle\frac{8Q^{3}}{27}+Q^{2}\Bigl{[}\frac{N-36}{27}+\frac{2(N+14)\zeta_{3}}{27}\Bigr{]}+Q\Bigl{[}\frac{646+23N-2N^{2}}{216}-\frac{2(N+14)\zeta_{3}}{9}\Bigr{]}$
$\displaystyle+\frac{3N^{2}-72N-860}{432}+\frac{4(N+14)\zeta_{3}}{27},\ $
$\displaystyle\gamma_{4}=$
$\displaystyle\frac{14Q^{4}}{27}+Q^{3}\Bigl{[}\frac{4(N-73)}{81}+\frac{2(6N+65)\zeta_{3}}{81}+\frac{5(N+30)\zeta_{5}}{81}\Bigr{]}$
(29)
$\displaystyle+Q^{2}\Bigl{[}\frac{(1634-15N-2N^{2})}{162}+\frac{(3N^{2}-14N-360)\zeta_{3}}{81}-\frac{(N+8)(N+14)\pi^{4}}{4860}$
$\displaystyle\ \ \ \ \ \ \ \ \ +\frac{70(N-2)\zeta_{5}}{81}\Bigr{]}$
$\displaystyle+Q\Bigl{[}\frac{3N^{3}+137N^{2}+2448N-16748}{1944}-\frac{(N^{3}+36N^{2}+168N-816)\zeta_{3}}{324}$
$\displaystyle\ \ \ \ \ \ \ \
+\frac{(3N^{2}+76N+380)\pi^{4}}{4860}-\frac{35(7N+18)\zeta_{5}}{81}\Bigr{]}$
$\displaystyle-\frac{9N^{3}+1144N^{2}+20952N-22120}{15552}+\frac{(N^{3}+24N^{2}+176N+104)\zeta_{3}}{324}$
$\displaystyle-\frac{(N^{2}+27N+134)\pi^{4}}{2430}+\frac{10(17N+62)\zeta_{5}}{81}\
.$

The expression of $\gamma_{5}$ can be found in Appendix C.

The Wilson-Fisher fixed point is define by

$\frac{\partial g}{\partial\ln\mu}\Bigr{|}_{g=g*}=\epsilon
g^{*}-\beta(g^{*})=0\ .$ (30)

At the Wilson-Fisher fixed point, the scaling dimension of $\varphi^{Q}$ can
be expanded into series of $\epsilon$,

$\triangle_{\varphi^{Q}}=Q(1-\frac{\epsilon}{2})+\gamma_{{\varphi^{Q}}}\Bigr{|}_{g=g^{*}}\equiv
Q+Q\sum_{l=1}^{\infty}\epsilon^{l}P_{l}(Q)\ ,$ (31)

in which $Q(1-\frac{\epsilon}{2})$ is the classical scaling dimension of
$\varphi^{Q}$ in the $4-\epsilon$ dimensional space.

Our result for $\triangle_{\varphi^{Q}}$ is in complete agreement with
Antipin:2020abu , in which $P_{1}(Q)$ to $P_{4}(Q)$ are fixed using both
semiclassical and known perturbative results at $Q=1,2,4$. The 5-loop
contribution $P_{5}(Q)\equiv\sum_{i=0}^{5}P_{5}^{i}Q^{i}$ is given by

$\displaystyle P_{5}^{5}=$ $\displaystyle\frac{256}{(8+N)^{5}},$ (32)
$\displaystyle P_{5}^{4}=$
$\displaystyle\frac{3(7N^{2}-1248N-8752)}{(8+N)^{6}}+\frac{28(28+3N)\zeta_{3}}{(8+N)^{5}}+\frac{40(22+N)\zeta_{5}}{(8+N)^{5}}+\frac{14(62+N)\zeta_{7}}{(8+N)^{5}},$
$\displaystyle P_{5}^{3}=$
$\displaystyle-\frac{6(-173104-48736N-3036N^{2}+73N^{3}+N^{4})}{(8+N)^{7}}+\frac{7(206+361N)\zeta_{7}}{4(8+N)^{5}}$
$\displaystyle+\frac{(18N^{3}-673N^{2}-15210N-60128)\zeta_{3}}{(8+N)^{6}}-\frac{(65+6N)\pi^{4}}{60(8+N)^{4}}-\frac{5(30+N)\pi^{6}}{756(8+N)^{4}}$
$\displaystyle+\frac{(1104+158N+13N^{2})\zeta_{3}^{2}}{2(8+N)^{5}}+\frac{5(-16736-1744N+166N^{2}+5N^{3})\zeta_{5}}{2(8+N)^{6}},$
$\displaystyle P_{5}^{2}=$
$\displaystyle\frac{-19621472-8157920N-1037502N^{2}-19068N^{3}+3419N^{4}+\frac{291N^{5}}{4}+\frac{N^{6}}{2}}{(8+N)^{8}}$
$\displaystyle-\frac{(-3372512-1321456N-134112N^{2}+1116N^{3}+362N^{4}+N^{5})\zeta_{3}}{2(8+N)^{7}}$
$\displaystyle-\frac{(-7464-1720N-44N^{2}+3N^{3})\pi^{4}}{120(8+N)^{5}}+\frac{49(-76+4N+9N^{2})\zeta_{7}}{4(8+N)^{5}}$
$\displaystyle-\frac{(37472+12328N+1258N^{2}+53N^{3})\zeta_{3}^{2}}{(8+N)^{6}}-\frac{5(-2+N)\pi^{6}}{54(8+N)^{4}}$
$\displaystyle+\frac{(86496-26060N-4324N^{2}+301N^{3}+2N^{4})\zeta_{5}}{2(8+N)^{6}},$

$\displaystyle P_{5}^{1}=$
$\displaystyle-\frac{1}{64(8+N)^{9}}\Bigl{(}-11161819136-6069470208N-1140742656N^{2}$
(33)
$\displaystyle-61478976N^{3}+5999960N^{4}+824504N^{5}+25772N^{6}+463N^{7}+4N^{8}\Bigr{)}$
$\displaystyle+\frac{(-514112-181184N-12224N^{2}+904N^{3}+52N^{4}+N^{5})\pi^{4}}{480(8+N)^{6}}$
$\displaystyle+\frac{5(2496+958N+65N^{2})\pi^{6}}{756(8+N)^{5}}$
$\displaystyle-\frac{\zeta_{3}}{8(8+N)^{8}}\Bigl{(}159301888+82561536N+12675296N^{2}+111872N^{3}$
$\displaystyle-94960N^{4}-3608N^{5}-48N^{6}+N^{7}\Bigr{)}$
$\displaystyle+\frac{(1491072+729440N+116880N^{2}+7554N^{3}+247N^{4})\zeta_{3}^{2}}{2(8+N)^{7}}$
$\displaystyle-\frac{(-1323968-1009544N-77280N^{2}+22856N^{3}+2021N^{4}+6N^{5})\zeta_{5}}{2(8+N)^{7}}$
$\displaystyle-\frac{7(190840+76762N+7771N^{2}+189N^{3})\zeta_{7}}{4(8+N)^{6}}\
,$

$\displaystyle P_{5}^{0}=$
$\displaystyle\frac{1}{256(8+N)^{9}}\Bigl{(}-18215837696-9842300928N-1755406336N^{2}$
(34)
$\displaystyle-53369344N^{3}+18477984N^{4}+1992280N^{5}+47688N^{6}+648N^{7}+3N^{8}\Bigr{)}$
$\displaystyle-\frac{(-300096-100480N-3496N^{2}+1056N^{3}+40N^{4}+N^{5})\pi^{4}}{480(8+N)^{6}}$
$\displaystyle+\frac{\zeta_{3}}{16(8+N)^{8}}\Bigl{(}155171328+78629632N+10305472N^{2}-485824N^{3}$
$\displaystyle-141784N^{4}-3208N^{5}-68N^{6}+3N^{7}\Bigr{)}-\frac{5(1240+418N+25N^{2})\pi^{6}}{378(8+N)^{5}}$
$\displaystyle-\frac{(481088+242512N+38280N^{2}+2278N^{3}+77N^{4})\zeta_{3}^{2}}{(8+N)^{7}}$
$\displaystyle+\frac{(-1340096-657800N+13132N^{2}+26662N^{3}+1759N^{4}+4N^{5})\zeta_{5}}{2(8+N)^{7}}$
$\displaystyle+\frac{21(31580+12236N+1145N^{2}+21N^{3})\zeta_{7}}{2(8+N)^{6}}\
.$

The $P_{55}$ and $P_{54}$ in (32) are also in agreement with the semiclassical
computation in Antipin:2020abu . The anomalous dimension and scaling dimension
at the Wilson-Fisher fixed point are also available in ancillary files to the
arXiv submission of this article.

In Badel:2019oxl , the scaling dimension in the $\phi^{4}$ theory with $U(1)$
symmetry was computed to 5-loop by combining semiclassical and perturbative
results. This scaling dimension can also be obtained by setting $N=2$ in our
result. We find that $P_{1}$ to $P_{3}$ in Badel:2019oxl  are consistent with
our result. $P_{4}$ is different, which stems from a typo in Calabrese:2002bm
, as was pointed out in Antipin:2020abu . $P_{5}$ was given in numerical form
in Badel:2019oxl , and it is in agreement with our result.

Recently, the complete $Q$ dependence of $P_{5}(Q)$ was also obtained in
Henriksson:2022rnm  by combining the results of semi-classical and
perturbative computations999The $Q=3$ results was extracted from the general
results in Bednyakov:2021ojn ., and it was in complete agreement with our
result.

##  4 Discussion

In this work we have computed the anomalous dimensions of the $\phi^{Q}$-type
operators to 5-loop. It is worth extending the computation to 6-loop, which
gives further check to the scaling dimensions computed using the semiclassical
method, and also allows us to test the efficiency of the UV decomposition
method beyond 5-loop. It should not be hard to generate the 6-loop form factor
integrand using the unitarity cut method, and the analytic expressions of
propagator type master integrals needed for the computation of 6-loop UV
divergence can be found in Georgoudis:2021onj .

The large-$Q$ scaling dimensions in a $\phi^{4}$ theory with $U(N)\times U(N)$
symmetry were computed using semiclassical methods in Antipin:2020rdw , and
the results were verified using direct Feynman diagram computation to 4-loop
in Jack:2021lja . With the help of unitarity cut and UV decomposition method,
one should be able to obtain the full $Q$-dependence of scaling dimensions to
5-loop.

The same methods may also be applied to other operators, like the
$\partial^{2m}\phi^{n}$-type operators considered in Cao:2021cdt , and the
operators in more complicated $O(N)$ representations considered in
Kehrein:1995ia .

###### Acknowledgements.

We would like to thank Bo Feng, Johan Henriksson, Hui Luo, Rui Yu, Gang Yang
for helpful discussions.

##  Appendix A The factorization of $O(N)$ tensor structure

In each loop Feynman diagram contributing to the $\phi^{Q}\rightarrow Q\phi$
form factor, the tree form factor can be factored out, and the integrand is
free of $O(N)$ indices. To understand this feature, let us rewrite the
Lagrangian with the help of an auxiliary field $\sigma$,

$\mathcal{L}\rightarrow\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{1}{2}\sigma^{2}+\sqrt{\frac{g}{12}}\sigma\phi^{2}\
.$ (35)

The original Lagrangian is recovered when the auxiliary field is integrated
out.

In the new Lagrangian, $\sigma$ is a $O(N)$ singlet, and the interaction
vertices are shown in Figure 8. In the Figure, dashed lines represents
$\sigma$, solid lines without arrows represents $\phi^{i}$, and solid lines
with arrows represents $\varphi\equiv q\cdot\phi$ and $\bar{\varphi}$.

Figure 8: The interaction vertices of the Lagrangian given by (35).

An example of Feynman diagrams contributing to $\phi^{Q}\rightarrow Q\phi$ is
shown in Figure 9(a). In these diagrams, $\phi$-fields either form a closed
$\phi$-loop, or a $\phi$-line connecting two external legs. The $\phi$-lines
attached to the operator will not form $\phi$-loops, because the operator is
traceless.

Figure 9: Feynman diagram contributing to (a) $\phi^{Q}\rightarrow Q\phi$ and
(b) $\varphi^{Q}\rightarrow Q\varphi$.

The $O(N)$ tensor structure of the diagram can be easily read from the
diagram. Each $\phi$-loop contributes a $N$ factor, and each $\phi$-line
contributes a $\delta_{ij}$ factor. So it is clear that the tensor structure
is proportional to the Feynman rule of the operator.

An example of Feynman diagrams contributing to $\varphi^{Q}\rightarrow
Q\varphi$ is shown in Figure 9(b). These diagrams only contain $\varphi$-lines
and $\phi$-loops.

##  Appendix B The superficial degree of divergence of $\phi^{Q}\rightarrow
Q\phi$ form factor

An important property of $\mathcal{F}_{Q}^{(L)}$ is that the superficial
degree of UV divergence is 0. To see this, let $V$ and $E$ be the number of
vertices and propagators of a Feynman diagram contributing to
$\mathcal{F}_{Q}^{(L)}$. The diagram has $(V-1)$ 4-vertices, a single
$(Q+1)$-vertex, and $Q+1$ external legs, so the number of propagators of the
diagram is

$E=\frac{4(V-1)+(Q+1)-(Q+1)}{2}=2(V-1)\ .$ (36)

A connected Feynman diagram also satisfies

$V-E+L=1\ .$ (37)

Combining (36) and (37), we derive $E=2L$. $E$ equals the number of loop
propagators in the integral, therefore the superficial degree of divergence in
4-D is given by

$\omega=4L-2E=0\ .$ (38)

The renormalization Z-factors only depend on the local UV divergences of
integrals, which have polynomial dependence on external momenta (see e.g.
Jin:2022zcj ). When $\omega=0$, the local UV divergence does not depend on the
external momenta:

$\mathbf{L}\mathcal{F}_{Q}^{(L)}(p_{1},\cdots,p_{Q})=\mathbf{L}\mathcal{F}_{Q}^{(L)}(0,\cdots,0)\
.$ (39)

##  Appendix C The 5-loop corrections to the anomalous dimension

The expression of $\gamma_{5}$ in (26) is given by:

$\displaystyle\gamma_{5}\equiv$ $\displaystyle\sum_{i=0}^{5}\gamma_{5i}Q^{i}\
,\ \ \ \ \gamma_{55}=\frac{256}{243},$ (40) $\displaystyle\gamma_{54}=$
$\displaystyle\frac{7N-800}{81}+\frac{28(3N+28)\zeta_{3}}{243}+\frac{40(N+22)\zeta_{5}}{243}+\frac{14(N+62)\zeta_{7}}{243}\
,$ $\displaystyle\gamma_{53}=$
$\displaystyle\frac{3074-33N-2N^{2}}{81}+\frac{(18N^{2}-385N-3994)\zeta_{3}}{243}-\frac{(N+8)(6N+65)\pi^{4}}{14580}$
$\displaystyle-\frac{5(N+8)(N+30)\pi^{6}}{183708}+\frac{(13N^{2}+158N+1104)\zeta_{3}^{2}}{486}$
$\displaystyle+\frac{5(5N^{2}+198N-832)\zeta_{5}}{486}+\frac{7(361N+206)\zeta_{7}}{972}\
,$ $\displaystyle\gamma_{52}=$
$\displaystyle\frac{12N^{3}+881N^{2}+10230N-407000}{5832}-\frac{(3N^{3}+28N^{2}-216N-2056)\pi^{4}}{29160}$
$\displaystyle+\frac{(41348+5776N-161N^{2}-2N^{3})\zeta_{3}}{972}+\frac{(12372-924N+285N^{2}+2N^{3})\zeta_{5}}{486}$
$\displaystyle-\frac{5(N-2)(N+8)\pi^{6}}{13122}-\frac{(53N^{2}+474N+1912)\zeta_{3}^{2}}{243}$
$\displaystyle+\frac{49(-76+4N+9N^{2})\zeta_{7}}{972}\ ,$
$\displaystyle\gamma_{51}=$
$\displaystyle\frac{14368320+669056N+37734N^{2}-2981N^{3}-48N^{4}}{186624}$
$\displaystyle+\frac{(-205528-37536N+868N^{2}-9N^{3}-2N^{4})\zeta_{3}}{3888}$
$\displaystyle+\frac{(10136+5960N+878N^{2}+44N^{3}+N^{4})\pi^{4}}{116640}$
$\displaystyle+\frac{5(1380+628N+53N^{2})\pi^{6}}{183708}-\frac{(34252+5740N+1155N^{2}+6N^{3})\zeta_{5}}{486}$
$\displaystyle+\frac{(2580+1376N+211N^{2})\zeta_{3}^{2}}{486}-\frac{7(7286+2731N+189N^{2})\zeta_{7}}{972}\
,$ $\displaystyle\gamma_{50}=$
$\displaystyle\frac{-6934336-1059072N-84662N^{2}+2223N^{3}+9N^{4}}{186624}$
(41)
$\displaystyle+\frac{(91864+19560N-460N^{2}+13N^{3}+3N^{4})\zeta_{3}}{3888}-\frac{5(682+253N+19N^{2})\pi^{6}}{91854}$
$\displaystyle+\frac{(-14552-6176N-758N^{2}-32N^{3}-N^{4})\pi^{4}}{116640}+\frac{(70-293N-59N^{2})\zeta_{3}^{2}}{243}$
$\displaystyle+\frac{(24280+5594N+845N^{2}+4N^{3})\zeta_{5}}{486}+\frac{7(1186+389N+21N^{2})\zeta_{7}}{162}\
.$

## References

  * (1) G. Badel, G. Cuomo, A. Monin and R. Rattazzi, _The Epsilon Expansion Meets Semiclassics_ , _JHEP_ 11 (2019) 110 [1909.01269].
  * (2) O. Antipin, J. Bersini, F. Sannino, Z.-W. Wang and C. Zhang, _Charging the $O(N)$ model_, _Phys. Rev. D_ 102 (2020) 045011 [2003.13121].
  * (3) K.G. Wilson, _Renormalization group and critical phenomena. 1. Renormalization group and the Kadanoff scaling picture_ , _Phys. Rev. B_ 4 (1971) 3174.
  * (4) K.G. Wilson, _Renormalization group and critical phenomena. 2. Phase space cell analysis of critical behavior_ , _Phys. Rev. B_ 4 (1971) 3184.
  * (5) K.G. Wilson and M.E. Fisher, _Critical exponents in 3.99 dimensions_ , _Phys. Rev. Lett._ 28 (1972) 240.
  * (6) J.C. Le Guillou and J. Zinn-Justin, _Critical Exponents from Field Theory_ , _Phys. Rev. B_ 21 (1980) 3976.
  * (7) D.T. Son, _Semiclassical approach for multiparticle production in scalar theories_ , _Nucl. Phys. B_ 477 (1996) 378 [hep-ph/9505338].
  * (8) S. Hellerman, D. Orlando, S. Reffert and M. Watanabe, _On the CFT Operator Spectrum at Large Global Charge_ , _JHEP_ 12 (2015) 071 [1505.01537].
  * (9) L. Alvarez-Gaume, O. Loukas, D. Orlando and S. Reffert, _Compensating strong coupling with large charge_ , _JHEP_ 04 (2017) 059 [1610.04495].
  * (10) L. Alvarez-Gaume, D. Orlando and S. Reffert, _Large charge at large N_ , _JHEP_ 12 (2019) 142 [1909.02571].
  * (11) G. Arias-Tamargo, D. Rodriguez-Gomez and J.G. Russo, _The large charge limit of scalar field theories and the Wilson-Fisher fixed point at $\epsilon=0$_, _JHEP_ 10 (2019) 201 [1908.11347].
  * (12) L.A. Gaumé, D. Orlando and S. Reffert, _Selected topics in the large quantum number expansion_ , _Phys. Rept._ 933 (2021) 2180 [2008.03308].
  * (13) I. Jack and D.R.T. Jones, _Anomalous dimensions at large charge in d=4 O(N) theory_ , _Phys. Rev. D_ 103 (2021) 085013 [2101.09820].
  * (14) O. Antipin, J. Bersini, F. Sannino, Z.-W. Wang and C. Zhang, _Charging non-Abelian Higgs theories_ , _Phys. Rev. D_ 102 (2020) 125033 [2006.10078].
  * (15) I. Jack and D.R.T. Jones, _Anomalous dimensions at large charge for U(N) $\times$U(N) theory in three and four dimensions_, _Phys. Rev. D_ 104 (2021) 105017 [2108.11161].
  * (16) O. Antipin, J. Bersini, F. Sannino, Z.-W. Wang and C. Zhang, _Untangling scaling dimensions of fixed charge operators in Higgs theories_ , _Phys. Rev. D_ 103 (2021) 125024 [2102.04390].
  * (17) Z. Bern, L.J. Dixon, D.C. Dunbar and D.A. Kosower, _One loop n point gauge theory amplitudes, unitarity and collinear limits_ , _Nucl.Phys._ B425 (1994) 217 [hep-ph/9403226].
  * (18) Z. Bern, J. Carrasco and H. Johansson, _New Relations for Gauge-Theory Amplitudes_ , _Phys.Rev._ D78 (2008) 085011 [0805.3993].
  * (19) R. Boels and C. Schwinn, _Deriving CSW rules for massive scalar legs and pure Yang-Mills loops_ , _JHEP_ 0807 (2008) 007 [0805.1197].
  * (20) Z. Bern, J. Carrasco, H. Johansson and R. Roiban, _The Five-Loop Four-Point Amplitude of N=4 super-Yang-Mills Theory_ , _Phys.Rev.Lett._ 109 (2012) 241602 [1207.6666].
  * (21) F. Cachazo, S. He and E.Y. Yuan, _Scattering of Massless Particles in Arbitrary Dimensions_ , _Phys.Rev.Lett._ 113 (2014) 171601 [1307.2199].
  * (22) F. Cachazo, S. He and E.Y. Yuan, _Scattering of Massless Particles: Scalars, Gluons and Gravitons_ , 1309.0885.
  * (23) Z. Bern, S. Davies and J. Nohle, _Double-Copy Constructions and Unitarity Cuts_ , _Phys. Rev._ D93 (2016) 105015 [1510.03448].
  * (24) G. Yang, _Color-kinematics duality and Sudakov form factor at five loops for N=4 supersymmetric Yang-Mills theory_ , _Phys. Rev. Lett._ 117 (2016) 271602 [1610.02394].
  * (25) Q. Jin, _The Decomposition of UV Divergence_ , 2203.11554.
  * (26) V.M. Braun and A.N. Manashov, _Evolution equations beyond one loop from conformal symmetry_ , _Eur. Phys. J. C_ 73 (2013) 2544 [1306.5644].
  * (27) M. Kompaniets and K.J. Wiese, _Fractal dimension of critical curves in the $O(n)$-symmetric $\phi^{4}$ model and crossover exponent at 6-loop order: Loop-erased random walks, self-avoiding walks, Ising, XY , and Heisenberg models_, _Phys. Rev. E_ 101 (2020) 012104 [1908.07502].
  * (28) P. Calabrese, A. Pelissetto and E. Vicari, _Multicritical phenomena in O(n(1)) + O(n(2)) symmetric theories_ , _Phys. Rev. B_ 67 (2003) 054505 [cond-mat/0209580].
  * (29) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279.
  * (30) H. Kleinert, J. Neu, V. Schulte-Frohlinde, K.G. Chetyrkin and S.A. Larin, _Five loop renormalization group functions of O(n) symmetric phi**4 theory and epsilon expansions of critical exponents up to epsilon**5_ , _Phys. Lett. B_ 272 (1991) 39 [hep-th/9503230].
  * (31) J. Henriksson, _The critical O(N) CFT: Methods and conformal data_ , 2201.09520.
  * (32) A. Bednyakov and A. Pikelner, _Six-loop beta functions in general scalar theory_ , _JHEP_ 04 (2021) 233 [2102.12832].
  * (33) A. Georgoudis, V. Gonçalves, E. Panzer, R. Pereira, A.V. Smirnov and V.A. Smirnov, _Glue-and-cut at five loops_ , _JHEP_ 09 (2021) 098 [2104.08272].
  * (34) W. Cao, F. Herzog, T. Melia and J.R. Nepveu, _Renormalization and non-renormalization of scalar EFTs at higher orders_ , _JHEP_ 09 (2021) 014 [2105.12742].
  * (35) S.K. Kehrein, _The Spectrum of critical exponents in phi**2 in two-dimensions theory in D = (4-epsilon)-dimensions: Resolution of degeneracies and hierarchical structures_ , _Nucl. Phys. B_ 453 (1995) 777 [hep-th/9507044].


# New links between PDE’s and Voronoi patterns

Yuriria Cortés-Poza, David Padilla-Garza, Pablo Padilla

###### Abstract

We present several results in PDE’s, both analytical and numerical, in which
Voronoi patterns arise:

  * a)

The solution to an elliptic equation and its probabilistic interpretation and
its relation to a stochastic colonization game.

  * b)

Parabolic equations in Riemannian manifolds, with applications to chemical
reactions and diffusive fronts.

  * c)

The result of wavefronts interactions from point light sources and Huygens
principle, as the singular set of weak solutions to the eikonal equation.

  * d)

A transport problem and the corresponding Monge-Ampère equation.

For these problems, the connection with Voronoi cells is presented, and in the
case of the colonization game, we also provide new agent-based numerical
simulations.

##  1 Introduction

Voronoi tessellations appear in many contexts and have important applications
in several fields, including archaeology, biology, computer science,
economics, engineering, epidemiology, medicine, physics, image processing, and
urban planning [1]. They go back as far as Descartes in 1644. He informally
drew a diagram of what is now called a Voronoi tessellation when explaining
his vortex theory about the universe’s structure. Peter Gustav Lejeune
Dirichlet used two-dimensional and three-dimensional Voronoi diagrams to study
quadratic forms in 1850. In Britain, the physician John Snow used a Voronoi
diagram in 1854 to establish that most people who died in the Broad Street
cholera outbreak lived closer to the infected Broad Street water pump than to
any other. Voronoi diagrams are named after Georgy Feodosievych Voronoi, who
defined and studied the general n-dimensional case in 1908 (see [1] for more
applications and details on the history of the subject).

Voronoi tessellations are defined in a straightforward way using the distance
function as follows. Consider a finite set of points in a metric space. Two
arbitrary points in the space are in the same cell if they have a common
closest point to the original set of given points. More precisely, we have the
following definition:

######  Definition 1.1.

Given $N$ points, $P=\\{x_{1},...x_{N}\\}$ on $M^{n}$ an $n$-dimensional
Riemannian manifold (in particular $\mathbf{R}^{n}$), the $i-th$ Voronoi cell,
$V_{i}$ is the set of points in the manifold for which

$x_{i}=\text{argmin }d(x,P).$ (1)

In other words,

$V_{i}=\\{x\in M^{n}|d(x,x_{i})<d(x,x_{j})\text{ for }i\neq j\\}.$ (2)

In [16] and references therein, an up-to-date description of the methods
available in computational geometry can be found for actually obtaining the
Voronoi cells. In [8] a detailed analytical and computational study is made of
generalized Voronoi regions. In this paper, we formalize a simple dynamical
principle that allows characterizing the boundary of Voronoi cells. Consider
$n$ point sources located at the points $x_{i}$ and suppose that a spherical
front starts from each of them at time $t=0$. We consider the set of points
where two fronts coming from different sources meet as time progresses (see
figure 1), i.e., the ”collision” set of the fronts. The case of two points is
easily pictured and determines the perpendicular bisector of the segment
joining $x_{1}$ and $x_{2}$. This set contains the boundary of the Voronoi
tessellation since, in order for two fronts to collide at a certain point,
they have to be equidistant from at least two of the originally given points.
From this set, we can obtain the boundary of the tessellation by removing
those points in the collision set that belong to the interior of a Voronoi
cell.

Figure 1: Two fronts coming from different sources (blue and green) that meet
as time progresses.

This principle is obvious but has several advantages. First, it can be used to
characterize the boundary of the Voronoi cells as the singular set of a weak
solution to a Hamilton-Jacobi equation, the eikonal equation, and then use the
results from [10] to conclude that it is rectifiable, even for the general
problem in higher dimensions. Second, it has a natural and straightforward
numerical implementation. Third, it can be interpreted in terms of diffusive
waves, and this provides a natural interpretation for some biological and
physical phenomena. Recently this has been verified experimentally (see [4]).
This last approach in terms of diffusion also allows to study this problem
from a stochastic process point of view and formulate a dynamic representation
of Voronoi diagrams using random walks. We present in section 2 a colonization
game and its relation to the solution of an elliptic equation based on the
probabilistic representation just mentioned. We also consider another related
problem, namely a harmonic function on a perforated domain, and use the
solution to define what we call a harmonic tessellation. In section 3, we
discuss the formulation in terms of diffusive waves motivated by some
applications to pattern formation in biological, chemical, and physical
phenomena and present some results based on short time estimates of the heat
kernel. We include a short discussion on nonlinear parabolic equations, e.g.,
the KPP equation or the porous medium equation exhibiting spherical fronts as
natural candidates to be studied analytically and numerically. We devote
section 4 to the characterization of Voronoi cells in terms of spherical waves
using Huygen’s principle. Finally, in section 5, we study a simple
transportation problem that gives rise to a Voronoi tessellation and its
connection to the Monge-Ampère equation. The last section is devoted to open
questions and concluding remarks.

##  2 Elliptic equations and Voronoi cells

In the first part of this section, we consider the elliptic problem associated
with the stationary concentration of a chemical diffusing on the plane
resulting from $n$ point sources of equal intensity. We discuss its relation
with Voronoi cells motivated by a colonization game, show numerical
simulations, and describe how it can model different systems, including
bacterial growth and urban dynamics. The second part introduces a tessellation
induced by the graph of a harmonic function on the perforated plane. In order
to fix ideas, we start by considering the stationary concentration of a
substance diffusing on the plane from $n$ sources. More precisely, let
$x_{1},...,x_{n}$ be as above and consider the elliptic problem

$\begin{split}\Delta u&=\sum_{i=1}^{n}\delta(x-x_{i})\quad{\rm
on}\quad\mathbb{R}^{2}\\\ u&\to 0\quad{\rm as}\quad x\to\infty.\end{split}$
(3)

If $\varphi(x)$ is the fundamental solution for Laplace’s equation, then the
solution to the above problem can be written as

$u(x)=\sum_{i=1}^{n}\varphi(x-x_{i}).$ (4)

This solution has a natural probabilistic interpretation: it is (up to a
multiplicative normalization constant) the stationary probability distribution
of many Brownian walks starting at each source. If each of the Brownian walks
is labeled with a different color, according to the source it started from, it
is natural to expect that when the number of Brownian paths tends to infinity,
on average, the solution will be colored as a Voronoi tessellation. This
phenomenon suggests the following colonization game: From each point, $x_{i}$
an equal number of Brownian explorers start each one claiming the territory it
moves around until it finds an explorer of a different color in its
$\varepsilon$-neighborhood, at which point they stop and are replaced by two
new explorers coming out from their respective sources.

The following algorithm formalizes the description above:

  * a)

Domain. We define a rectangle ${\bf R}=\\{(x,y):0\leq x\leq R,0\leq y\leq
U\\}$.

  * b)

Sources and particles: We randomly choose $n$ points in ${\bf R}$; these
points will be the sources. $m$ particles will emerge from each source and
start moving in a random walk with a fixed step size.

  * c)

Coalitions: When a particle $P$ from source $S_{i}$ collides with the
particle’s trajectory from a different source or with the boundary, $P$ will
be sent to its source $S_{i}$ and will resume moving from there.

  * d)

Iterations: The particles will move synchronously and will stop after a
certain number of iterations.

The parameters we used are shown in table (1).

Parameter Value
---
Number of sources. | $3,4,5$
Number of particles per source | $100$
Number of iterations | $300-500$
Steps of the particles. |  $\alpha=1/10$ where $\alpha\in(1-,1)$
Domain ${\bf R}$ | ${\bf R}=\\{(x,y):x,y\in[0,2]\\}$

Table 1: Algorithm parameters

.

We find that typically 300 to 500 iterations and 100 particles are sufficient
to recover Voronoi cells.

The following figures (2,3,4) show the output of several simulations of our
program varying number of sources and number of iterations. In all of them, we
draw with colored points the trajectories of the particles. The ones that
emerge from a given cell have the same color. On top of the colored points, we
draw Voronoi cells, computed using the Delaunay triangulation. The black
circle centers correspond to the sources and are drawn so that they can easily
be identified.

Figure 2: Simulation with 4 sources, 100 particles per source, 800 iterations.
Figure 3: Simulation with 5 sources, 100 particles per source, 500 iterations.
Figure 4: Simulation with 6 sources, 100 particles per source, 500 iterations.

As we can see in all figures, Voronoi cells are recovered from the random
walkers following the abovementioned rules. This colonization process takes
place in biological and social systems and chemical reactions, as can be seen
for instance in [4], where a chemical reaction is documented starting with a
diffusion front from $n$ points replicating the formation of Voronoi patterns
on the plane. It has also been observed in colonies of bacteria starting to
grow from different locations in a Petri dish with limited resources that the
end distribution of the population is a Voronoi tessellation. Finally, some
authors have considered urban growth from different commercial centers and
have found that the resulting neighborhoods correspond to Voronoi cells.
Voronoi diagrams have also been used in the context of social planning (for
more details on the subject, see [17], and [12]). Considering models with
nonlinear diffusive fronts in a more realistic scenario would be interesting,
such as the KPP or the porous medium equation.

In what follows, we consider a problem motivated by the previous setting that
will allow us to define a tessellation determined by a harmonic function.
Before doing so, we notice that given a point $x$ on the plane, another way of
determining to what Voronoi cell it belongs to consists in following the
steepest ascent path determined by the distance function to the set
$\bigcup_{i=1}^{n}\\{x_{i}\\}$. For almost all points on the plane, this path
will end at the $x_{i}$ whose cell $x$ belongs to. So let $u$ be a solution to
the problem.

$\begin{split}\Delta u&=0\quad{\rm in}\quad\Omega\\\ u&=1\quad{\rm
on}\quad\partial\Omega\\\
\Omega&=\mathbb{R}^{2}\setminus\bigcup_{i=1}^{n}B_{\varepsilon}(x_{i})\\\
u&\to 0\quad{\rm as}\quad x\to\infty.\end{split}$ (5)

By the maximum principle, the maximum of $u$ is attained at the boundary of
$\Omega$. Thus, if we follow the path of the steepest ascent starting at a
given arbitrary x, it will lead to a point in $B_{\varepsilon}(x_{i})$ for
some $i$, for almost every point on the plane. For $\varepsilon$ sufficiently
small, this procedure provides the same result for a given $x$, determining a
tessellation that, for simplicity, we call a harmonic tessellation.

##  3 Voronoi patterns and the heat equation

In this section, we consider a mathematical model for the experiment presented
in [4]. In this experiment, $n$ different diffusive sources are placed on a
Petri dish simultaneously. The diffusion front can be tracked visually, and it
stops whenever it meets a front coming from another source. Notice that it is
a very similar process as the one described in the previous section. As a
first approximation, we propose a linear heat equation to model this
phenomenon and prove that for sufficiently short times, the effective
concentration associated with a diffusion source is restricted to the
corresponding Voronoi cell (see theorem 3.1). As mentioned in the previous
section, we suggest improving our model, including a nonlinear term, and
consider, for instance, a Fisher-KPP type equation. The main result in this
section is the following theorem.

######  Theorem 3.1.

Let $M^{n}$ be a smooth Riemannian manifold of dimension $n$, and let
$\\{x_{i}\\}_{i=1}^{N}$ be $N$ distinct points in $M^{n},$ let
$\\{V_{i}\\}_{i=1}^{N}$ be the corresponding Voronoi cells.

Consider $u(x,t)$ the solution of the heat equation with initial condition

$u(x,0)=\sum_{i=1}^{N}w_{i}\delta_{x_{i}},$ (6)

where $w_{i}>0.$

Let $S^{n}$ be the $n-1$ unit sphere in $n$ dimensions and let $\theta\in
S^{n}.$ Assume that $V_{i}$ is compact and contained in the complement of the
cut locus of $x_{i}$ (see [13], section 9.2 for the definition of cut locus).
Then for all $\delta,\epsilon>0$ there exists $T>0$ such that for all $t<T,$
and for every $i,$

$u(\text{exp}_{x_{i}}(s\theta),t)$ (7)

is decreasing in $s$ as long as $s>\delta$ and

$d(\text{exp}_{x_{i}}(s\theta),\partial V_{i})>\epsilon.$ (8)

Furthermore, $T$ can be chosen as a continuous and monotone (in each
component) function of $\delta$ and $\epsilon.$

######  Corollary 3.2.

Let $\epsilon,\delta,T$ be as above. Let $\gamma(s)$ be a smooth path joining
$x_{i}$ and $x_{j}$ for $i\neq j.$ Assume that, whenever $\gamma(s)\in V_{i},$
that

$\frac{d}{ds}d(\gamma(s),x_{i})\neq 0.$ (9)

Then any local minimums of $u(\gamma(s),t)$ occur when

$d\left(\gamma(s),\bigcup_{i=1}^{N}\partial V_{i}\right)\leq\epsilon.$ (10)

We need three lemmas to give the proof of the Theorem 3.1.

######  Lemma 3.3.

Let $M^{n}$ be a smooth compact manifold, and let $\\{x_{i}\\}_{i=1}^{N}$ be
$N$ distinct points in $M^{n}.$ Let $\\{V_{i}\\}_{i=1}^{N}$ be the
corresponding Voronoi cells.

Let

$B_{i}^{\epsilon}=\bigcup_{x\in\partial V_{i}}B(x,\epsilon)$ (11)

and let

$V_{i}^{\epsilon}=V_{i}\setminus B_{i}^{\epsilon}.$ (12)

Assume that $V_{i}$ is bounded. Then there exists a strictly positive, non-
increasing continuous function $\varphi(\epsilon)$ such that for any $x\in
V_{i}^{\epsilon}$ and any $j\neq i,$

$d^{2}(x,x_{j})\geq d^{2}(x,x_{i})+\varphi(\epsilon).$ (13)

###### Proof.

Consider the function

$\Phi(x)=\min_{j\neq i}\\{d(x,x_{j})-d(x,x_{i})\\}.$ (14)

Since $V_{i}$ is bounded, $V_{i}^{\epsilon}$ is compact; and therefore $\Phi$
achieves a minimum in $V_{i}^{\epsilon},$ since $\Phi$ is clearly continuous.
By definition of Voronoi cell, we have that

$\Phi(x)\geq 0,$ (15)

and therefore

$\min_{y\in V_{i}^{\epsilon}}\Phi(y)\geq 0.$ (16)

We now claim that

$\min_{y\in V_{i}^{\epsilon}}\Phi(y)>0.$ (17)

To see this, proceed by contradiction, and assume that

$\min_{y\in V_{i}^{\epsilon}}\Phi(y)=0.$ (18)

Then there exists $x\in V_{i}^{\epsilon},$ and $j\neq i$ such that

$d(x,x_{j})=d(x,x_{i}).$ (19)

Since $x_{i}\neq x_{j}$ this implies that there exists $y\in B(x,\epsilon)$
such that

$d(y,x_{j})<d(y,x_{i}),$ (20)

but this is a contradiction because $B(x,\epsilon)\in V_{i}$ and therefore

$d(y,x_{j})\geq d(y,x_{i}).$ (21)

Let

$\alpha(\epsilon)=\min_{y\in V_{i}^{\epsilon}}\Phi(y),$ (22)

then for any $x\in V_{i}^{\epsilon},$

$\begin{split}d^{2}(x,x_{j})&\geq(d(x,x_{i})+\alpha(\epsilon))^{2}\\\
&=d^{2}(x,x_{i})+2\alpha(\epsilon)d(x,x_{i})+\alpha(\epsilon)^{2}\\\
&>d^{2}(x,x_{i})+\alpha(\epsilon)^{2}.\end{split}$ (23)

Then

$\varphi(\epsilon)=\alpha^{2}(\epsilon)$ (24)

satisfies that

$d^{2}(x,x_{j})\geq d^{2}(x,x_{i})+\varphi(\epsilon).$ (25)

Clearly $\varphi$ is non-increasing. We will show that it is continuous by
showing that it is continuous from the right and from the left. We will first
show continuity from the left . Let

$\epsilon_{n}\to\epsilon$ (26)

monotonically from the left. Clearly

$\varphi(\epsilon)\geq\lim_{n\to\infty}\varphi(\epsilon_{n})$ (27)

since

$V_{i}^{\epsilon}\subset V_{i}^{\epsilon_{n}}$ (28)

and therefore

$\varphi(\epsilon)\geq\varphi(\epsilon_{n}).$ (29)

Now let

$x_{n}=\text{argmin}_{y\in V_{i}^{\epsilon_{n}}}\Phi(y).$ (30)

Modulo a subsequence,

$x_{n}\to x\in V_{i}^{\epsilon},$ (31)

and by continuity of $\Phi,$

$\begin{split}\alpha(\epsilon)&\leq\Phi^{2}(x)\\\
&=\lim_{n\to\infty}\Phi^{2}(x_{n})\\\
&=\lim_{n\to\infty}\alpha(\epsilon).\end{split}$ (32)

Hence,

$\varphi(\epsilon)\leq\lim_{n\to\infty}\varphi(\epsilon_{n}),$ (33)

and therefore $\varphi$ is continuous from the left. To show that it is
continuous from the right, let

$\epsilon_{n}\to\epsilon$ (34)

from the right. Then

$V_{i}^{\epsilon_{n}}\subset V_{i}^{\epsilon}$ (35)

and

$\varphi(\epsilon)\leq\lim_{n\to\infty}\varphi(\epsilon_{n}).$ (36)

Let

$x=\text{argmin}_{y\in V_{i}^{\epsilon}}\Phi(y),$ (37)

and let $x_{n}\in V_{i}^{\epsilon_{n}}$ be such that

$\lim_{n\to\infty}x_{n}=x.$ (38)

Then

$\begin{split}\alpha(\epsilon)&=\Phi^{2}(x)\\\ &=\lim\Phi^{2}(x_{n})\\\
&\geq\lim_{n\to\infty}\min_{y\in V_{i}^{\epsilon_{n}}}\Phi^{2}(y)\\\
&=\lim_{n\to\infty}\alpha(\epsilon_{n}).\end{split}$ (39)

Therefore $\varphi$ is continuous from the right and the left and hence
continuous. ∎

Before giving the proof of the next lemma, we need the definition of the
fundamental solution to the heat equation:

######  Proposition 3.4.

([6], Proposition 4.1)

Let $M^{n}$ be an $n-$dimensional smooth manifold. There exists a unique
differentiable function

$p:(0,\infty)\times M^{n}\times M^{n}$ (40)

denoted by $p_{t}(x,y)$ such that

$\frac{\partial}{\partial t}p_{t}(x,y)=\Delta_{y}p_{t}(x,y)$ (41)

and

$\lim_{t\to 0}p_{t}(x,\cdot)=\delta_{x}$ (42)

for every $x,y\in M^{n}$ and $t>0$, where $\Delta$ is the Laplace-Beltrami
operator. Such a function will be called the heat kernel for $\Delta$. It is
non-negative, it is mass preserving, i.e.

$\int_{M^{n}}p_{t}(x,y)d\pi(y)=1$ (43)

for every $x\in M^{n}$ and $t>0$, it is symmetric, i.e.

$p_{t}(x,y)=p_{t}(y,x)$ (44)

for every $x,y\in M^{n}$ and $t>0$ and it satisfies the semi-group property
i.e.

$\int_{M^{n}}p_{t}(x,y)p_{s}(y,z)d\pi(y)=p_{t+s}(x,z)$ (45)

for every $x,y\in M^{n}$ and $t,s>0$.

######  Lemma 3.5.

Let $M$ be a manifold, and let $p_{t}(x,y)$ be the fundamental solution to the
heat equation. Then

$\log|\nabla_{y}p_{t}(x,y)|\leq
C+\left(n-1\right)\log(t)-\frac{d^{2}(x,y)}{4t}+\log\left(d(x,y)\right),$ (46)

where the constant $C$ depends only on $M^{n}.$

###### Proof.

Let

$L_{x}^{t}(y)=t\log(p_{t}(x,y)).$ (47)

Then by chain rule

$\begin{split}\nabla_{y}p_{t}&=\frac{1}{t}\nabla
L_{x}^{t}\exp\left(\frac{1}{t}L_{x}^{t}\right)\\\
&=\frac{1}{t}p_{t}(x,y)\nabla L_{x}^{t}.\end{split}$ (48)

Using lemma 4.7 of [6] we have that, for $t\in(0,1),$

$p_{t}(x,y)\leq\frac{C}{t^{n-\frac{1}{2}}}\exp\left(\frac{d^{2}(x,y)}{4t}\right),$
(49)

where $C$ depends only on $M^{n}.$

On the other hand, theorem 5.5.3 of [7] implies that

$\frac{1}{t}|\nabla L_{x}^{t}|\leq
C\left[\frac{d(x,y)}{t}+\frac{1}{\sqrt{t}}\right],$ (50)

where $C$ depends only on $M^{n}.$ Hence,

$|\nabla L_{x}^{t}|\leq C\left[d(x,y)+\sqrt{t}\right].$ (51)

Putting everything together, we have

$|\nabla p_{t}|\leq
C\left[d(x,y)+\sqrt{t}\right]\frac{1}{t^{n-\frac{1}{2}}}\exp\left(\frac{d^{2}(x,y)}{4t}\right),$
(52)

which implies

$\begin{split}\log|\nabla_{y}p_{t}(x,y)|&\leq
C+\left(n-\frac{1}{2}\right)\log(t)-\frac{d^{2}(x,y)}{4t}-\log(t)+\log\left(d(x,y)\right)+\frac{1}{2}\log(t)\\\
&\leq
C+\left(n-1\right)\log(t)-\frac{d^{2}(x,y)}{4t}+\log\left(d(x,y)\right).\end{split}$
(53)

∎

######  Lemma 3.6.

Let $M^{n}$ be a manifold of dimension $n$, and let $p_{t}(x,y)$ be the
fundamental solution to the heat equation. Let $x\in M^{n}$ and let $K\subset
M^{n}$ be a compact subset of $M^{n}$ containing $x,$ and contained in the
complement of the cut-locus of $x.$ Let $S^{n}$ be the $n-1$ dimensional unit
sphere in $n$ dimensions, and let $\theta\in S^{n}.$ Then

$\begin{split}&\log\left(\frac{d}{ds}\left(p_{t}(x,\exp_{x}(s\theta))\right)\right)\geq\\\
&C-\left(1+\frac{n}{2}\right)\log(t)-\frac{s}{4t}+\log(s)+\log(err),\end{split}$
(54)

where $\log(err)\to 0$ uniformly in $K$ and $C$ depends only on $M^{n}.$

###### Proof.

We start by quoting Corollary 2.29 from [9]:

Let $x,y\in M^{n}$ let

$E(x,y)=\frac{1}{2}dist^{2}(x,y).$ (55)

Then for fixed $y,$ the function $E(x,y)$ is smooth in $x$ for $x\in
M^{n}\setminus C_{M^{n}}(y),$ where $C_{M^{n}}(y)$ is the cut locus of $y.$
Moreover,

$\lim_{t\to 0}[\nabla L_{x}^{t}]=-\nabla_{y}E$ (56)

uniformly on compact subsets of $M^{n}\setminus C_{M^{n}}(y),$ where
$L_{x}^{t}$ is given by equation 47.

Proceeding as in the proof of lemma 3.5, we have that

$\nabla_{y}p_{t}=\frac{1}{t}p_{t}(x,y)\nabla L_{x}^{t}.$ (57)

Using lemma 4.7 of [6], we have that

$p_{t}(x,y)\geq\frac{C}{t^{\frac{n}{2}}}\exp\left(\frac{d^{2}(x,y)}{4t}\right),$
(58)

where $C$ depends only on $M^{n}.$

Note that

$\langle\nabla_{y}E(x,y),\frac{d}{ds}\exp_{x}(s\theta)\rangle=s$ (59)

if

$\exp_{x}(s\theta)=y.$ (60)

Therefore

$\begin{split}\frac{d}{ds}\left(p_{t}(x,\exp_{x}(s\theta))\right)&=\langle\nabla_{y}p_{t}(x,\exp_{x}(s\theta)),\frac{d}{ds}\exp_{x}(s\theta)\rangle\\\
&=\frac{1}{t}p_{t}(x,\exp_{x}(s\theta))\langle\nabla
L_{x}^{t},\frac{d}{ds}\exp_{x}(s\theta))\rangle\\\
&=\frac{1}{t}p_{t}(x,\exp_{x}(s\theta))\cdot s\cdot err\\\
&\geq\frac{1}{t}\frac{C}{t^{\frac{n}{2}}}\exp\left(\frac{d^{2}(x,\exp_{x}(s\theta)))}{4t}\right)\cdot
s\cdot err\\\
&=\frac{C}{t^{1+\frac{n}{2}}}\exp\left(\frac{s^{2}}{4t}\right)\cdot s\cdot
err,\end{split}$ (61)

where $err\to 1$ uniformly in $K$ in virtue of Corollary 2.29 of [9].

This implies that

$\begin{split}&\log\left(\frac{d}{ds}\left(p_{t}(x,\exp_{x}(s\theta))\right)\right)\geq\\\
&C-\left(1+\frac{n}{2}\right)\log(t)-\frac{s}{4t}+\log(s)+\log(err).\end{split}$
(62)

∎

We are now ready to give the proof of theorem 3.1.

###### Proof.

(Of theorem 3.1).

Note that $u_{t}$ can be written as

$u_{t}(y)=\sum_{j}^{N}w_{j}p_{t}(x_{j},y).$ (63)

Without loss of generality, assume $i=1.$ Then we have that

$\nabla
u_{t}(y)=w_{1}\nabla_{y}p_{t}(x_{1},y)+\sum_{j=2}^{N}w_{j}\nabla_{y}p_{t}(x_{j},y).$
(64)

Let

$M=\max d(x_{j},x_{k}).$ (65)

Using lemma 3.5, we have that

$\begin{split}&\left|\sum_{j=2}^{N}w_{j}\nabla_{y}p_{t}(x_{j},y)\right|\leq\\\
&N\max_{j}\left|w_{j}\nabla_{y}p_{t}(x_{j},y)\right|\leq\\\
&\exp\Bigg{(}C+\left(n-1\right)\log(t)-\frac{d^{2}(x,y)}{4t}+\\\ &\ \ \
\log\left(M\right)+\log(N)+\max_{j}\log(w_{j})\Bigg{)}.\end{split}$ (66)

Using lemma 3.3, we have that

$\begin{split}&\frac{d}{ds}\left(\sum_{j=2}^{N}w_{j}p_{t}(x_{j},\exp_{x_{1}}(s\theta))\right)\leq\\\
&\left|\sum_{j=2}^{N}w_{j}\nabla_{y}p_{t}(x_{j},\exp_{x_{1}}(s\theta))\right|\leq\\\
&\exp\Bigg{(}C+\left(n-1\right)\log(t)-\frac{s^{2}+\varphi(\epsilon)}{4t}\\\
&\ \ \ \ +\log\left(M\right)+\log(N)+\max_{j}\log(w_{j})\Bigg{)}.\end{split}$
(67)

Note that for any $\epsilon,\delta>0$ and $s>\delta$ there exists $T>0$ such
that, if $t<T$ then

$\begin{split}&C+\left(n-\frac{1}{2}\right)\log(t)-\frac{s^{2}+\varphi(\epsilon)}{4t}-\log(t)+\log\left(M\right)+\frac{1}{2}\log(t)+\log(N)+\max_{i}\log(w_{i})\leq\\\
&C-\left(1+\frac{n}{2}\right)\log(t)-\frac{\delta^{2}}{4t}+\log(\delta)+\min_{j}\log(w_{j})+\log(err).\end{split}$
(68)

Since $\varphi$ is continuous and monotone, $T$ can be chosen to be continuous
and monotone in each variable.

Also, by Lemma 3.6, we have that

$\begin{split}&C-\left(1+\frac{n}{2}\right)\log(t)-\frac{\delta^{2}}{4t}+\log(\delta)+\min_{j}\log(w_{j})+\log(err)\leq\\\
&C-\left(1+\frac{n}{2}\right)\log(t)-\frac{s^{2}}{4t}+\log(s)+\min_{j}\log(w_{j})+\log(err)\leq\\\
&\log\left(\frac{d}{ds}\left(w_{1}p_{t}(x,\exp_{x}(s\theta))\right)\right)\end{split}$
(69)

Therefore, for $\epsilon,\delta,t,T$ as above, we have that

$\left|\frac{d}{ds}\left(\sum_{j=2}^{N}w_{j}p_{t}(x_{j},\exp_{x_{1}}(s\theta))\right)\right|\leq\frac{d}{ds}w_{1}p_{t}(x_{1},\exp_{x_{1}}(s\theta)),$
(70)

and therefore

$\frac{d}{ds}\left(u_{t}(\exp_{x_{1}}(s\theta))\right)<0.$ (71)

∎

######  Remark 3.7.

The dependence of $T$ on the parameters is quite delicate. The choice of $T$
will depend on the number of points, the maximum distance between them, and
the weights. However, $T$ also depends on the function $\varphi,$ which is
closely linked to the geometry of the Voronoi cells.

##  4 Huygens principle and Voronoi cells in optics

For the sake of completeness, we also include some comments about hyperbolic
equations, although the following remarks are straightforward.

The same dynamic construction mentioned in the previous sections relative to
diffusive spherical fronts applies to the propagation of light in optics from
$n$ given sources utilizing the Huygens principle. This can be formalized
using the eikonal equation:

$|\nabla u|^{2}=1,$ (72)

on the plane. It is well known that the distance function $d$ constitutes a
solution. Given $n$ points on the plane, the distance function to this set,
$P=\\{x_{1},...x_{N}\\}$

$d(x,P),$ (73)

is a singular solution, which is smooth in the interior of a Voronoi cell and
whose singular set is precisely the boundary of the Voronoi cells. The same
observation holds in higher dimensions, and we could use the results in [10]
to conclude that the singular set is rectifiable.

##  5 Voronoi patterns, optimal transport, and the Monge-Ampere equation

This section is devoted to a link between Voronoi cells and the Monge-Ampere
equation. Based on a classical result in Optimal Transportation [2], we
construct a family of solutions to the Monge-Ampere equation. This family of
solutions is a counterexample to several results for the Monge Ampere equation
(particularly concerning regularity) when the hypotheses are slightly
modified.

######  Theorem 5.1.

Let $\Omega\subset\mathbf{R}^{n}$ be an open bounded set and let
$\\{w_{i}\\}_{i=1}^{N}\in\mathbf{R}.$ Let $\\{x_{i}\\}_{i=1}^{N}$ be $N$
distinct points in $\Omega.$ Fix $\lambda\in(0,1)$ and let

$\Phi_{\lambda}(x)=\frac{1}{2}|x|^{2}+\frac{\lambda-1}{2}\min_{i}\\{|x-x_{i}|^{2}+w_{i}\\}.$
(74)

Then $\Phi_{\lambda}$ is a strictly convex function that solves

$\det(\nabla^{2}\Phi_{\lambda})=\lambda^{n}$ (75)

in the Brennier sense (Definition 5.2).

######  Definition 5.2.

A convex function $\varphi$ is a solution of

$\det(\nabla^{2}\varphi)=\frac{f(x)}{g(\nabla\varphi(x))}$ (76)

in the Brennier sense if

$\int\xi(y)g(y)\,dy=\int\xi(\nabla\varphi(x))f(x)\,dx$ (77)

for all $\xi\in C^{\infty}$ and $g(\nabla\varphi(x))\neq 0.$

######  Lemma 5.3.

If $\varphi\in C^{\infty}$ is convex and satisfies 77 for all $\xi\in
C^{\infty}$, then $\varphi$ satisfies 76.

###### Proof.

The proof is essentially found in chapter $4$ of [15], but we include it here
for completeness.

If $\varphi\in C^{\infty}$ and is convex, we can perform a change of variables

$y=\nabla\varphi(x)$ (78)

on the left-hand-side and get

$\int\xi(\nabla\varphi(x))g(\nabla\varphi(x))\det(\nabla^{2}\varphi(x))\,dx=\int\xi(\nabla\varphi(x))f(x)\,dx$
(79)

for all $\xi\in C^{\infty}.$ Since this is true for all $\xi\in C^{\infty},$
we have that

$g(\nabla\varphi(x))\det(\nabla^{2}\varphi(x))=f(x).$ (80)

Since $g(\nabla\varphi(x))\neq 0,$ we have that

$\det(\nabla^{2}\varphi(x))=\frac{f(x)}{g(\nabla\varphi(x))}.$ (81)

∎

We will now prove Theorem 5.1

###### Proof.

(Of theorem 5.1)

For simplicity, we assume that $\Omega$ has volume $1,$ (otherwise, we could
apply this argument to a rescaled solution).

Let $V_{i}$ be the power set of $x_{i},$ defined as

$V_{i}=\\{x\in\Omega||x-x_{i}|^{2}+w_{i}<|x-x_{j}|^{2}+w_{j}\forall j\neq
i\\}.$ (82)

Assume without loss of generality that $V_{i}$ is nonempty for each $i.$

Consider the sets

$V_{i}^{\lambda}=\\{x_{i}+\lambda(x-x_{i})|x\in V_{i}\\}.$ (83)

Consider $\mu$ the uniform probability measure on $\Omega,$ and also the
probability measure $\nu,$ defined as

$d\nu=\frac{1}{\lambda^{n}}\sum_{i=1}^{N}\mathbf{1}_{V_{i}^{\lambda}}.$ (84)

We can easily verify that for each $x\in V_{i},$

$\nabla\Phi_{\lambda}=x+(\lambda-1)(x-x_{i}).$ (85)

and therefore $\nu$ is the push-forward of $\mu$ under $\nabla\Phi_{\lambda}:$

$\nu=\nabla\Phi_{\lambda}\\#\nu.$ (86)

Therefore $\Phi_{\lambda}$ solves 77 with

$\begin{split}f(x)&=1\\\
g(x)&=\frac{1}{\lambda^{n}}\sum_{i=1}^{N}\mathbf{1}_{V_{i}^{\lambda}}.\end{split}$
(87)

On the other hand, we also have that

$g(\nabla\Phi_{\lambda}(x))=\frac{1}{\lambda^{n}}$ (88)

for all $x.$ Therefore $\Phi_{\lambda}$ solves (75) in the sense of Brennier.
∎

######  Remark 5.4.

One way to get intuition on this result is to first solve for $\Phi_{\lambda}$
such that

$\nu=\nabla\Phi_{\lambda}\\#\nu.$ (89)

This leads to finding $\Phi_{\lambda}$ such that, for each $x\in V_{i},$

$\nabla\Phi_{\lambda}=x+(\lambda-1)(x-x_{i}).$ (90)

This leads to the formula for $\Phi_{\lambda}.$

As a corollary, we recover the solution to the corresponding problem in
Optimal Transportation.

######  Corollary 5.5.

Consider the Monge transportation problem

$\min_{T\\#\mu=\nu}\int_{\Omega}\left|x-T(x)\right|^{2}d\mu,$ (91)

with $\mu$ the uniform measure, and $\nu$ given by

$\nu=\sum_{i=1}^{N}\frac{1}{\lambda^{n}}\mathbf{1}_{V_{i}^{\lambda}}.$ (92)

Then by Brennier’s theorem, the solutions exists and is given by

$T=\nabla\Phi_{\lambda}.$ (93)

######  Remark 5.6.

Letting $\lambda$ tend to $0,$ we obtain that Corollary 5.5 remains true when

$\nu=\sum_{i=1}^{N}|V_{i}|\delta_{x_{i}}.$ (94)

This last remark recovers a classical result in Optimal Transportation [2].

$\Phi_{\lambda}$ is a counterexample to several results for Alexandrov
solutions of the Monge Ampere equation (see [11] for the definition of
Alexandrov solution), mainly concerning regularity. We list the results to
which $\Phi_{\lambda}$ is a counterexample and omit the proof.

######  Theorem 5.7.

([11] Theorem 1.2) Let $u$ be an Alexandrov solution to

$0<\lambda\leq Det(D^{2}u)\leq\Lambda<\infty$ (95)

in $B^{1}\subset\mathbf{R}^{n}.$ Then $u\in W^{2,1}_{loc}(B_{1})$

######  Theorem 5.8.

([11] Theorem 1.3) Assume that $u$ and $v$ are Alexandrov solutions to

$det(D^{2}u)=det(D^{2}v)=f$ (96)

in an open connected set $\Omega\subset\mathbf{R}^{n},$ with $f\in
C^{1,\alpha}(\overline{\Omega})$ strictly positive. If $u=v$ in an open subset
of $\Omega$ then $u=v$ in $\Omega.$

######  Theorem 5.9.

([3]) A strictly convex (Alexandrov) solution $u$ of

$det(D^{2}u)=d\,\mu,$ (97)

with $d\mu$ satisfying $0<\lambda\leq d\mu\leq\Lambda<\infty$ is
$C^{1,\alpha},$ with the $C^{1,\alpha}$ norm depending only on the Lipschitz
norm of $u$ and on its strict convexity. 111The exact statement of Theorem 2
is more general, see [3].

######  Theorem 5.10.

([14] Theorem 1.1) Suppose that $u:\Omega_{1}\to\mathbf{R}$ is a strictly
convex Alexandrov solution of

$\det(\nabla^{2}u)=f$ (98)

where $f$ satisfies

$0<\lambda\leq f\leq\Lambda<\infty.$ (99)

Then we have $u\in W^{2,1+\epsilon}_{loc}(\Omega)$ for some positive constant
$\epsilon$ which depends only on $n,\lambda,\Lambda.$

######  Corollary 5.11.

([5] Corollary 2.11) Let $\Omega$ be an open bounded set,
$g:\partial\Omega\to\mathbf{R}$ a continuous function, and $\nu$ a Borel
measure in $\Omega.$ Then there exists at most one convex function
$u:\Omega\to\mathbf{R}$ solving (in the Alexandriv sense) the problem

$\begin{split}&\det(\nabla^{2}u)=\nu\text{ in }\Omega\\\ &u=g\text{ in
}\partial\Omega\end{split}$ (100)

##  6 Conclusion

In this paper, we have established some new connections between partial
differential equations and Voronoi tessellations. An interesting direction for
future work is to consider general Riemannian manifolds and other nonlinear
equations, such as Fisher-KPP and the porous medium equations. Another
possible direction for future work is to extend the numerical simulations to a
broader class of colonization games.

Extensions to objects different from points, higher dimensions, manifolds, and
other functions(not only distance function).

## References

  * [1] F. Aurenhammer, R. Klein, and D.T. Lee.  Voronoi Diagrams and Delaunay Triangulations.  World Scientific, 2013.
  * [2] Franz Aurenhammer, Friedrich Hoffmann, and Boris Aronov.  Minkowski-type theorems and least-squares clustering.  Algorithmica, 20(1):61–76, 1998.
  * [3] Luis A Caffarelli.  Some regularity properties of solutions of monge ampere equation.  Communications on pure and applied mathematics, 44(8-9):965–969, 1991.
  * [4] Ben de Lacy Costello.  Calculating voronoi diagrams using simple chemical reactions.  Parallel Processing Letters, 25(01):1540003, 2015.
  * [5] Alessio Figalli.  The Monge–Ampère equation and its applications.  2017\.
  * [6] David García-Zelada et al.  Concentration for coulomb gases on compact manifolds.  Electronic Communications in Probability, 24, 2019.
  * [7] Elton P Hsu.  Stochastic analysis on manifolds.  Number 38. American Mathematical Soc., 2002.
  * [8] Lisa Larsson.  Generalized Voronoi Regions.  PhD thesis, McGill University Libraries, 2014.
  * [9] Paul Malliavin, Daniel W Stroock, et al.  Short time behavior of the heat kernel and its logarithmic derivatives.  Journal of Differential Geometry, 44(3):550–570, 1996.
  * [10] Carlo Mantegazza and Andrea Carlo Mennucci.  Hamilton-jacobi equations and distance functions on riemannian manifolds.  arXiv preprint math/0201296, 2002.
  * [11] Connor Mooney.  Partial regularity for singular solutions to the monge-ampère equation.  Communications on Pure and Applied Mathematics, 68(6):1066–1084, 2015.
  * [12] Mir Abolfazl Mostafavi, Leila Hashemi Beni, and Karine Hins-Mallet.  Representing dynamic spatial processes using voronoi diagrams: Recent developements.  In 2009 Sixth International Symposium on Voronoi Diagrams, pages 109–117, 2009.
  * [13] Peter Petersen.  Riemannian geometry, volume 171.  Springer, 2006.
  * [14] Thomas Schmidt.  W2, 1+ $\varepsilon$ estimates for the monge–ampère equation.  Advances in Mathematics, 240:672–689, 2013.
  * [15] Cédric Villani.  Topics in optimal transportation.  Number 58. American Mathematical Soc., 2003.
  * [16] Camille Wormser.  Generalized Voronoi Diagrams and Applications.  PhD thesis, Université Nice Sophia Antipolis, 2008.
  * [17] Zheng Xinqi, Wang Shuqing, Fu Meichen, Zhao Lu, and Yang Shujia.  Urban cluster layout based on voronoi diagram.  In Khaled Elleithy, editor, Advanced Techniques in Computing Sciences and Software Engineering, pages 13–18, Dordrecht, 2010. Springer Netherlands.