# Symmetry-Preserving Coupling Method for Topological Acoustic Metamaterials

Ssu-Ying Chen<EMAIL_ADDRESS>Department of Physics, New Jersey Institute of
Technology, Newark, NJ, USA  Camelia Prodan<EMAIL_ADDRESS>Department of
Physics, New Jersey Institute of Technology, Newark, NJ, USA

###### Abstract

In this paper we investigate different types of couplings used in acoustic
metamaterials requiring preservation of symmetries. For testing we use the SSH
model to test whether topologically edge and interface modes are supported
with the different types of connection. We observed that a modular platform
where the resonators are coupled through the bottom is the simplest method
that is accurate and flexible.

##  I I. INTRODUCTION

Topology opened a gate to realize new mechanical and acoustic systems. In
recent years, phononic crystals have grasped the attention for offering
inspiring opportunities to manipulate sound in new and unanticipated ways
based on topological concepts.

Periodic acoustic systems, such as the topological boundary states based on
the analogue of quantum Hall effect[5, 6], the analogue of quantum spin Hall
effect[7], the Floquet topological insulator[8], and the valley Hall effect[9]
have been successively proposed and experimentally verified. Besides periodic
systems, topological systems that lack periodicity can also achieve
topological edge states, such as topological quasi-crystals [10, 11, 12] where
the topological structure is caused by disorder. The idea of topological
acoustic insulators provides new schemes for designing devices with advanced
functionalities. For example, the potential improvement of leaky-wave acoustic
antenna[14], directional topological acoustic antenna controlling sound for
versatile applications[15], specific signal filtering achieved by adding
random disorder to clean structures. [16], various sound proof strategies[17,
18, 19], and for a growing discussion of different fractal geometries.[20, 21,
22] In essence, researchers have been and continue to intensely explored
various topological acoustic systems, and among these researches, topological
acoustic meta-material consisting of numerous phononic crystals is one of the
focuses.

Figure 1: Simulations showing that the bottom connection preserves the
symmetry. a Dimensions of the resonators and coupling bridges. H = 40 mm, r =
5 mm. The coupling bridge has the width = 5 mm, the length d = 26 mm, the
thickness t = 3 mm. b The height position of the coupling bridge h1 = 5 mm for
side connection, h1 = 0 mm for bottom connection. c Dimer structure and the
band spectrum of both side and bottom coupling methods. d Trimer structure and
the band spectrum of both side and bottom coupling methods. e Pentamer
structure and the band spectrum of both side and bottom coupling methods. It
is clear in the band spectrum that in all three cases, the symmetry is
preserved when the resonators are connected through bottom. Red dash lines
indicates the symmetrical axis.

In view of the fact that topological insulators are also called symmetry-
protected topological phases of matter, these gapped phases have topological
properties relying on the presence of symmetries. This result is from a global
property of the topological insulator’s band structure: local perturbations
cannot alter or damage this surface state. Since the properties of topological
insulators and their surface states highly depend on the dimension of the
material and its symmetries, they can be categorized using the periodic table
of topological insulators[13]. From the experimental perspective, it is
crucial to satisfy certain structures and to display the symmetry, moreover,
preserving the symmetry is equivalent to preserving the topological
properties. To preserve symmetry, coupling methods for topological acoustic
experiments must be carefully chosen. The following question then arises: in
discrete acoustic resonant models, how can the resonators be efficiently
coupled to form a desired structure while preserving the symmetry so that it
complies with the idea of topological metamaterial?

Some methods preserve the symmetry using double side connection, meaning they
added a coupling bridge that is also on the side and made the structure
symmetric[23, 24, 25, 26]. The resonators are 3D printed using photosensitive
resins or other types of materials as coupling bridges must be fabricated with
the resonators which makes as manufacture more complicated and time consuming.
On top of that, if the dimension of the coupling bridge has to be changed,
everything has to be made from scratch all over again.

In this article, we test an acoustic coupling method to connect the resonators
through the bottom to preserve the symmetry and create a simple, flexible and
Lego-like platform. We will compare two types of coupling methods throughout
the paper: side coupling and bottom coupling(Fig. 1b). We will start by
showing that bottom coupling preserves symmetry in simple periodic acoustic
models while single side coupling fails to do so. Via the simulations of
resonant modes of dimer, trimer, pentamer, and a classic Su–Schrieffer–Heeger
model (SSH model) consisting of 14 resonators, as well as another SSH model
with 28 resonators and a domain boundary. Furthermore, we address the
corresponding experimental results of the dimer, the SSH model, and the SSH
model with a domain boundary. All experimental results clearly show great
agreement to the simulations, proving the effectiveness and simplicity of
bottom coupling method.

##  II II. RESULTS and DISCUSSION

The dimensions of the resonators that were used throughout the simulations and
experiments are shown in (Fig. 1a,b).

###  II.1 Side Connection vs. Bottom Connection

To demonstrate and to compare the effect of side connection and bottom
connection for topological acoustic resonators, numerical simulations were
done using COMSOL Multiphysics software.

The simulated results is reported in Fig. 1 where the structures of dimer,
trimer and pentamer of resonators with both side bridge coupling and bottom
bridge coupling (Fig. 1b). The band spectrum (Fig. 1c,d,e) were generated by
sweeping the width of coupling bridge from 1 mm to 20 mm with 1 mm steps and
plotting corresponding eigenfrequencies versus widths. The wider coupling
bridge means stronger coupling strength. In the band spectrum, the second mode
split and the gap became larger when the coupling strength went up. With red
dash lines as symmetry axes, it is clear that with side connection, the
splitting modes are not symmetrical. Once we moved the coupling bridges to the
bottom, the symmetry was restored.

###  II.2 Dimer experiment

The experiment of dimer coupling was done to confirm the simulation results.
The assembly process is shown in Fig. 2a, and the experimental set up is shown
in Fig. 2c. The speaker and the microphone were placed on top of the same
resonator to give off and collect sound, respectively. The frequency of the
input signal sent to the speaker was swept from 4 kHz to 5 kHz in intervals of
10 Hz. Experimental results of dimer with bottom coupling bridges show a great
agreement to the simulations (Fig. 2b). The middle of 2 peaks is about 4.3
kHz, same as where the 2 modes start to split in the band spectrum (around
18.5 $kHz^{2}$). However, with the side connection, the peaks shifted from the
spectrum. Additionally, the height difference in peaks for bottom coupling is
noticeably smaller than that of the side coupling (Fig. 2b). Comparing with
the results from bottom coupling method, the symmetry was distinctly absent
for side coupling.

Figure 2: Dimer experimental setups as well as the comparison of simulation
and experiment results. a Assembly process of dimer connected through side.
The dimer was 3D-printed together with a side coupling bridge, the height
position of the top of the side coupling bridge h1 is equal to 7 mm. b
Assembly process of dimer connected through bottom. Two resonators were
3D-printed separately, and the bottom coupling bridge is grooved in the middle
acrylic sheet as depicted. Red boxes in both band spectrum indicate the
coupling width (5 mm) used in the experiment. c The experimental setup.

###  II.3 SSH acoustic model

To further explore how the position of connection between acoustic crystals
would influence topological gaps in band spectrum, we started by simulating 3
types of connection for SSH model of 14 resonators. The first type is to
couple them with bridges connected through side, the second one is double side
coupling bridges, and the third one is bottom connection. Fig. 3a shows the
top view of SSH coupling bridges. $r_{\textrm{1}}$ and $r_{\textrm{2}}$ are
the widths of coupling bridges, therefore they are related to the alternate
coupling strength. Light blue dash lines indicate where the resonators were
placed. The geometries of 3 coupling types are shown in Fig. 3b, c, d. The
resonant spectrum were generated by sweeping S – half of the difference in
widths between strong and weak coupling bridge – from -9 mm to 9 mm. When S
equals to zero, all coupling bridges have an identical width of 10 mm.

Figure 3: Band spectrum of different types of connection in SSH model of 14
resonators, sweeping S. S is half of the difference in widths between strong
and weak coupling bridge. a Top view of SSH coupling bridges. $r_{\textrm{1}}$
and $r_{\textrm{2}}$ are related to the alternate coupling strength. Light
blue dash lines are where the resonators were placed. b SSH model with bottom
connection. c SSH model with double connection from the sides. b and c has
very similar spectrum. Blue boxes includes edge resonant modes. d SSH model
with side connection, the edge band merges into bulk. e Band spectrum of SSH
model sweeping the height of the coupling bridges. Red parts represent when
the coupling bridges protrude from the bottom or top. Dark blue dash line
marks 18.5 $kHz^{2}$.

