# Comparing quantum fluctuations in the spin-$\frac{1}{2}$ and spin-$1$ XXZ
Heisenberg models
on square and honeycomb lattices

Masahiro Kadosawa  Department of Physics, Chiba University, Chiba 263-8522,
Japan  Masaaki Nakamura  Department of Physics, Ehime University, Ehime
790-8577, Japan  Yukinori Ohta  Department of Physics, Chiba University, Chiba
263-8522, Japan  Satoshi Nishimoto  Department of Physics, Technical
University Dresden, 01069 Dresden, Germany  Institute for Theoretical Solid
State Physics, IFW Dresden, 01069 Dresden, Germany

###### Abstract

We present a detailed investigation of the XXZ Heisenberg model for spin-$1/2$
and spin-$1$ systems on square and honeycomb lattices. Utilizing the density-
matrix renormalization group (DMRG) method, complemented by Spiral Boundary
Conditions (SBC) for mapping two-dimensional (2D) clusters onto one-
dimensional (1D) chains, we meticulously explore the evolution of staggered
magnetization and spin gaps across a broad spectrum of easy-axis anisotropies.
Our study reveals that, despite the lower site coordination number of
honeycomb lattice, which intuitively suggests increased quantum fluctuations
in its Néel phase compared to the square lattice, the staggered magnetization
in the honeycomb structure exhibits only a marginal reduction. Furthermore,
our analysis demonstrates that the dependence of staggered magnetization on
the XXZ anisotropy $\Delta$, except in close proximity to $\Delta=1$, aligns
with series expansion predictions up to the 12th order. Notably, for the
$S=1/2$ honeycomb lattice, deviations from the 10th order series expansion
predictions near the isotropic Heisenberg limit emphasize the critical
influence of quantum fluctuations on the spin excitation in its Néel state.
Additionally, our findings are numerically consistent with the singular
behavior of the spin gap near the isotropic Heisenberg limit as forecasted by
spin-wave theory. The successful implementation of SBC marks a methodological
advancement, streamlining the computational complexity involved in analyzing
2D models and paving the way for more precise determinations of physical
properties in complex lattice systems.

##  I introduction

In the complex world of quantum magnetism, the interplay between spin
interactions and lattice geometry crafts a fascinating landscape of ground
states and excitations (e.g., see Refs. [1, 2, 3]). A central subject for this
exploration is the XXZ Heisenberg model [4], a cornerstone that has deeply
enriched our comprehension of anisotropic magnetic systems [5, 6]. The model,
celebrated for its versatility in representing real materials, allows for the
examination of quantum fluctuations – the quintessential quantum mechanical
effect that destabilizes classical magnetic order, paving the way for the
emergence of novel quantum phases like spin liquids [7, 8].

This paper focuses on comparing the manifestations of quantum fluctuations
within the $S=1/2$ and $S=1$ XXZ models on two fundamentally distinct lattice
structures: the square and the honeycomb. These lattices, emblematic of
different coordination environments and geometric constraints, provide a
compelling backdrop against which the interplay of spin magnitude and lattice
topology can be meticulously dissected. The square lattice, with its direct
links to high-temperature superconductivity [9] and magnetic order [3] in
solid-state compounds, and the honeycomb lattice, notable for hosting exotic
phenomena such as the quantum spin Hall effect [10] and potential quantum spin
liquid states [11], are ideal platforms for this comparative study.

Numerical simulations of such systems often pose significant challenges. When
instantiated on a finite-size lattice, the total degrees of freedom
exponentially increase with lattice size. This constraint on the geometry and
size of the cluster becomes particularly notable for systems in more than two
dimensions. Consequently, an extrapolation of physical quantity to an infinite
system size is imperative to ascertain the bulk value. However, the execution
of such finite-size scaling is usually fraught with challenges due to the
presence of multiple scaling dimensions, such as the x and y directions in a
2D case. In our previous studies [12, 13], we introduced an efficient
numerical method for determining the local order parameter in 2D systems
through the use of spiral boundary conditions (SBC). This method provides a
promising approach to address the challenges associated with finite-size
scaling in numerical simulations.

Applying SBC allows for the exact projection of lattice models, even those
extending beyond 2D, onto 1D periodic chains that maintain translational
symmetry. Within this projected 1D chain, each lattice site is denoted by a
single coordinate, contrary to the dual coordinates used in the original 2D
cluster. This simplification means that only one finite-size scaling analysis
along the chain direction is necessary to ascertain a physical quantity in the
thermodynamic limit. We have demonstrated the capability of precisely
determining the magnitude of staggered magnetization for the XXZ Heisenberg
model on a square lattice ranging from $S=1/2$ to $6$ [13]. We here consider
the extension of this technique to studies of the honeycomb-lattice model and
further demonstrates the systematic calculation of excitation energy in the
bulk limit.

In this paper, we investigate the $S=1/2$ and $S=1$ XXZ Heisenberg models on
square and honeycomb lattices employing the DMRG method. We delve into the
evolution of staggered magnetization and the accompanying spin gap with easy-
axis anisotropy. It is shown that by applying SBC to both lattice types, a
finite-size scaling analysis towards the thermodynamic limit can be
effortlessly conducted for the studied physical quantities. The efficacy of
our approach is corroborated by comparing our findings with pre-existing
numerical and analytical results. Our analysis reveals that for most cases,
the staggered magnetization and spin gap within the easy-axis Néel phase can
be approximately accounted for by series expansions (SE) of the 10th to 12th
order in terms of $1/\Delta$. However, for the $S=1/2$ honeycomb lattice, the
results significantly deviate from those of the 10th order SE across a broad
range near the isotropic Heisenberg limit due to strong quantum fluctuations.
Furthermore, for all models considered, we obtain results that are numerically
consistent with the singular behavior of the spin gap near the isotropic
Heisenberg limit as predicted by spin-wave theory (SWT). Also, we find a
marked reduction in quantum fluctuations transitioning from $S=1/2$ to $S=1$
across all physical quantities assessed.

The paper is structured as follows: Sec. II provides a detailed description of
our spin model. In Sec. III, we elucidate the method of mapping 2D models to
1D using SBC, along with the procedures for calculating physical quantities
via the DMRG technique. Sec. IV presents our numerical findings, examining the
influence of lattice type, XXZ anisotropy magnitude, and spin size on the
stability of staggered magnetization and the magnitude of the spin gap.
Additionally, we incorporate a discussion on the specific behavior of the spin
gap for the $S=1/2$ honeycomb-lattice case. Finally, in Sec. V, we conclude
the paper with a summary and further insights into the observed phenomena.

##  II model

The Hamiltonian of the XXZ Heisenberg model is represented as follows:

$\displaystyle\mathcal{H}=\sum_{\braket{ij}}(S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}+\Delta
S^{z}_{i}S^{z}_{j})$ (1)

where $S^{\gamma}_{i}$ $(\gamma=x,y,z)$ are the spin-$S$ operators, $\Delta$
is the anisotropy parameter, and the sum $\braket{ij}$ runs over all nearest-
neighbor pairs. In this context, we consider two types of lattice structures:
the square lattice and the honeycomb lattice. The XXZ models associated with
these lattices have been significantly examined thus far.

There are three phases depending on $\Delta$ [14, 15, 16, 17, 18]: (i) For
$\Delta>1$ easy-axis Néel phase with antiferromagnetic (AFM) spin alignment
along the $z$-direction, (ii) for $-1<\Delta<1$ easy-plane Néel (XY) phase
with AFM spin alignment along some arbitrary direction in the $xy$-plane, and
(iii) for $\Delta<-1$ ferromagnetic (FM) phase with fully-polarized spins
along the $z$-direction. The phase transitions at $\Delta=\pm 1$ are both
first order. For $\Delta=-1$ this model can be exactly solved: the Néel and FM
states are degenerate at the ground state where the energies are
$E_{0}=-2NS^{2}$ and $E_{0}=-(3/2)NS^{2}$ for square- and honeycomb-lattice
models. At $\Delta=-1^{+}$, the wave function of Néel state is expressed as

$\displaystyle|\Psi_{0}({\rm XY})\rangle=\sum_{m}\lambda_{m}|\psi_{m}\rangle,$
(2)

where $|\psi_{m}\rangle$ are bases restricted to $S^{z}_{\rm
tot}=\sum_{i=1}^{L}\langle S^{z}_{i}\rangle=0$ subspace, $m$ is summed over
all possible combinations of the spin configurations, and $\lambda_{m}$ are
determined for each $S$ [13]. The magnitude of staggered magnetization is $S$
with the direction parallel to the $xy$-plane. At $\Delta=-1^{-}$, the
function of FM state is

$\displaystyle|\Psi_{0}({\rm
FM})\rangle=\frac{1}{\sqrt{2}}(|\Uparrow\rangle+|\Downarrow\rangle),$ (3)

where $|\Uparrow\rangle$ and $|\Downarrow\rangle$ denote fully-polarized
states toward $z$ and $-z$ directions, respectively. These states
$|\Psi_{0}({\rm XY})\rangle$ and $|\Psi_{0}({\rm FM})\rangle$ are orthogonal.

In regards to the square lattice model, it has been numerically confirmed that
for $S=1/2$ [19, 13], Néel long-range order (LRO) always exists for
$\Delta>-1$. On the other hand, for the honeycomb-lattice model, due to fewer
bonds between adjacent sites compared to the square lattice, quantum
fluctuations are larger, and there is not yet a complete consensus on which
$S$ and $\Delta$ regions stabilize Néel LRO [20, 21].

##  III method

Figure 1:  (a) 2D square-lattice $4\times 4$ cluster and (b) 2D honeycomb-
lattice $3\times 3$ cluster. The regions outlined by red lines represent the
original clusters. Through the application of SBC, the 2D clusters are
projected onto 1D periodic chains, where the sites are aligned along the green
lines. (c) Schematic depiction of the 1D open chain utilized for our DMRG
calculations, focusing on staggered spin moment $m^{z}_{\alpha}(S,\Delta)$ and
spin gap $\varepsilon_{\alpha}(S,\Delta)$. In the square lattice, pinning is
applied to $L$ sites at each end, while in the honeycomb lattice, pinning is
applied to $2L$ sites. Upward and downward arrows signify pinned spins with
$S^{z}=S$ and $S^{z}=-S$, respectively.

###  III.1 spiral boundary conditions

Applying the DMRG method to 2D systems introduces significant challenges,
primarily due to two factors. Firstly, the entanglement entropy, which
quantifies the quantum correlations within different parts of the system,
follows an “area law”. This law indicates that entanglement entropy scales
with the boundary area of a region, complicating the simulation of large
systems. Secondly, the DMRG’s sweeping process, which optimizes quantum states
site by site in a linear fashion in 1D systems, encounters difficulties in the
more complex geometries of 2D systems. This complexity can lead to
inaccuracies because a straightforward sweeping mechanism is harder to
implement across 2D lattices.

To address these challenges, careful management of boundary conditions is
essential for accurate DMRG simulations. Traditional approaches often employ
cylinder or torus configurations for 2D systems. Yet, these configurations can
create short bond loops and impose an unnatural periodicity on the wave
function, leading to inaccuracies such as an unexpected plaquette constraint
on particles or spins. An inappropriate choice of boundary conditions might
also skew the energy states observed in finite clusters away from those
relevant in the thermodynamic limit, instead of systematic errors due to the
finite-size effects.

A promising alternative that circumvents these limitations involves the
implementation of SBC [22, 12]. SBC enables the exact projection of lattice
models, including those extending beyond two dimensions, onto 1D periodic
chains that preserve translational symmetry. This projection effectively
transforms a 2D $L\times L$ cluster into a 1D chain, maintaining nearest-
neighbor and ($L-1$)th-neighbor bonds for square lattices, and nearest- and
($2L-1$)th-neighbor bonds for honeycomb lattices, as depicted in Fig. 1(a,b).
This innovative approach prevents the emergence of artificial short bond loops
and ensures an even distribution of quantum entanglement across the chain,
leveraging translational symmetry.

Notably, SBC minimizes the distance of the longest bonds, denoted as $d$, to
$L-1$ for square and $2L-1$ for honeycomb lattices, optimizing conditions for
DMRG calculations. In contrast, conventional periodic boundary conditions
would increase $d$ to $2L$ and $4L-2$, respectively, posing challenges for
DMRG analysis.

