#  Reconciling the theoretical and experimental electronic structure of NbO2

Samuel Berman<EMAIL_ADDRESS>School of Physics and Centre for Research on
Adaptive Nanostructures and Nanodevices (CRANN), Trinity College Dublin, The
University of Dublin, Dublin 2, D02 PN40 Ireland  Ainur Zhussupbekova  School
of Physics and Centre for Research on Adaptive Nanostructures and Nanodevices
(CRANN), Trinity College Dublin, The University of Dublin, Dublin 2, D02 PN40
Ireland  School of Chemistry, Trinity College Dublin, The University of
Dublin, Dublin 2, D02 PN40 Ireland  Jos E. Boschker  Leibniz-Institut für
Kristallzüchtung, Max-Born-Str. 2, 12489 Berlin, Germany  Jutta Schwarzkopf
Leibniz-Institut für Kristallzüchtung, Max-Born-Str. 2, 12489 Berlin, Germany
David D. O’Regan  School of Physics and Centre for Research on Adaptive
Nanostructures and Nanodevices (CRANN), Trinity College Dublin, The University
of Dublin, Dublin 2, D02 PN40 Ireland  Igor V. Shvets  School of Physics and
Centre for Research on Adaptive Nanostructures and Nanodevices (CRANN),
Trinity College Dublin, The University of Dublin, Dublin 2, D02 PN40 Ireland
Kuanysh Zhussupbekov<EMAIL_ADDRESS>School of Physics and Centre for
Research on Adaptive Nanostructures and Nanodevices (CRANN), Trinity College
Dublin, The University of Dublin, Dublin 2, D02 PN40 Ireland  School of
Chemistry, Trinity College Dublin, The University of Dublin, Dublin 2, D02
PN40 Ireland  Kh. Dosmukhamedov Atyrau University, Studenchesky Ave.,1,
Atyrau, Kazakhstan

###### Abstract

Metal-insulator transition materials such as NbO2 have generated much
excitement in recent years for their potential applications in computing and
sensing. NbO2 has generated considerable debate over the nature of the phase
transition, and the values for the band gap/band widths in the insulating
phase. We present a combined theoretical and experimental study of the band
gap and electronic structure of the insulating phase of NbO2. We carry out ab-
initio density functional theory plus $U$ calculations, directly determining
$U$ and $J$ parameters for both the Nb 4d and O 2p subspaces through the
recently introduced minimum-tracking linear response method. We find a
fundamental bulk band gap of 0.80 eV for the full DFT+$U$+$J$ theory. We also
perform calculations and measurements for a (100) oriented thin film. Scanning
tunnelling spectroscopy measurements show that the surface band gap varies
from 0.75 eV to 1.35 eV due to an excess of oxygen in and near the surface
region of the film. Slab calculations indicate metallicity localised at the
surface region caused by an energy level shift consistent with a reduction in
Coulomb repulsion. We demonstrate that this effect in combination with the
simple, low cost DFT+$U$+$J$ method can account for the band widths and p-d
gap observed in X-ray photoelectron spectroscopy experiments. Overall, our
results indicate the possible presence of a 2D anisotropic metallic layer at
the (100) surface of NbO2.

##  I I. Introduction

Materials that undergo a metal to insulator transition (MIT) attract
considerable attention due to their potential applications in areas such as
memory devices Kim  _et al._ (2013); Cha  _et al._ (2016), sensing Strelcov
_et al._ (2009); Chen  _et al._ (2019), and neural-inspired computing Kumar
_et al._ (2017); Pickett  _et al._ (2013). Niobium dioxide (NbO2) undergoes
such a transition at 1080 K, with a change in crystal structure from the
insulating distorted rutile, body centred tetragonal (BCT) phase to the
metallic rutile phase. Reconciling theoretical calculations with the
experimental measurements for this material has proven challenging, most
notably in regard to the fundamental (indirect) band gap of the insulating BCT
phase. Theoretical calculations of the band gap range from 0.058 eV Wahila
_et al._ (2019) to 1.48 eV O’Hara  _et al._ (2014). Experimental measurements
also show a broad range, from 0.7 eV O’Hara  _et al._ (2014) to 1.23 eV Sakai
_et al._ (1985). Similar discrepancies between theory and experiment exist for
other spectral properties, such as the widths of the Nb 4dxy and O 2p bands,
as well as the p-d gap. The present work seeks to resolve this discrepancy.

Experimentally, the large variation in band gap can be due to many factors,
such as crystalline quality and the technique used to measure the gap. The
smallest gaps are reported for room temperature optical measurements, usually
around 0.7 eV O’Hara  _et al._ (2014). In this case the temperature will
certainly effect the measured gap, as well as exciton binding effects. The
highest reported gap was measured on a polycrystalline sample, where a value
for the gap of 1.23 eV was extracted from resistance versus temperature
measurements Sakai  _et al._ (1985). The sample used in this measurement was
ceramic NbO2 obtained by reducing an Nb2O5 sample under a hydrogen atmosphere.
It is therefore highly likely that residual Nb2O5 in the sample would effect
the measured gap. Combined photoemission and inverse photoemission
spectroscopy (PES and IPES) experiments determined a lower bound on the band
gap to be 1.0 eV Posadas  _et al._ (2014). However, core level XPS
measurements also showed the presence of two oxidation states of Nb on the
surface, indicating a possible capping layer. Due to the extreme surface
sensitivity of IPES this could significantly impact the measured band gap. The
small thickness of the film used (4 nm) could also give rise to quantum size
effects, distorting the true gap of the bulk phase. For these reasons, it is
uncertain whether this PES/IPES measurement can be taken as a true
representation of the bulk NbO2 band structure. Recently, Stoever et al.
Stoever  _et al._ (2020) finally obtained strong agreement between the optical
gap measured at low temperature (0.85 eV) and the electronic gap determined
from resistance versus temperature measurements (0.88 eV), for high quality
bulk-like NbO2 thin films (thickness of 100 nm). They also obtained a room
temperature optical band gap of 0.76 eV, in line with other reports O’Hara
_et al._ (2014). Therefore, it would seem that experimentally there is still
some degree of uncertainty, but the most reliable measurements for the bulk
band gap indicate a value of around 0.85 to 0.88 eV.

With regard to theoretical calculations, the importance of electronic exchange
and correlation in NbO2 has become a point of some contention in the
literature. For the rutile phase, Brito et al. Brito  _et al._ (2017)
performed cluster dynamical mean field theory (cDMFT) calculations (with
assumed Nb 4d subspace parameters of $U=6$ eV, $J=1$ eV) and determined that
correlation plays a significant role in that phase. For the BCT phase, O’Hara
et al. O’Hara  _et al._ (2014) performed DFT+$U$ calculations varying $U$ from
0 to 5 eV (only on the Nb 4d subspace). They found a band gap of 0.83 eV with
$U=2$ eV, in line with the lower end of the reported experimental gaps,
however the band widths and p-d gap showed poor agreement with experiment.
With an emphasis more on exchange, calculations with the hybrid
Heyd–Scuseria–Ernzerhof (HSE) functional from both O’Hara et al. O’Hara  _et
al._ (2014) and Posades et al. Posadas  _et al._ (2014) show a band gap 1.48
eV, much higher than any reported experimental value. Other authors suggest
that correlations play at most a minor role in NbO2 Eyert (2002); Kulmus  _et
al._ (2021). Kulmus et al. Kulmus  _et al._ (2021) carried out GW plus Bethe-
Salpeter (GW-BSE) calculations and found a band gap of 0.98 eV, a slight
overestimation on the likely experimental value, yet much improved compared to
LDA or HSE. Additionally they found excellent agreement with experiment for
their calculated dielectric function. However, in addition to the band gap,
the band widths and p-d gap remain problematic, being underestimated and
overestimated, respectively.

In this paper we show that first principles DFT+$U$+$J$ calculations show
solid agreement with experiment in terms of the bulk band gap, and through
scanning tunnelling spectroscopy (STS) we reveal that excess of oxygen can
significantly influence the local surface band gap. Additionally we
demonstrate that slab calculations including DFT$+U$+$J$ (which are feasible
due to the relatively low computational expense of this method) can predict
the correct band widths and p-d gap observed in X-ray photoelectron
spectroscopy (XPS). Overall, this indicates that a relatively simple treatment
of exchange and correlation is sufficient for NbO2, when relevant surface
effects are included.

##  II II. Experimental and Computational details

###  II.1 A. Experimental Details

Experimental measurements were performed across several ultra-high vacuum
(UHV) systems. The NbO2 film was deposited on top of a rutile (110) TiO2
substrate via pulsed laser deposition (PLD) of a NbO2 target, the purity of
which was 99.9985$\%$. An aperture was employed to the central part of a KrF
laser (wavelength 248 nm), which was resulting in a power of 65 mJ/pulse and a
fluence of 2.1 cm2. For the synthesis, a pulse repetition rate of 5 Hz and in
total 12000 pulses was used. The distance between the target and substrate was
60 mm and the substrate temperature was 650∘C. A partial pressure of 0.1 mbar
of Ar gas was utilised during the synthesis.

To provide compositional and structural information of the thin NbO2 film,
XPS, X-ray diffraction (XRD), reflection high-energy electron diffraction
(RHEED) and low-energy electron diffraction (LEED) were carried out. RHEED was
conducted in situ during and after the film growth. XRD was utilised to
analyse the structure and orientation of the NbO2 crystal lattice (XRD data is
shown in the SI Fig. 1(a,b) SI (1)). The thickness of the NbO2 film was
estimated at approximately 80 nm by the evaluation of the thickness fringes
around the film peak in the $\theta$/2$\theta$ XRD scans. For the XRD
measurements, a Bruker D8 HRXRD was used with a Cu anode X-ray tube with a
solid-state 1D detector, a Göbel mirror for converting the diverging output of
the tube into a parallel beam, and a two bounce Ge monochromator for the
selection the Cu-Kα1 emission with wavelength of 1.5406 Å from the spectrum.
An Omicron MultiProbe XPS system with a monochromatized Al Kα source (XM 1000,
1486.7 eV) was used for compositional analysis and resolving a valence band
(VB) structure of NbO2 film. The base pressure of XPS is 2$\times$10-10 mbar,
and the instrumental resolution is 0.6 eV. Peak positions and elemental
components were selected from the elemental library Fairley  _et al._ (2021)
with Gaussian/Lorentzian line shapes and Shirley background. LEED measurements
were performed, in conjunction with scanning tunneling microscopy and
spectroscopy (STM/S), to analyse the surface structure of the NbO2 film grown
on top of rutile TiO2(110) substrate.

The scanning tunnelling microscope (STM) utilised in this study is a
commercial low-temperature system from Createc with a base pressure of
$5\times 10^{-11}$ mbar. All STM images were obtained at liquid nitrogen
temperature (77 K) in constant-current mode (CCM). For STS measurements, the
tip was moved across the grid with the tip height stabilised by the CCM
scanning parameters. Before the spectra was obtained at any point, the
feedback loop was turned off such that the current was allowed to be recorded
as a function of the voltage, while maintaining a constant tip height. The
feedback was subsequently turned back on and the tip was moved to the next
grid position Zhussupbekov  _et al._ (2020a, 2021a); Berman  _et al._ (2023).
The preparation chamber of the UHV system is fitted with a cooling/heating
stage, LEED and an ion gun for sputtering. The STM tips used were of
[001]-oriented single-crystalline tungsten, which electrochemically etched in
NaOH. The bias was applied to the sample with respect to the tip.

The NbO2 film was transferred under ambient conditions from the PLD camber to
the LEED and STM/S chambers. For surface preparation, the sample was subjected
to Ar+ sputtering and annealing cycles, with a voltage of 750 V, emission
current of 11 $\mu$A, chamber pressure of 10-5 mbar, for 10 minutes, followed
by annealing at 700 ∘C for 1 hour. The maximum temperature, which was measured
by the system’s K-type thermocouple, was up to 700 ∘C for 2 h in UHV. This
procedure was employed in order to remove contamination and any Nb2O5 layer
from the surface. A vacuum suitcase with a base pressure of 2$\times$10-10
mbar was employed for transferring sample between STM/LEED and XPS chambers
under UHV environment Zhussupbekov  _et al._ (2021b).

###  II.2 B. Computational Details

The low temperature, insulating phase of NbO2 takes on a BCT (distorted
rutile) crystal structure. The primitive cell for this structure contains 48
atoms and is shown in Fig. 1(a), alongside the 96 atom standard conventional
unit cell in Fig. 1(b). The crystallography of this material is most usually
described in terms of the conventional cell, which is the cell we utilise in
this study. It is worth noting that when discussing crystallographic
directions some papers use the distorted rutile (DR) cell as a reference, and
some use the almost equivalent rutile (R) cell. For clarity, the relations
between some important directions are as follows: [100]${}_{\textrm{DR}}$
$\simeq$ [110]${}_{\textrm{R}}$, [110]${}_{\textrm{DR}}$ $\simeq$
[100]${}_{\textrm{R}}$. Crystal structures were visualised using the VESTA
package Momma and Izumi (2011).

