# Activity gradients in two- and three-dimensional active nematics

Liam J. Ruske<EMAIL_ADDRESS>Julia M. Yeomans  Rudolf Peierls
Centre For Theoretical Physics, University of Oxford, UK

###### Abstract

We numerically investigate how spatial variations of extensile or contractile
active stress affect bulk active nematic systems in two and three dimensions.
In the absence of defects, activity gradients drive flows which re-orient the
nematic director field and thus act as an effective anchoring force. At high
activity, defects are created and the system transitions into active
turbulence, a chaotic flow state characterized by strong vorticity. We find
that in two-dimensional (2D) systems active torques robustly align $+1/2$
defects parallel to activity gradients, with defect heads pointing towards
contractile regions. In three-dimensional (3D) active nematics disclination
lines preferentially lie in the plane perpendicular to activity gradients due
to active torques acting on line segments. The average orientation of the
defect structures in the plane perpendicular to the line tangent depends on
the defect type, where wedge-like $+1/2$ defects align parallel to activity
gradients, while twist defects are aligned anti-parallel. Understanding the
response of active nematic fluids to activity gradients is an important step
towards applying physical theories to biology, where spatial variations of
active stress impact morphogenetic processes in developing embryos and affect
flows and deformations in growing cell aggregates, such as tumours.

##  I Introduction

Active systems continuously take energy from their environment and convert it
into mechanical stress which drives the system out of equilibrium. Active
matter has recently received considerable attention in the community because
of its potential as a way of interpreting biological mechanics and as examples
of non-equilibrium statistical physics [1, 2, 3].

Dense active systems, and in particular the continuum theory of active
nematics, have been used successfully to model the dynamics of a variety of
biological systems, such as myosin-driven actin-microtubule networks [4],
bacterial biofilms [5] and confluent cell monolayers [6, 7, 8, 9]. A key
property of active nematics, which distinguishes them from passive liquid
crystals, is active turbulence. This is a chaotic flow state characterised by
strong vorticity and topological defects which are continually created and
destroyed. Unlike in passive systems, where topological defects of opposite
charge tend to annihilate over time due to the relaxation of the large elastic
energy associated with them, in active systems the number of topological
defects is finite at steady-state, and motile $+1/2$ defects move through the
fluid like self-propelled, oriented particles [10]. Considerable experimental
and theoretical work has been devoted to understanding the statistical
properties of active turbulence and the motion of the associated topological
defects, in both two and three dimensions [11, 12, 13, 14].

A major research theme has been to investigate how to control the chaotic
flows and defect motion in active materials. Guiding active flow and defects
along selected paths is likely to be a prerequisite to using active matter for
powering micro-machines or targeting the delivery of drugs. It has been
proposed that this control can be achieved in two-dimensional (2D) systems by
patterns of activity on a substrate. This causes the orientation of $+1/2$
defects to align along activity gradients and hence the defects to move along
paths defined by activity patterns [15, 16]. Further work has investigated the
alignment of individual topological defects by activity gradients in 2D using
theoretical arguments and numerical simulations [17, 18]. More recently it has
also proved possible to design active materials that allow imaging of a three-
dimensional (3D) active nematic and, in particular, the associated motile
disclination loops and lines [11]. These articles have focused primarily on
how activity patterns affect flows and the alignment of individual topological
defects, and it is still poorly understood how activity gradients affect
systems in the absence of topological defects or how well the concept of 2D
topological defects as effective oriented particles translates to disclination
line segments in 3D materials.

In this paper we numerically investigate how active nematics respond to
spatial variations in activity in both two and three dimensions. We start by
outlining the well-established hydrodynamic theory of active nematics and the
numerical methods we use to solve the equations of motion. In section III we
show how non-uniform activity creates active forces which, in the absence of
defects, set up flows which rotate the nematic director field. This
effectively creates an active force which favours either normal or tangential
alignment of the nematic director with respect to activity gradients for
contractile or extensile active stress, respectively. In the following section
IV, we show that activity gradients in 2D active turbulence induce polar order
of $+1/2$ defects, which is consistent with previous theoretical predictions
considering active torques on isolated defects [17]. In section V we
investigate how activity gradients in 3D systems affect nematic alignment and
the polarisation of disclination lines. We find that disclination lines, which
locally resemble $+1/2$ or twist defects, are polarized in opposite
directions, aligning parallel or anti-parallel to activity gradients,
respectively. The last section of the paper summarizes the key results,
highlights similarities between activity gradients and active-passive
interfaces and points out possible connections to experiment.

##  II  Equations of motion

We numerically solve the continuum equations of motion of an active nematic
fluid [19]. The ordering of the nematic constituents in 3D is quantified by
the tensor order parameter $Q_{ij}=\frac{3}{2}S(n_{i}n_{j}-\delta_{ij}/3)$,
where the scalar $S$ quantifies the magnitude of the order and the director
field $\mathbf{n}$ the direction. We work in the uniaxial limit, ignoring
higher-order effects such as biaxial order close to disclination cores [20,
21]. The evolution of the nematic tensor $\mathbf{Q}$ is coupled to the
velocity field $\mathbf{u}$ of the fluid and follows [22]

$D_{t}\mathbf{Q}-\mathbf{\mathcal{W}}=\Gamma\mathbf{H}\>,$ (1)

where $D_{t}$ denotes the material derivative and $\Gamma$ the rotational
diffusivity which controls the relaxation of the order parameter towards
equilibrium through the molecular field
$\mathbf{H}=-\delta\mathcal{F}/\delta\mathbf{Q}+(\mathbf{I}/3)\text{Tr}(\delta\mathcal{F}/\delta\mathbf{Q})$.
The free energy $\mathcal{F}=\int fdV$ of the system includes a Landau-de
Gennes bulk term
$f_{bulk}=\frac{A}{2}\text{Tr}~{}\mathbf{Q}^{2}+\frac{B}{3}\text{Tr}~{}\mathbf{Q}^{3}+\frac{C}{4}(\text{Tr}~{}\mathbf{Q}^{2})^{2}$,
where we chose the coefficients $A,B,C$ such that nematic ordering $S_{eq}>0$
is favoured in equilibrium, and an elastic term
$f_{el}=K_{el}\left(\mathbf{\nabla}\mathbf{Q}\right)^{2}$, which penalizes
distortions in the director field $\mathbf{n}$ using the one-elastic-constant
approximation. This commonly used approximation has been shown to accurately
describe many experimental systems in 2D and 3D [9, 6, 23, 11, 24, 20]. It
should also be noted that active nematic dynamics are only weakly affected by
the bulk energy scale and mainly depend on the magnitude of $K_{el}$. The
second term on the left-hand side of eq. (1) is the co-rotational term
$\mathbf{\mathcal{W}}$ which describes the response of a flow-tumbling
director field to gradients in the flow field $\mathbf{u}$ [19],

$\mathbf{\mathcal{W}}_{ij}=\Omega_{ik}Q_{jk}+Q_{ik}\Omega_{jk}\>,$ (2)

where $\Omega_{ij}=(\partial_{j}u_{i}-\partial_{i}u_{j})/2$ is the anti-
symmetric part of the velocity gradient tensor. The time evolution of the
velocity field obeys the incompressible Navier-Stokes equation

$\rho\>D_{t}\mathbf{u}=\mathbf{\nabla}\cdot\left(\Pi^{visc}+\Pi^{el}+\Pi^{act}\right)\>,$
(3)

where $\rho$ is the density of the fluid and the stress tensor on the right-
hand side comprises the viscous stress $\Pi^{visc}$ of a Newtonian fluid,
elastic stress $\Pi^{el}$ caused by the relaxation of $\mathbf{Q}$ (for
details see Supplementary Information [25]), and an active stress component
$\Pi^{act}$. Activity locally induces dipolar stress, of magnitude $\zeta$,
along the axis of the director field, $\Pi^{act}_{ij}=-\zeta Q_{ij}$ which
drives the system out of equilibrium [26]. The parameter $\zeta$ quantifies
the type and magnitude of dipolar stress produced by the constituents, where
$\zeta>0$ describes extensile systems in which particles push out fluid along
their long axis $\mathbf{n}$ and pull in fluid in the plane perpendicular to
$\mathbf{n}$. $\zeta<0$ describes contractile activity where the flow
direction is reversed. It is apparent from eq. (3) that active forces are
induced by gradients in either $\mathbf{Q}$ or $\zeta$. They induce the well-
known hydrodynamic bend or splay instability of the director field
$\mathbf{n}$ in homogeneous active nematic bulk systems [26]. We solve the
equations of motion using a hybrid Lattice Boltzmann-finite difference method
[19] with parameter values stated in the Supplementary Information [25].

We consider periodic variations of activity along one coordinate axis with
alternating extensile ($\zeta>0$) and contractile activity ($\zeta<0$) in two
and three-dimensions,

$\zeta(\mathbf{x})=\zeta_{max}\cdot\cos(\mathbf{k}\cdot\mathbf{x})\>,$ (4)

with wave-vector $\mathbf{k}$ indicating the direction
$\mathbf{k}/|\mathbf{k}|$ and wavelength $\lambda=2\pi/|\mathbf{k}|$ of
activity patterns (Fig. 1 a,f). In the following we will denote the coordinate
axes parallel and perpendicular to the plane of constant activity as the
tangential and normal direction, respectively. We define a dimensionless
activity number $\mathcal{A}=\lambda\sqrt{|\zeta_{max}|/K_{el}}$ as the ratio
of the wavelength $\lambda$ to the active nematic length-scale
$\ell_{\zeta}=\sqrt{K_{el}/|\zeta_{max}|}$. For more complex, aperiodic
activity patterns $\mathcal{A}$ can be defined as the ratio of a
characteristic domain size $L$ to the active length-scale $\ell_{\zeta}$.

##  III  Active forces in the absence of defects

We start by summarising results for systems with heterogeneous activity in
which the magnitude of active stress is too small to create active defects.
Assuming gradients in the order parameter $\mathbf{Q}$ are small, the active
force density can be approximated by

$\mathbf{F}_{act}\approx
S\left(\mathbf{\nabla}\zeta-2\left(\mathbf{\nabla}\zeta\cdot\mathbf{n}\right)\mathbf{n}\right)\>,$
(5)