As one can see in the spectrum (Fig. 3b, c, d, there is no gap for a uniform
coupling connection (S = 0 mm, $r_{\textrm{1}}$ = $r_{\textrm{2}}$ = 10 mm). A
bulk spectral gap opens once the connecting channels are set in an alternating
strong/weak coupling strength ($S\neq 0~{}mm$). Furthermore, the bulk spectrum
remains symmetric with respect to the middle of the bulk spectral gap (the
small deviations are less than 5% when compared with the overall width of the
bulk spectrum).

In (Fig. 3b), SSH model with bottom connection, one can observe the expected
edge resonant modes, whose energies are pinned in the middle of the bulk
spectral gap around 18.5 $kHz^{2}$, where the second mode of resonant
frequency of a single resonator is. We can say that SSH model with bottom
connection, COMSOL-simulated spectra displays an exact chiral symmetry. In
Fig. 3c, SSH model with double side coupling bridges gives the same spectra,
meaning both types of coupling approaches achieve the goal of protecting the
chiral symmetry.

On the other hand, the spectra of SSH model connected through sides does not
display the same symmetry. In Fig. 3d, not only do the edge modes merge into
lower bulk bands, but the middle point of the spectrum shifted to around 2
$kHz^{2}$, demonstrating the same impracticality as in dimer experiments.

We also checked to see if the height position (h1 in Fig. 1b) of side coupling
plays a role in preserving the symmetry. The height of side coupling bridge
was swept from bottom to top to generate the spectra in Fig. 3e, while the
coupling width were fixed at $r_{\textrm{1}}=15$ and $r_{\textrm{2}}=5$ (S =
-5 mm). When h1 = -21.5 and 21.5 mm, the resonators are coupled from the
bottom or from the top which resulted in bulk modes symmetric with respect to
18.5 $kHz^{2}$, and the edge modes appear at middle. Red parts represent when
the coupling bridges protrude from the bottom or the top, meaning part of them
are outside of the resonators, there are still edge modes existing in the gap
but the symmetry shifts further and further away from the dash line. One can
see that once the coupling bridges enter the body of resonators, the edge
modes merge into bulk and the symmetry disappears.

Next, we did experiments with SSH model connected through bottom to confirm if
the symmetry is also seen in the experimental results. Fig. 4a shows how the
setup was made. The COMSOL-simulated acoustic pressure fields of the edge
resonant modes are shown in Fig. 4b. The simulated band structure in Fig. 4c
is added as a reference (same as Fig. 3b).

In these measurements, same as dimer experiments, the speaker and the
microphone were inserted in the same resonator via two holes open at the top
and the speaker’s frequency was swept from 4 kHz to 5 kHz and the microphone
picked up corresponding signals. The measurements were repeated for all
resonators and the collected data was assembled in the local density of states
plot shown in Fig. 4d. It is worth mentioning here that all resonators are
removable and interchangeable, so the one with holes can be placed at any
probe position desired for measurements. Panel e provides an alternative
depiction of the same data. The spectral gap as well as the expected edge and
interface modes can be clearly identified and they are well aligned with the
simulation in panel c.

The density of states reported in Fig. 4d was obtained by integrating the
local density of states acquired from resonators whose index are the same as
the position of the probe. The same instrumentation was used. The measurements
were repeated while moving the position of the probe. For each measurement,
the frequency was scanned from 4 kHz to 5 kHz in 20 Hz steps.

Figure 4: SSH experiments demonstrate the simplicity and efficiency of bottom
coupling method. a The assembling process of SSH model of 14 resonators
connected through bottom. b Acoustic pressure field distribution for the edge
modes marked as red dot in panel c when S = -5 mm. c COMSOL simulated SSH
model resonant spectrum. The red vertical box indicates S = -5 mm which makes
$r_{\textrm{1}}=15$ mm and $r_{\textrm{2}}=5$ mm, the parameters were used in
the experiments. d Experimentally measured local density of states, assembled
from normalized microphone readings from the top of the block resonators. The
bright dispersive modes indicates the bulk and edge modes. e Collapse on the
frequency axis of the intensity plot in d. The spectral gap is clearly
recognized and the edge modes that show up in the gap are marked with a red
star.

Finally, we test and compare two coupling methods in a SSH model with a domain
boundary(interface), both simulationally and experimentally. A SSH model with
a domain boundary separates two topologically distinct SSH insulating phases,
with one non-trivial edge(left) and one trivial(right). It is shown in Fig. 5a
where $r_{\textrm{1}}=15$ mm and $r_{\textrm{2}}=5$ mm. Fig. 5b and Fig. 5c
show the band spectrum of side connection and bottom connection, respectively.
Red vertical boxes include resonant modes when S = -5 mm. Red and blue stars
label where interface mode and edge modes appear, and the acoustic pressure
field maps are shown below. Similar to the previous results, with bottom
coupling, both simulation and experiment verified that the system contains
topological resonant modes at the non-trivial edge as well as the domain
boundary as expected, and they are located in the middle of the bulk band gap.
In opposition, with side coupling, the edge mode disappears and the interface
mode is close to the bulk, besides, it is noticeable that the energy is not as
concentrated at interface from the acoustic pressure field distribution.

We also did experiments for bottom-connected SSH model with a domain boundary,
the same protocol was applied and frequency was swept from 4 kHz to 5 kHz. The
measurements were repeated for all 28 resonators and the collected data was
then used to plot the local density of states which is shown in Fig. 5d. Panel
e shows the collapse of the data in panel d on the frequency axis. The
spectral gap as well as the expected edge modes can be clearly identified and
are well aligned with the simulation in panel c. The COMSOL-simulated acoustic
pressure field maps of the edge resonant modes are shown in Fig. 5b, c.

##  III III. CONCLUSION

In this paper we investigate acoustic coupling that preserve symmetries
through Comsol simulations and experiments. We tested couplings through the
bottom, and with bridges at different height through the side. The first test
was done on dimmers, then on single SSH set up, as well as one with interface.
We observed that the simplest way for such coupling is to connect the
resonators through the bottom. Coupling through the side requires a second
connection. The advantages of coupling through the bottom is the modular
structure of the set up. The resonators and coupling are printed individually,
allowing for an easy variation in the coupling strength when the experiment
requires. This platform has an exceptional flexibility since the resonators
can be stored for further use.

##  IV IV. MATERIAL and METHODS

###  IV.1 Simulation

The simulations reported in all figures were performed with the COMSOL
Multiphysics pressure acoustic module. The wave propagation domain shown in
Fig. 1 was filled with air with a mass density 1.3 kg/m3 and the sound’s speed
was set at 343 m/s, which is appropriate for room temperature. We shall
consider the 3D printing UV resin material as hard boundary because of the
huge acoustic impedance mismatch comparing with air.

###  IV.2 Experiment

The resonators were 3D-printed using Anycubic Photon 3D printer, which uses UV
resin and has 47 um XY-resolution and 10 um Z-resolution. The thickness of the
walls is 2 mm, which ensures a high Q factor and justifies the rigid
boundaries in the simulations. The inner dimensions of the resonators are
shown in Fig. 1a.

A dimer with a coupling bridge on the sides is 3D-printed as a whole. The
width, length and the position of the coupling bridge are as labeled in Fig.
1a. One side was left open for ethanol rinsing and UV-curing. The dimer was
then placed on a base of two layers of acrylic plates (top layer: 2 mm in
thickness, bottom layer: 3 mm in thickness) to create a closed space for wave
propagating(Fig. 2a). The reason that the top layer was 2 mm is to accommodate
the side bridge. The bottom connection is achieved by assembling the
supporting base, which consists of three layers of 3 mm thick acrylic
plates(Fig. 2b). The middle layer of a groove is to account for acoustic
coupling. The acrylic plates with patterns of the supporting bases were cut by
the Boss Laser-1630 Laser Engraver. The nominal tolerance of the laser-cutter
is 250 um.

Figure 5: The SSH model consists of 28 resonators with a domain boundary in
the middle, topologically non-trivial at one end and trivial at the other. a
Coupling bridges and its dimension. $r_{\textrm{1}}=5$ mm and
$r_{\textrm{2}}=15$ mm. b SSH model with side connection. The acoustic
pressure field map below shows the interface mode at the red star in the
spectrum when $S=-5$ mm. The model lacks edge modes. c SSH model with bottom
connection. Band structure and acoustic pressure field maps of both edge mode
and interface mode. Blue and Red stars in the panel mark edge mode and
interface mode, respectively. d Experimentally measured local density of
states, assembled from normalized microphone readings from the top of the
block resonators. The bright dispersive modes indicates the bulk and edge
modes. e Collapse on the frequency axis of the intensity plot in d. The
spectral gap can be clearly identified and the edge and interface modes show
up in the gap marked with a red star.

For the SSH model connected from bottom, the same method was utilized. 14
resonators were placed and coupled through the channels with alternating
widths grooved in the acrylic plates of the base. These resonators are
detachable and interchangeable so that they can be moved around, thus acoustic
crystals with different probe positions can be generated, and the resonators
can be taken apart, stored and reassembled for new projects or designs.

The protocol for the acoustic measurements shown in Fig. 2, Fig. 4 and Fig. 5
was as follows: Sinusoidal signals of duration 1 s and amplitude of 0.5 V
generated by a Rigol DG 1022 function generator were sent out to a speaker
placed in a porthole opened on top of a resonator. A dbx RTA-M Reference
Microphone with a Phantom Power was inserted in a porthole next to the
previous one and was used to acquire the acoustic signals (Fig. 2c). The
signals were then read by a custom LabVIEW code via National Instruments
USB-6122 data acquisition box and the data was stored for graphic renderings.

##  V ACKNOWLEDGMENT

The authors acknowledges support from the National Science Foundation, grant
CMMI-2131759.

## References

  * [1] Intro: topological acoustics
  * [2] Hasan, M Zahid and Kane, Charles L Colloquium: topological insulators, Reviews of modern physics 82, 3045 (2010).
  * [3] Ozawa, Tomoki et al. Topological photonics, Reviews of modern physics 91, 015006 (2019).
  * [4] Zhang, Xiujuan et al. Topological sound, Communications Physics 1, 1 (2018).
  * [5] Khanikaev, Alexander B et al. Topologically robust sound propagation in an angular-momentum-biased graphene-like resonator lattice, Nature communications 6, 1 (2015).
  * [6] Yang, Zhaoju et al. Topological acoustics, Physical review letters 114, 114301 (2015).
  * [7] He, Cheng et al. Acoustic topological insulator and robust one-way sound transport, Nature physics 12, 1124 (2016).
  * [8] Peng, Yu-Gui et al. Experimental demonstration of anomalous Floquet topological insulator for sound, Nature communications 7, 13368 (2016).
  * [9] Lu, Jiuyang et al. Valley vortex states in sonic crystals, Physical review letters 116, 093901 (2016).
  * [10] Ni, Xiang et al. Observation of Hofstadter butterfly and topological edge states in reconfigurable quasi-periodic acoustic crystals, Communications Physics 2, 55 (2019).
  * [11] Cheng, Wenting and Prodan, Emil and Prodan, Camelia Experimental Demonstration of Dynamic Topological Pumping across Incommensurate Bilayered Acoustic Metamaterials, Physical Review Letters 125, 1224301 (2020).
  * [12] Apigo, David J et al. Observation of topological edge modes in a quasiperiodic acoustic waveguide, Physical Review Letters 122, 095501 (2019).
  * [13] Chiu, Ching-Kai et al. Classification of topological quantum matter with symmetries, Reviews of Modern Physics 88. 035005 (2016).
  * [14] Fleury et al. Floquet topological insulators for sound, Nature Communications 7, 11744 (2016).
  * [15] Zhang et al. Directional Acoustic Antennas Based on Valley-Hall Topological Insulators, Advanced Material 30, 1803229 (2018).
  * [16] Zangeneh-Nejad et al. Disorder-Induced Signal Filtering with Topological Metamaterials, Advanced Material 32, 2001034 (2020).
  * [17] T. Yamamoto et al. Topology design of multi-material soundproof structures including poroelastic media to minimize sound pressure levels, Comput. Methods Appl. Mech. Engrg. 198. 1439–1455 (2009).
  * [18] Xu, Zi-xiang et al.Topology-Optimized Omnidirectional Broadband Acoustic Ventilation Barrier, Physical Review Applied 14. 054013 (2020).
  * [19] Liu, Chenkai et al. Three-Dimensional Soundproof Acoustic Metacage, Physical Review Letter 127 084301 (2021).
  * [20] Song, Gang Yong et al. Broadband Focusing Acoustic Lens Based on Fractal Metamaterials, Scientific Report 6 35929 (2016).
  * [21] Zhao, Xinzhe et al. Fractal acoustic metamaterials for transformer noise reduction, Applied Physics Letters. 113 074101 (2018).
  * [22] Singh, Sanjeet Kumar et al. Hybrid fractal acoustic metamaterials for low-frequency sound absorber based on cross mixed micro-perforated panel mounted over the fractals structure cavity, Scientific Report 12 20444 (2022).
  * [23] Zheng, Shengjie et al. Observation of fractal higher-order topological states in acoustic metamaterials, Science Bulletin 67, 2069-2075 (2022).
  * [24] Xiao, Yi-Xin et al. Topological Subspace-Induced Bound State in the Continuum, Physical Review Letters 118, 166803 (2017).
  * [25] Ma, Guancong et al. Topological phases in acoustic and mechanical systems, Nature Review Physics 1, 281 (2019).
  * [26] Qi, Yajuan et al. Acoustic Realization of Quadrupole Topological Insulators, Physical Review Letters 124, 206601 (2020).


# On The Continuous Coverage Problem for a Swarm of UAVs

Hazim Shakhatreh  Department of Electrical and Computer Engineering, New
Jersey Institute of Technology  Abdallah Khreishah  Department of Electrical
and Computer Engineering, New Jersey Institute of Technology  Jacob Chakareski
Department of Electrical and Computer Engineering, University of Alabama
Haythem Bany Salameh  Department of Telecommunications Engineering, Yarmouk
University  Issa Khalil  Qatar Computing Research Institute

###### Abstract

Unmanned aerial vehicles (UAVs) can be used to provide wireless network and
remote surveillance coverage for disaster-affected areas. During such a
situation, the UAVs need to return periodically to a charging station for
recharging, due to their limited battery capacity. We study the problem of
minimizing the number of UAVs required for a continuous coverage of a given
area, given the recharging requirement. We prove that this problem is NP-
complete. Due to its intractability, we study partitioning the coverage graph
into cycles that start at the charging station. We first characterize the
minimum number of UAVs to cover such a cycle based on the charging time, the
traveling time, and the number of subareas to be covered by the cycle. Based
on this analysis, we then develop an efficient algorithm, the cycles with
limited energy algorithm. The straightforward method to continuously cover a
given area is to split it into N subareas and cover it by N cycles using N
additional UAVs. Our simulation results examine the importance of critical
system parameters: the energy capacity of the UAVs, the number of subareas in
the covered area, and the UAV charging and traveling times. We demonstrate
that the cycles with limited energy algorithm requires 69%-94% fewer
additional UAVs relative to the straightforward method, as the energy capacity
of the UAVs is increased, and 67%-71% fewer additional UAVs, as the number of
subareas is increased.

###### Index Terms:

Unmanned aerial vehicles, charging, continuous coverage.

##  I Introduction

In 2005, Hurricane Katrina in the United States caused over 1,900 deaths, 3
million land-line phone interruptions, and more than 2,000 base stations going
out of service [1, 2, 3]. Another example of a large-scale interruption of
telecommunications service is the World Trade Center attack in 2001, when it
took just minutes for the nearby base stations to be overloaded. The attacks
caused the disturbance of a phone switch with over 200,000 lines, 20 cell
sites, and 9 TV broadcast stations [1, 4]. These incidents demonstrate the
need for a quick/efficient deployment network for emergency cases.

The authors in [5] proposed a UAV-based replacement network during disasters,
where the UAVs serve as aerial wireless base stations. However, this study did
not consider how the UAVs will guarantee a continuous coverage when they need
to return to the charging station for recharging. Though a UAV has limited
energy capacity and needs to recharge its battery before running out of power
during the coverage process, only few studies have considered this constraint
in the UAV coverage problem. Concretely, the author in [6] determined the
minimum number of UAVs that can provide continuous

Figure 1: The Continuous Coverage Problem

coverage for a single area using UAVs with uniform and non-uniform energy
capacity. However, no consideration has been made for the case when there are
multiple subareas that need to be covered, which is the typical scenario
during disasters. The authors in [7, 8] formulated the Mobile Charging
Problem, in which multiple mobile chargers collaborate to charge static
sensors with minimum number of mobile chargers subject to speed and energy
limits of the mobile chargers, such that the energy consumption of the mobile
chargers is minimized. Note that the mobile charging problem is different from
the problem we study. In the mobile charging problem, the chargers will not
cover the sensors continuously. The mobile charger will visit the sensor and
stay for a specific time to charge the sensor and after finishing the charging
process, it will visit other sensors.

In [9], the authors studied the continuous coverage problem for mobile
targets, where during the coverage process a UAV that runs out of energy is
replaced by a new one. Many studies [10, 11, 12, 13] focused on minimizing the
total transmission power of the UAVs during the coverage of a geographical
area, however, no limits on the UAV energy capacity and the need for
recharging have been considered. The work in [14] reported that the energy
consumption during data transmission and reception is much smaller than the
energy consumption during the UAV hovering, i.e., it only constitutes
10$\%$-20$\%$ of the UAV energy capacity. Thus, it is important to conduct
studies that take into account the energy consumption during the UAV hovering
rather than focusing on minimizing the energy consumption during data
transmission and reception.

Contrary to the related work above, we integrate the recharging requirements
into the coverage problem and examine the minimum number of required UAVs for
enabling continuous coverage under that setting (see Figure 1). To the best of
our knowledge, this is the first study that jointly considers the coverage and
recharging problems where multiple subareas are to be covered. Our main
contributions in this context are: 1) We formulate the problem of minimizing
the number of UAVs required to provide continuous coverage of a given area,
given the recharging requirement. 2) We prove that this problem is NP-
complete. 3) Due to the intractability of the problem, we study partitioning
the coverage graph into cycles that start at the charging station. 4) Based on
this analysis, we develop an efficient algorithm.