Norm-conserving pseudopotentials generated with the OPIUM OPI  code were used
throughout the study. DFT-LDA calculations were first carried out in the PWscf
code of the Quantum Espresso suite Giannozzi  _et al._ (2009, 2017) to
determine the cutoff energy and k-point sampling necessary to converge the
total energy to within $<$1 meV per atom, and the optimised ionic
geometry/unit cell parameters. A cutoff energy of 60 Ry/815 eV and k-point
sampling of 2$\times$2$\times$4 for the standard conventional cell were deemed
necessary. These parameters were then used to inform the parameters chosen for
calculations in the ONETEP Skylaris  _et al._ (2005); Prentice  _et al._
(2020) linear-scaling DFT code, which includes a DFT+$U$+$J$ implementation
O’Regan  _et al._ (2012). A 2$\times$2$\times$4 supercell consisting of 1536
atoms was used with psinc spacing equal to a/($3\cdot 5\cdot 7$)=0.4839 bohr
along the distorted rutile a axes and c/(7$\cdot$13)=0.4934 bohr along the
distorted rutile c axis. A nonorthogonal generalised Wannier function (NGWF)
cutoff radius of 10 bohr was utilised for both niobium and oxygen. For
niobium, there were 10 NGWFs and 13 valence electrons per atom, compared to 4
NGWFs and 6 valance electrons per atom for oxygen.

When using a supercell that is not simply a direct scaling of the primitive
cell (i.e., when the transformation matrix $\mathbf{S}$ is not diagonal), it
is important to verify that the high symmetry k-points in the Brillouin zone
are sampled within the supercell calculations. This is especially important
for determining the band gap in an indirect gap material like BCT NbO2.
Following the method in Lloyd-Williams and Monserrat Lloyd-Williams and
Monserrat (2015), in order for a k point $\mathbf{q}$ to be commensurate with
the supercell, the vector $\mathbf{Q}$ resulting from the transformation
matrix $\mathbf{S}$ acting on $\mathbf{q}$ must contain only integers (Eq.
1.). The necessary equations for determining if a high symmetry point is
sampled are

$\mathbf{S}\mathbf{q}=\mathbf{Q}$ (1)

$\mathbf{S}\begin{pmatrix}0\\\ {}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{2}\\\
0\end{pmatrix}=\begin{pmatrix}0\\\ 1\\\
2\end{pmatrix}\quad\mathbf{S}\begin{pmatrix}{}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{2}\\\
{}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{2}\\\
-{}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{2}\end{pmatrix}=\begin{pmatrix}2\\\ 0\\\
0\end{pmatrix}\quad\mathbf{S}\begin{pmatrix}{}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{4}\\\
{}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{4}\\\
-{}^{1}{\mskip-5.0mu/\mskip-3.0mu}_{4}\end{pmatrix}=\begin{pmatrix}0\\\ 0\\\
2\end{pmatrix}$

$\mathbf{q}=$ N $\mathbf{q}=$ Z $\mathbf{q}=$ P

In this case, we can see that the high symmetry points are correctly sampled,
in particular the N point where the valence band maximum is located.

In order to better account for the exchange and correlation effects present in
this system than is provided for in the underlying LDA functional, we utilise
the simplified, rotationally invariant DFT+$U$+$J$ method of the form
presented in Ref. Himmetoglu  _et al._ (2011). Functionals of this type are
widely used for calculations on metal oxides to account for the self
interaction error (SIE) present in approximate exchange correlation
functionals.  The Hund’s $J$ term, in particular, is associated with static
correlation error, and the DFT+$U$+$J$ functional used here was designed to
improve the description of spin-flip interactions. It has recently been
tested, for example in Refs. Lim  _et al._ (2016); Linscott  _et al._ (2018);
Orhan and O’Regan (2020); Winczewski  _et al._ (2022), and some of its
limitations have been explored in Ref. Burgess  _et al._ (2023).

Table 1:  Calculated Hubbard $U$ and Hund’s $J$ parameters for the Nb 4d and O 2p subspaces in units of eV.  | Niobium | Oxygen
---|---|---
$U$ | 2.48 | 9.02
$J$ | 0.23 | 0.90
$U_{\mathrm{eff}}=U-J$ | 2.25 | 8.12
$U_{\mathrm{full}}=U-2J$ | 2.02 | 7.22
Figure 1:  (a-b) Primitive and conventional units cells for insulating NbO2
(Nb atoms shown in green, O in red), (c-d) Linear response plots for
calculation of Hubbard $U$ and Hund’s $J$. (e-f) Band structure calculated by
spectral function unfolding for LDA and DFT+$U$+$J$ respectively, using the
implementation described in Refs. Constantinescu and Hine (2015); Popescu and
Zunger (2012) and recently demonstrated in Refs. Graham  _et al._ (2020); Loh
_et al._ (2021).

Ab-initio calculations of the LDA-appropriate Hubbard $U$ and Hund’s $J$
parameters were carried out within the minimum-tracking linear-response
approach described in Linscott et al. Linscott  _et al._ (2018), and recently
applied successfully to rutile and anatase TiO2 Orhan and O’Regan (2020). For
a closed shell system such as that of NbO2, the ‘scaled 2$\times$2’ approach
reduces to the ‘simple 2$\times$2’ approach as described in the aforementioned
papers. Within this method, $U$ and $J$ are determined by the matrix elements
$f^{\sigma\sigma^{\prime}}$ as in Eq. 2, with these matrix elements being
determined by Eq. 3. The necessary equations are

$U=\frac{1}{2}(f^{\uparrow\uparrow}+f^{\uparrow\downarrow})\qquad
J=-\frac{1}{2}(f^{\uparrow\uparrow}-f^{\uparrow\downarrow})$ (2)
$f^{\sigma\sigma^{\prime}}=\left[\left(\frac{\delta v_{\mathrm{KS}}}{\delta
v_{\mathrm{ext}}}-\mathds{1}\right)\left(\frac{\delta n}{\delta
v_{\mathrm{ext}}}\right)^{-1}\right]^{\sigma\sigma^{\prime}}$ (3)

Based on this procedure we obtain $U$ and $J$ values for niobium and oxygen,
as summarised in Table 1. The relatively large first-principles oxygen $U$
value and low metal $U$ value is consistent with literature findings,
generated using other DFT codes, on diverse closed-shell early transition-
metal oxides Kirchner-Hall  _et al._ (2021); Lambert and O’Regan (2023); Moore
_et al._ (2022). This trend reflects, among other factors, that the O 2p
orbitals are at least as localized, and sometimes considerably more so, than
the metal d orbitals in these systems (see Fig. 3 of Ref. Lambert and O’Regan
(2023)). Interestingly, and similarly to the findings on TiO2 in Ref. Orhan
and O’Regan (2020), the computed $J$ values for both species with LDA are
rather low in this system, at approximately one-tenth of the corresponding
$U$.

##  III III. Results

###  III.1 A. Band Gap

For the conventional 96 atom unit cell our LDA calculations give lattice
parameters of a=13.4445 Å (2.15% lower than experimental value 13.7020 Å) and
c=5.9394 Å (0.89% lower than the experimental value of 5.9850 Å) Bolzan  _et
al._ (1994); Dhamdhere  _et al._ (2016). The most crucial aspect of the NbO2
crystal structure is the Nb-Nb dimerization along the c-axis ([001]
direction). The lattice parameter we obtain along this direction is very close
to the experimental value, which is vital since the electronic structure of
NbO2 is known to be highly sensitive to the dimerization Páez Fajardo  _et
al._ (2021); O’Hara and Demkov (2015). The short and long Nb-Nb bond lengths
at room temperature are 2.71 Å and 3.30 Å respectively Bolzan  _et al._
(1994); Pynn  _et al._ (1976); Pynn and Axe (1976) compared to the values
obtained in our calculation of 2.66 Å and 3.31 Å, again close to the room
temperature values. Of course it should be noted that the bond length varies
as a function of temperature (roughly 0.1 Å between 100 ∘C and 800 ∘C Páez
Fajardo  _et al._ (2021); Wahila  _et al._ (2019)), but no low temperature
measurements of the bond length exist in the literature to the authors’
knowledge.

Our LDA calculated band gap of 0.39 eV is slightly higher than other LDA
calculations O’Hara  _et al._ (2014) (possibly due to details of the
pseudopotential construction), but is still underestimating the true gap by at
least a factor of 2. A similar trend emerges for the other spectral features,
with slight increase over previous LDA calculations yet still not meaningfully
approaching the experimental values. The LDA calculated band structure is
shown in Fig. 1(e), with the valance band maximum (VBM) at the N point, and
the conduction band minimum (CBM) at the $\Gamma$ point in line with previous
reports O’Hara  _et al._ (2014); Kulmus  _et al._ (2021).

Table 2:  Summary of experimental and theoretical values for the band gap. Theoretical values pertain to that of bulk NbO2. Method | Band Gap [eV]
---|---
LDA O’Hara  _et al._ (2014) | 0.35
SCAN Páez Fajardo  _et al._ (2021) | 0.48
HSE O’Hara  _et al._ (2014) | 1.48
DFT+cDMFT Brito  _et al._ (2017) | 0.73
LDA Wahila  _et al._ (2019) | 0.058
PBE Wahila  _et al._ (2019) | 0.127
mBJ Wahila  _et al._ (2019) | 0.844
GW-BSE Kulmus  _et al._ (2021) | 0.98
|
LDA (this work) | 0.39
DFT+$U^{d}$ (this work) | 0.83
DFT+$U^{d}_{\mathrm{eff}}$ (this work) | 0.78
DFT+$U^{d}$+$J^{d}$ (this work) | 0.74
DFT+$U^{d,p}$ (this work) | 0.92
DFT+$U^{d,p}_{\mathrm{eff}}$ (this work) | 0.86
DFT+$U^{d,p}$+$J^{d,p}$ (this work) | 0.80
|
Resistance vs Temperature Sakai  _et al._ (1985) | 1.23
Resistance vs Temperature Stoever  _et al._ (2020) | 0.88
Scanning Tunnelling Spectroscopy | 1.05 $\pm$ 0.29
(this work) |
Inverse Photoemission Posadas  _et al._ (2014) (4nm film) | 1.0 (lower bound)
Optical Absorbtion Stoever  _et al._ (2020) (T=4K) | 0.85
Film Absorption Edge Lee and Durand (1984) | 0.88
Optical Ellipsometry O’Hara  _et al._ (2014) (indirect gap) | 0.7

The linear response plots used to calculate the $U$ and $J$ for the Nb and O
are shown in Fig. 1(c-d), and reflect good linear response. There are a number
of different ways of implementing the Hubbard $U$ and Hund’s $J$ corrections.
In many published works, the correction is applied only the transition metal d
orbitals (such as in Ref. O’Hara  _et al._ (2014)). We consider three
possibilities for the specific form of the correction. The first option is to
only add the $+U$, not including explicit opposite-spin interaction terms in
the corrective functional (effectively treating $J$ as being zero).
Alternatively, the Dudarev functional Dudarev  _et al._ (1998) additionally
accounts for exchange-related corrections arising from interactions between
like spins, by taking $U_{\textrm{eff}}=U-J$. Finally, corrections arising
from interactions between both like and unlike spins can be obtained by
carrying out full spin polarised calculations with DFT+$U$+$J$, or
equivalently for non spin-polarised systems, as presented in Ref. Orhan and
O’Regan (2020) by using the Dudarev functional with $U_{\textrm{full}}=U-2J$
and $\alpha=J/2$. In all three of these cases, we can then also apply the same
form of the correction to the O 2p subspace, resulting in a total of 6
possible implementations of these $U$ and $J$ values.

Table 2 shows the calculated band gaps for the 6 different corrections tested,
along with experimental and theoretical values from the literature.
Interestingly we can see that even though the band gap exists between Nb 4d
states, adding the correction to the O 2p subspace does have a small effect on
the calculated gap. The gap increased by about 0.06-0.09 eV or roughly 10%
across the three different forms of the functional. The
DFT+$U^{d,p}$+$J^{d,p}$ calculated band structure is shown in Fig. 1(f). We
see that the VBM remains at the N point, however the CBM is shifted away from
the $\Gamma$ point along the $\Gamma$-X direction. This is in line with Ref.
O’Hara  _et al._ (2014), where they observed a similar change in the band
curvature resulting from applying a Hubbard $U$ correction.

In the supplemental information SI (1) (see also references Huang  _et al._
(2011); Wong  _et al._ (2014) therein) we detail the structural and chemical
characterisation of the thin film. Crucially, we see no evidence of
contamination from Nb2O5 from core level XPS (Fig. SI1), as has been present
in some previous studies Páez Fajardo  _et al._ (2021); Wahila  _et al._
(2019).

Figure 2:  (a) Constant current STM image of the (100) oriented film ($V$ =
-2.5 V and $I$ = 70 pA), (b) Band gap map from same area consisting of 35
$\times$ 35=1225 $I(V)$ curves, (c) statistical analysis of the 1225 measured
band gaps (red line represents mean of (1.05 $\pm$ 0.29) eV, blue line
represents electrical band gap measurement from ref. Stoever  _et al._ (2020)
of 0.88 eV), black line represents trimodal gaussian fit with distributions
centred around (0.53, 1.03, 1.37) eV (individual distributions shown in light
blue), inset shows the average $\mathrm{d}I/\mathrm{d}V$ curve.

