4D Spherically Symmetric Time-Dependent Quantum Gravity Amplitudes.

J. A. Rosabal<EMAIL_ADDRESS>
Departamento de Electromagnetismo y Electrónica, Universidad de Murcia, Campus
de Espinardo, 30100 Murcia, Spain.

In these short notes, we compute non-perturbatively the time-dependent quantum
gravity amplitudes for a four-dimensional spherically symmetric space-time
with space-like and time-like boundaries. We solve the 4D classical and
quantum constraints in a novel way. We identify the classical solution of the
constraints as a canonical transformation, where the integration constants are
the new variables. We apply this canonical transformation to the path integral
representation of the amplitudes we are interested in. This procedure allows
us to get the time-dependent amplitudes from the path integral without solving
it explicitly. From these amplitudes, we show that for most of the boundary
conditions time evolution in quantum gravity is non-unitary. There is however
a special case where unitary evolution could be achieved.

###### Contents

  1. 1 Introduction
  2. 2 Canonical formalism for spherically symmetric space-time with space-like and time-like boundaries
  3. 3 Spherically symmetric QG with space-like and time-like boundaries
  4. 4 Conclusions
  5. A The boundary measure
  6. B The bulk measure
  7. C Unitarity

##  1 Introduction

Since the seminal works Thiemann:1992jj ; Kastrup:1993br ; Thiemann:1993qh ;
Kuchar:1994zk  spherically symmetric quantum gravity (QG) has become a
paradigm. They showed that this system is fully quantizable in a non-
perturbative way. References Thiemann:1992jj ; Kastrup:1993br ;
Thiemann:1993qh  deal with the problem of gravity quantization using Asthekar
variables Ashtekar:1986yd , which provide a straightforward way for solving
the constraint. In reference Kuchar:1994zk  The main ingredient is a canonical
transformation which reduces the action to the true degree of freedom. In this
last reference, the GQ constraints are also solved, however, it is not an easy
task.

On the other hand in two dimensions there is a straightforward procedure to
solve the constraints Henneaux:1985nw . This method, as well as the previously
described, applies to open or compact spaces but not to spaces with spatial
boundaries. The inclusion of spatial boundaries opens the windows to
incorporate time in a QG description.

Recently Rosabal:2021fao ; Rosabal:2023kdg  there has been some advance in the
study of QG on manifold with spatial boundaries together with boundaries in
time. These kinds of manifolds are of special interest because they provide a
natural framework to compute time-dependent QG amplitudes. For early works on
this subject see Hayward:1992ix ; Hayward:1993my ; Lau:1995fr ; Kummer:1996si
; Hawking:1996ww ; Hawking:1995fd .

This work aims to extend the results of Rosabal:2023kdg  to four dimensions.
In particular, we are interested in the quantum gravity time-dependent
amplitudes for the spherically symmetric space-time. Our results are not
limited to gravity on these particular bounded manifolds. Our results could
also be applied to open and compact spaces with this symmetry. However, for
these cases, no physical time appears in the description, so no time evolution
is available.

The main motivation to compute these amplitudes is because through them we can
diagnose whether or not the time evolution of a QG system is unitary, without
having to compute the time-dependent entropy of the system Hawking:1976ra .

In Rosabal:2023kdg , a recipe to compute the time-dependent quantum gravity
amplitudes in two dimensions was used. Although in this reference it was not
explicitly given, here we list its steps.

  1. 1.

Find a classical solution to the Hamiltonian and the Momentum constraints in
the form (if possible)
$\text{P}_{\text{A}}=\text{M}\big{[}\text{Q}^{\text{A}}\big{]}$. Where
$\big{(}\text{Q}^{\text{A}},\text{P}_{\text{A}}\big{)}$, represent the
coordinates and the momenta of the theory.

  2. 2.

Find a quantum solution to the constraints
($\text{P}_{\text{A}}\rightarrow-\text{i}\frac{\delta{}}{\delta\text{Q}^{\text{A}}}$)
in the form $\Psi=\text{exp}[\text{i}\Omega]$.

  3. 3.

Promote the classical solution to a change of variable in the path integral,
the new variables will be the integration constants associated with this
solution. By construction, it is a canonical transformation.

  4. 4.

Identify $\Omega$ as the generating function of the transformation.

  5. 5.

Perform this transformation to the path integral.

  6. 6.

Solve for the new theory.

  7. 7.

Get the amplitudes.

Although this procedure applies to open or compact spaces, here we are
interested in space-times with spatial boundaries. To account for them we need
to complement the previous list with one more step before getting the
amplitudes. This is

  * •

Reduce the boundary action to be able to apply the canonical transformation.

This step can be done always, for any gravity theory on a manifold with
spatial boundaries. See for instance the form of the boundary action in
(Rosabal:2021fao, ) for more general settings beyond spherically symmetric
space-times. What we do not know is whether the reduced action is always
adapted to the canonical transformation, if it exists. For sure for the two
Rosabal:2023kdg , and four-dimensional (as we will see shortly) cases the
reduced boundary action is perfected and adapted to be transformed by the
canonical transformation associated with the classical solution to the
constraints.

In these short notes, we will apply this recipe to the calculation of the 4D
time-dependent QG amplitudes in the spherically symmetric space-time with
spatial boundaries. The key ingredient will be the use of the same method
proposed by Marc Henneaux in (Henneaux:1985nw, ), and refined by D. Louis-
Martinez et al. in Louis-Martinez:1993bge  to solve the Hamiltonian and the
Momentum constraints but in our case for the four-dimensional Hamiltonian,
which differs from the two dimensional one. Once this solution is found the
rest of the steps can be easily performed, almost mimicking the two-
dimensional case (Rosabal:2023kdg, ). The second fundamental ingredient is a
canonical transformation that allows to reduce the gravity action in the same
spirit of Kuchar:1994zk . The novelty of this transformation is that it can be
constructed directly from the classical solution of the constraints.

The paper is organized as follows, in section 2 we briefly review the
canonical formalism for spherically symmetric space-time with space-like and
time-like boundaries. We present a new procedure to solve the constraint for
this case. This approach is inspired by the method presented in
Henneaux:1985nw  for the two-dimensional case.

In section 3 we develop the QG theory of the spherically symmetric space-time
with spatial boundaries together with boundaries in time. We solve the quantum
version of the constraints. After, we identify the classical and quantum
solutions of the constraints with the canonical transformation and the
generating function of the transformation, respectively. Finally, we present
the time-dependent amplitude. Interestingly enough for most of the boundary
conditions, the amplitude develops a non-unitary behavior. Nonetheless, there
is a very special case in which the amplitude could evolve unitary.
Conclusions are presented in section 4 followed by three appendices. Two of
these appendices are devoted to working out the changes of variables in the
path integral. The third concerns the conditions for the wave function to
evolve unitary.

##  2 Canonical formalism for spherically symmetric space-time with space-like
and time-like boundaries

The canonical formalism for spherically symmetric (open) space-time was
developed in Kuchar:1994zk . The introduction of spatial boundaries together
with the proper boundary terms was presented in Lau:1995fr . The formalism for
more general space-times, beyond those with spherical symmetry, has been
developed in Hayward:1993my ; Hayward:1992ix ; Hawking:1996ww ;
Rosabal:2021fao ; Hawking:1995fd . Here we briefly review it using a more
modern notation. We also present a new procedure to find the classical
solution of the constraints in four dimensions. Although in two dimensions a
similar procedure has been used and extensively celebrated (Henneaux:1985nw, )
in four dimensions, up to the knowledge of the author, it has never been
explored.

The ADM metric for the spherically symmetric space-time reads like

$\text{ds}^{2}=-\text{N}^{2}\text{dt}^{2}+\Lambda^{2}(\text{dr}+\text{N}^{\text{r}}\text{dt})^{2}+4\phi\
\text{ds}^{2}_{\text{S}^{2}}.$ (1)

We have used $\phi$ to contrast more easily with the two-dimensional case
Rosabal:2023kdg . The Einstein Hilbert action in the Hamiltonian form on a
spherically symmetric space-time with spatial boundaries takes the form

$\displaystyle\text{S}=\text{S}_{\text{M}}+\text{S}_{\text{B}}=$
$\displaystyle\int\limits_{t_{i}}^{t_{f}}\text{d}t\int\limits_{0}^{\text{r}_{0}}\text{dr}\Big{[}\text{P}_{\phi}\partial_{t}\phi+\text{P}_{\Lambda}\partial_{t}\Lambda-\text{N}\cal{H}-\text{N}^{\text{r}}\cal{H}_{\text{r}}\Big{]}$
(2)
$\displaystyle+\int\limits_{\text{B}}\text{d}t\Big{(}\Lambda\text{N}^{\text{r}}\
\text{P}_{\Lambda}-2\text{N}\Lambda^{-1}\partial_{\text{r}}\phi+2\partial_{t}\phi\
\text{arcsinh}\big{(}\frac{\Lambda\text{N}^{\text{r}}}{\overline{\text{N}}}\big{)}\Big{)}.$

Where

${\cal
H}=-\frac{1}{2}\text{P}_{\phi}\text{P}_{\Lambda}+2\partial_{\text{r}}\big{(}\Lambda^{-1}\partial_{\text{r}}\phi\big{)}-\frac{1}{2}\Lambda\
+\frac{1}{8}\Lambda\phi^{-1}\Big{(}\text{P}_{\Lambda}^{2}-(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)},$
(3) ${\cal
H}_{\text{r}}=\partial_{\text{r}}\phi\text{P}_{\phi}-\Lambda\partial_{\text{r}}\text{P}_{\Lambda}.$
(4)

The fields in the action (2) are subjected to the following boundary
conditions. On the initial and final slices we fix $(\Lambda,\phi)_{i}$, and
$(\Lambda,\phi)_{f}$, respectively. The induced metric on the spatial
boundary, at a constant $\text{r}=\text{r}_{0}$, takes the form

$\text{ds}^{2}_{\big{|}_{\text{r}_{0}}}=-(\text{N}^{2}-(\Lambda\text{N}^{\text{r}})^{2})\text{dt}^{2}+4\phi\
\text{ds}^{2}_{\text{S}^{2}}=-\overline{\text{N}}^{2}\text{dt}^{2}+4\phi\
\text{ds}^{2}_{\text{S}^{2}}.$ (5)

On this boundary the combination

$\text{N}^{2}-(\Lambda\text{N}^{\text{r}})^{2}=\overline{\text{N}}^{2},$ (6)

is fixed, with $\overline{\text{N}}$, a given function depending only on time.
So, over B, we fix
$\big{(}\sqrt{\text{N}^{2}-(\Lambda\text{N}^{\text{r}})^{2}},\phi\big{)}_{\text{r}_{0}}$.

Classically, variations with respect to N, and $\text{N}^{\text{r}}$, in the
bulk, impose the constraint

${\cal H}={\cal H}_{\text{r}}=0.$ (7)

Note that over the spatial boundary N, and $\text{N}^{\text{r}}$, are not
independent (6), so extra care is needed with these variations over the time-
like boundary.

In the same spirit 111Note that in (Henneaux:1985nw, ), and in general for all
two dimensional dilatonic gravities Louis-Martinez:1993bge , the Hamiltonian
corresponds with the three first terms in (3). of (Henneaux:1985nw, ) we can
make the linear combination

${\cal H}\partial_{\text{r}}\phi+\frac{1}{2}{\cal
H}_{\text{r}}\text{P}_{\Lambda}=0,$

to get

$-\partial_{\text{r}}\Big{(}\text{P}_{\Lambda}^{2}-(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}+\frac{1}{2}\phi^{-1}\partial_{\text{r}}\phi\Big{(}\text{P}_{\Lambda}^{2}-(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}-2\partial_{\text{r}}\phi=0.$
(8)

Previous relation is a differential equation for the function

$\text{F}=\text{P}_{\Lambda}^{2}-(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2},$
(9)
$-\partial_{\text{r}}\text{F}+\frac{1}{2}\phi^{-1}\partial_{\text{r}}\phi\text{F}-2\partial_{\text{r}}\phi=0.$
(10)

The particular form of this equation allows us to solve it in terms of $\phi$,
solely

$\text{F}=\text{C}(t)\sqrt{\phi}-4\phi=\text{P}_{\Lambda}^{2}-(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}.$
(11)

To finally get the full solution of the constraints as

$\displaystyle\text{P}_{\Lambda}$ $\displaystyle=$
$\displaystyle\Big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}^{\frac{1}{2}},$
$\displaystyle\text{P}_{\phi}$ $\displaystyle=$
$\displaystyle\frac{\Lambda}{\partial_{r}\phi}\partial_{\text{r}}\text{P}_{\Lambda}=\frac{\Lambda}{\partial_{r}\phi}\frac{\partial_{r}\Big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}}{2\Big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}^{\frac{1}{2}}}.$
(12)

Note that (2) written in terms of F can be generalized through the inclusion
of more inhomogeneous terms depending on $\phi$ and its derivatives.

$-\partial_{\text{r}}\text{F}+\frac{1}{2}\phi^{-1}\partial_{\text{r}}\phi\text{F}-2\partial_{\text{r}}\phi=\text{W}[\phi].$
(13)

This generalization in four dimensions is equivalent to the statement that all
dilatonic gravities in two dimensions are solvable Louis-Martinez:1993bge . In
four dimensions, however, extra care is needed because $\phi$, is not a scalar
but a component of the metric. This means that $\text{W}[\phi]$, should come
from an invariant term in the action. For instance, a cosmological constant
term in four dimensions is allowed. In this case, in the action, we would have
an extra term

$\lambda\sqrt{-g}=\text{N}\Lambda(4\lambda\phi)\Rightarrow\text{W}[\phi]\propto\lambda\phi\partial_{\text{r}}\phi=\frac{1}{2}\lambda\
\partial_{\text{r}}\phi^{2}.$ (14)

Extensions of W, to include $\Lambda$, and their derivatives are possible but
they are out of the scope of this work.

##  3 Spherically symmetric QG with space-like and time-like boundaries

In the previous section, we have performed the step one of the recipe to get
the amplitudes we are interested in. In this section, we perform the rest of
the steps on the list.

To complete the second point we need to solve the functional differential
equations