This approximation is justified as long as
$|\mathbf{\nabla}\zeta|\gg|\zeta~{}\mathbf{n}(\mathbf{\nabla}\cdot\mathbf{n})|$,
which is satisfied in the transition area between extensile and contractile
domains $\left(\zeta\approx 0\right)$ or in the absence of topological defects
in active regions which create strong distortions in the order parameter. The
force component along the tangential direction, perpendicular to the activity
gradient, is [27]

$F^{||}_{act}=-2|\mathbf{\nabla}\zeta|\cdot S\cos\theta\sin\theta\>,$ (6)

where $\theta\in[-\pi/2,\pi/2]$ denotes the angle of director field
$\mathbf{n}$ with respect to the activity gradient vector
$\mathbf{\nabla}\zeta$. Whenever $\mathbf{n}$ is not aligned parallel or
perpendicular to a finite activity gradient, active forces will generate flows
along the tangential direction, which are maximal at the steepest point of the
activity gradient (Fig. 1 b,d). As a consequence of the co-rotation term
$\mathbf{\mathcal{W}}$ in eq. (1), the resulting shear flow rotates the
director field $\mathbf{n}$ until a stable equilibrium is reached where the
elastic energy associated with director deformations balances co-rotation by
the flow.

In the case of $0<\theta<\pi/2$ and considering a $\mathbf{\nabla}\zeta>0$,
the resulting shear flow rotates the director clockwise left of the maximum of
$|\mathbf{\nabla}\zeta|$ and anti-clockwise to its right. Hence the director
in the contractile domain points along the normal, parallel to the activity
gradient $\mathbf{k}$, while in the extensile domain the director aligns
tangentially, perpendicular to the activity gradient (Fig. 1 c). The magnitude
of nematic alignment depends on the relative strengths of the flow-induced
rotation of the director field through the co-rotational term $\mathcal{W}$
and the relaxation of nematic distortions through the molecular field
$\mathbf{H}$. Thus spatial variations of $|\theta|$ along the normal
coordinate become larger for increasing activity $\zeta_{max}$ (faster flow)
and decreasing elastic constant $K_{el}$ (slower relaxation) (Fig. 1 e).

The nematic alignment induced by active flows at points of large activity
gradients therefore acts as an effective anchoring force, favouring normal
anchoring on the contractile side and tangential anchoring on the extensile
side of transition regions. To investigate how the shape of activity profiles
modifies the dynamics, we also used activity profiles which feature sharper
interfaces between extensile and contractile domains,

$\zeta(\mathbf{x})=\begin{cases}\zeta_{max},&\text{if }0\leq x\leq L/3,\\\
-\zeta_{max},&\text{if }L/2\leq x\leq 5/6L,\\\
\zeta_{max}\cdot\cos\left(\frac{6\pi}{L}x\right),&\text{otherwise}.\end{cases}\>$
(7)

We find that this does not qualitatively change the dynamics in the laminar
regime and only slightly modifies the shape of flow and director profiles
(Fig. S1 a,b [25]).

In 3D systems where activity varies along one of the coordinate axes, the
director alignment with respect to the gradient direction $\mathbf{k}$ is
quantified by the director angle
$|\theta_{x}|=|\text{cos}^{-1}(\mathbf{k}/|\mathbf{k}|\cdot\mathbf{n})|$. As
for the 2D case, the director field follows variations in activity and aligns
normal to the activity gradient in contractile domains. In extensile domains,
however, the tangential director orientation is degenerate and points along a
random direction within the two-dimensional plane perpendicular to
$\mathbf{k}$ (Fig. 1 f,g).

Figure 1: (a) Cosinusoidal variation of active stress along the normal
direction. Extensile and contractile regions are labeled as E and C
respectively. (b) Active flows form along the tangential direction in the
transition regions where $\zeta\approx 0$. (c) Director field aligns
tangential (normal) to the gradient direction $\mathbf{k}$ in extensile
(contractile) regions. (d,e) Before the onset of active turbulence, the
magnitude of flows and the degree of tangential (normal) alignment increases
with activity number $\mathcal{A}=\lambda\sqrt{|\zeta_{max}|/K_{el}}$. (f,g)
Three-dimensional systems in which activity varies along one coordinate axis
show the same quasi-two-dimensional behaviour. However, the director alignment
in extensile regions is degenerate within the two-dimensional tangential
plane.

##  IV  Defect polarization in 2D active turbulence

Figure 2: For sufficiently large activity, the system transitions into an
active turbulent state, which is characterised by (a) chaotic flow fields and
(b) motile topological defects. Comet-like $+1/2$ defects are self-propelled
and are shown in red; trefoil-like $-1/2$ defects are only passively advected
and are highlighted in blue. (c) The onset of active turbulence is
characterised by chaotic flows with components along the normal coordinate, as
indicated by a sharp increase of the average flow component along the normal
coordinate $\mathbf{u}_{N}$ relative to the root-mean-squared velocity
$\mathbf{u}_{rms}$ in the system. The critical activity $\zeta_{max}$ at which
this transition occurs depends on the dimensionless activity number
$\mathcal{A}=\lambda\sqrt{|\zeta_{max}|/K_{el}}\sim 10$. (d) The average
alignment of the director field in the turbulent regime decreases with
$\mathcal{A}$ as shown by the time-averaged director angle
$\langle|\theta|\rangle=\langle|\text{cos}^{-1}(\mathbf{k}/|\mathbf{k}|\cdot\mathbf{n})|\rangle$.
(e) The magnitude of alignment $A_{\theta}$, which is defined as the amplitude
of variations of $\langle|\theta|\rangle$ along the normal direction,
increases monotonically with $\mathcal{A}$ in the laminar regime and decays
with increasing activity in the turbulent regime. The mean director angle and
normal flow component were calculated by averaging quantities over the
tangential coordinate and over time $t=[0,80000]$ ($80$ snapshots).

We will now consider systems with large active stress $\zeta_{max}$, where the
soft confinement introduced by passive boundaries between extensile and
contractile regions is too weak to prevent the system transitioning to active
turbulence. Active turbulence is characterised by strong vorticity and
topological $\pm 1/2$ defects which are continually created and destroyed
(Fig. 2 a,b). In 2D systems the transition occurs at a critical activity
number $\mathcal{A}\gtrsim 10$, thus depending on the magnitude of active
stress $\zeta_{max}$, the elastic constant $K_{el}$ and the wavelength of the
activity patterns (Fig. 2 c). Although topological defects and chaotic flow
fields tend to destroy any nematic ordering, the time-averaged director field
still shows some finite normal (tangential) alignment in contractile
(extensile) domains in the turbulent state as shown by the time-averaged
director angle $\langle|\theta|\rangle$ between director $\mathbf{n}$ and
gradient axis (Fig. 2 d). The magnitude of alignment $A_{\theta}$, which is
defined as the amplitude of variations of $\langle|\theta|\rangle$ along the
normal direction, reaches a maximum at the onset of active turbulence and
decays with increasing activity (Fig. 2 e). The director $\mathbf{n}$
eventually approaches a random distribution,
$\theta\sim\mathcal{U}\left(-\pi/2,\pi/2\right)$, which yields
$\langle|\theta|\rangle\rightarrow\pi/4$ throughout the system.

In the active turbulent regime, both $\pm 1/2$ defects move in the background
flow, but comet-shaped $+1/2$ defects are also self-motile, with their
propulsion speed being proportional to activity $|\zeta|$ [28]. Their
propulsion direction $\mathbf{v}$ is parallel (anti-parallel) to the defect
orientation vector
$\mathbf{p}=\mathbf{\nabla}\cdot\mathbf{Q}/|\mathbf{\nabla}\cdot\mathbf{Q}|$
(which points from head to tail), in contractile (extensile) systems. Activity
gradients across $+1/2$ defects give rise an active torque $\Gamma_{act}$
which acts to align $+1/2$ defects parallel to activity gradients,
$\mathbf{p}\propto\mathbf{\nabla}\zeta$ (Fig. S3 a,b [25]). On the other hand,
the flow-driven normal (tangential) director field alignment in contractile
(extensile) domains acts as a soft boundary condition on the director field
$\mathbf{n}$. $+1/2$ defects which align anti-parallel to the activity
gradient, $\mathbf{p}\propto-\mathbf{\nabla}\zeta$, minimize distortions in
$\mathbf{n}$, thus elastic forces create an elastic torque $\Gamma_{el}$ which
minimize the elastic energy of the system (Fig. S3 c,d [25]). The elastic
torque $\Gamma_{el}$ on $+1/2$ defects thus opposes the active torque
$\Gamma_{act}$ (Fig. 3 a). While the magnitude of active and elastic torques
can be easily measured for single $+1/2$ defects subject to external activity
gradients or nematic boundary conditions, it is not immediately clear which
contribution dominates in the regime of active turbulence, where a finite
alignment magnitude $A_{\theta}$ of the average director field
$\langle|\theta|\rangle$ gives rise to soft boundary conditions and
interactions between defects become important.

To investigate the behaviour of defects in the regime of active turbulence, we
track the position and orientation $\mathbf{p}$ of $+1/2$ defects in the
system. By considering active and elastic forces in the vicinity of defects,
we can calculate elastic and active torques on $+1/2$ defects as a function of
defect orientation (see Supplementary Information [25]). By averaging over
defects and time, we confirm that elastic torques indeed oppose active torques
(Fig. 3 b). However, active torques remain the dominant contribution to the
total torque $\Gamma=\Gamma_{el}+\Gamma_{act}$ in the regime of active
turbulence (Fig. 3 c).

