# Finite-size versus finite-temperature effects
in the critical long-range $O(N)$ model

Dario Benedetti  CPHT, CNRS, École polytechnique, Institut Polytechnique de
Paris, 91120 Palaiseau, France
Razvan Gurau  Heidelberg University, Institut für Theoretische Physik,
Philosophenweg 19, 69120 Heidelberg, Germany
Sabine Harribey  NORDITA, Stockholm University and KTH Royal Institute of
Technology,
Hannes Alfvéns väg 12, SE-106 91 Stockholm, Sweden

emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Davide Lettera  Heidelberg University, Institut für Theoretische Physik,
Philosophenweg 19, 69120 Heidelberg, Germany

###### Abstract

In this paper we consider classical and quantum versions of the critical long-
range $O(N)$ model, for which we study finite-size and finite-temperature
effects, respectively, at large $N$. First, we consider the classical
(isotropic) model, which is conformally invariant at criticality, and we
introduce one compact spatial direction. We show that the finite size
dynamically induces an effective mass and we compute the one-point functions
for bilinear primary operators with arbitrary spin and twist. Second, we study
the quantum model, mapped to a Euclidean anisotropic field theory, local in
Euclidean time and long-range in space, which we dub _fractional Lifshitz
field theory_. We show that this model admits a fixed point at zero
temperature, where it displays anisotropic Lifshitz scaling, and show that at
finite temperature a thermal mass is induced. We then compute the one-point
functions for an infinite family of bilinear scaling operators.

In both the classical and quantum model, we find that, as previously noted for
the short-range $O(N)$ model in [1], the large-$N$ two-point function contains
information about the one-point functions, not only of the bilinear operators,
but also of operators that appear in the operator product expansion of two
fundamental fields only at subleading order in $1/N$, namely powers of the
Hubbard-Stratonovich intermediate field.

###### Contents

  1. 1 Introduction
  2. 2 Generalized free field theory
    1. 2.1 Propagator on $S_{L}^{1}\times\mathbb{R}^{d-1}$
    2. 2.2 One-point functions of the minimal-twist operators
    3. 2.3 One-point functions of higher-twist operators
  3. 3 Classical long-range $O(N)$ model with a finite size
    1. 3.1 Schwinger-Dyson equation in the hyper-strip geometry
    2. 3.2 Long-range two-point function with mass
    3. 3.3 Finite-size mass solution
    4. 3.4 One-point functions of the minimal-twist operators
    5. 3.5 One-point functions of higher-twist operators
  4. 4 Quantum long-range $O(N)$ model at finite temperature
    1. 4.1 The massless free fractional Lifshitz field theory at zero temperature
    2. 4.2 The massless free fractional Lifshitz field theory at finite temperature
    3. 4.3 The interacting fractional Lifshitz field theory at zero temperature
    4. 4.4 The interacting fractional Lifshitz field theory at finite temperature
  5. 5 Conclusions
  6. A Bilinear operators
    1. A.1 Normalization of two-point functions
    2. A.2 OPE coefficient $f_{\phi\phi J}$
  7. References

##  1 Introduction

Quantum field theories in high energy physics are typically based on
Lagrangians which only depend on the fields and their first derivatives.
Higher-order derivatives in time are excluded in order to preserve unitarity,
and then by Lorentz invariance higher-order derivatives in space are excluded
as well. Fractional derivatives of order less than one are instead excluded by
locality, or the need of having a particle interpretation. Such requirements
might be dropped in Euclidean field theories, if the base space is interpreted
as representing only spatial directions and the theory itself as a model of
classical equilibrium statistical physics. It is within this perspective that
we can frame the study of non-local models containing fractional derivatives,
also known as long-range models. Such models have been realized experimentally
in recent years (see [2] for a review), but have been studied from a
theoretical standpoint since a long time because of several interesting
features: they can violate the Mermin-Wagner theorem and admit phase
transitions even in one dimension, as proved by Dyson for the long-range Ising
model [3]; they typically lead to one-parameter families of universality
classes (as found for example by two-loop [4, 5] and three-loop computations
[6]) and tuning such parameter one can study, at fixed dimension, phenomena
such as the crossover from a long-range to a short-range universality class
[7, 8, 9, 10]; one can rigorously construct fixed points and trajectories of
the renormalization group in three dimensions [11, 12, 13, 14, 15]; and so on.
Their realizations in the context of defect or boundary field theories have
also contributed to renewing theoretical interest in such models [16, 17].

In this paper, we want to address some aspects of finite-size and finite-
temperature effects, which in general are very important (see for example
[18]), in long-range models. In particular, we will study the long-range
$O(N)$ model with one compact direction. This can be seen as a special case of
finite-size effects in the corresponding classical statistical model [19].
However, by far the most common motivation for introducing only one compact
direction is the fact that the quantum statistical mechanics of a
$d$-dimensional system at finite temperature $T=1/\beta$ can be described as a
$(d+1)$-dimensional Euclidean field theory with compact Euclidean time of
length $\beta$, that is, on a spacetime $S^{1}_{\beta}\times\mathbb{R}^{d}$
with Euclidean signature (e.g. [20]). Such a description has two complementary
limits: the $\beta\to\infty$ limit describing a quantum system at zero-
temperature as a Euclidean field theory in $d+1$ dimensions, and the $\beta\to
0$ limit describing a classical system at equilibrium in $d$ space dimensions.
Therefore, in general with a given Euclidean field theory comes a choice of
interpreting it as representing either as classical or a quantum statistical
system. However, this is not always the case, as the quantum interpretation is
more constrained: we demand in particular a unitary and local time
evolution.111By local time evolution we mean one satisfying the Schrodinger
equation, or in a discrete setting, one admitting a transfer matrix
description. Unitarity, in this Euclidean setting, is equivalent to the
condition of reflection positivity, which implies positivity of the
Hamiltonian or transfer matrix. Such constraints imply in particular that
derivatives of order higher than two or fractional derivatives of any order
cannot be assigned to an imaginary-time direction. We stress this point
because of our focus on long-range models, where the fact that finite-
temperature and finite-size cases are structurally different seems to be
sometimes overlooked. As we will discuss below, the finite-size classical
system involves a compact direction among those appearing in the non-local
kernel, while the field theoretic description of the finite-temperature
quantum system requires the introduction of one extra compact direction with
local interaction. This observation is not really new, but with few exceptions
[21, 22] it is often only left implicit in the details of numerical
simulations [23, 24], and it has not been explored much at an analytical
level.

In order to clarify this point, let us consider a finite but continuous
system222See [25] for a similar story in the case of discrete spin systems.
with Hamiltonian:

$H(\Pi,Q)=\frac{1}{2}\sum_{i\in\Lambda}\Pi_{i}^{2}+U(Q)\;,\quad\qquad
U(Q)=-\sum_{i,j\in\Lambda}J_{ij}Q_{i}Q_{j}+\sum_{i\in\Lambda}V(Q_{i})\;,$
(1.1)

where $\Lambda\subset\mathbb{Z}^{d}$ is a finite $d$-dimensional lattice, and
$J_{ij}$ is not restricted to nearest neighbouring sites, but rather assumed
to be of the form:

$J_{ij}\sim\frac{1}{|i-j|^{d+2\zeta}}\;,\quad\qquad 0<\zeta<1\;,$ (1.2)

with the case $\zeta=1$ being equivalent to the nearest-neighbour (or short-
range) model.

In the classical canonical partition function, the configuration variables $Q$
and their conjugate momenta $\Pi$ are independent variables (e.g. taking
values in $\mathbb{R}$), and expectation values in the canonical ensemble are
computed as integrals over phase space with the Boltzmann weight. When
considering a $Q$-dependent observable $\mathcal{A}(Q)$, the integral over
$\Pi$ factors out:

$\langle\mathcal{A}(Q)\rangle_{\rm class.}\propto\int\,e^{-\beta
U(Q)}\mathcal{A}(Q)\,\prod_{i\in\Lambda}dQ_{i}\;.$ (1.3)

In the formal continuum limit $\Lambda\to\mathbb{R}^{d}$ and $Q_{i}\to\phi(x)$
this becomes a Euclidean functional integral, and $\beta$ is absorbed into the
definition of the couplings of the action $\beta U(Q)\to S[\phi]$:

$S[\phi]\sim-\int d^{d}xd^{d}y\,\frac{\phi(x)\phi(y)}{|x-y|^{d+2\zeta}}+\int
d^{d}x\,V(\phi)\;.$ (1.4)

In the quantum case, $\Pi$ and $Q$ are instead replaced by non-commuting
operators $\hat{\Pi}$ and $\hat{Q}$, and expectation values computed as
$\langle\mathcal{A}(\hat{Q})\rangle_{\rm
quant.}\propto\operatorname{Tr}[\mathcal{A}(Q)e^{-\beta\hat{H}}]$. Noticing
that $e^{-\beta\hat{H}}$ has the form of the time evolution operator
$e^{\imath t\hat{H}}$ with imaginary time $t=\imath\beta$, one then proceeds
using the Trotter formula to split the interval $[0,\beta]$ in $n$ steps, and
introducing resolutions of the identity in both the $|\Pi\rangle$ and
$|Q\rangle$ bases at each step, thus reducing again to commuting variables.
Integrating over $\Pi$ one reaches:

$\langle\mathcal{A}(\hat{Q})\rangle_{\rm
quant.}\propto\int\,e^{-\frac{1}{2\epsilon}\sum_{<k,l>}\sum_{i\in\Lambda}(Q_{i,k}-Q_{i,l})^{2}-\epsilon\sum_{k}U(Q_{k})}\mathcal{A}(Q)\,\prod_{k=1}^{n}\prod_{i\in\Lambda}dQ_{i,k}\;,$
(1.5)

where $\epsilon=\beta/n$ and the lattice now includes a (Euclidean) time
direction, with periodic boundary condition $Q_{i,n+1}=Q_{i,1}$, along which
the action is short-range ($<k,l>$ denotes nearest-neighbours). In the
continuum limit, this becomes again a Euclidean field theory, but on a
$(d+1)$-dimensional space, and with an anisotropy in the (Euclidean) time
direction, if $\zeta\neq 1$.333 Notice that such a construction goes through
in a very similar fashion for any Hamiltonian of the form
$H(\Pi,Q)=K(\Pi)+U(Q)$, with the only difference that if the integral in $\Pi$
is not Gaussian then the result will not be Gaussian in $(Q_{i,k}-Q_{i,l})$;
but it will still be short-ranged in time. For a case with mixing terms, see
for example the rotor model [26].

The main point of this example is to highlight that it would be naive and
problematic to interpret as a finite-temperature quantum system the model
(1.4) with a compactified direction. In particular, we would face a problem in
reversing the path integral construction and identifying a quantum
Hamiltonian.

In this paper, we will be interested in distinguishing and studying the same
kind of long-range model in two different regimes: as a classical statistical
field theory with one compact spatial direction and as a quantum statistical
field theory at finite temperature (i.e. with compact Euclidean time). In
other words, we will consider a continuous version of the model (1.5), defined
in general on $S^{1}_{\beta}\times S^{1}_{L}\times\mathbb{R}^{d-1}$, with
long-range interaction in the $d$ spatial directions and standard local
interaction in the imaginary-time direction, i.e. the coordinate on
$S^{1}_{\beta}$ (see equation (4.1) below for the explicit form of the
action). However, we will never consider explicitly the case with both $\beta$
and $L$ finite: we will start with the classical system at $\beta\to 0$ and
$0<L<\infty$ (i.e. a purely long-range model on
$S^{1}_{L}\times\mathbb{R}^{d-1}$), and then study the finite-temperature
quantum system at $0<\beta\leq\infty$ and $L\to\infty$ (i.e. a short-range in
time and long-range in space model on $S^{1}_{\beta}\times\mathbb{R}^{d}$).

In both cases, we will focus on the critical model. As usual, this will
correspond to tuning the non-compactified version of the model (i.e. on
$\mathbb{R}^{d}$ or $\mathbb{R}^{d+1}$) to the infrared fixed point of the
renormalization group, where the model is scale invariant. In the isotropic
case (the classical model), it is known that such invariance is promoted to
full conformal invariance [16], and as we review below, an interesting
question is what happens to the conformal structure of the theory once a
spatial direction is compactified. In the anisotropic case (the quantum
model), the situation is more complicated, as the $SO(d+1)$ group of rotations
is broken to $SO(d)$ even at zero temperature. The scaling invariance is in
this case of the Lifshitz type (see below), and at present it is not clear to
us whether this will be enhanced to some form of anisotropic conformal
invariance (see [27] for examples of that).

#### Finite-size effects in CFTs.

In general, a conformal field theory (CFT) in $\mathbb{R}^{d}$ is
unambiguously defined by its CFT data, that is, the spectrum of scaling
dimensions of its primary operators and the set of all the operator product
expansion (OPE) coefficients. One-point functions vanish identically, two-
point functions are fixed up to normalization by the spin and scaling
dimensions of the operators, and three-point functions are fixed by the same
data plus the OPE coefficients. For higher $n$-point functions no new data are
needed, as in principle they can be reconstructed by repeated use of the OPE.
Things change when one or more spatial directions are compactified. In fact,
such compactifications obviously introduce a scale (the size $L$ of the
compact direction) that breaks both scale and rotation invariance. Luckily,
not all the beauty of CFT is lost: if the system is still tuned to the same
critical point corresponding to the CFT of the non-compact case, we expect
that the compactification will not affect small-distance444Smaller than the
smallest size of the compact directions. properties, such as the OPE and the
short-scale divergence of the two-point functions (hence the scaling
dimensions). As a consequence, the reduction of $n$-point functions to lower
correlators is still possible, at least for configurations of the $n$ points
inside a ball of radius smaller than $L$. On the other hand, this reasoning
does not apply to one-point functions: since translations are preserved by
toroidal compactifications, one-point functions have to be constant, and thus
do not admit a short-distance limit. Therefore, as a consequence of the broken
scale invariance, the one-point function of an operator $\mathcal{O}$ can
acquire a non-vanishing value, proportional to $L^{-\Delta_{\mathcal{O}}}$,
where $\Delta_{\mathcal{O}}$ is the scaling dimension of the operator. The
non-vanishing coefficients of one-point functions thus become new CFT data,
that permeate through OPE to all the higher $n$-point functions.

Such aspects of finite-size/temperature effects in CFTs have been discussed in
[28, 29, 1, 30, 31, 32], and we follow in particular the presentation of [1].
On $S^{1}_{L}\times\mathbb{R}^{d-1}$, conformal symmetry constrains the one-
point functions to take the form:

$\langle\mathcal{O}^{\mu_{1}\cdots\mu_{J}}\rangle=\frac{b_{\mathcal{O}}}{L^{\Delta_{\mathcal{O}}}}\left(e^{\mu_{1}}\cdots
e^{\mu_{J}}-\text{traces}\right)\;,$ (1.6)

where $\Delta_{\mathcal{O}}$ and $J$ are the scaling dimension and spin (order
of the symmetric traceless tensor representation of $SO(d)$) of the operator
$\mathcal{O}$, and $e$ is the unit vector in the compact direction. The $L$
dependence in (1.6) is fixed but the coefficient $b_{\mathcal{O}}$ is not: it
is the new CFT data we discussed above. The non-vanishing of $b_{\mathcal{O}}$
affects the OPE of the two-point functions, which otherwise would only include
the identity operator. Indeed, starting with a two-point function in direct
space:

$G(y,\mathbf{x})=\langle\phi(y,\mathbf{x})\phi(0,0)\rangle\;,\qquad y\in
S^{1}_{L}\;,\;\mathbf{x}\in\mathbb{R}^{d-1}\;,$ (1.7)

and using the OPE for the two fields $\phi$ (which is convergent in the region
$x^{2}=y^{2}+\mathbf{x}^{2}<L^{2}$, as argued in [1]) one obtains:

$G(y,\mathbf{x})=\sum_{\mathcal{O}\in\phi\times\phi}\frac{f_{\phi\phi\mathcal{O}}}{c_{\mathcal{O}}}|x|^{\Delta_{\mathcal{O}}-2\Delta_{\phi}-J}x_{\mu_{1}}\cdots
x_{\mu_{J}}\langle\mathcal{O}^{\mu_{1}\cdots\mu_{J}}(0)\rangle\;,$ (1.8)

where $f_{\phi\phi\mathcal{O}}$ is the OPE coefficient and $c_{\mathcal{O}}$
is the normalization coefficient of the two-point function of $\mathcal{O}$.
As the tensor structure in (1.6) is symmetric and traceless one can use [33,
34]:

$|x|^{-J}x_{\mu_{1}}\cdots x_{\mu_{J}}\left(e^{\mu_{1}}\cdots
e^{\mu_{J}}-\text{traces}\right)=C^{\left(\frac{d-2}{2}\right)}_{J}(\eta)\frac{J!\Gamma(\frac{d-2}{2})}{2^{J}\Gamma(\frac{d-2}{2}+J)}\;,$
(1.9)

where $\eta=y/|x|$ and $C^{\left(\frac{d-2}{2}\right)}_{J}(\eta)$ is a
Gegenbauer Polynomial and rewrite the OPE of the two-point function as:

$G(y,\mathbf{x})=\sum_{\mathcal{O}\in\phi\times\phi}\frac{a_{\mathcal{O}}}{L^{\Delta_{\mathcal{O}}}}C^{\left(\frac{d-2}{2}\right)}_{J}(\eta)|x|^{\Delta_{\mathcal{O}}-2\Delta_{\phi}}\,,\qquad
a_{\mathcal{O}}\equiv\frac{b_{\mathcal{O}}f_{\phi\phi\mathcal{O}}}{c_{\mathcal{O}}}\frac{J!\Gamma(\frac{d-2}{2})}{2^{J}\Gamma(\frac{d-2}{2}+J)}\;.$
(1.10)

Note that in the infinite-size case one can go through the same steps and
write an OPE expansion for the two-point function, but since only the term
corresponding to the identity operator (with $\Delta_{\mathcal{O}}=0$) has a
non-vanishing expectation value in the limit $L\rightarrow\infty$ one recovers
the standard conformal form of the two-point function.

In the finite-size case one can in principle read off the coefficients
$a_{\mathcal{O}}$ by expanding the two-point function as a sum over terms of
the form
$C^{\left(\frac{d-2}{2}\right)}_{J}(\eta)|x|^{\Delta-2\Delta_{\phi}}$.
Conversely, once the $a_{\mathcal{O}}$ coefficients are known one can
reconstruct the two-point function. In practice, both of these reconstructions
are typically hard, as they require knowing exactly either the two-point
function or the CFT data. However, they can become useful in special
circumstances in which either of these pieces of information is available, as
in the large-$N$ limit that we will employ here.

The new data $b_{\mathcal{O}}$, for $\mathcal{O}$ in the OPE of two
fundamental fields of the short-range $O(N)$ model, have been computed in [1,
32] at leading and subleading order of the $1/N$ expansion. Here we will
tackle similar computations for the long-range case.

#### Fractional Lifshitz field theories.

As argued above, the quantum version of a long-range model, at zero or non-
zero temperature, is mapped to a Euclidean field theory with a strong
anisotropy. That is, it is not just the compactification of the Euclidean time
direction that breaks rotation invariance: the invariance is broken already by
the local dynamics. The model that arises in such case (see our action (4.1),
as well as [21, 22]) is a kind of Lifshitz field theory, that is, the type of
theory used as Landau free energy of a Lifshitz point [35]. Lifshitz points
are characterized by a distinctive phase diagram, including a modulated phase,
and an anisotropic scale invariance:

$\tau\to\Omega^{z}\tau,\quad x_{i}\to\Omega x_{i}\;,$ (1.11)