The rest of this paper is organized as follows. Section II presents our system
model including the problem formulation and the proof of the NP completeness.
In Section III, we show how to find the minimum number of additional UAVs that
are required to guarantee the continuous coverage. Then, we present our
proposed algorithm. In Section IV, we present our experimental results.
Finally, Section V concludes the paper.

##  II System Model

###  II-A Problem Statement

Consider a geographical area $G$={$g_{1}$,…,$g_{N}$}, where $g_{i}$ represents
a subarea $i$, the subarea $g_{1}$ $\in$ $G$ includes the charging station and
all subareas except subarea $g_{1}$ need to be covered $G_{0}$=$G$
$\backslash$ $g_{1}$. We want to find the minimum number of UAVs that can
provide a continuous coverage over $G_{0}$ by placing the UAVs at locations
where each UAV will provide full coverage for one subarea. In the continuous
coverage problem, we assume: (1) Time slotted system in which the slot
duration is 1 time unit and the total coverage duration is $T$. (2) All UAVs
start the coverage process from the charging station and they need to return
to the charging station after they complete the coverage process. (3) Each UAV
has limited energy capacity $E$ and it needs to return to the charging station
to recharge the battery before running out during the coverage process. (4)
Each UAV can move (from the charging station to location $i$), (from location
$i$ to location $j$) or (from location $j$ to the charging station) and this
process will take one time slot. (5) Each UAV covers a given subarea for one
or multiple time slots. (6) Each subarea will be covered by only one UAV. (7)
The UAV cannot travel to the charging station or to any other location until
the handoff process is completed in which another UAV arrives to cover the
subarea such that the continuous coverage is guaranteed. (8) The recharging
process takes $T_{charge}$ at the charging station.

###  II-B Problem Formulation

Now, we formulate the continuous coverage problem. In order to present the
problem formulation, we introduce the binary variable $x_{m}$ that takes the
value of 1 if the UAV $m$ visits any subarea from charging station during the
coverage duration $T$ and equals 0 otherwise; the binary variable
$y_{ij,m}^{t}$ that takes the value of 1 if the if the UAV $m$ moves through
edge $ij$ during the time slot $t$ and equals 0 otherwise; the binary variable
$z_{j,m}^{t}$ that takes the value 1 if the UAV $m$ covers the subarea $j$ at
time slot $t$ and equals 0 otherwise. Table I provides a list of the major
notations used in this paper.

TABLE I: A List of Notations $M$ | The set of fully charged UAVs available at the charging
---|---
| station.
$E$ |  The energy capacity of each UAV $m$ $\in$ $M$.
$T$ | The coverage duration.
$E_{ij}^{Travel}$ | The energy consumed when the UAV travels from subarea
|  $i$ to subarea $j$ where $i,j\in G$.
$E_{j}^{Cover}$ |  The energy consumed when the UAV covers the subarea $j$
|  for one time slot where $j\in G$ (constant).
$T_{charge}$ | The time that the UAV needs to recharge the battery at
| the charging station.

