# Antiferromagnetic magnons and local anisotropy: dynamical mean-field study

A. Niyazi  Institute of Solid State Physics, TU Wien, 1040 Vienna, Austria  D.
Geffroy  Department of Condensed Matter Physics, Faculty of Science, Masaryk
University, Kotlářská 2, 611 37 Brno, Czechia  Institute of Solid State
Physics, TU Wien, 1040 Vienna, Austria  J. Kuneš  Institute of Solid State
Physics, TU Wien, 1040 Vienna, Austria

(August 29, 2024)

###### Abstract

We present a dynamical mean-field study of antiferromagnetic magnons in one-,
two- and three-orbital Hubbard model of square and bcc cubic lattice at
intermediate coupling strength. We investigate the effect of anisotropy
introduced by an external magnetic field or single-ion anisotropy. For the
latter we tune continuously between the easy-axis and easy-plane models. We
also analyze a model with spin-orbit coupling in cubic site-symmetry setting.
The ordered states as well as the magnetic excitations are sensitive to even a
small breaking of $SU(2)$ symmetry of the model and follow the expectations of
spin-wave theory as well as general symmetry considerations.

##  I Introduction

The scalability to multi-orbital systems made the dynamical mean-field theory
(DMFT) Metzner and Vollhardt (1989); Georges and Kotliar (1992); Georges  _et
al._ (1996) a widely used tool for the investigation of electronic
correlations on material specific level Kotliar  _et al._ (2006); Held (2007).
The primary niche of DMFT are the one-particle correlation functions (1PCFs)
such as the generalized band structures. Nevertheless, DMFT allows the
calculation of higher order correlation functions as well. In particular, two-
particle correlation functions (2PCFs) play a crucial role in the description
of continuous phase transitions and, in the form of transport coefficients,
dynamical structure factors or susceptibilities are directly related to
numerous experiments. The solution of the Bethe-Salpeter equation (BSE) makes
the calculation of 2PCFs numerically very demanding. The DMFT calculations of
2PCFs have so far been limited to simple models and high symmetry Jarrell
(1992); Boehnke  _et al._ (2011); Hafermann  _et al._ (2014); van Loon  _et
al._ (2014); Stepanov  _et al._ (2018), or have involved substantial
approximations Park  _et al._ (2011). Numerical and analytical developments in
representation of 2PCFs Shinaoka  _et al._ (2020) and solution of BSE
Wallerberger  _et al._ (2020); Otsuki  _et al._ (2019); Krien (2019); Katanin
(2020) make calculations for realistic models with three or more orbitals and
several atoms in the unit cell feasible.

The collective modes described by 2PCFs play a particularly important role in
ordered phases with spontaneously broken symmetry. Presence or absence of a
gap in 2PCF spectrum related to the type, discrete vs continuous, of broken
symmetry is decisive for finite-temperature stability of the ordered state in
two dimensions Mermin and Wagner (1966). Therefore it is important to
understand not only the analytic properties of a given theory, for DMFT see
Refs. van Loon  _et al._ , 2015; Krien  _et al._ , 2017; Krien, 2018, but also
the properties of the actual numerical implementation. Recently, some of us
demonstrated Geffroy  _et al._ (2019); Niyazi  _et al._ (2020) that the
dynamical susceptibilities obtained using DMFT respect the Goldstone theorem
Goldston (1961); Nambu and Jona-Lasinio (1961); Goldston (1961) in the case of
$U(1)\times U(1)$ symmetry breaking in spinful excitonic condensate Kuneš
(2014); Kuneš and Augustinský (2014); Nasu  _et al._ (2016).

The most common and largely studied long-range order in strongly correlated
materials is the antiferromagnetic (AFM) one Anderson (1952). The AFM
Heisenberg model Manousakis (1991), which describes fluctuating spins in Mott
insulators, was studied in detail with analytical Chakravarty  _et al._ (1989)
as well as with numerically exact methods Sandvik (1997). Investigation of
antiferromagnetism in the fermionic Hubbard model, which allows for the
description of doped Mott insulators or AFM metals, relies on weak-coupling
methods such as the random phase approximation Rowe  _et al._ (2012); Del Re
and Toschi (2021), two-particle self-consistent approximation Vilk  _et al._
(1994) or fluctuation-exchange approximation Bickers  _et al._ (1989),
numerical simulations on finite systems Lüscher and Läuchli (2009); Varney
_et al._ (2009) or local approximations to interaction vertices such as DMFT
and its cluster Kent  _et al._ (2005); Fuchs  _et al._ (2011) or diagrammatic
extensions Hirschmeier  _et al._ (2015); Rohringer  _et al._ (2018).

In this Article, we use DMFT and 2PCFs to study AFM ordering in the half-
filled Hubbard model with one, two and three orbitals. While the magnetic
phase diagram of the single-orbital Hubbard model in three dimensions (3D) has
been investigated with a number of methods Kent  _et al._ (2005); Fuchs  _et
al._ (2011); Rohringer  _et al._ (2011) including DMFT Ulmke  _et al._ (1995);
Kent  _et al._ (2005); Rohringer  _et al._ (2011), we calculate the magnon
spectra and we introduce various symmetry breaking terms such as external
field or single-ion anisotropy, in order to analyze their effect on magnon
dispersions. The calculations are performed for 2D and 3D lattices. Our goal
is to demonstrate the utility of the present approach for the study of the
magnetically ordered phases of multi-orbital models. In the 2D case, where
spontaneous symmetry breaking is forbidden in the case of continuous symmetry
Mermin and Wagner (1966) but allowed in case of discrete symmetry, we
investigate a model, which can be continuously tuned between these cases.

##  II Computational Method

The studied Hamiltonian consists of the inter-site hopping, diagonal in
orbital and spin indices, and the on-site part $H_{i}$ consisting of the
electron-electron interaction, external Zeeman field and spin-orbit coupling
(the specific form will be given later for each studied case)

$H=t\sum_{<ij>,\sigma}\sum_{l=1}^{N}c_{il\sigma}^{\dagger}c_{jl\sigma}^{\phantom{\dagger}}+\sum_{i}H_{i}.$
(1)

Here, $c^{{\dagger}}_{il\sigma}$ and $c_{il\sigma}^{\phantom{\dagger}}$ are
the fermionic creation and annihilation operators for electrons with spin
$\sigma$ in orbital $l$ at site $i$ of a square or bcc cubic lattice. The
number of orbitals $N$ ranges from 1 to 3 in the studied models. For later use
we define the occupation number operator $n_{il\sigma}\equiv
c^{{\dagger}}_{il\sigma}c_{il\sigma}^{\phantom{\dagger}}$ and the local spin
operators
$S^{\alpha}_{i}\equiv\sum\limits_{l\nu\nu^{\prime}}\sigma^{\alpha}_{\nu\nu^{\prime}}c^{{\dagger}}_{il\nu}c_{il\nu^{\prime}}^{\phantom{\dagger}}$,
where $\sigma^{\alpha}$ are the Pauli matrices. For the sake of consistency
with previous work Kuneš (2011) we choose $t=1/8$ and use it for all studied
cases.

The calculations follow the standard DMFT procedure. The lattice model is
mapped onto an auxiliary Anderson impurity model with self-consistently
determined parameters Georges and Kotliar (1992); Jarrell (1992), for which
the 1PCFs are evaluated using the ALPS implementation Bauer  _et al._ (2011);
Shinaoka  _et al._ (2017); Gaenko  _et al._ (2017) of the strong-coupling
continuous-time quantum Monte-Carlo (CT-QMC) algorithm Werner  _et al._
(2006). The model hosts two competing phases: the normal paramagnetic one and
the AFM phase with a staggered spin configuration characterized by the Néel
vector $\mathbf{N}=\tfrac{1}{2}\langle\mathbf{S}_{A}-\mathbf{S}_{B}\rangle$.
We focus on the dynamical spin susceptibility
$\chi^{\alpha\alpha}(\mathbf{q},\omega)$, which is obtained by analytic
continuation from its Matsubara representation

$\frac{1}{N}\expectationvalue{\tilde{S}^{\alpha}_{-\mathbf{q}}\tilde{S}^{\alpha}_{\mathbf{q}}}_{\omega_{n}}=\frac{1}{N}\int_{0}^{\beta}\\!\\!\\!\\!\\!\mathrm{d}\tau
e^{i\omega_{n}\tau}\expectationvalue{\tilde{S}^{\alpha}_{-\mathbf{q}}(\tau)\tilde{S}^{\alpha}_{\mathbf{q}}(0)}.$
(2)

Here $\expectationvalue{X}=\frac{1}{Z}\Tr Xe^{-\beta H}$ denotes the thermal
average, $N$ is the number of lattice sites and
$\tilde{X}=X-\expectationvalue{X}$. We also evaluate the 1P observables such
as the static magnetization $\expectationvalue{\mathbf{S}_{\mathbf{R}}}$ and
the 1P spectral functions. The analytic continuation employs the maximum
entropy method Gubernatis  _et al._ (1991), for details of the present
implementation see Refs. Geffroy  _et al._ , 2019; Kaufmann and Held, 2021.

The calculations are performed in a two-site unit cell where the sites are
labeled by the sublattice index ${s=A,B}$. The reciprocal space operators in
Eq. 2 are then given by
${\tilde{S}^{\alpha}_{\mathbf{q}}=\sum_{\mathbf{R}}e^{-i\mathbf{q}\cdot\mathbf{R}}\left(\tilde{S}^{\alpha}_{\mathbf{R},A}+e^{-i\mathbf{q}\cdot\mathbf{s}}\tilde{S}^{\alpha}_{\mathbf{R},B}\right)}$,
where the sublattice vector $\mathbf{s}$ assumes the value
$(\tfrac{1}{2},\tfrac{1}{2})$ or $(\tfrac{1}{2},\tfrac{1}{2},\tfrac{1}{2})$
for the 2D and 3D models, respectively. The correlation function in Eq. 2 is
then obtained as a linear combination of the sublattice contributions

$\begin{split}\expectationvalue{\tilde{S}_{-\mathbf{q}}\tilde{S}_{\mathbf{q}}}_{\omega}=\expectationvalue{\tilde{S}_{-\mathbf{q},A}\tilde{S}_{\mathbf{q},A}}_{\omega}+\expectationvalue{\tilde{S}_{\mathbf{q},B}\tilde{S}_{-\mathbf{q},B}}_{\omega}\\\
+e^{i\mathbf{q}\cdot\mathbf{s}}\expectationvalue{\tilde{S}_{-\mathbf{q},B}\tilde{S}_{\mathbf{q},A}}_{\omega}+e^{-i\mathbf{q}\cdot\mathbf{s}}\expectationvalue{\tilde{S}_{-\mathbf{q},A}\tilde{S}_{\mathbf{q},B}}_{\omega}.\end{split}$