$\displaystyle-\text{i}\frac{\delta}{\delta\Lambda}\Psi$ $\displaystyle=$
$\displaystyle\Big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}^{\frac{1}{2}}\Psi,$
$\displaystyle-\text{i}\frac{\delta}{\delta\phi}\Psi$ $\displaystyle=$
$\displaystyle\frac{\Lambda}{\partial_{r}\phi}\frac{\partial_{r}\Big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}}{2\Big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\Big{)}^{\frac{1}{2}}}\Psi.$
(15)

These equations are equivalent to the Wheeler-DeWitt equation and the Momentum
constraint. Their integration is similar to the two-dimensional case
(Henneaux:1985nw, ; Louis-Martinez:1993bge, ; Rosabal:2023kdg, )

$\Psi=\text{exp}\Big{[}\text{i}\ \Omega[\Lambda,\phi;\text{C}]\Big{]},$ (16)

with

$\Omega[\Lambda,\phi;\text{C}]=\Omega^{\prime}[\Lambda,\phi;\text{C}]+\text{G}(\text{C}),$
(17) $\displaystyle\Omega^{\prime}[\Lambda,\phi;\text{C}]=$
$\displaystyle\int\limits_{\Sigma_{t}}\text{dr}\
\Lambda\Big{[}\sqrt{\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}}$
$\displaystyle+(2\Lambda^{-1}\partial_{\text{r}}\phi)\text{ln}\Big{(}\frac{(2\Lambda^{-1}\partial_{\text{r}}\phi)+\sqrt{\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}}}{(2\Lambda^{-1}\partial_{\text{r}}\phi)-\sqrt{\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}}}\Big{)}\Big{]},$
(18)

where $\text{G}(\text{C})$, is a real arbitrary function of the constant of
integration $\text{C}(t)$. This arbitrary function can be found by imposing an
initial condition for the amplitude (which is not (16) but (30)).

If we were computing the wave function on a compact or an open manifold (16)
would be the final solution. Note however that there is no room for
introducing time.

Let us now pose the problem on the corresponding space-time with boundaries.
We are interested in computing the transition amplitude between an initial
configuration $(\Lambda,\phi)_{i}$, at some initial time $t_{i}$, and a final
one $(\Lambda,\phi)_{f}$, at some final time $t_{f}$. While over B the
configuration is $(\overline{\text{N}},\phi)_{\text{B}}$.

In the path integral formulation, this transition amplitude is represented by

$\displaystyle\Psi\Big{[}(\Lambda,\phi)_{f},(\Lambda,\phi)_{i}\
;(\overline{\text{N}},\phi)_{\text{B}}\Big{]}$
$\displaystyle=\int\text{D}\big{[}\text{N},\text{N}^{\text{r}},\phi,\Lambda,\text{P}_{\phi},\text{P}_{\Lambda}\big{]}_{\big{|}_{\text{M}}}$
$\displaystyle\ \ \
\int\text{D}\big{[}\text{N},\text{N}^{\text{r}},\Lambda,\text{P}_{\Lambda}\big{]}_{\big{|}_{\text{B}}}$
$\displaystyle\delta\big{[}\text{N}^{2}-(\Lambda\text{N}^{\text{r}})^{2}-\overline{\text{N}}^{2}\big{]}_{\big{|}_{\text{B}}}\text{e}^{\text{i}\text{S}}.$
(19)

Where the measure $\text{D}\big{[}\ldots\big{]}$, might depend on
$(\phi,\Lambda)$, but it does not depend on the momenta
$(\text{P}_{\phi},\text{P}_{\Lambda})$, nor on the Lagrange multipliers
$(\text{N},\text{N}^{\text{r}})$, except for the functional Dirac delta. Note
that after integration in $(\text{P}_{\phi},\text{P}_{\Lambda})$, we can
recover the Einstein-Hilbert action in its Lagrangian form. See Kuchar:1994zk
, for the Lagrangian form of the action.

Performing the change of variables over the spatial boundaries only

$\displaystyle\Lambda\text{N}^{\text{r}}$ $\displaystyle=$ $\displaystyle{\cal
R}\ \text{sinh}(\eta),$ (20) N $\displaystyle=$ $\displaystyle{\cal R}\
\text{cosh}(\eta),$

where ${\cal R}$, and $\eta$, range in the interval $[0,\infty)$. The
transition amplitude can be written as, appendix A

$\Psi\Big{[}(\Lambda,\phi)_{f},(\Lambda,\phi)_{i}\
;(\overline{\text{N}},\phi)_{\text{B}}\Big{]}=\int\text{D}\big{[}\text{N},\text{N}^{\text{r}},\phi,\Lambda,\text{P}_{\phi},\text{P}_{\Lambda}\big{]}_{\big{|}_{\text{M}}}\int\text{D}\big{[}\eta,\Lambda,\text{P}_{\Lambda}\big{]}_{\big{|}_{\text{B}}}\text{e}^{\text{i}\text{S}}.$
(21)

Where now

$S_{\text{B}}=\int\limits_{\text{B}}\text{d}t\Big{(}\overline{\text{N}}\
\text{P}_{\Lambda}\
\text{sinh}(\eta)-2\overline{\text{N}}\Lambda^{-1}\partial_{\text{r}}\phi\
\text{cosh}(\eta)+2\partial_{t}\phi\eta\Big{)}.$ (22)

Proceeding with the last complementary point on the list we can reduce the
boundary action. To do that we should note first that the measure after the
change of variables does not depend on $\eta$, see appendix A. Second, the
equation of motion for $\eta$, does not contain derivatives. This fact allows
us to solve for $\eta$, in terms of the other boundary degree of freedom and
plug back the solution into the action, to get Rosabal:2023kdg

$\text{S}_{\text{B}}=-\int\limits_{\text{B}}\text{d}t\overline{\text{N}}\sqrt{(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}-\text{P}_{\Lambda}^{2}}.$
(23)

We have assumed that

$\phi_{|_{\text{B}}}=\text{const}\Rightarrow\partial_{t}\phi_{|_{\text{B}}}=0.$
(24)

See also Rosabal:2023kdg ; Lau:1995fr ; Kummer:1996si , where a similar
treatment of the boundary action is performed.

After checking points 3 and 4 on the list we are in a condition of performing
the change of variables (2) in the path integral (21). Note that (2) is the
transformation that maps ${\cal H}$, and ${\cal H}_{\text{r}}$, to zero. So,
we expect the action to be reduced considerably. In addition, it is a
canonical transformation, therefore the new action will contain two boundary
terms at the initial and final time.

After performing the canonical transformation we get

S
$\displaystyle=\int\limits_{t_{i}}^{t_{f}}\text{d}t\int\limits_{0}^{\text{r}_{0}}\text{dr}\Big{[}\text{P}_{\text{C}}\partial_{t}\text{C}+\tilde{\text{P}}_{\phi}\partial_{t}\phi-{\cal
K}[\text{C},\phi,\text{P}_{\text{C}},\tilde{\text{P}}_{\phi}]\Big{]}$
$\displaystyle+\int\limits_{t_{i}}^{t_{f}}\text{d}t\frac{\text{d}}{\text{d}t}\Omega[\Lambda,\phi;\text{C}]-\int\limits_{\text{B}}\text{d}t\
\overline{\text{N}}\sqrt{-\text{F}}.$ (25)

where

$\displaystyle\text{P}_{\Lambda}$ $\displaystyle=$
$\displaystyle\frac{\delta}{\delta\Lambda}\Omega[\Lambda,\phi;\text{C}],$ (26)
$\displaystyle\text{P}_{\phi}$ $\displaystyle=$
$\displaystyle\tilde{\text{P}}_{\phi}+\frac{\delta}{\delta\phi}\Omega[\Lambda,\phi;\text{C}],$
$\displaystyle\text{P}_{\text{C}}$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial\text{C}}\Omega[\Lambda,\phi;\text{C}],$

and the new Hamiltonian ${\cal
K}[\text{C},\phi,\text{P}_{\text{C}},\tilde{\text{P}}_{\phi}]$,

${\cal
K}[\text{C},\phi,\text{P}_{\text{C}},\tilde{\text{P}}_{\phi}]=-\frac{1}{2}\text{N}\tilde{\text{P}}_{\phi}\sqrt{\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}}+\text{N}^{\text{r}}\partial_{\text{r}}\phi\tilde{\text{P}}_{\phi}.$
(27)

From the equation defining $\text{P}_{\text{C}}$, in (26) we can immediately
infer that since $\text{P}_{\text{C}}\in\mathbb{R}$, the arbitrary function in
$\Omega$, $\text{G}(\text{C})$, has to be a real function. Now integration in
N, or $\text{N}^{\text{r}}$, imposes $\tilde{\text{P}}_{\phi}=0$, and the
action further reduces to

$\text{S}=\int\limits_{t_{i}}^{t_{f}}\text{d}t\ \Pi_{\text{C}}\
\partial_{t}\text{C}+\Omega[\Lambda,\phi;\text{C}]\Big{|}_{t_{i}}^{t_{f}}-\int\limits_{\text{B}}\text{d}t\
\overline{\text{N}}\sqrt{-\text{F}},$ (28)

where $\Pi_{\text{C}}(t)=\int\limits_{0}^{\text{r}_{0}}\text{dr}\
\text{P}_{\text{C}}(t)$.

At the classical level the boundary term
$\Omega[\Lambda,\phi;\text{C}]\Big{|}_{t_{i}}^{t_{f}}$, in (28) does not play
any role because it does not affect the equations of motions. However, at the
quantum level, as we are assuming we are performing the change of variables
within the path integral we can not discard this term.

The path integral in the remaining variables
$\big{(}\Pi_{\text{C}}(t),\text{C}(t)\big{)}$, and the new action can be
reduced even more. As the measure does not depend on $\Pi_{\text{C}}$, see
appendix B the path integration in $\Pi_{\text{C}}$, implies
$\text{C}(t)=\text{const}=\text{C}\ \ \forall\ \ t$, and the new action
reduces to

$\text{S}=\Omega[\Lambda,\phi;\text{C}]\Big{|}_{t_{i}}^{t_{f}}-\int\limits_{\text{B}}\text{d}t\
\overline{\text{N}}\sqrt{-\text{F}},$ (29)

where, now C, is a constant.

Putting all these results together we finally get the spherically symmetric
time-dependent amplitude 222See also Rosabal:2023kdg .

$\Psi\Big{[}(\phi,\Lambda)_{f},(\phi,\Lambda)_{i}\
;(\phi,\tau)_{\text{B}}\Big{]}=\int\limits_{0}^{\infty}\text{dC}\
\chi(\text{C})\text{exp}\Big{[}\text{i}\Omega^{\prime}[\Lambda,\phi;\text{C}]\Big{|}_{t_{i}}^{t_{f}}+\sqrt{F}_{|_{\text{B}}}\tau\Big{]},$
(30)

where

$\chi(\text{C})=\text{exp}[\text{i}\text{G}(\text{C})],$ (31)

is a complex function that can be determined by imposing initial conditions
for the functional $\Psi\Big{[}(\phi,\Lambda)_{f},(\phi,\Lambda)_{i}\
;(\phi,\tau)_{\text{B}}\Big{]}$, and

$\tau=\int\limits_{t_{i}}^{t_{f}}\text{d}t\ \overline{\text{N}}.$ (32)

is the proper time over the boundary. This means that in order to fully
determine the function $\chi(\text{C})$, we need to prescribe the state at
$\tau=0\Leftrightarrow t_{i}=t_{f}$, namely

$\Psi\Big{[}(\phi,\Lambda)_{f},(\phi,\Lambda)_{i}\
;(\phi,0)_{\text{B}}\Big{]}.$ (33)

As in the 2D case, Rosabal:2023kdg  in this example several boundary
conditions for the metric on the spatial boundary will lead to a non-unitary
evolution. This is determined by the sign of the function F. Note however that
we also have room for boundary conditions leading to unitary evolution, see
appendix C.

The amplitude (30) will potentially evolve unitarily only if
$\text{F}_{\big{|}_{\text{B}}}<0$, for a given value of $\phi$ over B. When
$\text{C}<4\sqrt{\phi}$, $\text{F}<0$, and vice versa for the other values of
C. This means that to have an entirely negative function for all values of C,
we need to push to infinite the zero of F. In other words,

$\text{F}_{\big{|}_{\text{B}}}<0\ \ \forall\ \ \text{C},\ \text{only when}\ \
\phi_{\big{|}_{\text{B}}}\rightarrow\infty.$ (34)

To have a well-defined limit in the action we need to complement it with a
boundary term (counterterm) evaluated on a classical reference background
Hawking:1995fd . At this point, it is worth pointing out that choosing a
classical solution is equivalent to fixing a value of C.

The boundary action on the reference background takes the same form as in
(23). After performing the canonical transformation, or using the classical
equation of motion for the reference background we have

$\text{S}_{\text{B}_{\text{reference}}}=-\sqrt{-\text{F}_{0}}_{\big{|}_{\text{B}}}\tau,$
(35)

where

$\text{F}_{0}=\text{C}_{0}\sqrt{\phi}-4\phi,$ (36)

with $\text{C}_{0}$, fixed.

It is straightforward to check that

$\lim_{\phi\to\infty}\Big{(}\text{S}_{\text{B}}-\text{S}_{\text{B}_{\text{reference}}}\Big{)}=\lim_{\phi\to\infty}-\text{i}\Big{(}\sqrt{\text{F}}-\sqrt{\text{F}}_{0}\Big{)}\tau=\Big{(}\frac{1}{4}\text{C}_{0}-\frac{1}{4}\text{C}\Big{)}\tau.$
(37)

So, $\text{C}_{0}$, enters in the wave function in the phase
$\text{exp}\Big{(}\text{i}\frac{1}{4}\text{C}_{0}\tau\Big{)}$, namely
$\text{C}_{0}$, does not contribute to the probability

$\Psi\Big{[}(\phi,\Lambda)_{f},(\phi,\Lambda)_{i}\
;(\infty,\tau)_{\text{B}}\Big{]}=\text{e}^{(\text{i}\frac{1}{4}\text{C}_{0}\tau)}\int\limits_{0}^{\infty}\text{dC}\
\chi(\text{C})\text{exp}\Big{[}\text{i}\Omega^{\prime}[\Lambda,\phi;\text{C}]\Big{|}_{t_{i}}^{t_{f}}-\text{i}\frac{1}{4}\text{C}\tau\Big{]}.$
(38)