In order to quantify to which extent defects are locally polarized by active
torques, we measure the angular distribution of defect orientations
$\mathbf{p}$, averaging over time (Fig. S2 a [25]). We define the defect
orientation order parameter $\mathcal{P}=\sqrt{\langle
p_{x}\rangle^{2}+\langle p_{y}\rangle^{2}}$, which is is zero if $+1/2$
defects are oriented randomly and equals one if the defects perfectly align
along one direction. The polarization direction, $\langle\mathbf{p}\rangle$,
is calculated as the circular mean of the distribution of $\mathbf{p}$. Even
though the motion of individual defects is highly chaotic in the regime of
active turbulence, we find that $+1/2$ defects are strongly polarized at the
boundaries between extensile and contractile domains (Fig. 3 d). The average
defect orientation aligns with and is proportional to the activity gradient,
$\langle\mathbf{p}\rangle\propto\mathbf{\nabla}\zeta$, in agreement with
previous theoretical predictions which considered active torques on individual
defects [15, 16]. A strong polarization of defects along activity gradients
thus confirms once more that active torques on $+1/2$ defects are dominating
elastic torques, $\Gamma_{act}>\Gamma_{el}$. It should be emphasised, however,
that both torques originate from activity gradients since the nematic boundary
condition giving rise to $\Gamma_{el}$ is driven by active tangential flows.

In systems with sharp activity interfaces, non-zero defect polarization
$\langle\mathbf{p}\rangle$ and average director alignment
$\langle|\theta|\rangle$ are limited to boundaries and vanish in contractile
and extensile bulk regions as expected for isotropic active turbulence (Fig.
S1 c,d [25]).

Figure 3: (a) $+1/2$ defects are subject to an elastic torque $\Gamma_{el}$,
which tries to minimize the elastic energy of the system (orange), and an
active torque $\Gamma_{act}$, which is caused by differential self-propulsion
$\mathbf{v}$ across the defect (red). (b) Activity gradients create an active
torque on $+1/2$ defects, which is strongest if the defect orientation
$\mathbf{p}$ is perpendicular to the activity gradient $\mathbf{\nabla}\zeta$
($\alpha=\pm\pi/2$) and vanishes at $\alpha=0,\pi$. Active torques align
$+1/2$ defects parallel to activity gradients, which is evident from the
stable fixed-point at $\alpha=0$. Elastic torques, on the other hand, oppose
active torques and favour $+1/2$ defects to align anti-parallel to activity
gradients (see Fig. S3 [25]). We find that the average torque on defects
$\Gamma=\Gamma_{el}+\Gamma_{act}$ is dominated by active contributions in the
regime of active turbulence. Average torques $\Gamma_{act}$ and $\Gamma_{el}$
were obtained by averaging over defects within an orientation window
$[\alpha-\pi/12,\alpha+\pi/12]$. (c) The relative magnitude of active and
elastic torques does not significantly change as function of activity number
$\mathcal{A}$. Error bars show the mean and standard deviation of
$\Gamma_{el}/\Gamma_{act}$ when averaged over defect orientations $\alpha$.
(d) $+1/2$ defects are polarized at positions of large activity gradients
between contractile and extensile domains with their tail pointing towards
extensile domains. The fluctuations in polarization $\mathcal{P}$ are caused
by the finite sample size of defect orientations. The defect polarization
$\mathcal{P}$ as a function of normal coordinate $x$ was obtained by averaging
over defect orientations $\mathbf{p}$ in intervals
$[x-\lambda/20,x+\lambda/20]$ along $x$ and averaging over time $t=[0,400000]$
($200$ snapshots).

##  V  Disclination line properties in 3D active turbulence

Figure 4: (a) Three-dimensional active turbulence is characterised by motile
disclination lines which form closed loops or span the entire system as open
lines. (b) Disclination lines can be locally classified by a continuous twist-
angle $0<\beta<\pi$ which quantifies how the director field winds around the
disclination core in the plane perpendicular to the line tangent. $\beta$ is
the angle between the line tangent $\mathbf{t}$ and the rotation axis
$\mathbf{\Omega}$ around which the director winds while moving around the
disclination core. Quasi-two-dimensional $\pm 1/2$ defects correspond to
$\beta=\approx 0,\pi$, respectively, and $\beta\approx\pi/2$ describes twist
disclinations, which are unique to three-dimensional systems (radial twist-
type shown here). (c) The relative density of $\pm 1/2$ and twist segments
depends on the type of activity (extensile/contractile) and follows spatial
variations in active stress (fluctuations are caused by finite sample size).
(d) The director field alignment in 3D systems in the turbulent regime decays
with $\mathcal{A}=\lambda\sqrt{|\zeta_{max}|/K_{el}}$, as shown by the time-
averaged director angle $\langle\cos(\theta_{x})\rangle$. (e) As for 2D
systems, the magnitude of alignment $A_{\theta}$, which in 3D is defined as
the amplitude of variations of $\langle\cos(\theta_{x})\rangle$ along the
normal direction, decays with increasing activity number $\mathcal{A}$. Black
diamonds show $A_{\theta}/A_{\theta}^{max}$ as a function of $\zeta_{max}$ for
different parameters $K_{el}=[0.04,0.08,0.12,0.2]$, where $A_{\theta}^{max}$
is the maximum alignment realized just before the onset of active turbulence
(dashed line). Pink diamonds show normalized data taken from Fig. 2 e. The
mean director angle and disclination line densities were calculated by
averaging quantities over both tangential coordinates and time $t=[0,500000]$
($250$ snapshots).

We next consider 3D bulk systems in which activity varies along one direction.
The point defects of 2D nematics are replaced by disclination lines which can
continuously transform from a local $-1/2$ configuration, in the plane
perpendicular to the line tangent, into a $+1/2$ configuration through an
intermediate twist defect, which is inherently three-dimensional (Fig. 4 a,b).
The director configuration of the defects along a disclination line can be
locally classified by the twist angle $\beta$, where $-1/2$ ($+1/2$) wedge-
type defects correspond to $\beta=0$ $(\pi)$ and line segments with twist
defects are indicated by $\beta=\pi/2$. Radial twist and $+1/2$ defects both
have a unique polar direction $\mathbf{p}$ in the plane perpendicular to the
line tangent (Fig. 4 b).

The flows and morphological dynamics of disclination lines in 3D active
nematics are governed by the local director profile [12, 13]. While
disclination lines in purely extensile systems are predominantly twist-like,
contractile activity favours the creation of the wedge-like line segments with
cross-sections that resemble $\pm 1/2$ defects [29, 14]. This is reflected by
variations in the spatial density of $\pm 1/2$ and twist line segments across
a system with extensile and contractile domains (Fig. 4 c).

In 3D active turbulence with homogeneous activity the orientation of the
director field $\mathbf{n}$ is, by symmetry, randomly distributed,
$\cos\theta_{x}\sim\mathcal{U}(0,1)$, where $\theta_{x}$ is the angle between
the director field and a random reference vector. Deviations from the average
director angle $\langle\cos(\theta_{x})\rangle=1/2$ thus indicate that the
time-averaged director field is either aligned parallel
$\left(\langle\cos\theta_{x}\rangle>1/2\right)$ or perpendicular
$\left(\langle\cos\theta_{x}\rangle<1/2\right)$ with respect to the reference
vector. For 3D active turbulence with heterogeneous activity we find that the
time-averaged director field is aligned normal (tangential) with respect to
the activity gradient in contractile (extensile) domains, as for 2D systems
(Fig. 4 d). The nematic alignment in 3D systems can be quantified by
$A_{\theta}$, the amplitude of variations of $\langle\cos(\theta_{x})\rangle$
along the normal direction. To compare the variation of the alignment
magnitude with activity in 2D ($A_{\theta}\in[0,\pi/4]$) and 3D systems
($A_{\theta}\in[0,0.5]$), we normalize $A_{\theta}$ with respect to
$A_{\theta}^{max}$, the maximum alignment realized just before the onset of
active turbulence. We find that with increasing activity, the average
alignment of the director field in 3D decays much more slowly compared to 2D
systems, where $A_{\theta}$ sharply drops at the onset of active turbulence
(Fig. 4 e).

To gain further insight, we measure the angle $\gamma_{x}$ between
disclination line segments and the gradient direction $\mathbf{k}$. We find
that disclination lines are preferentially oriented perpendicular to activity
gradients, $\gamma_{x}\approx\pi/2$ (Fig. 5 a). This preferred alignment is
caused by active torques $\Gamma_{\gamma}$ on line segments when they lie
along activity gradients (Fig. 5 b). This implies that the orientation of the
defect structures in the plane perpendicular to line segments can be either
parallel ($\mathbf{p}\propto\mathbf{\nabla}\zeta$), anti-parallel
($\mathbf{p}\propto-\mathbf{\nabla}\zeta$) or perpendicular
($\mathbf{p}\cdot\mathbf{\nabla}\zeta=0$) to the activity gradient.

We define the polarization order parameter for disclination lines as
$\mathcal{P}=\sqrt{\langle p_{x}\rangle^{2}+\langle p_{y}\rangle^{2}+\langle
p_{z}\rangle^{2}}$, which is is zero if defects are oriented randomly and
equals one if the defects perfectly align along one direction. The
polarization order parameter $\mathcal{P}$ and direction of $+1/2$ and twist
segments are shown in Fig. 5 c and Fig. 5 d, respectively. While $+1/2$
disclination are polarized parallel to activity gradients,
$\langle\mathbf{p}\rangle\propto\mathbf{\nabla}\zeta$, like in 2D systems,
twist disclinations are polarized anti-parallel to the activity gradient
$\langle\mathbf{p}\rangle\propto-\mathbf{\nabla}\zeta$.

We argue that this is caused by a different balance between elastic torques
$\Gamma_{el}$ and active torques $\Gamma_{act}$ on different defect types.
While active torques $\Gamma_{act}$ dominate elastic torques $\Gamma_{el}$ for
$+1/2$ disclinations (see Fig. 3 b), active forces around twist disclinations
feature out-of-plane components, yielding an active torque which does not
align with the line tangent (Fig. S6 [25]). Thus the orientation $\mathbf{p}$
of twist defects is less affected by active torques than the orientation of
$+1/2$ defects (Fig. S6 c), causing twist defects to align in a way which
minimizes the elastic energy of the system. Indeed, since active torque is
generated by differential velocities across defects, we would expect a smaller
active torque on twist defects because they have a lower self-propulsion speed
than $+1/2$ defects [12, 14].

Sharp activity interfaces in 3D systems give rise to similar dynamics, where
defect polarization $\langle\mathbf{p}\rangle$ and director alignment
$\langle\cos\theta_{x}\rangle$ are limited to boundary regions and vanish in
contractile and extensile bulk regions (Fig. S4 [25]).

