Quanto: Optimizing Quantum Circuits with Automatic Generation of Circuit Identities

Jessica Pointing
Department of Physics, University of Oxford, Oxford, OX1 3PU, United Kingdom
Department of Computer Science, Stanford University, Stanford, California, 94305, United States

Oded Padon
VMware Research, Palo Alto, California, 94304, United States
Department of Computer Science, Stanford University, Stanford, California, 94305, United States

Zhihao Jia
Department of Computer Science, Stanford University, Stanford, California, 94305, United States
Department of Computer Science, Carnegie Mellon University, Pittsburgh, Pennsylvania, 15213, United States

Henry Ma
Department of Computer Science, University of California, Los Angeles, California, 90095, United States

Auguste Hirth
Department of Computer Science, University of California, Los Angeles, California, 90095, United States

Jens Palsberg
Department of Computer Science, University of California, Los Angeles, California, 90095, United States

Alex Aiken
Department of Computer Science, Stanford University, Stanford, California, 94305, United States

Existing quantum compilers focus on mapping a logical quantum circuit to a quantum device and its native quantum gates. Only simple circuit identities are used to optimize the quantum circuit during the compilation process. This approach misses more complex circuit identities, which could be used to optimize the quantum circuit further.

We propose Quanto, the first quantum optimizer that automatically generates circuit identities. Quanto takes as input a gate set and generates provably correct circuit identities for the gate set.

Quanto's automatic generation of circuit identities includes single-qubit and two-qubit gates, which leads to a new database of circuit identities, some of which are novel to the best of our knowledge. In addition to the generation of new circuit identities, Quanto's optimizer applies such circuit identities to quantum circuits and finds optimized quantum circuits that have not been discovered by other quantum compilers, including IBM Qiskit and Cambridge Quantum Computing Tket. Quanto's database of circuit identities could be applied to improve existing quantum compilers and Quanto can be used to generate identity databases for new gate sets.

§ INTRODUCTION

Quantum programs are expressed as circuits of primitive gates; the specific gate set available varies depending on the quantum computer. Compilers are used to optimize the circuit in an effort to reduce the execution time, the noise that naturally arises in quantum computations [17], or both.

A quantum algorithm is implemented on a quantum computer with a quantum circuit, which consists of a sequence of quantum gates. A quantum gate applies a transformation to a qubit and can be represented as a matrix. Generally, the larger the circuit depth (length of the longest path from the input to the output), the noisier the computation, leading to more errors in the results. Because of the possible circuit of incorrect results and the inherent statistical nature of quantum computation, quantum programs must be executed many times to obtain an accurate picture of their results. Therefore, reducing the depth of a quantum circuit not only makes a quantum program faster, but also reduces the potential for noise and reduces the number of program evaluations required.  Thus, finding equivalent but smaller depth circuits is a core problem in the implementation of quantum programs.

Quantum hardware imposes many constraints on quantum compilers, including the layout and connectivity of the qubits. Such aspects are addressed by the compilers built by the hardware vendor (such as IBM or Google). Less attention has been given to optimizing the logical circuit itself – finding circuit identities for a given gate set.

Circuit identities have been found manually [11] and automatically for single-qubit gates [14], but as far as we know there is no system today that automatically finds identities for full gate sets and no compiler that optimizes quantum circuits based on such circuit identities.

We present Quanto, the first quantum circuit optimizer that automatically generates circuit identities. Figure <ref> shows an overview of Quanto. In a  first, offline phase, Quanto generates a database of circuit identities. The Quanto compiler takes a quantum circuit and identity database as inputs and uses a cost-based search algorithm to replace the quantum circuit with an optimized quantum circuit.

Overview of Quanto

Figure <ref> shows an example of a quantum circuit with two qubits and a depth of five. The circuit contains the two-qubit gate $CX$ and the single-qubit gates $H$ and $Z$. The single-qubit gates are represented by $2 \times 2$ unitary complex matrices and the two-qubit gates are represented by $4 \times 4$ unitary complex matrices. This quantum circuit is a logical quantum circuit, which is a quantum circuit defined in terms of unitary matrices and is agnostic to the quantum hardware.

$q_0$         1         1

Example of a quantum circuit that can be optimized with Quanto

If we wish to run this logical quantum circuit on the quantum hardware, we would have to compile the quantum circuit into the hardware's native quantum gates. These are the set of quantum gates that have been implemented on particular quantum hardware. For example, IBM's native quantum gate set was [IBM recently changed their native gate set to {CX, ID, RZ, SX, X}] {U1,U2,U3,CX} where

U1(λ) =
1     0

0     e^iλ

U2(ψ,λ) = 1/√(2)
1     -e^iλ

e^iψ     e^i(λ+ψ)

U3(θ,ψ,λ) =
cos(θ/2)     -e^iλ sin(θ/2)

e^iψ sin(θ/2)     e^i(λ+ ψ) cos(θ/2)

CX =
1     0     0     0

0     1     0     0

0     0     0     1

0     0     1     0

In this example, the source gates are the gates present in the logical quantum circuit {X, H, CZ} and the target gates are the native quantum gates {U1, U2, U3, CX}. Therefore, in order to run this logical quantum circuit on the quantum hardware, we must find circuit identities for the source gates in terms of the target gates. In this case $Z = U1(\pi)$ and $H = U2(\pi,0)$.

Another step for running a logical quantum circuit on the hardware is to map the qubit and gate connections to those of the physical constraints of the hardware. Most compilers focus on these two aspects: (1) mapping source gates to target gates and (2) mapping the quantum circuit to the hardware constraints. Then, generally, a set of simple manual rules are used to optimize the logical quantum circuit itself. For example, voqc (Verified Optimizer for Quantum Circuits) [12] includes five rules for optimizations – (0) not propagation, (1) Hadamard reduction, (2) single-qubit gate cancellation, (3) two-qubit gate cancellation, and (4) rotation merging. The voqc optimizer applies these rules in the following order in order to optimize a circuit: 0, 1, 3, 2, 3, 1, 2, 4, 3, 2.

Applying a limited number of optimization rules limits the opportunity to find other optimizations. Instead of defining a set of optimization rules, our compiler is the first to automatically generate circuit identities. Quanto finds an optimized logical quantum circuit that is equivalent to the input logical quantum circuit.

The first phase in the compiler is generating circuit identities. Circuit identities are quantum circuits that are mathematically equivalent. A quantum circuit can be represented as a matrix and two quantum circuits are equivalent when their matrices are equal (for a certain floating point precision).

Quanto's circuit identity generator enumerates all possible quantum circuits over a given set of gates up to a fixed depth and number of qubits. In the example in Figure <ref>, the circuit identity generator would find all possible circuits with two qubits and a depth of five containing the gates {X, H, CZ} (or another gate set specified). To efficiently find all such circuit identities, Quanto constructs a hash table where quantum circuits are stored based on the hash of their unitary matrix.
The operation of a quantum gate can be represented using a unitary complex matrix. The operation of a quantum circuit can be represented by a unitary matrix obtained from the matrices of individual gates using matrix multiplication and tensor products. In the example in Figure <ref>, the unitary matrix is $U = (I \otimes H) \cdot CX \cdot (Z \otimes Z) \cdot CX \cdot (I \otimes H)$.

The second phase is the optimizer. When given an input quantum circuit, our optimizer searches for quantum circuits equivalent to the input quantum circuit that minimise the cost function. In this paper we use circuit depth as the cost; that is, a circuit will be replaced by the lowest-depth equivalent circuit known.  In general, however, the cost function can be chosen to optimize any quantifiable criteria of a circuit. For example, a cost function could include the total number of gates, total number of specific gates, such as $T$ gates or two-qubit gates, and/or gate errors. Because we generate all possible circuits, Quanto returns the optimal quantum circuit with respect to the cost function. In the example in Figure <ref>, our optimizer returns an equivalent quantum circuit shown in Figure <ref>. Quanto has reduced the depth of the circuit from five to one. This circuit has the minimum depth possible for that given gate set.

The optimal circuit of the quantum circuit in Figure <ref> for the given set

The circuit identities are obtained by simulating the quantum circuits. One problem, however, is that simulating a quantum circuit is expensive if the quantum circuit has a large number of qubits and circuit depth. The promise of quantum computing is dependent on the fact that simulating quantum circuits is classically inefficient. Therefore, the larger the quantum circuit becomes (i.e. more qubits and greater depth), and the more gates that are added, the more configurations there are to generate, which eventually can become infeasible. In order to solve this issue in our optimizer, we introduce tiles. A tile is a subset of a quantum circuit with a specified length and width. A tile is the source quantum sub-circuit that is matched to equivalent quantum circuits. In the example in Figure <ref>, if our tile width is equal to three and our tile length is equal to two, then the tiles in Figure <ref> would be generated.

i) scale=1

$q_0$     [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]     1         1

ii) scale=1

$q_0$         1[2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]         1

iii) scale=1

$q_0$         1     [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]     1

Splitting the circuit into tiles of width 3 and height 2

Quanto searches for quantum circuits equivalent to the tiles in the circuit identities database. A cost-based search algorithm is used to apply the substitution, which is a quantum circuit in the circuit identities database that replaces a tile in the quantum circuit. A substitution is the target quantum circuit that defines an equivalent new quantum circuit to replace the matched quantum circuit (i.e. the tile). Figure <ref> show how the circuit gets optimized when it is split into tiles.

Step 1) scale=0.65

$q_0$         1[2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]         1

$q_0$         [2,steps=1,style=dashed,
rounded corners,fill=orange!20, inner sep=4pt,

Step 2) scale=1

