# Covariance-Based Joint Device Activity and Delay Detection in Asynchronous
mMTC

Zhaorui Wang, Ya-Feng Liu, and Liang Liu The work of Y.-F. Liu was supported
in part by the National Natural Science Foundation of China (NSFC) under Grant
12022116 and Grant 12021001. The work of L. Liu was supported by the Research
Grants Council, Hong Kong, China, under Grant 25215020. _(Corresponding
author: Ya-Feng Liu.)_ Z. Wang is with the Department of Information
Engineering, The Chinese University of Hong Kong, Hong Kong (e-mail:
zrwang2009@gmail.com).Y.-F. Liu is with the State Key Laboratory of Scientific
and Engineering Computing, Institute of Computational Mathematics and
Scientific/Engineering Computing, Academy of Mathematics and Systems Science,
Chinese Academy of Sciences, Beijing 100190, China (e-mail:
yafliu@lsec.cc.ac.cn).L. Liu is with the Department of Electronic and
Information Engineering, The Hong Kong Polytechnic University, Hong Kong
(e-mail: liang-eie.liu@polyu.edu.hk).

###### Abstract

In this letter, we study the joint device activity and delay detection problem
in asynchronous massive machine-type communications (mMTC), where all active
devices asynchronously transmit their preassigned preamble sequences to the
base station (BS) for device identification and delay detection. We first
formulate this joint detection problem as a maximum likelihood estimation
problem, which depends on the received signal only through its sample
covariance, and then propose efficient coordinate descent type of algorithms
to solve the formulated problem. Our proposed covariance-based approach is
sharply different from the existing compressed sensing (CS) approach for the
same problem. Numerical results show that our proposed covariance-based
approach significantly outperforms the CS approach in terms of the detection
performance since our proposed approach can make better use of the BS antennas
than the CS approach.

###### Index Terms:

Asynchronous mMTC, coordinate descent, joint activity and delay detection,
random access

##  I Introduction

This letter studies massive machine-type communications (mMTC) which provides
efficient random access communications for a large number of devices, out of
which only a small number of them are active[1, 2]. To reduce communication
latency, a grant-free random access scheme was proposed in [3, 4], where the
active devices can directly transmit the data signals after the preamble
signals without permissions from the base station (BS). Note that in practice,
since the IoT devices are equipped with low-cost local oscillators, they are
not perfectly synchronized. This letter studies the problem of joint device
activity and delay detection in the asynchronous mMTC system.

Most of the existing works on massive random access are based on the perfect
synchronization assumption, see [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18] and the references therein. There are generally two lines of
researches for massive random access. One line of research is to apply the
compressed sensing (CS) approach to solve the joint activity detection and
channel estimation problem by taking advantage of the sporadical nature of the
device activity pattern [4, 5, 6, 7, 8, 9, 10, 11]. This line of approaches
generally suffer from poor channel estimation performance (especially in
massive MIMO) as shown in [5, 19] since the channels are estimated through
nonorthogonal preamble sequences. The second line of research is to apply the
covariance-based approach[12, 13, 14, 15, 16, 17, 18] to solve the problem.
Specifically, if only device activity is detected, but channels are not
estimated (e.g., a sensor to send a fire alarm whose activity itself is an
alarming signal), the activity can be detected merely based on the covariance
of the received preamble signals. The minimum preamble sequence length
required in the covariance-based approach is much smaller than that of the CS
approach for the same performance[12, 14]. However, in many applications, data
needs to be transmitted after the preamble, and it is important to estimate
the channels based on the preamble. In this case, [20] proposed a new
protocol. Specifically, at the first phase, the BS applies the covariance-
based approach to detect device activity; at the second phase, the BS
transmits a common feedback message to all the active devices to select
orthogonal preambles; at the third phase, the BS estimates channels and
detects data from the active devices based on the orthogonal preambles. In
short, since the channels are estimated through orthogonal preambles in the
three-phase protocol[20], channel estimation performance should be much better
than that from the grant-free protocol[3, 4] using the CS approach.

All the above works are based on the perfect synchronization assumption. This
is impractical since it is hard to perfectly synchronize all devices in mMTC
in general. In this case, the unknown delays of the active devices should be
estimated as well. However, the current detection framework based on the
synchronous assumption is inapplicable since the sensing/signature matrix in
all of these works is not consistent with the actual one with delays.
Recently, [21] considered the more practical _asynchronous_ massive random
access scheme, and applied the CS approach to solve the joint device activity,
delay detection, and channel estimation problem.

In this letter, we consider the asynchronous massive random access scheme, and
follow the covariance-based three-phase communication protocol in [20]. We
particularly focus on the joint activity and delay detection problem in the
first phase of the three-phase protocol in [20], and propose to use the
covariance-based approach to solve it. We emphasize that the delay detection
is necessary under the protocol proposed in [20] since the knowledge of delay
can be utilized in the rest of phases. In this case, a key challenge is how to
efficiently estimate device activity and delay at the BS based on the
preassigned preambles. To solve the above challenge, we formulate the
detection problem as a maximum likelihood estimation problem. Based on the
sample covariance of the received signals, we first adopt the coordinate
descent (CD) algorithm[12], then propose a block coordinate descent (BCD)
algorithm to solve the detection problem, respectively. Numerical results show
that our proposed covariance-based approaches have significant performance
improvement compared with the CS approach in [21]. A key intuition is that,
the covariance-based approach can make better use of the BS antennas such that
its solution is more accurate than that of the CS approach. It is worth
mentioning that the techniques developed in this paper can be further extended
to the cell-free massive connectivity[22] in which the devices are more
difficult to synchronize than the single-cell setup studied in this paper.

##  II System Model and Problem Formulation

###  II-A System Model

This letter studies a narrowband uplink system which consists of a BS equipped
with $M$ antennas and $N$ single-antenna devices. Due to the sporadic traffic,
only $K\ll N$ devices are active at each coherence block. Let $\alpha_{n}$ be
an activity indicator of device $n$, and