Figure 5: (a) The angle $\gamma_{x}$ between disclination line segments and
the normal axis along which activity varies is measured for all disclination
lines in the system. The time-averaged distribution of $|\cos(\gamma_{x})|$
shows that disclination line segments preferentially point perpendicular to
activity gradients ($\gamma_{x}=\pi/2$), independent of twist-angle $\beta$.
The distribution of orientation $\gamma_{x}$ was obtained by averaging over
all disclination line segments and time $t=[0,500000]$. (b) When a $+1/2$
disclination line points along an activity gradients $\mathbf{\nabla}\zeta$
($\gamma_{x}=0$), active forces (black arrows) create a torque
$\Gamma_{\gamma}$ (pink arrow), which re-orients the disclination line
perpendicular the gradient direction ($\gamma_{x}=\pi/2$). The magnitude of
the line torque $\Gamma_{\gamma}$ depends on the self-propulsion speed of the
disclination type, thus it is strongest for $+1/2$ line segments and vanishes
for $-1/2$ segments. Since disclination line segments of different types
$\beta$ are connected with each other, line torques $\Gamma_{\gamma}$ acting
on $+1/2$ and twist segments also align passive $-1/2$ segments. (c) As in 2D
systems, active torques align the orientation $\mathbf{p}$ of $+1/2$ defects
parallel to activity gradients $\mathbf{\nabla}\zeta$, causing a large
polarization of $+1/2$ disclinations ($\beta\in[7\pi/8,\pi]$) at points of
large activity gradients. (d) Twist disclinations ($\beta\in[3\pi/8,5\pi/8]$)
are dominated by elastic torque rather than active torque, thus they are
polarized anti-parallel to activity gradients. The fluctuations in
polarization $\mathcal{P}$ are caused by the finite sample size of $+1/2$ and
twist defect orientations. The polarization $\mathcal{P}$ of disclination
lines was calculated by averaging over defect orientations $\mathbf{p}$ in
intervals $[x-\lambda/20,x+\lambda/20]$ along the normal coordinate $x$ and
averaging over time $t=[0,500000]$ ($250$ snapshots).

##  VI  Conclusion

In this paper we have investigated how two- and three-dimensional active
nematic bulk systems respond to one-dimensional spatial variations of active
stress. We have shown that, in the absence of defects, activity gradients
induce an effective anchoring force which aligns the director field along the
normal direction in contractile domains and tangentially in extensile domains.
In line with recent analytical work [16, 17], we found that the orientation of
motile $+1/2$ defects in 2D active turbulence are dominated by active torque
and polarized parallel to activity gradients,
$\mathbf{p}\propto\mathbf{\nabla}\zeta$. Since the director field in the
vicinity of polarized $+1/2$ defects opposes the director alignment set up by
active flows perpendicular to activity gradients, the average alignment of the
director field in 2D active turbulence remains weak.

In 3D systems we find that disclination line segments are oriented
perpendicular to activity gradients, thereby allowing the orientation of the
defect structures in the plane perpendicular to line segments to point
parallel, anti-parallel or perpendicular to activity gradients. Like in 2D
systems, active torque robustly aligns $+1/2$ disclination lines parallel to
activity gradients, even in the turbulent regime. In contrast, the alignment
of twist disclinations is dominated by elastic interactions rather than active
torques, causing them to be polarized anti-parallel to activity gradients.
Thereby the director field in the vicinity of twist defects matches with the
alignment set up by active flows in the plane perpendicular to activity
gradients and the average alignment of the director field decays more slowly
with activity number $\mathcal{A}$ than in 2D systems. Deriving exact
expressions for the elastic torque on disclination line segments in the
turbulent regime remains challenging since the local torque on line segments
depends not only on the nematic boundary conditions in the plane perpendicular
to the line segment, but also on elastic forces acting along disclination
lines and their local curvature. Active-passive systems which are purely
extensile or contractile with spatial activity profiles
$\zeta(\mathbf{x})=\zeta_{max}(\cos(\mathbf{k}\cdot\mathbf{x})\pm 1)$ are
qualitatively similar to extensile-contractile systems and are presented in
the Supplementary Information [25] (Fig. S2,S5).

We note that the alignment of the director field and defects due to active
force gradients is not unique to bulk systems. In systems with isotropic-
nematic interfaces ($\mathbf{\nabla}S>0$), active forces create flows which
align the director field either perpendicular or planar to the interface in
contractile or extensile systems, respectively [30]. This effect, termed
active anchoring, also affects the orientation of defects close to an
interface and can be observed both in experiments [31] and in simulations in
two and three dimensions [32, 14].

Experiments with spatially structured activity have been realised by combining
actin filaments with light-activated gear-shifting myosin motors, where light
locally induces extensile active stress in an otherwise passive fluid [15].
Extending the concept of light-modulated activity to contractile systems, it
might be possible to create artificial systems with distinct extensile and
contractile domains in the future by selectively activating molecular motors
mediating either extensile or contractile stress. Spatial variations of active
stress are also ubiquitous in biological systems, ranging from structure
formation during embryonic development driven by a heterogeneous distribution
of myosin motors [33, 34, 35] to flows and deformations of cell aggregates,
such as tumors, caused by differential cell growth and death as a consequence
of limited access to nutrients in the core of the aggregate [36, 37, 38, 39].
Understanding the response of two- and three-dimensional active nematic
systems to activity gradients is thus an important step towards applying
physical theories in the life sciences.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

This project was funded by the European Commission’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No
812780.

## References

  * Jülicher _et al._ [2018] F. Jülicher, S. W. Grill and G. Salbreux, _Reports on Progress in Physics_ , 2018, 81, 076601.
  * Marchetti _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, _Reviews of Modern Physics_ , 2013, 85, 1143\.
  * Chen _et al._ [2018] T. Chen, T. B. Saw, R.-M. Mège and B. Ladoux, _Journal of Cell Science_ , 2018, 131, jcs218156.
  * Lee _et al._ [2021] G. Lee, G. Leech, M. J. Rust, M. Das, R. J. McGorty, J. L. Ross and R. M. Robertson-Anderson, _Science Advances_ , 2021, 7, eabe4334.
  * Yaman _et al._ [2019] Y. I. Yaman, E. Demir, R. Vetter and A. Kocabas, _Nature Communications_ , 2019, 10, 1–9.
  * Saw _et al._ [2017] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans and B. Ladoux, _Nature_ , 2017, 544, 212–216.
  * Mueller _et al._ [2019] R. Mueller, J. M. Yeomans and A. Doostmohammadi, _Physical Review Letters_ , 2019, 122, 048004.
  * Duclos _et al._ [2017] G. Duclos, C. Erlenkämper, J.-F. Joanny and P. Silberzan, _Nature Physics_ , 2017, 13, 58–62.
  * Saw _et al._ [2018] T. B. Saw, W. Xi, B. Ladoux and C. T. Lim, _Advanced Materials_ , 2018, 30, 1802579.
  * Vromans and Giomi [2016] A. J. Vromans and L. Giomi, _Soft Matter_ , 2016, 12, 6490–6495.
  * Duclos _et al._ [2020] G. Duclos, R. Adkins, D. Banerjee, M. S. Peterson, M. Varghese, I. Kolvin, A. Baskaran, R. A. Pelcovits, T. R. Powers, A. Baskaran _et al._ , _Science_ , 2020, 367, 1120–1124.
  * Binysh _et al._ [2020] J. Binysh, Ž. Kos, S. Čopar, M. Ravnik and G. P. Alexander, _Physical Review Letters_ , 2020, 124, 088001.
  * Čopar _et al._ [2019] S. Čopar, J. Aplinc, Ž. Kos, S. Žumer and M. Ravnik, _Physical Review X_ , 2019, 9, 031051.
  * Ruske and Yeomans [2021] L. J. Ruske and J. M. Yeomans, _Physical Review X_ , 2021, 11, 021001\.
  * Zhang _et al._ [2021] R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant, M. L. Gardel _et al._ , _Nature Materials_ , 2021, 20, 875–882.
  * Shankar and Marchetti [2019] S. Shankar and M. C. Marchetti, _Physical Review X_ , 2019, 9, 041047\.
  * Tang and Selinger [2021] X. Tang and J. V. Selinger, _Physical Review E_ , 2021, 103, 022703\.
  * Mozaffari _et al._ [2021] A. Mozaffari, R. Zhang, N. Atzin and J. J. de Pablo, _Physical Review Letters_ , 2021, 126, 227801.
  * Marenduzzo _et al._ [2007] D. Marenduzzo, E. Orlandini and J. M. Yeomans, _Physical Review Letters_ , 2007, 98, 118102.
  * Schimming and Viñals [2020] C. D. Schimming and J. Viñals, _Physical Review E_ , 2020, 102, 010701.
  * Long _et al._ [2021] C. Long, X. Tang, R. L. Selinger and J. V. Selinger, _Soft Matter_ , 2021, 17, 2265–2278.
  * Beris and Edwards [1994] A. N. Beris and B. J. Edwards, _Thermodynamics of flowing systems with internal microstructure_ , Oxford University Press, New York ; Oxford, 1994\.
  * Norton _et al._ [2018] M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran and M. F. Hagan, _Physical Review E_ , 2018, 97, 012702.
  * Zhang _et al._ [2018] R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel and J. J. De Pablo, _Proceedings of the National Academy of Sciences_ , 2018, 115, E124–E133.
  * [25] _See Supplementary Information for supplementary figures and details on numerical implementation._
  * Simha and Ramaswamy [2002] R. A. Simha and S. Ramaswamy, _Physical Review Letters_ , 2002, 89, 058101.
  * Bhattacharyya and Yeomans [2021] S. Bhattacharyya and J. M. Yeomans, _Soft Matter_ , 2021, 17, 10716–10722.
  * Giomi [2015] L. Giomi, _Physical Review X_ , 2015, 5, 031003.
  * Nejad and Yeomans [2022] M. R. Nejad and J. M. Yeomans, _Physical Review Letters_ , 2022, 128, 048001.
  * Blow _et al._ [2014] M. L. Blow, S. P. Thampi and J. M. Yeomans, _Physical Review Letters_ , 2014, 113, 248303.
  * Dell’Arciprete _et al._ [2018] D. Dell’Arciprete, M. Blow, A. Brown, F. Farrell, J. S. Lintuvuori, A. McVey, D. Marenduzzo and W. C. Poon, _Nature Communications_ , 2018, 9, 1–9.
  * Doostmohammadi _et al._ [2016] A. Doostmohammadi, S. P. Thampi and J. M. Yeomans, _Physical Review Letters_ , 2016, 117, 048102.
  * Streichan _et al._ [2018] S. J. Streichan, M. F. Lefebvre, N. Noll, E. F. Wieschaus and B. I. Shraiman, _Elife_ , 2018, 7, e27454.
  * He _et al._ [2014] B. He, K. Doubrovinski, O. Polyakov and E. Wieschaus, _Nature_ , 2014, 508, 392–396.
  * Mason _et al._ [2013] F. M. Mason, M. Tworoger and A. C. Martin, _Nature Cell Biology_ , 2013, 15, 926–936.
  * Jagiella _et al._ [2016] N. Jagiella, B. Müller, M. Müller, I. E. Vignon-Clementel and D. Drasdo, _PLoS Computational Biology_ , 2016, 12, e1004412.
  * Delarue _et al._ [2013] M. Delarue, F. Montel, O. Caen, J. Elgeti, J.-M. Siaugue, D. Vignjevic, J. Prost, J.-F. Joanny and G. Cappello, _Physical Review Letters_ , 2013, 110, 138103.
  * Delarue _et al._ [2014] M. Delarue, J.-F. Joanny, F. Jülicher and J. Prost, _Interface Focus_ , 2014, 4, 20140033.
  * Ranft _et al._ [2010] J. Ranft, M. Basan, J. Elgeti, J.-F. Joanny, J. Prost and F. Jülicher, _Proceedings of the National Academy of Sciences_ , 2010, 107, 20863–20868.


