# Storage and retrieval of von Neumann measurements via indefinite causal
order structures

Paulina Lewandowska  IT4Innovations, VSB - Technical University of Ostrava,
17. listopadu 2172/15, 708 33 Ostrava, Czech Republic  Ryszard Kukulski∗
Faculty of Physics, Astronomy and Applied Computer Science,
ul. Łojasiewicza 11, Jagiellonian University, 30-348 Kraków, Poland
<EMAIL_ADDRESS>
###### Abstract

This work presents the problem of learning an unknown von Neumann measurement
of dimension $d$ using indefinite causal structures. In the considered
scenario, we have access to $N$ copies of the measurement. We use formalism of
process matrices to store information about the given measurement, that later
will be used to reproduce its best possible approximation. Our goal is to
compute the maximum value of the average fidelity function $F_{d}(N)$ of our
procedure. We prove that $F_{d}(N)=1-\Theta\left(\frac{1}{N^{2}}\right)$ for
arbitrary but fixed dimension $d$. Furthermore, we present the SDP program for
computing $F_{d}(N)$. Basing on the numerical investigation, we show that for
the qubit von Neumann measurements using indefinite causal learning structures
provide better approximation than quantum networks, starting from $N\geq 3$.

##  I Introduction

In the setup of storage and retrieval (SAR) of quantum measurements, we would
like to to approximate a given, unknown measurement, which we were able to
perform $N$ times experimentally. Such strategy usually consists of the notion
of a quantum networks known also as quantum combs [1]. Then, the scheme is
created by preparing some initial quantum state, applying the unknown
operation $N$ times, which stores information about the unknown operation for
later use, and finally, a retrieval operation that produces an approximation
of the black box on some arbitrary quantum state. The scheme is optimal when
it achieves the highest possible fidelity of the approximation [2, 3, 4].

The seminal work in this field was the paper [5]. The Authors have shown that,
in general, the optimal algorithm for quantum measurement learning cannot be
parallel. In [2], whereas, the Authors have discovered the asymptotic
behaviour of the maximum value of the average fidelity function over all
possible learning schemes is given by $1-\Theta\left(\frac{1}{N^{2}}\right)$.

In this work, we introduce a new aspect of von Neumann measurement learning by
using indefinite causal structure theory. The topic of indefinite causal
structures has recently gained traction in quantum information research. This
more general model of computation has the potential to outperform algorithms
based on quantum networks in specific tasks, such as learning or
discriminating between two quantum channels [6, 7, 8]. In the problem of von
Neumann measurements learning, indefinite causal structures find a place in
the storage part of the procedure. Their mathematical description is
formalized in the language of process matrices [9, 10].

As for the results of our work, we will prove that for $2\rightarrow 1$
learning scheme, using indefinite causal structures does not improve the
average fidelity function $F_{d}(2)$ for any dimension $d$. Next, however, we
will show the numerical advantage of using causal structure theory in the
$N\rightarrow 1$ learning scheme for $N\geq 3$. Finally, we determine the
asymptotic behavior of $F_{d}(N)$ for $N\rightarrow\infty$.

This paper is organized as follows. In Section II we formulate $N\rightarrow
1$ learning scheme of von Neumann measurement and necessary mathematical tools
needed to describe this problem. Section III presents a theoretical results of
this paper. In Subsection III.1, we show that the indefinite causal structures
do not improve the average fidelity function of learning for two copies of von
Neumann measurement, whereas Subsection III.2 presents the semidefinite
programs for computing the maximum value of the fidelity function for a finite
number of copies. In Subsection III.3, we also prove the asymptotic behavior
of the fidelity function. Section IV shows a numerical advantage of using the
indefinite causal structure of von Neumann measurements learning over any
strategies based on quantum combs. Finally, the concluding remarks are
presented in Section V.

##  II Problem formulation

This section presents the formulation of learning scheme of an unknown von
Neumann measurement via indefinite causal structures and necessary
mathematical tools needed to describe this problem.

###  II.1 Mathematical framework

Let us introduce the following notation. Consider two complex Euclidean spaces
and denote them by $\mathcal{X},\mathcal{Y}$. By
$\mathrm{L}(\mathcal{X},\mathcal{Y})$ we denote the collection of all linear
mappings of the form $M:\mathcal{X}\rightarrow\mathcal{Y}$. As a shorthand we
put $\mathrm{L}(\mathcal{X})\coloneqq\mathrm{L}(\mathcal{X},\mathcal{X})$. By
$\mathrm{Herm}(\mathcal{X})$ we denote the set of Hermitian operators while
the subset of $\mathrm{Herm(\mathcal{X})}$ consisting of positive semidefinite
operators will be denoted by $\mathrm{Pos}(\mathcal{X})$. The set of quantum
states defined on space $\mathcal{X}$, that is the set of positive
semidefinite operators having unit trace, will be denoted by
$\Omega(\mathcal{X})$. We will also need a linear mapping transforming
$\mathrm{L}(\mathcal{X})$ into $\mathrm{L}(\mathcal{Y})$ as
$\mathcal{T}:\mathrm{L}(\mathcal{X})\mapsto\mathrm{L}(\mathcal{Y}).$ There
exists a bijection between introduced linear mappings $\mathcal{T}$ and set of
matrices $\mathrm{L}(\mathcal{Y}\otimes\mathcal{X})$, known as the Choi-
Jamiołkowski isomorphism [11, 12]. Its explicit form is
${T}=\sum_{i,j}\mathcal{T}(\outerproduct{i}{j})\otimes\outerproduct{i}{j}$. We
will denote linear mappings with calligraphic font
$\mathcal{L},\SS,\mathcal{T}$ etc., whereas the corresponding Choi-
Jamiołkowski matrices as plain symbols: $L,S,T$ etc.

A general quantum measurement (POVM) $\mathcal{Q}$ can be viewed as a set of
positive semidefinite operators $\mathcal{Q}=\\{Q_{i}\\}_{i}$ such that
$\sum_{i}Q_{i}={\rm 1l}$. These operators are usually called effects. The von
Neumman measurements, $\mathcal{P}_{U}$, are a special subclass of
measurements whose all effects are rank-one projections given by
$\mathcal{P}_{U}=\\{P_{U,i}\\}_{i}=\\{U\outerproduct{i}{i}U^{\dagger}\\}_{i}$
for some unitary matrix $U\in\mathrm{L}(\mathcal{X})$. The Choi matrix of
$\mathcal{P}_{U}$ is
$P_{U}=\sum_{i}\outerproduct{i}{i}\otimes\overline{P_{U,i}},$ which will be
utilized throughout this work.

Let us consider a composition of mappings
$\mathcal{R}=\mathcal{N}\circ\mathcal{M}$, where
$\mathcal{N}:\mathrm{L}(\mathcal{Z})\rightarrow\mathrm{L}(\mathcal{Y})$ and
$\mathcal{M}:\mathrm{L}(\mathcal{X})\rightarrow\mathrm{L}(\mathcal{Z})$ with
Choi matrices $N\in\mathrm{L}(\mathcal{Z}\otimes\mathcal{Y})$ and
$M\in\mathrm{L}(\mathcal{X}\otimes\mathcal{Z})$, respectively. Then, the Choi
matrix of $\mathcal{R}$ is given by [13] $R=\tr_{\mathcal{Z}}\left[\left({\rm
1l}_{\mathcal{Y}}\otimes M^{T_{\mathcal{Z}}}\right)\left(N\otimes{\rm
1l}_{\mathcal{X}}\right)\right],$ where $M^{T_{\mathcal{Z}}}$ denotes the
partial transposition of $M$ on the subspace $\mathcal{Z}$. The above result
can be expressed by introducing the notation of the link product of the
operators $N$ and $M$ as

$N*M\coloneqq\tr_{\mathcal{Z}}\left[\left({\rm 1l}_{\mathcal{Y}}\otimes
M^{T_{\mathcal{Z}}}\right)\left(N\otimes{\rm 1l}_{\mathcal{X}}\right)\right].$
(1)

Finally, we define the operator $\prescript{}{\mathcal{X}}{Y}$ as
$\prescript{}{\mathcal{X}}{Y}=\frac{{\rm
1l}_{\mathcal{X}}}{\dim(\mathcal{X})}\otimes\tr_{\mathcal{X}}Y,$ for every
$Y\in\mathrm{L}(\mathcal{X}\otimes\mathcal{Z})$, where $\mathcal{Z}$ is an
arbitrary complex Euclidean space and the projective operator

$L_{V}(W)=\prescript{}{\left[1-\prod_{i}\left(1-\mathcal{A}_{O}^{i}+\mathcal{A}_{I}^{i}\mathcal{A}_{O}^{i}\right)+\prod_{i}\mathcal{A}_{I}^{i}\mathcal{A}_{O}^{i}\right]}{W}.$
(2)

We say that
$W\in\mathrm{Pos}(\mathcal{A}_{I}^{1}\otimes\mathcal{A}_{O}^{1}\otimes\ldots\otimes\mathcal{A}_{I}^{N}\otimes\mathcal{A}_{O}^{N})$
is $N$-partite process matrix if it fulfills the following conditions [10]

$W=L_{V}(W)\,\,\,\text{and}\,\,\tr(W)=\dim(\mathcal{A}_{O}^{1})\cdot\ldots\cdot\dim(\mathcal{A}_{O}^{N}),$
(3)

where the projection operator $L_{V}$ is defined by Eq. (2).

###  II.2 Learning setup

Imagine we are given a black box with the promise that it contains some von
Neumann measurement, $\mathcal{P}_{U}$, which is parameterized by a unitary
matrix $U$. The exact value of $U$ is unknown to us. We are allowed to use the
black box $N$ times. Our goal is to prepare a storage strategy $\SS$ and a
retrieval measurement $\mathcal{R}$ such that we are able to approximate
$\mathcal{P}_{U}$ on an arbitrary state $\rho\in\Omega(\mathcal{X}_{in})$.
This approximation will be denoted throughout this work as $\mathcal{Q}_{U}$.
We would like to point out that, generally, $\mathcal{Q}_{U}$ will _not_ be a
von Neumann measurement. The learning scheme will be denoted by $\mathcal{L}$
with Choi matrix $L$ being a concatenation $L=R*S$. Additionally, we assume
that the Choi matrix of the storage $S$ has an indefinite causal order. More
precisely, the storage $S$ is described by $N-$partite process matrix $W$,
that is $W=\tr_{\mathcal{X}_{a}}S$. We provide an overview of the learning
scheme in Fig. 1.

Figure 1: A schematic representation of the setup for $N\rightarrow 1$
learning scheme of von Neumann measurements $\mathcal{P}_{U}$ with the usage
of indefinite causal order structures.

As a measure of quality of approximating a von Neumann measurement
$\mathcal{P}_{U}=\\{P_{U,i}\\}_{i}$ with a POVM
$\mathcal{Q}_{U}=\\{Q_{U,i}\\}_{i}$ we choose the fidelity function [3], which
is defined as follows

$\mathcal{F}(\mathcal{P}_{U},\mathcal{Q}_{U})\coloneqq\frac{1}{d}\sum_{i}\tr(P_{U,i}Q_{U,i}),$
(4)

where $d$ is the dimension of the measured system. Note that in the case when
$\mathcal{P}_{U}$ is a von Neumann measurement we obtain the value of fidelity
function $\mathcal{F}$ belongs to the interval $[0,1]$ and equals to one if
and only if $P_{U,i}=Q_{U,i}$ for all $i$. As there is no prior information
about $\mathcal{P}_{U}$ provided, we assume that $U$ is sampled from a
distribution pertaining to the Haar measure. Therefore, considering a von
Neumann measurement $\mathcal{P}_{U}$ and its approximation $\mathcal{Q}_{U}$
we introduce the _average fidelity function_ [14] with respect to Haar measure
as

$\mathcal{F}^{\text{
avg}}_{d}\coloneqq\int_{U}dU\mathcal{F}(\mathcal{P}_{U},\mathcal{Q}_{U}).$ (5)

Our main goal is to maximize $\mathcal{F}_{\text{avg}}$ over all possible
learning schemes $\mathcal{L}$. We introduce the notation of the maximum value
of the average fidelity function

$F_{d}(N)\coloneqq\max_{\mathcal{L}}\mathcal{F}_{\text{avg}}.$ (6)

##  III $N\rightarrow 1$ learning scheme of von Neumann measurements

This section presents the various results for $N\rightarrow 1$ learning scheme
of von Neumann measurements. The solution for one copy of von Neumann
measurements was provided in [14]. The Authors then have proved that
$F_{d}(1)=\frac{d+1}{d^{2}}$. But what if we have access to more copies of von
Neumann measurements?

###  III.1  $2\rightarrow 1$ learning scheme

Let us consider a learning scheme $\mathcal{L}$ in which we learn a von
Neumann measurement $\mathcal{P}_{U}$ by using two copies of it. Then, its
Choi operator
$L\in\mathrm{L}(\mathcal{A}_{I}\otimes\mathcal{A}_{O}\otimes\mathcal{B}_{I}\otimes\mathcal{B}_{O}\otimes\mathcal{X}_{in}\otimes\mathcal{X}_{out})$
satisfies the condition

$\tr_{\mathcal{X}_{out}}L={\rm 1l}_{{\mathcal{X}_{in}}}\otimes W,$ (7)

where $W$ is a bipartite process matrix $(N=2)$. The following proposition
says that no matter which approach we use, either parallel or indefinite
causal learning strategies, the maximum value of the average fidelity for
$2\rightarrow 1$ learning scheme of von Neumann measurements is the same. The
solution for the parallel case was provided in [14].

######  Proposition 1.

The usage of indefinite causal order structure does not improve the maximum
value of the average fidelity function in $2\rightarrow 1$ learning scheme of
von Neumann measurements $\mathcal{P}_{U}$ of dimension $d$.

###### Proof.

Due to the fact that the quantum network $\mathcal{L}$ has classical labels on
spaces $\mathcal{A}_{O},\mathcal{B}_{O}$ and $\mathcal{X}_{out}$ then its Choi
matrix
$L\in\mathrm{Herm}(\mathcal{A}_{I}\otimes\mathcal{A}_{O}\otimes\mathcal{B}_{I}\otimes\mathcal{B}_{O}\otimes\mathcal{X}_{in}\otimes\mathcal{X}_{out})$
has the following form
$L=\sum_{i}\outerproduct{i}{i}_{\mathcal{X}_{out}}\otimes L_{i},$ where
$L_{i}=\sum_{j,k}\outerproduct{j}{j}_{\mathcal{B}_{O}}\otimes\outerproduct{k}{k}_{\mathcal{A}_{O}}\otimes
L_{ijk}.$ Simultaneously, we know that
$\sum_{i}L_{i}=\sum_{i,j,k}\outerproduct{j}{j}_{\mathcal{B}_{O}}\otimes\outerproduct{k}{k}_{\mathcal{A}_{O}}\otimes
L_{ijk}={\rm 1l}_{\mathcal{X}_{in}}\otimes W.$ Without loss of generality we
can assume that $L_{ijk}$ satisfies the commutative relation (see Lemma 9.3 in
[14])