Each term is a contraction of the generalized susceptibility
$\chi_{ijs,kls^{\prime}}(\mathbf{q},\omega)$ over the spin-orbital indices

$\begin{split}\frac{1}{N}\expectationvalue{\tilde{S}^{\alpha}_{-\mathbf{q},s}\tilde{S}^{\alpha}_{\mathbf{q},s^{\prime}}}_{\omega}=M^{\alpha}_{\underline{ij}}M^{\alpha}_{\underline{kl}}\chi_{\underline{ij}s,\underline{kl}s^{\prime}}(\mathbf{q},\omega).\end{split}$

with matrix elements $M^{\alpha}_{ij}$ following from the definition of the
spin operators above. To simplify the notation we use underline to indicate
summation over the corresponding indices 111In case of summation over the
reciprocal lattice vectors this amounts to $\frac{1}{N}\sum_{\mathbf{k}}$, in
case of summation over Matsubara frequencies to $T\sum_{\nu}$ where $T$ is the
temperature..

The calculation of the generalized susceptibility within DMFT Jarrell (1992);
Georges  _et al._ (1996); Kuneš (2011); Boehnke  _et al._ (2011) requires a
more general object – the 2PCF
$X_{ijs\underline{\nu},kls^{\prime}\underline{\nu^{\prime}}}(\mathbf{q},\omega)$,
where $i,j,k,l$ are the spin-orbital indices, $s,s^{\prime}$ are the
sublattice indices, and $\nu,\nu^{\prime}$ are fermionic indices, which
represent the imaginary time evolution, e.g. Matsubara frequencies or
imaginary time. Since the fermionic indices appear only as dummy variables in
the BSE, the equation is invariant under their unitary transformation. In the
present calculations we use the Legendre basis Boehnke  _et al._ (2011) for
the fermionic indices. The susceptibility
$\chi_{\ldots,\ldots}(\mathbf{q},\omega)$ is obtained by the contraction of
$X_{\dots\nu,\ldots\nu^{\prime}}(\mathbf{q},\omega)$ with the basis dependent
structure factor $F_{\nu}$ Boehnke  _et al._ (2011) 222$F_{\nu}=1$ in the
Matsubara frequency basis.

$\chi_{ijs,kls^{\prime}}(\mathbf{q},\omega)=F_{\underline{\nu}}F_{\underline{\nu}^{\prime}}X_{ijs\underline{\nu},kls^{\prime}\underline{\nu}^{\prime}}(\mathbf{q},\omega).$
(3)

The 2PCF $X_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)$ is the
solution of the lattice BSE (4) using the local 2P-irreducible vertices
$\Gamma^{s}_{ij\nu,kl\nu^{\prime}}(\omega)$ and the lattice bubbles
$X^{0}_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)$. The vertices are
related to the impurity 2PCF
$\operatorname{x}_{ij\nu,kl\nu^{\prime}}(s;\omega)$ via the impurity BSE (5)
for each sublattice $s$.

$\displaystyle
X_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)=X^{0}_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)+X^{0}_{ijs\nu,\underline{mns_{1}\nu_{1}}}(\mathbf{q},\omega)\Gamma^{\underline{s_{1}}}_{\underline{mn\nu_{1}},\underline{pq\nu_{2}}}(\omega)X_{\underline{pqs_{1}\nu_{2}},kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)$
(4)
$\displaystyle\operatorname{x}_{ij\nu,kl\nu^{\prime}}(s;\omega)=\operatorname{x}^{0}_{ij\nu,kl\nu^{\prime}}(s;\omega)+\operatorname{x}^{0}_{ij\nu,\underline{mn\nu_{1}}}(s;\omega)\Gamma^{s}_{\underline{mn\nu_{1}},\underline{pq\nu_{2}}}(\omega)\operatorname{x}_{\underline{pq\nu_{2}},kl\nu^{\prime}}(s;\omega).$
(5)

The lattice and local bubbles

$\displaystyle
X^{0}_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)=-\delta_{\nu\nu^{\prime}}G_{is,ks^{\prime}}(\underline{\mathbf{k}}\\!+\\!\mathbf{q},\nu\\!+\\!\omega)G_{ls^{\prime},js}(\underline{\mathbf{k}},\nu)$
$\displaystyle\operatorname{x}^{0}_{ij\nu,kl\nu^{\prime}}(s;\omega)=-\delta_{\nu\nu^{\prime}}G_{is,ks}(\underline{\mathbf{k}},\nu\\!+\\!\omega)G_{ls,js}(\underline{\mathbf{k}^{\prime}},\nu).$

are obtained from the 1P propagator

$G_{is,js^{\prime}}(\mathbf{k},\nu)=\left[i\nu-
h_{\mathbf{k}}-\Sigma(\nu)\right]^{-1}_{is,js^{\prime}}.$

Note that in the 2P quantities, such as
$X_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)$, the spin-orbital
indices $i$ and $j$ ($k$ and $l$) point to the same lattice site and thus
share the sublattice index $s$ ($s^{\prime}$). This is because we are
interested in correlators of local operators, i.e., products of the type
$c_{\mathbf{R}s\ldots}^{\dagger}c_{\mathbf{R}s\ldots}^{\phantom{\dagger}}$,
and because the DMFT vertex $\Gamma^{s}_{jk\nu,kl\nu^{\prime}}(\omega)$ is
local. As a result the corresponding matrices scale with the square of the
number of spin-orbitals per site $2N$, but only linearly with the number of
sites per unit cell.

We used the 1P propagators at 300 Matsubara frequencies and a uniform
$55\times 55$ ($25\times 25\times 25$) $k$-mesh in the 2D (3D) case to compute
the lattice and local bubbles
$X_{ijs\nu,kls^{\prime}\nu^{\prime}}(\mathbf{q},\omega)$ and
$\operatorname{x}^{0}_{ij\nu,kl\nu^{\prime}}(s;\omega)$. These are then
transformed into the Legendre representation Boehnke  _et al._ (2011). The 2P
correlation function $\operatorname{x}_{ij\nu,kl\nu^{\prime}}(s;\omega)$ is
sampled using the CT-QMC directly in the Legendre basis. The local
2P-irreducible vertex $\Gamma^{s}_{ij\nu,kl\nu^{\prime}}(\omega)$ is obtained
from the impurity BSE (5). Next we solve the lattice BSE (4) independently for
each bosonic Matsubara frequency and $\mathbf{q}$-point. We use from 22 (for
the zeroth bosonic frequency) to 45 Legendre coefficients (for the fourteenth
bosonic frequency). A sizable reduction of the computational and storage cost
can be achieved with the procedure of Refs. Otsuki  _et al._ , 2017; Shinaoka
_et al._ , 2020.

We found that 15 bosonic Matsubara frequencies allow for a stable analytic
continuation in the low-energy region, which for the chosen interaction
strengths dominates over the high-energy particle-hole continuum, which we do
not attempt to extract. In particular, it allows to detect the opening of spin
gaps and the suppression of otherwise divergent spectral weights caused by
minor lowering of the Hamiltonian symmetry. Given the insufficient data
statistics we have used a conservative estimate of the uncorrelated standard
deviations of 0.02 at each Matsubara frequency point Kaufmann and Held (2021).
Similar calculations for weaker interaction strength, for which the particle-
hole continuum co-exists with magnon modes, would require more Matsubara
frequencies.

The spectral functions
$B^{\alpha\alpha}_{\mathbf{k}}(\omega)=-\tfrac{1}{\pi}\imaginary\chi^{\alpha\alpha}(\mathbf{q},\omega)$
are our final product. Given the divergence of $B_{\mathbf{k}}(\omega)$ at the
ordering wave vector, we introduce an effective cutoff by plotting
$\tilde{B}_{\mathbf{k}}(\omega)=\tfrac{B_{\mathbf{k}}(\omega)}{C+B_{\mathbf{k}}(\omega)}$
instead of $B_{\mathbf{k}}(\omega)$ itself. To represent the amplitude we also
plot the spectral weight
$\Omega_{\mathbf{k}}=\int_{0}^{0.5}B_{\mathbf{k}}(\omega)d\omega$.

Figure 1: The AFM ordering in the single-band Hubbard model. Left: Imaginary
part of the dynamical susceptibility $\chi(\mathbf{q},\omega)$
($\tilde{B}_{\mathbf{k}}(\omega)$ with $C=5.5$ is plotted) along the
$\Gamma(0,0)$–$X(\pi,\pi)$–$M(2\pi,0)$ path in the extended Brillouin zone for
the 2-site unit cell plotted as color plot for various temperatures. The top
row corresponds to the normal state at $T\approx 1.12T_{c}$. The lower rows
are obtained in the AFM state with $\mathbf{N}$ along the $x$-axis. The white
line (right axis) shows the integrated spectral weight $\Omega_{\mathbf{k}}$.
Right: The corresponding 1P spectral functions. Figure 2: (a) The phase
boundary between AFM and normal (N) phase adapted from Ref. Kuneš, 2011. (b)
The schematic phase diagram in the $T-h$ plane ($U=1$). Symbols mark the
points of actual calculations.

Figure 3:  The same as in Fig. 1 for various amplitudes of the external field
$\bm{h}$ (along the $z$-axis) at $T=1/30$. The top row corresponds to the AFM
state without field, the bottom row to a fully polarized normal state. The
arrows show the size and tilt of the sublattice spin polarization. The Néel
vector $\mathbf{N}$ points along the $x$-axis.

##  III Results and Discussion

###  A S=1/2 and magnetic field

First, we discuss the antiferromagnetism of the single orbital 2D model. The
local term in Eq. 1 adopts the form

$H_{i}=Un_{i\uparrow}n_{i\downarrow}+h\left(n_{i\uparrow}-n_{i\downarrow}\right).$
(6)

At $U$=1 the model is close to the maximal transition temperature between the
weak-coupling RPA and strong-coupling Heisenberg regimes Kuneš (2011);
Rohringer and Toschi (2016); Schäfer  _et al._ (2021). In Fig. 1 we show the
evolution of the electron spectral density and the dynamical spin
susceptibility across the AFM transition along with the local 1P spectral
densities. The location of the studied temperatures in the phase diagram is
shown in Fig. 2a. The direction of the staggered magnetization,
$\expectationvalue{\mathbf{S}_{\mathbf{R}}}=(-1)^{|\mathbf{R}|}\mathbf{N}$ is
chosen along the $x$-axis. The elements $\chi^{yy}$ and $\chi^{zz}$ reflect
the two linear Goldstone modes Watanabe and Murayama (2012) arising from
breaking of $SU(2)$ symmetry in an antiferromagnet. Note, that their spectral
weight diverges at the ordering wave vector ($M$-point). While the Hamiltonian
(1) for $h=0$ is isotropic in spin space, the numerical treatment of
$\chi^{yy}$ is kept independent of that of $\chi^{zz}$. Nevertheless, the
results reflect the symmetry quite accurately.