$q_0$     [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,

$q_0$     [1,steps=1,style=dashed,
rounded corners,fill=orange!20, inner sep=4pt,

Steps to optimize the circuit when the circuit is split into tiles. The boxes shaded in blue are the tiles. The boxes shaded in orange are the substitutions.

Our main contribution is presenting Quanto, which is the first quantum optimizer to automatically generate quantum circuit identities. Our automatic generator finds identities used by existing compilers, and in addition finds novel identities (to the best of our knowledge). Quanto can generate novel identities for any given gate set, which is especially useful for quantum hardware with native gate sets that are not standard (for example, the two-qubit iSwap gate is less common than the two-qubit CX gate). In our evaluation we demonstrate that the IBM Qiskit compiler [1] can not find these novel identities.

Quanto also includes a tile-based method to optimize a large quantum circuit using the generated identities. We demonstrate that Quanto's optimizer finds equivalent circuits with smaller depth compared to IBM's Qiskit and Cambridge Quantum Computing's Tket [18] compilers.

§ RESULTS

§.§ Generating Existing Circuit Identities

Quanto can automatically discover identities that are known from previous work. [11, 7] We ran the compilers (Quanto, IBM's Qiskit, and Cambridge Quantum Computing's Tket) on a set of quantum circuits with known optimizations. The results show that Quanto is able to find the known optimizations for all of the input quantum circuits, whereas other compilers fail to find some of the identities. The other compilers could find simple identities, such as a circuit with the same gate in a row, but they failed to find more complicated identities.

These results are shown in Tables <ref>, <ref>, <ref>, <ref>, <ref>. The second column of the Tables show the input (uncompiled) circuit and the other columns show the compiled circuit after running on the compiler specified in that column. The last columns in Table <ref> and Table <ref> show the circuits that are generated when the uncompiled circuits are run through the Quanto compiler and then the IBM compiler. The reason we show this column is because our optimizer is at the logical level and does not map the qubits to the hardware, a task we leave to the IBM compiler. We show that for the circuits in Table <ref> and some of the circuits in Table <ref> the circuits are unchanged by the mapping to the hardware. In particular, it is interesting to look at Circuit A and Circuit C in Table <ref>. One might say that the IBM compiler circuit increased the depth of the uncompiled circuit from depth three to six in Circuit A and from depth three to seven in Circuit C because it needs to map to the hardware, but when we run our compiled circuit (which reduced the depth of the circuit to two in both circuits), the IBM compiler left the optimization unchanged. These results demonstrate that optimizing at the logical level is a useful task that can be done before mapping to the specific quantum hardware.

The Quanto compiled circuits of Circuits B and D in Table <ref> are changed when they are mapped to the hardware. In particular, two-qubit gates that are not applied to neighbouring qubits are changed. In order to address this, we implemented a two-qubit gate neighbouring constraint in Quanto that can be turned on or off. If this constraint is on, it only returns optimized circuits that include two-qubit gates on neighbouring qubits. This is a simple hardware constraint. Quanto could be improved by implementing more specific hardware constraints. The last column in Table <ref> shows the compiled circuits when the constraint is on.

Circuit     Uncompiled
circuit     IBM compiled
circuit     Tket compiled
circuit     Quanto compiled
circuit     Quanto-IBM
compiled circuit

Depth     2     0     0     0     0

$q_0$         T^†

Depth     2     0     0     0     0

$q_0$         1

$q_0$     U_2(0,π)     1     U_2(0,π)

$q_1$     U_2(0,π)          U_2(0,π)

$q_0$     U_3(π/2,0,π)     1     U_3(π/2,0,π)

$q_1$     U_3(π/2,0,π)          U_3(π/2,0,π)

$q_1$     -1

$q_1$     -1

Depth     3     3     3     1     1

$q_0$         Z     1

$q_0$     U_2(π,π)     1     U_2(0,π)

$q_1$     U_2(0,π)          U_2(0,π)

$q_0$     U_3(3π/2,0,0)     1     U_3(π/2,0,π)

$q_1$     U_3(π/2,0,π)          U_3(π/2,0,π)

$q_1$         -1

$q_0$     U_3(π,0,π)

$q_1$         -1

Depth     4     3     3     2     2

Comparing Quanto to other compilers for optimizing quantum circuits with known optimizations. IBM and Tket backend is ibmq essex.

Circuit     Uncompiled
circuit     IBM compiled
circuit     Tket compiled
circuit     Quanto compiled

$q_0$     1         1

$q_0$     1     U3(π,0,π)     1

$q_0$     1     U3(π,0,π)     1

Depth     3     3     3     1

$q_0$     1         1

$q_0$     1         1

$q_1$          U3(π,0,π)

$q_1$     U3(π,0,π)

Depth     3     3     1     1

$q_0$     1         1

$q_0$     1     U1(π)     1

$q_0$     U1(π)

Depth     3     3     1     1

$q_0$     1         1

$q_0$     1         1

$q_1$          U1(π)

$q_0$     1         1

$q_1$          U1(π)

Depth     3     3     3     1

$q_0$     1         1

$q_0$         1     U3(π,0,π)     1

$q_1$     U2(0,π)                   U2(0,π)

$q_0$         1     U3(π,0,π)     1

$q_1$     U3(π/2,0,π)                   U3(π/2,0,π)

Depth     3     5     5     1

$q_0$     1         1

$q_0$         1         1

$q_1$     U2(0,π)          U1(3π)          U2(0,π)

$q_0$         1         1

$q_1$     U3(π/2,0,π)          U1(π)            U3(π/2,0,π)

Depth     3     5     5     1

Comparing Quanto to other compilers for optimizing quantum circuits with known optimizations. IBM and Tket backend is Aerbackend.

Circuit     Uncompiled
circuit     IBM compiled
circuit     Tket compiled
circuit     Quanto compiled

$q_0$     1         1

$q_0$         1     U1(π)     1

$q_1$     U2(0,π)                   U2(0,π)

$q_0$     U1(π)

Depth     3     3     1     1

$q_0$     1         1

$q_0$         1         1

$q_1$     U2(0,π)          U3(-π,π/2,3π/2)          U2(0,π)

$q_0$         1         1

$q_1$     U3(π/2,0,π)          U3(3π,0,π)          U3(π/2,0,π)

Depth     3     5     5     1

$q_0$         1         1

$q_0$         1         1

$q_1$     U2(0,π)          U3(π,0,π)          U2(0,π)

$q_1$     U1(π)

Depth     5     5     1     1

$q_0$         1         1

$q_0$         1     U1(π)     1

$q_1$     U2(0,π)          U1(π)          U2(0,π)

$q_0$     U1(π)     1         1

$q_1$     U3(π/2,0,π)          U1(π)          U3(π/2,0,π)

Depth     5     5     5     1

Comparing Quanto to other compilers for optimizing quantum circuits with known optimizations. IBM and Tket backend is Aerbackend.

Circuit     Uncompiled
circuit     IBM compiled
circuit     Tket compiled
circuit     Quanto compiled
circuit     Quanto-IBM
compiled circuit     Quanto compiled
(with constraint)

$q_0$         1     2

$q_1$     1

$q_0$              1          1

$q_1$     -1               -1

$q_2$         -1                 -1

$q_0$         1     2

$q_1$     1

$q_0$     1

$q_1$          1

$q_0$     1

$q_1$          1

$q_0$     1

$q_1$          1

Depth     3     6     3     2     2     2

$q_0$     1         1

$q_1$          1          1

$q_0$     1         1

$q_1$          1          1

$q_0$     1         1

$q_1$          1          1

$q_0$     2

$q_1$     1            1      -1

$q_2$           -1

$q_0$     1         1

$q_1$          1          1

Depth     4     4     4     1     4     4

$q_0$     2         1

$q_1$         1

$q_0$     1             1                1

$q_1$            -1           1     -1           1

$q_0$     2         1

$q_1$         1

$q_0$     1

$q_1$          1

$q_0$     1

$q_1$          1

$q_0$     1

$q_1$          1

Depth     3     7     3     2     2     2

$q_0$         1

$q_1$     1          1

$q_0$         1

$q_1$     1          1

$q_0$         1

$q_1$     1          1

$q_0$     1     2

$q_1$          1             1      -1

$q_2$     -1           -1

$q_0$         1

$q_1$     1          1

Depth     3     3     3     2     5     3

$q_0$         1     2             1     2             1     2

$q_1$     1                  1                  1

$q_3$         1     2             1     2             1     2

$q_4$     1                  1                  1

Circuit too large
for displaying image

$q_0$         1     2             1     2             1     2

$q_1$     1                  1                  1

$q_3$         1     2             1     2             1     2

$q_4$     1                  1                  1

$q_0$         1

$q_1$     1

$q_3$         1

$q_4$     1

$q_0$         1

$q_1$     1

$q_3$         1

$q_4$     1

$q_0$         1

$q_1$     1

$q_3$         1

$q_4$     1

Depth     9     31     9     2     2     2

Comparing Quanto to other compilers for optimizing quantum circuits with known optimizations. IBM backend is ibmq lima and Tket backend is Aerbackend.

§.§ Generating Novel Circuit Identities

Quanto can also generate new identities, especially when given custom gates, as it is more likely that identities for custom gates have not been previously discovered. Existing compilers have a fixed repertoire of identities using common gates. The advantage of Quanto is that it can find identities for any custom gate, leading to optimizations that may be missed by other compilers. In Table <ref> we show some of these novel identities (these identities were found from our identity database and we are not aware of any previous reports of these identities). We compare our novel identity optimizations to the IBM Qiskit compiler, as it allows us to specify basis gates for compiling. IBM's compiler is unable to find these circuit identities.

§.§ Optimizing Quantum Circuits

Other than generating existing and novel identities, Quanto can also optimize any quantum circuit as shown in Table <ref>.

Circuit     Uncompiled
circuit     IBM compiled
circuit     Quanto compiled

$q_0$         1

$q_1$     S^†

$q_0$         1

$q_1$     U1(-π/2)          U1(π/2)

$q_0$     1

$q_1$     Y

Depth     3     3     1

$q_0$         1

$q_0$     1

$q_0$     1

Depth     3     1     1

$q_0$             1

$q_1$                  -1

$q_0$     U2(0,3π/2)     1

$q_1$     U1(π/2)          -1     U2(0,π)

$q_0$     [2]iSWAP

Depth     5     4     1

Comparing Quanto to other compilers for optimizing quantum circuits for unique gate sets

Circuit     Uncompiled
circuit     IBM compiled
circuit     Quanto compiled

$q_0$              1     S^†

$q_1$         -1

$q_0$              1     U1(-π/2)

$q_1$     U1(π/2)     -1          U1(π/2)

$q_0$         1          1

$q_1$              -1

Depth     4     4     3

$q_0$         1

$q_1$                  -1

$q_0$     U_3(π,π/2,π/2)     1     U_3(π,0,π)

$q_1$              U_3(π,π/2,π/2)     -1

$q_0$     1

$q_1$          -1

Depth     4     4     3

$q_0$         [2]iSWAP

$q_0$     U1(π)     [2]iSWAP

$q_1$     U3(π,0,π)         U2(0,π)

$q_0$     [2]iSWAP

Depth     3     3     2

$q_0$         1     [2]iSWAP

$q_1$     S^†

$q_0$         1                         1

$q_1$     S^†                       -1

$q_0$         1                 U2(0,3π/2)     1

$q_1$     U1(-π/2)          U1(π/2)          -1     U2(0,π)

$q_0$          1

$q_1$     -1

Depth     7     6     2

Comparing Quanto to other compilers for optimizing quantum circuits with unknown optimizations

§ DISCUSSION

§.§ Related Work

§.§.§ Automatically generating identities

To the best of our knowledge, previous work on automatically discovering quantum circuit identities is limited to single-qubit gates [14]. In addition, these automatically generated identities have not been applied to optimizing quantum programs.

§.§.§ Superoptimization

Superoptimization is a compiler technique that was originally designed to find the optimal machine code for a sequence of instructions [15] [5] [4]. Our approach is inspired by TASO [13]. TASO applies superoptimization to deep learning by optimizing tensor computation graphs. For a given set of tensor operators (e.g., addition, matrix convolution) TASO automatically finds identities and uses these as substitutions to optimize computation graphs of deep learning applications. There are some differences: In the TASO paper, they generate a graph and then go through the list of substitutions and try to see whether the substitution is a sub-graph in the original graph. They have a limited number of substitutions, in the thousands, so they can do that fairly efficiently. In our case, however, we have can have hundreds of thousands or millions of substitutions and the time complexity of that technique would grow large. Instead we split our circuit into sub-circuits, tiles, and then look for the hash of that tile in the hash table, which is a constant lookup. Therefore, our technique is more efficient for our case and only depends on the number of tiles in the circuit, which grows modestly.

§.§.§ Quantum Optimization

There are other compilers, [1, 18, 12] but their focus is mapping the logical circuit to the quantum device. Therefore, they do not pay much attention to the potential logical transformations. This can be shown in our results – other optimizers fail to recognize certain circuit identities. These optimizers have manual rules for logical substitutions and therefore are not finding the full range of optimization.

§.§ Limitations and Future Work

One limitation of Quanto is the scalability of generating a large number of circuit identities, which limits the size (depth and number of qubits) of the circuit identities. This issue could be mitigated by applying pruning techniques during the generation of substitutions to remove redundant substitutions or by using heuristics to produce substitutions. Another limitation is missing out on potential substitutions, due to invalid tiles in the quantum circuit as shown in Figure <ref>.

To use gates with parameters in our optimizer, we can set the parameter to a small set of specific values when generating identities. Since many parameters are used for qubit rotations, and certain angles of rotation are very common in quantum algorithms, this approach gives some reasonable coverage of potential optimizations involving parameters.  However, our optimizer would be improved if it included symbolic identities, where substitutions for parametrized gates are generated regardless of the value of the parameter.

§.§ Conclusion

Quanto is the first quantum optimizer that automatically generates general circuit identities. Quanto finds novel optimizations and finds better optimizations for quantum circuits, compared to other compilers. This work contributes to reducing the circuit depth of quantum circuits, which could lead to less noise and more accurate execution of quantum algorithms.

§ METHODS

§.§ Circuit Identity Generator

This section describes the Quanto circuit identity generator, which automatically generates circuit identities given a gate set, as well as the number of qubits and depth bounds. This generation algorithm finds all circuit identities given the number of qubits, circuit depth, and gate set.

§.§.§ Generation Algorithm

Table <ref> shows the inputs to the generation algorithm.

Symbol     Definition     Description

n     number of qubits     maximum number of qubits of the circuit identities generated

d     depth of the circuit     the maximum depth of the circuit identities generated

dp     decimal point accuracy     the decimal points of the  floating-point  numbers  in  the  matrices  of  the gates

gs     gate set     a list of gates to use in the circuit identities, including the identity gate

Inputs of the generation algorithm

Step 1: Represent the quantum circuit as a grid.
We represent a quantum circuit as a grid, where the number of rows $m$ is the number of layers in the quantum circuit, also known as the circuit depth. The number of columns $n$ is the number of qubits in the quantum circuit. If there is no gate present at a particular row for a particular column, we place the Identity gate in that position.

Step 2: For a single layer (depth = 1), generate all the possible circuits of quantum circuits using the gate set.

a) Single-qubit gates only

If there are only single-qubit gates in the gate set, then each gate is applied to a single qubit. If $n$ is the number of qubits and $g$ is the number of single-qubit gates, then there are $g^n$ possible circuits.

Example: Gate Set = [I, X, H] and $n = 2$

I = \begin{pmatrix}
1 & 0 \\
0 & 1
\end{pmatrix},
X = \begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix},
H = \begin{pmatrix}
\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\
\frac{1}{\sqrt{2}} & -\frac{1}{\sqrt{2}}
\end{pmatrix}