$\begin{array}[]{ll}\displaystyle\min\displaystyle\sum_{m\in M}x_{m}\\\
\textrm{subject to}\\\ y_{ij,m}^{t}\leq x_{m}~{}~{}~{}~{}~{}~{}\forall i,j\in
G,\forall t\in[0,T],\forall m\in M&(1)\\\ z_{j,m}^{t}\leq
x_{m}~{}~{}~{}~{}~{}~{}~{}\forall j\in G_{0},\forall t\in(0,T),\forall m\in
M&(2)\\\ \displaystyle\sum_{m\in M}y_{1j,m}^{0}=1~{}~{}\forall j\in
G_{0}&(3)\\\ \displaystyle\sum_{m\in M}z_{j,m}^{t}=1~{}~{}~{}\forall j\in
G_{0},\forall t\in(0,T)&(4)\\\ \displaystyle\sum_{i\in G,i\neq
j}\displaystyle\sum_{m\in M}y_{ij,m}^{t}\leq 1~{}~{}~{}~{}~{}~{}\forall j\in
G_{0},\forall t\in[0,T]&(5)\\\ \displaystyle\sum_{i_{1}\in
G}y_{i_{1}j,m_{1}}^{t}=\displaystyle\sum_{i_{2}\in
G}y_{ji_{2},m_{2}}^{t+1}~{}\forall j\in G_{0},\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
t\in(0,T),m_{1}\neq m_{2}&(6)\\\ \displaystyle\sum_{m\in
M}\displaystyle\sum_{t\in[0,T)}\displaystyle\sum_{i\in
G}y_{ij,m}^{t}\leq\displaystyle\sum_{m\in
M}\displaystyle\sum_{t\in(0,T)}z_{j,m}^{t}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in G_{0}&(7)\\\ \displaystyle\sum_{j\in G_{0}}\displaystyle\sum_{\tau\in
T_{charge}}[y_{j1,m}^{t}+y_{1j,m}^{t+\tau}]\leq 1\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
m\in M,\forall t\in(0,T)&(8)\\\ \displaystyle\sum_{m\in
M}\displaystyle\sum_{t\in[0,T]}\displaystyle\sum_{j\in
G_{0}}y_{1j,m}^{t}=\displaystyle\sum_{m\in
M}\displaystyle\sum_{t\in[0,T]}\displaystyle\sum_{i\in
G_{0}}y_{i1,m}^{t}&(9)\\\
\displaystyle\sum_{t\in[t_{1},t_{2}]}\displaystyle\sum_{i,j\in
G}E_{ij}^{Travel}y_{ij,m}^{t}+\displaystyle\sum_{t\in[t_{1},t_{2}]}\displaystyle\sum_{j\in
G}E_{j}^{Cover}z_{j,m}^{t}\\\ \leq E~{}~{}~{}~{}~{}\forall m\in
M,\forall[t_{1},t_{2}]\in[0,T],t_{1}=arg~{}y_{1j,m}^{t}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}t_{2}=arg~{}y_{i1,m}^{t},t_{2}>t_{1},\forall
t_{3}\in(t_{1},t_{2})\\\ ~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}t_{3}\neq
arg~{}y_{1j,m}^{t}\neq arg~{}y_{i1,m}^{t}&(10)\end{array}$

The objective is to minimize the number of UAVs that are needed to provide a
continuous coverage during the coverage duration $T$ subject to various design
constraints. Constraints (1) and (2) ensure that the UAV can travel and cover
the subareas only if we select it to participate in the coverage process.
Constraint (3) ensures that all subareas will be covered at the first time
slot. Constraint (4) guarantees the continuous coverage for each subarea.
Constraint (5) allows the UAV to visit a new subarea (when $y_{ij,m}^{t}$=1)
or to continue covering the current subarea (when $y_{ij,m}^{t}$=0).
Constraint (6) characterizes the handoff process between the UAVs, when the
UAV $m_{1}$ wants to visit the subarea $j$ from subarea $i_{1}$ at time $t$
($y_{i_{1}j,m_{1}}^{t}$=1), the UAV $m_{2}$ that covers the subarea $j$ will
travel to subarea $i_{2}$ at time $t+1$ ($y_{ji_{2},m_{2}}^{t+1}$=1).
Constraint (7) describes the relation between the traveling process and the
covering process, where the number of times that the subarea $j$ is covered
will be greater than or equal the number of times that it is visited.
Constraint (8) shows that the recharging process will take $T_{charge}$ at the
charging station. Constraint (9) ensures that the number of UAVs outgoing from
the charging station and the number of UAVs incoming to charging station are
the same after we complete the coverage process. Constraint (10) shows that
the energy capacity of the UAV can cover the wasted energy during the
traveling and the covering processes in each cycle where $t_{1}$ represents
the time that the UAV travels from the charging station with full energy
capacity and $t_{2}$ represents the time that the UAV arrives to the charging
station to charge the battery. Now we will prove that the continuous coverage
problem is an NP-complete.

###  II-C NP completeness

######  Theorem 1.

The Continuous Coverage Problem is NP-complete.

###### Proof.

The number of constraints is polynomial in terms of the number of subareas,
the number of UAVs and the number of time slots. Given any solution for our
problem, we can check the solution’s feasibility in polynomial time, then the
problem is NP.
To prove that the problem is NP-hard, we reduce the Bin Packing Problem which
is NP-hard [15] to a special case of our problem. The special case will be the
discrete coverage problem in which each subarea will be visited one time by
one UAV during the coverage process. In the Bin Packing Problem, we have $p$
items where each item has volume $z_{p}$. All items must be packed into a
finite number of bins ($b_{1}$, $b$,…,$b_{B}$), each of volume $V$ in a way
that minimizes the number of bins used. The reduction steps are: 1) The $b$-th
bin in the Bin Packing Problem is mapped to the $m$-th UAV in our problem
(where the volume $V$ for each bin is mapped to the energy capacity of the UAV
$E$). 2) The $p$-th item is mapped to the $n$-th subarea, (where the volume
for each item $p$ is mapped to the energy consumed when the UAV (visits and
covers) subarea $n$. 3) All UAVs have the same energy capacity $E$. 4) The
energy consumed (during the traveling and the covering processes) when the UAV
visits subarea $j$ from any subarea ($i$ $\in$ $G$ $\backslash$ {j}) will be
constant. 5) The energy required for the UAV to return to the charging station
from any subarea $i$ will be zero ($E_{i1}^{Travel}$=0). 6) The time that the
UAV needs to recharge the battery at the charging station will be infinity. 7)
Each subarea will be visited one time by one UAV during the coverage process
(discrete coverage). If there exists a solution to the bin packing problem
with cost $C$, then the selected bins will represent the UAVs that are
selected and the items in each bin will represent the subareas that the UAV
must visit and the total cost of our problem is $C$. If there exists a
solution to our problem with cost $C$, then the selected UAVs will represent
the bins that are selected and the subareas that the UAV must visit will
represent the items in each bin and the total cost of the bin packing problem
is $C$. We prove that there exists a solution to the bin packing problem with
cost $C$ iff there exists a solution to our problem with cost $C$.∎

##  III Heuristic algorithm

Due to the intractability of the problem, we study partitioning the coverage
graph into cycles that start at the charging station. We first characterize
the minimum number of UAVs to cover each cycle based on the charging time, the
traveling time, and the number of subareas to be covered by the cycle. Our
analysis based on the uniform coverage in which the UAV covers each subarea in
a given cycle for a constant time. Based on this analysis, we then develop an
efficient algorithm, the cycles with limited energy algorithm, that minimizes
the required number of UAVs that guarantees a continues coverage.

###  III-A Analysis

It is obvious that we need $N$ UAVs to cover $N$ subareas at any given time,
but the question here is how many additional UAVs are needed to guarantee a
continuous coverage. In this subsection, we assume that the UAV visits the
subareas based on a cycle that starts from the charging station and ends at
the charging station for charging process. We also assume that a given UAV
covers the subareas in the cycle uniformly, in which the UAV covers each
subarea in a given cycle for a constant time. In Theorem 2, we find the
minimum number of additional UAVs that are needed to guarantee a continuous
coverage for a cycle, which will help us while developing Algorithm 1.

######  Theorem 2.

The minimum number of additional UAVs $k$ that are required to provide
continuous and uniform coverage for a cycle that contains $n$ subareas must
satisfy this inequality:

$k\frac{T_{Coverage}}{n}\geq(n+1)T+T_{Charge}$ (1)

where $T_{Coverage}$ is the time that the UAV allocates to cover all subareas
in the cycle, $T$ is the time that the UAV needs to travel from subarea $i$ to
subarea $j$ and $T_{charge}$ is the time that the UAV needs to recharge the
battery at the charging station.

###### Proof.

Consider that all $n$ subareas in the cycle are covered by $n$ UAVs and the
UAV that covers the last subarea want to return to the charging station to
recharge its battery. The handoff process needs to begin between one of the
additional UAVs from the charging station and the UAV that covers the first
subarea in the cycle.
The UAV that covers the last subarea needs to wait $(n-1)$ $T$ to do the
handoff process, during this time the additional UAVs are covering the first
subarea. After the handoff process is completed, the UAV needs $T$ time units
to return to the charging station, $T_{charge}$ to recharge the battery and T
to visit the first subarea in the cycle again. Then, we have:

$k\frac{T_{Coverage}}{n}\geq(n-1)T+T+T_{Charge}+T$ (2)

∎

###  III-B The cycles with limited energy algorithm

The straightforward method (SM) to continuously cover $N$ subareas is to
allocate two UAVs for each subarea. At the first time slot, $N$ UAVs cover the
$N$ subareas. Then, any UAV wants to return to the charging station to
recharge the battery will do the handoff process with one of the additional
UAVs that are available at the charging station. By applying SM, we need $N$
additional UAVs and we have $N$ cycles to cover all the subareas.
Our proposed algorithm, the cycles with limited energy algorithm (CLE), is
inspired by the nearest neighbor algorithm, the nearest neighbor algorithm is
used to solve the Traveling Salesman Problem [16], in which the salesman keeps
visiting the nearest unvisited vertex until all the vertices are visited. In
our algorithm, the UAV (salesman) has limited energy capacity and before
visiting any new subarea, we must check if the remaining energy is enough to
return to the charging station from the new location or not. In the previous
subsection, we show how to find the minimum number of additional UAVs that are
required to guarantee the continuous coverage for a given cycle, we use the
Theorem 2 to find the minimum number of additional UAVs that are required to
provide the continuous coverage for a given area, by finding the cycles that
need only one additional UAV. The pseudo code of this algorithm is shown in
Algorithm 1.