$\left[L_{ijk},U_{\mathcal{A}}\otimes U_{\mathcal{B}}\otimes
U^{\dagger}\right]=0,$ (8)

where $U_{\mathcal{A}}\in\mathrm{L}(\mathcal{A}_{I}\otimes\mathcal{A}_{O})$,
$U_{\mathcal{B}}\in\mathrm{L}(\mathcal{B}_{I}\otimes\mathcal{B}_{O})$ and
$U^{\dagger}\in\mathrm{L}\left({\mathcal{X}_{in}\otimes\mathcal{X}_{out}}\right)$.
So, we have $\sum_{i}L_{ijk}={\rm
1l}_{\mathcal{X}_{in}}\otimes\bra{j}_{\mathcal{B}_{O}}\bra{k}_{\mathcal{A}_{O}}W\ket{j}_{\mathcal{B}_{O}}\ket{k}_{\mathcal{A}_{O}}.$
for all $j,k$. From relabeling symmetry property (see Lemma 9.4 in [14]) given
by $L_{ijk}=L_{\sigma(i)\sigma(j)\sigma(k)}$, we have

$\begin{split}\sum_{i}L_{ijk}&=\sum_{i}L_{\sigma(i)\sigma(j)\sigma(k)}\\\
&={\rm
1l}_{\mathcal{X}_{in}}\otimes\bra{\sigma(j)}\bra{\sigma(k)}W\ket{\sigma(j)}\ket{\sigma(k)},\end{split}$
(9)

for any permutation $\sigma$. Therefore, we have

$\begin{split}W&=\sum_{j,k}\outerproduct{j}{j}_{\mathcal{B}_{O}}\otimes\outerproduct{k}{k}_{\mathcal{A}_{O}}\otimes\frac{1}{d}\sum_{i}\tr_{\mathcal{X}_{in}}L_{ijk}\\\
&=\sum_{j,k}\outerproduct{j}{j}_{\mathcal{B}_{O}}\otimes\outerproduct{k}{k}_{\mathcal{A}_{O}}\otimes
W_{jk}.\end{split}$ (10)

Hence, $\forall j,k,\sigma\,\,\ W_{jk}=W_{\sigma(j)\sigma(k)}.$ It implies
that $W_{11}=\ldots=W_{dd}$ and $W_{12}=W_{ab}$ for all $a\neq b$. These
properties together with Eq. (10) imply that $W={\rm
1l}_{\mathcal{B}_{O}}\otimes{\rm 1l}_{\mathcal{A}_{O}}\otimes
P+J(\Delta)\otimes\left(Q-P\right),$ where $P=W_{12}$ and $Q=W_{11}$ and
$J(\Delta)$ denotes the Choi matrix of of the completely dephasing channel
$\Delta$. From the definition of the process matrix (more precisely from the
condition
$W+\prescript{}{\mathcal{A}_{O}\mathcal{B}_{O}}{W}=\prescript{}{\mathcal{A}_{O}}{W}+\prescript{}{\mathcal{B}_{O}}{W}$)
we obtain that $P=Q$, which completes the proof. ∎

###  III.2 Semidefinite program for calculating the maximum value of the
average fidelity

In the general approach, to compute the maximum value of the average fidelity
$F_{d}(N)$ we use the semidefinite programming (SDP). We will present the
original primal problem (Program 1) for computing $F_{d}(N)$ for $N\rightarrow
1$ learning scheme of von Neumann measurement $\mathcal{P}_{U}$ of dimension
$d$. Next, we will describe a simplified version of the primal problem
presented in Program 2.

To optimize this problem, we used the Julia programming language along with
quantum package QuantumInformation.jl [15] and SDP optimization via SCS solver
[16, 17] with absolute convergence tolerance $10^{-5}$. The code is available
on GitHub [18].

Original problem

$\begin{split}\text{maximize:}\quad&\int_{U}dU\frac{1}{d}\sum_{i=1}^{d}\tr\left[L_{i}^{\top}\left(P_{U,i}\otimes
P_{U}^{\otimes N}\right)\right]\\\\[5.69054pt] \text{subject
to:}\quad&L\in\mathrm{Pos}\left(\bigotimes_{i=1}^{N}\mathcal{A}_{I,O}^{i}\otimes\mathcal{X}_{in}\otimes\mathcal{X}_{out}\right),\\\
&L=\sum_{i=1}^{d}\outerproduct{i}{i}_{\mathcal{X}_{out}}\otimes L_{i},\\\
&\prescript{}{\mathcal{X}_{out}}{L}=\prescript{}{\mathcal{X}_{out},\mathcal{X}_{in}}{L},\\\
&\prescript{}{\mathcal{X}_{out}}{L}=\prescript{}{[{\rm 1l}-\Pi_{i}({\rm
1l}-\mathcal{A}_{O}^{i}+\mathcal{A}_{I,O}^{i})+\Pi_{i}\mathcal{A}_{I,O}^{i}]}{\left[\prescript{}{\mathcal{X}_{out}}{L}\right]},\\\
&\tr(L)=d^{N+1}.\end{split}$

Table 1: Semidefinite program for maximizing the value of the average fidelity
function $F$ for $N\rightarrow 1$ learning scheme of von Neumann measurement
$\mathcal{P}_{U}$ of dimension $d$.

Here, we present a simplified description of the primal problem associated
with Program 1. Let
$\mathcal{Y}=\sum_{j\in\mathcal{X}}\left(\sum_{i=1}^{d}\sum_{j\not\in\mathcal{X}}\tr_{\mathcal{X}_{in},\Pi_{\not\mathcal{X}}\mathcal{A}_{I}^{i}}L_{i,j}\right)\otimes\outerproduct{j}{j}_{\mathcal{X}}$
such that $\mathcal{X}\neq\emptyset$. Then, we have:

Simplified problem

$\begin{split}\text{maximize:}\quad&\frac{1}{d}\sum_{i=1}^{d}\sum_{j_{1},\ldots,j_{N}=1}^{d}\tr\left[L_{i,j}\outerproduct{i}{i}\otimes\outerproduct{j}{j}\right]\\\\[5.69054pt]
\text{subject
to:}\quad&L_{i,j}\in\mathrm{Pos}\left(\mathcal{X}_{in}\otimes\mathcal{A}_{I}\right),\\\
&\sum_{i}L_{i,j}=\prescript{}{\mathcal{X}_{in}}{\left[\sum_{i}L_{i,j}\right]}\quad\forall_{j},\\\
&[L_{i,j},\bar{U}\otimes U^{\otimes N}]=0,\\\
&\prescript{}{[\Pi_{\mathcal{X}}({\rm
1l}-\mathcal{A}_{O}^{i})]}{\left[\mathcal{Y}\right]}=0,\\\
&\sum_{i,j}\tr(L_{i,j})=d^{N+1}.\\\ &\end{split}$

Table 2: A simplified description of the SDP Program 1.

######  Remark 1.

We would like to point out that the commutation relation
$[L_{i,j},\bar{U}\otimes U^{\otimes N}]=0$ can be equivalently exchanged with
$L_{i,j}=\bigoplus_{\mu\in\mathrm{irrepS(U\otimes U^{\otimes N})}}{\rm
1l}_{d_{\mu}}\otimes Q_{i,j,\mu}$, where the summation goes over the
irreducible representation of $L_{i,j}$ [14].

###  III.3  $N\rightarrow 1$ learning scheme

In this section, we analyze the asymptotic behavior of $F_{d}(N)$ for
$N\rightarrow\infty$. Our main result can be summarized as the following
theorem.

######  Theorem 1.

Let $F_{d}(N)$ be the maximum value of the average fidelity function, defined
in Eq. (6) for $N\rightarrow 1$ learning scheme of von Neumann measurements.
Then, for arbitrary but fixed dimension $d$ we obtain

$F_{d}(N)=1-\Theta\left(\frac{1}{N^{2}}\right).$ (11)

###### Proof.

We will follow the approach from [2, Lemma 3]. Based on the results from [2]
it is enough to show $F_{d}(N)\leq 1-\Theta\left(\frac{1}{N^{2}}\right).$

A learning network $L$ can be described as a concatenation of a storage $S$
and a retrieval $R$, that is $L=R*S$. What is more, we can assume that storage
is given as a purification, so the Choi-Jamiołkowski isomorphism of $S$ is
pure,
$S=\left|\left.{{X}}\right\rangle\\!\right\rangle\\!\\!\left\langle\\!\left\langle{{X}}\right.\right|\in\mathcal{A}_{I,O}\otimes\mathcal{X}_{a}$
and $X\geq 0$ [2, Lemma 3], see Fig. 1. Then,
$W=\tr_{\mathcal{X}_{a}}(S)=X^{2}.$ From the SDP program, we have $[W,{\rm
1l}_{\mathcal{A}_{O}}\otimes U^{\otimes N}]=0$ for each unitary matrix $U$.
Therefore, $[X,{\rm 1l}_{\mathcal{A}_{O}}\otimes U^{\otimes N}]=0$ and the
memory state that keeps the information of $P_{U}$ has the form

$\begin{split}&\tr_{\mathcal{A}_{I,O}}\left(S\left({\rm
1l}_{{\mathcal{X}_{a}}}\otimes\left({\rm 1l}_{\mathcal{A}_{O}}\otimes
U^{\otimes N}J(\Delta){\rm 1l}_{\mathcal{A}_{O}}\otimes{U^{\dagger}}^{\otimes
N}\right)\right)\right)=\\\ &\left({\rm
1l}_{\mathcal{A}_{O}}\otimes\overline{U}^{\otimes N}\right)\rho\left({\rm
1l}_{\mathcal{A}_{O}}\otimes{U^{\top}}^{\otimes N}\right),\end{split}$ (12)

where $\rho$ is some state. That means, we can upper bound the value of the
fidelity within the new scheme, where we are given in parallel $N$ copies of
unitary channel $\Phi_{\overline{U}}$ and we try to learn $\mathcal{P}_{U}$.
According, to [2, Lemma 7] we get $F_{d}(N)\leq
1-\Theta\left(\frac{1}{N^{2}}\right).$ ∎

##  IV Numerical results for qubit von Neumann measurements

Although the maximum value of the average fidelity function behaves
asymptotically the same using quantum combs or indefinite causal order, we
will show here a numerical advantage of SAR for qubit von Neumann measurements
($d=2$). To show that, we compare the results for learning scheme for $N\geq
3$ with the parallel and adaptive learning schemes introduced in [2].

In the qubit case, we make two simplifications of SDP Program 2. First, the
relation $\forall_{U}\,\,[L_{i,j},\bar{U}\otimes U^{\otimes N}]=0$ is
equivalent with $\forall_{U}\,\,[L_{i,j},U^{\otimes N+1}]=0$. Second, as
$L=\sum_{i,j}\outerproduct{i}{i}_{\mathcal{X}_{out}}\otimes\outerproduct{j}{j}_{\mathcal{A}_{O}}\otimes
L_{i,j}$, then
$W=\sum_{j}\outerproduct{j}{j}_{\mathcal{A}_{O}}\otimes\frac{1}{2}\tr_{\mathcal{X}_{in}}\left(\sum_{i}L_{i,j}\right)$
is a $N$-partite block-diagonal process matrix. Describing $W$ in the Pauli
basis [6] $\left\\{{\rm
1l}_{2},\sigma_{X},\sigma_{Y},\sigma_{Z}\right\\}^{\otimes 2N}$, we get that
$W$ belongs to the subspace spanned by $\left\\{{\rm
1l}_{2},\sigma_{Z}\right\\}^{\otimes N}\bigotimes\left\\{{\rm
1l}_{2},\sigma_{X},\sigma_{Y},\sigma_{Z}\right\\}^{\otimes N}$ and within this
subspace it is orthogonal to

$\mathrm{R}=\\{\sigma_{Z}^{k_{1}}\otimes\cdots\otimes\sigma_{Z}^{k_{N}}\otimes
M_{1}^{k_{1}}\otimes\cdots\otimes M_{N}^{k_{N}}:$ $0\neq
k=(k_{1},\ldots,k_{N})\in\\{0,1\\}^{N},M_{n}^{k}\in\\{{\rm
1l}_{2},\sigma_{X},\sigma_{Y},\sigma_{Z}\\}\\}.$

Verifying if $W$ is orthogonal to $\mathrm{R}$ is easier than verifying
$\prescript{}{[\Pi_{\mathcal{X}}({\rm
1l}-\mathcal{A}_{O}^{i})]}{\left[\mathcal{Y}\right]}=0$. It can be done simply
by checking if

$\sum_{j}(-1)^{k\cdot
j}\tr_{\mathcal{X}_{in},\Pi_{l:k_{l}=0}\mathcal{A}_{I}^{l}}\left(\sum_{i}L_{i,j}\right)=0,$

for each $0\neq k\in\\{0,1\\}^{N}$.

Below we present numerical results for $d=2$ and $N=1,\ldots,5$ that compares
different storing strategies ($\mathcal{L}_{\text{Causal}}$ \- indefinite
causal order learning strategies, $\mathcal{L}_{\text{Adaptive}}$ \- adaptive
learning strategies, $\mathcal{L}_{\text{Parallel}}$ \- parallel learning
strategies) and show the advantage of using indefinite causal order
strategies.

$\begin{array}[]{l|c|c|c|c|c}N&1&2&3&4&5\\\
\hline\cr\mathcal{F}_{2}^{\text{avg}}(\mathcal{L}_{\text{Causal}})&0.7500&0.8114&0.8698&0.8981&0.9204\\\
\hline\cr\mathcal{F}_{2}^{\text{avg}}(\mathcal{L}_{\text{Adaptive}})&0.7499&0.8114&0.8684&0.8968&0.9189\\\
\hline\cr\mathcal{F}_{2}^{\text{avg}}(\mathcal{L}_{\text{Parallel}})&0.7499&0.8114&0.8676&0.8955&0.9187\\\
\hline\cr\end{array}$

The above results are presented also in Figure 2.

Figure 2: The average fidelity function $\mathcal{F}_{2}^{\text{avg}}$ for
$N\rightarrow 1$ learning scheme of qubit von Neumann measurements, where
$N=3,4,5$ – optimal indefinite causal order learning strategy
$\mathcal{L}_{\text{Causal}}$ (orange circles); optimal adaptive learning
strategy $\mathcal{L}_{\text{Adaptive}}$ (purple squares); optimal parallel
learning scheme $\mathcal{L}_{\text{Parallel}}$ (green triangles).

##  V Conclusions and discussion

