# On-the-Fly Ensemble Pruning in Evolving Data Streams

Sanem Elbasi  Bilkent Information Retrieval GroupBilkent University
<EMAIL_ADDRESS>,  Alican Büyükçakır  Bilkent Information
Retrieval GroupBilkent University<EMAIL_ADDRESS>,  Hamed
Bonab  College of Information and Computer SciencesUniversity of Massachusetts
Amherst<EMAIL_ADDRESS>and  Fazli Can  Bilkent Information Retrieval
GroupBilkent University<EMAIL_ADDRESS>
###### Abstract.

Ensemble pruning is the process of selecting a subset of component classifiers
from an ensemble which performs at least as well as the original ensemble
while reducing storage and computational costs. Ensemble pruning in data
streams is a largely unexplored area of research. It requires analysis of
ensemble components as they are running on the stream, and differentiation of
useful classifiers from redundant ones. We present CCRP, an on-the-fly
ensemble pruning method for multi-class data stream classification empowered
by an imbalance-aware fusion of class-wise component rankings. CCRP aims that
the resulting pruned ensemble contains the best performing classifier for each
target class and hence, reduces the effects of class imbalance. The conducted
experiments on real-world and synthetic data streams demonstrate that
different types of ensembles that integrate CCRP as their pruning scheme
consistently yield on par or superior performance with 20% to 90% less average
memory consumption. Lastly, we validate the proposed pruning scheme by
comparing our approach against pruning schemes based on ensemble weights and
basic rank fusion methods.

Ensemble pruning, Ensemble efficiency, Multiple classifiers systems, Thinning,
Concept drift

##  1\. Introduction

Data streams are environments such as network event logs, video streams, call
records, and transaction records where vast amounts of data is generated at
high speed (Aggarwal, 2007). Classification models in data streams have to
work under strict time and memory constraints, and be able to adapt changes in
the distribution of data over time (Bifet et al., 2009). Such changes in the
relation between data points and their corresponding classes over time are
called concept drifts (Gama, 2010). Robustness against concept drifts is a
topic that is extensively studied in the literature (Bonab and Can, 2018;
Brzezinski and Stefanowski, 2014). As data continuously arrive, learned
information from the past instances become irrelevant under concept drifts.
Classification models need to adapt to the new concepts by not only learning
the new concepts but also forgetting the now-obsolete ones (Brzezinski and
Stefanowski, 2014).

Ensembles are common choices for models in stream environments, as they often
improve predictive performance while providing robustness against concept-
drifting and time-evolving data (Bifet et al., 2009). Ensembles employ
different approaches for dealing with concept drifts, such as replacing the
weakest / lowest weight classifier with a new one that is trained with the
more recent data (Bonab and Can, 2018), or updating ensemble weights regularly
(Wang et al., 2003; Bonab and Can, 2018). Despite the popularity and
prevalence of ensemble learning in data streams, ensemble pruning is still a
vastly undiscovered area of research in stream mining community.

Figure 1. Pipeline of Proposed Method. Prediction history of each component in
the ensemble are recorded. At each row of records, colored cells indicate
relevance scores for corresponding classes. Then, these records are used in
calculating class-wise losses for each component, and class-wise rankings of
those components. The acquired rankings are combined using a rank fusion
algorithm (Modified Borda Count). The resulting top $\varphi$ components
($\varphi<K$) are selected, and the rest is pruned.

In this study, we tackle the task of ensemble pruning in data streams for
multi-class classification. Our contributions are as follows. We propose,

  * •

an explicit pruning technique that can be integrated into any type of
streaming ensemble, and called whenever pruning is requested. Our method
results in pruned ensembles that are better performing in terms of accuracy
and more efficient in terms of memory consumption (as demonstrated in Section
5).

  * •

to the best of our knowledge, the first large-scale 111The previous studies
(Krawczyk, 2015; Cruz et al., 2019) experimented on datasets with number of
instances in the scales of hundreds or thousands on-the-fly ensemble pruning
method for streaming data.

  * •

a class imbalance-aware pruning method, where a pruned ensemble does not lose
its ability to classify rare or less-frequent classes, as a result of our
class-wise component prediction analysis and ranking (as described in Section
4).

##  2\. Problem Definition and Notation

We consider the problem of ensemble pruning within the context of supervised
classification in data streams. Data stream $\mathcal{D}$ is a (possibly
infinite) sequence of time-ordered data that consists of pairs
$(X^{(t)},y^{(t)})$ where $t$ denotes the arrival time. In multi-class
classification, the target $y^{(t)}$ has $L$-many classes where $L>2$. An
ensemble model $\xi$ is a group of $K$ component classifiers, i.e.
$\xi=\\{C_{1},C_{2},\dots,C_{K}\\}$ where the prediction of the model for an
instance is a combination of individual hypotheses ($h_{k}(X)$) of its
component classifiers where $k$ denotes the classifier index.

Let $\hat{\rho}_{k,l}$ be the record of the predictions of $k$th component for
the $l$th class on a sliding window containing the latest $N$ instances.
Similarly, let $\rho_{l}$ be the record of the ground truth information for
the $l$th class on the sliding window.

The aim of ensemble pruning is selecting a subset of classifiers
$\xi^{\prime}\subset\xi$ such that the ensemble’s both efficiency and
predictive performance are improved. Here, let us denote the size of the
pruned ensemble as $|\xi^{\prime}|=\varphi$ where $\varphi<K$.

##  3\. Related Work

Online Ensembles. The literature on the online ensembles is abundant (Krawczyk
et al., 2017). Here, we only mention the ensembles that are used in our
experiments. One of the most commonly used and well-known ensembles in data
streams is AWE (Accuracy Weighted Ensemble) (Wang et al., 2003). It is a
chunk-based ensemble where each classifier are weighted according to their
expected accuracy in the latest chunk. Similarly, GOOWE (Bonab and Can, 2018)
is another chunk-based stacking ensemble that assigns optimal weights to its
components using a geometric interpretation of target space and solving an
optimization problem as the chunks are processed. Both of these ensembles are
dynamic and they replace one of their components at the end of each chunk.

Pruning in Data Streams. A related topic to on-the-fly ensemble pruning is
Dynamic Ensemble Selection (DES) (Ko et al., 2008). In DES, a subset of
components from a pool is selected for each test instance. DES assigns
estimated levels of competence for each classifier in the pool, so that a
selection method selects the classifiers to be assigned to incoming instances.
Note that, unlike ensemble pruning, selection of a subset of classifiers for
every instance does not reduce the size of the pool of classifiers. Firefly
(Krawczyk, 2015) introduces a swarm intelligence approach to pruning with one-
class classifiers where input classifiers are members of a firefly population.
Interactions between fireflies describe the effectiveness of the classifiers,
and best representatives of diverse groups in the population are selected as
the pruning step. Lastly, FIRE-DES++ (Cruz et al., 2019) incorporates a pre-
selection stage where only the classifiers that can correctly classify at
least a pair of instances with different classes can proceed to the next stage
of classification. Again, the pre-selection stage in FIRE-DES++ does not
change the size of the pool of classifiers.

##  4\. Proposed Method

Algorithm 1 CCRP: Class-wise Component Ranking-based Pruner

1:$\mathcal{D}$: data stream, $\xi$: ensemble, $\varphi$: pruned ensemble size

2:$\xi^{\prime}$: pruned ensemble

3:Initialize $\rho$ as an empty FIFO buffer of size $N$.

4:Initialize $\hat{\rho}$ as an empty FIFO buffer of size $K\times N$.

5:for $(X,y)\in\mathcal{D}$ do

6: $\hat{\rho}_{k}.$ append$(h_{k}(X))\;$ for $\forall k$ $\triangleright$[
Preliminary Phase ]

7: $\rho.$ append$(y)$

8: if prune then $\triangleright$[ Phase I ]

9: for $l\in L$ do $\triangleright$Start CCRP

10: for $k\in K$ do

11: $scores.$append$($MSE$(\rho_{l},\hat{\rho}_{k,l}))$ $\triangleright$Class-
wise MSEs

12: end for

13: $ranks[l]=$ argsort$(scores)$ $\triangleright$Class-wise order of
components

14: end for

15: $ccrp\\_ranks=$ ModifiedBorda$(ranks)$ $\triangleright$[ Phase II ]

16: $\xi^{\prime}\leftarrow$ top ranked $\varphi$ components based on
$ccrp\\_ranks$

17: $\xi\leftarrow\xi^{\prime}$ $\triangleright$Pruned ensemble in effect

18: end if

19:end for

We introduce, CCRP, Class-wise Component Rankings based Pruner for multi-class
online ensembles. CCRP can be integrated to work with both dynamic and static
multi-class ensembles to improve the performance and memory consumption. CCRP
can be performed at any point in time on the stream. The proposed method
consists of three phases: Preliminary Phase: Recording Component Predictions
on the Sliding Window, Phase I: On-the-Fly Performance Analysis and Class-wise
Ranking of Components, and Phase II: Fusion of Rankings and Component
Selection (see Figure 1).

#### Preliminary Phase: Recording Component Predictions on the Sliding Window.

This phase only applies to the ensembles that are not already recording the
component predictions on the sliding window (e.g. OzaBagging (Oza, 2005)). In
this phase, for the latest $N$ data instances, predictions of classifiers and
the ground truth are recorded as $\rho$ and $\hat{\rho}$ respectively (Alg.1
line 6). Since only the latest $N$ instances are kept in the records, the
required memory for this process is fixed, despite being proportional to the
number of classes.

#### Phase I: On-the-Fly Performance Analysis and Class-wise Ranking of
Components.

This is the initial phase of CCRP in which per class performances of
classifiers are measured. In this phase, for each class $l$, predictions of
classifiers and the ground truth are extracted as $\hat{\rho}_{k,l}$ and
$\rho_{l}$ respectively (Alg.1 line 8). Then, scores of classifiers are
calculated using Mean Square Error (Eqn. 1) for each class (Alg.1 line 11).

(1)
$\mathcal{L}_{k,l}(X)=\sum^{N}_{i=1}\Big{(}\hat{\rho}_{k,l}^{(i)}-\rho_{l}^{(i)}\Big{)}^{2},\qquad
1\leq k\leq K,\;1\leq l\leq L$

#### Phase II: Fusion of Rankings and Component Selection.

In the final phase, CCRP generates the overall ranking of the classifiers,
using a modified version Borda Count (Nuray and Can, 2006) rank fusion method
(Alg.1 line 15). In Modified Borda Count (MBC, hereafter), CCRP assigns
$K\times L$ points (instead of $K$ points in the regular Borda Count (Nuray
and Can, 2006)) to the highest ranking components in the class-wise rankings.
This ensures the winning component for each class to appear at the top $L$
places in the overall ranking. Afterwards, the second classifier gets $K-1$
points and each proceeding classifier gets 1 point less than its successor.
Then, the top $\varphi$ classifiers from the overall ranking are selected as
the members of the pruned ensemble. It is recommended (Bonab and Can, 2019)
that the pruned ensemble size should be at least equal to the number of
classes. Taking this into account, CCRP guarantees that the best performing
classifier for each class is included in the pruned ensemble when $\varphi\geq
L$.

##  5\. Experimental Setup and Results

###  5.1. Setup

The proposed method and modifications to existing ensembles are integrated
into scikit-multiflow (Montiel et al., 2018) library. The experiments are
evaluated prequentially (Gama et al., 2009). We report prequential and overall
accuracy as performance metrics. In addition, for memory efficiency analysis
in the experiments, we define Mean Memory Consumption Ratio ($\mu$) which
indicates the percentage of average memory the pruned ensemble uses with
respect to the original ensemble (Eqn. 2).

