#  Exact solution for quantum strong long-range models via a generalized
Hubbard-Stratonovich transformation

Juan Román-Roche  Instituto de Nanociencia y Materiales de Aragón (INMA),
CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain  Departamento de Física de
la Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain  Víctor
Herráiz-López  Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-
Universidad de Zaragoza, Zaragoza 50009, Spain  Departamento de Física de la
Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain  David
Zueco  Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-
Universidad de Zaragoza, Zaragoza 50009, Spain  Departamento de Física de la
Materia Condensada, Universidad de Zaragoza, Zaragoza 50009, Spain

###### Abstract

We present an exact analytical solution for quantum _strong_ long-range models
in the canonical ensemble by extending the classical solution proposed in
[Campa _et al._ , J. Phys. A 36, 6897 (2003)]. Specifically, we utilize the
equivalence between generalized Dicke models and interacting quantum models as
a generalization of the Hubbard-Stratonovich transformation. To demonstrate
our method, we apply it to the Ising chain in transverse field and discuss its
potential application to other models, such as the Fermi-Hubbard model,
combined short- and long-range models and models with antiferromagnetic
interactions. Our findings indicate that the critical behaviour of a model is
independent of the range of interactions, within the strong long-range regime,
and the dimensionality of the model. Moreover, we show that the order
parameter expression is equivalent to that provided by mean-field theory, thus
confirming the exactness of the latter. Finally, we examine the algebraic
decay of correlations and characterize its dependence on the range of
interactions in the full phase diagram.

##  I Introduction

Long-range systems are those in which two-body interactions decay as a power-
law at large distances. They are ubiquitous in nature, with some examples
given by dipolar, Coulomb or Wan-der-Walls interactions. Recent experimental
advances in atomic, molecular and optical systems have lead to a resurgence of
interest in long-range models [1, 2, 3, 4]. In these experiments, the
effective interactions between spins are often long-ranged and tunable,
renewing the need for a comprehensive understanding of long-range systems.
Although less studied than their short-ranged counterparts, there are already
some rigorous and numerical results available [5, 6, 7, 8, 9]. Some
equilibrium and dynamical properties have been discussed in comparison with
short-range systems. Notable examples are the existence (or absence) of an
area law of entanglement [10, 11, 12, 13], the algebraic decay of two-point
correlators out of criticality [14, 15, 16], the spreading of correlations
[17], the existence of Majorana modes [18] and topological properties [19].

In these examples, the phenomenology can be understood within a
(sub)classification in terms of the range of interactions they exhibit. To fix
notation and ideas, let us introduce this classification with the models
considered in this paper: quantum long-range models in an $N$-site lattice
with a coupling of the form

$\mathcal{H}_{\rm c}=-\sum_{ij}^{N}J_{ij}\mathcal{C}_{i}\mathcal{C}_{j}\,,$
(1)

where $\mathcal{C}_{i}$ is a local hermitian operator acting on site $i$. We
consider models with power-law decaying interactions
$J_{ij}=\Gamma\tilde{J}(\mathbf{r}_{ij})/\tilde{N}$,

$\tilde{J}(\mathbf{r}_{ij})=\begin{cases}b&\text{ if
}\quad\mathbf{r}_{ij}=0\\\ |\mathbf{r}_{ij}|^{-\alpha}&\text{
otherwise}\end{cases}$ (2)

and periodic boundary conditions (PBC).

Figure 1: Classification of long-range models following reviews [5] and [8]
valid for both the classical and quantum models. This work presents a solution
for the quantum strong long-range regime.

The distance between sites $\mathbf{r}_{ij}$ is then given by the nearest
image convention. Through this work we will focus on the case of attractive or
ferromagnetic interaction, so the interaction strength is $\Gamma>0$, although
the extension to antiferromagnetic or repulsive models will be discussed. $b$
is a parameter that can be tuned to shift the spectrum of $J$. The decay rate,
$\alpha$, sets the range of the interactions. For $\alpha<d$, where $d$ is the
dimensionality of the lattice, the interactions decay slowly enough that the
sum in the coupling term (1) depends superlinearly on $N$, breaking the
extensivity of the model (See App. A.1.). Kac’s renormalization factor
$1/\tilde{N}$ restores extensivity, ensuring a well-defined thermodynamic
limit. Here $\tilde{N}=\sum_{i}\tilde{J}_{ij}$, note that PBC make the model
translation invariant and thus $\sum_{i}\tilde{J}_{ij}$ is independent of $j$.
Regardless, the model remains non-additive in this regime. Non-additivity
brings about particular statistical and dynamical phenomena that differ from
the commonly studied short-range models, such as ensemble inequivalence,
negative specific heat and quasistationary states [8]. Accordingly, the regime
$\alpha<d$ is identified as _(strong) long-range_. In the regime $\alpha>d$,
the model is naturally extensive, and Kac’s renomarlization factor amounts to
a rescaling of the interaction strength. Within the regime $\alpha>d$ two
further subregimes can be identified: for $\alpha>\alpha_{*}$ the critical
exponents of the model match those of the nearest-neighbours model
($\alpha\to\infty$), this is the _short-range regime_ ; for
$d<\alpha<\alpha_{*}$ the model presents critical exponents that differ from
the short-range ones, the effects of long-range interactions are felt but the
model is additive, this is the _weak long-range regime_ [5, 8]. For
convenience, we summarize this classification in Fig. 1. _This work deals with
the strong long-range regime_.

Strong long-range models are commonly disregarded in many analytical and
numerical studies on the grounds of the ill-defined thermodynamic limit
brought about by the non-extensivity. Kac’s rescaling eliminates this barrier,
making their study possible. For quantum models, seminal numerical studies
solving the transverse-field Ising model in the strong long-range regime are
found in Refs. [20, 21]. They confirm that in this regime the model is within
the mean-field universality class. This is in agreement with the claim that
mean field is exact for quantum spin models in the strong long-range regime
[22] that generalizes similar findings in classical systems [23, 24, 25, 26].
These works are crucial for the rigorous understanding of the physics of long-
range systems. On the one hand, they provide an exact way to solve them, on
the other hand, they provide a starting point for approximations that tackle
the weak long-range regime.

This work provides a recipe to analytically solve, in the canonical ensemble,
quantum strong long-range models. Therefore, it complements the work of Mori
[22] and confirms that _in the strong long-range regime mean field is exact_.
Besides, it extends the work of Campa and coworkers for classical strong long-
range models to the quantum case [27]. Our work introduces a generalized
Hubbard-Stratonovich transformation (HST) and provides a closed expression for
the free energy at any temperature. Technically, we show how to use the
equivalence between generalized Dicke models and interacting quantum models as
a quantum HST. We show that only _strong_ long-range models admit this mapping
and formulate their canonical solution in terms of the associated Dicke model,
which is then tackled following the prescription of Wang and Hioe [28, 29]. We
illustrate the method on the Ising chain in transverse field. We find that the
critical behaviour is universal for $\alpha<d$ and any lattice dimensionality.
The expression for the magnetization (the order parameter) is shown to be
equivalent to the mean-field solution, thus proving the exactness of the
latter. Finally, we study the algebraic decay of correlations as a function of
the decay rate of interactions $\alpha$.

The rest of the paper is organized as follows. In Section II, we provide a
brief overview of the HST as a tool to solve classical models, which forms the
basis for our further development. In Section III, we establish the
relationship between generalized Dicke and long-range models and introduce the
generalized HST. Section IV presents the solution for strong long-range models
and a discussion of which models can be treated with this method. We perform
the calculations for the long-range transverse field Ising model, including
the full phase diagram and the decay of two-point correlations in section V.
Finally, we conclude the article with some general remarks and relegate more
technical details to the Appendices.

##  II Sketch of the solution for classical systems

To warm up, it is convenient to understand how to solve classical strong long
range models, mainly following the works of Campa and coworkers [23, 27]. For
simplicity consider the Ising model,

$\mathcal{H}_{\rm cl}=h\sum_{i}^{N}s_{i}-\sum_{ij}^{N}J_{ij}s_{i}s_{j}\;.$ (3)

Here, $s_{i}$ is a discrete variable. The solution is based on two main
observations. First, diagonalizing the interaction matrix $J=\Lambda
D\Lambda^{T}$, which allows one to write the coupling term as
$\sum_{k}D_{k}\left(\sum_{i}\Lambda_{ik}s_{i}\right)^{2}$, where $\\{D_{k}\\}$
are the eigenvalues of the interaction matrix. Note that, in this form, the
coupling is written as a sum of interaction terms that are quadratic in
$\sum_{i}^{N}\Lambda_{ik}s_{i}$. Second, eliminating these quadratic
interactions by use of the Hubbard-Stratonovich transformation which is based
on the equality,

$\begin{split}Z&=\sum_{s_{i}}e^{-\beta\mathcal{H}_{cl}}\\\ &\propto\int
du_{k}\sum_{s_{i}}e^{-\beta\big{(}h\sum_{i}s_{i}+\sum_{k}u_{k}^{2}/D_{k}-2\sum_{ik}\Lambda_{ik}s_{i}u_{k}\big{)}}\end{split}$
(4)

here, $u_{k}$ are real auxiliary variables. This equality for the partition
function follows from Gaussian integral formulas.

Notice that we have decoupled the interaction $J_{ij}s_{i}s_{j}$, therefore
the sum over $s_{i}$-configurations are trivial. Finally, the integral over
the real variables $u_{k}$ can be done within the saddle point approximation
in the $N$-large limit. This is true if some conditions are met on the
eigenspectrum of the $J_{ij}$ matrix, see below and the discussion in App. B.

##  III From the Dicke model to quantum long-range models and back

###  III.1 Effective theory of the Dicke model

Figure 2: Analysis of the eigenvalues of the coupling matrix $J$ (1) for
$d=1$. Left: For $\alpha=0.2<1$, a plot of the eigenvalues for $N=100$ on top
and a histogram of the eigenvalues as a function of $N$ on the bottom. Right:
Same but for $\alpha=1.8>1$.