$\alpha_{n}=\left\\{\begin{array}[]{ll}1,~{}{\rm if~{}device}~{}n~{}{\rm
is~{}active},\\\ 0,~{}{\rm otherwise},\end{array}~{}~{}n=1,\dots,N.\right.$
(1)

The set containing the indices of $K$ active devices is denoted as
$\mathcal{K}=\\{n:\alpha_{n}=1,\forall n\\}$. We assume block-fading channels,
i.e., the channels remain roughly constant within each coherence block, but
may vary among different coherence blocks. Let
$\sqrt{\beta_{n}}\mbox{\boldmath{$h$}}_{n}\in\mathbb{C}^{M\times 1}$ be the
channel between device $n$ and the BS, where $\beta_{n}$ is the large-scale
fading component, and $\mbox{\boldmath{$h$}}_{n}$ follows the independent and
identically distributed (i.i.d.) Rayleigh fading model, i.e.,
$\mbox{\boldmath{$h$}}_{n}\sim\mathcal{CN}(0,\mbox{\boldmath{$I$}})$, $\forall
n$. In addition, the devices are not perfectly synchronized since the devices
are equipped with low-cost local oscillators. We assume that the transmitted
packet of each active device $n\in\mathcal{K}$ is delayed by
$\tau_{n}\in\\{0,\dots,\tau_{{\rm max}}\\}$ symbols, where $\tau_{{\rm max}}$
is the maximum delay for all the devices and is assumed to be known at the BS.
The delay of each active device is unknown and needs to be estimated.

Each device is preassigned a unique preamble sequence
$\mbox{\boldmath{$s$}}_{n}\in\mathbb{C}^{L\times 1}$, $n=1,\dots,N$, where
$L<N$ is the preamble sequence length. Then, for active device
$n\in\mathcal{K}$, it starts to transmit sequence $\mbox{\boldmath{$s$}}_{n}$
at the beginning of the $\tau_{n}+1$ time slot. Note that in this letter, one
time slot equals one symbol duration. We use the two terms interchangeably.
Given the delay $\tau_{n}$, define the effective preamble sequence of device
$n$ as

$\displaystyle\bar{\mbox{\boldmath{$s$}}}_{n,\tau_{n}}=[\underbrace{0,\dots,0}_{\tau_{n}},\mbox{\boldmath{$s$}}_{n},\underbrace{0,\dots,0}_{\tau_{{\rm
max}}-\tau_{n}}]^{T},~{}n=1,\dots,N.$ (2)

In this case, the received signal
$\mbox{\boldmath{$Y$}}\in\mathbb{C}^{(L+\tau_{{\rm max}})\times M}$ from time
slot 1 to time slot $L+\tau_{{\rm max}}$ is expressed as

$\displaystyle\mbox{\boldmath{$Y$}}=\sum_{n=1}^{N}\alpha_{n}\sqrt{\beta_{n}}\bar{\mbox{\boldmath{$s$}}}_{n,\tau_{n}}\mbox{\boldmath{$h$}}^{T}_{n}+\mbox{\boldmath{$Z$}},$
(3)

where
$\mbox{\boldmath{$Z$}}=[\mbox{\boldmath{$z$}}_{1},\dots,\mbox{\boldmath{$z$}}_{M}]$
is the additive white Gaussian noise with
$\mbox{\boldmath{$z$}}_{m}\sim\mathcal{CN}(0,\sigma^{2}\mbox{\boldmath{$I$}})$,
$\forall m$. To express the received signal at the BS in a more compact form,
denote

$\theta_{n,\tau}=\left\\{\begin{array}[]{ll}1,~{}{\rm
if}~{}\alpha_{n}=1~{}{\rm and}~{}\tau=\tau_{n},\\\ 0,~{}{\rm
otherwise},\end{array}~{}\forall n,\tau.\right.$ (4)

In other words, $\theta_{n,\tau}=1$ only if device $n$ is active, and its
delay is of $\tau_{n}$ symbol durations. In addition, denote
$\mbox{\boldmath{$S$}}=\left[\bar{\mbox{\boldmath{$S$}}}_{1},\bar{\mbox{\boldmath{$S$}}}_{2},\dots,\bar{\mbox{\boldmath{$S$}}}_{N}\right]$,
where
$\bar{\mbox{\boldmath{$S$}}}_{n}=\left[\bar{\mbox{\boldmath{$s$}}}_{n,0},\bar{\mbox{\boldmath{$s$}}}_{n,1},\dots,\bar{\mbox{\boldmath{$s$}}}_{n,\tau_{{\rm
max}}}\right]$; $\mbox{\boldmath{$\gamma$}}={\rm
diag}\left(\mbox{\boldmath{$\gamma$}}_{1},\mbox{\boldmath{$\gamma$}}_{2},\dots,\mbox{\boldmath{$\gamma$}}_{N}\right)$,
where $\mbox{\boldmath{$\gamma$}}_{n}={\rm
diag}\left(\gamma_{n,0},\dots,\gamma_{n,\tau_{{\rm max}}}\right)$ and
$\gamma_{n,\tau}=\beta_{n}\theta_{n,\tau}$;
$\mbox{\boldmath{$H$}}=\left[\mbox{\boldmath{$H$}}_{1},\mbox{\boldmath{$H$}}_{2},\dots,\mbox{\boldmath{$H$}}_{N}\right]^{T}$,
where
$\mbox{\boldmath{$H$}}_{n}=[\underbrace{\mbox{\boldmath{$h$}}_{n},\mbox{\boldmath{$h$}}_{n},\dots,\mbox{\boldmath{$h$}}_{n}}_{\tau_{{\rm
max}}+1}]$. In this case, $Y$ in (3) is rewritten as

$\displaystyle\mbox{\boldmath{$Y$}}=\sum_{n=1}^{N}\sum_{\tau=0}^{\tau_{{\rm
max}}}\theta_{n,\tau}\sqrt{\beta_{n}}\bar{\mbox{\boldmath{$s$}}}_{n,\tau}\mbox{\boldmath{$h$}}^{T}_{n}+\mbox{\boldmath{$Z$}}=\mbox{\boldmath{$S$}}\mbox{\boldmath{$\gamma$}}^{\frac{1}{2}}\mbox{\boldmath{$H$}}+\mbox{\boldmath{$Z$}}.$
(5)

From the above received signal expression in (5), we model the joint activity
and delay detection problem as an effective sequence detection problem.
Specifically, we detect the effective sequence through the detection of
$\gamma$, which contains the device activity and delay information in its
diagonal elements $\theta_{n,\tau}$’s. For example, a non-zero diagonal
element $\gamma_{3,1}$ in $\gamma$ indicates that the effective sequence
$\bar{\mbox{\boldmath{$s$}}}_{3,1}\in\mbox{\boldmath{$S$}}$ is detected, which
means device $3$ is active, and its delay is of one symbol duration. In this
case, to estimate the device activity and delay information, we can estimate
$\gamma$ in (5) directly, instead of estimating $\alpha_{n}$’s and
$\tau_{n}$’s separately.

###  II-B Problem Formulation

In the following, we show the approach to estimate $\gamma$ based on the
received signal $Y$ and the effective sequence $S$ in (5). Once $\gamma$ is
estimated, we can determine the active devices and their corresponding delays
based on the diagonal values in the estimated $\gamma$. We formulate the
estimation of $\gamma$ as a maximum likelihood estimation problem.
Specifically, we first derive the likelihood of $Y$ given $\gamma$ and $S$ as

$\displaystyle
p(\mbox{\boldmath{$Y$}}|\mbox{\boldmath{$\gamma$}},\mbox{\boldmath{$S$}})$
$\displaystyle=\prod_{m=1}^{M}p(\mbox{\boldmath{$y$}}_{m}|\mbox{\boldmath{$\gamma$}},\mbox{\boldmath{$S$}})$
$\displaystyle\overset{(a)}{=}\frac{1}{|\mbox{\boldmath{$\Sigma$}}\pi|^{M}}\prod_{m=1}^{M}\exp{(-\mbox{\boldmath{$y$}}_{m}^{H}\mbox{\boldmath{$\Sigma$}}^{-1}\mbox{\boldmath{$y$}}_{m})}$
$\displaystyle\overset{(b)}{=}\frac{1}{|\mbox{\boldmath{$\Sigma$}}\pi|^{M}}\exp{\left(-{\rm
Tr}\left(\mbox{\boldmath{$\Sigma$}}^{-1}\mbox{\boldmath{$Y$}}\mbox{\boldmath{$Y$}}^{H}\right)\right)},$
(6)

where $\mbox{\boldmath{$y$}}_{m}$ is the $m$-th column of $Y$. In (6), (a) is
because given $\gamma$ and $S$, $\mbox{\boldmath{$y$}}_{m}$’s are independent
of each other, and (b) is because given $\gamma$ and $S$,
$\mbox{\boldmath{$y$}}_{m}$ is Gaussian distributed with zero mean and
covariance matrix

$\displaystyle\mbox{\boldmath{$\Sigma$}}=\mathbb{E}(\mbox{\boldmath{$y$}}_{m}\mbox{\boldmath{$y$}}_{m}^{H})$
$\displaystyle=\mbox{\boldmath{$S$}}\mbox{\boldmath{$\gamma$}}^{\frac{1}{2}}{\rm
diag}\left(\mbox{\boldmath{$E$}}_{1},\dots,\mbox{\boldmath{$E$}}_{N}\right)\mbox{\boldmath{$\gamma$}}^{\frac{1}{2}}\mbox{\boldmath{$S$}}^{H}+\sigma^{2}\mbox{\boldmath{$I$}}$
$\displaystyle=\mbox{\boldmath{$S$}}\mbox{\boldmath{$\gamma$}}\mbox{\boldmath{$S$}}^{H}+\sigma^{2}\mbox{\boldmath{$I$}},$
(7)

where $\mbox{\boldmath{$E$}}_{n}\in\mathbb{C}^{(\tau_{{\rm
max}}+1)\times(\tau_{{\rm max}}+1)}$ is a all-one matrix. From (7), we know
that the covariance matrix $\Sigma$ does not depend on the index $m$. Note
that maximizing
$p(\mbox{\boldmath{$Y$}}|\mbox{\boldmath{$\gamma$}},\mbox{\boldmath{$S$}})$ in
(6) is equivalent to minimizing
$-\log{(p(\mbox{\boldmath{$Y$}}|\mbox{\boldmath{$\gamma$}},\mbox{\boldmath{$S$}}))}$.
In this case, the optimization problem is formulated as

$\displaystyle\underset{\mbox{\boldmath{$\gamma$}}}{\rm
minimize}~{}~{}~{}\log{|\mbox{\boldmath{$\Sigma$}}|}+{\rm
Tr}\left(\mbox{\boldmath{$\Sigma$}}^{-1}\tilde{\mbox{\boldmath{$\Sigma$}}}\right)$
(8) $\displaystyle{\rm subject~{}to}~{}~{}~{}\mbox{\boldmath{$\gamma$}}\geq
0,$ (9)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}|\mbox{\boldmath{$\gamma$}}_{n}|_{0}\leq
1,~{}n=1,\dots,N,$ (10)

where
$\tilde{\mbox{\boldmath{$\Sigma$}}}=\frac{1}{M}\mbox{\boldmath{$Y$}}\mbox{\boldmath{$Y$}}^{H}$
is the sample covariance of the received signals;
$|\mbox{\boldmath{$\gamma$}}_{n}|_{0}$ denotes the number of nonzero elements
in the vector $\mbox{\boldmath{$\gamma$}}_{n}$. The constraint in (9) is due
to the non-negative path loss coefficients in the diagonal elements of
$\gamma$, and the constraint in (10) is because there is at most one possible
delay for each device. In addition, the estimation problem depends on the
received signals only through the sample covariance
$\tilde{\mbox{\boldmath{$\Sigma$}}}$. Note that, the sample covariance
$\tilde{\mbox{\boldmath{$\Sigma$}}}$ approaches the real covariance $\Sigma$
as the increase of $M$, which in turn leads to a better solution.

Algorithm 1 CD Algorithm followed by Constraint Enforcement for Solving
Problem (8)

1:Initialization: $\hat{\mbox{\boldmath{$\gamma$}}}=0$,
$\hat{\mbox{\boldmath{$\Sigma$}}}=\sigma^{2}\mbox{\boldmath{$I$}}$, $f_{0}$,
$\epsilon=1$, $i=0$.

2:while $\epsilon>\delta$ do

3: $i\leftarrow i+1$;

4: for $n=1,\dots,N$ do

5: for $\tau=0,\dots,\tau_{{\rm max}}$ do

6: Compute $\eta$ according to (18) on the next page;

7: $\hat{\gamma}_{n,\tau}\leftarrow\hat{\gamma}_{n,\tau}+\eta$;

8:
$\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}\leftarrow\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}-\eta\frac{\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}\mbox{\boldmath{$s$}}_{n,\tau}\mbox{\boldmath{$s$}}_{n,\tau}^{H}\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}}{1+\eta\mbox{\boldmath{$s$}}_{n,\tau}^{H}\mbox{\boldmath{$\Sigma$}}^{-1}\mbox{\boldmath{$s$}}_{n,\tau}}$;

9: end for

10: end for

11: Compute $f_{i}$;

12: $\epsilon=f_{i-1}-f_{i}$;

13:end while$\triangleright$ _End of the CD algorithm_.

14:Enforce constraint (10) on $\hat{\mbox{\boldmath{$\gamma$}}}$ according to
(13);

15:Output $\hat{\mbox{\boldmath{$\gamma$}}}$.

##  III Proposed Algorithms for Joint device Activity and Delay Detection

In this section, we propose two algorithms for solving the optimization
problem in (8). Note that the optimization problem in (8) is not a convex
problem. We first adopt the standard CD algorithm[12], and then propose a BCD
algorithm for solving the problem, respectively.

###  III-A CD Algorithm Followed by Constraint Enforcement

The optimization problem (8) can be efficiently solved by adopting the
standard CD algorithm if we drop the constraint (10) in the problem (8). In
this subsection, we first apply the standard CD algorithm to get an estimate
$\hat{\mbox{\boldmath{$\gamma$}}}$ of $\gamma$. Then, we enforce the block-
wise constraint in (10) on the estimate $\hat{\mbox{\boldmath{$\gamma$}}}$ as
follows:

$\displaystyle\hat{\gamma}_{n,\tau}=\left\\{\begin{array}[]{ll}\\!\\!\\!\\!\hat{\gamma}_{n,\tau},~{}{\rm
if}~{}\hat{\gamma}_{n,\tau}=\hat{\gamma}_{n}^{{\rm max}},\\\ 0,~{}~{}~{}{\rm
otherwise},\end{array}\\!\forall n,\tau,\right.$ (13)

where $\hat{\gamma}_{n}^{{\rm max}}=\max_{\tau=0}^{\tau=\tau_{{\rm
max}}}\\{\hat{\gamma}_{n,\tau}\\}$.

The CD followed by constraint enforcement (CD-E) algorithm is summarized in
Algorithm 1, where $i$ denotes the index of iteration, and $f_{i}$ denotes the
objective value of problem (8) at the $i$-th iteration. In particular, $f_{0}$
denotes the objective value when the initial estimation
$\hat{\mbox{\boldmath{$\gamma$}}}$ is set to be zero; $\delta$ is a small
number to control the convergence of the algorithm. After executing the CD-E
algorithm, we get an estimate $\hat{\mbox{\boldmath{$\gamma$}}}$. Then, we
need to do a thresholding procedure as follows:

$\displaystyle\hat{\gamma}_{n,\tau}=\left\\{\begin{array}[]{ll}\\!\\!\\!\\!\hat{\gamma}_{n,\tau},~{}{\rm
if}~{}\hat{\gamma}_{n,\tau}\geq t_{{\rm CD}},\\\ 0,~{}~{}~{}{\rm
otherwise},\end{array}\\!\forall n,\tau,\right.$ (16)

where $t_{{\rm CD}}$ is a threshold to determine the nonzero diagonal elements
in $\hat{\mbox{\boldmath{$\gamma$}}}$ since the estimated
$\hat{\mbox{\boldmath{$\gamma$}}}$ has been corrupted by noise. Last, the
activity and delay indicator $\theta_{n,\tau}$ in (4) is estimated as

$\hat{\theta}_{n,\tau}=\left\\{\begin{array}[]{ll}1,~{}{\rm
if}~{}\hat{\gamma}_{n,\tau}>0,\\\ 0,~{}{\rm
otherwise},\end{array}~{}~{}\forall n,\tau.\right.$ (17)

The indices of nonzero $\hat{\theta}_{n,\tau}$ are the active devices and
their corresponding delays.

$\eta=\max\left\\{\frac{\mbox{\boldmath{$s$}}_{n,\tau}^{H}\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}\tilde{\mbox{\boldmath{$\Sigma$}}}\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}\mbox{\boldmath{$s$}}_{n,\tau}-\mbox{\boldmath{$s$}}_{n,\tau}^{H}\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}\mbox{\boldmath{$s$}}_{n,\tau}}{\left(\mbox{\boldmath{$s$}}_{n,\tau}^{H}\hat{\mbox{\boldmath{$\Sigma$}}}^{-1}\mbox{\boldmath{$s$}}_{n,\tau}\right)^{2}},-\hat{\gamma}_{n,\tau}\right\\}$
(18)

###  III-B Block Coordinate Descent Algorithm

The CD-E algorithm presented in the last subsection applies the CD algorithm
to optimize the problem (8) by relaxing the constraint (10), and then enforces
the obtained solution to satisfy the constraint in (10). In this subsection,
we introduce a BCD algorithm which takes the constraint (10) into account
within its optimization process. Thus, the BCD algorithm may lead to a better
solution compared with the CD-E algorithm.

For the BCD algorithm, at each time, we only optimize
$\mbox{\boldmath{$\gamma$}}_{\bar{n}}$, while keeping
$\mbox{\boldmath{$\gamma$}}_{n}$ for $n\neq\bar{n}$ unchanged. In this case,
the optimization problem in (8) reduces to the following problem:

$\displaystyle\underset{\mbox{\boldmath{$\gamma$}}_{\bar{n}}}{\rm
minimize}~{}~{}~{}\log{|\mbox{\boldmath{$\Sigma$}}|}+{\rm
Tr}\left(\mbox{\boldmath{$\Sigma$}}^{-1}\tilde{\mbox{\boldmath{$\Sigma$}}}\right)$
(19) $\displaystyle{\rm
subject~{}to}~{}~{}~{}\mbox{\boldmath{$\gamma$}}_{\bar{n}}\geq 0,$ (20)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}|\mbox{\boldmath{$\gamma$}}_{\bar{n}}|_{0}\leq
1.$ (21)

We propose to solve the problem (19) by decomposing it into $\tau_{{\rm
max}}+1$ subproblems. Specifically, for the $\bar{\tau}$-th subproblem, we fix
$\gamma_{\bar{n},\tau}=0$ for $\tau\neq\bar{\tau}$, and optimize
$\gamma_{\bar{n},\bar{\tau}}$ by solving the following problem:

$\displaystyle\underset{\gamma_{\bar{n},\bar{\tau}}}{\rm
minimize}~{}~{}~{}\log{|\mbox{\boldmath{$\Sigma$}}|}+{\rm
Tr}\left(\mbox{\boldmath{$\Sigma$}}^{-1}\tilde{\mbox{\boldmath{$\Sigma$}}}\right)$
(22) $\displaystyle{\rm subject~{}to}~{}~{}\gamma_{\bar{n},\bar{\tau}}\geq 0.$
(23)

The above problem (22) is exactly a subproblem within the CD-E algorithm
(i.e., Line 6 to Line 8 in Algorithm 1) and has a closed-form solution. For
the $\bar{\tau}$-th subproblem, let $\hat{\gamma}_{\bar{n},\bar{\tau}}$ be the
solution to the subproblem (22),
$\mbox{\boldmath{$\gamma$}}_{\bar{n}}^{(\bar{\tau})}\in\mathbb{R}^{(\tau_{{\rm
max}}+1)\times 1}$ be a all-zero vector except for the $\bar{\tau}$-th element
being $\hat{\gamma}_{\bar{n},\bar{\tau}}$, and $f_{\bar{n},\bar{\tau}}$ be the
objective value of problem (22) at $\hat{\gamma}_{\bar{n},\bar{\tau}}$. After
solving subproblem (22) for $\bar{\tau}=0,\dots,\tau_{{\rm max}}$, the
solution for problem (19) is

$\displaystyle\mbox{\boldmath{$\gamma$}}_{\bar{n}}=\mbox{\boldmath{$\gamma$}}_{\bar{n}}^{(\hat{\tau})},~{}~{}{\rm
for}~{}~{}\hat{\tau}=\underset{\tau}{\arg\min}~{}f_{\bar{n},\tau}.$ (24)

The problem (19) is solved from $\bar{n}=1$ to $\bar{n}=N$. The BCD algorithm
is summarized in Algorithm 2. After executing the BCD algorithm, we get an
estimate $\hat{\mbox{\boldmath{$\gamma$}}}$. Then, we determine the nonzero
diagonal elements in $\hat{\mbox{\boldmath{$\gamma$}}}$ according to (16)
except that the threshold now is $t_{{\rm BCD}}$. Last, the activity and delay
indicator $\theta_{n,\tau}$ in (4) is similarly estimated as in (17).

In Section IV, we show that the BCD algorithm outperforms the CD-E algorithm.
The possible reasons are as follows. First, BCD solves the original problem
while CD-E solves the relaxation problem albeit with a constraint enforcement
procedure. Second, the monotonicity of the objective values in BCD is
guaranteed. However, the constraint enforcement procedure in CD-E may destroy
the monotonicity of the objective function, and possibly increases the
objective value without any control.

Algorithm 2 BCD Algorithm for Solving Problem (8)

1:Initialization: $\hat{\mbox{\boldmath{$\gamma$}}}=0$, $f_{0}$, $\epsilon=1$,
$i=0$.

2:while $\epsilon>\delta$ do

3: $i\leftarrow i+1$;

4: for $\bar{n}=1,\dots,N$ do

5: for $\bar{\tau}=0,\dots,\tau_{{\rm max}}$ do

6: Set $\hat{\gamma}_{\bar{n},\tau}=0$ for $\tau\neq\bar{\tau}$;

7: Compute $\eta$ according to (18);

8:
$\hat{\gamma}_{\bar{n},\bar{\tau}}\leftarrow\hat{\gamma}_{\bar{n},\bar{\tau}}+\eta$;

9: Compute $f_{\bar{n},\bar{\tau}}$;

10: end for

11: Determine $\mbox{\boldmath{$\gamma$}}_{\bar{n}}$ according to (24);

12: end for

13: Compute $f_{i}$;

14: $\epsilon=f_{i-1}-f_{i}$;

15:end while

16:Output $\hat{\mbox{\boldmath{$\gamma$}}}$.

##  IV Numerical Results

In this section, we present the detection performance of the proposed CD-E and
BCD algorithms in terms of missed detection probability (MDP) and false alarm
probability (FAP). The MDP is the probability that an active device is
detected to be inactive, or an active device is correctly detected but its
delay is not correctly detected; the FAP is the probability that an inactive
device is detected to be active. The number of potential devices is $N=200$,
the number of active devices is $K=90$, the preamble sequence length is
$L=100$, and the maximum number of delays is $\tau_{{\rm max}}=4$ symbol
durations. All preamble sequences follow the i.i.d. complex Gaussian
distribution with zero mean and unit variance. The noise power spectrum
density is -169dBm/Hz, the bandwidth is 10 MHz, the transmit power of each
device is set to be 23dBm, and the path loss is modeled as
$128.1+37.6\log_{10}(d)$, where $d$ is the distance between the BS and devices
in km[23]. For ease of demonstration, we consider a worst case in which all
devices are located in the cell edge such that the large-scale fading
coefficients are the same for all devices. In the simulations, we set $d=1$
km. Moreover, we set $\delta=10^{-3}$, $t_{{\rm CD}}=0.1$, and $t_{{\rm
BCD}}=0.12$. Note that $t_{{\rm BCD}}$ is slightly larger than $t_{{\rm CD}}$
since we observed that the estimated $\hat{\mbox{\boldmath{$\gamma$}}}$ by the
BCD algorithm is slightly larger than that by the CD-E algorithm. Last, we
have two benchmarks to evaluate the performance of the proposed algorithms as
follows:

  * •

CS approach in [21]. Ref. [21] also studied the joint device activity and
delay detection problem. Different from us, [21] solved this problem by
applying a CS approach.

  * •

