# Influence of Micro-turbulence on Neoclassical Tearing Mode Onset

Tonghui Shi,1 L. Wei,2 H.H. Wang,1 E. Li,1 B. Shen,1 J.P. Qian,1 Y.M. Wang,1
T. Zhang,1 H.L. Zhao,1 L. Zeng,1 Y. Zhang,1 H.Q. Liu,1 Q. Ma,1 D.L. Chen,1
Z.P. Luo,1 Y.Y. Li,1 Z.C. Shen,1 L.Q. Xu,1 B. Zhang,1 M.H. Li,1 Z.X. Wang,2
B.L. Ling,1 X.Z. Gong,1 Y. Sun1,∗ and B. Wan1,† 1Institute of Plasma Physics,
Chinese Academy of Sciences, Hefei 230031, China
2School of Physics, Dalian University of Technology, Dalian 116024, China

###### Abstract

Direct evidence of micro-turbulence effect on the onset of neoclassical
tearing mode (NTM) is reported for the first time in this letter. A puzzling
positive correlation between critical width of seed island of NTM and
normalized plasma pressure $\beta_{p}$ is first observed employing a novel
method for clearly separating the processes of seed island and the onset of
NTM in the EAST tokamak. Different from the methods developed before, the
width of the seed island is well controlled by slowly ramping up the current
in resonant magnetic perturbation coils. It is revealed that the positive
correlation is mainly attributed to the enhancement of perpendicular transport
by micro-turbulence, which overcomes the destabilizing effect of $\beta_{p}$
on the onset of NTM. Reduced magnetohydrodynamics (MHD) modeling well
reproduced the two states of nonlinear bifurcations observed in this
experiment by including the finite transport effect. This result provides a
new route for understanding multi-scale interaction in plasma physics.

Magnetic island formed by tearing mode instability is a focus issue in space
and laboratory plasmas, which provides an excellent medium to study
bifurcation phenomena and multi-scale interaction between micro-turbulence and
macro-instability in magnetized plasmas Magnetic recnection RMP 2010 ; Ida.K
2020 . In tokamak plasmas, Neo-classical Tearing Mode (NTM), namely nonlinear
neoclassical pressure driven tearing mode, is one of the significant threats
to a fusion reactor TFTR 1995 PRL . It has been observed that small-scale
micro-turbulence can be influenced by large-scale magnetohydrodynamics (MHD)
instabilities Bardoczi PRL 2016 ; Sunpj NF 2018 , and in turn its enhanced
cross-field transport can affect the stability of large magnetic island
Bardoczi pop 2017a . As a nonlinear MHD instability, NTM requires a large
enough seed island to be triggered LaHaye 2006 . Plasma theory predicted
small-scale micro-turbulence can affect the threshold width of seed island for
triggering NTM Sanae PRL 2003 ; CaiHS NF 2019 . In this letter, we report for
the first time a direct evidence of the influence of micro-turbulence on the
onset of NTM by a new method to uncover the entire growth process of NTM, i.e.
separating the seed island phase and the nonlinear triggering of NTM phase by
well controlling the ramping rate of the seed island width using resonant
magnetic perturbations (RMPs), and present a distinct dependency between
critical width of seed island and normalized plasma pressure $\beta_{p}$.

Many achievements have been reached in NTM seed island research. In some
experiments, NTM was found to be triggered by some MHD instabilities LaHaye
2006 ; NTM-HL2A ; NTM-EAST ; NTM-KSTAR , in which its triggering process was
very fast-about tens of milliseconds. In other experiments, using classical
tearing mode and external magnetic perturbation as seed island and seed island
driving force respectively TCV TM to NTM 2002 ; AUG 2019 , the seed island
phase was demonstrated to sustain nearly $100ms$, and the bifurcation related
to NTM triggering threshold was distinguished. Besides, marginal width of NTM
magnetic island, indirectly reflecting the characteristics of NTM seed island,
was obtained by ramping-down input power, presenting a proportional relation
to ion banana orbit width LaHaye NF 2006 w_mar and rho_i , which could be well
explained by NTM theory including finite orbit width effect Poli PRL 2002 ;
Imada PRL 2018 ; Wang Feng 2019 . Although NTM theory, which contains a
variety of effects Sauter 1997 , has successfully explained many observations
LaHaye 2006 , uncertainties still remain in the physics describing small
magnetic island, e.g. NTM seed island, because the growth process of NTM is
difficult to be controlled and the plasma is unstable before NTMs saturated.

Figure 1: (color online). Temporal evolutions of (a) $n=1$ perturbed magnetic
field $\delta{B}_{r,n=1}$ measured by saddle loops, (b) Mirnov signal
$\dot{\delta{B}_{\theta}}$, (c) and (d) spectra of Mirnov signals for toroidal
mode number $n$ and RMP currents for shots $91970$ (red dashed lines) and
$91972$ (black solid lines), respectively.

In this letter, we introduce a new method that the width of seed island is
well controlled by slowly ramping up the current in RMP coils, seed island
phase sustains hundreds of milliseconds, and the plasma is more stable before
NTM triggering. EAST tokamak EAST 2019  has a flexible RMPs coil system
consisting of $(2\times 8)$ coils installed in the low side EAST RMP PRL 2016
; EAST RMP NF 2017 . NTMs are often stable in EAST high confinement discharges
due to low intrinsic error field EAST error field NF 2016 . All shots shown in
this letter have been set as even coil connection between the upper and lower
RMP coils, which can easily generate mode penetration in EAST EAST error field
NF 2016 .

Fig. 1 shows an example of the triggering of $m/n=2/1$ NTM by $n=1$ RMP in two
shots $\sharp 91970$ and $\sharp 91972$, in which mode penetrations ($1th$-BF)
are both observed when the RMP currents exceed some critical values, where
magnetic island amplitude $\delta{B}_{r,n=1}$ (Fig. 1(a)) and poloidal and
toroidal mode numbers $(m,n)$ are measured by $(5\times 8)$ saddle coils
installed in different toroidal and poloidal angles. The two shots have nearly
same plasma parameters except the RMP current, i.e. magnetic safety factor
$q_{95}\simeq 5.4$, $\beta_{N}\simeq 0.72$, $\langle{N}_{e0}\rangle\simeq
2.6(\times 10^{19}/m^{3})$ measured by POlarimeter-INTerferometer (POINT)
system EAST POINT 2014 , $T_{i0}\simeq 1.0kev$ measured by charge exchange
recombination spectroscopy (CXRS) EAST CXRS 2014 , lower hybrid wave heating
(LHW) EAST LHCD NF 2015  $P_{LHW}\simeq 1.0MW$, and co-current neutral beam
injection (NBI) EAST NBI 2015  $P_{NBI}=2.0MW$. In shot $\sharp 91972$ (Fig. 1
black solid line), a bifurcation ($2nd$-BF) related to NTM triggering
threshold ($t=5.9s$) is observed when RMP current grows to a critical value.
Besides, the magnetic island unlocks when RMP current is turned off as shown
in Mirnov signal (Fig. 1(b) black solid line) and its spectrum (Fig. 1(c)).
For comparison, in shot $\sharp 91970$ (Fig. 1 red dashed line), if magnetic
island width induced by RMP doesn’t exceed some critical value, NTM can not be
triggered and magnetic island decays immediately as RMP current is turned off
(Fig. 1(a) and 1(b) red dashed lines and 1(d) spectrum). Therefore, by slowly
ramping up the current in RMP coils, seed island width is well controlled, and
the growth of seed island and the onset of NTM are clearly separated in this
experiment.

Figure 2: (color online). Temporal evolutions of (a) $2/1$ magnetic island
width $\mathcal{W}$ from experiment (black doted line), the simulation results
of mode-1 (blue line) and model-2 (red line) by using eq. (1), (b) toroidal
phase of $n=1$ perturbations in plasma response (black doted line) and vacuum
field (red solid line), respectively, (c) normalized variations of line
averaged density (doted line) and perturbed electron temperature (solid line)
near the $q=2$ surface.

Applying this method, the dynamics of NTM triggering are clearly displayed. As
shown in Fig. 2, the rapid growth of seed island width (Fig. 2(a) black dotted
line) during $t_{a}<t<t_{b}$ is resulted from the fast-change of $\phi_{re}$
(Fig. 2(b) black dotted line), where plasma response phase $\phi_{res}$ and
vacuum field phase $\phi_{vac}$ are measured by saddle coils, island width
$\mathcal{W}\simeq{a}_{eq}\sqrt{\delta{B}_{r,n=1}}$, and $a_{eq}$ is related
to plasma equilibrium parameters island formula PRL 2002 . The width of seed
island during $t_{b}<t<t_{c}$ grows slowly and is well controlled by slowly
ramping up the current in RMP coils. After phase transition (Fig. 2(b),
$t>t_{b}$), the plasma response phase is nearly locked on vacuum field, even
when $2nd$-BF ($t=t_{c}$) happens. The normalized perturbations near the $q=2$
surface, i.e. electron temperature $\delta{T_{e}}/{T_{e}}$ (Fig. 2(c) colored
solid line) and electron density
$\delta{\langle{N_{e}}\rangle}/\langle{N_{e}}\rangle$ (Fig. 2(c) black doted
line), decline near $1th$-BF $(t=t_{a})$ and $2nd$-BF $(t=t_{c})$, where
$T_{e}$ and $\langle{N_{e}}\rangle$ are measured by electron cyclotron
emission (ECE) system EAST ECE 2018  and POINT respectively. But the declines
near $2nd$-BF are sharper than that near $1th$-BF. This property of NTM
triggering threshold has been disclosed in NTM theory and simulations Poli PRL
2002 ; Imada PRL 2018 ; Wang Feng 2019 . Using Eq. (1) with the parameters in
table I, Model-2 (Fig. 2(a) red solid line) reconciles the experimental result
(Fig. 2(a) black doted line), nevertheless, in Model-1 (Fig. 2(a) blue solid
line), only one bifurcation is observed. In essence, the neoclassical driving
force $\Delta_{bs}$ (Eq. (2)) becomes significant only when
$\mathcal{W}>\mathcal{W}_{cri}$ in Model-2, which is consistent with
experimental observations as shown in Fig. 2(c). Therefore, the dynamics of
NTM triggering is clearly presented by this method, and the growth process of
NTM is well simulated by MRE model as follows.

A Modified Rutherford Equation (MRE) Sauter 1997 ; Fitzpatrick NF 1993  with a
switch $f_{bs}(\mathcal{W})$ about neoclassical effect (bootstrap current) is
introduced as

$\frac{\tau_{R_{nc}}}{r_{s}}\frac{d\mathcal{W}}{dt}=r_{s}\Delta^{\prime}+\Delta_{rmp}+\Delta_{bs},$
(1)

with