Input:

The geografical area $G$={$g_{1}$,…,$g_{N}$},

$T$: The required time to travel between two subareas in the area,

$E$: The energy capacity of the UAV,

$T_{Charge}$: The time that the UAV needs to recharge the battery at charging
station,

$e$: The energy consumed by the UAV when it covers the subarea for 1 sec.

$i$=1 // $i$ is the index of the cycle

START:

while $G$ not empty do

$c_{i}$={$g_{1}$}

Do:

1 $v$= most recently added subarea to cycle $c_{i}$

2 Find {$g$}= $argmin_{b\in G-\\{v\\}}$ $distance(v,b)$

3 Calculate $E_{Coverage}$=$E$-$E_{Travel}$-$E_{ReturntoBS}$

4 Calculate $T_{Coverage}$ $=$ $\frac{E_{Coverage}}{e}$

5 If $\frac{T_{Coverage}}{|c_{i}|}$ $\geq$ $(|c_{i}|+1)T$ $+$ $T_{Charge}$
then

6 $c_{i}$ $\longleftarrow$ $c_{i}$ $\cup$ $\\{g\\}$

7 $G$ $\longleftarrow$ $G$ $\backslash$ $\\{g\\}$

while ($\frac{T_{Coverage}}{|c_{i}|}$ $\geq$ $(|c_{i}|+1)T$ $+$ $T_{Charge}$ )

8 $c_{i}$ $\longleftarrow$ $c_{i}$ $\cup$ $\\{g_{1}\\}$

9 C $\longleftarrow$ C $\cup$ $c_{i}$

10 $i$=$i$+1

end while

Output: C and $M^{{}^{\prime}}$ $=$( (
$\displaystyle\sum_{c_{i}\in\textbf{C}}|c_{i}|-2))+|\textbf{C}|$

Algorithm 1 The cycles with limited energy algorithm

##  IV Performance Evaluation

###  IV-A Power Consumption Models

In this section, we quantify the power consumption by each UAV when it is
hovering, traveling and transmitting data.

####  IV-A1 Power consumption during hovering

The power consumption in watt by the UAV during hovering can be given by [17]:

$P=4\frac{T_{h}^{3/2}}{\sqrt{2QS}}+p$ (3)

where $T_{h}$ is the fourth of the quadcopter total weight in $N$, $Q$ is the
density of the air in kg/$m^{3}$, $S$ is the rotor swept area in $m^{2}$ and
$p$ is the power consumption of electronics in watt.

####  IV-A2 Power consumption during traveling

The power consumption in kW by the UAV during traveling can be given by [18]:

$P=\frac{(m_{p}+m_{v})v}{370\eta r}+p$ (4)

where $m_{p}$ is the payload mass in kg, $m_{v}$ is the vehicle mass in kg,
$r$ is the lift-to-drag ratio (equals 3 for vehicle that is capable of
vertical takeoff and landing), $\eta$ is the power transfer efficiency for
motor and propeller, $p$ is the power consumption of electronics in kW and $v$
is the velocity in km/h.

####  IV-A3 Power consumption during data transmission

The power consumption in dB by the UAV during data transmission can be given
by [10]:

$P_{t}(dB)=P_{r}(dB)+\bar{L}(R,h)$ (5) $\bar{L}(R,h)=P(LOS)\times
L_{LOS}+P(NLOS)\times L_{NLOS}$ (6)
$P(LOS)=\dfrac{1}{1+\alpha.exp(-\beta[\frac{180}{\pi}\theta-\alpha])}$ (7)
$L_{LOS}(dB)=20log(\dfrac{4\pi f_{c}d}{c})+\xi_{LOS}$ (8)
$L_{NLOS}(dB)=20log(\dfrac{4\pi f_{c}d}{c})+\xi_{NLOS}$ (9)

In equation (5), $P_{t}$ is the transmit power, $P_{r}$ is the required
received power to achieve a SNR greater than threshold $\gamma_{th}$,
$\bar{L}(R,h)$ is the average path loss as a function of the altitude $h$ and
coverage radius $R$.

In equation (6), $P(LOS)$ is the probability of having line of sight (LOS)
connection at an evaluation angle of $\theta$, $P(NLOS)$ is the probability of
having non LOS connection and equal (1-$P(LOS)$), $L_{LOS}$ and $L_{NLOS}$ are
the average path loss for LOS and NLOS paths.

In equations (7), (8) and (9), $\alpha$ and $\beta$ are constant values which
depend on the environment, $f_{c}$ is the carrier frequency, $d$ is the
distance between the UAV and the user, $c$ is the speed of the light ,
$\xi_{LOS}$ and $\xi_{NLOS}$ are the average additional loss which depend on
the environment.

Actually, the power consumed by the UAV during data transmission and reception
is much smaller than the power consumed during the UAV hovering or traveling
[14]. In this paper, we assume that the power wasted during data transmission
is constant.

###  IV-B Simulation Setup

Given a geographical area $G$ , the number of the subareas that we need to
cover and the density of the users, the question here is how to find the
optimal boundaries of the subareas that to be covered by the UAVs. To answer
this question, the authors of [11] used transport theory to find the optimal
boundaries of the subareas. Unfortunately, this approach needs to solve
${N\choose 2}$ non-linear equations at each iteration, where $N$ is the number
of subareas. In this paper, we divide the geographical area uniformly and
apply the SM and CLE algorithm to find the minimum number of additional UAVs
that provides the continuous coverage. We study the effect of the UAV energy
capacity, the grid size of the geographical area, the charging time and the
traveling time on the number of the additional UAVs. Table II lists the
parameters used in the numerical analysis [19].

TABLE II: Parameters in numerical analysis UAV energy capacity | 0.88kW.h
---|---
Power consumption by the electronics | 0.15kW
Grid size | 4x4
Area of the graph | 1kmx1km
Traveling time through edge | 2.5 min
Charging station location (x,y) | (0,0)
Charging time | 5 min
UAV weight with battery | 8.5 k.g
Maximum payload weight | 2 k.g
MAX forward speed | 12 m/s
Figure 2: The cycles that cover the subareas using the CLE algorithm Figure 3:
Energy capacity vs. The number of addtional UAVs Figure 4: Grid size vs. The
number of additional UAVs

In Figure 2, we uniformly divide the geographical area into 16 subareas and
apply the CLE algorithm to find the cycles with minimum number of additional
UAVs. From the figure, we notice that 5 cycles are needed to cover all
subareas with 5 additional UAVs. Also, we note that the paths of the cycles
are intersected in many locations. To avoid the collisions between the UAVs,
we operate the paths (cycles) at different altitudes with small altitude
differences.

In Figure 3, we study the effect of the UAV energy capacity on the number of
the additional UAVs needed to cover the subareas. Actually, when increasing
the energy capacity of the UAV and apply SM, the number of additional UAVs
needed will not change because each subarea is covered by one cycle and two
UAVs. When increasing the energy capacity of the UAVs, only the coverage time
of each UAV increases. on the other hand, increasing the energy capacity of
each UAV results in minimizing the number of additional UAVs that needed using
CLE. This is because increasing the energy capacity of each UAV gives the UAV
a chance to visit and to cover more subareas, which minimizes the number of
the cycles that are needed to cover the subareas.

Figure 5: Charging time vs. The number of additional UAVs Figure 6: Traveling
time vs. The number of additional UAVs

In Figure 4, the slope of the line produced by SM is greater than the curve of
CLE. When applying SM, the number of additional UAVs increases linearly with
the grid size. This is because the number of additional UAVs equals the grid
size. Also, when applying the CLE, the number of additional UAVs increases
with the grid size. This is because more cycles are needed to cover more
subareas and each cycle will need one additional UAV.

In Figure 5, we study the effect of the charging time on the number of
additional UAVs needed. Changing the charging time will not affect the number
of additional UAVs needed when applying SM. This is because the coverage time
of each UAV will cover the time that the UAV needs to return to the charging
station to recharge the battery and to visit the subarea again. On the other
hand, when applying CLE, it will be a critical issue (see Theorem 2).
Actually, charging the battery of the UAV takes long time. For this reason,
each UAV has a replacement battery [19]. In this paper, we assume the time
needed to replace the battery for each UAV is 5 minutes.

In Figure 6, we study the effect of the traveling time on the number of
additional UAVs. Changing the traveling time will not affect the number of
additional UAVs when applying SM. On the other hand, it will be a critical
issue to choose the appropriate traveling time when applying CLE. When
increasing the traveling time, the wasted energy during traveling will
increase and the coverage time will decrease. Hence, the chance to visit other
subareas will decrease.

##  V Conclusion

In this paper, we study the continuous coverage problem with minimum number of
replacement UAVs and prove that it is NP complete. We design an efficient
algorithm to solve it, the cycles with limited energy algorithm. The proposed
algorithm covers the $N$ subareas by cycles, in which each cycle needs one
additional UAV to ensure continuous coverage. We showed that the energy
capacity of the UAVs, the number of subareas in the affected area, and the UAV
charging and traveling times will all impact the required number of UAVs. Our
simulation results showed that applying the cycles with limited energy
algorithm, can efficiently reduce the number of additional UAVs needed
relative to the straightforward method. As future work, we will study the
continuous coverage problem using UAVs with non-uniform energy capacities and
the use of green energy.

## Acknowledgment

The work of Jacob Chakareski has been partially supported by the NSF under
award CCF-1528030.