where we denoted $\tau$ the Euclidean time and the exponent $z$ is known as
_dynamic critical exponent_ , because of the similarity with the scale
invariance of dynamic critical phenomena [36]. Typical Lifshitz field theories
are scalar field theories whose kinetic term has two time derivatives555In the
classical interpretation of Lifshitz points, “time” is actually one spatial
direction of a spatially anisotropic classical system [35]. However, an
interpretation as time has also appeared in the context of quantum Lifshitz
points [37] or of Lorentz-violating quantum field theories [38]. and $2z$
spatial derivatives, with $z$ typically being an integer (most commonly
$z=2$). Our action will be of a similar form, but with non-integer $z=\zeta$,
fixed by the corresponding classical long-range model (1.4), and thus we will
call it _fractional Lifshitz field theory_ (FLFT).

Given that most of the literature on Lifshitz field theories is focused on
actions with $z=2$, or other integer values, there are several open questions
for our FLFT, even at zero temperature. For example, as we mentioned above, we
do not currently know whether at the interacting fixed point scale invariance
is enhanced to some form of anisotropic conformal invariance. Moreover, we do
not know yet if at such fixed point the dynamic critical exponent takes a
value different from $\zeta$,666As a reminder, in the isotropic long-range
model there is no anomalous dimension. and we also do not know if above some
$\zeta^{*}<1$ there is a crossover to the short-range universality class
($\zeta=1$), as in the isotropic long-range model. We hope to address these
questions in a separate paper (the last two require in particular to go beyond
the leading order of the large-$N$ approximation, to which we stick in the
present work), while here we wish to use this model for a comparative study of
some finite-temperature versus finite-size effects.

#### Plan of the paper.

The paper is organised as follows. In section 2, we start by studying the case
of a generalized free field theory on $S_{L}^{1}\times\mathbb{R}^{d-1}$, as a
warm up. In the process we need to deal with higher twist operators that are
absent in the short-range case. We reproduce the results of [1] for the
$a_{\cal O}$ coefficients in the expansion of the two-point function for
operators with arbitrary spin and twist by a different technique. In section
3, we study the long-range $O(N)$ model with finite size. We prove the
existence of a finite-size mass and compute the coefficients $b_{\mathcal{O}}$
for spinning operators with minimal twist. We also compute the one-point
functions for operators with arbitrary twist, and deduce their $a_{\cal O}$
coefficients. In section 4, we switch to the case of the quantum long-range
model at finite temperature, described by a FLFT on
$S^{1}_{\beta}\times\mathbb{R}^{d}$. We show that at zero temperature this
model admits an interacting infrared fixed point. We then show that a thermal
mass is generated when putting the model on
$S^{1}_{\beta}\times\mathbb{R}^{d}$ and we derive the one-point functions for
an infinite family of scaling operators at finite temperature.

In appendix A, we give some details on the computation of the zero temperature
CFT data for bilinear operators, that is the normalizations of the two-point
functions and the OPE coefficients $f_{\phi\phi J}$.

##  2 Generalized free field theory

As a warm up, we first consider the case of a conformal generalized free field
theory (GFFT), i.e. a massless Gaussian theory with a non-local inverse
covariance, also known as mean field theory [39] or fractional Gaussian field
[40]. It is worth discussing it in some detail because it is the simplest case
of (non-local) CFT, and because typical long-range models can be defined as
perturbations of a GFFT.

We begin by recalling some basic facts in the non-compact case, viewed as a
classical statistical field theory at criticality, that is, the GFFT defined
on $\mathbb{R}^{d}$, with Euclidean signature. By definition this is a CFT
whose only non-vanishing connected $n$-point function is the two-point
function, given by the covariance $C(x)$ of a Gaussian measure, which however
has a non-canonical scaling exponent $\Delta_{\phi}>d/2-1$:

$C(x)=\frac{c(\Delta_{\phi})}{|x|^{2\Delta_{\phi}}}\;,\qquad
c(\Delta_{\phi})=\frac{\Gamma(\Delta_{\phi})}{2^{d-2\Delta_{\phi}}\pi^{d/2}\Gamma(\frac{d}{2}-\Delta_{\phi})}\;.$
(2.1)

Writing $\Delta_{\phi}=d/2-\zeta$, such GFFT can be obtained from a functional
integral with the action:777We typically assume $0<\zeta<1$. The restriction
to $\zeta<1$ is imposed to preserve reflection positivity (unitarity in
Lorentzian signature), but also because $\zeta>1$ would correspond to a strong
short-range rather than long-range action, and moreover the operator with
$\zeta=1$ would in that case be a relevant perturbation. The restriction to
$\zeta>0$ is instead chosen to avoid a strong long-range action, with its
associated unusual thermodynamic features [41].

$S_{\rm GFFT}[\phi]=\frac{1}{2}\int{\rm
d}^{d}x\,\phi(x)(-\partial^{2})^{\zeta}\phi(x)\;,$ (2.2)

where $\partial^{2}=\partial_{\mu}\partial^{\mu}$ is the Laplace operator. The
fractional power of the Laplacian can be defined in many equivalent ways [42],
and in particular as the inverse of the operator displayed in (2.1), which has
a kernel of the same form, but with the replacement $\Delta_{\phi}\to
d-\Delta_{\phi}$, leading to the kinetic term of (1.4).888Notice that
$c(d-\Delta_{\phi})=c(d/2+\zeta)<0$ for $0<\zeta<1$, thus matching the minus
sign in (1.4). Notice also that despite such minus sign, the GFFT action is
positive definite, as best understood in momentum space: this is because the
position space representation of the fractional Laplacian is singular and
needs a subtraction, which is not written explicitly in (1.4) (see footnote 18
of [43]). The easiest definition is of course in Fourier space, where
$(-\partial^{2})^{\zeta}$ is defined as the multiplication operator
$(p^{2})^{\zeta}$, which is the inverse of the Fourier transform of (2.1) (the
coefficient $c(\Delta_{\phi})$ is chosen precisely to have a unit coefficient
in momentum space, $\tilde{C}(p)=1/(p^{2})^{\zeta}$).

Since the theory is Gaussian, the scaling dimensions of local composite
operators are simply given by $n_{\phi}\Delta_{\phi}+n_{\partial}$, where
$n_{\phi}$ is the number of $\phi$ fields and $n_{\partial}$ the number of
derivatives. This is the same as in the local free field theory, but the
operator content is slightly different in the two cases: in the GFFT case,
objects like the energy-momentum tensor or the field equations are non-local,
hence they do not appear in the algebra of local operators. As a consequence,
we cannot use the Schwinger-Dyson equation:

$\langle\frac{\delta
S}{\delta\phi(x_{0})}\phi(x_{1})\cdots\phi(x_{n})\rangle=\sum_{i=1}^{n}\delta(x_{0}-x_{i})\langle\prod_{j\neq
i}^{1,\ldots,n}\phi(x_{j})\rangle\;,$ (2.3)

to conclude that insertions of $(-\partial^{2})\phi$ in any $n$-point function
(at separate points) vanish.999In fact it is straightforward to check that
they do not vanish: using the fact that the theory is Gaussian we reduce
$\langle[\partial^{2}\phi](x_{0})\phi(x_{1})\cdots\phi(x_{n})\rangle$ to
products of two-point functions, and we observe that
$\langle[\partial^{2}\phi](x_{0})\phi(x_{i})\rangle=\partial^{2}_{x_{0}}\langle\phi(x_{0})\phi(x_{i})\rangle\neq\delta(x_{0}-x_{i})$
for $\Delta_{\phi}\neq d/2-1$.

Another consequence of Gaussianity is that in the OPE of two fundamental
fields, $\phi\times\phi$, only bilinear operators in $\phi$ appear, that is,
operators of the schematic form:

$[\phi\phi]_{k,\mu_{1}\cdots\mu_{J}}(x)=\,:\phi(x)(\partial^{2})^{k}\left(\partial_{\mu_{1}}...\partial_{\mu_{J}}-\text{traces}\right)\phi(x):\;,$
(2.4)

where $J\in 2\mathbb{N}_{0}$, $k\in\mathbb{N}_{0}$, and the colon notation
implies a renormalization of the composite operator (normal ordering in this
case). These operators have dimension:

$\Delta_{k,J}=2\Delta_{\phi}+J+2k\;.$ (2.5)

Introducing the twist $\tau_{k,J}=\Delta_{k,J}-J$, we see that the operators
with $k=0$ are the minimal-twist operators, also known as double-twist
operators [44]. Once more, the appearance of the higher-twist operators, with
$k>0$, marks a difference with respect to the local free theory. Another way
to see this difference is by using the conformal partial wave expansion [45,
46, 47, 48] of the four-point function. A conformal four-point function of
identical scalars can be written as follows:

$\begin{split}\langle\phi(x_{1})\phi(x_{2})\phi(x_{3})\phi(x_{4})\rangle=&\langle\phi(x_{1})\phi(x_{2})\rangle\langle\phi(x_{3})\phi(x_{4})\rangle+\\\
&\ \
+\sum_{J}\int_{\frac{d}{2}-\imath\infty}^{\frac{d}{2}+\imath\infty}\frac{dh}{2\pi\imath}\
\frac{1}{1-k(h,J)}\,\mu^{d}_{\Delta_{\phi},J}(h,J)\,G^{\Delta_{\phi}}_{h,J}(x_{i})\;,\end{split}$
(2.6)

where $G^{\Delta_{\phi}}_{h,J}(x_{i})$ is the conformal block,
$\mu^{d}_{\Delta_{\phi},J}(h,J)$ is a known
function101010$\mu^{d}_{\Delta_{J}}(h,J)$ is the Plancherel weight $\rho(h,J)$
multiplied by the right normalization of three-point functions, necessary to
make them an orthonormal basis of the appropriate space of bilocal functions.
and $k(h,J)$ is the eigenvalue of the Bethe-Salpeter four-point kernel. In the
case of a free theory, we have $k(h,J)=0$ and one can close the contour of
integration of $h$ on the right of the complex plane and pick the poles of
$\mu^{d}_{\Delta_{\phi},J}(h,J)\,G^{\Delta_{\phi}}_{h,J}(x_{i})$ on the right
of $\frac{d}{2}$. The location of such poles are the scaling dimensions of
bilinear primaries and the residues are their OPE coefficients with two
$\phi$. In [49] it is shown that for $\Delta_{\phi}>\Delta_{SR}=\frac{d-2}{2}$
the location of poles are at $h=\Delta_{k,J}=2\Delta_{\phi}+J+2k$, but in the
limit $\Delta_{\phi}\rightarrow\Delta_{SR}$ all the residues with $k\geq 1$ go
to zero, hence those poles disappear.

As written in (2.4), the bilinear operators are in a schematic form because we
have not taken into account the primary constraint, ensuring that they
transform appropriately under the full conformal group. In the short-range
case the precise form of bilinear primaries has been extensively studied [45,
50, 51] and we review it in appendix A. However, as we already emphasized
above, in that case $k$ is forced to be zero by the equation of motion. As a
consequence, the operators with $k>0$, relevant for the GFFT, have not been
studied as much in the literature, and to the best of our knowledge their
general form is not known, although a recurrence relation was derived in [52].
Below we give the explicit expressions for some primary operators with $k=1$
and $k=2$:

$\begin{split}&\frac{[\phi\phi]_{1,0}(x)}{\mathcal{N}_{1,0}}=:\phi\partial^{2}\phi(x):+\left(\frac{\Delta_{SR}}{\Delta_{\phi}}-1\right):\partial_{\alpha}\phi\partial_{\alpha}\phi(x):\;,\\\
&\frac{[\phi\phi]_{2,0}(x)}{\mathcal{N}_{2,0}}=:\phi\partial^{2}\partial^{2}\phi(x):+\left(4\frac{\left(\Delta_{SR}-\Delta_{\phi}-1\right)}{\Delta_{\phi}}\right):\partial_{\alpha}\phi\partial_{\alpha}\partial^{2}\phi(x):\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(\frac{\left(\Delta_{\phi}-\Delta_{SR}+1\right)\left(\Delta_{\phi}(\Delta_{\phi}+4)+2-\Delta_{SR}\right)}{\Delta_{\phi}(\Delta_{\phi}+1)\left(\Delta_{\phi}-\Delta_{SR}\right)}\right):\partial^{2}\phi\partial^{2}\phi(x):\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(\frac{2\left(\Delta_{\phi}-\Delta_{SR}\right)\left(\Delta_{\phi}-\Delta_{SR}+1\right)}{\Delta_{\phi}(\Delta_{\phi}+1)}\right):\partial_{\alpha}\partial_{\beta}\phi\partial_{\alpha}\partial_{\beta}\phi(x):\;,\\\
&\frac{[\phi\phi]_{1,\mu_{1}\mu_{2}}(x)\xi^{\mu_{1}}\xi^{\mu_{2}}}{\mathcal{N}_{1,2}}=:\phi\partial^{2}(\imath\xi\cdot\partial)^{2}\phi(x):+\left(\frac{(\Delta_{\phi}+2)(\Delta_{\phi}+3)}{\Delta_{\phi}(\Delta_{\phi}+1)}\right):\partial^{2}\phi(\imath\xi\cdot\partial)^{2}\phi(x):\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(\frac{2\left(\Delta_{\text{SR}}-\Delta_{\phi}\right)}{\Delta_{\phi}}\right):\partial_{\alpha}\phi(\imath\xi\cdot\partial)^{2}\partial_{\alpha}\phi(x):-\left(\frac{2(\Delta_{\phi}+3)}{\Delta_{\phi}}\right):\phi\partial^{2}(\imath\xi\cdot\partial)\phi(\imath\xi\cdot\partial)\phi(x):\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\left(\frac{2(\Delta_{\phi}+2)\left(\Delta_{\phi}-\Delta_{\text{SR}}\right)}{\Delta_{\phi}(\Delta_{\phi}+1)}\right):\partial_{\alpha}(\imath\xi\cdot\partial)\phi\partial_{\alpha}(\imath\xi\cdot\partial)\phi(x):\;.\end{split}$
(2.7)

where we recall that $\Delta_{SR}=\frac{d-2}{2}$ is the scaling dimension of
the field in the short-range case, $\xi$ is a complex vector with $\xi^{2}=0$
and $\mathcal{N}_{k,J}$ is a normalization that we are free to fix at
convenience.

We notice that for $\Delta_{\phi}=\Delta_{SR}$ only terms with at least one
$\partial^{2}\phi$ survive, but, as explained above, such terms vanish in the
short-range case due to the field equations. Also, we notice that the
operators with $k>0$ do not coincide with the trace part subtracted in the
spinning operators. For example, from the formulas of appendix A we have:

$\begin{split}\frac{[\phi\phi]_{0,\mu\nu}(x)}{\mathcal{N}_{0,2}}=&:\phi\partial_{\mu}\partial_{\nu}\phi(x):-\frac{\Delta_{\phi}+1}{\Delta_{\phi}}:\partial_{\mu}\phi\partial_{\nu}\phi(x):-\frac{\delta_{\mu\nu}}{d}\left(:\phi\partial^{2}\phi(x):-\frac{\Delta_{\phi}+1}{\Delta_{\phi}}:\partial_{\alpha}\phi\partial_{\alpha}\phi(x):\right)\;,\end{split}$
(2.8)

and we see that the trace part is a different linear combination than the one
appearing in $[\phi\phi]_{1,0}(x)$.

###  2.1 Propagator on $S_{L}^{1}\times\mathbb{R}^{d-1}$

As anticipated in the introduction, we now place the model in a hyper-strip
geometry with one compactified dimension of length $L$ and periodic boundary
conditions. The propagator in momentum space is just:

$\tilde{C}_{L}(q_{n},\mathbf{p})=\frac{1}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}}\;,$
(2.9)

with $q_{n}=2\pi n/L$, $n\in\mathbb{Z}$, and $\mathbf{p}\in\mathbb{R}^{d-1}$.
The subscript $L$ signifies that the covariance corresponds to a compactified
strip geometry, and in the non-compact case we recover the GFFT covariance,
$C_{\infty}\equiv C$. In order to Fourier transform back to direct space, we
will introduce a fictitious integration in $q$ via $1=\int
dq\,\delta(q-q_{n})$ and use the Poisson summation formula for distributions:

$\frac{1}{L}\sum_{n}\delta(q-q_{n})=\frac{1}{2\pi}\sum_{m}e^{\imath Lmq}\;.$
(2.10)

Considering $q$ as the first component of a $d$-dimensional continuous
momentum $p=(q,\mathbf{p})$, we get:

$C_{L}(y,\mathbf{x})=\frac{1}{L}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\frac{e^{\imath\mathbf{p}\cdot\mathbf{x}+\imath
q_{n}y}}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}}=\sum_{m=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot x_{m}}}{(p^{2})^{\zeta}}\;,$ (2.11)

where $y\in S_{L}^{1}$, $\mathbf{x}\in\mathbb{R}^{d-1}$, and
$x_{m}=(y+mL,\mathbf{x})$. At this stage we could proceed by using a Schwinger
parametrization for the denominator (as for example in [6]). However this does
not generalize nicely to the massive case of the next section, hence we
evaluate here the two-point function in position space by a different method.

We write the momentum space propagator as a Stieltjes transform (or Källén-
Lehmann representation):

$\frac{1}{(p^{2})^{\zeta}}=\int_{0}^{+\infty}ds\,\frac{\rho(s)}{p^{2}+s}\;,\qquad\rho(s)=\frac{\sin(\pi\zeta)}{\pi}s^{-\zeta}\;.$
(2.12)

The Fourier transform in (2.11) is then reduced to the standard one for the
short-range massive scalar free propagator, and we thus obtain:

$\begin{split}C_{L}(y,\mathbf{x})&=\frac{\sin(\pi\zeta)}{(2\pi)^{d/2}\pi}\sum_{m=-\infty}^{+\infty}\int_{0}^{+\infty}ds\,s^{-\zeta}\,\left(\frac{\sqrt{s}}{|x_{m}|}\right)^{d/2-1}K_{\frac{d}{2}-1}(\sqrt{s}|x_{m}|)\\\
&=\frac{\sin(\pi\zeta)}{(2\pi)^{d/2}\pi}\sum_{m=-\infty}^{+\infty}\frac{1}{|x_{m}|^{d-2\zeta}}\int_{0}^{+\infty}ds\,s^{-\zeta+d/4-1/2}\,K_{\frac{d}{2}-1}(\sqrt{s})\\\
&=\frac{\sin(\pi\zeta)\Gamma(1-\zeta)\Gamma(d/2-\zeta)}{2^{2\zeta}\pi^{d/2}\pi}\sum_{m=-\infty}^{+\infty}\frac{1}{|x_{m}|^{d-2\zeta}}\\\
&=c(\Delta_{\phi})\sum_{m=-\infty}^{+\infty}\frac{1}{|x_{m}|^{2\Delta_{\phi}}}\;,\end{split}$
(2.13)

where $K_{\frac{d}{2}-1}(x)$ is the modified Bessel function of second kind.
This might look like a useless exercise, as we knew already the result of the
Fourier transform in (2.11), but its usefulness lies in highlighting a
strategy that will be fruitful in the massive case.

###  2.2 One-point functions of the minimal-twist operators

We are now interested in computing the one-point function of bilinear primary
operators introduced in (4.12). The most convenient way to deal with the
tracelessness condition of higher-spin operators is by means of the index-free
formalism of [45, 50]: define an operator $\mathcal{O}_{k,J}(x,\xi)$ by
contracting each derivative with a complex null vector $\xi=(\chi,\bm{\xi})$
such that $\bm{\xi}^{2}=-\chi^{2}$ (see appendix A for more details). For now,
we will consider $k=0$.