$\begin{split}\left\\{\begin{aligned}
&\Delta_{rmp}=A_{rmp}\frac{\mathcal{W}_{vac}^{2}}{\mathcal{W}^{2}}\cos(\Delta\phi(t)),\\\
&\Delta_{bs}=A_{bs}\beta_{p}f_{bs}\frac{{r_{s}}\mathcal{W}}{\mathcal{W}^{2}+\mathcal{W}_{cri}^{2}},\end{aligned}\right.\end{split}$
(2)

where $f_{bs}(\mathcal{W})={1}/[1+(\mathcal{W}_{cri}/\mathcal{W})^{n_{bs}}]$,
$\mathcal{A}_{rmp}$ and $\mathcal{A}_{bs}$ are related to plasma equilibrium
parameters, $\tau_{R_{nc}}$ is the effective resistive time Sauter 1999 ,
vacuum magnetic island width $\mathcal{W}_{vac}$ calculated by MAPS code Sun
MAPS code , and $\Delta\phi(t)$ is phase difference. If setting $n_{bs}=0$,
Eq. (1) becomes the well-known model (Model-1) Fitzpatrick finite transport
1995 . And if setting $n_{bs}\gg 1$, Eq. (1) is Model-2, where
$f(\mathcal{W})\simeq 1$ when $\mathcal{W}>\mathcal{W}_{cri}$, and
$f(\mathcal{W})\simeq 0$ when $\mathcal{W}<\mathcal{W}_{cri}$.

Table 1: Modeling parameters in model-1 and model-2 $\tau_{R_{nc}}(s)$ |  $r_{s}\Delta^{\prime}$ |  $A_{rmp}$ |  $A_{bs}\beta_{p}$ |  $\mathcal{W}_{cri}(m)$ |  $r_{s}(m)$
---|---|---|---|---|---
0.3 |  -2.61 |  6.78 |  0.24 |  0.03 |  0.23

Figure 3: Dependencies of threshold RMP current $Ir_{cri}$ (a) and seed island
critical width $\mathcal{W}_{cri}$ (b) on $\beta_{p,cri}$ at the NTM onset
time.

Based on the statistical results, a puzzling positive correlation between seed
island critical width and $\beta_{p}$ is discovered for the first time. As
illustrated in Fig. 3, both RMP critical current value $Ir_{cri}$ and seed
island critical width $\mathcal{W}_{cri}$ increase with $\beta_{p,cri}$, where
the plasma parameters are similar except for $q_{95}\subset[4.5,5.0,5.4]$,
i.e. $\langle{N_{e0}}\rangle=2.57\pm 0.12(\times 10^{19}/m^{3})$,
$T_{i0}=1.08\pm 0.14kev$. As is well known, the fraction of bootstrap current
is nearly proportional to plasma poloidal beta $\beta_{p}$. Then, dose it mean
the bootstrap current has a stabilizing effect on NTM onset?

In fact, from NTM theory Sauter 1997 , neoclassical effect (bootstrap current)
always plays a destabilizing role in NTM triggering. From a recent simulation
Q Yu NF 2020 , bootstrap current shows a stabilizing effect on NTM onset
because of particle pumping-in effect caused by RMP when magnetic island
rotates in co-current direction before mode penetration. Nevertheless, in our
case, seed island has been locked on vacuum field and island rotation remains
nearly zero for a long time before $2nd$-BF happening. Moreover, curvature
effect is very small in conventional high aspect ratio tokamaks and can’t
change the property of bootstrap current destabilization on NTM onset Sauter
1997 . And, we don’t discuss neo-classical polarization current effect or
finite orbit width effect Poli PRL 2002 ; Imada PRL 2018 ; Wang Feng 2019
because these shots have nearly same ion temperatures as shown in Fig. 3. In
addition, in other NTM experiments NTM-COMPASSD ; NTM-DIIID ; NTM-EAST , none
has observed an obvious dependency between seed island threshold width and the
friction of bootstrap current (or $\beta_{p}$), although bootstrap current is
the main driving force of NTM.

Figure 4: (color online). Local perturbations near the $q=2$ surface. Cross-
powers spectra of (a) electron density $\widetilde{N}_{e1}\widetilde{N}_{e2}$
and (b) electron temperature $\widetilde{T}_{e1}\widetilde{T}_{e2}$. (c)
Dependencies of $\langle\widetilde{N}_{e1}\widetilde{N}_{e2}\rangle$ and
$\langle\widetilde{T}_{e1}\widetilde{T}_{e2}\rangle$ on $\beta_{p}$, where
$\langle\widetilde{f}_{1}\widetilde{f}_{2}\rangle_{f=T_{e},N_{e}}$ presents
the sum of cross-powers $\widetilde{f}_{1}\widetilde{f}_{2}$ in the
$0.1MHz-0.5MHz$ frequency range.

From the measurements of turbulence powers, it is found that the cross-field
transport enhanced by micro-turbulence can affect the onset of NTM. As shown
in Fig. 4 (a) and (b), the amplitudes of cross-powers of electron density
$\widetilde{N}_{e1}\widetilde{N}_{e2}$ and electron temperature
$\widetilde{T}_{e1}\widetilde{T}_{e2}$ increase with $\beta_{p}$ in full
$0.1MHz-0.5MHz$ frequency range, where fluctuations of $\widetilde{N}_{e1,2}$
and $\widetilde{T}_{e1,2}$ are measured by reflectometer EAST refectometry
and ECE respectively. These three similar shots have nearly same plasma
equilibrium parameters except $\beta_{p}$ caused by different LHW heating
powers. And the fluctuated amplitudes
$\langle\widetilde{T}_{e1}\widetilde{T}_{e2}\rangle$ and
$\langle\widetilde{N}_{e1}\widetilde{N}_{e2}\rangle$ have positive dependency
on $\beta_{p}$ (Fig. 4(c)), which are associated with broadband low-$k$
turbulence powers and calculated by nearest neighbor cross-powers just before
$2nd$-BF happening in the $0.1MHz-0.5MHz$ frequency range. Besides, in these
shots (Figs. 3 and 4), the increase of $\beta_{p}$ mainly arises from the
increase of electron temperature caused by LHW heating. Consequently, the
enhancement of turbulence power may be related to the gradient of electron
temperature, and thus the cross-field transport enhanced by micro-turbulence
plays a stabilizing role in the onset of NTM.

Figure 5: (color online). Temporal evolutions of (a) normalized magnetic
island width, (b) input normalized RMP amplitude. (c) NTM marginally stable
boundary: $\beta_{p,mar}\propto\mathcal{W}_{d}/r_{a}$.

These observations are well explained by NTM theory with the finite transport
effect. Using a reduced MHD code including the finite transport effect L Wei
NF 2016 , the two states of nonlinear bifurcations in this experiment (Fig. 1)
are well reproduced as shown in Fig. 5 (a) and (b). When the plasma is close
to marginally stable state of NTM (i.e. the fraction of bootstrap current is
not large enough), this positive correlation between seed island critical
width and $\beta_{p}$ is also simulated as illustrated in Fig. 5(c), in which
the stable and unstable regions of NTM are clearly separated by the marginal
stable boundary of NTM, i.e. $\beta_{p,mar}\propto{\mathcal{W}_{d}}/r_{a}$,
where the scale length of finite transport is calculated by a formula
$\mathcal{W}_{d}=[5.1r_{s}/(\epsilon_{s}ns_{s})^{1/2}](\chi_{\perp}/\chi_{\parallel})^{1/4}$
Fitzpatrick finite transport 1995 , $\epsilon_{s}=r_{s}/R$ and local magnetic
shear $s_{s}=rq^{\prime}/q|_{r=r_{s}}$. And, $\chi_{\perp}/\chi_{\parallel}$
and $\beta_{p,mar}$ are the output of this reduced MHD code. In addition,
applying MRE model Eq. (1), we also obtain a proportional dependency, i.e.
$\beta_{p,cri}\sim\beta_{p,mar}=-(r_{s}\Delta^{\prime}/\mathcal{A}_{bs})(\mathcal{W}_{cri}/r_{s})$
NTM-COMPASSD , where $\mathcal{W}_{cri}\sim\mathcal{W}_{d}$ related to
perpendicular transport strength enhanced by micro-turbulence. Therefore, when
plasma is close to the marginal state of NTM, micro-turbulence enhanced cross-
field transport overcomes the destabilization of bootstrap current and shows a
distinct stabilizing effect on NTM onset.

In summary, direct evidence of the influence of micro-turbulence on the onset
of NTM is reported for the first time in the EAST tokamak. Applying a novel
method, a positive dependency between critical width of seed island of NTM and
plasma poloidal beta is first discovered. Different from the methods developed
before, the process of seed island and the onset of NTM is distinguished more
clearly, and the width of the seed island is well controlled by slowly ramping
up the current in RMPs. From the analysis of fluctuated amplitudes associated
with turbulence power, it is revealed that the positive dependency is mainly
resulted from the enhancement of cross-field transport by micro-turbulence,
which overcomes the destabilizing effect of $\beta_{p}$ and shows a
stabilizing effect on the onset of NTM. Finally, using a reduced MHD model
including the finite transport effect, the two bifurcations related to mode
penetration and NTM triggering in this experiment are well reproduced, and the
stabilizing effect of perpendicular transport enhanced by micro-turbulence on
the NTM onset is also simulated. In brief, this method offers remarkable
potentials for studying magnetic island dynamics and multi-scale interaction
between micro-turbulence and small-scale magnetic island.

One of authors, Tonghui Shi, would like to thank Q. Yu (Max-Planck-Institut
für Plasmaphysik, Germany) for helpful discussion. This work is supported by
the National Key R$\&$D Program of China (Grants No: 2017YFE0301100,
2018YFE0302100, 2019YFE03010002), and the National Natural Science Foundation
of China (Grants No: 11905250 and 11875292).

<EMAIL_ADDRESS><EMAIL_ADDRESS>
## References

  * (1) M. Yamada, R. Kulsrud, and H. Ji, Rev. Mod. Phys. 82, 603 (2010).
  * (2) K. Ida, ADVANCES IN PHYSICS:X 5, 1801354 (2020).
  * (3) Z. Chang, J.D. Callen, E.D. Fredrickson, R.V. Budny, C.C. Hegna, K.M. McGuire, M.C. Zarnstorff, and TFTRgroup, Phys. Rev. Lett. 74, 4663 (1995).
  * (4) L. Bardóczi, T.L. Rhodes, T.A. Carter, A. Bañón Navarro, W.A. Peebles, F. Jenko and G. McKee, Phys. Rev. Lett. 116, 215001 (2016).
  * (5) P.J. Sun et al, Nucl. Fusion 58, 016003 (2018).
  * (6) L. Bardóczi, T.A. Carter, R.J. La Haye, T.L. Rhodes, and G.R. McKee, Phys. Plasmas 24, 122503 (2017).
  * (7) R.J. La Haye, Phys. Plasmas 13, 055501 (2006), and references therein.
  * (8) Sanae-I. Itoh, Kimitaka Itoh, and Masatoshi Yagi, Phys. Rev. Lett. 91, 045003 (2003).
  * (9) Huishan Cai, Nucl. Fusion 59, 026009 (2019).
  * (10) X.Q. Ji, Q.W. Yang, Y. Liu, J. Zhou, B.B. Feng, B.S. Yuan, Chin. Phys. Lett. 27, 065202 (2010).
  * (11) T.H. Shi et al, Plasma Phys. Control. Fusion 55, 055007 (2013).
  * (12) M. Joung et al, AIP Conference Proceedings 1580, 506 (2014).
  * (13) H. Reimerdes, O. Sauter, T. Goodman, and A. Pochelon, Phys. Rev. Lett. 88, 105005 (2002).
  * (14) Dmytro Meshcheriakov et al, Phys. Plasmas 26, 042504 (2019).
  * (15) R.J. La Haye, R. Prater, R.J. Buttery, N. Hayashi2, A. Isayama, M.E. Maraschek, L. Urso and H. Zohm, Nucl. Fusion 46, 451-461 (2006).
  * (16) E. Poli, A.G. Peeters, A. Bergmann, S. Günter, and S.D. Pinches, Phys. Rev. Lett. 88, 075001 (2002).
  * (17) K. Imada, H.R. Wilson, J.W. Connor, A.V. Dudkovskaia, and P. Hill, Phys. Rev. Lett. 121, 175001 (2018).
  * (18) F. Wang, J.Q. Li, H.P. Qu, X.D. Peng, and Y. Xiao, Phys. Plasmas 26, 052516 (2019).
  * (19) O. Sauter et al, Phys. Plasmas 4, 1654 (1997), and references therein.
  * (20) B. Wan et al, Nucl. Fusion 59, 112003 (2019).
  * (21) Y. Sun, Y. Liang, Y.Q. Liu, S. Gu, X. Yang, W. Guo et al, Phys. Rev. Lett. 117, 115001 (2016).
  * (22) Y. Sun et al, Nucl. Fusion 57, 036007 (2017).
  * (23) H.H. Wang et al, Nucl. Fusion 56, 066011 (2016).
  * (24) H.Q. Liu et al, Rev. Sci. Instrum. 85, 11D405 (2014).
  * (25) Y.Y. Li et al, Rev. Sci. Instrum. 85, 11E428 (2014).
  * (26) H. Zhao, T.F. Zhou, Y. Liu, A. Ti, B.L. Ling, M.E. Austin, Saeid Houshmandyar, H. Huang, W.L. Rowan, and L.Q. Hu, Rev. Sci. Instrum. 89, 10H111 (2018).
  * (27) F.K. Liu et al, Nucl. Fusion 55, 123022 (2015).
  * (28) C.D. Hu et al, Plasma Sci. Technol. 17, 817 (2015).
  * (29) F. Salzedas, F.C. Schüller, A.A.M. Oomens, and the RTP Team, Phys. Rev. Lett. 88, 075002 (2002).
  * (30) R. Fitzpatrick, Nucl. Fusion 33, 1049-1084 (1993).
  * (31) O. Sauter, C. Angioni and Y.R. Lin-Liu, Phys. Plasmas 6, 2834 (1999).
  * (32) Y. Sun, Y. Liang, J. Qian, B. Shen and B. Wan, Plasma Phys. Control. Fusion 57, 045003 (2015).
  * (33) R. Fitzpatrick, Phys. Plasmas 2, 825 (1994).
  * (34) L. Wei, Z.X. Wang, J.L. Wang and X.F. Yang, Nucl. Fusion 56, 106015 (2016).
  * (35) Q. Yu, Nucl. Fusion 60, 084001 (2020).
  * (36) D.A. Gates et al, Nucl. Fusion 37, 1593 (1997).
  * (37) R.J. La Haye, O. Sauter, Nucl. Fusion 38, 987 (1998).
  * (38) H.M. Xiang et al, Rev. Sci. Instrum. 89, 10H103 (2018).


# A Study of Multihop mmW Aerial Backhaul Links

Mohammad T. Dabiri1, Mazen O. Hasna1, Tamer Khattab1, and Khalid Qaraqe2
This publication was made possible by NPRP13S-0130-200200 from the Qatar
National Research Fund (a member of The Qatar Foundation). The statements made
herein are solely the responsibility of the author[s]. 1Department of
Electrical Engineering,Qatar University, Doha, Qatar,
E-mails<EMAIL_ADDRESS> 2Department of Electrical and
Computer Engineering,Texas A&M University at Qatar, Doha, Qatar, E-mail:
<EMAIL_ADDRESS>
###### Abstract

The main contribution of this paper is to analyze a long networked flying
platform (NFP)-based millimeter wave (mmWave) backhaul link that is offered as
a cost effective and easy to deploy solution to connect a disaster or remote
area to the nearest core network. For this aim, we characterize the backhaul
channel as a function of realistic physical parameters such as heights and
distances of obstacles along the route, flight altitude and the intensity of
NFPs’ vibrations, the real 3D antenna pattern provided by 3GPP, etc. For the
characterized channel, we derive an analytical closed-form expression for the
outage probability. Finally, using the obtained results, we provide a fast
algorithm for the optimal parameter design of the considered system that
minimizes the cost.

###### Index Terms:

Antenna pattern, backhaul/frounthaul links, positioning, mmWave communication,
unmanned aerial vehicles (UAVs).

##  I Introduction

Tornadoes, hurricanes, earthquakes, tsunami and other major natural disasters
all have the potential to cut or entirely destroy fibre infrastructure to the
disaster area. Any disruption to the fragile fibre causes data outages which
take days or weeks to locate and repair. However, providing wireless
alternative connectivity in the immediate moments after a disaster event is a
key feature to facilitate rescue operations. Implementing a long terrestrial
wireless backhaul link faces major challenges, including providing a line of
sight (LoS) between the nearest core network to the disaster area, especially
for mountainous or forested areas. Due to their unique capabilities such as
flexibility, maneuverability, and adaptive altitude adjustment, unmanned
aerial vehicles (UAVs) acting as networked flying platforms (NFPs) can be
considered as a promising solution to provide a temporary wireless backhaul
connectivity while improving flexibility and reliability of backhaul
operations [1, 2]. More recently, millimeter wave (mmWave) backhauling has
been proposed as a promising approach for aerial communications because of
three reasons [3]. First, unlike terrestrial mmWave communication links that
suffer from blockage, the flying nature of UAVs offers a higher probability of
LoS. Second, the large available bandwidth at mmWave frequencies can provide
high data rate. Third, to compensate the negative effects of the high path-
loss at the mmWave bands, the small wavelength enables the realization of a
compact form of highly directive antenna arrays which is suitable for small
UAVs with limited payload. Although NFP-based mmWave backhaul link has been
studied in recent works [3, 4, 5, 6, 7, 8, 9, 10, 11], the results of these
studies are limited for short-distance NFP-based backhaul communications.

Figure 1: An illustration of a long aerial backhaul link in order to transfer
data from the nearest core network to the disaster or remote area by using
four NFPs.

In this study, we consider a long NFP-assisted backhaul link as shown in Fig.
1 that is offered as a cost effective and easy to deploy solution to connect a
disaster or remote area to the nearest core network in a short time. 111The
results of this study are provided for the types of NFPs that hover in the sky
such as multi-rotor UAVs or tethered balloons and thus, for the optimal design
of systems based on fixed-wing UAVs that are moving at high speed, a separate
study should be done, which is beyond the scope of this work. In particular,
we characterize a more practical scenario by tacking into account the effects
of realistic physical parameters such as heights and distances of obstacles
along the route, flight altitude and the severity of UAVs’ vibrations, the
real 3D antenna pattern provided by 3GPP, atmospheric channel loss as a
function of temperature and air pressure, etc. In this work, the QoS is
evaluated using the outage probability metric, which is a very important
criterion in wireless communication. For this aim, we characterize the
backhaul channel as a function of all aforementioned channel parameters and we
derive an analytical closed-form expression for the outage probability. By
using Monte Carlo simulations, the accuracy of the derived analytical
expression is verified. Then, we investigate the effects of key channel
parameters such as antenna pattern gain, strength of UAV’s vibrations, UAVs’
positions in the sky, frequency bands, and physical parameters of a given
region on the performance of the considered long NFP-based mmWave backhaul
link in terms of outage probability. By providing sufficient simulation
results, we carefully study the relationships between these parameters in
order to reduce the negative effect of UAV vibrations by adjusting the 3D
antenna pattern and at the same time, decrease channel loss by adjusting
optimal positions of UAVs in the sky compared to the obstacles. Finally, using
the obtained results, we provide a fast algorithm for the optimal design of
system parameters that minimizes the cost.

##  II The System Model

As illustrated in Fig. 1, we consider a multi-hop NFP-based mmWave backhaul
link to transfer data from the nearest core network (source) to the disaster
area (destination) where the terrestrial infrastructure is damaged due to
natural events. Let $M$ be the number of NFPs where $U_{1}$ and $U_{M}$ denote
the first and last NFPs, respectively, and $U_{i}$ for $i\in\\{2,...,M-1\\}$
denotes intermediate relays as shown in Fig. 1. The core network sends signals
towards $U_{1}$ and after passing through the intermediate relays, finally by
the last relay, $U_{M}$, the signals will be directed to the destination in
the remote area.

Each NFP is equipped with two mmWave directional antennas, one for the
transmitter and one for the receiver. Let $A_{t,i}$ and $A_{r,i}$ denote
respectively the transmitter and receiver antennas mounted on $U_{i}$ for
$i\in\\{1,...,M\\}$. Also, let $A_{s}$ and $A_{d}$ be the core network and
destination antennas, respectively. Using location information obtained from
the GPS, $U_{1}$ adjusts the direction of $A_{r,1}$ and $A_{t,1}$ toward the
$A_{s}$ and $A_{r,2}$, respectively. Similarly, for $i\in\\{2,...,M-1\\}$,
$U_{i}$ adjusts the direction of $A_{r,i}$ and $A_{t,i}$ toward the
$A_{t,i-1}$ and $A_{r,i+1}$, respectively. $U_{M}$ adjusts the direction of
$A_{r,M}$ and $A_{t,M}$ toward the $A_{t,M-1}$ and $A_{d}$, respectively.
Assuming that each link is assigned a separate frequency band, the received
signal at the first NFP, $i$th NFPs for $i\in\\{2,...,M\\}$, and destination
are obtained respectively as

$\displaystyle\left\\{\begin{array}[]{rl}&\\!\\!\\!\\!\\!\\!\\!\\!P_{r,1}=P_{t,s}h_{L_{s}}G_{s}(\theta_{tx,s},\theta_{ty,s})G_{r,1}(\theta_{rx,1},\theta_{ry,1}),\\\
&\\!\\!\\!\\!\\!\\!\\!\\!P_{r,i}=P_{t,i-1}h_{L_{i-1}}G_{t,i-1}(\theta_{tx,i-1},\theta_{ty,i-1})\\\
&~{}~{}~{}~{}\times G_{r,i}(\theta_{rx,i},\theta_{ry,i}),\\\
&\\!\\!\\!\\!\\!\\!\\!\\!P_{r,d}=P_{t,M}h_{L_{M}}G_{t,M}(\theta_{tx,M},\theta_{ty,M})G_{d}(\theta_{rx,d},\theta_{ry,d}),\end{array}\right.$
(5)

where $P_{t,s}$ and $G_{s}(\theta_{tx,s},\theta_{ty,s})$ are respectively the
transmitted power and the antenna pattern gain of ground core network
transmitter, $P_{t,i}$ and $G_{t,i}(\theta_{tx,i},\theta_{ty,i})$ are
respectively the transmitted power and the antenna pattern gain of $A_{t,i}$,
$G_{r,i}(\theta_{rx,i},\theta_{ry,i})$ is the antenna pattern gain of
$A_{r,i}$, and $G_{d}(\theta_{rx,d},\theta_{ry,d})$ is antenna pattern gain of
the ground destination. The parameters $h_{L_{0}}$, $h_{L_{M}}$, and
$h_{L_{i}}$ for $i\in\\{1,...,M-1\\}$ respectively denote the channel loss of
source to $U_{1}$ link, $U_{M}$ to destination link, and $U_{i}$ to $U_{i+1}$
link. The UAVs’ vibrations are characterized by random variables
$\theta_{tx,i}$s, $\theta_{ty,i}$s, $\theta_{rx,i}$s, $\theta_{ry,i}$s.

As shown in Fig. 1, minimum elevation angels of the source is denoted by
$\psi_{s,\text{min}}$ which is a function of source height $H_{s}$, and height
of the nearest obstacle to the source $H_{bs}$. Also, minimum elevation angels
of the destination is denoted by $\psi_{d,\text{min}}$ which is a function of
destination height $H_{d}$, and height of the nearest obstacle to the
destination $H_{bd}$. The height difference between the source and the
destination is also denoted by $H_{sd}$. $L_{sd}$ is the horizon link length
between the core network and the remote area and $H_{b,\text{max}}$ is the
height of highest obstacle.

###  II-1 Channel Propagation Loss

In normal atmospheric conditions, water vapor (H2O) and oxygen (O2) molecules
are strongly absorptive of radio signals, especially at mmWave frequencies and
higher. Channel loss is usually expressed in dB and it can be calculated using
the formula

$\displaystyle h_{L,\text{dB}}^{\text{tot}}(f_{c})=20\log\left(\frac{4\pi
L}{\lambda}\right)+h_{L,\text{dB}}^{o,w}(f_{c}),$ (6)

where $L$ is the link length (in m), $\lambda$ is the wavelength (in m),
$f_{c}$ is mmWave frequency (in GHz),
$h_{L,\text{dB}}^{o,w}(f_{c})=\frac{h_{L,\text{dB/km}}^{o,w}(f_{c})L}{1000}$
is the attenuation due to oxygen and water (in dB),
$h_{L,\text{dB/km}}^{o,w}(f_{c})=h_{L,\text{dB/km}}^{o}(f_{c})+h_{L,\text{dB/km}}^{w}(f_{c})$
is the attenuation due to oxygen and water (in dB/km). At 20°C surface
temperature and at sea level, approximate expressions for the attenuation
constants of oxygen and water vapor (in dB/km) as defined by the International
Telecommunications Union (ITU) are [12]:

$\displaystyle h_{L,\text{dB/km}}^{o,0}(f_{c})=0.001\times f_{c}^{2}$ (7)
$\displaystyle\times\left\\{\begin{array}[]{rl}\frac{6.09}{f_{c}^{2}+0.227}+\frac{4.81}{(f_{c}-57)^{2}+1.5}{}{}{}{}{}{}{}{}{}{}{}&~{}~{}~{}f_{c}<57\\\
h_{L,\text{dB/km}}^{o,0}(f_{c}=57)+1.5(f_{c}-57)&~{}~{}~{}57<f_{c}<63\\\
\frac{4.13}{(f_{c}-63)^{2}+1.1}+\frac{0.19}{(f_{c}-118.7)^{2}+2}{}{}{}{}{}{}&~{}~{}~{}63<f_{c}<350\end{array}\right.$
(11)

and

$\displaystyle h_{L,\text{dB/km}}^{w,0}(f_{c})=0.0001\times
f_{c}^{2}\rho_{0}\left(0.05+\frac{3.6}{(f_{c}-22.2)^{2}+8.5}\right.$
$\displaystyle\left.+\frac{10.6}{(f_{c}-183.3)^{2}+9}+\frac{8.9}{(f_{c}-325.4)^{2}+26.3}\right),~{}~{}f_{c}<350,$
(12)

where $\rho_{0}=7.5~{}\text{g/m}^{3}$ is the water vapor density at sea level,
and $h_{L,\text{dB/km}}^{o,0}(f_{c}=57)$ is the value of the first expression
at $f_{c}=57$. In general, the attenuation constants of oxygen and water vapor
are functions of altitude, since they depend on factors such as temperature
and pressure. These quantities are often assumed to vary exponentially with
height $H$, as $\rho(H)=\rho_{0}\exp\left(-H/H_{\text{scale}}\right)$ where
$H_{\text{scale}}$ is known as the scale height, which is typically 1-2 km.
From this, the specific attenuation as a function of height can be
approximately modeled as

$\displaystyle
h_{L,\text{dB/km}}^{o,w}(f_{c},H)=h_{L,\text{dB/km}}^{o,w,0}(f_{c})\exp\left(-H/H_{\text{scale}}\right).$
(13)

For a slant atmospheric path from height $H_{1}$ to $H_{2}$ at an angle
$\psi$, the total atmospheric attenuation is determined by integration of the
specific attenuation as

$\displaystyle
h_{L,\text{dB/km}}^{o,w}(f_{c})\simeq\frac{h_{L,\text{dB/km}}^{o,w,0}(f_{c})\left(e^{-H_{1}/H_{\text{scale}}}-e^{-H_{2}/H_{\text{scale}}}\right)H_{s}}{\sin(\psi)}.$
(14)

In our system model, both core network to first NFP (CU) link and last NFP
node to destination (UD) link are slant and the inter-NFP links are
approximately horizontal.

###  II-2 3D Antenna Pattern

The array radiation gain is mainly formulated in the direction of $\theta$ and
$\phi$. In our model, $\theta$ and $\phi$ can be defined as functions of
random variables (RVs) $\theta_{x}$ and $\theta_{y}$ as follows:
$\theta=\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{x})+\tan^{2}(\theta_{y})}\right)$,
$\phi=\tan^{-1}\left({\tan(\theta_{y})}\big{/}{\tan(\theta_{x})}\right)$. By
taking into account the effect of all elements, the array radiation gain in
the direction of angles $\theta_{x}$ and $\theta_{y}$ will be:

$\displaystyle
G(\theta_{x},\theta_{y})=G_{0}(N)\,\underbrace{G_{e}(\theta_{x},\theta_{y})\,G_{a}(\theta_{x},\theta_{y})}_{G^{\prime}(\theta_{x},\theta_{y^{\prime}})},$
(15)

where $G_{a}$ is an array factor, $G_{e}$ is single element radiation pattern
and $G_{0}$ is a constant defined in the sequel. From the 3GPP single element
radiation pattern, $G_{e,\textrm{3dB}}=10\times\log_{10}(G_{e})$ of each
single antenna element is obtained as [13].

If the amplitude excitation of the entire array is uniform, then the array
factor $G_{a}(\theta_{x},\theta_{y})$ for a square array of $N\times N$
elements can be obtained as [14, eqs. (6.89) and (6.91)]

$\displaystyle G_{a}(\theta_{x},\theta_{y})=$ (16)
$\displaystyle\frac{\sin^{2}\left(\frac{N(kd_{x}\sin(\theta)\cos(\phi)+\beta_{x})}{2}\right)}{N^{2}\sin^{2}\left(\frac{kd_{x}\sin(\theta)\cos(\phi)+\beta_{x}}{2}\right)}\frac{\sin^{2}\left(\frac{N(kd_{y}\sin(\theta)\sin(\phi)+\beta_{y})}{2}\right)}{N^{2}\sin^{2}\left(\frac{kd_{y}\sin(\theta)\sin(\phi)+\beta_{y}}{2}\right)}$