##  4 Conclusions

In this work, we have computed for the first time the time-dependent
transition amplitudes for a spherically symmetric space-time (30) with space-
like and time-like boundaries. From it, one can show that for most of the
boundary conditions, the quantum evolution is non-unitary. There is however a
case where unitary evolution could be achieved. It is when the
$\phi_{|_{\text{B}}}$, component of the metric goes to infinite on the spatial
boundary (38). In addition, there is one more condition to be satisfied. It is
given in equation (58). This is a very strong condition on the measure of the
internal product between states.

We do not discard the possibility that such a measure exists. Nonetheless, it
is worth emphasizing that only when these two conditions are met this quantum
gravity system will evolve unitary. This is a quite revealing fact since we
have been struggling for decades trying to show that QG is unitary. It turns
out that for most physical configurations it is not, and the windows for which
it could be is very small. Certainly, finding such a measure, or showing it
does not exist will end this battle for unitarity. We leave this study for
future work.

Besides the results commented on in the previous paragraphs, over the paper we
have derived others results. Perhaps the most important among them is the
solutions of the classical (2) and quantum constraints (3). We have used a
method proposed in Henneaux:1985nw  for the case of two-dimensional gravity.
Its extension to four-dimensional gravity, as far as the author is concerned,
is a new contribution never before explored in the literature. The power of
this procedure in four dimensions and its beauty get condensed in equation
(10) (or its generalization (13)), and its solution (11).

In connection with previous results is the canonical transformation described
in section 3. It resembles the one presented in (Kuchar:1994zk, ). The novelty
of this transformation is that it can be identified directly with the
classical solution of the constraints (2). Remarkably, the generating
function, (17) and (18), of this transformation can be derived from the
quantum solution of the constraints (3).

The fact that we can derive a canonical transformation from the classical
solution of the constraints opens new avenues to address some QG systems. As
long as we can find a classical solution to the constraints, as commented in
the introduction, the system will be quantizable.

Certainly, an urgent point should be addressed is the inclusion of matter for
the spherically symmetric space-time on this bounded manifold. With its
inclusion, we will be able to incorporate in the Hilbert space states
describing the black hole formation and evaporation. We leave the inclusion of
matter for future works.

Before ending we would like to comment again and emphasize some points about
the (non) unitary issue. Since the work of Hawking Hawking:1976ra  the physics
community became aware that QG might evolve non-unitary. Since then a lot of
effort has been put into showing that QG is unitary.

In this work we have conclusively shown using the time-dependent amplitudes
and not the entropy that for most of the boundary configurations, pure
spherically symmetric QG in four dimensions is indeed non-unitary. The only
case where it could be unitary is when $\phi_{|_{\text{B}}}\rightarrow\infty$.
However, in addition, one more condition has to be satisfied to ensure
unitarity (58). This condition is related to the measure of the internal
product of the Hilbert space.

The message one should take from this work is that for most of the boundary
configuration, QG evolves non-unitarily, and the case where it could be
unitary is so restrictive that seems not to exist (58).

## Acknowledgments

We thank Jose J. Fernandez-Melgarejo for comments and support and Marc
Henneaux for the correspondence.

##  Appendix A The boundary measure

In this appendix, we develop the path integral boundary measure. Although it
is not needed to fully know it, we need to know certain information at some
points to perform some steps. We are assuming that the measure might depend on
the fields degree of freedom $\text{h}=(\Lambda,\phi)$, but it does not depend
on the momenta $(\text{P}_{\phi},\text{P}_{\Lambda})$, nor on the Lagrange
multipliers $(\text{N},\text{N}^{\text{r}})$.

We start with the boundary measure in (3)

$\text{D}\big{[}\text{N},\text{N}^{\text{r}},\Lambda,\text{P}_{\Lambda}\big{]}_{\big{|}_{\text{B}}}\\\
\delta\big{[}\text{N}^{2}-(\Lambda\text{N}^{\text{r}})^{2}-\overline{\text{N}}^{2}\big{]}_{\big{|}_{\text{B}}}.$

Given what is stated in the previous paragraph, this measure can be
interpreted as

$\prod_{t\in[t_{i},t_{f}]}\text{H}[\text{h}]\
\text{d}\text{N}(t)\wedge\text{d}\text{N}^{\text{r}}(t)\wedge\text{d}\Lambda(t)\wedge\text{d}\text{P}_{\Lambda}(t)\delta\Big{(}\text{N}(t)^{2}-(\Lambda(t)\text{N}^{\text{r}}(t))^{2}-\overline{\text{N}}(t)^{2}\Big{)}_{\big{|}_{\text{B}}}.$
(39)

Now, performing the change of variables, the first we should note is that we
can differentiate both expressions (at a fixed time) in (20) to compute the
wedge product in (39). The first step yields to

$\text{d}\Lambda\text{N}^{\text{r}}\wedge\text{d}\text{N}+\Lambda\text{d}\text{N}^{\text{r}}\wedge\text{d}\text{N}={\cal
R}\text{d}\eta\wedge\text{d}{\cal R}.$ (40)

While a second step brings the product to the desirable form

$\Lambda\text{d}\text{N}\wedge\text{d}\text{N}^{\text{r}}\wedge\text{d}\Lambda={\cal
R}\text{d}\eta\wedge\text{d}\Lambda\wedge\text{d}{\cal R}.$ (41)

The Dirac delta can be straightforwardly massaged to get

$\delta\Big{(}{\cal
R}^{2}-\overline{\text{N}}^{2}\Big{)}=\frac{\delta\Big{(}{\cal
R}-\overline{\text{N}}\Big{)}}{2\overline{\text{N}}}+\frac{\delta\Big{(}{\cal
R}+\overline{\text{N}}\Big{)}}{2\overline{\text{N}}}.$ (42)

As $\overline{\text{N}}>0$, and the integration in ${\cal R}$, is on the
interval $[0,\infty)$, only the first term will contribute to the path
integral.

Putting all the results together we can see that (A) can be expressed as

$\prod_{t\in[t_{i},t_{f}]}\text{H}[\text{h}]_{\big{|}_{\text{B}}}\
\Lambda^{-1}(t){\cal R}(t)\
\text{d}\eta(t)\wedge\text{d}\Lambda(t)\wedge\text{d}{\cal R}(t)\
\frac{\delta\Big{(}{\cal
R}(t)-\overline{\text{N}}(t)\Big{)}}{2\overline{\text{N}}(t)}.$ (43)

Integration in the variable ${\cal R}$, yields to

$\prod_{t\in[t_{i},t_{f}]}\frac{1}{2}\text{H}[\text{h}]_{\big{|}_{\text{B}}}\
\Lambda^{-1}(t)\ \text{d}\eta(t)\wedge\text{d}\Lambda(t).$ (44)

Where now the ${\cal R}$, integration can be straightforwardly performed.

What is important to note from (44) is that it does not depend on $\eta$. This
justifies the manipulations to get (23).

##  Appendix B The bulk measure

In this appendix, we perform the canonical transformation (2) to the bulk
measure

$\text{D}\big{[}\text{N},\text{N}^{\text{r}},\phi,\Lambda,\text{P}_{\phi},\text{P}_{\Lambda}\big{]}_{\big{|}_{\text{M}}}.$
(45)

The interpretation of this measure is similar to the previous one but of
course, the transformation is much more involved.

$\prod_{t\in[t_{i},t_{f}]}\prod_{\text{r}\in[0,\text{r}_{0})}\text{H}[\text{h}]\
\text{d}\text{N}\wedge\text{d}\text{N}^{\text{r}}\wedge\text{d}\phi\wedge\text{d}\Lambda\wedge\text{d}\text{P}_{\phi}\wedge\text{d}\text{P}_{\Lambda}.$
(46)

Under the assumptions about the measure stated in the previous section and to
solve the last integration in $\big{(}\Pi_{\text{C}}(t),\text{C}(t)\big{)}$,
in section 3, it is important to show that the measure does not depend on
$\Pi_{\text{C}}(t)$.

First, we can differentiate (at a fixed (t,r)) the expressions in (2). The
relevant parts are

$\displaystyle\text{d}\text{P}_{\Lambda}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\phi}}{2\big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\big{)}^{\frac{1}{2}}}\text{d}\text{C}+\ldots\
,$ (47) $\displaystyle\text{d}\text{P}_{\phi}$ $\displaystyle=$
$\displaystyle\text{d}\tilde{\text{P}}_{\phi}+\ldots\ ,$ (48)

where the ellipsis represents the terms with the differentials that will
cancel in the wedge product (46). Plugging these expression in (46) we arrive
at

$\prod_{t\in[t_{i},t_{f}]}\prod_{\text{r}\in[0,\text{r}_{0})}\text{H}[\text{h}]\
\text{d}\text{N}\wedge\text{d}\text{N}^{\text{r}}\wedge\text{d}\phi\wedge\text{d}\Lambda\wedge\text{d}\tilde{\text{P}}_{\phi}\wedge\frac{\sqrt{\phi}\
\text{d}\text{C}}{2\big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\big{)}^{\frac{1}{2}}}.$
(49)

Now we need to incorporate $\text{P}_{\text{C}}$, in the measure. From the
last equation, in (26) we can get explicitly $\text{P}_{\text{C}}$, in terms
of the other degrees of freedom.

$\text{P}_{\text{C}}=-\int\limits_{0}^{\text{r}_{0}}\text{dr}\frac{\Big{(}\Lambda^{2}\text{F}+4(\partial_{\text{r}}\phi)^{2}\Big{)}^{\frac{1}{2}}}{2(\text{C}-4\sqrt{\phi})}.$
(50)

Its differential is given by

$\text{dP}_{\text{C}}=-\int\limits_{0}^{\text{r}_{0}}\text{dr}\frac{\sqrt{\phi}\
\text{d}\Lambda}{2\big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\big{)}^{\frac{1}{2}}}+\ldots\
,$ (51)

where the ellipsis represent the term containing $\text{d}\phi(t)$.

Now we can form the infinite wedge product
$\prod_{t\in[t_{i},t_{f}]}\text{dP}_{\text{C}}(t)$. Note that this product
will contain a sum of infinite terms 333For handling the right-hand side of
(51) we need to discretize the integral. but there is one term containing the
product over all the points in the bulk of $\text{d}\Lambda(t,r)$. The other
terms contain at least one $\text{d}\phi(t,\text{r})$. That is 444Strictly
speaking, every term also contains an infinite product of the discretized
differentials $\Delta\text{r}$, coming from the discretization of (51). This
infinite product can be omitted without loss of generality.

$\prod_{t\in[t_{i},t_{f}]}\text{dP}_{\text{C}}=-\text{i}\prod_{t\in[t_{i},t_{f}]}\prod_{\text{r}\in[0,\text{r}_{0})}\frac{\sqrt{\phi}\
\text{d}\Lambda}{2\big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\big{)}^{\frac{1}{2}}}+\ldots\
,$ (52)

where the $-\text{i}$ factor comes from $(-1)^{\sum 1}=(-1)^{-\frac{1}{2}}$.

The observation here is that upon wedge-multiplicating (52) with the infinite
wedge product of the differentials $\text{d}\phi$, the ellipsis contribution
cancels out and we are left with

$\prod_{t\in[t_{i},t_{f}]}\prod_{\text{r}\in[0,\text{r}_{0})}\text{dP}_{\text{C}}\wedge\text{d}\phi=-\text{i}\prod_{t\in[t_{i},t_{f}]}\prod_{\text{r}\in[0,\text{r}_{0})}\frac{\sqrt{\phi}\
\text{d}\Lambda\wedge\text{d}\phi}{2\big{(}\text{F}+(2\Lambda^{-1}\partial_{\text{r}}\phi)^{2}\big{)}^{\frac{1}{2}}}+\ldots\
,$ (53)

From here it is clear that the measure (49) can be written in terms of
$\text{d}\text{P}_{\text{C}}$, as

$\prod_{t\in[t_{i},t_{f}]}\prod_{\text{r}\in[0,\text{r}_{0})}\text{H}[\text{h}]\
\text{d}\text{N}\wedge\text{d}\text{N}^{\text{r}}\wedge\text{d}\phi\wedge\text{d}\tilde{\text{P}}_{\phi}\wedge\text{d}\text{C}\wedge\text{d}\text{P}_{\text{C}}.$
(54)

Note that (54) does not depend on $\text{P}_{\text{C}}$. This fact is what
justifies the manipulations to go from equation (28) to (29).

##  Appendix C Unitarity

One of the axioms of quantum mechanics (QM) or quantum field theory (QFT) is:
The evolution of a closed system is unitary. The state
vector$|\psi(t)\rangle$, at time $t$, is derived from the state vector
$|\psi(t_{0})\rangle$, at time $t_{0}$, by applying a unitary operator
$U(t,t_{0})$, called the evolution operator $|\psi(t)\rangle=U(t,t_{0})\
|\psi(t_{0})\rangle$.

A direct consequence of unitary evolution is the probability $\text{P}(t)$,
conservation

$\text{P}(t)=\langle\psi(t)|\psi(t)\rangle=\langle\psi(t_{0})|\psi(t_{0})\rangle=\text{P}(t_{0}).$
(55)

In other words, unitary evolution is equivalent to say the probability does
not depend on time.

For defining probability in a quantum mechanic’s settings there should exist a
notion of internal product between states $\langle\cdot|\cdot\rangle$. In the
coordinate representation of the wave function, it is defined as

$\langle\psi(t)|\psi(t)\rangle=\int\text{d}x^{\text{d}}\text{M(x)}\psi(t,x)\psi^{*}(t,x),$
(56)

where M(x), is a measure.

Over this paper we have been working with the concept of (non) unitary
evolution, however, we do not have a clear notion of the product between
states. In what follows we will establish the conditions under which the wave
function (38) evolves unitary even without explicitly knowing the internal
product between states.

As seen, there is a possibility of achieving unitary evolution when
$\phi_{|_{\text{B}}}\rightarrow\infty$. Let us pose the calculation of the
probability for a generic state of the form (38)