When computing in momentum space the one-point function of the primary
operator defined in (A.1) we just have to substitute
$\overset{\scriptscriptstyle\leftarrow}{\partial}\rightarrow-\imath p$ and
$\overset{\scriptscriptstyle\rightarrow}{\partial}\rightarrow\imath p$. As a
result, only the $n=0$ term in the sum (A.2) survives, the prefactor
simplifies to one and we get:

$\begin{split}\langle\mathcal{O}_{0,J}(x,\xi)\rangle&\equiv\xi^{\mu_{1}}\cdots\xi^{\mu_{J}}\langle[\phi_{a}\phi_{a}]_{0,\mu_{1}\cdots\mu_{J}}(x)\rangle\\\
&=\frac{1}{L}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\frac{(\imath\mathbf{p}\cdot\bm{\xi}+\imath
q_{n}\chi)^{J}}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}}-\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath
p\cdot\xi)^{J}}{(p^{2})^{\zeta}}\;,\end{split}$ (2.14)

where the subtraction term comes from the normal ordering of the operator, as
implemented before compactifying one direction, and an ultraviolet
regularization is implicitly assumed. As written, it is not obvious that the
normal ordering is sufficient in making the one-point function finite when we
remove the regularization. In order to see that indeed it is finite, it is
convenient to use the Poisson formula, and rewrite:

$\begin{split}\langle\mathcal{O}_{0,J}(x,\xi)\rangle&=\sum_{m=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath\mathbf{p}\cdot\bm{\xi}+\imath
q\chi)^{J}\,e^{\imath
Lmq}}{(\mathbf{p}^{2}+q^{2})^{\zeta}}-\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath
p\cdot\xi)^{J}}{(p^{2})^{\zeta}}\\\ &=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath\mathbf{p}\cdot\bm{\xi}+\imath
q\chi)^{J}\,e^{\imath Lmq}}{(\mathbf{p}^{2}+q^{2})^{\zeta}}\;.\end{split}$
(2.15)

The oscillating factor renders each term in the sum finite (as can be shown by
deforming the contour of integration of $q$ in the complex plane), and the
normal ordering removes the only term without oscillating factor. Actually,
the subtracted term vanishes for $J>0$ in any regularization preserving
rotation invariance, such as a cutoff on $p^{2}$ (via the insertion in the
integrand of a cutoff function $f(p^{2}/\Lambda^{2})$ decaying sufficiently
fast at infinity): by symmetry reasons the resulting integral over $p$ must be
proportional to a linear combination of products of $J/2$ Kronecker delta
functions, which vanish when contracted with the null vector $\xi$. The $m\neq
0$ terms are instead not rotation invariant, due to the phase factor depending
only on one momentum component, and thus are non-vanishing.

Following an idea from [32], we define a generating function for these
expectations:

$\sum_{J\geq 0}\frac{\langle\mathcal{O}_{0,J}(x,\xi)\rangle}{J!}=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot\xi_{m}}}{(p^{2})^{\zeta}}\equiv C^{(0)}(\xi)\;,$ (2.16)

where $\xi_{m}=(\chi+mL,\bm{\xi})$. We recognize that such a generating
function is the two-point function, with the $m=0$ term subtracted, and
evaluated at the null vector $\xi$, $C^{(0)}(\xi)=C_{L}(\xi)-C(\xi)$.
Therefore, it is enough to use our previous results, evaluated at $x=\xi$ and
read off the coefficients of the one-point function from the expansion in
powers of $\chi$.

Note that the general OPE formula for the two-point function in (1.10), when
evaluated at a null $\xi$ (and subtracting the divergent contribution from the
identity operator, $\Delta_{\mathcal{O}}=J=0$), becomes:

$C^{(0)}(\xi)=\sum_{J\in
2\mathbb{N}_{0}}\frac{1}{L^{2\Delta_{\phi}+J}}\frac{b^{(\text{free})}_{0,J}f_{\phi\phi
J}}{c_{0,J}}\chi^{J}\;.$ (2.17)

On the other hand, using the general formula (1.6), and contracting with our
null vector $\xi$, we find:

$\langle\mathcal{O}_{0,J}(x,\xi)\rangle=\frac{b^{(\text{free})}_{0,J}}{L^{2\Delta_{\phi}+J}}\chi^{J}\;.$
(2.18)

Therefore, by comparison to (2.16), in our conventions it must be $f_{\phi\phi
J}/c_{0,J}=1/J!$. In appendix A we prove that this is indeed the case.

After subtracting the divergent $m=0$ term (corresponding to the contribution
from the identity operator) from (2.13), and replacing $x\to\xi$, we have:

$C^{(0)}(\xi)=c(\Delta_{\phi})\sum_{m\neq
0}\frac{1}{|Lm|^{2\Delta_{\phi}}(1+\frac{2\chi}{mL})^{\Delta_{\phi}}}\;.$
(2.19)

Expanding in $\chi$, and combining the terms with $m>0$ and $m<0$, the odd
powers of $\chi$ cancel while the even ones reconstitute a Riemann zeta
function, which we denote $\zeta_{R}$, to avoid confusion with our long-range
parameter:

$C^{(0)}(\xi)=c(\Delta_{\phi})\sum_{J\in
2\mathbb{N}_{0}}\frac{\chi^{J}}{L^{J+2\Delta_{\phi}}}\frac{\Gamma(J+\Delta_{\phi})}{J!\Gamma(\Delta_{\phi})}2^{J+1}\zeta_{R}(J+2\Delta_{\phi})\;.$
(2.20)

From such expansion we can extract the one-point function coefficients:

$b^{(\text{free})}_{0,J}=c(\Delta_{\phi})\frac{\Gamma(J+\Delta_{\phi})}{\Gamma(\Delta_{\phi})}2^{J+1}\zeta_{R}(J+2\Delta_{\phi})\;.$
(2.21)

Lastly, using the relation $f_{\phi\phi J}/c_{0,J}=1/J!$ we find:

$\begin{split}a^{(\text{free})}_{0,J}&\equiv\frac{J!\Gamma(\Delta_{SR})}{2^{J}\Gamma(\Delta_{SR}+J)}\frac{b^{(\text{free})}_{0,J}f_{\phi\phi
J}}{c_{0,J}}=c(\Delta_{\phi})\frac{2\Gamma(\Delta_{\phi}+J)\Gamma(\Delta_{SR})}{\Gamma(\Delta_{\phi})\Gamma(\Delta_{SR}+J)}\zeta_{R}(J+2\Delta_{\phi})\;,\end{split}$
(2.22)

which matches with the result in [1], up to the different normalization of the
two-point function $c(\Delta_{\phi})$.

###  2.3 One-point functions of higher-twist operators

As we already discussed, in the case of long-range models the OPE of
$\phi\times\phi$ contains also higher-twist operators, of the form (2.4) with
$k>0$. As we do not have a general formula for the primaries with $k>0$, we
cannot explicitly compute their two-point functions, and the associated
coefficient, or their OPE coefficient. However, the precise distribution of
derivatives is irrelevant for the tadpole computation, as the derivatives
acting on the left and right field $\phi$ in Fourier space both reduce to
(plus or minus) the same loop momentum. Therefore, we can implicitly define
the normalization to be the one such that the one-point function reduces
to:111111In the cases of the operators that we wrote in (2.7), the explicit
normalizations are
$\mathcal{N}_{1,0}=\frac{\Delta_{\phi}}{\Delta_{SR}}\,,\quad\mathcal{N}_{2,0}=\frac{\Delta_{\phi}(\Delta_{\phi}+1)\left(\Delta_{\phi}-\Delta_{SR}\right)}{2\Delta_{\phi}\left(\Delta_{SR}^{2}+\Delta_{SR}+1\right)-\Delta_{SR}^{2}\left(2\Delta_{SR}+1\right)+\Delta_{SR}+2}\,,\quad\mathcal{N}_{1,2}=-\frac{\Delta_{\phi}\left(\Delta_{\phi}+1\right)}{2\Delta_{SR}}\,.$

$\langle\mathcal{O}_{k,J}(x,\xi)\rangle=\frac{1}{L}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\frac{(\imath\mathbf{p}\cdot\bm{\xi}+\imath
q_{n}\chi)^{J}}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}}(\mathbf{p}^{2}+q_{n}^{2})^{k}-\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath
x\cdot\xi)^{J}}{(p^{2})^{\zeta}}(p^{2})^{k}\;,$ (2.23)

which is finite, by the same argument used for $k=0$, and again the
subtraction is really needed only for $J=0$. As before, we will aim at a
generating function for the corresponding coefficient $a_{k,J}$, defined as in
(1.10), which is independent of the normalization of $\mathcal{O}_{k,J}$.

We begin by defining a generating function of such one-point functions. That
is, after using again the Poisson formula, we write:

$\begin{split}C^{(k)}(\chi,\bm{\xi})&\equiv\sum_{J\geq
0}\frac{\langle\mathcal{O}_{k,J}(x,\xi)\rangle}{J!}=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot\xi_{m}}}{(p^{2})^{\zeta}}(p^{2})^{k}\\\
&=(-\partial^{2})^{k}C^{(0)}(x)_{|_{x=(\chi,\bm{\xi})}}\;,\end{split}$ (2.24)

which unsurprisingly is the result of applying $k$ Laplacians to the
(subtracted) covariance and evaluating it at $x=\xi$. Then we can find an
alternative representation of $C^{(k)}(\chi,\bm{\xi})$ just by acting with $k$
Laplacians on (1.8) and evaluating at $x=\xi$, we find:

$C^{(k)}(\chi,\bm{\xi})=(-1)^{k}\sum_{J\in
2\mathbb{N}_{0}}\frac{\chi^{J}}{J!}\frac{a_{k,J}}{L^{2\Delta_{\phi}+2k+J}}\frac{2^{J+2k}k!\Gamma(k+J+\Delta_{SR}+1)}{(J+\Delta_{SR})\Gamma(\Delta_{SR})}\;.$
(2.25)

On the other hand, from the Fourier representation in (2.24) it is clear that
$C^{(k)}(\chi,\bm{\xi})$ corresponds to the replacement
$\zeta\rightarrow\zeta-k$ (i.e. $\Delta_{\phi}\rightarrow\Delta_{\phi}+k$) in
the subtracted covariance $C^{(0)}(\chi,\bm{\xi})$, and thus we obtain:

$C^{(k)}(\chi,\bm{\xi})=c(\Delta_{\phi}+k)\sum_{m\neq
0}\frac{1}{|\xi_{m}|^{2\Delta_{\phi}+2k}}\;,$ (2.26)

which, after repeating the same steps as before, becomes:

$\begin{split}C^{(k)}(\chi,\bm{\xi})&=c(\Delta_{\phi}+k)\sum_{J\in
2\mathbb{N}_{0}}\frac{\chi^{J}}{L^{J+2k+2\Delta_{\phi}}}\frac{\Gamma(J+\Delta_{\phi}+k)}{J!\Gamma(\Delta_{\phi}+k)}2^{J+1}\zeta_{R}(J+2\Delta_{\phi}+2k)\;.\end{split}$
(2.27)

Lastly, by comparison with (2.25), we extract the $a_{k,J}$ coefficients:

$\begin{split}a_{k,J}&=c(\Delta_{\phi})\frac{2(J+\Delta_{SR})\Gamma(\Delta_{SR})\Gamma(\Delta_{\phi}+J+k)\Gamma(\Delta_{\phi}-\Delta_{SR}+k)}{k!\Gamma(\Delta_{\phi})\Gamma(\Delta_{SR}+J+k+1)\Gamma(\Delta_{\phi}-\Delta_{SR})}\zeta_{R}(J+2k+2\Delta_{\phi})\;,\end{split}$
(2.28)

that again match with the result in [1] up to the $c(\Delta_{\phi})$ factor.

##  3 Classical long-range $O(N)$ model with a finite size

The long-range $O(N)$ model is defined by the action:

$S[\phi]=\int_{\mathbb{R}^{d}}d^{d}x\left(\frac{1}{2}\phi_{a}(x)(-\partial^{2})^{\zeta}\phi_{a}(x)+\frac{1}{2}M_{\text{bare}}^{2\zeta}\phi^{2}(x)+\frac{\lambda}{4N}(\phi^{2}(x))^{2}\right)\;,$
(3.1)

where $a=1,\ldots,N$, $\phi^{2}(x)\equiv\sum_{a}\phi_{a}(x)^{2}$, and
$M_{\text{bare}}$ is the bare mass, which we treat as a pure counterterm
because we are interested in the critical theory. It is convenient to
introduce also the intermediate field representation, by the usual Hubbard-
Stratonovich transformation:

$S[\phi,\sigma]=\int
d^{d}x\left(\frac{1}{2}\phi_{a}(x)(-\partial^{2})^{\zeta}\phi_{a}(x)+\frac{1}{2}M_{\text{bare}}^{2\zeta}\phi^{2}(x)+N\frac{\sigma(x)^{2}}{4\lambda}+\frac{\operatorname{\mathrm{i}}}{2}\,\sigma(x)\phi^{2}(x)\right)\;.$
(3.2)

For any $\zeta\in(d/4,1)$ one finds a non-trivial IR fixed point, most easily
accessible in the large-$N$ expansion, or in the small-$\epsilon$ expansion
where $\zeta=(d+\epsilon)/4$ [4]. Characteristic features of such fixed point
are the absence of anomalous dimension (that is, the two-point function scales
as in the corresponding GFFT, as proved rigorously in [14]), and
$\zeta$-dependent critical exponents. Such $\zeta$-dependent fixed points are
then associated to a one-parameter family of universality classes,
collectively called _long-range universality class_. An interesting crossover
phenomenon actually occurs before reaching $\zeta=1$: there exists a
$\zeta^{\star}\in(d/4,1)$ such that for $\zeta>\zeta^{\star}$ the operator
$\phi_{a}\partial^{2}\phi_{a}$ (i.e. the local kinetic term) becomes relevant
at the IR fixed point, thus destabilizing it and driving the system to the
short-range universality class [7]. The value of such crossover exponent for
$2\leq d<4$ is $\zeta^{\star}=1-\eta_{SR}/2$, where $\eta_{SR}$ is the
anomalous dimension of the short-range model.121212For $d=1$ the upper limit
is $\zeta^{\star}=1/2$, beyond which there can be no spontaneous symmetry
breaking [3]. However, since $\eta_{SR}=O(1/N)$ [53], at large $N$ there is no
crossover, and our analysis will apply all the way to $\zeta=1$.

In the large-$N$ limit, the propagator only receives mass corrections, via
tadpoles, and we thus have the following Schwinger-Dyson equation:

$\tilde{G}(p)^{-1}=\tilde{C}(p)^{-1}+M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\int\frac{d^{d}p^{\prime}}{(2\pi)^{d}}\tilde{G}(p^{\prime})\;,$
(3.3)

where $\tilde{C}(p)$ is as in the GFFT, and we defined
$\langle\phi_{a}(x)\phi_{a}(0)\rangle=\delta_{ab}G(x)$ and $\tilde{G}(p)$ the
Fourier transform of $G(x)$. Since the tadpole is momentum-independent, the
solution must be of the form $\tilde{G}(p)^{-1}=\tilde{C}(p)^{-1}+M^{2\zeta}$.
Tuning the model to criticality requires choosing a bare mass
$M_{\text{bare}}$ such that $M=0$, which amounts to canceling the tadpole
contribution, hence we have:

$M_{\text{bare}}^{2\zeta}=-\frac{\lambda}{2}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{1}{(p^{2})^{\zeta}}\equiv-\frac{\lambda}{2}\mathcal{T}\;.$
(3.4)

In the same limit, the connected four-point function is given by a sum of
chains of bubble diagrams, with full leading-order propagator on each edge,
see Figure 1.

Figure 1: Chain of bubbles with $n\geq 1$ vertices. The connected four-point
function is given by the sum of all such chains.

Therefore, the proper vertex in momentum space reads:

$\Gamma^{(4)}(p)=\frac{2\lambda}{1+\lambda B(p)}\;,$ (3.5)

where $p$ is the total transferred momentum and we defined the bubble
integral:

$B(p)=\int\frac{d^{d}q}{(2\pi)^{d}}\frac{1}{q^{2\zeta}(p-q)^{2\zeta}}=\frac{c(d/2-\zeta)^{2}}{c(d-2\zeta)}|p|^{d-4\zeta}\;.$
(3.6)

We define the renormalized dimensionless coupling $g$ at some renormalization
group scale $|p|=\mu$:

$g=\mu^{4\Delta_{\phi}-d}\Gamma^{(4)}(\mu)/2\;.$ (3.7)

It follows that the beta function is:131313Notice that for
$d-4\Delta_{\phi}=4\zeta-d=\epsilon$ (i.e. for $\zeta=(d+\epsilon)/4$) the
second term is finite (and positive) in the limit $\epsilon\to 0$, because
$c(2\Delta_{\phi})\sim 1/\Gamma(\epsilon/2)$, and we recover the standard one-
loop result $\beta(g)=-\epsilon g+2g^{2}/((4\pi)^{d/2}\Gamma(d/2))$ (see [6]).

$\beta(g)=\mu\frac{dg}{d\mu}=-(d-4\Delta_{\phi})g+\frac{c(\Delta_{\phi})^{2}}{c(2\Delta_{\phi})}(d-4\Delta_{\phi})g^{2}\;.$
(3.8)

Therefore, for $\Delta_{\phi}<d/4$ (i.e. $\zeta>d/4$) the beta function has a
zero (i.e. a fixed point) at
$g^{\star}=c(2\Delta_{\phi})/c(\Delta_{\phi})^{2}$. It is worth noticing that
since we can trivially invert the relation between bare and renormalized
coupling to get:

$\mu^{4\Delta_{\phi}-d}\lambda=\frac{g}{1-\frac{c(\Delta_{\phi})^{2}}{c(2\Delta_{\phi})}g}\;,$
(3.9)

the fixed point corresponds to the limit $\lambda\to\infty$.

Another well known large-$N$ result for the short-range $O(N)$ model, which
generalizes straightforwardly to the long-range case, is the fact that at the
IR fixed point the composite operator $\phi^{2}(x)$ (of dimension
$\Delta_{0,0}=2\Delta_{\phi}=d-2\zeta$ in the GFFT) is mapped to its shadow
operator (of dimension $\Delta_{\sigma}=d-\Delta_{0,0}=2\zeta$), which can be
identified with the intermediate field $\sigma(x)$. This is a particular
realization of the general result by Gubser and Klebanov [54] on large-$N$
flows driven by a double-trace operator (in this case, the $(\phi^{2})^{2}$
interaction). There are many other ways to see the same result, for example by
means of the conformal partial wave representation of the four-point function
[43], or diagrammatically, as we will see below.

The other bilinear primaries (2.4) remain instead in the spectrum. In
addition, in the interacting model we also have other operators, with more
powers of $\sigma$ or $\phi$. Naively, at large $N$ such operators are
suppressed in the OPE of $\phi\times\phi$. However, the situation is more
subtle: as argued in [1], in the expansion (1.10) it can happen that for
certain $\mathcal{O}$’s while $f_{\phi\phi\mathcal{O}}/\sqrt{c_{\mathcal{O}}}$
is suppressed, $b_{\mathcal{O}}/\sqrt{c_{\mathcal{O}}}$ is boosted, thus
resulting in a finite contribution to $a_{\mathcal{O}}$ even at large $N$.
This in fact does happen for the operators $\sigma^{n}(x)$ which yield non-
trivial contributions to the OPE, as we will discuss in 3.5.