STM of this film shows a relatively smooth surface with an RMS roughness of
$\simeq$ 2 Å over 10 nm2 (Fig. 2(a)), with a granular structure, often
observed for oxide thin films Kurnosikov  _et al._ (2001); Stavale  _et al._
(2012). Due to the high reactivity of Nb with oxygen, preparing a smooth
enough surface for obtaining atomic resolution is likely not possible with a
thin film (as compared with the less reactive Ti, for which atomic resolution
STM images can be obtained on thin films Ashworth and Thornton (2001); Shi
_et al._ (2020). Much more aggressive sputter anneal cycles would be
necessary, which could only be performed on a single crystal. Grid
spectroscopy measurements were performed on the same area as the STM image in
Fig. 2(a). Fig. 2(b) shows the resulting band gap map from this grid
measurement. 35 $\times$ 35 individual $I(V)$ measurements were performed on
the grid. The CCM scanning parameters utilised to move tip between points were
$V$ = -2.5 V and $I$ = 70 pA. The voltage has been swept between $\pm$ 2 V.
Fig. 2(c) shows the statistical analysis of the grid STS. The inset of Fig.
2(c) illustrate averaged dI/dV spectrum of the whole grid. These STS
measurements demonstrate that the majority of the measured local band gaps
vary from roughly 0.75 eV to 1.35 eV, with a statistical mean of (1.05 $\pm$
0.29) eV. We can fit this with three Gaussian distributions as shown in Fig.
2(c). There are a number of possible effects that could contribute to the
width of the distribution. Excess oxygen could cause an increase in the local
surface band gap as measured by STS, and we do indeed observe a slightly lower
Nb:O ratio of 31%:69% from in-situ XPS measurements. Therefore, we suggest
that the areas of higher band gap likely arise due to an excess of oxygen on
and near the surface. There are a number of ways this excess oxygen could be
incorporated at the surface. Some reports of thin film growth show both Nb4+
and Nb5+ present at the surface Posadas  _et al._ (2014); Páez Fajardo  _et
al._ (2021), indicating the presence of some higher oxides in parts of the
surface Zhussupbekov  _et al._ (2020b). However, for our surface as scanned,
core level XPS on the Nb 3d state shows only an Nb4+ component, indicating no
higher valence oxide layer on the surface. Analysis of the oxygen 1s core
level state shows two chemical states of oxygen. Further sputter anneal cycles
reduce the intensity of the secondary oxygen peak (see SI Fig. 1(d-e)),
suggesting that this oxygen is more weakly bound to the surface. Therefore we
believe that this component likely comes from oxygen adatoms/hydroxide groups
at the surface, which have been observed on other rutile (110) systems such as
TiO2 Epling  _et al._ (1998); Zhang  _et al._ (2019); Skelton and Walsh
(2007); Wen  _et al._ (2020). LEED also supports this viewpoint (Fig. SI1(a)),
as we see no evidence of a major surface reconstruction caused by the
formation of additional oxides. Instead, LEED measurements confirm the long
range crystalline order with the expected distorted rutile (100)/rutile (110)
like symmetry.

###  III.2 B. Band widths

The initially measured valence band XPS (Fig. 3(a), black curve) shows
characteristics similar to that observed in Ref. Posadas  _et al._ (2014) with
low counts per second, and a diminished intensity of the Nb 4$d$ peak. We
believe that this is due to oxygen contamination at the surface, likely in the
form of oxygen adatoms due to the lack of any other Nb components observed in
core-level XPS. Further sputter anneal cycles on our film dramatically improve
the counts per second, and we see an increase in the intensity of the Nb 4$d$
peak (accompanied by a decrease in the secondary oxygen component in core-
level XPS). This type of spectra (Fig. 3(a), green curve) is consistent with
that observed in Refs. Wahila  _et al._ (2019); Wong  _et al._ (2014); O’Hara
_et al._ (2014). Across both of the observed spectra we see values for the
band widths and $p-d$ gap that are generally consistent with range of values
from previous reports (see Table 3). While in some works a sharp $p-d$ gap is
observed, we do not see this in our spectra. This could be due to weak
disorder and/or the remaining oxygen contamination not removed by sputter
anneal cycles. Ultra-violet photoelectron spectroscopy (UPS) measurements
(Fig. 3(b)) show that the surface exhibits metallic behaviour.

Figure 3:  (a) Measured valence band XPS for the NbO2(100) film before and
after further sputter anneal cycles, (b) UPS measured after further sputter
anneal cycles, (c) LDA and DFT+$U$+$J$ calculated DOS for the bulk, (d) LDA
and DFT+$U$+$J$ calculated DOS for the NbO2(100) slab.

Since XPS is a surface sensitive technique, it is possible that modifications
of the electronic structure at the surface could contribute to the lack of
agreement between the experiment and theory. To this end, in addition to the
bulk calculations, we also consider a (100)${}_{\textrm{DR}}$ surface, using a
slab with 8 NbO2 layers (24 atomic layers, symmetric to ensure no spurious
dipole moment). Fig. 3(c-d) shows the calculated DOS for the bulk and surface
respectively, for LDA and O 2p. Table 3 shows the calculated band widths and
$p$-$d$ gaps for our calculations, along with experimental and theoretical
values from the literature.

For the bulk calculations, we see that when $U$ and $J$ corrections are
applied to both the Nb 4$d$ and O 2$p$ subspaces, we get values for the $p$
bandwidth that are comparable with experiment. The $p$-$d$ gap is also
improved when applying the correction to both atoms, though overshooting the
experimental values. This reinforces the findings of Ref. Orhan and O’Regan
(2020), namely that for a charge-transfer gap like this $p$-$d$ sub gap,
appropriate corrections for both subspaces must be included to obtain the
correct value for that gap. For the 4$d$ bandwidth, there is marginal
improvement when including the $U$ and $J$ correction, but still well shy of
the experimentally measured values.

Moving to the slab calculations, we see that there is a clear non-zero DOS at
the Fermi level. This leads to an increase the Nb 4dxy band width, and as a
result slightly decreases the p-d gap. Unsurprisingly, the O 2p band width is
largely unaffected by this effect. We see that overall, the
DFT+$U^{d,p}$+$J^{d,p}$ slab calculation presents solid agreement with the
range of experimentally measured values. The Nb 4dxy width is increased to 1.5
eV, in line with experiment. The p-d gap is now 2.2 eV, much closer to the
experimental value. The O 2p width is unchanged from the bulk calculation of
the same functional, remaining in good agreement with experiment. It therefore
seems that the underestimation of the calculated bulk Nb 4dxy band width
compared to the measured value is not purely related to the to the exchange-
correlation functional, but also due to metallisation at the surface.

Table 3:  Summary of experimental and theoretical values for the valence band XPS spectral features. Results from the most elaborate theoretical treatment, including both surface effects, as well as Hubbard $U$ and Hund $J$ corrections, are highlighted in bold. Method |  p-d Gap [eV] |  p bandwidth [eV] |  d bandwidth [eV]
---|---|---|---
LDA O’Hara  _et al._ (2014) | 1.58 | 5.55 | 0.70
HSE O’Hara  _et al._ (2014) | 1.79 | 5.76 | 0.71
GW-BSE Kulmus  _et al._ (2021) | 2.4 | 5.7 | 0.79
|  |  |
LDA bulk (this work) | 1.51 | 5.65 | 0.80
DFT+$U^{d}$ bulk (this work) | 1.16 | 5.70 | 0.91
DFT+$U^{d}_{\mathrm{eff}}$ bulk (this work) | 1.17 | 5.70 | 0.88
DFT+$U^{d}$+$J^{d}$ bulk (this work) | 1.20 | 5.71 | 0.86
DFT+$U^{d,p}$ bulk (this work) | 2.80 | 6.90 | 0.87
DFT+$U^{d,p}_{\mathrm{eff}}$ bulk (this work) | 2.64 | 6.65 | 0.88
DFT+$U^{d,p}$+$J^{d,p}$ bulk (this work) | 2.31 | 6.40 | 0.85
LDA (100) slab (this work) | 1.48 | 5.65 | 0.95
DFT+$\boldsymbol{U^{d,p}}$+$\boldsymbol{J^{d,p}}$ (100) slab (this work) | 2.20 | 6.40 | 1.50
|  |  |
XPS (this work, initial) |  $\approx$ 2.0 |  $\approx$ 6.2 |  $\approx$ 1.5
XPS (this work, after further sputter anneal cycles) |  $\approx$ 1.7 |  $\approx$ 6.4 |  $\approx$ 1.9
XPS Posadas  _et al._ (2014) |  $\approx$ 1.8 |  $\approx$ 6.3 |  $\approx$ 1.5
XPS O’Hara  _et al._ (2014) |  $\approx$ 2.0 |  $\approx$ 6.2 |  $\approx$ 1.7
XPS Wong  _et al._ (2014) |  $\approx$ 2.0 |  $\approx$ 6.8 |  $\approx$ 1.2
HAXPES Wahila  _et al._ (2019) |  $\approx$ 1.9 |  $\approx$ 6.0 |  $\approx$ 1.5

The optimised geometry of the slab shows that the dimerization collapses
almost entirely for the Nb atoms under the bridging oxygens (blue highlighted
rows in Fig. 4(a)). The dimerization is also reduced relative to the bulk for
the 5-fold Nb atoms and those in the 2${}^{\textrm{nd}}$ layer of the slab.
Fig. 4(d) shows the dimerization plotted as a function of depth into the
surface. It may be tempting then to ascribe the metallisation of the surface
to this breakdown in dimerization, however there is another competing effect
which can change the surface band structure, namely the reduced Coulomb
repulsion at the surface. Fig. 4(b-c) show the layer resolved DOS for the slab
without and with ionic relaxation. We see that even when the ionic geometry is
kept fixed (no breakdown in dimerization allowed) metallicity still emerges at
the surface.

In order to further quantify these two effects, we plot the DOS, decomposed by
magnetic quantum number $m_{l}$ for the Nb atoms in the slab. We utilise the
method developed in Ref. Aarons  _et al._ (2019), and use pseudo-atomic
orbitals as projectors. The coordinate system chosen for the basis functions
is such that the dimerization lies in the $xy$ plane. We first note that for
the bulk Nb atoms (Fig. 4(e) in the slab our ml projected DOS resembles that
of Ref. Kulmus  _et al._ (2021), with the $d_{xy}$/$d_{xy}^{*}$ splitting
caused by dimerization, with degenerate $d_{yz}$ and $d_{xz}$ orbitals. Now
looking at the plots for the surface atoms, we see that for the Nb atoms under
the bridging oxygens (Fig. 4(f)), (dimerization $\simeq$ 0) the
$d_{xy}$/$d_{xy}^{*}$ splitting has disappeared, as is the expected behaviour
from a Peierls type splitting. We also see that the energy of the $d_{yz}$ and
$d_{xz}$ orbitals has decreased, as expected due to the reduction in Coulomb
repulsion at the surface. For the 5-fold surface Nb atoms (Fig. 4(g)), we see
the same reduction in energy of $d_{yz}$ and $d_{xz}$ orbitals, but with the
$d_{xy}$/$d_{xy}^{*}$ splitting still present. Looking at both plots, we
notice that the non-zero DOS at the Fermi level is coming mostly from the
$d_{yz}$ and $d_{xz}$ orbitals, indicating that the lowering in energy of
these orbitals is the underlying cause of the metallicity at the surface.
Finally, we note some of the outstanding discrepancies between experiment and
theory in the p-d gap region, due to neglected effects in the latter. In
particular, referring to the XPS data in Fig. 3(a-b), there is non-zero
spectral weight between -3 and -1 eV, where in the calculated (generalized
Kohn-Sham) spectra there is a clear p-d gap. There are certainly surface
disorder and finite-temperature effects that are neglected in the simulations.
Referring to Fig. 4(b), it is clear that surface ionic displacements can
readily populate this region with spectral features. Furthermore, there are
many-electron effects that cannot be seen at the level of theory used here.
These include plasmon satellite bands, Hubbard bands, and enhanced spectral
splitting due to strong and orientation-dependent intra-orbital interactions
associated exchange and correlation effects. Noting the significant first-
principles Hubbard U for O 2p orbitals, exceeding their band-width, we cannot
rule out the possibility of unconventional strong 2p electron correlation
effects, as well as the more usually anticipated Nb 4d ones.

Figure 4:  (a) Surface crystal structure showing side view (blue panel) and
top down view (black panel) (b-c) Layer resolved local density of states for
the (100) slab with shifted bulk DOS superimposed without and with ionic
relaxation (d) dimerization (=$\frac{\textrm{long bond}}{\textrm{short
bond}}-1)$ versus layer number (e-g) $m_{l}$ projected DOS for bulk atom, and
the two types of surface atom.

##  IV IV. Conclusions

Overall, we find that the $U$ and $J$ values for both Nb 4d and O 2p orbitals
in NbO2 can be straightforwardly calculated due to well behaved linear
response. The inclusion of a correction based on these values can
significantly improve agreement with experimentally measured spectral
quantities, demonstrating the efficiency of this first-principles method in
capturing the essential physics of the material. The band gap and band widths
are relatively robust against the exact form of the correction used, due in
part to the relatively small Hund’s $J$ value for both Nb 4d and O 2p. At
approximately one-tenth of the corresponding Hubbard $U$ values, the low
Hund’s $J$ values indicate an already good treatment of static correlation
effects at the LDA level in this system.

Additionally, in order to explain the measured XPS spectra we must account for
surface effects, due to the emergence of metallic behaviour. It seems that
although the dimerization does break down at the surface, it is not the
driving force behind the transition from insulating to metallic behaviour at
the surface. Therefore it cannot be considered a purely Peierls transition.
Our results show that a transition from insulating to metallic behaviour can
be induced in NbO2 from a change in the Coulomb potential at the surface. This
is the case irrespective of whether corrections for correlation effects at the
DFT+$U$+$J$ level are incorporated. Indeed, due to the relatively small Hund’s
$J$ parameters, it is sufficient to include corrections beyond LDA at the
Hubbard $U$ level alone in order to achieve good agreement with spectroscopic
measurements. Such corrections do not explicitly address correlation within
this formalism, in the strict sense of the word, but only Coulomb repulsion
and exchange. Thus, our results give only partial credence to the reports that
electronic correlation play an important role in this system Rana  _et al._
(2019); Brito  _et al._ (2017), and we cannot rule out agreement with other
reports which claim correlations to be unimportant Wahila  _et al._ (2019);
Kulmus  _et al._ (2021).

We emphasize, however, that our findings pertain primarily to surface effects,
as distinct from bulk effects such as the temperature dependent metal-
insulator transition. We have not, for example, explored the relationship
between bulk dimerization or finite temperature effects with the calculated
$U$ and $J$ parameters, which remains an avenue for future research. Our
measurements also show the presence of excess oxygen on and near the surface,
which could also play a role in modifying the observed spectra. However, since
we only see evidence of oxygen adatoms and not a continuous layer of a higher
valence oxide, we do not expect that these adatoms could account for the
change we predict between the surface and bulk properties (e.g., the 0.5 eV
widening of the Nb dxy band width observed in Fig. 3(c-d)).

In order to further confirm that the metallicity is confined to the surface,
the surface and bulk contributions to the valence band photoelectron signal
would need to be decoupled, using a technique sensitive to both the bulk and
surface. Hard X-ray photoelectron Spectroscopy measurements were reported in
Ref. Wahila  _et al._ (2019), but only for a single photon energy. Multiple
measurements varying the photon energy would be needed to decouple the surface
and bulk contributions, and to determine whether or not the observed Nb dxy
band width comes mostly from the surface effect calculated here, or a
combination of bulk and surface. Angle resolved XPS/UPS could also shed light
on this effect.

This metallisation at the surface could give rise to new applications and new
physics. Anisotropic 2D electron gases have attracted significant attention
recent years due to their novel properties such as 2D superconductivity,
unique magnetic phases, and nonlinear transport characteristics Annadi  _et
al._ (2013); Lev  _et al._ (2018); Lee  _et al._ (2012); Wang  _et al._
(2014). With the top two layers near the surface becoming metallic, while the
bulk remains insulating, a 2D electron gas may form. Additionally, NbO2 is
known to have significant anisotropy in its conductivity, with the lowest
resistance along the dimerization axis Sakai  _et al._ (1985); Belanger  _et
al._ (1974); Stoever  _et al._ (2020). Due to the partial breakdown of
dimerization at the surface, we would expect the anisotropy of the
conductivity to be even higher for these metallic surface electrons, giving
rise to a highly anisotropic 2D electron gas.

##  V Acknowledgements

This work was supported by the Irish Research Council (IRC) Laureate Award
(IRCLA/2019/171), Science Foundation Ireland. The authors wish to acknowledge
the Irish Centre for High-End Computing (ICHEC) for the provision of
computational facilities and support, as well as the Trinity Centre for High
Performance Computing (Boyle cluster). K.Z would like to acknowledge support
from the Ministry of Science and Higher Education of the Republic of
Kazakhstan (No. AP19174589). This work was further supported (D.D.O’R) by
Science Foundation Ireland (No. 19/EPSRC/3605) and the Engineering and
Physical Sciences Research Council (No. EP/S030263/1). A.Z. and K.Z. would
also like to acknowledge funding from IRC through No. GOIPD/2022/443 and No.
GOIPD/2022/774 awards. J.E.B. and J.S. acknowledge support from the Leibniz
Competition project with the title “Physics and Control of Defects in Oxide
Films for Adaptive Electronics”. D.D.O’R thanks N.D.M Hine and E.B. Linscott
for helpful discussions. We acknowledge and thank the ONETEP Developers Group.

##  VI Bibliography

## References

  * Kim  _et al._ (2013) S. Kim, J. Park, J. Woo, C. Cho, W. Lee, J. Shin, G. Choi, S. Park, D. Lee, B. H. Lee, and H. Hwang, Microelectronic Engineering 107, 33 (2013).
  * Cha  _et al._ (2016) E. Cha, J. Park, J. Woo, D. Lee, A. Prakash, and H. Hwang, Applied Physics Letters 108, 153502 (2016).
  * Strelcov  _et al._ (2009) E. Strelcov, Y. Lilach, and A. Kolmakov, Nano Letters 9, 2322 (2009).
  * Chen  _et al._ (2019) J. Chen, W. Mao, B. Ge, J. Wang, X. Ke, V. Wang, Y. Wang, M. Döbeli, W. Geng, H. Matsuzaki, J. Shi, and Y. Jiang, Nature Communications 10, 1 (2019).
  * Kumar  _et al._ (2017) S. Kumar, J. P. Strachan, and R. S. Williams, Nature 548, 318 (2017).
  * Pickett  _et al._ (2013) M. D. Pickett, G. Medeiros-Ribeiro, and R. S. Williams, Nature Materials 12, 114 (2013).
  * Wahila  _et al._ (2019) M. J. Wahila, G. Paez, C. N. Singh, A. Regoutz, S. Sallis, M. J. Zuba, J. Rana, M. B. Tellekamp, J. E. Boschker, T. Markurt, J. E. N. Swallow, L. A. H. Jones, T. D. Veal, W. Yang, T.-L. Lee, F. Rodolakis, J. T. Sadowski, D. Prendergast, W.-C. Lee, W. A. Doolittle, and L. F. J. Piper, Phys. Rev. Materials 3, 074602 (2019).
  * O’Hara  _et al._ (2014) A. O’Hara, T. N. Nunley, A. B. Posadas, S. Zollner, and A. A. Demkov, J. Appl. Phys. 116, 213705 (2014).
  * Sakai  _et al._ (1985) Y. Sakai, N. Tsuda, and T. Sakata, J. Phys. Soc. Jpn. 54, 1514 (1985).
  * Posadas  _et al._ (2014) A. B. Posadas, A. O’Hara, S. Rangan, R. A. Bartynski, and A. A. Demkov, Appl. Phys. Lett. 104, 092901 (2014).
  * Stoever  _et al._ (2020) J. Stoever, J. E. Boschker, S. Bin Anooz, M. Schmidbauer, P. Petrik, J. Schwarzkopf, M. Albrecht, and K. Irmscher, Appl. Phys. Lett. 116, 182103 (2020).
  * Brito  _et al._ (2017) W. H. Brito, M. C. O. Aguiar, K. Haule, and G. Kotliar, Phys. Rev. B 96, 195102 (2017).
  * Eyert (2002) V. Eyert, EPL 58, 851 (2002).
  * Kulmus  _et al._ (2021) K. Kulmus, S. Gemming, M. Schreiber, D. Pashov, and S. Acharya, Phys. Rev. B 104, 035128 (2021).
  * SI (1) See supplementary information .
  * Fairley  _et al._ (2021) N. Fairley, V. Fernandez, M. Richard‐Plouet, C. Guillot-Deudon, J. Walton, E. Smith, D. Flahaut, M. Greiner, M. Biesinger, S. Tougaard, D. Morgan, and J. Baltrusaitis, Applied Surface Science Advances 5, 100112 (2021).
  * Zhussupbekov  _et al._ (2020a) K. Zhussupbekov, C. P. Cullen, A. Zhussupbekova, I. V. Shvets, G. S. Duesberg, N. McEvoy, and C. Ó. Coileáin, RSC Advances 10, 42001 (2020a).
  * Zhussupbekov  _et al._ (2021a) K. Zhussupbekov, L. Ansari, J. B. McManus, A. Zhussupbekova, I. V. Shvets, G. S. Duesberg, P. K. Hurley, F. Gity, C. Ó. Coileáin, and N. McEvoy, npj 2D Materials and Applications 5 (2021a), 10.1038/s41699-020-00196-8.
  * Berman  _et al._ (2023) S. Berman, A. Zhussupbekova, B. Walls, K. Walshe, S. I. Bozhko, A. Ionov, D. D. O’Regan, I. V. Shvets, and K. Zhussupbekov, Phys. Rev. B 107, 165425 (2023).
  * Zhussupbekov  _et al._ (2021b) K. Zhussupbekov, K. Walshe, B. Walls, A. Ionov, S. I. Bozhko, A. Ksenz, R. N. Mozhchil, A. Zhussupbekova, K. Fleischer, S. Berman, I. Zhilyaev, D. D. O’Regan, and I. V. Shvets, The Journal of Physical Chemistry C 125, 5549 (2021b).
  * Momma and Izumi (2011) K. Momma and F. Izumi, Journal of Applied Crystallography 44, 1272 (2011).
  * (22) “Opium : The optimized pseudopotential interface unification module, http://opium.sourceforge.net/.” .
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, and et al, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, and et al, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Skylaris  _et al._ (2005) C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, The Journal of Chemical Physics 122, 084119 (2005).
  * Prentice  _et al._ (2020) J. C. A. Prentice, J. Aarons, J. C. Womack, A. E. A. Allen, L. Andrinopoulos, L. Anton, R. A. Bell, A. Bhandari, G. A. Bramley, R. J. Charlton, R. J. Clements, D. J. Cole, G. Constantinescu, F. Corsetti, S. M.-M. Dubois, K. K. B. Duff, J. M. Escartín, A. Greco, Q. Hill, L. P. Lee, E. Linscott, D. D. O’Regan, M. J. S. Phipps, L. E. Ratcliff, Álvaro Ruiz Serrano, E. W. Tait, G. Teobaldi, V. Vitale, N. Yeung, T. J. Zuehlsdorff, J. Dziedzic, P. D. Haynes, N. D. M. Hine, A. A. Mostofi, M. C. Payne, and C.-K. Skylaris, The Journal of Chemical Physics 152, 174111 (2020).
  * O’Regan  _et al._ (2012) D. D. O’Regan, N. D. M. Hine, M. C. Payne, and A. A. Mostofi, Phys. Rev. B 85, 085107 (2012).
  * Lloyd-Williams and Monserrat (2015) J. H. Lloyd-Williams and B. Monserrat, Physical Review B 92, 184301 (2015).
  * Himmetoglu  _et al._ (2011) B. Himmetoglu, R. M. Wentzcovitch, and M. Cococcioni, Phys. Rev. B 84, 115108 (2011).
  * Lim  _et al._ (2016) J. S. Lim, D. Saldana-Greco, and A. M. Rappe, Phys. Rev. B 94, 165151 (2016).
  * Linscott  _et al._ (2018) E. B. Linscott, D. J. Cole, M. C. Payne, and D. D. O’Regan, Physical Review B 98, 235157 (2018).
  * Orhan and O’Regan (2020) O. K. Orhan and D. D. O’Regan, Physical Review B 101, 245137 (2020).
  * Winczewski  _et al._ (2022) S. Winczewski, J. Dziedzic, T. Miruszewski, J. Rybicki, and M. Gazda, The Journal of Physical Chemistry C 126, 18439 (2022).
  * Burgess  _et al._ (2023) A. C. Burgess, E. Linscott, and D. D. O’Regan, Phys. Rev. B 107, L121115 (2023).
  * Constantinescu and Hine (2015) G. C. Constantinescu and N. D. M. Hine, Phys. Rev. B 91, 195416 (2015).
  * Popescu and Zunger (2012) V. Popescu and A. Zunger, Phys. Rev. B 85, 085201 (2012).
  * Graham  _et al._ (2020) A. J. Graham, J. Zultak, M. J. Hamer, V. Zolyomi, S. Magorrian, A. Barinov, V. Kandyba, A. Giampietri, A. Locatelli, F. Genuzio, N. C. Teutsch, C. Salazar, N. D. M. Hine, V. I. Fal’ko, R. V. Gorbachev, and N. R. Wilson, 2D Materials 8, 015016 (2020).
  * Loh  _et al._ (2021) S. M. Loh, X. Xia, N. R. Wilson, and N. D. M. Hine, Phys. Rev. B 103, 245410 (2021).
  * Kirchner-Hall  _et al._ (2021) N. E. Kirchner-Hall, W. Zhao, Y. Xiong, I. Timrov, and I. Dabo, Applied Sciences 11 (2021), 10.3390/app11052395.
  * Lambert and O’Regan (2023) D. S. Lambert and D. D. O’Regan, Phys. Rev. Res. 5, 013160 (2023).
  * Moore  _et al._ (2022) G. C. Moore, M. K. Horton, A. M. Ganose, M. Siron, E. Linscott, D. D. O’Regan, and K. A. Persson, “High-throughput determination of hubbard u and hund j values for transition metal oxides via linear response formalism,” (2022), arXiv:2201.04213 .
  * Bolzan  _et al._ (1994) A. A. Bolzan, C. Fong, B. J. Kennedy, and C. J. Howard, Journal of Solid State Chemistry 113, 9 (1994).
  * Dhamdhere  _et al._ (2016) A. R. Dhamdhere, T. Hadamek, A. B. Posadas, A. A. Demkov, and D. J. Smith, Journal of Applied Physics 120, 245302 (2016).
  * Páez Fajardo  _et al._ (2021) G. J. Páez Fajardo, S. A. Howard, E. Evlyukhin, M. J. Wahila, W. R. Mondal, M. Zuba, J. E. Boschker, H. Paik, D. G. Schlom, J. T. Sadowski, S. A. Tenney, B. Reinhart, W.-C. Lee, and L. F. J. Piper, Chemistry of Materials 33, 1416 (2021).
  * O’Hara and Demkov (2015) A. O’Hara and A. A. Demkov, Phys. Rev. B 91, 094305 (2015).
  * Pynn  _et al._ (1976) R. Pynn, J. D. Axe, and R. Thomas, Phys. Rev. B 13, 2965 (1976).
  * Pynn and Axe (1976) R. Pynn and J. D. Axe, Journal of Physics C: Solid State Physics 9, L199 (1976).
  * Lee and Durand (1984) J. C. Lee and W. W. Durand, J. Appl. Phys. 56, 3350 (1984).
  * Dudarev  _et al._ (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * Huang  _et al._ (2011) H. Huang, J.-T. Sun, Y. P. Feng, W. Chen, and A. T. S. Wee, Physical Chemistry Chemical Physics 13, 20933 (2011).
  * Wong  _et al._ (2014) F. J. Wong, N. Hong, and S. Ramanathan, Phys. Rev. B 90, 115135 (2014).
  * Kurnosikov  _et al._ (2001) O. Kurnosikov, F. C. de Nooij, P. LeClair, J. T. Kohlhepp, B. Koopmans, H. J. M. Swagten, and W. J. M. de Jonge, Phys. Rev. B 64, 153407 (2001).
  * Stavale  _et al._ (2012) F. Stavale, N. Nilius, and H.-J. Freund, New Journal of Physics 14, 033006 (2012).
  * Ashworth and Thornton (2001) T. Ashworth and G. Thornton, Thin Solid Films 400, 43 (2001).
  * Shi  _et al._ (2020) H. Shi, H. Liu, Z. Li, W. Wang, and X. Shao, Physical Chemistry Chemical Physics 22, 27077 (2020).
  * Zhussupbekov  _et al._ (2020b) K. Zhussupbekov, K. Walshe, S. I. Bozhko, A. Ionov, K. Fleischer, E. Norton, A. Zhussupbekova, V. Semenov, I. V. Shvets, and B. Walls, Scientific Reports 10 (2020b).
  * Epling  _et al._ (1998) W. S. Epling, C. H. Peden, M. A. Henderson, and U. Diebold, Surface Science 412-413, 333 (1998).
  * Zhang  _et al._ (2019) Q. Zhang, H. F. Wen, Y. Adachi, M. Miyazaki, Y. Sugawara, R. Xu, Z. H. Cheng, and Y. J. Li, The Journal of Physical Chemistry C (2019).
  * Skelton and Walsh (2007) A. A. Skelton and T. R. Walsh (Taylor & Francis Group, 2007) pp. 379–389.
  * Wen  _et al._ (2020) H. F. Wen, Y. Sugawara, and Y. J. Li, Nanomaterials 10, 1506 (2020).
  * Aarons  _et al._ (2019) J. Aarons, L. G. Verga, N. D. M. Hine, and C.-K. Skylaris, Electronic Structure 1, 035002 (2019).
  * Rana  _et al._ (2019) R. Rana, J. M. Klopf, J. Grenzer, H. Schneider, M. Helm, and A. Pashkin, Phys. Rev. B 99, 041102 (2019).
  * Annadi  _et al._ (2013) A. Annadi, Q. Zhang, X. R. Wang, N. Tuzla, K. Gopinadhan, W. M. Lü, A. R. Barman, Z. Q. Liu, A. Srivastava, S. Saha, Y. L. Zhao, S. W. Zeng, S. Dhar, E. Olsson, B. Gu, S. Yunoki, S. Maekawa, H. Hilgenkamp, T. Venkatesan, and A. Ariando, Nature Communications 4, 1 (2013).
  * Lev  _et al._ (2018) L. L. Lev, I. O. Maiboroda, M. A. Husanu, E. S. Grichuk, N. K. Chumakov, I. S. Ezubchenko, I. A. Chernykh, X. Wang, B. Tobler, T. Schmitt, M. L. Zanaveskin, V. G. Valeyev, and V. N. Strocov, Nature Communications 9, 1 (2018).
  * Lee  _et al._ (2012) S. W. Lee, Y. Liu, J. Heo, and R. G. Gordon, Nano Letters 12, 4775 (2012).
  * Wang  _et al._ (2014) Z. Wang, Z. Zhong, X. Hao, S. Gerhold, B. Stöger, M. Schmid, J. Sánchez-Barriga, A. Varykhalov, C. Franchini, K. Held, and U. Diebold, Proceedings of the National Academy of Sciences of the United States of America 111, 3933 (2014).
  * Belanger  _et al._ (1974) G. Belanger, J. Destry, G. Perluzzo, and P. M. Raccah, Canadian Journal of Physics 52, 2272 (1974).

Supporting Information
Reconciling the theoretical and experimental electronic structure of NbO2

The NbO2 film was introduced to the UHV chamber and went through Ar ion
sputtering and annealing cycles in a UHV environment at
$700\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 2 h. The LEED measurement of
the synthesised NbO2 film reveals a long-range crystalline order as presented
in Fig.SI1(a). The LEED pattern in Fig.SI1(a) was obtained with a primary beam
energy of 90 eV, at room temperature. It shows a number of fairly sharp and
well-resolved diffraction reflexes with two-fold symmetry. The directions of
the reciprocal $\overrightarrow{a}^{*}$ and $\overrightarrow{b}^{*}$ vectors
are indicated with blue colour on top of the LEED pattern Huang  _et al._
(2011).

Fig.SI1(b) shows the side and top views of the relaxed crystal geometry of the
NbO2 surface. The top-view of the surface is displaying the (100) plane of the
NbO2. The real-space vectors of the NbO2 ($\overrightarrow{a}$ and
$\overrightarrow{b}$) are shown in the schematic crystal structure model in
blue colour for a better illustration.

Fig.SI1(c-e) present the core level XPS of Nb 3d and O 1s levels,
respectively. We can see for the Niobium core level in (c) that only one
component corresponding to Nb4+ is present Wong  _et al._ (2014). The lack of
any Nb5+ contribution (as has been reported in some other studies O’Hara  _et
al._ (2014); Páez Fajardo  _et al._ (2021); Wahila  _et al._ (2019)) confirms
that we are probing solely the NbO2 surface structure. The initially measured
oxygen core level shows two components (d), with the secondary component
decreasing after further sputter anneal cycles (e). This suggests that the
secondary component is likely arising from hydroxide groups or oxygen adatoms
weakly bound to the surface.

Fig.SI2 shows the RHEED and AFM images of the (100) film. The sharp reflexes
and streak pattern demonstrates the high crystalline quality of the film,
while the AFM images show that the film displays low roughness. Fig.SI3 shows
the XRD $\theta-2\theta$ scan and rocking curve, demonstrating the high
quality single crystalline nature of the film. All in all, these measurements
of RHEED, AFM, LEED, and XRD confirm the high crystalline quality of the
epitaxial NbO2 film, with a (100) distorted rutile (or (110) rutile) surface
orientation.

Figure SI1:  (a) Measured LEED pattern for NbO2(100) film at E=90eV. (b-c)
Side and top view of fully relaxed crystal structure. (d-e) Core level XPS of
Nb 3d and O 1s levels respectively.  Figure SI2:  (a) RHEED image after
deposition, (b-d) AFM images of (100) film.  Figure SI3:  (a) $\theta-2\theta$
scan, (b) rocking curve.


# Delay Minimization for Federated Learning Over
Wireless Communication Networks

Zhaohui Yang  Mingzhe Chen  Walid Saad  Choong Seon Hong  Mohammad Shikh-
Bahaei  H. Vincent Poor  Shuguang Cui

###### Abstract

In this paper, the problem of delay minimization for federated learning (FL)
over wireless communication networks is investigated. In the considered model,
each user exploits limited local computational resources to train a local FL
model with its collected data and, then, sends the trained FL model parameters
to a base station (BS) which aggregates the local FL models and broadcasts the
aggregated FL model back to all the users. Since FL involves learning model
exchanges between the users and the BS, both computation and communication
latencies are determined by the required learning accuracy level, which
affects the convergence rate of the FL algorithm. This joint learning and
communication problem is formulated as a delay minimization problem, where it
is proved that the objective function is a convex function of the learning
accuracy. Then, a bisection search algorithm is proposed to obtain the optimal
solution. Simulation results show that the proposed algorithm can reduce delay
by up to 27.3% compared to conventional FL methods.

Machine Learning, ICML

##  1 Introduction

In future wireless systems, due to privacy constraints and limited
communication resources for data transmission, it is impractical for all
wireless devices to transmit all of their collected data to a data center that
can implement centralized machine learning algorithms for data analysis (Wang
et al., 2018; Chen et al., 2019a; Huang et al., 2020; Dong et al., 2019; Gao
et al., 2020). To this end, distributed edge learning approaches, such as
federated learning (FL), were proposed (Saad et al., to appear, 2020; Park et
al., 2019; Chen et al., 2020; Samarakoon et al., 2018; Gündüz et al., 2019;
Chen et al., 2019b). In FL, the wireless devices individually establish local
learning models and cooperatively build a global learning model by uploading
the local learning model parameters to a base station (BS) instead of sharing
training data(McMahan et al., 2016; Yang et al., 2020; Wang et al., 2019). To
implement FL over wireless networks, the wireless devices must transmit their
local training results over wireless links (Zhu et al., 2018a), which can
affect the FL performance, because both local training and wireless
transmission introduce delay. Hence, it is necessary to optimize the delay for
wireless FL implementation.

Some of the challenges of FL over wireless networks have been studied in (Zhu
et al., 2018b; Ahn et al., 2019; Yang et al., 2018; Zeng et al., 2019; Chen et
al., 2019; Tran et al., 2019). To minimize latency, a broadband analog
aggregation multi-access scheme for FL was designed in (Zhu et al., 2018b).
The authors in (Ahn et al., 2019) proposed an FL implementation scheme between
devices and access point over Gaussian multiple-access channels. To improve
the statistical learning performance for on-device distributed training, the
authors in (Yang et al., 2018) developed a sparse and low-rank modeling
approach. The work in in (Zeng et al., 2019) proposed an energy-efficient
strategy for bandwidth allocation with the goal of reducing devices’ sum
energy consumption while meeting the required learning performance. However,
the prior works (Konečnỳ et al., 2016; Zhu et al., 2018b; Ahn et al., 2019;
Yang et al., 2018; Zeng et al., 2019) focused on the delay/energy consumption
for wireless consumption without considering the delay/energy tradeoff between
learning and transmission. Recently, in (Chen et al., 2019) and (Tran et al.,
2019), the authors considered both local learning and wireless transmission
energy. In (Chen et al., 2019), the authors investigated the FL loss function
minimization problem with taking into account packet errors over wireless
links. However, this prior work ignored the computation delay of local FL
model. The authors in (Tran et al., 2019) considered the sum learning and
transmission energy minimization problem for FL, where all users transmit
learning results to the BS. However, the solution in (Tran et al., 2019)
requires all users to upload their learning model synchronously.

The main contribution of this paper is a framework for optimizing FL over
wireless networks. In particular, we consider a wireless-powered FL algorithm
in which each user locally computes its FL model parameters under a given
learning accuracy and the BS broadcasts the aggregated FL model parameters to
all users. Considering the tradeoff between local computation delay and
wireless transmission delay, we formulate a joint transmission and computation
optimization problem aiming to minimize the delay for FL. We theoretically
show that the delay is a convex function of the learning accuracy. Based on
the theoretical finding, we propose a bisection-based algorithm to obtain the
optimal solution.

##  2 System Model and Problem Formulation

Consider a cellular network that consists of one BS serving a set
$\mathcal{K}$ of $K$ users, as shown in Fig. 1. Each user $k$ has a local
dataset $\mathcal{D}_{k}$ with $D_{k}$ data samples. For each dataset
$\mathcal{D}_{k}=\\{\boldsymbol{x}_{kl},y_{kl}\\}_{l=1}^{D_{k}}$,
$\boldsymbol{x}_{kl}\in\mathbb{R}^{d}$ is an input vector of user $k$ and
$y_{kl}$ is its corresponding output111For simplicity, this paper only
considers an FL algorithm with a single output. Our approach can be extended
to the case with multiple outputs (Konečnỳ et al., 2016)..

Figure 1: FL over wireless communication networks.

###  2.1 FL Model

For FL, we define a vector $\boldsymbol{w}$ to capture the parameters related
to the global FL model that is trained by all datasets. Hereinafter, the FL
model that is trained by all users’ data set is called _global FL model_ ,
while the FL model that is trained by each user’s dataset is called _local FL
model_. We introduce the loss function
$f(\boldsymbol{w},\boldsymbol{x}_{kl},y_{kl})$, that captures the FL
performance over input vector $\boldsymbol{x}_{kl}$ and output $y_{kl}$. For
different learning tasks, the loss function will be different. Since the
dataset of user $k$ is $\mathcal{D}_{k}$, the total loss function of user $k$
will be:

$\vspace{-1em}F_{k}(\boldsymbol{w})=\frac{1}{D_{k}}\sum_{l=1}^{D_{k}}f(\boldsymbol{w},\boldsymbol{x}_{kl},y_{kl}).$
(1)

In order to deploy FL, it is necessary to train the underlying model. Training
is done in order to compute the global FL model for all users without sharing
their local datasets due to privacy and communication issue. The FL training
problem can be formulated as follows (Wang et al., 2018):

$\vspace{-0.5em}\min_{\boldsymbol{w}}F(\boldsymbol{w})=\sum_{k=1}^{K}\frac{D_{k}}{D}F_{k}(\boldsymbol{w})=\frac{1}{D}\sum_{k=1}^{K}\sum_{l=1}^{D_{k}}f(\boldsymbol{w},\boldsymbol{x}_{kl},y_{kl}),$
(2)

where $D=\sum_{k=1}^{K}D_{k}$ is the total data samples of all users.

To solve problem (2), we adopt the FL algorithm in (Konečnỳ et al., 2016),
which is summarized in Algorithm 1.

Algorithm 1 FL Algorithm

1: Initialize global regression vector $\boldsymbol{w}^{0}$ and iteration
number $n=0$.

2: repeat

3: Each user $k$ computes $\nabla F_{k}(\boldsymbol{w}^{(n)})$ and sends it to
the BS.

4: The BS computes $\nabla
F(\boldsymbol{w}^{(n)})=\frac{1}{K}\sum_{k=1}^{K}\nabla
F_{k}(\boldsymbol{w}^{(n)})$, which is broadcast to all users.

5: parallel for user $k\in\mathcal{K}$

6: Solve local FL problem (2.1) with a given learning accuracy $\eta$ and the
solution is $\boldsymbol{h}_{k}^{(n)}$.

7: Each user sends $\boldsymbol{h}_{k}^{(n)}$ to the BS.

8: end for

9: The BS computes
$\boldsymbol{w}^{(n+1)}=\boldsymbol{w}^{(n)}+\frac{1}{K}\sum_{k=1}^{K}\boldsymbol{h}_{k}^{(n)}$
and broadcasts the value to all users.

10: Set $n=n+1$.

11: until the accuracy $\epsilon_{0}$ of problem (2) is obtained.

In Algorithm 1, at each iteration of the FL algorithm, each user downloads the
global FL model parameters from the BS for local computing, while the BS
periodically gathers the local FL model parameters from all users and sends
the updated global FL model parameters back to all users. We define
$\boldsymbol{w}^{(n)}$ as the global FL parameter at a given iteration $n$.
Each user computes the local FL problem:

$\displaystyle\vspace{-.5em}\min_{\boldsymbol{h}_{k}\in\mathbb{R}^{d}}\quad
G_{k}$ $\displaystyle(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k})\triangleq
F_{k}(\boldsymbol{w}^{(n)}+\boldsymbol{h}_{k})$ $\displaystyle-(\nabla
F_{k}(\boldsymbol{w}^{(n)})-\xi\nabla
F(\boldsymbol{w}^{(n)}))^{T}\boldsymbol{h}_{k},$ (3)

by using the gradient method with a given accuracy. In problem (2.1), $\xi$ is
a constant value. The solution $\boldsymbol{h}_{k}$ in problem (2.1) means the
updated value of local FL parameter for user $k$ in each iteration, i.e.,
$\boldsymbol{w}^{(n)}+\boldsymbol{h}_{k}$ denotes user $k$’ local FL parameter
at the $n$-th iteration. Since it is hard to obtain the actual optimal
solution of problem (2.1), we obtain a solution of (2.1) with some accuracy.
The solution $\boldsymbol{h}_{k}^{(n)}$ of problem (2.1) at the $n$-th
iteration with accuracy $\eta$ means that

$\displaystyle\vspace{-1em}G_{k}$
$\displaystyle(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k}^{(n)})-G_{k}(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k}^{(n)*})$
$\displaystyle\leq\eta(G_{k}(\boldsymbol{w}^{(n)},\boldsymbol{0})-G_{k}(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k}^{(n)*})),$
(4)