Furthermore, SBC provides a significant benefit for finite-size scaling
analysis. By projecting the original 2D lattice onto a 1D chain, SBC enables
the indexing of each lattice site with a singular coordinate, rather than the
dual-coordinate system inherent to 2D clusters. This transition to a single-
coordinate framework simplifies the analytical process to a unidimensional
scale. It facilitates a more direct method for extrapolating physical
quantities to the thermodynamic limit, enhancing the accuracy and efficiency
of our simulations.

###  III.2 density-matrix renormalization group

The investigation of the ground state of the 1D chain, transformed via SBC, is
conducted using the DMRG method [23]. For this purpose, we implement open
boundary conditions on the 1D chain, a choice that significantly enhances the
precision of our DMRG calculations. Our study encompasses open chains with
lengths up to $N=L^{2}=196$ sites for the square lattice and up to
$N=2L^{2}=162$ sites for the honeycomb lattice. To ensure the robustness of
our calculations, we retain up to $m=8000$ density-matrix eigenstates, with
all calculated values subsequently extrapolated to the limit of
$m\rightarrow\infty$. The maximum discarded weight observed is on the order of
$10^{-6}$.

Furthermore, we intentionally break the spin-rotation symmetry by employing a
spin pinning technique. This approach effectively lifts the degeneracy of the
ground state, thereby efficiently reducing the dimensionality of the Hilbert
space required for our calculations. As a result, even for computations of 2D
systems, the discarded weight remains minimal, enhancing the accuracy and
feasibility of our analysis.

Figure 2:  Comparative profiles of the local spin moment for $12\times 12$
square (left panel) and $8\times 8$ honeycomb (right panel) lattices at spin-
isotropic case $\Delta=1$. The profiles depict $\braket{S^{z}_{i}}$, the
expectation value of the $z$-component of spin at site $i$, for two distinct
total spin sectors: $S^{z}=0$ (top row) and $S^{z}=1$ (middle row). The bottom
row illustrates the differential profiles, showing the variance in
$\braket{S^{z}_{i}}$ between the $S^{z}=0$ and $S^{z}=1$ sectors.

###  III.3 physical quantities

####  III.3.1 staggered magnetization

We examine the Néel state for $\Delta\geq 1$, where the magnitude of staggered
magnetization serves as the order parameter in both square and honeycomb
lattice XXZ models. In the original 2D clusters, Néel order characterized by
$\boldsymbol{k}=(\pi,\pi)$ is transformed into a magnetic order with $k=\pi$
along the projected 1D chain through our implementation of SBC. For the
analysis, we utilize open chains and determine the magnitude of staggered
magnetization by measuring half the amplitude of the Friedel oscillation of
$\braket{S^{z}_{i}}$ from the system edges. Given the condition $\Delta\geq
1$, it is sufficient to focus on the $z$-component of the spin moment.

To establish such an open chain configuration, we sever $L$ ($L+1$) bonds
between adjacent sites on the projected 1D periodic chain for square
(honeycomb) lattices, as described in recent works [13, 24]. We specifically
examine the local moments of the two central spins, $\braket{S^{z}_{N/2}}$ and
$\braket{S^{z}_{N/2+1}}$, employing spin pinning near the system edges, such
as setting $\braket{S^{z}_{i}}=(-1)^{i-1}S$ at sites $i=1,\dots,L$ and
$i=L^{2}-(L-1),\dots,L^{2}$ for the square lattice; at sites $i=1,\dots,2L$
and $i=2L^{2}-(2L-1),\dots,2L^{2}$ for honeycomb lattice.

While pinning is typically positioned at the system edges, i.e., at $i=1$,
$i=N$, the last $L-1$ ($2L-1$) sites at both ends of projected 1D chain
corresponding to square (honeycomb) lattice are left without their original
bonds, necessitating the placement of pinnings at these outer sites to
accurately estimate staggered magnetization [see Fig.1(c)]. Therefore, we
define the magnitude of staggered magnetization for a given spin $S$ and XXZ
anisotropy $\Delta$ as

$\displaystyle
m_{\alpha}^{z}(S,\Delta)=\lim_{N\to\infty}|\braket{S^{z}_{N/2}}-\braket{S^{z}_{N/2+1}}|/2.$
(4)

Examples of local spin moment profiles, $\braket{S^{z}_{i}}$, for the
$S^{z}=0$ sector at $\Delta=1$ are depicted in the top panels of Fig.2. Here,
the parameter $\alpha$ adopts the value ‘sq’ for square lattices and ‘hon’ for
honeycomb lattices. In both square and honeycomb lattices, the oscillation of
$\braket{S^{z}_{i}}$ smoothly decays towards the center of the system,
validating the approach of defining the Néel order parameter via the local
moments of the two central spins, $\braket{S^{z}_{N/2}}$ and
$\braket{S^{z}_{N/2+1}}$.

####  III.3.2 spin gap

The spin gap offers insight into phenomena such as quantum fluctuations and
the stability of the Néel state. It also provides us with understanding of how
classical behavior emerges from quantum systems. The spin gap (singlet-triplet
gap) for given $S$ and $\Delta$ is defined as follows:

$\displaystyle\varepsilon_{\alpha}(S,\Delta)=\lim_{N\to\infty}[E_{0}(N,1)-E_{0}(N,0)],$
(5)

where $E_{0}(N,S^{z})$ is the total ground-state energy of the system with $N$
sites and the $z$-component of the total spin, $S^{z}$. Here, the parameter
$\alpha$ adopts the value ‘sq’ for square lattices and ‘hon’ for honeycomb
lattices.

To verify the validity of our spin gap calculations under the imposed pinning
distribution, we examine the spatial distribution of a spinon as the spin
sector transitions from $S^{z}=0$ to $S^{z}=1$. The spatial distribution of
the spinon is visualized by seeing the variance in $\braket{S^{z}_{i}}$
between the $S^{z}=0$ and $S^{z}=1$ sectors, $\delta\braket{S^{z}_{i}}$. We
plot the profile of the spinon distribution for $\Delta=1$ in the bottom
panels of Fig.2. In both square and honeycomb lattices, it is observed that
the probability of spinon presence is maximized near the center of the system,
with minimal presence near the edges. This observation confirms that the spin
gap in the bulk limit can be correctly estimated using Eq. (5). We note that
the distribution of $\braket{S^{z}_{i}}$ for the $S^{z}=1$ sector and
$\delta\braket{S^{z}_{i}}$ are aymmetric because the rotational symmetry as
well as mirror symmetry of the system is broken.

##  IV Results

###  IV.1 staggered magnetization

We start by looking at what we found about staggered magnetization. This
measure adeptly quantifies the AFM alignment of magnetic moments throughout
the lattice, offering profound insights into the stability and resilience of
the Néel state against external perturbations. Our investigation spans two
lattice configurations, i.e., square and honeycomb structures, and encompasses
systems with spin magnitudes of $S=1/2$ and $S=1$. The primary objective of
this analysis is to delineate the degree to which staggered magnetization is
influenced by the lattice geometry and spin values.

####  IV.1.1 square lattice

Figure 3:  Extrapolated values of $m_{\mathrm{sq}}^{z}(S,\Delta)$ to the
thermodynamic limit as a function of $\Delta$ for the (a) $S=1/2$ and (b)
$S=1$ square-lattice XXZ models. Dashed lines represent the results from the
SE up to the 12th order $1/\Delta$ [25, 26, 27]. Red crosses denote the
results obtained via the CCM [28] for the $S=1/2$ case. Insets show
$m_{\mathrm{sq}}^{z}(S,\Delta)$ versus $1/\Delta^{2}$ in the large-$\Delta$
region, with solid lines indicating the SE results up to the 6th order of
$1/\Delta$.

In our preceding study [13], the magnitude of staggered magnetization for the
$S=1/2$ and $S=1$ square-lattice XXZ Heisenberg models was quantified as a
function of $\Delta$. We here revisit these results to consider the extent of
quantum fluctuations in the Néel state at $\Delta\geq 1$. They are plotted in
Fig. 3. For any $S$, the quantum fluctuations are maximized at $\Delta=1$,
correspondingly minimizing the staggered magnetization magnitude.
Specifically, for $S=1/2$, we previously estimated $m^{z}_{\rm
sq}(\frac{1}{2},1.0)=0.3065$ [13], aligning closely with extant numerical
estimates: $0.3067$ by DMRG [29], $0.3093$ via the coupled cluster method
(CCM) [28], and $0.30743$ through quantum Monte Carlo (QMC) [30]. These values
approximate $60\%$ of the classical one $m^{z}_{\rm
sq}(\frac{1}{2},\infty)=0.5$. For $S=1$, the staggered magnetization rises to
about $80\%$ of its classical counterpart $m_{z}^{\rm st}(1,\infty)=1$, as
numerically determined to be $m^{z}_{\rm sq}(1,1.0)=0.8017$ [13] via DMRG and
$m^{z}_{\rm sq}(1,1.0)=0.802$ [31] by infinite projected entangled pair states
(iPEPS). These results underscore the substantial suppression of quantum
fluctuations with an increase in $S$. In fact, an expansion in terms of $1/S$
yields $m^{z}_{\rm sq}(S,1.0)=1-0.1966019S^{-1}+0.00087S^{-3}+O(S^{-4})$ [32,
33, 34]. Since the coefficients of higher order terms than $1/S$ are very
small, a rapid convergence to the classical (or Ising) limit $m^{z}_{\rm
sq}(S,1.0)/S\to 1$ with increasing $S$ is naively expected. This trend has
received numerical validation [13].

To explore the evolution of quantum fluctuations with $\Delta$, we compare our
DMRG results to SE predictions. The SE results are plotted by dotted lines in
Fig. 3. At $\Delta=1$, our DMRG analysis yields $m^{z}_{\rm
sq}(\frac{1}{2},1.0)=0.3065$ and $m^{z}_{\rm sq}(1,1.0)=0.8017$, whereas SE up
to the 12th order offers $m^{z}_{\rm sq}(\frac{1}{2},1.0)=0.3462$ and
$m^{z}_{\rm sq}(1,1.0)=0.8579$. Despite incorporating terms up to the 12th
order, a discrepancy remains due to unaccounted quantum fluctuations, with
errors around $7.9\%$ for $S=1/2$ and $5.6\%$ for $S=1$. Increasing $\Delta$
slightly to $1.05$ reduces these deviations significantly to $0.76\%$ and
$0.44\%$ for $S=1/2$ and $S=1$, respectively, illustrating a marked decrease
in quantum fluctuations with enhanced XXZ anisotropy. In the $\Delta=\infty$
limit, SE approaches exactness, nullifying quantum fluctuations. These
discrepancies highlight the challenge of fully accounting for quantum
fluctuations, with a deviation of approximately $7.9\%$ for $S=1/2$ and
$5.6\%$ for $S=1$. However, a slight increase in $\Delta$ to $1.05$ yields
$m^{z}_{\rm sq}(\frac{1}{2},1.05)=0.3640$ (DMRG) and $m^{z}_{\rm
sq}(\frac{1}{2},1.05)=0.3678$ (SE12) for $S=1/2$, and $m^{z}_{\rm
sq}(1,1.05)=0.8767$ (DMRG) and $m^{z}_{\rm sq}(1,1.05)=0.8811$ (SE12) for
$S=1$, significantly reducing deviations to $0.76\%$ and $0.44\%$,
respectively. (Hereinafter, when considering up to the $n$th order in the SE,
we denote it as SE$n$ if needed.) This definitely suggests a pronounced
reduction in quantum fluctuations attributable to XXZ anisotropy.

In the limit of $\Delta=\infty$, the SE analysis approaches exact since
quantum fluctuations disappear. Utilizing SE up to the 6th order in
$1/\Delta$, we express the staggered magnetization for $S=1/2$ systems as
$2m^{z}_{\rm
sq}(\frac{1}{2},\Delta)=1+m_{2}/\Delta^{2}+m_{4}/\Delta^{4}+m_{6}/\Delta^{6}$.
The coefficients are calculated to be $m_{2}=-2/9=-0.222222\dots$,
$m_{4}=-8/255=-0.0355555\dots$, and $m_{6}=-0.01894258$ for $S=1/2$ [35, 36,
37, 25, 27]. In the case of $S=1$ systems, staggered magnetization is
similarly formulated as $m^{z}_{\rm
sq}(1,\Delta)=1+m_{2}/\Delta^{2}+m_{4}/\Delta^{4}+m_{6}/\Delta^{6}$, with
coefficients $m_{2}=-4/49=-0.081632653\dots$, $m_{4}=-0.026959099$, and
$m_{6}=-0.0136997515$ [26, 27]. The magnitude of the lowest-order term, i.e.,
$m_{2}$, for $S=1/2$ is approximately 2.7 times that for $S=1$, implying a
significant difference in quantum fluctuations between the two. Fitting our
data for $0\leq 1/\Delta\leq 0.05$, we obtain $m_{2}=-0.222222225$,
$m_{4}=-0.0355542736$, and $m_{6}=-0.0189663810$ for $S=1/2$, and
$m_{2}=-0.081632653$, $m_{4}=-0.026959397$, and $m_{6}=-0.013867377$ for
$S=1$. These coefficients are in almost perfect agreement with the SE results.