#### Two-point functions of spinning operators.

When evaluating at leading order in $1/N$ the two-point functions of bilinear
operators of the $O(N)$ model, we need to sum chains of bubble diagrams, as
for the four-point function. We now show that for the spinning operators at
the fixed point the sum over the chain of bubbles reduces to the contribution
of a single bubble, and thus their two-point functions (and hence scaling
dimensions) are the same as in the GFFT.

First, we check once again that the $\phi^{2}$ operator drops out of the CFT.
The sum of bubbles becomes a simple geometric sum in Fourier space, hence we
have:

$\langle[\phi_{a}\phi_{a}]_{0}(x_{1})[\phi_{a}\phi_{a}]_{0}(x_{2})\rangle=\int\frac{d^{d}p}{(2\pi)^{d}}e^{\imath
p\cdot x_{12}}\frac{B(p)}{1+\lambda B(p)}\;.$ (3.10)

For $\lambda\to\infty$ the result vanishes as expected. The two-point function
of $\sigma(x)$ is similar, but with the replacement $B(p)\to\lambda$ in the
numerator (it coincides with the $\Gamma^{(4)}(p)$ above), hence it remains
finite in the limit. We infer that $\phi^{2}$ is replaced by its shadow
operator $\sigma$ at the interacting fixed point.

For higher-spin operators we have instead to distinguish the first and last
bubbles:

$\langle[\phi_{a}\phi_{a}]_{0,\mu_{1}\cdots\mu_{J}}(x_{1})[\phi_{a}\phi_{a}]_{0,\nu_{1}\cdots\nu_{J}}(x_{2})\rangle=\int\frac{d^{d}p}{(2\pi)^{d}}e^{\imath
p\cdot
x_{12}}\left(B_{\mu_{1}\cdots\mu_{J},\nu_{1}\cdots\nu_{J}}-B_{\mu_{1}\cdots\mu_{J}}\frac{\lambda}{1+\lambda
B(p)}B_{\nu_{1}\cdots\nu_{J}}\right)\;,$ (3.11)

where $B_{\mu_{1}\cdots\mu_{J},\nu_{1}\cdots\nu_{J}}$ is the Fourier transform
of the same two-point function in the free theory, i.e. the diagram without
internal vertices, while $B_{\mu_{1}\cdots\mu_{J}}$ is the Fourier transform
of the two-point function of the spin-$J$ operator with the spin-0 one, also
in the free theory. However, the latter vanishes: technically because of
orthogonality of the Gegenbauer polynomials [50], and conceptually because in
a CFT the two-point functions of primaries of different dimensions must
vanish.

The explicit computation of these two-point functions is reported in appendix
A.

###  3.1 Schwinger-Dyson equation in the hyper-strip geometry

We now place the model in a hyper-strip geometry with one finite dimension of
length $L$ and periodic boundary conditions. The action is as in (3.1) with
the replacement $\mathbb{R}^{d}\to S^{1}_{L}\times\mathbb{R}^{d-1}$. Such
setting has been extensively studied in the short-range case, and in
particular in the large-$N$ limit (e.g. [55, 56, 28, 29, 1, 30, 32]).

Like in the short-range model, also in the long-range case the propagator will
acquire a mass due to the introduction of the compact direction, as we now
show by reducing the Schwinger-Dyson equation to a mass gap equation.

At large $N$, the Schwinger-Dyson equation in momentum space is:

$\tilde{G}(q_{n_{0}},\mathbf{p}_{0})^{-1}=\tilde{C}_{L}(q_{n_{0}},\mathbf{p}_{0})^{-1}+M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\frac{1}{L}\sum_{n}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\tilde{G}(q_{n},\mathbf{p})\;.$
(3.12)

The inverse free propagator is
$\tilde{C}_{L}(q_{n},\mathbf{p})^{-1}=(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}$ as
in the GFFT, $M_{\text{bare}}^{2\zeta}$ is the bare mass (3.4), and
$q_{n}=2\pi n/L$, as before. As the tadpole term in (3.12) does not depend on
the external momentum, the full inverse two-point function is
$\tilde{G}(q_{n},\mathbf{p}_{0})^{-1}=(\mathbf{p}_{0}^{2}+q_{n}^{2})^{\zeta}+M_{L}^{2\zeta}$
for some $M_{L}$ to be determined self-consistently.

The self-consistent gap equation for the finite-size mass is:

$\begin{split}M_{L}^{2\zeta}&=M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\frac{1}{L}\sum_{n}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\frac{1}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}+M_{L}^{2\zeta}}\\\
&=\frac{\lambda}{2}\left(-M_{L}^{2\zeta}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{1}{((p^{2})^{\zeta}+M_{L}^{2\zeta})(p^{2})^{\zeta}}+\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
Lmq}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}\right)\\\
&\equiv\frac{\lambda}{2}\,\left(\mathcal{T}(M_{L})-\mathcal{T}\right)\;,\end{split}$
(3.13)

where in the second equality we used the Poisson summation formula (2.10), we
have introduced again a $d$-dimensional momentum $p=(q,\mathbf{p})$, and we
have combined the $m=0$ term with the explicit expression of the bare mass. In
the last step we have defined the subtracted tadpole
$\mathcal{T}(M_{L})-\mathcal{T}$.

The first integral in the subtracted tadpole can be computed explicitly and we
find:

$\mathcal{T}(M_{L})-\mathcal{T}=\frac{\pi
M_{L}^{d-2\zeta}}{(4\pi)^{d/2}\Gamma(d/2)\zeta\sin\left(\frac{d\pi}{2\zeta}\right)}+\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
Lmq}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}\;.$ (3.14)

The second term is not manifestly UV convergent, but its convergence can be
inferred from the convergence of the $M_{L}=0$ case of the previous section
(see (2.13) without the $m=0$ term, evaluated at $x_{m}=(Lm,\mathbf{0})$,
giving $b^{(\text{free})}_{0,0}$), as the inclusion of the mass does not make
the UV behavior worse.141414We can also explicitly show its convergence
deforming the contour of integration of $q$ to a Hankel contour in the upper
(lower) complex plane for positive (negative) $m$. Its explicit evaluation is
however rather complicated, and we will take a different route to the gap
equation below.

As we recalled above, at the IR fixed point the bare coupling $\lambda$ flows
to infinity, hence the mass gap equation corresponds to finding the value of
$M_{L}$ such that the subtracted tadpole (3.14) vanishes. We emphasize that
the subtraction is essential to the existence of a solution, as otherwise the
tadpole would seem positive, which of course is naive as it is in fact
divergent: for $d/4<\zeta<1$ and $2<d<4$ (and for
$1/4<\zeta<\zeta^{\star}=1/2$ if $d=1$) the sine in the first term is
negative, and thus a cancellation between the two terms can occur.

We will come back to this equation below, from a different perspective, and
prove that it admits a solution.

###  3.2 Long-range two-point function with mass

The two-point function in direct space reads:

$G(y,\mathbf{x})=\frac{1}{L}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\frac{e^{\imath\mathbf{p}\cdot\mathbf{x}+\imath
yq_{n}}}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}+M_{L}^{2\zeta}}\;.$ (3.15)

As before, we introduce $1=\int dq\delta(q-q_{n})$ and we use the Poisson
summation formula to rewrite it as:

$G(y,\mathbf{x})=\sum_{m=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot x_{m}}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}\;,$ (3.16)

where $p$ is a $d$-dimensional momentum and $x_{m}=(y+mL,\mathbf{x})$.

In order to proceed further, we rewrite again the momentum-space propagator
$\tilde{G}(p)=\left((p^{2})^{\zeta}+M_{L}^{2\zeta}\right)^{-1}$ as a Stieltjes
transform. While such a rewriting was trivial in the massless case because the
spectral density had to be scaling invariant, it is slightly more involved in
the presence of a mass. First, as the momentum-space propagator is only a
function of the squared momentum, we can define $g(p^{2})=\tilde{G}(p)$, and
we use Cauchy’s theorem to write:

$g(p^{2})=\frac{1}{2\pi\imath}\oint_{\gamma}\frac{g(z)}{z-p^{2}}dz\;,$ (3.17)

with $\gamma$ a counterclockwise contour encircling $z=p^{2}$. Then, we deform
the contour to infinity, and in the process we pick up poles and branch cuts
of $g(z)$. Now, for $\zeta\in(0,1)$, the function
$g(z)=1/(z^{\zeta}+M_{L}^{2\zeta})$ has a branch cut on the negative axis, and
no poles in the principal sheet.151515The poles are at
$z=M_{L}^{2}e^{\frac{\imath\pi(1+2k)}{\zeta}}$, for integer $k$, and it can be
seen that for $0<\zeta<1$ none of them is in the principal sheet, which by
definition corresponds to ${\rm Arg}(z)\in(-\pi,\pi)$. Therefore, the
deformation of the contour $\gamma$ leads to a Hankel contour around the cut,
which in turn leads to the standard formula for the inversion of the Stieltjes
transform
$\rho(s)=\frac{1}{2\pi\imath}\left(g(-s-\imath\epsilon)-g(-s+\imath\epsilon)\right)$,
and we find:

$\tilde{G}(p)=\frac{1}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}=\int_{0}^{+\infty}ds\,\frac{\rho(s)}{p^{2}+s}\;,\qquad\rho(s)=\frac{\sin(\pi\zeta)}{\pi}\frac{s^{\zeta}}{s^{2\zeta}+2M_{L}^{2\zeta}s^{\zeta}\cos(\pi\zeta)+M_{L}^{4\zeta}}\;.$
(3.18)

Going back to position space, and setting $s=t^{2}$, we have:

$G(y,\mathbf{x})=\frac{2\sin(\pi\zeta)}{(2\pi)^{d/2}\pi}\sum_{m=-\infty}^{+\infty}\int_{0}^{+\infty}dt\;\left(\frac{t}{|x_{m}|}\right)^{d/2-1}K_{\frac{d}{2}-1}(t|x_{m}|)\frac{t^{1+2\zeta}}{t^{4\zeta}+2M_{L}^{2\zeta}t^{2\zeta}\cos(\pi\zeta)+M_{L}^{4\zeta}}\;.$
(3.19)

Because of the mass, in this case we cannot factor out of the integral the
$x$-dependence.

At this step one might worry that for $\zeta\to 1$ everything vanishes because
of the $\sin(\pi\zeta)$ factor. This is not the case because in this limit the
denominator in the integrand becomes $(t^{2}-M_{L}^{2})^{2}$ and the integral
develops a singularity (this is the analogue of the $\Gamma(1-\zeta)$ factor
in the massless case of (2.13)). In fact in the $\zeta\to 1$ limit we recover
the short-range case (2.13), as it can be shown by using a variation of the
Plemelj–Sokhotski formula.

For later reference we write explicitly the case $d=3$, in which the Bessel
function takes a simple form:

$G(y,x)=\frac{\sin(\pi\zeta)}{2\pi^{2}}\sum_{m=-\infty}^{+\infty}\int_{0}^{+\infty}dt\;\frac{e^{-t|x_{m}|}}{|x_{m}|}\frac{t^{1+2\zeta}}{t^{4\zeta}+2M_{L}^{2\zeta}t^{2\zeta}\cos(\pi\zeta)+M_{L}^{4\zeta}}\;.$
(3.20)

###  3.3 Finite-size mass solution

The $m\neq 0$ part of the subtracted tadpole (3.14) is the $m\neq 0$ part of
the two-point function, evaluated at $x=0$. Therefore, we can use the
representation (3.20) to rewrite the mass gap equation in $d=3$ as:

$\frac{LM_{L}}{4\pi\zeta\sin\left(\frac{3\pi}{2\zeta}\right)}-\frac{\sin(\pi\zeta)}{\pi^{2}}\int_{0}^{+\infty}\frac{t^{2\zeta+1}}{t^{4\zeta}+2\cos(\pi\zeta)t^{2\zeta}+1}\ln\left(1-e^{-LM_{L}t}\right)=0\;,$
(3.21)

which is more suitable for a numerical evaluation, having eliminated the sum
in favor of a single integral that is manifestly convergent. We find that the
equation always admits a real positive solution, which we plot as a function
of $\zeta$ in Figure 2. The solution goes to zero at $\zeta=3/4$, where the IR
fixed point merges with the Gaussian one, and it goes to the short-range value
(e.g. [56]) for $\zeta\to 1$.

Figure 2: Numerical solutions for the mass gap equation as a function of
$\zeta$, at $d=3$. The red dot is the value of the short-range thermal mass
and it is approached in the $\zeta\rightarrow 1$ limit as expected.

###  3.4 One-point functions of the minimal-twist operators

The one-point functions of the minimal-twist operators in the hyper-strip
geometry are defined similarly to the GFFT case, except for the presence of
the mass and the restriction to $J>0$, as the $\phi^{2}$ operator is not in
the spectrum:

$\begin{split}\langle\mathcal{O}_{0,J}(x,\xi)\rangle&=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath\mathbf{p}\cdot\bm{\xi}+\imath
q\chi)^{J}\,e^{\imath
Lmq}}{(\mathbf{p}^{2}+q^{2})^{\zeta}+M_{L}^{2\zeta}}\;.\end{split}$ (3.22)

Notice that the restriction to $m\neq 0$ comes from the fact that the $m=0$
term vanishes in any rotationally invariant regularization, due to the
contraction with $\xi$, and as before the other terms do not need any
subtraction.

As in the previous section, following [32], we can evaluate the subtracted
two-point function at a null vector $\xi=(\chi,\bm{\xi})$ and interpret it as
the generating function of the one-point functions:

$G^{(0)}(\chi,\bm{\xi})=\sum_{J>0}\frac{\langle\mathcal{O}_{0,J}(x,\xi)\rangle}{J!}=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot\xi_{m}}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}-M_{L}^{2\zeta}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{1}{((p^{2})^{\zeta}+M_{L}^{2\zeta})(p^{2})^{\zeta}}\;.$
(3.23)

Notice that in order to get an exponential in the last step we had to include
a spurious $J=0$ term, which however we cancel with the last term, thanks to
the mass gap equation imposing the vanishing of the subtracted tadpole (3.14).

An alternative explanation for the $J=0$ subtraction in (3.23) is obtained by
considering the general expansion of the two-point function in (1.10). As we
know, at large $N$ the spectrum of the $\phi\times\phi$ OPE contains the
identity operator (with $\Delta_{\rm id}=0$), the intermediate field (with
$\Delta_{\sigma}=2\zeta$), and the bilinear operators with
$\Delta_{k,J}=2\Delta_{\phi}+J+2k$, the $J=k=0$ operator being
excluded.161616As discussed earlier in the text (and originally pointed out in
[1]), one-point functions can have a positive scaling in $N$ that boosts
operators whose OPE is suppressed in $1/N$, thus introducing new operators in
(1.10) even at large $N$. However, such operators have larger scaling
dimension (at fixed spin) than those we are considering here, so they do not
affect this argument.  When setting $x=\xi$, since $|\xi|=0$ and the
Gegenbauer polynomial is of order $J$, all the bilinear operators with $k>0$
drop from the expansion, and so does also $\sigma$, since
$2\zeta>2\Delta_{\phi}$ for $\zeta>d/4$. Therefore, we should find only the
$\Delta_{0,J>0}$ contributions plus a divergent one from the identity
operator. The latter corresponds to the massless propagator in the non-compact
case, but the $m=0$ term we have removed in the first term of (3.23) is the
massive propagator. In order to get the right operator content of (1.10), we
should thus add back this massive propagator, and then subtract the massless
one, which is precisely what the last term in (3.23) does, since:

$\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath p\cdot
x}}{((p^{2})^{\zeta}+M_{L}^{2\zeta}}-\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot
x}}{(p^{2})^{\zeta}}=-M_{L}^{2\zeta}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot x}}{((p^{2})^{\zeta}+M_{L}^{2\zeta})(p^{2})^{\zeta}}\;,$ (3.24)

is a function of $|x|$, convergent at $|x|=0$, and thus its value at $x=\xi$
is the same as at $x=0$.

Subsequently we can extract the $b_{k,J}$ coefficients from the formal series
expansion in $\chi$. In $d=3$, we have:

$G^{(0)}(\chi,\bm{\xi})=\frac{\pi
M_{L}^{3-2\zeta}}{4\pi\zeta\sin\left(\frac{3\pi}{2\zeta}\right)}+\frac{\sin(\pi\zeta)}{2\pi^{2}}\sum_{m\neq
0}\int_{0}^{+\infty}dt\;\frac{e^{-t|\xi_{m}|}}{|\xi_{m}|}\frac{t^{1+2\zeta}}{t^{4\zeta}+2M_{L}^{2\zeta}t^{2\zeta}\cos(\pi\zeta)+M_{L}^{4\zeta}}\;,$
(3.25)

and the expansion in $\chi$ inside the integral takes the same form as in
[32], but with $t$ replacing the mass of the short-range case. We thus get an
integral representation for the one-point functions of the $J>0$ spinning
operators:171717Again we can use $f_{\phi\phi J}/c_{0,J}=1/J!$ and easily
compute $a_{0,J}$ from $b_{0,J}$

$b_{0,J}=\frac{\sin(\pi\zeta)}{2\pi^{2}}\sum_{n=0}^{J}\frac{2^{n-J+1}(2J-n)!}{n!(J-n)!}\int_{0}^{+\infty}dt\
\text{Li}_{J+1-n}\left(e^{-tLM_{L}}\right)\frac{t^{1+2\zeta+n}(LM_{L})^{n+2-2\zeta}}{t^{4\zeta}+2t^{2\zeta}\cos(\pi\zeta)+1}\;.$
(3.26)

Unfortunately the integral above can not be evaluated in terms of usual
special functions and needs to be evaluated numerically. We compare results to
the short-range case by means of a plot in Figure 3.

Figure 3: Plot of the ratio between $b_{0,J}$ for $\zeta=0.8$ and the
corresponding short-range value.

###  3.5 One-point functions of higher-twist operators

As in the GFFT, we define the one-point functions of the higher-twist bilinear
operators with an implicit choice of normalization such that:

$\langle\mathcal{O}_{k,J}(x,\xi)\rangle=\frac{1}{L}\sum_{n}\int\frac{d^{d-1}\mathbf{p}}{(2\pi)^{d-1}}\frac{(\imath\mathbf{p}\cdot\bm{\xi}+\imath
q_{n}\chi)^{J}}{(\mathbf{p}^{2}+q_{n}^{2})^{\zeta}+M_{L}^{2\zeta}}(\mathbf{p}^{2}+q_{n}^{2})^{k}-\text{counterterms}\;.$
(3.27)

In order to make explicit the structure of divergences and related
counterterms we make use again of the Poisson formula, and we write:

$\langle\mathcal{O}_{k,J}(x,\xi)\rangle=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(\imath
p\cdot\xi)^{J}\,(p^{2})^{k}\,e^{\imath
Lmq}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}+\delta_{J,0}\mathcal{C}_{k}\;,$ (3.28)

where:

$\mathcal{C}_{k}\equiv\int\frac{d^{d}p}{(2\pi)^{d}}\frac{(p^{2})^{k}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}-\sum_{n=0}^{\lfloor\frac{d+2k-2\zeta}{2\zeta}\rfloor}(-M_{L}^{2\zeta})^{n}\int\frac{d^{d}p}{(2\pi)^{d}}(p^{2})^{k-(n+1)\zeta}\;.$
(3.29)