where $\boldsymbol{h}_{k}^{(n)*}$ is the actual optimal solution of problem
(2.1).

In Algorithm 1, the iterative method involves a number of global iterations
(i.e., the value of $n$ in Algorithm 1) to achieve a global accuracy
$\epsilon_{0}$ of global FL model. The solution $\boldsymbol{w}^{(n)}$ of
problem (2) with accuracy $\epsilon_{0}$ means that

$\displaystyle\vspace{-1em}F(\boldsymbol{w}^{(n)})-F(\boldsymbol{w}^{*})\leq\epsilon_{0}(F(\boldsymbol{w}^{(0)})-F(\boldsymbol{w}^{*})),$
(5)

where $\boldsymbol{w}^{*}$ is the actual optimal solution of problem (2).

To analyze the convergence of Algorithm 1, we assume that
$F_{k}(\boldsymbol{w})$ is $L$-Lipschitz continuous and $\gamma$-strongly
convex, i.e.,

$\vspace{-0.5em}\gamma\boldsymbol{I}\preceq\nabla^{2}F_{k}(\boldsymbol{w})\preceq
L\boldsymbol{I},\quad\forall k\in\mathcal{K}.$ (6)

Under assumption (6), we provide the following lemma about convergence rate of
Algorithm 1.

###### Lemma 1