The method described above and utilized in Ref. [27] cannot be
straightforwardly applied to quantum models. The application of the HST
requires splitting the exponential that constitutes the kernel of the
partition function into a product of exponentials, which in the quantum case
is prevented by the non-commutativity of the long-range interaction term and
other terms in the Hamiltonian. There are ways in which the HST can be applied
to solve quantum systems, but it requires a reframing of the partition
function in terms of commuting quantities. A field-theory formulation or
imaginary-time Trotterization are examples of this. Here we present an
alternative which is the closest to the classical formulation.

Our method utilizes some results from quantum optics in order to draw an
equivalence between some quantum long-range models and a cavity QED model.
Specifically, we utilize the generalized Dicke model as our starting point to
develop this equivalence:

$\mathcal{H}=\sum_{k=0}^{M-1}\omega_{k}a_{k}^{\dagger}a_{k}+\mathcal{H}_{0}-\sum_{k,i}\left(a_{k}+a_{k}^{\dagger}\right)\frac{\lambda_{ik}}{\sqrt{N}}\mathcal{C}_{i}\,.$
(5)

Here $\mathcal{H}_{0}$ is an exactly solvable Hamiltonian of the “mater”
degrees of freedom and $\mathcal{C}_{i}$ is the local hermitian operator that
couples site $i$ to the bosonic modes with
$[a_{k},a_{k^{\prime}}^{\dagger}]=\delta_{k,k^{\prime}}$, finally
$\lambda_{ik}$ are real coupling constants. In two previous publications [30,
31] we show that the generalized Dicke model, after integrating out the
electromagnetic modes, yields an exact effective Hamiltonian description for
the matter degrees of freedom alone in the limit $N\to\infty$ (thermodynamic
limit). The resulting Hamiltonian,

$\mathcal{H}_{\rm m}^{\rm
eff}=\mathcal{H}_{0}-\sum_{ij}^{N}\sum_{k=0}^{M-1}\frac{\lambda_{ik}\lambda_{jk}}{N\omega_{k}}\mathcal{C}_{i}\mathcal{C}_{j}\,,$
(6)

corresponds to a quantum model with interactions given by $(J_{\rm
eff})_{ij}=\sum_{k}^{M-1}\lambda_{ik}\lambda_{jk}/(N\omega_{k})$. The mode
structure of the cavity determines the resulting effective model. However, it
is important to note that the exact mapping between Hamiltonians (5) and (6)
is limited to the thermodynamic limit, $N\to\infty$, and a number of modes $M$
such that $\lim_{N\to\infty}M/N=0$. Below, we demonstrate how we can reverse
the effective theory to solve a quantum model. The first question that arises
is which family of quantum models, with interaction given by Eq. (1), can be
solved this way, _i.e._ which can be cast in the form of Eq. (6). Below, we
show that this is the case for strong long range models, $\alpha<d$, this is
the first result of this paper.

###  III.2 Mapping a quantum model to the Dicke model

If we start from an arbitrary extensive 111meaning that Kac’s prescription is
used to ensure extensivity if the model is strong-long-ranged model of the
form $\mathcal{H}_{\rm
m}=\mathcal{H}_{0}-\sum_{ij}J_{ij}\mathcal{C}_{i}\mathcal{C}_{j}$, Cf. Eq.
(1), the first step is to diagonalize the interaction matrix $J=\Lambda
D\Lambda^{T}$, where $D$ is a diagonal matrix, $D_{kp}\equiv
D_{k}\delta_{kp}$. Note that $\Lambda$ is orthogonal because $J$ is symmetric.
The matrix elements are then given by

$J_{ij}=\sum_{k=0}^{N-1}\Lambda_{ik}D_{k}\Lambda_{jk}$ (7)

Assuming that $J_{ij}>0$, the smallest eigenvalue of $J$ can always be set to
zero by adjusting its diagonal elements, which we denote $b$. Fixing $b\neq 0$
introduces, generally, non-trivial diagonal terms of the form $\Gamma
b/\tilde{N}\mathcal{C}_{i}^{2}$. These can be shown to be negligible in the
thermodynamic limit, so the freedom to set $b$ remains (See App. A.2.). For a
general interaction matrix the number of non zero eigenvalues, $M$, scales
with the size of the matrix, $N$. Conveniently, it can be shown that for a
model with power-law decaying interactions and PBC such as the one considered
in this work (1), the number of non-zero modes in the thermodynamic limit
($N\to\infty$) depends on the decay rate of the interaction [27]. For a model
in the strong long-range regime, $\alpha<d$, only a small number of modes have
a non zero eigenvalue, such that $\lim_{N\to\infty}M/N=0$. This can be seen
analytically in models with a translation-invariant interaction matrix, which
can be diagonalized in Fourier space, obtaining a closed expression for its
eigenvalues:

$D(\mathbf{q})=\frac{\Gamma}{\tilde{N}}\sum_{\mathbf{r}}\tilde{J}(\mathbf{r})\exp[-i\mathbf{q}\mathbf{r}]\,.$
(8)

Here $\mathbf{q}$ denotes any of the $N$ reciprocal-space vectors in the first
Brillouin zone and the sum runs over all lattice points. The large-$N$
behaviour of $D(\mathbf{q})$ can then be estimated by replacing the sum with
an integral [27].

Figure 3: Schematic of the generalized Hubbard-Stratonovich transformation
mapping a quantum long-range model to (and from) a generalized Dicke model.
Blue dots represent “matter” degrees of freedom and red squares represent
cavity modes. Cf. Hamiltonians (5) and (6) in the main text.

Complementarily, we provide in Fig. 2 a graphical analysis of this phenomenon
by showing the typical distribution of eigenvalues depending on $\alpha$ for
$d=1$ ( the same behaviour is observed in other dimensions, not shown). This
graphical analysis can be useful for models without translation invariance. In
Fig. 2 we show that for a strong long-range model the eigenvalues bunch around
zero as $N$ increases, whereas they remain more uniformly distributed in the
weak long-range regime. This can be condensed into a criterion for determining
whether arbitrary models are tractable: knowing that the eigenvalues of $J$
are non-negative and bounded by construction, if only a vanishingly small
fraction $M/N$ are non-zero for $N\to\infty$, then their average will tend to
zero and vice versa. Thus, for an arbitrary interaction matrix $J$, if

$\lim_{N\to\infty}\frac{1}{N}\sum_{k=0}^{M-1}D_{k}=\lim_{N\to\infty}\frac{1}{N}{\rm
Tr}(J)=0$ (9)

the model is tractable, i.e. the number $M$ of non-zero eigenvalues scales as
$\lim_{N\to\infty}M/N=0$. If we apply this criterion to translation invariant
models we find $\lim_{N\to\infty}1/N\sum_{k}D_{k}=\lim_{N\to\infty}\Gamma
b/\tilde{N}$, which is zero for $\alpha<1$ and non-zero otherwise (See Apps.
A.1 and A.2.).

Once it is established that a given model has a sufficiently small number of
non-zero eigenvalues, one can sort them by decreasing value and truncate the
sum in Eq. (7) to consider only the first $M$ terms for which $D_{k}\neq 0$.
For these remaining non-zero eigenvalues, we can identify
$\omega_{k}=1/D_{k}$. This, together with the rescaled elements of the change
of basis matrix, $\lambda_{ik}=\sqrt{N}\Lambda_{ik}$, leads to Eq. (6) and
effectively defines the mode structure for the associated Dicke model (5).

In summary, models for which $\lim_{N\to\infty}M/N=0$, in particular strong
long-range models, can be mapped to generalized Dicke models via the effective
theory described in [30, 31]. To gain further insight, it is useful to compare
Hamiltonians (5) and (6) with the system depicted in Figure 3, which outlines
our procedure. The interacting model, shown in the left-hand side of the
figure, where constituents are depicted as blue nodes and interactions as
black edges, is mapped to a larger system where the _physical_ (matter)
degrees of freedom are uncoupled and interact with _auxiliary_ bosonic modes
represented as red squares 222Fig. 3 is actually an oversimplification, as it
only depicts long-range interactions, which are the ones replaced by the
auxiliary bosonic modes of the effective theory. The theory is also applicable
to models containing a combination of short- and long-range interactions. For
a discussion about the applicability of the method see Sec. IV. Integrating
out these bosonic modes would lead back to the desired interactions $J_{ij}$.
The parallelism between the auxiliary bosonic modes in the effective theory
and the auxiliary classical fields in the standard HST motivates _the claim of
a generalized Hubbard-Stratonovich transformation_.

##  IV Exact solution of strong long-range models

At this point we have shown how to map a strong long-range quantum model (1)
to a generalized Dicke model (5) as illustrated in Fig. 3. To solve the
latter, we will follow the steps outlined in the original solution of the
Dicke model by Wang and Hioe [28, 29]. In the thermodynamic limit, the trace
over the photonic degrees of freedom is replaced by a collection of complex
Gaussian integrals and the bosonic creation and annihilation operators,
$a_{k}^{\dagger}$, $a_{k}$, are replaced by complex fields, $\alpha_{k}^{*}$,
$\alpha_{k}$

$Z=\int\prod_{k=0}^{M-1}\frac{d^{2}\alpha_{k}}{\pi}{\rm Tr}_{\rm
m}\left\\{\exp\left[-\beta\left(\sum_{k=0}^{M-1}\omega_{k}|\alpha_{k}|^{2}+\mathcal{H}_{0}+\sum_{k,i}\frac{2\lambda_{ik}x_{k}}{\sqrt{N}}\mathcal{C}_{i}\right)\right]\right\\}\,,$
(10)

where $\alpha_{k}=x_{k}+iy_{k}$. At this point, the parallelism with the
standard Hubbard-Stratonovich transformation used in the classical model is
even more explicit. The Gaussian integral over the imaginary parts
$\\{y_{k}\\}$ yields an unimportant constant. To tackle the integration over
the real parts, we perform a change of variables $u_{k}^{2}=x_{k}^{2}/N$ and
define

$\begin{split}Z_{\rm m}[u_{k}]&\equiv Z_{\rm m}(u_{0},\ldots,u_{M-1})\\\
&={\rm Tr}_{\rm
m}\left\\{e^{-\beta\left(\mathcal{H}_{0}+\sum_{k,i}2\lambda_{ik}u_{k}\mathcal{C}_{i}\right)}\right\\}\end{split}$
(11)