We have used again the fact that for $J>0$ the $m=0$ term vanishes in any
rotationally invariant regularization. For $J=0$ we have subtracted instead
the purely divergent parts of the $m=0$ term. Alternatively, we could drop
such explicit subtractions and evaluate the $m=0$ term by analytic
continuation in $k$ from a region of convergence (i.e. $-d/2<k<-d/2+\zeta$).
It is however interesting, and actually useful, to keep them explicit as we
will see that they have a direct interpretation in terms of (1.10).

We are going to define again a generating function for these one-point
functions, as:

$\begin{split}G^{(k)}(\chi,\bm{\xi})&\equiv\sum_{J\geq
0}\frac{\langle\mathcal{O}_{k,J}(x,\xi)\rangle}{J!}=\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{\imath
p\cdot\xi_{m}}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}(p^{2})^{k}+\mathcal{C}_{k}\\\
&=\left((-\partial^{2})^{k}G(x)-\mathcal{L}_{k}(x)\right)_{x=\xi}\;,\end{split}$
(3.30)

where:

$\begin{split}\mathcal{L}_{k}(x)&=\sum_{n=0}^{\lfloor\frac{d+2k-2\zeta}{2\zeta}\rfloor}(-M_{L}^{2\zeta})^{n}\int\frac{d^{d}p}{(2\pi)^{d}}(p^{2})^{k-(n+1)\zeta}\,e^{\imath
p\cdot x}\\\
&=\sum_{n=0}^{\lfloor\frac{d+2k-2\zeta}{2\zeta}\rfloor}\frac{(-M_{L}^{2\zeta})^{n}\,c(d/2+k-(n+1)\zeta)}{|x|^{d+2k-2(n+1)\zeta}}\;.\end{split}$
(3.31)

The interpretation of $\mathcal{L}_{k}(x)$ is the following. When acting $k$
times with a Laplacian on (1.10) and setting $x=\xi$, we have a similar
structure as in the GFFT from the bilinear operators and the identity
operator. In addition, we have other operators that at fixed $k$ give
vanishing contributions (see footnote 16), but we also have operators that
give new divergent contributions from some $k$ on: these are operators
$\mathcal{O}$ with $\Delta_{\mathcal{O}}-2\Delta_{\phi}<2k$ not even. For
example, the intermediate field $\sigma$, having $J=0$ and non-integer
dimension $\Delta_{\sigma}=2\zeta$, gives vanishing contribution at $k=0$, but
a diverging one as soon as $k\geq 1$ if $d\geq 2$: indeed
$\partial^{2}|x|^{\Delta_{\sigma}-2\Delta_{\phi}}=\partial^{2}|x|^{4\zeta-d}=(4\zeta-d)(4\zeta-2)|x|^{4\zeta-d-2}$,
and $4\zeta-d-2<0$ for $\zeta<1$ and $d\geq 2$. More in general, we see by
comparing (3.31) with (1.10), that the operators leading to singularities in
$(-\partial^{2})^{k}G(\xi)$ have dimension $\Delta_{\sigma^{n}}=2n\zeta$, and
thus correspond to the composite fields $\sigma^{n}$:181818It can be shown
that spinning operators involving $\sigma$ or higher powers of $\phi^{2}$ are
subleading in the $1/N$ expansion (remember that we compute $b_{\mathcal{O}}$
at finite $L$, but $f_{\phi\phi\mathcal{O}}$ and $c_{\mathcal{O}}$ at infinite
$L$). This explains why we only find contributions from $\sigma^{n}$, which
has $a_{\sigma^{n}}\sim O(N^{0})$.

$(-\partial^{2})^{k}\frac{a_{\sigma^{n}}}{L^{2n\zeta}}|x|^{2n\zeta-2\Delta_{\phi}}=\frac{a_{\sigma^{n}}}{L^{2n\zeta}}(-4)^{k}\frac{\Gamma(n\zeta-\Delta_{\phi}+1)\Gamma(d/2+n\zeta-\Delta_{\phi})}{\Gamma(n\zeta-\Delta_{\phi}+1-k)\Gamma(d/2+n\zeta-\Delta_{\phi}-k)}|x|^{2n\zeta-2\Delta_{\phi}-2k}\;.$
(3.32)

Therefore, we need to subtract these dangerous terms from
$(-\partial^{2})^{k}G(x)$ before setting $x=\xi$, which is what we are doing
with $\mathcal{L}_{k}(x)$. Equivalently, we can introduce a regularization,
subtract $\mathcal{L}_{k}^{\rm reg.}(\xi)$, and then remove the regulator,
which is what we are doing with the subtractions of $\mathcal{C}_{k}$.

The above discussion has an important consequence: evaluating
$(-\partial^{2})^{k}G(x)$ allows us to extract not only the $a_{k,J}$
coefficients (from the finite part at $x=\xi$, as we did in the GFFT), but
also the $a_{\sigma^{m}}$ coefficients (from the divergent part). The latter
is obtained straightforwardly by comparing (3.31) and (3.32):

$\begin{split}a_{\sigma^{n}}&=(-1)^{n+k}(LM_{L})^{2n\zeta}\,\frac{c(\Delta_{\phi}+k-n\zeta)}{4^{k}}\frac{\Gamma(n\zeta-\Delta_{\phi}+1-k)\Gamma(d/2+n\zeta-\Delta_{\phi}-k)}{\Gamma(n\zeta-\Delta_{\phi}+1)\Gamma(d/2+n\zeta-\Delta_{\phi})}\\\
&=(-1)^{n}\frac{(LM_{L})^{2n\zeta}}{4^{(n+1)\zeta}\pi^{d/2}}\,\frac{\pi}{\sin(\pi(\Delta_{\phi}-n\zeta))\Gamma(n\zeta-\Delta_{\phi}+1)\Gamma(d/2+n\zeta-\Delta_{\phi})}\;.\end{split}$
(3.33)

Notice that any dependence on $k$ has dropped out, and the result is valid for
all $n\geq 1$. In the limit $\zeta\to 1$ this yields the short-range result:

$a_{\sigma^{n}}^{SR}=(LM_{L})^{2n}\,c(\Delta_{SR})\frac{\Gamma(1-\Delta_{SR})}{4^{n}\Gamma(n+1)\Gamma(n-\Delta_{SR}+1)}\;,$
(3.34)

which in $d=3$ matches the result of [1]. However, notice that while in the
short-range case in integer dimensions there are degeneracies (essentially all
operators have integer canonical dimension), this is not the case for general
$\zeta<1$, and therefore we can safely identify (3.33) with the one-point
function coefficient of $\sigma^{n}$.

Finally, we are ready to set $x=\xi$ in (3.30) and extract the $a_{k,J}$
coefficients. This can be easily done as we did in section (2.3), that is, by
comparing the Taylor expansion of (3.30) in $\chi$ with equation (2.25) and
matching order by order the coefficients of the powers of $\chi$. The reason
is that equation (2.25) was derived starting from (1.10) assuming that the
dimension of operators appearing in the OPE of $\phi\times\phi$ are
$\Delta_{k,J}=2\Delta_{\phi}+2k+J$, and this is also the case for the
interacting large-$N$ theory up to the $\sigma^{n}$ operators that we
subtracted explicitly in (3.30). To this purpose it is convenient to write
$\frac{(p^{2})^{k}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}$ as a Stieltjes transform:

$\frac{p^{2k}}{(p^{2})^{\zeta}+M_{L}^{2\zeta}}=\int_{0}^{+\infty}ds\,\frac{\rho(s)}{p^{2}+s}\;,\qquad\rho(s)=\frac{\sin(\pi\zeta)}{\pi}\frac{s^{\zeta+k}}{s^{2\zeta}+2M_{L}^{2\zeta}s^{\zeta}\cos(\pi\zeta)+M_{L}^{4\zeta}}\
.$ (3.35)

Then specializing to $d=3$ we find:

$G^{(k)}(\chi,\bm{\xi})=\frac{\sin(\pi\zeta)}{2\pi^{2}}\sum_{m\neq
0}\int_{0}^{+\infty}dt\;\frac{e^{-t|x_{m}|}}{|x_{m}|}\frac{t^{1+2\zeta+2k}}{t^{4\zeta}+2M_{L}^{2\zeta}t^{2\zeta}\cos(\pi\zeta)+M_{L}^{4\zeta}}+\mathcal{C}_{k}\
,$ (3.36)

where $\mathcal{C}_{k}$ contributes to the $J=0$ term only, while the sum over
$m$ will contribute to any $J$. After performing the small $\chi$ expansion we
get:

$\begin{split}G^{(k)}(\chi,\bm{\xi})=&-\frac{1}{L^{3-2\zeta+2k}}\left(\frac{\sin(\pi\zeta)}{\pi^{2}}\int_{0}^{+\infty}dt\ln\left(1-e^{-LM_{L}t}\right)\frac{t^{1+2\zeta+2k}(LM_{L})^{2-2\zeta+2k}}{t^{4\zeta}+2t^{2\zeta}\cos(\pi\zeta)+1}\right)+\mathcal{C}_{k}\\\
&+\sum_{J>0}\frac{\chi^{J}}{J!L^{3-2\zeta+2k+J}}\left(\frac{\sin(\pi\zeta)}{2\pi^{2}}\sum_{n=0}^{J}\frac{2^{n-J+1}(2J-n)!}{n!(J-n)!}\right.\crcr&\qquad\qquad\left.\times\int_{0}^{+\infty}dt\
\text{Li}_{J+1-n}\left(e^{-tLM_{L}}\right)\frac{t^{1+2\zeta+n+2k}(LM_{L})^{n+2-2\zeta+2k}}{t^{4\zeta}+2t^{2\zeta}\cos(\pi\zeta)+1}\right)\;,\end{split}$
(3.37)

where the first line corresponds to spin-zero operator and the sum to all the
spinning operators. Comparing with (2.25) we finally find:

$a_{k,0}=\frac{-\sin(\pi\zeta)}{(2k+1)!\pi^{2}}\int_{0}^{+\infty}dt\ln\left(1-e^{-tLM_{L}}\right)\frac{t^{1+2\zeta+2k}(LM_{L})^{2-2\zeta+2k}}{t^{4\zeta}+2t^{2\zeta}\cos(\pi\zeta)+1}+\frac{(LM_{L})^{3-2\zeta+2k}}{(2k+1)!4\pi\zeta\sin\left(\pi\frac{2k+3}{2\zeta}\right)}\;,$
(3.38)
$\begin{split}a_{k,J}=&\frac{4^{-k-J}(J+\tfrac{1}{2})\sin(\pi\zeta)}{\pi^{3/2}\Gamma(J+k+\frac{3}{2})k!}\sum_{n=0}^{J}\frac{2^{n}(2J-n)!}{n!(J-n)!}\\\
&\qquad\qquad\times\int_{0}^{+\infty}dt\
\text{Li}_{J+1-n}\left(e^{-tLM_{L}}\right)\frac{t^{1+2\zeta+n+2k}(LM_{L})^{n+2-2\zeta+2k}}{t^{4\zeta}+2t^{2\zeta}\cos(\pi\zeta)+1}\;,\end{split}$
(3.39)

where we used the fact that $\mathcal{C}_{k}$ in $d=3$ evaluates to:

$\mathcal{C}_{k}=\frac{M_{L}^{3-2\zeta+2k}}{4\pi\zeta\sin\left(\pi\frac{2k+3}{2\zeta}\right)}\
.$ (3.40)

##  4 Quantum long-range $O(N)$ model at finite temperature

The simplest example of an interacting fractional Lifshitz field theory
(FLFT), of the type discussed in the introduction is described by the action:

$\begin{split}S[\phi]=\int d\tau\int
d^{d}x\bigg{(}\frac{1}{2}&\phi_{a}(\tau,x)(-\partial_{\tau}^{2})\phi_{a}(\tau,x)+\frac{1}{2}\phi_{a}(\tau,x)(-\partial_{i}\partial_{i})^{\zeta}\phi_{a}(\tau,x)\crcr&\qquad+\frac{1}{2}M_{\text{bare}}^{2\zeta}\phi^{2}(\tau,x)+\frac{\lambda}{4N}(\phi^{2}(\tau,x))^{2}\bigg{)}\;,\end{split}$
(4.1)

where $\partial_{\tau}$ is the ordinary derivative with respect to the
Euclidean time $\tau$, $(-\partial_{i}\partial_{i})^{\zeta}$ is the fractional
Laplacian for the $d$ spatial directions and, as before, $a=1,\ldots,N$,
$\phi^{2}(\tau,x)\equiv\sum_{a}\phi_{a}(\tau,x)^{2}$, and $M_{\text{bare}}$ is
the bare mass. As argued in the introduction, we can view this model, when
defined on $\mathbb{R}^{d+1}$ or $S_{\beta}^{1}\times\mathbb{R}^{d}$, as
describing the quantum long-range $O(N)$ model at zero or finite temperature,
respectively. With the same perspective, a similar model has been considered
in [21, 22], where it has been studied at finite $N$ by other means.
Alternatively, one could view it as a classical anisotropic model with short-
range interactions in one direction (compact or not) and long-range
interactions in the remaining (non-compact) ones (see also [57] for a previous
study of anisotropic long-range models).

We notice first of all that dimensional analysis of the quadratic part of the
action implies that we have the mass dimensions $[x]=-1$, $[\tau]=-\zeta$, and
therefore $[\phi]=\Delta_{\phi}$ given by:

$2\Delta_{\phi}+2\zeta-\zeta-d=0\quad\Rightarrow\quad\Delta_{\phi}=\frac{d-\zeta}{2}\;.$
(4.2)

As a consequence, the quartic interaction in (4.1) is marginal (i.e.
$[\lambda]=0$) for $4\Delta_{\phi}=d+\zeta$, that is, $\zeta=d/3$. The model
we are interested in has $\zeta<1$ and a relevant quartic interaction,
therefore we will concentrate on $d=2$ and $2/3<\zeta<1$, but for now we keep
$d$ and $\zeta$ unspecified.

The action (4.1) is invariant under global $O(N)$ transformations with $\phi$
transforming in the fundamental representation. Moreover, it is clearly
translation invariant, both in time and space, while for $\zeta\neq 1$ it has
rotation invariance only in space. We will be interested in the limiting cases
in which it has also a scaling symmetry, namely, the ultraviolet Gaussian
fixed point and an infrared interacting fixed point. Due to the anisotropic
nature of the model we expect a Lifshitz type of scaling invariance, that is:

$\tau\to\Omega^{\zeta}\tau,\qquad x\to\Omega
x,\qquad\phi(\tau,x)\to\phi(\Omega^{\zeta}\tau,\Omega
x)=\Omega^{-\Delta_{\phi}}\phi(\tau,x)\;,$ (4.3)

for any constant $\Omega>0$. This symmetry is in general not promoted to a
generalized form of conformal invariance.191919In the case $\zeta=2$, when
formulated as a dynamic critical model (or stochastic partial differential
equation problem), it has been shown in [58] that there is a larger symmetry,
namely the Schrödinger group (see [27] and references therein), but only for
the response function, not the correlator. This result has been generalized in
[59] to the case $\zeta=2/n$, for $n\in\mathbb{N}$, but this does not include
our interesting $d=2$ case with $\zeta\in(2/3,1)$. The anisotropy and lack of
conformal invariance imply that even at fixed points the correlators are much
less restricted than in the isotropic case. One-point functions are still
vanishing (for non-compact spacetime), because of translation and scale
invariance. However, for the two-point function we have:

$\langle\phi(\tau,x)\phi(0,0)\rangle=|x|^{-2\Delta_{\phi}}f_{x}\left(\frac{|\tau|^{1/\zeta}}{|x|}\right)=|\tau|^{-\frac{2\Delta_{\phi}}{\zeta}}f_{\tau}\left(\frac{|x|}{|\tau|^{1/\zeta}}\right)\;,$
(4.4)

where we have expressed it in two equivalent ways, and we expect the
$f_{x}(u)$ and $f_{\tau}(u)$ to be two functions that are regular in zero and
vanishing at infinity as $u^{-2\Delta_{\phi}}$, but otherwise completely
unconstrained.

Even if the fixed-point theories are not conformally invariant, we can still
make use of the OPE. Indeed, assuming that a fixed-point Lifshitz field theory
admits a basis in the space of operators consisting in scaling operators
${\cal O}_{i}$ with dimensions $\Delta_{i}$, such that ${\cal
O}_{i}(\Omega^{\zeta}\tau,\Omega x)=\Omega^{-\Delta_{i}}{\cal O}_{i}(\tau,x)$,
the OPE expansion writes:

${\cal O}_{i_{1}}(\tau_{1},x_{1}){\cal
O}_{i_{2}}(\tau_{2},x_{2})=\sum_{k}f_{i_{1}i_{2}k}(\tau_{12},x_{12})\;{\cal
O}_{k}(\tau_{2},x_{2})\;,$ (4.5)

where $\tau_{12}=\tau_{1}-\tau_{2}$ and $x_{12}=x_{1}-x_{2}$. Applying a
scaling transformation and the OPE in either order we have:

$\begin{split}{\cal O}_{i_{1}}(\Omega^{\zeta}\tau_{1},\Omega x_{1}){\cal
O}_{i_{2}}(\Omega^{\zeta}\tau_{2},\Omega
x_{2})&=\Omega^{-\Delta_{i_{1}}-\Delta_{i_{2}}}{\cal
O}_{i_{1}}(\tau_{1},x_{1}){\cal O}_{i_{2}}(\tau_{2},x_{2})\\\
&=\Omega^{-\Delta_{i_{1}}-\Delta_{i_{2}}}\sum_{k}f_{i_{1}i_{2}i}(\tau_{12},x_{12})\;{\cal
O}_{k}(\tau_{2},x_{2})\;,\crcr{\cal O}_{i_{1}}(\Omega^{\zeta}\tau_{1},\Omega
x_{1}){\cal O}_{i_{2}}(\Omega^{\zeta}\tau_{2},\Omega
x_{2})&=\sum_{k}f_{i_{1}i_{2}k}(\Omega^{\zeta}\tau_{12},\Omega x_{12})\;{\cal
O}_{k}(\Omega^{\zeta}\tau_{2},\Omega x_{2})\\\
&=\sum_{k}f_{i_{1}i_{2}k}(\Omega^{\zeta}\tau_{12},\Omega
x_{12})\Omega^{-\Delta_{k}}{\cal O}_{k}(\tau_{2},x_{2})\;,\end{split}$ (4.6)

hence the OPE coefficient is a scaling function
$f_{i_{1}i_{2}k}(\Omega^{\zeta}\tau,\Omega
x)=\Omega^{\Delta_{k}-\Delta_{i_{1}}-\Delta_{i_{2}}}f_{i_{1}i_{2}k}(\tau,x)$,
that is, we can always write:

$f_{ijk}(\tau,x)=|\tau|^{\frac{\Delta_{k}-\Delta_{i}-\Delta_{j}}{\zeta}}f_{ijk}\left(\frac{|x|}{|\tau|^{1/\zeta}}\right)\;.$
(4.7)

This will be useful in the finite-temperature case (compact time), in order to
relate the two-point function to the one-point functions via the zero-
temperature OPE coefficients.

We will first discuss the free theory as a warm up, and then move on to the
interacting case.

###  4.1 The massless free fractional Lifshitz field theory at zero
temperature

The action of the massless free FLFT at zero temperature is (we suppress the
$O(N)$ indices here as the $N$ fields are decoupled in this case):

$S[\phi]=\int_{-\infty}^{+\infty}d\tau\int_{\mathbb{R}^{d}}d^{d}x\left(\frac{1}{2}\phi(\tau,x)(-\partial_{0}^{2})\phi(\tau,x)+\frac{1}{2}\phi(\tau,x)(-\partial_{i}\partial_{i})^{\zeta}\phi(\tau,x)\right)\;.$
(4.8)