## References

  * [1] P. Bupe, R. Haddad, and F. Rios-Gutierrez, “Relief and emergency communication network based on an autonomous decentralized uav clustering network,” in _SoutheastCon 2015_. IEEE, 2015, pp. 1–8.
  * [2] R. Miller, “Hurricane katrina: Communications & infrastructure impacts,” DTIC Document, Tech. Rep., 2006.
  * [3] F. F. Townsend _et al._ , “The federal response to hurricane katrina: Lessons learned,” _Washington, DC: The White House_ , 2006.
  * [4] E. M. Noam, “What the world trade center attack has shown us about our communications networks,” _Global Economy and Digital Society. Amsterdam: Elsevier Science_ , pp. 375–378, 2004.
  * [5] I. Dalmasso, I. Galletti, R. Giuliano, and F. Mazzenga, “Wimax networks for emergency management based on uavs,” in _IEEE First AESS European Conference on Satellite Telecommunications (ESTEL)_. IEEE, 2012, pp. 1–6.
  * [6] T. Ha, “The uav continuous coverage problem,” DTIC Document, Tech. Rep., 2010\.
  * [7] J. Wu, “Collaborative mobile charging and coverage,” _Journal of Computer Science and Technology_ , vol. 29, no. 4, pp. 550–561, 2014.
  * [8] S. Zhang, J. Wu, and S. Lu, “Collaborative mobile charging,” _IEEE Transactions on Computers_ , vol. 64, no. 3, pp. 654–667, 2015.
  * [9] L. D. P. Pugliese, F. Guerriero, D. Zorbas, and T. Razafindralambo, “Modelling the mobile target covering problem using flying drones,” _Optimization Letters_ , pp. 1–32, 2015.
  * [10] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Drone small cells in the clouds: Design, deployment and performance analysis,” in _IEEE Global Communications Conference (GLOBECOM)_ , 2015, pp. 1–6.
  * [11] ——, “Optimal transport theory for power-efficient deployment of unmanned aerial vehicles,” _IEEE International Conference on Communications (ICC), Kuala Lumpur, Malaysia,_ , 2016.
  * [12] D. Zorbas, T. Razafindralambo, F. Guerriero _et al._ , “Energy efficient mobile target tracking using flying drones,” _Procedia Computer Science_ , vol. 19, pp. 80–87, 2013.
  * [13] E. T. Ceran, T. Erkilic, E. Uysal-Biyikoglu, T. Girici, and K. Leblebicioglu, “Optimal energy allocation policies for a high altitude flying wireless access point,” _Transactions on Emerging Telecommunications Technologies_ , 2016.
  * [14] L. Gupta, R. Jain, and G. Vaszkun, “Survey of important issues in uav communication networks,” _IEEE Communications Surveys and Tutorials_ , 2015\.
  * [15] R. E. Korf, “A new algorithm for optimal bin packing,” in _AAAI/IAAI_ , 2002, pp. 731–736.
  * [16] D. S. Johnson and L. A. McGeoch, “The traveling salesman problem: A case study in local optimization,” _Local search in combinatorial optimization_ , vol. 1, pp. 215–310, 1997.
  * [17] M. Wannberg, “The quadrotor platform: from a military point of view,” 2012.
  * [18] R. D’Andrea, “Guest editorial can drones deliver?” _IEEE Transactions on Automation Science and Engineering_ , vol. 11, no. 3, pp. 647–648, 2014.
  * [19] “Sentinel+ drone, http://www.airbornedrones.co/,” 2016.


# Magnetic fluctuations in Pb9Cu(PO4)6O

Makoto Shimizu  Department of Physics, Graduate School of Science, Kyoto
University, Kyoto 606-8502, Japan  Junya Otsuki  Research Institute for
Interdisciplinary Science, Okayama University, Okayama 700-8530, Japan  Harald
O. Jeschke  Research Institute for Interdisciplinary Science, Okayama
University, Okayama 700-8530, Japan

###### Abstract

The hope that copper doped lead apatite Pb9Cu(PO4)6O is a room-temperature
superconductor has largely been dashed by global research efforts.
Nevertheless, according to the current state of knowledge, the material has
interesting magnetic properties, and research groups around the world have
prepared high quality samples. We use a fluctuation exchange approximation
(FLEX) approach to study the magnetic tendencies in Pb9Cu(PO4)6O. We find
ferromagnetic fluctuations very close to the filling of the stoichiometric
compound which can be understood from Fermi surface nesting at the $M$ point.
This is similar to the one-band triangular lattice Hamiltonian at three-
quarter filling. Interestingly, the special $k_{z}$ dependence of the
Pb9Cu(PO4)6O band structure makes it very sensitive to doping. Only slight
charge doping switches between antiferromagnetic and ferromagnetic
fluctuations. If the material could become superconducting, it might be easily
switchable between singlet and triplet superconductivity.

##  I Introduction

Copper doped lead apatite Pb9Cu(PO4)6O has recently been put forward as a
candidate for room temperature superconductivity at ambient pressure [1, 2],
causing substantial excitement in the scientific community. Efforts at
independent synthesis of the material and at reproduction of the experimental
observations have been undertaken by many groups [3, 4, 5]. Every failure to
reproduce superconductivity [6] has added to the doubts concerning the
interpretation of the original experiments. However, the suggestion that a
well known phase transition in Cu2S [7] which is a known impurity [1] is
responsible for the sudden change in resistivity around 104 Kelvin in the
mixed-phase Pb9Cu(PO4)6O samples is rather convincing [8]. Nevertheless, it is
still an interesting question what the properties of the impurity free
Pb9Cu(PO4)6O material are. In particular, the half levitation [2] that has
been reproduced by some groups [5, 9] points to interesting magnetic
properties. The fact that excellent single crystals of Pb9Cu(PO4)6O have been
grown [10] means that this question is now accessible to precise experiments.

The experimental claims for Pb9Cu(PO4)6O have been met with a concerted
response of the theoretical materials science community. Very fast density
functional theory studies [11, 12, 13, 14, 15] have discussed the importance
of flat bands and the shape of the Fermi surface. Several different tight
binding models have been proposed [13, 16, 17] and analyzed with respect to
topological properties [18]. Three DFT+DMFT studies show that at the integer
filling of the stoichiometric Pb9Cu(PO4)6O compound, correlations open a gap
[19, 20, 21] in agreement with the transparent single crystals [10].
Superconductivity has been theoretically studied using the two-dimensional
$t$-$J$ model [22] and spin fluctuation theory [23]. First statements about
the magnetism of Pb9Cu(PO4)6O have been obtained [24].

In this study, we will study the magnetic fluctuations in Cu doped Pb10(PO4)6O
using the fluctuation exchange approximation (FLEX) [25, 26]. The
stoichiometric Pb9Cu(PO4)6O compound is described by a three-quarter filled
two-orbital model. Its Fermi surface has been described as rugby ball shaped.
We will take the three-dimensional nature of the compound seriously and work
out similarity and differences between the electronic structure of
Pb9Cu(PO4)6O and a one-band model at different fillings on the triangular
lattice.

Figure 1: (a) Predicted simplified crystal structure of Pb9Cu(PO4)6O. (b)
Electronic band structure near the Fermi level at the GGA level with Cu $3d$
orbital character highlighted. Figure 2: GGA Fermi surfaces of Pb9Cu(PO4)6O
for several doping levels $\Delta n=n-7$. (c) is the Fermi surface of the
stoichiometric compound, (a) and (b) correspond to hole doping of 0.2 and 0.1
electrons per formula unit, and (d) shows electron doping of 0.1 electrons.
The color red (blue) indicates high (low) Fermi velocity. Figure 3:  Cuts
trough the static spin susceptibility $\chi^{\mathrm{s}}(\bm{q})$ for (a)
$q_{z}=0$, (b) $q_{z}=\pi/2c$ and (c) $q_{z}=\pi/c$ at $T=290\,\mathrm{\>K}$
calculated within FLEX. Interaction parameters are $U=200\mathrm{\>meV}$,
$V=U/2$ and $J=J^{\prime}=U/4$. High-symmetry points are
$\bm{Q}_{\Gamma}\equiv(0,0,0)$, $\bm{Q}_{\mathrm{K}}\equiv(1/\sqrt{3},1/3,0)$,
$\bm{Q}_{\mathrm{M}}\equiv(1/\sqrt{3},0,0)$,
$\bm{Q}_{\mathrm{A}}\equiv(0,0,1/2)$,
$\bm{Q}_{\mathrm{H}}\equiv(1/\sqrt{3},1/3,1/2)$,
$\bm{Q}_{\mathrm{L}}\equiv(1/\sqrt{3},0,1/2)$, in units of $2\pi/a$ for the
$x$ and $y$ components and of $2\pi/c$ for the $z$ component. The high
symmetry points are shown in (a1) and (a3).  Figure 4:  The static spin
susceptibility $\chi^{\mathrm{s}}(\bm{q})$ of the triangular lattice at
$T=|t|/4$ calculated within FLEX. The interaction is $U=9|t|=W$ where $t$ is
the nearest neighbor hopping parameter and $W$ is the band width. We use a
$32\times 32$ $\bm{k}$-mesh.  Figure 5:  Momentum-resolved density of states
$\rho_{0}(\bm{k})$ at the Fermi level in a non-interacting triangular lattice
model at $n=0.8$ and $n=1.5$ An artificial broadening of width $0.1t$ was
applied.  Figure 6:  Band structure of Pb9Cu(PO4)6O along a triangular $k$
path for different values of $k_{z}$. The $x$ and $y$ components are in units
of $2\pi/a$, and the $z$ component is in units of $2\pi/c$. The path at
$k_{z}=0$ is equivalent to $\Gamma$-M-K-$\Gamma$.

##  II Electronic structure

We perform all electron density functional theory (DFT) calculations for
Pb9Cu(PO4)6O using the full potential local orbital (FPLO) basis [27] and the
generalized gradient approximation exchange correlation functional [28]. We
derive a precise four-band tight binding model using the symmetry conserving
projected Wannier functions implemented within FPLO [29].