In this work, we studied the problem of learning an unknown von Neumann
measurement from a finite number of copies $N$ using indefinite causal order
structures. To do so, we have introduced a notion of $N$-partite process
matrices. Our main goal was to compute the maximum value of the average
fidelity function $F_{d}(N)$ of the approximation, having access to $N$ copies
of the given measurement. We have proved that
$F_{d}(N)=1-\Theta\left(\frac{1}{N^{2}}\right)$ for arbitrary but fixed
dimension $d$. Next, we have considered various learning schemes for different
numbers of accessible copies of von Neumann measurements. For $N=2$, we proved
that using an indefinite causal structures do not improve the average fidelity
function $F_{d}(2)$. Next, however, we show numerically advantage of using
indefinite causal order structures for $d=2$ and $N\geq 3$. For this purpose,
we have stated a SDP program and provided its simplified version.

Our results give additional confirmation of potential benefits of using
indefinite causal order structures in quantum information theory. Previously,
applications were observed for example in quantum channel discrimination [6],
quantum communication [19] or unitary channels transformation [20, 21]. Here,
we showed the usage of the theory of indefinite causal order in the problem of
learning of quantum channels, in particular by using storage and retrieval
scheme.

## Acknowledgments

PL is supported by the Ministry of Education, Youth and Sports of the Czech
Republic through the e-INFRA CZ (ID:90254), with the financial support of the
European Union under the REFRESH – Research Excellence For
REgionSustainability and High-tech Industries project number
CZ.10.03.01/00/22_003/0000048 via the Operational Programme Just Transition.

RK is supported by the National Science Centre, Poland, under the contract
number 2021/03/Y/ST2/00193 within the QuantERA II Programme that has received
funding from the European Union’s Horizon 2020 research and innovation
programme under Grant Agreement No 101017733.

## References

  * [1] A. Bisio, G. Chiribella, G. M. D’Ariano, S. Facchini, and P. Perinotti, “Optimal quantum learning of a unitary transformation,” Physical Review A, vol. 81, no. 3, p. 032324, 2010.
  * [2] P. Lewandowska, R. Kukulski, Ł. Pawela, and Z. Puchała, “Storage and retrieval of von neumann measurements,” Physical Review A, vol. 106, no. 5, p. 052423, 2022.
  * [3] M. Raginsky, “A fidelity measure for quantum channels,” Physics Letters A, vol. 290, no. 1-2, pp. 11–18, 2001.
  * [4] V. P. Belavkin, G. M. D’Ariano, and M. Raginsky, “Operational distance and fidelity for quantum channels,” Journal of Mathematical Physics, vol. 46, no. 6, p. 062106, 2005.
  * [5] A. Bisio, G. M. D’Ariano, P. Perinotti, and M. Sedlák, “Quantum learning algorithms for quantum measurements,” Physics Letters A, vol. 375, no. 39, pp. 3425–3434, 2011.
  * [6] J. Bavaresco, M. Murao, and M. T. Quintino, “Strict hierarchy between parallel, sequential, and indefinite-causal-order strategies for channel discrimination,” Physical Review Letters, vol. 127, no. 20, p. 200504, 2021\.
  * [7] M. T. Quintino and D. Ebler, “Deterministic transformations between unitary operations: Exponential advantage with adaptive quantum circuits and the power of indefinite causality,” Quantum, vol. 6, p. 679, 2022.
  * [8] J. Bavaresco, M. Murao, and M. T. Quintino, “Unitary channel discrimination beyond group structures: Advantages of sequential and indefinite-causal-order strategies,” Journal of Mathematical Physics, vol. 63, no. 4, p. 042203, 2022.
  * [9] P. Lewandowska, Ł. Pawela, and Z. Puchała, “Strategies for single-shot discrimination of process matrices,” Scientific Reports, vol. 13, no. 1, p. 3046, 2023.
  * [10] M. Araújo, C. Branciard, F. Costa, A. Feix, C. Giarmatzi, and Č. Brukner, “Witnessing causal nonseparability,” New Journal of Physics, vol. 17, no. 10, p. 102001, 2015.
  * [11] M.-D. Choi, “Completely positive linear maps on complex matrices,” Linear Algebra and its Applications, vol. 10, no. 3, pp. 285–290, 1975.
  * [12] A. Jamiołkowski, “Linear transformations which preserve trace and positive semidefiniteness of operators,” Reports on Mathematical Physics, vol. 3, no. 4, pp. 275–278, 1972.
  * [13] G. Chiribella, G. M. D’Ariano, and P. Perinotti, “Theoretical framework for quantum networks,” Physical Review A, vol. 80, no. 2, p. 022339, 2009.
  * [14] A. Bisio, G. Chiribella, G. D’Ariano, and P. Perinotti, “Quantum networks: general theory and applications,” Acta Physica Slovaca, vol. 61, no. 3, pp. 273–390, 2011.
  * [15] P. Gawron, D. Kurzyk, and Ł. Pawela, “QuantumInformation.jl—a julia package for numerical computation in quantum information theory,” PLOS ONE, vol. 13, p. e0209358, dec 2018.
  * [16] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, “Conic optimization via operator splitting and homogeneous self-dual embedding,” Journal of Optimization Theory and Applications, vol. 169, pp. 1042–1068, June 2016.
  * [17] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, “SCS: Splitting conic solver, version 3.2.1.” https://github.com/cvxgrp/scs, Nov. 2021.
  * [18] https://https://github.com/rkukulski/
sar-measurements-causal.  Permanent link to code/repository, Accessed:
2024-05-13.

  * [19] D. Ebler, S. Salek, and G. Chiribella, “Enhanced communication with the assistance of indefinite causal order,” Physical review letters, vol. 120, no. 12, p. 120502, 2018.
  * [20] M. T. Quintino, Q. Dong, A. Shimbo, A. Soeda, and M. Murao, “Probabilistic exact universal quantum circuits for transforming unitary operations,” Physical Review A, vol. 100, no. 6, p. 062339, 2019.
  * [21] M. T. Quintino, Q. Dong, A. Shimbo, A. Soeda, and M. Murao, “Reversing unknown quantum transformations: Universal quantum circuit for inverting general unitary operations,” Physical Review Letters, vol. 123, no. 21, p. 210502, 2019.


# Integer Linear Programming Formulations for Triple and Quadruple Roman
Domination Problems.

Sanath Kumar Vengaldas<EMAIL_ADDRESS>Computer Science and
Engineering, National Institute of Technology, Warangal
Telangana, India
Adarsh Reddy Muthyala<EMAIL_ADDRESS>Computer Science and
Engineering, National Institute of Technology, Warangal
Telangana, India
Bharath Chaitanya Konkati<EMAIL_ADDRESS>Computer Science
and Engineering, National Institute of Technology, Warangal
Telangana, India
P. Venkata Subba Reddy
<EMAIL_ADDRESS>Computer Science and Engineering, National Institute of
Technology, Warangal
Telangana, India

(13 March 2023)

###### Abstract

Roman domination is a well researched topic in graph theory. Recently two new
variants of Roman domination, namely triple Roman domination and quadruple
Roman domination problems have been introduced, to provide better defense
strategies. However, triple Roman domination and quadruple Roman domination
problems are NP-hard. In this paper, we have provided genetic algorithm for
solving triple and quadruple Roman domination problems. Programming (ILP)
formulations for triple Roman domination and quadruple Roman domination
problems have been proposed. The proposed models are implemented using IBM
CPLEX 22.1 optimization solvers and obtained results for random graphs
generated using NetworkX Erdős-Rényi model.

###### keywords:

Dominating Set; Roman Domination Number; Triple Roman Domination; Quadruple
Roman Domination; Integer Linear Programming.

Mathematics Subject Classification: 05C69, 68Q25

##  1 Introduction

Graphs $G(V,E)$ considered are undirected, connected and simple. Let $V(G)$
represent the vertex set and $E(G)$ represent the edge set of $G$. $N(v)$
represents the open neighbourhood of vertex $v$, $N(v)=\\{u\mid(u,v)\in E\\}$
and $N(v)\cup\\{v\\}$ represents the closed neighbourhood $N[v]$ of vertex
$v$.

A function $f:V\rightarrow\\{0,1\\}$ on a graph $G$ satisfying the condition
that for every vertex $v\in V$ if $f(v)=0$ then $v$ has a neighbour $u$ such
that $f(u)=1$ is called a dominating function of $G$. Domination number
($\gamma(G)$) is the minimum possible sum of $f$ values assigned to all
vertices of $G$. Optimization version of domination problem is NP-hard and
decision version is NP-complete [7].

Based on the defense tactics to defend the Roman empire, the concept of Roman
domination was introduced in 2004 [5],[12],[13]. A legion from a powerful
neighbour might repel an unexpected invasion on any neighbouring undefended
city and if such moving would leave its existing city defenceless, no legion
could relocate. A function $f:V\rightarrow\\{0,1,2\\}$ on $G$ is called a
Roman dominating function if for every vertex $v\in V$, either
$f(v)\in\\{1,2\\}$ or if $f(v)=0$ then it should have at least 1 neighbour $u$
such that $f(u)=2$. Roman domination number ($\gamma_{R}(G)$) is the minimum
possible sum of $f$ values assigned to all vertices of $G$. Optimization
version of Roman domination problem is NP-hard [11].

Double Roman domination was introduced as a stronger version of weak Roman
domination and Roman domination in [2]. A function
$f:V\rightarrow\\{0,1,2,3\\}$ on a graph $G$ is called a double Roman
dominating function if following conditions are satisfied, for every vertex
$v\in V$, if $f(v)=0$, then vertex $v$ must have at least one neighbour $u$
with $f(u)=3$ or at least 2 neighbours assigned value 2 under $f$ , and if
$f(v)=1$, then $v$must have at least 1 neighbour $w$ with $f(w)\geq 2$. A
Double Roman domination number $\gamma_{dR}(G)$ is the minimum possible sum of
$f$ values of all vertices of $G$.

In [1], triple Roman domination problem is introduced because double Roman
domination problem is not efficient in some cases as discussed in [2].A
function $f:V\rightarrow\\{0,1,2,3,4\\}$ on a graph $G=(V,E)$ is called a
triple Roman dominating function $(3RDF)$ if the following conditions are
satisfied,for every vertex $v\in V$, if $f(v)=0$ then $v$ must have at least 3
neighbours assigned value 2 under $f$ or two neighbours $u$,$w$ with $f(u)\geq
2$ and $f(w)=3$ or one neighbour $z$ with $f(z)=4$, and if $f(v)=1$ then
vertex $v$ must have at least 2 neighbours assigned value 2 under $f$ or at
least 1 neighbour $w$ such that $f(w)\geq 3$ , and if $f(v)=2$, then vertex
$v$ must have at least 1 neighbour $w$ with $f(w)\geq 2$. Weight of a 3RDF is
the sum of function $f$ values of all the vertices. Triple Roman domination
number of a graph $G$ ($\gamma_{3R}(G)$) is minimum weight of all the possible
triple Roman dominating functions. Determining $\gamma_{3R}(G)$ for a graph
$G$ is known as triple Roman domination problem (M3RDP). M3RDP for bipartite
graphs and chordal graphs, triple Roman domination problem is NP-hard [1].

A function $f:V\rightarrow\\{0,1,2,3,4,5\\}$ on a graph $G=(V,E)$ is called
quadruple Roman dominating function (4RDF) if the following conditions are
satisfied, for every vertex $v\in V$, if $f(v)=0$ then $v$ must have at least
1 neighbour vertex $u$ with $f(u)=5$ or two vertices $u,w$ such that $f(u)=4$
and $f(w)\geq 2$ or at least 2 vertices such that $f(u)\geq 3$ or at least 4
vertices such that $f(u)=2$ or two vertices with $f(u)=2$ and one vertex with
$f(u)=3$, and if $f(v)=1$ then vertex $v$ must have at least 1 neighbour $u$
with $f(u)\geq 4$ or at least 2 vertices $u,w$ such that $f(u)=3$ and
$f(w)\geq 2$ or at least 3 vertices with $f(u)=2$, and if $f(v)=2$ then vertex
$v$ must have at least 1 neighbour $u$ with $f(u)\geq 3$ or at least 2
vertices with $f(u)=2$, and if $f(v)=3$ then it must have at least 1 neighbour
$u$ with $f(u)\geq 2$. Weight of a 4RDF is the sum of function $f$ values of
all the vertices. Quadruple Roman domination number of a graph $G$
($\gamma_{4R}(G)$) is the minimum weight of all the possible 4RDFs.
Determining $\gamma_{4R}(G)$ for a graph $G$ is known as quadruple Roman
domination problem (M4RDP). Throughout the paper, by $\gamma_{3R}$-function we
mean any 3RDF of graph $G$ with weight $\gamma_{3R}(G)$. Similarly
$\gamma_{4R}$-function is defined.

In general $[k]$-Roman dominating function ($k$RDF) of a graph $G$ is also
discussed in [1] and it is a mapping $f:V\rightarrow\\{0,1,2,3,\ldots,k+1\\}$
such that for every vertex $u\in V\ for\ which\ f(u)\leq k$ following
condition should satisfy.

${h(AN[u])\geq\;\lvert AN(u)\rvert+k}$

Here, $AN(u)=\\{z:z\in N(u)\;\text{and}\;f(v)\neq 0\\},$
$AN[u]=AN(u)\cap\\{u\\}$ and $\lvert AN(u)\rvert$ represents is the
cardinality of $AN(u)$. The weight of a $k$RDF is $\sum\limits_{{u\in
V}}f(u)$. The $[k]$-Roman domination number $\gamma_{[kR]}(G)$ is the minimum
weight of all the $k$RDFs on $G$.
Motivated by the Integer Linear Programming (ILP) formulations given for
different domination problems [4] [8] [9] [10], in this paper we present
different ILP formulations for M3RDP and M4RDP problems. The proposed models
are implemented using IBM CPLEX 22.1 optimization solvers and obtained results
for random graphs generated using NetworkX Erdős-Rényi model, a widely used
model for generating random graphs [3] [6].

##  2 ILP Formulations for Triple Roman Domination Problem

Proposed ILP formulations for M3RDP are discussed in this section.

###  2.1 M3RDP-1

This is very basic ILP formulation for M3RDP where we assigned four decision
variables ($p_{v},q_{v},r_{v},s_{v}$) for each vertex $v$ and their
definitions are as follows

$p_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 1}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$q_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 2}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$r_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 3}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$s_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 4}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$t_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{\begin{math}\sum\limits_{{u\in
N(v)}}q_{u}\geq 1\end{math}}\\\ 0,&\mbox{otherwise}.\end{array}\right.$
$x_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{\begin{math}\sum\limits_{{u\in
N(v)}}r_{u}\geq 1\end{math}}\\\ 0,&\mbox{otherwise}.\end{array}\right.$

An ILP formulation for M3RDP is given below:

$Minimize\hskip 28.45274pt\sum\limits_{{v\in V}}p_{v}+\sum\limits_{{v\in
V}}2q_{v}+\sum\limits_{{v\in V}}3r_{v}+\sum\limits_{{v\in V}}4s_{v}$ (1a) Such
that $p_{v}+q_{v}+r_{v}+s_{v}+\sum_{u\in N(v)}s_{u}+\frac{1}{3}\sum_{u\in
N(v)}q_{u}+\frac{1}{2}\left(t_{v}+x_{v}\right)\geq 1$ (1b)

$\frac{1}{2}\sum_{u\in N(v)}q_{u}+\sum_{u\in N(v)}\left(s_{u}+r_{u}\right)\geq
p_{v}$ (1c) $\sum_{u\in N(v)}\left(q_{u}+r_{u}+s_{u}\right)\geq q_{v}$ (1d)
$p_{v}+q_{v}+r_{v}+s_{v}\leq 1$ (1e)
$p_{v},q_{v},r_{v},s_{v},t_{v},x_{v}\in\\{0,1\\}$ (1f)

As we can see number of decision variables are 4$\lvert V\rvert$
($p_{v},q_{v},r_{v},s_{v}$) and number of constraints are 5$\lvert V\rvert$.

Constraint (1a) gives the objective function value which is also the value of
triple Roman domination number $\gamma_{3R}(G)$. Constraint (1b) is to ensure
that every vertex with $f(v)=0$ should have at least 1 neighbour with $f(v)=4$
or at least 3 neighbours with $f(v)=2$ or at least 2 vertices $v$,$w$ such
that $f(v)$=2 and $f(w)=3$. Constraint (1c) is to ensure that for every vertex
$v$ with $f(v)=1$ there should be at least 2 neighbours with $f(v)=2$ or at
least 1 neighbour with $f(v)\geq 3$. Constraint (1d) ensures that for every
vertex $v$ with $f(v)=2$ there should be at least 1 neighbour with $f(v)\geq
2$. Constraint (1e) is to ensure that every vertex gets only one label from
the set {0,1,2,3,4}. Constraint (1f) ensures that variables
$p_{v},q_{v},r_{v},s_{v}$ does not take any other integer values other than 0
or 1, restricting them to be decision variables.

###  2.2 M3RDP-2

In this section, we propose an ILP formulation for M3RDP with only three
binary variables. Substituting $p_{v}$ = 0 in M3RDP-1 gives us a new model
M3RDP-2. After reducing one variable the formulation will become

$Minimize\hskip 28.45274pt\sum\limits_{{v\in V}}2q_{v}+\sum\limits_{{v\in
V}}3r_{v}+\sum\limits_{{v\in V}}4s_{v}$ (2a) Such that

$q_{v}+r_{v}+s_{v}+\sum_{u\in N(v)}s_{u}+\frac{1}{3}\sum_{u\in
N(v)}q_{u}+\frac{1}{2}\left(t_{v}+x_{v}\right)\geq 1$ (2b) $\sum_{u\in
N(v)}\left(q_{u}+r_{u}+s_{u}\right)\geq q_{v}$ (2c) $q_{v}+r_{v}+s_{v}\leq 1$
(2d) $q_{v},r_{v},s_{v},t_{v},x_{v}\in\\{0,1\\}$ (2e)

Here the total number of decision variables are 3$\lvert V\rvert$
($q_{v},r_{v},s_{v}$) and the number of constraints are 4$\lvert V\rvert$.
Let’s say $\gamma_{3R}^{{}^{\prime}}(G)$ is the weight of triple Roman
domination function we get from M3RDP-2.
Next we show that ILP formulations M3RDP-1 and M3RDP-2 are equivalent except
for the number of variables and constraints.

######  Lemma 2.1.

Let $h$ be a $\gamma_{3R}$-function of $G$. Then there doesn’t exist two
adjacent vertices $v$,$u$ with labels as $h(v)=1$ and $h(u)=4$.

######  Proof 2.2.

By contradiction, let $s$ and $t$ be two adjacent vertices of $G$ with labels
1 and 4 under $h$. Now, we obtain a labelling $h^{\prime}$ of $G$ from $h$ as
follows.

$h^{\prime}(v)=\begin{cases}0,&\text{ if }\;v=s\\\ 4,&\text{ if }\;v=t\\\
h(v),&\text{otherwise}\end{cases}$ (3)

Clearly $h^{\prime}$ is a 3RDF of $G$, because no other label depends on
vertex with label 1. Hence $w(h^{\prime})<w(h)$, a contradiction.

######  Theorem 2.3.

M3RDP-1 is equivalent to M3RDP-2.

######  Proof 2.4.

Clearly $\gamma_{3R}(G)$ is the optimal value of M3RDP-1 for a given graph
$G$. Let $\gamma_{3R}^{\prime}(G)$ is the optimal value of M3RDP-2. Next, we
show that $\gamma_{3R}(G)=\gamma_{3R}^{{}^{\prime}}(G)$. In M3RDP-2 we have
removed the variable $p_{v}$ which represents the vertices with $f(v)=1$. In
fact we don’t need $p_{v}$ that is we can always rearrange labels of vertices
in the closed neighbourhood such the weight of the closed neighbourhood does
not change. Let’s see how we can rearrange.

Let’s say vertex $v$ has $f(v)=1$ then the conditions are such that it should
have at least 2 adjacent vertices whose $f(u)=2$ or at least 1 adjacent vertex
whose $f(u)\geq 3$. Now lets consider the possible scenarios of neighbourhood
of vertex $v$.

1) Vertex $v$ has two such adjacent vertices $u,w$ where $f(u)=2$ and
$f(w)=2$. Now we rearrange labels of $v,u,w$ as below

$f(v)=0,f(u)=3,f(w)=2$

we can create a 3RDF $g$ defined as follows: $g(x)=f(x)$ for all
$x\notin\\{v,u\\}$ and $g(v)=0$, $g(u)=3$. $g$ is a 3RDF with same weight as
$f$.

2) Our second scenario is that vertex $v$ has at least 1 adjacent vertex $u$
such that $f(u)=3$ then we can rearrange labels of $v$ and $u$ as below

$f(v)=0,f(u)=4$

Modified $f$ is also a 3RDF.

3) As discussed in lemma 1, there will be no optimal labelling of adjacent
vertices $v$,$u$ with $f(v)=1$ and $f(u)=4$.

###  2.3 M3RDP-3

Further optimising the above ILP formulation M3RDP-2 we can remove constraint
(2d). So our new optimised formulation will be

$Minimize\hskip 28.45274pt\sum\limits_{{v\in V}}2q_{v}+\sum\limits_{{v\in
V}}3r_{v}+\sum\limits_{{v\in V}}4s_{v}$ (4a)

Such that

$q_{v}+r_{v}+s_{v}+\sum_{u\in N(v)}s_{u}+\frac{1}{3}\sum_{u\in
N(v)}q_{u}+\frac{1}{2}\left(t_{v}+x_{v}\right)\geq 1$ (4b) $\sum_{u\in
N(v)}\left(q_{u}+r_{u}+s_{u}\right)\leq q_{v}$ (4c)
$q_{v},r_{v},s_{v},t_{v},x_{v}\in\\{0,1\\}$ (4d)

Here the number of decision variables are 3$\lvert V\rvert$
($q_{v},r_{v},s_{v}$) and number of constraints are 3$\lvert V\rvert$.

######  Theorem 2.5.

Let $\gamma_{3R}^{{}^{\prime}}(G)$ and $\gamma_{3R}^{{}^{\prime\prime}}(G)$
represent the optimal weight of 3RDF obtained using M3RDP-2 and M3RDP-3. Then
$\gamma_{3R}^{{}^{\prime}}(G)=\gamma_{3R}^{{}^{\prime\prime}}(G)$

######  Proof 2.6.

Since the absence of constraint (2d) in M3RDP-3 allows
$p_{v},q_{v},r_{v},s_{v}$ to take 1 simultaneously i.e for example both
$r_{v},s_{v}$ can be 1 which has no significance while calculating triple
Roman domination number using those decision variables for vertex $v$. If
$S_{2}$ is the set of all possible 3RDFs using M3RDP-2 formulation and
similarly $S_{3}$ is the set of all possible 3RDFs for M3RDP-3 then it is
clear that $S_{2}\subseteq S_{3}$ which is shown in figure 1.
$S_{3}$$S_{2}$

From the above Venn diagram it is clear that any 3RDF given by M3RDP-3 is
either from region $S_{2}$ or from region $S_{32}(=S_{3}-S_{2}$). If we choose
from $S_{2}$ there is no chance that we get optimal value lesser than M3RDP-2
formulation so we need to check from the set $S_{32}$. In $S_{32}$ at least 2
variables from $q_{v},r_{v},s_{v}$ will be equal to 1. Lets assume
$H(v)=2q_{v}+3r_{v}+4s_{v}$. For a vertex $v\in V$ and $f\in S_{32}$, the
maximum value of $H(v)$ will be 2+3+4=9 and minimum value of $H(v)$ will be 5
with M3RDP-3 formulation where as if $f\in S_{2}$ then maximum and minimum
values of $H(v)$ are 4 and 2 respectively. We can see that maximum value
$H(v)$ when $f\in S_{2}$ is less than that of minimum value $H(v)$ of when
$f\in S_{32}$. So definitely the optimal value with function from $S_{32}$ set
will be strictly greater than the optimal value with function from $S_{2}$ set
with M3RDP-3 formulation so we will be choosing function from $S_{2}$ set
without the constraint (2d) also. From Theorem 1, we know that M3RDP-2 gives
optimal answer. Therefore from the fact that M3RDP-3 chooses a function from
set $S_{2}$ it follows that both the models are give the same result.

##  3 ILP Formulations for Quadruple Roman Domination Problem

In the similar way of construction of ILP models for M3RDP we have also
provided ILP models for M4RDP.

###  3.1 M4RDP-1

$p_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 1}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$q_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 2}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$r_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 3}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$s_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 4}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$t_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{$f(v)$ = 5}.\\\
0,&\mbox{otherwise}.\end{array}\right.$
$x_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{\begin{math}\sum\limits_{{u\in
N(v)}}q_{u}\geq 1\end{math}}.\\\ 0,&\mbox{otherwise}.\end{array}\right.$
$y_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{\begin{math}\sum\limits_{{u\in
N(v)}}r_{u}\geq 1\end{math}}.\\\ 0,&\mbox{otherwise}.\end{array}\right.$
$z_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{\begin{math}\sum\limits_{{u\in
N(v)}}s_{u}\geq 1\end{math}}.\\\ 0,&\mbox{otherwise}.\end{array}\right.$
$a_{v}=\left\\{\begin{array}[]{ll}1,&\mbox{\begin{math}\sum\limits_{{u\in
N(v)}}q_{u}\geq 2\end{math}}.\\\ 0,&\mbox{otherwise}.\end{array}\right.$

$Minimize\hskip 28.45274pt\sum\limits_{{v\in V}}p_{v}+2\sum\limits_{{v\in
V}}q_{v}+3\sum\limits_{{v\in V}}r_{v}+4\sum\limits_{{v\in
V}}s_{v}+5\sum\limits_{{v\in V}}t_{v}\\\ $ (5a) Such that
$\begin{split}p_{v}+q_{v}+r_{v}+s_{v}+t_{v}+\sum\limits_{{u\in
N(V)}}t_{u}+\frac{1}{2}(x_{v}+z_{v})+\frac{1}{2}(z_{v}+y_{v})+\frac{1}{2}(a_{v}+y_{v})\\\
+\frac{1}{2}\sum\limits_{{u\in N(v)}}s_{u}+\frac{1}{2}\sum\limits_{{u\in
N(v)}}r_{u}+\frac{1}{4}\sum\limits_{{u\in N(v)}}q_{u}\geq 1\end{split}$ (5b)
$\sum\limits_{{u\in N(v)}}t_{u}+\sum\limits_{{u\in
N(v)}}s_{u}+\frac{1}{2}\left(y_{v}+x_{v}\right)+\frac{1}{2}\sum\limits_{{u\in
N(v)}}r_{u}+\frac{1}{3}\sum\limits_{{u\in N(v)}}q_{u}\geq p_{v}$ (5c)
$\sum\limits_{{u\in
N(v)}}\left(t_{u}+s_{u}+r_{u}\right)+\frac{1}{2}\sum\limits_{{u\in
N(v)}}q_{u}\geq q_{v}$ (5d) $\sum\limits_{{u\in
N(v)}}\left(t_{u}+s_{u}+q_{u}+r_{u}\right)\geq r_{v}$ (5e)
$p_{v}+q_{v}+r_{v}+s_{v}+t_{v}\leq 1$ (5f)
$p_{v},q_{v},r_{v},s_{v},t_{v},x_{v}\in\\{0,1\\}$ (5g)

As we can see number of decision variables are 5$\lvert V\rvert$
($p_{v},q_{v},r_{v},s_{v},t_{v}$) and number of constraints are 6$\lvert
V\rvert$.

Constraint (4a) gives the objective function value which is the quadruple
Roman domination number $\gamma_{4R}(G)$. Constraint (4b) is to ensure that
every vertex with $f(v)=0$ should have at least 1 neighbour with $f(v)=5$ or
at least 2 vertices $v$,$w$ such that $f(v)=4$ and $f(w)\geq 2$ or at least 2
vertices $v$,$w$ with $f(v)=3$ and $f(w)\geq 3$ or at least four vertices with
$f(v)=2$. Constraint (4c) is to ensure that for every vertex $v$ with $f(v)=1$
there should be at least 1 vertex with $f(v)\geq 4$ or at least 2 vertices
$v$,$w$ with $f(v)=3$ and $f(w)\geq 2$. Constraint (4d) ensures that for every
vertex $v$ with $f(v)=2$ there should be at least 1 neighbour with $f(v)\geq
3$ or at least 2 vertices with $f(v)\geq 2$. Constraint (4e) is to ensure that
for every vertex $v$ with $f(v)=3$ there should be at least 1 neighbour with
$f(v)\geq 2$. Constraint (4f) is to ensure that every vertex gets only one
label from the set {0,1,2,3,4,5}. Constraint (4g) ensures that variables
$p_{v},q_{v},r_{v},s_{v},t_{v},x_{v}$ does not take any other integer values
other than 0 or 1, restricting them to be decision variables.

###  3.2 M4RDP-2

In this section, we proposed an ILP formulation of M4RDP with only 4 decision
variables. Substituting $p_{v}$ = 0 in M4RDP-1 gives us a new model M4RDP-2.
After reducing one variable the formulation will become

$Minimize\hskip 14.22636pt2\sum\limits_{{v\in V}}q_{v}+3\sum\limits_{{v\in
V}}r_{v}+4\sum\limits_{{v\in V}}s_{v}+5\sum\limits_{{v\in V}}t_{v}\\\ $ (6a)
Such that $\begin{split}q_{v}+r_{v}+s_{v}+t_{v}+\sum\limits_{{u\in
N(v)}}t_{u}+\frac{1}{2}(x_{v}+z_{v})+\frac{1}{2}(z_{v}+y_{v})+\frac{1}{2}(a_{v}+y_{v})\\\
+\frac{1}{2}\sum\limits_{{u\in N(v)}}s_{u}+\frac{1}{2}\sum\limits_{{u\in
N(v)}}r_{u}+\frac{1}{4}\sum\limits_{{u\in N(v)}}q_{u}\geq 1\end{split}$ (6b)
$\sum\limits_{{u\in
N(v)}}\left(t_{u}+s_{u}+r_{u}\right)+\frac{1}{2}\sum\limits_{{u\in
N(v)}}q_{u}\geq q_{v}$ (6c) $\sum\limits_{{u\in
N(v)}}\left(t_{u}+s_{u}+q_{u}+r_{u}\right)\geq r_{v}$ (6d)
$q_{v}+r_{v}+s_{v}+t_{v}\leq 1$ (6e)
$q_{v},r_{v},s_{v},t_{v},x_{v}\in\\{0,1\\}$ (6f)

######  Lemma 3.1.

For optimal labelling of vertices for quadruple Roman domination function,
there doesn’t exist a labelling where two adjacent vertices $v$,$u$ labels as
$f(v)=1$ and $f(u)=5$.

######  Proof 3.2.

Similar to lemma 1, there will be no such optimal combination since we can
change vertex $v$ label to 0 because no other label depends on vertex with
label 1.Changing label of vertex $v$ from 1 to 0 decreases weight of the
function $f$ and resulting in more optimal value.

######  Theorem 3.3.

M4RDP-1 is equivalent to M4RDP-2.

######  Proof 3.4.

Similar to theorem 1, in M4RDP-2 we have removed the variable $p_{v}$ which
represents the vertices with $f(v)=1$. Rearranging labels can be done as
follows.

Let’s say vertex $v$ has $f(v)=1$ then the conditions are such that it should
have at least 1 adjacent vertex $u$ whose $f(u)=5$ or at least 1 adjacent
vertex $u$ whose $f(u)=4$ or at least one adjacent vertex $u$ with $f(u)=3$
and at least adjacent vertex $w$ with $f(w)=2$ or at least 3 adjacent vertices
with label. $f(u)\geq 3$. Now lets consider the possible scenarios of
neighbourhood of vertex $v$.
1) From lemma 2, there will no optimal labelling of adjacent vertices $v,u$
with $f(v)=1$ and $f(u)=5$.
2) Vertex $v$ has at least adjacent vertex $u$ where $f(u)=4$ . Now we
rearrange labels of $v,u,w$ as below