(2)
$\mu=\frac{\sum_{\forall\text{chunk}}\;Size(\xi^{\prime})}{\sum_{\forall\text{chunk}}\;Size(\xi)}$

Throughout the experiments, we use two typical dynamic ensembles (AWE (Wang et
al., 2003) and GOOWE (Bonab and Can, 2018)) with Hoeffding Trees (Domingos and
Hulten, 2000) as their base classifiers. CCRP is performed whenever the
ensembles are full, i.e. when the ensemble size reaches $K$. After pruning,
the ensembles continue to grow until they reach the maximum size $K$, where
pruning takes place again. This process continues indefinitely.

Table 1. Datasets | Name | # Features | # Classes | # Instances
---|---|---|---|---
Real | COVTYPE (Bifet et al., 2013) | 54 | 7 | 581.012
Poker Hand (Bifet et al., 2013) | 10 | 10 | 829.201
Rialto (Losing et al., 2016) | 27 | 10 | 82.250
Synth | Mov. Squares (Losing et al., 2016) | 2 | 4 | 200.000
Mov. RBF (Losing et al., 2016) | 10 | 5 | 200.000
Trn. Chssbrd (Losing et al., 2016) | 2 | 8 | 200.000

The models are ran on three real-world and three synthetic well-known datasets
with concept drifts (Bifet et al., 2013; Losing et al., 2016). A summary of
the dataset information is provided in Table 1.

Figure 2. The impact of CCRP ($\varphi=L$) on the prequential accuracy of two
ensemble models on three datasets over the first few thousand instances. The
first occurrence of pruning is denoted with a vertical red line.

###  5.2. Effect of CCRP on Predictive Performance and Memory Efficiency

Does CCRP yield more effective and efficient ensemble models? We respond to
this question by taking both accuracy and memory consumption ($\mu$) into
account. We perform experiments with AWE (where $K=20$) and GOOWE ($K=30$) to
investigate the effect of CCRP on the ensembles.

To examine the change in the behavior of ensembles, ceteris paribus, we
consider the case when CCRP is performed for the first time on the stream. In
Figure 2, it can be observed that the overall accuracy of the pruned ensemble
is always higher than the original ensemble after pruning with CCRP. In
addition to the improvement in accuracy, it should be noted that the ensembles
after pruning are relatively more robust to the changes in the data stream.
This can be observed from the steepness of the declining accuracies of the
original ensembles compared to the pruned ensembles.

Table 2. Effect of CCRP on Accuracy and Mean Memory Consumption Ratio ($\mu$) for Different Ensembles |  AWE (Wang et al., 2003) |  GOOWE (Bonab and Can, 2018)
---|---|---
| Default | CCRP | $\mu$ | Default | CCRP | $\mu$
COVTYPE |  0.661 | 0.659 | 75% |  0.852 | 0.829 | 10%
PokerHand | 0.525 |  0.585 | 83% |  0.712 | 0.651 | 14%
Rialto | 0.377 |  0.379 | 80% |  0.428 | 0.420 | 55%
Mov. Squares | 0.478 |  0.532 | 71% | 0.354 |  0.359 | 52%
Mov. RBF |  0.499 |  0.499 | 72% |  0.433 | 0.396 | 19%
Trn. Chssbrd |  0.146 | 0.141 | 71% | 0.685 |  0.791 | 70%

Figure 2 indicates that the initial calls of CCRP locally improves the
ensemble performance. Yet, is that the case throughout the stream? To inspect
the effect of CCRP on both accuracy and memory consumption throughout the
stream, we experiment with AWE ($K=20$) and GOOWE ($K=30$) on all of the
datasets, and report overall accuracies, as well as $\mu$ (Table 2).

In Table 2, notice that CCRP yields $10\%\leq\mu\leq 83\%$, which means the
pruned ensembles consume at most $90\%$ and at least $17\%$ less memory over
the course of their run. In 5 out of 12 cases, CCRP improves the performance
of ensembles while reducing the memory consumption by $20$ to $25\%$. In a
small number of cases (3 out of 12: GOOWE on COVTYPE, PokerHand and
MovingRBF), ensembles without pruning outperforms the ones with pruning. These
are also the cases where CCRP yields the most drastic reductions in memory
consumption, around $80$ to $90\%$. This behaviour may occur due to removing
old and strong components of the ensemble in the presence of drifts. For the
remaining 4 cases, the overall accuracy for the pruned models are on par with
the ones that are not pruned. A decrease in memory consumption around $30\%$
can be observed in these cases, as well. These memory savings imply similar
savings in CPU time, since memory consumption of an ensemble is proportional
to the complexity of its components.

###  5.3. Superiority of CCRP over Different Baseline Selection Techniques

Is using CCRP’s pruning scheme indeed better than selecting the highest
weighted $\varphi$ components from a weighted ensemble? What difference does
using MBC make with respect to the regular Borda Count? To investigate these
questions, we design the following experiment: We choose a weighted ensemble
(GOOWE (Bonab and Can, 2018)) which normally employs a replacement policy
where the lowest weight component is replaced at the end of each chunk. We,
then, modify GOOWE’s component replacement policy so that it removes its
component using CCRP, and Borda Count with $\varphi=K-1$, and places the newly
trained component to the freed-up spot. Overall accuracy values for several
datasets are reported in Table 3. The winning method(s) for each dataset is
underlined.

Table 3. Accuracy of GOOWE (with $K=10$) using Different Pruning Schemes (with $\varphi=9$) for the Selected Datasets | CCRP | Weight-based | Regular Borda
---|---|---|---
COVTYPE |  0.839 | 0.796 | 0.703
Poker |  0.726 | 0.667 | 0.701
Mov. RBF |  0.543 | 0.537 |  0.543
Trn. Chessboard |  0.452 | 0.450 | 0.449

This experiment shows that the internal component weights of an ensemble do
not necessarily indicate their contribution to the overall performance of the
ensemble. In addition, CCRP’s superiority over using Regular Borda (Nuray and
Can, 2006) shows MBC’s capability of handling class imbalance. By assigning
$K\times L$ points to the best performing classifier for a less frequent
class, MBC ensures that that classifier is not ranked poorly and selected for
the resulting subset of classifiers.

##  6\. Conclusion and Future Work

In this work, we present an on-the-fly ensemble pruner for evolving data
streams for the first time in the literature. Our approach utilizes
predictions of ensembles to generate class-wise rankings and combine these
rankings with a rank fusion algorithm. We demonstrate that the proposed
pruning scheme yields smaller and more efficient ensembles with on par or
better predictive performance consistently. Additionally, we show that the
proposed rank fusion algorithm works better than its alternatives. Our
experiments on different ensembles provide evidence that our method can be
integrated into any streaming ensemble.

For future work, we plan to do a rigorous analysis of diversity of the
ensemble components before and after pruning, considering various diversity
measures (Krawczyk et al., 2017). The proposed approach can be extended by
incorporating diversity of components into the rank fusion process.
Additionally, ”how to select the optimal prune size ($\varphi$)” is still an
open question. The effect of ensemble parameters and dataset properties on
$\varphi$ can be examined. Lastly, we observed that the time of the pruning is
crucial to obtain high predictive performance. Therefore, integrating a
concept drift detector into the proposed method and pruning the ensemble
whenever a drift is detected may further improve the performance of models.

## References

  * (1)
  * Aggarwal (2007) Charu C Aggarwal. 2007\.  _Data streams: models and algorithms_. Vol. 31.  Springer Science & Business Media.
  * Bifet et al. (2009) Albert Bifet, Geoff Holmes, Bernhard Pfahringer, Richard Kirkby, and Ricard Gavaldà. 2009\.  New ensemble methods for evolving data streams. In _ACM SIGKDD_. ACM, 139–148.
  * Bifet et al. (2013) Albert Bifet, Bernhard Pfahringer, Jesse Read, and Geoff Holmes. 2013.  Efficient data stream classification via probabilistic adaptive windows. In _Proceedings of the 28th Annual ACM Symposium on Applied Computing_. ACM, 801–806.
  * Bonab and Can (2018) Hamed Bonab and Fazli Can. 2018.  GOOWE: Geometrically Optimum and Online-Weighted Ensemble classifier for evolving data streams.  _ACM TKDD_ 12, 2 (2018), 25.
  * Bonab and Can (2019) Hamed Bonab and Fazli Can. 2019.  Less is more: a comprehensive framework for the number of components of ensemble classifiers.  _IEEE TNNLS_ (2019).
  * Brzezinski and Stefanowski (2014) Dariusz Brzezinski and Jerzy Stefanowski. 2014.  Reacting to different types of concept drift: The accuracy updated ensemble algorithm.  _IEEE TNNLS_ 25, 1 (2014), 81–94.
  * Cruz et al. (2019) Rafael MO Cruz, Dayvid VR Oliveira, George DC Cavalcanti, and Robert Sabourin. 2019.  FIRE-DES++: Enhanced online pruning of base classifiers for dynamic ensemble selection.  _Pattern Recognition_ 85 (2019), 149–160.
  * Domingos and Hulten (2000) Pedro Domingos and Geoff Hulten. 2000.  Mining High-speed Data Streams. In _ACM SIGKDD_. ACM, 71–80.
  * Gama (2010) Joao Gama. 2010\.  _Knowledge discovery from data streams_.  Chapman and Hall/CRC.
  * Gama et al. (2009) João Gama, Raquel Sebastião, and Pedro Pereira Rodrigues. 2009\.  Issues in evaluation of stream learning algorithms. In _ACM SIGKDD_. ACM, 329–338.
  * Ko et al. (2008) Albert HR Ko, Robert Sabourin, and Alceu Souza Britto Jr. 2008\.  From dynamic classifier selection to dynamic ensemble selection.  _Pattern Recognition_ 41, 5 (2008), 1718–1731.
  * Krawczyk (2015) Bartosz Krawczyk. 2015\.  One-class classifier ensemble pruning and weighting with firefly algorithm.  _Neurocomputing_ 150 (2015), 490–500.
  * Krawczyk et al. (2017) Bartosz Krawczyk, Leandro L Minku, João Gama, Jerzy Stefanowski, and Michał Woźniak. 2017.  Ensemble learning for data stream analysis: A survey.  _Information Fusion_ 37 (2017), 132–156.
  * Losing et al. (2016) Viktor Losing, Barbara Hammer, and Heiko Wersing. 2016\.  KNN classifier with self adjusting memory for heterogeneous concept drift. In _IEEE ICDM_. 291–300.
  * Montiel et al. (2018) Jacob Montiel, Jesse Read, Albert Bifet, and Talel Abdessalem. 2018.  Scikit-multiflow: a multi-output streaming framework.  _The Journal of Machine Learning Research_ 19, 1 (2018), 2915–2914.
  * Nuray and Can (2006) Rabia Nuray and Fazli Can. 2006.  Automatic ranking of information retrieval systems using data fusion.  _IPM_ 42, 3 (2006), 595–614.
  * Oza (2005) Nikunj C Oza. 2005\.  Online Bagging and Boosting. In _2005 IEEE International Conference on Systems, Man and Cybernetics_ , Vol. 3. IEEE, 2340–2345.
  * Wang et al. (2003) Haixun Wang, Wei Fan, Philip S Yu, and Jiawei Han. 2003\.  Mining Concept-Drifting Data Streams using Ensemble Classifiers. In _ACM SIGKDD_. ACM, 226–235.


# Beyond Diagonal RIS: Key to Next-Generation Integrated Sensing and
Communications?