$\displaystyle\text{P}(\tau)=\int\text{D}[\phi_{f},\Lambda_{f}]\Psi\Big{[}(\phi,\Lambda)_{f},(\phi,\Lambda)_{i}\
;(\infty,\tau)_{\text{B}}\Big{]}\Psi^{*}\Big{[}(\phi,\Lambda)_{f},(\phi,\Lambda)_{i}\
;(\infty,\tau)_{\text{B}}\Big{]}=$
$\displaystyle\int\text{D}[\phi_{f},\Lambda_{f}]\int\limits_{0}^{\infty}\text{dC}\
\int\limits_{0}^{\infty}\text{dC}^{\prime}\
\chi(\text{C})\chi^{*}(\text{C}^{\prime})\text{exp}\Big{[}\text{i}\Big{(}\Omega^{\prime}[\Lambda,\phi;\text{C}]-\Omega^{\prime}[\Lambda,\phi;\text{C}^{\prime}]\Big{)}\Big{|}_{t_{i}}^{t_{f}}\Big{]}$
$\displaystyle\text{e}^{-\text{i}\frac{1}{4}\big{(}\text{C}-\text{C}^{\prime}\big{)}\tau}.$
(57)

Although we do not know the measure $\text{D}[\phi_{f},\Lambda_{f}]$, we can
make some conclusions about the time evolution of the probability.

From (C) we can see that the condition for the probability does not depend on
$\tau$, is

$\int\text{D}\phi_{f}\text{D}\Lambda_{f}\ \text{M}[\phi_{f},\Lambda_{f}]\
\text{exp}\Big{[}\text{i}\Big{(}\Omega^{\prime}[\Lambda_{t_{f}},\phi_{t_{f}};\text{C}]-\Omega^{\prime}[\Lambda_{t_{f}},\phi_{t_{f}};\text{C}^{\prime}]\Big{)}\Big{]}=\text{f}(\text{C})\delta\big{(}\text{C}-\text{C}^{\prime}\big{)},$
(58)

where $\text{f}(\text{C})$, could be any function. In this case, the
probability is

$\text{P}(\tau)=\int\limits_{0}^{\infty}\text{dC}\
\chi(\text{C})\chi^{*}(\text{C})\ \text{f}(\text{C})=\text{P}(0).$ (59)

Equation (58) is extremely difficult to satisfy for some
$\text{M}[\phi_{f},\Lambda_{f}]$. Nonetheless, we do not discard the
possibility that there exists a measure such that (58) holds. We leave this
window open and the search for this measure for future works.

On the one hand, if this measure exists there is only one case where the wave
function evolves unitary. It is when $\phi_{|_{\text{B}}}\rightarrow\infty$.
On the other hand, if this measure does not exist we inevitably must conclude
that time evolution in QG is non-unitary.

## References

  * (1) T. Thiemann and H. A. Kastrup, “Canonical quantization of spherically symmetric gravity in Ashtekar’s selfdual representation,” Nucl. Phys. B 399, 211-258 (1993).
  * (2) H. A. Kastrup and T. Thiemann, “Spherically symmetric gravity as a completely integrable system,” Nucl. Phys. B 425, 665-686 (1994).
  * (3) T. Thiemann, “The Reduced phase space of spherically symmetric Einstein-Maxwell theory including a cosmological constant,” Nucl. Phys. B 436, 681-720 (1995).
  * (4) K. V. Kuchar, “Geometrodynamics of Schwarzschild black holes,” Phys. Rev. D 50, 3961-3981 (1994).
  * (5) A. Ashtekar, “New Variables for Classical and Quantum Gravity,” Phys. Rev. Lett. 57, 2244-2247 (1986).
  * (6) M. Henneaux, “QUANTUM GRAVITY IN TWO-DIMENSIONS: EXACT SOLUTION OF THE JACKIW MODEL,” Phys. Rev. Lett. 54, 959-962 (1985).
  * (7) J. A. Rosabal, “JT Gravity on a Finite Lorentzian Strip: Time dependent Quantum Gravity Amplitudes,” arXiv:2302.11863.
  * (8) J. A. Rosabal, “Quantum gravity on a manifold with boundaries: Schrödinger evolution and constraints,” Eur. Phys. J. C 82, no.7, 589 (2022).
  * (9) D. Louis-Martinez, J. Gegenberg and G. Kunstatter, “Exact Dirac quantization of all 2-D dilaton gravity theories,” Phys. Lett. B 321, 193-198 (1994).
  * (10) G. Hayward and K. Wong, “Boundary Schrodinger equation in quantum geometrodynamics,” Phys. Rev. D 46, 620-626 (1992), Addendum Phys. Rev. D 47, (1993).
  * (11) G. Hayward, “Gravitational action for space-times with nonsmooth boundaries,” Phys. Rev. D 47, 3275-3280 (1993).
  * (12) S. R. Lau, “On the canonical reduction of spherically symmetric gravity,” Class. Quant. Grav. 13, 1541-1570 (1996).
  * (13) W. Kummer and S. R. Lau, “Boundary conditions and quasilocal energy in the canonical formulation of all (1+1) models of gravity,” Annals Phys. 258, 37-80 (1997).
  * (14) S. W. Hawking and C. J. Hunter, “The Gravitational Hamiltonian in the presence of nonorthogonal boundaries,” Class. Quant. Grav. 13, 2735-2752 (1996).
  * (15) S. W. Hawking and G. T. Horowitz, “The Gravitational Hamiltonian, action, entropy and surface terms,” Class. Quant. Grav. 13, 1487-1498 (1996).
  * (16) S. W. Hawking, “Breakdown of Predictability in Gravitational Collapse,” Phys. Rev. D 14, 2460-2473 (1976).


$B_{c}^{+}\to B_{1}^{0{\prime}}\mu^{+}\nu_{\mu}$ | ${-0.17}^{+0.02+0.02}_{-0.00-0.01}$ | $-0.27$ | $0.82$ | $0.52$ | $-0.36$
$B_{c}^{+}\to B_{s1}^{0}e^{+}\nu_{e}$ | ${0.10}^{+0.00+0.00}_{-0.06-0.06}$ | $-0.06$ | $1.00$ | $0.36$ | $-0.06$
$B_{c}^{+}\to B_{s1}^{0{\prime}}e^{+}\nu_{e}$ | ${-0.24}^{+0.07+0.10}_{-0.00-0.00}$ | $-0.32$ | $1.00$ | $0.48$ | $-0.29$
$B_{c}^{+}\to B_{s1}^{0}\mu^{+}\nu_{\mu}$ | ${0.14}^{+0.02+0.00}_{-0.04-0.02}$ | $-0.02$ | $0.73$ | $0.40$ | $-0.15$
$B_{c}^{+}\to B_{s1}^{0{\prime}}\mu^{+}\nu_{\mu}$ | ${-0.14}^{+0.06+0.08}_{-0.00-0.00}$ | $-0.15$ | $0.73$ | $0.48$ | $-0.30$
$B_{c}^{-}\to{D}_{1}^{0}e^{+}\nu_{e}$ | ${-0.05}^{+0.00+0.00}_{-0.00-0.34}$ | $0.14$ | $1.00$ | $0.27$ | $0.15$
$B_{c}^{+}\to{D}_{1}^{0{\prime}}e^{+}\nu_{e}$ | ${-0.01}^{+0.00+0.00}_{-0.00-0.19}$ | $-0.72$ | $1.00$ | $0.65$ | $-0.58$
$B_{c}^{+}\to{D}_{1}^{0}\mu^{+}\nu_{\mu}$ | ${-0.05}^{+0.00+0.00}_{-0.00-0.34}$ | $0.14$ | $1.00$ | $0.27$ | $0.15$
$B_{c}^{+}\to{D}_{1}^{0{\prime}}\mu^{+}\nu_{\mu}$ | ${-0.002}^{+0.001+0.000}_{-0.001-0.188}$ | $-0.71$ | $0.99$ | $0.65$ | $-0.58$
$B_{c}^{+}\to{D}_{1}^{0}\tau^{+}\nu_{\tau}$ | ${-0.05}^{+0.01+0.00}_{-0.01-0.25}$ | $0.04$ | $0.72$ | $0.30$ | $0.07$
$B_{c}^{+}\to{D}_{1}^{0{\prime}}\tau^{+}\nu_{\tau}$ | ${0.17}^{+0.01+0.00}_{-0.00-0.16}$ | $-0.20$ | $0.45$ | $0.62$ | $-0.52$
$B_{c}^{+}\to\chi_{c1}e^{+}\nu_{e}$ | ${-0.58}^{+0.01+0.00}_{-0.01-0.00}$ | $0.35$ | $1.00$ | $0.18$ | $0.40$
$B_{c}^{+}\to\chi_{c1}^{{\prime}}e^{+}\nu_{e}$ | ${-0.04}^{+0.01+0.04}_{-0.01-0.00}$ | $-1.25$ | $1.00$ | $0.89$ | $-0.88$
$B_{c}^{+}\to\chi_{c1}\mu^{+}\nu_{\mu}$ | ${-0.58}^{+0.01+0.00}_{-0.01-0.00}$ | $0.36$ | $0.99$ | $0.18$ | $0.40$
$B_{c}^{+}\to\chi_{c1}^{{\prime}}\mu^{+}\nu_{\mu}$ | ${-0.02}^{+0.01+0.04}_{-0.01-0.00}$ | $-1.20$ | $0.95$ | $0.89$ | $-0.88$
$B_{c}^{+}\to\chi_{c1}\tau^{+}\nu_{\tau}$ | ${-0.41}^{+0.01+0.00}_{-0.01-0.00}$ | $0.11$ | $0.53$ | $0.20$ | $0.34$
$B_{c}^{+}\to\chi_{c1}^{{\prime}}\tau^{+}\nu_{\tau}$ | ${0.28}^{+0.01+0.05}_{-0.02-0.01}$ | $-0.13$ | $-0.04$ | $0.79$ | $-0.76$
Table 6: Predictions of the branching ratios for $B_{c}\to A^{({\prime})}\ell\nu_{\ell}$ decays from other theoretical works. Decays |  RQM Ebert:2010zu |  MGI-CLFQM Li:2023wgq |  CLFQM Shi:2016gqt
---|---|---|---
$B_{c}^{+}\to B_{1}^{0}e^{+}\nu_{e}$ | $3.30\times 10^{-6}$ | $1.636\times 10^{-6}$ | $1.52\times 10^{-5}$
$B_{c}^{+}\to B_{1}^{0{\prime}}e^{+}\nu_{e}$ | $7.20\times 10^{-6}$ | $1.743\times 10^{-5}$ | $7.70\times 10^{-5}$
$B_{c}^{+}\to B_{1}^{0}\mu^{+}\nu_{\mu}$ | $3.00\times 10^{-6}$ | $1.244\times 10^{-6}$ | $1.28\times 10^{-5}$
$B_{c}^{+}\to B_{1}^{0{\prime}}\mu^{+}\nu_{\mu}$ | $5.70\times 10^{-6}$ | $1.461\times 10^{-5}$ | $6.28\times 10^{-5}$
$B_{c}^{+}\to B_{s1}^{0}e^{+}\nu_{e}$ | $2.0\times 10^{-5}$ | $3.946\times 10^{-5}$ | $8.31\times 10^{-5}$
$B_{c}^{+}\to B_{s1}^{0{\prime}}e^{+}\nu_{e}$ | $4.5\times 10^{-5}$ | $5.824\times 10^{-5}$ | $5.38\times 10^{-4}$
$B_{c}^{+}\to B_{s1}^{0}\mu^{+}\nu_{\mu}$ | $1.8\times 10^{-5}$ | $2.793\times 10^{-5}$ | $6.33\times 10^{-5}$
$B_{c}^{+}\to B_{s1}^{0{\prime}}\mu^{+}\nu_{\mu}$ | $3.1\times 10^{-5}$ | $4.230\times 10^{-5}$ | $3.98\times 10^{-4}$
$B_{c}^{-}\to{D}_{1}^{0}e^{+}\nu_{e}$ | $1.1\times 10^{-5}$ | $3.510\times 10^{-5}$ | -
$B_{c}^{+}\to{D}_{1}^{0{\prime}}e^{+}\nu_{e}$ | $1.9\times 10^{-5}$ | $9.84\times 10^{-6}$ | -
$B_{c}^{-}\to{D}_{1}^{0}\mu^{+}\nu_{\mu}$ | - | $3.491\times 10^{-5}$ | -
$B_{c}^{+}\to{D}_{1}^{0{\prime}}\mu^{+}\nu_{\mu}$ | - | $9.79\times 10^{-6}$ | -
$B_{c}^{+}\to{D}_{1}^{0}\tau^{+}\nu_{\tau}$ | $3.9\times 10^{-6}$ | $1.026\times 10^{-5}$ | -
$B_{c}^{+}\to{D}_{1}^{0{\prime}}\tau^{+}\nu_{\tau}$ | $1.1\times 10^{-5}$ | $5.19\times 10^{-6}$ | -
Table 7: Predictions of the branching ratios ($A_{FB}$) of $B_{c}\to\chi_{c1}^{({\prime})}\ell\nu_{\ell}$ decays from other theoretical works. Decays | $B_{c}^{+}\to\chi_{c1}e^{+}\nu_{e}$ | $B_{c}^{+}\to\chi_{c1}^{{\prime}}e^{+}\nu_{e}$ | $B_{c}^{+}\to\chi_{c1}\tau^{+}\nu_{\tau}$ | $B_{c}^{+}\to\chi_{c1}^{{\prime}}\tau^{+}\nu_{\tau}$
---|---|---|---|---
RQM Ebert:2010zu | $8.2\times 10^{-4}$ | $9.6\times 10^{-4}$ | $9.2\times 10^{-5}$ | $7.7\times 10^{-5}$
pQCD Rui:2018kqr | $1.53\times 10^{-3}$ | $1.06\times 10^{-3}$ | $2.0\times 10^{-4}$ | $1.3\times 10^{-4}$
BS Wang:2011jt | $(1.1\pm 0.3)\times 10^{-3}$ | $(2.8\pm 0.8)\times 10^{-3}$ | $(9.7\pm 6.5)\times 10^{-5}$ | $(1.9\pm 1.3)\times 10^{-4}$
NRQM Hernandez:2006gt | $6.6\times 10^{-4}$ | $1.70\times 10^{-3}$ | $7.2\times 10^{-5}$ | $1.5\times 10^{-4}$
$A_{FB}$ | $-0.60_{-0.01}$ | $-0.83_{-0.05}\times 10^{-2}$ | $-0.46$ | $0.35$
RQM Ivanov:2005fd | $9.8\times 10^{-4}$ | $3.10\times 10^{-3}$ | $1.2\times 10^{-4}$ | $2.7\times 10^{-4}$
$A_{FB}$ | $0.19$ | $-3.6\times 10^{-2}$ | $0.34$ | $0.31$
RCQM Ivanov:2006ni | $9.2\times 10^{-4}$ | $2.70\times 10^{-3}$ | $8.9\times 10^{-5}$ | $1.7\times 10^{-4}$
QCDSR Azizi:2009ny | $(1.46\pm 0.42)\times 10^{-3}$ | $(1.42\pm 0.40)\times 10^{-3}$ | $(1.47\pm 0.44)\times 10^{-4}$ | $(1.37\pm 0.38)\times 10^{-4}$
CLFQM Wang:2009mi | $(1.4^{+0.0}_{-0.1})\times 10^{-3}$ | $(3.1^{+0.5}_{-0.8})\times 10^{-3}$ | $(1.5^{+0.10}_{-0.20})\times 10^{-4}$ | $(2.2^{+0.2}_{-0.4})\times 10^{-4}$
MGI-CLFQM Li:2023wgq | $2.953\times 10^{-4}$ | $1.155\times 10^{-3}$ | $0.256\times 10^{-4}$ | $0.051\times 10^{-3}$