Figure 4: (top) The spectral weight of the transverse ($xx$ or $yy$) magnons
of Fig. 3 compared to the induced moment $\langle S_{z}\rangle$. (bottom) The
static uniform transverse susceptibility $\chi^{xx}(0,0)=\chi^{yy}(0,0)$
compared to its exact value $\tfrac{\langle S_{z}\rangle}{h}$. The dotted line
in the top panel is piecewise: a quadratic fit to the data and a constant
(reflecting the spin saturation). In the lower panel the same curve is divided
by $h$. Note that at $h=0$ the equivalence of $\chi^{xx}$ and $\chi^{yy}$ is
broken. The empty circle marks the value of $\chi^{xx}(0,0)$, while the full
circle that of $\chi^{yy}(0,0)$.

Application of a uniform magnetic field is known to polarize an isotropic
antiferromagnet perpendicularly to the field, $\mathbf{N}\perp\bm{h}$, while
the spin density acquires a uniform component parallel to the field:
$\expectationvalue{\mathbf{S}_{\mathbf{R}}}=\mathbf{M}+(-1)^{|\mathbf{R}|}\mathbf{N}$,
with $\mathbf{M}\parallel\bm{h}$. The symmetry of the Hamiltonian is reduced
from $SU(2)$ to $U(1)$ 333We mention only the spin symmetries, which are
broken or varied in our models and ignore the ones, which are not changed such
the global $U(1)$ due to charge conservation.. Its breaking leads to a single
linear Goldstone mode ($\chi^{yy}$) with polarization perpendicular to both
$\mathbf{N}$ and $\mathbf{M}$, which corresponds to spin rotation in the
$xy$-plane. Indeed, the ($\chi^{zz}$) mode polarized along the field direction
is gapped, as shown in Fig. 3. At the same time its spectral weight at the
$M$-point becomes finite. $\chi^{yy}$ acquires a finite spectral weight at the
$\Gamma$-point at energy equal to $2h$ Krien (2018). With increasing $h$ and
progressing tilt of the local moments, $\chi^{zz}$ loses its spectral weight
while the weight of $\chi^{xx}$ grows from the $\Gamma$-point outwards. These
observations agree with zero temperature exact diagonalization of Ref. Lüscher
and Läuchli, 2009, although multi-magnon satellites are as expected absent in
the DMFT spectra. Eventually, the AFM order is lost, as shown in Fig. 2b, the
system becomes fully spin-polarized and recovers the $U(1)$ symmetry of its
Hamiltonian. The perpendicular susceptibilities ($\chi^{xx}$ and $\chi^{yy}$)
describe gapped magnons (perhaps better called magnetic excitons - a spin flip
propagating in a spin-polarized background). The field-driven transition can
be viewed as the Bose-Einstein condensation of these magnons, which takes
place when the magnon gap is closed. The direction of $\mathbf{N}$ within the
$xy$-plane plays the role of the condensate phase.

The uniform susceptibilities (at the $\Gamma$-point) reflect the conservation
of the total spin moment along the field direction ($S^{z}$), which leads to
(i) ${\chi^{zz}(0,\omega)=\beta\expectationvalue{(S^{z})^{2}}\delta(\omega)}$
having no dynamics and therefore vanishing imaginary part, (ii)
${\chi^{xx}(0,\omega)=\chi^{yy}(0,\omega)}$. For a simple proof see the
Appendix A. Our empirical observations showed that the numerical noise in
$\chi^{\alpha\alpha}(\mathbf{q},\omega_{n})$ is essentially independent of
$\mathbf{q}$. This is consistent with the fact that such noise mostly
originates from the QMC calculation of the local 2PCF, which is used in the
determination of the local irreducible vertex. As a consequence, one may
subtract
$\chi^{\alpha\alpha}(\mathbf{q},\omega_{n})-\chi^{\alpha\alpha}(0,\omega_{n})$
for $\omega_{n}>0$ in order to reduce the noise in cases where (i) holds.

The uniform susceptibility in external field offers further simple consistency
tests. First, its static part is exactly $\chi^{xx}(0,0)=\tfrac{\langle
S_{z}\rangle}{h}$, reflecting the fact that application of an infinitesimal
transverse field simply rotates the net moment in the new field direction.
Second, the spectral weight is equal to the net moment
$\Omega^{xx}_{0}=\langle S_{z}\rangle$ (see Eq. D.4 of Ref. Krien, 2018). Fig.
4 shows that our numerical results respcect these properties with great
accuracy across all field values.

###  B S=1 and single-ion anisotropy

Next, we investigate the effect of single-ion anisotropy on the magnon
dispersion. To this end we study a two-orbital model (at half filling $n=2$)

$\begin{split}H_{i}=&U\sum_{l=1,2}n_{il\uparrow}n_{il\downarrow}+U^{\prime}\sum_{\sigma,\sigma^{\prime}}n_{i1\sigma}n_{i2\sigma^{\prime}}\\\
-&J\sum_{\sigma}\left(n_{i1\sigma}n_{i2\sigma}+\gamma
c_{i1\sigma}^{\dagger}c_{i1\bar{\sigma}}^{\phantom{\dagger}}c_{i2\bar{\sigma}}^{\dagger}c_{i2\sigma}^{\phantom{\dagger}}\right)\end{split}$
(7)

with $U=1$, $J=0.25$ and $U^{\prime}=0.5$ 444For the sake of simplicity we
have omitted the pair-hopping term
$c_{i1\uparrow}^{\dagger}c_{i1\downarrow}^{\dagger}c_{i2\downarrow}c_{i2\uparrow}+H.c.$
This term does not affect the broken spin symmetry and is expected to play a
minor role when the local high-spin state dominates as in the AFM phases
studied here. The single ion anisotropy is introduced by unequal weights of
the Ising and spin-flip terms in the interaction Hamiltonian. This way the
$SU(2)$ symmetry ($\gamma=1$) is reduced to a $Z_{2}\times U(1)$ for
$\gamma\neq 1$. The residual symmetry of the AFM state depends on $\gamma$.
For $\gamma>1$ the atomic ground state corresponds to
$\ket{S,S_{z}}=\ket{1,0}$, i.e., a state with no spin dipole moment. However,
for moderate deviations $\gamma\gtrsim 1$ the inter-site exchange, which
favors (dipole) magnetic order, enforces AFM order within the $xy$-plane. The
in-plane order breaks the $U(1)$ symmetry and thus one linear Goldstone mode
is expected. For $\gamma<1$ the atomic ground state corresponds to
$\ket{S,S_{z}}=\ket{1,\pm 1}$, i.e., an Ising ground state. The inter-site
exchange leads to the formation of AFM order with moments along the $z$-axis.
The residual symmetry of the ordered state is $U(1)$ and only the discrete
$Z_{2}$ symmetry is broken at the transition. Therefore no gapless Goldstone
mode is expected. The numerical results are presented in Fig. 5. In the
$SU(2)$ symmetric case we have numerically tested the stability of the AFM
state with arbitrary orientation of the Néel vector $\mathbf{N}$. We observe
two linear Goldstone modes in the response to a field perpendicular to the
Néel vector, i.e., $\chi^{yy}$ and $\chi^{zz}$ for $\mathbf{N}\parallel x$.
For $\gamma\neq 1$ the system self-consistently picks the expected Néel
vector. For $\gamma>1$ with the in-plane order, we choose $\mathbf{N}\parallel
x$. As in the $SU(2)$ symmetric case we find minor longitudinal response and
progressive gapping of the out-of-plane $\chi^{zz}$ component with increasing
$\gamma$. For $\gamma<1$ the system picks $\mathbf{N}\parallel z$. In this
case, the residual $U(1)$ symmetry is reflected in the equivalence of
$\chi^{xx}=\chi^{yy}$, both of which are progressively gapped when lowering
$\gamma<1$.

The presence of a spin gap is known to stabilize the long-range order in 2D.
On the other hand, linear gapless magnon mode is detrimental to the long-range
order at any finite temperature Mermin and Wagner (1966) if the thermal
population of this mode is properly taken into account. This is not the case
of the DMFT treatment, in which the low-energy long-range spin fluctuations do
not feed back to the calculation of the 1P and 2P vertices. Nevertheless, the
present results show that DMFT accurately captures the behavior of the spin
gap. Therefore it provides a useful reference and a starting point for more
sophisticated approaches such as D$\Gamma$A Toschi  _et al._ (2007); Kaufmann
_et al._ (2021) or dual fermions Rubtsov  _et al._ (2008); Rohringer  _et al._
(2018). The model with tunable $\gamma$ not only hosts both states with gapped
and gappless magnons, but allows a continuous tuning between them.
Investigation of $T_{c}(\gamma)$ dependence may thus provide a useful test of
the above methods as well as a quantitative measure of their accuracy.

Figure 5:  Imaginary part of the dynamical susceptibility
$\chi(\mathbf{q},\omega)$ ($\tilde{B}_{\mathbf{k}}(\omega)$ with $C=8.5$ is
plotted) at $T=1/30$ in the AFM state of two-orbital Hubbard model for various
single-ion anisotropy. The rows 1-2 correspond to easy ($xy$) plane, row 3 to
$SU(2)$ symmetry, and rows 4-5 to easy ($z$) axis anisotropy. The white line
(right axis) shows the the integrated spectral weight $\Omega_{\mathbf{k}}$.

###  C S=3/2 and spin-orbit coupling

Finally, we study the more realistic 3D case of a three orbital model on a bcc
cubic lattice and introduce spin-orbit coupling as the source of spin
anisotropy

$\begin{split}H_{i}=&\sum_{l,l^{\prime}}\sum_{\sigma,\sigma^{\prime}}h^{\text{soc}}_{l\sigma,l^{\prime}\sigma^{\prime}}c_{il\sigma}^{\dagger}c_{il^{\prime}\sigma^{\prime}}^{\phantom{\dagger}}\\\
+&U\sum_{l}n_{il\uparrow}n_{il\downarrow}+U^{\prime}\sum_{l>l^{\prime}}\sum_{\sigma,\sigma^{\prime}}n_{il\sigma}n_{il^{\prime}\sigma^{\prime}}\\\
-&J\sum_{l>l^{\prime}}\sum_{\sigma}\left(n_{il\sigma}n_{il^{\prime}\sigma}+c_{il\sigma}^{\dagger}c_{il\bar{\sigma}}^{\phantom{\dagger}}c_{il^{\prime}\bar{\sigma}}^{\dagger}c_{il^{\prime}\sigma}^{\phantom{\dagger}}\right),\end{split}$
(8)