It is trivial to check that this action is invariant under Lifshitz scaling
(4.3). This implies that the two-point function, or covariance $C(\tau,x)$ of
the Gaussian measure, must have the form (4.4). Indeed, writing the direct
space covariance as a Fourier transform:

$C(\tau,x)=\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{2\pi}\;\frac{e^{\imath
p\cdot x+\imath\omega\tau}}{(p^{2})^{\zeta}+\omega^{2}}\;,$ (4.9)

the integral over $\omega$ can be performed easily, by deforming the contour
and picking up a pole:

$C(\tau,x)=\int\frac{d^{d}p}{(2\pi)^{d}}e^{\imath p\cdot
x}\frac{e^{-|\tau|(p^{2})^{\zeta/2}}}{2(p^{2})^{\zeta/2}}\;,$ (4.10)

where the integral is convergent for $\zeta<d$. Simple rescalings reproduce
(4.4) and, recalling that $2\Delta_{\phi}=d-\zeta$, we have the limit
cases:202020Keeping both $\tau$ and $x$ different from zero, we can reduce to
single integral representations in various ways. For instance, we can use
spherical coordinates in (4.10), perform the angular integral and remain with
an integral over the radial part, involving an hypergeometric function.
Alternatively, and perhaps more interestingly, we can use again a spectral
representation of the integrand, that is, perform a Stieltjes transform with
respect to $p^{2}$ and then use the standard result for the short-range
massive propagator, to obtain:
$C(\tau,x)=\int_{0}^{+\infty}\frac{ds}{(2\pi)^{d/2+1}}\left(\frac{\sqrt{s}}{|x|}\right)^{\frac{d-2}{2}}K_{\frac{d-2}{2}}(\sqrt{s}|x|)\;s^{-\frac{\zeta}{2}}\operatorname{Im}\left(e^{\frac{\imath\pi}{2}\zeta}\;e^{-|\tau|s^{\frac{\zeta}{2}}e^{-\frac{\imath\pi}{2}\zeta}}\right)\;.$

$\begin{split}&C(\tau,0)=|\tau|^{-\frac{2\Delta_{\phi}}{\zeta}}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{e^{-(p^{2})^{\zeta/2}}}{2(p^{2})^{\zeta/2}}=\frac{\Gamma\left(\frac{d-\zeta}{\zeta}\right)}{(4\pi)^{d/2}\zeta\Gamma\left(\frac{d}{2}\right)}|\tau|^{-\frac{2\Delta_{\phi}}{\zeta}}\;,\\\
&C(0,x)=\int\frac{d^{d}p}{(2\pi)^{d}}e^{\imath p\cdot
x}\frac{1}{2(p^{2})^{\zeta/2}}=\frac{c(\Delta_{\phi})}{2}|x|^{-2\Delta_{\phi}}\;,\end{split}$
(4.11)

and we notice that $C(0,x)$ is (half) the GFFT covariance of section 2, with
the replacement $\zeta\to\zeta/2$.

#### Bilinear operators.

As the theory is Gaussian and massless, we expect again that in the (non-
conformal) OPE of two fundamental fields, $\phi\times\phi$, only bilinear
operators in $\phi$ appear, that is, operators of the schematic form:

$[\phi\phi]_{k,l,\mu_{1}\cdots\mu_{J}}(\tau,x)=\,:\phi(\tau,x)(\partial_{\tau})^{k}(\partial_{i}\partial_{i})^{l}\left(\partial_{i_{1}}...\partial_{i_{J}}-\text{traces}\right)\phi(\tau,x):\;,$
(4.12)

whose scaling dimension is

$\Delta_{k,l,J}=2\Delta_{\phi}+J+2l+k\zeta\;.$ (4.13)

Notice that while we have the usual degeneracy for a fixed integer $n$ such
that $J+2l=2n$, the time derivatives are in general (for non-rational $\zeta$)
unambiguously identifiable. On the other hand, lacking special conformal
transformations to select primary operators, the distribution of derivatives
on the two fields is rather arbitrary.

###  4.2 The massless free fractional Lifshitz field theory at finite
temperature

The massless free FLFT at finite temperature is defined by the action:

$S[\phi]=\int_{0}^{\beta}d\tau\int_{\mathbb{R}^{d}}d^{d}x\left(\frac{1}{2}\phi(\tau,x)(-\partial_{\tau}^{2})\phi(\tau,x)+\frac{1}{2}\phi(\tau,x)(-\partial_{i}\partial_{i})^{\zeta}\phi(\tau,x)\right)\;.$
(4.14)

In this case, the covariance, or two-point function, is:

$C_{\beta}(\tau,x)=\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\,\frac{e^{\imath
p\cdot x+\imath\omega_{n}\tau}}{p^{2\zeta}+\omega_{n}^{2}}\;,$ (4.15)

where $\omega_{n}=2\pi n/\beta$ and the subscript $\beta$ tracks the inverse
temperature. In this notation, the covariance at zero temperature of the
previous section is denoted $C_{\infty}(\tau,x)$. After using the Poisson
summation formula, $C_{\beta}$ becomes:

$C_{\beta}(\tau,x)=\sum_{m=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{(2\pi)}\,\frac{e^{\imath
p\cdot x+\imath\omega\tau_{m}}}{p^{2\zeta}+\omega^{2}}\;,$ (4.16)

where $\tau_{m}=\tau+\beta m$. Following the same steps as in the previous
section we get the massless covariance at finite temperature:

$C_{\beta}(\tau,x)=\sum_{m=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}e^{\imath
p\cdot x}\frac{e^{-|\tau_{m}|(p^{2})^{\zeta/2}}}{2(p^{2})^{\zeta/2}}\;,$
(4.17)

which becomes, after summing over $m$ (restricting to the fundamental domain
$0\leq\tau<\beta$):212121 For $0\leq\tau<\beta$, we have
$\sum_{m}e^{-a|\tau+\beta m|}=e^{-a\tau}+\sum_{m\geq 1}e^{-a\tau-a\beta
m}+\sum_{m\geq 1}e^{a\tau-a\beta
m}=e^{-a\tau}+\frac{e^{a\tau}+e^{-a\tau}}{e^{a\beta}-1}$ for any $a$ with
positive real part.

$C_{\beta}(\tau,x)=\int\frac{d^{d}p}{(2\pi)^{d}}\,e^{\imath p\cdot
x}\frac{1}{2(p^{2})^{\zeta/2}}\,\left(e^{-\tau(p^{2})^{\frac{\zeta}{2}}}+\frac{e^{\tau(p^{2})^{\frac{\zeta}{2}}}+e^{-\tau(p^{2})^{\frac{\zeta}{2}}}}{e^{\beta(p^{2})^{\frac{\zeta}{2}}}-1}\right)\;.$
(4.18)

The first term in parentheses is the $m=0$ mode and reproduces the zero-
temperature result. This is the only surviving term in the $\beta\to\infty$
limit. The second term captures the finite temperature effects.

As in the case of the GFFT at finite size, the two-point function can be seen
as the generating function for the one-point functions of scaling operators.
For example if we set $x=0$ and Taylor expand in $\tau$ in (4.15) we find:

$C_{\beta}(\tau,0)=\sum_{k\geq
0}\frac{(-1)^{k}\tau^{2k}}{(2k)!}\left(\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{\omega_{n}^{2k}}{p^{2\zeta}+\omega_{n}^{2}}\right)\;,$
(4.19)

and the coefficient of $\tau^{2k}$ is the one-point function of the bilinear
scaling operator ${\cal
O}_{k}(\tau^{\prime},x)=:\phi\partial_{\tau}^{2k}\phi(\tau^{\prime},x):$ at
finite temperature.222222Operators with odd number of time derivatives vanish
even at finite temperature, due to time-reversal symmetry. As the integral in
(4.19) is divergent, one needs to regularize it by an explicit subtraction.
This is fairly straightforward in this case: in order to cure all the
ultraviolet divergences it is enough to subtract the zero-temperature two-
point function,

$\begin{split}&\sum_{k\geq 0}\frac{(-1)^{k}\tau^{2k}}{(2k)!}\braket{{\cal
O}_{k}(\tau^{\prime},x)}=C_{\beta}(\tau,0)-C_{\infty}(\tau,0)\crcr&\qquad=\sum_{k\geq
0}\frac{(-1)^{k}\tau^{2k}}{(2k)!}\left(\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{\omega_{n}^{2k}}{p^{2\zeta}+\omega_{n}^{2}}-\int\frac{d\omega}{2\pi}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{\omega^{2k}}{p^{2\zeta}+\omega^{2}}\right)\;,\end{split}$
(4.20)

which comes to subtracting the $m=0$ mode after Poisson summation.

This subtraction can be understood directly from the representation (4.18),
which we now rewrite in $d=2$ and at $x=0$, with the change of variables
$t=p^{2}$:

$C_{\beta}(\tau,0)=\frac{1}{8\pi}\int_{0}^{\infty}dt\;\;t^{-\frac{\zeta}{2}}\;\bigg{(}e^{-\tau
t^{\frac{\zeta}{2}}}+\frac{e^{\tau t^{\frac{\zeta}{2}}}+e^{-\tau
t^{\frac{\zeta}{2}}}}{e^{\beta t^{\frac{\zeta}{2}}}-1}\bigg{)}\;.$ (4.21)

A naive Taylor expansion of the integrand yields ultraviolet (large $t$)
divergent terms at every order, coming form the first term (the zero-
temperature $m=0$ mode):

$\frac{1}{8\pi}\int_{0}^{\infty}dt\;\sum_{k\geq
0}\frac{(-1)^{k}}{k!}\tau^{k}\;t^{\frac{\zeta}{2}(k-1)}\;.$ (4.22)

However, once they are summed over $k$ they reconstitute $C_{\infty}(\tau,0)$
(compare to (4.11)):

$\frac{1}{8\pi}\int_{0}^{\infty}dt\;\;t^{-\frac{\zeta}{2}}\;e^{-\tau
t^{\frac{\zeta}{2}}}=\tau^{-\frac{2\Delta_{\phi}}{\zeta}}\;\frac{\Gamma(2/\zeta-1)}{4\pi\zeta}\equiv\tau^{-\frac{2\Delta_{\phi}}{\zeta}}\;d_{-2\Delta_{\phi}/\zeta}\;,$
(4.23)

which is unsurprising as the divergence came from expanding in $\tau$ the
$m=0$ mode.

Subtracting these terms in the spectral representation and Taylor expanding
the second term in $\tau$ we identify the one-point functions of the scaling
operators ${\cal O}_{k}$:

$(-1)^{k}\braket{{\cal
O}_{k}(\tau^{\prime},x)}=\frac{1}{4\pi}\int_{0}^{\infty}dt\;\;\frac{t^{\zeta(k-\frac{1}{2})}}{e^{\beta
t^{\frac{\zeta}{2}}}-1}=\frac{\Gamma(2\Delta_{\phi}/\zeta+2k)\text{Li}_{2\Delta_{\phi}/\zeta+2k}(1)}{2\pi\zeta\,\beta^{2\Delta_{\phi}/\zeta+2k}}\equiv\beta^{-\frac{2\Delta_{\phi}}{\zeta}-2k}\;d_{2k}\;.$
(4.24)

Observe that the counterterm part
$C_{\infty}(\tau,0)\sim\tau^{-2\Delta_{\phi}/\zeta}$ is non-analytic at
$\tau=0$, that is, the two-point function is:

$C_{\beta}(\tau,0)=\tau^{-\frac{2\Delta_{\phi}}{\zeta}}\;d_{-2\Delta_{\phi}/\zeta}+\sum_{k\geq
0}\frac{\tau^{2k}}{(2k)!}\;\beta^{-\frac{2\Delta_{\phi}}{\zeta}-2k}d_{2k}\;,$
(4.25)

and the subtraction is designed to take out the non-analytic piece of the two-
point function, corresponding to the identity operator’s contribution to the
OPE.

###  4.3 The interacting fractional Lifshitz field theory at zero temperature

We now turn on the quartic interaction in the zero-temperature theory, that
is, we consider the full model (4.1) on $\mathbb{R}^{d+1}$, and we are going
to show that at large $N$ the theory exhibits an interacting infrared fixed
point, much like its isotropic counterpart.

The Schwinger-Dyson equation at zero temperature relates the free covariance
$\tilde{C}^{-1}(\omega,p)=\omega^{2}+(p^{2})^{\zeta}$ and the full interacting
propagator (two-point function). As the combinatorics is the same as in the
usual $O(N)$ model, at large $N$ we have the usual dominance of cactus
diagrams, and thus the Schwinger-Dyson equation reads:

$\tilde{G}(\omega_{0},p_{0})^{-1}=\tilde{C}(\omega_{0},p_{0})^{-1}+M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{2\pi}\;\tilde{G}(\omega,p)\,,$
(4.26)

which as usual gives only a mass correction to the free propagator. Like in
the isotropic case, we tune the bare mass in order to cancel such mass
correction and obtain a massless propagator
$\tilde{G}^{-1}(\omega,p)=\omega^{2}+(p^{2})^{\zeta}$, that is, we impose the
condition:

$0=M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{2\pi}\;\frac{1}{\omega^{2}+(p^{2})^{\zeta}}\equiv
M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\mathcal{T}\;,$ (4.27)

where of course a ultraviolet regularization is implicit for the massless
tadpole $\mathcal{T}=C(0,0)$.

In order to show that the model exhibits an infrared fixed point in the
large-$N$ limit, we derive the beta function of the quartic coupling. As usual
in the vector large-$N$ limit, the beta function is one-loop exact and at one
loop only the bubble graph contributes. We regularize the infrared by
performing a subtraction at a non-zero external momentum $(\omega_{0},0)$ such
that $\omega_{0}=\mu^{\zeta}$.232323Alternatively we can subtract at an
external momentum $(0,p)$ such that $p^{2}=\mu^{2}$. The computation in this
case is slightly more involved but leads to the same result, as expected. The
bubble integral with massless propagator is then:

$\begin{split}B&=\int\frac{d\omega}{2\pi}\int\frac{d^{d}q}{(2\pi)^{d}}\frac{1}{(\omega^{2}+(q^{2})^{\zeta})((\omega+\omega_{0})^{2}+(q^{2})^{\zeta})}\\\
&=\int\frac{d^{d}q}{(2\pi)^{d}}\frac{(q^{2})^{-\zeta/2}}{4(q^{2})^{\zeta}+\omega_{0}^{2}}=-\frac{2^{1-\frac{d}{\zeta}}\pi\omega_{0}^{\frac{d}{\zeta}-3}}{(4\pi)^{\frac{d}{2}}\Gamma(\frac{d}{2})\zeta\cos(\frac{\pi
d}{2\zeta})}\,,\end{split}$ (4.28)

where we integrated over $\omega$ by deforming the contour and using the
residue theorem. In the weakly relevant case $\zeta=\frac{d+\epsilon}{3}$,
this integral exhibits a pole
$B=\frac{\mu^{-\epsilon}}{2(4\pi)^{d/2}\Gamma(d/2)\epsilon}$. The large-$N$
proper vertex at the subtraction point is:

$\Gamma^{(4)}(\omega_{0},0)=\frac{2\lambda}{1+\lambda B}\,,$ (4.29)

and defining the running coupling
$g=\mu^{-\epsilon}\Gamma^{(4)}(\omega_{0},0)/2$, we obtain the large-$N$ beta
function:

$\beta(g)=\mu\partial_{\mu}g=-\epsilon
g+\frac{g^{2}}{2(4\pi)^{d/2}\Gamma(d/2)}\,.$ (4.30)

This beta function displays a non-trivial fixed point
$g^{\star}=2(4\pi)^{d/2}\Gamma(d/2)\epsilon$, corresponding to a bare coupling
$\lambda$ flowing to infinity.

###  4.4 The interacting fractional Lifshitz field theory at finite
temperature

At last, we turn on also the temperature, that is, we compactify the Euclidean
time direction and study the full model (4.1) on
$S_{\beta}^{1}\times\mathbb{R}^{d}$, and we are going to discuss the thermal
mass and one-point functions of the scaling operators.

The model is tuned to the same fixed point, still at large $N$, but the
compactification of time implies that the tadpole is no longer canceled by the
bare mass in (4.27), and therefore we have a solution of the Schwinger-Dyson
equation of the form
$\tilde{G}_{\beta}^{-1}(\omega_{n},p)=\omega_{n}^{2}+(p^{2})^{\zeta}+M_{th}^{2\zeta}$,
with $\omega_{n}=2\pi n/\beta$ and non-vanishing thermal mass $M_{th}$. The
thermal mass is fixed self-consistently by the mass gap equation:

$\begin{split}M_{th}^{2\zeta}&=M_{\text{bare}}^{2\zeta}+\frac{\lambda}{2}\;\frac{1}{\beta}\sum_{n}\int\frac{d^{d}p}{(2\pi)^{d}}\frac{1}{\omega_{n}^{2}+(p^{2})^{\zeta}+M_{th}^{2\zeta}}\\\
&=\frac{\lambda}{2}\Bigg{(}\sum_{m\neq
0}\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{2\pi}\;\frac{e^{\imath\beta
m\omega}}{\omega^{2}+(p^{2})^{\zeta}+M_{th}^{2\zeta}}\\\ &\qquad\qquad-
M_{th}^{2\zeta}\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{2\pi}\;\frac{1}{(\omega^{2}+(p^{2})^{\zeta}+M_{th}^{2\zeta})(\omega^{2}+(p^{2})^{\zeta})}\Bigg{)}\\\
&\equiv\frac{\lambda}{2}\,\left(\mathcal{T}(M_{th})-\mathcal{T}\right)\;,\end{split}$
(4.31)

where in the second equality we used the Poisson summation formula (2.10), and
we have combined the $m=0$ term with the explicit expression of the bare mass.
In the last step we have defined the subtracted thermal tadpole
$\mathcal{T}(M_{th})-\mathcal{T}=G_{\beta}(0,0)-C_{\infty}(0,0)$.

At the fixed point, the bare coupling $\lambda$ diverges, and thus the mass
gap equation becomes:

$\mathcal{T}(M_{th})-\mathcal{T}=0\;.$ (4.32)

We are going to show that it admits a unique non-vanishing solution for
$\zeta\in(d/3,1)$.

First, we write the two-point function at finite temperature in position space
in the Fourier representation:

$G_{\beta}(\tau,x)=\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\;\frac{e^{\imath
p\cdot
x+\imath\omega_{n}\tau}}{p^{2\zeta}+\omega_{n}^{2}+M_{th}^{2\zeta}}=\sum_{m=-\infty}^{+\infty}\int\frac{d^{d}p}{(2\pi)^{d}}\int\frac{d\omega}{2\pi}\frac{e^{\imath
p\cdot x+\imath\omega\tau_{m}}}{p^{2\zeta}+\omega^{2}+M_{th}^{2\zeta}}\;,$
(4.33)

with $\tau_{m}=\tau+\beta m$. The integral over $\omega$ is evaluated again by
the residue theorem and restricting to the fundamental domain
$\tau\in[0,\beta)$ the sum over $m$ is just a geometric series (see footnote
21), leading to:

$G_{\beta}(\tau,x)=\int\frac{d^{d}p}{(2\pi)^{d}}e^{\imath p\cdot
x}\frac{1}{2\sqrt{M_{th}^{2\zeta}+(p^{2})^{\zeta}}}\bigg{(}e^{-\tau\sqrt{M_{th}^{2\zeta}+(p^{2})^{\zeta}}}+\frac{e^{\tau\sqrt{M_{th}^{2\zeta}+(p^{2})^{\zeta}}}+e^{-\tau\sqrt{M_{th}^{2\zeta}+(p^{2})^{\zeta}}}}{e^{\beta\sqrt{M_{th}^{2\zeta}+(p^{2})^{\zeta}}}-1}\bigg{)}\;,$
(4.34)

where the first term in parentheses is the $m=0$ term.

Setting $d=2$, taking $\tau=x=0$, and changing variables to $p^{2}=t$, we
arrive at the following integral form of the mass gap equation:

$\mathcal{T}(M_{th})-\mathcal{T}=\frac{1}{8\pi}\int_{0}^{\infty}dt\;\left[\frac{1}{\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}\left(1+\frac{2}{e^{\beta\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}-1}\right)-\frac{1}{\sqrt{t^{\zeta}}}\right]=0\;,$
(4.35)

that always admits a unique solution.242424We rewrite the mass gap equation as
$f_{1}(M_{th})=f_{2}(M_{th})$ with:
$f_{1}(M_{th})=\int_{0}^{\infty}dt\;\frac{1}{\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}\left(\frac{2}{e^{\beta\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}-1}\right)\;,\qquad
f_{2}(M_{th})=\int_{0}^{\infty}dt\;\left[\frac{1}{\sqrt{t^{\zeta}}}-\frac{1}{\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}\right]\;,$
(4.36) and $f_{2}$ is convergent for $\zeta>\frac{2}{3}$ as
$f_{2}\sim\int^{\infty}dt\;t^{-3\zeta/2}$. We note that $f_{1}(0)>0$ and
$f_{1}$ is strictly decreasing with $M_{th}$ while $f_{2}(0)=0$ and $f_{2}$ is
strictly increasing with $M_{th}$. We have plotted a numerical solution as a
function of $\zeta$ in Figure 4.

Figure 4: Numerical solutions for the thermal mass as a function of $\zeta$ in
$2+1$ dimensions. The red dot is the value of the short-range thermal mass,
and it is approached in the $\zeta\rightarrow 1$ limit, while the solution
goes to zero for $\zeta\to 2/3$.

Keeping $\tau\neq 0$, we obtain instead:

$G_{\beta}(\tau,0)=\frac{1}{8\pi}\int_{0}^{\infty}dt\;\frac{1}{\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}\left(e^{-\tau\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}+\frac{e^{\tau\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}+e^{-\tau\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}}{e^{\beta\sqrt{M_{th}^{2\zeta}+t^{\zeta}}}-1}\right)\;,$
(4.37)

and we are interested in its small $\tau$ expansion, which again we relate to
the one-point function of the bilinear scaling operator ${\cal
O}_{k}(\tau^{\prime},x)=:\phi\partial_{\tau}^{2k}\phi(\tau^{\prime},x):$ at
finite temperature:

$\sum_{k\geq 0}\frac{(-1)^{k}\tau^{2k}}{(2k)!}\braket{{\cal
O}_{k}(\tau^{\prime},x)}=G_{\beta}(\tau,0)-\text{counterterms}\;.$ (4.38)

As in the isotropic case, the $k=0$ operator will actually drop from this
equation, as it is not in the spectrum of the interacting fixed point, where
it is replaced by the intermediate field $\sigma$. Moreover, as in section
3.5, the counterterms will contain useful information about other operators in
the spectrum.

The need for counterterms comes from the first term in parentheses in (4.37),
i.e. the $m=0$ mode: its naive Taylor expansion in $\tau$ exhibits ultraviolet
divergent terms at each order. Such divergences are renormalized by
subtracting the beginning of the Taylor expansion in $M^{2\zeta}$ up to the
appropriate order:

$\frac{1}{8\pi}\int_{0}^{\infty}dt\;\sum_{q\geq
0}\frac{(-\tau)^{q}}{q!}\bigg{[}(M_{th}^{2\zeta}+t^{\zeta})^{\frac{q-1}{2}}-\sum_{n=0}^{\lfloor\frac{q}{2}\rfloor}\frac{\Gamma(\frac{q+1}{2})}{n!\Gamma(\frac{q+1}{2}-n)}(M_{th}^{2\zeta})^{n}t^{\zeta(\frac{q-1}{2}-n)}\bigg{]}\;.$
(4.39)

The odd $q=2k+1$ terms are entirely subtracted, since
$(M_{th}^{2\zeta}+t^{\zeta})^{k}$ is a polynomial in $M^{2\zeta}$, while for
even $q=2k$ only the terms containing a power of $t$ larger than $-1$ are
subtracted:

$\frac{1}{8\pi}\frac{\tau^{2k}}{(2k)!}\int_{0}^{\infty}dt\;\bigg{[}(M_{th}^{2\zeta}+t^{\zeta})^{k-\frac{1}{2}}-\sum_{n=0}^{k}\frac{\Gamma(k+\frac{1}{2})}{n!\Gamma(k-n+\frac{1}{2})}(M_{th}^{2\zeta})^{n}t^{\zeta(k-n-\frac{1}{2})}\bigg{]}\;,$
(4.40)

resulting in an expression that is both ultraviolet ($t\to\infty$) and
infrared ($t\to 0$) convergent for $\zeta>2/3$. Indeed, convergence in the
infrared is obvious term by term; convergence in the ultraviolet follows by
noting that the subtracted integrand behaves like $t^{-\frac{3}{2}\zeta}$ at
large $t$. Grouping the counterterms by their power in $M_{th}$, we obtain the
following small-$\tau$ expansion of the massive two-point function in finite
temperature (recall that $2\Delta_{\phi}=2-\zeta$):

$\displaystyle
G_{\beta}(\tau,0)=\tau^{-\frac{2\Delta_{\phi}}{\zeta}}\left(\sum_{n=0}^{\infty}(\tau
M_{th}^{\zeta})^{2n}d_{-2\Delta_{\phi}/\zeta+2n}+\sum_{k\geq 0}\frac{(\tau
M_{th}^{\zeta})^{2k+\frac{2\Delta_{\phi}}{\zeta}}}{(2k)!}d_{2k}^{(\beta
M_{th}^{\zeta})}\right)\;,$ (4.41)

exhibiting the finite-temperature scaling behavior
$G_{\beta}(\tau,0)=\tau^{-2\Delta_{\phi}/\zeta}\,h\left(\tau
M_{th}^{\zeta}\right)$ for a function $h(u)$ that is finite, but non-analytic,
at $u=0$. Taking into account the prefactor, the non-analytic part in $\tau$
of $G_{\beta}(\tau,0)$ is the first term in parentheses in (4.41), while the
second term gives the analytic part.

The coefficients of the non-analytic part of $G_{\beta}(\tau,0)$ are:

$\begin{split}d_{-2\Delta_{\phi}/\zeta+2n}&=\int_{0}^{\infty}\frac{dt}{8\pi}\;\sum_{q\geq
2n}\frac{(-1)^{q}}{q!}\;\frac{\Gamma(\frac{q+1}{2})}{n!\Gamma(\frac{q+1}{2}-n)}\,t^{\zeta(\frac{q-1}{2}-n)}\\\
&=\int_{0}^{\infty}\frac{dt}{8\pi}\;\left(\frac{{}_{1}F_{2}\left(1;\frac{1}{2},n+1;\frac{t^{\zeta}}{4}\right)}{4^{n}\,t^{\zeta/2}(n!)^{2}}-\frac{\,{}_{0}F_{1}\left(;n+\frac{3}{2};\frac{t^{\zeta}}{4}\right)}{(2n+1)!}\right)\;,\end{split}$
(4.42)

and, although each term at $q$ fixed is ultraviolet divergent, performing the
sum over $q$ first we obtain a convergent integral expression due to the
alternating sign (the linear combination of hypergeometric functions decays
like $t^{-\zeta(n+1)/2}e^{-t^{\zeta/2}}$ for $t\to+\infty$).

The analytic part of $G(\tau,0)$ comprises effects due to both the finite
temperature and the presence of a non-zero thermal mass:

$\begin{split}d_{2k}^{(\beta
M_{th}^{\zeta})}=&\int_{0}^{\infty}\frac{dt}{8\pi}\;\bigg{[}(1+t^{\zeta})^{k-\frac{1}{2}}-\sum_{n=0}^{k}\frac{\Gamma(k+\frac{1}{2})}{n!\Gamma(k-n+\frac{1}{2})}t^{\zeta(k-n-\frac{1}{2})}+\frac{2\left(1+t^{\zeta}\right)^{k-\frac{1}{2}}}{e^{\beta
M_{th}^{\zeta}\sqrt{1+t^{\zeta}}}-1}\bigg{]}\\\
=&\int_{0}^{\infty}\frac{dt}{8\pi}\;\bigg{[}-\frac{\Gamma\left(k+\frac{1}{2}\right)\,_{2}{F}_{1}\left(1,\frac{3}{2};k+2;-t^{-\zeta}\right)}{2\sqrt{\pi}\,(k+1)!\,t^{3\zeta/2}}+\frac{2\left(1+t^{\zeta}\right)^{k-\frac{1}{2}}}{e^{\beta
M_{th}^{\zeta}\sqrt{1+t^{\zeta}}}-1}\bigg{]}\;,\end{split}$ (4.43)

and we stress again that this is a convergent integral for $\zeta>2/3$.

As in the CFT case, the coefficients of the small-$\tau$ expansion can be
related to the one-point functions of composite operators. Indeed, from (4.5)
and (4.7), we have:

$\phi(\tau,0)\phi(0,0)=\sum_{\cal O}f_{\phi\phi{\cal
O}}(0)\;|\tau|^{\frac{\Delta_{\cal O}-2\Delta_{\phi}}{\zeta}}\;{\cal
O}(0,0)\;,$ (4.44)

where ${\cal O}$ runs over all the scaling operators in the theory. Taking the
expectation value of (4.44), we conclude that the two-point function is a
generating function for the one-point functions of the operators ${\cal O}$
having $f_{\phi\phi{\cal O}}(0)\neq 0$. Notice that the operators in (4.44)
are not primaries, as we do not have conformal invariance here, hence they
include also total derivatives of composite operators. However, the latter
have vanishing one-point functions even at finite temperature, due to
translation invariance.

Seen from such an OPE perspective, we recognize that (4.41) is indexed by two
classes of operators. First, we have operators with momentum dimensions
$\frac{\Delta_{\cal O}-2\Delta_{\phi}}{\zeta}=2k$, i.e. $\Delta_{\cal
O}=2\Delta_{\phi}+2k\zeta$. These are identified with the operators ${\cal
O}_{k}(\tau,x)=:\phi\partial_{\tau}^{2k}\phi(\tau,x):$. By consistency with
(4.38), we must have that with our choices of normalizations $f_{\phi\phi{\cal
O}}(0)=(-1)^{k}/(2k)!$, and thus we conclude that:

$(-1)^{k}\braket{{\cal
O}_{k}(\tau^{\prime},x)}=M_{th}^{2\Delta_{\phi}+2k\zeta}d_{2k}^{(\beta
M_{th}^{\zeta})}\;.$ (4.45)

Notice that the mass gap equation corresponds to $d_{0}^{(\beta M_{th})}=0$.

Second, we have a family of operators with momentum dimension
$\frac{\Delta_{\cal
O}-2\Delta_{\phi}}{\zeta}=2n-\frac{2\Delta_{\phi}}{\zeta}$, i.e. $\Delta_{\cal
O}=2n\zeta$. As in the isotropic case, the intermediate field $\sigma(x)$,
which can be introduced as in (3.2), is identified with the shadow of
$\phi^{2}(x)$ and has momentum dimension
$\Delta_{\sigma}=d+\zeta-2\Delta_{\phi}=2\zeta$. Therefore, the second class
of operators contributing to the OPE of the two-point function correspond to
composite $\sigma^{n}$ operators:

$f_{\phi\phi\sigma^{n}}(0)\braket{\sigma^{n}(\tau^{\prime},x)}=M_{th}^{2n\zeta}d_{-2\Delta_{\phi}/\zeta+2n}\;.$
(4.46)

##  5 Conclusions

We have studied classical and quantum versions of the long-range $O(N)$ model
at large $N$, focusing on the effects of having one compact spatial direction
in the former or finite temperature in the latter. The quantum model is mapped
as usual to a Euclidean field theory with one extra (time) dimension, along
which the action is local, thus resulting in an anisotropic action, which at
criticality has a Lifshitz scale invariance, but no conformal invariance. This
should be contrasted to the short-range case, where the classical and quantum
model are described by the same field theory, in $d$ and $d+1$ dimensions,
respectively, and therefore they are both conformal at criticality. Therefore,
in the long-range model also the finite-size and finite-temperature
compactifications are quite different: the finite-size geometry corresponds to
compactifying a direction along which the model has long-range interaction,
while the finite temperature corresponds to compactifying the local Euclidean
time direction.

In both cases, the non-compact versions of the model exhibit an ultraviolet
Gaussian fixed point and an interacting infrared fixed point. Upon
compactification, the fixed-point data (conformal data in the classical case,
or simply fixed-point data in the quantum one) of the non-compact theories
needs to be supplemented with the non-zero values of one-point functions of
primary (respectively scaling) operators.

We have focused on the operators arising in the OPE of two fundamental fields.
The corresponding new data can be inferred from the small argument behaviour
of the two-point function, and consequently our main technical tool consists
in finding suitable expansions of the latter. Our results are:

  * •

for the free (ultraviolet) compactified theories, we obtained the one-point
functions of bilinear operators. In the finite-size case we obtained such data
for all the bilinear primaries, with arbitrary spin and twist. In the finite-
temperature case we computed this data only for scalar bilinears with time
derivatives.

  * •

in the case of the compactified interacting (infrared) theories, at large $N$,
the effect of the interaction is encoded in the appearance of a dynamically
generated mass scale which we determined self-consistently. We obtained the
fixed-point data for the same set of bilinear operators as in the free case.
Moreover, in the interacting cases, a second infinite family of operators
contributes to the OPE, consisting in arbitrary powers of the
Hubbard–Stratonovich intermediate field. We computed also the one-point
functions of such operators.

A natural extension of this work would be to consider other configurations of
the background geometry, for example considering the classical model with more
or all directions compactified, or the quantum model with compact spatial
directions, either at zero or finite temperature.

Another important future direction would be the inclusion of subleading
corrections in $1/N$. In the short-range case, the effect of next-to-leading
corrections on the new finite-temperature data has been studied in [32]. In
the long-range case, we expect further complications. In the classical
isotropic case, we have seen that even in the leading-order we considered
here, the non-standard propagator makes analytic computations much harder, and
certainly this will affect the loop integrals that are needed for subleading
orders. In the anisotropic (quantum) case, computations are harder even at
zero temperature; in fact, as we mentioned in the introduction, the fractional
Lifshitz field theory we considered here has not been studied much before, and
there are several open directions even at zero temperature. In particular, it
would be interesting to determine the deviation of the dynamic critical
exponent $z$, characterizing the Lifshitz scaling at the fixed point, from its
leading order value $z=\zeta$, and to ascertain the existence of a crossover
(or transition) to isotropic short-range behavior above some critical value of
$\zeta$. Such aspects of the zero-temperature FLFT require at least the next-
to-leading order in $1/N$, therefore we postpone them to future work.

## Acknowledgments

Nordita is supported in part by NordForsk. R. G. and D. L. have been supported
by the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreement No 818066) and by the
Deutsche Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy EXC–
2181/1 – 390900948 (the Heidelberg STRUCTURES Cluster of Excellence).

##  Appendix A Bilinear operators

The most general bilinear primary operator with spin $J$ is generated by
acting with a precise combination of derivatives on the two fields. The
combination is fixed by requiring that the operator is annihilated by the
generator of special conformal transformations. The result of imposing such
constraint was obtained in [45, 50] (see also [51, 60]), and it is most easily
expressed by introducing a complex polarization vector $\xi^{\mu}$ such that
$\xi^{2}=0$ with which to form a scalar:

$\begin{split}\mathcal{O}_{0,J}(x,\xi)&\equiv[\phi_{a}\phi_{a}]_{0,\mu_{1}\cdots\mu_{J}}(x)\xi^{\mu_{1}}\cdots\xi^{\mu_{J}}=\frac{1}{N}:\phi_{a}(x)f_{J}\left(\xi\cdot\overset{\scriptscriptstyle\leftarrow}{\partial},\xi\cdot\overset{\scriptscriptstyle\rightarrow}{\partial}\right)\phi_{a}(x):\;,\end{split}$
(A.1)

where the colon notation stands as usual for normal (or Wick) ordering (e.g.
[61]), which serves to cancel the divergence in the one-point function. The
function $f_{J}(u,v)$ is found to be

$\begin{split}f_{J}(u,v)&=\frac{J!\,\Gamma\left(\nu^{(\zeta)}\right)}{4^{J}\,\Gamma\left(\nu^{(\zeta)}+J\right)}\left(u+v\right)^{J}C^{\nu^{(\zeta)}}_{J}\left(\frac{u-v}{u+v}\right)=\sum_{n=0}^{J/2}\gamma_{n,J}\left(u+v\right)^{2n}\left(u-v\right)^{J-2n}\;,\end{split}$
(A.2)

with $C^{\nu^{(\zeta)}}_{J}(z)$ being a Gegenbauer polynomial with
$\nu^{(\zeta)}=\frac{d-2\zeta-1}{2}$, and

$\gamma_{n,J}=(-1)^{n}\frac{J!\,}{2^{J}\,\Gamma\left(\nu^{(\zeta)}+J\right)}\frac{\Gamma\left(J+\nu^{(\zeta)}-n\right)}{2^{2n}\,n!(J-2n)!}\;.$
(A.3)

The normalization is chosen for convenience as in [50] (in particular it gives
finite coefficients for $J>0$ in $d=3$).

###  A.1 Normalization of two-point functions

The two-point function at zero temperature can be computed in the free theory
by applying the above formula to the (formal) spin-0 case [50]

$\langle\mathcal{O}_{0,J}(x_{1},\xi_{1})\mathcal{O}_{0,J}(x_{2},\xi_{2})\rangle=\frac{2}{N}C(x_{12})f_{J}\left(\xi_{1}\cdot\overset{\scriptscriptstyle\leftarrow}{\partial}_{1},\xi_{1}\cdot\overset{\scriptscriptstyle\rightarrow}{\partial}_{1}\right)f_{J}\left(\xi_{2}\cdot\overset{\scriptscriptstyle\leftarrow}{\partial}_{2},\xi_{2}\cdot\overset{\scriptscriptstyle\rightarrow}{\partial}_{2}\right)C(x_{12})\;.$
(A.4)

In order to extract the normalization constant $c_{0,J}$ in

$\langle\mathcal{O}_{0,J}(x,\xi_{1})\mathcal{O}_{0,J}(0,\xi_{2})\rangle=c_{0,J}\frac{\left(\xi_{1}^{\mu}\left(\delta_{\mu\nu}-2\frac{x_{\mu}x_{\nu}}{x^{2}}\right)\xi_{2}^{\nu}\right)^{J}}{(x^{2})^{d-2+J}}\;,$
(A.5)

it is enough to choose $\xi_{1}=\xi_{2}=\xi$, in which case the Kronecker
delta terms drop out. The derivatives can be evaluated explicitly using the
representation

$C(x)=\frac{\Gamma(d/2-\zeta)}{2^{2\zeta}\pi^{d/2}\Gamma(\zeta)}\frac{1}{\left(x^{2}\right)^{d/2-\zeta}}=\frac{1}{2^{2\zeta}\pi^{d/2}\Gamma(\zeta)}\int_{0}^{+\infty}d\alpha\,\alpha^{d/2-\zeta-1}e^{-\alpha
x^{2}}\;,$ (A.6)

and the differential formula