Tara Esmaeilbeig, Kumar Vijay Mishra and Mojtaba Soltanalian Tara Esmaeilbeig
and Mojtaba Soltanalian are with the ECE Departement, University of Illinois
at Chicago, Chicago, IL 60607 USA. Email: {zesmae2<EMAIL_ADDRESS>Kumar Vijay
Mishra is with the United States DEVCOM Army Research Laboratory, Adelphi, MD
20783 USA. E-mail<EMAIL_ADDRESS>work was sponsored by the Army Research
Office, accomplished under Grant Number W911NF-22-1-0263. The views and
conclusions contained in this document are those of the authors and should not
be interpreted as representing the official policies, either expressed or
implied, of the Army Research Office or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Government
purposes notwithstanding any copyright notation herein.

###### Abstract

Reconfigurable intelligent surface (RIS) have introduced unprecedented
flexibility and adaptability toward smart wireless channels. Recent research
on integrated sensing and communication (ISAC) systems has demonstrated that
RIS platforms enable enhanced signal quality, coverage, and link capacity. In
this paper, we explore the application of fully-connected beyond diagonal RIS
(BD-RIS) to ISAC systems. BD-RIS introduces additional degrees of freedom by
allowing non-zero off-diagonal elements for the scattering matrix, potentially
enabling further functionalities and performance enhancements. In particular,
we consider the joint design objective of maximizing the weighted sum of the
signal-to-noise ratio (SNR) at the radar receiver and communication users by
leveraging the extra degrees-of-freedom offered in the BD-RIS setting. These
degrees-of-freedom are unleashed by formulating an alternating optimization
process over known and auxiliary (latent) variables of such systems. Our
numerical results reveal the advantages of deploying BD-RIS in the context of
ISAC and the effectiveness of the proposed algorithm by improving the SNR
values for both radar and communication users by several orders of magnitude.

###### Index Terms:

Fully-connected RIS, ISAC, non-convex optimization, NLoS sensing, smart radio
environment.

##  I Introduction

Reconfigurable intelligent surfaces (RIS) have demonstrated great potential in
realizing smart programmable environments by modifying wireless channels [1,
2] and establishing reliable links for communication and sensing [3]. In
particular, RIS have been designed to enhance energy efficiency [4], physical
security [5, 6], radar target estimation and detection [7], and facilitate
integrated sensing and communications (ISAC) [8, 9]. In this paper, we focus
on RIS-enabled ISAC.

ISAC improves spectrum access and usage while reducing both hardware, power,
and signaling costs [10, 11, 12]. The literature indicates that RIS benefits
ISAC in both line-of-sight (LoS) and non-LoS (NLoS) scenarios [13]. When LoS
links exist between base station (BS) to communications users or radar
targets, RIS is used to compensate for the propagation loss [14, 15]. When the
LoS to either users or targets is blocked, RIS is utilized to establish a
virtual LoS or NLoS link to bypass the obstructions [16, 17].

All of the above-mentioned works employ conventional RIS architecture
characterized by a diagonal scattering matrix (D-RIS), also known as phase
shift matrix. In D-RIS each element is controlled by a tunable impedance
connected to ground. Whereas in [18], beyond diagonal RIS (BD-RIS) is
introduced by connecting the elements of RIS to each other through a
reconfigurable impedance network. When all or a group of the elements are
connected to each other the BD-RIS is called fully-connected or group-
connected BD-RIS, respectively. The maximum degrees-of-freedom offered by the
fully-connected BD-RIS has shown great potential in meeting the upper bound of
received signal power in [19]. BD-RIS extends the concept of conventional RIS
by actively manipulating the cross-polarization components, inter-element
coupling, and coherent scattering effects [20, 21]. By dynamically adjusting
these off-diagonal elements, beyond-diagonal RIS reshapes the wavefront in
ways previously unattainable, empowering ISAC systems with enhanced
resolution, interference suppression, and target discrimination capabilities
[22, 23, 24]. In this paper, we explore the benefit of fully-connected BD-RIS
in ISAC applications. In general, ISAC systems may follow communications-
centric design, radar-centric or integrated design [25]. Prior works,
including [23], have focused on radar-centric design of BD-RIS for ISAC by
maximizing the target detection probability subject to a minimum quality of
service at communications users.

Contrary to these studies, we employ a joint communications and sensing design
that offers a tunable trade-off between communications and sensing
performance. In particular, our design criterion is a weighted sum of the
signal-to-noise ratio (SNR) at the communications and radar (sensing)
receivers. This approach of weighted joint criterion has previously been
considered in ISAC system design studies [26, 27]. We show that optimization
of this objective with respect to the hardware constraints of fully-connected
BD-RIS is nonconvex. In order to tackle the resulting complicated non-convex
optimization problem, we develop an efficient alternating optimization
algorithm following a penalty-based framework for unimodular quartic
optimization over a finite alphabet set. Our numerical simulations, study the
performance of BD-RIS with different number of resolution bits. Additionally,
we illustrate that the BD-RIS outperforms the D-RIS with the same number of
reflecting elements in terms of enhancing the communications SNR, radar SNR,
and, hence, their weighted sum.

Throughout this paper, we use bold lowercase and bold uppercase letters for
vectors and matrices, respectively. The identity matrix of size $s\times s$ is
denoted by $\mathbf{I}_{s}$. The sets of complex and real numbers are
$\mathbb{C}$ and $\mathbb{R}$, respectively; $(\cdot)^{\top}$ and
$(\cdot)^{\mathrm{H}}$ are the vector/matrix transpose and the Hermitian
transpose, respectively; trace of a matrix is $\operatorname{Tr}(.)$. The
Kronecker product is denoted by $\otimes$. The vectorized form of a matrix
$\mathbf{B}$ is written as $\mathrm{vec}\left(\mathbf{B}\right)$. The maximum
eigenvalue of $\mathbf{B}$ is denoted by $\lambda_{max}(\mathbf{B})$. The
angle/phase and real part of a complex number is denoted by
$\operatorname{arg}\left(\cdot\right)$ and
$\operatorname{Re}\left(\cdot\right)$, respectively. The operation
$\mathrm{vec}_{{}_{K,L}}^{-1}\left(\mathbf{c}\right)$ reshapes the input
vector $\mathbf{c}\in\mathbb{C}^{KL}$ into a matrix
$\mathbf{C}\in\mathbb{C}^{K\times L}$ such that
$\mathrm{vec}\left(\mathbf{C}\right)=c$. For a symmetric $n\times n$ matrix
$\mathbf{A}$, the half-vectorization, denoted by
$\operatorname{vech}(\mathbf{A})$, is the $(n(n+1)/2)\times 1$ column vector
obtained by vectorizing only the lower triangular part of $\mathbf{A}$ i.e.
$\operatorname{vech}(\mathbf{A})=[\mathbf{A}_{1,1},\ldots,\mathbf{A}_{n,1},\mathbf{A}_{2,2},\ldots,\mathbf{A}_{n,2},\ldots,\mathbf{A}_{n-1,n-1},\mathbf{A}_{n,n-1},$
$\mathbf{A}_{n,n}]^{\mathrm{T}}$. Moreover, the unique _duplication_ matrix
$\mathbf{D}_{n}$ transforms $\operatorname{vech}(\mathbf{A})$ into
$\operatorname{vec}(\mathbf{A})$ as
$\operatorname{vec}(\mathbf{A})=\mathbf{D}_{n}\operatorname{vech}(\mathbf{A})$.

##  II System Model

Consider a BD-RIS-enabled ISAC system (Fig. 1) comprising a dual-function base
station (DFBS) with $N$ elements each in transmit (Tx) and receive (Rx)
antenna arrays. The BD-RIS comprises $L=L_{x}\times L_{y}$ reflecting elements
arranged as a uniform planar array (UPA) with $L_{x}$ ($L_{y}$) elements along
the x- (y-) axes in the Cartesian coordinate plane. Define the RIS steering
vector
$\mathbf{a}(\theta_{h},\theta_{v})=\mathbf{a}_{x}(\theta_{h},\theta_{v})\otimes\mathbf{a}_{y}(\theta_{h},\theta_{v})$,
where $\theta_{h}$ ($\theta_{v}$) is the azimuth (elevation) angle,
$\mathbf{a}_{x}(\theta_{h},\theta_{v})=[1,e^{\textrm{j}\frac{2\pi
d}{\lambda}cos\theta_{h}sin\theta_{v}},\ldots,e^{\textrm{j}\frac{2\pi
d(L_{x}-1)}{\lambda}cos\theta_{h}sin\theta_{v}}]^{\top}$ and
$\mathbf{a}_{y}(\theta_{h},\theta_{v})=[1,e^{\textrm{j}\frac{2\pi
d}{\lambda}cos\theta_{h}sin\theta_{v}},\ldots,e^{\textrm{j}\frac{2\pi
d(L_{y}-1)}{\lambda}cos\theta_{h}sin\theta_{v}}]^{\top}$, $\lambda=c/f_{c}$ is
the carrier wavelength, $c=3\times 10^{8}$ m/s is the speed of light, $f_{c}$
is the carrier frequency, and $d=0.5\lambda$ is the inter-element spacing. The
fully-connected BD-RIS operation is characterized by the scattering or phase
matrix $\boldsymbol{\Phi}$ which is a complex symmetric unitary, i.e.,

$\displaystyle\boldsymbol{\Phi}=\boldsymbol{\Phi}^{\top}$
$\displaystyle\text{and}\;\;\boldsymbol{\Phi}^{\mathrm{H}}\boldsymbol{\Phi}=\mathbf{I}.$
(1)

The phases may be assumed continuous or discrete/quantized. In general,
discrete phase settings simplify the hardware implementation and
control/programming of the RIS. Since discrete phase states are achieved using
digital switches or phase shifters, power consumption is lower than
continuously variable analog counterparts. For discrete-valued phases with
resolution $M$, we have
$[\boldsymbol{\Phi}]_{{}_{kl}}\in\boldsymbol{\Omega}_{{}_{M}}$, where

$\displaystyle\boldsymbol{\Omega}_{{}_{M}}$
$\displaystyle=\left\\{s=e^{\textrm{j}\omega},\omega\in\Pi_{{}_{M}}\right\\},$
$\displaystyle\Pi_{{}_{M}}$
$\displaystyle=\left\\{1,\frac{2\pi}{M},\cdots,\frac{2\pi(M-1)}{M}\right\\}.$
(2)

The continuous BD-RIS is equivalent to the discrete RIS with infinite
resolution i.e. $M\rightarrow\infty$ [28]. In other words, we have
$[\boldsymbol{\Phi}]_{{}_{kl}}\in\boldsymbol{\Omega}$ where
$\boldsymbol{\Omega}=\left\\{s=e^{\textrm{j}\omega},\omega\in[0,2\pi)\right\\}$.
In this paper, to achieve a discrete-valued BD-RIS phase design, we obtain the
optimal phases on the set $\boldsymbol{\Omega}$ and then project the solution
onto the set $\boldsymbol{\Omega}_{M}$.

DFBSCommunications usersTargetBD-RISBlocked LoS

Figure 1: A illustration of BD-RIS-enabled ISAC system. When the LoS is
blocked, the NLoS paths via the BD-RIS allow for establishing the link between
the targets/users with the DFBS.

We now describe the received communications signal at the users and the
sensing signal back-scattered from the target and received at the DFBS.
Communications Rx signal: Denote the direct or LoS channel state information
(CSI) and RIS-reflected non-line-of-sight (NLoS) CSI matrices by
$\mathbf{F}\in\mathds{C}^{K\times N}$ and $\mathbf{H}\in\mathds{C}^{K\times
L}$, respectively. Then, at each communications receiver after sampling, the
discrete-time received signal is $\mathbf{y}_{c,k}$. Concatenating the
received signals from all users, we obtain the $K\times 1$ vector:

$\mathbf{y}_{c}=(\mathbf{F}+\mathbf{H}\boldsymbol{\Phi}\mathbf{G})\mathbf{P}\mathbf{s}+\mathbf{n}_{c},$
(3)

where $\mathbf{P}\in\mathds{C}^{N\times K}$ is the DFBS precoder, the Tx-RIS
CSI matrix $\mathbf{G}\in\mathds{C}^{L\times N}$ is assumed to be estimated a
priori through suitable channel estimation techniques and
$\mathbf{n}_{c}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}_{c}\mathbf{I}_{K})$
is the spatially and temporally white noise at communications receivers.
Sensing Rx signal at DFBS: Consider a Swerling-0 [29] radar target located at
a NLoS location with respect to the DFBS. The direction-of-arrival (DoA) is
$(\theta_{h_{t}},\theta_{v_{t}})$ with respect to the RIS and radar cross-
section (RCS) $\alpha_{{}_{T}}$. Define
$\mathbf{R}=\alpha_{{}_{T}}\mathbf{G}^{T}\boldsymbol{\Phi}\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})^{\top}\boldsymbol{\Phi}\mathbf{G}$.