If we run Algorithm 1 with $0<\xi\leq\frac{\gamma}{L}$ for

$\vspace{-.5em}n\geq\frac{a}{1-\eta}\triangleq I_{0},$ (7)

iterations with $a=\frac{2L^{2}}{\gamma^{2}\xi}\ln\frac{1}{\epsilon_{0}}$, we
have
$F(\boldsymbol{w}^{(n)})-F(\boldsymbol{w}^{*})\leq\epsilon_{0}(F(\boldsymbol{w}^{(0)})-F(\boldsymbol{w}^{*}))$.

The proof of Lemma 1 can be found in (Yang et al., 2019). From Lemma 1, we can
find that the number of global iterations $n$ increases with the local
accuracy. This is because more iterations are needed if the local computation
has a low accuracy.

###  2.2 Computation and Transmission Model

The FL procedure between the users and their serving BS consists of three
steps in each iteration: Local computation at each user (using several local
iterations), local FL parameter transmission for each user, and result
aggregation and broadcast at the BS. During the local computation step, each
user calculates its local FL parameters by using its local dataset and the
received global FL parameters.

####  2.2.1 Local Computation

We solve the local learning problem (2.1) by using the gradient method. In
particular, the gradient procedure in the $(i+1)$-th iteration is given by:

$\vspace{-.5em}\boldsymbol{h}_{k}^{(n),(i+1)}=\boldsymbol{h}_{k}^{(n),(i)}-\delta\nabla
G_{k}(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k}^{(n),(i)}),$ (8)