where $\beta_{x}$ and $\beta_{y}$ are progressive phase shift between the
elements along the $x$ and $y$ axes, respectively. For a fair comparison
between antennas with different $N$, we assume that the total radiated power
of antennas with different $N$ are the same. From this, we have

$\displaystyle
G_{0}(N)=\left(\int_{0}^{\pi}\int_{0}^{2\pi}G^{\prime}(\theta,\phi)\sin(\theta)d\theta
d\phi\right)^{-1}.$ (17)

More details on the elements and array radiation pattern is provided in [13,
14].

Figure 2: An illustration of UAVs’ orientation fluctuations in the $x-z$
plane.

###  II-3 The Effect of UAV’s Instabilities

In practical situations, an error in the mechanical control system of UAVs,
mechanical noise, position estimation errors, air pressure, and wind speed can
affect the UAV’s angular and position stability [15, 16]. The instantaneous
orientation of a UAV can randomly deviate from its means denoted by $\theta$.
Let RVs $\theta_{tx,i}\sim\mathcal{N}(0,\sigma^{2}_{\theta})$ and
$\theta_{rx,i}\sim\mathcal{N}(0,\sigma^{2}_{\theta})$ denote respectively the
instantaneous angle of departure (AoD) of Tx and angle of arrival (AoA) of Rx
antennas mounted on $U_{i}$ in the $x-z$ Cartesian coordinates as depicted in
Fig. 2, and RVs $\theta_{ty,i}\sim\mathcal{N}(0,\sigma^{2}_{\theta})$ and
$\theta_{ry,i}\sim\mathcal{N}(0,\sigma^{2}_{\theta})$ denote the instantaneous
AoD of Tx and AoA of Rx antennas mounted on $U_{i}$ in the $y-z$.

##  III Optimal System Design

For a given region with physical parameters such as air pressure, temperature,
$\psi_{s,\text{min}}$, $\psi_{d,\text{min}}$, $H_{s}$, $H_{bs}$, $H_{d}$,
$H_{bd}$, $H_{sd}$ $L_{sd}$, etc, our aim is to adjust the tunable system
parameters such as NFPs’ positions and antenna patterns in such a way that the
cost is minimized. Here, the cost mainly depends on the number of NFPs,
denoted by $M$. From (16), the antenna pattern is adjusted by the number of
antenna elements. Let $N_{t,i}\times N_{t,i}$, $N_{r,i}\times N_{r,i}$,
$N_{s}\times N_{s}$ and $N_{d}\times N_{d}$ denote the number of antenna
elements of $A_{t,i}$, $A_{r,i}$, $A_{s}$, and $A_{d}$, respectively. The
position of $U_{1}$ is characterized by $L_{s}$ and elevation angle $\psi_{s}$
where $\psi_{s,\text{min}}\leq\psi_{s}<\pi/2$. Moreover, the position of
$U_{M}$ is characterized by $L_{d}$ and elevation angle $\psi_{d}$ where
$\psi_{d,\text{min}}\leq\psi_{d}<\pi/2$. Our optimization problem can be
formulated as follows:

$\displaystyle\displaystyle\min_{\begin{subarray}{c}\text{position
of}~{}U_{i},~{}i\in\\{1,...,M\\}\\\ N_{r,i},N_{t,i},~{}i\in\\{1,...,M\\}\\\
N_{s},N_{d}\end{subarray}}~{}~{}~{}~{}~{}~{}~{}{M}$ (18a)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{s.t.}~{}~{}~{}\mathbb{P}_{\text{out}}<\mathbb{P}_{\text{out,tar}}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(18b)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\psi_{s,\text{min}}\leq\psi_{s}<\frac{\pi}{2},\psi_{d,\text{min}}\leq\psi_{d}<\frac{\pi}{2},$
(18c)

where $\mathbb{P}_{\text{out}}$ is the total outage probability which is
calculated below, and $\mathbb{P}_{\text{out,tar}}$ is the target outage
probability which is determined based on the requested QoS.

Lemma 1. Outage probability of the considered system is derived as:

$\displaystyle\mathbb{P}_{\text{out}}\simeq\mathbb{P}_{\text{out,s1}}+\mathbb{P}_{\text{out,Md}}+\sum_{i=2}^{M}\mathbb{P}_{\text{out,i}}$
(19)

where

$\displaystyle\mathbb{P}_{\text{out,s1}}=\sum_{j=1}^{JK}\left(e^{-\frac{2(j-1)^{2}}{J^{2}N_{r,1}^{2}\sigma_{\theta}^{2}}}-e^{-\frac{2j^{2}}{J^{2}N_{r,1}^{2}\sigma_{\theta}^{2}}}\right)\mathbb{Y}(P_{r,\text{th}}-P^{\prime}_{r,1}(j)),$
$\displaystyle\mathbb{P}_{\text{out,Md}}=\sum_{j=1}^{JK}\left(e^{-\frac{2(j-1)^{2}}{J^{2}N_{t,M}^{2}\sigma_{\theta}^{2}}}-e^{-\frac{2j^{2}}{J^{2}N_{t,M}^{2}\sigma_{\theta}^{2}}}\right)\mathbb{Y}(P_{r,\text{th}}-P^{\prime}_{r,d}(j)),$
(20)

$\displaystyle\mathbb{P}_{\text{out,i}}=\sum_{j=1}^{JK}\sum_{j^{\prime}=1}^{JK}\left(e^{-\frac{2(j-1)^{2}}{J^{2}N_{t,i-1}^{2}\sigma_{\theta}^{2}}}-e^{-\frac{2j^{2}}{J^{2}N_{t,i-1}^{2}\sigma_{\theta}^{2}}}\right)$
$\displaystyle~{}\times\left(e^{-\frac{2(j^{\prime}-1)^{2}}{J^{2}N_{r,i}^{2}\sigma_{\theta}^{2}}}-e^{-\frac{2j^{\prime
2}}{J^{2}N_{r,i}^{2}\sigma_{\theta}^{2}}}\right)\mathbb{Y}(P_{r,\text{th}}-P^{\prime}_{r,i}(j,j^{\prime})),$
(21)