Stacking the echoes for all receiver antennas, the $N\times 1$ signal at radar
receiver after down conversion and sampling is

$\mathbf{y}_{{}_{r}}=\mathbf{R}\mathbf{P}\mathbf{s}+\mathbf{n}_{r},$ (4)

where
$\mathbf{n}_{r}\sim\mathcal{N}(\boldsymbol{0},\sigma^{2}_{r}\mathbf{I}_{N})$
is the noise at radar receiver.

Denote $\mathbf{C}=\mathbf{F}+\mathbf{H}\boldsymbol{\Phi}\mathbf{G}$. The
output SNR at communications and DFBS Rx are, respectively,

$\displaystyle\textrm{SNR}_{c}$
$\displaystyle=\frac{1}{\sigma^{2}_{c}}\mathrm{Tr}\left(\mathbf{C}\mathbf{P}\mathbf{P}^{\mathrm{H}}\mathbf{C}^{\mathrm{H}}\right),\textrm{and
}$ (5) $\displaystyle\textrm{SNR}_{r}$
$\displaystyle=\frac{1}{\sigma^{2}_{r}}\mathrm{Tr}\left(\mathbf{R}\mathbf{P}\mathbf{P}^{\mathrm{H}}\mathbf{R}^{\mathrm{H}}\right).$
(6)

In the next section, we use (5)-(6) to formalize the design problem for BD-
RIS-enabled ISAC system.

##  III BD-RIS design for ISAC

As the design criterion under BD-RIS hardware constraints (1) and (II), we use
the weighted sum of the radar and communications SNRs, with a weight factor
$\beta$ as
$\textrm{SNR}_{{}_{T}}=\beta\textrm{SNR}_{r}+(1-\beta)\textrm{SNR}_{c}$. Our
design problem is

$\displaystyle\mathcal{P}_{1}:$
$\displaystyle\;\underset{\boldsymbol{\Phi}\in\boldsymbol{\Omega}_{M}^{{}^{L\times
L}}}{\textrm{maximize}}\;\;\textrm{SNR}_{{}_{T}}$
$\displaystyle\;\textrm{subject
to}\;\;\boldsymbol{\Phi}^{\mathrm{H}}\boldsymbol{\Phi}=\mathbf{I}\;\;\text{and}\;\;\;\boldsymbol{\Phi}=\boldsymbol{\Phi}^{\top}.$
(7)

Since $\boldsymbol{\Phi}$ is symmetric, we solve $\mathcal{P}_{1}$ w.r.t
$\boldsymbol{\phi}=\operatorname{vech}(\boldsymbol{\Phi})\in\boldsymbol{\Omega}_{M}^{{}^{L(L+1)/2}}$
by means of a penalty-based optimization method. Following Lemmata 1 and 2,
state that the SNR at communications and the radar receivers are quadratic and
quartic with respect to $\boldsymbol{\phi}$, respectively.

######  Lemma 1.

For a symmetric scattering matrix $\boldsymbol{\Phi}$,

$\displaystyle\text{SNR}_{c}$
$\displaystyle=\boldsymbol{\phi}^{\mathrm{H}}\widetilde{\mathbf{Q}}\boldsymbol{\phi}+2\operatorname{Re}\left(\mathbf{q}^{\mathrm{H}}\boldsymbol{\phi}\right)+\mathbf{f}^{H}\widehat{\mathbf{P}}\mathbf{f},$
(8)

where
$\widetilde{\mathbf{Q}}=\mathbf{D}_{{}_{L}}^{\top}(\mathbf{G}^{H}\otimes\mathbf{H})^{\mathrm{H}}\widehat{\mathbf{P}}(\mathbf{G}^{H}\otimes\mathbf{H})\mathbf{D}_{{}_{L}}$,
with $\mathbf{D}_{{}_{L}}$ denoting the duplication matrix for
$\boldsymbol{\Phi}$, and
$\widehat{\mathbf{P}}=\frac{1}{\sigma^{2}_{c}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{K}\right)^{\mathrm{H}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{K}\right)$.
Also,
$\mathbf{q}=\mathbf{D}_{{}_{L}}^{\top}(\mathbf{G}^{H}\otimes\mathbf{H})^{\mathrm{H}}\widehat{\mathbf{P}}\mathbf{f}$.

###### Proof:

From (5) we have

$\displaystyle\textrm{SNR}_{c}$
$\displaystyle=\frac{1}{\sigma^{2}_{c}}\mathrm{Tr}\left(\mathbf{C}\mathbf{P}\mathbf{P}^{\mathrm{H}}\mathbf{C}^{\mathrm{H}}\right)=\frac{1}{\sigma^{2}_{c}}\vec{\mathbf{C}\mathbf{P}}^{\mathrm{H}}\vec{\mathbf{C}\mathbf{P}}$
$\displaystyle=\frac{1}{\sigma^{2}_{c}}\vec{\mathbf{C}}^{\mathrm{H}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{K}\right)^{\mathrm{H}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{K}\right)\vec{\mathbf{C}}$
$\displaystyle=\vec{\mathbf{C}}^{\mathrm{H}}\widehat{\mathbf{P}}\vec{\mathbf{C}}$
$\displaystyle=\mathbf{f}^{H}\mathbf{Q}\mathbf{f}+\mathbf{f}^{H}\widehat{\mathbf{P}}(\mathbf{G}^{H}\otimes\mathbf{H})\boldsymbol{\phi}+\boldsymbol{\phi}^{\mathrm{H}}(\mathbf{G}^{H}\otimes\mathbf{H})^{\mathrm{H}}\widehat{\mathbf{P}}\mathbf{f}$
$\displaystyle+\boldsymbol{\phi}^{\mathrm{H}}(\mathbf{G}^{H}\otimes\mathbf{H})^{\mathrm{H}}\widehat{\mathbf{P}}(\mathbf{G}^{H}\otimes\mathbf{H})\boldsymbol{\phi},$
(9)

where the last equality used
$\vec{\mathbf{C}}=\vec{\mathbf{F}}+(\mathbf{G}^{H}\otimes\mathbf{H})\vec{\boldsymbol{\Phi}}=\mathbf{f}+(\mathbf{G}^{H}\otimes\mathbf{H})\mathbf{D}_{{}_{L}}\boldsymbol{\phi}$.
Substituting $\widetilde{\mathbf{Q}}$ and $\mathbf{q}$ in (III) yields (8). ∎

######  Lemma 2.

For a symmetric scattering matrix $\boldsymbol{\Phi}$,

$\displaystyle\text{SNR}_{r}$
$\displaystyle=\boldsymbol{\phi}^{\mathrm{H}}\bar{\boldsymbol{Q}}(\boldsymbol{\Phi})\boldsymbol{\phi},$
(10)

where

$\displaystyle\bar{\boldsymbol{Q}}(\boldsymbol{\Phi})$
$\displaystyle=\mathbf{D}_{{}_{L}}^{\top}\left(\mathbf{G}^{\top}\otimes\mathbf{G}^{\top}\boldsymbol{\Phi}\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})\mathbf{a}^{\top}(\theta_{h_{t}},\theta_{v_{t}})\right)^{\mathrm{H}}\bar{\mathbf{P}}$
$\displaystyle\left(\mathbf{G}^{\top}\otimes\mathbf{G}^{\top}\boldsymbol{\Phi}\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})\mathbf{a}^{\top}(\theta_{h_{t}},\theta_{v_{t}})\right)\mathbf{D}_{{}_{L}},$
(11) $\displaystyle\bar{\mathbf{P}}$
$\displaystyle=\frac{|\alpha_{T}|^{2}}{\sigma^{2}_{r}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{N}\right)^{\mathrm{H}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{N}\right).$
(12)

###### Proof:

From (6) we have

$\displaystyle\textrm{SNR}_{r}$
$\displaystyle=\frac{1}{\sigma^{2}_{r}}\mathrm{Tr}\left(\mathbf{R}\mathbf{P}\mathbf{P}^{\mathrm{H}}\mathbf{R}^{\mathrm{H}}\right)=\frac{1}{\sigma^{2}_{r}}\vec{\mathbf{R}\mathbf{P}}^{\mathrm{H}}\vec{\mathbf{R}\mathbf{P}}$
$\displaystyle=\frac{1}{\sigma^{2}_{r}}\vec{\mathbf{R}}^{\mathrm{H}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{N}\right)^{\mathrm{H}}\left(\mathbf{P}^{\top}\otimes\mathbf{I}_{N}\right)\vec{\mathbf{R}}$
$\displaystyle=\boldsymbol{\phi}^{\mathrm{H}}\bar{\boldsymbol{Q}}(\boldsymbol{\Phi})\boldsymbol{\phi},$
(13)

where the last equality used
$\vec{\mathbf{R}}=\alpha_{T}\vec{\mathbf{G}^{\top}\boldsymbol{\Phi}\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})\mathbf{a}^{\top}(\theta_{h_{t}},\theta_{v_{t}})\boldsymbol{\Phi}\mathbf{G}}=\alpha_{T}\left(\mathbf{G}^{\top}\otimes\mathbf{G}^{\top}\boldsymbol{\Phi}\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})\mathbf{a}^{\top}(\theta_{h_{t}},\theta_{v_{t}})\right)\mathbf{D}_{{}_{L}}\boldsymbol{\phi}$.
∎

From Lemma 1-2, we conclude that the objective function in $\mathcal{P}_{1}$
is quartic w.r.t $\boldsymbol{\Phi}$. Define
$f(\boldsymbol{\phi})=\textrm{SNR}_{{}_{T}}=\beta\boldsymbol{\phi}^{\mathrm{H}}\bar{\boldsymbol{Q}}(\boldsymbol{\phi})\boldsymbol{\phi}+(1-\beta)\boldsymbol{\phi}^{\mathrm{H}}\widetilde{\mathbf{Q}}\boldsymbol{\phi}+2(1-\beta)\operatorname{Re}\left(\mathbf{q}^{\mathrm{H}}\boldsymbol{\phi}\right)$.
In order to reduce the problem from quartic to quadratic optimization we
introduce a latent (auxiliary) variable
$\boldsymbol{\Phi}_{0}\in\boldsymbol{\Omega}^{L\times L}$ as a symmetric
matrix. Define the new objective as

$\displaystyle f(\boldsymbol{\phi},\boldsymbol{\phi}_{0})$
$\displaystyle=\beta\boldsymbol{\phi}^{\mathrm{H}}\bar{\boldsymbol{Q}}(\boldsymbol{\Phi}_{0})\boldsymbol{\phi}+(1-\beta)\boldsymbol{\phi}^{\mathrm{H}}\widetilde{\mathbf{Q}}\boldsymbol{\phi}$
$\displaystyle+2(1-\beta)\operatorname{Re}\left(\mathbf{q}^{\mathrm{H}}\boldsymbol{\phi}\right).$
(14)