$f(v)=0,f(u)=5$

we can create a triple Roman domination function $g$ defined as follows:
$g(x)=f(x)$ for all $x\notin\\{v,u\\}$ and $g(v)=0$, $g(u)=5$. $g$ is a 4RDF
with same weight as $f$.
3) Vertex $v$ has at least 1 adjacent vertex $u$ such that $f(u)=3$ and at
least 1 adjacent vertex $w$ such that $f(w)=2$ then we can rearrange labels of
$v$,$u$ and $w$ as below

$f(v)=0,f(u)=3,f(w)=3$

Modified $f$ also 4RDF of $G$.
4) Vertex $v$ has at least 3 adjacent vertex $u,w,y$ such that $f(u)=2$ then
we can rearrange the labels of $v$,$u$,$w$ and $y$ as below

$f(v)=0,f(u)=2,f(w)=2,f(y)=3$

Here number of decision variables are 4$\lvert V\rvert$
($q_{v},r_{v},s_{v},t_{v}$) and number of constraints are 5$\lvert V\rvert$.

###  3.3 M4RDP-3

Further optimising the above ILP formulation M4RDP-2, we can remove constraint
(5e). So our new optimised formulation will be

$Minimize\hskip 14.22636pt2\sum\limits_{{v\in V}}q_{v}+3\sum\limits_{{v\in
V}}r_{v}+4\sum\limits_{{v\in V}}s_{v}+5\sum\limits_{{v\in V}}t_{v}\\\ $ (7a)
Such that $\begin{split}q_{v}+r_{v}+s_{v}+t_{v}+\sum\limits_{{u\in
N(v)}}t_{u}+\frac{1}{2}(x_{v}+z_{v})+\frac{1}{2}(z_{v}+y_{v})+\frac{1}{2}(a_{v}+y_{v})\\\
+\frac{1}{2}\sum\limits_{{u\in N(v)}}s_{u}+\frac{1}{2}\sum\limits_{{u\in
N(v)}}r_{u}+\frac{1}{4}\sum\limits_{{u\in N(v)}}q_{u}\geq 1\end{split}$ (7b)
$\sum\limits_{{u\in
N(v)}}\left(t_{u}+s_{u}+r_{u}\right)+\frac{1}{2}\sum\limits_{{u\in
N(v)}}q_{u}\geq q_{v}$ (7c) $\sum\limits_{{u\in
N(v)}}\left(t_{u}+s_{u}+q_{u}+r_{u}\right)\geq r_{v}$ (7d)
$q_{v},r_{v},s_{v},t_{v},x_{v}\in\\{0,1\\}$ (7e)

Here number of decision variables are 4$|V|$ ($q_{v},r_{v},s_{v},t_{v}$) and
number of constraints are 4 $|V|$. Let’s say
$\gamma_{4R}^{{}^{\prime\prime}}(G)$ is the optimal value of quadruple Roman
domination number we get from M4RDP-3. Similar to theorem 2, we can also prove
that $\gamma_{4R}^{{}^{\prime\prime}}(G)$ is same as $\gamma_{4R}(G)$

##  4 Computational Results

In this section, we showed the results of all the models we introduced in
section 3. All the models are solved using IBM CPLEX 22.1 optimization
solvers. NetworkX, a python package is used to generate random graphs. Whole
computaion is performed on MacOS Monterey version 12.6, M1 chip 2020, 8GB RAM.
Random graphs are generated using NetworkX Erdős–Rényi model. We have 2 tables
representing the running time of our models on the random graphs. Table 1
contains 4 columns instance refers to properties of graph $\lvert V\rvert$
means number of vertices, $\lvert E\rvert$ number of edges, $p$ is the
probability of each edge whether to take into generation of graph or not while
generating random graphs. Second column represent our first model in section 2
and it contains 2 sub columns result and time, result is $\gamma_{3R}(G)$ and
time is the time taken to get solution using M3RDP-1. Similarly next columns
represent our remaining model performances. We can see that all graphs are
giving same results for same graph which represents correctness of our models.
Similarly Table-2 represents our ILP model performances for minimum quadruple
Roman domination problem. We are using ’-’ sign to represent that optimal
solution is not being found in 1800 seconds for that particular model. We can
see that for most of the cases M3RDP-3 model performs better than the other
two models for triple Roman domination problem. For quadruple Roman domination
problem, M4RDP-3 performs better than remaining two.

Table 1: Performances of ILP Model for triple Roman domination problem

Instance | M3RDP-1 | M3RDP-2 | M3RDP-3 |
---|---|---|---|---
$\mid V\mid$ | $\mid E\mid$ | p | result | time | result | time | result | time |
25 | 50 | 0.2 | 26 | 0.44 | 26 | 0.38 | 26 | 0.39 |
| 141 | 0.5 | 11 | 0.45 | 11 | 0.43 | 11 | 0.81 |
| 248 | 0.8 | 7 | 0.37 | 7 | 0.43 | 7 | 0.35 |
50 | 250 | 0.2 | 24 | 14.05 | 24 | 5.82 | 24 | 3.95 |
| 631 | 0.5 | 12 | 174.49 | 12 | 330.18 | 12 | 89.37 |
| 985 | 0.8 | 8 | 3.51 | 8 | 2.04 | 8 | 2.42 |
75 | 556 | 0.2 | - | - | - | - | - | - |
| 1275 | 0.5 | - | - | 12 | 51.85 | 12 | 156.68 |
| 2238 | 0.8 | 8 | 23.79 | 8 | 6.16 | 8 | 11.75 |
100 | 1015 | 0.2 | - | - | - | - | - | - |
| 2502 | 0.5 | - | - | - | - | - | - |
| 3961 | 0.8 | 8 | 17.02 | 8 | 16.72 | 8 | 19.36 |
150 | 2257 | 0.2 | - | - | - | - | - | - |
| 5530 | 0.5 | - | - | - | - | - | - |
| 8973 | 0.8 | 8 | 16.52 | 8 | 14.22 | 8 | 15.45 |

Table 2: Performances of ILP Model for quadruple Roman domination problem

Instance | M4RDP-1 | M4RDP-2 | M4RDP-3
---|---|---|---
$\mid V\mid$ | $\mid E\mid$ | p | result | time | result | time | result | time
10 | 10 | 0.2 | 18 | 1.15 | 18 | 1.42 | 18 | 0.91
| 26 | 0.5 | 11 | 4.38 | 11 | 1.29 | 11 | 0.94
| 31 | 0.8 | 9 | 3.09 | 9 | 0.89 | 9 | 0.81
25 | 56 | 0.2 | - | - | 25 | 760.52 | 25 | 642.75
| 145 | 0.5 | - | - | 14 | 1247.1 | - | -
| 244 | 0.8 | - | - | 9 | 76.41 | 9 | 61.31
50 | 242 | 0.2 | - | - | - | - | - | -
| 601 | 0.5 | - | - | - | - | - | -
| 999 | 0.8 | - | - | 9 | 1416 | - | -

##  5 Conclusion

It is evident from both the tables that M3RDP-3 model is outperforming in most
of the cases for M3RDP and M4RDP-3 model is outperforming for M4RDP. As can be
seen from Table-1, models are becoming infeasible from 100 vertices for M3RDP
where as in Table-2 models are becoming infeasible from 50 vertices. From both
the tables we can conclude that if the number of vertices are increasing,
models are becoming infeasible within time limit. Hence, providing ILP
formulations for M3RDP and M4RDP with lesser number of constraints so that
CPLEX solvers can successfully find solutions of larger graph instances is an
interesting direction for future work

## References

  * [1] H. Abdollahzadeh Ahangar, M.P. Álvarez, M. Chellali, S.M. Sheikholeslami, and J.C. Valenzuela- Tripodoro. Triple roman domination in graphs. Applied Mathematics and Computation, 391:125444, 2021\.
  * [2] Robert A. Beeler, Teresa W. Haynes, and Stephen T. Hedetniemi. Double roman domination. Discrete Applied Mathematics, 211:23–29, 2016.
  * [3] Anthony Bonato and Changping Wang. A note on domination parameters in random graphs. Discus- siones Mathematicae Graph Theory, 28:335–343, 01 2008.
  * [4] Qingqiong Cai, Neng Fan, Yongtang Shi, and Shunyu Yao. Integer linear programming formulations for double roman domination problem. Optimization Methods and Software, 37(1):1–22, 2022.
  * [5] Ernie J Cockayne, Paul A Dreyer, Sandra M Hedetniemi, and Stephen T Hedetniemi. Roman domination in graphs. Discrete Mathematics, 278(1):11–22, 2004.
  * [6] Roman Glebov, Anita Liebenau, and Tibor Szabó. On the concentration of the domination number of the random graph, 2012.
  * [7] Teresa W Haynes, Stephen Hedetniemi, and Peter Slater. Fundamentals of domination in graphs. CRC press, 2013.
  * [8] Marija Ivanovic. A mixed integer linear programming formulation for restrained roman domination problem. Theory and Applications of Mathematics Computer Science TAMCS, 5, 11 2015.
  * [9] Marija Ivanovic. Improved mixed integer linear programing formulations for roman domination problem. Publications de l’Institut Mathematique, 99:51–58, 01 2016.
  * [10] Marija Ivanović. Improved integer linear programming formulation for weak roman domination problem. Soft Comput., 22(19):6583–6593, oct 2018.
  * [11] Mathieu Liedloff, Ton Kloks, J. Q. Liu, and Sheng-Hueng Peng. Roman domination in some special classes of graphs. 2008.
  * [12] Charles S ReVelle and Kenneth E Rosing. Defendens imperium romanum: a classical problem in military strategy. The American Mathematical Monthly, 107(7):585–594, 2000.
  * [13] Ian Stewart. Defend the roman empire! Scientific American, 281(6):136–138, 1999. 10


# Greybody Factors Imprinted on Black Hole Ringdowns:
an alternative to superposed quasi-normal modes

Naritaka Oshita  Center for Gravitational Physics and Quantum Information
(CGPQI), Yukawa Institute for Theoretical Physics (YITP), Kyoto University,
606-8502, Kyoto, Japan  The Hakubi Center for Advanced Research, Kyoto
University, Yoshida Ushinomiyacho, Sakyo-ku, Kyoto 606-8501, Japan  RIKEN
iTHEMS, Wako, Saitama, 351-0198, Japan

###### Abstract

It is shown that the spectral amplitude of gravitational-wave (GW) ringdown of
a Kerr black hole sourced by an extreme mass ratio merger can be modeled by
the greybody factor, which quantifies the scattering nature of the black hole
geometry. The estimation of the mass and spin of the remnant is demonstrated
by fitting the greybody factor to GW data without using black hole quasi-
normal modes. We propose that the ringdown modeling with the greybody factor
may strengthen the test of gravity as one can avoid the possible overfitting
issue and the start time problem in the ringdown modeling with superposed
quasi-normal modes.

††preprint: YITP-23-111††preprint: RIKEN-iTHEMS-Report-23

##  I Introduction

The Kerr solution, describing a spinning black hole, is one of the most
simplest solutions to the Einstein equation. Based on the black hole no-hair
theorem Israel (1967, 1968); Carter (1971), the spacetime structure near an
astrophysical Kerr black hole is characterized by two parameters only, i.e,
the mass $M$ and angular momentum $J$ of the black hole. Therefore, a black
hole is a suitable site to test gravity in strong gravity regimes. In the
context of the test of the no-hair theorem, the black hole spectroscopy
Echeverria (1989); Finn (1992); Berti et al. (2006) has been actively studied
so far. The black hole spectroscopy is an extraction of each black hole quasi-
normal (QN) mode Leaver (1986); Sun and Price (1988); Andersson (1995);
Glampedakis and Andersson (2001, 2003); Nollert and Schmidt (1992); Andersson
(1997); Nollert and Price (1999) from gravitational wave (GW) ringdown which
is a superposition of multiple QN modes. There are infinite number of QN modes
and each mode has a complex frequency $\omega=\omega_{lmn}\in{\mathbb{C}}$
labeled by the overtone number $n$ for each angular and azimuthal mode
$(l,m)$. The real and imaginary part of $\omega_{lmn}$ represent the frequency
and damping rate of the mode, respectively.