CD-E with $\tau_{{\rm max}}=0$ and $L=104$. We consider a case in which the
delay is zero, and the sequence length is $L=104$ which is the same as the
effective sequence length of the case with $\tau_{{\rm max}}=4$ and $L=100$.
In this case, without delay, the CD-E algorithm only needs to detect the
device activity, and the CD-E algorithm does not need to execute the
constraint enforcement procedure in (13). Thus, CD-E with $\tau_{{\rm max}}=0$
and $L=104$ provides a performance upper bound of CD-E with $\tau_{{\rm
max}}=4$ and $L=100$. In addition, the BCD algorithm is the same as the CD-E
algorithm when $\tau_{{\rm max}}=0$.

Figure 1: MDP comparison over different numbers of antennas.
Figure 2: FAP comparison over different numbers of antennas.

We compare their performance in terms of MDP and FAP over different numbers of
BS antennas in Fig. 1 and Fig. 2, respectively. First, when $\tau_{{\rm
max}}=4$ and $L=100$, CD-E and BCD outperform CS significantly both in MDP and
FAP. A key intuition is that, the covariance-based approach can make better
use of the BS antennas such that the sample covariance of the received signals
approaches the real covariance as the increase of the number of BS antennas,
which in turn leads to a better solution. Second, when $\tau_{{\rm max}}=4$
and $L=100$, BCD performs slightly better than CD-E. Specifically, from Fig.
1, they have the same MDP performance, while from Fig. 2, BCD performs better
than CD-E in FAP when $M$ is larger than 8. Third, in terms of MDP, “CD-E/BCD:
$\tau_{{\rm max}}=4$ and $L=100$” perform roughly the same as “CD-E:
$\tau_{{\rm max}}=0$ and $L=104$”, which is the performance limit of
“CD-E/BCD: $\tau_{{\rm max}}=4$ and $L=100$”; in terms of FDP, the performance
gap between “CD-E/BCD: $\tau_{{\rm max}}=4$ and $L=100$” and “CD-E:
$\tau_{{\rm max}}=0$ and $L=104$” is small.

## References

  * [1] C. Bockelmann, N. Pratas, H. Nikopour, K. Au, T. Svensson, C. Stefanovic, P. Popovski, and A. Dekorsy, “Massive machine-type communications in 5G: Physical and MAC-layer solutions,” _IEEE Commun. Mag._ , vol. 54, no. 9, pp. 59–65, Sept. 2016.
  * [2] X. Chen, D. W. K. Ng, W. Yu, E. G. Larsson, N. Al-Dhahir, and R. Schober, “Massive access for 5G and beyond,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 3, pp. 615–637, Mar. 2020.
  * [3] L. Liu, E. G. Larsson, W. Yu, P. Popovski, C. Stefanovic, and E. De Carvalho, “Sparse signal processing for grant-free massive connectivity: A future paradigm for random access protocols in the Internet of Things,” _IEEE Signal Process. Mag._ , vol. 35, no. 5, pp. 88–99, Sept. 2018.
  * [4] K. Senel and E. G. Larsson, “Grant-free massive MTC-enabled massive MIMO: A compressive sensing approach,” _IEEE Trans. Commun._ , vol. 66, no. 12, pp. 6164–6175, Dec. 2018.
  * [5] L. Liu and W. Yu, “Massive connectivity with massive MIMO—Part I: Device activity detection and channel estimation,” _IEEE Trans. Signal Process._ , vol. 66, no. 11, pp. 2933–2946, Jun. 2018.
  * [6] Z. Chen, F. Sohrabi, and W. Yu, “Sparse activity detection for massive connectivity,” _IEEE Trans. Signal Process._ , vol. 66, no. 7, pp. 1890–1904, Apr. 2018.
  * [7] T. Jiang, Y. Shi, J. Zhang, and K. B. Letaief, “Joint activity detection and channel estimation for IoT networks: Phase transition and computation-estimation tradeoff,” _IEEE Internet of Things J._ , vol. 6, no. 4, pp. 6212–6225, Aug. 2018.
  * [8] M. Ke, Z. Gao, Y. Wu, X. Gao, and R. Schober, “Compressive sensing-based adaptive active user detection and channel estimation: Massive access meets massive MIMO,” _IEEE Trans. Signal Process._ , vol. 68, pp. 764–779, Jan. 2020.
  * [9] T. Ding, X. Yuan, and S. C. Liew, “Sparsity learning-based multiuser detection in grant-free massive-device multiple access,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 7, pp. 3569–3582, Jul. 2019.
  * [10] Z. Sun, Z. Wei, L. Yang, J. Yuan, X. Cheng, and L. Wan, “Exploiting transmission control for joint user identification and channel estimation in massive connectivity,” _IEEE Trans. Commun._ , vol. 67, no. 9, pp. 6311–6326, Sept. 2019.
  * [11] Y. Mei, Z. Gao, Y. Wu, W. Chen, J. Zhang, D. W. K. Ng, and M. Di Renzo, “Compressive sensing based joint activity and data detection for grant-free massive IoT access,” _IEEE Trans. Wireless Commun._ , Sept. 2021.
  * [12] S. Haghighatshoar, P. Jung, and G. Caire, “Improved scaling law for activity detection in massive MIMO systems,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , June 2018, pp. 381–385.
  * [13] Z. Chen, F. Sohrabi, Y.-F. Liu, and W. Yu, “Covariance based joint activity and data detection for massive random access with massive MIMO,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , May 2019, pp. 1–6.
  * [14] ——, “Phase transition analysis for covariance based massive random access with massive MIMO,” to appear in _IEEE Trans. Inf. Theory_ , 2021.
  * [15] Z. Wang, Z. Chen, Y.-F. Liu, F. Sohrab, and W. Yu, “An efficient active set algorithm for covariance based joint data and activity detection for massive random access with massive MIMO,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Process. (ICASSP)_ , Jun. 2021.
  * [16] U. K. Ganesan, E. Björnson, and E. G. Larsson, “An algorithm for grant-free random access in cell-free massive MIMO,” in _Proc. IEEE Workshop Signal Process. Adv. Wireless Commun. (SPAWC)_ , May 2020, pp. 1–5.
  * [17] Z. Chen, F. Sohrabi, and W. Yu, “Sparse activity detection in multi-cell massive MIMO exploiting channel large-scale fading,” _IEEE Trans. Signal Process._ , vol. 69, pp. 3768 – 3781, Jun. 2021.
  * [18] Z. Wang, Y.-F. Liu, Z. Chen, and W. Yu, “Accelerating coordinate descent via active set selection for device activity detection for multi-cell massive random access,” in _Proc. IEEE Workshop Signal Process. Adv. Wireless Commun. (SPAWC)_ , Sept. 2021.
  * [19] L. Liu and W. Yu, “Massive connectivity with massive MIMO—Part II: Achievable rate characterization,” _IEEE Trans. Signal Process._ , vol. 66, no. 11, pp. 2947–2959, Jun. 2018.
  * [20] J. Kang and W. Yu, “Minimum feedback for collision-free scheduling in massive random access,” _IEEE Trans. Inf. Theory_ , vol. 67, no. 12, pp. 8094 – 8108, Dec. 2021.
  * [21] L. Liu and Y.-F. Liu, “An efficient algorithm for device detection and channel estimation in asynchronous IoT systems,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Process. (ICASSP)_ , Jun. 2021.
  * [22] M. Ke, Z. Gao, Y. Wu, X. Gao, and K.-K. Wong, “Massive access in cell-free massive MIMO-based Internet of Things: Cloud computing and edge computing paradigms,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 3, pp. 756–772, Aug. 2020.
  * [23] 3GPP, “Evolved universal terrestrial radio access (E-UTRA); radio frequency (RF) system scenarios,” Technical Specification (TS) 36.942.


# Doping phase diagram of a Hubbard model for twisted bilayer cuprates

Xiancong Lu  Department of Physics, Xiamen University, Xiamen 361005, China
D. Sénéchal  Département de physique and Institut quantique, Université de
Sherbrooke, Sherbrooke, Québec, Canada J1K 2R1

###### Abstract

We study the twisted Hubbard model of a cuprate bilayer at a fixed twist angle
$\theta=53.13^{\circ}$ using the variational cluster approximation, a method
that treats short-range dynamical correlations exactly. At intermediate
interlayer tunneling, the phase difference $\phi$ between the $d$-wave order
parameters of two layers is $\pi$ in the overdoped regime, while it is zero in
the underdoped regime, close to the Mott phase. At strong interlayer
tunneling, we observe a clear time-reversal symmetry breaking phase near
optimal doping, in which the phase difference $\phi$ changes continuously from
0 to $\pi$. However, this phase has trivial topology. We also apply a cluster
extension of dynamical mean field theory to the same problem, but fail to
detect a time-reversal breaking phase with that method.

##  I Introduction

The experimental discovery of correlated insulators and unconventional
superconductivity in twisted bilayer graphene (TBG) Cao  _et al._ (2018a, b)
has opened up the new field of twistronics Andrei and MacDonald (2020); Carr
_et al._ (2017). By twisting two graphene sheets by a small relative angle, a
long-period moiré pattern forms in the bilayer. At special magic angles, the
moiré band structure of TBG exhibits isolated flat bands near charge
neutrality Bistritzer and MacDonald (2011); Lopes dos Santos  _et al._ (2007);
Tarnopolsky  _et al._ (2019), which lead to a variety of strongly correlated
phenomena. Following this discovery, various twisted van der Waals
heterostructures have been constructed and investigated Kennes  _et al._
(2021), including transition metal dichalcogenides Regan  _et al._ (2020);
Tang  _et al._ (2020); Wang  _et al._ (2020); Huang  _et al._ (2021), double
bilayer graphene Liu  _et al._ (2020); Cao  _et al._ (2020); Zhang  _et al._
(2021), and trilayer graphene Chen  _et al._ (2019a, b, 2020); Park  _et al._
(2021).