In order to have the surrogate objective equal to the original objective, i.e.
$f(\boldsymbol{\phi},\boldsymbol{\phi}_{0})=f(\boldsymbol{\phi})$, we propose
to optimize $f(\boldsymbol{\phi},\boldsymbol{\phi}_{0})$ under the constraint
$\boldsymbol{\phi}=\boldsymbol{\phi}_{0}$, where
$\boldsymbol{\phi}_{0}=\mathrm{vech}\left(\boldsymbol{\Phi}_{0}\right)$.

Additionally, in order to impose the constraint
$\boldsymbol{\Phi}\in\boldsymbol{\Omega}_{{}_{M}}^{{}^{L\times L}}$ in (III),
we introduce a second latent variable,
$\boldsymbol{\Phi}_{1}\in\boldsymbol{\Omega}^{{}^{L\times L}}_{{}_{M}}$ as a
symmetric matrix and impose $\boldsymbol{\Phi}=\boldsymbol{\Phi}_{1}$. Denote
$\boldsymbol{\phi}_{1}=\mathrm{vech}\left(\boldsymbol{\Phi}_{1}\right)$. In
order to impose the unitarity constraint on $\boldsymbol{\Phi}$ in
$\mathcal{P}_{1}$, we introduce a third auxiliary variable $\mathbf{U}$ with
the constraints $\mathbf{U}=\boldsymbol{\Phi}$ and
$\mathbf{U}^{\mathrm{H}}\mathbf{U}=\mathbf{I}$. Introducing these auxiliary
variables into problem $\mathcal{P}_{1}$, we obtain the equivalent problem

$\displaystyle\mathcal{P}_{2}:\;\underset{\boldsymbol{\phi},\boldsymbol{\phi}_{0},\boldsymbol{\phi}_{1},\mathbf{U}}{\textrm{maximize}}$
$\displaystyle\;\;f(\boldsymbol{\phi},\boldsymbol{\phi}_{0})$ subject to
$\displaystyle\;\boldsymbol{\phi}=\boldsymbol{\phi}_{0}=\boldsymbol{\phi}_{1}=\mathrm{vech}\left(\mathbf{U}\right),$
(15) $\displaystyle\mathbf{U}^{\mathrm{H}}\mathbf{U}=\mathbf{I},\;$ (16)
$\displaystyle\;\boldsymbol{\phi},\boldsymbol{\phi}_{0}\in\boldsymbol{\Omega}^{{}^{L(L+1)/2}}\text{and}\;\boldsymbol{\phi}_{1}\in\boldsymbol{\Omega}^{{}^{L(L+1)/2}}_{M}.$
(17)

We solve this non-convex problem using the framework proposed in [30] for
unimodular quadratic optimization.

##  IV Optimization Framework

The Lagrangian of $\mathcal{P}_{2}$ is

$\displaystyle\mathcal{L}(\boldsymbol{\phi},\boldsymbol{\phi}_{0},\boldsymbol{\phi}_{1},\mathbf{U}$
$\displaystyle)=-f(\boldsymbol{\phi},\boldsymbol{\phi}_{0})+\sum_{i=0}^{1}\frac{\rho_{i}}{2}\|\boldsymbol{\phi}_{i}-\boldsymbol{\phi}\|^{2}_{2}$
(18)
$\displaystyle+\frac{\rho_{{}_{2}}}{2}\|\mathrm{vec}\left(\mathbf{U}\right)-\boldsymbol{\phi}\|^{2}_{2},$

where $\rho_{i}>0$, $i\in\\{0,\ldots,3\\}$ are penalty parameters. Next, we
employ a penalty based framework and determine
$(\boldsymbol{\phi},\boldsymbol{\phi}_{0},\boldsymbol{\phi}_{1},\mathbf{U})$
in an alternating iterative fashion to minimize
$\mathcal{L}(\boldsymbol{\phi},\boldsymbol{\phi}_{0},\boldsymbol{\phi}_{1},\mathbf{U})$
under the constraints in (16)-(17).

Update of $\boldsymbol{\phi}$: With fixed
$\boldsymbol{\phi}^{(r)}_{0},\boldsymbol{\phi}^{(r)}_{1},\mathbf{U}^{(r)}$ at
iteration r we need to minimize
$\mathcal{L}(\boldsymbol{\phi},\boldsymbol{\phi}^{(r)}_{0},\boldsymbol{\phi}^{(r)}_{1},\mathbf{U}^{(r)})$
w.r.t $\boldsymbol{\phi}$, which from (III)-(18), is equivalent to

$\displaystyle\underset{\boldsymbol{\phi}\in\boldsymbol{\Omega}^{{}^{L(L+1)/2}}}{\textrm{minimize}}\;\;\boldsymbol{\phi}^{\mathrm{H}}\mathbf{R}^{(r)}\boldsymbol{\phi}+\operatorname{Re}\left(\mathbf{c}^{(r)\mathrm{H}}\boldsymbol{\phi}\right)$
(19)

where
$\mathbf{R}^{(r)}=-\beta\bar{\mathbf{Q}}(\boldsymbol{\phi}_{0}^{(r)})+(\beta-1)\widetilde{\mathbf{Q}}$
and
$\mathbf{c}^{(r)}=2(\beta-1)\mathbf{q}-\sum_{i=0}^{1}\rho_{i}\boldsymbol{\phi}_{i}^{(r-1)}-\rho_{{}_{2}}\mathrm{vec}\left(\mathbf{U}\right)$.
(19) is equivalent to

$\underset{\bar{\boldsymbol{\phi}}\in\boldsymbol{\Omega}^{{}^{L(L+1)/2+1}}}{\textrm{minimize}}\quad\bar{\boldsymbol{\phi}}^{\mathrm{H}}\mathbf{A}^{(r)}\bar{\boldsymbol{\phi}},$
(20)

where $\mathbf{A}^{(r)}=\begin{bmatrix}\mathbf{R}^{(r)}&\mathbf{c}^{(r)}/2\\\
\mathbf{c}^{(r)\mathrm{H}}/2&0\end{bmatrix}$ and
$\bar{\boldsymbol{\phi}}=[\boldsymbol{\phi}^{\top}\;1]^{\top}$. This problem
with feasible set $\boldsymbol{\Omega}^{{}^{L(L+1)/2+1}}$ is a unimodular
quadratic program for which the following iterations will lead to non-
increasing objective

$\bar{\boldsymbol{\phi}}^{(r,s+1)}=e^{\mathrm{j}\operatorname{arg}\left(\bar{\mathbf{A}}^{(r)}\bar{\boldsymbol{\phi}}^{(r,s)}\right)}.$
(21)

Note that $s$ denotes the iteration number over (21) and we used the diagonal
loading technique, i.e.,
$\bar{\mathbf{A}}^{(r)}=\lambda_{m}\mathbf{I}-\mathbf{A}^{(r)}$, where the
loading parameter $\lambda_{m}\geq\lambda_{max}(\mathbf{A}^{(r)})$ [31].

Update of $\boldsymbol{\phi}_{0}$: The Lagrangian function
$\mathcal{L}(\boldsymbol{\phi}^{(r)},\boldsymbol{\phi}_{0},\boldsymbol{\phi}_{1}^{(r-1)},\mathbf{U}^{(r-1)})$
is quadratic w.r.t $\boldsymbol{\phi}_{0}$. To show this, from (11) we write
$\boldsymbol{\phi}^{(r)\mathrm{H}}\bar{\boldsymbol{Q}}(\boldsymbol{\Phi}_{0})\boldsymbol{\phi}^{(r)}=\boldsymbol{\phi}_{0}^{\mathrm{H}}\bar{\bar{\mathbf{Q}}}(\boldsymbol{\Phi}^{(r)})\boldsymbol{\phi}_{0}$,
where

$\bar{\bar{\mathbf{Q}}}(\boldsymbol{\phi}^{(r)})=\frac{|\alpha_{T}|^{2}}{\sigma^{2}_{r}}\left(\mathbf{G}^{\top}\boldsymbol{\Phi}^{(r)\top}\mathbf{a}\mathbf{a}^{\top}\otimes\mathbf{G}^{\top}\right)^{\mathrm{H}}\bar{\mathbf{P}}\left(\mathbf{G}^{\top}\boldsymbol{\Phi}^{(r)\top}\mathbf{a}\mathbf{a}^{\top}\otimes\mathbf{G}^{\top}\right)$
(22)

and
$\boldsymbol{\Phi}^{(r)}=\mathrm{vec}_{{}_{L,L}}^{-1}\left(\mathbf{D}_{{}_{L}}\boldsymbol{\phi}^{(r)}\right)$.
For the sake of brevity, we denoted
$\mathbf{a}(\theta_{h_{t}},\theta_{v_{t}})$ by $\mathbf{a}$. By substituting
(22) in (18), the minimization of the Lagrangian w.r.t $\boldsymbol{\phi}_{0}$
is equivalent to

$\displaystyle\underset{\boldsymbol{\phi}_{0}\in\boldsymbol{\Omega}^{{}^{L(L+1)/2}}}{\textrm{minimize}}\;\;\boldsymbol{\phi}_{0}^{\mathrm{H}}\bar{\mathbf{R}}^{(r)}\boldsymbol{\phi}_{0}+\operatorname{Re}\left(\bar{\mathbf{c}}^{(r)\mathrm{H}}\boldsymbol{\phi}_{0}\right),$
(23)

where
$\bar{\mathbf{R}}^{(r)}=-\beta\bar{\bar{\mathbf{Q}}}(\boldsymbol{\phi}^{(r)})$
and $\bar{\mathbf{c}}^{(r)}=-\rho_{0}\boldsymbol{\phi}^{(r)}$. Similar to the
update for $\boldsymbol{\phi}$, by letting
$\hat{\mathbf{A}}^{(r)}=\begin{bmatrix}\bar{\mathbf{R}}^{(r)}&\bar{\mathbf{c}}^{(r)}/2\\\
\bar{\mathbf{c}}^{(r)\mathrm{H}}/2&0\end{bmatrix}$ and
$\bar{\boldsymbol{\phi}}_{0}=[\boldsymbol{\phi}_{0}^{\top}\;1]^{\top}$,
similar power method like iterations as in (21), will be applicable to the
problem. We have

$\bar{\boldsymbol{\phi}}_{0}^{(r,t+1)}=e^{\mathrm{j}\operatorname{arg}\left(\widehat{\widehat{\mathbf{A}}}^{(r)}\bar{\boldsymbol{\phi}}_{0}^{(r,t)}\right)}.$
(24)

where $t$ is the number for power method like iterations and
$\widehat{\widehat{\mathbf{A}}}^{(r)}=\lambda_{m}\mathbf{I}-\widehat{\mathbf{A}}^{(r)}$,
with the loading parameter is chosen such that
$\lambda_{m}\geq\lambda_{max}(\widehat{\mathbf{A}}^{(r)})$.
Update of $\boldsymbol{\phi}_{1}$: Given
$\boldsymbol{\phi}^{(r)},\boldsymbol{\phi}^{(r)}_{0},\mathbf{U}^{(r-1)}$, we
minimize
$\mathcal{L}(\boldsymbol{\phi}^{(r)},\boldsymbol{\phi}^{(r)}_{0},\boldsymbol{\phi}_{1},\mathbf{U}^{(r-1)})$
under the finite alphabet constraint
$\boldsymbol{\phi}_{1}\in\boldsymbol{\Omega}_{M}^{L^{2}\times 1}$.
Analogously, we solve