and $f_{\rm m}[u_{k}]=\ln(Z_{\rm m}[u_{k}])/N$. In the resulting integral

$Z=\int\prod_{k=0}^{M-1}\sqrt{\frac{N}{\pi\omega_{k}}}du_{k}\exp\left(N\phi[u_{k}]\right)\,,$
(12)

where

$\phi[u_{k}]=-\beta\sum_{k=0}^{M-1}\omega_{k}u_{k}^{2}+f_{\rm m}[u_{k}]\,,$
(13)

the exponent depends explicitly linearly on $N$, allowing one to use the
saddle-point method (exactly for $N\to\infty$) to express the partition
function as the value of the integrand at the maximum of $\phi[u_{k}]$

$Z=\prod_{k=0}^{M-1}\sqrt{\frac{N}{\pi\omega_{k}}}\exp\left(N\phi[\bar{u}_{k}]\right)\,,$
(14)

with,

$\phi[\bar{u}_{k}]=\max_{\\{u_{k}\\}}\phi[u_{k}]\,.$ (15)

Computing the partition function is thus reduced to a multivariate
maximization problem. In order for the zero-order saddle-point approximation
to be exact, one has to verify that there exists a maximum
$\\{\bar{u}_{k}\\}$, i.e. that $\phi$ admits a stationary point
$\\{\bar{u}_{k}\\}$ and the eigenvalues of the Hessian of $\phi$ at the
stationary point, $H_{\phi}[\bar{u}_{k}]$, are all negative. In the presence
of several maxima, one has to find the global maximum. Finding global extrema
of a multivariate scalar function is normally a complex task, without
guarantee or provability of success, but in the present case it is greatly
facilitated for homogeneous or near-homogeneous systems (See Sec. V).
Additionally, the second-order corrections to the partition function in the
form of a factor $(\det H_{\phi}[\bar{u}_{k}])^{-1/2}$ must be negligible with
respect to the zero-order term, $\exp\left\\{N\phi[\bar{u}_{k}]\right\\}$, but
this is generally true, see App. B. _This is the main result of this paper,
i.e. the exact expression for the partition function of strong long-range
models ( 14)_.

In the next section and in order to give concrete formulas, we particularize
for the case of the Ising model in transverse field. However, the ideas
presented here can be applied to other models. For instance, our next section
generalizes easily to a spin-$s$ system where $s>1/2$ and also to the
inclusion of a longitudinal field, such that
${\mathcal{H}}_{0}=\omega_{z}\sum_{i}S^{z}_{i}+\omega_{x}\sum_{i}S^{x}_{i}$
and $\mathcal{C}_{i}=2S_{i}^{x}$ with
$[S_{i}^{\alpha},S_{j}^{\beta}]=i\epsilon_{\alpha\beta\gamma}S^{\gamma}\delta_{ij}$
spin-$s$ operators. The Fermi-Hubbard model with long-range interactions could
also be treated with our method. Here
${\mathcal{H}}_{0}=t_{ij}c_{i}^{\dagger}c_{j}$ and
${\mathcal{C}}_{i}=c_{i}^{\dagger}c_{i}$ with
$\\{c_{i},c_{j}^{\dagger}\\}=\delta_{ij}$ fermionic operators. Finally, we
could consider for $\mathcal{H}_{0}$ any model such that
$\mathcal{H}_{0}+\sum_{i}\xi_{i}\mathcal{C}_{i}$ is solvable, where the
$\xi_{i}$ are constants. In doing so, we could combine short-range models (as
the one-dimensional short-range Ising model in transverse field, the XY model,
and so on) with strong long-range interactions. This is because our method
requires knowledge of the eigenstates of
$\mathcal{H}_{0}+\sum_{k,i}2\lambda_{ik}u_{k}\mathcal{C}_{i}$, Cf. Eq. (11).

We have, thus far in the paper, focused only on ferromagnetic (attractive)
models. However, the discussion of the applicability and generalizations of
our method demands that we consider antiferromagnetic (repulsive) models [34,
35]. Frustrated antiferromagnetic long-range models cannot be tackled with our
method. To see why, it suffices to look at Fig. 2, frustrated
antiferromagnetic models arise from changing the global sign of the
interaction in Eq. (1), which in turn results in a change of sign of the
eigenvalues of the coupling matrix. For a general model, a shift to render the
smallest eigenvalue equal to zero is not possible as it would require a $b$ of
the order of $\tilde{N}$, leading to non-vanishing diagonal elements even in
the $N\to\infty$ limit. For a model in which $\mathcal{C}_{i}^{2}=1$, the
shift is possible, but after a shift to render the smallest eigenvalue equal
to zero, we find that the majority of the eigenvalues are non-zero, regardless
of the range of interactions $\alpha$. In contrast, it is possible to define
unfrustrated long-range antiferromagnetic models,
$J_{ij}=\Gamma(-1)^{i+j}\tilde{J}(\boldsymbol{r}_{ij})/\tilde{N}$, as an
extension of unfrustrated nearest-neighbour antiferromagnetic interactions
[36]. Here, the sign change is alternating, rather than global, effectively
defining two sublattices. The interaction matrix defined this way shares the
eigenvalues of its ferromagnetic counterpart and the corresponding model can
thus be tackled with our method. However, a number of interesting subtleties
arise later on in the solution that deserve a detailed discussion, we reserve
this for a future publication.

##  V Solution of the long-range Ising model in a transverse field.

To showcase the effectiveness of the formalism presented in the previous
sections, we particularize now to an Ising chain in transverse field

$\mathcal{H}=\frac{\omega_{z}}{2}\sum_{i}^{N}\sigma^{z}_{i}-\sum_{ij}^{N}J_{ij}\sigma_{i}^{x}\sigma_{j}^{x}\,,$
(16)

where $\sigma^{x,z}$ are the usual Pauli matrices and $J_{ij}$ is given by Eq.
(2). This corresponds to setting $d=1$,
$\mathcal{H}_{0}=\frac{\omega_{z}}{2}\sum_{i}\sigma^{z}_{i}$ and
${\mathcal{C}}_{i}=\sigma^{x}_{i}$. In App. C.1 we show that this leads to

$\phi[u_{k}]=-\beta\sum_{k=0}^{M-1}\omega_{k}u_{k}^{2}+\frac{1}{N}\sum_{i}^{N}\ln\left[2\cosh(\beta\epsilon_{i})\right]\,,$
(17)

with
$2\epsilon_{i}[u_{k}]=\sqrt{\omega_{z}^{2}+4\left(2\sum_{k}\lambda_{ik}u_{k}\right)^{2}}$.

In the case of a homogeneous ${\mathcal{H}}_{0}$ we show in Apps. C.2 and C.3
that the global maximum is homogeneous in the lattice. In terms of the
minimization variables, homogeneity implies that $u_{0}=u\neq 0$ and $u_{k\neq
0}=0$, see App. C.2. This means that only the zero mode, which is constant on
the lattice, $\lambda_{i0}=1\;\forall i$, is relevant in determining the
thermodynamic properties of the model. In turn, one finds that the critical
properties of the model are independent of the decay rate of interactions
$\alpha$, since the latter only determines the degree to which higher-
frequency modes ($k>0$) have to be considered in the diagonalization of $J$.
In more intuitive terms, homogeneity is revealed in the fact that
$2\epsilon_{i}=2\epsilon=\sqrt{\omega_{z}^{2}+16u^{2}},\;\forall i$. In any
case, the multivariate maximization problem simplifies to a single variable
maximization problem $\max_{u}\phi(u)$. Taking the derivative of $\phi$ with
respect to $u$ yields the condition

$\bar{u}\bar{\epsilon}=2\Gamma\bar{u}\tanh\left(\beta\bar{\epsilon}\right)\,,$
(18)

which is manifestly $\alpha$-independent. Note that
$\bar{\epsilon}=\epsilon(\bar{u})$. For $\Gamma<\omega_{z}/4$, $\bar{u}=0$ is
the only solution. For $\Gamma>\omega_{z}/4$, the solution depends on $\beta$,
for $\beta>\beta_{c}$, with $\beta_{c}$ given by
$\omega_{z}=4\Gamma\tanh\left(\beta_{c}\omega_{z}/2\right)$, there is another
solution to Eq. (18) given by
$\bar{\epsilon}=2\Gamma\tanh\left(\beta\bar{\epsilon}\right)$. The solution
$\bar{u}=0$ corresponds to a maximum in the regime where it is the only
solution and becomes a minimum for $\beta>\beta_{c}$ with the maximum given by
the other solution [28]. This marks the paramagnetic-ferromagnetic transition
point. This is the well-known mean-field critical behaviour of the standard
(single-mode homogeneous coupling) Dicke model [37, 28], which is shared by
the Lipkin-Meshkov-Glick (LMG) model (all-to-all homogeneous Ising) [38, 39]
and, as we just showed, is also universal to all strong-long-range Ising
models and their associated Dicke models, _i.e._ we have demonstrated that the
critical point is independent of $\alpha$. This can be visualized in Fig. 4
where the vertical line marks the phase transition, located at the maximum for
the susceptibility (see below), and is independent of $\alpha$. Besides, in
Fig. 5 we compute the critical line, in red, in the $(\Gamma,1/\beta)$-plane
and compare it against the simulations in Ref. [21]. We find excellent
agreement with their numerical results and showcase that the critical point is
independent of $\alpha$ and coincides with the mean-field value.

In terms of observables, we focus now on the calculation of magnetization. In
order to do so from the partition function, we introduce a perturbative
longitudinal field to the Hamiltonian, such that
$\mathcal{H}\to\mathcal{H}-\sum_{i}h_{i}\sigma^{x}_{i}$. Then one can compute
the order parameter $\beta\langle\sigma^{x}_{i}\rangle=\partial\ln Z/\partial
h_{i}$ and the susceptibilities 333The magnetization must be kept
$\\{h_{n}\\}$-dependent in order to compute the susceptibility, the
magnetization of the $\\{h_{n}\\}$-independent model is defined as
$\beta\langle\sigma^{x}_{i}\rangle=\lim_{\\{h_{n}\\}\to 0}\partial\ln
Z/\partial h_{i}$