where $\delta$ is the step size, $\boldsymbol{h}_{k}^{(n),(i)}$ is the value
of $\boldsymbol{h}_{k}$ at the $i$-th local iteration with given vector
$\boldsymbol{w}^{(n)}$, and $\nabla
G_{k}(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k}^{(n),(i)})$ is the gradient of
function $G_{k}(\boldsymbol{w}^{(n)},\boldsymbol{h}_{k})$ at point
$\boldsymbol{h}_{k}=\boldsymbol{h}_{k}^{(n),(i)}$. We set the initial solution
$\boldsymbol{h}_{k}^{(n),(0)}=\boldsymbol{0}$.

Next, we provide the number of local iterations needed to achieve a local
accuracy $\eta$ in (2.1). We set $v=\frac{2}{(2-L\delta)\delta\gamma}$.

###### Lemma 2

If we set step $\delta<\frac{2}{L}$ and run the gradient method for $i\geq
v\log_{2}(1/\eta)$ iterations at each user, we can solve local FL problem
(2.1) with an accuracy $\eta$.

The proof of Lemma 2 can be found in paper (Yang et al., 2019). Let $f_{k}$ be
the computation capacity of user $k$, which is measured by the number of CPU
cycles per second. The computation time at user $k$ needed for data processing
is:

$\vspace{-0.5em}\tau_{k}=\frac{vC_{k}D_{k}\log_{2}(1/\eta)}{f_{k}}=\frac{A_{k}\log_{2}(1/\eta)}{f_{k}},\quad\forall
k\in\mathcal{K},$ (9)

where $C_{k}$ (cycles/bit) is the number of CPU cycles required for computing
one sample data at user $k$, $v\log_{2}(1/\eta)$ is the number of local
iterations for each user as given by Lemma 2, and $A_{k}=vC_{k}D_{k}$.

####  2.2.2 Wireless Transmission

After local computation, all users upload their local FL parameters to the BS
via frequency domain multiple access (FDMA). The achievable rate of user $k$
can be given by:

$\vspace{-.5em}r_{k}=b_{k}\log_{2}\left(1+\frac{g_{k}p_{k}}{N_{0}b_{k}}\right),\quad\forall
k\in\mathcal{K},$ (10)

where $b_{k}$ is the bandwidth allocated to user $k$, $p_{k}$ is the transmit
power of user $k$, $g_{k}$ is the channel gain between user $k$ and the BS,
and $N_{0}$ is the power spectral density of the Gaussian noise. Due to the
limited bandwidth, we have $\sum_{k=1}^{K}b_{k}\leq B,$ where $B$ is the total
bandwidth.

In this step, user $k$ needs to upload the local FL parameters to the BS.
Since the dimensions of the vector $\boldsymbol{h}_{k}^{(n)}$ are fixed for
all users, the data size that each user needs to upload is constant, and can
be denoted by $s$. To upload data of size $s$ within transmit time $t_{k}$, we
must have: $t_{k}{r_{k}}\geq s.$

####  2.2.3 Information Broadcast

In this step, the BS aggregates the global prediction model parameters. The BS
broadcasts the global prediction model parameters to all users in the
downlink. Due to the high power of the BS and large downlink bandwidth, we
ignore the downlink time. Note that the local data $\mathcal{D}_{k}$ is not
accessed by the BS, so as to protect the privacy of users, as is required by
FL. The delay of each user includes the local computation time and transmit
time. Based on (7) and (9), the delay $T_{k}$ of user $k$ will be:

$\displaystyle\vspace{-1em}T_{k}=I_{0}(\tau_{k}+t_{k})=\frac{a}{1-\eta}\left(\frac{A_{k}\log_{2}(1/\eta)}{f_{k}}+t_{k}\right).$
(11)

We define $T=\max_{k\in\mathcal{K}}T_{k}$ as the delay for training the whole
FL algorithm.

###  2.3 Problem Formulation

We now pose the delay minimization problem:

$\displaystyle\mathop{\min}_{T,\boldsymbol{t},\boldsymbol{b},\boldsymbol{f},\boldsymbol{p},\eta}\>$
$\displaystyle T$ (12) s.t.
$\displaystyle\frac{a}{1-\eta}\left(\frac{A_{k}\log_{2}(1/\eta)}{f_{k}}+t_{k}\right)\leq
T,\quad\forall k\in\mathcal{K},$ (12a) $\displaystyle
t_{k}b_{k}\log_{2}\left(1+\frac{g_{k}p_{k}}{N_{0}b_{k}}\right)\geq
s,\quad\forall k\in\mathcal{K},$ (12b) $\displaystyle\sum_{k=1}^{K}b_{k}\leq
B,$ (12c) $\displaystyle
0\\!\leq\\!f_{k}\\!\leq\\!f_{k}^{\max},0\\!\leq\\!p_{k}\\!\leq\\!p_{k}^{\max},\>\forall
k\in\mathcal{K},$ (12d) $\displaystyle 0\leq\eta\leq 1,$ (12e) $\displaystyle
t_{k}\geq 0,b_{k}\geq 0,\quad\forall k\in\mathcal{K},$ (12f)

where $\boldsymbol{t}=[t_{1},\cdots,t_{K}]^{T}$,
$\boldsymbol{b}=[b_{1},\cdots,b_{K}]^{T}$,
$\boldsymbol{f}=[f_{1},\cdots,f_{K}]^{T}$, and
$\boldsymbol{p}=[p_{1},\cdots,p_{K}]^{T}$. $f_{k}^{\max}$ and $p_{k}^{\max}$
are, respectively, the maximum local computation capacity and maximum transmit
power of user $k$. (12a) indicates that the execution time of the local tasks
and the transmit time for all users should not exceed the delay of the whole
FL algorithm. The data transmission constraint is given by (12b), while the
bandwidth constraint is given by (12c). (12d) represents the maximum local
computation capacity and transmit power limits of all users. The accuracy
constraint is given by (12e).