$\displaystyle\underset{\boldsymbol{\phi}_{1}\in\boldsymbol{\Omega}^{{}^{L(L+1)/2}}_{M}}{\textrm{maximize}}\;\;\operatorname{Re}\left(\hat{\mathbf{c}}^{(r)\mathrm{H}}\boldsymbol{\phi}_{1}\right),$
(25)

where $\hat{\mathbf{c}}^{(r)}=\rho_{1}\boldsymbol{\phi}^{(r)}$. Rewrite this
problem as

$\displaystyle\underset{\gamma_{l}}{\textrm{maximize}}\;\;\cos{(\gamma_{l}-\boldsymbol{\psi}_{l})}$
$\displaystyle\;\textrm{subject
to}\;\;\gamma_{l}\in\left\\{1,\frac{2\pi}{M},\cdots,\frac{2\pi(M-1)}{M}\right\\}\;\;,$
(26)

where $\gamma_{l}$ is the phase of $\boldsymbol{\phi}_{1,l}$ and
$\boldsymbol{\psi}_{l}$ is the phase of $\hat{\mathbf{c}}^{(r)}_{l}$. The
solution is obtained by an exhaustive search for $m$ such that
$m^{(r)}_{l}=\underset{m}{\textrm{argmax}}\;\;\cos{(2\pi\frac{m}{M}-\boldsymbol{\psi}_{l})}$
$m\in\\{0,\ldots,M-1\\}$. Then, for $l\in\\{1,\ldots,L^{2}\\}$, the update of
the $l$-th element in $\boldsymbol{\phi}_{1}^{(r)}$ is

$\boldsymbol{\phi}_{1_{l}}^{(r)}=e^{\mathrm{j}2\pi\frac{m^{(r)}_{l}}{M}}.$
(27)

Update of $\mathbf{U}$: Consider the following problem

$\displaystyle\mathcal{P}_{2}:$
$\displaystyle\;\underset{\mathbf{U}}{\textrm{minimize}}\;\;\|\boldsymbol{\Phi}^{(r)}-\mathbf{U}\|_{{}_{2}}$
$\displaystyle\;\textrm{subject
to}\;\;\mathbf{U}^{\mathrm{H}}\mathbf{U}=\mathbf{I}.$ (28)

Consider the singular value decomposition (SVD) of $\boldsymbol{\Phi}^{(r)}$
as $\boldsymbol{\Phi}^{(r)}=\mathbf{U}_{1}\Sigma\mathbf{U}_{2}^{\mathrm{H}}$,
then the solution of $\mathcal{P}_{2}$ is given by [32]

$\mathbf{U}^{(r)}=\mathbf{U}_{2}\mathbf{U}_{1}^{\mathrm{H}}.$ (29)

Algorithm 1 summarizes the overall design procedure of BD-RIS scattering
matrix.

Algorithm 1 Optimal design of BD-RIS scattering matrix for ISAC.

1:Input Initializations $\boldsymbol{\phi}^{(0)}$,
$\boldsymbol{\phi}_{0}^{(0)}$, $\boldsymbol{\phi}_{1}^{(0)}$,
$\mathbf{U}^{(0)}$, $r=1$.

2:Set $s=1$, while
$\|\bar{\boldsymbol{\phi}}^{(r,s)}-\bar{\boldsymbol{\phi}}^{(r,s+1)}\|>\epsilon$
do update $\bar{\boldsymbol{\phi}}^{(r,s)}$ according to (21) and $s=s+1$.

3:$\boldsymbol{\phi}^{(r)}=[\mathbf{I}_{{}_{L(L+1)/2}}\;\;\boldsymbol{0}]\bar{\boldsymbol{\phi}}^{(r,s)}$.

4:Set $t=1$, while
$\|\bar{\boldsymbol{\phi}}_{0}^{(r,t)}-\bar{\boldsymbol{\phi}}_{0}^{(r,t+1)}\|>\epsilon_{1}$
do update $\bar{\boldsymbol{\phi}}_{0}^{(r,t)}$ according to (24) and $t=t+1$.

5:Update $\boldsymbol{\phi}_{1}^{(r)}$ and $\mathbf{U}^{(r)}$ as in (27) and
(29), respectively.

6:$r=r+1$.

7:Repeat steps $1$-$5$ until a prescribed stop criterion is satisfied, e.g.,
$\|\boldsymbol{\phi}^{(r)}-\boldsymbol{\phi}^{(r+1)}\|\leq\epsilon_{2}$.

8:Output
$\boldsymbol{\Phi}^{*}=\mathrm{vec}_{{}_{L,L}}^{-1}\left(\mathbf{D}_{{}_{L}}\boldsymbol{\phi}^{(r)}\right)$.

######  Remark 1.

For continuous-valued scattering matrices i.e.
$\boldsymbol{\phi}\in\boldsymbol{\Omega}^{L^{2}}$, it is sufficient to remove
the latent variable $\boldsymbol{\phi}_{1}$ from $\mathcal{P}_{2}$ and solve
it with respect to $\boldsymbol{\phi}$, $\boldsymbol{\phi}_{0}$ and
$\mathbf{U}$.

######  Remark 2.

Contrary to BD-RIS above, D-RIS design [16] is less complex because it
introduces a single latent variable to simplify quartic optimization into a
quadratic alternative.

##  V Numerical Experiments

(a)  (b)

Figure 2: (a) SNR${}_{{}_{T}}$ is plotted against the increasing number of
elements in BD-RIS. (b) Performance of D-RIS and BD-RIS in terms of SNRc and
SNRr vs. different values of $\beta$.

To evaluate our proposed algorithm, we assume the DFBS located at $(0,0,0)$m,
in the 3-D Cartesian plane, is equipped with $N=8$ Tx and Rx antennas with
half-wavelength antenna spacing. The ISAC system is serving K = 4 single-
antenna users located at $(20,20,0)$ m and simultaneously detects a Swerling 0
point-like target with complex reflectivity coefficient
$\alpha_{T}\sim\mathcal{CN}(0,1)$ and located at $(20,20,0)$ m. The ISAC
system is being assisted by a BD-RIS located at $(30,30,0)$m. The variance of
noise at communication users and sensing Rx is set at
$\sigma^{2}_{c}=\sigma^{2}_{r}=-100$ dBm. We assumed an i.i.d Rayleigh fading
channel in the DFBS to BD-RIS, BD-RIS to users and DFBS to users paths which
we abbreviate by RD,UR and UD, respectively. The path loss is modelled as
$p_{ij}=p_{0}d^{-\alpha_{ij}}$, where $p_{0}$ is the path loss at distance
$1$m and $ij\in\\{\text{RD, UR, UD}\\}$. We set $p_{0}=-30$ dB,
$\alpha_{{}_{\text{RD}}}=2$, $\alpha_{{}_{\text{UD}}}=4$ and
$\alpha_{{}_{\text{UR}}}=2.8$ for the paths that signal propagates through.
The CSI matrices are generated according to
$\mathbf{G}\sim\mathcal{CN}(0,p_{{}_{\text{RD}}}\mathbf{I})$,
$\mathbf{H}\sim\mathcal{CN}(0,p_{{}_{\text{UR}}}\mathbf{I})$ and
$\mathbf{F}\sim\mathcal{CN}(0,p_{{}_{\text{UD}}}\mathbf{I})$.

Fig. 2a illustrates the performance of Algorithm 1 for different resolutions
of BD-RIS with discrete valued optimized phases in comparison with the
continuous phase and no RIS. We set $\beta=0.5$, and manually tuned the
penalty parameters $\rho_{i}$, $i\in\\{0,\ldots,3\\}$ until convergence. To
design BD-RIS with continuous phases, we adapt Algorithm 1 as per Remark 1.
Fig. 2b shows that fully-connected BD-RIS outperforms conventional D-RIS in
terms of the optimized SNRr and SNRc. Here, D-RIS is employed with phases
optimized according to Remark 2. Moreover, changing $\beta$ in the interval
$[10^{-3},1]$ tunes our design criterion. Increasing $\beta$ shifts the system
from communications-centric to radar-centric as evidenced from increasing
(decreasing) SNRr (SNRc).