where $U=2$, $J=0.5$ and $U^{\prime}=1$ The form of
$h^{\text{soc}}_{l\sigma,l^{\prime}\sigma^{\prime}}$ is given in Appendix B.
We point out that the present CT-QMC calculation may become impossible at low
temperatures due to the sign problem associated with the spin-orbit coupling
Kim  _et al._ (2020). Nevertheless, in the half-filled case the effect of
spin-orbit coupling is rather moderate Kim  _et al._ (2020); Ahn  _et al._
(2017) and we are able to reach the AFM phase without problems. The Néel
vector picks an orientation along a cubic axis ($\mathbf{N}\parallel z$). In
Fig. 6 we show the calculated transverse susceptibilities. Without spin-orbit
coupling we observe a linear Goldstone mode. The spin-orbit coupling leads to
the opening of a finite spin gap. In the 3D case, the AFM order is physically
realistic. The DMFT transition temperature $T_{c}$ provides a realistic
estimate, which misses corrections due to long-range spin fluctuations, while
in 2D the corrections dominate.

Figure 6: Imaginary part of the transverse dynamical susceptibility
$\chi(\mathbf{q},\omega)$ ($\tilde{B}_{\mathbf{k}}(\omega)$ with $C=20$ is
plotted) at $T=1/40$ in the AFM state of 3-orbital Hubbard model
($\mathbf{N}\parallel z$-axis) with (right) and without (left) the spin-orbit
coupling. The data are plotted along the
$\Gamma(0,0,0)$–$H(2\pi,0,0)$–$N(\pi,\pi,0)$–$P(\pi,\pi,\pi)$ path in the
extended Brillouin zone for the 2-site unit cell. The white line (right axis)
shows the the integrated spectral weight $\Omega_{\mathbf{k}}$.

##  IV Conclusions

We have presented DMFT calculations of the AFM phase of the half-filled
Hubbard model with one, two and three orbitals in the intermediate coupling
regime in two and three dimensions. We find that the expected behavior of
magnons in response to external magnetic field or single-ion anisotropy is
well captured. While the 3D description of AFM ordering is physically
relevant, the ordering behavior in 2D is not correct since DMFT is not
sensitive to dimensionality and violates the Mermin-Wagner theorem, which
prohibits spontaneous breaking of continuous symmetry at finite temperature.
Nevertheless, the fact that DMFT magnons properly describe the opening of the
spin gap suggests that DMFT is a good starting point of theories, which
properly include the long-range spin fluctuations. The two-orbital model with
variable single-ion anisotropy provides an ideal test case for such theories
as it allows one to continuously tune between the breaking of a discrete and
of a continuous symmetry. On the computational level we have shown that
3-orbital calculations, which cover for example models of ruthenates or
iridates, are numerically feasible.

###### Acknowledgements.

The authors thank A. Kauch, F. Krien, J. Kaufmann, A. Toschi and M.
Wallerberger for comments and critical reading of the manuscript. This work
was supported by the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No.
646807-EXMAG). D.G was supported by the Czech Science Foundation (GAČR) under
Project No. GA19-16937S. The authors acknowledge support by the Czech Ministry
of Education, Youth and Sports from the Large Infrastructures for Research,
Experimental Development and Innovations project “IT4Innovations National
Supercomputing Center – LM2015070”. Part of the calculations were performed at
the Vienna Scientific Cluster.

## Appendix A

The susceptibilities at the $\Gamma$-point correspond to correlators of total
spin momenta $S^{\alpha}=\sum_{\mathbf{R}}S^{\alpha}_{\mathbf{R}}$. The
Hamiltonian (6) commutes with $S^{z}$, $[S^{z},H]=0$. First, we show that the
$zz$ correlator does not depend on the imaginary time $\tau$

$\expectationvalue{S^{z}(\tau)S^{z}}\equiv\expectationvalue{e^{\tau
H}S^{z}e^{-\tau H}S^{z}}=\expectationvalue{(S^{z})^{2}}.$

This implies that only $\chi^{zz}(0,\nu_{n}=0)$ is finite and equal to
$\beta\expectationvalue{(S^{z})^{2}}$.

To prove the equality of $\chi^{xx}$ and $\chi^{yy}$ we write the
corresponding spin-spin correlation functions with the help of ladder
operators ${S^{x}=S^{-}+S^{+}}$ and ${S^{y}=i(S^{-}-S^{+})}$:

$\begin{split}\expectationvalue{S^{x}(\tau)S^{x}}=&\expectationvalue{e^{\tau
H}S^{+}e^{-\tau H}S^{-}}+\expectationvalue{e^{\tau H}S^{-}e^{-\tau H}S^{+}}\\\
&+\expectationvalue{e^{\tau H}S^{-}e^{-\tau H}S^{-}}+\expectationvalue{e^{\tau
H}S^{+}e^{-\tau H}S^{+}}\end{split}$
$\begin{split}\expectationvalue{S^{y}(\tau)S^{y}}=&\expectationvalue{e^{\tau
H}S^{+}e^{-\tau H}S^{-}}+\expectationvalue{e^{\tau H}S^{-}e^{-\tau H}S^{+}}\\\
&-\expectationvalue{e^{\tau H}S^{-}e^{-\tau H}S^{-}}-\expectationvalue{e^{\tau
H}S^{+}e^{-\tau H}S^{+}}.\end{split}$

Since $H$ commutes with $S^{z}$ only the $+-$ and $-+$ contributions are non-
zero, while the $++$ and $--$ contributions are equal to zero, which implies
the $\Gamma$-point equality $\chi^{xx}(0,\omega)=\chi^{yy}(0,\omega)$.

## Appendix B

The spin-orbit coupling mimics that in the $t_{2g}$ subspace of atomic
$d$-shell Du  _et al._ (2013) commonly found in real materials (with cubic
site-symmetry). In the on-site basis
${\\{1\\!\\!\uparrow,1\\!\\!\downarrow,2\\!\\!\uparrow,2\\!\\!\downarrow,3\\!\\!\uparrow,3\\!\\!\downarrow\\}}$
it reads

$h^{\text{SOC}}=-\dfrac{\zeta}{2}\begin{pmatrix}0&0&-i&0&0&1\\\
0&0&0&i&-1&0\\\ i&0&0&0&0&-i\\\ 0&-i&0&0&-i&0\\\ 0&-1&0&i&0&0\\\
1&0&i&0&0&0\\\ \end{pmatrix}.$

The calculations presented in this work use $\zeta=0.1$.