Additionally, in contrast to the singular behavior near $\Delta=1$ predicted
by SWT, which posits that $m^{z}_{\rm
sq}(\frac{1}{2},\Delta)=\sum_{n=0}^{\infty}m_{n}(\Delta-1)^{n/2}$ [27, 15,
36], our observations reveal that $m^{z}_{\rm sq}(\frac{1}{2},\Delta)$ is
almost linearly proportional to $\Delta-1$ within the range
$1\leq\Delta\lesssim 1.01$. This behavior aligns with results from the CCM)
[28]. However, we further ascertain that, except in the immediate vicinity of
$\Delta=1$, the staggered magnetization for the square-lattice XXZ Heisenberg
model at $\Delta\gtrsim 1$ can be approximately and quantitatively captured by
SE, provided the expansion includes up to the 12th order of $1/\Delta$.

####  IV.1.2 honeycomb lattice

Figure 4:  Magnitude of staggered magnetization for the $S=1/2$ and $S=1$
honeycomb-lattice XXZ models as a function of $\Delta$. Finite-size scaling
analyses of (a) $m^{z}_{\rm hon}(\frac{1}{2},\Delta)$ and (b) $m^{z}_{\rm
hon}(1,\Delta)$. Extrapolated values of (c) $m^{z}_{\rm
hon}(\frac{1}{2},\Delta)$ and (d) $m^{z}_{\rm hon}(1,\Delta)$ to the
thermodynamic limit as a function of $\Delta$. The dashed lines show the SE
results up to the 12th order of $1/\Delta$ [38]. Insets show $m^{z}_{\rm
hon}(S,\Delta)$ versus $1/\Delta^{2}$ in the large-$\Delta$ region. Solid
lines represent the SE results up to the 6th order of $1/\Delta$.

In exploring the honeycomb-lattice XXZ model, we note that each site is
connected to three neighboring sites, a contrast to the four neighbors in the
square-lattice model. This difference leads us to hypothesize that the Néel
state in the honeycomb-lattice model might exhibit less stability compared to
its counterpart in the square model due to the potential for increased quantum
fluctuations. To investigate this, we delve into the dependence of staggered
magnetization on the anisotropy parameter $\Delta$ for both $S=1/2$ and $S=1$
within the honeycomb lattice, drawing comparisons with our previous findings
for the square-lattice model.

In Figs. 4(a,b), we perform finite-size scaling analyses for
$m^{z}_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and
$m^{z}_{\mathrm{hon}}(1,\Delta)$ over various values of $\Delta$. Broadly,
this scaling examines the decay of Friedel oscillations with distance from a
spin pinned at the system edge. Given that the Néel state features a staggered
(commensurate) arrangement of spins in our projected 1D chain, we anticipate a
smooth decay of the Friedel oscillations with distance. This has been verified
in Sec. III.3.1. Consequently, $m^{z}_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and
$m^{z}_{\mathrm{hon}}(1,\Delta)$ extrapolate smoothly to the thermodynamic
limit as functions of inverse system size. Furthermore, the observation that
convergence with respect to size occurs more rapidly as $\Delta$ increases
reflects the diminishing quantum fluctuations.

The extrapolated values of $m^{z}_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and
$m^{z}_{\mathrm{hon}}(1,\Delta)$ are plotted as a function of $\Delta$ in
Figs. 4(c,d), respectively. In the case where the quantum fluctuations are
largest at $\Delta=1$, we obtain
$m^{z}_{\mathrm{hon}}(\frac{1}{2},1.0)=0.2764$ for $S=1/2$ and
$m^{z}_{\mathrm{hon}}(1,1.0)=0.7646$ for $S=1$. These values are in close
agreement with various numerical methods from previous studies: $0.2857$ [39],
$0.2720$ [40], $0.2611$ [41] by DMRG, $0.2677$ [42] by QMC, $0.262$ [43] by
exact diagonalization, $0.2730$ [19] by CCM for $S=1/2$, and $0.7412$ [19] by
CCM for $S=1$. These are $\sim 55\%$ and $\sim 76\%$ of their respective
classical values, only slightly smaller despite of the greater quantum
fluctuations compared to the square lattice case. On the other hand, the SE
analyses up to the 12th order of $1/\Delta$ lead to $m^{z}_{\rm
hon}(\frac{1}{2},1)=0.3409$ and $m^{z}_{\rm hon}(1,1)=0.8139$ for $S=1/2$ and
$S=1$, respectively. These values, when compared with those obtained from our
DMRG simulations, exhibit discrepancies of $6.5\%$ for $S=1/2$ and $4.9\%$ for
$S=1$. Interestingly, these deviations are rather smaller than those observed
for the square lattice, where the discrepancies are notably lower at $7.9\%$
for $S=1/2$ and $5.6\%$ for $S=1$. Furthermore, the increase in magnetization
when $S$ changes from $1/2$ to $1$ is similar to that in the square lattice.
This may be expected from the coefficients of the $1/S$ series expansion for
the honeycomb lattice $m^{z}_{\mathrm{hon}}(S,1.0)=1-0.258193/S+\cdots$ [15],
which are close to those for the square lattice.

Let us then see the $\Delta$-dependence. As illustrated in Figs. 4(c,d), the
magnetization rapidly approaches classical values with increasing $\Delta$. It
is also evident that, apart from the immediate vicinity of $\Delta=1$, the
magnetization is well captured by the SE12 predictions. Indeed, a slight
increase in $\Delta$ from $1$ to $1.05$ yields
$m^{z}_{\text{hon}}(\frac{1}{2},1.05)=0.3188$ (DMRG) versus
$m^{z}_{\text{hon}}(\frac{1}{2},1.05)=0.3478$ (SE12) for $S=1/2$, and
$m^{z}_{\text{hon}}(1,1.05)=0.8242$ (DMRG) versus
$m^{z}_{\text{hon}}(1,1.05)=0.8427$ (SE12) for $S=1$. The discrepancies
between the DMRG and SE12 values significantly decrease from $6.5\%$ and
$4.9\%$ at $\Delta=1$ to $2.9\%$ and $1.9\%$ at $\Delta=1.05$ for $S=1/2$ and
$S=1$, respectively. Nevertheless, the quantum fluctuations remain
comparatively large, and thus the reduction is not as dramatic as in the case
of the square lattice, where the discrepancies decreased from $7.9\%$ and
$5.6\%$ to $0.76\%$ and $0.44\%$.

By employing expansions up to the 4th order in $1/\Delta$, the staggered
magnetization for both spin-$1/2$ and $S=1$ systems can be expressed as
$m_{z}(S=1/2)=0.5-m_{2}/\Delta^{2}+m_{4}/\Delta^{4}+o(1/\Delta^{6})$. The
coefficients are found to be $m_{2}=-3/16=-0.1975$,
$m_{4}=31/768=0.0403645833\dots$ for $S=1/2$ and $m_{2}=-3/25=-0.12$,
$m_{4}=-17977/54000=-0.3329074\dots$ for $S=1$ [38]. Fitting our data for
$0\leq 1/\Delta\leq 0.02$, we obtain $m_{2}=-0.187499964$,
$m_{4}=0.0394523682$ for $S=1/2$ and $m_{2}=-0.120000004$,
$m_{4}=-0.0332430099$ for $S=1$. These values are in good agreement with those
from SE.

Thus, our examination of the honeycomb-lattice XXZ model reveals that, akin to
the case for square lattice, staggered magnetization can be reasonably
approximated by SE12 for $\Delta>1$.

###  IV.2 Spin gap

Next, we delve into the investigation of the spin gap, which serves as an
indicator of the stability of Néel LRO when spin rotation symmetry about the
z-axis is explicitly broken by staggered magnetization. This parameter is
essential for understanding the energy required to excite the system from its
Néel ground state, a facet less explored compared to direct magnetization
measurements. Employing methodologies analogous to those used in our
magnetization studies, we extend our analysis to both square and honeycomb
lattice configurations for $S=1/2$ and $S=1$ systems. Our aim is to elucidate
the behavior of the spin gap across varying lattice geometries and spin
magnitudes, offering insights into the intricacies associated with excitations
from the Néel state.

####  IV.2.1 square lattice

Figure 5:  Spin gap for the $S=1/2$ and $S=1$ square-lattice XXZ models as a
function of $\Delta$. Finite-size scaling analyses of (a)
$\varepsilon_{\mathrm{sq}}(\frac{1}{2},\Delta)$ and (b)
$\varepsilon_{\mathrm{sq}}(1,\Delta)$. Extrapolated values of (c)
$\varepsilon_{\mathrm{sq}}(\frac{1}{2},\Delta)$ and (d)
$\varepsilon_{\mathrm{sq}}(1,\Delta)$. to the thermodynamic limit as a
function of $\Delta$. The dashed lines show the SE results up to the 10th
order of $1/\Delta$ [27]. Insets show $\varepsilon_{\mathrm{sq}}(S,\Delta)$ in
the large $\Delta$ region.Solid lines represent the SE results up to the 6th
order of $1/\Delta$.

Let us first examine the spin gap in the case of a square lattice. In Figs.
5(a,b), finite-size scaling analyses for
$\varepsilon_{\mathrm{sq}}(\frac{1}{2},\Delta)$ and
$\varepsilon_{\mathrm{sq}}(1,\Delta)$ are shown across various values of
$\Delta$. Our scaling analysis reveals that a smooth extrapolation of the spin
gap often suggests that the scaling function resembles the contour of magnon
band near the Fermi level. For both spin magnitudes at $\Delta=1$, the spin
gap data aligns closely with a linear fit, extrapolating towards zero in the
thermodynamic limit, albeit with minor numerical uncertainties inherent to the
extrapolation process. This linear fit aligns with expectations for gapless
systems where the linear magnon band structure near the Fermi points
dominates. The actual values of extrapolated spin gap are
$\varepsilon_{\mathrm{sq}}(\frac{1}{2},1)=-0.00170006$ and
$\varepsilon_{\mathrm{sq}}(1,1)=0.00684318$. This hints at a Néel LRO that,
despite being stable, exhibits a spin gap of zero due to the arbitrary
direction of symmetry breaking. Moreover, at slightly increased values of
$\Delta$, namely $1.1$ and $1.2$, the observed small gaps corroborate the
quadratic dispersion expected near the Fermi points. As the gap widens,
quantum fluctuations wane, leading to a narrower bandwidth and, thus, a
diminished size dependence of the gap. Utilizing SBC to transform the 2D
lattice into an effective 1D system allows for the original 2D Fermi surface
to be conceptualized as a ’Fermi line,’ aiding in the smooth scaling of the
gap.

Figs. 5(c,d) display the extrapolated spin gaps
$\varepsilon_{\mathrm{sq}}(\frac{1}{2},\Delta)$ and
$\varepsilon_{\mathrm{sq}}(1,\Delta)$, showcasing a trend similar to that of
magnetization with increasing $\Delta$, rapidly approaching the classical
values $4\Delta S$. Excluding the region immediately around $\Delta=1$, the
behavior of the spin gap correlates well with SE predictions up to the 10th
order in $1/\Delta$. Specifically, a SE up to the 2nd order in $1/\Delta$
formulates $\varepsilon_{\mathrm{sq}}(S,\Delta)=2+m_{2}/\Delta^{2}+{\cal
O}(1/\Delta^{4})$, where $m_{2}$ equals $-5/3=1.6666666\cdots$ for $S=1/2$ and
$-50/21=2.3809523\cdots$ for $S=1$. Our data fitting for $0\leq 1/\Delta\leq
0.02$ yields $m_{2}=-1.66803135$ for $S=1/2$ and $m_{2}=-2.37996804$ for
$S=1$, in close agreement with these SE coefficients, demonstrating the
validity of our method.