# A Deployable Quantum Access Points Selection Algorithm for Large-Scale
Localization

Ahmed Shokry  Computer Science and Engineering
Pennsylvania State University
PA, USA
<EMAIL_ADDRESS>Moustafa youssef  Computer Science and Engineering
American University in Cairo
Cairo, Egypt
<EMAIL_ADDRESS>
###### Abstract

Effective access points (APs) selection is a crucial step in localization
systems. It directly affects both localization accuracy and computational
efficiency. Classical APs selection algorithms are usually computationally
expensive, hindering the deployment of localization systems in a large
worldwide scale.

In this paper, we introduce a quantum APs selection algorithm for large-scale
localization systems. The proposed algorithm leverages quantum annealing to
eliminate redundant and noisy APs. We explain how to formulate the APs
selection problem as a quadratic unconstrained binary optimization (QUBO)
problem, suitable for quantum annealing, and how to select the minimum number
of APs that maintain the same overall localization system accuracy as the
complete APs set. Based on this, we further propose a logarithmic-complexity
algorithm to select the optimal number of APs.

We implement our quantum algorithm on a real D-Wave Systems quantum machine
and assess its performance in a real test environment for a floor localization
problem. Our findings reveal that by selecting fewer than 14% of the available
APs in the environment, our quantum algorithm achieves the same floor
localization accuracy as utilizing the entire set of APs and a superior
accuracy over utilizing the reduced dataset by classical APs selection
counterparts. Moreover, the proposed quantum algorithm achieves more than an
order of magnitude speedup over the corresponding classical APs selection
algorithms, emphasizing the efficiency of the proposed quantum algorithm for
large-scale localization systems.

###### Index Terms:

quantum location determination, access point selection, quantum annealing,
QUBO, next-generation location tracking systems

##  I Introduction

Large-scale fingerprinting localization systems are commonly employed both
indoors [1, 2] and outdoors [3, 4] due to their high accuracy. These systems
utilize the Radio Frequency (RF) technology and follow a two-stage process. In
the offline phase, they create an “RF fingerprint” by capturing the received
signal strength (RSS) signature from different access points (APs) at
different locations in the area of interest and store this information in a
database. These APs can include WiFi access points, cell towers, or Bluetooth
beacons. In the online tracking phase, the RSS received from the APs at the
user’s device with unknown location is compared with the fingerprint samples
in the database to identify the location closest, in the RSS space, to the
unknown location.

As the fingerprinting localization techniques need to compare online RSS
measurements with offline ones at every fingerprint location, their time and
space complexity increase with the increase in number of APs. Classical
algorithms try to infer a promising reduced subset of APs by adding or
removing APs and comparing the results [5]. But ultimately, when the number of
APs is large (especially in a worldwide localization system), this exhaustive
exploration of all combinations of APs becomes computationally prohibitive and
a trade-off remains between the accuracy of selecting the effective APs and
the effort spent finding the best subset of APs. This complexity does not
allow for scaling of classical APs selection techniques to support large-scale
indoor/outdoor localization systems on a global scale, e.g. in IoT
environments where there may be a significant number of APs with wireless
interfaces.

Recently, there have been successful proposals of quantum fingerprint matching
algorithms in the literature [6, 7, 8, 9, 10, 11, 12]. These techniques
leverage universal gate-based quantum machines to achieve exponential saving
in matching time and space in terms of the number of APs and in terms of the
number of fingerprint locations [13]. However, due to quantum decoherence
[14], the current gate-based quantum machines can only handle a limited number
of qubits accurately, which restricts the deployment of large-scale
localization systems in the near future on them.

In this paper, we propose a novel approach for APs selection that can help in
reducing the computational power of fingerprinting-based localization
algorithm. Specifically, we propose a quantum APs selection algorithm that
leverages quantum annealing to eliminate redundant and noisy APs which results
in improving the overall localization accuracy. We present a formulation of
the APs selection problem as a Quadratic Unconstrained Binary Optimization
(QUBO) problem [15, 16], suitable for processing on quantum annealers, with
the aim of selecting the powerful APs that improves overall localization
system efficiency. The proposed quantum algorithm utilizes binary variables to
denote APs selection and formulates an objective function that tries to select
APs that have a strong effect on the location information (i.e. maximize
important APs) and have a weak correlation between each others (i.e. minimize
redundancy). Furthermore, we present an algorithm designed to identify the
minimum set of powerful APs necessary to enhance localization accuracy. Our
proposed algorithms leverages the real D-Wave Systems quantum machines, which
are specialized practical quantum computers for solving quantum annealing
problems [17].

Results from deploying the proposed quantum algorithm and evaluating its
performance in a real testbed for a floor identification task reveal that by
selecting fewer than 14% of the available APs in the environment, our quantum
algorithm achieves the same floor localization accuracy as utilizing the
entire set of APs and a superior accuracy over utilizing the reduced dataset
by classical APs selection counterparts. Moreover, the proposed quantum
algorithm achieves a speedup of more than an order of magnitude over the
corresponding classical APs selection algorithms.

The rest of the paper is organized as follows: Section II gives a background
on quantum computing and the floor localization problem. Section III provides
the details of our quantum APs selection algorithm. We present the details of
the implementation and evaluation results in Section IV. Finally, Section V
concludes the paper.

##  II Background

In this section, we start by a background on quantum computing and quantum
annealing. Then, we give a brief introduction to the floor estimation problem,
which is the case study we apply our APs selection algorithm to.

###  II-A Quantum Computing and Quantum Annealing

Quantum computing represents a paradigm shift in computational theory and
practice, harnessing the principles of quantum mechanics to enable the
solution of complex problems that classical computers struggle with. At its
core, quantum computing leverages quantum bits (qubits) that, unlike classical
bits, can exist in a superposition state of 0 and 1 at the same time,
exponentially increasing computational capacity [18].

Quantum entanglement is another phenomenon where two or more qubits become
correlated and the state of one qubit instantaneously influences the state of
another, regardless of the distance between them [18]. This entanglement
property enables the creation of highly-correlated quantum states, a
phenomenon that is fundamental to the efficiency of quantum algorithms. The
power of quantum computing lies in its ability to perform parallel
computations due to superposition and to exploit quantum entanglement for
correlations between qubits.

Quantum annealing is a specialized approach within quantum computing [19],
tailored for solving optimization problems which involve finding the best
solution from a vast number of possibilities. Quantum annealing has been used
in various fields including finance [20] and artificial intelligence [21]. It
follows the principles of adiabatic quantum computation [19], where a system
is slowly evolved from a simple Hamiltonian (i.e. energy function) to a
complex one, reaching the ground state that encodes the solution. Adiabatic
quantum computation involves initializing a quantum system in the ground state
of a simple Hamiltonian ($H_{0}$), and slowly transforming it into the ground
state of a more complex final Hamiltonian ($H_{f}$). This adiabatic evolution
allows the system to settle into the state representing the optimal solution
to an optimization problem. The simple Hamiltonian $H_{0}$ represents the
initial state of the quantum system. It is carefully chosen such that its
ground state is easy to prepare, typically a state where each qubit is in a
simple state. The final Hamiltonian $H_{f}$ encodes the optimization problem.
Its ground state contains information about the optimal solution to the
problem. The adiabatic quantum computation ensures that the system evolves
into this final ground state, capturing the optimal solution.

The evolution of a quantum system during annealing is described by the
Schrödinger equation [22]. For a quantum state $|\psi(t)\rangle$, the time-
dependent Schrödinger equation is given by:

$i\hbar\frac{\partial}{\partial t}|\psi(t)\rangle=H(t)|\psi(t)\rangle$ (1)

Here, $i$ is the imaginary unit, $\hbar$ is the reduced Planck constant,
$H(t)$ is the time-dependent Hamiltonian operator representing the system’s
energy, and $\ket{\psi(t)}$ is the state of the quantum system at time $t$.
The Hamiltonian operator is typically expressed as a combination of the simple
Hamiltonian ($H_{0}$) and the final Hamiltonian ($H_{f}$). The annealing
schedule is controlled by a parameter $s(t)$, which varies from 0 to 1 as the
system evolves. The time-dependent Hamiltonian is then defined as:

$H(t)=(1-s(t))H_{0}+s(t)H_{f}$ (2)

The objective is to start with the simple Hamiltonian ($H_{0}$) and gradually
transition to the final Hamiltonian ($H_{f}$) to guide the system towards the
optimal solution. This process exploits quantum phenomena such as
superposition and entanglement to explore a vast solution space more
efficiently than classical methods.

The previous construction plays an important role not only in the realm of
quantum mechanics, but also in the field of optimization. Nowadays, many real-
world optimization problems can be expressed as QUBO problems [15, 16], a
class of problems that involves finding the minimum (or maximum) of a
quadratic objective function subject to binary constraints. Formally, a QUBO
problem can be defined as:

$\text{Minimize }f(\mathbf{x})=\sum_{i,j}Q_{ij}x_{i}x_{j}+\sum_{i}P_{i}x_{i}$
(3)

where $\mathbf{x}$ is a binary vector, $Q_{ij}$ are coefficients representing
the quadratic terms, and $P_{i}$ are coefficients representing the linear
terms. QUBO problems can be reformulated as an Ising Hamiltonian, which is a
mathematical model that provides a simplified representation of magnetic
interactions in a system of atomic spins [23]. The basic idea is to map the
QUBO problem onto a physical system described by the Ising model, where the
spins of magnetic particles represent variables and their interactions encode
the problem’s constraints and objectives. In this mapping, binary variables
correspond to spins, and the QUBO objective function maps onto the energy
function of the Ising model (i.e. Hamiltonian). This transformation allows
quantum annealers, which naturally operate with Ising models, to be employed
in solving QUBO problems. Quantum annealers offer a practical avenue for
solving QUBO problems. The encoding of the QUBO problem into the Ising model
allows quantum annealers to explore the solution space through quantum
fluctuations, seeking the optimal configuration.

###  II-B The Floor Estimation Problem

The floor estimation problem [24, 25, 26, 27] relies on utilizing pre-
installed signal sources (e.g., WiFi APs or cellular base stations) to
determine users’ floor. Typically, the floor estimation algorithms operate in
two phases: (a) Offline calibration phase: During this phase, signal
information such as RSS and signal source IDs are collected at various known
locations on each floor to construct a fingerprint. Each fingerprint is stored
as vectors, where each entry represents the RSS from one AP, along with the
corresponding floor number. (b) Online tracking phase: In this phase, the
information from online signal sources is compared to the fingerprint vectors,
and the closest floor in the signal space is reported as the estimated floor.

##  III The Proposed Quantum APs Selection Algorithm

TABLE I:  Table of notations. Symbol | Description
---|---
$m$ | Number of fingerprint samples
$n$ | Number of APs
$r$ | A fingerprint RSS sample vector
$b$ | Number of bins in the RSS range
$l$ | A fingerprint floor label
$f$ | Number of distinct locations in the fingerprint
$\mathcal{D}=\\{(r,l)\\}$ | Original fingerprint
$\mathcal{D}_{S}=\\{(r_{S},l)\\}$ | Reduced fingerprint after selecting effective S APs
$I_{i}$ | Importance of AP $i$
$R_{ij}$ | Mutual correlation between two APs (redundancy)
$\alpha$ | Balancing parameter (importance and redundancy)
$\bm{x}^{*}$ | A binary vector indicating the selected APs
$k=\lVert\bm{x}^{*}\rVert_{1}$ | Number of selected APs

In this section, we describe our proposed APs selection algorithm. We start
with formulating APs selection problem as a QUBO problem. Then, we describe
how to select the minimum number of APs that achieve the best localization
accuracy. Based on this, we present our APs selection algorithm. Table I
summarizes the notations used in this section.

###  III-A QUBO Formulation

(a) Importance

(b) Redundancy

Figure 1: Importance and Redundancy in the full dataset. The importance
reflects the influence of each AP on estimating the location. The Redundancy
reflects APs inter-dependency.

Assume a localization problem on a fingerprint with $m$ fingerprint samples
and $n$ APs, $\mathcal{D}=\\{(r^{i},l^{i})\\}_{i\in\\{1,2,\dots,m\\}}$ with
$n$-dimensional APs received signal strength (RSS) vector
$r^{i}\in\mathbb{R}^{n}$ and floor information labels $l^{i}\in\mathbb{N}$ for
all $i\in\\{1,2,3,\dots,m\\}$. The APs selection problem tries to find a
subset of APs $S\subset\\{1,2,3,\dots,n\\}$, such that the reduced dataset
$\mathcal{D}_{S}=\\{(r^{i}_{S},l^{i})\\}_{i\in\\{1,2,3,\dots,m\\}}$ with the
reduced dimension vector $r_{S}$ referring to the RSS from the selected APs
only; leads to comparable localization accuracy as the full original RSS
vectors $r$ with all APs. Typically, this subset is found by searching for the
“optimal” subset of APs that lead to the best localization accuracy.

Our APs selection problem goal is to select APs that have the highest
influence on estimating the location (i.e., highest importance). It also tries
to reduce the APs that depend on each others (i.e., lowest redundancy) as
shown in Figure 1. We propose a quantum algorithm that targets maximizing the
importance and minimize the redundancy of APs. Specifically, we represent the
APs subset $S$ as a binary indicator vector
$\bm{x}=(x_{1},\dots,x_{n})\in\\{0,1\\}^{n}$, such that AP $i\in S$ if and
only if $x_{i}=1$ for all $i\in\\{1,2,3,\dots,n\\}$. Therefore, the optimal
subset of APs, $x^{*}$ can be obtained as

$\bm{x}^{*}=\underset{\bm{x}\in\\{0,1\\}^{n}}{\arg\min}\,Q(\bm{x},\alpha)$ (4)

Where our QUBO objective function, $Q(\bm{x},\alpha)$, is defined as

$Q(\bm{x},\alpha)=-\alpha\sum_{i=1}^{n}I_{i}x_{i}+(1-\alpha)\sum_{i,j=1}^{n}R_{ij}x_{i}x_{j},$
(5)

where the first sum of linear terms represents the individual contributions
from APs (i.e., importance, $I_{i}$) and the second sum contains quadratic
interaction terms between APs pairs (i.e., redundancy, $R_{ij}$). The user-
defined parameter $\alpha\in[0,1]$ balances the influence of the two terms:
the importance and the redundancy. Finding the set of binary variables
$\bm{x}_{i}$ that minimizes this objective function is equivalent to finding
the optimal set of APs.

In the balance of this section, we will give the details of how Importance and
Redundancy are calculated in our algorithm. Finally, we will provide an
algorithm for controlling the number of selected APs.

####  III-A1 Importance

We define the importance of each AP $i$ as the Cramér’s Variation (Cramér’s V)
measure of the signal strength association with the floor labels [28]. This
measure is an extension of chi-squared test, which is used to determine if
there is a significant association between two categorical variables. The
Cramer’s V measure further standardizes the chi-squared statistic to a value
between 0 and 1, making it easier to interpret; regardless of the specific
variables distributions or scales.

Hence, our Importance measure ranges from 0 to 1, where 0 indicates no
association between the AP RSS and floor numbers, and 1 indicates a strong
association.

Specifically, the Importance between an AP $i$ RSS values and the fingerprint
location labels $l$ can then be calculated as

$I_{i}=\sqrt{\frac{\chi^{2}}{m\times\min(\text{{b-1}},\text{{f-1}})}}$ (6)

Where $\chi^{2}$ is the chi-square statistic which measures the discrepancy
between the observed frequencies of AP $i$ RSS with the output labels, and the
frequencies that would be expected if the two variables were independent. $m$
is the number of fingerprint samples, $f$ is the number of distinct labels,
and $b$ is the number of possible RSS values.

Figure 1(a) shows the importance of all APs within the environment. The figure
shows that the fingerprint contains many APs with low importance. These APs
can be removed to reduce computational complexity without affecting the
accuracy of the localization system.

####  III-A2 Redundancy

We define the redundancy between each pair of APs, $i$ and $j$, as the
strength of the absolute correlation between them.

$R_{ij}=\text{$|$Corr(}r_{i},r_{j}\text{)$|$}\geq 0$ (7)

where $\bm{R}\in\mathbb{R}^{n\times n}$ is symmetric and positive semi-
definite matrix. This matrix represents the absolute mutual correlation
$|Corr(r_{i},r_{j})|$ among the APs RSS and therefore measures their
redundancy. For $i=j$, we set $R_{ii}=0$. Without loss of generality, we use
Pearson correlation as the measure of redundancy (Corr(.)) [29].

Figure 1(b) shows the correlation between each pair of APs within the
environment. The figure shows that the fingerprint contains many APs pairs
that depend on each other (color towards yellow). This also adds to the
computational cost on the localization system without a real gain on accuracy.

###  III-B Objective Function in Matrix Form

The QUBO objective function in Eq 5 can be written a matrix form
$\mathbf{P}(\alpha)\in\mathbb{R}^{n\times n}$ for better mathematical
manipulation, where a diagonal item $\mathbf{P}_{ii}$ presents the importance
of each AP $i$ (i.e. $I_{i}$), and a non-diagonal item $\mathbf{P}_{ij}$ is
the redundancy (i.e. $R_{ij}$). Hence, the optimal subset of APs in Eq 4 can
be expressed in matrix form as

$\bm{x}^{*}=\underset{\bm{x}\in\\{0,1\\}^{n}}{\arg\min}\
\bm{x}^{T}\mathbf{P}(\alpha)\bm{x}$ (8)
$P_{ij}(\alpha)=R_{ij}-\alpha(R_{ij}+\delta_{ij}I_{i}),$ (9)

where $\delta_{ij}=1_{\\{i=j\\}}$ denotes the Kronecker delta.

###  III-C Controlling the Number of Selected APs

In addressing our APs selection problem, the objective is to choose the
minimum number of APs that yield a certain localization accuracy. This
requirement can be incorporated by introducing a constraint to Eq 4, ensuring
the selection of precisely $k$ APs. One common strategy is to append a
Lagrangian penalty term [30], such as
$\lambda\left(\left(\sum_{i}x_{i}\right)-k\right)^{2}$, to $Q(\bm{x},\alpha)$.
This penalty term is only minimized when exactly $k$ APs are selected, where
$\lambda$ controls how important it is to exactly match the constraint [30].
Solving the problem yields the desired solution. However, two primary
challenges emerge from this approach. First, determining an appropriate value
for $\lambda$ is not a straightforward. Second, it is uncertain whether a
manually-selected $k$ will lead to the best localization accuracy [31].