## References

  * Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
  * Georges and Kotliar (1992) A. Georges and G. Kotliar, Phys. Rev. B 45, 6479 (1992).
  * Georges  _et al._ (1996) A. Georges, G. Kotliar, W. Krauth, and M. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ (2006) G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  * Held (2007) K. Held, Advances in Physics 56, 829 (2007).
  * Jarrell (1992) M. Jarrell, Phys. Rev. Lett. 69, 168 (1992).
  * Boehnke  _et al._ (2011) L. Boehnke, H. Hafermann, M. Ferrero, F. Lechermann, and O. Parcollet, Phys. Rev. B 84, 075145 (2011).
  * Hafermann  _et al._ (2014) H. Hafermann, E. G. C. P. van Loon, M. I. Katsnelson, A. I. Lichtenstein, and O. Parcollet, Phys. Rev. B 90, 235105 (2014).
  * van Loon  _et al._ (2014) E. G. C. P. van Loon, H. Hafermann, A. I. Lichtenstein, A. N. Rubtsov, and M. I. Katsnelson, Phys. Rev. Lett. 113, 246407 (2014).
  * Stepanov  _et al._ (2018) E. A. Stepanov, S. Brener, F. Krien, M. Harland, A. I. Lichtenstein, and M. I. Katsnelson, Phys. Rev. Lett. 121, 037204 (2018).
  * Park  _et al._ (2011) H. Park, K. Haule, and G. Kotliar, Phys. Rev. Lett. 107, 137007 (2011).
  * Shinaoka  _et al._ (2020) H. Shinaoka, D. Geffroy, M. Wallerberger, J. Otsuki, K. Yoshimi, E. Gull, and J. Kuneš, SciPost Phys. 8, 012 (2020).
  * Wallerberger  _et al._ (2020) M. Wallerberger, H. Shinaoka, and A. Kauch, (2020), arXiv:2012.05557 .
  * Otsuki  _et al._ (2019) J. Otsuki, K. Yoshimi, H. Shinaoka, and Y. Nomura, Phys. Rev. B 99, 165134 (2019).
  * Krien (2019) F. Krien, Phys. Rev. B 99, 235106 (2019).
  * Katanin (2020) A. Katanin, Phys. Rev. B 101, 035110 (2020).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
  * van Loon  _et al._ (2015) E. G. C. P. van Loon, H. Hafermann, A. I. Lichtenstein, and M. I. Katsnelson, Phys. Rev. B 92, 085106 (2015).
  * Krien  _et al._ (2017) F. Krien, E. G. C. P. van Loon, H. Hafermann, J. Otsuki, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 96, 075155 (2017).
  * Krien (2018) F. Krien,  _Conserving dynamical mean field approaches to strongly correlated systems_ , Ph.D. thesis, University of Hamburg (2018).
  * Geffroy  _et al._ (2019) D. Geffroy, J. Kaufmann, A. Hariki, P. Gunacker, A. Hausoel, and J. Kuneš, Phys. Rev. Lett. 122, 127601 (2019).
  * Niyazi  _et al._ (2020) A. Niyazi, D. Geffroy, and J. Kuneš, Phys. Rev. B 102, 085159 (2020).
  * Goldston (1961) J. Goldston, Nuovo Cim. 19, 154 (1961).
  * Nambu and Jona-Lasinio (1961) Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961).
  * Kuneš (2014) J. Kuneš, Phys. Rev. B 90, 235140 (2014).
  * Kuneš and Augustinský (2014) J. Kuneš and P. Augustinský, Phys. Rev. B 90, 235112 (2014).
  * Nasu  _et al._ (2016) J. Nasu, T. Watanabe, M. Naka, and S. Ishihara, Phys. Rev. B 93, 205136 (2016).
  * Anderson (1952) P. W. Anderson, Phys. Rev. 86, 694 (1952).
  * Manousakis (1991) E. Manousakis, Rev. Mod. Phys. 63, 1 (1991).
  * Chakravarty  _et al._ (1989) S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys. Rev. B 39, 2344 (1989).
  * Sandvik (1997) A. W. Sandvik, Phys. Rev. B 56, 11678 (1997).
  * Rowe  _et al._ (2012) W. Rowe, J. Knolle, I. Eremin, and P. J. Hirschfeld, Phys. Rev. B 86, 134513 (2012).
  * Del Re and Toschi (2021) L. Del Re and A. Toschi, (2021), arXiv:2011.04080 .
  * Vilk  _et al._ (1994) Y. M. Vilk, L. Chen, and A.-M. S. Tremblay, Phys. Rev. B 49, 13267 (1994).
  * Bickers  _et al._ (1989) N. E. Bickers, D. J. Scalapino, and S. R. White, Phys. Rev. Lett. 62, 961 (1989).
  * Lüscher and Läuchli (2009) A. Lüscher and A. M. Läuchli, Phys. Rev. B 79, 195102 (2009).
  * Varney  _et al._ (2009) C. N. Varney, C.-R. Lee, Z. J. Bai, S. Chiesa, M. Jarrell, and R. T. Scalettar, Phys. Rev. B 80, 075116 (2009).
  * Kent  _et al._ (2005) P. R. C. Kent, M. Jarrell, T. A. Maier, and T. Pruschke, Phys. Rev. B 72, 060411 (2005).
  * Fuchs  _et al._ (2011) S. Fuchs, E. Gull, L. Pollet, E. Burovski, E. Kozik, T. Pruschke, and M. Troyer, Phys. Rev. Lett. 106, 030401 (2011).
  * Hirschmeier  _et al._ (2015) D. Hirschmeier, H. Hafermann, E. Gull, A. I. Lichtenstein, and A. E. Antipov, Phys. Rev. B 92, 144409 (2015).
  * Rohringer  _et al._ (2018) G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Rev. Mod. Phys. 90, 025003 (2018).
  * Rohringer  _et al._ (2011) G. Rohringer, A. Toschi, A. Katanin, and K. Held, Phys. Rev. Lett. 107, 256402 (2011).
  * Ulmke  _et al._ (1995) M. Ulmke, V. Janiš, and D. Vollhardt, Phys. Rev. B 51, 10411 (1995).
  * Kuneš (2011) J. Kuneš, Phys. Rev. B 83, 085102 (2011).
  * Bauer  _et al._ (2011) B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner, and S. Wessel, J. Stat. Mech. Theory Exp. 2011, P05001 (2011).
  * Shinaoka  _et al._ (2017) H. Shinaoka, E. Gull, and P. Werner, Comput. Phys. Commun. 215, 128 (2017).
  * Gaenko  _et al._ (2017) A. Gaenko, A. E. Antipov, G. Carcassi, T. Chen, X. Chen, Q. Dong, L. Gamper, J. Gukelberger, R. Igarashi, S. Iskakov, M. Könz, J. P. LeBlanc, R. Levy, P. N. Ma, J. E. Paki, H. Shinaoka, S. Todo, M. Troyer, and E. Gull, Comput. Phys. Commun. 213, 235 (2017).
  * Werner  _et al._ (2006) P. Werner, A. Comanac, L. de’ Medici, M. Troyer, and A. J. Millis, Phys. Rev. Lett. 97, 076405 (2006).
  * Gubernatis  _et al._ (1991) J. E. Gubernatis, M. Jarrell, R. N. Silver, and D. S. Sivia, Phys. Rev. B 44, 6011 (1991).
  * Kaufmann and Held (2021) J. Kaufmann and K. Held, (2021), arXiv:2105.11211 .
  * Note (1) In case of summation over the reciprocal lattice vectors this amounts to $\frac{1}{N}\sum@\slimits@_{\mathbf{k}}$, in case of summation over Matsubara frequencies to $T\sum@\slimits@_{\nu}$ where $T$ is the temperature.
  * Note (2) $F_{\nu}=1$ in the Matsubara frequency basis.
  * Otsuki  _et al._ (2017) J. Otsuki, M. Ohzeki, H. Shinaoka, and K. Yoshimi, Phys. Rev. E 95, 061302 (2017).
  * Rohringer and Toschi (2016) G. Rohringer and A. Toschi, Phys. Rev. B 94, 125144 (2016).
  * Schäfer  _et al._ (2021) T. Schäfer, N. Wentzell, F. Šimkovic, Y.-Y. He, C. Hille, M. Klett, C. J. Eckhardt, B. Arzhang, V. Harkov, F.-M. Le Régent, A. Kirsch, Y. Wang, A. J. Kim, E. Kozik, E. A. Stepanov, A. Kauch, S. Andergassen, P. Hansmann, D. Rohe, Y. M. Vilk, J. P. F. LeBlanc, S. Zhang, A.-M. S. Tremblay, M. Ferrero, O. Parcollet, and A. Georges, Phys. Rev. X 11, 011058 (2021).
  * Watanabe and Murayama (2012) H. Watanabe and H. Murayama, Phys. Rev. Lett. 108, 251602 (2012).
  * Note (3) We mention only the spin symmetries, which are broken or varied in our models and ignore the ones, which are not changed such the global $U(1)$ due to charge conservation.
  * Note (4) For the sake of simplicity we have omitted the pair-hopping term $c_{i1\uparrow}^{\dagger}c_{i1\downarrow}^{\dagger}c_{i2\downarrow}c_{i2\uparrow}+H.c.$ This term does not affect the broken spin symmetry and is expected to play a minor role when the local high-spin state dominates as in the AFM phases studied here.
  * Toschi  _et al._ (2007) A. Toschi, A. A. Katanin, and K. Held, Phys. Rev. B 75, 045118 (2007).
  * Kaufmann  _et al._ (2021) J. Kaufmann, C. Eckhardt, M. Pickem, M. Kitatani, A. Kauch, and K. Held, Phys. Rev. B 103, 035120 (2021).
  * Rubtsov  _et al._ (2008) A. N. Rubtsov, M. I. Katsnelson, and A. I. Lichtenstein, Phys. Rev. B 77, 033101 (2008).
  * Kim  _et al._ (2020) A. J. Kim, P. Werner, and R. Valentí, Phys. Rev. B 101, 045108 (2020).
  * Ahn  _et al._ (2017) K.-H. Ahn, K. Pajskr, K.-W. Lee, and J. Kuneš, Phys. Rev. B 95, 064416 (2017).
  * Du  _et al._ (2013) L. Du, L. Huang, and X. Dai, The European Physical Journal B 86 (2013).


# Single Image Super-Resolution via a Dual Interactive Implicit Neural Network

Quan H. Nguyen and William J. Beksi
The University of Texas at Arlington
Arlington, TX, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we introduce a novel implicit neural network for the task of
single image super-resolution at arbitrary scale factors. To do this, we
represent an image as a decoding function that maps locations in the image
along with their associated features to their reciprocal pixel attributes.
Since the pixel locations are continuous in this representation, our method
can refer to any location in an image of varying resolution. To retrieve an
image of a particular resolution, we apply a decoding function to a grid of
locations each of which refers to the center of a pixel in the output image.
In contrast to other techniques, our dual interactive neural network decouples
content and positional features. As a result, we obtain a fully implicit
representation of the image that solves the super-resolution problem at (real-
valued) elective scales using a single model. We demonstrate the efficacy and
flexibility of our approach against the state of the art on publicly available
benchmark datasets.

##  1 Introduction

Single image super-resolution (SISR) is a fundamental low-level computer
vision problem that aims to recover a high-resolution (HR) image from its low-
resolution (LR) counterpart. There are two main reasons for performing SISR:
(i) to enhance the visual quality of an image for human consumption, and (ii)
to improve the representation of an image for machine perception. SISR has
many practical applications including robotics, remote sensing, satellite
imaging, thermal imaging, medical imaging, and much more [40, 39]. Despite
being a challenging and ill-posed subject, SISR has remained a crucial area of
study in the research community.

Recent deep learning approaches have provided high-quality SISR results [3,
37]. In perception systems, images are represented as 2D arrays of pixels
whose quality, sharpness, and memory footprint are controlled by the
resolution of the image. Consequently, the scale of the generated HR image is
fixed depending on the training data. For example, if a neural network is
trained to recover HR images of $\times 2$ scale then it only performs well on
what it is trained for (_i.e_., performance will be poor on $\times 3$,
$\times 4$, or other scales). Thus, instead of training multiple models for
various resolutions, it can be extremely useful in terms of practicality to
have a single SISR architecture that handles arbitrary scale factors. This is
especially true for embedded vision platforms (e.g., unmanned ground/aerial
vehicles) with multiple on-board cameras that must execute difficult tasks
using limited computational resources.

Figure 1: Our proposed dual interactive implicit neural network (DIINN) is
capable of producing images of arbitrary resolution, using a single trained
model, by capturing the underlying implicit representation of the input image.

The notion of an implicit neural representation, also known as coordinate-
based representation, is an active field of research that has yielded
substantial results in modeling 3D shapes [2, 5, 18, 30, 31]. Inspired by
these successes, learning implicit neural representations of 2D images is a
natural solution to the SISR problem since an implicit system can produce
output at arbitrary resolutions. While this idea has been touched upon in
several works [4, 33, 6, 26], in this paper we propose a more expressive
neural network for SISR with significant improvements over the existing state
of the art, Figure 1. Our contributions are summarized as follows.

  * •

We develop a novel dual interactive implicit neural network (DIINN) for SISR
that handles image content features in a modulation branch and positional
features in a synthesis branch, while allowing for interactions between the
two.

  * •

We learn an implicit neural network with a pixel-level representation, which
allows for locally continuous super-resolution synthesis with respect to the
nearest LR pixel.

  * •

We demonstrate the effectiveness of our proposed network by setting new
benchmarks on public datasets.

Our source code is available at [10]. The remainder of this paper is organized
as follows. Related research is discussed in Section 2. In Section 3, we
present a detailed description of our model for SISR at arbitrary scales using
an implicit representation. Experimental results are presented in Section 4.
The paper concludes in Section 5 and discusses future work.

Figure 2: An overview of our DIINN architecture. Note that each layer in the
modulation network is immediately followed by a ReLU activation, (2) and (4).
Similarly, each layer in the synthesis network is followed by a sine
activation, (3) and (5).

##  2 Related Work

This section highlights pertinent literature on the task of SISR. First, we
discuss deep learning techniques for SISR. Then, we provide an overview of
implicit neural representations. Lastly, we cite the nascent domain of
implicit representations for images. The SISR problem is an ill-defined
problem in the sense that there are many possible HR images that can be
downsampled to a single LR image. In this work, we focus on learning
deterministic mappings, rather than stochastic mappings (_i.e_., generative
models). In general, the input to an SISR system is an LR image, and the
output is a super-resolved (SR) image that may or may not have the same
resolution as a target HR image.

###  2.1 Deep Learning for SISR

Existing work on SISR typically utilizes convolutional neural networks (CNNs)
coupled with upsampling operators to increase the resolution of the input
image.

####  2.1.1 Upscaling $+$ Refinement

SRCNN [11], VDSR [20], and DRCN [21] first interpolate an LR image to a
desired resolution using bicubic interpolation, followed by a CNN-based neural
network to enhance the interpolated image and produce an SR image. The
refining network acts as a nonlinear mapping, which aims to improve the
quality of the interpolation. These methods can produce SR images at arbitrary
scales, but the performance is severely affected by the noise introduced
during the interpolation process. The refining CNNs also have to operate at
the desired resolution, thus leading to a longer runtime.

####  2.1.2 Learning Features $+$ Upscaling