Figure 5: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to B_{1}\ell\nu_{\ell}$ decays.

Figure 6: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to B_{1}^{\prime}\ell\nu_{\ell}$ decays.

Figure 7: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to B_{s1}\ell\nu_{\ell}$ decays.

Figure 8: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to B_{s1}^{\prime}\ell\nu_{\ell}$ decays.

Figure 9: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to D_{1}\ell\nu_{\ell}$ decays.

Figure 10: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to D_{1}^{\prime}\ell\nu_{\ell}$ decays.

Figure 11: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to\chi_{c1}\ell\nu_{\ell}$ decays.

Figure 12: $q^{2}$ variation of branching ratios and physical observables in
$B_{c}\to\chi_{c1}^{\prime}\ell\nu_{\ell}$ decays.


# Topological Corner States in Bilayer and Trilayer Systems with
Vertically Stacked Topologically Distinct Layers

Natsuko Ishida<EMAIL_ADDRESS>Research Center for Advanced
Science and Technology, The University of Tokyo, 4-6-1 Komaba, Meguro-ku,
Tokyo 153-8505, Japan  Motohiko Ezawa  Department of Applied Physics, The
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Guangtai
Lu  Research Center for Advanced Science and Technology, The University of
Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan  Wenbo Lin  Institute of
Innovative Research, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-
ku, Tokyo 152-8550, Japan  Yasutomo Ota  Research Center for Department of
Applied Physics and Physico-Informatics, Keio University, 3-14-1 Hiyoshi,
Kohoku-ku, Yokohama, Kanagawa 223-8522, Japan  Yasuhiko Arakawa  Institute for
Nano Quantum Information Electronics, The University of Tokyo, 4-6-1 Komaba,
Meguro-ku, Tokyo 153-8505, Japan  Satoshi Iwamoto  Research Center for
Advanced Science and Technology, The University of Tokyo, 4-6-1 Komaba,
Meguro-ku, Tokyo 153-8505, Japan  Institute of Industrial Science, The
University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan

###### Abstract

We investigate bilayer and trilayer systems composed of topologically
distinct, vertically stacked layers based on the Benalcazar-Bernevig-Hughes
model. We have identified a topological phase transition that significantly
alters the number of the topological corner states in these systems.
Additionally, we find that traditional nested Wilson loop analysis
inaccurately classifies certain phases, leading us to evaluate multipole
chiral numbers (MCNs) as a more appropriate topological invariant. This
approach accurately identifies topological phases and the number of MCNs
coincides with the number of corner modes, enhancing our understanding of
topological insulators and opening avenues for further applications.

##  I Introduction

Topological insulators are fascinating emerging materials found in condensed
matter physics[1, 2, 3]. These materials are distinguished by their unique
combination of an insulating bulk and metallic surfaces, i.e., edge states.
The distinctive behavior of the edge states originates from the topological
features in the bulk band structure giving a non-zero topological invariant,
demonstrating the essential concept of bulk-boundary correspondence[4, 5]. The
basic concept of topological materials has been expanded into classical wave
phenomena, including photonics[6, 7, 8, 9, 10], acoustic waves[11, 12, 13,
14], mechanical waves[15, 16, 17, 18], thermal waves[19], as well as electric
circuits[20, 21, 22, 23, 24]. This expansion has led to attractive
applications such as topological lasers[25, 26].

Expanding beyond the conventional topological insulators, the concept of a
higher-order topological insulators (HOTIs [27, 28]) has been proposed and is
gaining significant attention. This class of materials is characterized by
topological states that are at least two dimensions lower than the system
itself. For instance, in a two-dimensional (2D) HOTI, one can observe zero-
dimensional corner states.

Building on this framework, the majority of research has been directed towards
the study of 2D monolayers[29, 30, 31, 32] and their vertically stacked three-
dimensional structures[33, 34, 35, 36, 37, 38]. These investigations have a
broad scope, aiming to reveal the distinct topological properties intrinsic to
these structures[39, 40, 41], as well as to examine their potential
applications in diverse fields, including nanocavity lasers [42, 43, 44],
nonlinear optics[45, 46], cavity quantum electrodynamics (QED) experiments[47,
48] and interferometer[49, 50]. There are some studies on vertically stacking
of identical models in multilayer topological systems [51, 52, 53]. However,
multilayer systems composed of topologically distinct layers have not been
investigated.

In this paper, we investigate multilayer systems composed of topologically
distinct layers. Our work is based on the Benalcazar-Bernevig-Hughes (BBH)
model [27], a simplest model for studying quadrupole topological insulators.
This model consists of $\pi$-flux square lattices and it supports topological
corner modes. The realization of quadrupole insulators has been demonstrated
through various experimental setups, including optical ring resonator arrays
[54], waveguide arrays[55], electric circuits[21], acoustic system[56], and
microwave resonator[57].

Our findings indicate a phase transition induced by the interlayer coupling,
which significantly influences the emergence and annihilation of topological
corner states. Interestingly, we identified a phase where the conventional
method, using the nested Wilson loop, incorrectly classifies the phase as
trivial, despite the presence of corner states. Therefore, we calculated a
recently developed alternative topological invariant known as multipole chiral
numbers (MCNs) [58] to evaluate the topological phases. Due to its capability
to handle $\mathbb{Z}$-class higher-order topological systems, the MCN has
begun to be applied across various models[59, 60, 61, 62]. We found that the
MCN precisely identifies the phases present in our model and accurately
matches the number of corner modes at zero-energy. By adopting the MCN as a
topological invariant, we not only identify changes in the number of corner
modes but also uncover the phase transition induced by interlayer coupling.

##  II Topological corner states

The BBH lattice consists of four sites in a unit cell where the $\pi$-flux is
inserted by introducing the negative coupling in a plaquette, as shown with
the red bond in Fig. 1 (a). There are two topologically distinct phases,
determined by the ratio of the intra-cell coupling to the inter-cell coupling.
When the intra-cell coupling is weaker than inter-cell coupling, there exists
four degenerate states that form a single set of corner modes, localized at
each corner of the structure within a band gap at energy $E=0$. The bulk
Hamiltonian for the system shown in Fig. 1 (a) is given by

$\mathcal{H}_{\mathrm{BBH}}(k_{x},k_{y})=\begin{pmatrix}0&0&-(\kappa_{1}+\kappa_{2}e^{-ik_{x}})&\kappa_{1}+\kappa_{2}e^{-ik_{y}}\\\
0&0&\kappa_{1}+\kappa_{2}e^{ik_{y}}&\kappa_{1}+\kappa_{2}e^{ik_{x}}\\\
-(\kappa_{1}+\kappa_{2}e^{ik_{x}})&\kappa_{1}+\kappa_{2}e^{-ik_{y}}&0&0\\\
\kappa_{1}+\kappa_{2}e^{ik_{y}}&\kappa_{1}+\kappa_{2}e^{-ik_{x}}&0&0\\\
\end{pmatrix}.$ (1)

Throughout this paper, we use $\kappa_{1}$ ($\kappa_{2}$) and $\kappa_{2}$
($\kappa_{1}$) as the intra-cell and inter-cell couplings for the topological
(trivial) system, respectively, with its corresponding Hamiltonian
$\mathcal{H}_{\text{BBH}}^{\text{Topo}}\left(k_{x},k_{y}\right)$
($\mathcal{H}_{\text{BBH}}^{\text{Tri}}\left(k_{x},k_{y}\right)$). Here, we
suppose $\kappa_{1}<\kappa_{2}$, and both are positive.

The BBH system has two occupied degenerate bands that possess two cancelling
dipole moments and a non-vanishing quadrupole moment. By adapting the Wilson
loop calculation, the degeneracy is lifted, resulting in two distinct single-
band subspaces where Wannier bands become non-degenerate. Consequently, each
Wannier band carries its own topological invariant, which can be obtained
through the nested Wilson loop calculation [27]. Following the calculation, a
quadrupole moment $q_{xy}=\frac{1}{2}$ signifies that the system is
topological, while $q_{xy}=0$ signifies that the system is trivial.

In our study, we consider a vertically stacked configuration composed of
topologically distinct layers, taking into account the vertical nearest-
neighbor site coupling denoted by $t_{\mathrm{v}}$. Fig. 1 (b) illustrates the
considered BBH bilayer system and its unit cell, where the top and bottom
layers are topological and trivial, respectively. For trilayer systems, we
consider two cases as shown in Figs. 1 (d) and (e), where the arrangements of
the topological and trivial layers are flipped.

###  II.1 Topological corner states in BBH bilayer lattice

The Hamiltonian describing the bilayer system, shown in Fig. 1 (b), is given
by

$\mathcal{H}_{\text{bilayer}}=\left(\begin{array}[]{cc}\mathcal{H}_{\text{BBH}}^{\text{Topo}}&\mathcal{H}_{\text{v}}\\\
\mathcal{H}_{\text{v}}&\mathcal{H}_{\text{BBH}}^{\text{Tri}}\end{array}\right),$
(2)

Figure 1:  (a) Blue and red bonds represent the positive and negative
couplings, respectively. Corner states, shown in red at each corner, appear
when $\kappa_{1}<\kappa_{2}$. (b) An illustration of the BBH bilayer system
composed of topological and trivial layers. $\kappa_{1}$ ($\kappa_{2}$) and
$\kappa_{2}$ ($\kappa_{1}$) denote the intra- and inter-cell couplings for
topological (trivial) system, respectively.

(c) The labeling lattice sites a, b, c, and d represent the Hamiltonian basis
for each layer. (d,e) Two distinct BBH trilayer systems where topologically
distinct layers are stacked on top of each other.

where

$\mathcal{H}_{\text{v}}=t_{\text{v}}\left(\begin{array}[]{cccc}0&0&1&0\\\
0&0&0&1\\\ 1&0&0&0\\\
0&1&0&0\end{array}\right)=t_{\text{v}}\sigma_{x}\otimes\mathbb{I}_{2}.$ (3)

Here, we use the site index shown in Fig. 1(c) in which chiral symmetry is
preserved. This ordering is essential for the theoretical analysis presented
later in this paper. $\sigma_{x}$ and $\mathbb{I}_{2}$ are Pauli matrix and
two by two identical matrix, respectively. The band structure has four pairs
of doubly degenerate bands, and exhibits both opening and closing of a band
gap around the $\Gamma$ point as $t_{\mathrm{v}}$ changes. (Fig. 2). The
critical value $t_{\mathrm{c}}^{\rm{bilayer}}$ corresponds to the point at
which the band gap closes, defined as
$t_{\mathrm{c}}^{\rm{bilayer}}=\sqrt{2}(\kappa_{1}+\kappa_{2})$. See appendix
A for details.

Figure 2:  Interlayer coupling $t_{\mathrm{v}}$ dependence of band structures
for BBH bilayer system where top and bottom layers are composed of two
topologically distinct layers. Parameters are chosen as (a)
$t_{\mathrm{v}}$=0.5, (b) $t_{\mathrm{v}}=t_{\mathrm{c}}^{\rm{bilayer}}-1$,
(c) $t_{\mathrm{v}}$=$t_{\mathrm{c}}^{\rm{bilayer}}$ , (d)
$t_{\mathrm{v}}=t_{\mathrm{c}}^{\rm{bilayer}}+1$, $\kappa_{1}$=1, and
$\kappa_{2}$=3. $t_{\mathrm{c}}^{\rm{bilayer}}(=4\sqrt{2}$) denotes the
critical value where a band gap closes. Figure 3:  (a) Energy spectrum vs
$t_{\mathrm{v}}$ for BBH bilayer system. Blue dots represent the bulk
spectrum, with or without gapped edges, whereas the red line corresponds to a
single set of corner modes. The gapped edge states appear when corner modes
exist. (b) Field amplitude summed over the four corner modes at $E=0$ for
$t_{\mathrm{v}}$=$t_{\mathrm{c}}^{\rm{bilayer}}-1$. Parameters are chosen as
$\kappa_{1}$=1, $\kappa_{2}=3$ and $N=20$.

We now consider a finite BBH model with square geometry, with each layer
forming a BBH lattice consisting of $2N\times 2N$ sites. Fig. 3(a) shows the
energy spectrum as a function of varying $t_{\mathrm{v}}$. In the weak
interlayer coupling regime ($t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{bilayer}}$),
the system possesses a single set of corner modes, consisting of a total of
four degenerate states, each predominantly localized at the corners of the
layers. These are depicted by the red line in the figure. Fig. 3(b) shows the
field amplitude of the corner state, summed over for all four states. Notably,
these corner states disappear when the interlayer coupling exceeds the
critical value $t_{\mathrm{c}}^{\rm{bilayer}}$. Sudden disappearance of the
corner states, accompanied by a gap closing, suggests a topological phase
transition associated with increasing interlayer coupling strengths.