##  3 Optimal Resource Allocation

Although the delay minimization problem (12) is nonconvex due to constraints
(12a)-(12b), the globally optimal solution is shown to be obtained by using
the bisection method.

###  3.1 Optimal Resource Allocation

Let
$(T^{*},\boldsymbol{t}^{*},\boldsymbol{b}^{*},\boldsymbol{f}^{*},\boldsymbol{p}^{*},\eta^{*})$
be the optimal solution of problem (12). We provide the following lemma about
the feasibility conditions of problem (12).

###### Lemma 3

Problem (12) with fixed $T<T^{*}$ is always feasible, while problem (12) with
fixed $T>T^{*}$ is infeasible.

Proof: Assume that
$(\bar{T},\bar{\boldsymbol{t}},\bar{\boldsymbol{b}},\bar{\boldsymbol{f}},\bar{\boldsymbol{p}},\bar{\eta})$
is a feasible solution of problem (12) with $T=\bar{T}<T^{*}$. Then, solution
$(\bar{T},\bar{\boldsymbol{t}},\bar{\boldsymbol{b}},\bar{\boldsymbol{f}},\bar{\boldsymbol{p}},\bar{\eta})$
is feasible with lower value of the objective function than solution
$(T^{*},\boldsymbol{t}^{*},\boldsymbol{b}^{*},\boldsymbol{f}^{*},\boldsymbol{p}^{*},\eta^{*})$,
which contradicts the fact that
$(T^{*},\boldsymbol{t}^{*},\boldsymbol{b}^{*},\boldsymbol{f}^{*},\boldsymbol{p}^{*},\eta^{*})$
is the optimal solution. For problem (12) with $T=\bar{T}>T^{*}$, we can
always construct a feasible solution
$(\bar{T},\boldsymbol{t}^{*},\boldsymbol{b}^{*},\boldsymbol{f}^{*},\boldsymbol{p}^{*},\eta^{*})$
to problem (12) by checking all constraints. $\Box$

According to Lemma 3, we can use the bisection method to obtain the optimal
solution of problem (12). Denote

$\vspace{-.5em}T_{\min}=0,T_{\max}=\max_{k\in\mathcal{K}}\frac{2aA_{k}}{f_{k}^{\max}}+\frac{2aKs}{B\log_{2}\left(1+\frac{g_{k}p_{k}^{\max}K}{N_{0}B}\right)}.$
(13)

If $T>T_{\max}$, problem (12) is always feasible by setting
$f_{k}=f_{k}^{\max}$, $p_{k}=p_{k}^{\max}$, $b_{k}=\frac{B}{K}$,
$\eta=\frac{1}{2}$, and

$\vspace{-.5em}t_{k}=\frac{Ks}{B\log_{2}\left(1+\frac{g_{k}p_{k}^{\max}K}{N_{0}B}\right)}.$
(14)

Hence, the optimal $T^{*}$ of problem (12) must lie in the interval
$(T_{\min},T_{\max})$. At each step, the bisection method divides the interval
in two by computing the midpoint $T_{\text{mid}}=(T_{\min}+T_{\max})/2$. There
are now only two possibilities: 1) if problem (12) with $T=T_{\text{mid}}$ is
feasible, we have $T^{*}\in(T_{\min},T_{\text{mid}}]$ and 2) if problem (12)
with $T=T_{\text{mid}}$ is infeasible, we have
$T^{*}\in(T_{\text{mid}},T_{\max})$. The bisection method selects the
subinterval that is guaranteed to be a bracket as the new interval to be used
in the next step. As such an interval that contains the optimal $T^{*}$ is
reduced in width by 50% at each step. The process continues until the interval
is sufficiently small.

With a fixed $T$, we still need to check whether there exists a feasible
solution satisfying constraints (12a)-(12g). From constraints (12a) and (12c),
we can see that it is always efficient to utilize the maximum computation
capacity, i.e., $f_{k}^{*}=f_{k}^{\max},\forall k\in\mathcal{K}$. In addition,
from (12b) and (12d), we can see that minimizing the delay can be done by
having: $p_{k}^{*}=p_{k}^{\max},\forall k\in\mathcal{K}$. Substituting the
maximum computation capacity and maximum transmission power into (12), delay
minimization problem becomes:

$\displaystyle\min_{T,\boldsymbol{t},\boldsymbol{b},\eta}\quad\>$
$\displaystyle T$ (15) s.t. $\displaystyle
t_{k}\leq\frac{(1-\eta)T}{a}+\frac{A_{k}\log_{2}\eta}{f_{k}^{\max}},\quad\forall
k\in\mathcal{K},$ (15a) $\displaystyle\frac{s}{t_{k}}\leq
b_{k}\log_{2}\left(1+\frac{g_{k}p_{k}^{\max}}{N_{0}b_{k}}\right),\quad\forall
k\in\mathcal{K},$ (15b) $\displaystyle\sum_{k=1}^{K}b_{k}\leq B,$ (15c)
$\displaystyle 0\leq\eta\leq 1,$ (15d) $\displaystyle t_{k}\geq 0,b_{k}\geq
0,\quad\forall k\in\mathcal{K}.$ (15e)

We provide the sufficient and necessary condition for the feasibility of set
(15a)-(15e) using the following lemma.

###### Lemma 4

With a fixed $T$, set (15a)-(15e) is nonempty if an only if

$\vspace{-.5em}B\geq\min_{0\leq\eta\leq
1}\quad\sum_{k=1}^{K}u_{k}(v_{k}(\eta)),$ (16)

where

$\vspace{-0.5em}u_{k}(\eta)=-\frac{(\ln 2)\eta}{W\left(-\frac{(\ln
2)N_{0}\eta}{g_{k}p_{k}^{\max}}\text{e}^{-\frac{(\ln
2)N_{0}\eta}{g_{k}p_{k}^{\max}}}\right)+\frac{(\ln
2)N_{0}\eta}{g_{k}p_{k}^{\max}}},$ (17)

and

$\vspace{-1em}v_{k}(\eta)=\frac{s}{\frac{(1-\eta)T}{a}+\frac{A_{k}\log_{2}\eta}{f_{k}^{\max}}}.$
(18)

Proof: To prove this, we first define a function
$y=x\ln\left(1+\frac{1}{x}\right)$ with $x>0$. Then, we have

$\vspace{-.5em}y^{\prime}=\ln\left(1+\frac{1}{x}\right)-\frac{1}{x+1},y^{\prime\prime}=-\frac{1}{x(x+1)^{2}}<0.$
(19)

According to (19), $y^{\prime}$ is a decreasing function. Since
$\lim_{t_{i}\rightarrow+\infty}y^{\prime}=0$, we have $y^{\prime}>0$ for all
$0<x<+\infty$. Hence, $y$ is an increasing function, i.e., the right hand side
of (15b) is an increasing function of bandwidth $b_{k}$. To ensure that the
maximum bandwidth constraint (15c) can be satisfied, the left hand side of
(15b) should be as small as possible, i.e., $t_{k}$ should be as long as
possible. Based on (15a), the optimal time allocation should be:

$\vspace{-.5em}t_{k}^{*}=\frac{(1-\eta)T}{a}+\frac{A_{k}\log_{2}\eta}{f_{k}^{\max}},\quad\forall
k\in\mathcal{K}.$ (20)

Substituting (20) into (15b), we can construct the following problem:

$\displaystyle\min_{\boldsymbol{b},\eta}\>$ $\displaystyle\sum_{k=1}^{K}b_{k}$
(21) s.t. $\displaystyle v_{k}(\eta)\leq
b_{k}\log_{2}\left(1+\frac{g_{k}p_{k}^{\max}}{N_{0}b_{k}}\right),\>\forall
k\in\mathcal{K},$ (21a) $\displaystyle 0\leq\eta\leq 1,$ (21b) $\displaystyle
b_{k}\geq 0,\quad\forall k\in\mathcal{K},$ (21c)

where $v_{k}(\eta)$ is defined in (18). We can observe that set (15a)-(15e) is
nonempty if an only if the optimal objective value of (21) is less than $B$.
Since the right hand side of (15b) is an increasing function, (15b) should
hold with equality for the optimal solution of problem (21). Setting (15b)
with equality, problem (21) reduces to (16). $\Box$

To effectively solve (16) in Lemma 4, we provide the following lemma.

###### Lemma 5

In (17), $u_{k}(v_{k}(\eta))$ is a convex function.

Proof: We first prove that $v_{k}(\eta)$ is a convex function. To show this,
we define:

$\vspace{-.5em}\phi(\eta)=\frac{s}{\eta},\quad 0\leq\eta\leq 1,$ (22)

and

$\vspace{-.5em}\varphi_{k}(\eta)={\frac{(1-\eta)T}{a}+\frac{A_{k}\log_{2}\eta}{f_{k}^{\max}}},\quad
0\leq\eta\leq 1.$ (23)

According to (18), we have: $v_{k}(\eta)=\phi(\varphi_{k}(\eta)).$ Then, the
second-order derivative of $v_{k}(\eta)$ can be given by:

$\displaystyle\vspace{-0.5em}v_{k}^{\prime\prime}(\eta)$
$\displaystyle=\phi^{\prime\prime}(\varphi_{k}(\eta))(\varphi_{k}^{\prime}(\eta))^{2}+\phi^{\prime}(\varphi_{k}(\eta))\varphi_{k}^{\prime\prime}(\eta).$
(24)

According to (22) and (23), we have:

$\vspace{-.5em}\phi^{\prime}(\eta)=-\frac{s}{\eta^{2}}\leq
0,\quad\phi^{\prime\prime}(\eta)=\frac{2s}{\eta^{3}}\geq 0,$ (25)

and

$\vspace{-.5em}\varphi_{k}^{\prime\prime}(\eta)=-\frac{A_{k}}{(\ln
2)f_{k}^{\max}\eta^{2}}\leq 0.$ (26)

Combining (24)-(26), we can find that $v_{k}^{\prime\prime}(\eta)\geq 0$,
i.e., $v_{k}(\eta)$ is a convex function.

Then, we can show that $u_{k}(\eta)$ is an increasing and convex function.
According to the proof of Lemma 4, $u_{k}(\eta)$ is the inverse function of
the right hand side of (15b). If we further define function:

$\vspace{-.5em}z_{k}(\eta)=\eta\log_{2}\left(1+\frac{g_{k}p_{k}^{\max}}{N_{0}\eta}\right),\quad\eta\geq
0,$ (27)

$u_{k}(\eta)$ is the inverse function of $z_{k}(\eta)$, which gives
$u_{k}(z_{k}(\eta))=\eta$.

According to (19), function $z_{k}(\eta)$ is an increasing and concave
function, i.e., $z_{k}^{\prime}(\eta)\geq 0$ and
$z_{k}^{\prime\prime}(\eta)\leq 0$. Since $z_{k}(\eta)$ is an increasing
function, its inverse function $u_{k}(\eta)$ is also an increasing function.

Based on the definition of concave function, for any $\eta_{1}\geq 0$,
$\eta_{2}\geq 0$ and $0\leq\theta\leq 1$, we have:

$\vspace{-.5em}z_{k}(\theta\eta_{1}+(1-\theta)\eta_{2})\geq\theta
z_{k}(\eta_{1})+(1-\theta)z_{k}(\eta_{2}).$ (28)

Applying the increasing function $u_{k}(\eta)$ on both sides of (28) yields:

$\vspace{-.5em}\theta\eta_{1}+(1-\theta)\eta_{2}\geq u_{k}(\theta
z_{k}(\eta_{1})+(1-\theta)z_{k}(\eta_{2})).$ (29)

Denote $\bar{\eta}_{1}=z_{k}(\eta_{1})$ and $\bar{\eta}_{2}=z_{k}(\eta_{2})$,
i.e., we have $\eta_{1}=u_{k}(\bar{\eta}_{1})$ and
$\eta_{2}=u_{k}(\bar{\eta}_{2})$. Thus, (29) can be rewritten as:

$\theta u_{k}(\bar{\eta}_{1})+(1-\theta)u_{k}(\bar{\eta}_{1})\geq
u_{k}(\theta\bar{\eta}_{1}+(1-\theta)\bar{\eta}_{2}),$ (30)

which indicates that $u_{k}(\eta)$ is a convex function. As a result, we have
proven that $u_{k}(\eta)$ is an increasing and convex function, which shows:

$\vspace{-.5em}u_{k}^{\prime}(\eta)\geq 0,\quad u_{k}^{\prime\prime}(\eta)\geq
0.$ (31)

To show the convexity of $u_{k}(v_{k}(\eta))$, we have:

$\displaystyle\vspace{-.5em}u_{k}^{\prime\prime}(v_{k}(\eta))=u_{k}^{\prime\prime}(v_{k}(\eta))(v_{k}^{\prime}(\eta))^{2}+u_{k}^{\prime}(v_{k}(\eta))v_{k}^{\prime\prime}(\eta)\geq
0,$