GW ringdown appears after the inspiral phase of a binary black hole system. If
the ringdown starts around the strain peak, it would be possible to measure
several QN modes included in a ringdown signal by truncating GW data before
the assumed start time of ringdown and by fitting several QN modes to the
truncated data Giesler et al. (2019). However, some issues in the black hole
spectroscopy with (superposed) QN modes have been pointed out like the start
time problem Sun and Price (1988); Nollert and Price (1999); Andersson (1997);
Berti and Cardoso (2006) and overfitting problem Baibhav et al. (2023). Then,
it would be natural to ask if there is another nice quantity being suitable to
test gravity other than QN modes.

In this paper, we propose that the black hole greybody factor,
$\Gamma_{lm}(\omega)$, would be an important quantity in the test of the no-
hair theorem and the estimation of the remnant parameters from GW ringdown. We
here consider a particle plunging into a massive black hole as a source of GW.
Then we show that for $(l,m)=(2,2)$, $\Gamma_{lm}$ can be imprinted on the GW
spectral amplitude $|\tilde{h}_{lm}(\omega)|$ in $\omega\gtrsim
f_{lm}\equiv\text{Re}(\omega_{lm0})$ with the form of

$|\tilde{h}_{lm}(\omega)|\simeq c_{lm}\times\gamma_{lm}(\omega)\equiv
c_{lm}\times\sqrt{1-\Gamma_{lm}(\omega)}/\omega^{3}\ \text{for}\ \omega\gtrsim
f_{lm},$ (1)

where $\omega$ is a frequency of GWs and $c_{lm}$ is a constant corresponding
to the GW amplitude. The frequency dependence of the greybody factor is
determined by the two remnant parameters only, i.e., the mass and spin of the
black hole. It means that if (1) holds, one can detect the greybody factor
from the ringdown to test the no-hair theorem as the spectral amplitude in
$\omega\gtrsim f_{lm}$ corresponds to the ringdown signal. The reflectivity
${\cal R}_{lm}\equiv 1-\Gamma_{lm}$ has an exponential damping at high
frequencies ($\omega\gtrsim f_{lm}$) and the strength of the damping in the
frequency domain is unique for the remnant mass and spin like a complex QN
mode frequency. As the damping in ${\cal R}_{lm}$ is strong for rapid spins,
the reflectivity well govern the dependence of $|\tilde{h}_{lm}|$ on frequency
and our model works well especially for rapidly spinning remnant black holes.
One of the important difference in the ringdown modeling with the greybody
factor and QN modes is that given the remnant mass $M$ and spin $j(\equiv
J/M^{2})$, we know where the universal damping of ${\cal
R}_{lm}=1-\Gamma_{lm}$ appears in the frequency space, i.e., $\omega\gtrsim
f_{lm}(M,j)$, but it is unknown when the excitation of superposed QN modes
appears in the time domain, which is recognized as the start time problem or
the time-shift problem.

The original idea of the modeling of ringdown with the greybody factor was
introduced in the previous paper by the author Oshita (2023). In this paper,
we investigate the importance of the greybody factors in the ringdown sourced
by an extreme mass ratio merger in more detail. In Sec. II.1, we explain our
methodology to compute GW waveform in the linear perturbation regime. The
definition and the property of the greybody factor is provided in Sec. II.2.
In Sec. III.1, we study why the greybody factor can be imprinted on the GW
ringdown by carefully considering the effect of the source term. In Sec.
III.2, we investigate the exponential damping in the greybody factor and how
it is consistent with the exponential damping in the spectral amplitude of the
GW ringdown. In Sec. III.3, we perform the measurement of the remnant mass and
spin only with the fit of the greybody factor. In Sec. IV, our conclusion is
provided and we discuss the pros and cons of using the greybody factor and QN
modes in the test of the no-hair theorem, measurement of the remnant
quantities, and the modeling of GW ringdown. Throughout the manuscript, we use
the natural unit of $c=\hbar=1$ and $G=1$.

##  II Formalism

In this section, we describe how we compute GW spectral amplitude for an
extreme mass ratio merger and the greybody factor of a spinning black hole.
Here we concentrate on a particle plunging into the hole and its trajectory is
restricted on the equatorial plane.

###  II.1 extreme mass ratio merger and gravitational waveform

The background geometry is approximated by the Kerr spacetime when we consider
an extreme mass ratio merger with a massive black hole. Therefore, the
background geometry can be covered by the Boyer-Lindquist coordinates
$(t,r,\theta,\phi)$ and one can compute the GW spectrum
$\tilde{h}_{lm}(\omega)$ sourced by the merger event in a linear manner. Let
us begin with solving the Sasaki-Nakamura equation Sasaki and Nakamura (1982):

$\left[\frac{d^{2}}{dr^{\ast}{}^{2}}-F_{lm}\frac{d}{dr^{\ast}}-U_{lm}\right]X_{lm}=\rho_{lm},$
(2)

where the explicit forms of $F_{lm}$ and $U_{lm}$ are given in the original
paper by Sasaki and Nakamura Sasaki and Nakamura (1982), and the spectrum
$\tilde{h}_{lm}$ is obtained from the perturbation variable $X_{lm}$ as is
shown later explicitly. The source term $\rho_{lm}$ depends on the plunging
orbit of a particle with mass $\mu$. The form of $\rho_{lm}$ for the plunging
particle on the equatorial plane ($\theta=\pi/2$) is Kojima and Nakamura
(1984)

$\rho_{lm}=\frac{\gamma_{0}\Delta}{(r^{2}+a^{2})^{3/2}r^{2}}W\exp{\left(-i\int^{r}\frac{K(r^{\prime})}{\Delta(r^{\prime})}dr^{\prime}\right)},$
(3)

where $\Delta(r)\equiv r^{2}-2Ma+a^{2}$, $K(r)\equiv(r^{2}+a^{2})\omega-am$,
and $a\equiv J/M$. The functions $\gamma_{0}$ and $W$ are shown in Appendix A
and B in Ref. Kojima and Nakamura (1984), respectively.

Figure 1: Spectral amplitude $|\tilde{h}_{22}|$ for $j=0.7$ (black dash-
dotted), $0.9$ (red solid), and $0.99$ (blue dashed). The source term is
obtained with $L_{z}=0.5$ and the observation angle is $\theta=\pi/2$. The
black vertical lines show the real part of the fundamental QN mode frequency
$f_{22}(j)$. The exponential damping of the spectral amplitudes appear in
$\omega\gtrsim f_{22}$.

The trajectory of a particle is determined by the following differential
equations Carter (1968); Kojima and Nakamura (1984):

$\displaystyle r^{2}\frac{dt}{d\tau}$
$\displaystyle=-a(a-L_{z})+\frac{r^{2}+a^{2}}{\Delta}P,$ (4) $\displaystyle
r^{2}\frac{d\phi}{d\tau}$ $\displaystyle=-a(a-L_{z})+\frac{a}{\Delta}P,$ (5)
$\displaystyle r^{2}\frac{dr}{d\tau}$ $\displaystyle=-\sqrt{R},$ (6)
$\displaystyle\theta$ $\displaystyle=\pi/2,$ (7)

where $P\equiv r^{2}+a^{2}-L_{z}a$, $R\equiv
2Mr^{3}-L_{z}^{2}r^{2}+2Mr(L_{z}-a)^{2}$, $\mu L_{z}$ is the orbital angular
momentum, and $\tau$ is the proper time of the particle. We obtain the source
term $\rho_{lm}$ by substituting the trajectory of the particle,
$(t(\tau),r(\tau),\theta=\pi/2,\phi(\tau))$, into (3). We then numerically
compute the Sasaki-Nakamura equation with the source term for a plunging orbit
on the equatorial plane $\rho_{lm}$. The GW spectral amplitude we obtained for
$l=m=2$ are shown in Figure 1.

###  II.2 greybody factors

spin parameter $(j)$ | 0.001 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6
---|---|---|---|---|---|---|---
decay frequency $(T_{22})$ | 0.067 | 0.066 | 0.065 | 0.064 | 0.063 | 0.062 | 0.060
spin parameter $(j)$ | 0.7 | 0.8 | 0.9 | 0.95 | 0.99 | 0.995 | 0.998
---|---|---|---|---|---|---|---
decay frequency $(T_{22})$ | 0.057 | 0.053 | 0.045 | 0.036 | 0.019 | 0.014 | 0.0096
Table 1: The value of $T_{lm}$ $(l=m=2)$ with respect to the spin parameter
$j$. We read the value of $T_{22}$ from $1-\Gamma_{22}(\omega)$ in the range
of $\alpha\times f_{22}\leq\omega\leq 1.99/(2M)$ where we choose a constant
$\alpha$ in the range of $1\leq\alpha\leq 1.2$. The constant $\alpha$ should
be larger for a lower spin in order for $1-\Gamma_{22}$ to be well
approximated with $e^{-(\omega-f_{22})/T_{22}}$ in the frequency range. The
value of $\alpha$ we set is shown in (23) in Appendix A.

The greybody factor quantifies the absorptive nature of a black hole geometry
and is independent of the source term. It is determined only by the no-hair
parameters of a black hole, i.e., the mass and spin of a black hole, like the
black hole quasi-normal modes. We obtain the greybody factor by computing a
homogeneous solution to the Sasaki-Nakamura equation
$X_{lm}=X_{lm}^{\rm(hom)}$ with the boundary condition of

$X_{lm}^{\rm(hom)}=e^{-ik_{\rm H}r^{\ast}}\ \text{for}\ r^{\ast}\to-\infty,$
(8)

where $k_{\rm H}\equiv\omega-m\Omega_{\rm H}$ with $\Omega_{\rm H}\equiv
j/(2r_{+})$. We then read the asymptotic ingoing and outgoing amplitudes at a
distant region as

$X_{lm}^{\rm(hom)}=A_{\rm in}e^{-i\omega r^{\ast}}+A_{\rm out}e^{i\omega
r^{\ast}}\ \text{for}\ r^{\ast}\to\infty.$ (9)

The reflectivity of the angular momentum barrier is given by Brito et al.
(2015); Nakano et al. (2017)

${\cal R}_{lm}\equiv\left|\frac{C}{c_{0}}\right|^{2}\left|\frac{A_{\rm
out}}{A_{\rm in}}\right|^{2}\equiv 1-\Gamma_{lm},$ (10)

where ${\cal R}_{lm}$ and $\Gamma_{lm}$ are the reflectivity and the greybody
factor (i.e., transmissivity), respectively. The factors $|C|^{2}$ and $c_{0}$
are Starobinsky (1973); Starobinsky and Churilov (1974)

$\displaystyle|C|^{2}$
$\displaystyle\equiv\lambda^{4}+4\lambda^{3}+\lambda^{2}(-40a^{2}\omega^{2}+40am\omega+4)+48a\lambda\omega(a\omega+m)+144\omega^{2}(a^{4}\omega^{2}-2a^{3}m\omega+a^{2}m^{2}+M^{2}),$
(11) $\displaystyle c_{0}$
$\displaystyle\equiv\lambda(\lambda+2)-12a\omega(a\omega-m)-i12M\omega,$ (12)

respectively, and $\lambda$ is the separation constant of the spin-weighted
spheroidal harmonics. We numerically compute the greybody factor by solving
the Sasaki-Nakamura equation111For more details of our numerical computation,
see Appendix A. Recently, the greybody factor of Kerr black holes was
analytically computed by solving the connection problem of the confluent Heun
equation and its analytic form was obtained in Ref. Bonelli et al. (2022)..
Our computation reproduce the exponential decay of ${\cal R}_{22}$ at high
frequencies ($\omega\gtrsim f_{22}$) and the superradiant amplification at
$\omega<m\Omega_{\rm H}$ as is shown in Figure 2. The exponential damping of
${1-\Gamma_{lm}}$ at high-frequency region ($\omega\gtrsim f_{lm}$) can be
approximated as

$1-\Gamma_{lm}\simeq e^{-(\omega-f_{lm})/T_{lm}},$ (13)

where $T_{lm}$ quantifies the strength of the exponential damping of ${\cal
R}_{lm}$ in the frequency domain. The factor $T_{lm}$ is a no-hair quantity
which depends only on the mass and spin of the remnant black hole like the QN
mode frequency. The values of $T_{lm}$ extracted from our numerical data are
shown in Table 1 and the fitting methodology we used is provided in Appendix
A.

Figure 2: Reflectivity ${\cal R}_{lm}=1-\Gamma_{lm}$ for $(l,m)=(2,2)$ with
$j=0.7$ (black dot-dashed), $0.9$ (red solid), and $0.99$ (blue dashed). The
black vertical lines show the real part of the fundamental QN mode frequency
$\omega=f_{22}$. The exponential damping of the reflectivity appear in
$\omega\gtrsim f_{22}$.

##  III Greybody Factors in Ringdown

In this Section, we study the greybody factor imprinted on GW ringdown. As is
shown in Figure 3, we find that the greybody factor can model the spectral
amplitude of GW ringdown in $\omega\gtrsim f_{lm}$ when the GW is sourced by
an extreme mass ratio merger. This still holds even for some higher harmonic
modes (Figure 4) and for various values of the orbital angular momentum $\mu
L_{z}$ and the spin parameter $j$ of the massive black hole. We also confirm
that the frequency dependence of GW spectrum is insensitive to the observation
angle $\theta$ as shown in Appendix B. This universal nature in the black hole
ringdown is important to test the no-hair theorem with high precision by
combining the black hole spectroscopy.

Figure 3: The comparison of the GW spectral amplitude $|\tilde{h}_{22}|$ and a
function $\gamma_{22}(\omega)=\sqrt{1-\Gamma_{22}}/\omega^{3}$ determined by
the greybody factor $\Gamma_{22}$. The observation angle $\theta$ is set to
$\pi/2$. For other values of $\theta$, see Appendix B. The match between the
two quantities for $\omega\gtrsim f_{22}$ is explained in Sec. III.1. The
vertical black lines indicate $\omega=f_{22}$. Figure 4: The comparison of the
GW spectral amplitude $|\tilde{h}_{lm}|$ and a function
$\gamma_{lm}(\omega)=\sqrt{1-\Gamma_{lm}}/\omega^{3}$ for $(l,m)=(3,3)$,
$(4,4)$ and $(5,5)$. The vertical black lines indicate $\omega=f_{lm}$. We set
$j=0.9$, $L_{z}=0.5$, and $\theta=\pi/2$.

We here study how the greybody factor can be imprinted on the ringdown signal.
We also demonstrate the estimation of the remnant mass and spin by using the
greybody factor only. We then find that it works well, which implies that the
greybody factor is important to test the black hole no-hair theorem. Not only
using the QN modes but also using the greybody factor would enhance the
accuracy of the test of gravity and measurability of the remnant quantities at
least for extreme mass ratio mergers.

###  III.1 Why are the greybody factors imprinted on ringdown?