Methods following this approach first feed an LR image through a CNN to obtain
a deep feature map at the same resolution. In this way, the CNNs are of
cheaper cost since they are applied at LR, which allows for deeper
architectures. Next, an upscaling operator is used to produce an SR image. The
most common upscaling operators are deconvolution (FSRCNN [12], DBPN [14]),
and sub-pixel convolution (ESPCN [32], EDSR [23]). It is also possible to
perform many iterations of learning features $+$ upscaling and explicitly
exploit the relationship between intermediate representations [14]. These
methods only work with integer scale factors and produce fixed-sized outputs.

EDSR [23] attempts to mitigate these problems by training a separate upscaling
head for each scale factor. On the other hand, Meta-SR [15] is among the first
attempts to solve SISR at arbitrary real-valued scale factors via a soft
version of the sub-pixel convolution. To predict the signal at each pixel in
the SR image, Meta-SR uses a meta-network to determine the weights for
features of a $(3\times 3)$ window around the nearest pixel in the LR image.
Effectively, each channel of the predicted pixel in the SR image is a weighted
sum of a $(C\times 3\times 3)$ volume, where $C$ is the number of channels in
the deep feature map. While Meta-SR has a limited generalization capability to
scale factors larger than its training scales, it can be viewed as a hybrid
implicit/explicit model.

###  2.2 Implicit Neural Representations

Implicit neural representations are an elegant way to parameterize signals
continuously in comparison to conventional representations, which are usually
discrete. Chen _et al_. [7], Mescheder _et al_. [27], and Park _et al_. [28]
are among the first to show that implicit neural representations outperform 3D
representations (_e.g_., meshes, voxels, and point clouds) in 3D modeling.
Many works that achieve state-of-the-art results in 3D computer vision have
followed. For example, Chabra _et al_. [5] learned local shape priors for the
reconstruction of 3D surfaces coupled with a deep signed distance function. A
new implicit representation for 3D shape learning called a neural distance
field was proposed by Chibane _et al_. [9]. Jiang _et al_. [18] leveraged
voxel representations to enable implicit functions to fit large 3D scenes, and
Peng _et al_. [30] increased the expressiveness of 3D scenes with various
convolutional models. It also is possible to condition the implicit neural
representations on the input signals [5, 8, 18, 30], which can be considered
as a hybrid implicit/explicit model.

####  2.2.1 Implicit Neural Representations of 2D Images

Despite there being many uses of implicit neural representations in 3D
computer vision, to the extent of our knowledge, it has been under-explored in
the domain of 2D imaging. Parameterizing images implicitly via neural networks
can be traced back to Stanley _et al_. [34] in 2007. Various types of signals,
including images, are implicitly represented with neural networks using
periodic activation functions by Sitzmann _et al_. [33]. Bemana _et al_. [4]
learned representations of changes in view, time, and light using Jacobians of
pixel positions to naturally interpolate images. Tancik _et al_. [35] showed
that a coordinate-based multilayer perceptron (MLP) can implicitly synthesize
2D images with sharpness by transforming the coordinate inputs with a Fourier
feature mapping. These works have limited generalizability to different
instances, especially as the complexity of the data increases. Recently, Mehta
_et al_. [26] proposed to modulate a synthesis network with a separate network
whose inputs are local features, which enables generalization while
maintaining high fidelity.

####  2.2.2 Implicit Neural Representations for SISR

Learning implicit neural representations of 2D images is immediately useful
for SISR since it can sample the pixel signals at any location in the spatial
domain. Chen _et al_. [6] proposed a local implicit image representation in
which 2D deep feature maps of the input images coupled with an ensemble of
local predictions ensure a smooth (_i.e_., continuous) transition between
different locations. Ma _et al_. [24] extended [6] by enforcing sharpness
constraints and additionally imposing multiple loss functions (including L1,
perceptual [19], and generative adversarial [13] losses) versus the L1 loss in
[6] to generate perceptually-pleasant details. In contrast to [6], we improve
the implicit decoding function by decoupling the content and the positional
features as proposed by [26] at pixel-level representations.

Figure 3: An image at  $(2,2)$ and  $(4,4)$ resolutions. We consider pixels as
squares rather than points, the global and local coordinates reference the
centers of pixel squares with respect to the center of the image and the
center of the nearest LR pixel, respectively. Both are scaled to $[-1,1]^{2}$.

##  3 Dual Interactive Implicit Neural Network

In this section, we present the core components of our proposed model for
SISR. The network consists of the following two parts: an encoder and a dual
interactive implicit decoder, Figure 2. The encoder learns the content of the
LR image and produces a deep feature map. The implicit decoder predicts the
signal (_i.e_., the $(r,g,b)$ values) at any query location within the image
space, conditioned on the associated features provided by the encoder. For a
target resolution, we query the signal of every pixel in the SR output image.
Our network can be expressed as

$S_{(x,y)}=f_{\theta}(g_{\gamma}(L)_{(x,y)},p),$ (1)

where $S_{(x,y)}$ is the SR predicted signal at $(x,y)$, $L$ is the LR image,
$g_{\gamma}(\cdot)$ is the encoder function parameterized by $\gamma$,
$g_{\gamma}(L)_{(x,y)}$ is the extracted features at $(x,y)$ (referred to as
content features), $p$ is the positional features, and
$f_{\theta}(\cdot,\cdot)$ is the implicit decoder function parameterized by
$\theta$. Note that our method considers pixels of an image as squares rather
than points and the location of a pixel references its center.

###  3.1 Encoder

The encoder supplies the decoder with a content feature representation
associated with each query location within the image. The encoder used in our
method could be a CNN, similar to those used in previous works [6, 15],
without any upscaling module. The input is an LR image and its output is a
deep feature map that preserves the spatial content of the image. With proper
settings (_e.g_., kernel size, padding, _etc_.), the output has the same
spatial size as the input image. Thus, given a set of coordinates for query
locations, we can sample its corresponding features via interpolation.

####  3.1.1 Feature Unfolding

Following previous works [6, 15], before interpolation we enrich the feature
map by applying feature unfolding with a kernel size of $3$. That is, we
concatenate the features of a $3\times 3$ neighborhood around each pixel of
the deep feature map, increasing the number of channels by 9. Since our
implicit decoder handles continuous queries and the predictions are locally
continuous, we choose nearest interpolation to avoid the additional smoothness
constraint as well as artifacts that come with bilinear or bicubic
interpolation. For any query location, the encoder effectively supplies the
decoder with the nearest features in the deep feature map. As expressed in
(1), $g_{\gamma}(L)_{(x,y)}$ learns a feature map of $L$ and allows for
querying content features (via nearest interpolation) associated with a
coordinates $(x,y)$.

###  3.2 Decoder

The decoder predicts the signal at every pixel in the SR image at a target
resolution. Our implicit decoder uses both the content features (from the
encoder) and the positional features for prediction. For SISR, the positional
features of a query location typically contain information about its relative
position to the nearest LR pixel (center) and some information about the scale
factor [15]. We refer to the coordinates that reference the position of a
pixel with respect to the center of the image as the global coordinates.
Conversely, we refer to the coordinates that reference the position of a pixel
with respect to its nearest LR pixel as the local coordinates. In our method,
the global coordinates allow for uniquely identifying each pixel and computing
its nearest LR pixel while the local ones are used as input to the decoder.
Figure 3 shows examples of global and local coordinates.

####  3.2.1 Modulated Periodic Activations Network

Recently, neural networks with periodic activation functions [33, 26] exhibit
excellent performance in reconstructing high-fidelity images and videos. We
use a similar dual MLP architecture, also known as a modulated periodic
activations neural network [26], for our decoder with adjustments to address
the SISR task. Our decoder of $N$ layers (_i.e_., $N$ layers of modulation and
$N$ layers of synthesis) can be written recursively as

$\displaystyle m_{0}$ $\displaystyle=ReLU(w_{0}z+b_{0}),$ (2) $\displaystyle
s_{0}$ $\displaystyle=m_{0}\odot sin(w^{\prime}_{0}p+b^{\prime}_{0}),$ (3)
$\displaystyle m_{i}$ $\displaystyle=ReLU(w_{i}[s_{i-1}\;z]+b_{i}),$ (4)
$\displaystyle s_{i}$ $\displaystyle=m_{i}\odot
sin(w^{\prime}_{i}s_{i-1}+b^{\prime}_{i}),$ (5)

where $w$, $w^{\prime}$, $b$, and $b^{\prime}$ are the weights and biases, $z$
is the content features, and $p$ is the positional features. The last output
of $s_{i}$ is then passed through a final dense layer (without activation) to
output the predicted SR signal.

The two main differences between our method and [26] are the following.

  1. (i)

The $z$ vector in the image reconstruction experiment in [26] represents a
patch while ours represents a pixel in the LR image. Typically, coarse-grain
features are associated with higher levels of semantic information and thus
may miss low-level details (_e.g_., edges, corners, _etc_.), which are crucial
for SISR. Additionally, pixel-level representations have consistently
performed well in the SISR literature [23, 15, 6]. Therefore, we opt for the
finest-grain representation (_i.e_., pixel-level). The architecture of our
encoder reflects this choice.

  2. (ii)

In (4), we use a concatenation of the output of the previous synthesis layer
(instead of the previous modulation layer in [26]) and the latent feature
vector as inputs for the latter modulation layer. We argue that the outputs of
the synthesis network are progressively more refined toward different query
locations and therefore provide better information to the modulation network
while the latent feature vector serves as residual feedback. We show the
benefits of our modification over [26] in Section 4.

Dataset Scale $\times 2$ $\times 3$ $\times 4$ Method PSNR SSIM LR-PSNR PSNR
SSIM LR-PSNR PSNR SSIM LR-PSNR DIV2K Bicubic 31.06 0.8937 40.38 28.26 0.8138
39.45 26.7 0.753 38.68 Meta-SR [15] 34.64 0.938 57.13 30.91 0.8736 54.93 28.89
0.8173 53.19 LIIF [6] 34.46 0.9367 53.63 30.81 0.8724 52.23 28.88 0.8163 50.95
DIINN (ours) 34.63 0.938 55.96 30.93 0.8741 53.65 28.98 0.8193 52 B100 Bicubic
28.27 0.8316 38.54 25.88 0.7171 38.19 24.64 0.6411 37.93 Meta-SR [15] 29.23
0.8651 41.5 27.51 0.7823 48.6 25.94 0.701 47.36 LIIF [6] 30.66 0.8891 52.27
27.68 0.7862 51.67 26.17 0.7097 50.61 DIINN (ours) 30.69 0.8896 52.62 27.73
0.7873 52.36 26.22 0.7119 50.82 Set5 Bicubic 31.81 0.9097 40.77 28.63 0.8385
38.79 26.7 0.7739 37.29 Meta-SR [15] 35.68 0.9439 57.21 31.6 0.8878 48.82
29.95 0.8605 53.79 LIIF [6] 35.5 0.9427 54.06 32.15 0.9008 52.61 29.92 0.8601
50.52 DIINN (ours) 35.67 0.9438 56.38 32.26 0.902 54.2 30.06 0.8631 51.89
Set14 Bicubic 28.33 0.8437 38.25 25.74 0.7413 37.18 24.24 0.6648 36.57 Meta-SR
[15] 30.9 0.8897 52.93 27.65 0.8005 46.51 26.25 0.7383 50.61 LIIF [6] 31.15
0.8919 51.57 28.04 0.8083 50.43 26.34 0.7407 49.22 DIINN (ours) 31.29 0.8937
53.92 28.14 0.8101 51.53 26.43 0.7437 50.33 Urban100 Bicubic 25.44 0.8284
35.54 23.01 0.7151 35.02 21.69 0.6334 34.64 Meta-SR [15] 29.32 0.907 47.53
25.53 0.8152 42.79 24.12 0.7507 48.35 LIIF [6] 30.02 0.9147 48.92 26.29 0.8296
48.34 24.29 0.7555 47.33 DIINN (ours) 30.29 0.9176 51.2 26.46 0.8337 49.67
24.49 0.7624 48.75