Facing these challenges, we present an alternative approach for determining
the number of selected APs. Instead of relying on a penalty term, we exploit
the intrinsic connection between the parameter $\alpha$, controls the balance
between importance and redundancy, and the number of APs in the solution
($\lVert\bm{x}^{*}\rVert_{1}$). This premise is substantiated by considering
the extreme values of $\alpha\in[0,1]$. Specifically, setting $\alpha=0$
renders all diagonal entries of $\mathbf{P}(0)$ as $0$, emphasizing
redundancy. Naturally, both the empty set of features and any single AP are
minimally redundant, resulting in the optimal selection of APs being either
$\\{\emptyset\\}$ or $\\{\\{i\\};i\in\\{1,2,3,\dots,n\\}\\}$. Conversely, when
$\alpha=1$, the problem becomes linear with solely negative coefficients,
leading to the selection of all APs as the optimal solution. Consequently, by
iteratively adjusting $\alpha$ from $0$ to $1$ in sufficiently small
increments, we observe a monotonic increase of the number of selected APs in
the optimal solution in steps of one from $0$ to $n$. This insight implies
that we can modify $\alpha$ to obtain a subset of any desired size
$k\in\\{0,1,\dots,n\\}$ as we quantify in Section IV. The following lemma
proves that we can control the number of selected APs using $\alpha$.

######  Lemma 1.

For all $Q(\cdot,\alpha)$ defined as in Eq 5 and $k\in\\{0,1,\dots,n\\}$,
there exist an $\alpha\in[0,1]$ such that
$\bm{x}^{*}\in\arg\min_{\bm{x}}Q(\bm{x},\alpha)$ and
$\norm{\bm{x}^{*}}_{1}=k$. Moreover, the number of selected APs $k$
monotonically increases with $\alpha$.

###### Proof.

The proof can be found in Appendix A. ∎

This lemma shows that we do not need additional constraints on the QUBO
instance to control the number of features present in the global optimum.
Hence, we can determine a suitable value for $\alpha$ ($\alpha^{*}$) that
results in the minimum number of selected APs ($k$) in $\mathcal{O}(\log n)$
steps using binary search. It is important to note that the value of
$\alpha^{*}$ may not be singular. The outlined algorithm is provided in Alg 1.

Algorithm 1 Binary Search for $k$ Optimal APs

Input :

Redundancy Matrix $\bm{R}$.

Importance vector $\bm{I}$.

Number of APs $n$.

Termination threshold $\epsilon$.

Required Number of APs $k$.

Output :

Value of $\alpha^{*}$, Optimal value for $\alpha$.

Selected APs $\bm{x}^{*}$ with $\lVert\bm{x}^{*}\rVert_{1}=k$.

// Initialization

$\bm{x}\leftarrow ones(n)$

base_acc $\leftarrow$ Localizer($\bm{x}$) $\triangleright$ Localization with
full dataset

acc $\leftarrow$ base_acc

$a\leftarrow 0.0$ $\triangleright$ Initialize $\alpha$ range min value

$b\leftarrow 1.0$ $\triangleright$ Initialize $\alpha$ range max value

// Search for Optimal APs

do

prev_acc $\leftarrow$ acc $\alpha\leftarrow\frac{(a+b)}{2}$ $\triangleright$
$\alpha$ is a floating point number. $\bm{x}^{*}\leftarrow
Q^{*}(\bm{x},\alpha)$ $\triangleright$ Using quantum/classical annealing
$k^{\prime}\leftarrow\norm{\bm{x}^{*}}_{1}$ $\triangleright$ Number of
selected APs acc $\leftarrow$ Localizer($\bm{x}^{*}$) $\triangleright$
Localization accuracy if  _acc $<$ base_acc_ then

$a\leftarrow\alpha$ $\triangleright$ Increase selected APs

else

$b\leftarrow\alpha$

end if

while  _acc - prev_acc $\geq\epsilon$_

$\alpha^{*}\leftarrow\alpha$ return _$\alpha^{*}$ , $\bm{x}^{*}$_

##  IV Evaluation

In this section, we evaluate our quantum APs selection algorithm and implement
it on a real D-wave quantum hybrid solver which leverages the D-Wave Advantage
quantum processing unit (QPU) [32, 17] with 5000 qubits [33]. We start by
describing our real testbed. Then, we show the effect of the balancing
parameter $\alpha$. Finally, we compare the proposed quantum algorithm with
the classical counterpart.

###  IV-A Experiment Setup

The collected data cover a building with five floors of a public dataset [34].
The floors area is around $110\textrm{m}^{2}$ [34]. The total number of APs is
520. We used the already deployed WiFi APs in the building as well as the
overheard APs from the neighbors buildings. The total number of samples used
in the fingerprint is 10000 samples. Without loss of generality, we used the
Random Forest Classifier as a classification model for floor localization
[35]. The classical machine has a 2.2GHz CPU and 12GB RAM.

###  IV-B Effect of Parameter $\alpha$

Figure 1 shows the effect of the balancing parameter $\alpha$ on the number of
selected APs and the floor localization accuracy. The figure confirms that the
$\alpha$ parameter can control the number of selected APs as we proved in
Section III-C.

Figure 2: Effect of parameter $\alpha$ on the number of selected APs (red
line) and floor localization accuracy (blue bars). Figure 3:  Selected number
of APs. Figure 4:  Floor localization accuracy using the classical simulated
annealing algorithm and the quantum annealing algorithm. Figure 5: Running
time in log-scale of classical and quantum APs selection algorithms.

###  IV-C Comparison with Classical APs Selection

In this section, we compare our proposed quantum annealing APs selection
algorithm to the classical simulated annealing algorithm [36] under different
metrics: number of selected APs, accuracy, and running time. We end the
section by showing the importance and redundancy of the selected APs for both
algorithms.

####  IV-C1 Selected number of APs

Figure 3 shows the optimal number of APs obtained by running Algorithm 1 using
the classical simulated annealing algorithm and the proposed quantum annealing
algorithm. The figure shows that we can achieve the same accuracy as the
complete set of APs (520 APs) using 119 APs with simulated annealing and only
70 APs with the proposed quantum annealing algorithm. The figure also confirms
that the proposed quantum annealing algorithm can reach better effective APs
due to quantum tunneling [37], where the quantum annealing process can
penetrate through a potential barrier, enabling exploring the solution space
efficiently.

####  IV-C2 Accuracy

Figure 4 illustrates the classical simulated annealing floor accuracy (using
the selected 119 APs), the proposed quantum annealing accuracy (using the
selected 70 APs), and the accuracy of using the complete set of APs (520 APs).
Although the three sets of used APs achieve a comparable accuracy, the
proposed quantum algorithm achieves a superior floor localization accuracy
compared to the simulated annealing accuracy and even the accuracy while using
all APs. This highlight the ability of the proposed quantum algorithm to
select the set of effective APs that maximize the importance and minimize the
redundancy. This helps also in removing noisy APs, which is one reason why the
set of APs selected by the quantum annealing algorithm performs better than
the full set of APs.

####  IV-C3 Running time

Figure 5 further shows the running time of the classical and the quantum
version of Algorithm 1. The figure highlights that the quantum APs selection
algorithm can provide more than an order of magnitude saving in running time
than the classical counterpart.

####  IV-C4 Importance and Redundancy

Figure 6 shows the full set of APs importance and redundancy. The figure shows
that the dataset has many redundant (color more towards yellow at off-diagonal
elements) and irrelevant APs (low importance). Figures 7 and 8 further show
the quality of the selected APs using classical simulated annealing and
quantum annealing, respectively. The figures confirms that the proposed
quantum annealing algorithm capture more of the best APs. This is reflected by
showing higher Importance of the selected APs and less yellowish color on the
off-diagonal elements in the redundancy matrix.

(a) Importance

(b) Redundancy

Figure 6: Importance and Redundancy in all APs dataset (redrawn for
convenience).

(a) Importance

(b) Redundancy

Figure 7: Importance and Redundancy of the 120 selected APs by the simulated
annealing.

(a) Importance

(b) Redundancy

Figure 8: Importance and Redundancy of the 70 selected APs by the proposed
quantum algorithm.

##  V Conclusion

In this paper, we introduced a quantum APs selection algorithm for large-scale
localization systems. The proposed quantum algorithm leverages quantum
annealing to eliminate redundant and noisy APs from the fingerprint. We
illustrated how to formulate the APs selection problem as a QUBO problem and
how to encode the number of powerful APs constraint in the QUBO objective
function. We deployed our quantum algorithm on a real D-Wave Systems quantum
machine and evaluated its performance in a real testbed for floor
localization. Our results showed that by selecting fewer than 14% of the
available APs, our quantum algorithm achieves the same floor localization
accuracy compared to utilizing the entire set of APs and a better accuracy
over utilizing the reduced dataset by the classical APs selection counterpart.
Moreover, the proposed quantum algorithm achieved a speedup of more than an
order of magnitude over the corresponding classical APs selection algorithms
which highlight the promise of the proposed quantum APs selection algorithm
for large-scale localization.