We here discuss why the greybody factors can be imprinted on the ringdown for
extreme mass ratio mergers. For higher mass ratios, the background geometry is
governed only by a massive black hole and the perturbation theory in the Kerr
spacetime works to compute GW waveform. The GW strain, $h$, is given by

$\displaystyle h=\sum_{l,m}\frac{e^{im\phi}}{\sqrt{2\pi}r}\int
d\omega\tilde{h}_{lm}(\omega)e^{-i\omega t}$
$\displaystyle=\sum_{l,m}\frac{e^{im\phi}}{\sqrt{2\pi}r}\int
d\omega\frac{-2}{\omega^{2}}{}_{-2}S_{lm}(a\omega,\theta)R_{lm}(\omega)e^{-i\omega
t},$ (14) $\displaystyle=\sum_{l,m}\frac{e^{im\phi}}{\sqrt{2\pi}r}\int
d\omega\frac{-2}{\omega^{2}}{}_{-2}S_{lm}(a\omega,\theta)\frac{A_{\rm
out}(\omega)}{2i\omega A_{\rm in}(\omega)}\tilde{\rho}_{lm}(\omega)e^{-i\omega
t},$ (15)

where ${}_{-2}S_{lm}(a\omega,\theta)$ is the spin-weighted spheroidal
harmonics, $R_{lm}$ is the radial Teukolsky variable, and $\tilde{\rho}_{lm}$
is

$\tilde{\rho}_{lm}(\omega)=\frac{-4\omega^{2}}{\lambda(\lambda+2)-12iM\omega-12a^{2}\omega^{2}}\int_{r_{+}}^{\infty}dr^{\prime}\frac{\rho_{lm}(\omega,r^{\prime})X_{lm}^{\rm(hom)}(\omega,r^{\prime})}{A_{\rm
out}(\omega)},$ (16)

Then we find

$|\tilde{h}_{lm}|=\frac{\sqrt{1-\Gamma_{lm}}}{\omega^{3}}t_{lm}=\gamma_{lm}t_{lm},$
(17)

where $\gamma_{lm}$ is a universal quantity that depends only on the two
remnant quantities ($M,j$) and another factor $t_{lm}$ includes the source
term and the spheroidal harmonics, depending on the external information like
the GW source and the observation angle, respectively:

$t_{lm}\equiv\left|\frac{c_{0}}{C}\right||\tilde{\rho}_{lm}||{}_{-2}S_{lm}|.$
(18)

The factor of $t_{lm}$ is hereinafter referred to as the renormalized source
term. The frequency dependence of the GW spectral amplitude,
$|\tilde{h}_{lm}|$, is governed by $\gamma_{lm}$ at higher frequencies
$\omega\gtrsim f_{lm}$ provided that $t_{lm}$ has the small dependence in
$\omega$ for $\omega\gtrsim f_{lm}$. The value of the source term $t_{22}$ is
shown in Figure 5. One can see that $t_{22}(\omega)$ is indeed nearly constant
and $\gamma_{22}$ governs the frequency dependence of the GW spectrum. It
might imply that a compact object plunging into a large black hole can be
regarded as an instantaneous source of GW ringdown and the associated source
term can be nearly constant in the frequency domain222Remember that an
instantaneous pulse like a delta function or a sharp Gaussian distribution in
the time domain has a (nearly) constant distribution in the frequency domain..
We leave a more detailed study of our ringdown model for other harmonic modes,
e.g., the sensitivity of our model for $(l,m)\neq(2,2)$ to external parameters
like the orbital angular momentum, for a future work.

Figure 5: Plot of the GW spectral amplitude $|\tilde{h}_{22}|$ (black solid),
$\gamma_{22}$ (blue dot-dashed), and $t_{22}$ (red dashed) for $j=0.8$,
$L_{z}=0$, and $\theta=\pi/2$.

###  III.2 Exponential decay in GW spectral amplitudes and in the greybody
factors

We find the spectral amplitude in high frequency region ($\omega\gtrsim
f_{22}$) can be modeled by the greybody factor as is shown in Figures 3 and 5.
This holds for various values of the orbital angular momentum of the plunging
particle $\mu L_{z}$. Fitting the Boltzmann factor $\exp[-(\omega-
f_{22})/T^{\rm(GW)}_{22}]$ to the simulated GW data333The detailed methodology
of our fitting analysis is provided in Appendix A., we read the damping
exponent of the GW spectral amplitude $T^{\rm(GW)}_{22}$ in the frequency
domain. The result is shown in Figure 6 and the best fit values of
$T_{22}^{\rm(GW)}$ (dots) is consistent with $T_{22}$ (solid line) especially
for $j\gtrsim 0.8$. The value of $T_{22}$ is sensitive to the spin parameter
$j$ for rapid spins, but is insensitive to $j$ for lower spins (see table 1 as
well).

In addition to $T_{lm}$, another quantity $f_{lm}$ is also important to model
$|\tilde{h}_{lm}|$ as the exponential damping in $|\tilde{h}_{lm}|$ appears at
$\omega\gtrsim f_{lm}$ (see Figures 1 and 3). In the next section, we show
that the two remnant values, i.e., $M$ and $j$, can be extracted from the GW
spectral amplitude by fitting $\gamma_{lm}$ characterized by $T_{lm}(M,j)$ and
$f_{lm}(M,j)$.

Figure 6: The dots show the best fit values of $T^{\rm(GW)}_{22}$ extracted
from the numerical GW waveform data with $\theta=\pi/2$. The black solid line
is the best fit value of the damping frequency $T_{22}$ obtained from the
numerical computation of the greybody factor.

###  III.3 estimation of the remnant quantities

The two no-hair quantities $(M,j)$ can be extracted by fitting the function of
$\gamma_{lm}(M,j)\equiv\sqrt{1-\Gamma_{lm}(M,j)}/\omega^{3}$ to the spectral
amplitude of GW data $|\tilde{h}_{lm}|$ as the greybody factor is
characterized by the two remnant quantities. The fitting parameters are $M$,
$j$, and an amplitude $c_{lm}$. Here we demonstrate the extraction of the two
no-hair parameters $(M,j)$ from our clean numerical GW waveform by fitting
$\gamma_{lm}$ with $(l,m)=(2,2)$. We here use the analytic model function that
models $\gamma_{22}$, whose explicit form is shown in Appendix C. The
estimation of $(M,j)$ with noise is important to quantify the feasibility for
a specific detector, and it will be studied elsewhere.

We estimate the mismatch ${\cal M}$ between the GW spectral amplitude
$|\tilde{h}_{22}|$ and $c_{22}\times\gamma_{22}$ on the mass-spin space with

${\cal
M}(M,a)=\left|1-\frac{\braket{}{\tilde{h}_{22}}{c_{22}\times\gamma_{22}}}{\sqrt{\braket{}{\tilde{h}_{22}}{|\tilde{h}_{22}|}\braket{c_{22}\times\gamma_{22}}{c_{22}\times\gamma_{22}}}}\right|=\left|1-\frac{\braket{}{\tilde{h}_{22}}{\gamma_{22}}}{\sqrt{\braket{}{\tilde{h}_{22}}{|\tilde{h}_{22}|}\braket{\gamma_{22}}{\gamma_{22}}}}\right|,$
(19)

where $\braket{a(\omega)}{b(\omega)}$ is

$\braket{a(\omega)}{b(\omega)}=\int_{\omega_{i}}^{\omega_{f}}d\omega
a(\omega)b^{\ast}(\omega).$ (20)

Note that the mismatch ${\cal M}$ is independent of the scale $c_{22}$ and
depends only on the other two fitting parameters $(M,j)$ only. This makes the
fit and extraction of the remnant quantities quite simpler than the case in
the multiple QN mode fitting. Also, we could avoid the overfitting issue. For
the fit of multiple overtones, on the other hand, there are many fitting
parameters, i.e., an amplitude and phase for each QN mode. It was pointed out
Baibhav et al. (2023) that the inclusion of many QN modes in the ringdown
model may cause overfitting when we use a GW waveform beginning with the
strain peak Baibhav et al. (2023)444On the other hand, the previous work of
Ref. Giesler et al. (2019) fit multiple QN modes to the numerical relativity
GW waveform beginning from the strain peak. Then they reproduced the injected
remnant mass and spin values. This implies that the fit of multiple QN modes
may work at least when GW data has no contamination by noise..

We estimate the mismatch ${\cal M}$ by computing the inner product (20) with
the range of the integral of $\omega_{i}=2M$ and $\omega_{f}=2M\times 1.99$.
Note that the $M$ in $\omega_{i/f}$ is not the true value but the fitting
parameter of the black hole mass. The mismatch is computed in the mass-spin
domain and the result is shown in Figure 7. We find that the mass-spin
estimation works well even though we here use the greybody factor without the
fit of multiple QN modes. We also find that the best fit mass and spin are not
sensitive to an artificial choice of the range of the data we use
$\omega\in[\omega_{i},\omega_{f}]$ as is shown in Figure 8. On the other hand,
in the fit of QN modes, the mass-spin measurement is sensitive to the assumed
start time of ringdown Giesler et al. (2019). Although the feasibility of the
extraction of the greybody factor depends on noise, combining this with the
black hole spectroscopy may strengthen not only the measurability of the
remnant quantities but also the precision of the test of gravity. We will come
back to this in the future.

Figure 7: The mismatch ${\cal M}$ between $|\tilde{h}_{22}|$ and $\gamma_{22}$
for (a) $j=0.7$, (b) $0.9$, and (c,d) $0.99$. The injected (true) values of
the remnant quantities are indicated with the white solid lines. The source
particle has the orbital angular momentum of $\mu L_{z}=0.5\mu$ and the
observation angle is set to $\theta=\pi/2$. The frequency range used in the
estimation of the mismatch is $[2M\omega_{i},2M\omega_{f}]=[1,1.99]$. Figure
8: The mismatch ${\cal M}$ between $|\tilde{h}_{22}|$ and $\gamma_{22}$ for
$j=0.9$. The injected (true) values of the remnant quantities are indicated
with the white solid lines. The source particle has the orbital angular
momentum of $\mu L_{z}=0.5\mu$ and we set $\theta=\pi/2$. We change the
frequency range of the data used in the computation of ${\cal M}$ as
$2M\omega_{i}=1$, $1.1$, $1.25$ and $2M\omega_{f}$ is fixed to $1.99$.

##  IV Discussions

The superposed QN modes is one of the most established model of the black hole
ringdown. In this paper, we discussed another universal nature of ringdown
that is described by the black hole greybody factor $\Gamma_{lm}$. We
considered how GW ringdown can be modeled by the greybody factor, which is
another no-hair quantity that depends only on the mass and spin of the remnant
black hole like the black hole QN modes. We found that the spectral amplitude
of GW ringdown $|\tilde{h}_{lm}|$ with $(l,m)=(2,2)$ sourced by an extreme
mass ratio merger can be modeled by $|\tilde{h}_{lm}|\sim
c\times\gamma_{lm}(\omega)$ for $\omega\gtrsim f_{lm}$ where
$\gamma_{lm}(\omega)\equiv\sqrt{1-\Gamma_{lm}}/\omega^{3}$ and $c$ is an
amplitude. In order for the greybody factor to be imprinted on the ringdown,
the renormalized source term $t_{lm}(\omega)$, depending on the GW source,
should be nearly constant with respect to $\omega$ for $\omega\gtrsim f_{lm}$.
We confirmed that $t_{22}(\omega)$ satisfies the condition when GW is sourced
by a compact object plunging into a massive black hole in the extreme mass
ratio regime (Figure 5). We may expect that this is the case as long as a
particle plunging into a massive black hole can be regarded as an
instantaneous source of GW ringdown. We numerically computed GW waveforms with
several values of the orbital angular momentum $\mu L_{z}$. We then confirmed
that the GW spectral amplitude is well modeled by $\gamma_{22}$, determined by
the greybody factor, at higher frequencies $\omega\gtrsim f_{22}$ for various
values of $|L_{z}|\lesssim 0.5$ (see Figures 3 and 6). Also, this model works
well especially for a rapidly spinning remnant black hole. Indeed, the
measurement of the innermost stable circular orbit of supermassive black holes
(SMBHs) based on the X-ray observation puts the lower bound on the spin $j$ of
SMBHs, and some of them take $j>0.9$ and can be even near extremal as $j>0.99$
Reynolds (2019, 2021).

As the greybody factor $\Gamma_{lm}$ is another no-hair quantity, the
extraction of not only the QN modes but also the greybody factor from GW
ringdown would improve the accuracy of the measurement of the remnant mass and
spin and strengthens the test of gravity (Figure 7). The pros and cons in the
modeling of GW ringdown with QN modes and greybody factors are summarized
below.

  1. 1.

For the ringdown modeling with QN modes, the relevant data range in the time
domain $t\gtrsim t_{\rm start}$ is difficult to identify, where $t_{\rm
start}$ is the start time of ringdown. On the other hand, for the ringdown
modeling with the greybody factor, the relevant data range $\omega\gtrsim
f_{lm}(M,j)$ is uniquely determined once we fix the remnant quantities $M$ and
$j$.

  2. 2.

Many fitting parameters are needed to extract QN modes from GW ringdown
especially when several QN modes are excited simultaneously as
$\sum_{n}C_{lmn}\exp[-i\omega_{lmn}t+\varphi_{lmn}]$ for a dominant angular
mode of $(l,m)$. For the extraction of the greybody factor, on the other hand,
the spectral amplitude of the black hole ringdown at $\omega\gtrsim f_{lm}$ is
modeled by $c\times\gamma_{lm}(M,j)$. The scale $c$ is irrelevant for the
minimization of the mismatch ${\cal M}$. As such, one can search the least
value of ${\cal M}$ with the only two fitting parameters $(M,j)$ while
avoiding the overfitting issue. It is much simpler than the QN mode fitting
which involves many fitting parameters, i.e., amplitude $C_{lmn}$ and phase
$\varphi_{lmn}$ for each QN mode.

  3. 3.

GW ringdown can be modeled by the superposition of QN modes regardless of the
frequency-dependence of the source term. However, the modeling of ringdown
with the greybody factor does not always work due to the contamination from
the source term. Note that the greybody factor can be extracted only when the
normalized source term $t_{lm}(\omega)$ is nearly constant in $\omega\gtrsim
f_{lm}$ (Figure 5).

Given the pros and cons in the modeling of ringdown with the greybody factor
and in that with QN modes, combining those two models may improve the test of
the no-hair theorem and the estimation of the remnant quantities. We could
also relate the excitation of overtones with the greybody factors as the
residue of $\gamma_{lm}$ at QN modes can be regarded as the excitation factor,
which quantifies the excitability of each QN mode Leaver (1986); Berti and
Cardoso (2006); Zhang et al. (2013); Oshita (2021). It would be important to
understand the relation between the greybody factor and excitation factor to
reveal the universality in the black hole ringdown.