Table 1: A comparison against state-of-the-art SISR methods that allow for
arbitrary scale upsampling at trained scales.

Scale $\times 2.5$ $\times 3.5$ Method PSNR SSIM LR-PSNR PSNR SSIM LR-PSNR
Bicubic 29.41 0.8514 39.89 27.39 0.781 39.02 Meta-SR [15] 31.36 0.8898 46.1
29.21 0.8305 47.18 LIIF [6] 32.29 0.9036 52.93 29.73 0.8427 51.62 DIINN (ours)
32.34 0.9046 53.43 29.82 0.845 52.39

Table 2: A comparison against state-of-the-art SISR methods that allow for
arbitrary scale upsampling on the DIV2K dataset at “inter”-scales.

Compared to the closest competing framework, LIIF [6], our method uses a more
expressive neural network for the decoder, which leads to better performance.
Unlike LIIF where the content features ($z$) and positional features ($p$) are
concatenated and fed to a single-branch decoder, we decouple these features
and use two separate branches while allowing for interactions between the two
branches. We show in Section 4 that our method is both better and faster.

###  3.3 Architecture Details

Note that while we present our method in a point-wise manner, in practice we
implement the model to process the whole image. We use RDN [41], without their
upsampling module, as the encoder. Our decoder consists of two 4-layer MLPs
(each with 256 hidden units) as described in Section 3.2. We implement the
MLPs as convolutional layers with 256 kernels of size 1. As discussed in
Section 3.1, after feeding the encoder with the LR image, we obtain a deep
feature map of the same spatial size with 64 channels. Then, we perform the
nearest interpolation thus effectively increasing the spatial size of the deep
feature map to be the same as the target resolution.

The right-hand side of (2) takes this upsampled deep feature map (denoted $z$)
as input. Then, we construct a 2D grid of global coordinates and compute the
respective local coordinates, which results in a tensor of the same spatial
size as the target resolution with 2 channels for $x$ and $y$. We additionally
attach to it $\frac{1}{upscale\\_ratio}$ and denote the concatenation as $p$.
(3) takes this concatenated tensor as input. Finally, we pass the output of
the synthesis network ($s_{N-1}$) through a convolutional layer with a kernel
size of 1 to produce the predicted SR image.

Scale $\times 6$ $\times 8$ $\times 10$ $\times 15$ Method PSNR SSIM LR-PSNR
PSNR SSIM LR-PSNR PSNR SSIM LR-PSNR PSNR SSIM LR-PSNR Bicubic 24.87 0.6761
37.75 23.75 0.6326 37.14 22.95 0.6056 36.72 21.63 0.57 36.02 Meta-SR [15]
26.53 0.7294 49.51 25.13 0.6704 46.64 24.18 0.633 45.16 22.63 0.5819 43.27
LIIF [6] 26.65 0.7368 49.3 25.32 0.6865 48.3 24.38 0.6527 47.55 22.82 0.6038
46.49 DIINN (ours) 26.74 0.7404 50.18 25.41 0.6898 49.17 24.46 0.6557 48.48
22.89 0.6063 47.55

Table 3: A comparison against state-of-the-art SISR methods that allow for
arbitrary scale upsampling on the DIV2K dataset at “outer”-scales.

Output size $(128\times 128)$ $(256\times 256)$ $(512\times 512)$ Method
Runtime (ms) Bicubic 0.11 0.15 0.49 Meta-SR [15] 50.27 69.27 162.22 LIIF [6]
73.36 185.42 693.27 DIINN (ours) 56.6 95.51 243.59

Table 4: A runtime comparison in milliseconds. The input is a single RGB image
of size $(48\times 48)$. We report the average runtime of the forward pass for
each method and SR size over 100 runs on a single NVIDIA Quadro RTX 3000.
Bicubic is for reference only.

##  4 Experiments

In this section, we present an experimental evaluation of DIINN for SISR. We
give an overview of the datasets and metrics used in Section 4.1, along with
the training details in Section 4.2. In Section 4.3, we highlight our
benchmark results on popular datasets and provide a comparison with state-of-
the-art methods. We discuss and illustrate our qualitative results on various
image scales in Section 4.4. Finally, we present an ablation study in Section
4.5.

###  4.1 Datasets and Metrics

The DIV2K dataset [1], released for the NTIRE 2017 challenge on SISR [36],
consists of 1000 HR images each of which has a height or width equal to 2040.
Our model is trained using a split of 800 HR images. We do not retrain and we
preserve all hyperparameters when testing the model on the following four
standard benchmark datasets: DIV2K validation set (100 HR images), Set5,
Set14, BSDS100 [25], and Urban100 [16]. The results are evaluated using the
peak signal-to-noise ratio (PSNR) [17] and the structural similarity index
measure (SSIM) [38]. Additionally, we evaluate the LR consistency by measuring
the LR-PSNR. LR-PSNR is computed as the PSNR between the downsampled SR image
and the downsampled ground-truth HR image with the same bicubic kernel.

###  4.2 Training Details

Our training procedure is similar to [15, 6]. For each scale factor
$s\in\\{2,3,4\\}$ and HR image in the minibatch, we randomly cropped a
$48s\times 48s$ patch and downsample it using bicubic interpolation via the
resize function available in Torchvision [29]. We randomly applied horizontal,
vertical, and/or diagonal flips, each with a probability of 0.5. A minibatch
of 4 HR images was used, which resulted in 12 pairs of LR and HR images across
the three scales. We trained our model for 1000 epochs, where each epoch is a
full pass through 800 HR images in the DIV2K training set, using the Adam [22]
optimizer with the default hyperparameters provided by PyTorch [29]. The
learning rate was initialized at $10^{-4}$ and halved every 200 epochs. To
ensure a fair comparison, we followed previous works [15, 6] and used the L1
loss to train our network.

###  4.3 Benchmark Results

We compared DIINN against other methods that allow for arbitrary scale
upsampling, namely, Meta-SR [15] and LIIF [6]. We retrained Meta-SR and LIIF
under the same settings. In Tables 1-3, we show the results on the upsampling
scales that the models were trained with (_e.g_., $\times 2,\times 3,$ and
$\times 4$), “inter”-scales (_e.g_., $\times 2.5$ and $\times 3.5$), and the
“outer”-scales (_e.g_., $\times 6,\times 8,\times 10,$ and $\times 15$). The
HR images were downsampled by each scale to be used as inputs and the SR
images are of the same size as the HR images.

For the scales that the models were trained with (Table 1), we observe that in
terms of PSNR and SSIM, DIINN only underperforms Meta-SR in two settings
($\times 2$ scale for DIV2K and Set5) while it outperforms both Meta-SR and
LIIF in all other settings. For both “inter”-scales (Table 2) and
“outer”-scales (Table 3), DIINN consistently beats both Meta-SR and LIIF
across all settings. In addition, we report the inference time of all the
methods with the same input and output sizes in Table 4. The results confirm
that our framework is competitive with the state of the art in both
performance and generalizability.

###  4.4 Qualitative Results

Figure 4: Qualitative results on an image patch from the Urban100 dataset.

In Figure 4, we show a qualitative comparison between bicubic interpolation
(first column), Meta-SR (second column), LIIF (third column), and DIINN (last
column). SR predictions are obtained across all methods from the same LR input
(downsampled via bicubic interpolation by a factor of $2\pi\approx 6.28$) at
increasing scale factors. We denote the size of the images for each row and
the method for each column. As expected, bicubic interpolation only smooths
the LR image and is unable to recover sharp details. SR images predicted by
Meta-SR have artifacts, especially around the boundary locations corresponding
to the edges of LR pixels. We observe significantly better results from LIIF
and DIINN, which shows the benefits of implicit decoders for SISR.

Compared to LIIF, DIINN produces sharper edges and details. Note that if we
zoom in on the results of our method, we can see similar artifacts to those of
Meta-SR, but much less distinguishable. LIIF overcomes these types of
artifacts by averaging their predictions over a small window, which is termed
“local ensemble”. However, as we observe in the third column in Figure 4, the
“local ensemble” introduces blunter edges at a higher computational cost. We
also note that the aforementioned artifacts do not appear visually with DIINN
when the LR input is of higher quality (_e.g_., downsampled HR image at a
factor of $4$ or less).

Scale $\times 3.14$ $\times 4$ $\times 8$ Model MI IP PSNR SSIM LR-PSNR PSNR
SSIM LR-PSNR PSNR SSIM LR-PSNR (a) $[m]$ Yes 25.67 0.8105 47.5 24.07 0.7465
47.25 20.68 0.5643 45.75 (b) $[m]$ No 25.78 0.814 47.93 24.17 0.7506 48.02
20.74 0.5673 46.42 (c) $[m\;z]$ Yes 25.89 0.8168 48.73 24.26 0.7538 48.26
20.77 0.5704 46.26 (d) $[m\;z]$ No 26.07 0.8218 49.11 24.42 0.7598 48.88 20.87
0.5762 47.05 (e) $[s\;z]$ Yes 25.93 0.8179 48.58 24.3 0.7555 48.12 20.8 0.5725
45.95 (f) $[s\;z]$ No 26.14 0.8233 49.4 24.49 0.7624 48.75 20.91 0.5798 46.23

Table 5: A quantitative ablation study on the Urban100 dataset. “MI” stands
for the modulation inputs (_i.e_., the inputs to the subsequent modulation
layers) and “IP” stands for the initialization of the positional features
($p$).

###  4.5 Ablation Study

To more closely examine the design choices of our architecture, we performed
an ablation study on variations of the dual implicit decoder. In Table 5, we
summarize the quantitative results of six variations on the Urban100 dataset.
Note that model (f) is our final model, which we used to report the results in
the previous tables.