## References

  * [1] Moustafa Youssef.  Towards truly ubiquitous indoor localization on a worldwide scale.  In ACM SIGSPATIAL., pages 1–4, 2015.
  * [2] Ahmed Shokry, Moustafa Elhamshary, and Moustafa Youssef.  The tale of two localization technologies: Enabling accurate low-overhead wifi-based localization for low-end phones.  In ACM SIGSPATIAL, pages 1–10, 2017.
  * [3] Ahmed Shokry, Marwan Torki, and Moustafa Youssef.  Deeploc: a ubiquitous accurate and low-overhead outdoor cellular localization system.  In ACM SIGSPATIAL., pages 339–348, 2018.
  * [4] Mohamed Ibrahim and Moustafa Youssef.  Cellsense: An accurate energy-efficient gsm positioning system.  IEEE Transactions on Vehicular Technology, 61(1):286–296, 2011.
  * [5] Alan Jović, Karla Brkić, and Nikola Bogunović.  A review of feature selection methods with applications.  In 2015 38th international convention on information and communication technology, electronics and microelectronics (MIPRO), pages 1200–1205. Ieee, 2015.
  * [6] Ahmed Shokry and Moustafa Youssef.  Challenge: Quantum Computing for Location Determination.  arXiv e-prints, pages arXiv–2106, 2021.
  * [7] Ahmed Shokry and Moustafa Youssef.  Towards Quantum Computing for Location Tracking and Spatial Systems.  In Proceedings of the 29th International Conference on Advances in Geographic Information Systems, pages 278–281, 2021.
  * [8] Ahmed Shokry and Moustafa Youssef.  Device-independent Quantum Fingerprinting for Large Scale Localization.  In 2022 MedComNet, pages 208–215. IEEE, 2022.
  * [9] Ahmed Shokry and Moustafa Youssef.  Quantum fingerprinting for heterogeneous devices localization.  Computer Communications, 2023.
  * [10] Ahmed Shokry and Moustafa Youssef.  A Quantum Algorithm for RF-based Fingerprinting Localization Systems.  IEEE Conference on Local Computer Networks (LCN), 2022.
  * [11] Ahmed Shokry and Moustafa Youssef.  QLoc: A Realistic Quantum Fingerprint-based Algorithm for Large Scale Localization.  IEEE International Conference on Quantum Computing and Engineering (QCE), 2022.
  * [12] Ahmed Shokry and Moustafa Youssef.  Qradar: A deployable quantum euclidean similarity large-scale localization system.  In 2023 IEEE 48th Conference on Local Computer Networks (LCN), pages 1–8. IEEE, 2023.
  * [13] Yousef Zook, Ahmed Shokry, and Moustafa Youssef.  A quantum fingerprinting algorithm for next generation cellular positioning.  arXiv preprint arXiv:2306.08108, 2023.
  * [14] Maximilian Schlosshauer.  Decoherence, the measurement problem, and interpretations of quantum mechanics.  Reviews of Modern physics, 76(4):1267, 2005.
  * [15] Gary Kochenberger, Jin-Kao Hao, Fred Glover, Mark Lewis, Zhipeng Lü, Haibo Wang, and Yang Wang.  The unconstrained binary quadratic programming problem: a survey.  Journal of combinatorial optimization, 28:58–81, 2014.
  * [16] Fred Glover, Gary Kochenberger, and Yu Du.  A tutorial on formulating and using qubo models.  arXiv preprint arXiv:1811.11538, 2018.
  * [17] Sergio Boixo, Troels F Rønnow, Sergei V Isakov, Zhihui Wang, David Wecker, Daniel A Lidar, John M Martinis, and Matthias Troyer.  Evidence for quantum annealing with more than one hundred qubits.  Nature physics, 10(3):218–224, 2014.
  * [18] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * [19] Catherine C McGeoch.  Adiabatic quantum computation and quantum annealing: Theory and practice.  Springer Nature, 2022.
  * [20] Román Orús, Samuel Mugel, and Enrique Lizaso.  Quantum computing for finance: Overview and prospects.  Reviews in Physics, 4:100028, 2019.
  * [21] Rajdeep Kumar Nath, Himanshu Thapliyal, and Travis S Humble.  A review of machine learning classification using quantum annealing for real-world applications.  SN Computer science, 2:1–11, 2021.
  * [22] David J Griffiths and Darrell F Schroeter.  Introduction to quantum mechanics.  Cambridge university press, 2018.
  * [23] Fred Glover, Gary Kochenberger, and Yu Du.  A tutorial on formulating and using qubo models.  arXiv preprint arXiv:1811.11538, 2018.
  * [24] Litao Han, Li Jiang, Qiaoli Kong, Ji Wang, Aiguo Zhang, and Shiming Song.  Indoor Localization within Multi-Story Buildings Using MAC and RSSI Fingerprint Vectors.  Sensors, 19(11):2433, May 2019.
  * [25] K. Kim, S. Lee, and K. Huang.  A scalable deep neural network architecture for multi-building and multi-floor indoor localization based on Wi-Fi fingerprinting.  Big Data Analytics, 3(4), 2018.
  * [26] Khaled Alkiek, Aya Othman, Hamada Rizk, and Moustafa Youssef.  Deep Learning-Based Floor Prediction Using Cell Network Information.  In Proceedings of the 28th International Conference on Advances in Geographic Information Systems, SIGSPATIAL ’20, page 663–664, New York, NY, USA, 2020. Association for Computing Machinery.
  * [27] Giuseppe Caso, Luca De Nardis, Filip Lemic, Vlado Handziski, Adam Wolisz, and Maria-Gabriella Di Benedetto.  ViFi: Virtual Fingerprinting WiFi-Based Indoor Positioning via Multi-Wall Multi-Floor Propagation Model.  IEEE Transactions on Mobile Computing, 19(6):1478–1491, 2020.
  * [28] Harald Cramér.  Mathematical methods of statistics, volume 26.  Princeton university press, 1999.
  * [29] Philip Sedgwick.  Pearson’s correlation coefficient.  Bmj, 345, 2012.
  * [30] Brian Beavis and Ian Dobbs.  Optimisation and stability theory for economic analysis.  Cambridge university press, 1990.
  * [31] Sascha Mücke, Raoul Heese, Sabine Müller, Moritz Wolter, and Nico Piatkowski.  Feature selection on quantum computers.  Quantum Machine Intelligence, 5(1):11, 2023.
  * [32] D-Wave Systems.  Hybrid Solvers for Quadratic Optimization.  2022\.
  * [33] D-Wave Systems.  D-Wave Hybrid Solver Service + Advantage: Technology Update.  2022\.
  * [34] Joaquín Torres-Sospedra, Raúl Montoliu, Adolfo Martínez-Usó, Joan P. Avariento, Tomás J. Arnau, Mauri Benedito-Bordonau, and Joaquín Huerta.  UJIIndoorLoc: A new multi-building and multi-floor database for WLAN fingerprint-based indoor localization problems.  In 2014 International Conference on Indoor Positioning and Indoor Navigation (IPIN), pages 261–270, 2014.
  * [35] Tin Kam Ho.  Random decision forests.  In Proceedings of 3rd international conference on document analysis and recognition, volume 1, pages 278–282. IEEE, 1995.
  * [36] Peter JM Van Laarhoven, Emile HL Aarts, Peter JM van Laarhoven, and Emile HL Aarts.  Simulated annealing.  Springer, 1987.
  * [37] Siddharth Muthukrishnan, Tameem Albash, and Daniel A Lidar.  Tunneling and speedup in quantum optimization for permutation-symmetric problems.  Physical Review X, 6(3):031010, 2016.

##  Appendix A Proof of lemma 1

Recall that the both the redundancy $R(\bm{x})$ and the importance $I(\bm{x})$
are $\geq 0$ due to non-negativity of the absolute value of the Pearson
correlation and the Cramer’s V measurements, respectively.

When $\alpha=0$, the objective function in Eq 5, $Q(\cdot,0)$, tries to
minimize the redundancy, leading to either the zero vector $\bar{\bm{0}}$ as
the optimal solution or one of one-hot encoding vectors with the index of AP
$i$ is one and zero otherwise, $\forall i\in\\{1,2,3,\dots,n\\}$ is optimal.
This covers the cases of not selecting any AP ($k=0$) or selecting a single AP
($k=1$).

When $\alpha=1$, the objective function in Eq 5, becomes

$Q(\bm{x},1)=-I(\bm{x})=-\sum_{i\in\\{1,2,\dots,n\\}}I_{i}x_{i}$ (10)

which is minimized by a vector of ones $\bar{\bm{1}}$, covering the case of
selecting all APs ($k=n$).

Now, for all $k\in\\{0,\dots,n\\}$, consider the functions,

$Q^{*}_{\leq
k}(\alpha):=\min_{\bm{x}\in\\{0,1\\}^{n}}Q_{\alpha}(\bm{x})~{}\text{s.t.
}\norm{\bm{x}}_{1}\leq k\;.$ (11)

These functions are piece-wise linear and strictly decreasing in $\alpha$
[31], due to

$\frac{\partial Q_{\alpha}(\bm{x})}{\partial\alpha}=-(R(\bm{x})+I(\bm{x}))\leq
0$ (12)

This implies further for all $\alpha\in[0,1]$ and $k\in\\{1,2,3,\dots,n\\}$
that

$\displaystyle\min_{\begin{subarray}{c}\bm{x}\in\\{0,1\\}^{n}\\\
\norm{\bm{x}}_{1}\leq
k-1\end{subarray}}R(x)\leq\min_{\begin{subarray}{c}\bm{x}\in\\{0,1\\}^{n}\\\
\norm{\bm{x}}_{1}\leq k\end{subarray}}R(x)$ (13)
$\displaystyle\max_{\begin{subarray}{c}\bm{x}\in\\{0,1\\}^{n}\\\
\norm{\bm{x}}_{1}\leq
k-1\end{subarray}}I(x)\leq\max_{\begin{subarray}{c}\bm{x}\in\\{0,1\\}^{n}\\\
\norm{\bm{x}}_{1}\leq k\end{subarray}}I(x)\;.$ (14)

Hence, for any $k<k^{\prime}$, unless $Q^{*}_{\leq k}$ = $Q^{*}_{\leq
k^{\prime}}$,

$Q^{*}_{\leq k^{\prime}}(\alpha)\leq Q^{*}_{\leq k}(\alpha)$ (15)

as a consequence of $Q^{*}_{\leq k}$ and $Q^{*}_{\leq k^{\prime}}$ being non-
increasing, from which follows the proof. If indeed $Q^{*}_{\leq k}$ =
$Q^{*}_{\leq k^{\prime}}$, both binary vectors $\bm{x}$ and $\bm{x}^{\prime}$
with $\norm{\bm{x}}_{1}=k$ and $\norm{\bm{x}^{\prime}}_{1}=k^{\prime}$ would
be optimal, from which the proof still follows. ∎