and $\mathbb{Y}(x)=\left\\{\begin{array}[]{rl}1&~{}~{}~{}{\rm
for}~{}~{}~{}x\geq 0\\\ 0&~{}~{}~{}{\rm for}~{}~{}~{}x<0\\\
\end{array}\right.$ is the sign function,
$P^{\prime}_{r,1}(j)=P_{t,s}h_{L_{s}}(\psi_{s},L_{s})G_{0}(N_{r,1})G_{0}(N_{s})N_{s}^{2}10^{\frac{G_{\text{max}}}{5}}\mathbb{G}(j,N_{q,i})$
and
$P^{\prime}_{r,i}(j,j^{\prime})=P_{t,i-1}h_{L_{i-1}}G_{0}(N_{t,i-1})G_{0}(N_{r,i})10^{\frac{G_{\text{max}}}{5}}\mathbb{G}(j,N_{t,i-1})\\\
\mathbb{G}(j^{\prime},N_{r,i})\mathbb{G}(j,N_{q,i})$.

###### Proof:

Please refer to Appendix A. ∎

Through the simulations, we will show that the results obtained from (19) are
very close to the results obtained from the Monte-Carlo simulations. More
importantly, the runtime of (19) is much shorter than the Monte-Carlo
simulation, especially for lower values of outage probability.

Without loss of generality and for notation simplicity, we assume that
$H_{sd}\simeq 0$, $H_{s}<<H_{u1}$ and $H_{d}<<H_{ud}$. The effective value of
CU link length is $L_{s,\text{ef}}=L_{s}\cos(\psi_{s})$. Therefore $\psi_{s}$
must be decreased to increase $L_{s,\text{ef}}$. On the other hand, based on
(14), $\psi_{s}$ must be increased to reduce the attenuation loss. As a
results, for any values of $L_{s}$, optimizing $\psi_{s}$ requires balancing
an inherent tradeoff between decreasing $\psi_{s}$ to increase
$L_{s,\text{ef}}$ and increasing it to decrease the attenuation loss. Since
the transmitted power of ground station is higher than the transmitted power
of NFP (due to the limitations of the NFP’s transmitted power), the length of
the CU link is greater than the UD link, so we expect $U_{1}$ to be at a
higher altitude than $U_{M}$. Therefore, in order to provide LoS between
$U_{1}$ and $U_{M}$, it is necessary that $U_{1}$ has a height higher than the
highest obstacle, i.e., $H_{u_{1}}>H_{b,\text{max}}$.

Remark 1. At 70 GHz, the optimal value for $\psi_{s}$ and $\psi_{d}$ is placed
in the following interval:

$\displaystyle H_{\text{scale}}<L_{s}\sin(\psi_{w})<2H_{\text{scale}},$ (22)

where $w\in\\{s,d\\}$. Moreover, when $20^{o}<\psi_{q,\text{min}}<40^{o}$, the
optimal value for $\psi_{q}w$ is close to $\psi_{w,\text{min}}$ and for
$40^{o}<\psi_{w,\text{min}}$, the optimal value for $\psi_{w}$ is equal
$\psi_{w,\text{min}}$.

Remark 2. For $H_{u_{i}}>H_{\text{scale}}$, height changes in the order of a
few hundred meters do not have a significant effect on the amount of channel
attenuation at 70 GHz frequency.

In the next section, by providing simulations for conventional values of
channel parameters at 70 GHz, we will confirm the results of Remarks 1 and 2.
According to the results of Remark 2, for each given $M$, the best place for
NFPs is to be on a straight line as shown in Fig. 3. Assuming that the NFPs
are of the same type, the length of the links between the $U_{i}$s are equal
and is calculated as follows:

$\displaystyle
L_{i}=\frac{\sqrt{(L_{SD}-L_{s,\text{ef}}-L_{d,\text{ef}})^{2}+(H_{u_{1}}-H_{u_{M}})^{2}}}{M-1},$
(23)

where $i\in\\{1,...,M-1\\}$. Moreover, due to the symmetry, the optimal values
for $N_{t,i}$s and $N_{r,i}$s must be the same as

$\displaystyle N_{t,1}=...=N_{t,M-1}=N_{r,2}=...=N_{r,M}=N_{u}.$ (24)

As a summary of the aforementioned results, the tunable parameters of the
considered long aerial backhaul link are listed as

$\displaystyle\left\\{\begin{array}[]{rl}&\\!\\!\\!\\!\\!\text{CU
link:}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\psi_{s,\text{min}}<\psi_{s}<\pi/2,L_{s},N_{r,1},\\\
&\\!\\!\\!\\!\\!\text{UD
link:}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\psi_{d,\text{min}}<\psi_{d}<\pi/2,L_{d},N_{t,M},\\\
&\\!\\!\\!\\!\\!\text{Inter NFPs'
links:}~{}~{}L_{i},N_{u},M.\end{array}\right.$ (28)

Equation (28) is made up of nine variables, and it is clear that finding the
optimal values for the nine variables simultaneously takes a lot of time.
Based on (23), $L_{i}$ is a function of tunable parameters $\psi_{s}$,
$\psi_{d}$, $L_{s}$, $L_{d}$, and $M$, and once these parameters are known,
the exact value of $L_{i}$ is calculated from (23). Moreover, based on the
results of Remark 1, when $\psi_{s\text{max}}>20^{o}$ and
$\psi_{d\text{max}}>20^{o}$, we can set $\psi_{s}\simeq\psi_{s\text{max}}$ and
$\psi_{d}\simeq\psi_{d\text{max}}$. Therefore, the tunable parameters of our
optimization problem can be reduced to six parameters as

$\displaystyle\left\\{L_{s},N_{r,1},~{}~{}~{}L_{d},N_{t,M},~{}~{}~{},N_{u},M\right\\}.$
(29)

Figure 3: As illustrated for each given $M$, assuming that separate frequency
bands are assigned to each link, the best place for NFPs is to be on a
straight line.

Lemma 2. For any given $\psi_{s}$ and $L_{s}$, the optimal value for $N_{r,1}$
can be obtained easily as follows:

$\displaystyle\min_{N_{r,1}}\left|C_{1}h_{L_{s}}(\psi_{s},L_{s})\left(\frac{\sin\left(\frac{N_{r,1}kd_{x}\sin(\frac{2j}{JN_{r,1}})}{2}\right)}{N_{r,1}\sin\left(\frac{kd_{x}\sin(\frac{2j}{JN_{r,1}})}{2}\right)}\right)^{2}-P_{r,\text{th}}\right|,$
$\displaystyle\textrm{s.t.}~{}~{}~{}1-\mathbb{P}_{\text{out,tar}}<\sum_{j=1}^{J^{\prime}}\left(e^{-\frac{2(j-1)^{2}}{J^{2}N_{r,1}^{2}\sigma_{\theta}^{2}}}-e^{-\frac{2j^{2}}{J^{2}N_{r,1}^{2}\sigma_{\theta}^{2}}}\right).$
(30)

Similarly, for any given $\psi_{d}$ and $L_{d}$, we can compute the optimal
value for $N_{t,M}$ from (III) by substituting the parameters $N_{t,M}$,
$N_{d}$, $h_{L_{d}}(\psi_{d},L_{d})$, and $P_{t,M}$ instead of the parameters
$N_{r,1}$, $N_{s}$, $h_{L_{s}}(\psi_{s},L_{s})$, and $P_{t,s}$, respectively.

###### Proof:

Please refer to Appendix B. ∎

Using the results of Lemma 2, it is observed that the 6-dimensional search
space in (29) is simplified to a 4-dimensional space containing variables
$L_{s}$, $L_{d}$, $N_{u}$, and $M$, which significantly reduces the time to
find the optimal parameters.

From the obtained results, we provide Algorithm 1, which reduces the
8-dimensional search space to a 4-dimensional one, which significantly reduces
the optimization time. As discussed, for any given $\sigma_{\theta}$, the
optimal values of $N_{r,1}$ and $N_{t,M}$, independent of the other tunable
parameters, are only a function of $L_{s}$ and $L_{d}$, respectively.
Therefore, in the first two for loops of Algorithm 1, we calculate and store
the optimal vectors of $\bar{N}_{r,1}$ and $\bar{N}_{t,M}$ based on the given
$\sigma_{\theta}$ and intervals $L_{s,\text{min}}\leq L_{s}\leq
L_{s,\text{max}}$, and $L_{d,\text{min}}\leq L_{d}\leq L_{d,\text{max}}$.
Then, each time we increase $M$ and calculate the optimal values for tunable
parameters according to the three-dimensional for loops. The first value of
$M$ for which the end-to-end outage probability is lower than the threshold is
the optimal value for $M$ and Algorithm 1 is terminated. In addition, to
reduce the execution time of the algorithm, instead of the time consuming
Monte Carlo simulations, the equation presented in Lemma 1 is used to
calculate the outage probability.

Algorithm 1 Optimal design of parameters.

Input: $\sigma_{to}$, $L_{sd}$, $\psi_{\textrm{s,\text{min}}}$,
$\psi_{\textrm{d,\text{min}}}$, $f_{c}$, $P_{t,s}$, $P_{t,i}$, $N_{s}$,
$N_{d}$, $H_{b,\text{max}}$, $H_{\text{scale}}$, $H_{bd}$, $G_{\text{max}}$,
$\mathbb{P}_{\text{out,tr}}$;
Output: optimal pattern
$\big{(}N_{r,1}^{\text{opt}},N_{t,M}^{\text{opt}},N_{i}^{\text{opt}}\big{)}$,
optimal NFPs’ positions
$(L_{s}^{\text{opt}},L_{d}^{\text{opt}},L_{i}^{\text{opt}})$, optimal number
of NFPs $M_{\text{opt}}$;
Initialize: $M=1$, $\mathbb{P}_{\text{out,min}}=1$;
for: $L_{s}=H_{b,\text{max}}/\sin({\psi_{s,\text{min}}})$ to
$L_{s,\text{max}}$ do
$|$ find $N_{r,1}$ based on Lemma 2 and store in $\bar{N}_{r,1}$,
end for
for: $L_{d}=H_{bd}/\sin({\psi_{d,\text{min}}})$ to $L_{d,\text{max}}$ do
$|$ find $N_{t,M}$ based on Lemma 2 and store in $\bar{N}_{t,M}$,
end for
while: $\mathbb{P}_{\text{out,min}}>\mathbb{P}_{\text{out,tr}}$ do $M=M+1$,
$|$ for: $L_{s}=H_{b,\text{max}}/\sin({\psi_{s,\text{min}}})$ to
$L_{s,\text{max}}$ do
$|$ find $N_{r,1}$ from vector $\bar{N}_{r,1}$,
$|$ for: $L_{d}=H_{bd}/\sin({\psi_{d,\text{min}}})$ to $L_{d,\text{max}}$ do
$|$ find $N_{t,M}$ from vector $\bar{N}_{t,M}$,
$|$ compute $L_{i}$ based on (23),
$|$ for: $N_{u}=N_{u,\text{min}}$ to $N_{u,\text{max}}$ do
$|$ compute $\mathbb{P}_{\text{out}}$ based on Lemma 1,
$|$ if $\mathbb{P}_{\text{out}}<\mathbb{P}_{\text{out,min}}$ then do
$\mathbb{P}_{\text{out,min}}=\mathbb{P}_{\text{out}}$,
$|$
$N_{r,1}=N_{r,1}^{\text{opt}},N_{t,M}=N_{t,M}^{\text{opt}},N_{i}=N_{i}^{\text{opt}}$,
$|$
$L_{s}=L_{s}^{\text{opt}},L_{d}=L_{d}^{\text{opt}},L_{i}=L_{i}^{\text{opt}}$,
$M=M_{\text{opt}}$
$|$ end if
$|$ end for
$|$ end for
$|$ end for
end while

##  IV Simulation Results

Figure 4: Outage probability of inter NFPs links versus link length for
$H_{u_{i}}=2.5$ km, $f_{c}=70$ GHz, and $\sigma_{\theta}=2^{o}$.

In Fig. 4, outage probability of inter NFPs links is depicted versus link
length for $f_{c}=70$ GHz and different antenna patterns characterized by
$N_{q,i}$. The results of Fig. 4 are obtained for $P_{t,i}=200$ mW,
$\sigma_{\theta}=1.5^{o}$ and $\sigma_{\theta}=2^{o}$. For shorter link
lengths, smaller values of $N_{q,i}$ that create a larger beam width perform
better than larger ones. But as the link length increases, the channel loss
increases and as a result, the average received signal power at the receiver
antenna decreases. Therefore, we expect that the higher antenna pattern gain
performs better for longer link lengths. The results of Fig. 4 confirm the
accuracy of this point.

Oxygen and water vapor can be strongly absorptive of radio signals, especially
at mmWave frequencies and higher. The 60 GHz band is considered as unlicensed
communication bands, however atmospheric attenuation peaks at 60 GHz, with a
value of over 15 dB/km. On the other hand, some frequency bands higher than 60
GHz such as 71–76, 81–86 and 92–95 GHz bands do not suffer from oxygen
absorption, but require a transmitting license in the US from the Federal
Communications Commission (FCC). In Fig. 5, we plot the outage probability of
inter-NFP links in Fig. 5 for different NFP heights and two frequency bands
$f_{c}=60$ and 70 GHz. For $f_{c}=60$ GHz, as the height of the NFPs
increases, it is observed that the maximum achievable link length increases
and thus, this confirms that for the design of the considered system in the 60
GHz frequency band, the optimal choice of flight altitude of the NFPs is of
great importance. However, for $f_{c}=70$ GHz, it is observed that with
increasing the altitude of the NFPs from 2 to 3 km, the performance of the
system does not change significantly.

Figure 5: Comparison of the effect of $H_{u_{i}}$ on the outage probability of
inter NFP links for frequency bands $f_{c}=60$ GHz and $f_{c}=70$ GHz.

Finally, in Table I, for a given region with physical parameters $L_{sd}=40$
km, $\psi_{s,\text{min}}=40^{o}$, $\psi_{d,\text{min}}=20^{o}$,
$H_{b\text{max}}=2$ km, $H_{\text{scale}}=1.5$ km, and $T=20^{o}$C, we find
the optimal values for the channel parameters such as optimal antenna pattern
characterized by
$\big{(}N_{r,1}^{\text{opt}},N_{t,M}^{\text{opt}},N_{i}^{\text{opt}}\big{)}$,
optimal NFPs’ positions characterized by
$(L_{s}^{\text{opt}},L_{d}^{\text{opt}},L_{i}^{\text{opt}})$, optimal number
of NFPs, denoted by $M_{\text{opt}}$. The results of Table I are obtained for
two different NFPs’ instabilities $\sigma_{\theta}=2^{o}$ and
$\sigma_{\theta}=1.5^{o}$. From the results of Table I, for
$\sigma_{\theta}=2^{o}$, the minimum required number of NFPs is $M=10$ to
guarantee $\mathbb{P}_{\text{out}}<10^{-3}$ in the considered 40 km distance
between the core network and the remote area, while for NFPs with higher
stability $\sigma_{\theta}=1.5^{o}$, the minimum number of NFPs is reduced to
$M=7$.

TABLE I: Optimal parameter design for $L_{sd}=40$ km, and two different NFPs’ stability $\sigma_{\theta}=2^{o}$ and $1.5^{o}$. $\sigma_{\theta}$ | $N_{r,1}^{\text{opt}}$ | $N_{t,M}^{\text{opt}}$ | $N_{i}^{\text{opt}}$ | $L_{s}^{\text{opt}}$ | $L_{d}^{\text{opt}}$ | $L_{i}^{\text{opt}}$ | $M_{\text{opt}}$
---|---|---|---|---|---|---|---
|  |  |  | (km) | (km) | (km) |
$2^{o}$ | 6 | 6 | 8 | 9.6 | 6.2 | 3.1 | 10
$1.5^{o}$ | 8 | 8 | 10 | 7.4 | 11.4 | 4.3 | 7

## Appendix A

For the considered multihop DF relaying system, the outage probability of the
end-to-end system is given by [17]

$\displaystyle\mathbb{P}_{\text{out}}$
$\displaystyle=\text{Prob}\left\\{\text{min}(P_{r,1},...,P_{r,M},P_{r,d})<P_{r,\text{th}}\right\\}$
$\displaystyle=1-(1-\mathbb{P}_{\text{out,s1}})(1-\mathbb{P}_{\text{out,Md}})\prod_{i=2}^{M}(1-\mathbb{P}_{\text{out,i}}),$
(31)

where $\mathbb{P}_{\text{out,s1}}$, $\mathbb{P}_{\text{out,Md}}$, and
$\mathbb{P}_{\text{out,i}}$ denote the outage probability of $A_{s}$ to
$A_{r,1}$ link, $A_{t,M}$ to $A_{d}$ link, and $A_{t,i-1}$ to $A_{r,i}$ link,
respectively, and $P_{r,\text{th}}$ is the received power threshold. For lower
values of outage probability, (A) can be approximated as (19).

Considering the approximate symmetry of the antenna pattern along the $\psi$
axis, (15) can be approximated as

$\displaystyle G(\theta_{q_{i}})$
$\displaystyle=G_{0}(N_{q,i})10^{\frac{G_{\text{max}}}{10}}\left(\frac{\sin\left(\frac{N_{q,i}kd_{x}\sin(\theta_{qU_{i}})}{2}\right)}{N_{q,i}\sin\left(\frac{kd_{x}\sin(\theta_{qU_{i}})}{2}\right)}\right)^{2},$
(32)

where $q\in{t,r}$,
$\theta_{q_{i}}=\sqrt{\theta_{qx_{i}}^{2}+\theta_{qy_{i}}^{2}}$. Since RVs
$\theta_{qx_{i}}$ and $\theta_{qy_{i}}$ have Gaussian distributions, the
random variable $\theta_{q_{i}}$ follows the Rayleigh distribution as

$\displaystyle
f_{\theta_{q_{i}}}(\theta_{q_{i}})=\frac{\theta_{q_{i}}}{\sigma_{\theta}^{2}}\exp\left(-\theta_{q_{i}}^{2}/2\sigma_{\theta}^{2}\right).$
(33)

Let us approximate (32) as

$\displaystyle G(\theta_{q_{i}})\simeq
G_{0}(N_{q,i})10^{\frac{G_{\text{max}}}{10}}\sum_{j=1}^{KJ}\mathbb{G}(j,N_{q,i})$
$\displaystyle~{}~{}~{}\times\left[\mathbb{Y}\left(\theta_{q,i}-\frac{2(j-1)}{JN_{q,i}}\right)-\mathbb{Y}\left(\theta_{q,i}-\frac{2j}{JN_{q,i}}\right)\right],$
(34)

where
$\mathbb{G}(j,N_{q,i})=\left(\frac{\sin\left(\frac{N_{q,i}kd_{x}\sin(\frac{2j}{JN_{q,i}})}{2}\right)}{N_{q,i}\sin\left(\frac{kd_{x}\sin(\frac{2j}{JN_{q,i}})}{2}\right)}\right)^{2}$,
$\mathbb{Y}(x)=\left\\{\begin{array}[]{rl}1&~{}~{}~{}{\rm for}~{}~{}~{}x\geq
0\\\ 0&~{}~{}~{}{\rm for}~{}~{}~{}x<0\\\ \end{array}\right.$ is the sign
function, and the parameters $J$ and $K$ are the natural numbers that for
large values of $J$, (A) tends to (32). Also, $K=1$ refers to the main lobe of
the antenna pattern and $K>1$ refers to the number of sidelobes. Note that the
problem of orientation fluctuations is related to the antennas mounted on the
NFPs. For the ground antennas $A_{s}$ and $A_{d}$, it is assumed that the
ground station does not face weight and power limitations and uses a
stabilizer with high accuracy and response speed to track the first and second
drones. Based on this assumption, and using (5), (33), and (A), the outage
probability of $A_{s}$ to $A_{1}$ link is obtained in Lemma 1. From (5) and
(A), we have

$\displaystyle P_{r,i}\simeq
P_{t,i-1}h_{L_{i-1}}G_{0}(N_{t,i-1})G_{0}(N_{r,i})10^{\frac{G_{\text{max}}}{5}}\sum_{j=1}^{KJ}\sum_{j^{\prime}=1}^{KJ}\mathbb{G}(j,N_{t,i-1})$
$\displaystyle\mathbb{G}(j^{\prime},N_{r,i})\left[\mathbb{Y}\left(\theta_{t,i-1}-\frac{2(j-1)}{JN_{t,i-1}}\right)-\right.\left.\mathbb{Y}\left(\theta_{t,i-1}-\frac{2j}{JN_{t,i-1}}\right)\right]$
$\displaystyle\left[\mathbb{Y}\left(\theta_{r,i}-\frac{2(j^{\prime}-1)}{JN_{r,i}}\right)-\mathbb{Y}\left(\theta_{r,i}-\frac{2j^{\prime}}{JN_{r,i}}\right)\right].$
(35)

As we see, $P_{r,i}$ is a function of two independent RVs $\theta_{t,i-1}$ and
$\theta_{r,i}$. From (33) and (A) and using [18], after some manipulations,
the outage probability of $A_{t,i-1}$ to $A_{r,i}$ link is obtained in Lemma
1.

## Appendix B

Next, we show that for any given values of $L_{w}$ and $\psi_{w}$ where
$w\in\\{s,d\\}$, $h_{L_{w}}$, we can obtain optimal values for $N_{r,1}$ and
$N_{t,M}$. In the sequel, we prove this for $N_{r,1}$ and similarly we can
prove this for $N_{t,M}$. For lower values of NFP’s orientation fluctuations,
we can approximate $P^{\prime}_{r,1}(j)$ with a good accuracy as

$\displaystyle P^{\prime}_{r,1}(j)\simeq$ (36)
$\displaystyle\left\\{\begin{array}[]{rl}&\\!\\!\\!\\!\\!\\!\\!C_{1}h_{L_{s}}(\psi_{s},L_{s})\left(\frac{\sin\left(\frac{N_{r,1}kd_{x}\sin(\frac{2j}{JN_{r,1}})}{2}\right)}{N_{r,1}\sin\left(\frac{kd_{x}\sin(\frac{2j}{JN_{r,1}})}{2}\right)}\right)^{2},~{}~{}~{}~{}j<J,\\\
&\\!\\!\\!\\!\\!\\!\\!0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}j>0,\end{array}\right.$
(39)

where $C_{1}$ is a constant. Based on (36), it can easily be shown that
$\frac{\text{d}P^{\prime}_{r,1}(j)}{\text{d}j}\leq 0$ and, therefore,
$P^{\prime}_{r,1}(j)$ is a descending function of $j$. Using this and after
some manipulations, we obtain

$\displaystyle\mathbb{P}^{\prime}_{\text{out,s1}}(N_{r,1})\simeq
1-\sum_{j=1}^{J^{\prime}}\left(e^{-\frac{2(j-1)^{2}}{J^{2}N_{r,1}^{2}\sigma_{\theta}^{2}}}-e^{-\frac{2j^{2}}{J^{2}N_{r,1}^{2}\sigma_{\theta}^{2}}}\right),$
(40)

where $J^{\prime}$ can be obtained easily from (36) when
$P^{\prime}_{r,1}(J^{\prime})=P_{r,\text{th}}$. Based on these, for any given
$h_{L_{s}}(\psi_{s},L_{s})$, independent of other tunable parameters, any
$N_{r,1}$ that satisfies the conditions of (III) is the optimal value.

## References

  * [1] M. Alzenad, M. Z. Shakir, H. Yanikomeroglu, and M.-S. Alouini, “Fso-based vertical backhaul/fronthaul framework for 5g+ wireless networks,” _IEEE Communications Magazine_ , vol. 56, no. 1, pp. 218–224, 2018.
  * [2] W. Khawaja, I. Guvenc, D. W. Matolak, U.-C. Fiebig, and N. Schneckenberger, “A survey of air-to-ground propagation channel modeling for unmanned aerial vehicles,” _IEEE Commun. Surveys Tuts._ , vol. 21, no. 3, pp. 2361–2391, May 2019.
  * [3] M. T. Dabiri, H. Safi, S. Parsaeefard, and W. Saad, “Analytical channel models for millimeter wave UAV networks under hovering fluctuations,” _IEEE Trans. Wireless Commun._ , Jan. 2020.
  * [4] M. T. Dabiri, M. Rezaee, V. Yazdanian, B. Maham, W. Saad, and C. S. Hong, “3D channel characterization and performance analysis of UAV-assisted millimeter wave links,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 110–125, 2020.
  * [5] M. Gapeyenko, V. Petrov, D. Moltchanov, S. Andreev, N. Himayat, and Y. Koucheryavy, “Flexible and reliable UAV-assisted backhaul operation in 5G mmWave cellular networks,” _IEEE J. Sel, Areas Commun._ , vol. 36, no. 11, pp. 2486–2496, 2018.
  * [6] B. Galkin, J. Kibilda, and L. A. DaSilva, “Backhaul for low-altitude UAVs in urban environments,” in _Proc. IEEE International Conference on Communications (ICC), Kansas City, MO, USA_ , May 2018.
  * [7] N. Tafintsev, D. Moltchanov, M. Gerasimenko, M. Gapeyenko, J. Zhu, S.-p. Yeh, N. Himayat, S. Andreev, Y. Koucheryavy, and M. Valkama, “Aerial access and backhaul in mmWave B5G systems: Performance dynamics and optimization,” _IEEE Commun. Mag._ , vol. 58, no. 2, pp. 93–99, 2020.
  * [8] W. Wang, N. Cheng, Y. Liu, H. Zhou, X. Lin, and X. Shen, “Content delivery analysis in cellular networks with aerial caching and mmWave backhaul,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 5, pp. 4809–4822, 2021.
  * [9] C. T. Cicek, H. Gultekin, B. Tavli, and H. Yanikomeroglu, “Backhaul-aware optimization of UAV base station location and bandwidth allocation for profit maximization,” _IEEE Access_ , vol. 8, pp. 154 573–154 588, 2020.
  * [10] Z. Feng, L. Ji, Q. Zhang, and W. Li, “Spectrum management for mmwave enabled UAV swarm networks: Challenges and opportunities,” _IEEE Commun. Mag._ , vol. 57, no. 1, pp. 146–153, 2018.
  * [11] Y. Yu, X. Bu, K. Yang, H. Yang, and Z. Han, “UAV-aided low latency mobile edge computing with mmWave backhaul,” in _ICC 2019-2019 IEEE International Conference on Communications (ICC)_. IEEE, 2019, pp. 1–7.
  * [12] “Attenuation by atmospheric gases,” Recommendation ITU-R P.676-3.
  * [13] 3GPP TR 37.840 v12.1.0, “Technical specification group radio access network; study of radio frequency (RF) and electromagnetic compatibility (EMC) requirements for active antenna array system (AAS) base station,” _Tech. Rep._ , 2013.
  * [14] C. A. Balanis, _Antenna theory: analysis and design_. John wiley & sons, 2016.
  * [15] M. T. Dabiri, S. M. S. Sadough, and M. A. Khalighi, “Channel modeling and parameter optimization for hovering UAV-based free-space optical links,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 2104–2113, 2018.
  * [16] M. Dabiri, S. Khankalantary, M. Piran, I. Ansari, M. Uysal, W. Saad, and C. Hong, “UAV-assisted free space optical communication system with amplify-and-forward relaying,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 9, pp. 8926–8936, 2021.
  * [17] M. O. Hasna and M.-S. Alouini, “Outage probability of multihop transmission over Nakagami fading channels,” _IEEE Commun. Let._ , vol. 7, no. 5, pp. 216–218, 2003.
  * [18] A. Jeffrey and D. Zwillinger, _Table of Integrals, Series, and Products_. Elsevier, 2007.


# Designer spin models in tunable two-dimensional nanographene lattices

J. C. G. Henriques1,2, Mar Ferri-Cortés3, J. Fernández-Rossier2 1International
Iberian Nanotechnology Laboratory (INL), Av. Mestre José Veiga, 4715-330
Braga, Portugal  2 Universidade de Santiago de Compostela, 15782 Santiago de
Compostela, Spain  3 Departamento de Física Aplicada, Universidad de Alicante,
03690 San Vicente del Raspeig, Spain  On permanent leave from Departamento de
Física Aplicada, Universidad de Alicante, 03690 San Vicente del Raspeig,
Spain.

###### Abstract

Motivated by recent experimental breakthroughs, we propose a strategy to
design two-dimensional spin lattices with competing interactions that lead to
non-trivial emergent quantum states. We consider $S=1/2$ nanographenes with
$C_{3}$ symmetry as building blocks, and we leverage the potential to control
both the sign and the strength of exchange with first neighbours to build a
family of spin models. Specifically, we consider the case of a Heisenberg
model in a triangle-decorated honeycomb lattice with competing ferromagnetic
and antiferromagnetic interactions whose ratio can be varied in a wide range.
Based on exact diagonalization of both fermionic and spin models we predict a
quantum phase transition between a valence bond crystal of spin singlets with
triplon excitations living in a Kagomé lattice and a Néel phase of effective
$S=3/2$ in the limit of dominant ferromagnetic interactions.

Artificial lattices, created through precision engineering at the nanoscale,
can serve as powerful tools for investigating emergent quantum phenomena
across various physical platforms, such as ultra-cold atoms in optical
latticesGross and Bloch (2017), trapped ionsMonroe  _et al._ (2021), Rydberg
atomsScholl  _et al._ (2021), superconducting circuitsHouck  _et al._ (2012),
semiconductor quantum dotsHensgens  _et al._ (2017) and magnetic
adatomsKhajetoorians  _et al._ (2019). These lattices mimic the periodic
structures found in condensed matter systems and can help us to emulate and
understand exotic quantum states, predicted by theoretical modelsKitaev
(2006); Wei  _et al._ (2011), that may be challenging to observe directly in
natural materials and can play a crucial role in the development of quantum
coherent nanotechnologyHeinrich  _et al._ (2021).

Here, we propose to leverage the gigantic know-how of organic chemistry using
open-shell nanographenes as building blocks to create artificial quantum
matterde Oteyza and Frederiksen (2022). Open shell nanographenes are planar
carbon molecules with a finite spin $S$ in the ground state. The formation of
supramolecular structures, such as dimersMishra  _et al._ (2020); Krane  _et
al._ (2023), trimersDu  _et al._ (2023), ringsHieulle  _et al._ (2021),
chainsMishra  _et al._ (2021) and two-dimensional arraysDelgado  _et al._
(2023), that preserve the open shell nature of the building blocksCatarina
_et al._ (2023), has recently been demonstrated. This bottom-up approach
permits one to designOrtiz (2023); Ortiz  _et al._ (2023); Catarina  _et al._
(2023) two-dimensional carbon crystals with very narrow bands, linear
combination of carbon $p_{z}$ orbitals, and strong electron correlations,
complementing the top-bottom approach of magic angle twisted bilayer graphene
Cao  _et al._ (2018a, b).

In the case of one-dimensional nanographene chains, quantum fluctuations
destroy long-range order, and outstanding quantum magnetism phenomena, such as
the valence bond solid (VBS) Affleck  _et al._ (1987) with fractional edge
states in $S=1$ chainsMishra  _et al._ (2021), and the topological dimerized
phase of the alternate bond $S=1/2$ chainZhao , have been reported. In two
dimensions, spin lattices can display either broken symmetry phases with
magnetic order and gapless spin excitations as well as gapped quantum
disordered ground states, some of which can have emergent properties, and have
potential for applications in quantum information Kitaev (2006); Wei  _et al._
(2011).

The properties of both ground state and spin excitations depend crucially on
the lattice as well as the sign, symmetry and strength of the interactions
therein. For instance, experiments show that $S=3/2$ triangulenes naturally
form honeycomb latticesDelgado  _et al._ (2023) and theory showsCatarina  _et
al._ (2023) that these are described with non-linear Heisenberg spin models
that belong to the same class of the AKLT Hamiltonian Affleck  _et al._ (1987)
with a VBS ground state. However, calculations Catarina  _et al._ (2023) show
that the non-linear exchange is way too small in this system when compared
with the AKLT model, which very likely will drive the system towards a Néel
ground state Ganesh  _et al._ (2011) with gapless magnon excitations.

Here we pose the challenge of engineering a planar nanographene lattice that,
depending on structural design parameters, can display either a broken
symmetry ground state with gapless spin excitations, or quantum disordered
ground state with a gapped spectrum. For that matter, we propose a system that
realizes the so called star lattice, a triangle decorate honeycomb graph, with
$S=1/2$ spins (See Figure 1). We choose as our $S=1/2$ units phenalenyl
molecules Reid (1965), known to form open-shell dimersKrane  _et al._ (2023)
and trimersDu  _et al._ (2023). In each site of the honeycomb lattice we place
three phenalenyl, coupled via a benzene linker that, as we show, induces a
weak ferromagnetic exchange. Each phenalenyl in a trimer is coupled with one
additional molecule from a different site of the honeycomb lattice, in a
configuration known to give antiferromagnetic exchange Jacob and Fernández-
Rossier (2022); Krane  _et al._ (2023).

We start our discussion with the introduction of the Hubbard model, which we
shall use to describe the phenalenyl dimers and trimers considered in this
work. The use of Hubbard model to model planar open-shell nanographenes is
known to yield good agreement both with quantum chemistryFernández-Rossier and
Palacios (2007); Ortiz  _et al._ (2019) and experimental resultsMishra  _et
al._ (2020, 2021).

Explicitly, the model readsArovas  _et al._ (2022):

$\displaystyle{\cal
H}=\sum_{i,j,\sigma}t_{i,j}c^{\dagger}_{i\sigma}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow},$
(1)

where the indices $i,j$ run over carbon atoms (which are first and third
neighbors, only), $t_{i,j}$ stands for the hopping between sites $i$ and $j$,
and $U$ is the on-site Hubbard repulsion. In the rest of the paper we shall
consider the first neighbor hopping to be $t=-2.7$ eV, the third neighbor
hopping to be $t_{3}=t/10$ and we set the Hubbard repulsion to $U=|t|$, which
has been shown to agree well with _ab initio_ results in previous works Ortiz
_et al._ (2022); Catarina  _et al._ (2023). The operators
$c^{\dagger}_{i\sigma}$ ($c_{i\sigma}$) represent the creation (annihilation)
of an electron in site $i$ with spin projection along a quantization axis
$\sigma=\uparrow,\downarrow$, and
$n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the corresponding number
operator.

For a phenalenyl molecule, due to its sublattice imbalance of 1 (one
sublattice contains one additional atom when compared with the other), a
single mode at zero energy is expected in its single-particle spectrum
Sutherland (1986). Accounting for the Hubbard repulsion, Lieb’s theorem Lieb
(1989) predicts a many-body ground state with spin $S=1/2$, in agreement with
the Ovchinnikov rule Ovchinnikov (1978). This has been verified experimentally
Turco  _et al._ (2023) and theoretically Su  _et al._ (2019); Ahmed and Mandal
(2022) in previous works.

When phenalenyl molecules are combined to form supra-molecular structures,
such as dimers or trimers (see Fig. 1a and b), Lieb’s theoremLieb (1989)
predicts that the many-body ground state will have a spin which is larger or
smaller than $S=1/2$ depending on how the molecules are connected. In
particular, if three phenalenyl molecules are linked together via a central
benzene ring, preserving $C_{3}$ symmetry, Lieb’s theorem states that the
three $S=1/2$ spins couple ferromagnetically, leading to a many-body ground
state with $S=3/2$; this is a consequence of the sublattice imbalance
increasing to 3, and can be microscopically explained via Hund’s ferromagnetic
exchange between the zero modes hosted by the the three phenalenyl building
blocks Ortiz  _et al._ (2019); Jacob and Fernández-Rossier (2022).

Figure 1: Schematic representation of a phenalenyl trimer a), and dimer b)
with different number of benzene rings, $n$, between the two molecules; the
effective $S=1/2$ spin models, and the predicted spin of the many-body ground
state according to Lieb’s theorem are also shown; c) Ratio of the trimer
ferromagnetic and the dimer antiferromagnetic exchange for different number of
benzene rings in the dimer (results obtained from a CAS calculation with the 4
least energetic orbitals in the Hilbert space); d) Depiction of a fraction of
a 2D crystal obtained by combining phenalenyl trimers; the unit-cell is marked
in orange.