In order to determine the topological phases, we calculate the quadrupole
moment $q_{xy}$ [27] via the nested Wilson loop. As a result, we find that
$q_{xy}=\frac{1}{2}$ for $t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{bilayer}}$,
identifying the system as topological, while $q_{xy}=0$ for
$t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{bilayer}}$, indicating a trivial phase.
These results indicate that the interlayer coupling induces a topological
phase transition in the BBH bilayer system. The presence and absence of
topological corner states change at the critical value
$t_{\mathrm{c}}^{\rm{bilayer}}$ bilayer as the band gap closes and re-opens.
This critical value marks the transition point where the number of topological
corner states changes. The detailed calculation of the quadrupole moment is
shown in the Appendix B.

###  II.2 Topological corner states in BBH trilayer lattice

We proceed to study a BBH trilayer system depicted in Fig. 1(d). This
arrangement consists of a trivial layer positioned between two topological
layers, where the layers are coupled by an interlayer coupling
$t_{\mathrm{v}}$.

The system Hamiltonian is given by

$\mathcal{H}_{\text{trilayerI}}=\left(\begin{array}[]{ccc}\mathcal{H}_{\text{BBH}}^{\text{Topo}}&\mathcal{H}_{\text{v}}&0\\\
\mathcal{H}_{\text{v}}&\mathcal{H}_{\text{BBH}}^{\text{Tri}}&\mathcal{H}_{\text{v}}\\\
0&\mathcal{H}_{\text{v}}&\mathcal{H}_{\text{BBH}}^{\text{Topo}}\end{array}\right),$
(4)

with the interlayer coupling matrix $\mathcal{H}_{\text{v}}$. Similar to the
bilayer system, the labeling of sites in the Hamiltonian varies between odd
and even layers to maintain the chiral symmetry, as shown in Fig. 1(c). This
labeling is crucial for applying the methods we will discuss in the next
section. The trilayer system possesses six pairs of doubly degenerate bands.
Analogous to the bilayer system, band gap closes at a critical value
$t_{\mathrm{c}}^{\rm{trilayer}}=\kappa_{1}+\kappa_{2}$ (See appendix A). The
band gap re-opens again when $t_{\mathrm{v}}$ exceeds the critical value.

Figure 4: (a) Energy spectrum vs $t_{\mathrm{v}}$ for BBH trilayer system
(topological-trivial-topological structure). Blue dots represent the gapped
edge states and bulk spectrum. The green and red lines correspond to four and
a single set of corner modes, respectively. (b) Field amplitude summed over
the four corner modes at $E=0$ for
$t_{\mathrm{v}}$=$t_{\mathrm{c}}^{\rm{trilayer}}-1$ and
$t_{\mathrm{v}}$=$t_{\mathrm{c}}^{\rm{trilayer}}+1$. ’Corner 1’ persists even
after exceeding the critical value, while ’corner 2’ vanishes at the critical
value $t_{\mathrm{c}}^{\rm{trilayer}}$. Parameters are chosen as
$\kappa_{1}$=1, $\kappa_{2}=3$, and $N=20$.

For a finite BBH trilayer system with square geometry, with each layer
consisting of $2N\times 2N$ sites, we calculate the eigenenergies as a
function of the interlayer coupling strength(Fig 4(a)). When
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$, which is below the critical
value, there are two sets of corner states, with a total number of eight
(illustrated by the green line in the figure). In contrast, when
$t_{\mathrm{v}}$ exceeds $t_{\mathrm{c}}^{\rm{trilayer}}$, only a single set
of corner states survives, with a total number of four (indicated by the red
line). Each of these sets consists of four degenerate modes located at $E=0$.
These two sets of corner states are distinguished by their unique field
distributions(Fig. 4 (b)). Specifically, the corner states that disappear at
the critical interlayer coupling show non-zero field amplitudes across all
layers, which corresponds to symmetric coupling between the corners (corner
2). Conversely, the corner states that remain after exceeding the critical
interlayer coupling exhibit a complete absence of field amplitude in the
middle layer, which corresponds to anti-symmetric coupling between the corners
(corner 1). The presence of topological corner states, ’Corner 1’, at large
interlayer couplings, is understood by an effective model in the vicinity of
the zero-energy, derived based on the perturbation theory. See Appendix C.

Figure 5:  (a) Energy spectrum vs $t_{\mathrm{v}}$ for BBH trilayer system
(trivial-topological-trivial structure). Blue dots correspond to the bulk
spectrum, with or without gapped edges, whereas the red line represents a
single set of corner modes. Gapped edges are present when the system exhibits
corner states. (b) Field amplitude summed over the four corner modes at $E=0$
for $t_{\mathrm{v}}$=$t_{\mathrm{c}}^{\rm{trilayer}}-1$. Parameters are chosen
as $\kappa_{1}$=3, $\kappa_{2}=1$ and $N=20$.

As discussed in the subsection $\mathrm{II.A}$ of the bilayer system, we
calculate the quadrupole moment $q_{xy}$ to identify the topological phases of
the trilayer system. Consequently, our analysis shows that $q_{xy}=0$ for
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$, classifying the system as
trivial, despite the presence of two sets of corner states at $E=0$. In the
case where $t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{trilayer}}$, we have
$q_{xy}=\frac{1}{2}$, indicating the system as topological. This leads to the
issue where the quadrupole moment fails to fully capture the topological
phases in certain instances. Therefore, the introduction of another invariant
is essential to accurately describe the topological phases of both bilayer and
trilayer BBH systems. We will calculate and discuss this alternative invariant
in the next section.

Next, we study a trilayer system which involves a topological layer sandwiched
between two trivial layers, with each layer coupled by an interlayer coupling
$t_{\mathrm{v}}$. This arrangement is depicted in Fig. 1 (e). In this case,
the Hamiltonian is given by

$\mathcal{H}_{\text{trilayerII}}=\left(\begin{array}[]{ccc}\mathcal{H}_{\text{BBH}}^{\text{Tri}}&\mathcal{H}_{\text{v}}&0\\\
\mathcal{H}_{\text{v}}&\mathcal{H}_{\text{BBH}}^{\text{Topo}}&\mathcal{H}_{\text{v}}\\\
0&\mathcal{H}_{\text{v}}&\mathcal{H}_{\text{BBH}}^{\text{Tri}}\end{array}\right).$
(5)

Its critical value is determined by
$t_{\mathrm{c}}^{\rm{trilayer}}=\kappa_{1}+\kappa_{2}$, which is identical to
that in the previous trilayer case.

Fig. 5(a) shows the energy spectrum as a function of $t_{\mathrm{v}}$. There
is a single set of corner modes that includes four degenerate corner states at
$E=0$ (depicted by the red line) when
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$; however, all corner modes
disappear above the critical value $t_{\mathrm{c}}^{\rm{trilayer}}$. The field
amplitudes of the corner modes, shown in Fig. 5 (b), are distributed across
each layer. The quadrupole moments are obtained as $q_{xy}=\frac{1}{2}$ for
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$, identifying the phase as
topological, while $q_{xy}=0$ for
$t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{trilayer}}$, classifying it as trivial. We
note that in BBH trilayer systems with large interlayer coupling
$t_{\mathrm{v}}$, the size of the band-gap varies between two distinct
systems, despite them having identical band structures, as shown in Figs 4(a)
and 5(a). This discrepancy arises due to the presence of gapped edge states
within the band gap as illustrated in Fig. 4(a).

##  III 2D winding number in real space

As we observed a discrepancy between the quadrupole moment and the presence of
corner states when $t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$, in the
trilayer case (’topological-trivial-topological’ structure), we need to
calculate another topological invariant to discuss its topological phase in
detail. We choose the MCN as the additional invariant, which is capable of
accurately evaluating the topological phases of $\mathbb{Z}$-class chiral-
symmetric higher-order topological systems, as recently introduced by [58].
The MCNs are analogous to the 2D winding number in real space, based on the
work described in Ref.[63]. This development is especially crucial in
situations where existing theoretical frameworks, such as the nested Wilson
loop, may erroneously categorize certain phases as trivial. The MCN provides
an innovative approach for analyzing these systems. Notably, the value of MCN
directly correlates with the number of degenerate zero-energy corner states
present in a finite system, thereby offering a more precise invariant for
examining these unique topological states. In this study, we use the MCNs to
identify the topological phases of both BBH bilayer and trilayer systems in
their chiral-symmetric forms. This is particularly relevant in cases where the
conventional nested Wilson approach inaccurately classified the phase as
trivial, despite the presence of zero-energy corner states.

To compute the topological invariant for BBH bilayer system, we modify the
Hamiltonian in its chiral-symmetric form for its finite-size system as
follows:

$\mathcal{H}_{\text{bilayer}}=\begin{pmatrix}0&h_{\text{bilayer}}\\\
h^{\dagger}_{\text{bilayer}}&0\end{pmatrix},$ (6)

Note that the labeling of sites in the Hamiltonian differs between layer 1 and
layer 2 in order to preserve the chiral symmetry, as shown in Fig. 1(c). Also,
sites a and b belong to sublattice A, while sites c and d belong to sublattice
B. When a 2D square-lattice has $N$ unit cells in the $x$ and $y$-directions
(i.e. each layer consists of $N\times N$ unit cells), $h_{\text{bilayer}}$ is
a $4N^{2}\times 4N^{2}$ matrix. Here, the eigenstates of
$\mathcal{H}_{\text{bilayer}}$ in equation (6) are denoted as
$\ket{\psi}^{T}=(a_{1}^{1},\cdots,a_{N^{2}}^{1},b_{1}^{1},\cdots,b_{N^{2}}^{1},a_{1}^{2},\cdots,a_{N^{2}}^{2},b_{1}^{2},\cdots,b_{N^{2}}^{2},c_{1}^{1},\cdots,c_{N^{2}}^{1},d_{1}^{1},\cdots,d_{N^{2}}^{1},c_{1}^{2},\cdots,c_{N^{2}}^{2},d_{1}^{2},\cdots,d_{N^{2}}^{2})$.
Superscripts represent the layer index, where layer 1 and 2 correspond to
upper and lower layer, respectively. The quadrupole moment operator for the
BBH bilayer system is represented by a following sublattice quadrupole moment
operator[58, 64]

$\mathcal{Q}_{xy}^{S}=\sum_{\mathbf{R},\alpha\in\mathcal{S}}e^{-i\frac{2\pi
xy}{N^{2}}}\ket{\mathbf{R},\alpha}\bra{\mathbf{R},\alpha}$ (7)

for $\mathcal{S}=A,B$. Each layer consists of $N\times N$ unit cells, each
labeled as $\mathbf{R}=(x,y)$. We diagonalize the Hamiltonian
$\mathcal{H}_{\text{bilayer}}$ by using the singular value decomposition (SVD)
of $h_{\text{bilayer}}=U_{A}\Sigma U_{B}^{\dagger}$, where $U_{A}$ ($U_{B}$)
denotes a $4N^{2}\times 4N^{2}$ unitary matrix while $\Sigma$ is a diagonal
matrix of singular values. Using
$\mathcal{\tilde{Q}}_{xy}^{S}=U_{\mathcal{S}}^{\dagger}\mathcal{Q}_{xy}^{S}U_{\mathcal{S}}$,
the MCN is defined as [58]

$N_{xy}=\frac{1}{2\pi
i}\mathrm{Tr}\left[\mathrm{log}(\mathcal{\tilde{Q}}_{xy}^{A}\mathcal{\tilde{Q}}_{xy}^{B\dagger})\right].$
(8)

Our calculations reveal that the MCN invariant is $N_{xy}=1$ for
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{bilayer}}$ while $N_{xy}=0$ for
$t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{bilayer}}$, where the MCN number matches
the number of corner states at $E=0$. This result is consistent with the
Wilson loop calculation. The comprehensive calculations for MCNs are presented
in the Appendix D. It is important to note that to calculate the topological
invariant, sufficiently large systems are required.

This calculation can readily be adapted for trilayer cases. First, we rewrite
the trilayer Hamiltonian in equation (4) and (5) into a form that preserves
chiral symmetry, given by

$\mathcal{H}_{\text{trilayer}}=\begin{pmatrix}0&h_{\text{trilayer}}\\\
h^{\dagger}_{\text{trilayer}}&0\end{pmatrix}.$ (9)

The Hamiltonian $h_{\text{trilayer}}$ is represented by a matrix of size
$6N^{2}\times 6N^{2}$. By using the quadrupole moment operators in equation
(7), we obtain the MCN, given by equation (8). For a trilayer system depicted
in Fig. 1(d), the invariant is obtained as $N_{xy}=2$ for
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$ while $N_{xy}=1$ for
$t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{trilayer}}$. The MCN number perfectly
matches the number of corner states at $E=0$. Similarly, for the alternative
trilayer system illustrated in Fig. 1(e), we obtain $N_{xy}=1$ for
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$ and $N_{xy}=0$ for
$t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{trilayer}}$, which again coincides with
the number of corner states at $E=0$. Therefore, we conclude that the MCN
accurately captures the topological phases that even the quadrupole moment may
fail to identify.

##  IV Discussion and Conclusions

We investigated bilayer and trilayer BBH structures, constructed by stacking
topologically distinct layers. By adjusting the interlayer coupling strengths,
we identified changes in the number of corner modes, indicating a topological
phase transition. This topological phase transition, accompanied by band gap
closing, is induced by the interlayer coupling.

Figure 6:  Topological phase diagram of the BBH trilayer system as a function
of the the coupling strengths ratio and the interlayer coupling, highlighting
areas defined by distinct multipole chiral numbers $N_{xy}$ and quadrupole
moments $q_{xy}$. The closing of band gaps is marked by red solid lines.

To identify the topological phases, we first calculated the quadrupole moments
using the nested Wilson loop approach, a well-known method for analyzing
quadrupole insulators. The quadrupole moments explained the topological phases
of BBH bilayer and ’trivial-topological-trivial’ trilayer structures. However,
we observed a mismatch between the quadrupole moments and the presence of
corner modes in ’topological-trivial-topological’ trilayer structure when the
interlayer coupling was below the critical value. In contrast, MCNs, which are
a topological invariant for $\mathbb{Z}$-class chiral-symmetric higher-order
topological systems recently reported, can systematically explain the results.
We also confirmed that the number of corner states coincides with the MCNs
across all interlayer coupling regimes for all three structures we considered.
As we identified a topological phase that may have been misclassified as
trivial when analyzed through the conventional nested Wilson loop approach,
this finding highlights the need for the topological invariant, the MCN, as a
more appropriate method of evaluating $\mathbb{Z}$-class higher-order
topological systems.