Table <ref> shows all the possible circuits that would be generated for the gate set $[I, X, H]$

$q_0$     I

$q_1$     I

[I, I]     [X, X]     [H, H]

$q_0$     I

$q_0$     I

$q_1$     I

[I, X]     [I, H]     [X, I]

$q_1$     I

[X, H]     [H, I]     [H, X]

possible circuits for single-qubit gates given gate set $[I, X, H]$
. In the case of the possible circuit $[H, X]$ the $H$ gate is applied to the first qubit and the $X$ gate is applied to the second qubit as shown in Table <ref>.

b) Two-qubit and single-qubit gates

If there are two-qubit gates, then we split the two-qubit gate into the first and second qubits it operates on.

Example: Gate Set = [I, CX] and $n=3$

CX = \begin{pmatrix}
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 0 & 1\\
0 & 0 & 1 & 0
\end{pmatrix}

In the notation below, CXC denotes applying the CX gate to the $control$ qubit and CXT denotes applying the CX gate to the $target$ qubit.

$q_0$     I

$q_1$     I

$q_2$     I

$q_0$     1

$q_2$     I

$q_0$     2

$q_1$     I

$q_1$     -1

$q_2$     I

[I, I, I]     [CXC, CXT, I]     [CXC, I, CXT]     [CXT, CXC, I]

$q_0$     I

$q_1$     1

$q_1$     I

$q_2$     -2

$q_0$     I

$q_2$     -1

[I, CXC, CXT]     [CXT, I, CXC]     [I, CXT, CXC]

Possible circuits for two-qubit gates given gate set $[I, CX]$

Table <ref> shows all the possible circuits that would be generated for the gate set $[I, CX]$. In the case of the possible circuit $[CXC, CXT, I]$ the $CX$ gate is applied to the first and second qubit where the first qubit is the $control$ qubit and the second qubit is the $target$ qubit. The $I$ gate is applied to the third qubit.

Verifying the validity of the two-qubit gates

$CXC$ and $CXT$ must both be present in the possible circuits; it is not possible to have one present without the other. For example, $[CXC, I, I]$ would be an invalid configuration. Quanto verifies that $CXC$ and $CXT$ are present in each possible circuit.

Step 3: For multiple layers (depth $> 1$), generate all the possible circuits of the quantum circuits using the gate set. For an additional layer in the circuit, we add all the possible circuits for a single layer onto the previous layer.

Example: For the gate set presented in Table <ref>, we take each possible circuit and append it to the end of each possible circuit. Table <ref> illustrates appending the possible circuit $[X, X]$ to each possible circuit. This is repeated for $[H, H], [I, X], etc...$

$q_0$     I

$q_1$     I

$[I, I], [X, X]$     $[X, X], [X, X]$     $[H, H], [X, X]$

$q_0$     I

$q_0$     I

$q_1$     I

$[I, X], [X, X]$     $[I, H], [X, X]$     $[X, I], [X, X]$

$q_1$     I

$[X, H], [X, X]$     $[H, I], [X, X]$     $[H, X], [X, X]$

Possible circuits for multiple layers

Step 4: Calculate the unitary matrices and generate their fingerprints. For each possible circuit, we find the unitary matrix of the possible circuit. We adopt an idea from compiler superoptimization [4] and compute a fingerprint for each possible circuit, which is a hash of the unitary matrix. Figure <ref> shows a possible circuit and how to calculate the unitary matrix. The matrices in part (c) correspond to the gates in part (b) and (a) of Figure <ref>.

\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \circuitX & \qw \\
\lstick{$q_1$} & \circuitX & \circuitX & \qw \\
\end{tikzcd}

$$ (H \otimes X) \times (X \otimes X) $$

\begin{equation*}
\resizebox{1\hsize}{!}{$%
\Bigg(\frac{1}{\sqrt{2}}
\begin{pmatrix}
1 & 1 \\
1 & -1
\end{pmatrix} \otimes
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix}
\Bigg)
\times
\Bigg(
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix} \otimes
\begin{pmatrix}
0 & 1 \\
1 & 0
\end{pmatrix} \Bigg) \\
= \frac{1}{\sqrt{2}} \begin{pmatrix}
1 & 0 & 1 & 0 \\
0 & 1 & 0 & 1 \\
-1 & 0 & 1 & 0\\
0 & -1 & 0 & 1
\end{pmatrix}
\end{equation*}