In order to estimate the magnitude of the effective ferromagnetic interaction
that leads to a $S=3/2$ ground state, we solve the Hubbard model for the
trimer of Fig. 1a using the Configuration Interaction (CI) method in the
Complete Active Space (CAS) approximation, where the Hilbert space is
restricted to the three single particle states at zero energy Ortiz  _et al._
(2019); sup . Doing so, we find that, in agreement with Lieb’s theorem, the
molecule presents a quartet ground state corresponding to a $S=3/2$ spin; the
first excited state is also a quartet, made of two degenerate doublets. Based
on the low energy spectrum of the CAS calculation, we postulate a spin model
to describe this system, consisting of three $S=1/2$ spins coupled
ferromagnetically all-to-all, with an exchange coupling given approximately by
$J_{\textrm{FM}}\approx-3$ meV.

Contrarily to the trimer, the phenalenyl dimer has no sublattice imbalance,
and as a consequence its many-body ground state is a singlet. Microscopically,
this implies that intermolecular exchange is antiferromagnetic Jacob and
Fernández-Rossier (2022); Krane  _et al._ (2023), and favors low spin ground
states. Performing a CAS calculation for the dimer (where only the zero modes
of the two phenalenyl are accounted for in the Hilbert space), one indeed
finds a singlet ground state split by approximately $63$ meV from a triplet
excited state sup , in qualitative agreement with recent experimentsKrane  _et
al._ (2023); at higher energies, well separated from these low energy states,
other excitations appear. The results of the multiconfigurational calculation
lead us to describe this molecule as two antiferromagnetically coupled spin
$1/2$ with an antiferromagnetic exchange of $J_{\textrm{AFM}}\approx 63$ meV.
Importantly, by introducing benzene rings between the molecules, our
calculations show that the magnitude of the antiferromagnetic intermolecular
exchange is reduced dramatically (see Fig. 1c). This reduction has also been
reported in the case of $S=1$ triangulene dimersMishra  _et al._ (2020).