according to $v_{k}^{\prime\prime}(\eta)\geq 0$ and (31). As a result,
$u_{k}(v_{k}(\eta))$ is a convex function. $\Box$

Lemma 5 implies that the optimization problem in (16) is a convex problem,
which can be effectively solved. By finding the optimal solution of (16), the
sufficient and necessary condition for the feasibility of set (15a)-(15e) can
be simplified using the following theorem.

###### Theorem 1

With a fixed $T$, set (15a)-(15e) is nonempty if and only if

$\vspace{-0.5em}B\geq\sum_{k=1}^{K}u_{k}(v_{k}(\eta^{*})),$ (32)

where $\eta^{*}$ is the solution to
$\sum_{k=1}^{K}u_{k}^{\prime}(v_{k}(\eta^{*}))v_{k}^{\prime}(\eta^{*})=0.$

Theorem 1 directly follows from Lemmas 4 and 5. Due to the convexity of
function $u_{k}(v_{k}(\eta))$,
$\sum_{k=1}^{K}u_{k}^{\prime}(v_{k}(\eta^{*}))v_{k}^{\prime}(\eta^{*})$ is an
increasing function of $\eta^{*}$. As a result, the unique solution of
$\eta^{*}$ to
$\sum_{k=1}^{K}u_{k}^{\prime}(v_{k}(\eta^{*}))v_{k}^{\prime}(\eta^{*})=0$ can
be effectively solved via the bisection method.

Algorithm 2 Delay Minimization

1: Initialize $T_{\min}$, $T_{\max}$, and the tolerance $\epsilon_{0}$.

2: repeat

3: Set $T=\frac{T_{\min}+T_{\max}}{2}$.

4: Check the feasibility condition (32).

5: If set (15a)-(15e) has a feasible solution, set $T_{\max}=T$. Otherwise,
set $T_{\min}=T$.

6: until $(T_{\max}-T_{\min})/T_{\max}\leq\epsilon_{0}$.

Based on Theorem 1, the algorithm for obtaining the minimal delay is
summarized in Algorithm 2.

##  4 Simulation Results

For our simulations, we deploy $K=50$ users uniformly in a square area of size
$500$ m $\times$ $500$ m with the BS located at its center. The path loss
model is $128.1+37.6\log_{10}d$ ($d$ is in km) and the standard deviation of
shadow fading is $8$ dB (Yang et al., 2020). In addition, the noise power
spectral density is $N_{0}=-174$ dBm/Hz. We use the real open blog feedback
dataset in (Buza, 2014). This dataset with a total number of 60,021 data
samples originates from blog posts and the dimensional of each data sample is
281. The prediction task associated with the data is the prediction of the
number of comments in the upcoming 24 hours. Parameter $C_{k}$ is uniformly
distributed in $[1,3]\times 10^{4}$ cycles/sample. The effective switched
capacitance in local computation is $\kappa=10^{-28}$. In Algorithm 1, we set
$\xi=1/10$, $\delta=1/10$, and $\epsilon_{0}=10^{-3}$. Unless specified
otherwise, we choose an equal maximum average transmit power
$p_{1}^{\max}=\cdots=p_{K}^{\max}=p^{\max}=10$ dBm, an equal maximum
computation capacity $f_{1}^{\max}=\cdots=f_{K}^{\max}=f^{\max}=2$ GHz, a
transmit data size $s=28.1$ kbits, and a bandwidth $B=20$ MHz. Each user has
$D_{k}=500$ data samples, which are randomly selected from the dataset with
equal probability. All statistical results are averaged over 1000 independent
runs.

Figure 2: Value of the loss function as the number of iterations varies for
convex and nonconvex loss functions.

In Fig. 2, we show the value of the loss function as the number of iterations
varies for convex and nonconvex loss functions. For this feedback prediction
problem, we consider two different loss functions: convex loss function
$f_{1}(\boldsymbol{w},\boldsymbol{x},y)=\frac{1}{2}(\boldsymbol{x}^{T}\boldsymbol{w}-y)^{2}$,
and nonconvex loss function
$f_{2}(\boldsymbol{w},\boldsymbol{x},y)=\frac{1}{2}(\max\\{\boldsymbol{x}^{T}\boldsymbol{w},0\\}-y)^{2}$.
From this figure, we can see that, as the number of iterations increases, the
value of the loss function first decreases rapidly and then decreases slowly
for both convex and nonconvex loss functions. According to Fig. 2, the initial
value of the loss function is $F(\boldsymbol{w}^{(0)})=10^{6}$ and the value
of the loss function decreases to $F(\boldsymbol{w}^{(500)})=1$ for convex
loss function after 500 iterations. For our prediction problem, the optimal
model $\boldsymbol{w}^{*}$ is the one that predicts the output without any
error, i.e., the value of the loss function value should be
$F(\boldsymbol{w}^{*})=0$. Thus, the actual accuracy of the proposed algorithm
is
$\frac{F(\boldsymbol{w}^{(500)})-F(\boldsymbol{w}^{*})}{F(\boldsymbol{w}^{(0)})-F(\boldsymbol{w}^{*})}=10^{-6}$
after 500 iterations. Meanwhile, Fig. 2 clearly shows that the FL algorithm
with a convex loss function can converge faster than that the one having a
nonconvex loss function. According to Fig. 2, the loss function monotonically
decreases as the number of iterations varies for even nonconvex loss function,
which indicates that the proposed FL scheme can also be applied to the
nonconvex loss function.

Figure 3: Delay versus maximum average transmit power of each user.

We compare the proposed FL scheme with the FL FDMA scheme with equal bandwidth
$b_{1}=\cdots=b_{K}$ (labelled as ‘EB-FDMA’), the FL FDMA scheme with fixed
local accuracy $\eta=1/2$ (labelled as ‘FE-FDMA’), and the FL time division
multiple access (TDMA) scheme in (Tran et al., 2019) (labelled as ‘TDMA’).
Fig. 3 shows how the delay changes as the maximum average transmit power of
each user varies. We can see that the delay of all schemes decreases with the
maximum average transmit power of each user. This is because a large maximum
average transmit power can decrease the transmission time between users and
the BS. We can clearly see that the proposed FL scheme achieves the best
performance among all schemes. This is because the proposed approach jointly
optimizes bandwidth and local accuracy $\eta$, while the bandwidth is fixed in
EB-FDMA and $\eta$ is not optimized in FE-FDMA. Compared to TDMA, the proposed
approach can reduce the delay by up to 27.3%.

##  5 Conclusions

In this paper, we have investigated the delay minimization problem of FL over
wireless communication networks. The tradeoff between computation delay and
transmission delay is determined by the learning accuracy. To solve this
problem, we first proved that the total delay is a convex function of the
learning accuracy. Then, we have obtained the optimal solution by using the
bisection method. Simulation results show the various properties of the
proposed solution.

## References

  * Ahn et al. (2019) Ahn, J.-H., Simeone, O., and Kang, J.  Wireless federated distillation for distributed edge learning with heterogeneous data.  _arXiv preprint arXiv:1907.02745_ , 2019.
  * Buza (2014) Buza, K.  Feedback prediction for blogs.  In _Data analysis, machine learning and knowledge discovery_ , pp. 145–152. Springer, 2014.
  * Chen et al. (2019a) Chen, M., Challita, U., Saad, W., Yin, C., and Debbah, M.  Artificial neural networks-based machine learning for wireless networks: A tutorial.  _IEEE Commun. Surveys Tut._ , pp. 1–1, 2019a.  ISSN 1553-877X.  doi: 10.1109/COMST.2019.2926625.
  * Chen et al. (2019b) Chen, M., Yang, Z., Saad, W., Yin, C., Poor, H. V., and Cui, S.  Performance optimization of federated learning over wireless networks.  In _Proc. IEEE Global Commun. Conf._ , pp. 1–6, Waikoloa, HI, USA, Dec. 2019b.
  * Chen et al. (2019) Chen, M., Yang, Z., Saad, W., Yin, C., Poor, H. V., and Cui, S.  A joint learning and communications framework for federated learning over wireless networks.  _arXiv preprint arXiv:1909.07972_ , 2019.
  * Chen et al. (2020) Chen, M., Semiari, O., Saad, W., Liu, X., and Yin, C.  Federated echo state learning for minimizing breaks in presence in wireless virtual reality networks.  _IEEE Trans. Wireless Commun._ , to appear, 2020.
  * Dong et al. (2019) Dong, P., Zhang, H., Li, G. Y., Gaspar, I. S., and NaderiAlizadeh, N.  Deep CNN-based channel estimation for mmWave massive MIMO systems.  _IEEE J. Sel. Topics Signal Process._ , 13(5):989–1000, Sept. 2019.
  * Gao et al. (2020) Gao, S., Dong, P., Pan, Z., and Li, G. Y.  Reinforcement learning based cooperative coded caching under dynamic popularities in ultra-dense networks.  _IEEE Trans. Veh. Technol._ , 69(5):5442–5456, 2020.
  * Gündüz et al. (2019) Gündüz, D., de Kerret, P., Sidiropoulos, N. D., Gesbert, D., Murthy, C. R., and van der Schaar, M.  Machine learning in the air.  _IEEE J. Sel. Areas Commun._ , 37(10):2184–2199, Oct. 2019.
  * Huang et al. (2020) Huang, C., Mo, R., and Yuen, C.  Reconfigurable intelligent surface assisted multiuser MISO systems exploiting deep reinforcement learning.  _IEEE J. Sel. Areas Commun._ , pp. 1–1, 2020.
  * Konečnỳ et al. (2016) Konečnỳ, J., McMahan, H. B., Ramage, D., and Richtárik, P.  Federated optimization: Distributed machine learning for on-device intelligence.  _arXiv preprint arXiv:1610.02527_ , 2016.
  * McMahan et al. (2016) McMahan, H. B., Moore, E., Ramage, D., Hampson, S., and Arcas, B. A. y.  Communication-efficient learning of deep networks from decentralized data.  _arXiv preprint arXiv:1602.05629_ , 2016.
  * Park et al. (2019) Park, J., Samarakoon, S., Bennis, M., and Debbah, M.  Wireless network intelligence at the edge.  _Proceedings of the IEEE_ , 107(11):2204–2239, Nov. 2019.
  * Saad et al. (to appear, 2020) Saad, W., Bennis, M., and Chen, M.  A vision of 6G wireless systems: Applications, trends, technologies, and open research problems.  _IEEE Network_ , to appear, 2020.
  * Samarakoon et al. (2018) Samarakoon, S., Bennis, M., Saad, W., and Debbah, M.  Distributed federated learning for ultra-reliable low-latency vehicular communications.  _arXiv preprint arXiv:1807.08127_ , 2018.
  * Tran et al. (2019) Tran, N. H., Bao, W., Zomaya, A., and Hong, C. S.  Federated learning over wireless networks: Optimization model design and analysis.  In _Proc. IEEE Conf. Computer Commun._ , pp. 1387–1395, Paris, France, June 2019.
  * Wang et al. (2018) Wang, S., Tuor, T., Salonidis, T., Leung, K. K., Makaya, C., He, T., and Chan, K.  When edge meets learning: Adaptive control for resource-constrained distributed machine learning.  In _IEEE Conf. Computer Commun._ , pp. 63–71, Honolulu, HI, USA, Apr. 2018.
  * Wang et al. (2019) Wang, S., Tuor, T., Salonidis, T., Leung, K. K., Makaya, C., He, T., and Chan, K.  Adaptive federated learning in resource constrained edge computing systems.  _IEEE J. Sel. Areas Commun._ , 37(6):1205–1221, June 2019.
  * Yang et al. (2020) Yang, H. H., Liu, Z., Quek, T. Q. S., and Poor, H. V.  Scheduling policies for federated learning in wireless networks.  _IEEE Trans. Commun._ , to appear, 2020.
  * Yang et al. (2018) Yang, K., Jiang, T., Shi, Y., and Ding, Z.  Federated learning via over-the-air computation.  _arXiv preprint arXiv:1812.11750_ , 2018.
  * Yang et al. (2019) Yang, Z., Chen, M., Saad, W., Hong, C. S., and Shikh-Bahaei, M.  Energy efficient federated learning over wireless communication networks.  _arXiv preprint arXiv:1911.02417_ , 2019.
  * Yang et al. (2020) Yang, Z., Chen, M., Saad, W., Xu, W., Shikh-Bahaei, M., Poor, H. V., and Cui, S.  Energy-efficient wireless communications with distributed reconfigurable intelligent surfaces, 2020.
  * Zeng et al. (2019) Zeng, Q., Du, Y., Leung, K. K., and Huang, K.  Energy-efficient radio resource allocation for federated edge learning.  _arXiv preprint arXiv:1907.06040_ , 2019.
  * Zhu et al. (2018a) Zhu, G., Liu, D., Du, Y., You, C., Zhang, J., and Huang, K.  Towards an intelligent edge: Wireless communication meets machine learning.  _arXiv preprint arXiv:1809.00343_ , 2018a.
  * Zhu et al. (2018b) Zhu, G., Wang, Y., and Huang, K.  Low-latency broadband analog aggregation for federated edge learning.  _arXiv preprint arXiv:1812.11494_ , 2018b.