Recently, twistronics concepts have been extended to high-temperature
superconductors Can  _et al._ (2021); Volkov  _et al._ (2020), which are
strongly correlated materials by themselves. This was motivated by the
experimental realization of two-dimensional (2D) monolayer Bi2Sr2CaCu2O8+δ
(Bi2212), whose transition temperature is shown to be very close to that of
bulk samples Yu  _et al._ (2019); Zhao  _et al._ (2019). It is theoretically
predicted that, at large twist angles (close to $45^{\circ}$), a fully gapped
$d+id$ superconducting phase emerges, which spontaneously breaks time-reversal
symmetry (TRS) and is topologically nontrivial Can  _et al._ (2021). This TRS
breaking superconducting phase is also predicted to be stable at small twist
angle, due to the strong renormalization of Bogoliubov-de Gennes (BdG)
quasiparticles near the nodes Volkov  _et al._ (2020). In order to determine
the pairing symmetry of cuprate superconductors, $c$-axis twisted Josephson
junctions, formed by stacking two Bi2212 crystals along the $c$-axis, have
been realized Li  _et al._ (1999); Takano  _et al._ (2002); Latyshev  _et al._
(2004); Klemm (2005); Yang  _et al._ (2018); Zhu  _et al._ (2021). However,
most experimental works did not observe the angular dependence of the
Josephson current Li  _et al._ (1999); Takano  _et al._ (2002); Zhu  _et al._
(2021). Owed to the novel technique of van der Waals stacking, high-quality
twisted Bi2212 Josephson junctions with an atomically sharp interface have
been successfully fabricated recently Zhu  _et al._ (2021); Zhao  _et al._
(2021).

Previous theoretical work on twisted bilayer cuprates are mainly based on
Bogoliubov-de-Gennes mean-field theory Yang  _et al._ (2018); Can  _et al._
(2021); Volkov  _et al._ (2020); Tummuru  _et al._ (2021); Volkov  _et al._
(2021), which does not take into account the effects of strong correlations.
To overcome this, a twisted $t$-$J$ model of cuprates has been proposed and
studied within slave-boson mean-field theory Song  _et al._ (2021), in which a
topological-trivial time-reversal symmetry breaking superconductor is also
found, but within a small range of twist angles around $45^{\circ}$,
questioning the possibility of topological superconductors in this region. In
spite of this work, the stability of the novel superconducting phases against
doping has not been fully addressed before in the literature. In this paper,
we will numerically study the twisted Hubbard model of bilayer cuprates using
the variational cluster approach (VCA) and cluster dynamical mean field theory
(CDMFT). These approaches have been successfully used in the past to study
high-temperature superconductors and the Hubbard model at intermediate
coupling is arguably a better representation of these materials. We will focus
on a fixed twisted angle $\theta=53.13^{\circ}$, at which these cluster
methods are easily applicable, and investigate the superconducting phase
diagram as a function of doping for two different sets of interlayer
tunneling.

This paper is organized as follows. In Sec. II, we introduce the Hubbard model
for the twisted bilayer. In Sec. III, we review the variational cluster
approximation (VCA) and present our main results obtained from this method,
e.g., the phase diagram of bilayer as a function of hole doping. In Sect. IV,
we present the corresponding results from cluster dynamical mean field theory
(CDMFT).

##  II Model

###  II.1 Hamiltonian

We assume that each of the two layers of the system can be described by the
one-band Hubbard model (the sites correspond to the location of copper atoms).
The bilayer is then described by the following tight-binding Hubbard model:Can
_et al._ (2021); Song  _et al._ (2021)

$H=H^{(1)}+H^{(2)}+H_{\perp}.$ (1)

where the intra-layer Hamiltonian $H^{(\ell)}$ is

$H^{(\ell)}=\sum_{\mathbf{r},\mathbf{r}^{\prime}\in\ell,\sigma}t_{\mathbf{r}\mathbf{r}^{\prime}}c_{\mathbf{r},\ell,\sigma}^{\dagger}c_{\mathbf{r}^{\prime},\ell,\sigma}+U\sum_{\mathbf{r}}n_{\mathbf{r},\ell,\uparrow}n_{\mathbf{r},\ell,\downarrow}-\mu\sum_{\mathbf{r},\sigma}n_{\mathbf{r},\ell,\sigma},$
(2)

where $c_{\mathbf{r},\ell,\sigma}(c_{\mathbf{r},\ell,\sigma}^{\dagger})$ is
the annihilation (creation) operator of an electron at site $\mathbf{r}$ on
layer $\ell$ ($\ell=1,2$) with spin $\sigma=\uparrow,\downarrow$, and
$n_{\mathbf{r},\ell,\sigma}$ is the associated number density operator. The
labels $\mathbf{r},\mathbf{r}^{\prime}$ run over the possible sites of a
square lattice (each layer has it own). We will keep nearest-neighbor ($t$)
and next-nearest-neighbor ($t^{\prime}$) hopping terms only, so that the
dispersion relation on a square lattice is $\varepsilon(\mathbf{k})=-2t(\cos
k_{x}+\cos k_{y})+4t^{\prime}\cos k_{x}\cos k_{y}-\mu$. Only on-site
interactions are considered here. For Bi2212, the nearest-neighbor hopping is
$t=126$meV Markiewicz  _et al._ (2005). In the remainder of this paper, we set
$t$ as the energy unit, and choose the other parameters to be
$t^{\prime}=-0.3$ and $U=8$.

set |  $V_{1}$ |  $V_{2}$ |  $V_{3}$
---|---|---|---
I |  0.1 |  0.05 |  0.03
II |  0.4 |  0.2 |  0.012
Table 1: The two sets of inter-layer hopping terms used in this work.

The interlayer tunneling is represented by

$H_{\perp}=\sum_{n=1}^{3}V_{n}\sum_{\langle\mathbf{r},\mathbf{r}^{\prime}\rangle_{\bot,n},\sigma}\Big{[}c_{\mathbf{r},1,\sigma}^{\dagger}c_{\mathbf{r}^{\prime},2,\sigma}+\mbox{H.c.}\Big{]}$
(3)

where the notation $\langle\mathbf{r},\mathbf{r}^{\prime}\rangle_{\bot,n}$
($n=1,2,3$) stands for the set of square lattice sites $\mathbf{r}$ on layer 1
and $\mathbf{r}^{\prime}$ on layer 2 such that their projection on the plane
are $n^{\rm th}$ neighbors. This is illustrated on Fig. 1 for $V_{1}$, $V_{2}$
and $V_{3}$. For instance, $V_{1}$ is the interlayer tunneling between sites
located exactly on top of each other, $V_{2}$ for sites that are first
neighbors when projected on a common plane, etc. Such an interlayer tunneling
model is obviously oversimplified, as it ignores the complexity of the CuO2
layers and of the rare-earth layers that will intervene between the twisted
CuO2 layers. In this work we will use two sets of values for $V_{n}$, shown in
Table 1. These values have been chosen heuristically, those of set II being
four times larger than those of set I and certainly unrealistic, but necessary
in order to unravel TRS breaking, as we will see below. Given the accepted
values of hopping along the $c$-axis in bulk cuprates, even the values of set
I are large, and will be referred to as intermediate tunneling, whereas those
of set II will constitute strong tunneling.

Figure 1: (Color online). Unit cell of the bilayer twisted by an angle $\theta=2\arctan\frac{1}{2}=53.13^{\circ}$. The ten sites within the unit cell are labeled and their color (orange or blue) indicates the layer. The three most important interlayer tunneling terms ($V_{1,2,3}$) are illustrated in gray. (the red dashed enclosures are the clusters used in VCA; see below).  Figure 2: (Color online) Symmetries on the bilayer system. The bases of the two layers are indicated in blue and red, respectively. The rotations $C_{x}$, $C_{y}$ and $C_{d}$ are indicated; the rotations $C_{2}$ and $C_{4}$, within each plane, are not. Table 2:  Character table of $D_{4}$, with a list of the simplest gap basis functions. The rotation $C_{d}^{\prime}$ is defined about the other diagonal axis, at right angle from $C_{d}$.  | $\displaystyle e$ | $\displaystyle 2C_{4}$ | $\displaystyle C_{2}$ | $\displaystyle C_{x,y}$ | $\displaystyle C_{d,d^{\prime}}$ | gap functions
---|---|---|---|---|---|---
$\displaystyle A_{1}$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle 1$
$\displaystyle A_{2}$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle-1$ | $\displaystyle-1$ | $\displaystyle\sin k_{x}\;\sin k_{y}(\cos k_{x}-\cos k_{y})$
$\displaystyle B_{1}$ | $\displaystyle 1$ | $\displaystyle-1$ | $\displaystyle 1$ | $\displaystyle 1$ | $\displaystyle-1$ | $\displaystyle\cos k_{x}-\cos k_{y}$
$\displaystyle B_{2}$ | $\displaystyle 1$ | $\displaystyle-1$ | $\displaystyle 1$ | $\displaystyle-1$ | $\displaystyle 1$ | $\displaystyle\sin k_{x}\;\sin k_{y}$
$\displaystyle E$ | $\displaystyle 2$ | $\displaystyle 0$ | $\displaystyle-2$ | $\displaystyle 0$ | $\displaystyle 0$ | $\displaystyle(\sin k_{x},\sin k_{y})$

In order to simplify as much as possible our numerical work, we will restrict
our analysis to a twist angle of $\theta=2\arctan\frac{1}{2}=53.13^{\circ}$.
The unit cell of the twisted bilayer at that angle is illustrated on Fig. 1
and contains ten sites (five per layer).

###  II.2 Symmetries

The bilayer system is invariant under a $\pi/4$ rotation around the $z$ axis
(perpendicular to the bilayer plane) and under the $\pi$ rotations $C_{x}$,
$C_{y}$ and $C_{d}$ illustrated on Fig. 2, which make up the $D_{4}$ point
group, the same as for an isolated layer. Possible superconducting gap
functions for this system should in principle be classified according to the
irreducible representations of $D_{4}$. Table 2 shows the character table and
the simplest gap functions associated with each irreducible representation.
Representations $B_{1}$ and $B_{2}$ correspond to what is usually called
$d_{x^{2}-y^{2}}$ and $d_{xy}$, respectively. Representations $A_{1}$ and
$A_{2}$ correspond respectively to $s$-wave (or extended $s$-wave) and
$f$-wave, and the two-dimensional representation $E$ would correspond to
(triplet) $p$-wave, with basis $(p_{x},p_{y})$. Thus, the only possibility of
a (pure) chiral representation is $p_{x}+ip_{y}$, a triplet state that will
not occur in this cuprate system. We rather expect representations $B_{1}$ and
$B_{2}$ to be realized here, owing to the $d$-wave character of
superconductivity in single layers. In principle, according to the Landau
theory of phase transitions, one of those two should prevail just below
$T_{c}$, but there is always the possibility that, the two states ($B_{1}$ and
$B_{2}$) being very close in energy, a second phase transition occurs below
$T_{c}$ and a complex combination $d_{x^{2}-y^{2}}+id_{xy}$ is present at zero
temperature. This is the scenario anticipated in Ref. Can  _et al._ (2021) and
investigated here.