$\chi_{ij}=\lim_{\\{h_{n}\\}\to
0}\frac{\partial\langle\sigma^{x}_{i}\rangle}{\partial h_{j}}\,.$ (19)

The introduction of longitudinal fields leads to the substitution
$2\sum_{k}\lambda_{ik}u_{k}\to 2\sum_{k}\lambda_{ik}u_{k}+h_{i}$ in Eq. (17).
The magnetization is then

$\langle\sigma^{x}_{i}\rangle=\tanh\left(\beta\epsilon_{i}[\bar{u}_{k}]\right)\frac{2\sum_{k=0}^{M-1}\lambda_{ik}\bar{u}_{k}+h_{i}}{\epsilon_{i}[\bar{u}_{k}]}\,.$
(20)

Here, the magnetization appears as a function of the maximization variables
$\\{\bar{u}_{k}\\}$. However, it is possible to show that
$\sum_{k}\lambda_{ik}\bar{u}_{k}=\sum_{j}J_{ij}\langle\sigma_{j}^{x}\rangle$,
rewriting Eq. (20) as a self-consistent equation on
$\\{\langle\sigma_{i}^{x}\rangle\\}$ which is precisely the self-consistent
equation that arises in a mean-field solution. Our exact analytical method is
thus equivalent to mean-field theory, proving that mean-field theory is exact
for strong long-range models and any lattice dimensionality $d$. Anecdotally,
our theory evidences that the self-consistent solution from mean-field theory
is redundant, in the sense that the solution involves a transcendental
equation of $N$ variables (the magnetizations
$\\{\langle\sigma_{i}^{x}\rangle\\}$), whereas the same problem can be
rewritten in terms of $M$ variables (the $\\{u_{k}\\}$), with
$\lim_{N\to\infty}M/N=0$.

###  V.1 Decay of correlations

Figure 4: Susceptibility for maximally separated spins $\chi_{N/2}$ for strong
long-range quantum Ising models as a function of the interaction strength
$\Gamma$. The parameters used were $\omega_{z}=1$, $\beta\omega_{z}=10$ and
$N=100$.

Our final result concerns the decay of correlations. In weak long-range and
short-range systems correlations decay exponentially at long distances. Only
at the critical point do these systems exhibit power law decay of correlations
[14, 15, 16]. Conversely, strong long-range systems exhibit power law decay of
correlations at all distances. In the absence of exponential decay, the
concept of correlation length cannot be straightforwardly defined, although
there have been some attempts [41]. Here we study the susceptibility
$\chi_{ij}$ (19) as a measure of correlations between spins, as it is
proportional to the Kubo correlator [42][Chap. 4]

$\chi_{ij}=\beta\left(\frac{1}{Z}\int_{0}^{\beta}{\rm
Tr}\left(e^{-(\beta-s)\mathcal{H}}\sigma_{i}^{x}e^{-s\mathcal{H}}\sigma_{j}^{x}\right)-\langle\sigma^{x}_{i}\rangle\langle\sigma^{x}_{j}\rangle\right)\,.$
(21)

The susceptibility can be computed analytically from Eq. (19) for a
translation invariant model [43][Chap. 6] (See. App. D for a derivation) or
numerically otherwise. The analytical derivation yields

$\chi_{ij}=Y\delta_{ij}+\frac{1}{N}\sum_{k=0}^{M-1}\lambda_{ij}\left(\chi_{k}-Y\right)\lambda_{jk}\,,$
(22)

where $Y$ is a quantity that depends on $\bar{u}$ (47) and $\\{\chi_{k}\\}$
are the Fourier modes of the susceptibility (50). Eq. (22) evidences that
$\chi_{i\neq j}$ goes to zero in the thermodynamic limit with a speed that is
determined by the ratio $M/N$ and thus ultimately by $\alpha$ (by its relation
to $d$).

For a numerical calculation, the introduction of a site-dependent field
$h_{i}$ breaks the homogeneity of the model and the multivariate maximization
of $\phi$ is carried out numerically, $\langle\sigma^{x}_{i}\rangle$ is then
computed according to Eq. (20) and $\chi_{ij}$ is computed as a finite
difference. We have verified that both methods yield the same results for the
current model. This is noteworthy because the numerical calculation relies on
a multivariate optimization which could, a priori, converge to an incorrect
result corresponding to a local maxima. We believe the success is due to the
fact that the only deviation from homogeneity stems from the introduction of a
perturbative field and is thus small. Hence, although the optimization is
strictly multivariate, the landscape does not differ much from the univariate
case.

Despite the fact that the analytical results have been obtained under the
assumption that we worked in the thermodynamic limit $N\to\infty$. The
computation of the susceptibility, whether numerically or according to Eq.
(22), requires us to fix a finite value of $N$ and $M$. For each value of
$\alpha$ and $N$, we increase the value of $M$ until convergence is reached
while enforcing the constraint that $\lim_{N\to\infty}M(N)/N=0$.

Because the model is translation invariant, the susceptibility is only a
function of distance, allowing us to define $\chi_{ii+r}\equiv\chi_{r}$. In
Fig. 4 we study the susceptibility at a fixed distance: we plot the half-chain
susceptibility $\chi_{N/2}$ as a function of the interaction strength $\Gamma$
for different decay rates $\alpha<1$ at zero temperature $\beta\to\infty$. The
half-chain susceptibility displays $\alpha$-independent divergence at the
critical point and some dependence on $\alpha$ away from it. Intuitively, the
correlations remain larger for longer-ranged models. We now turn to the
spatial dependence of the susceptibility. In the absence of a correlation
length, we study the susceptibility decay rate, $\alpha_{\chi}$, defined from
the relation $\chi_{r}=A\cdot r^{-\alpha_{\chi}}$. Interestingly, one finds
that $\alpha_{\chi}$ depends linearly on $\alpha$, $\alpha_{\chi}=a\alpha+b$.
In Fig. 5 we plot the slope $a$ as a function of interaction strength $\Gamma$
and inverse temperature $1/\beta$. Close to the critical line, the
susceptibility decay rate $\alpha_{\chi}$ becomes independent of the
interaction decay rate $\alpha$, in agreement with Fig. 4. As one moves
further from the critical line, $a\to 1$, varying continuously from $0$ to $1$
in intermediate regions. In all cases we find $b\approx 0$. Similar algebraic
decays have been described previously for the connected correlator
$\langle\sigma_{1}^{x}\sigma_{r}^{x}\rangle-\langle\sigma_{1}^{x}\rangle\langle\sigma_{r}^{x}\rangle$
in the paramagnetic phase [14]. There, a linear relation between
$\alpha_{\chi}$ and $\alpha$ is also reported. Here we extend those findings
to the full phase diagram.

Figure 5: Phase diagram of the strong long-range Ising chain in transverse
field. The red line corresponds to the universal critical line for $\alpha<1$.
The numerical data is taken from Ref. [21]. The colormap shows the slope of
the linear dependence $\alpha_{\chi}=a\alpha+b$, where $\alpha_{\chi}$ is
susceptibility decay rate and $\alpha$ is the rate of decay of interactions,
computed with parameters $\omega_{z}=1$ and $N=100$.

##  VI Conclusions

In this paper, we have presented a method for solving strong long-range models
in the quantum domain based on the Hubbard-Stratonovich transformation for
classical systems. Our method is a physics-based solution rooted in light-
matter interaction models in which, in the thermodynamic limit, light can be
integrated out leaving an effective long-range model. Solutions of the former,
_i.e._ Dicke models, are due to Hepp and Lieb [44, 45], and Wang and Hioe [28,
29], which we have recently generalized [30, 31].

We have shown that our method can be applied in the strong long-range regime
and confirmed that mean field theory is exact in this regime. Cf. Fig. 1 Eqs.
(2), (5) and (6). In doing so, this paper complements the work of Mori [22].
Besides, it extends the work of Campa and coworkers for classical strong long-
range models to the quantum case [27]. It is worth noting that neither our
method nor the equivalent mean field theory can be used to compute non-local
quantities such as the entanglement entropy, which can be non-trivial in
strong long-range systems [46].

Our method is flexible and could be applied e.g. to spin-$s$ systems where
$s>1/2$, to models with a longitudinal field, such as the long-range XXZ
model, or to the Fermi-Hubbard model with long-range interactions.
Additionally, many exactly solvable models could be complemented with long-
range interactions and solved with our method, since it relies on knowing the
eigenvalues of the system without long-range interactions and with a “field”
term proportional to the long-range coupling operator. Unfrustrated
antiferromagnetic systems are also within the scope of the method. In
conclusion, our work provides a new and powerful tool for solving quantum
long-range models.

## Acknowledgements

We thank Alessandro Campa and Takashi Mori for discussing and explaining
details of their previous results, they were of great help in completing this
work. The authors acknowledge funding from the EU (QUANTERA SUMO and FET-OPEN
Grant 862893 FATMOLS), the Spanish Government Grants
PID2020-115221GB-C41/AEI/10.13039/501100011033 and TED2021-131447B-C21 funded
by MCIN/AEI/10.13039/501100011033 and the EU “NextGenerationEU”/PRTR, the
Gobierno de Aragón (Grant E09-17R Q-MAD) and the CSIC Quantum Technologies
Platform PTI-001. This work has been financially supported by the Ministry of
Economic Affairs and Digital Transformation of the Spanish Government through
the QUANTUM ENIA project call - Quantum Spain project, and by the European
Union through the Recovery, Transformation and Resilience Plan -
NextGenerationEU within the framework of the “Digital Spain 2026 Agenda”. J
R-R acknowledges support from the Ministry of Universities of the Spanish
Government through the grant FPU2020-07231.

##  Appendix A Properties of the long-range interaction matrix

###  A.1 Loss of extensivity