We base our model of the structure of Pb9Cu(PO4)6O on the crystal structure of
lead apatite Pb10(PO4)6O as determined by Krivovichev and Burns [30]. Its
$P6_{3}/m$ (no. 176) space group has an O(4) position that is only occupied by
one quarter. We first simplify by choosing one of four symmetry equivalent
O(4) positions which reduces the symmetry to $P3$ (no. 143). According to Ref.
13, the Pb(2) position with four symmetry equivalent sites is energetically
more favorable for Cu substitution. We thus replace one out of four Pb by Cu
and perform an internal relaxation of this structure. As DFT structure
prediction has well known issues for correlated oxides, we use experimental
[1] rather than DFT predicted lattice parameters. The resulting crystal
structure is shown in Fig. 1 (a). It is similar to the structures considered
in other works [13, 31].

We show the DFT band structure and density of states of Pb9Cu(PO4)6O in Fig. 1
(b) and (c). Two almost degenerate bands of Cu $3d_{xz}$ and $3d_{yz}$
character are crossing the Fermi level. They are slightly hybridizing with two
occupied oxygen bands arising from the extra O in the formula, the O(4) in the
Pb10(PO4)6O structure; this has also been observed in Ref. [31]. Therefore, in
this study we work with a four-band tight binding model. The upper of the two
bands at the Fermi level leads to the three-dimensional Fermi surface shown in
Fig. 2 (c). It has mount Fuji shape because the finite $k_{z}$ dispersion
leads to different filling levels of the underlying triangular lattice at
every value of $k_{z}$. Slight hole doping (Fig. 2 (b) and (a)) lead to the
rugby ball shape that has been found previously [11] and tiny Fermi surfaces
from the second Cu $3d$ band. The shape remains qualitatively the same between
undoped and moderate electron doping (Fig. 2 (d)). This $k_{z}$ dispersion of
the electronic structure will constitute the main difference to a fully two-
dimensional triangular lattice Hamiltonian.

We determine [29] a four-band tight-binding model including $3d_{xz}$ and
$3d_{yz}$ orbitals of Cu and $2p_{x}$ and $2p_{y}$ orbitals of O. The
Hamiltonian $\mathcal{H}_{0}$ is given by

$\mathcal{H}_{0}=\sum_{ij\mu\sigma\sigma}t_{ij}^{\mu\nu}c^{{\dagger}}_{i\mu\sigma}c_{j\nu\sigma},$
(1)

where $i$ and $j$ are site indices, $\mu$ and $\nu$ are orbital indices, and
$\sigma$ is the spin index. The occupation number $n$ is $n=7$ for the
stoichiometric Pb9Cu(PO4)6O. Three bands are fully occupied, and the top band
is half-filled without doping and interactions. This means that at
stoichiometry, one band forms the Fermi surface. However, the two upper bands
deriving almost entirely from Cu $3d$ are nearly degenerate, and the orbital
character of the bands and Fermi surface changes between $3d_{xz}$ and
$3d_{yz}$ depending on $k$ direction.

##  III Spin fluctuations

In order to discuss spin fluctuations, we consider interactions between $3d$
electrons on Cu. We adopt the usual local interactions, which are represented
by the Hamiltonian

$\displaystyle\mathcal{H}_{\mathrm{int}}$
$\displaystyle=\sum_{i}\Big{[}\frac{U}{2}\sum_{l\sigma}n_{il\sigma}n_{il\bar{\sigma}}+\frac{V}{2}\sum_{l\neq
m}n_{il}n_{im}$ $\displaystyle-\frac{J}{2}\sum_{l\neq
m}\bm{S}_{il}\cdot\bm{S}_{im}+\frac{J^{\prime}}{2}\sum_{l\neq
m}\sum_{\sigma}c^{\dagger}_{il\sigma}c^{\dagger}_{il\sigma}c_{im\bar{\sigma}}c_{im\sigma}\Big{]},$
(2)

where the orbital indices $l$ and $m$ run over two orbitals ($3d_{xz}$ and
$3d_{yz}$) out of four. $U$, $V$, $J$, and $J^{\prime}$ are intra-orbital
Coulomb interaction, inter-orbital Coulomb interaction, Hund’s coupling, and
pair hopping, respectively.

We apply the FLEX approximation to the four-band Hubbard model. We change the
electron number in the range of $6.5\leq n\leq 7.5$. We choose spin rotation-
invariant interaction parameters [32] $V=U/2$ and $J=J^{\prime}=U/4$ for the
Cu $3d$ orbitals, and we assume $U=200\mathrm{\>meV}$. We use a $16\times
16\times 4$ $\bm{k}$-mesh.

We discuss the momentum-dependent static spin susceptibility,
$\chi^{\mathrm{s}}(\bm{q})=\sum_{lm}\chi^{\mathrm{s}}_{lm}(\bm{q})$, in which
contributions from $3d_{xz}$ and $3d_{yz}$ orbitals are summed up. Figure 3
shows cuts of $\chi^{\mathrm{s}}(\bm{q})$ at $q_{z}=0$, $\pi/2c$, and $\pi/c$
(zone boundary) for $6.7\leq n\leq 7.4$. At the formal electron filling of
Pb9Cu(PO4)6O ($n=7.0$), $\chi^{\mathrm{s}}(\bm{q})$ exhibits peaks at
$\bm{q}=\bm{Q}_{\mathrm{H}}$. This fluctuation corresponds to an
antiferromagnetic (AFM) stacking of the 120∘ state on the triangular lattice.
The AFM fluctuation persists for hole doping. On the electron doping side, on
the other hand, the magnetic properties change drastically. At $n=7.1$, a peak
is present at $\bm{q}=\bm{Q}_{\mathrm{\Gamma}}$, which corresponds to the
ferromagnetic (FM) fluctuation. The FM fluctuation becomes dominant for
$n\gtrsim 7.2$.

Let us consider the origin of the AFM and FM fluctuations in Pb9Cu(PO4)6O in
terms of a single-band Hubbard model on a triangular lattice. This model
approximately describes the energy dispersion of the top band through
$\Gamma$-M-K-$\Gamma$ in Fig. 1. Figure 4 shows the spin susceptibility
$\chi^{\mathrm{s}}(\bm{q})$ for several values of the occupation number $n$.
Here, $n$ in this model corresponds to $n-6$ in the four-band model of
Pb9Cu(PO4)6O. $\chi^{\mathrm{s}}(\bm{q})$ exhibits peaks at
$\bm{Q}_{\mathrm{K}}$ around half-filling ($0.6\leq n\leq 1.0$) while it has a
peak at $\bm{Q}_{\mathrm{\Gamma}}$ around 3/4-filling ($1.4\leq n\leq 1.6$).
These magnetic fluctuations can be understood well with the electronic
structure. Around half-filling, the Fermi surface in the corresponding non-
interacting system is a circle around the $\Gamma$ point as shown in Fig. 5
(a). The nesting vector agrees with $\bm{Q}_{\mathrm{K}}$ as indicated by the
arrow in Fig. 5 (a). On the other hand, the Fermi surface around three-quarter
filling is a hexagon inscribed in the first Brillouin zone as shown in Fig. 5
(b). The density of states is large at the M point because of the van Hove
singularity. Therefore, dominant particle-hole excitations occur in the
$\bm{q}$-vector that connects two M points, that is, $(2\pi/a,0,0)$, which
corresponds to $\bm{Q}_{\mathrm{\Gamma}}$. The Fermi surface thus explains AFM
fluctuations corresponding to the 120∘ state for $n\lesssim 1$ and the strong
FM fluctuations for $n\simeq 1.5$.

The variations of the magnetic properties between AFM and FM can be understood
with the triangular lattice. However, Pb9Cu(PO4)6O shows two unique properties
that the triangular lattice does not show. First, the change from
antiferromagnetic to ferromagnetic is very sensitive to electron doping. It is
antiferromagnetic at $n=7.0$, and even with very small electron doping, it
starts to be ferromagnetic at $n=7.1$. 0.1 electrons are needed to cause
ferromagnetic fluctuations when we start from stoichiometric integer filling.
In the triangular lattice, on the other hand, 0.5 electrons are needed. Doping
0.5 electrons into a material is usually very difficult in reality. The second
unique property of Pb9Cu(PO4)6O is the fact that ferromagnetism persists in
the whole region of electron doping while it disappears in the triangular
lattice.

These two properties of Pb9Cu(PO4)6O result from the three-dimensional Fermi
surface. Even though the highest energy band is almost half-filled in total,
various situations of filling are realized if we look into the two-dimensional
$k$-path for different $k_{z}$: almost empty at $k_{z}=\pi/c$ and almost
fully-filled at $k_{z}=0$ as shown in Fig. 6. Therefore, the van Hove
singularity at ($1/\sqrt{3}$, 0, $k_{z}$) appears on the Fermi level even at
the stoichiometric filling, and the ferromagnetic fluctuation is enhanced by
the mechanism discussed in Fig. 5 (b). Moreover, when electrons are doped, the
Fermi level passes the singularity point at a different $k_{z}$. This explains
the result that the ferromagnetic fluctuation is dominant in a wide range of
electron doping levels in Pb9Cu(PO4)6O.

##  IV Conclusions

We have studied the consequences of the electronic structure of Pb9Cu(PO4)6O
for magnetic fluctuations. We find that electronically, Pb9Cu(PO4)6O provides
an interesting three-quarter filled, two orbital triangular lattice
Hamiltonian. As such, it can be partly understood from the one-band triangular
lattice Hamiltonian that has previously been studied in the context of organic
charge transfer salts and of sodium cobaltate. However, we have worked out two
important differences that arise from the three-dimensionality of
Pb9Cu(PO4)6O, i.e. from the $k_{z}$ dispersion of the electronic structure.
One consequence is that the change of doping that triggers the transition
between antiferromagnetic and ferromagnetic spin fluctuations is much smaller
in Pb9Cu(PO4)6O compared to the simple triangular lattice. The other
consequence is that the ferromagnetic fluctuations in the electron doped
region is significantly stabilized in the three-dimensional two-orbital model.