Based on the analysis we have just performed, one sees that if two phenalenyl
trimers are combined in the way depicted in Fig. 1d, then it is possible to
build a 2-dimensional nanographene crystal with competing ferromagnetic and
antiferromagnetic interactions in the unit cell. Since we have argued that the
building blocks of such a crystal are well described with spin models, it is
only natural to extend the same reasoning to the crystal. Hence, we propose
the following spin $1/2$ Heisenberg Hamiltonian in the star lattice to
describe the structure depicted in Fig. 1d:

$H=\sum_{i,i^{\prime}}J_{\textrm{FM}}\vec{S}_{i}\cdot\vec{S}_{i^{\prime}}+\sum_{i,j}J_{\textrm{AFM}}\vec{S}_{i}\cdot\vec{S}_{j}$
(2)

where $i$ and $i^{\prime}$ are first neighbours belonging to the same trimer,
while $i$ and $j$ are first neighbours which belong to adjacent trimers. The
values of $J_{\textrm{FM}}$ and $J_{\textrm{AFM}}$ are obtained from our
Hubbard-CAS calculations for the trimer and dimer, respectively. In Fig. 2a we
depict the spin model for the crystal, where the geometry of a star lattice is
evident.

In order to validate the proposed spin model, we compared the results obtained
from exact diagonalization for the unit cell of the crystal with periodic
boundary conditions, to the ones obtained with CAS for the same structure for
the case of $n=0$ benzene spacers in the dimer. An excellent agreement between
the two approaches is seen in Fig. 2b, thus validating the proposed spin
model. In the CAS calculation only the six single particle states closer to
zero energy were considered, as these are expected to give the biggest
contribution to the low energy part of the many-body spectrum.

Figure 2: a) Schematic representation of the spin model for the crystal,
forming a star lattice. The dashed rectangle highlights the unit cell of the
crystal b) Comparison between the results obtained from CAS(6,6) and exact
diagonalization of the spin model for the unit cell (u.c.) of the crystal with
periodic boundary conditions (p.b.c.)

The proposed star lattice model has two competing interactions,
$J_{\textrm{FM}}$ and $J_{\textrm{AFM}}$, whose ratio can be controlled, in
the proposed system, by changing the number of benzene spacers, $n$, in the
phenalenyl dimers. We now discuss two limiting cases of the model,
$J_{\textrm{FM}}=0$ and $J_{\textrm{AFM}}=0$. In the former, the star lattice
becomes a Kagomé lattice of independent pairs of antiferromagnetically coupled
phenalenyl molecules. The ground state wave function in this cases is a
valence bond crystal (VBC), with an excited manifold of $3N$ degenerate $S=1$
excitations with energy $J_{\textrm{AFM}}$, where $N$ is the number of unit
cells in the crystal. When $J_{\textrm{FM}}$ is turned on, but remains much
smaller than $J_{\textrm{AFM}}$, the ground state remains close to a VBC and
the $S=1$ excited states acquire a dispersion, forming triplon bands (more on
this below); this explains the energy splitting seen in the excited manifold
in Fig. 2b. Solving the model analytically, one finds the splitting of the
first excited manifold to be given by $3J_{\textrm{FM}}/2$, in good agreement
with the numerical result. In the opposite limit, i.e. $J_{\textrm{AFM}}=0$,
we find $2N$ independent trimers with $S=3/2$ spin. When $J_{\textrm{AFM}}$ is
turned on, we are left with a honeycomb lattice of $S=3/2$ spins with
antiferromagnetic interactions. This broken symmetry limit has been discussed
by Owerre Owerre (2017), who found topological magnons.

The star lattice Hamiltonian has been studied before in the case where both
interactions are antiferromagnetic Yang  _et al._ (2010); Jahromi and Orús
(2018). In our case, however, the interactions within each triangle are
ferromagnetic instead. To infer the properties of this system we now exactly
diagonalize the spin model for a ring made of six trimers, containing a total
of 18 spins, with periodic boundary conditions (see Fig. 2a). We focus on two
aspects, the nature of the ground state wave function and the low energy
excitation spectrum.

To gain insight about the nature of the ground state of this system, we
compute the overlap between its wave function and another two particularly
relevant states: i) the wave function of the VBC; ii) the ground state of the
$S=3/2$ antiferromagnetic Heisenberg model. For the first case, we simply
build the VBC wave function by solving the spin model with
$J_{\textrm{FM}}=0$, and then compute the overlap with the solution for
several values of $J_{\textrm{FM}}/J_{\textrm{AFM}}$. Regarding the second
case, we solve the $S=3/2$ Heisenberg model which we then project into a
spin-1/2 basis, thus allowing the computation of the overlap with the wave
function of our system; the results are displayed in Fig. 3b.

Figure 3:  a) Overlap between the wave function of the ground state of a ring
of 6 trimers, with a total of 18 spins, with p.b.c.,
$|\psi_{\textrm{GS}}\rangle$, and the wave functions of the valence bond
crystal (VBC) and the $S=3/2$ Heisenberg model ground state for the same
geometry. The dashed vertical line marks the point where the two overlaps are
identical, and the dotted lines refer to the values of
$J_{\textrm{FM}}/J_{\textrm{AFM}}$ for different number number of benzene
spacers in the phenalenyl dimers (c.f. Fig. 1c); b) Energy of the first
excitation, $\Delta$, as a function of $J_{\textrm{FM}}/J_{\textrm{AFM}}$; c)
Exact diagonalization of spin model the 6 trimer ring with
$J_{\textrm{FM}}=-10J_{\textrm{AFM}}$ for the $S=1/2$ star lattice and
$\mathcal{J}^{\textrm{eff.}}_{S=3/2}=J_{\textrm{AFM}}/9$ for the $S=3/2$
honeycomb geometry.

There, we see that, as expected, the overlap between the exact ground state
wave function and the one from the VBC starts at 1 for $J_{\textrm{FM}}=0$,
and quickly decreases as the ratio of ferro- and antiferromagnetic
interactions increases. At odds, with this, the overlap with the $S=3/2$
Heisenberg model wave function starts at approximately 0.14 and increases as
$|J_{\textrm{FM}}/J_{\textrm{AFM}}|$ grows. In the limit $|J_{\textrm{FM}}|\gg
J_{\textrm{AFM}}$ the overlap approaches 1, as expected given the physical
picture presented before. Crucially, by controlling the molecular linking
groups it seems possible to induce a quantum phase transition in the system,
as its ground states changes from a quantum disordered one in the VBC limit,
to a broken symmetry one in the Heisenberg limit.

In Fig. 3b and c we focus on the excitation spectrum obtained via exact
diagonalization of the spin model for the 18 spin cluster we have just
considered. For all considered values of $J_{\textrm{FM}}/J_{\textrm{AFM}}$ we
find a singlet ground state followed by a triplet excitation. In panel b we
plot the energy of the first excited manifold, $\Delta$, as a function of the
ratio of ferro- and antiferromagnetic interactions. In the limit
$|J_{\textrm{FM}}|/J_{\textrm{AFM}}\ll 1$, where the system is close to a VBC
with a quantum disordered ground state, the gap of the system is essentially
$J_{\textrm{AFM}}$ and the first set of excited states corresponds to a set of
almost degenerate triplets, slightly split due to the finite
$J_{\textrm{FM}}$. As $|J_{\textrm{FM}}|/J_{\textrm{AFM}}$ increases, and the
system transitions from a VBC to a $S=3/2$ Heisenberg model, the gap of the
first excitation decreases sharply, and saturates at ${\cal
J}_{S=3/2}^{\textrm{eff.}}=J_{\textrm{AFM}}/9$. The finite value of this gap
is due to the finite size of our simulation lattice. The relation between
${\cal J}_{S=3/2}^{\textrm{eff.}}$ and $J_{\textrm{AFM}}$ is identical to the
one expected from a classical model sup . In Figure 3c we compare the low
energy spectrum of the $S=1/2$ star lattice, using
$J_{\textrm{FM}}=-10J_{\textrm{AFM}}$, with the energy spectrum of a $S=3/2$
honeycomb Heisenberg model with exchange ${\cal J}_{S=3/2}^{\textrm{eff.}}$.
There one finds that the two data set are in good agreement; we have verified
that as we approach the limit
$|J_{\textrm{FM}}|/J_{\textrm{AFM}}\rightarrow\infty$ an exact agreement
appears sup .

We now focus on the excitations in the small
$|J_{\textrm{FM}}|/J_{\textrm{AFM}}$ limit. These excitations are triplons,
which correspond to $S=1$ excitations that propagate in the Kagomé lattice
defined by the bonds of the honeycomb lattice of antiferromangetically coupled
phenalenyls. In this limit, we can employ the mean-field self-consistent bond-
operator formalism introduced by Sachdev and Bhatt Sachdev and Bhatt (1990);
Gopalan  _et al._ (1994); Normand and Rice (1996). Within this framework, we
assign a singlet and three triplet creation operators to each
antiferromagnetically coupled dimer of the unit cell. Expressing the
Hamiltonian of the spin model in terms of these bosonic bond operators, and
diagonalizing it via a paraunitary transformation Pires (2021), we find the
triplon dispersion depicted in Fig. 4 sup . These bands follow the usual
dispersion of Kagomé systems, with two dispersive bands, and a single flat
one. This set of three bands is centered approximately at the energy of the
antiferromagnetic exchange, and the bandwidth is controlled by the
ferromagnetic coupling. Reducing the strength of the antiferromagnetic
coupling, while still keeping it larger than the ferromagnetic interaction,
simply leads to a downwards shift of the band structure; further decreasing
$J_{\textrm{AFM}}$ would originate an instability in the triplon spectrum,
indicating the transition from a quantum disordered ground state to a broken
symmetry one.