Figure 6:  $\Delta$-normalized spin gap for the (a) $S=1/2$ and (b) $S=1$ square-lattice XXZ models as a function of $1-\Delta^{-2}$. Circles denote DMRG results in the thermodynamic limit, while the solid line represents a fit using Eq. (6).  Table 1:  Comparison of spin gaps as a function of $\Delta$ for the $S=1/2$ square-lattice XXZ Heisenberg model obtained via DMRG and CCM calculations.  $\Delta$ | $\varepsilon_{\mathrm{sq}}(\frac{1}{2},\Delta)$ |  $\Delta$ | $\varepsilon_{\mathrm{sq}}(\frac{1}{2},\Delta)$
---|---|---|---
DMRG |  CCM |  DMRG |  CCM
1.00 |  -0.0017 |  -0.0086 |  1.60 |  2.1995 |  2.2279
1.10 |  0.6810 |  0.5601 |  1.70 |  2.4578 |  2.4921
1.15 |  0.8722 |  0.7811 |  1.80 |  2.7083 |  2.7465
1.20 |  1.0461 |  0.9805 |  1.90 |  2.9535 |  2.9934
1.25 |  1.2088 |  1.1646 |  2.00 |  3.1942 |  3.2344
1.30 |  1.3627 |  1.3371 |  2.50 |  4.3496 |  4.3828
1.35 |  1.5118 |  1.5004 |  3.00 |  5.4546 |  5.4790
1.40 |  1.6560 |  1.6563 |  3.50 |  6.5306 |  6.5481
1.50 |  1.9334 |  1.9509 |  4.00 |  7.5880 |  7.6008

Investigating the behavior of spin gap near $\Delta=1$, SWT anticipates a
singularity, encapsulated by the following relation:

$\displaystyle\frac{\varepsilon}{\Delta}=\eta_{1}(1-\Delta^{-2})^{1/2}+\eta_{2}(1-\Delta^{-2})+\eta_{1}(1-\Delta^{-2})^{3/2}+\cdots,$
(6)

with constants $\eta_{1}=4S-0.431436$ and $\eta_{2}=1.2732$ derived from SWT
predictions [44]. To scrutinize this predicted behavior, we present DMRG
results for the normalized spin gap, $\frac{\varepsilon}{\Delta}$, as a
function of $1-\Delta^{-2}$ around $\Delta=1$ in Figs. 6(a,b). For both
$S=1/2$ and $S=1$, we can reasonably fit our DMRG data by Eq. (6), yielding
$\eta_{1}=1.36394987$, $\eta_{2}=0.0667245772$, $\eta_{3}=0.572388620$ for
$S=1/2$ and $\eta_{1}=3.23450360$, $\eta_{2}=0.754467400$,
$\eta_{3}=0.0226881660$ for $S=1$. These leading coefficients are close to the
SWT prediction of $\eta_{1}=1.568564$ for $S=1/2$ and $\eta_{1}=3.568564$ for
$S=1$, indicating an increase in $\eta_{1}$ with $S$, consistent with the fact
that in the large $S$ limit, the gap jumps to
$\varepsilon_{\mathrm{sq}}(S=\infty,\Delta=1^{+})=4\Delta S$ as soon as XXZ
anisotropy is introduced. Our numerical investigation thus substantiates the
SWT-predicted singularity in the spin gap near $\Delta=1$. Nonetheless, for
$S=1/2$, contrasting viewpoints emerge, such as those from CCM analyses, which
suggest a near-linear relation, $\varepsilon\propto\Delta$, diverging from the
expected singular behavior. This discrepancy could arise from our
extrapolations to the thermodynamic limit, particularly near $\Delta=1$, which
are marginally higher than those deduced from CCM. As reference, Table 1
compares the spin gap values derived from DMRG and CCM for various $\Delta$
settings in the $S=1/2$ case. A more in-depth examination will be imperative
in future studies to resolve these discrepancies and fully delineate the
characteristics of spin gap near $\Delta=1$.

####  IV.2.2 honeycomb lattice

Figure 7:  Spin gap for the $S=1/2$ and $S=1$ honeycomb-lattice XXZ models as
a function of $\Delta$. Finite-size scaling analyses of (a)
$\varepsilon_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and (b)
$\varepsilon_{\mathrm{hon}}(1,\Delta)$. Extrapolated values of (c)
$\varepsilon_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and (d)
$\varepsilon_{\mathrm{hon}}(1,\Delta)$. to the thermodynamic limit as a
function of $\Delta$. The dashed lines show the SE results up to the 10th
order of $1/\Delta$ [38]. Insets show $\varepsilon_{\mathrm{sq}}(S,\Delta)$ in
the large $\Delta$ region.Solid lines represent the SE results up to the 6th
order of $1/\Delta$.

Finally, we examine the spin gap for the honeycomb-lattice model. In Figs.
7(a,b), we conduct finite-size scaling analyses for
$\varepsilon_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and
$\varepsilon_{\mathrm{hon}}(1,\Delta)$ over various values of $\Delta$. We see
that smooth scaling is possible for all values of $\Delta$, similar to the
case of the square lattice. Consistently, at $\Delta=1$, the scaling function
is almost linear, approaching to nearly zero as $1/L$ decreases; while for
$\Delta>1$, a quadratic behavior indicative of gap opening is observed. The
actual extrapolated values are
$\varepsilon_{\mathrm{hon}}(\frac{1}{2},1)=-0.01916648$ and
$\varepsilon_{\mathrm{hon}}(1,1)=-0.01447616$, indicating slight but larger
deviations from zero than those observed in the square-lattice case. The
square lattice, having one site per structural unit cell, allowed calculations
up to $L=14$, whereas the honeycomb lattice, with two sites per structural
unit cell, limits computations to $L=9$ for a comparable computational cost,
potentially leading to relatively larger scaling errors towards the
thermodynamic limit in the honeycomb-lattice case.

The extrapolated values of $\varepsilon_{\mathrm{hon}}(\frac{1}{2},\Delta)$
and $\varepsilon_{\mathrm{sq}}(1,\Delta)$ are plotted as a function of
$\Delta$ in Figs. 7(c,d), showing a behavior broadly similar to the
magnetization versus $\Delta$ relationship. Interestingly, the saturation
towards the classical value $3\Delta S$ seems to be more gradual than observed
in the square-lattice case, possibly a reflection of heightened quantum
fluctuations within the honeycomb structure.

A notable finding is that SE analyses up to the 10th order in $1/\Delta$ are
almost inapplicable in the region $\Delta\lesssim 2$ for S=1/2. Due to the
large coefficients of higher-order terms, adding each successive term causes
the values to oscillate significantly near $\Delta=1$, making approximation by
SE very hard in the vicinity of $\Delta=1$ (see Appendix A). This emphasizes
the critical influence of quantum fluctuations on the spin excitation in its
Néel state. Conversely, for $S=1$, the SE results up to the 10th order in
$1/\Delta$ can reasonably describe the general gap behavior, except in the
immediate vicinity of $\Delta=1$.

To substantiate the accuracy of our DMRG calculations, we undertake a
comparison with SE predictions, particularly in the regime of large $\Delta$.
The SE up to the 2nd order in $1/\Delta$ is expressed as
$\varepsilon_{\mathrm{hon}}(S,\Delta)=2+m_{2}/\Delta^{2}+{\cal
O}(1/\Delta^{4})$, where the coefficients are determined to be
$m_{2}=-15/8=1.875$ for $S=1/2$ and $m_{2}=-39/20=1.95$ for $S=1$ [38].
Fitting our DMRG data within the range $0\leq 1/\Delta\leq 0.02$ yields
coefficients $m_{2}=-1.85284825$ for $S=1/2$ and $m_{2}=1.94992536$ for $S=1$.
These findings are in remarkable concordance with the established SE
coefficients, further affirming the reliability of our DMRG computational
approach.

Figure 8:  $\Delta$-normalized spin gap for the (a) $S=1/2$ and (b) $S=1$
honeycomb-lattice XXZ models as a function of $1-\Delta^{-2}$. Circles denote
DMRG results in the thermodynamic limit, while the solid line represents a fit
using Eq. (6).

In alignment with observations for the square lattice, SWT also forecasts
singular behavior near $\Delta=1$ for the honeycomb lattice. The asymptotic
form of this behavior is encapsulated by Eq. (6), with SWT providing the
coefficients as $\eta_{1}=3S-0.423239$ and $\eta_{2}=1.2405$[15]. To verify if
our DMRG data align with these predictions, we analyze the normalized spin
gap, $\frac{\varepsilon}{\Delta}$, as a function of $1-\Delta^{-2}$
approaching $\Delta=1$, as depicted in Figs. 8(a,b).

Our fits to Eq. (6) for both $S=1/2$ and $S=1$ seem to be reasonable and yield
$\eta_{1}=1.03093109$, $\eta_{2}=-0.614589918$, $\eta_{3}=1.06200153$ for
$S=1/2$, and $\eta_{1}=2.12914213$, $\eta_{2}=0.923738066$,
$\eta_{3}=-0.0550192661$ for $S=1$. These leading coefficients exhibit a
notable correspondence with SWT-anticipated $\eta_{1}=1.076761$ for $S=1/2$
and $\eta_{1}=2.576761$ for $S=1$, underscoring our numerical validation of
the predicted singularity in the spin gap as $\Delta$ approaches $1$.
Incidentally, the observation that $\eta_{1}$ for $S=1$ is larger than that
for $S=1/2$ suggests a trend towards the discontinuity where the gap leaps to
$\varepsilon_{\mathrm{sq}}(S=\infty,\Delta=1^{+})=3\Delta S$ with the
introduction of XXZ anisotropy in the limit of $S=\infty$.

##  V summary and discussion

Figure 9:  (a) Ratio of the spin gap to the magnitude of magnetization, which
is normalized to approach unity in the classical limit, as a function of
$\Delta$ for the the $S=1/2$ and $S=1$ XXZ Heisenberg models on square and
honeycomb lattices. (b) The values plotted in (a) subtracted by
$r_{\mathrm{hon}}(\frac{1}{2},\Delta)$.

We have achieved a comprehensive study of the $S=1/2$ and $S=1$ XXZ Heisenberg
model on square and honeycomb lattices. By employing the DMRG method, we
systematically analyzed the evolution of staggered magnetization and the
associated spin gap across a wide range of easy-axis anisotropies. A key to
enhancing DMRG performance was the implementation of SBC, which enabled an
exact mapping of the original 2D clusters onto a 1D chain. This technique
significantly improved our ability to perform efficient finite-size scaling
analysis, thereby facilitating the extrapolation of physical quantities to the
thermodynamic limit with high accuracy.

Given the difference in the number of adjacent sites – four for the square
lattice versus three for the honeycomb lattice – it is reasonable to
anticipate greater quantum fluctuations in the Néel phase for the honeycomb
structure. This implies a potentially lower stability of the Néel LRO in the
honeycomb lattice as compared to the square lattice. Contrary to what might be
expected from the increased quantum fluctuations, we found that the magnitude
of staggered magnetization in the honeycomb lattice is only marginally smaller
than that in the square lattice. Furthermore, across all models investigated,
the dependence of staggered magnetization on $\Delta$, except very close to
$\Delta=1$, is well captured by SE up to the 12th order. The $\Delta$
dependence of the spin gap closely mirrors that of staggered magnetization,
with most cases being approximately describable by SE up to the 10th order.
However, for the $S=1/2$ honeycomb lattice, significant deviations from the
10th order SE predictions are observed near the isotropic Heisenberg limit,
underscoring the pivotal role of quantum fluctuations on the spin gap in its
Néel state. Moreover, for all models considered, our results align numerically
with the singular behavior of the spin gap near the isotropic Heisenberg limit
as predicted by SWT.

We here delve into why the $\Delta$-dependence of the spin gap for the $S=1/2$
honeycomb lattice case significantly deviates from the SE predictions. The
spin gap in this context serves as an indicator of the stability of staggered
magnetization in the z-direction within the Néel LRO. To analyze this, we
consider the ratio of the spin gap to the magnitude of magnetization,
normalized to approach unity in the classical limit, defined for the square
and honeycomb lattices respectively as:

$\displaystyle
r_{\mathrm{sq}}(S,\Delta)=\frac{\varepsilon_{\mathrm{sq}}(S,\Delta)}{4\Delta
m^{z}_{\mathrm{sq}}(S,\Delta)},\
r_{\mathrm{hon}}(S,\Delta)=\frac{\varepsilon_{\mathrm{hon}}(S,\Delta)}{3\Delta
m^{z}_{\mathrm{hon}}(S,\Delta)}.$ (7)

We plot the DMRG results for all four models as a function of $\Delta$ in
Fig.9(a). In all cases, as $\Delta$ approaches 1, increased quantum
fluctuations reduce the stability of staggered magnetization, driving the
ratio towards zero. Notably, at the same $\Delta$,
$r_{\mathrm{hon}}(\frac{1}{2},\Delta)$ is particularly lower than the other
three cases, suggesting that the $S=1/2$ honeycomb lattice experiences larger
quantum fluctuations relative to its magnetization size, implying a relatively
unstable Néel LRO. This observation is consistent with the larger coefficients
for higher-order terms in SE. The differences between
$r_{\mathrm{hon}}(\frac{1}{2},\Delta)$ and the other three ratios are plotted
in Fig.9(b), revealing the largest quantum fluctuations around $\Delta\approx
1.5$ for the $S=1/2$ honeycomb lattice case, with the $S=1/2$ square lattice
next in line, as intuitively expected. Conversely, for $S=1$, this physical
quantity appears less dependent on the lattice geometry. The data plotted in
Fig.9 are shown in Appendix B.

Additionally, we have tested the application of a recently proposed novel
method for numerical calculations in 2D systems, demonstrating its
effectiveness. The application of SBC has proven to be particularly effective
in mapping 2D models onto 1D periodic chains, significantly simplifying the
computational complexity of our analyses and enabling the precise
determination of physical quantities in the thermodynamic limit. This
methodological innovation opens new avenues for the study of quantum phenomena
in complex lattice systems, providing a robust framework for exploring the
effects of lattice geometry and spin interactions on magnetic order and
excitations.

## Acknowledgements

We thank Ulrike Nitzsche for technical support. This work was supported by
Grants-in-Aid for Scientific Research from JSPS (Projects No. JP20H01849, No.
JP20K03769, and No. JP21J20604). M.K. acknowledges support from the JSPS
Research Fellowship for Young Scientists. M.N. acknowledges the research
fellow position of the Institute of Industrial Science, The University of
Tokyo. S.N. acknowledges support from SFB 1143 project A05 (project-id
247310070) of the Deutsche Forschungsgemeinschaft.

##  Appendix A Series expansion for spin gap in the $S=1/2$ honeycomb-lattice
XXZ Heisenberg model

Figure 10:  $\Delta$-normalized spin gap for the $S=1/2$ honeycomb lattice
case obtained via DMRG and SE, plotted as a function of $\Delta$. Here, the SE
results incorporating terms up to the $n$th order are denoted as SE$n$.

Around the classical limit, perturbative expansion up to the 10th order in
$1/\Delta$ for the $S=1/2$ honeycomb-lattice XXZ Heisenberg model yields the
spin gap as:

$\displaystyle\varepsilon_{\mathrm{hon}}(\frac{1}{2},\Delta)$
$\displaystyle=\frac{3}{2}-\frac{15}{8\Delta^{2}}+\frac{2.3046875}{\Delta^{4}}-\frac{7.05102539062}{\Delta^{6}}$
$\displaystyle+\frac{26.3766856347}{\Delta^{8}}-\frac{111.596182008}{\Delta^{10}}.$
(8)

We denote the results obtained by considering terms up to the $n$th order as
SE$n$ and plot them as a function of $\Delta$ in Fig. 10, alongside a
comparison with our DMRG results. As indicated by Eq.(8), the inclusion of
higher-order terms introduces significant oscillations near $\Delta=1$, due to
the increasing coefficients of these terms. Consequently, within the region
$\Delta\lesssim 2$, there is a noticeable deviation from the DMRG results. To
approximate the DMRG findings even closely, it is anticipated that a
considerably higher order of terms must be accounted for, suggesting that the
impact of quantum fluctuations on spin excitations is significant in this Néel
phase.

##  Appendix B DMRG data for staggered magnetization and spin gap

Table 2:  Extrapolated values of staggered magnetization, spin gap, and their ratio defined by Eq. (7) to the thermodynamic limit as a function of $\Delta$ for the $S=1/2$ square-lattice XXZ Heisenberg model obtained by our DMRG calculations.  $\Delta$ | $S=1/2$ | $S=1$
---|---|---
$\varepsilon_{\mathrm{sq}}\\!(\frac{1}{2},\Delta)$ |  $m^{z}_{\mathrm{sq}}\\!(\frac{1}{2},\Delta)$ |  $r_{\mathrm{sq}}\\!(\frac{1}{2},\Delta)$ |  $\varepsilon_{\mathrm{sq}}\\!(1,\Delta)$ |  $m^{z}_{\mathrm{sq}}\\!(1,\Delta)$ |  $r_{\mathrm{sq}}\\!(1,\Delta)$
1.00 |  -0.00170 |  0.30651 |  -0.00139 |  0.00684 |  0.80170 |  0.00213
1.05 |  0.46158 |  0.36402 |  0.30191 |  1.11193 |  0.87667 |  0.30199
1.10 |  0.68102 |  0.38422 |  0.40284 |  1.62512 |  0.89616 |  0.41214
1.15 |  0.87217 |  0.39752 |  0.47696 |  2.04938 |  0.91139 |  0.48883
1.20 |  1.04607 |  0.40789 |  0.53429 |  2.42961 |  0.92172 |  0.54915
1.30 |  1.36273 |  0.42488 |  0.61679 |  3.09541 |  0.93763 |  0.63487
1.40 |  1.65602 |  0.43681 |  0.67699 |  3.69616 |  0.94862 |  0.69578
1.50 |  1.93344 |  0.44593 |  0.72262 |  4.25626 |  0.95670 |  0.74149
1.60 |  2.19954 |  0.45312 |  0.75847 |  4.78831 |  0.96291 |  0.77699
1.70 |  2.45781 |  0.45891 |  0.78761 |  5.29980 |  0.96780 |  0.80532
1.80 |  2.70835 |  0.46366 |  0.81128 |  5.79559 |  0.97174 |  0.82836
1.90 |  2.95351 |  0.46761 |  0.83108 |  6.27901 |  0.97497 |  0.84740
2.00 |  3.19416 |  0.47094 |  0.84781 |  6.75247 |  0.97765 |  0.86335
2.50 |  4.34959 |  0.48173 |  0.90292 |  9.02033 |  0.98619 |  0.91467
3.00 |  5.45464 |  0.48742 |  0.93257 |  11.19124 |  0.99058 |  0.94148
3.50 |  6.53058 |  0.49081 |  0.95042 |  13.31055 |  0.99315 |  0.95731
4.00 |  7.58805 |  0.49298 |  0.96200 |  15.39881 |  0.99479 |  0.96747
Table 3:  Extrapolated values of staggered magnetization, spin gap, and their ratio defined by Eq. (7) to the thermodynamic limit as a function of $\Delta$ for the $S=1/2$ honeycomb-lattice XXZ Heisenberg model obtained by our DMRG calculations.  $\Delta$ | $S=1/2$ | $S=1$
---|---|---
$\varepsilon_{\mathrm{hon}}\\!(\frac{1}{2},\Delta)$ |  $m^{z}_{\mathrm{hon}}\\!(\frac{1}{2},\Delta)$ |  $r_{\mathrm{hon}}\\!(\frac{1}{2},\Delta)$ |  $\varepsilon_{\mathrm{hon}}\\!(1,\Delta)$ |  $m^{z}_{\mathrm{hon}}\\!(1,\Delta)$ |  $r_{\mathrm{hon}}\\!(1,\Delta)$
1.00 |  -0.01917 |  0.27637 |  -0.02312 |  -0.01448 |  0.75925 |  -0.00636
1.05 |  0.30264 |  0.31882 |  0.30135 |  0.77160 |  0.82632 |  0.29644
1.10 |  0.43473 |  0.34256 |  0.38456 |  1.14549 |  0.85580 |  0.40561
1.15 |  0.55835 |  0.35957 |  0.45010 |  1.46041 |  0.87559 |  0.48345
1.20 |  0.67549 |  0.37295 |  0.50312 |  1.73880 |  0.89047 |  0.54241
1.30 |  0.89595 |  0.39335 |  0.58403 |  2.24251 |  0.91197 |  0.63051
1.40 |  1.10423 |  0.40856 |  0.64351 |  2.69467 |  0.92701 |  0.69210
1.50 |  1.30397 |  0.42045 |  0.68919 |  3.11703 |  0.93819 |  0.73831
1.60 |  1.49774 |  0.43004 |  0.72558 |  3.51818 |  0.94683 |  0.77412
1.70 |  1.68690 |  0.43792 |  0.75531 |  3.90403 |  0.95368 |  0.80267
1.80 |  1.87243 |  0.44451 |  0.78007 |  4.27815 |  0.95924 |  0.82592
1.90 |  2.05486 |  0.45008 |  0.80097 |  4.64300 |  0.96381 |  0.84514
2.00 |  2.23472 |  0.45484 |  0.81887 |  5.00033 |  0.96764 |  0.86126
2.50 |  3.10562 |  0.47085 |  0.87944 |  6.71114 |  0.97988 |  0.91319
3.00 |  3.94431 |  0.47960 |  0.91379 |  8.34704 |  0.98624 |  0.94039
3.50 |  4.76235 |  0.48494 |  0.93529 |  9.94248 |  0.98997 |  0.95649
4.00 |  5.56612 |  0.48843 |  0.94967 |  11.51335 |  0.99237 |  0.96683

For additional context, we present the values of staggered magnetization and
spin gap obtained via the DMRG method at the thermodynamic limit for various
values of $\Delta$. Tables 2 and 3 detail these quantities for the square and
honeycomb lattices, respectively. Moreover, the values of the ratio of
staggered magnetization to the spin gap, as plotted in Fig. 9 of the main
text, are also provided for reference.