$(\xi\cdot\partial)^{k}e^{-\alpha x^{2}}=(-2\,\alpha\,\xi\cdot
x)^{k}e^{-\alpha x^{2}}\;.$ (A.7)

Noticing that $\partial_{1}C(x_{12})=-\partial_{2}C(x_{12})$ and remembering
that $J$ is even, we find

$\begin{split}&\langle\mathcal{O}_{0,J}(x,\xi)\mathcal{O}_{0,J}(0,\xi)\rangle=\frac{2}{N}\left(\frac{1}{2^{2\zeta}\pi^{d/2}\Gamma(\zeta)}\right)^{2}\int_{0}^{+\infty}d\alpha_{1}\,\alpha_{1}^{d/2-\zeta-1}\int_{0}^{+\infty}d\alpha_{2}\,\alpha_{2}^{d/2-\zeta-1}\\\
&\qquad\quad\qquad\qquad\times f_{J}(-2\,\alpha_{1}\,\xi\cdot
x,-2\,\alpha_{2}\,\xi\cdot x)^{2}e^{-(\alpha_{1}+\alpha_{2})x^{2}}\\\
&\qquad\quad=\frac{2}{N}\left(\frac{1}{2^{2\zeta}\pi^{d/2}\Gamma(\zeta)}\frac{J!\,\Gamma\left(\nu^{(\zeta)}\right)}{4^{J}\,\Gamma\left(\nu^{(\zeta)}+J\right)}\right)^{2}\frac{(-2\,\xi\cdot
x)^{2J}}{(x^{2})^{d-2\zeta+2J}}\int_{0}^{+\infty}d\alpha_{1}\,\alpha_{1}^{d/2-\zeta-1}\int_{0}^{+\infty}d\alpha_{2}\,\alpha_{2}^{d/2-\zeta-1}\\\
&\qquad\quad\qquad\qquad\times\left(\alpha_{1}+\alpha_{2}\right)^{2J}C^{\nu^{(\zeta)}}_{J}\left(\frac{\alpha_{1}-\alpha_{2}}{\alpha_{1}+\alpha_{2}}\right)^{2}e^{-(\alpha_{1}+\alpha_{2})}\\\
&\qquad\quad=\frac{2}{N}\left(\frac{1}{2^{2\zeta}\pi^{d/2}\Gamma(\zeta)}\frac{J!\,\Gamma\left(\nu^{(\zeta)}\right)}{4^{J}\,\Gamma\left(\nu^{(\zeta)}+J\right)}\right)^{2}\frac{(-2\,\xi\cdot
x)^{2J}}{(x^{2})^{d-2\zeta+2J}}2^{2\zeta+1-d}\int_{0}^{+\infty}ds\,s^{d-2\zeta+2J-1}\,e^{-s}\\\
&\qquad\quad\qquad\qquad\times\int_{-1}^{+1}dt\,(1-t^{2})^{\nu^{(\zeta)}-\frac{1}{2}}C^{\nu^{(\zeta)}}_{J}\left(t\right)^{2}\\\
&\qquad\quad=\frac{2}{N}\frac{\Gamma(d-2\zeta+2J)}{2^{4\zeta}\pi^{d}\Gamma(\zeta)^{2}}\left(\frac{J!\,\Gamma\left(\nu^{(\zeta)}\right)}{4^{J}\,\Gamma\left(\nu^{\zeta}+J\right)}\right)^{2}\,2^{J}\frac{\pi
2^{1-4\nu^{(\zeta)}}\Gamma(J+2\nu^{(\zeta)})}{J!(J+\nu^{(\zeta)})\Gamma(\nu^{(\zeta)})^{2}}\,\frac{(2\,(\xi\cdot
x)^{2})^{J}}{(x^{2})^{d-2\zeta+2J}}\;,\end{split}$ (A.8)

where we used the change of variables $\alpha_{1}=s(1+t)/2$,
$\alpha_{2}=s(1-t)/2$, and the orthogonality of Gegenbauer polynomials

$\int_{-1}^{+1}dt\,(1-t^{2})^{\alpha-1/2}C^{\alpha}_{J}\left(t\right)^{2}=\frac{\pi
2^{1-2\alpha}\Gamma(J+2\alpha)}{J!(J+\alpha)\Gamma(\alpha)^{2}}\;.$ (A.9)

We thus find

$\begin{split}c_{0,J}&=\frac{2}{N}\frac{\Gamma(d-2\zeta+2J)}{2^{4\zeta}\pi^{d}\Gamma(\zeta)^{2}}\left(\frac{J!\,\Gamma\left(\nu^{(\zeta)}\right)}{4^{J}\,\Gamma\left(\nu^{(\zeta)}+J\right)}\right)^{2}\,\frac{\pi
2^{1-4\nu^{(\zeta)}+J}\Gamma(J+2\nu^{(\zeta)})}{J!(J+\nu^{(\zeta)})\Gamma(\nu^{(\zeta)})^{2}}\\\
&=\frac{2}{N}\frac{\pi^{\frac{1}{2}-d}J!2^{-d-2\zeta-J+2}\Gamma(d+J-2\zeta-1)\Gamma\left(\frac{d}{2}+J-\zeta\right)}{\Gamma(\zeta)^{2}\Gamma\left(\frac{d}{2}+J-\zeta-\frac{1}{2}\right)}\;.\end{split}$
(A.10)

This result agrees in the limit $\zeta=1$ with [51] up to the
$\left(\frac{J!\,\Gamma\left(\nu^{(1)}\right)}{4^{J}\,\Gamma\left(\nu^{(1)}+J\right)}\right)^{2}$
factor from the normalization of [50].

###  A.2 OPE coefficient $f_{\phi\phi J}$

Now we compute the OPE coefficient between a spinning operators of the kind
$\mathcal{O}_{0,J}$ and two fundamental fields. The functional form of the
three-point function is constrained by conformal symmetry to be

$\langle\mathcal{O}_{0,J}(x_{0},\xi)\phi_{i}(x_{1})\phi_{j}(x_{2})\rangle=\delta_{ij}f_{\phi\phi
J}\left(\frac{\xi\cdot x_{20}}{x_{20}^{2}}-\frac{\xi\cdot
x_{10}}{x_{10}^{2}}\right)^{J}\frac{1}{(x_{02}^{2})^{\frac{d}{2}-\zeta}(x_{01}^{2})^{\frac{d}{2}-\zeta}}\;.$
(A.11)

We can compute the three-point function explicitly by Wick theorem as we did
for the two-point function and extract the OPE coefficient by comparing with
(A.11). We have

$\begin{split}\langle\mathcal{O}_{0,J}(x_{0},\xi)\phi_{i}(x_{1})\phi_{i}(x_{2})\rangle=&\frac{2}{N}\left(\frac{1}{2^{2\zeta}\pi^{d/2}\Gamma(\zeta)}\right)^{2}\int_{0}^{+\infty}d\alpha_{1}\,\alpha_{1}^{d/2-\zeta-1}\int_{0}^{+\infty}d\alpha_{2}\,\alpha_{2}^{d/2-\zeta-1}\\\
&\qquad\times
f_{J}(-2\,\alpha_{1}\,\xi\cdot(x_{0}-x_{1}),-2\,\alpha_{2}\,\xi\cdot(x_{0}-x_{2}))e^{-\alpha_{1}x_{01}^{2}-\alpha_{2}x_{02}^{2}}\;.\end{split}$
(A.12)

Without loss of generality we can set $x_{0}=0$ and we get to

$\begin{split}&\langle\mathcal{O}_{0,J}(x_{0},\xi)\phi_{i}(x_{1})\phi_{i}(x_{2})\rangle=\frac{2}{N}\frac{J!\,\Gamma(\nu^{(\zeta)})}{4^{2\zeta+J}\pi^{d}\Gamma(\zeta)^{2}\Gamma(\nu^{(\zeta)}+J)}\int_{0}^{+\infty}d\alpha_{1}\,\alpha_{1}^{d/2-\zeta-1}\int_{0}^{+\infty}d\alpha_{2}\,\alpha_{2}^{d/2-\zeta-1}\\\
&\qquad\qquad\times\left(2\alpha_{1}\frac{\xi\cdot
x_{1}}{x_{1}^{2}}+2\alpha_{2}\frac{\xi\cdot
x_{2}}{x_{2}^{2}}\right)^{J}C^{(\nu^{(\zeta)})}_{J}\left(\frac{2\alpha_{1}\frac{\xi\cdot
x_{1}}{x_{1}^{2}}+2\alpha_{2}\frac{\xi\cdot
x_{2}}{x_{2}^{2}}}{2\alpha_{1}\frac{\xi\cdot
x_{1}}{x_{1}^{2}}-2\alpha_{2}\frac{\xi\cdot
x_{2}}{x_{2}^{2}}}\right)\frac{e^{-\alpha_{1}-\alpha_{2}}}{(x_{1}^{2})^{d/2-\zeta}(x_{2}^{2})^{d/2-\zeta}}\;.\end{split}$
(A.13)

Considering the region $x_{2}^{2}\gg x_{1}^{2}$ the argument the argument of
the Gegenbauer polynomials become $1$ and it is easy to extract

$\begin{split}f_{\phi\phi
J}&=\frac{2}{N}\frac{J!\,\Gamma(\nu^{(\zeta)})}{4^{2\zeta+J}\pi^{d}\Gamma(\zeta)^{2}\Gamma(\nu^{(\zeta)}+J)}2^{J}\Gamma\left(\frac{d}{2}-\zeta\right)\Gamma\left(\frac{d}{2}-\zeta+J\right)\frac{\Gamma(2\nu^{(\zeta)}+J)}{\Gamma(2\nu^{(\zeta)})J!}\\\
&=\frac{2}{N}\frac{\Gamma\left(\frac{d-1}{2}-\zeta\right)\Gamma\left(\frac{d}{2}-\zeta\right)\Gamma(d+J-2\zeta-1)\Gamma\left(\frac{d}{2}+J-\zeta\right)}{\pi^{d}2^{4\zeta+J}\Gamma(\zeta)^{2}\Gamma(d-2\zeta-1)\Gamma\left(J+\frac{d-1}{2}-\zeta\right)}\;.\end{split}$
(A.14)

Now it is simple to check by a direct computation that in our conventions

$\frac{f_{\phi\phi J}}{c_{0,J}}=\frac{1}{J!}\;.$ (A.15)

## References

  * [1] L. Iliesiu, M. Koloğlu, R. Mahajan, E. Perlmutter and D. Simmons-Duffin, _The Conformal Bootstrap at Finite Temperature_ , _JHEP_ 10 (2018) 070 [1802.10266].
  * [2] N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo and A. Trombettoni, _Long-range interacting quantum systems_ , 2109.01063.
  * [3] F.J. Dyson, _Existence of a phase transition in a one-dimensional Ising ferromagnet_ , _Commun. Math. Phys._ 12 (1969) 91.
  * [4] M.E. Fisher, S.-k. Ma and B. Nickel, _Critical exponents for long-range interactions_ , _Phys.Rev.Lett._ 29 (1972) 917.
  * [5] Y. Yamazaki and M. Suzuki, _Critical behavior of isotropic systems with long range interactions_ , _Prog. Theor. Phys._ 57 (1977) 1886.
  * [6] D. Benedetti, R. Gurau, S. Harribey and K. Suzuki, _Long-range multi-scalar models at three loops_ , _J. Phys. A_ 53 (2020) 445008 [2007.04603].
  * [7] J. Sak, _Recursion relations and fixed points for ferromagnets with long-range interactions_ , _Phys. Rev. B_ 8 (1973) 281.
  * [8] E. Brezin, G. Parisi and F. Ricci-Tersenghi, _The crossover region between long-range and short-range interactions for the critical exponents_ , _J. Stat. Phys._ 157 (2014) 855 [1407.3358].
  * [9] C. Behan, L. Rastelli, S. Rychkov and B. Zan, _Long-range critical exponents near the short-range crossover_ , _Phys. Rev. Lett._ 118 (2017) 241601 [1703.03430].
  * [10] C. Behan, L. Rastelli, S. Rychkov and B. Zan, _A scaling theory for the long-range to short-range crossover and an infrared duality_ , _J.Phys.A_ 50 (2017) 354002 [1703.05325].
  * [11] D.C. Brydges, P.K. Mitter and B. Scoppola, _Critical (Phi**4)(3, epsilon)_ , _Commun. Math. Phys._ 240 (2003) 281 [hep-th/0206040].
  * [12] A. Abdesselam, _A complete renormalization group trajectory between two fixed points_ , _Commun.Math.Phys._ 276 (2007) 727 [math-ph/0610018].
  * [13] G. Slade, _Critical exponents for long-range ${O(n)}$ models below the upper critical dimension_, _Communications in Mathematical Physics_ 358 (2017) 343?436 [1611.06169].
  * [14] M. Lohmann, G. Slade and B.C. Wallace, _Critical two-point function for long-range o(n) models below the upper critical dimension_ , _Journal of Statistical Physics_ 169 (2017) 1132 [1705.08540].
  * [15] A. Giuliani, V. Mastropietro and S. Rychkov, _Gentle introduction to rigorous Renormalization Group: a worked fermionic example_ , _JHEP_ 01 (2021) 026 [2008.04361].
  * [16] M.F. Paulos, S. Rychkov, B.C. van Rees and B. Zan, _Conformal Invariance in the Long-Range Ising Model_ , _Nucl. Phys. B_ 902 (2016) 246 [1509.00008].
  * [17] S. Giombi and H. Khanchandani, _$O(N)$ models with boundary interactions and their long range generalizations_, _JHEP_ 08 (2020) 010 [1912.08169].
  * [18] J. Cardy, ed., _Finite-size Scaling_ , Current physics, North-Holland (1988).
  * [19] C. Linhares, A. Malbouisson and I. Roditi, _Ginzburg–Landau Theory of Phase Transitions in Compactified Spaces_ , _Advances in Quantum Field Theory, Prof. Sergey Ketov (Ed.), InTech_ (2012) .
  * [20] M. Le Bellac, _Thermal Field Theory_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (2011), 10.1017/CBO9780511721700.
  * [21] A. Dutta and J.K. Bhattacharjee, _Phase transitions in the quantum ising and rotor models with a long-range interaction_ , _Phys. Rev. B_ 64 (2001) 184106.
  * [22] N. Defenu, A. Trombettoni and S. Ruffo, _Criticality and phase diagram of quantum long-range o(n) models_ , _Physical Review B_ 96 (2017)  [1704.00528].
  * [23] E. Gonzalez-Lazo, M. Heyl, M. Dalmonte and A. Angelone, _Finite-temperature critical behavior of long-range quantum Ising models_ , _SciPost Phys._ 11 (2021) 076 [2104.15070].
  * [24] J. Zhao, M. Song, Y. Qi, J. Rong and Z.Y. Meng, _Finite-temperature critical behaviors in 2d long-range quantum heisenberg model_ , 2306.01044.
  * [25] M. Suzuki, _Relationship between d-Dimensional Quantal Spin Systems and (d+1)-Dimensional Ising Systems: Equivalence, Critical Exponents and Systematic Approximants of the Partition Function and Spin Correlations_ , _Progress of Theoretical Physics_ 56 (1976) 1454 [https://academic.oup.com/ptp/article-pdf/56/5/1454/5264429/56-5-1454.pdf].
  * [26] S. Sachdev, _Quantum Phase Transitions_ , Cambridge University Press (4, 2011), 10.1017/cbo9780511973765.
  * [27] M. Henkel and M. Pleimling, _Non-equilibrium phase transitions. Vol. 2: ageing and dynamical scaling far from equilibrium_ , Theoretical and Mathematical Physics, Springer Dordrecht (2010), 10.1007/978-90-481-2869-3.
  * [28] A.C. Petkou and N.D. Vlachos, _Finite size effects and operator product expansions in a CFT for d $>$ 2_, _Phys. Lett. B_ 446 (1999) 306 [hep-th/9803149].
  * [29] A.C. Petkou and N.D. Vlachos, _Finite size and finite temperature effects in the conformally invariant O(N) vector model for 2 less than d less than 4_ , in _5th International Workshop on Thermal Field Theories and Their Applications_ , 9, 1998 [hep-th/9809096].
  * [30] A.C. Petkou and A. Stergiou, _Dynamics of Finite-Temperature Conformal Field Theories from Operator Product Expansion Inversion Formulas_ , _Phys. Rev. Lett._ 121 (2018) 071602 [1806.02340].
  * [31] J.R. David and S. Kumar, _Thermal one-point functions: CFT’s with fermions, large d and large spin_ , _JHEP_ 10 (2023) 143 [2307.14847].
  * [32] O. Diatlyk, F.K. Popov and Y. Wang, _Beyond $N=\infty$ in Large $N$ Conformal Vector Models at Finite Temperature_, 2309.02347.
  * [33] F.A. Dolan and H. Osborn, _Conformal Partial Waves: Further Mathematical Results_ , 1108.6194.
  * [34] M.S. Costa, T. Hansen, J. Penedones and E. Trevisani, _Projectors and seed conformal blocks for traceless mixed-symmetry tensors_ , _JHEP_ 07 (2016) 018 [1603.05551].
  * [35] P. Chaikin and T. Lubensky, _Principles of Condensed Matter Physics_ , Cambridge University Press (2000).
  * [36] P.C. Hohenberg and B.I. Halperin, _Theory of Dynamic Critical Phenomena_ , _Rev. Mod. Phys._ 49 (1977) 435.
  * [37] E. Ardonne, P. Fendley and E. Fradkin, _Topological order and conformal quantum critical points_ , _Annals Phys._ 310 (2004) 493 [cond-mat/0311466].
  * [38] D. Anselmi and M. Halat, _Renormalization of Lorentz violating theories_ , _Phys. Rev. D_ 76 (2007) 125011 [0707.2480].
  * [39] D. Karateev, P. Kravchuk and D. Simmons-Duffin, _Harmonic Analysis and Mean Field Theory_ , _JHEP_ 10 (2019) 217 [1809.05111].
  * [40] A. Lodhia, S. Sheffield, X. Sun and S.S. Watson, _Fractional Gaussian fields: A survey_ , _Probability Surveys_ 13 (2016) 1  [1407.5598].
  * [41] A. Campa, T. Dauxois and S. Ruffo, _Statistical mechanics and dynamics of solvable models with long-range interactions_ , _Physics Reports_ 480 (2009) 57 [0907.0323].
  * [42] M. Kwasnicki, _Ten equivalent definitions of the fractional Laplace operator_ , _Fractional Calculus and Applied Analysis_ 20 (2017) 51 [1507.07356].
  * [43] D. Benedetti, R. Gurau, S. Harribey and D. Lettera, _The F-theorem in the melonic limit_ , _JHEP_ 02 (2022) 147 [2111.11792].
  * [44] Z. Komargodski and A. Zhiboedov, _Convexity and Liberation at Large Spin_ , _JHEP_ 11 (2013) 140 [1212.4103].
  * [45] V.K. Dobrev, V.B. Petkova, S.G. Petrova and I.T. Todorov, _Dynamical Derivation of Vacuum Operator Product Expansion in Euclidean Conformal Quantum Field Theory_ , _Phys. Rev. D_ 13 (1976) 887.
  * [46] D. Simmons-Duffin, D. Stanford and E. Witten, _A spacetime derivation of the lorentzian OPE inversion formula_ , _Journal of High Energy Physics_ 2018 (2018) .
  * [47] J. Liu, E. Perlmutter, V. Rosenhaus and D. Simmons-Duffin, _$d$ -dimensional SYK, AdS Loops, and $6j$ Symbols_, _JHEP_ 03 (2019) 052 [1808.00612].