In Fig. 6 we illustrate the extensivity (or lack thereof) of a model with
power-law decaying interactions in $d=1$. We compute
$\tilde{N}=\sum_{j}\tilde{J}_{ij}$ as a measure of the coupling energy per
spin. In the absence of Kac’s rescaling, this quantity must not scale with the
number of spins $N$ to keep the total coupling energy extensive. Fig. 6 shows
that this is not the case for $\alpha<1$. The threshold case $\alpha=1$ is
highlighted for clarity and corresponds to a logarithmic dependence of
$\tilde{N}$ on $N$. For $\alpha>1$ the dependence is sublogarithmic, i.e.
$\tilde{N}$ becomes independent of $N$ at large $N$. As discussed in the main
text, loss of extensivity is prevented with Kac’s rescaling, which we can now
understand as a renormalization of the total coupling energy by the energy per
spin.

In fact, for $d=1$ and $N\to\infty$ the scaling of $\tilde{N}$ can be computed
analytically, since

$\tilde{N}=\sum_{i}^{N}\tilde{J}_{ij}=b+2\sum_{r=1}^{\infty}r^{-\alpha}\;.$
(23)

As $b$ converges to a constant value when $N\to\infty$ (will be shown in App.
A.2), the convergence of $\tilde{N}$ will be ruled by the convergence of the
infinite series. For $\alpha>1$ the series is convergent, so $\tilde{N}$
becomes independent of $N$ at large $N$. For $\alpha\leq 1$ the series
diverges.

###  A.2 Diagonal terms can be neglected in strong-long range models

Setting $b\neq 0$ introduces a new diagonal term in the Hamiltonian

$-\sum_{i}^{N}\Gamma b/\tilde{N}\mathcal{C}_{i}^{2}\,.$ (24)

Importantly, this term contains a factor $b/\tilde{N}$. We know from App. A.1
that $\lim_{N\to\infty}\tilde{N}=\infty$ for $\alpha<d$. So if $b$ is
independent of $N$ for $N\to\infty$, the diagonal term vanishes in the
thermodynamic limit.

It can be shown analytically that this is the case for $d=1$ and $\alpha>0$.
From Eq. (8) we see that the smallest eigenvalue of $J$ when $N\to\infty$ is
given by

$D_{\rm
min}=\frac{\Gamma}{\tilde{N}}\sum_{\mathbf{r}}\tilde{J}\left(\mathbf{r}\right)\left(-1\right)^{r}=\frac{\Gamma}{\tilde{N}}\left(b+2\sum_{r=1}^{\infty}\left(-1\right)^{r}r^{-\alpha}\right)\,.$
(25)

Here $r=|\mathbf{r}|$. Hence, setting $D_{\rm min}=0$ fixes

$b=-2\sum_{r=1}^{\infty}(-1)^{r}r^{-\alpha}\,.$ (26)

The convergence of this series is proven by means of the alternating series
test, since the absolute value of its terms monotonically decrease to 0. For
$\alpha=0$ and finite $N$, as $\tilde{J}_{i\neq j}=1$, $b$ must be fixed to
$1$ to ensure that the smallest eigenvalue of $\tilde{J}$ is zero. This is
manifestly independent of $N$.

In other dimensions or other models, this test can be done graphically. In
Fig. 6 we show the value of $b$, for a model with power-law decaying
interactions in $d=1$, computed numerically for different values of $N$ when
$b$ is chosen so that the smallest eigenvalue of $J$ is zero. One can see that
in this case the numerical results converge to the analytical prediction. The
same behaviour is observed in other dimensions.

Figure 6: Scaling of $\tilde{N}=\sum_{j}^{N}\tilde{J}_{ij}$ (left) and $b$
(right) for a one dimensional model ($d=1$) with power-law decaying
interactions (1). Here $b$ is fixed such that the smallest eigenvalue of
$\tilde{J}$ is zero. The dotted lines indicate the analytical asymptote.

##  Appendix B Negligible second order corrections to the saddle point method

The second order term of the saddle-point expansion is proportional to $(\det
H_{\phi}[\bar{u}_{k}])^{-1/2}$. Accordingly, it corresponds to a correction to
the free energy per particle of the form

$\frac{1}{N}\ln(\det
H_{\phi}[\bar{u}_{k}])=\frac{1}{N}\sum_{k=0}^{M-1}\ln\nu_{k}\propto\frac{M}{N}\,,$
(27)

where the $\\{\nu_{k}\\}$ are the eigenvalues of $\det H_{\phi}[\bar{u}_{k}]$.
This correction scales as $M/N$ and thus vanishes in the thermodynamic limit.
Notably, if the applicability of the effective theory to map a long-range
interacting model to a generalized Dicke model constitutes the first
appearance of the restriction $\lim_{N\to\infty}M/N=0$, the argument contained
in this Appendix constitutes a second independent one. In fact, this second
occurrence of the restriction also appears in classical systems, where it is
actually the only restriction to $\lim_{N\to\infty}M/N=0$, as in classical
systems an unrestricted standard HST can be used, as outlined in Sec. II.

##  Appendix C Solving the long-range Ising model in transverse field

###  C.1 Solving the associated Dicke model

Particularizing Eq. (11) for the Ising model, we have to compute

$\begin{split}Z_{\rm m}[u_{k}]={\rm Tr}_{\rm
m}\left\\{\exp\left[-\beta\sum_{i}^{N}\Biggl{(}\right.\right.&\frac{\omega_{z}}{2}\sigma_{i}^{z}\\\
&\left.\left.\left.+\sum_{k=0}^{M-1}2\lambda_{ik}u_{k}\sigma_{i}^{x}\right)\right]\right\\}\,.\end{split}$
(28)

Because the spins are now decoupled, the trace over spins factorizes. The
resulting single spin Hamiltonian has eigenvalues

$\epsilon_{i}^{\pm}=\pm\epsilon_{i}=\pm\frac{1}{2}\sqrt{\omega_{z}^{2}+4\left(\sum_{k=0}^{M-1}2\lambda_{ik}u_{k}\right)^{2}}\,.$
(29)

Accordingly,

$\displaystyle Z_{\rm m}[u_{k}]=\prod_{i}2\cosh(\beta\epsilon_{i})\,,$ (30)
$\displaystyle f_{\rm
m}[u_{k}]=\frac{1}{N}\sum_{i}\ln(2\cosh\beta\epsilon_{i})\,,$ (31)
$\displaystyle\phi_{\rm
m}[u_{k}]=-\beta\sum_{k}\omega_{k}u_{k}^{2}+\frac{1}{N}\sum_{i}\ln(2\cosh\beta\epsilon_{i})\,.$
(32)

###  C.2 Existence of a homogeneous maximum of $\phi$

To find the maximum of $\phi[u_{k}]$ we impose a vanishing gradient:
$\nabla\phi=0$, which translates to the following condition for the
maximization variables

$\bar{u}_{k}\omega_{k}=\frac{1}{N}\sum_{i}^{N}\tanh(\beta\bar{\epsilon}_{i})\frac{\lambda_{ik}2\sum_{l=0}^{M-1}\lambda_{il}\bar{u}_{l}}{\bar{\epsilon}_{i}}\,.$
(33)

From here, let us consider a solution that is homogeneous in the lattice, we
will later prove that possible inhomogeneous maxima are only local maxima C.3.
Let us define $\mu_{i}=2\sum_{k}\lambda_{ik}u_{k}$ and consider it as an
alternate optimization variable. Homogeneity implies that
$\bar{\mu}_{i}\equiv\bar{\mu}$, to see how this affects the variables
$\\{\bar{u}_{k}\\}$ it is useful to invert the relation and write $u_{k}$ in
terms of the $\\{\mu_{i}\\}$, yielding

$u_{k}=\frac{1}{2N}\sum_{i}^{N}\lambda_{ik}\mu_{i}\,.$ (34)

Now, homogeneity implies

$\bar{u}_{k}=\frac{\bar{\mu}}{2N}\sum_{i}^{N}\lambda_{ik}\,.$ (35)

Since the $\\{\lambda_{ik}\\}$ are the Fourier modes resulting from the
diagonalization of $J$, we have $\sum_{i}\lambda_{ik}=N\delta_{k0}$.
Accordingly, we find $\bar{u}_{k\neq 0}=0$ and $2\bar{u}_{0}=\bar{\mu}$. So,
if the solution is homogeneous, the only relevant mode is the zero mode
$\bar{u}_{0}\equiv\bar{u}$ and the rest of the maximization variables are
zero, with $\bar{u}$ satisfying the condition

$\bar{u}\omega_{0}=\frac{2\bar{u}}{\bar{\epsilon}}\tanh(\beta\bar{\epsilon})\,.$
(36)

From Eq. 8 we have $\omega_{0}=1/\Gamma$, which when replaced in Eq. (36)
yields Eq. (18).

At this point we can compute the Hessian of $\phi$, $H_{\phi}$. As we have
shown, for a homogeneous solution, the optimization problem becomes single-
valued such that

$\begin{split}H_{\phi}=\frac{\partial^{2}\phi}{\partial
u^{2}}=&-2\beta\omega_{0}+\left[1-\tanh^{2}(\beta\epsilon)\right]\left(\frac{4\beta
u}{\epsilon}\right)^{2}\\\
&+\beta\tanh(\beta\epsilon)\left(\frac{4}{\epsilon}-\frac{16u^{2}}{\epsilon^{3}}\right)\,.\end{split}$
(37)

If we evaluate the Hessian at $\bar{u}=0$, which is always a solution of Eq.
36 we obtain

$H_{\phi}(\bar{u}=0)=-2\beta\omega_{0}+\beta\tanh\left(\beta\frac{\omega_{z}}{2}\right)\frac{8}{\omega_{z}}\,,$
(38)

which is always negative for $\Gamma<\omega_{z}/4$, i.e. for
$\omega_{0}=1/\Gamma>4/\omega_{z}$. For $\omega_{0}>4/\omega_{z}$, the sign
depends on $\beta$, being negative for $\beta<\beta_{c}$ with
$\omega_{0}\omega_{z}=4\tanh\left(\beta_{c}\omega_{z}/2\right)$. So in the
regime where $\bar{u}=0$ is the only solution to Eq. (36), it is a maximum.
For $\omega_{0}<4/\omega_{z}$ and $\beta>\beta_{c}$ a non-trivial solution
given by $\omega_{0}\bar{\epsilon}=2\tanh(\beta\bar{\epsilon})$ appears and
can be shown graphically to be the maximum. Therefore, Eq. (36) always has a
solution that is a maximum of $\phi$.