Figure 3: (Color online) Schematics of cluster methods. Left panel: the
lattice is tiled into identical $2\times 2$ clusters with Hamiltonian
$H^{\prime}$. Right panel: the reduced Brillouin zone (RBZ) is then four times
smaller than the original Brillouin zone.

##  III Results from the Variational Cluster Approximation

High-temperature superconductors have strong correlations. There is a limited
number of numerical methods that can tackle such systems, and methods based on
small clusters of sites embedded into an effective medium are amongst the most
successful. These so-called quantum cluster methods are approximation
strategies for the electron Green function $G(\mathbf{k},\omega)$, by which
the electron self-energy $\Sigma$ on the infinite lattice is approximated by
that of a small cluster. In this work we will apply two of these methods to
the bilayer Hamiltonian (1).

The first of these methods is the variational cluster approximation (VCA)
Potthoff (2003); Potthoff  _et al._ (2003); Potthoff (2014). It is based on a
variational principle proposed by Potthoff Potthoff (2003) and can be seen as
a variational extension of cluster perturbation theory (CPT) Sénéchal  _et
al._ (2000); Gros and Valenti (1993). Let us start by briefly summarizing the
latter. In CPT, the lattice is tiled into identical clusters, and the
Hamiltonian is written as $H=H^{\prime}+V$, where $H^{\prime}$ is the
restriction of $H$ to the clusters and $V$ only contains hopping terms between
different clusters. If the model contains $N_{b}$ bands and each cluster
contains $L$ lattice sites, then $LN_{b}$ must be small enough to allow for an
exact numerical solution of $H^{\prime}$, and the associated one-particle
Green function $\mathbf{G}_{c}(\omega)$ on the cluster is a $2LN_{b}\times
2LN_{b}$ matrix (the factor of 2 because of spin). The tiling into clusters
defines a superlattice, and the corresponding Brillouin zone is $L$ times
smaller than the original Brillouin zone (see Fig. 3). We call it the reduced
Brillouin zone and its wave vectors are noted $\mathbf{\tilde{k}}$. The
hopping matrix in $H$ can be expressed as a $2LN_{b}\times 2LN_{b}$ matrix
$\mathbf{t}(\mathbf{\tilde{k}})$, a function of $\mathbf{\tilde{k}}$, which is
the sum of a $\mathbf{\tilde{k}}$-independent part $\mathbf{t}_{c}$ and of the
inter-cluster part $\mathbf{V}(\mathbf{\tilde{k}})$ :
$\mathbf{t}(\mathbf{\tilde{k}})=\mathbf{t}_{c}+\mathbf{V}(\mathbf{\tilde{k}})$.
The self-energy $\bm{\Sigma}_{c}(\omega)$ associated with the cluster Green
function $\mathbf{G}_{c}(\omega)$ is thus defined by Dyson’s equation on the
cluster:

$\mathbf{G}_{c}^{-1}=\omega-\mathbf{t}_{c}-\bm{\Sigma}_{c}(\omega)$ (4)

In CPT, the electron self-energy is approximated by that of the restriction
$H^{\prime}$ of the Hamiltonian to the cluster. In the mixed momentum-cluster
site basis, the electron Green function is then given by the following
relation:

$\mathbf{G}^{-1}(\mathbf{\tilde{k}},\omega)=\omega-\mathbf{t}(\mathbf{\tilde{k}})-\bm{\Sigma}_{c}(\omega)=\mathbf{G}_{c}^{-1}(\omega)-\mathbf{V}(\mathbf{\tilde{k}})$
(5)

We assume here that the chemical potential $\mu$ is included in the hopping
matrix $\mathbf{t}(\mathbf{\tilde{k}})$.

CPT is unable to describe broken symmetry states: it is not a self-consistent
approach, nor is it based on a variational principle. The VCA adds a
variational aspect to CPT: the cluster Hamiltonian $H^{\prime}$ is augmented
by a certain number of Weiss fields:

$H^{\prime}\to H^{\prime}+\sum_{a}h_{a}{\hat{O}}_{a}$ (6)

where the operators ${\hat{O}}_{a}$ are defined on the cluster only, and
possibly represent broken symmetries. These additional terms are in turn
subtracted from $V$, so that the original Hamiltonian $H$ in unaffected. The
values $h_{a}$ of these Weiss fields are not arbitrary, but set by Potthoff’s
variational principle: The following function:

$\Omega(h_{a})=\Omega^{\prime}-\int\frac{d\omega}{2\pi}\sum_{\mathbf{\tilde{k}}}\ln\det\left[\mathbf{1}-\mathbf{V}(\mathbf{\tilde{k}})\mathbf{G}_{c}(\omega)\right]$
(7)

should be stationary with respect to these fields $h_{a}$. In that expression,
$\Omega^{\prime}$ is the ground state energy of the cluster Hamiltonian
$H^{\prime}$ and $\mathbf{G}_{c}(\omega)$ is the electron Green function
derived from the cluster Hamiltonian $H^{\prime}$ that includes the Weiss
fields $h_{a}{\hat{O}}_{a}$.

In the problem at hand, it might seem natural to use the 10-site unit cell
shown in Fig. 1 as the repeated cluster, especially since 10 sites is an
easily manageable size for an exact-diagonalization solver. However, the set
of numbered sites in Fig. 1 does not have the $D_{4}$ symmetry of the full
Hamiltonian, and this complicates the VCA computations. We will rather use a
slight refinement of the method described above, assuming that the repeated
unit is a supercluster of 10 sites obtained by assembling an octagonal cluster
of 8 sites and a point-like cluster of 2 sites, each delimited by a red dashed
line in Fig. 1. The self-energy of the supercluster is then a direct sum of
the self-energies of an 8-site and of a 2-site cluster. Otherwise, the method
is unchanged from the general approach described above.

Figure 4: (Color online) Top panel: Order parameter as a function of electron
density $n$, as obtained in VCA, for both inter-layer hopping sets. Lower
panel: the corresponding relative phase $\phi$ of the order parameters on the
two layers. At intermediate inter-layer hopping, the system jumps from
$\phi=0$ to $\phi=\pi$ at $n\approx 0.93$. At strong inter-layer hopping, the
switch is gradual.

On the octagonal cluster, we will define Weiss fields associated with $d$-wave
superconductivity on each layer. On each layer of the lattice, we can define
an operator field that describes $d$-wave superconductivity:

$\displaystyle\Delta^{(\ell)}=\sum_{\mathbf{r}\in\ell}$
$\displaystyle\Big{\\{}c_{\mathbf{r},\ell,{\uparrow}}c_{\mathbf{r}+\mathbf{x}^{(\ell)},\ell,\downarrow}-c_{\mathbf{r},\ell,\downarrow}c_{\mathbf{r}+\mathbf{x}^{(\ell)},\ell,{\uparrow}}$
$\displaystyle-
c_{\mathbf{r},\ell,{\uparrow}}c_{\mathbf{r}+\mathbf{y}^{(\ell)},\ell,\downarrow}+c_{\mathbf{r},\ell,\downarrow}c_{\mathbf{r}+\mathbf{y}^{(\ell)},\ell,{\uparrow}}\Big{\\}}$
(8)

where $\mathbf{x}^{(\ell)}$ and $\mathbf{y}^{(\ell)}$ are the orthogonal
lattice vectors on layer $\ell$. We can then add the following combinations to
the cluster Hamiltonian:

$H^{\prime}\to
H^{\prime}+\sum_{\ell=1,2}d^{(\ell)}\Delta^{(\ell)}_{c}+\mbox{H.c}$ (9)

where $\Delta^{(\ell)}_{c}$ is a restriction to the cluster of the lattice
operator (III) and $d^{(\ell)}$ is a complex amplitude. The real and imaginary
parts of $d^{(\ell)}$ are then Weiss fields in the sense of the coefficients
$h_{a}$ of Eq. (6). Because of overall phase symmetry, one can always assume
that $d^{(1)}$ is real, but we must assume in all generality that $d^{(2)}$ is
complex. The complex phase of $d^{(2)}$ is then the relative phase $\phi$ of
the superconducting order parameters of the two layers, and a value other than
zero or $\pi$ would signal a spontaneous TRS breaking and possible topological
properties.

Figure 5: (Color online) Potthoff functional as a function of inter-layer
phase $\phi$, for different values of the chemical potential, for intermediate
(top) and strong (bottom) inter-layer tunneling.

We applied the VCA method on this cluster system, using the two sets of
interlayer tunneling defined in Table 1. In practice, this means computing the
cluster Green function $\mathbf{G}_{c}(\omega)$ repeatedly while adjusting the
Weiss fields $d^{(\ell)}$ so as to make the Potthoff functional stationary (in
fact, minimum). Once the stationary values are found, the Green function (5)
can be used to compute the ground state average of any one-body operator, in
particular the order parameter $\Psi^{(\ell)}=\langle\Delta^{(\ell)}\rangle/N$
($N$ is the number of sites) on each layer. The electron density $n$ can be
likewise computed from the Green function (the chemical potential $\mu$ is the
actual control parameter that is varied).

Fig. 4 shows the order parameter $\Psi^{(\ell)}$ as a function of electron
density $n$ for hole doping and the two sets of interlayer tunneling
(intermediate and strong). We note the characteristic dome shape that is
typically obtained in quantum cluster methods, qualitatively agreeing with the
known properties of cuprates. The electron density computed from the Green
function (5) has some systematic error, as can be seen from the fact that the
order parameter vanishes not at $n=1$, as it should from Mott physics, but at
$n=1.006$. The bottom layer of the figure shows the relative phase of the
order parameters $\Psi^{(2)}$ and $\Psi^{(1)}$ (on the two layers). At
intermediate inter-layer tunneling, this phase is 0 at low doping, which is
the signature of the $B_{1}$ representation of Table 2. Beyond about 7%
doping, this phase switches to $\pi$, a signature of the $B_{2}$
representation. There is thus a doping-induced transition of the bilayer
superconducting state, which coincides with the passage from underdoped to
overdoped, judging by the location of optimal doping on the upper panel of the
figure.