While our FLEX result for the undoped material indicates antiferromagnetic
fluctuations, we also find that ferromagnetic fluctuations are reached by very
little electron doping. Thus, slightly inhomogeneous samples Pb10-xCux(PO4)6O
with $x\approx 1$ would have ferromagnetic as well as antiferromagnetic
regions, and this could explain the levitation observed in some of the samples
[2, 5, 9] of Pb9Cu(PO4)6O.

By focusing on a single site for the position of the doped Cu2+ ion, we have
introduced a significant simplification, just like in many other studies [20,
23, 21, 23]. It is possible that disorder effects will play an important role
in Pb9Cu(PO4)6O [24]. Nevertheless, our study shows that the electronic
structure that is present in slightly simplified Pb9Cu(PO4)6O has very
interesting properties. It is characterized by two degenerate orbitals on a
triangular lattice, three-quarter filling and a finite $k_{z}$ dispersion. As
this is a rather generic situation, it could already be realized in other
existing materials, or it would not be too hard to design materials with this
property. As the hope for superconductivity in Pb9Cu(PO4)6O is probably in
vain, it makes sense to understand why and possibly find or design a similar
material that does support a superconducting ground state. A likely
disadvantage of Pb9Cu(PO4)6O is the large ratio between interaction strength
on the $d$ orbitals and the band width [23]; in the search for a better
Pb9Cu(PO4)6O, wider bands or weaker local interactions would be desirable. If
superconductivity could be realized in a dopable material with the electronic
structure of Pb9Cu(PO4)6O, then we would expect that the order parameter could
be tuned between singlet and triplet.

###### Acknowledgements.

We acknowledge useful discussions with S. Kasahara and K. Kobayashi. The
computation in this work has been done using the facilities of the
Supercomputer Center, the Institute for Solid State Physics, the University of
Tokyo. M.S. was supported by JSPS KAKENHI grant No. 22H01181. J.O. was
supported by JSPS KAKENHI grant No. 21H01041, No. 21H01003, No. 23H04869.

## References

  * Lee  _et al._ [2023a] S. Lee, J.-H. Kim, and Y.-W. Kwon, The first room-temperature ambient-pressure superconductor (2023a), arXiv:2307.12008 [cond-mat.supr-con] .
  * Lee  _et al._ [2023b] S. Lee, J. Kim, H.-T. Kim, S. Im, S. An, and K. H. Auh, Superconductor Pb10-xCux(PO4)6O showing levitation at room temperature and atmospheric pressure and mechanism (2023b), arXiv:2307.12037 [cond-mat.supr-con] .
  * Kumar  _et al._ [2023a] K. Kumar, N. K. Karn, and V. P. S. Awana, Synthesis of possible room temperature superconductor LK-99: Pb9Cu(PO4)6O (2023a), arXiv:2307.16402 [cond-mat.supr-con] .
  * Hou  _et al._ [2023] Q. Hou, W. Wei, X. Zhou, Y. Sun, and Z. Shi, Observation of zero resistance above 100∘ K in Pb10-xCux(PO4)6O (2023), arXiv:2308.01192 [cond-mat.supr-con] .
  * Guo  _et al._ [2023] K. Guo, Y. Li, and S. Jia, Ferromagnetic half levitation of LK-99-like synthetic samples (2023), arXiv:2308.03110 [cond-mat.supr-con] .
  * Kumar  _et al._ [2023b] K. Kumar, N. K. Karn, Y. Kumar, and V. P. S. Awana, Absence of superconductivity in LK-99 at ambient conditions (2023b), arXiv:2308.03544 [cond-mat.supr-con] .
  * Jain [2023] P. K. Jain, Phase transition of copper (i) sulfide and its implication for purported superconductivity of LK-99 (2023), arXiv:2308.05222 [cond-mat.supr-con] .
  * Zhu  _et al._ [2023] S. Zhu, W. Wu, Z. Li, and J. Luo, First order transition in Pb10-xCux(PO4)6O ($0.9<x<1.1$) containing Cu2S (2023), arXiv:2308.04353 [cond-mat.supr-con] .
  * Wu  _et al._ [2023] H. Wu, L. Yang, B. Xiao, and H. Chang, Successful growth and room temperature ambient-pressure magnetic levitation of LK-99 (2023), arXiv:2308.01516 [cond-mat.supr-con] .
  * Puphal  _et al._ [2023] P. Puphal, M. Y. P. Akbar, M. Hepting, E. Goering, M. Isobe, A. A. Nugroho, and B. Keimer, Single crystal synthesis, structure, and magnetism of Pb10-xCux(PO4)6O (2023), arXiv:2308.06256 [cond-mat.supr-con] .
  * Lai  _et al._ [2024] J. Lai, J. Li, P. Liu, Y. Sun, and X.-Q. Chen, First-principles study on the electronic structure of Pb10-xCux(PO4)6O ($x=0$, 1), J. Mater. Sci. Technol. 171, 66 (2024).
  * Griffin [2023] S. M. Griffin, Origin of correlated isolated flat bands in copper-substituted lead phosphate apatite (2023), arXiv:2307.16892 [cond-mat.supr-con] .
  * Si and Held [2023] L. Si and K. Held, Electronic structure of the putative room-temperature superconductor Pb9Cu(PO4)6O (2023), arXiv:2308.00676 [cond-mat.supr-con] .
  * Kurleto  _et al._ [2023] R. Kurleto, S. Lany, D. Pashov, S. Acharya, M. van Schilfgaarde, and D. S. Dessau, Pb-apatite framework as a generator of novel flat-band CuO based physics, including possible room temperature superconductivity (2023), arXiv:2308.00698 [cond-mat.supr-con] .
  * Cabezas-Escares  _et al._ [2023] J. Cabezas-Escares, N. F. Barrera, C. Cardenas, and F. Munoz, Theoretical insight on the LK-99 material (2023), arXiv:2308.01135 [cond-mat.supr-con] .
  * Tavakol and Scaffidi [2023] O. Tavakol and T. Scaffidi, Minimal model for the flat bands in copper-substituted lead phosphate apatite (2023), arXiv:2308.01315 [cond-mat.supr-con] .
  * Lee and Dai [2023] P. A. Lee and Z. Dai, Effective model for pb9cu(po4)6o (2023), arXiv:2308.04480 [cond-mat.supr-con] .
  * Hirschmann and Mitscherling [2023] M. M. Hirschmann and J. Mitscherling, Tight-binding models for SG 143 (P3) and application to recent DFT results on copper-doped lead apatite (2023), arXiv:2308.03751 [cond-mat.mes-hall] .
  * Korotin  _et al._ [2023] D. M. Korotin, D. Y. Novoselov, A. O. Shorikov, V. I. Anisimov, and A. R. Oganov, Electronic correlations in promising room-temperature superconductor Pb9Cu(PO4)6O: a DFT+DMFT study (2023), arXiv:2308.04301 [cond-mat.supr-con] .
  * Si  _et al._ [2023] L. Si, M. Wallerberger, A. Smolyanyuk, S. di Cataldo, J. M. Tomczak, and K. Held, Pb10-xCux(PO4)6O: a mott or charge transfer insulator in need of further doping for (super)conductivity (2023), arXiv:2308.04427 [cond-mat.supr-con] .
  * Yue  _et al._ [2023] C. Yue, V. Christiansson, and P. Werner, Correlated electronic structure of Pb10-xCux(PO4)6O (2023), arXiv:2308.04976 [cond-mat.str-el] .
  * Oh and Zhang [2023] H. Oh and Y.-H. Zhang, S-wave pairing in a two-orbital t-J model on triangular lattice: possible application to Pb10-xCux(PO4)6O (2023), arXiv:2308.02469 [cond-mat.supr-con] .
  * Witt  _et al._ [2023] N. Witt, L. Si, J. M. Tomczak, K. Held, and T. Wehling, No superconductivity in Pb9Cu1(PO4)6O found in orbital and spin fluctuation exchange calculations (2023), arXiv:2308.07261 [cond-mat.supr-con] .
  * Sun  _et al._ [2023] Y. Sun, K.-M. Ho, and V. Antropov, Metallization and spin fluctuations in Cu-doped lead apatite (2023), arXiv:2308.03454 [cond-mat.supr-con] .
  * Bickers  _et al._ [1989] N. E. Bickers, D. J. Scalapino, and S. R. White, Conserving Approximations for Strongly Correlated Electron Systems: Bethe-Salpeter Equation and Dynamics for the Two-Dimensional Hubbard Model, Phys. Rev. Lett. 62, 961 (1989).
  * Ikeda  _et al._ [2010] H. Ikeda, R. Arita, and J. Kuneš, Phase diagram and gap anisotropy in iron-pnictide superconductors, Phys. Rev. B Condens. Matter 81, 054502 (2010).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59, 1743 (1999).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Koepernik  _et al._ [2023] K. Koepernik, O. Janson, Y. Sun, and J. van den Brink, Symmetry-conserving maximally projected Wannier functions, Phys. Rev. B 107, 235135 (2023).
  * Krivovichev and Burns [2003] S. V. Krivovichev and P. C. Burns, Crystal chemistry of lead oxide phosphates: crystal structures of Pb4O(PO4)2, Pb8O5(PO4)2 and Pb10(PO4)6O, Z. Kristallogr. 218, 357 (2003).
  * Jiang  _et al._ [2023] Y. Jiang, S. B. Lee, J. Herzog-Arbeitman, J. Yu, X. Feng, H. Hu, D. Călugăru, P. S. Brodale, E. L. Gormley, M. G. Vergniory, C. Felser, S. Blanco-Canosa, C. H. Hendon, L. M. Schoop, and B. A. Bernevig, Pb9Cu(PO4)6(OH)2: Phonon bands, localized flat band magnetism, models, and chemical analysis (2023), arXiv:2308.05143 [cond-mat.supr-con] .
  * Guterding  _et al._ [2015] D. Guterding, H. O. Jeschke, P. J. Hirschfeld, and R. Valentí, Unified picture of the doping dependence of superconducting transition temperatures in alkali metal/ammonia intercalated FeSe, Phys. Rev. B 91, 041112 (2015).