###  C.3 Proof that the homogeneous solution is the global maximum

We cannot rule out the existence of inhomogenoeus extrema of $\phi$, i.e.
inhomogeneous solutions of Eq. (33). Instead, we show that if there exists an
inhomogeneous solution and it is a maximum, it is a local maximum, with the
global maximum given by the homogeneous solution.

Let us express the self-consistent condition for the extrema of $\phi$ given
in Eq. (33) in terms of $\\{\mu_{i}\\}$ and $\\{\epsilon_{i}\\}$

$\frac{1}{2N}\sum_{i}^{N}\lambda_{ik}\bar{\mu}_{i}\omega_{k}=\frac{1}{N}\sum_{i}^{N}\tanh(\beta\bar{\epsilon}_{i})\frac{\lambda_{ik}\bar{\mu}_{i}}{\epsilon_{i}}\,.$
(39)

Isolating the $\\{\bar{\mu}_{i}\\}$ yields the self-consistent condition

$\bar{\mu}_{i}=2\sum_{j}^{N}\tanh(\beta\bar{\epsilon}_{j})\frac{\bar{\mu}_{j}}{\bar{\epsilon}_{j}}J_{ij}\,,$
(40)

which is simply a reformulation of the maximization problem in terms of new
variables. Accordingly, $\phi$ reads

$\phi[\mu_{i}]=-\frac{\beta}{4N}\sum_{ij}^{N}\mu_{i}J^{+}_{ij}\mu_{j}+\frac{1}{N}\sum_{i}^{N}\ln\left[2\cosh(\beta\epsilon_{i})\right]\,,$
(41)

with $NJ^{+}_{ij}=\sum_{k}\lambda_{ik}\omega_{k}\lambda_{jk}$. Note that
$JJ^{+}J=J$. Substituting Eq. (40) in Eq. (41) yields

$\begin{split}\phi[\bar{\mu}_{i}]&=-\frac{\beta}{2N}\sum_{i}^{N}\frac{\tanh(\beta\bar{\epsilon}_{i})}{\bar{\epsilon}_{i}}\bar{\mu}_{i}^{2}+\frac{1}{N}\sum_{i}^{N}\ln\left[2\cosh(\beta\bar{\epsilon}_{i})\right]\\\
&=\frac{1}{N}\sum_{i}^{N}\xi(\bar{\mu}_{i})\,.\end{split}$ (42)

We can particularize this expression for the homogenous solution,
$\bar{\mu}_{i}\equiv\bar{\mu}$,

$\begin{split}\phi(\bar{\mu})&=-\frac{\beta}{2}\frac{\tanh(\beta\bar{\epsilon})}{\bar{\epsilon}}\bar{\mu}^{2}+\ln\left[2\cosh(\beta\bar{\epsilon})\right]\\\
&=\frac{1}{N}\sum_{i}^{N}\xi(\bar{\mu})\,.\end{split}$ (43)

Note that because $\bar{\mu}$ maximizes $\phi$, it also maximizes $\xi$.
Therefore, an inhomogeneous maximum of $\phi$ given by $\\{\bar{\mu}_{i}\\}$
cannot maximize $\xi$ for all $\bar{\mu}_{i}$ (to the extent that some
$\bar{\mu}_{i}$ must deviate from $\bar{\mu}$ in order for the configuration
to be inhomogeneous) and thus $\phi(\bar{\mu})\geq\phi[\bar{\mu}_{i}]$. The
global maximum of $\phi$ is given by the homogeneous solution.

##  Appendix D Analytical calculation of susceptibilities

From Eqs. (20) and (33) we realize that

$\bar{u}_{k}\omega_{k}=\frac{1}{N}\sum_{i}^{N}\langle\sigma_{i}^{x}\rangle\lambda_{ik}$
(44)

and thus

$\frac{\partial\bar{u}_{k}}{\partial
h_{j}}\omega_{k}=\frac{1}{N}\sum_{i}^{N}\chi_{ij}\lambda_{ik}\,.$ (45)

From Eq. (19) we have

$\chi_{ij}=\lim_{\\{h_{n}\\}\to
0}\left(2\sum_{k=0}^{M-1}\lambda_{ik}\frac{\partial\bar{u}_{k}}{\partial
h_{j}}+\delta_{ij}\right)Y_{i}\,,$ (46)

with

$\begin{split}Y_{i}=&\left(1-\tanh^{2}(\beta\bar{\epsilon}_{i})\right)\beta\left(\frac{\bar{\mu}_{i}+h_{i}}{\bar{\epsilon}_{i}}\right)^{2}\\\
&+\tanh(\beta\bar{\epsilon}_{i})\frac{1}{\bar{\epsilon}_{i}}\left[1-\left(\frac{\bar{\mu}_{i}+h_{i}}{\bar{\epsilon}_{i}}\right)^{2}\right]\,.\end{split}$
(47)

From Eqs. (45) and (46) and after some manipulation, we find

$\chi_{ij}=\left(2\sum_{r}^{N}J_{ir}\chi_{rj}+\delta_{ij}\right)Y\,,$ (48)

where $Y=\lim_{\\{h_{n}\\}\to 0}Y_{i}$. For a translation-invariant model, Eq.
(48) can be solved in Fourier space. We define

$\chi_{ij}=\frac{1}{N}\sum_{k=0}^{N-1}\lambda_{ik}\chi_{k}\lambda_{jk}$ (49)

and find

$\chi_{k}=\frac{Y}{1-2YD_{k}}\,.$ (50)

Here $D_{k}$ are the eigenvalues of $J$, Cf. Eq. (7). The susceptibilities in
real space are thus given by

$\begin{split}\chi_{ij}&=\frac{1}{N}\left(\sum_{k=0}^{M-1}\lambda_{ik}\chi_{k}\lambda_{jk}+Y\sum_{k=M}^{N-1}\lambda_{ik}\lambda_{jk}\right)\\\
&=\frac{1}{N}\left[\sum_{k=0}^{M-1}\lambda_{ik}\chi_{k}\lambda_{jk}+Y\left(N\delta_{ij}-\sum_{k=0}^{M-1}\lambda_{ik}\lambda_{jk}\right)\right]\,.\end{split}$
(51)

leading to Eq. (22).