For strong interlayer tunneling, the situation is different: an intermediate
phase appears in which the relative complex phase of the two order parameters
changes continuously from 0 to $\pi$. This intermediate phase breaks time
reversal symmetry and correponds roughly to a $B_{1}+iB_{2}$ state, except
that the two components do not have the same amplitude.

A more detailed view of how this is happening from the VCA perspective is
shown on Fig. 5. On the top half of the figure, we show the profile of the
Potthoff functional (7) as a function of the relative phase $\phi$ of the
order parameters on the two layers. For $n=0.89$ and $n=0.91$, the minimum is
at $\phi=\pm\pi$ ($B_{2}$ representation). Near $n=0.93$, the profile changes
suddenly to one where the minimum is at $\phi=0$ ($B_{1}$ representation).
Note that the vertical scale is tiny ($10^{-5}$), in multiples of $t$, which
defines the energy unit here. This means that the energy difference between
the two representations $B_{1}$ and $B_{2}$ might just be too small to be of
consequence experimentally ($\sim 10^{-2}$meV or $\sim 10^{-1}$K in terms of
temperature), at an intermediate interlayer tunneling of $V_{1}=0.1$.

On the bottom half of Fig. 5, the same type of data is shown at strong
interlayer tunneling ($V_{1}=0.4$). There the transition between $B_{2}$ and
$B_{1}$ is gradual as the position of the minimum moves continuously from
$\phi=\pm\pi$ to $\phi=0$, with a spontaneous breaking of the $\phi\to-\phi$
symmetry. Even though this TRS breaking state is what we are looking for, such
a strong value of interlayer tunneling is unrealistic.

Figure 6: (Color online) Impurity models used in CDMFT. On the left: the
4-site cluster used for sites (2,3,4,5) and (7,8,9,10) of the unit cell, as
labeled in Fig. 1. On the right, the 2-site cluster used for sites (1,6).

Does this TRS breaking state have nontrivial topology? In a strongly
correlated system, this question may be answered through the properties of the
approximate interacting Green function (5) Wang and Zhang (2012a, b). The key
idea is to define a “topological Hamiltonian”
$h_{t}(\mathbf{k})=-G^{-1}(\mathbf{k},\omega=0)$, which can be diagonalized:

$h_{t}(\mathbf{k})|\alpha,\mathbf{k}\rangle=\mu_{\alpha}(\mathbf{k})|\alpha,\mathbf{k}\rangle$
(10)

One can then define a generalized Chern number just like in noninteracting
systems:

$C_{1}=\int\frac{d^{2}k}{2\pi}\;\mathcal{F}_{xy}(\mathbf{k})\quad\quad\mathcal{F}_{xy}(\mathbf{k})=\frac{\partial\mathcal{A}_{y}}{\partial
k_{x}}-\frac{\partial\mathcal{A}_{x}}{\partial k_{y}}$ (11)

with the Berry connection

$\mathcal{A}_{j}(\mathbf{k})=-i\sum_{\mu_{\alpha}(\mathbf{k})<0}\langle\alpha,\mathbf{k}|\partial_{k_{j}}|\alpha,\mathbf{k}\rangle,\qquad(j=x,y)$
(12)

When applying this formula to the TRS states found by VCA, we find the
topology to be trivial (the Chern number vanishes). This results from a
compensation between different regions of the Brillouin zone, with opposite
Berry curvature.

##  IV Results from Cluster Dynamical Mean Field Theory

In order to test the robustness of our predictions, we have also studied the
same system using cluster dynamical mean field theory (CDMFT) Lichtenstein and
Katsnelson (2000); Kotliar  _et al._ (2001); Liebsch  _et al._ (2008);
Sénéchal (2015) with an exact diagonalization solver at zero temperature (or
ED-CDMFT). Here the Weiss fields of VCA are replaced by a bath of uncorrelated
orbitals whose parameters are determined self-consistently. Because the
presence of this bath increases the size of the problem, the cluster cannot be
as large as in VCA and typically contains no more than 4 sites.

Each cluster, together with the associated bath, defines an Anderson impurity
model (AIM):

$H_{\rm
imp}=H_{c}+\sum_{\mu,\alpha}\theta_{\mu,\alpha}\left(c_{\mu}^{\dagger}a_{\alpha}^{\phantom{\dagger}}+\mbox{H.c.}\right)+\sum_{\alpha\beta}\epsilon_{\alpha\beta}a_{\alpha}^{\dagger}a_{\beta}^{\phantom{\dagger}}\leavevmode\nobreak\
\leavevmode\nobreak\ ,$ (13)

where $a_{\alpha}$ annihilates an electron in the bath orbital labeled
$\alpha$. The Nambu formalism must be used to incorporate pairing between bath
sites, within the matrix $\epsilon_{\alpha\beta}$, or within the hybridization
$\theta_{\mu,\alpha}$, depending on the impurity model. The index $\mu$ then
labels different sites of the cluster, together with the Nambu index, and
takes $2L$ values in a cluster with $L$ sites.

Figure 7: (Color online) Order parameter as a function of electron density
$n$, as found in CDMFT, for the two sets of inter-layer hopping. The anomalous
bath parameters were initially set to have opposite (top) or identical
(bottom) signs.

The bath parameters $\theta_{\mu,\alpha}$ and $\epsilon_{\alpha\beta}$ are
determined by an approximate self-consistent procedure, as proposed initially
in Caffarel and Krauth (1994), that goes as follows: (i) initial values of
these parameters are chosen on the first iteration. (ii) For each iteration,
the cluster Hamiltonian (13) is solved, i.e., the cluster Green function
$\mathbf{G}_{c}(\omega)$ is computed. The latter can be expressed as

$\mathbf{G}_{c}(\omega)^{-1}=\omega-\mathbf{t}_{c}-\bm{\Gamma}(\omega)-\bm{\Sigma}_{c}(\omega)$
(14)

where $\bm{\Gamma}(\omega)$ is the bath hybridization matrix:

$\Gamma_{ij}(\omega)=\sum_{\alpha,\alpha^{\prime}}\theta_{i\alpha}\left(\frac{1}{\omega-\bm{\epsilon}}\right)_{\alpha\alpha^{\prime}}\theta_{j\alpha^{\prime}}^{*}$
(15)

(iii) The bath parameters are updated, by minimizing the distance function:

$d(\bm{\epsilon},\bm{\theta})=\sum_{i\omega_{n}}W(i\omega_{n})\left[\mathbf{G}_{c}(i\omega_{n})^{-1}-\bar{\mathbf{G}}(i\omega_{n})^{-1}\right]$
(16)

where $\bar{\mathbf{G}}(\omega)$, the projected Green function, is defined as

$\bar{\mathbf{G}}(\omega)=\frac{1}{N}\sum_{\mathbf{k}}\mathbf{G}(\mathbf{k},\omega)\quad,\quad\mathbf{G}(\mathbf{k},\omega)=\frac{1}{\omega-\mathbf{t}_{\mathbf{k}}-\bm{\Sigma}_{c}(\omega)}\leavevmode\nobreak\
\leavevmode\nobreak\ .$ (17)

Ideally, $\bar{\mathbf{G}}(\omega)$ should coincide with the impurity Green
function $\mathbf{G}_{c}(\omega)$, but the finite number of bath parameters
does not allow for this correspondence at all frequencies, and so a distance
function $d(\epsilon_{r},\theta_{ir})$ is defined, with emphasis on low
frequencies along the imaginary axis. The weight function $W(i\omega_{n})$ is
where the method has some arbitrariness; in this work $W(i\omega_{n})$ is
taken to be a constant for all Matsubara frequencies lower than a cutoff
$\omega_{c}=2t$, with a fictitious temperature $\beta^{-1}=t/50$. (iv) We go
back to step (ii) and iterate until the bath parameters or the bath
hybridization function $\bm{\Gamma}(\omega)$ stop varying within some preset
tolerance.

In the current problem, the 10-site unit cell was separated in three impurity
problems: a four-site cluster on each layer (which together are equivalent to
the 8-site cluster used in VCA in the last section), made respectively of the
orbitals (2,3,4,5) and (7,8,9,10) as labeled on Fig. 1, and a two-site cluster
made of orbitals (1,6). These clusters are illustrated on Fig. 6. The 4-site
cluster is connected to 8 uncorrelated bath orbitals, and contains 6
independent parameters: Two bath energies $\epsilon_{1,2}$, two hybridization
$\theta_{1,2}$ and two pairing amplitudes $\Delta_{1,2}$ between bath
orbitals, with signs appropriate for describing $d$-wave superconductivity.
This way of parametrizing the bath is not the most general possible, but has
been successfully used in the past Kancharla  _et al._ (2008); Kyung  _et al._
(2009); Foley  _et al._ (2019). The two-site cluster connects the two layers
and also contains 6 bath parameters, except that the anomalous part is
contained in the hybridization, i.e., it connects the bath sites to the
cluster sites, not the bath sites themselves. In order to allow for a relative
phase between the pairing on the two layers, the pairing bath parameters
$\Delta_{1,2}$ on the square cluster of the second layer are allowed to take
complex values, whereas those on the first layer are assumed to be real. Once
a converged CDMFT solution is found, the same order parameters
$\Delta^{(\ell)}$ as in the previous section are computed.

Fig. 7 shows the results of CDMFT applied to this system, for both
intermediate (left) and strong (right) interlayer tunneling. The results
depend on the initial set of bath parameters. On the top panels, the bath
pairing parameters were initialized with opposite values on the two layers,
whereas on the bottom panels, they were initialized with the same values. At
intermediate interlayer tunneling ($V_{1}=0.1$), the order parameters stay
opposite throughout the doping range if the bath pairings are initialized this
way; in other words, if the system is primed in the $B_{2}$ representation, it
will stay in that representation. At strong interlayer tunneling
($V_{1}=0.4$), this only occurs if doping is large enough. In other words, for
doping 12% or less, the system primed in the $B_{2}$ representation will
either not converge, or converge to a normal solution, indicating its
incompatibility with the $B_{2}$ initial conditions. On the other hand, if the
system is primed in the $B_{1}$ representation, then it stays in the $B_{1}$
representation, except that, at strong interlayer tunneling, it converges for
larger values of doping, and converges to a normal solution at very small
doping.