## References

  * [1] J. A. Hodge, K. V. Mishra, B. M. Sadler, and A. I. Zaghloul, “Index-modulated metasurface transceiver design using reconfigurable intelligent surfaces for 6G wireless networks,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 17, no. 6, pp. 1248–1263, 2023.
  * [2] J. A. Hodge, K. V. Mishra, and A. I. Zaghloul, “Intelligent time-varying metasurface transceiver for index modulation in 6G wireless networks,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 19, no. 11, pp. 1891–1895, 2020.
  * [3] S. P. Chepuri, N. Shlezinger, F. Liu, G. C. Alexandropoulos, S. Buzzi, and Y. C. Eldar, “Integrated sensing and communications with reconfigurable intelligent surfaces: From signal modeling to processing,” _IEEE Signal Processing Magazine_ , vol. 40, no. 6, pp. 41–62, 2023.
  * [4] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 8, pp. 4157–4170, 2019.
  * [5] W. Shi, J. Xu, W. Xu, C. Yuen, A. L. Swindlehurst, and C. Zhao, “On secrecy performance of RIS-assisted MISO systems over rician channels with spatially random eavesdroppers,” _IEEE Transactions on Wireless Communications_ , 2024, in press.
  * [6] K. V. Mishra, A. Chattopadhyay, S. S. Acharjee, and A. P. Petropulu, “OptM3Sec: Optimizing multicast IRS-aided multiantenna DFRC secrecy channel with multiple eavesdroppers,” in _IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2022, pp. 9037–9041.
  * [7] Z. Esmaeilbeig, K. V. Mishra, and M. Soltanalian, “IRS-aided radar: Enhanced target parameter estimation via intelligent reflecting surfaces,” in _IEEE Sensor Array and Multichannel Signal Processing Workshop_ , 2022, pp. 286–290.
  * [8] E. Shtaiwi, H. Zhang, A. Abdelhadi, A. L. Swindlehurst, Z. Han, and H. V. Poor, “Sum-rate maximization for RIS-assisted integrated sensing and communication systems with manifold optimization,” _IEEE Transactions on Communications_ , vol. 71, no. 8, pp. 4909–4923, 2023.
  * [9] R. Liu, M. Li, Y. Liu, Q. Wu, and Q. Liu, “Joint transmit waveform and passive beamforming design for RIS-aided DFRC systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 5, pp. 995–1010, 2022.
  * [10] K. V. Mishra, M. R. B. Shankar, V. Koivunen, B. Ottersten, and S. A. Vorobyov, “Toward millimeter-wave joint radar communications: A signal processing perspective,” _IEEE Signal Process. Mag._ , vol. 36, no. 5, pp. 100–114, 2019.
  * [11] A. Liu, Z. Huang, M. Li, Y. Wan, W. Li, T. X. Han, C. Liu, R. Du, D. K. P. Tan, J. Lu, Y. Shen, F. Colone, and K. Chetty, “A survey on fundamental limits of integrated sensing and communication,” _IEEE Communications Surveys & Tutorials_, vol. 24, no. 2, pp. 994–1034, 2022.
  * [12] F. Liu, Y. Cui, C. Masouros, J. Xu, T. X. Han, Y. C. Eldar, and S. Buzzi, “Integrated sensing and communications: Toward dual-functional wireless networks for 6G and beyond,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 6, pp. 1728–1767, 2022.
  * [13] X. Hu, C. Liu, M. Peng, and C. Zhong, _Reconfigurable intelligent surface-enabled integrated sensing and communication in 6G_. Springer, 2024.
  * [14] R. Liu, M. Li, Q. Liu, and A. L. Swindlehurst, “SNR/CRB-constrained joint beamforming and reflection designs for RIS-ISAC systems,” _arXiv preprint arXiv:2301.11134_ , 2023.
  * [15] T. Wei, L. Wu, K. V. Mishra, and M. B. Shankar, “Multi-IRS-aided Doppler-tolerant wideband DFRC system,” _IEEE Transactions on Communications_ , vol. 71, no. 11, pp. 6561–6577, 2023.
  * [16] Z. Esmaeilbeig, A. Eamaz, K. V. Mishra, and M. Soltanalian, “Quantized phase-shift design of active IRS for integrated sensing and communications,” in _IEEE International Conference on Acoustics, Speech, and Signal Processing Workshops_ , 2023, pp. 1–5.
  * [17] M. Hua, Q. Wu, C. He, S. Ma, and W. Chen, “Joint active and passive beamforming design for IRS-aided radar-communication,” _IEEE Transactions on Wireless Communications_ , vol. 22, no. 4, pp. 2278–2294, 2022\.
  * [18] S. Shen, B. Clerckx, and R. Murch, “Modeling and architecture design of reconfigurable intelligent surfaces using scattering parameter network analysis,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 2, pp. 1229–1243, 2022.
  * [19] M. Nerini, S. Shen, and B. Clerckx, “Closed-form global optimization of beyond diagonal reconfigurable intelligent surfaces,” _IEEE Transactions on Wireless Communications_ , vol. 23, no. 2, pp. 1037–1051, 2024.
  * [20] H. Li, S. Shen, and B. Clerckx, “Beyond diagonal reconfigurable intelligent surfaces: From transmitting and reflecting modes to single-, group-, and fully-connected architectures,” _IEEE Transactions on Wireless Communications_ , vol. 22, no. 4, pp. 2311–2324, 2023.
  * [21] Q. Wu, T. Q. Duong, D. W. K. Ng, R. Schober, and R. Zhang, _Intelligent surfaces empowered 6G wireless network_. John Wiley & Sons, 2023.
  * [22] T. Fang and Y. Mao, “A low-complexity beamforming design for beyond-diagonal RIS aided multi-user networks,” _IEEE Communications Letters_ , vol. 28, no. 1, pp. 203–207, 2024.
  * [23] B. Wang, H. Li, S. Shen, Z. Cheng, and B. Clerckx, “A dual-function radar-communication system empowered by beyond diagonal reconfigurable intelligent surface,” _arXiv preprint arXiv:2301.03286_ , 2023.
  * [24] I. Santamaria, M. Soleymani, E. Jorswieck, and J. Gutiérrez, “SNR maximization in beyond diagonal RIS-assisted single and multiple antenna links,” _IEEE Signal Processing Letters_ , vol. 30, pp. 923–926, 2023.
  * [25] J. A. Zhang, F. Liu, C. Masouros, R. W. Heath, Z. Feng, L. Zheng, and A. Petropulu, “An overview of signal processing techniques for joint communication and radar sensing,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 15, no. 6, pp. 1295–1315, 2021.
  * [26] J. Liu, K. V. Mishra, and M. Saquib, “Co-designing statistical MIMO radar and in-band full-duplex multi-user MIMO communications,” _arXiv preprint arXiv:2006.14774_ , 2020.
  * [27] A. R. Chiriyath, B. Paul, and D. W. Bliss, “Radar-communications convergence: Coexistence, cooperation, and co-design,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 3, no. 1, pp. 1–12, 2017.
  * [28] M. Nerini, S. Shen, and B. Clerckx, “Discrete-value group and fully connected architectures for beyond diagonal reconfigurable intelligent surfaces,” _IEEE Transactions on Vehicular Technology_ , vol. 72, no. 12, pp. 16 354–16 368, 2023.
  * [29] M. I. Skolnik, _Radar handbook_ , 3rd ed. McGraw-Hill, 2008.
  * [30] X. Yu, G. Cui, J. Yang, J. Li, and L. Kong, “Quadratic optimization for unimodular sequence design via an ADPM framework,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 3619–3634, 2020.
  * [31] M. Soltanalian and P. Stoica, “Designing unimodular codes via quadratic optimization,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 5, pp. 1221–1234, 2014.
  * [32] H. He, P. Stoica, and J. Li, “Designing unimodular sequence sets with good correlations – Including an application to MIMO radar,” _IEEE Transactions on Signal Processing_ , vol. 57, no. 11, pp. 4391–4405, 2009.


# Can the Blandford-Znajek mechanism power steady jets?

A.R. King School of Physics and Astronomy, University of Leicester, Leicester,
LE1 7RH, UK  Astronomical Institute Anton Pannekoek, University of Amsterdam,
Science Park 904, NL-1098 XH Amsterdam, Netherlands  Leiden Observatory,
Leiden University, Niels Bohrweg 2, NL-2333 CA Leiden, Netherlands  J.E.
Pringle Institute of Astronomy, University of Cambridge, Madingley Road,
Cambridge CB3 OHA, UK

###### Abstract

We consider the Blandford-Znajek (BZ) mechanism for extracting black hole spin
energy to drive astrophysical jets. Analyses of the BZ mechanism generally
take no account of any electric charge on the black hole. But, as noted by
Wald and others, if the medium surrounding the black hole is an ionised plasma
with mobile charges, then a spinning hole quickly acquires an electric charge.
The effect of this charge is to nullify the electric field structures which
drive the BZ mechanism. Since jets are now observed in a wide variety of
classes of accreting objects, most of which do not contain a central black
hole, it seems likely that the jet driving mechanism in all astrophysical
objects uses energy directly from the accretion disc, rather than black hole
spin.

astrophysical jets – accretion discs – black hole physics

##  1 Introduction

Early maps of double–lobe radio galaxies (e.g. Mitton & Ryle, 1969) showed
amorphous blobs of radio emission symmetrically placed each side of the
central galaxy. Rees (1971) suggested that an unknown object in the galactic
nucleus channels energy to the radio lobes through jets. There is now almost
universal agreement that the central object in radio galaxies is a
supermassive black hole (Rees 1984), and that the high energy activity in the
nucleus is powered by accretion (Salpeter, 1964), most likely through an
accretion disc (Lynden–Bell, 1969). Thus the relativistic jets seen to emanate
from the nucleus are ultimately powered by accretion of matter on to the
central black hole (see the reviews by Begelman, Blandford & Rees, 1984;
Heckman & Best, 2014; Blandford, Meier & Readhead, 2019).

Blandford & Znajek (1977, hereafter BZ77) proposed a radical new mechanism in
which the jets from galactic nuclei are powered by direct electromagnetic
extraction of the spin energy from the central black hole. This is the
Blandford–Znajek (BZ) mechanism.

In Section 2.1 we provide a brief overview of this mechanism, which invokes a
spinning black hole situated in an aligned magnetic field. We note that the
central black hole is tacitly assumed always to have negligible electric
charge, that is, any current entering the hole is assumed to be balanced by a
precisely opposing current entering elsewhere, so that there is an effective
current through the hole. In Section 2.2 we draw attention to the analysis by
Wald (1974; see also King et al., 1975; Petterson, 1975; Gibbons et al., 2013)
in which he finds that if the black hole is surrounded by a plasma which
contains mobile charges, it will acquire a specific electric charge. The
charge is exactly such as to render the BZ mechanism inoperative. We provide a
brief discussion in Section 3.

##  2 The Blandford–Znajek mechanism

###  2.1 The basic mechanism; the uncharged black hole

Black holes are generally assumed to have zero net electric charge. The reason
is that if for example a Schwarzschild (non–rotating) black hole is given a
charge, the neighbourhood of the hole then acquires an electric field. If (and
only if) charge separation is allowed, charged particles in the surrounding
astrophysical plasma move parallel or antiparallel to the electric field. In
this way the black hole selectively acquires charges so that it moves quickly
towards zero net charge. Charge separation is well known to occur in
electrical media in which the charge carriers are able to move independently.
It occurs for example in electrolytes (Debye & Hückel, 1923) and in ionized
astrophysical plasmas (Salpeter, 1954), and it affects the rates of nuclear
burning in stars (e.g. Clayton, 1968). The timescale to reach zero net charge
is typically quite short in realistic astrophysical environments (see, for
example, the discussion by Cardoso et al., 2016). Thus, Blandford (1987)
asserts that “charged, Kerr–Newman black holes are irrelevant to astronomy”.

Soon after black holes were recognised as a realistic astrophysical
possibility, Wald (1974) devised an elegant method to compute the effect on
the electromagnetic field structure when a black hole, with zero net charge,
is placed in a uniform magnetic field. King et al. (1975) extended this result
to consider more general aligned magnetic field structures. Wald (1974) showed
that for a spinning (Kerr) black hole, with zero net electric charge, and with
spin aligned with an external uniform magnetic field ${\bf B}$ in a vacuum,
the spin induces an electric field with ${\bf E\cdot B\neq 0}$ (for a
classical analogy, see Ruffini & Treves, 1973). Wald assumed a field uniform
at infinity, but King et al. (1975) showed that the field structure near the
hole is essentially the same for any realistic aligned field. Wald also noted
that an electric field with $\bf E\cdot B\neq 0$ would lead to movement of any
external charged particles. King et al. (1975) showed that the induced
electric field structure is such that $\bf E\cdot B$ has opposite signs near
the poles and near the equator of the spinning hole. This implies that charges
of one sign would be attracted towards the poles, whereas charges of the
opposite sign would be attracted to a band near the equator.

This result led BZ77 to consider a spinning black hole at the centre of an
accretion disc with the disc providing the required currents to give rise to
an aligned magnetic field. The Bardeen–Petterson effect (Bardeen & Petterson,
1975) ensures that the central disc and black hole spin are generally aligned,
and so by symmetry the magnetic field should also be aligned with the spin.
BZ77 then noted that if the black hole is surrounded by a conducting medium
the structure of the induced electric field found by Wald (1974) and by King
et al. (1975) must produce an effective electric current through the black
hole, returning through the surrounding medium (see, for example, Thorne &
Blandford, 1982; Blandford, 1987). Any dissipation of this current’s energy in
the surrounding material then taps the spin energy of the black hole.

This is the BZ–mechanism. It provides, in principle, a mechanism for the
continuous extraction of energy from a spinning black hole. BZ77 speculated
that it could be used to power the astrophysical jets seen to emanate from
active galactic nuclei (see also Blandford, Meier & Readhead, 2019). The
luminosity $L_{\rm BZ}$ produced by this process is, on dimensional grounds,

$L_{\rm BZ}\sim a^{2}\times\frac{B^{2}}{8\pi}\times\frac{4\pi R_{\rm
s}^{3}}{3}\times\frac{c}{R_{s}}.$ (1)

for a black hole of mass $M$, Schwarzschild radius $R_{\rm s}\sim 2GM/c^{2}$,
dimensionless spin parameter $a$, where $0\leq a^{2}\leq 1$, placed in a
magnetic field of strength $B$, i.e. approximately the magnetic energy
contained by the formal ‘volume’ of the hole, emitted every light crossing
time of the hole, and moderated by the amount of spin energy available.

For a ‘fiducial’ field strength (Begelman, Blandford & Rees, 1984; but see
Ghosh & Abramowicz, 1997) of

$B\sim 2\times 10^{4}M_{8}^{-1/2}\>{\rm G},$ (2)

this gives a luminosity

$L_{\rm BZ}\sim 2\times 10^{45}a^{2}M_{8}\>{\rm erg\,s^{-1}},$ (3)

where $M_{8}$ is the mass of the black hole in units of $10^{8}M_{\odot}$.
This is comparable to the Eddington luminosity