(a) A possible circuit generated from a gate set (b) the gate representation of the circuit (c) the matrix calculation of the circuit to obtain the unitary matrix

\adjustbox{scale=1}{%
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \circuitX & \ctrl{1} & \qw \\
\lstick{$q_1$} & \gate{I} & \circuitX & \targ{} & \qw \\
\end{tikzcd}

$$[[H, I], [X, X], [CXC, CXT]]$$

(a) Quantum circuit which can be represented as (b) a matrix, which can be represented computationally as (c)

Step 5: Create hash table of fingerprints.
The hash table of fingerprints includes the quantum circuit as the key and the fingerprint of the quantum circuit as the value. This allows us to look up a quantum circuit's fingerprint in constant time during the optimizer stage, without having to recalculate the unitary matrix and the fingerprint of the quantum circuit. Figure <ref> shows the sequence of steps to generate a hash table of fingerprints for each possible circuit. Table <ref> is an example hash table with one entry which shows the quantum circuit representation as the key and its corresponding fingerprint as the value.

A quantum circuit can be represented by an unitary matrix. The hash of this unitary matrix is its fingerprint.

Key     Value

$[[H, X], [X, X]]$     $x7ffx9exa0xe6?x00…$

Example of a hash table with one entry demonstrating the quantum circuit representation as the key and the fingerprint as the value. Note that the fingerprint in this table is only for illustrative purposes (as the actual fingerprint is too long to be displayed).

Step 6: Create hash table of circuit identities.
The hash table of substitutions includes the fingerprints as keys and an array of matching circuit identities as values. For each possible circuit we generate, we check if its fingerprint is in the hash table of circuit identities. If it is, we add the possible circuit to the matching fingerprint. If it is not, we add the fingerprint to the hash table and add the possible circuit to the array of matching circuit identities. The identities are stored implicitly in the hashtable - any set of matrices associated with the same key are equal and can be substituted for each other. The identities are up to some (configurable) floating point precision. Table <ref> is an example hash table with one entry which shows the fingerprint as the key and all the circuit identities that match to the fingerprint as the value.

Key     Value

$x7ffx9exa0xe6?x00…$     $[[X, X], [X, X]]$

$[[H, H], [H, H]]$

Hash table with one entry demonstrating the fingerprint as the key and the array of substitutions as the value.

§.§ Scaling

To estimate the computational resources to generate circuit identities, we derived a formula for the total number of possible circuits to be generated. The symbols are defined in Table <ref>.

\begin{equation}
S_{l} =\sum_{r=0}^{\lfloor n/2 \rfloor}\frac{n!}{r! (n-2r)!} g^{n-2r} t^r
\end{equation}

\begin{equation}
S = S_{l}^{d}
\end{equation}

Symbol     Definition

$S$     total number of possible circuits

$S_{l}$     number of possible circuits per depth

$g$     number of single-qubit gates

$n$     number of qubits

$t$     number of two-qubit gates

$S$     number of possible circuits for the depth

$d$     depth

Symbols for scaling formula

As one can see, the number of possible circuits scales quickly with $n$ which is to be expected. We always use the optimizer, and if the size of the tile is equivalent to the size of the input circuit, there is only one tile, meaning we enumerated this circuit during generation.

§.§ Optimizer

This section describes the Quanto optimizer, which applies the automatically generated circuit identities to optimize a quantum circuit.

§.§.§ Algorithm

Table <ref> shows the inputs.

Symbol     Definition

quantum circuit     quantum circuit to optimize

$i$     the maximum number of qubits of the tiles

$j$     the maximum depth of the tiles

Possible circuits for multiple layers

Step 1: Generate tiles of the input quantum circuit.
If the input circuit, which has depth (or width) $m$ and number of qubits (or length) $n$, is larger than the size of the circuit identities generated, we split a circuit into sub-circuits, called $tiles$. The tile width, $j$, is the maximum depth of the circuit identities and the tile length, $i$, is the maximum number of qubits of the circuit identities. We iterate through the rows and columns in the quantum circuit matrix and generate the indices of all possible tiles. The number of possible $i \times j$ tiles for an $n \times m$ circuit is $(n-i+1) \cdot (m-j+1)$.

Figure <ref> shows the dimensions of the tiles and quantum circuit. Figure <ref> shows all the possible tiles generated if we chose a tile with $j=2$ and $i=2$ in a quantum circuit with $m=3$ and $n=2$. Using the formula for the number of possible tiles, we see that there are four possible tiles, which corresponds with Figure <ref>.

Dimensions of a tile compared to the dimensions of a quantum circuit

$q_0$     [2,steps=2,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]         I

$q_2$             I

$q_0$         [2,steps=2,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]     I

$q_2$             I

$q_0$             I

$q_1$     [2,steps=2,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,

$q_2$             I

$q_0$              I

$q_1$         [2,steps=2,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,

$q_2$             I

All the possible tiles (shaded in blue) generated for a tile width and length = 2 for a quantum circuit with width and length = 3

Step 2: Check for valid tiles.
Not all tiles generated from the circuit will be valid. If a two-qubit gate is not placed at the boundaries of the tile, then the tile is invalid if the two-qubit gate is cut off by the tile. In this case, we discard such tiles. We can keep a tile if the two-qubit gate is placed at the boundary of the tile even it has been cut off by the tile. Figure <ref> shows an example of an invalid tile, whereas Figure <ref> shows an example of a valid tile.

$q_0$         1     I

$q_1$     [2,steps=3,style=dashed,
rounded corners,fill=red!20, inner sep=4pt,

$q_2$              I     H

An example of an invalid tile (shaded in red), because the $CX$ gate is cut off and is not first or last on its wire".

$q_0$         1     I

$q_1$          [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,

$q_2$              I     H

An example of a valid tile (shaded in blue), because even though the $CX$ gate is cut off, it is first or last on its wire.

Step 3: Search for tile in hash table of circuit identities.
For each tile we have generated, we first search in the hash table of fingerprints.  If a tile contains a two-qubit gate that has been cut off only at the boundary of the tile, then we first have to replace the two-qubit gate with an Identity gate in order to make the tile valid and searchable in the hash table of circuit identities. Part (a) in Figure <ref> shows this transformation. The now-valid tile is represented as a matrix, as shown in Part (b) in Figure <ref> and we can search for this matrix in the hash table of fingerprints. This gives us the fingerprint of the tile, without needing to calculate the unitary matrix of the tile. After we have obtained the fingerprint, we can search for circuit identities in the hash table of circuit identities as shown in Part (c) in Figure <ref>. Part (d) in Figure <ref> shows circuit identities for the tile.

Once we have the fingerprint, we then look up the possible substitutions by its fingerprint in the hash table of circuit identities.

\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \targ{} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \circuitX & \circuitX & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \circuitH & \circuitH & \qw
\end{tikzcd}
}  \longrightarrow
\adjustbox{scale=0.8}{%
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \gate{I} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \circuitX & \circuitX & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \circuitH & \circuitH & \qw
\end{tikzcd}

I & Y \\
X & H \\
X & H
\end{pmatrix}
\longrightarrow \begin{tabular}{ |c | c | }
\hline
\thead{Key} & \thead{Value} \\
\hline
$[[I, Y], [X, H], [X, H]]$ & $x7ffx9exa0xe6?x00…$ \\
\hline
\end{tabular}

$$\begin{tabular}{ |c | l | }
\hline
\thead{Key} & \thead{Value} \\
\hline
$x7ffx9exa0xe6?x00…$ & \makecell{$[[I, Y], [I, H], [I, H]]$ \\
$[[I, Y], [X, I], [X, I]]$ \\
$[[I, Y], [I, I], [I, I]]$}
\\
\hline
\end{tabular}$$

\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \gate{I} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY & \circuitH & \circuitH & \qw
\end{tikzcd}
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_2$} & \circuitY & \gate{I} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \circuitX & \circuitX & \qw \\
\lstick{$q_1$} & \circuitS & \circuitY & \gate{I} & \gate{I} & \qw
\end{tikzcd}
} \adjustbox{scale=0.6}{%
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \gate{I} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY & \gate{I} & \gate{I} & \qw
\end{tikzcd}

(a) the invalid $CX$ gate in the tile is replaced with an Identity gate to make it a valid tile (b) the matrix representation of the tile is searched for in the hash table of fingerprints, which returns its fingerprint (c) the fingerprint is used to find the circuit identities for the tile in the hash table of circuit identities (d) the circuit identities for the tile

Step 4: Apply substitutions via cost-based search.
We then iterate through the circuit identities for each tile and evaluate the cost of each circuit identity. In this paper, we simply use the depth of the circuit (ignoring identity gates) as the circuit's cost; the cost function can be customized if desired. Table <ref> shows the cost for the circuit identities presented in Figure <ref>. In this case, circuit (iii) in Table <ref> has the lowest cost and will be used to substitute the tile in the quantum circuit.

Circuit Identity     Cost

$q_0$         1     I     I

$q_1$         I [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]     I     I

$q_0$         1     I     I

$q_2$         I [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,

$q_1$             I     I

$q_0$         1     I     I

$q_1$         I [2,steps=3,style=dashed,
rounded corners,fill=blue!20, inner sep=4pt,
background]     I     I

$q_2$             I     I

Circuit identities and their cost

Step 5: Check for valid substitutions.
Adding a substitution might make the entire quantum circuit invalid. For example, two-qubit gates may be cut off or not aligned properly. Additional steps, explained below, are needed to ensure the final circuit has valid two-qubit gates.

Step 6: Apply the substitution to the circuit.
We apply the substitution to the quantum circuit. If the original tile contained an invalid $CX$, then we have to do a couple more transformations before applying the substitution. Figure <ref> shows the steps. Step (a) is equivalent to step (a) in Figure <ref>. In step (b), we apply the substitution with the lowest cost, which was the circuit identity (iii) in Table <ref>. We now need to transform this new circuit so that it is consistent with the gates in the original circuit. In particular, this new circuit has the $CX$ target gate replaced with an Identity gate. Step (b) replaces the Identity gate with the original $CX$ target gate. Step (c) shows the total transformation from the original circuit to an optimized circuit, after only one tile substitution. The Identity gates have been omitted to show that the depth of the circuit has decreased and therefore the circuit has been improved. We then continue to the next tile with the newly replaced substitution.

\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \targ{} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \circuitX & \circuitX & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \circuitH & \circuitH & \qw
\end{tikzcd}
}  \longrightarrow
\adjustbox{scale=0.8}{%
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \gate{I} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \circuitX & \circuitX & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \circuitH & \circuitH & \qw
\end{tikzcd}

\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \gate{I} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \gate{I} & \gate{I} & \qw
\end{tikzcd}
}  \longrightarrow
\adjustbox{scale=0.8}{%
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \targ{} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \gate{I} & \gate{I} & \qw
\end{tikzcd}
} $$

\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \gate{I} & \gate{I} & \qw \\
\lstick{$q_1$} & \circuitY & \targ{} \gategroup[2,steps=3,style={dashed,
rounded corners,fill=blue!20, inner sep=4pt},
background]{{\sc}} & \circuitX & \circuitX & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \circuitH & \circuitH & \qw
\end{tikzcd}
}  \longrightarrow \adjustbox{scale=0.8}{%
\begin{tikzcd}
\lstick{$q_0$} & \circuitH & \ctrl{1} & \qw \\
\lstick{$q_1$} & \circuitY & \targ{} & \qw \\
\lstick{$q_2$} & \circuitS & \circuitY{} & \qw
\end{tikzcd}

Applying the substitution

Step 7: Repeat process.
Steps 1 - 6 are repeated for a fixed number of times. In our implementation, ten times were sufficient.

§.§ Implementation

Quanto is implemented in Python and can be used as an API. In particular, Quanto has the following modules:

§.§.§ Circuit module

Quanto can take qasm (quantum assembly language) files [9] as input, [19, 6, 8, 3, 10], which is a popular way to construct circuits and makes our optimizer compatible with other quantum circuit frameworks, such as Qiskit [1] and Cirq [2], which can input/output qasm files. Our code converts the qasm file into our object representation of a circuit as shown in Figure <ref>. We represent a circuit as a grid, where the rows represent the layers and the columns represent the gates applied to the qubits.

§.§.§ Subs module

The user can specify a gate set to generate circuit identities. Quanto also has several pre-set gate set options apart from customized gate sets. If no gate set is specified, then our code automatically detects all the gates used in the input circuit and uses that as the gate set. This module generates the circuit identities and builds the circuit identities database is used for the optimizer module.

§.§.§ Optimize module

The optimize module takes the generated circuit identities database and optimizes the input circuit. It outputs an optimized circuit, along with metrics, such as the final circuit depth. The circuit can be saved as a qasm file or a qiskit circuit (applicable for certain gates).

[1]
<https://qiskit.org/>, 2017.

[2]
<https://quantumai.google/cirq>, 2018.

[3]
S. Balensiefer, Lucas Kreger-Stickles, and M. Oskin.
Quale: quantum architecture layout evaluator.

[4]
Sorav Bansal and Alex Aiken.
Automatic Generation of Peephole Superoptimizers.

[5]
Sorav Bansal and Alex Aiken.
Binary translation using peephole superoptimizers.
In Proceedings of the 8th USENIX Conference on Operating Systems
Design and Implementation, OSDI'08, page 177–192, USA, 2008. USENIX

[6]
I. Chuang.
<http://www.media.mit.edu/quanta/qasm2circ/>, 2005.

[7]
Gavin E. Crooks.
Gates, states, and circuits.
<https://threeplusone.com/pubs/on_gates.pdf>, 2021.

[8]
A. Cross.
<http://www.media.mit.edu/quanta/quanta-web/projects>, 2005.

[9]
A. Cross, L. Bishop, J. Smolin, and J. Gambetta.
Open quantum assembly language.
arXiv: Quantum Physics, 2017.

[10]
M. Dousti, A. Shafaei, and Massoud Pedram.
Squash 2: a hierarchical scalable quantum mapper considering ancilla
Quantum Inf. Comput., 16:332–356, 2016.

[11]
Juan Carlos Garcia-Escartin and Pedro Chamorro-Posada.
Equivalent Quantum Circuits.
arXiv e-prints, page arXiv:1110.2998, October 2011.

[12]
Kesha Hietala, Robert Rand, Shih-Han Hung, Xiaodi Wu, and Michael Hicks.
A verified optimizer for quantum circuits.
Proc. ACM Program. Lang., 5(POPL):1–29, 2021.

[13]
Zhihao Jia, Oded Padon, James J. Thomas, Todd Warszawski, M. Zaharia, and
A. Aiken.
Taso: optimizing deep learning computation with automatic generation
of graph substitutions.
Proceedings of the 27th ACM Symposium on Operating Systems
Principles, 2019.

[14]
Chris Lomont.
Quantum Circuit Identities.
arXiv e-prints, pages quant–ph/0307111, July 2003.

[15]
Henry Massalin.
Superoptimizer: a look at the smallest program.
ACM SIGARCH Computer Architecture News, 15(5):122–126,
November 1987.

[16]
IBM recently changed their native gate set to {CX, ID, RZ, SX, X}.

[17]
John Preskill.
Quantum Computing in the NISQ era and beyond.
2019:R62.001, January 2019.

[18]
Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Edgington, and Ross
Tket: A retargetable compiler for nisq devices.
03 2020.

[19]
K. Svore, A. Aho, A. Cross, I. Chuang, and I. Markov.
A layered software architecture for quantum computing design tools.
Computer, 39:74–83, 2006.


# Restoring Strong Cosmic Censorship in Reissner-Nordström-de Sitter Black
Holes via Non-Minimal Electromagnetic-Scalar Couplings

Jie Jiang2<EMAIL_ADDRESS>Jia Tan1 Corresponding author.
<EMAIL_ADDRESS>1Jiangsu Key Laboratory of Micro and Nano Heat Fluid Flow
Technology and Energy Application, School of Physical Science and Technology,
Suzhou University of Science and Technology, Suzhou, 215009, China  2College
of Education for the Future, Beijing Normal University, Zhuhai, 519087, China

###### Abstract

We investigate whether the Strong Cosmic Censorship (SCC) Conjecture can be
reinstated in Reissner-Nordström-de Sitter (RNdS) black holes by introducing
non-minimal couplings between the electromagnetic and scalar fields in
Einstein-Maxwell-scalar (EMS) theory. By conducting numerical calculations, we
find that the SCC can be restored within a specific range of the coupling
constant. Notably, for a given value of the cosmological constant, there
exists a critical coupling constant above which the EMS theory satisfies the
SCC. These findings suggest that the non-minimal couplings between the
electromagnetic and scalar fields may play a crucial role in the restoration
of the SCC in RNdS spacetime.

##  I Introduction

The Strong Cosmic Censorship Conjecture (SCC) represents a critical and
fascinating issue in classical general relativity (GR). First proposed by
Roger Penrose in 1969 Penrose:1969pc , the conjecture asserts that the Cauchy
horizon (CH) of a black hole must be singular and, consequently, inextendible
into the future. Based on the principle of determinism in classical GR, the
SCC contends that our universe’s future is dictated by initial conditions
specified on a Cauchy surface. A non-singular CH would result in an
undetermined future and a non-deterministic spacetime, which is why the SCC is
regarded as a foundational principle in classical GR.

Though there are multiple formulations of the SCC, the most contemporary
version was introduced by Demetrios Christodoulou in 2008 Christodoulou:2008nj
. This formulation posits that the spacetime metric cannot extend beyond the
CH in the form of weak solutions to the field equations. In essence, the CH
must be a mass-inflation singularity, signifying that field perturbations must
evolve to become divergent on the CH. As a highly complex issue, proving or
disproving the SCC remains one of the most vital and demanding challenges in
the realm of GR.

The validity of the Strong Cosmic Censorship Conjecture (SCC) has been
confirmed in specific black hole solutions, including asymptotically flat
Reissner-Nordström and Kerr black holes Simpson:1973ua ; Poisson:1990eh ;
Dafermos:2003wr . However, recent findings have identified violations of the
SCC in the nearly extremal region of the Reissner-Nordström-de Sitter (RNdS)
black hole Cardoso:2017soq . These violations are attributed to the presence
of cosmological horizons, which can suppress the blueshift effect responsible
for the mass inflation singularity at the Cauchy horizon (CH) Brady:1996za ;
Costa:2014aia ; Hintz:2015jkj ; Costa:2016afl . Following this, further
violations have been discovered in the RNdS background by massless charged
scalar fields, massless Dirac fields, and others Mo:2018nnu ; Ge:2018vjq ;
Liu:2019rbq ; Dias:2018ufh ; Cardoso:2018nvb ; Destounis:2018qnb ; Liu:2019lon
; Destounis:2020yav . Hollands et al. Hollands:2019whz  made progress in
understanding the SCC in RNdS black holes, showing that non-singular initial
quantum field data can restore the SCC. The review Kontou:2020bta  provides a
comprehensive discussion on energy conditions in GR and quantum field theory.
Moreover, investigations into the validity of the SCC for the RNdS black hole
under classical perturbations of non-minimally coupled scalar fields in
certain modified gravitational theories have been conducted Destounis:2019omd
; Guo:2019tjy . Within the framework of the Einstein-Maxwell-Scalar-Gauss-
Bonnet (EMSG) theory, Ref. Sang:2022frw  has demonstrated that the SCC can be
restored for RNdS black holes under specific parameter ranges when
incorporating non-minimally coupled scalar field perturbations. This finding
offers a fresh perspective for studying the SCC in modified gravitational
theories. The understanding of the SCC continues to be an active research
area, bearing significant implications for the fundamental principles of
classical GR and the ultimate fate of our universe.

The Einstein-Maxwell-Scalar (EMS) theory has attracted significant interest
due to its non-minimal couplings between the electromagnetic and scalar fields
and the resulting black hole solutions. These couplings have long been
considered in the context of theories such as Kaluza-Klein and supergravity
Gibbons:1987ps ; Garfinkle:1990qj , and more recently, in the context of black
hole spontaneous scalarization, a strong gravity phase transition
Damour:1993hw ; Herdeiro:2018wub . EMS models use a non-minimal coupling
between the scalar field and Maxwell invariant to induce scalarization,
requiring the presence of electric or magnetic charge. Studying the EMS models
has led to the discovery that spontaneous scalarization of charged black holes
occurs dynamically, resulting in scalarized, perturbatively stable black holes
Herdeiro:2018wub ; Fernandes:2019rez ; Myung:2018jvi ; Myung:2018vug ;
Myung:2019oua .

Given the similarities between EMS and ESGB, it is worth exploring whether SCC
can also be restored in EMS. Studying SCC in EMS is highly significant as it
can provide valuable insights into the interplay between non-minimal
couplings, spontaneous scalarization, and the fundamental principles of
classical GR Brady:1996za ; Cardoso:2017soq . Understanding whether SCC can be
restored in EMS may ultimately deepen our knowledge of black holes and the
underlying principles of gravitational theories Cardoso:2018nvb ;
Destounis:2018qnb .

The structure of this paper is outlined as follows. Sec. II presents the field
equation of the scalar perturbation in the EMS theory, reviews the violation
condition of the SCC, and discusses the non-minimal couplings between the
electromagnetic and scalar fields. In Sec. III, we describe the numerical
methods we used to compute the quasinormal modes (QNM) frequency, and we
present the corresponding results. Finally, we summarize and discuss our
findings in Sec. IV, which includes the identification of a critical coupling
constant above which the SCC is reinstated for all black holes in RNdS
spacetime.

| $l=0$ | $l=1$ | $l=2$ | $l=10$ | $l=20$
---|---|---|---|---|---
Pseudospectral method | $0.07967913$ | $0.06393574$ | $0.06182984$ | $0.06075811$ | $0.06071149$
Direct integration method | $0.07967913$ | $0.06393574$ | $0.06182984$ | $0.06075811$ | $0.06071149$
WKB approximation |  |  |  | $0.06073984$ | $0.06070675$

Table 1: The lowest-lying QNMs $\beta=-{\text{Im}(\omega)}/{\kappa_{-}}$ with
different $l$ calculated by different numerical methods for $\Lambda
M^{2}=0.06,\alpha=0.1$ and $Q/Q_{\text{ext}}=0.9$.

| $l=0$ | $l=1$ | $l=2$ | $l=10$ | $l=20$
---|---|---|---|---|---
Pseudospectral method | $0.57893700$ | $0.45043291$ | $0.43826722$ | $0.43203320$ | $0.43124841$
Direct integration method | $0.57893700$ | $0.45043291$ | $0.43826723$ | $0.43203451$ | $0.43176229$
WKB approximation |  |  |  | $0.43209053$ | $0.43177729$

Table 2: The lowest-lying QNMs $\beta=-{\text{Im}(\omega)}/{\kappa_{-}}$ with
different $l$ calculated by different numerical methods for $\Lambda
M^{2}=0.06,\alpha=0.1$ and $Q/Q_{\text{ext}}=0.99$.

##  II Strong cosmic censorship and quasinormal mode in Einstein-Maxwell-
Scalar theory

In this paper, we investigate the Einstein-Maxwell-Scalar (EMS) theory, which
involves a non-minimal coupling between a massless scalar field $\Phi$ and a
Maxwell invariant term. The full action of this theory, given by

$\displaystyle S=\frac{1}{16\pi}\int$ $\displaystyle
d^{4}x\sqrt{-g}\left[R-2\Lambda\right.$ (1)
$\displaystyle\left.-2\nabla_{a}\Phi\nabla^{a}\Phi-f(\Phi)F_{ab}F^{ab}\right],$

includes the cosmological constant $\Lambda$, the Ricci scalar $R$, and the
electromagnetic tensor $F_{ab}$. The coupling function $f(\Phi)$ determines
the strength of the non-minimal coupling of $\Phi$ to the Maxwell term. To
study the scalar field perturbation under first-order approximation, we
require $f^{\prime}(0)=0$, which can be achieved using a quadratic coupling
function $1+\alpha\Phi^{2}$, an exponential coupling function
$e^{\alpha\Phi^{2}}$, or other forms that satisfy

$f(0)=1\,\quad\text{and}\quad f^{\prime\prime}(0)=2\alpha\,,$ (2)

where $\alpha$ is the coupling constant. The equations of motion, given by

$\displaystyle G_{ab}+\Lambda g_{ab}=T_{ab}^{\text{sc}}+T_{ab}^{\text{EM}}\,,$
(3) $\displaystyle\nabla^{2}\Phi-\frac{f^{\prime}(\Phi)}{4}F_{ab}F^{ab}=0\,,$
$\displaystyle\nabla^{a}[f(\Phi)F_{ab}]=0\,,$

are obtained by varying the action and consist of the Einstein equation, the
scalar field equation, and the electromagnetic field equation. The energy-
momentum tensor of the scalar field and the electromagnetic field are denoted
by $T_{ab}^{\text{sc}}$ and $T_{ab}^{\text{EM}}$, respectively, and are
defined by

$\displaystyle T_{ab}^{\text{sc}}$
$\displaystyle=2\nabla_{a}\Phi\nabla_{b}\Phi-
g_{ab}\nabla_{c}\Phi\nabla^{c}\Phi\,,$ (4) $\displaystyle T_{ab}^{\text{EM}}$
$\displaystyle=2f(\Phi)F_{a}{}^{c}F_{bc}-\frac{f(\Phi)}{2}F_{cd}F^{cd}g_{ab}\,.$

The EMS theory allows for the formation of scalarized black hole solutions
through the non-minimal coupling between the scalar field and Maxwell
invariant. By choosing an appropriate coupling function in the EMS theory, the
electromagnetic vacuum solution, namely the RNdS solution, can become unstable
under perturbations, resulting in the formation of scalarized black hole
solutions. However, recent literatures Cai:2020wrp ; An:2021plu  have shown
that nonlinear electrodynamics black holes with scalar hairs, which lack a
Cauchy horizon, effectively satisfy SCC. Therefore, the SCC problem only needs
to be tested in the perturbatively stable black hole solutions, namely the
RNdS solutions, which is given by

$\displaystyle ds^{2}$
$\displaystyle=-f(r)dt^{2}+\frac{1}{f(r)}dr^{2}+r^{2}d\Omega^{2}\,,$ (5)
$\displaystyle A_{a}$ $\displaystyle=-\frac{Q}{r}(dt)_{a}\,,$

with the blackening factor

$\displaystyle f(r)=1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}-\frac{\Lambda
r^{2}}{3}\,,$ (6)

describes a special solution of the EMS theory with vanishing scalar field.

Assuming $r_{c}$, $r_{+}$ and $r_{-}$ are the cosmological horizon, event
horizon and Cauchy horizon respectively, we can rewrite the blackening factor
as:

$f(r)=\frac{\Lambda}{3r^{2}}(r_{c}-r)(r-r_{+})(r-r_{-})(r-r_{o}),$ (7)

where $r_{o}$ is the minimum root of $f(r)=0$ and can be found to be
$r_{o}=-r_{c}-r_{+}-r_{-}$. The surface gravity of each horizon can be defined
as:

$\kappa_{i}=\frac{1}{2}|f^{\prime}(r_{i})|\quad\text{with}\quad i={c,+,-,o}.$
(8)

We can expand the coupled scalar field $\Phi(t,r,\theta,\phi)$ on the RNdS
background spacetime as a perturbation Kokkotas:1999bd ; Berti:2009kk ;
Konoplya:2011qq . Using the symmetries of the spacetime, we can express the
scalar field in terms of spherical harmonics as

$\Phi(t,r,\theta,\phi)=\sum_{lm}e^{-i\omega
t}Y_{lm}(\theta,\phi)\frac{\psi(r)}{r}\,,$ (9)

where $Y_{lm}(\theta,\phi)$ is the spherical harmonics. By substituting Eqs.
(5) and (9) into the equation of motion (LABEL:eomall) satisfied by the scalar
field, we obtain a one-dimensional Schr$\ddot{\text{o}}$dinger-like equation

$\frac{d^{2}\psi(r)}{dr^{2}_{\ast}}+[\omega^{2}-V(r)]\psi(r)=0\,,$ (10)

where the effective potential is given by

$V(r)=\frac{f(r)}{r^{2}}\left[l(l+1)+rf^{\prime}(r)-\frac{\alpha
Q^{2}}{r^{2}}\right]\,,$ (11)

with the tortoise coordinate

$dr_{\ast}=\frac{dr}{f(r)}\,.$ (12)

In the physical region between the event horizon $r_{+}$ and cosmological
horizon $r_{c}$, the tortoise coordinate $r_{\ast}$ can be expressed as a sum
of logarithmic terms involving the surface gravity of each horizon $i$, i.e.,

$\displaystyle r_{\ast}=$
$\displaystyle-\frac{1}{2\kappa_{c}}\ln\left(1-\frac{r}{r_{c}}\right)+\frac{1}{2\kappa_{+}}\ln\left(\frac{r}{r_{+}}-1\right)$
(13)
$\displaystyle-\frac{1}{2\kappa_{-}}\ln\left(\frac{r}{r_{-}}-1\right)+\frac{1}{2\kappa_{o}}\ln\left(1-\frac{r}{r_{o}}\right)\,.$

Notably, the effective potential (11) vanishes on every horizon, implying that
the asymptotic solution of the Schr$\ddot{\text{o}}$dinger-like equation (10)
near each horizon $i$ takes the form:

$\psi\sim e^{\pm i\omega r_{\ast}},,\quad r\rightarrow r_{i}\,,$ (14)

where $e^{i\omega r_{\ast}}$ represents the outgoing wave and the other one
represents the ingoing wave. Physical considerations Konoplya:2011qq  dictate
that there is only an ingoing wave near the event horizon $r_{+}$ and only an
outgoing wave near the cosmological horizon $r_{c}$, which can be expressed as
the following boundary conditions:

$\displaystyle\psi\sim e^{-i\omega r_{\ast}}\,,\quad r\rightarrow r_{+}\,,$
(15) $\displaystyle\psi\sim e^{i\omega r_{\ast}}\,,\quad\,\,\,\,r\rightarrow
r_{c}\,.$

The QNMs are the solutions of Eq. (10) that satisfy the boundary condition
(LABEL:boundarycon), and their frequencies $\omega$ are discrete. The validity
of the SCC for an asymptotically de-Sitter black hole is determined by whether
its weak solution can be extended beyond the Cauchy horizon. It was shown in
Refs. Dias:2018ynt ; Cardoso:2017soq  that this is possible if and only if the
solution is locally square integrable derivative, meaning that it satisfies

$\displaystyle\beta\equiv-\frac{\text{Im}(\omega)}{\kappa_{-}}\geq 1/2\,.$
(16)

As long as there is a QNM having $\beta$ less than $1/2$, the Cauchy horizon
is unable to be extended, which ensures the validity of SCC. Hence, to examine
whether SCC is valid, we only need to investigate the lowest-lying QNM.

Figure 1: When $\Lambda M^{2}=0.06$, the lowest-lying QNMs with the frequency
$\beta=-{\text{Im}(\omega)}/{\kappa_{-}}$ are depicted for various coupling
constants $\alpha$ as a function of the black hole charge ratio
$Q/Q_{\text{ext}}$, for a given $l$. The SCC critical value of the charge
ratio is indicated by the horizontal dashed line at $\beta=1/2$.

##  III Numerical methods and results

In this section, we will present two novel numerical techniques to accurately
determine the QNM frequencies.

While several numerical approaches have been developed for computing QNM
frequencies with high precision Konoplya:2011qq , in this study, we will
introduce a pseudospectral method Jansen:2017oag ; Miguel:2020uln  for
calculating the QNM frequencies. To ensure the accuracy of our results, we
will also validate them with the direct integration method
Chandrasekhar:1975zza ; Molina:2010fb . Additionally, we will adopt the WKB
approximation Konoplya:2003ii  to calculate the photo-sphere modes, which
correspond to the QNMs in the large-$l$ limit.

We notice that the scalar field $\psi(r)$ oscillates significantly near the
two horizons, namely the event horizon and the cosmological horizon, as given
in equations (LABEL:boundarycon). To apply the pseudospectral method to the
field equation (10), we introduce a new variable $y(x)$ defined by

$\displaystyle\begin{aligned}
\psi(r)=\left(x+1\right)^{\frac{i\omega}{2\kappa_{+}}}\left(x-1\right)^{-\frac{i\omega}{2\kappa_{c}}}y(x)\,.\end{aligned}$
(17)

to transform the field equation into a regular form in the interval $[-1,1]$.
In particular, the variable $x$ is defined as a function of $r$ using the
relationship between $r$ and $x$ given in the equation below,

$r=\frac{r_{c}-r_{+}}{2}x+\frac{r_{c}+r_{+}}{2},.$ (18)

After this transformation, the field equation (10) can be written as

$\displaystyle\begin{aligned}
a_{0}(\omega,x)y(x)+a_{1}(\omega,x)y^{\prime}(x)+a_{2}(\omega,x)y^{\prime\prime}(x)=0\,,\end{aligned}$
(19)

which is a second-order differential equation with variable coefficients.

To apply the pseudospectral method, we first expand the field equation and
variable $y(x)$ by the cardinal function $C_{i}(x)$, which satisfies
$C_{i}(x_{j})=\delta_{ij}$, where $x_{j}$ denotes the $j$-th Gauss-Lobatto
point. Substituting this expansion into equation (19) and multiplying each
term by $C_{k}(x)$, we obtain a set of algebraic equations. By imposing the
boundary conditions given in equation (LABEL:boundarycon), we can obtain a
matrix equation of the form $(M_{0}+\omega M_{1})Y=0$, where $M_{0}$ and
$M_{1}$ are matrices and $Y$ is a vector containing the coefficients of the
expansion. Then, the QNM frequency can be obtained by solving the eigenvalue
problem of the matrix $(-M_{1}^{-1}M_{0})$. This method allows us to calculate
the QNM frequencies for scalar perturbations on the RNdS black hole background
with high precision.

Figure 2: The red shaded region in the $\alpha-Q/Q_{\text{ext}}$ diagram
represents the region where the SCC is violated for $\Lambda M^{2}=0.06$ and
$\Lambda M^{2}=0.14$, respectively. The critical value of
$\alpha=\alpha_{\text{crit}}$, represented by the black dot-dashed vertical
line, indicates the point where all black holes satisfy the SCC for EMS theory
with a coupling constant $\alpha>\alpha_{\text{crit}}$.

Tables 2 and 2 display the QNM frequencies obtained from both the
pseudospectral method and the direct integration method. For the direct
integration method, we used the series expansion of $y(x)$ near $x=\pm 1$ as
the boundary condition and solved equation (19) in the intervals $(-1,0]$ and
$[0,1)$. This was done for a given $\omega$ using _Mathematica_ , and the
acceptable frequency $\omega$ was determined by ensuring the two solutions
were smooth at $x=0$.

We have utilized both the pseudospectral method and direct integration method
to calculate the lowest-lying QNMs $\beta=-\text{Im}(\omega)/\kappa_{-}$ for
different $l$ and various black hole parameters. The results are presented in
Tables 2 and 2, which demonstrate the reliability of our numerical
computations. Here, $Q_{\text{ext}}$ represents the electric charge of the
extremal black hole solution, which can be derived from the condition
$r_{+}=r_{-}$. Furthermore, we have employed the WKB approximation to evaluate
the large-$l$ lowest-lying modes, and our results are consistent with other
methods. Therefore, the lowest-lying modes for large $l$ are determined by the
WKB approximation in our study. It has been previously established that the
SCC is hardly violated for the RNdS black holes that are not near extremal.
Our results for the EMS theory are consistent with previous research. Hence,
we only present the QNMs results in the nearly extremal region in this paper.

Fig. 1 displays the lowest-lying QNMs with frequency
$\beta=-{\text{Im}(\omega)}/{\kappa_{-}}$ for different coupling constants
$\alpha$ as a function of the black hole charge ratio $Q/Q_{\text{ext}}$, at a
given $l$. The graphs reveal that, for most coupling constants, there is a
range of charge parameters near the black hole’s extremal limit where all QNMs
violate the SCC by having $\beta>1/2$. We can determine the violation region
of the SCC by identifying the intersection point of the $\beta=1/2$ curve and
the curve of all lowest-lying modes for all possible values of $l$. However,
at specific values of the coupling constant, such as $\alpha=0.5$, all black
holes satisfy the SCC, meaning that all lowest-lying QNMs for all $l$ are less
than $1/2$. Moreover, the lowest-lying QNMs with $l=0$ (nearly extremal modes)
contribute the most to the preservation of the SCC. The plots reveal that the
lowest-lying modes correspond to either $l=0$ or large $l$. The top three
plots represent cases where the coupling constant is negative: $\alpha=-0.1$,
$\alpha=-0.2$, and $\alpha=-0.3$. As the coupling constant becomes more
negative, the lowest-lying QNMs with $l=0$ increase, indicating that restoring
the SCC in these cases requires increasing the coupling constant in the
positive direction. This is expected because a larger coupling constant
$\alpha$ leads to a smaller effective potential, which can even become
negative and destabilize the spacetime under scalar perturbations.

Finally, we investigate the impact of the non-minimal coupling constant
$\alpha$ on the validity of the SCC by plotting the violation region of the
SCC in the $\alpha-Q/Q_{\text{ext}}$ diagram (see Figs. 2) for two different
values of $\Lambda M^{2}$: 0.06 and 0.14. We exclude the region where $\alpha$
is negative since it extends to minus infinity and can be inferred from the
positive $\alpha$ region. Our analysis demonstrates that for small positive
coupling constants $\alpha$, the SCC is violated as the black hole approaches
the extremal limit. However, as $\alpha$ increases, the SCC can be partially
or fully restored, depending on the value of $\alpha$. We identify a critical
value $\alpha_{\text{crit}}$ of the coupling constant, above which the SCC is
reinstated for all black holes in RNdS spacetime, regardless of their charge-
to-mass ratio. In other words, the validity of the SCC no longer depends on
the charge-to-mass ratio of the black hole in these situations. Thus, we can
say that this theory satisfies the SCC. Our results provide valuable insights
into the restoration of the SCC in RNdS spacetime and suggest that the non-
minimal couplings between the electromagnetic and scalar fields may have an
important role in the process.

##  IV conclusion and discussion

In this paper, we have investigated the validity of the SCC in the context of
RNdS black holes within the framework of the EMS theory. We have examined the
impact of scalar field perturbations on the SCC in RNdS spacetime with non-
minimal couplings between the electromagnetic and scalar fields.

Our numerical analysis of the QNM frequencies of the non-minimally coupled
scalar field reveals that the SCC can be reinstated within a specific range of
the coupling constant $\alpha$. By plotting the violation region of the SCC in
the $\alpha-Q/Q_{\text{ext}}$ diagram for two different values of $\Lambda
M^{2}$, we have identified a critical coupling constant
$\alpha_{\text{crit}}$, which enables the EMS theory with
$\alpha>\alpha_{\text{crit}}$ to satisfy the SCC. Our results provide new
insights into the restoration of the SCC in RNdS spacetime and suggest that
non-minimal couplings between the electromagnetic and scalar fields may play
an important role in the restoration of the SCC.

Our findings have important implications for other modified gravitational
theories, and it would be valuable to explore whether similar restoration of
the SCC can be achieved in those theories. Additionally, our analysis has
focused on the linear level of the EMS theory, and it would be worthwhile to
extend our investigation to the non-linear regime. Finally, it would be
interesting to explore the physical implications of the restoration of the SCC
in the context of astrophysical black holes and the evolution of the universe.

## Acknowledgement

Jie Jiang is supported by the National Natural Science Foundation of China
with Grant No. 12205014, the Guangdong Basic and Applied Research Foundation
with Grant No. 217200003 and the Talents Introduction Foundation of Beijing
Normal University with Grant No. 310432102. Jia Tan is supported by the
starting funding of Suzhou University of Science and Technology with Grant no.
332114702, Jiangsu Key Disciplines of the Fourteenth Five-Year Plan with Grant
No. 2021135, and Natural Science Foundation of JiangSu Province (BK20220633).

## References

  * (1) R. Penrose, Riv. Nuovo Cim. 1, 252-276 (1969) doi:10.1023/A:1016578408204
  * (2) D. Christodoulou, “The Formation of Black Holes in General Relativity,” doi:10.1142/9789814374552_0002 [arXiv:0805.3880 [gr-qc]].
  * (3) M. Dafermos, Commun. Pure Appl. Math. 58, 0445-0504 (2005) [arXiv:gr-qc/0307013 [gr-qc]].
  * (4) M. Simpson and R. Penrose, Int. J. Theor. Phys. 7, 183-197 (1973)
  * (5) E. Poisson and W. Israel, Phys. Rev. D 41, 1796-1809 (1990) doi:10.1103/PhysRevD.41.1796
  * (6) V. Cardoso, J. L. Costa, K. Destounis, P. Hintz and A. Jansen, Phys. Rev. Lett. 120, no.3, 031103 (2018) doi:10.1103/PhysRevLett.120.031103 [arXiv:1711.10502 [gr-qc]].
  * (7) P. R. Brady, C. M. Chambers, W. Krivan and P. Laguna, Phys. Rev. D 55, 7538-7545 (1997) doi:10.1103/PhysRevD.55.7538 [arXiv:gr-qc/9611056 [gr-qc]].
  * (8) J. L. Costa, P. M. Girão, J. Natário and J. D. Silva, doi:10.1007/s40818-017-0028-6 [arXiv:1406.7261 [gr-qc]].
  * (9) P. Hintz and A. Vasy, J. Math. Phys. 58, no.8, 081509 (2017) doi:10.1063/1.4996575 [arXiv:1512.08004 [math.AP]].
  * (10) J. L. Costa and A. T. Franzen, Annales Henri Poincare 18, no.10, 3371-3398 (2017) doi:10.1007/s00023-017-0592-z [arXiv:1607.01018 [gr-qc]].
  * (11) Y. Mo, Y. Tian, B. Wang, H. Zhang and Z. Zhong, Phys. Rev. D 98, no.12, 124025 (2018) doi:10.1103/PhysRevD.98.124025 [arXiv:1808.03635 [gr-qc]].
  * (12) B. Ge, J. Jiang, B. Wang, H. Zhang and Z. Zhong, JHEP 01, 123 (2019) doi:10.1007/JHEP01(2019)123 [arXiv:1810.12128 [gr-qc]].
  * (13) X. Liu, S. Van Vooren, H. Zhang and Z. Zhong, JHEP 10, 186 (2019) doi:10.1007/JHEP10(2019)186 [arXiv:1909.07904 [hep-th]].
  * (14) O. J. C. Dias, H. S. Reall and J. E. Santos, Class. Quant. Grav. 36, no.4, 045005 (2019) doi:10.1088/1361-6382/aafcf2 [arXiv:1808.04832 [gr-qc]].
  * (15) V. Cardoso, J. L. Costa, K. Destounis, P. Hintz and A. Jansen, Phys. Rev. D 98, no.10, 104007 (2018) doi:10.1103/PhysRevD.98.104007 [arXiv:1808.03631 [gr-qc]].
  * (16) K. Destounis, Phys. Lett. B 795, 211-219 (2019) doi:10.1016/j.physletb.2019.06.015 [arXiv:1811.10629 [gr-qc]].
  * (17) H. Liu, Z. Tang, K. Destounis, B. Wang, E. Papantonopoulos and H. Zhang, JHEP 03, 187 (2019) doi:10.1007/JHEP03(2019)187 [arXiv:1902.01865 [gr-qc]].
  * (18) K. Destounis, R. D. B. Fontana and F. C. Mena, Phys. Rev. D 102, no.10, 104037 (2020) doi:10.1103/PhysRevD.102.104037 [arXiv:2006.01152 [gr-qc]].
  * (19) S. Hollands, R. M. Wald and J. Zahn, Class. Quant. Grav. 37, no.11, 115009 (2020)
  * (20) E. A. Kontou and K. Sanders, Class. Quant. Grav. 37, no.19, 193001 (2020) doi:10.1088/1361-6382/ab8fcf
  * (21) K. Destounis, R. D. B. Fontana, F. C. Mena and E. Papantonopoulos, JHEP 10, 280 (2019) doi:10.1007/JHEP10(2019)280 [arXiv:1908.09842 [gr-qc]].
  * (22) H. Guo, H. Liu, X. M. Kuang and B. Wang, Eur. Phys. J. C 79, no.11, 891 (2019) doi:10.1140/epjc/s10052-019-7416-x [arXiv:1905.09461 [gr-qc]].
  * (23) A. Sang and J. Jiang, Phys. Rev. D 105, no.8, 084047 (2022) doi:10.1103/PhysRevD.105.084047 [arXiv:2201.00664 [gr-qc]].
  * (24) G. W. Gibbons and K. i. Maeda, Nucl. Phys. B 298, 741-775 (1988)
  * (25) D. Garfinkle, G. T. Horowitz and A. Strominger, Phys. Rev. D 43, 3140 (1991) [erratum: Phys. Rev. D 45, 3888 (1992)] doi:10.1103/PhysRevD.43.3140
  * (26) T. Damour and G. Esposito-Farese, Phys. Rev. Lett. 70, 2220-2223 (1993) doi:10.1103/PhysRevLett.70.2220
  * (27) C. A. R. Herdeiro, E. Radu, N. Sanchis-Gual and J. A. Font, Phys. Rev. Lett. 121, no.10, 101102 (2018) doi:10.1103/PhysRevLett.121.101102 [arXiv:1806.05190 [gr-qc]].
  * (28) P. G. S. Fernandes, C. A. R. Herdeiro, A. M. Pombo, E. Radu and N. Sanchis-Gual, Class. Quant. Grav. 36, no.13, 134002 (2019) [erratum: Class. Quant. Grav. 37, no.4, 049501 (2020)] doi:10.1088/1361-6382/ab23a1 [arXiv:1902.05079 [gr-qc]].
  * (29) Y. S. Myung and D. C. Zou, Phys. Lett. B 790, 400-407 (2019) doi:10.1016/j.physletb.2019.01.046 [arXiv:1812.03604 [gr-qc]].
  * (30) Y. S. Myung and D. C. Zou, Eur. Phys. J. C 79, no.3, 273 (2019) doi:10.1140/epjc/s10052-019-6792-6 [arXiv:1808.02609 [gr-qc]].
  * (31) Y. S. Myung and D. C. Zou, Eur. Phys. J. C 79, no.8, 641 (2019) doi:10.1140/epjc/s10052-019-7176-7 [arXiv:1904.09864 [gr-qc]].
  * (32) Y. S. An, L. Li and F. G. Yang, Phys. Rev. D 104, no.2, 024040 (2021) doi:10.1103/PhysRevD.104.024040 [arXiv:2106.01069 [gr-qc]].
  * (33) R. G. Cai, L. Li and R. Q. Yang, JHEP 03, 263 (2021) doi:10.1007/JHEP03(2021)263 [arXiv:2009.05520 [gr-qc]].
  * (34) K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999) doi:10.12942/lrr-1999-2 [arXiv:gr-qc/9909058 [gr-qc]].
  * (35) E. Berti, V. Cardoso and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009) doi:10.1088/0264-9381/26/16/163001 [arXiv:0905.2975 [gr-qc]].
  * (36) R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793-836 (2011) doi:10.1103/RevModPhys.83.793 [arXiv:1102.4014 [gr-qc]].
  * (37) O. J. C. Dias, F. C. Eperon, H. S. Reall and J. E. Santos, Phys. Rev. D 97, no.10, 104060 (2018).
  * (38) A. Jansen, Eur. Phys. J. Plus 132, no.12, 546 (2017) doi:10.1140/epjp/i2017-11825-9 [arXiv:1709.09178 [gr-qc]].
  * (39) F. S. Miguel, Phys. Rev. D 103, no.6, 064077 (2021) doi:10.1103/PhysRevD.103.064077 [arXiv:2012.10455 [gr-qc]].
  * (40) S. Chandrasekhar and S. L. Detweiler, Proc. Roy. Soc. Lond. A 344, 441-452 (1975) doi:10.1098/rspa.1975.0112
  * (41) C. Molina, P. Pani, V. Cardoso and L. Gualtieri, Phys. Rev. D 81, 124021 (2010) doi:10.1103/PhysRevD.81.124021 [arXiv:1004.4007 [gr-qc]].
  * (42) R. A. Konoplya, Phys. Rev. D 68, 024018 (2003) doi:10.1103/PhysRevD.68.024018 [arXiv:gr-qc/0303052 [gr-qc]].