## References

  * Britton  _et al._ [2012] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. J. Wang, J. K. Freericks, H. Uys, M. J. Biercuk, and J. J. Bollinger, Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of spins, Nature 484, 489 (2012).
  * Knap  _et al._ [2013] M. Knap, A. Kantian, T. Giamarchi, I. Bloch, M. D. Lukin, and E. Demler, Probing real-space and time-resolved correlation functions with many-body ramsey interferometry, Phys. Rev. Lett. 111, 147205 (2013).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L. M. Duan, Z. X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 025001 (2021).
  * Browaeys and Lahaye [2020] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nat. Phys. 16, 132 (2020).
  * Mukamel [2008] D. Mukamel, Statistical mechanics of systems with long range interactions, AIP Conf Proc. 970, 22 (2008).
  * Campa  _et al._ [2009] A. Campa, T. Dauxois, and S. Ruffo, Statistical mechanics and dynamics of solvable models with long-range interactions, Phys. Rep. 480, 57 (2009).
  * Fey and Schmidt [2016] S. Fey and K. P. Schmidt, Critical behavior of quantum magnets with long-range interactions in the thermodynamic limit, Phys. Rev. B 94, 075156 (2016).
  * Defenu  _et al._ [2023a] N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni, Long-range interacting quantum systems, Rev. Mod. Phys. 95, 035002 (2023a).
  * Defenu  _et al._ [2023b] N. Defenu, A. Lerose, and S. Pappalardi, Out-of-equilibrium dynamics of quantum many-body systems with long-range interactions (2023b), arXiv:2307.04802 [cond-mat.quant-gas] .
  * Koffel  _et al._ [2012] T. Koffel, M. Lewenstein, and L. Tagliacozzo, Entanglement entropy for the long-range ising chain in a transverse field, Phys. Rev. Lett. 109, 267203 (2012).
  * Kuwahara and Saito [2020] T. Kuwahara and K. Saito, Area law of noncritical ground states in 1d long-range interacting systems, Nat. Commun. 11, 4478 (2020).
  * Vodola  _et al._ [2014] D. Vodola, L. Lepori, E. Ercolessi, A. V. Gorshkov, and G. Pupillo, Kitaev chains with long-range pairing, Phys. Rev. Lett. 113, 156402 (2014).
  * Ares  _et al._ [2018] F. Ares, J. G. Esteve, F. Falceto, and A. R. de Queiroz, Entanglement entropy in the long-range kitaev chain, Phys. Rev. A 97, 062301 (2018).
  * Vodola  _et al._ [2015] D. Vodola, L. Lepori, E. Ercolessi, and G. Pupillo, Long-range ising and kitaev models: phases, correlations and edge modes, New J. Phys. 18, 015001 (2015).
  * Vanderstraeten  _et al._ [2018] L. Vanderstraeten, M. VanDamme, H. P. Büchler, and F. Verstraete, Quasiparticles in quantum spin chains with long-range interactions, Phys. Rev. Lett. 121, 090603 (2018).
  * Francica and Dell’Anna [2022] G. Francica and L. Dell’Anna, Correlations, long-range entanglement, and dynamics in long-range kitaev chains, Phys. Rev. B 106, 155126 (2022).
  * Schneider  _et al._ [2021] J. T. Schneider, J. Despres, S. J. Thomson, L. Tagliacozzo, and L. Sanchez-Palencia, Spreading of correlations and entanglement in the long-range transverse ising chain, Phys. Rev. Res. 3, L012022 (2021).
  * Jäger  _et al._ [2020] S. B. Jäger, L. Dell'Anna, and G. Morigi, Edge states of the long-range kitaev chain: An analytical study, Phys. Rev. B 102, 035152 (2020).
  * Viyuela  _et al._ [2016] O. Viyuela, D. Vodola, G. Pupillo, and M. A. Martin-Delgado, Topological massive dirac edge modes and long-range superconducting hamiltonians, Phys. Rev.B 94, 125121 (2016).
  * Koziol  _et al._ [2021] J. A. Koziol, A. Langheld, S. C. Kapfer, and K. P. Schmidt, Quantum-critical properties of the long-range transverse-field ising model from quantum monte carlo simulations, Phys. Rev. B 103, 245135 (2021).
  * Lazo  _et al._ [2021] E. G. Lazo, M. Heyl, M. Dalmonte, and A. Angelone, Finite-temperature critical behavior of long-range quantum ising models, SciPost Physics 11, 076 (2021).
  * Mori [2012a] T. Mori, Equilibrium properties of quantum spin systems with nonadditive long-range interactions, Phys. Rev.E 86, 021132 (2012a).
  * Campa  _et al._ [2000] A. Campa, A. Giansanti, and D. Moroni, Canonical solution of a system of long-range interacting rotators on a lattice, Phys. Rev. E 62, 303 (2000).
  * Mori [2011] T. Mori, Instability of the mean-field states and generalization of phase separation in long-range interacting systems, Phys. Rev. E 84, 031128 (2011).
  * Mori [2010] T. Mori, Analysis of the exactness of mean-field theory in long-range interacting systems, Phys. Rev. E 82, 060103(R) (2010).
  * Mori [2012b] T. Mori, Microcanonical analysis of exactness of the mean-field theory in long-range interacting systems, J. Stat. Phys. 147, 1020 (2012b).
  * Campa  _et al._ [2003] A. Campa, A. Giansanti, and D. Moroni, Canonical solution of classical magnetic models with long-range couplings, J. Phys. A Math. Gen. 36, 6897 (2003).
  * Wang and Hioe [1973] Y. K. Wang and F. T. Hioe, Phase transition in the dicke model of superradiance, Phys. Rev. A 7, 831 (1973).
  * Hioe [1973] F. T. Hioe, Phase transitions in some generalized dicke models of superradiance, Phys. Rev. A 8, 1440 (1973).
  * Román-Roche  _et al._ [2021] J. Román-Roche, F. Luis, and D. Zueco, Photon condensation and enhanced magnetism in cavity qed, Phys. Rev. Lett. 127, 167201 (2021).
  * Román-Roche and Zueco [2022] J. Román-Roche and D. Zueco, Effective theory for matter in non-perturbative cavity QED, SciPost Phys. Lect. Notes , 50 (2022).
  * Note [1] Meaning that Kac’s prescription is used to ensure extensivity if the model is strong-long-ranged.
  * Note [2] Fig. 3 is actually an oversimplification, as it only depicts long-range interactions, which are the ones replaced by the auxiliary bosonic modes of the effective theory. The theory is also applicable to models containing a combination of short- and long-range interactions. For a discussion about the applicability of the method see Sec. IV.
  * Simon  _et al._ [2011] J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner, Quantum simulation of antiferromagnetic spin chains in an optical lattice, Nature 472, 307 (2011).
  * Kaicher  _et al._ [2023] M. P. Kaicher, D. Vodola, and S. B. Jäger, Study of the long-range transverse field ising model with fermionic gaussian states (2023), arXiv:2301.02939 [quant-ph] .
  * Gong  _et al._ [2016] Z.-X. Gong, M. F. Maghrebi, A. Hu, M. L. Wall, M. Foss-Feig, and A. V. Gorshkov, Topological phases with long-range interactions, Phys. Rev. B 93, 041102(R) (2016).
  * Lieb [1973] E. H. Lieb, The classical limit of quantum spin systems, Commun. Math. Phys 31, 327 (1973).
  * Brankov  _et al._ [1975] I. G. Brankov, V. A. Zagrebnov, and I. S. Tonchev, Asymptotically exact solution of the generalized dicke model, Theor. Math. Phys. 22, 13 (1975).
  * Gibberd [1974] R. Gibberd, Equivalence of the dicke maser model and the ising model at equilibrium, Aust. J. Phys. 27, 241 (1974).
  * Note [3] The magnetization must be kept $\\{h_{n}\\}$-dependent in order to compute the susceptibility, the magnetization of the $\\{h_{n}\\}$-independent model is defined as $\beta\langle\sigma^{x}_{i}\rangle=\mathop{lim}\displaylimits_{\\{h_{n}\\}\to 0}\partial\mathop{ln}\nolimits Z/\partial h_{i}$.
  * Sadhukhan and Dziarmaga [2021] D. Sadhukhan and J. Dziarmaga, Is there a correlation length in a model with long-range interactions? (2021), arXiv:2107.02508 [cond-mat.str-el] .
  * Kubo  _et al._ [1991] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical Physics II_ (Springer Berlin Heidelberg, Heidelberg, Germany, 1991).
  * Schwabl [2006] F. Schwabl,  _Statistical Mechanics_ (Springer Berlin Heidelberg, Heidelberg, Germany, 2006).
  * Hepp and Lieb [1973a] K. Hepp and E. H. Lieb, On the superradiant phase transition for molecules in a quantized radiation field: the dicke maser model, Ann. Phys 76, 360 (1973a).
  * Hepp and Lieb [1973b] K. Hepp and E. H. Lieb, Equilibrium statistical mechanics of matter interacting with the quantized radiation field, Phys. Rev. A 8, 2517 (1973b).
  * Latorre  _et al._ [2005] J. I. Latorre, R. Orús, E. Rico, and J. Vidal, Entanglement entropy in the lipkin-meshkov-glick model, Phys. Rev. A 71, 064101 (2005).


11institutetext: LaTIM, INSERM, UMR 1101, University Brest, Brest, France

# Squeeze-and-Excitation Normalization for Automated Delineation of Head and
Neck Primary Tumors in Combined PET and CT Images††thanks: Cite this paper as:
Iantsen A., Visvikis D., Hatt M. (2021) Squeeze-and-Excitation Normalization
for Automated Delineation of Head and Neck Primary Tumors in Combined PET and
CT Images. In: Andrearczyk V., Oreiller V., Depeursinge A. (eds) Head and Neck
Tumor Segmentation. HECKTOR 2020. Lecture Notes in Computer Science, vol
12603. Springer, Cham. https://doi.org/10.1007/978-3-030-67194-5_4

Andrei Iantsen  0000-0002-6690-6070 Dimitris Visvikis  0000-0003-0831-3637
Mathieu Hatt  0000-0002-8938-8667

###### Abstract

Development of robust and accurate fully automated methods for medical image
segmentation is crucial in clinical practice and radiomics studies. In this
work, we contributed an automated approach for Head and Neck (H&N) primary
tumor segmentation in combined positron emission tomography / computed
tomography (PET/CT) images in the context of the MICCAI 2020 Head and Neck
Tumor segmentation challenge (HECKTOR). Our model was designed on the U-Net
architecture with residual layers and supplemented with Squeeze-and-Excitation
Normalization. The described method achieved competitive results in cross-
validation (DSC 0.745, precision 0.760, recall 0.789) performed on different
centers, as well as on the test set (DSC 0.759, precision 0.833, recall 0.740)
that allowed us to win first prize in the HECKTOR challenge among 21
participating teams. The full implementation based on PyTorch and the trained
models are available at https://github.com/iantsen/hecktor.

###### Keywords:

Medical imaging Segmentation Head and Neck cancer U-Net SE Normalization.

##  1 Introduction

Combined positron emission tomography / computed tomography (PET/CT) imaging
is broadly used in clinical practice for radiotherapy treatment planning,
initial staging and response assessment. In radiomics analyses, quantitative
evaluation of radiotracer uptake in PET images and tissues density in CT
images, aims at extracting clinically relevant features and building
diagnostic, prognostic and predictive models. The segmentation step of the
radiomics workflow is the most time-consuming bottleneck and variability in
usual semi-automatic segmentation methods can significantly affect the
extracted features, especially in case of manual segmentation, which is
affected by the highest magnitude of inter- and intra-observer variability.
Under these circumstances, a fully automated segmentation is highly desirable
to automate the whole process and facilitate its clinical routine usage.

The MICCAI 2020 Head and Neck Tumor segmentation challenge (HECKTOR) [1] aims
at evaluating automatic algorithms for segmentation of Head and Neck (H&N)
tumors in combined PET and CT images. A dataset of 201 patients from four
medical centers in Québec (CHGJ, CHMR, CHUM and CHUS) with histologically
proven H&N cancer in the oropharynx is provided for a model development. A
test set comprised of 53 patients from a different center in Switzerland
(CHUV) is used for evaluation. All images were re-annotated by an expert for
the purpose of the challenge in order to determine primary gross tumor volumes
(GTV) on which the methods are evaluated using the Dice score (DSC), precision
and recall.

This paper describes our approach based on convolutional neural networks
supplemented with Squeeze-and-Excitation Normalization (SE Normalization or SE
Norm) layers to address the goal of the HECKTOR challenge.

##  2 Materials & Methods

###  2.1 SE Normalization

The key element of our model is SE Normalization layers [6] that we recently
proposed in the context of the Brain Tumor Segmentation Challenge (BraTS 2020)
[3]. Similarly to Instance Normalization [4], for an input
${\mathrm{X}=(x_{1},x_{2},\dots,x_{\mathrm{N}})}$ with $\mathrm{N}$ channels,
SE Norm layer first normalizes all channels of each example in a batch using
the mean and standard deviation:

$x^{\prime}_{i}=\frac{1}{\sigma_{i}}(x_{i}-\mu_{i})$ (1)

where $\mu_{i}=\textrm{E}[x_{i}]$ and
$\sigma_{i}=\sqrt{\textrm{Var}[x_{i}]+\epsilon}$ with $\epsilon$ as a small
constant to prevent division by zero. After, a pair of parameters
$\gamma_{i},\beta_{i}$ are applied to each channel to scale and shift the
normalized values:

$y_{i}=\gamma_{i}x^{\prime}_{i}+\beta_{i}$ (2)