Figure 4: (a) Triplon dispersion obtained with the mean-filed self-consistent
bond operator formalism. Solid black lines show the bands for the case with
$J_{\textrm{AFM}}=63$ meV and dashed orange lines refer to the case with
$J_{\textrm{AFM}}=20$ meV; in both cases $J_{\textrm{FM}}=-3$ meV; (b)
Schematic representation of the effective Kagomé lattice formed by the
phenalenyl pairs which link different trimers.

In summary, we propose a strategy to create artificial two-dimensional spin
lattices with tunable competing interactions using nanographenes.
Specifically, we have shown that phenalenyls can be used as building blocks to
realize the $S=1/2$ star lattice model with competing ferro-
($J_{\textrm{FM}}$) and antiferromagnetic ($J_{\textrm{AFM}}$) couplings. We
have shown that the figure of merit of this model, the ratio
$\frac{J_{\textrm{AFM}}}{J_{\textrm{FM}}}$, can be tuned in a range that
permits one to cross the transition between two very different phases, a
$S=1/2$ valence bond crystal state with gapped triplons, and a $S=3/2$ Néel
ordered phase. Our results illustrate the potential of nanographene based
artificial lattices to apply the gigantic power of synthetic organic chemistry
to engineer non-trivial quantum spin states.

Acknowledgements.  We acknowledge António Costa and Bruno Murta for fruitful
discussions and support with numerical calculations. MFC acknowledges funding
from Generalitat Valenciana (CIACIF/2021/434). J.F.-R., J.C.G.H. acknowledge
financial support from FCT (Grant No. PTDC/FIS-MAC/2045/2021), SNF Sinergia
(Grant Pimag, CRSII5-205987) and the European Union (Grant FUNLAYERS-
101079184). J.F.-R. acknowledges funding from Generalitat Valenciana
(Prometeo2021/017 and MFA/2022/045) and MICIN-Spain (Grants No.
PID2019-109539GB-C41 and PRTR-C1y.I1)

## References

  * Gross and Bloch (2017) Christian Gross and Immanuel Bloch, “Quantum simulations with ultracold atoms in optical lattices,” Science 357, 995–1001 (2017).
  * Monroe  _et al._ (2021) Christopher Monroe, Wes C Campbell, L-M Duan, Z-X Gong, Alexey V Gorshkov, Paul W Hess, Rajibul Islam, Kihwan Kim, Norbert M Linke, Guido Pagano, _et al._ , “Programmable quantum simulations of spin systems with trapped ions,” Reviews of Modern Physics 93, 025001 (2021).
  * Scholl  _et al._ (2021) Pascal Scholl, Michael Schuler, Hannah J Williams, Alexander A Eberharter, Daniel Barredo, Kai-Niklas Schymik, Vincent Lienhard, Louis-Paul Henry, Thomas C Lang, Thierry Lahaye, _et al._ , “Quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms,” Nature 595, 233–238 (2021).
  * Houck  _et al._ (2012) Andrew A Houck, Hakan E Türeci, and Jens Koch, “On-chip quantum simulation with superconducting circuits,” Nature Physics 8, 292–299 (2012).
  * Hensgens  _et al._ (2017) Toivo Hensgens, Takafumi Fujita, Laurens Janssen, Xiao Li, CJ Van Diepen, Christian Reichl, Werner Wegscheider, Sankar Das Sarma, and Lieven MK Vandersypen, “Quantum simulation of a fermi–hubbard model using a semiconductor quantum dot array,” Nature 548, 70–73 (2017).
  * Khajetoorians  _et al._ (2019) Alexander A Khajetoorians, Daniel Wegner, Alexander F Otte, and Ingmar Swart, “Creating designer quantum states of matter atom-by-atom,” Nature Reviews Physics 1, 703–715 (2019).
  * Kitaev (2006) Alexei Kitaev, “Anyons in an exactly solved model and beyond,” Annals of Physics 321, 2–111 (2006).
  * Wei  _et al._ (2011) Tzu-Chieh Wei, Ian Affleck, and Robert Raussendorf, “Affleck-kennedy-lieb-tasaki state on a honeycomb lattice is a universal quantum computational resource,” Phys. Rev. Lett. 106, 070501 (2011).
  * Heinrich  _et al._ (2021) Andreas J Heinrich, William D Oliver, Lieven MK Vandersypen, Arzhang Ardavan, Roberta Sessoli, Daniel Loss, Ania Bleszynski Jayich, Joaquin Fernandez-Rossier, Arne Laucht, and Andrea Morello, “Quantum-coherent nanoscience,” Nature nanotechnology 16, 1318–1329 (2021).
  * de Oteyza and Frederiksen (2022) Dimas G de Oteyza and Thomas Frederiksen, “Carbon-based nanostructures as a versatile platform for tunable $\pi$-magnetism,” Journal of Physics: Condensed Matter 34, 443001 (2022).
  * Mishra  _et al._ (2020) Shantanu Mishra, Doreen Beyer, Kristjan Eimre, Ricardo Ortiz, Joaquín Fernández-Rossier, Reinhard Berger, Oliver Gröning, Carlo A Pignedoli, Roman Fasel, Xinliang Feng, _et al._ , “Collective all-carbon magnetism in triangulene dimers,” Angewandte Chemie International Edition (2020).
  * Krane  _et al._ (2023) N Krane, E Turco, A Bernhardt, D Jacob, G Gandus, D Passerone, M Luisier, M Juríček, Roman Fasel, J Fernández-Rossier, _et al._ , “Exchange interactions and intermolecular hybridization in a spin-1/2 nanographene dimer,” Nano letters (2023).
  * Du  _et al._ (2023) Qingyang Du, Xuelei Su, Yufeng Liu, Yashi Jiang, Can Li, KaKing Yan, Ricardo Ortiz, Thomas Frederiksen, Shiyong Wang, and Ping Yu, “Orbital-symmetry effects on magnetic exchange in open-shell nanographenes,” Nature Communications 14, 4802 (2023).
  * Hieulle  _et al._ (2021) Jeremy Hieulle, Silvia Castro, Niklas Friedrich, Alessio Vegliante, Francisco Romero Lara, Sofía Sanz, Dulce Rey, Martina Corso, Thomas Frederiksen, Jose Ignacio Pascual, _et al._ , “On-surface synthesis and collective spin excitations of a triangulene-based nanostar,” Angewandte Chemie International Edition 60, 25224–25229 (2021).
  * Mishra  _et al._ (2021) Shantanu Mishra, Gonçalo Catarina, Fupeng Wu, Ricardo Ortiz, David Jacob, Kristjan Eimre, Ji Ma, Carlo A Pignedoli, Xinliang Feng, Pascal Ruffieux, Joaquin Fernandez-Rossier, and Roman Fasel, “Observation of fractional edge excitations in nanographene spin chains,” Nature 598, 287–292 (2021).
  * Delgado  _et al._ (2023) Aidan Delgado, Carolin Dusold, Jingwei Jiang, Adam Cronin, Steven G Louie, and Felix R Fischer, “Evidence for excitonic insulator ground state in triangulene kagome lattice,” arXiv preprint arXiv:2301.06171 (2023).
  * Catarina  _et al._ (2023) G Catarina, JCG Henriques, A Molina-Sánchez, AT Costa, and J Fernández-Rossier, “Broken-symmetry magnetic phases in two-dimensional triangulene crystals,” arXiv preprint arXiv:2306.17153 (2023).
  * Ortiz (2023) Ricardo Ortiz, “Olympicene radicals as building blocks of two-dimensional anisotropic networks,” arXiv preprint arXiv:2306.05346 (2023).
  * Ortiz  _et al._ (2023) Ricardo Ortiz, Geza Giedke, and Thomas Frederiksen, “Magnetic frustration and fractionalization in oligo (indenoindenes),” Physical Review B 107, L100416 (2023).
  * Cao  _et al._ (2018a) Yuan Cao, Valla Fatemi, Ahmet Demir, Shiang Fang, Spencer L Tomarken, Jason Y Luo, Javier D Sanchez-Yamagishi, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, _et al._ , “Correlated insulator behaviour at half-filling in magic-angle graphene superlattices,” Nature 556, 80–84 (2018a).
  * Cao  _et al._ (2018b) Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watanabe, Takashi Taniguchi, Efthimios Kaxiras, and Pablo Jarillo-Herrero, “Unconventional superconductivity in magic-angle graphene superlattices,” Nature 556, 43–50 (2018b).
  * Affleck  _et al._ (1987) Ian Affleck, Tom Kennedy, Elliott H. Lieb, and Hal Tasaki, “Rigorous results on valence-bond ground states in antiferromagnets,” Phys. Rev. Lett. 59, 799–802 (1987).
  * (23) Chenxiao Zhao, “Engineering interacting topological phases in nanographene-based spin-1/2 alternating heisenberg chains,” In preparation .
  * Ganesh  _et al._ (2011) R Ganesh, DN Sheng, Young-June Kim, and A Paramekanti, “Quantum paramagnetic ground states on the honeycomb lattice and field-induced néel order,” Physical Review B 83, 144414 (2011).
  * Reid (1965) D-H_ Reid, “The chemistry of the phenalenes,” Quarterly Reviews, Chemical Society 19, 274–302 (1965).
  * Jacob and Fernández-Rossier (2022) David Jacob and Joaquín Fernández-Rossier, “Theory of intermolecular exchange in coupled spin-1 2 nanographenes,” Physical Review B 106, 205405 (2022).
  * Fernández-Rossier and Palacios (2007) Joaquín Fernández-Rossier and J_ J Palacios, “Magnetism in graphene nanoislands,” Physical Review Letters 99, 177204 (2007).
  * Ortiz  _et al._ (2019) Ricardo Ortiz, Roberto Álvarez Boto, Noel García-Martínez, Juan Carlos Sancho-García, Manuel Melle-Franco, and Joaquín Fernández-Rossier, “Exchange rules for diradical $\pi$-conjugated hydrocarbons,” Nano Lett. 19, 5991–5997 (2019).
  * Arovas  _et al._ (2022) Daniel P. Arovas, Erez Berg, Steven A. Kivelson, and Srinivas Raghu, “The Hubbard model,” Annu. Rev. Condens. Matter Phys. 13, 239–274 (2022).
  * Ortiz  _et al._ (2022) R Ortiz, G Catarina, and J Fernández-Rossier, “Theory of triangulene two-dimensional crystals,” 2D Materials 10, 015015 (2022).
  * Sutherland (1986) Bill Sutherland, “Localization of electronic wave functions due to local topology,” Physical Review B 34, 5208 (1986).
  * Lieb (1989) Elliott H Lieb, “Two theorems on the hubbard model,” Physical review letters 62, 1201 (1989).
  * Ovchinnikov (1978) Alexandr A Ovchinnikov, “Multiplicity of the ground state of large alternant organic molecules with conjugated bonds,” Theoretica Chimica Acta 47, 297–304 (1978).
  * Turco  _et al._ (2023) Elia Turco, Annika Bernhardt, Nils Krane, Leosš Valenta, Roman Fasel, Michal Juriícček, and Pascal Ruffieux, “Observation of the magnetic ground state of the two smallest triangular nanographenes,” JACS Au (2023).
  * Su  _et al._ (2019) Jie Su, Mykola Telychko, Pan Hu, Gennevieve Macam, Pingo Mutombo, Hejian Zhang, Yang Bao, Fang Cheng, Zhi-Quan Huang, Zhizhan Qiu, _et al._ , “Atomically precise bottom-up synthesis of $\pi$-extended [5] triangulene,” Science advances 5, eaav7717 (2019).
  * Ahmed and Mandal (2022) Jasimuddin Ahmed and Swadhin K Mandal, “Phenalenyl radical: Smallest polycyclic odd alternant hydrocarbon present in the graphene sheet,” Chemical Reviews 122, 11369–11431 (2022).
  * (37) See Supplemental Material at [URL will be inserted by publisher].
  * Owerre (2017) SA Owerre, “Topological magnon bands in ferromagnetic star lattice,” Journal of Physics: Condensed Matter 29, 185801 (2017).
  * Yang  _et al._ (2010) Bohm-Jung Yang, Arun Paramekanti, and Yong Baek Kim, “Competing quantum paramagnetic ground states of the heisenberg antiferromagnet on the star lattice,” Physical Review B 81, 134418 (2010).
  * Jahromi and Orús (2018) Saeed S Jahromi and Román Orús, “Spin-1 2 heisenberg antiferromagnet on the star lattice: Competing valence-bond-solid phases studied by means of tensor networks,” Physical Review B 98, 155108 (2018).
  * Sachdev and Bhatt (1990) Subir Sachdev and RN Bhatt, “Bond-operator representation of quantum spins: Mean-field theory of frustrated quantum heisenberg antiferromagnets,” Physical Review B 41, 9323 (1990).
  * Gopalan  _et al._ (1994) Sudha Gopalan, TM Rice, and M Sigrist, “Spin ladders with spin gaps: A description of a class of cuprates,” Physical Review B 49, 8901 (1994).
  * Normand and Rice (1996) B Normand and TM Rice, “Electronic and magnetic structure of lacuo 2.5,” Physical Review B 54, 7180 (1996).
  * Pires (2021) Antonio Sergio Teixeira Pires,  _Theoretical tools for spin models in magnetic systems_ (IOP Publishing, 2021).