To further confirm the importance of the greybody factor in the modeling of GW
ringdown, we have to check the detectability of the greybody factor from GW
ringdown with the future detectors such as LISA. Also, it would be important
to take into account some higher harmonic modes, which would affect the
extraction of the greybody factor and increases the fitting parameters if
higher harmonic modes are significantly excited. We will come back to these
points in the future. It is interesting to note that as another different
direction, the authors in Ref. Völkel et al. (2019) studied an inverse problem
to read the greybody factor from quantum Hawking radiation. An interesting
aspect of the greybody factor is that it can be important in both quantum and
classical radiation of black holes, i.e., Hawking radiation and GW ringdown,
respectively.

###### Acknowledgements.

The author appreciate Niayesh Afshordi, Kazumasa Okabayashi, and Hidetoshi
Omiya for valuable comments on this work. The author thanks Daiki Watarai for
carefully reading an earlier version of this paper and for valuable comments.
The author also thanks Giulio Bonelli and Sebastian Völkel for sharing their
recent works and for helpful comments on an earlier version of this paper. The
author is supported by the Grant-in-Aid for Scientific Research (KAKENHI)
project for FY2023 (23K13111).

##  Appendix A numerical methodology and accuracy

We numerically solve the Sasaki-Nakamura equation (2) with the 4th Runge-Kutta
method. The source term $\rho_{lm}$ for the plunging particle is numerically
computed in the range of $r^{\ast}_{\rm min}\leq r^{\ast}\leq r^{\ast}_{\rm
max}=400M$. The minimum radius of the range of integral $r^{\ast}_{\rm min}$
is set to

$\displaystyle r^{\ast}_{\rm min}=\begin{cases}-40M,\ &\text{for}\ j\leq
0.9,\\\ -50M,\ &\text{for}\ 0.9<j<0.97,\\\ -80M,\ &\text{for}\ 0.97\leq
j<0.99,\\\ -100M,\ &\text{for}\ j=0.99,\\\ -120M,\ &\text{for}\ j=0.995,\\\
-160M,\ &\text{for}\ j=0.998.\end{cases}$ (21)

For the Sasaki-Nakamura equation, the exponential tail of the potential
$U_{lm}$ near the horizon becomes long range as $j\to 1$. As such,
$r^{\ast}_{\rm min}$ should be a larger negative value for rapid spins so that
one can impose the boundary condition of $e^{-ik_{\rm H}r^{\ast}}$ at the end
point of $r^{\ast}=r^{\ast}_{\rm min}$. The numerical integration of the
Sasaki-Nakamura equation is done in the range of $r^{\ast}_{\rm min}\leq
r^{\ast}\leq r^{\ast}_{\rm SNmax}=300+30/\omega$ for each frequency mode of
$\omega$.

Figure 9: GW spectral amplitude for $j=0.8$, $(l,m)=(2,2)$, $L_{z}=0.5$ and
$\theta=\pi/2$. We compute it for high, medium, and low resolutions with
$(N_{\rm source},N_{\rm SN})=(10^{4},2\times 10^{5})$, $(5\times
10^{3},10^{5})$, and $(2.5\times 10^{3},5\times 10^{4})$, respectively.

The source term $\rho_{lm}(\omega,r^{\ast})$ is obtained in the resolution of
$\Delta r^{\ast}=(r^{\ast}_{\rm max}-r^{\ast}_{\rm min})/N_{\rm source}$ with
$N_{\rm source}=5000$. The Sasaki-Nakamura equation is integrated with the
step size of $\Delta r^{\ast}=(r^{\ast}_{\rm SNmax}-r^{\ast}_{\rm min})/N_{\rm
SN}$ with $N_{\rm SN}=10^{5}$. The greybody factor is computed by reading the
asymptotic amplitude at $r^{\ast}=r^{\ast}_{\rm SNmax}$ by using the
Wronskian. We checked that our resolution is high enough to obtain high-
accuracy GW waveform and greybody factor (see Figure 9).

The damping frequency $T_{lm}$ in $1-\Gamma_{lm}(\omega)$ and
$T_{lm}^{\rm(GW)}$ in the GW spectral amplitude $|\tilde{h}_{lm}(\omega)|$ are
extracted at higher frequencies $\omega\gtrsim f_{lm}$ by using a
Mathematica’s function NonlinearModelFit for the log-scaled data,
$\log\Gamma_{lm}$ and $\log(|\tilde{h}_{lm}|^{2}\omega^{6})$, with the fitting
function of

$B(\omega-\omega_{i})+\log A,$ (22)

where $A$ and $B$ are the fitting parameters and $B$ is associated with
$T_{lm}$ or $T^{\rm(GW)}_{lm}$. The results for $(l,m)=(2,2)$ are shown in
Table 1 and Figure 6. The extraction of $T_{22}$ is done by fitting the
Boltzmann factor $e^{-(\omega-f_{22})/T_{22}}$ to the data in the frequency
range of $\omega_{i}\leq\omega\leq\omega_{f}=1.99/(2M)$ with

$\displaystyle\omega_{i}=\begin{cases}1.20\times f_{22},\ &\text{for}\
0.001\leq j\leq 0.75,\\\ 1.10\times f_{22},\ &\text{for}\ 0.8\leq j\leq
0.9,\\\ 1.05\times f_{22},\ &\text{for}\ 0.93\leq j\leq 0.99,\\\ 1.02\times
f_{22},\ &\text{for}\ j=0.995,\\\ 1.00\times f_{22},\ &\text{for}\
j=0.998.\end{cases}$ (23)

For the extraction of $T_{22}^{\rm(GW)}$, we fit the Boltzmann factor to the
numerical data in the range of $\omega_{i}\leq\omega\leq\omega_{f}$ with

$\displaystyle\omega_{i}=\begin{cases}1.20\times f_{22},\ &\text{for}\
0.001\leq j\leq 0.75,\\\ 1.10\times f_{22},\ &\text{for}\ 0.8\leq j\leq
0.9,\\\ 1.05\times f_{22},\ &\text{for}\ 0.93\leq j<0.98,\\\ 1.00\times
f_{22},\ &\text{for}\ 0.98\leq j\leq 0.998,\end{cases}$ (24)

and $\omega_{f}$ is set to a value at which
$(\omega_{f}^{3}|\tilde{h}_{22}(\omega_{f})|^{2})/(\omega_{i}^{3}|\tilde{h}_{22}(\omega_{i})|^{2})\simeq
0.01$. Also, the best fit value and error of $T_{22}^{\rm(GW)}$ in Figure 6
was estimated by Mathematica’s commands BestFitParameters and ParameterErrors
in a Mathematica’s function of NonlinearModelFit.

Figure 10: Comparison of the spectral amplitude of GW $|\tilde{h}_{22}|$
(black solid) and $\gamma_{22}$ (red dashed). We set $j=0.9$ and $L_{z}=0.5$.
We also set the observation angle as $\theta=\pi/6$, $\pi/4$, $\pi/2$,
$3\pi/4$, and $5\pi/6$. The mismatch ${\cal M}$ is evaluated for data in
$f_{22}\leq\omega\leq 1.99/(2M)$.

##  Appendix B greybody factor in GW ringdown and observation angle $\theta$

Our ringdown modeling for a harmonic mode $(l,m)$ is given by the product of
$\gamma_{lm}$ and the renormalized source term (17). The renormalized source
term is determined by a source of GW emission and the observation angle as it
includes the spin-weighted spheroidal harmonics
${}_{-2}S_{lm}(a\omega,\theta)$. We confirmed that the ringdown modeling with
the greybody factor works for a wide range of the observation angle $\theta$.
Indeed, the mismatch ${\cal M}$ defined in (19) is less than $10^{-3}$ at
least for $\pi/6\leq\theta\leq 5\pi/6$ as is shown in Figure 10. The mismatch
${\cal M}$ is evaluated for data in $f_{22}\leq\omega\leq 1.99/(2M)$.

##  Appendix C analytic model of the greybody factor

Our proposal in this paper is that the greybody factor $\Gamma_{lm}$ is
imprinted on the spectral amplitude of GW ringdown with the form of

$\gamma_{lm}=\sqrt{1-\Gamma_{lm}}/\omega^{3}.$ (25)

As the greybody factor is a universal quantity which depends only on the
remnant mass and spin like the black hole QN modes, the extraction of the
greybody factor from the signal is applicable to test the no-hair theorem and
the measurement of the remnant mass and spin. To demonstrate that in Sec.
III.3, we compute the mismatch ${\cal M}$ between GW spectral amplitude and
the function $\gamma_{lm}$. As the computation of the greybody factor involves
the numerical integration of the Sasaki-Nakamura equation in our approach, we
shorten the computation time of ${\cal M}$ by using an analytic model function
$\tilde{\Gamma}_{22}$ that models the greybody factor $\Gamma_{22}$ for
$\omega>0$555Another fitting function of the reflectivity for $0.6<j<0.8$ is
provided in Ref. Nakano et al. (2017).:

$1-\tilde{\Gamma}_{22}(\omega)=\frac{1+a_{1}Z[-2,2,2,M,j,\omega](1-\tanh{[(\omega-
f_{22})/a_{2}]})}{(1+\exp[(\omega-f_{22})/T_{22}])}\ \text{for}\ \omega>0,$
(26)

where $a_{1}=0.325$, $a_{2}=0.02$, and

$\displaystyle T_{22}$ $\displaystyle\simeq
0.223\sqrt{1-j}-0.33(1-j)+0.249(1-j)^{3/2}-0.0748(1-j)^{2},$ (27)
$\displaystyle f_{22}$ $\displaystyle\simeq
2-2.85\sqrt{1-j}+3.01(1-j)-2.01(1-j)^{3/2}+0.597(1-j)^{2},$ (28)
$\displaystyle Z[s,l,m,M,j,\omega]$ $\displaystyle\equiv 4m\Omega_{\rm
H}\frac{r_{+}}{\sqrt{1-j^{2}}}\left(\frac{(l-s)!(l+s)!}{(2l)!(2l+1)!!}\right)^{2}[2M\omega(1-j^{2})]^{2l+1}\prod_{k=1}^{l}\left(1+\frac{4}{k^{2}}\left(m\Omega_{\rm
H}\frac{r_{+}}{\sqrt{1-j^{2}}}\right)^{2}\right),$ (29)

where $s$ is the spin of the relevant field, e.g., $|s|=2$ for gravitational
field. This fitting model matches with the exact greybody factor within
$\tilde{{\cal M}}\lesssim 0.01$ as is shown in Figure 11. This fitting
function is applicable to the broad range of spin parameter $0.001\leq j\leq
0.998$ as is partially shown in the Figure.

Figure 11: Reflectivity ${\cal R}_{22}=1-\Gamma_{22}$ is shown in the linear
and log scale for several spin parameters in $0.001\leq j\leq 0.998$. The
mismatch $\tilde{\cal M}$ between the reflectivity obtained by our numerical
computation $1-\Gamma_{22}$ and the one modeled by the analytic function in
(26).

## References

  * Israel (1967) W. Israel, Phys. Rev. 164, 1776 (1967).
  * Israel (1968) W. Israel, Commun. Math. Phys. 8, 245 (1968).
  * Carter (1971) B. Carter, Phys. Rev. Lett. 26, 331 (1971).
  * Echeverria (1989) F. Echeverria, Phys. Rev. D 40, 3194 (1989).
  * Finn (1992) L. S. Finn, Phys. Rev. D 46, 5236 (1992), eprint gr-qc/9209010.
  * Berti et al. (2006) E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D 73, 064030 (2006), eprint gr-qc/0512160.
  * Leaver (1986) E. W. Leaver, Phys. Rev. D 34, 384 (1986).
  * Sun and Price (1988) Y. Sun and R. H. Price, Phys. Rev. D 38, 1040 (1988).
  * Andersson (1995) N. Andersson, Phys. Rev. D 51, 353 (1995).
  * Glampedakis and Andersson (2001) K. Glampedakis and N. Andersson, Phys. Rev. D 64, 104021 (2001), eprint gr-qc/0103054.
  * Glampedakis and Andersson (2003) K. Glampedakis and N. Andersson, Class. Quant. Grav. 20, 3441 (2003), eprint gr-qc/0304030.
  * Nollert and Schmidt (1992) H.-P. Nollert and B. G. Schmidt, Phys. Rev. D 45, 2617 (1992).
  * Andersson (1997) N. Andersson, Phys. Rev. D 55, 468 (1997), eprint gr-qc/9607064.
  * Nollert and Price (1999) H.-P. Nollert and R. H. Price, J. Math. Phys. 40, 980 (1999), eprint gr-qc/9810074.
  * Giesler et al. (2019) M. Giesler, M. Isi, M. A. Scheel, and S. Teukolsky, Phys. Rev. X 9, 041060 (2019), eprint 1903.08284.
  * Berti and Cardoso (2006) E. Berti and V. Cardoso, Phys. Rev. D 74, 104020 (2006), eprint gr-qc/0605118.
  * Baibhav et al. (2023) V. Baibhav, M. H.-Y. Cheung, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, and F. Duque (2023), eprint 2302.03050.
  * Oshita (2023) N. Oshita, JCAP 04, 013 (2023), eprint 2208.02923.
  * Sasaki and Nakamura (1982) M. Sasaki and T. Nakamura, Prog. Theor. Phys. 67, 1788 (1982).
  * Kojima and Nakamura (1984) Y. Kojima and T. Nakamura, Prog. Theor. Phys. 71, 79 (1984).
  * Carter (1968) B. Carter, Phys. Rev. 174, 1559 (1968).
  * Brito et al. (2015) R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015), eprint 1501.06570.
  * Nakano et al. (2017) H. Nakano, N. Sago, H. Tagoshi, and T. Tanaka, PTEP 2017, 071E01 (2017), eprint 1704.07175.
  * Starobinsky (1973) A. A. Starobinsky, Sov. Phys. JETP 37, 28 (1973).
  * Starobinsky and Churilov (1974) A. A. Starobinsky and S. M. Churilov, Sov.Phys.-JETP 38, 1 (1974).
  * Bonelli et al. (2022) G. Bonelli, C. Iossa, D. P. Lichtig, and A. Tanzini, Phys. Rev. D 105, 044047 (2022), eprint 2105.04483.
  * Reynolds (2019) C. S. Reynolds, Nature Astron. 3, 41 (2019), eprint 1903.11704.
  * Reynolds (2021) C. S. Reynolds, Ann. Rev. Astron. Astrophys. 59, 117 (2021), eprint 2011.08948.
  * Zhang et al. (2013) Z. Zhang, E. Berti, and V. Cardoso, Phys. Rev. D 88, 044018 (2013), eprint 1305.4306.
  * Oshita (2021) N. Oshita, Phys. Rev. D 104, 124032 (2021), eprint 2109.09757.
  * Völkel et al. (2019) S. H. Völkel, R. Konoplya, and K. D. Kokkotas, Phys. Rev. D 99, 104025 (2019), eprint 1902.07611.