The categorizing the phase transition induced by the interlayer coupling is
below:

Bilayer (topological-trivial)

from the topological phase with $N_{xy}=1$ to the trivial phase with
$N_{xy}$=0.

Trilayer (topological-trivial-topological)

from the topological phase with $N_{xy}=2$ to the topological phase with
$N_{xy}=1$.

Trilayer (trivial-topological-trivial)

from the topological phase with $N_{xy}=1$ to the trivial phase with
$N_{xy}$=0.

Figure 6 displays the topological phase diagram, which illustrates both the
MCNs and the quadrupole moments in BBH trilayer systems.

This study demonstrates that adjusting the interlayer coupling manipulates the
topological phase transition, which is accompanied by the appearance and
disappearance of corner modes, potentially expanding the scope of
applications. The configuration of vertically stacking layers, as discussed in
our study, can be implemented in existing platforms, such as electrical
circuits. Additionally, the field distribution of corner modes suggests that
the intensity of corner modes can be accumulated through stacked layers, which
can be adapted to develop a novel topological laser. We believe that our
findings provide valuable insights into the manipulation of localized modes
within multilayer systems, opening avenues for innovative device engineering.

## Appendix A: bulk energy for the BBH bilayer and trilayer system

($\mathrm{i}$) The bulk spectrum of the bilayer system
$\mathcal{H}_{\text{bilayer}}$ is obtained as

$E=\pm\sqrt{2F+t_{\text{v}}^{2}\pm\sqrt{2G}},$ (10)

with

$\displaystyle F$ $\displaystyle=$
$\displaystyle\kappa_{1}^{2}+\kappa_{2}^{2}+\kappa_{1}\kappa_{2}\left(\cos
k_{x}+\cos k_{y}\right),$ (11) $\displaystyle G$ $\displaystyle=$
$\displaystyle\left(\kappa_{1}+\kappa_{2}\right)^{2}t_{\text{v}}^{2}\left(\cos
k_{x}+\cos k_{y}+2\right),$ (12)

where each energy is two-fold degenerated. At the $\Gamma$ point, the energy
is obtained as

$\displaystyle E=\pm\sqrt{2}|\kappa_{1}+\kappa_{2}|\pm t_{\text{v}}.$ (13)

The gap closes at

$\displaystyle t_{\text{v}}=\pm\sqrt{2}(\kappa_{1}+\kappa_{2}).$ (14)

($\mathrm{ii}$) The bulk spectrum of the trilayer systems
$\mathcal{H}_{\text{trilayerI}}$ and $\mathcal{H}_{\text{trilayerII}}$ is
obtained as

$E=\pm\sqrt{2}\sqrt{F},\qquad\pm\sqrt{2}\sqrt{F+t_{\text{v}}^{2}\pm\sqrt{G}},$
(15)

where each energy is two-fold degenerated. It is remarkable that the bulk
spectrum is identical between $\mathcal{H}_{\text{trilayerI}}$ and
$\mathcal{H}_{\text{trilayerII}}$. At the $\Gamma$ point, the energy is
obtained as

$\displaystyle
E=\pm\sqrt{2}|\kappa_{1}+\kappa_{2}|,\qquad\pm\sqrt{2}(\kappa_{1}+\kappa_{2}\pm
t_{\text{v}}).$ (16)

The gap closes at

$\displaystyle t_{\text{v}}=\pm(\kappa_{1}+\kappa_{2}).$ (17)

## Appendix B: Quadrupole moment $q_{xy}$ in BBH system

In the BBH model, the bulk quadrupole moment $q_{xy}$, its boundary
polarizations $p_{x,y}$ and corner charges $Q$ have the following
relationship: $|p^{\mathrm{edge}\pm y}_{x}|=|p^{\mathrm{edge}\pm
x}_{y}|=|Q^{\mathrm{corner}\pm x,\pm y}|=|q_{xy}|$ [27]. The requirements for
quantizing the quadrupole moment $q_{xy}$ are existence of three symmetries,
which are inversion symmetry and $x,y$-reflection symmetry under a gauge
transformation. These constraints lead to quantized values of the Wannier band
polarizations
$p_{y}^{v_{x}^{\pm}},p_{x}^{v_{y}^{\pm}}\overset{\mathcal{I},\mathcal{M}_{x},\mathcal{M}_{y}}{=}0$
or $1/2$. Furthermore, when the system preserves $C_{4}$ symmetry, the
classification of the Wannier bands is $\mathbb{Z}_{2}$ [33]. Therefore, the
quadrupole moment $q_{xy}$ can be described by
$q_{xy}=p_{y}^{v_{x}^{\pm}}=p_{x}^{v_{y}^{\pm}}=0$ or $1/2$ and the
topological invariant is defined by the polarizations which are obtained as
follows:

$p_{y}^{v_{x}^{\pm}}=-\frac{i}{2\pi}\frac{1}{N_{x}}\sum_{k_{x}}\log[\tilde{W}_{y,k_{x}}^{\pm}]=0\
\mathrm{or}\ \frac{1}{2}$ (18)

where $\tilde{W}_{y,k_{x}}^{\pm}$ is a nested Wilson loop along $k_{y}$ [27].
The quadrupole moment of $\frac{1}{2}$ and 0 imply that the system is
topological and trivial, respectively.

Figure 7:  Wannier bands of (a) BBH bilayer, (b) ’Topological-Trivial-
Topological’ BBH trilayer systems and (c) ’Trivial-Topological-Trivial’ BBH
trilayer systems. (a1)-(c1) Below the critical interlayer coupling strength
($t_{\mathrm{v}}$=$t_{\mathrm{c}}$-1). (a2)-(c2) Above the critical interlayer
coupling strength ($t_{\mathrm{v}}$=$t_{\mathrm{c}}$+1). Each band is assigned
a number that represents its polarization. Parameters are chosen as
$\kappa_{1}$=1 and $\kappa_{2}=3$ for (a) and (b), whereas the values are
reversed ($\kappa_{1}$=3 and $\kappa_{2}=1$) for (c).

By adapting the Wilson loop calculation to the occupied bands of the BBH
bilayer shown in Fig. 2, we obtain four non-degenerate Wannier bands, as
depicted in Fig. 7(a1,a2). The gapping of the Wannier bands results from the
non-commutativity of the two reflection symmetries,
$\mathcal{M}_{x}^{\mathrm{bi}}$ and $\mathcal{M}_{y}^{\mathrm{bi}}$, as
discussed in [27]. These two reflection symmetries in the bilayer system are
derived from the direct product of the monolayer symmetries,
$\mathcal{M}_{x,y}^{\mathrm{mono}}$, expressed as
$\mathcal{M}_{x,y}^{\mathrm{bi}}=\mathbb{I}_{2}\otimes\mathcal{M}_{x,y}^{\mathrm{mono}}$.

The Wannier bands exhibit distinct structures below and above the critical
interlayer coupling strength, $t_{\mathrm{c}}^{\rm{bilayer}}$. The nested
Wilson loop calculation of the Wannier bands yields band polarizations of
$p_{y}=\\{1/2,0,0,1/2\\}$ for $t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{bilayer}}$
and $p_{y}=\\{0,0,0,0\\}$ for $t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{bilayer}}$.
This change in polarization, corresponding with the disappearance of corner
states in finite systems, indicates a transition to a strong interlayer
coupling regime, reducing the quantized quadrupole moment from one to zero.

The Wilson loop calculation applied to the occupied bands of the trilayer
structure in Fig. 1(d), yields six non-degenerate Wannier bands (Fig.
7(b1,b2)). As $t_{\mathrm{v}}$ increases, the Wannier bands experience
changes, leading to different polarizations. For
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$, we obtain polarizations of
$p_{y}=\\{1/2,1/2,0,0,1/2,1/2\\}$, identifying the phase as trivial even with
the existence of corner states. In contrast, for
$t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{trilayer}}$, the polarizations are
$p_{y}=\\{1/2,0,0,0,0,1/2\\}$, indicating a topological phase.

For the alternative trilayer structure depicted in Fig. 1(e), the Wilson loop
calculation generates six non-degenerate Wannier bands (Fig. 7(c1,c2)). The
nested Wilson loop computation yields polarizations of
$p_{y}=\\{1/2,0,0,0,0,1/2\\}$ for
$t_{\mathrm{v}}<t_{\mathrm{c}}^{\rm{trilayer}}$ and $p_{y}=\\{0,0,0,0,0,0\\}$
for $t_{\mathrm{v}}>t_{\mathrm{c}}^{\rm{trilayer}}$. These polarization
changes correspond with the number of corners obtained in finite-sized system
calculations.

## Appendix C: Perturbation theory

A numerical simulation shows that there are four corner states at large
interlayer couplings $t_{\text{v}}$ in the ’topological-trivial-topological’
trilayer structure. In order to understand it, we derive the effective
Hamiltonian for large interlayer couplings $t_{\text{v}}$. We first
diagonalize the interlayer coupling Hamiltonian $H_{0}$ defined by

$\mathcal{H}_{0}=\left(\begin{array}[]{ccc}0&\mathcal{H}_{\text{v}}&0\\\
\mathcal{H}_{\text{v}}&0&\mathcal{H}_{\text{v}}\\\
0&\mathcal{H}_{\text{v}}&0\end{array}\right),$ (19)

where only the interlayer coupling exist. It is diagonalized as

$U_{0}^{-1}\mathcal{H}_{0}U_{0}=\left(\begin{array}[]{ccc}-\sqrt{2}\mathbb{I}_{4}&0&0\\\
0&0&0\\\ 0&0&\sqrt{2}\mathbb{I}_{4}\end{array}\right),$ (20)

where $\mathbb{I}_{4}$ is the four by four identical matrix and the
transformation matrix $U_{0}$ is given by

$U_{0}=\left(\begin{array}[]{ccc}\mathbb{I}_{4}&-\sqrt{2}\mathbb{I}_{2}\otimes\sigma_{x}&\mathbb{I}_{4}\\\
-\mathbb{I}_{4}&0&\mathbb{I}_{4}\\\
\mathbb{I}_{4}&\sqrt{2}\mathbb{I}_{2}\otimes\sigma_{x}&\mathbb{I}_{4}\end{array}\right).$
(21)

We are interested in the middle four by four matrix, where there are three
degenerate zero energy states in (20). Next, we transform the full Hamiltonian
and take the middle four by four matrix, which is given by

$P_{0}^{-1}U_{0}^{-1}\mathcal{H}_{\text{trilayerI}}U_{0}P_{0}=\left(\begin{array}[]{cccc}0&0&\kappa_{1}+\kappa_{2}e^{-ik_{x}}&\kappa_{1}+\kappa_{2}e^{ik_{y}}\\\
0&0&\kappa_{1}+\kappa_{2}e^{-ik_{y}}&-\left(\kappa_{1}+\kappa_{2}e^{ik_{x}}\right)\\\
\kappa_{1}+\kappa_{2}e^{ik_{x}}&\kappa_{1}+\kappa_{2}e^{ik_{y}}&0&0\\\
\kappa_{1}+\kappa_{2}e^{-ik_{y}}&-\left(\kappa_{1}+\kappa_{2}e^{-ik_{x}}\right)&0&0\end{array}\right),$
(22)

where $P_{0}$ is a projection operator taking the middle four by four matrix.
It is the effective Hamiltonian, which is also the BBH model. Hence, there are
four corner modes for large interlayer couplings $t_{\text{v}}$.

## Appendix D: Chiral symmetric Hamiltonian for BBH bilayer system

As an illustrative example, we provide a detailed calculation for $N=2$ for
bilayer systems. The submatrix of BBH bilayer Hamiltonian in equation (6) is
composed of four $2N^{2}\times 2N^{2}$ submatrices, described as

$h_{\text{bilayer}}=\begin{pmatrix}h_{\text{Topo}}&\mathcal{H}_{\text{vp}}\\\
\mathcal{H}_{\text{vp}}&h_{\text{Tri}}\end{pmatrix},$ (23)

where

$h_{\text{Topo}}=\begin{pmatrix}-\kappa_{1}&0&0&0&\kappa_{1}&0&0&0\\\
-\kappa_{2}&-\kappa_{1}&0&0&0&\kappa_{1}&0&0\\\
0&0&-\kappa_{1}&0&\kappa_{2}&0&\kappa_{1}&0\\\
0&0&-\kappa_{2}&-\kappa_{1}&0&\kappa_{2}&0&\kappa_{1}\\\
\kappa_{1}&0&\kappa_{2}&0&\kappa_{1}&\kappa_{2}&0&0\\\
0&\kappa_{1}&0&\kappa_{2}&0&\kappa_{1}&0&0\\\
0&0&\kappa_{1}&0&0&0&\kappa_{1}&\kappa_{2}\\\
0&0&0&\kappa_{1}&0&0&0&\kappa_{1}\end{pmatrix},$ (24)

and $\mathcal{H}_{\text{vp}}=t_{\text{v}}\times\mathbb{I}_{2N^{2}}$. Note that
$h_{\text{Topo}}$ and $h_{\text{Tri}}$ have an inverse relationship with
respect to the ratio between $\kappa_{1}$ and $\kappa_{2}$.

For $N=2$, the sublattice quadrupole moment operator $\mathcal{Q}_{xy}^{S}$ of
the BBH bilayer system is a diagonal matrix where the diagonal elements form a
repeating sequence of ($-i,-1,-1,1$), this sequence repeats $2N$ times. Thus,
$\mathcal{Q}_{xy}^{S}$ = $\mathbb{I}_{4}\otimes\text{diag.}(-i,-1,-1,1)$.

The calculation process for bilayer systems can be readily extended to
trilayer systems. The submatrix of BBH trilayer Hamiltonian is

$h_{\text{trilayer}}=\begin{pmatrix}h_{\text{Topo}}&\mathcal{H}_{\text{vp}}&0\\\
\mathcal{H}_{\text{vp}}&h_{\text{Tri}}&\mathcal{H}_{\text{vp}}\\\
0&\mathcal{H}_{\text{vp}}&h_{\text{Topo}}\end{pmatrix},$ (25)