It is thus difficult to discriminate between the $B_{1}$ and $B_{2}$
representations within CDMFT, which does not have the fine energy resolution
that VCA has. Nevertheless, we sense from the above results that the $B_{1}$
representation is preferred at low doping and the $B_{2}$ representation at
higher doping, but a strong interlayer tunneling is needed for that. Also,
despite allowing in principle for an arbitrary complex phase between the
anomalous bath parameters of the two layers, only the phases 0 and $\pi$ are
found: no state with spontaneous breaking of time reversal is found in CDMFT.
This may be related to the fact that the main 4-site impurity model in CDMFT
is confined to each layer, i.e., the complex, twisted inter-layer structure
has an impact only through the self-consistency relation. In studying such
systems, it seems that the VCA is a better choice.

##  V Conclusion

In a one-band Hubbard model for a cuprate bilayer twisted by an angle of
53.13∘, the relative phase of the superconducting order parameter in the two
layers depends on hole doping away from half-filling. In the underdoped
regime, the relative phase vanishes, whereas it is $\pi$ in the overdoped
regime. If the interlayer tunneling is strong, then there is an intermediate
phase between those two in which this phase varies continuously from 0 to
$\pi$. Time reversal symmetry is broken in that intermediate phase, but the
topology is trivial, at least as computed from the electron Green function. At
intermediate interlayer tunneling, this TRS breaking phase does not exist.

It is possible that this TRS breaking phase survives at weaker interlayer
tunneling if the twist angle is closer to 45∘. A twist angle of
$43.60^{\circ}$ corresponds to a unit cell of 58 copper sites Can  _et al._
(2021) and might be amenable to a similar VCA study, albeit markedly more
complex numerically. Work in this direction will be necessary in order to
assess whether this putative phase is realistic in strong-coupling
superconductivity.

###### Acknowledgements.

Computing resources were provided by Compute Canada and Calcul Québec. X.L. is
supported by the National Natural Science Foundation of China (Grant No.
11974293) and the Fundamental Research Funds for Central Universities (Grant
No. 20720180015). D.S. acknowledges NSERC (Canada) under grant
RGPIN-2020-05060.

## References

  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018b).
  * Andrei and MacDonald (2020) E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020).
  * Carr  _et al._ (2017) S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin, and E. Kaxiras, Phys. Rev. B 95, 075420 (2017).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), https://www.pnas.org/content/108/30/12233.full.pdf .
  * Lopes dos Santos  _et al._ (2007) J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Phys. Rev. Lett. 99, 256802 (2007).
  * Tarnopolsky  _et al._ (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Phys. Rev. Lett. 122, 106405 (2019).
  * Kennes  _et al._ (2021) D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. N. Basov, A. N. Pasupathy, and A. Rubio, Nature Physics 17, 155 (2021).
  * Regan  _et al._ (2020) E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl, and F. Wang, Nature 579, 359 (2020).
  * Tang  _et al._ (2020) Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan, and K. F. Mak, Nature 579, 353 (2020).
  * Wang  _et al._ (2020) L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai, B. Kim, K. Watanabe, T. Taniguchi, X. Zhu, J. Hone, A. Rubio, A. N. Pasupathy, and C. R. Dean, Nature Materials 19, 861 (2020).
  * Huang  _et al._ (2021) X. Huang, T. Wang, S. Miao, C. Wang, Z. Li, Z. Lian, T. Taniguchi, K. Watanabe, S. Okamoto, D. Xiao, S.-F. Shi, and Y.-T. Cui, Nature Physics 17, 715 (2021).
  * Liu  _et al._ (2020) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. Haei Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Nature 583, 221 (2020).
  * Cao  _et al._ (2020) Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 583, 215 (2020).
  * Zhang  _et al._ (2021) C. Zhang, T. Zhu, S. Kahn, S. Li, B. Yang, C. Herbig, X. Wu, H. Li, K. Watanabe, T. Taniguchi, S. Cabrini, A. Zettl, M. P. Zaletel, F. Wang, and M. F. Crommie, Nature Communications 12, 2516 (2021).
  * Chen  _et al._ (2019a) G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi, J. Jung, Y. Zhang, and F. Wang, Nature Physics 15, 237 (2019a).
  * Chen  _et al._ (2019b) G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature 572, 215 (2019b).
  * Chen  _et al._ (2020) G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, Z. Shi, T. Senthil, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Nature 579, 56 (2020).
  * Park  _et al._ (2021) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021).
  * Can  _et al._ (2021) O. Can, T. Tummuru, R. P. Day, I. Elfimov, A. Damascelli, and M. Franz, Nature Physics 17, 519 (2021).
  * Volkov  _et al._ (2020) P. A. Volkov, J. H. Wilson, and J. H. Pixley, Magic angles and current-induced topology in twisted nodal superconductors (2020), arXiv:2012.07860 [cond-mat.supr-con] .
  * Yu  _et al._ (2019) Y. Yu, L. Ma, P. Cai, R. Zhong, C. Ye, J. Shen, G. D. Gu, X. H. Chen, and Y. Zhang, Nature 575, 156 (2019).
  * Zhao  _et al._ (2019) S. Y. F. Zhao, N. Poccia, M. G. Panetta, C. Yu, J. W. Johnson, H. Yoo, R. Zhong, G. D. Gu, K. Watanabe, T. Taniguchi, S. V. Postolova, V. M. Vinokur, and P. Kim, Phys. Rev. Lett. 122, 247001 (2019).
  * Li  _et al._ (1999) Q. Li, Y. N. Tsay, M. Suenaga, R. A. Klemm, G. D. Gu, and N. Koshizuka, Phys. Rev. Lett. 83, 4160 (1999).
  * Takano  _et al._ (2002) Y. Takano, T. Hatano, A. Fukuyo, A. Ishii, M. Ohmori, S. Arisawa, K. Togano, and M. Tachiki, Phys. Rev. B 65, 140513 (2002).
  * Latyshev  _et al._ (2004) Y. I. Latyshev, A. P. Orlov, A. M. Nikitina, P. Monceau, and R. A. Klemm, Phys. Rev. B 70, 094517 (2004).
  * Klemm (2005) R. A. Klemm, Philosophical Magazine 85, 801 (2005), https://doi.org/10.1080/14786430412331314573 .
  * Yang  _et al._ (2018) Z. Yang, S. Qin, Q. Zhang, C. Fang, and J. Hu, Phys. Rev. B 98, 104515 (2018).
  * Zhu  _et al._ (2021) Y. Zhu, M. Liao, Q. Zhang, H.-Y. Xie, F. Meng, Y. Liu, Z. Bai, S. Ji, J. Zhang, K. Jiang, R. Zhong, J. Schneeloch, G. Gu, L. Gu, X. Ma, D. Zhang, and Q.-K. Xue, Phys. Rev. X 11, 031011 (2021).
  * Zhao  _et al._ (2021) S. Y. F. Zhao, N. Poccia, X. Cui, P. A. Volkov, H. Yoo, R. Engelke, Y. Ronen, R. Zhong, G. Gu, S. Plugge, T. Tummuru, M. Franz, J. H. Pixley, and P. Kim, Emergent interfacial superconductivity between twisted cuprate superconductors (2021), arXiv:2108.13455 [cond-mat.supr-con] .
  * Tummuru  _et al._ (2021) T. Tummuru, S. Plugge, and M. Franz, Josephson effects in twisted cuprate bilayers (2021), arXiv:2108.13457 [cond-mat.supr-con] .
  * Volkov  _et al._ (2021) P. A. Volkov, S. Y. F. Zhao, N. Poccia, X. Cui, P. Kim, and J. H. Pixley, Josephson effects in twisted nodal superconductors (2021), arXiv:2108.13456 [cond-mat.supr-con] .
  * Song  _et al._ (2021) X.-Y. Song, Y.-H. Zhang, and A. Vishwanath, Doping a moiré mott insulator: A t-j model study of twisted cuprates (2021), arXiv:2109.08142 [cond-mat.supr-con] .
  * Markiewicz  _et al._ (2005) R. S. Markiewicz, S. Sahrakorpi, M. Lindroos, H. Lin, and A. Bansil, Phys. Rev. B 72, 054519 (2005).
  * Potthoff (2003) M. Potthoff, European Physical Journal B 32, 429 (2003).
  * Potthoff  _et al._ (2003) M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev. Lett. 91, 206402 (2003).
  * Potthoff (2014) M. Potthoff, in  _DMFT at 25: Infinite dimensions, Lecture Notes of the Autumn School on Correlated Electrons 2014_ , edited by E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein (Forschungszentrum J{ü}lich, 2014).
  * Sénéchal  _et al._ (2000) Sénéchal, D. Perez, and M. Pioro-Ladrière, Physical Review Letters 84, 522 (2000).
  * Gros and Valenti (1993) C. Gros and R. Valenti, Phys. Rev. B 48, 418 (1993).
  * Wang and Zhang (2012a) Z. Wang and S.-C. Zhang, Phys. Rev. X 2, 031008 (2012a).
  * Wang and Zhang (2012b) Z. Wang and S.-C. Zhang, Phys. Rev. B 86, 165116 (2012b).
  * Lichtenstein and Katsnelson (2000) A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 62, R9283 (2000).
  * Kotliar  _et al._ (2001) G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Phys. Rev. Lett. 87, 186401 (2001).
  * Liebsch  _et al._ (2008) A. Liebsch, H. Ishida, and J. Merino, Phys. Rev. B 78, 165123 (2008).
  * Sénéchal (2015) D. Sénéchal, in  _Many-Body Physics: From Kondo to Hubbard_ , Vol. 5, edited by E. Pavarini, E. Koch, and P. Coleman (Forschungszentrum Jülich, 2015) pp. 13.1–13.22.
  * Caffarel and Krauth (1994) M. Caffarel and W. Krauth, Phys. Rev. Lett. 72, 1545 (1994).
  * Kancharla  _et al._ (2008) S. S. Kancharla, B. Kyung, D. Sénéchal, M. Civelli, M. Capone, G. Kotliar, and A.-M. S. Tremblay, Physical Review B 77, 184516 (2008).
  * Kyung  _et al._ (2009) B. Kyung, D. Sénéchal, and A.-M. S. Tremblay, Physical Review B 80, 205109 (2009).
  * Foley  _et al._ (2019) A. Foley, S. Verret, A.-M. S. Tremblay, and D. Sénéchal, Phys. Rev. B 99, 184510 (2019).