####  4.5.1 Modulation Inputs

First, we adjust the inputs to the modulation network in following three ways.

  1. 1.

For models (a) and (b), each subsequent modulation layer takes only the output
of its previous modulation layer as input and the residual connections to the
content features ($z$) are removed. The following equation is used in place of
(4):

$m_{i}=ReLU(w_{i}m_{i-1}+b_{i}).$ (6)

  2. 2.

For models (c) and (d), each subsequent modulation layer takes a concatenation
of the output of its previous modulation layer and the content features as
input, as proposed in [26]. The following equation is used in place of (4):

$m_{i}=ReLU(w_{i}[m_{i-1}\;z]+b_{i}).$ (7)

  3. 3.

For models (e) and (f), each subsequent modulation layer takes a concatenation
of the output of the previous synthesis layer and the content features as
input, as expressed in (4).

Comparing models (a) and (c), or models (b) and (d), we observe that the
addition of the skip connections to the content features ($z$) increases the
SISR performance similar to the observation in [26] on image reconstruction.
This performance boost comes with additional parameters in the modulation
networks. Nonetheless, comparing models (c) and (e), or models (d) and (f), we
find that by simply using the output of the previous synthesis layer instead
of the previous modulation layer, we can increase the PSNR and SSIM
performance at no additional cost. We note that the consistency of the SR
images with respect to the LR inputs, measured by the LR-PSNR, is better with
models (c) and (d) compared to models (e) and (f) at larger scale factors.

####  4.5.2 Positional Features

Next, we explore a way to connect the positional features with the
neighborhood around the corresponding LR pixels introduced by the feature
unfolding operation (Section 3.1). We attempt to do this by adding an
initializing dense layer with a sine activation at the beginning of the
decoder (_i.e_., before we perform (2)). The following equations are used to
transform both the positional and the content features:

$\displaystyle p$ $\displaystyle\leftarrow\sin(w_{i}p+b_{i}),$ (8)
$\displaystyle z$ $\displaystyle\leftarrow p\odot z.$ (9)

Intuitively, we want this initializing layer to decide the weights for each LR
pixel in the corresponding neighborhood, which could be seen as a learnable
distance weighting function. As shown in Table 5, we observe no improvement
across the three variations (Section 4.5.1) despite allowing for more
flexibility with additional parameters in our model.

##  5 Conclusion

In this paper, we leveraged the learning of implicit representations to
develop a dual interactive implicit network for SISR. Our approach allows for
arbitrary scale factors without the need to train multiple models. Through
extensive experimental evaluations across many settings and datasets, we
empirically showed that DIINN achieves state-of-the-art results. Exploring
ways to transform our architecture to a larger receptive field along with
improved interactions between the encoder and the decoder will be left for
future work. We will also consider adapting our method to learn a space of
predictions to address the ill-posed nature of SISR.

## Acknowledgments

The authors acknowledge the Texas Advanced Computing Center (TACC) at the
University of Texas at Austin for providing software, computational, and
storage resources that have contributed to the research results reported
within this paper.

## References

  * [1] Eirikur Agustsson and Radu Timofte.  Ntire 2017 challenge on single image super-resolution: Dataset and study.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pages 126–135, 2017.
  * [2] Matan Atzmon and Yaron Lipman.  Sal: Sign agnostic learning of shapes from raw data.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 2565–2574, 2020.
  * [3] K Bai, X Liao, Q Zhang, X Jia, and S Liu.  Survey of learning based single image super-resolution reconstruction technology.  Pattern Recognition and Image Analysis, 30(4):567–577, 2020.
  * [4] Mojtaba Bemana, Karol Myszkowski, Hans-Peter Seidel, and Tobias Ritschel.  X-fields: Implicit neural view-, light-and time-image interpolation.  ACM Transactions on Graphics, 39(6):1–15, 2020.
  * [5] Rohan Chabra, Jan E Lenssen, Eddy Ilg, Tanner Schmidt, Julian Straub, Steven Lovegrove, and Richard Newcombe.  Deep local shapes: Learning local sdf priors for detailed 3d reconstruction.  In Proceedings of the European Conference on Computer Vision, pages 608–625. Springer, 2020.
  * [6] Yinbo Chen, Sifei Liu, and Xiaolong Wang.  Learning continuous image representation with local implicit image function.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 8628–8638, 2021.
  * [7] Zhiqin Chen and Hao Zhang.  Learning implicit fields for generative shape modeling.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5939–5948, 2019.
  * [8] Julian Chibane, Thiemo Alldieck, and Gerard Pons-Moll.  Implicit functions in feature space for 3d shape reconstruction and completion.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6970–6981, 2020.
  * [9] Julian Chibane, Aymen Mir, and Gerard Pons-Moll.  Neural unsigned distance fields for implicit function learning.  In Proceedings of the Advances in Neural Information Processing Systems, pages 21638–21652, 2020.
  * [10] https://github.com/robotic-vision-lab/Dual-Interactive-Implicit-Neural-Network.
  * [11] Chao Dong, Chen Change Loy, Kaiming He, and Xiaoou Tang.  Image super-resolution using deep convolutional networks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 38(2):295–307, 2015.
  * [12] Chao Dong, Chen Change Loy, and Xiaoou Tang.  Accelerating the super-resolution convolutional neural network.  In Proceedings of the European Conference on Computer Vision, pages 391–407. Springer, 2016.
  * [13] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Proceedings of the Advances in Neural Information Processing Systems, 27, 2014.
  * [14] Muhammad Haris, Gregory Shakhnarovich, and Norimichi Ukita.  Deep back-projection networks for super-resolution.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1664–1673, 2018.
  * [15] Xuecai Hu, Haoyuan Mu, Xiangyu Zhang, Zilei Wang, Tieniu Tan, and Jian Sun.  Meta-sr: A magnification-arbitrary network for super-resolution.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1575–1584, 2019.
  * [16] Jia-Bin Huang, Abhishek Singh, and Narendra Ahuja.  Single image super-resolution from transformed self-exemplars.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 5197–5206, 2015.
  * [17] Michal Irani and Shmuel Peleg.  Motion analysis for image enhancement: Resolution, occlusion, and transparency.  Journal of Visual Communication and Image Representation, 4(4):324–335, 1993.
  * [18] Chiyu Jiang, Avneesh Sud, Ameesh Makadia, Jingwei Huang, Matthias Nießner, Thomas Funkhouser, et al.  Local implicit grid representations for 3d scenes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6001–6010, 2020.
  * [19] Justin Johnson, Alexandre Alahi, and Li Fei-Fei.  Perceptual losses for real-time style transfer and super-resolution.  In Proceedings of the European Conference on Computer Vision, pages 694–711. Springer, 2016.
  * [20] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee.  Accurate image super-resolution using very deep convolutional networks.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1646–1654, 2016.
  * [21] Jiwon Kim, Jung Kwon Lee, and Kyoung Mu Lee.  Deeply-recursive convolutional network for image super-resolution.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1637–1645, 2016.
  * [22] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In Proceedings of the International Conference on Learning Representations, 2015.
  * [23] Bee Lim, Sanghyun Son, Heewon Kim, Seungjun Nah, and Kyoung Mu Lee.  Enhanced deep residual networks for single image super-resolution.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 136–144, 2017.
  * [24] Cheng Ma, Peiqi Yu, Jiwen Lu, and Jie Zhou.  Recovering realistic details for magnification-arbitrary image super-resolution.  IEEE Transactions on Image Processing, 2022.
  * [25] David Martin, Charless Fowlkes, Doron Tal, and Jitendra Malik.  A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, volume 2, pages 416–423, 2001.
  * [26] Ishit Mehta, Michaël Gharbi, Connelly Barnes, Eli Shechtman, Ravi Ramamoorthi, and Manmohan Chandraker.  Modulated periodic activations for generalizable local functional representations.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 14214–14223, 2021.
  * [27] Lars Mescheder, Michael Oechsle, Michael Niemeyer, Sebastian Nowozin, and Andreas Geiger.  Occupancy networks: Learning 3d reconstruction in function space.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 4460–4470, 2019.
  * [28] Jeong Joon Park, Peter Florence, Julian Straub, Richard Newcombe, and Steven Lovegrove.  Deepsdf: Learning continuous signed distance functions for shape representation.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 165–174, 2019.
  * [29] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  Proceedings of the Advances in Neural Information Processing Systems, 32:8024–8035, 2019.
  * [30] Songyou Peng, Michael Niemeyer, Lars Mescheder, Marc Pollefeys, and Andreas Geiger.  Convolutional occupancy networks.  In Proceedings of the European Conference on Computer Vision, pages 523–540. Springer, 2020.
  * [31] Eduard Ramon, Gil Triginer, Janna Escur, Albert Pumarola, Jaime Garcia, Xavier Giro-i Nieto, and Francesc Moreno-Noguer.  H3d-net: Few-shot high-fidelity 3d head reconstruction.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5620–5629, 2021.
  * [32] Wenzhe Shi, Jose Caballero, Ferenc Huszár, Johannes Totz, Andrew P Aitken, Rob Bishop, Daniel Rueckert, and Zehan Wang.  Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1874–1883, 2016.
  * [33] Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  Proceedings of the Advances in Neural Information Processing Systems, 33:7462–7473, 2020.
  * [34] Kenneth O Stanley.  Compositional pattern producing networks: A novel abstraction of development.  Genetic Programming and Evolvable Machines, 8(2):131–162, 2007\.
  * [35] Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  Proceedings of the Advances in Neural Information Processing Systems, 2020.
  * [36] Radu Timofte, Eirikur Agustsson, Luc Van Gool, Ming-Hsuan Yang, and Lei Zhang.  Ntire 2017 challenge on single image super-resolution: Methods and results.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, pages 114–125, 2017.
  * [37] Wei Wang, Yihui Hu, Yanhong Luo, and Tong Zhang.  Brief survey of single image super-resolution reconstruction based on deep learning approaches.  Sensing and Imaging, 21(1):1–20, 2020.
  * [38] Zhou Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli.  Image quality assessment: From error visibility to structural similarity.  IEEE Transactions on Image Processing, 13(4):600–612, 2004.
  * [39] Jianchao Yang and Thomas Huang.  Image super-resolution: Historical overview and future challenges.  In Super-Resolution Imaging, pages 1–34. CRC Press, 2017.
  * [40] Linwei Yue, Huanfeng Shen, Jie Li, Qiangqiang Yuan, Hongyan Zhang, and Liangpei Zhang.  Image super-resolution: The techniques, applications, and future.  Signal Processing, 128:389–408, 2016.
  * [41] Yulun Zhang, Yapeng Tian, Yu Kong, Bineng Zhong, and Yun Fu.  Residual dense network for image super-resolution.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2018.