## References

  * Schollwöck  _et al._ [2004] U. Schollwöck, J. Richter, D. J. Farnell, and R. F. Bishop,  _Quantum Magnetism_, Lecture Notes in Physics, Vol. 645 (Springer, 2004).
  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila,  _Introduction to frustrated magnetism: materials, experiments, theory_ , Vol. 164 (Springer Science & Business Media, 2011).
  * Sachdev [2011] S. Sachdev,  _Quantum Phase Transitions_ (Cambridge University Press, 2011).
  * Heisenberg [1928] W. Heisenberg, Zur theorie des ferromagnetismus, Z. Phys. 49, 619 (1928).
  * Lines [1963] M. E. Lines, Magnetic properties of co${\mathrm{cl}}_{2}$ and ni${\mathrm{cl}}_{2}$, Phys. Rev. 131, 546 (1963).
  * Achiwa [1969] N. Achiwa, Linear antiferromagnetic chains in hexagonal abcl3-type compounds (a; cs, or rb, b; cu, ni, co, or fe), Journal of the Physical Society of Japan 27, 561 (1969), https://doi.org/10.1143/JPSJ.27.561 .
  * Savary and Balents [2016] L. Savary and L. Balents, Quantum spin liquids: a review, Reports on Progress in Physics 80, 016502 (2016).
  * Zhou  _et al._ [2017a] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Reviews of Modern Physics 89, 025003 (2017a).
  * Schrieffer and Brooks [2007] J. R. Schrieffer and J. S. Brooks,  _Handbook of high-temperature superconductivity: theory and experiment_ (Springer, 2007).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Quantum spin hall effect and topological phase transition in hgte quantum wells, science 314, 1757 (2006).
  * Zhou  _et al._ [2017b] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017b).
  * Kadosawa  _et al._ [2022] M. Kadosawa, M. Nakamura, Y. Ohta, and S. Nishimoto, One-dimensional projection of two-dimensional systems using spiral boundary conditions, arXiv:2205.15775 (2022).
  * Kadosawa  _et al._ [2023a] M. Kadosawa, M. Nakamura, Y. Ohta, and S. Nishimoto, Study of staggered magnetization in the spin-s square-lattice heisenberg model using spiral boundary conditions, Journal of the Physical Society of Japan 92, 023701 (2023a).
  * Kubo and Kishi [1988] K. Kubo and T. Kishi, Existence of long-range order in the $\mathrm{XXZ}$ model, Phys. Rev. Lett. 61, 2585 (1988).
  * Weihong  _et al._ [1991a] Z. Weihong, J. Oitmaa, and C. J. Hamer, Second-order spin-wave results for the quantum xxz and xy models with anisotropy, Phys. Rev. B 44, 11869 (1991a).
  * Viswanath  _et al._ [1994] V. S. Viswanath, S. Zhang, J. Stolze, and G. Müller, Ordering and fluctuations in the ground state of the one-dimensional and two-dimensional s=1/2 xxz antiferromagnets: A study of dynamical properties based on the recursion method, Phys. Rev. B 49, 9702 (1994).
  * Yunoki [2002] S. Yunoki, Numerical study of the spin-flop transition in anisotropic spin-$\frac{1}{2}$ antiferromagnets, Phys. Rev. B 65, 092402 (2002).
  * Braiorr-Orrs  _et al._ [2019] B. Braiorr-Orrs, M. Weyrauch, and M. V. Rakov, Numerical studies of entanglement properties in one- and two-dimensional quantum ising and xxz models, Ukrainian Journal of Physics 61, 613 (2019).
  * Bishop and Li [2016] R. Bishop and P. Li, Large-s expansions for the low-energy parameters of the honeycomb-lattice heisenberg antiferromagnet with spin quantum number s, Journal of Magnetism and Magnetic Materials 407, 348 (2016).
  * Affleck  _et al._ [1988] I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Valence bond ground states in isotropic quantum antiferromagnets, Communications in Mathematical Physics 115, 477 (1988).
  * Wojtkiewicz  _et al._ [2023] J. Wojtkiewicz, K. Wohlfeld, and A. M. Oleś, Long-range order in the xy model on the honeycomb lattice, Phys. Rev. B 107, 064409 (2023).
  * Nakamura  _et al._ [2021] M. Nakamura, S. Masuda, and S. Nishimoto, Characterization of topological insulators based on the electronic polarization with spiral boundary conditions, Phys. Rev. B 104, L121114 (2021).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * Kadosawa  _et al._ [2023b] M. Kadosawa, M. Nakamura, Y. Ohta, and S. Nishimoto, Phase diagram of the kitaev-heisenberg model using various finite-size clusters, Journal of the Physical Society of Japan 92, 055001 (2023b).
  * Singh [1989] R. R. P. Singh, Thermodynamic parameters of the $t=0$, spin-1/2 square-lattice heisenberg antiferromagnet, Phys. Rev. B 39, 9760 (1989).
  * Singh [1990] R. R. P. Singh, Quantum renormalizations in the spin-1 heisenberg antiferromagnet on the square lattice, Phys. Rev. B 41, 4873 (1990).
  * Weihong  _et al._ [1991b] Z. Weihong, J. Oitmaa, and C. J. Hamer, Square-lattice heisenberg antiferromagnet at $t=0$, Phys. Rev. B 43, 8321 (1991b).
  * Bishop  _et al._ [2017] R. Bishop, P. Li, R. Zinke, R. Darradi, J. Richter, D. Farnell, and J. Schulenburg, The spin-half xxz antiferromagnet on the square lattice revisited: A high-order coupled cluster treatment, Journal of Magnetism and Magnetic Materials 428, 178 (2017).
  * White and Chernyshev [2007] S. R. White and A. L. Chernyshev, Neél order in square and triangular lattice heisenberg models, Phys. Rev. Lett. 99, 127004 (2007).
  * Sandvik and Evertz [2010] A. W. Sandvik and H. G. Evertz, Loop updates for variational and projector quantum monte carlo simulations in the valence-bond basis, Phys. Rev. B 82, 024407 (2010).
  * Niesen and Corboz [2017] I. Niesen and P. Corboz, Emergent haldane phase in the $s=1$ bilinear-biquadratic heisenberg model on the square lattice, Phys. Rev. B 95, 180404 (2017).
  * Hamer  _et al._ [1992] C. J. Hamer, Z. Weihong, and P. Arndt, Third-order spin-wave theory for the heisenberg antiferromagnet, Phys. Rev. B 46, 6276 (1992).
  * Igarashi [1992] J.-i. Igarashi, 1/s expansion for thermodynamic quantities in a two-dimensional heisenberg antiferromagnet at zero temperature, Phys. Rev. B 46, 10763 (1992).
  * Canali and Wallin [1993] C. M. Canali and M. Wallin, Spin-spin correlation functions for the square-lattice heisenberg antiferromagnet at zero temperature, Phys. Rev. B 48, 3264 (1993).
  * Davis [1960] H. L. Davis, New method for treating the antiferromagnetic ground state, Phys. Rev. 120, 789 (1960).
  * Huse [1988] D. A. Huse, Ground-state staggered magnetization of two-dimensional quantum heisenberg antiferromagnets, Phys. Rev. B 37, 2380 (1988).
  * Parrinello and Arai [1974] M. Parrinello and T. Arai, Infinite-order cumulant expansion for spins, Phys. Rev. B 10, 265 (1974).
  * Oitmaa  _et al._ [1992] J. Oitmaa, C. J. Hamer, and Z. Weihong, Quantum magnets on the honeycomb and triangular lattices at t=0, Phys. Rev. B 45, 9834 (1992).
  * Ganesh  _et al._ [2013] R. Ganesh, J. van den Brink, and S. Nishimoto, Deconfined criticality in the frustrated heisenberg honeycomb antiferromagnet, Phys. Rev. Lett. 110, 127203 (2013).
  * Zhu  _et al._ [2013] Z. Zhu, D. A. Huse, and S. R. White, Weak plaquette valence bond order in the $s\mathbf{=}1/2$ honeycomb ${J}_{1}\mathbf{-}{J}_{2}$ heisenberg model, Phys. Rev. Lett. 110, 127205 (2013).
  * Gong  _et al._ [2013] S.-S. Gong, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Phase diagram of the spin-$\frac{1}{2}$ ${J}_{1}$-${J}_{2}$ heisenberg model on a honeycomb lattice, Phys. Rev. B 88, 165138 (2013).
  * Castro  _et al._ [2006] E. V. Castro, N. M. R. Peres, K. S. D. Beach, and A. W. Sandvik, Site dilution of quantum spins in the honeycomb lattice, Phys. Rev. B 73, 054422 (2006).
  * Albuquerque  _et al._ [2011] A. F. Albuquerque, D. Schwandt, B. Hetényi, S. Capponi, M. Mambrini, and A. M. Läuchli, Phase diagram of a frustrated quantum antiferromagnet on the honeycomb lattice: Magnetic order versus valence-bond crystal formation, Phys. Rev. B 84, 024406 (2011).
  * Hamer  _et al._ [1991] C. J. Hamer, J. Oitmaa, and Z. Weihong, Zero-temperature properties of the quantum xy model with anisotropy, Phys. Rev. B 43, 10789 (1991).


# Federated Learning for IoUT: Concepts, Applications, Challenges and
Opportunities

Nancy Victor, Rajeswari. C, Mamoun Alazab, Sweta Bhattacharya, Sindri
Magnusson, Praveen Kumar Reddy Maddikunta, Kadiyala Ramana, Thippa Reddy
Gadekallu Corresponding Author: Thippa Reddy GadekalluNancy Victor, Rajeswari.
C, Sweta Bhattacharya, Praveen Kumar Reddy Maddikunta, Thippa Reddy Gadekallu
are with the School of Information Technology, Vellore Institute of
Technology, Tamil Nadu-632014, India (Email: {nancyvictor, rajeswari.c,
sweta.b, praveenkumarreddy, thippareddy.g}@vit.ac.in). Mamoun Alazab is with
College of Engineering, IT and Environment, and is the Director of the NT
Academic Centre for Cyber Security and Innovation (ACCI) at Charles Darwin
University, Australia. Email<EMAIL_ADDRESS>Magnusson is with
Department of Computer and Systems Science (DSV), Stockholm University,
Borgarfjordsgatan, Sweden Email<EMAIL_ADDRESS>Ramana is
with Department of Information Technology, Chaitanya Bharathi Institute of
Technology, Hyderabad - 500075, Telangana, India. Email<EMAIL_ADDRESS>
###### Abstract

Internet of Underwater Things (IoUT) have gained rapid momentum over the past
decade with applications spanning from environmental monitoring and
exploration, defence applications, etc. The traditional IoUT systems use
machine learning (ML) approaches which cater the needs of reliability,
efficiency and timeliness. However, an extensive review of the various studies
conducted highlight the significance of data privacy and security in IoUT
frameworks as a predominant factor in achieving desired outcomes in mission
critical applications. Federated learning (FL) is a secured, decentralized
framework which is a recent development in machine learning, that will help in
fulfilling the challenges faced by conventional ML approaches in IoUT. This
paper presents an overview of the various applications of FL in IoUT, its
challenges, open issues and indicates direction of future research prospects.

###### Index Terms:

Federated Learning (FL), Internet of Underwater Things (IoUT), Privacy,
Underwater Sensor Networks.

##  I Introduction

Internet of Underwater Things (IoUT) is a network of smart, interconnected
underwater devices that help in monitoring the underwater environment and
thereby using the data for various underwater applications such as
environmental monitoring, disaster prediction and prevention, underwater
exploration and so on. IoUT has gained great momentum with the advances in
underwater sensor-based and underwater wireless communication-based
technologies. The underwater devices are embedded with sensors and extremely
sophisticated tracking technologies, which enable capturing of data from the
underwater environment and then reacts accordingly. The terrestrial objects
like smart phones or similar devices get connected with underwater objects,
which are further connected with virtual objects storing information
pertaining to the features, origin and sensory aspect of the underwater
objects. This information is extremely important for further analysis,
tracking, monitoring and prediction, and thus is accessed by users with the
help of Human-to-Thing and Thing-to-Thing form of communication. As an
example, the underwater autonomous vehicles collect data, sense information,
communicate with each other and then transmits information to control centers
located on land for performing various types of predictive analytics. The
information collected with the help of IoUT act as a source of conducting
tasks relevant to crash surveying, searching of ship wrecks, prediction of
tsunami, marine life health monitoring, aquatic data collection and various
forms of archaeological expeditions. Fig. 1 depicts the network architecture
of IoUT.

Figure 1: Network architecture of IoUT.

In the recent past, significant efforts have been dedicated in the development
of practical IoUT applications catering to the need of fishing resource
detection, environmental monitoring, defense and various others which all
intend to add convenience and economic advantage to life on earth. The scope
of marine exploration is ever evolving, which triggers the need of advanced
mission critical IoUT systems. These systems help to fulfil non-functional
requirements like latency, reliability, security, and timeliness. The
inability to meet these requirements lead to inaccurate decision making that
may cause disasters. Although IoUT based mission critical systems seem to have
promising prospects, they have associated challenges as well. Firstly,
information sensing often acts as a challenge in case of IoUT due to the vast
span and complexity of underwater environment which damage or interfere with
the data collection devices. This leads to inaccuracy and insufficiency in
data collection affecting subsequent data analysis and the generated results.
Secondly, transmission of information becomes extremely difficult due to
variability in underwater temperature, hydrological characteristics,
underwater pressure, varying water currents and topography interfering with
the traditional communication models. Thirdly, the data collected from the
underwater devices are often erroneous which includes large amount of
irrelevant information creating delay in data transmission and generation of
colossal waste of communication systems [1]. The traditional network and
sensor based technologies employ deterministic methods which often fail to
serve the purpose of IoUT based systems. Traditional network and sensor based
technologies are classified into communication technology, transmission
technology and energy harvesting technology. The communication technologies
include IoUT based acoustic links. Transmission technologies include various
RFID tags such as acoustic, radio and passive integrated transponder tags
(PIT), and energy harvesting technologies such as microbial fuel cell (MFC),
that play an important role in IoUT devices. The use of machine learning (ML)
in IoUT based systems have revolutionized the traditional network
methodologies fulfilling all the requirements of efficiency, fault tolerance,
reliability, low-latency, information processing and information sensing.

Role of Federated Learning in IoUT Systems:

The cloud computing based frameworks are immensely predominant in IoUT based
applications. The data collected from the various devices are aggregated and
sent to the cloud server. The required valuable information is retrieved from
the large scale data set in the cloud server by applying ML algorithms. But
they may fail to generate accurate results in case of mission critical IoUT
systems due to high transmission latency. The reason behind high latency is
due to the longer distance between the cloud server and the IoUT devices,
compelling it to be forwarded by multiple intermediate nodes. Also, the
collected data include lot of irrelevant information increasing wastage in
resource consumption and enhanced latency. There is huge possibility of data
getting compromised as well in this process impacting the security of the IoUT
frameworks. Edge computing could be implemented in this regard. The autonomous
underwater vehicles and similar floating devices which are in close proximity
to the data sources could be considered as the edge nodes which would perform
data processing at network edge eliminating useless data, reduce transmission
data size and optimize utilization of resources. In addition to maximizing
resource utilization and adding value to the data, the need of security is
non-compromising. Federated Learning (FL) is a secured decentralized
collaborative ML framework which has immense potential in developing efficient
yet secured IoUT framework. As part of the FL based approach, each node in the
network trains its subsequent sub-model independently. Each of these nodes
interact only with the encrypted sub-model parameters in association with the
fusion node to achieve an integrated global ML model [2]. The FL based IoUT
framework reduces communication overhead, and effectively mines from the large
scale data set ensuring optimal security.

Even though there are several surveys on FL and IoUT separately, there is no
survey on FL for IoUT to the best of our knowledge, which is a motivating
factor for this survey. Some of the recent surveys on IoUT and FL are
highlighted below and also summarized in Table I.

In [3], the authors have presented a comprehensive survey on the fusion of
IoUT and big marine data analytics. In another interesting study, the authors
in [4] surveyed the various networking techniques and communication
technologies for IoUT.

TABLE I: Summary of state-of-the-art implementation in IoUT.

Ref.No. |  Contributions |  Limitations
---|---|---
[3] |  A review of IoUT network architecture, IoUT communication protocols, IoUT network topologies were presented. |  This work focused mainly on IoUT but not FL and applications.
[4] |  Applications of IoT and ML for incentive designs in IoUT are studied. |  This work did not focus on the potential of FL for IoUT and the use of IoT, AI techniques for IoUT.
[5] |  This work highlighted machine learning, deep learning and blockchain techniques for FL in secure IIoT. |  This review paper was only limited to the IIoT with FL for achieving privacy and on device learning capability, applications of IIoT with FL in healthcare and automobile industry, while FL for IoUT were not studied.
[6] |  Applications of FL for aerial communications is presented. |  This work focused mainly on FL for aerial communications but not IoUT and applications.

The authors in [5] have reviewed the applications of FL for industrial
internet of things. In another interesting study, the authors in [6] presented
a survey of the applications of FL for aerial communications.

The rest of the paper is organized as follows. Section II gives an overview of
IoUT and FL techniques. Section III discusses about the various applications
of FL in IoUT and section IV explains about the various challenges, open
issues and possible future directions for integration of FL with IoUT. Section
V concludes the paper.

##  II Background

In this section, a brief introduction about IoUT, FL and the motivation behind
integration of FL with IoUT is discussed.

###  II-A Internet of Underwater Things

IoUT enable several autonomous vehicles in underwater to collaborate and
communicate with each other. The sensors used in these vehicles can collect
the data from the underwater environment that will further be transmitted to
the storage such as cloud through the Internet. The information collected from
the IoUT can be used for several purposes such as effective management of
resources in the ocean, surveying shipwrecks and crashes, early detection of
tsunamis/cyclones, monitoring the health of living species in the oceans,
surveillance of marine environment, defense, and so on. The IoUT can influence
several applications from small scientific observatory, to a medium sized
harbor to covering the ocean trade at global scale. Even though IoUT has many
similarities with IoT with respect to the functioning and the structure, due
to several factors like communication environment, computation/energy
constrained resources, there are several differences between traditional IoT
and IoUT [3].

ML/deep learning techniques can be used to uncover the hidden patterns from
the large volumes of data acquired from IoUT to realize several aforementioned
applications. Some of the AI/ML based approaches in IoUT and their associated
challenges are discussed below:

AI in Information Sensing: IoUT are usually deployed in harsh underwater
environments wherein interruption and inaccurate information sensing is a
prevalent issue. Also considering the energy consumption of these devices, it
becomes extremely difficult to charge or replace the depleted battery, leading
to inability of the systems to function properly. There are additional issues
such as erroneous antenna signals, signal interferences, malicious attacks,
and obstacles causing system failures and malfunction. The need for fault
tolerant, reliable and effective information sensing in IoUT systems is thus a
dire necessity. ML and AI approaches help to establish relationship between
the factors that contribute to device failures, which is difficult to be
computed using a deterministic mathematical model. Deep learning techniques
help in this regard by analysing large amount of historical data thereby
generating interesting inferences and predictions. It also establishes non-
linear and dynamic relationship between probability of system failures and
relevant underlying reasons in the complex underwater environment. Multi
resource data fusioning is a technique wherein diverse sensor data are
integrated which eliminates lack of accuracy and robustness of collected data
from single IoUT devices. The traditional data fusion strategies fail to
provide the desired results due to its high-dimensional and non-linear nature.
AI and machine learning approaches provide enhanced multi-sourced data fusion
especially in mission critical IoUT systems.

Challenges: The application of ML approaches in IoUT is dependent on data
which often get compromised due to security concerns. The limited energy of
the under water sensor networks (USWN) make the frameworks susceptible to
malicious attacks reducing their lifespan. The nodes can be hacked by malicous
attackers and the data collected by the nodes get read or manipulated
resulting in generation of inaccurate decisions when AI algorithms are
applied.

AI in Data Transmission: In IoUT, there exist various challenges such as
narrow availability of bandwidth, doppler effect, multi-path inferences in
underwater wireless communication hindering efficient data transmission
services. There are primarily four underwater wireless communication methods
namely - Under water Acoustic Communication (UAC), radio communication,
wireless optical communication and magnetic - inductive communication. Instead
of implementing these methods solely, the recent trend involves using multi-
modal integrative communication methods which enables highly reliable
transmission. The opto-acoustic communication method is one such approach
which includes optical and acoustic communication to ensure higher rate and
long distant data communication. The use of reinforcement learning play a
significant role in this regard providing dynamic adaptive learning capability
ensuring multi-dimensional complex channel perception and intelligent decision
making.

Challenges: Although the AI based multi-modal integrative communication
methods have promising prospects, their high bit error rates, low bandwidth
and propagation delays in the acoustic channel often invite malicious attacks.
There are possibilities of routing attacks which prevent or divert the
delivery of data packets in the network. Denial of service attacks are also
possible which stop the legitimate nodes from performing their activities.

To summarize, conventional ML algorithms face several challenges in IoUT. Some
of them are listed below:

  * •

Privacy Preservation: The data collected from several autonomous vehicles such
as submarines is very sensitive and should not be leaked into the wrong hands.
In traditional ML approach, the collected data has to be transferred to the
central cloud for further processing. This may result in exposing sensitive
information and the malicious users may take advantage and get hold of this
information while transmission or during storage.

  * •

Increased Latency: In conventional machine learning algorithms, data from all
the devices/sensors has to be aggregated in the central cloud for training the
machine learning algorithms. This incurs substantial communication cost.

  * •

Handling Big Data: In IoUT, huge volumes of data will be generated from the
devices/sensors at a rapid pace. Labeling of these data at real time and
handling the volume and velocity of the data generated is a huge challenge.

  * •

Quality of the Data: Several data may be lost during the data transmission due
to network issues. Also, external noise may be added to the data which poses a
significant challenge regarding the quality of the data that has to be trained
by the ML algorithms.

###  II-B Federated Learning

FL is a recent development of ML, where the global ML model is distributed
across the individual devices. In FL, the data need not be transferred from
individual devices to the central cloud storage, Instead, the ML algorithm
itself will be executed in the local devices, and only the model updates will
be transferred across the central storage for global ML training. As the raw
data is not sent across the central server, the privacy of the data at local
devices will not be exposed to the potential malicious users/devices.

###  II-C Motivation for Integration of Federated Learning for Internet of
Underwater Things

The underlying principle of FL, where the computation is offloaded to the
local devices, has a great potential in solving several aforementioned
challenges posed by training the data from IoUT through conventional ML
approach.

The motivation behind the integration of FL with IoUT are summarized as
follows:

  * •

Preserving privacy of sensitive data generated from IoUT devices.

  * •

Reduced latency and communication costs while transmitting the data.

  * •

FL can handle the noisy data better when compared to traditional ML
algorithms.

  * •

Reduce the possibility of data quality issues during transmission of big data
generated from the IoUT devices/sensors to the cloud.

The advantages of FL for IoUT are depicted in Fig. 2.

Figure 2: Federated Learning for IoUT.

##  III Applications

In this section several applications of FL in different IoUT verticals are
discussed.

###  III-A Environmental Monitoring

The water quality is monitored using a variety of parameters, including pH
values, oxygen levels, temperature, and dissolved metals. A variety of drones
and underwater robots are deployed to collect data and transmit it to the FL
server. FL server provides timely and appropriate solutions, assisting
engineers in making faster decisions and keeping polluted water away from the
public. Several underwater sensors are deployed to monitor and detect oil
spills at pipelines. The sensors installed outside and inside the pipelines
measure various parameters of the oil, such as pressure, speed, and
vibrations, and send the data to the sink node. The data from all sink nodes
will be aggregated and sent to the base station (monitoring center). Based on
the data gathered at the base station, a decision must be made as soon as
possible to avoid disaster. Data transmission via an underwater physical
transmission medium has a high propagation delay and error rate. To address
these issues, the authors in [7] formed clusters from monitoring stations.
Each cluster acquires local model training, and the central server
consolidates all of the local model weight sets to create a global weight set
with common features. The local cluster downloads the final global model for a
new round of training resulting in faster decisions with higher accuracy. FL
allows the devices to store and build the model in a decentralised fashion
which is not possible in ML algorithms. The experimentation was carried out on
i7-5930 processor with 12 cores at 3.50 GHz and 64 Gb of memory. In the
proposed model, more number of blocks could be transferred per second.
Also,the encoding has been decreased by 9% for 20 MB of files.

###  III-B Underwater Exploration

Underwater exploration aids in investigating the biological, chemical and
physical conditions of the underwater environment so as to utilize it for
scientific and commercial purposes. Technologies such as gears, satellites,
underwater vehicles and buoys are generally used for exploring the underwater
environment. Underwater exploration helps in discovering lost treasures and
other natural resources, tracking underwater objects and also the fish
density. Preserving data privacy, location privacy and device privacy is
essential in such scenarios. FL approaches can significantly help in
preserving the privacy of data collected by locally training the data, instead
of sharing the data to a centralized server. This is a very unique feature for
which FL is preferred over ML.

A ”federated meta learning enhanced acoustic radio cooperative framework”,
also known as ARC/FML was proposed by H. Zhao et. al. [8] to wisely use the
data collected from distributed sources such as buoy nodes. This technique
helps in sharing data across air and water. The ARC/FML technique can greatly
help in sharing underwater exploration related sensitive data. The experiment
was conducted using DeepSink and RF channel dataset. The proposed model
acheives 97% of accuracy. Similarly, an edge computing framework is proposed
using blockchain and FL mechanisms by Z. Qin et. al. [9]. This system can
efficiently deal with the security issues in the marine IoT systems. The
implementation suggests how malicious worker in FL can be simulated. During
this process, initially the experimentation was carried with attack intensity
of 20% and achieved the loss function of FL to 0.092. Later the intensity rate
was increased to 80% and achieved loss function to 0.0195. As exploring
underwater environment and gathering data regarding the same will be
beneficial for the economy of the society, huge research prospects are
awaiting in this field.

###  III-C Disaster Prevention and Mitigation

The underwater environment is always prone to disasters, both man-made and
natural. Man-made disasters usually comprise of oil spills, poisonous gas and
substance leakage and illegal dumping. One of the major man-made disasters
happened in 2010, the Deepwater Horizon incident [10] on the ”Oil Drilling Rig
Deepwater Horizon”. Concrete core was used to seal the well and nitrogen gas
was used for curing. However, a gush of natural gas leaked through the core
and it couldn’t withstand the pressure, resulting in the gas to rise to the
platform and catch fire. The explosion killed 11 workers and 94 members in the
crew were cleared from the site. The rig was sunk and resulted in a
catastrophic oil leakage from the well, and has has affected nearly 70000
square miles of ocean in the Gulf of Mexico. Natural disasters in ocean often
fall into the categories such as tsunamis, hurricanes, and tropical storms.
The Indian ocean earthquake and tsunami happened in 2004 is considered as one
of the deadliest natural disasters [11].