The matrix $\mathcal{Q}_{xy}^{S}$ of the BBH trilayer system is given by
$\mathcal{Q}_{xy}^{S}$ = $\mathbb{I}_{6}\otimes\text{diag.}(-i,-1,-1,1)$ for
$N=2$.

## Acknowledgement

N.I. was supported by the JSPS KAKENHI (Grant No. JP21J40088). M.E. was
supported by CREST, JST (Grants No. JPMJCR20T2) and by the Grants-in-Aid for
Scientific Research from MEXT KAKENHI (Grant No.23H00171). Y.O. was supported
by the Grants-in-Aid for Scientific Research from MEXT KAKENHI (Grants
No.22H01994 and No. 22H00298). S.I. was supported by CREST, JST (Grant No.
JPMJCR19T1) and the Grants-in-Aid for Scientific Research from MEXT KAKENHI
(Grants No 22H00298 and No. 22H01994).

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Bansil  _et al._ [2016] A. Bansil, H. Lin, and T. Das, Colloquium: Topological band theory, Rev. Mod. Phys. 88, 021004 (2016).
  * Hatsugai [1993a] Y. Hatsugai, Chern number and edge states in the integer quantum hall effect, Phys. Rev. Lett. 71, 3697 (1993a).
  * Hatsugai [1993b] Y. Hatsugai, Edge states in the integer quantum hall effect and the riemann surface of the bloch function, Phys. Rev. B 48, 11851 (1993b).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * Khanikaev  _et al._ [2013] A. B. Khanikaev, S. Hossein Mousavi, W.-K. Tse, M. Kargarian, A. H. MacDonald, and G. Shvets, Photonic topological insulators, Nature Materials 12, 233 (2013).
  * Lu  _et al._ [2014] L. Lu, J. D. Joannopoulos, and M. Soljačić, Topological photonics, Nature Photonics 8, 821 (2014).
  * Ota  _et al._ [2020] Y. Ota, K. Takata, T. Ozawa, A. Amo, Z. Jia, B. Kante, M. Notomi, Y. Arakawa, and S. Iwamoto, Active topological photonics, Nanophotonics 9, 547 (2020).
  * Iwamoto  _et al._ [2021] S. Iwamoto, Y. Ota, and Y. Arakawa, Recent progress in topological waveguides and nanocavities in a semiconductor photonic crystal platform [invited], Opt. Mater. Express 11, 319 (2021).
  * Yang  _et al._ [2015] Z. Yang, F. Gao, X. Shi, X. Lin, Z. Gao, Y. Chong, and B. Zhang, Topological acoustics, Phys. Rev. Lett. 114, 114301 (2015).
  * He  _et al._ [2016] C. He, X. Ni, H. Ge, X.-C. Sun, Y.-B. Chen, M.-H. Lu, X.-P. Liu, and Y.-F. Chen, Acoustic topological insulator and robust one-way sound transport, Nature Physics 12, 1124 (2016).
  * Peng  _et al._ [2016] Y.-G. Peng, C.-Z. Qin, D.-G. Zhao, Y.-X. Shen, X.-Y. Xu, M. Bao, H. Jia, and X.-F. Zhu, Experimental demonstration of anomalous floquet topological insulator for sound, Nature Communications 7, 13368 (2016).
  * Xue  _et al._ [2022] H. Xue, Y. Yang, and B. Zhang, Topological acoustics, Nature Reviews Materials 7, 974 (2022).
  * Kane and Lubensky [2014] C. L. Kane and T. C. Lubensky, Topological boundary modes in isostatic lattices, Nature Physics 10, 39 (2014).
  * Süsstrunk and Huber [2015] R. Süsstrunk and S. D. Huber, Observation of phononic helical edge states in a mechanical topological insulator, Science 349, 47 (2015).
  * Huber [2016] S. D. Huber, Topological mechanics, Nature Physics 12, 621 (2016).
  * Chen  _et al._ [2021] C.-W. Chen, R. Chaunsali, J. Christensen, G. Theocharis, and J. Yang, Corner states in a second-order mechanical topological insulator, Communications Materials 2, 62 (2021).
  * Hu  _et al._ [2022] H. Hu, S. Han, Y. Yang, D. Liu, H. Xue, G.-G. Liu, Z. Cheng, Q. J. Wang, S. Zhang, B. Zhang, and Y. Luo, Observation of topological edge states in thermal diffusion, Advanced Materials 34, 2202257 (2022).
  * Ezawa [2018a] M. Ezawa, Higher-order topological electric circuits and topological corner resonance on the breathing kagome and pyrochlore lattices, Phys. Rev. B 98, 201402(R) (2018a).
  * Imhof  _et al._ [2018] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Topolectrical-circuit realization of topological corner modes, Nature Physics 14, 925 (2018).
  * Lee  _et al._ [2018] C. H. Lee, S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, and R. Thomale, Topolectrical circuits, Communications Physics 1, 39 (2018).
  * Li  _et al._ [2018] Y. Li, Y. Sun, W. Zhu, Z. Guo, J. Jiang, T. Kariyado, H. Chen, and X. Hu, Topological lc-circuits based on microstrips and observation of electromagnetic modes with orbital angular momentum, Nature Communications 9, 4598 (2018).
  * Ezawa [2019] M. Ezawa, Non-hermitian higher-order topological states in nonreciprocal and reciprocal systems with their electric-circuit realization, Phys. Rev. B 99, 201411(R) (2019).
  * Harari  _et al._ [2018] G. Harari, M. A. Bandres, Y. Lumer, M. C. Rechtsman, Y. D. Chong, M. Khajavikhan, D. N. Christodoulides, and M. Segev, Topological insulator laser: Theory, Science 359, eaar4003 (2018).
  * Bandres  _et al._ [2018] M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides, and M. Khajavikhan, Topological insulator laser: Experiments, Science 359, eaar4005 (2018).
  * Benalcazar  _et al._ [2017a] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Quantized electric multipole insulators, Science 357, 61 (2017a).
  * Benalcazar  _et al._ [2017b] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators, Phys. Rev. B 96, 245115 (2017b).
  * Ezawa [2018b] M. Ezawa, Higher-order topological insulators and semimetals on the breathing kagome and pyrochlore lattices, Phys. Rev. Lett. 120, 026801 (2018b).
  * Ezawa [2018c] M. Ezawa, Minimal models for wannier-type higher-order topological insulators and phosphorene, Phys. Rev. B 98, 045125 (2018c).
  * Xie  _et al._ [2018] B.-Y. Xie, H.-F. Wang, H.-X. Wang, X.-Y. Zhu, J.-H. Jiang, M.-H. Lu, and Y.-F. Chen, Second-order photonic topological insulator with corner states, Phys. Rev. B 98, 205147 (2018).
  * Wang  _et al._ [2021] H.-X. Wang, L. Liang, B. Jiang, J. Hu, X. Lu, and J.-H. Jiang, Higher-order topological phases in tunable c3 symmetric photonic crystals, Photon. Res. 9, 1854 (2021).
  * Lin and Hughes [2018] M. Lin and T. L. Hughes, Topological quadrupolar semimetals, Phys. Rev. B 98, 241103(R) (2018).
  * Schindler  _et al._ [2018] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Higher-order topological insulators, Science Advances 4, eaat0346 (2018).
  * Kooi  _et al._ [2018] S. H. Kooi, G. van Miert, and C. Ortix, Inversion-symmetry protected chiral hinge states in stacks of doped quantum hall layers, Phys. Rev. B 98, 245102 (2018).
  * Ghorashi  _et al._ [2020] S. A. A. Ghorashi, T. Li, and T. L. Hughes, Higher-order weyl semimetals, Phys. Rev. Lett. 125, 266804 (2020).
  * Wei  _et al._ [2021] Q. Wei, X. Zhang, W. Deng, J. Lu, X. Huang, M. Yan, G. Chen, Z. Liu, and S. Jia, 3d hinge transport in acoustic higher-order topological insulators, Phys. Rev. Lett. 127, 255501 (2021).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Ota  _et al._ [2019] Y. Ota, F. Liu, R. Katsumi, K. Watanabe, K. Wakabayashi, Y. Arakawa, and S. Iwamoto, Photonic crystal nanocavity based on a topological corner state, Optica 6, 786 (2019).
  * Chen  _et al._ [2019] X.-D. Chen, W.-M. Deng, F.-L. Shi, F.-L. Zhao, M. Chen, and J.-W. Dong, Direct observation of corner states in second-order topological photonic crystal slabs, Phys. Rev. Lett. 122, 233902 (2019).
  * Li  _et al._ [2020] M. Li, D. Zhirihin, M. Gorlach, X. Ni, D. Filonov, A. Slobozhanyuk, A. Alù, and A. B. Khanikaev, Higher-order topological states in photonic kagome crystals with long-range interactions, Nature Photonics 14, 89 (2020).
  * Smirnova  _et al._ [2020] D. Smirnova, A. Tripathi, S. Kruk, M.-S. Hwang, H.-R. Kim, H.-G. Park, and Y. Kivshar, Room-temperature lasing from nanophotonic topological cavities, Light: Science & Applications 9, 127 (2020).
  * Zhang  _et al._ [2020] W. Zhang, X. Xie, H. Hao, J. Dang, S. Xiao, S. Shi, H. Ni, Z. Niu, C. Wang, K. Jin, X. Zhang, and X. Xu, Low-threshold topological nanolasers based on the second-order corner state, Light: Science & Applications 9, 109 (2020).
  * Kim  _et al._ [2020] H.-R. Kim, M.-S. Hwang, D. Smirnova, K.-Y. Jeong, Y. Kivshar, and H.-G. Park, Multipolar lasing modes from topological corner states, Nature Communications 11, 5758 (2020).
  * Kirsch  _et al._ [2021] M. S. Kirsch, Y. Zhang, M. Kremer, L. J. Maczewsky, S. K. Ivanov, Y. V. Kartashov, L. Torner, D. Bauer, A. Szameit, and M. Heinrich, Nonlinear second-order photonic topological insulators, Nature Physics 17, 995 (2021).
  * Ezawa [2022] M. Ezawa, Nonlinear non-hermitian higher-order topological laser, Phys. Rev. Res. 4, 013195 (2022).
  * Xie  _et al._ [2020] X. Xie, W. Zhang, X. He, S. Wu, J. Dang, K. Peng, F. Song, L. Yang, H. Ni, Z. Niu, C. Wang, K. Jin, X. Zhang, and X. Xu, Cavity quantum electrodynamics with second-order topological corner state, Laser & Photonics Reviews 14, 1900425 (2020).
  * Kim  _et al._ [2021] E. Kim, X. Zhang, V. S. Ferreira, J. Banker, J. K. Iverson, A. Sipahigil, M. Bello, A. González-Tudela, M. Mirhosseini, and O. Painter, Quantum electrodynamics in a topological waveguide, Phys. Rev. X 11, 011015 (2021).
  * Li  _et al._ [2021] C.-A. Li, S.-B. Zhang, J. Li, and B. Trauzettel, Higher-order Fabry-Pérot interferometer from topological hinge states, Phys. Rev. Lett. 127, 026803 (2021).
  * Chaou  _et al._ [2023] A. Y. Chaou, P. W. Brouwer, and N. Sedlmayr, Hinge states of second-order topological insulators as a mach-zehnder interferometer, Phys. Rev. B 107, 035430 (2023).
  * Cook and Nielsen [2023] A. M. Cook and A. E. B. Nielsen, Finite-size topology, Phys. Rev. B 108, 045144 (2023).
  * Flores-Calderon  _et al._ [2023] R. Flores-Calderon, R. Moessner, and A. M. Cook, Time-reversal invariant finite-size topology, Phys. Rev. B 108, 125410 (2023).
  * Michał J. Pacholski [2023] A. M. C. Michał J. Pacholski, Crystalline finite-size topology, arXiv:2312.08552 (2023).
  * Mittal  _et al._ [2019] S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, Photonic quadrupole topological phases, Nature Photonics 13, 692 (2019).
  * Schulz  _et al._ [2022] J. Schulz, J. Noh, W. A. Benalcazar, G. Bahl, and G. von Freymann, Photonic quadrupole topological insulator using orbital-induced synthetic flux, Nature Communications 13, 6597 (2022).
  * Serra-Garcia  _et al._ [2018] M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Observation of a phononic quadrupole topological insulator, Nature 555, 342 (2018).
  * Peterson  _et al._ [2018] C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, A quantized microwave quadrupole insulator with topologically protected corner states, Nature 555, 346 (2018).
  * Benalcazar and Cerjan [2022] W. A. Benalcazar and A. Cerjan, Chiral-symmetric higher-order topological phases of matter, Phys. Rev. Lett. 128, 127601 (2022).
  * Wang  _et al._ [2023] D. Wang, Y. Deng, J. Ji, M. Oudich, W. A. Benalcazar, G. Ma, and Y. Jing, Realization of a $\mathbb{Z}$-classified chiral-symmetric higher-order topological insulator in a coupling-inverted acoustic crystal, Phys. Rev. Lett. 131, 157201 (2023).
  * Li  _et al._ [2023] Y. Li, J.-H. Zhang, F. Mei, B. Xie, M.-H. Lu, J. Ma, L. Xiao, and S. Jia, Large-chiral-number corner modes in $\mathbb{Z}$-class higher-order topolectrical circuits, Phys. Rev. Appl. 20, 064042 (2023).
  * Qi  _et al._ [2024] Y. Qi, Z. He, K. Deng, J. Li, and Y. Wang, Multipole higher-order topological semimetals, Phys. Rev. B 109, L060101 (2024).
  * Yang  _et al._ [2024] W.-J. Yang, S.-F. Li, X.-Y. Zou, and J.-C. Cheng, Characterization of zero-energy corner states in higher-order topological systems with chiral symmetry, Phys. Rev. B 109, 024114 (2024).
  * Lin  _et al._ [2021] L. Lin, Y. Ke, and C. Lee, Real-space representation of the winding number for a one-dimensional chiral-symmetric topological insulator, Phys. Rev. B 103, 224208 (2021).
  * Wheeler  _et al._ [2019] W. A. Wheeler, L. K. Wagner, and T. L. Hughes, Many-body electric multipole operators in extended systems, Phys. Rev. B 100, 245135 (2019).