In case of Instance Normalization, both parameters $\gamma_{i},\beta_{i}$,
fitted in the course of training, stay fixed and independent on the input
$\mathrm{X}$ during inference. By contrast, we propose to model the parameters
$\gamma_{i},\beta_{i}$ as functions of the input $\mathrm{X}$ by means of
Squeeze-and-Excitation (SE) blocks [5], i.e

$\displaystyle\gamma$ $\displaystyle=f_{\gamma}(X)$ (3) $\displaystyle\beta$
$\displaystyle=f_{\beta}(X)$ (4)

where ${\gamma=(\gamma_{1},\gamma_{2},\dots,\gamma_{\mathrm{N}})}$ and
${\beta=(\beta_{1},\beta_{2},\dots,\beta_{\mathrm{N}})}$ \- the scale and
shift parameters for all channels, $f_{\gamma}$ \- the original SE block with
the sigmoid, and $f_{\beta}$ is modeled as the SE block with the tanh
activation function to enable the negative shift (see Fig. 1). Both of SE
blocks first apply global average pooling (GAP) to squeeze each channel into a
single descriptor. Then, two fully connected (FC) layers aim at capturing non-
linear cross-channel dependencies. The first FC layer is implemented with the
reduction ratio $r$ to form a bottleneck for controlling model complexity.
Throughout this paper, we apply SE Norm layers with the fixed reduction ration
$r=2$.

Figure 1: Layers with SE Normalization: (1) SE Norm layer, (1) residual layer
with the shortcut connection, and (1) residual layer with the non-linear
projection. Output dimensions are depicted in italics.

###  2.2 Network Architecture

Our model is built upon a seminal U-Net architecture [7, 8] with the use of SE
Norm layers [6]. Convolutional blocks, that form the model decoder, are stacks
of $3\times 3\times 3$ convolutions and ReLU activations followed by SE Norm
layers. Residual blocks in the encoder consist of convolutional blocks with
shortcut connections (see Fig. 1). If the number of input/output channels in a
residual block is different, a non-linear projection is performed by adding
the $1\times 1\times 1$ convolutional block to the shortcut in order to match
the dimensions (see Fig. 1).

In the encoder, downsampling is done by applying max pooling with the kernel
size of $2\times 2\times 2$. To linearly upsample feature maps in the decoder,
$3\times 3\times 3$ transposed convolutions are used. In addition, we
supplement the decoder with three upsampling paths to transfer low-resolution
features further in the model by applying the $1\times 1\times 1$
convolutional block to reduce the number of channels, and utilizing trilinear
interpolation to increase the spatial size of the feature maps (see Fig. 2,
yellow blocks).

The first residual block placed after the input is implemented with the kernel
size of $7\times 7\times 7$ to increase the receptive field of the model
without significant computational overhead. The sigmoid function is applied to
output probabilities for the target class.

Figure 2: The model architecture with SE Norm layers. The input consists of
PET/CT patches of the size of ${144\times 144\times 144}$ voxels. The encoder
consists of residual blocks with identity (solid arrows) and projection
(dashed arrows) shortcuts. The decoder is formed by convolutional blocks.
Additional upsampling paths are added to transfer low-resolution features
further in the decoder. Kernel sizes and numbers of output channels are
depicted in each block.

###  2.3 Data Preprocessing & Sampling

Both PET and CT images were first resampled to a common resolution of
${1\times 1\times 1}$ $mm^{3}$ with trilinear interpolation. Each training
example was a patch of ${144\times 144\times 144}$ voxels randomly extracted
from a whole PET/CT image, whereas validation examples were received from
bounding boxes provided by organizers. Training patches were extracted to
include the tumor class with the probability of 0.9 to facilitate model
training.

CT intensities were clipped in the range of ${\left[-1024,1024\right]}$
Hounsfield Units and then mapped to ${\left[-1,1\right]}$. PET images were
transformed independently with the use of Z-score normalization, performed on
each patch.

###  2.4 Training Procedure

The model was trained for 800 epochs using Adam optimizer on two GPUs NVIDIA
GeForce GTX 1080 Ti (11 GB) with a batch size of 2 (one sample per worker).
The cosine annealing schedule was applied to reduce the learning rate from
$10^{-3}$ to $10^{-6}$ within every 25 epochs.

###  2.5 Loss Function

The unweighted sum of the Soft Dice Loss [9] and the Focal Loss [10] is
utilized to train the model. Based on [9], the Soft Dice Loss for one training
example can be written as

$L_{Dice}(y,\hat{y})=1-\frac{2\sum_{i}^{\mathrm{N}}y_{i}\hat{y}_{i}+1}{\sum_{i}^{\mathrm{N}}y_{i}+\sum_{i}^{\mathrm{N}}\hat{y}_{i}+1}$
(5)

The Focal Loss is defined as

$L_{Focal}(y,\hat{y})=-\frac{1}{\mathrm{N}}\sum_{i}^{\mathrm{N}}y_{i}(1-\hat{y}_{i})^{\gamma}\ln(\hat{y}_{i})$
(6)

In both definitions, ${y_{i}\in\\{0,1\\}}$ \- the label for the i-th voxel,
${\hat{y}_{i}\in\left[0,1\right]}$ \- the predicted probability for the i-th
voxel, and $\mathrm{N}$ \- the total numbers of voxels. Additionally we add +1
to the numerator and denominator in the Soft Dice Loss to avoid the zero
division in cases when the tumor class is not present in training patches. The
parameter $\gamma$ in the Focal Loss is set at 2.

Table 1: The performance results on different cross-validation splits. Average results (the row ’Average’) are provided for each evaluation metric across all centers in the leave-one-center-out cross-validation (first four rows). The mean and standard deviation of each metric are computed across all data samples in the corresponding validation center. The row ’Average (rs)’ indicates the average results on the four random data splits. Center |  DSC |  Precision |  Recall
---|---|---|---
CHUS (${n=72}$) |  $0.744\pm 0.206$ |  $0.763\pm 0.248$ |  $0.788\pm 0.226$
CHUM (${n=56}$) |  $0.739\pm 0.190$ |  $0.748\pm 0.224$ |  $0.819\pm 0.216$
CHGJ (${n=55}$) |  $0.801\pm 0.180$ |  $0.791\pm 0.208$ |  $0.839\pm 0.200$
CHMR (${n=18}$) |  $0.696\pm 0.232$ |  $0.739\pm 0.286$ |  $0.712\pm 0.228$
Average |  0.745 |  0.760 |  0.789
Average (rs) |  0.757 |  0.762 |  0.820

###  2.6 Ensembling

Our results on the test set were produced with the use of an ensemble of eight
models trained and validated on different splits of the training set. Four
models were built using a leave-one-center-out cross-validation, i.e, the data
from three centers was used for training and the data from the fourth center
was held out for validation. Four other models were fitted on random training
/ validation splits of the whole dataset. Predictions on the test set were
produced by averaging predictions of the individual models and applying a
threshold operation with a value equal to ${0.5}$.

Table 2: The test set results of the ensemble of eight models. Center |  DSC |  Precision |  Recall
---|---|---|---
CHUV (${n=53}$) |  0.759 |  0.833 |  0.740

##  3 Results & Discussion

Our validation results in the context of the HECKTOR challenge are summarized
in Table 1. The best outcome in terms of all evaluation metrics was received
for the ’CHGJ’ center with 55 patients. The model demonstrated the poorest
performance for the ’CHMR’ center that is least represented in the whole
dataset. The differences with the two other centers was minor for all
evaluation metrics. The small spread between all centers and the average
results implies that the model predictions were robust and any center-specific
data standardization was not required. This finding is supported by the lack
of significant difference in the average results between the leave-one-center-
out and random split cross-validations.

The ensemble results on the test set consisting of 53 patients from the ’CHUV’
center are presented in Table 2. On the previously unseen data, the ensemble
of eight models achieved the highest results among 21 participating teams with
the Dice score of 75.9%, precision 83.3% and recall 74%.

## References

  * [1] V. Andrearczyk, V. Oreiller, M. Jreige, M. Vallières, J. Castelli, H. Elhalawani, S. Boughdad, J. O. Prior, A. Depeursinge, ”Overview of the HECKTOR challenge at MICCAI 2020: Automatic Head and Neck Tumor Segmentation in PET/CT”, 2021.
  * [2] V. Andrearczyk, V. Oreiller, M. Vallières, J. Castelli, H. Elhalawani, M. Jreige, S. Boughdad, J. O. Prior, A. Depeursinge, ”Automatic Segmentation of Head and Neck Tumors and Nodal Metastases in PET-CT scans”, Medical Imaging with Deep Learning (MIDL), 2020.
  * [3] B. H. Menze et al., ”The Multimodal Brain Tumor Image Segmentation Benchmark (BRATS)” in IEEE Transactions on Medical Imaging, vol. 34, no. 10, pp. 1993-2024, 2015.
  * [4] D. Ulyanov, A. Vedaldi, and V. Lempitsky, ”Instance normalization: The missing ingredient for fast stylization”, arXiv preprint arXiv:1607.08022, 2016.
  * [5] J. Hu, L. Shen, and G. Sun, “Squeeze-and-excitation networks”, CoRR, vol. abs/1709.01507, 2017. [Online]. Available: http://arxiv.org/abs/1709.01507
  * [6] A. Iantsen, V. Jaouen, D. Visvikis, M. Hatt, ”Squeeze-and-Excitation Normalization for Brain Tumor Segmentation”, International MICCAI Brainlesion Workshop, 2020.
  * [7] O. Ronneberger, P. Fischer, T. Brox, ”U-net: Convolutional networks for biomedical image segmentation”, in MICCAI. Springer, 2015, pp. 234–241.
  * [8] Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, O. Ronneberger, ”3D U-Net: Learning dense volumetric segmentation from sparse annotation”, in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2016, pp. 424–432
  * [9] F. Milletari, N. Navab, S.-A. Ahmadi, ”V-net: Fully convolutional neural networks for volumetric medical image segmentation”, in International Conference on 3D Vision. IEEE, 2016, pp. 565–571.
  * [10] T.-Y. Lin, P. Goyal, R. Girshick, K. He, P. Dollár, ”Focal Loss for Dense Object Detection”, arXiv preprint arXiv:1708.02002, 2017.