FL enabled IoUT solutions significantly help in disaster prevention by
monitoring the underwater environment in real-time. Data from the underwater
environment can be collected using devices such as sensors, cameras, and
underwater vehicles. These devices can be deployed at various locations, thus
enabling continuous data collection in the underwater environment. Various
sensors such as seismic pressure sensors can be used for disaster prevention.
FL allows devices to learn shared prediction model rather storing it on the
central server. Barrier systems are one of the effective ways in mitigating
the effect of disasters. FL approaches are beneficial in disaster prevention
not just because of preserving the data privacy, but also due to the improved
model performance and scalability it offers.

###  III-D Defence/ Military

The naval activities in defence include submarine detection, mine warfare,
recovery operations, and surveillance. In the conventional model, these
activities are carried out by humans in under water vessels. With advancement
in technology, underwater wireless sensors (UWSN), a branch of wireless sensor
networks, have proven to be a leading –edge communication infrastructure that
enables under water activities without human intervention. UWSN serves to
detect and classify subjects of interest in the underwater environment.

The US navy emphasize their ships, submarines, drones and on-shore
intelligence analysts to be able to share data in real time [12]. Lack of
strong data leveraging tools is a key challenge in the deployment of this
Multi Domain Operations (MDO). AI and ML based approaches can address the
overwhelming amount of data, improve data analytics for better decision
making, increase the speed of action, reduced manpower and accommodate
uncertainties also. Federated machine learning is a collaborative training
approach where training data is not exchanged to the edge device so as to
overcome constraints on training data sharing (policy, security, and coalition
constraints) and insufficient network capacity. To address the aforementioned
issues, G. Cirincione and D. Verma [13] proposed a federated machine learning
approach for deployment of MDO.

###  III-E Navigation and Localization

Exploration of underwater environment requires modern assisted navigation and
localization technologies. Some of the underwater sensors are anchored at the
ocean bottom and some are floated at different positions or freely floated in
the oceans current. The data collected from the various sensors are
potentially useful only if the sensors’ location is identified exactly. Also
the underwater sensors location act as reference points for other smart
objects, swimmers, divers and explorers. Location identification is the key
point for source detection and tracking applications. Due to the
characteristics of underwater channel, localization protocols cannot work in
underwater application. Poor underwater positioning and navigation are the two
predominant issues of IoUT. Sybil, black hole, and worm hole are significant
attacks that utilizes or modifies the localization information generated by
UWSNs. Improvements in the treatment of databases and model building could be
adapted so as to improve the privacy in localization and navigation. Federated
machine learning learn collaboratively to improve performance and overcome
coalition sharing restrictions /network constraints that restrict sharing of
training data. Thus, a continuous FL process could be adapted for localization
and navigation, as the input is streaming data and the decision cycles are in
seconds and milliseconds. This is another predominant reason why FL is
preferred over ML in IoUT environments.

The applications of FL in IoUT are depicted in Fig. 3.

Figure 3: Application verticals of IoUT.

##  IV Challenges, Open Issues and Future Directions

Integrating FL techniques in an IoUT environment proves to be advantageous due
to the various reasons discussed in the previous sections. However, various
challenges are also identified in such a setting. This section deals with the
various challenges in integrating FL with IoUT, open issues for the research
community and possible solutions that can be adapted.

Device/ Network Configuration: Deploying FL over the network edge requires the
configuration of devices or the network itself. Network configuration issues
prominently occur in underwater networks due to the fact that the nodes are
mostly mobile, thus resulting in disturbing the connectivity between nodes.
This occurs due to various factors such as path loss, noise, multipath effect
and doppler shift. Edge devices suffer from issues such as limited battery
power, storage capacity, and computational capacity. In order to enable FL on
IoUT, lightweight models with self/auto configuration methods need to be
devised. Also, a hardware and algorithm co-design can significantly help.

Data Transmission: Data transmission in underwater environment is different
when compared with the terrestrial environment. One of the major issues is
with respect to the low frequency range with which the signals need to be
transmitted, so as to avoid the same being hindered by water. Due to the long
transmission range, there is a high chance for interference and collision.
When considering the deployment of FL in IoUT, this proves to be a significant
challenge, as synchronous update is done at the server side, and this requires
the data to be sent to the server even from the slowest client, in order to
enable aggregation. However, due to the prevailing underwater environmental
characteristics and the connectivity, availability and transmission issues,
efficient solutions are required for enabling FL in IoUT. One of the promising
research directions would be to use an asynchronous update at the server side
[14].

Unreliable Channel Condition: Unreliable channel conditions can occur in
underwater environment due to various factors such as limited bandwidth,
channel noise, node mobility and transmission delays. Communication bottleneck
is one of the significant challenges in deploying FL in an IoUT environment,
that makes it difficult for the clients to share the local updates to a
centralized server. Another bottleneck is in dealing with the dynamic changes
that occur in IoUT networks since the topology itself can change over time due
to the presence of mobile nodes. Optimization techniques can be devised for
deploying FL in an IoUT environment. Also, gradient compression schemes and
aggregation frameworks that focus on computation, communication and resource
efficiency can be used for dealing with bandwidth related issues.

TABLE II: Challenges and the possible solutions for integrating FL in an IoUT
environment.

Challenge Type |  Challenge Description |  Possible Solutions and Effective Methods
---|---|---
Device and Network Configuration |  Configuring the devices and the network itself is challenging in an UIoT environment due to the factors such as limited on-device resources, storage capacity and battery power |  Light weight models, auto or self reconfiguration methods, hardware-algorithm co-design
Data Transmission |  Transmitting data in UIoT suffers from significant issues such as network connectivity, data availability, propagation speed, transmission range and rate, due to which data aggregation is impeded at the server end |  Asynchronous update for FL
Unreliable Channel Condition |  Limited network bandwidth, transmission delays, and noise are some of the factors that result in an unreliable channel condition, that will result in significant communication delays with respect to local updates |  Optimization techniques, computation, communication and resource efficient aggregation frameworks
System Heterogeneity |  The different types of devices, data formats and software systems in an UIoT environment makes it difficult for the FL technique to be integrated, since it can lead to discrepancy in the local data structure in the clients’ side |  Frameworks that support heterogeneity, global inference models
Privacy |  Data, device and location privacy need to be guaranteed to the clients in a FL-enabled UIoT environment. Model inversion attacks result in the input data being accessed by the attackers, thus questioning the client’s privacy |  Secure, light-weight protocols
Real-time generation of labels |  For generating decisions in real-time, learning mechanisms need to applied, where real-time generation of labels is required, which is quite tedious in an underwater setting |  Unsupervised learning mechanisms, Automated label generation

System Heterogeneity: The underlying hardware and software systems in an IoUT
environment differ in various aspects when compared to the terrestrial IoT
environment. This result in system design issues due to asynchronous
communication, diverse data formats and dimensions. Enabling FL approaches in
IoUT requires the handling of heterogenous devices and its capabilities, with
the help of which load can be distributed among devices based on the
availability. Frameworks that support system heterogeneity can be employed
here, that will further enable in having a global reference model.

Privacy: The sensor nodes in an IoUT environment are sparsely deployed, by
making it quite difficult to manage. Hence, this is regarded as a complex and
sensitive environment. Privacy, here deals with data, device and location
privacy. Robust solutions are essential for preserving the privacy in IoUT
networks. FL techniques help in dealing with these trust-related issues.
However, FL also has privacy limitations that include reconstruction of user
data by using the gradient information [15]. The images collected from the
IoUT setting also can be reconstructed using attacks such as model inversion.
This is indeed an open challenge to researchers to deal with the privacy
related issues in a privacy-enabled setting. One possible solution could be
the implementation of secure, light-weight protocols.

Real-time generation of labels: Generation of class labels is essential for
applying supervised learning mechanisms in a dataset. However, in IoUT,
generating labels in real-time is a time-consuming and tedious task. In order
to address this issue, it is suggested to apply unsupervised learning
mechanisms, where generation of class labels is not required. Another method
would be to use automated tools for generating labels in real-time.

Table.II summarizes the challenges and the possible solutions for deploying FL
in an IoUT environment.

##  V Conclusion

This paper presents a review of the IouT systems, its applications, challenges
and scope of future direction of research. At the outset, a brief overview of
IoUT and the underlying technologies have been discussed. The applications of
AI/ML in information sensing and data transmission is presented. However the
review of the various studies highlight the need of FL in achieving desired
characteristics in an IoUT environment. The various application areas namely
environment monitoring, defence, underwater exploration and disaster
prevention is explored emphasizing on the integration of FL with IoUT. This
reveals the potential underlying challenges of deploying FL in underwater
environment and also highlights the possible solutions, effective methods
indicating scope of future research.

## References

  * [1] X. Hou, J. Wang, Z. Fang, X. Zhang, S. Song, X. Zhang, and Y. Ren, “Machine-learning-aided mission-critical internet of underwater things,” _IEEE Network_ , vol. 35, no. 4, pp. 160–166, 2021.
  * [2] M. Alazab, S. P. RM, M. Parimala, P. K. R. Maddikunta, T. R. Gadekallu, and Q.-V. Pham, “Federated learning for cybersecurity: Concepts, challenges, and future directions,” _IEEE Transactions on Industrial Informatics_ , vol. 18, no. 5, pp. 3501–3509, 2021.
  * [3] M. Jahanbakht, W. Xiang, L. Hanzo, and M. R. Azghadi, “Internet of underwater things and big marine data analytics—a comprehensive survey,” _IEEE Communications Surveys & Tutorials_, vol. 23, no. 2, pp. 904–956, 2021.
  * [4] M. Jouhari, K. Ibrahimi, H. Tembine, and J. Ben-Othman, “Underwater wireless sensor networks: A survey on enabling technologies, localization protocols, and internet of underwater things,” _IEEE Access_ , vol. 7, pp. 96 879–96 899, 2019.
  * [5] P. Boopalan, S. P. Ramu, Q.-V. Pham, K. Dev, P. K. R. Maddikunta, T. R. Gadekallu, T. Huynh-The _et al._ , “Fusion of federated learning and industrial internet of things: A survey,” _Computer Networks_ , p. 109048, 2022.
  * [6] Q. Pham, M. Zeng, Z. Han, W. Hwang _et al._ , “Aerial access networks for federated learning: Applications and challenges,” _IEEE Netw._ , 2022.
  * [7] Y. Gao, L. Liu, B. Hu, T. Lei, and H. Ma, “Federated region-learning for environment sensing in edge computing system,” _IEEE Transactions on Network Science and Engineering_ , vol. 7, no. 4, pp. 2192–2204, 2020.
  * [8] H. Zhao, F. Ji, Q. Guan, Q. Li, S. Wang, H. Dong, and M. Wen, “Federated meta learning enhanced acoustic radio cooperative framework for ocean of things underwater acoustic communications,” _arXiv preprint arXiv:2105.13296_ , 2021\.
  * [9] Z. Qin, J. Ye, J. Meng, B. Lu, and L. Wang, “Privacy-preserving blockchain-based federated learning for marine internet of things,” _IEEE Transactions on Computational Social Systems_ , 2021.
  * [10] C. L. Hagerty, _Deepwater Horizon oil spill: Selected issues for Congress_. Diane Publishing, 2010.
  * [11] J. Telford and J. Cosgrave, “The international humanitarian system and the 2004 indian ocean earthquake and tsunamis,” _Disasters_ , vol. 31, no. 1, pp. 1–28, 2007.
  * [12] M. Karagoz, “Maritime security operations and the combined joint operations from the sea center of excellence,” in _Maritime Security and Defence Against Terrorism_. IOS Press, 2012, pp. 87–91.
  * [13] G. Cirincione and D. Verma, “Federated machine learning for multi-domain operations at the tactical edge,” in _Artificial Intelligence and Machine Learning for Multi-Domain Operations Applications_ , vol. 11006. International Society for Optics and Photonics, 2019, p. 1100606.
  * [14] M. R. Sprague, A. Jalalirad, M. Scavuzzo, C. Capota, M. Neun, L. Do, and M. Kopp, “Asynchronous federated learning for geospatial applications,” in _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer, 2018, pp. 21–28.
  * [15] J. Geiping, H. Bauermeister, H. Dröge, and M. Moeller, “Inverting gradients-how easy is it to break privacy in federated learning?” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 16 937–16 947, 2020.