$L_{E}\sim 1.3\times 10^{46}M_{8}\>{\rm erg\,s^{-1}}.$ (4)

BZ77 investigated this process using the force–free approximation in a
charge–separated plasma in which particle inertia and interparticle collision
terms can be ignored (see also Komissarov, 2004). Since then a number of
authors have investigated this mechanism assuming that the surrounding medium
could be modelled using the MHD approximation in which collision terms
dominate, often by numerical means (see for example the reviews by Davis &
Tchekhovskoy, 2020; Komissarov & Porth, 2021). However, in all these
investigations, it is implicitly assumed that the movements of charges in any
surrounding plasma do not permit a change in the net electrical charge of the
black hole.

###  2.2 The acquisition of charge and its implication

The discussion of Section 2.1 above assumes, as is usually the case in
astronomy, that the spinning black hole has negligible net electric charge.
However, in his seminal paper, Wald (1974) reasoned that a black hole would be
surrounded by a standard astrophysical plasma in which the possibility of
charge separation would exist. In this case, it is to be expected that the
electric field drives the charges in such a way as to lead to a drop in
electric potential along the magnetic field lines (cf. Komissarov, 2004). He
argued further that for a spinning black hole in an aligned magnetic field
${\bf B}$ the movements of charge carriers (i.e. currents) induced by the
electric field with ${\bf E\cdot B\neq 0}$, together with the fact that the
charge carriers are individually mobile, would lead to the hole selectively
accreting net charge in such a way as to nullify the effects of the electric
field.

Specifically, Wald showed that the net charge reaches the value

$Q=2BJ$ (5)

in geometrized units (where $J=Ma$ is its total angular momentum). Then the
charge on the hole remains constant, removing the need for currents, which
might then drive the BZ effect.

This result was confirmed and generalised by Petterson (1975), who showed that
the precise value $Q^{\prime}$ of the critical charge in units of $Q$ depends
on the distribution of the source currents of the magnetic field.

The charge $Q^{\prime}$ is utterly negligible ($Q^{\prime 2}\ll M^{2}$) in
gravitational terms (Wald, 1974; Zajacek et al, 2018; Zajacek & Tursunov,
2019), so the spacetime metric is still to a high approximation uncharged
Kerr, in agreement with the remark by Blandford (1987) quoted above. In
particular the motion of uncharged particles is effectively identical to the
case $Q^{\prime}=0$. But charged particle motion is very different, and
strongly influenced by the charge $Q^{\prime}$. This is another illustration
of how extremely weak gravity is by comparison with electromagnetism. Thus, if
the surrounding conducting medium is treated as a realistic space plasma which
permits a net flow of charge into the black hole, then rather than providing a
continuous process for removing spin energy from the hole, the induced
electric currents are an initial transient effect which continues only until
the black hole acquires the charge $Q^{\prime}\simeq 2BJ$ (cf Zajacek et al.,
2018; Zajacek & Tursunov 2019)111For all conceivable boundary conditions far
from the hole its net charge tends monotonically to $Q^{\prime}$..

The energy released in this transient is, using the above numbers,

$E_{Q}\sim L_{\rm BZ}\times\frac{R_{\rm s}}{c}\sim 2\times
10^{48}a^{2}M_{8}^{2}\>{\rm erg},$ (6)

i.e. the transient emits the luminosity $L_{\rm BZ}$ for a time $\sim R_{\rm
s}/c\sim 10^{3}M_{8}\,{\rm s}$. Once the black hole acquires this charge, the
torque on it vanishes and no more spin energy can be extracted.

More recently numerical modelling of the BZ–mechanism has been undertaken
using Particle–In–Cell (PIC) plasma methods, which permit independent mobility
of individual charges. In principle these techniques should be able to test
Wald’s fundamental hypothesis that a spinning black hole immersed in a
magnetic field should acquire a net charge. The same technique has also been
applied to ionised plasma surrounding rotating neutron stars (e.g.
Kalapotharakos et al., 2018). However, in contrast to the neutron star case,
where Kalapotharakos et al. (2018) treat the inner boundary with some care,
noting that they ensure “current closure of charge carriers that reach the
stellar surface”, modelling of the black hole case is often done with inner
boundary conditions which either prevent or ignore the acquisition of charge
by the black hole. For example, Parfrey et al (2019, see also Crinquand et
al., 2019) do not comment on charge acquisition, and Hirotani et al. (2021)
impose that ${\bf E\cdot B}=0$ and that both the radial component of the
electric field and the meridional component of the magnetic field vanish at
the inner boundary222Note added in proof: Parfrey (private communication)
informs us that, during the timespan of the simulations presented in Parfrey
et al (2019), the black hole both acquires an electric charge and exhibits an
outflow of electromagnetic energy..

##  3 Discussion

The BZ mechanism is a standard cited mechanism for producing steady jets in
objects that contain accreting black holes – galactic nuclei and some X–ray
binaries. And indeed the process is often cited in papers which concern
numerical MHD simulations of jets and outflows produced by magnetic accretion
discs (see the review by Davis & Tchekhovskoy, 2020). We have argued above, in
line with the original ideas of Wald (1974; see also Gibbons et al., 2013)
that in a realistic space plasma which permits a net flow of charge into the
black hole the BZ mechanism cannot tap the spin energy of a black hole
continuously, and is therefore not a viable mechanism for powering continuous
astrophysical jets.

This is primarily because the conducting medium surrounding the black hole
should be treated as an astrophysical plasma with mobile charge carriers. When
charge separation is allowed, along with the possibility of a net flux of
charge into the black hole, any spinning black hole quickly acquires the net
electrical charge $Q^{\prime}$, and the electric fields which drive the
currents required for the BZ mechanism are nullified. The same process of
charge separation which ensures that a non–rotating black hole has zero charge
also ensures that a rotating hole acquires the charge $Q^{\prime}$ that makes
the BZ mechanism inoperative. We have noted that these ideas need to be
tested, for example using PIC plasma simulation techniques. For example, it
might be that collective plasma effects serve to counteract the tendency of
the black hole to acquire charge333We thank the referee for stressing this
possibility..

There are, of course, many other kinds of astrophysical objects which do not
contain black holes and nevertheless produce jets (Burgarella et al., 1993;
Smith, 2012). The jets emitted by young stellar objects are particularly
spectacular (see the review by Ray & Ferreira, 2021). Thus the application of
Occam’s razor444‘Do not multiply hypotheses’, or put simply, ‘don’t invent two
theories for the same thing’. has long suggested that the BZ mechanism, even
if it were viable, is in fact not required for the production of astrophysical
jets (Livio, 1997; see also Pringle, 1993; Price, Pringle & King, 2003). In
addition, Russell, Gallo & Fender (2013) have shown that the jet production
mechanism in binary X–ray sources is not consistent with the prediction of the
BZ mechanism that the jet power should depend on the the square of the
dimensionless jet spin parameter (see equation 3).

Thus, following Occam, if one is forced to choose a single mechanism capable
of producing all astrophysical jets which emanate from accreting objects, then
the most likely choice would be some form of MHD process resulting from
poloidal magnetic field threading the accretion disc (Livio, 1997; Livio et
al., 1999). A mechanism like this is already discussed by Blandford & Znajek
(1977), and early ideas on this process are given by Blandford & Payne (1982)
and, in the protostellar case, by Pudritz & Norman (1983, 1986).

## Acknowledgments

We thank Bob Carswell, Gary Gibbons, Chris Nixon, Colin Norman, Kyle Parfrey,
Roger Blandford and Roman Znajek for helpful comments, and the referee for a
thoughtful report.

##  REFERENCES

Bardeen, J.A. & Petterson, J.A., 1975, ApJL 195, L65

Begelman, M.C., Blandford, R.D., & Rees, M.J., 1984, Rev. Mod. Phys., 56, 255

Blandford, R.D., 1987, in Three Hundred Years of Gravitation, Cambridge
University Press, eds S. Hawking & W. Israel, pp. 277 – 329

Blandford, R.D., Meier, D. & Readhead, A. 2019, ARA&A, 57, 467

Blandford, R.D., & Payne, D.G., 1982, MNRAS, 199, 883

Blandford, R.D., & Znajek, R. L. 1977, MNRAS 179, 433

Burgarella, D., Livio, M., & O’Dea, C.P. (eds), 1993, Astrophysical Jets,
Cambridge University Press

Cardoso, V., Macedo C. F. B., Pani, P., & Ferrari, V., 2016, JCAP, 054

Clayton, D.D., 1968, Principles of Stellar Evolution and Nucleosynthesis,
Univ. of Chicago Press

Crinquand, B., Cerutti, B., Dubus, G., Parfrey, K., Philippov, A., 2021, A&A,
650, A163

Debye, P., & Hückel, E., 1923, Physikalische Zeitschrift, 24, 185

Ghosh, P, & Abramowicz, M.A., 1997, MNRAS, 292, 887

Gibbons, G. W., Mujtaba, A. H., & Pope, C. N. 2013, Classical and Quantum
Gravity, 30, 125008

Heckman, T.M., & Best, P. N., 2014, ARA&A, 52, 589

Hirotani, K., Krasnopolsky, R., Shang, H., Nishikawa, K., Watson, M., 2021,
ApJ, 908, 88

Kalapotharakos, C., Brambilla, G. Timokhin, A. Harding, A. K. & Kazanas, D.,
2018, ApJ, 857, 44

King, A.R., Lasota, J.P., & Kundt, W., 1975, Phys. Rev. D,12, 3037

Komissarov, S.S., 2004, MNRAS, 350, 427

Komissarov, S.S., & Porth, O., 2021, NewAR, 92, 101610

Livio, M., 1997, in Accretion Phenomena and Related Outflows eds D.T.
Wickramasinghe, L. Ferrario, & G.V. Bicknell, ASP Conference Series 121, 845

Livio., M., Ogilvie, G,I., & Pringle, J.E., 1999, ApJ, 512, 100

Lynden–Bell, D., 1969, Nature, 223, 690L

Mitton., S., & Ryle, M., 1969, MNRAS, 146, 221

Parfrey, K., Philippov, A., Cerutti, B., 2019, Phys. Rev. Lett. 122c5101

Petterson, J. A. 1975, Phys Rev D,12, 2218

Price, D.J., Pringle, J.E., & King, A.R., 2003, MNRAS 339, 1223

Pringle, J.E., 1993, in Astrophysical Jets, Cambridge University Press, eds
Burgarella, D., Livio, M., & O’Dea, C.P.

Pudritz, R.E., & Norman, C.A., 1983, ApJ, 274, 677

Pudritz, R.E., & Norman, C.A., 1986, ApJ, 301, 571

Ray, T.P, & Ferreira, J. 2021, NewAR, 93, 101615

Rees, M.J., 1971, Nature, 229, 312

Rees, M.J., 1984, ARA&A, 22, 471

Ruffini, R. & Treves, A. 1973, Astrophys. Lett., 13, 109

Russell, D.M., Gallo, E., & Fender, R.P., 2013, MNRAS, 431, 405

Salpeter, E.E., 1954, Aus J. Phys., 7, 373

Salpeter, E.E., 1964, ApJ, 140, 796

Smith, M.D., 2012, Astrophysical Jets and Beams Cambridge University Press

Thorne, K.S., & Blandford, R.D., 1982, in Proc IAU Symp 97 eds. Heeschen, D.S.
& Wade, C.M., 255

Wald, R.M., 1974, Phys Rev D, 10, 1680

Zajacek, M., Tursunov, A., Eckhart, A., Britze, S., 2018, MNRAS, 480, 4408

Zajacek, M., & Tursunov, A., 2019, The Observatory, 139, 231