# Band Structure and Superconductivity in Twisted Trilayer Graphene

Võ Tiến Phong Department of Physics and Astronomy, University of Pennsylvania,
Philadelphia PA 19104  Pierre A. Pantaleón  Imdea Nanoscience, Faraday 9,
28015 Madrid, Spain  Tommaso Cea  Imdea Nanoscience, Faraday 9, 28015 Madrid,
Spain  Francisco Guinea  Imdea Nanoscience, Faraday 9, 28015 Madrid, Spain
Donostia International Physics Center, Paseo Manuel de Lardizábal 4, 20018 San
Sebastián, Spain

(August 28, 2024)

###### Abstract

We study the symmetries of twisted trilayer graphene’s band structure under
various extrinsic perturbations, and analyze the role of long-range electron-
electron interactions near the first magic angle. The electronic structure is
modified by these interactions in a similar way to twisted bilayer graphene.
We analyze electron pairing due to long-wavelength charge fluctuations, which
are coupled among themselves via the Coulomb interaction and additionally
mediated by longitudinal acoustic phonons. We find superconducting phases with
either spin singlet/valley triplet or spin triplet/valley singlet symmetry,
with critical temperatures of up to a few Kelvin for realistic choices of
parameters.

Introduction: Recently, it has been shown that superconductivity and exotic
correlated phases can emerge when two or three monolayers of graphene are
laterally stacked with a small relative twist angle between successive layers
Cao  _et al._ (2018a, b); Park  _et al._ (2021); Hao  _et al._ (2021).
Observation of insulating behavior and/or superconductivity have also been
reported for graphene bilayers on hBN substrates Moriyama  _et al._ (2019);
Shen  _et al._ (2020), ABC-stacked trilayers on hBN substrates Chen  _et al._
(2019a, b, c); Chittari  _et al._ (2019), pairs of graphene bilayers twisted
with respect to each other He  _et al._ (2020); Liu  _et al._ (2020); Cao  _et
al._ (2020); Tsai  _et al._ (2019), rhombohedral tetralayers Kerelsky  _et
al._ (2019), rhombohedral trilayers (Zhou  _et al._ , 2021a, b) and twisted
transition-metal-dichalcogenide layers Wang  _et al._ (2020). However, twisted
bilayer graphene (TBG) Cao  _et al._ (2018a, b) and its recent cousin, the
alternating-angle twisted trilayer graphene (TTG) (Park  _et al._ , 2021; Hao
_et al._ , 2021), are the only systems to date where superconductivity has
been unambiguously established Lu  _et al._ (2019); Arora  _et al._ (2020);
Yankowitz  _et al._ (2019).

Twisted trilayer graphene with alternating angles shows a series of magic
angles with remarkably flat bands, where, under certain conditions, the
trilayer system behaves as a TBG with renormalized interlayer hoppings plus a
monolayer Khalaf  _et al._ (2019). This follows from the mirror symmetry
around the central layer of the system. This symmetry can be broken by various
experimentally-relevant perturbations. In this case, calculations of the
electronic structure in the absence of this mirror symmetry shows hybridized
monolayer and twisted bilayer features Park  _et al._ (2021); Wu  _et al._
(2021); Ramires and Lado (2021); Călugăru  _et al._ (2021); Shin  _et al._
(2021a); Fischer  _et al._ (2021). The effect of electron-electron
interactions have also been considered Ramires and Lado (2021); Călugăru  _et
al._ (2021); Shin  _et al._ (2021a); Fischer  _et al._ (2021); Lake and
Senthil (2021); Goodwin  _et al._ (2021); Christos  _et al._ (2021), as well
as the role of external magnetic fields Qin and MacDonald (2021).

In this work, we study the symmetries of the continuum model, the effects of
symmetry-breaking perturbations on the band structure, and the renormalization
of electronic energies by electron-electron interactions within a self-
consistent Hartree formalism, complemented by the analysis of the exchange
potential at half filling. Then by using an extended Kohn-Luttinger mechanism
that includes electron-hole pairs, plasmons, and phonons, we analyze how the
screened Coulomb interaction can induce pairing in TTG.

Figure 1: Single-particle band structure with various perturbations: (a)
Intrinsic TTG with a twist angle of $\theta=1.59^{\circ}$ where the monolayer
band is decoupled from the narrow bilayer bands. (b) A perpendicular
displacement field $\Delta V=50$ meV between the outermost layers hybridizes
these bands. (c) A non-zero staggered sublattice potential,
$\delta_{1}=\delta_{3}=10$ meV, generically breaks both $\mathcal{T}C_{2z}$
and $M_{z}$ and gaps out the Dirac cones near charge neutrality, endowing the
resulting bands with possible nonzero Chern numbers. (d) For the special case,
$\delta_{1}=-\delta_{3},$ while both $\mathcal{T}C_{2z}$ and $M_{z}$ are
broken individually, $\mathcal{T}C_{2z}M_{z}$ is preserved. In this case, we
find that the bands are mixed, but the Dirac cones are preserved, though they
can be pushed to different energies.

Band structure and symmetries of the continuum model: To form TTG, we consider
three monolayers stacked one on top of the other in perfect atomic registry,
i.e. in the $AA$ stacking configuration. Labeling the layers consecutively, we
twist layers $\ell=1$ and $\ell=3$ by $-\theta/2$ and layer $\ell=2$ by
$+\theta/2$ about a fixed hexagon center. To describe the low-energy physics
at small angles, we adopt a valley-projected continuum Hamiltonian Bistritzer
and MacDonald (2011); Koshino  _et al._ (2018); Khalaf  _et al._ (2019); Li
_et al._ (2019)

$\begin{split}\mathcal{H}^{\theta}_{\nu}(\mathbf{k},\mathbf{r})&=\begin{pmatrix}\mathfrak{h}_{\nu}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)+V_{1}&\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&\mathfrak{t}_{\text{AA}}\\\
\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})&\mathfrak{h}_{\nu}\left(\mathbf{k}^{-\frac{\theta}{2}}\right)&\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\\\
\mathfrak{t}_{\text{AA}}^{\dagger}&\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&\mathfrak{h}_{\nu}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)+V_{3}\end{pmatrix},\\\
\mathfrak{h}_{\nu}\left(\mathbf{k}\right)&=-\hbar
v_{F}\left(\mathbf{k}-\mathbf{K}_{\nu}\right)\left(\nu\sigma_{x},\sigma_{y}\right),\\\
\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})&=\begin{pmatrix}w_{0}&w_{1}\\\
w_{1}&w_{0}\end{pmatrix}+\begin{pmatrix}w_{0}&w_{1}e^{-i\phi}\\\
w_{1}e^{i\phi}&w_{0}\end{pmatrix}e^{i\nu\mathbf{G}_{1}^{\text{M}}\cdot\mathbf{r}}\\\
&+\begin{pmatrix}w_{0}&w_{1}e^{i\phi}\\\
w_{1}e^{-i\phi}&w_{0}\end{pmatrix}e^{i\nu\left(\mathbf{G}_{1}^{\text{M}}+\mathbf{G}_{2}^{\text{M}}\right)\cdot\mathbf{r}},\\\
\mathfrak{t}_{\text{AA}}&=\begin{pmatrix}\gamma_{2}&0\\\
0&\gamma_{2}\end{pmatrix}\end{split}$ (1)

where $\phi={2\pi\nu/3}$, $\hbar v_{F}/a=2.135$ eV, $w_{0}=79.7$ meV and
$w_{1}=97.5$ meV are the interlayer tunneling amplitudes whose imbalance
accounts for some lattice relaxation Koshino  _et al._ (2018), $\nu=\pm$ is
the valley index, $\mathbf{k}^{\theta}=R(\theta)\mathbf{k}$,
$\mathbf{K}_{\pm}=\frac{4\pi}{3a}\left(\mp 1,0\right)$ are the microscopic
unrotated zone corners, and $\gamma_{2}$ is a small $AA$ coupling between the
first and third layers. If we use the conventional Slater-Koster
parametrization Koshino  _et al._ (2018), $\gamma_{2}\approx 0.3$ meV, which
justifies its neglect hereafter. It is included here to clarify some symmetry
properties to be discussed. We note in passing that this approximation
warrants further inspection in future work since it plays an important role in
the modeling of graphite’s band structure McClure (1957); Slonczewski and
Weiss (1958). To study various perturbations (Cea  _et al._ , 2020; Shi  _et
al._ , 2021; Shin  _et al._ , 2021b; Mao and Senthil, 2021; Lin  _et al._ ,
2021), we also add on-site couplings to the first and third layers. For
substrate induced masses,
$V_{1}=\delta_{1}\sigma_{z},V_{3}=\delta_{3}\sigma_{z},$ while for a
perpendicular displacement field, $V_{1},V_{3}\propto\mathbb{1}.$ Direct
coupling to $\ell=2$ is ignored since this layer is encapsulated and therefore
difficult to access.

We now consider the symmetries of Hamiltonian (1) with $V_{1}=V_{3}=0$. Our
choice of twist center ensures that the structure is maximally symmetric. A
different twist center, for instance, about a registered site, retains only a
subgroup of the maximal group. The vertical rotations about the $z$ axis,
$C_{6z},$ are respected by each plane individually. The horizontal rotations
about the in-plane axes, $C_{2x}$ and $C_{2y},$ preserved by TBG, are
explicitly broken due to the presence of the middle layer in TTG. However,
reflection symmetry about the plane of the middle layer $M_{z}$ that exchange
the top and bottom layers is preserved in TTG. Alternatively, one can
substitute $M_{z}$ for inversion symmetry $\mathcal{I}$ that is equivalent to
$M_{z}$ followed by a rotation. Therefore, the point symmetry group of the
lattice is $C_{6h}\otimes\mathcal{T},$ where $\mathcal{T}$ is local time-
reversal symmetry. In addition to these microscopic point symmetries, we can
also now impose valley projection symmetry which becomes exact in the limit of
large moiré wavelength. In this case, only a subset of
$C_{6h}\otimes\mathcal{T}$ is preserved at a single valley, namely $C_{2z}$ is
broken because it takes $\mathbf{K}_{+}$ to $\mathbf{K}_{-}$ in momentum
space. The symmetry group of the Hamiltonian thus factorizes into a direct
product of valley symmetry and a magnetic point symmetry group that consists
of $\\{C_{3z},\mathcal{T}C_{2z},M_{z}\\}.$

The presence of $z$-mirror symmetry allows us to simplify Eq. (1)
significantly. Mirror symmetry leaves the middle layer, $\ell=2,$ invariant,
while it interchanges layers $\ell=1$ and $\ell=3.$ So we can form linear
combinations of the wavefunctions of these two exchanged layers that are odd
and even under mirror symmetry. Such a transformation is executed by the
unitary operator $\mathcal{U}$ Li  _et al._ (2019); Khalaf  _et al._ (2019)

$\begin{pmatrix}\ket{\mathcal{M}_{\text{o}}^{1}}\\\
\ket{\mathcal{M}_{\text{o}}^{2}}\\\ \ket{\mathcal{M}_{\text{e}}}\\\
\end{pmatrix}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&0&1\\\ 0&\sqrt{2}&0\\\
-1&0&1\\\ \end{pmatrix}\begin{pmatrix}\ket{\ell=1}\\\ \ket{\ell=2}\\\
\ket{\ell=3}\\\ \end{pmatrix}.$ (2)

Under this transformation, the Hamiltonian factorizes into

$\begin{split}&\tilde{\mathcal{H}}_{\nu}^{\theta}(\mathbf{k},\mathbf{r})=\mathcal{U}\mathcal{H}_{\nu}^{\theta}(\mathbf{k},\mathbf{r})\mathcal{U}^{\dagger}=\tilde{\mathcal{H}}_{\text{TBG}}^{\theta,\nu}\left(\mathbf{k},\mathbf{r}\right)\oplus\tilde{\mathcal{H}}_{\text{MLG}}^{\theta,\nu}\left(\mathbf{k}\right)\\\
&=\begin{pmatrix}\mathfrak{h}_{\nu}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)+\mathfrak{t}_{\text{AA}}&\sqrt{2}\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&0\\\
\sqrt{2}\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})&\mathfrak{h}_{\nu}\left(\mathbf{k}^{-\frac{\theta}{2}}\right)&0\\\
0&0&\mathfrak{h}_{\nu}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)-\mathfrak{t}_{\text{AA}}\end{pmatrix},\end{split}$
(3)

where the mirror-odd sector,
$\tilde{\mathcal{H}}_{\text{TBG}}^{\theta,\nu}\left(\mathbf{k},\mathbf{r}\right),$
behaves as though it were a TBG system with enhanced interlayer couplings by a
factor of $\sqrt{2}$ and the mirror-even sector,
$\tilde{\mathcal{H}}_{\text{MLG}}^{\theta,\nu}\left(\mathbf{k},\mathbf{r}\right),$
mimics a single rotated layer of graphene. This transformation is convenient
because it casts the Hamiltonian into a direct sum of two Hamiltonians that
are well-understood. Denoting operators in the new basis with an overline
tilde, we see that $\tilde{C}_{3z}=C_{3z}$ and
$\tilde{\mathcal{T}}\tilde{C}_{2z}=\mathcal{T}C_{2z}$ since these preserve
layer index. In the mirror-odd sector, there are additional symmetries thus
far neglected that exchanges basis states $\ket{\mathcal{M}_{\text{o}}^{1}}$
and $\ket{\mathcal{M}_{\text{o}}^{2}}$ when $\gamma_{2}=0.$ We can call them
$\tilde{C}_{2x}$ and $\mathcal{T}\tilde{C}_{2y},$ although they do not
correspond naturally to rotations in the original layer basis. When
$\gamma_{2}\neq 0,$ these symmetries are broken, but all the other symmetries
of the lattice in the layer basis are still preserved.

Figure 2: Hartree-Fock renormalization of band structure: Low energy band
structure and DOS obtained for $\theta=1.59^{\circ}$ with (a) $\Delta V=0$ meV
and (b) $\Delta V=60$ meV. Different colors denote different fillings, as
indicated in the DOS. The horizontal lines denote the Fermi level. The Fock
bands at half filling are shown in for (c) $\Delta V=0$ and (d) $\Delta V=50$
meV. Dotted lines are the bands in the opposite valley.

Within a single valley, we have one unbounded Dirac cone coming from the
monolayer spectrum located at either $K$ or $K^{\prime}$; in addition, we also
have two Dirac cones at $K$ and $K^{\prime}$ from the TBG’s flat bands near
charge neutrality. These Dirac cones are protected by $\mathcal{T}C_{2z}$ as
in TBG Zou  _et al._ (2018); Song  _et al._ (2019). The two Dirac cones from
the TBG spectrum are related to each other by $\tilde{C}_{2x},$ and
consequently, are pinned to the same energy. However, in the absence of
particle-hole symmetry, there is no symmetry which pins the monolayer Dirac
cone to the same energy. Because of this energy offset, unlike in TBG, TTG in
the single-particle limit is formally metallic for all energies in the
continuum model. This important difference alters the Fermi surface geometry
of TTG when compared to TBG, and may be important in the analysis of
correlated insulating phases. For instance, there is technically no magic
angle around which isolated energy bands become spectrally flatten like in
TBG. However, the presence of $z$-mirror symmetry allows us to still define a
quasi-flat band by first projecting into the mirror-odd sector and then
imposing band-flatting conditions on just the two bands nearest to charge
neutrality within this sector. The highly-dispersive monolayer band becomes a
spectator in this scheme. For our choice of parameters, this scheme produces a
magic angle at $\theta\approx 1.57^{\circ}$. This magic angle occurs at a
larger angle than in TBG because of the enhancement in interlayer coupling
strengths. In addition, we can also add an interlayer bias that couples
directly to the $\ell=1$ and $\ell=3$ layers. The presence of this potential
preserves all in-plane symmetries, but it breaks $M_{z}.$ Because of this, the
monolayer and bilayer bands are hybridized. However, in this case, the number
of Dirac cones (with the same chirality) remains preserved because the in-
plane symmetries which protect them are maintained. To liberate the Dirac
cones, we can add staggered sublattice potentials to the outer layers.
Generically, this breaks both $M_{z}$ and $C_{2z}$ symmetries, and therefore,
can gap out the Dirac cones as well as hybridizing the mirror-odd and mirror-
even sectors.

Band structures for different sets of parameters are shown in Fig. 1. In Fig.
1a, we observe that pristine TTG features decoupled bands from the mirror-odd
and mirror-even sectors. A perpendicular displacement field $\Delta
V=V_{1}-V_{3}$ couples these two sectors, as shown in Fig. 1b. The presence of
a substrate induces an staggered sublattice potential that can couple the two
mirror sectors as well as gapping out all the Dirac cones in the system, as
shown in Fig. 1c. In this situation, the narrow bands are isolated from all
the other bands and can carry a finite Berry curvature. Depending on the value
of the staggered sublattice potential, the Chern topology of these bands can
be tuned, as detailed in Ref. (SI, ).

Long-range Coulomb interaction: At fractional fillings, interaction-driven
band renormalization can significantly distort the bands represented in Fig.
1. We now consider the effect of the long range Coulomb interaction treated
within the self-consistent Hartree-Fock approximation in Eq (1). As emphasized
in Refs. (Guinea and Walet, 2018; Cea  _et al._ , 2019; Cea and Guinea, 2020)
the inhomogeneous charge distribution leads to an electrostatic potential of
the order of $V_{0}\sim e^{2}/\epsilon L$, where $e$ is the electron charge
and $L$ the moiré length. In TBG near the magic angle, it has been shown that
the electrostatic Hartree interaction strongly distorts the band structure
(Cea  _et al._ , 2019). In TTG, the situation is quite similar; we find that
the self-consistent order parameter of the Hartree Hamiltonian also depends
linearly on the filling fraction owing to the dominance of the narrow bands
over the monolayer Dirac cones. Figure 2 displays the band structure of TTG as
a function of filling for a screening of $\epsilon=10$. The effect of the
Hartree potential varies significantly as a function of filling. This leads to
an approximate pinning of the van Hove singularities at the Fermi energy at
certain fillings, even in the presence of a displacement field.

We also consider the Fock potential. In TBG, this term leads to gaps and
broken symmetry phases Cea and Guinea (2020). We analyze the system at half
filling, where the Hartree potential is expected to vanish. The calculations
have been carried out by projecting the Fock potential on six central bands HF
. Results are shown in Fig. 2. The non-interacting flat bands are shifted
downwards, and significantly widened and separated. The monolayer band is less
affected, and the system remains metallic.

Superconductivity: Next, we study the onset of the superconducting (SC) order
in TTG, as induced by the combined effects of the 2D plasmon and the
longitudinal acoustic phonons of the three constituting layers. The
calculation represents a straightforward generalization of the one developed
by two of us in the recent paper (Cea and Guinea, 2021) for the case of the
TBG. It consists in solving the linearized gap equation in which the pairing
interaction is mediated by the Coulomb potential screened by the particle-hole
excitations and the acoustic phonons within the random phase approximation
(RPA). The idea strictly recalls the Kohn-Luttinger mechanism of the
superconductivity induced by repulsive interactions (Kohn and Luttinger,
1965).

The linearized gap equation can be written as:

$\displaystyle\tilde{\Delta}^{m_{1}m_{2}}_{\alpha\beta}(\bm{k})=\sum_{n_{1}n_{2}}\sum_{\bm{q}}\Gamma^{m_{1}m_{2}}_{n_{1}n_{2};\alpha\beta}(\bm{k},\bm{q})\tilde{\Delta}^{n_{1}n_{2}}_{\alpha\beta}(\bm{q}),$
(4)

where $\tilde{\Delta}$ is the order parameter, $m_{1,2},n_{1,2}$ are band
indices, $\alpha,\beta$ are the spin/valley flavors, and $\Gamma$ is the
pairing interaction kernel, which encodes the screened potential as well the
temperature dependence of the Eq. (4). A detailed description of the method
including the derivation of the Eq. (4) is provided in Ref. SI . In what
follows, we consider inter-valley superconductivity, meaning that $\alpha$ and
$\beta$ have opposite valley indices in the Eq. (4). In contrast, we do not
impose any constraint to the spin texture of the SC ground state, meaning that
our model can be compatible with both spin singlet or spin triplet
superconductivity. However, it is worth noting that the violation of the Pauli
limit reported by the recent experiment (Cao  _et al._ , 2021) strongly
suggests that the SC ground state of the TTG is a spin triplet.

Given the temperature dependence of the kernel $\Gamma$, Eq. (4) defines the
SC critical temperature, $T_{c}$, as the one at which the largest eigenvalue
of $\Gamma$ is equal to $1$. Fig. 3 shows the critical temperature as a
function of the filling, $\nu$, obtained for the twist angle
$\theta=1.59^{\circ}$ and inter-layer potential: $\Delta V=0$eV (red line) and
$\Delta V=0.06$eV (blue line). The corresponding band structure and density of
states (DOS) of the central bands are shown in the Fig. 2. Here, the DOS is
expressed in units of meV${}^{-1}A_{C}^{-1}$, $A_{C}$ being the area of the
moiré unit cell. As is evident, the behavior of $T_{c}$ upon varying the
filling follows the evolution of the DOS at the Fermi level, the latter being
marked by the horizontal lines.

Figure 3: Superconducting order parameter and transition temperature.
Amplitude (a) and phase (b) of the SC order parameter in the valence $v$ and
conduction $c$ band, as obtained for: $\theta=1.59^{\circ}$, $\Delta V=60$ meV
and $\nu=-1$. The black lines identify the Fermi surface. c) Critical
temperature as a function of the filling, $\nu$, obtained for the twist angle
$\theta=1.59^{\circ}$ and inter-layer potential: $\Delta V=0$ meV (red line)
and $\Delta V=60$ meV (blue line).

Finally, we study the symmetry of the SC order parameter, as given by the
eigenvector of the kernel $\Gamma$, Eq. (4), corresponding to the eigenvalue
$1$, at $T=T_{c}$. As a representative case, Fig. 3 shows the amplitude (a)
and the phase (b) of the order parameter, in the valence $v$ and conduction
$c$ band, as obtained for: $\theta=1.59^{\circ}$, $\Delta V=0.06$eV and
$n=-1$. The black lines identify the Fermi surface. Remarkably, also the
conduction band contributes to the superconductivity, albeit marginally. The
reason for that can be understood by noting that the conduction band crosses
the Fermi surface close to the Dirac points, see Fig. 2(b), where the bands
flatten. The results of Fig. 3 suggest an extended $s$-wave symmetry without
sign changes along the Fermi surface, which is analogous to what obtained in
the Ref. (Cea and Guinea, 2021) for the case of the TBG.

Conclusions: We have analyzed the symmetries, band structure, role of
interactions, and superconductivity of twisted trilayer graphene near the
first magic angle. Our results are similar to those found for twisted bilayer
graphene. The presence of an additional layer, and the accompanying electronic
bands with monolayer dispersion, does not modify qualitatively the features
found in twisted bilayer graphene.

To summarize, our main results are:

\- The long-range electron-electron interaction modifies significantly the
central bands. The Hartree potential widens the bands, leading to a bandwidth
$\sim 20-40$ meV. The exchange potential at half filling separates the
occupied valence band from the empty conduction band, which become $20-40$ meV
apart. The presence of an additional linearly dispersive band prevents the
formation of a gap and an insulating state.

\- Long-wavelength charge fluctuations, renormalized by their coupling by the
electrostatic interaction, and also by their interaction with longitudinal
acoustic phonons, are sufficient to induce superconductivity, with critical
temperatures of up to a few Kelvin. The order parameter has a significant
structure near the Fermi surface pockets. Pairing by long-wavelength
excitations, which do not mix valleys, leads to degenerate spin singlet/valley
triplet and spin triplet/valley singlet superconducting phases.

Acknowledgements. This work was supported by funding from the European
Commision, under the Graphene Flagship, Core 3, grant number 881603, and by
the grants NMAT2D (Comunidad de Madrid, Spain), SprQuMat and SEV-2016-0686,
(Ministerio de Ciencia e Innovación, Spain). VTP acknowledges support from the
NSF Graduate Research Fellowships Program and the P.D. Soros Fellowship for
New Americans.

## References

  * Cao  _et al._ (2018a) 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 (2018a), arXiv:1802.00553 .
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * Park  _et al._ (2021) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021).
  * Hao  _et al._ (2021) Z. Hao, A. M. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * Moriyama  _et al._ (2019) S. Moriyama, Y. Morita, K. Komatsu, K. Endo, T. Iwasaki, S. Nakaharai, Y. Noguchi, Y. Wakayama, E. Watanabe, D. Tsuya, K. Watanabe, and T. Taniguchi,  (2019), arXiv:1901.09356 .
  * Shen  _et al._ (2020) C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, T. Taniguchi, R. Yang, Z. Y. Meng, D. Shi, O. V. Yazyev, and G. Zhang, Nature Physics 16, 520 (2020).
  * Chen  _et al._ (2019a) 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 (2019a), arXiv:1901.04621 .
  * Chen  _et al._ (2019b) 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 (2019b).
  * Chen  _et al._ (2019c) G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Y. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, Z. Shi, T. Senthil, D. Goldhaber-Gordon, Y.-H. Y. Zhang, and F. Wang, Nature 579, 56 (2019c), arXiv:1905.06535 .
  * Chittari  _et al._ (2019) B. L. Chittari, G. Chen, Y. Zhang, F. Wang, and J. Jung, Physical Review Letters 122, 016401 (2019), arXiv:1806.00462 .
  * He  _et al._ (2020) M. He, Y. Li, J. Cai, Y. Liu, K. Watanabe, T. Taniguchi, X. Xu, and M. Yankowitz, Nature Physics (2020), 10.1038/s41567-020-1030-6, arXiv:2002.08904 .
  * 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), arXiv:1903.08130 .
  * 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), arXiv:1903.08596 .
  * Tsai  _et al._ (2019) K.-T. Tsai, X. Zhang, Z. Zhu, Y. Luo, S. Carr, M. Luskin, E. Kaxiras, and K. Wang,  (2019), arXiv:1912.03375 .
  * Kerelsky  _et al._ (2019) A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, D. Basov, A. Rubio, and A. N. Pasupathy,  (2019), arXiv:1911.00007 .
  * Zhou  _et al._ (2021a) H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, and A. F. Young, (2021a), arXiv:2104.00653 [cond-mat.mes-hall] .
  * Zhou  _et al._ (2021b) H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, (2021b), arXiv:2106.07640 [cond-mat.mes-hall] .
  * 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), arXiv:1910.12147 .
  * Lu  _et al._ (2019) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, and D. K. Efetov, Nature 574, 653 (2019), arXiv:1903.06513 .
  * Arora  _et al._ (2020) H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, K. Watanabe, T. Taniguchi, J. Alicea, S. Nadj-Perge, H. Singh Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Zaky Wilson, X. Xu, J.-H. Chu, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, arXiv e-prints , arXiv:2002.03003 (2020), arXiv:2002.03003 .
  * Yankowitz  _et al._ (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019).
  * Khalaf  _et al._ (2019) E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019).
  * Wu  _et al._ (2021) Z. Wu, Z. Zhan, and S. Yuan, Sci. China-Phys. Mech. Astron. 64 (2021).
  * Ramires and Lado (2021) A. Ramires and J. L. Lado, (2021), arXiv:2102.03312 [cond-mat.mes-hall] .
  * Călugăru  _et al._ (2021) D. Călugăru, F. Xie, Z.-D. Song, B. Lian, N. Regnault, and B. A. Bernevig, Phys. Rev. B 103, 195411 (2021).
  * Shin  _et al._ (2021a) J. Shin, B. L. Chittari, and J. Jung, (2021a), arXiv:2104.01570 [cond-mat.mes-hall] .
  * Fischer  _et al._ (2021) A. Fischer, Z. A. H. Goodwin, A. A. Mostofi, J. Lischner, D. M. Kennes, and L. Klebl, (2021), arXiv:2104.10176 [cond-mat.mes-hall] .
  * Lake and Senthil (2021) E. Lake and T. Senthil, (2021), arXiv:2104.13920 [cond-mat.mes-hall] .
  * Goodwin  _et al._ (2021) Z. A. H. Goodwin, L. Klebl, V. Vitale, X. Liang, V. Gogtay, D. M. K. Xavier van Gorp, A. A. Mostofi, and J. Lischner, (2021), arXiv:2105.12641 [cond-mat.mes-hall] .
  * Christos  _et al._ (2021) M. Christos, S. Sachdev, and M. S. Scheurer, (2021), arXiv:2106.02063 [cond-mat.mes-hall] .
  * Qin and MacDonald (2021) W. Qin and A. H. MacDonald, (2021), arXiv:2104.14026 [cond-mat.mes-hall] .
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * Koshino  _et al._ (2018) M. Koshino, N. F. Q. Yuan, T. Koretsune, M. Ochi, K. Kuroki, and L. Fu, Phys. Rev. X 8, 031087 (2018).
  * Li  _et al._ (2019) X. Li, F. Wu, and A. H. MacDonald, arXiv preprint arXiv:1907.12338 (2019).
  * McClure (1957) J. W. McClure, Phys. Rev. 108, 612 (1957).
  * Slonczewski and Weiss (1958) J. C. Slonczewski and P. R. Weiss, Phys. Rev. 109, 272 (1958).
  * Cea  _et al._ (2020) T. Cea, P. A. Pantaleón, and F. Guinea, Phys. Rev. B 102, 155136 (2020).
  * Shi  _et al._ (2021) J. Shi, J. Zhu, and A. H. MacDonald, Phys. Rev. B 103, 075122 (2021).
  * Shin  _et al._ (2021b) J. Shin, Y. Park, B. L. Chittari, J.-H. Sun, and J. Jung, Phys. Rev. B 103, 075423 (2021b).
  * Mao and Senthil (2021) D. Mao and T. Senthil, Phys. Rev. B 103, 115110 (2021).
  * Lin  _et al._ (2021) X. Lin, K. Su, and J. Ni, 2D Materials 8, 025025 (2021).
  * Zou  _et al._ (2018) L. Zou, H. C. Po, A. Vishwanath, and T. Senthil, Phys. Rev. B 98, 085435 (2018).
  * Song  _et al._ (2019) Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, and B. A. Bernevig, Phys. Rev. Lett. 123, 036401 (2019).
  * (44) See Supplementary material.
  * Guinea and Walet (2018) F. Guinea and N. R. Walet, Proceedings of the National Academy of Sciences of the United States of America 115, 13174 (2018), arXiv:1806.05990 .
  * Cea  _et al._ (2019) T. Cea, N. R. Walet, and F. Guinea, Physical Review B 100, 205113 (2019), arXiv:1906.10570 .
  * Cea and Guinea (2020) T. Cea and F. Guinea, Physical Review B 102, 045107 (2020), arXiv:2004.01577 .
  * (48) We assume that the initial hoppings include the renormalization due to the exchange potential from the bands other than the sisx central bands, so that no substraction of the exchange potential is required.
  * Cea and Guinea (2021) T. Cea and F. Guinea, arXiv e-prints , arXiv:2103.01815 (2021), arXiv:2103.01815 [cond-mat.str-el] .
  * Kohn and Luttinger (1965) W. Kohn and J. M. Luttinger, Phys. Rev. Lett. 15, 524 (1965).
  * Cao  _et al._ (2021) Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, arXiv e-prints , arXiv:2103.12083 (2021), arXiv:2103.12083 [cond-mat.mes-hall] .

Supplementary Material

##  .1 Construction of Hamiltonian

Here, we study twisted trilayer graphene in which all three layers are rotated
with respect to adjacent layers. Before proceeding to the analysis of
multilayer systems, we first establish the convention of monolayer graphene.
We choose to align the $p_{z}$ orbitals on both sublattices so that they have
the same sign convention. In this case, $-\gamma_{0},$ the hopping integral
between in-plane nearest neighbors, is a negative constant (we have chosen a
minus sign so that $\gamma_{0}$ is positive). Let $a$ denote the lattice
constant, $\mathbf{a}_{1}=a\left(1,0\right),$ and
$\mathbf{a}_{2}=a\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right).$ The primitive
reciprocal lattice vectors
$\mathbf{g}_{1}=\frac{2\pi}{a}\left(1,-\frac{1}{\sqrt{3}}\right)$ and
$\mathbf{g}_{2}=\frac{2\pi}{a}\left(0,\frac{2}{\sqrt{3}}\right).$ The zone
corners are located at $\mathbf{K}_{\pm}=\frac{4\pi}{3a}\left(\mp 1,0\right).$
The minimal tight-binding model in the basis of $\ket{A}$ and $\ket{B}$ is

$\mathfrak{h}(\mathbf{k})=\begin{pmatrix}0&-\gamma_{0}\phi(\mathbf{k})\\\
-\gamma_{0}\phi(\mathbf{k})^{*}&0\end{pmatrix},$ (S1)

where $\phi(\mathbf{k})=\sum_{j=1}^{3}e^{i\mathbf{k}\cdot\bm{\delta}_{j}}$ and
$\bm{\delta}_{j}=\frac{a}{\sqrt{3}}\left(\sin\left(\frac{2\pi
j}{3}\right),-\cos\left(\frac{2\pi j}{3}\right)\right).$ Expanding about the
zone corners, we find

$\begin{split}\mathfrak{h}\left(\mathbf{p}=\mathbf{k}-\mathbf{K}_{-}\right)&=\begin{pmatrix}0&\frac{\sqrt{3}}{2}a\gamma_{0}\left(p_{x}+ip_{y}\right)\\\
\frac{\sqrt{3}}{2}a\gamma_{0}\left(p_{x}-ip_{y}\right)&0\end{pmatrix}=-\frac{\sqrt{3}}{2}a\gamma_{0}\mathbf{p}\cdot\left(-\sigma_{x},\sigma_{y}\right),\\\
\mathfrak{h}\left(\mathbf{p}=\mathbf{k}-\mathbf{K}_{+}\right)&=\begin{pmatrix}0&-\frac{\sqrt{3}}{2}a\gamma_{0}\left(p_{x}-ip_{y}\right)\\\
-\frac{\sqrt{3}}{2}a\gamma_{0}\left(p_{x}+ip_{y}\right)&0\end{pmatrix}=-\frac{\sqrt{3}}{2}a\gamma_{0}\mathbf{p}\cdot\left(+\sigma_{x},\sigma_{y}\right).\end{split}$
(S2)

Defining $\hbar v_{F}=\sqrt{3}a\gamma_{0}/2,$ we can summarize as follows:
$\mathfrak{h}\left(\mathbf{k}\right)=-\hbar
v_{F}\left(\mathbf{k}-\mathbf{K}_{\nu}\right)\cdot\left(\nu\sigma_{x},\sigma_{y}\right).$

Next, we consider stacking up three layers with a relative twist between them.
We denote a rotation matrix of the form

$R(\theta)=\begin{pmatrix}\cos\left(\theta\right)&-\sin\left(\theta\right)\\\
\sin\left(\theta\right)&\cos\left(\theta\right)\end{pmatrix}.$ (S3)

This performs a counter-clockwise rotation. We consider three layers
originally-aligned in $AA$ configuration. We then rotate layer $\ell=1$ by
$-\theta/2,$ layer $\ell=2$ by $\theta/2,$ and layer $3$ by $-\theta/2.$ The
primitive lattice vectors are rotated by
$\mathbf{a}^{\ell}_{i}=R(\pm\theta/2)\mathbf{a}_{i}$ and the primitive
reciprocal lattice vectors by
$\mathbf{g}^{\ell}_{i}=R(\pm\theta/2)\mathbf{g}_{i}.$ The primitive moiré
reciprocal lattice vectors are then given by

$\mathbf{G}_{1}^{\text{M}}=\mathbf{g}^{1}_{1}-\mathbf{g}^{2}_{1}=G_{\text{M}}\left(-\frac{1}{2},-\frac{\sqrt{3}}{2}\right),\quad\mathbf{G}_{2}^{\text{M}}=\mathbf{g}^{1}_{2}-\mathbf{g}^{2}_{2}=G_{\text{M}}\left(-1,0\right).$
(S4)

where $G_{\text{M}}=8\pi\sin\left(\theta/2\right)/(\sqrt{3}a).$ The primitive
moiré lattice vectors

$\mathbf{L}_{1}^{\text{M}}=L_{\text{M}}\left(0,-1\right),\quad\mathbf{L}_{2}=L_{\text{M}}\left(-\frac{\sqrt{3}}{2},\frac{1}{2}\right),$
(S5)

$L_{\text{M}}=a/(2\sin\left(\theta/2\right)).$ The basic valley-projected
Hamiltonian in mixed real-space and reciprocal-space representation is

$\mathcal{H}^{\theta}_{\nu}(\mathbf{k})=\begin{pmatrix}\mathfrak{h}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)&\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&0\\\
\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})&\mathfrak{h}\left(\mathbf{k}^{-\frac{\theta}{2}}\right)&\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\\\
0&\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&\mathfrak{h}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)\end{pmatrix},$
(S6)

where $\mathbf{k}^{(\pm\theta/2)}=R\left(\pm\theta/2\right)\mathbf{k}$ (note
the convention that layer 1 is rotated by $-\theta/2$ so the wavevector must
be rotated by $+\theta/2$ in the opposite direction to compensate), and

$\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})=\begin{pmatrix}w_{0}&w_{1}\\\
w_{1}&w_{0}\end{pmatrix}+\begin{pmatrix}w_{0}&w_{1}e^{-2\pi\nu i/3}\\\
w_{1}e^{2\pi\nu
i/3}&w_{0}\end{pmatrix}e^{i\nu\mathbf{G}_{1}^{\text{M}}\cdot\mathbf{r}}+\begin{pmatrix}w_{0}&w_{1}e^{2\pi\nu
i/3}\\\ w_{1}e^{-2\pi\nu
i/3}&w_{0}\end{pmatrix}e^{i\nu\left(\mathbf{G}_{1}^{\text{M}}+\mathbf{G}_{2}^{\text{M}}\right)\cdot\mathbf{r}}.$
(S7)

The Hamiltonian in Eq. (S6) is in the “natural” basis of sublattices and
layers. As shown in Ref. Khalaf  _et al._ (2019), we can perform a unitary
transformation to simplify the Hamiltonian

$\begin{split}\tilde{\mathcal{H}}_{\nu}^{\theta}(\mathbf{k},\mathbf{r})=&\begin{pmatrix}\frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\\\
0&1&0\\\
-\frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\end{pmatrix}\begin{pmatrix}\mathfrak{h}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)&\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&0\\\
\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})&\mathfrak{h}\left(\mathbf{k}^{-\frac{\theta}{2}}\right)&\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\\\
0&\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&\mathfrak{h}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)\end{pmatrix}\begin{pmatrix}\frac{1}{\sqrt{2}}&0&-\frac{1}{\sqrt{2}}\\\
0&1&0\\\ \frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\end{pmatrix}\\\
&=\begin{pmatrix}\mathfrak{h}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)&\sqrt{2}\left[\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})\right]^{\dagger}&0\\\
\sqrt{2}\mathfrak{t}_{\nu}^{\theta}(\mathbf{r})&\mathfrak{h}\left(\mathbf{k}^{-\frac{\theta}{2}}\right)&0\\\
0&0&\mathfrak{h}\left(\mathbf{k}^{+\frac{\theta}{2}}\right)\end{pmatrix}.\end{split}$
(S8)

The new basis states mix layers as follows

$\begin{pmatrix}\ket{\mathcal{M}_{\text{o}}^{1}}\\\
\ket{\mathcal{M}_{\text{o}}^{2}}\\\
\ket{\mathcal{M}_{\text{e}}}\end{pmatrix}=\begin{pmatrix}\frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\\\
0&1&0\\\
-\frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{2}}\end{pmatrix}\begin{pmatrix}\ket{\ell=1}\\\
\ket{\ell=2}\\\
\ket{\ell=3}\end{pmatrix}=\begin{pmatrix}\frac{1}{\sqrt{2}}\ket{\ell=1}+\frac{1}{\sqrt{2}}\ket{\ell=3}\\\
\ket{\ell=2}\\\
\frac{1}{\sqrt{2}}\ket{\ell=3}-\frac{1}{\sqrt{2}}\ket{\ell=1}\end{pmatrix}.$
(S9)

Figure S1: Representation of twisted trilayer graphene in real and reciprocal
space. (a) Carbon atoms from layers 1, 2, and 3 are shaded red, blue, and
green respectively. We note that the red and green circles often coincide
because if the middle layer were not present, the top and bottom layers would
be right on top of one another. (b) Moiré Brillouin zone is shown in green
with high-symmetry points indicated. The microsopic graphene Brillouin zones
are shown in blue, for layer 1, and red, for layer 2. For clarity, we do not
show explicitly the Brillouin zone for layer 3 since it is identical to that
of layer 1. The high-symmetry points in the moiré Brillouin zone are denoted
with bars for emphasis here to distinguish them from the microscopic high-
symmetry points; these bars are dropped in the main text for simplicity.

##  .2 Symmetry-Breaking Perturbations

The appeal of TTG is its robust tunability. As such, we now consider some ways
through which such tunability can be controlled. The most generic spatially-
uniform perturbation we consider takes the form, in the two different bases,

$\mathcal{V}=\begin{pmatrix}V_{1}&0&0\\\ 0&0&0\\\
0&0&V_{3}\end{pmatrix}\quad\text{and}\quad\tilde{\mathcal{V}}=\frac{1}{2}\begin{pmatrix}V_{1}+V_{3}&0&V_{3}-V_{1}\\\
0&0&0\\\ V_{3}-V_{1}&0&V_{1}+V_{3}\end{pmatrix},$ (S10)

where $V_{1}$ and $V_{3}$ are $2\times 2$ matrices that act on the sublattice
degrees of freedom. We do not consider perturbations that directly couple to
the middle layer since it is encapsulated and is presumably difficult to
practically access. For $V_{1}=V_{3},$ $z$-mirror symmetry is preserved;
otherwise, the mirror-odd and mirror-even sectors are mixed by $V_{3}-V_{1}.$

There are two limits to consider: $V_{1}=+V_{3}$ and $V_{1}=-V_{3}$. The
former retains the decoupled bilayer-monolayer form because the off-diagonal
terms vanish. In this case, we can still consider the composite system as two
separate systems: one monolayer and one twisted bilayer graphene. However,
both the bilayer system and the monolayer system may have additional
perturbations that may lower their symmetries. In the latter situation, the
bilayer and monolayer systems are mixed via the off-diagonal terms. An example
of the second limit is the displacement field, $\pm\Delta V/2$, between the
outermost layers. The generic system varies between this two limits. In
addition to the extrinsic perturbations discussed above, there is also an
intrinsic perturbation to Eq. (S6) in the original layer basis of the form

$\begin{split}\mathcal{V}_{AA}=\begin{pmatrix}0&0&\mathfrak{t}_{\text{AA}}\\\
0&0&0\\\
\mathfrak{t}_{\text{AA}}^{\dagger}&0&0\end{pmatrix},\quad\mathfrak{t}_{\text{AA}}=\begin{pmatrix}\gamma_{2}&0\\\
0&\gamma_{2}\end{pmatrix},\end{split}$ (S11)

where $\gamma_{2}$ is the $AA$ coupling constant of the first and third
layers. Because these two layers are far from each other, $\gamma_{2}$ is
presumably a very small parameter compared to the other energy scales in the
problem. To estimate, we use the conventional Slater-Koster parameterization
Koshino  _et al._ (2018)

$\gamma_{2}(r)=V_{pp\sigma}^{0}\exp\left(\frac{d_{0}-r}{r_{0}}\right).$ (S12)

We take $V_{pp\sigma}^{0}=0.48$ eV, $d_{0}=0.335$ nm, and $r_{0}=0.045$ nm. We
take $r=2d_{0}$ to find $\gamma_{2}\approx 0.3$ meV, indeed a very small value
as expected. However, as noted in the main text, this may be an underestimate
of $\gamma_{2}$ since it assumes that this parameter is environment-
independent.

Figure S2: Symmetry-Breaking Perturbations: a) Plot of TTG with a twist angle
of $\theta=1.59^{\circ}$ and $\gamma_{2}=10$ meV. b) For
$V_{1}=V_{3}=2\Delta\tilde{V}=10$ meV the Dirac cone is shifted but is still
decoupled from the narrow bands. c) For $V_{1}=-V_{3}$, both Dirac cones are
coupled to the narrow bands. This situation is equivalent to a displacement
field, $\Delta V$, between the outermost layers. In the figure we use $2\Delta
V=10$ meV.

##  .3 Band Structure

We now study the band structure of TTG. In the intrinsic case, $z$-mirror
symmetry allows us to partition the Hamiltonian into a direct sum of two
Hamiltonians: one describing a twisted graphene bilayer with enhanced
interlayer coupling by a factor of $\sqrt{2}$ and one describing a simple
graphene monolayer. Therefore, much of the band structure analysis of this
case is the same as the analysis of TBG. The monolayer Dirac cone is simply a
spectator here. The system as a whole is metallic for all energies in this
continuum construction without particle-hole symmetry. However, once we
project onto just the mirror-odd sector, we can still speak of flat bands
within just this sector of the Hamiltonian. This is useful for later because
once we turn on symmetry-breaking perturbations that hybridize the mirror-odd
and mirror-even sectors, some degeneracies will be gapped out, and true
isolated flat bands can emerge. Importantly, the flatness of these bands is
due primarily to the flatness of the parent bands that descend from the
mirror-odd sector. Thus, focusing for a moment on just the mirror-odd sector,
$\ket{\mathcal{M}_{\text{o}}^{1}}$ and $\ket{\mathcal{M}_{\text{o}}^{2}},$ we
first establish regions of parameter space where flat bands emerge.

By a magic angle, we mean an angle at which a pair of entangled flat bands
emerges near charge neutrality; these flat bands are isolated from the other
energy bands. This definition does not produce unique magic angles because the
criteria for flatness remain ambiguous. One common definition is to pin magic
angles at those at which the Dirac velocity at $\bar{K}$ and
$\bar{K}^{\prime}$ vanishes in all directions. However, this definition does
not always coincide with small bandwidth. Indeed, vanishing of the Dirac
velocity does not even imply band isolation from the other high-energy bands.
For us, we instead use the local minimalization of the bandwidth as the
definition of magic angle. To find the magic angle, we find regions of
$\theta$ where isolated bands can occur. This can be found numerically by
computing the four lowest energy eigenvalues near charge neutrality at
$\bar{\Gamma}.$ When the top two or bottom two eigenvalues coincide, this
means that a gap cannot exist. Then, for regions where a gap could exist, we
compute the corresponding bandwidths, defined as the maxima of the energies
computed at high-symmetry points $\bar{\Gamma},$ $\bar{K},$ $\bar{M},$ and
$\bar{M}^{\prime}$ minus the minima. These may, of course, not be the true
bandwidths, but they are sufficient to identify the magic angle. Using this
condition we find a magic angle around $\theta\approx 1.57^{\circ}$, as shown
in Fig. S4. We note that due to the enhancement of the interlayer coupling,
the magic angle, within the window of parameters in Fig. S4b, occurs at a
larger value than that known for TBG. In addition, this magic angle is
persistent even in the presence of a displacement field, as shown in Fig. S4a.

##  .4 Sublattice potential and topological phases

We now consider the presence of an staggered sublattice potential. As in the
case of TBG, a mass gap can be induced by the presence of a substrate (Cea
_et al._ , 2020; Shi  _et al._ , 2021; Shin  _et al._ , 2021b; Mao and
Senthil, 2021; Lin  _et al._ , 2021). For TTG, we consider a mass gap of the
form $\delta_{l}\sigma_{z}$ with $\ell=\\{1,3\\}$ labelling the corresponding
term acting on the bottom and top layer. As shown in Fig. 1c, a staggered
sublattice potential $\delta_{1}=\delta_{3}=10$ meV opens a gap between the
narrow bands and the Dirac cones. Fig. S3 shows a density plot of the
magnitude of the gap between the narrow bands and the Dirac cones as a
function of the mass terms $\delta_{1/3}$. Fig. S3a is the gap between the
upper Dirac cone and the top narrow band, Fig. S3b is the gap between the
narrow bands and Fig. S3c the gap between the lower narrow band and the lower
Dirac cone. Some effects are in order here; in Fig. S3a and Fig. S3c, we
notice that if $\delta_{1}=-\delta_{3}$, each narrow band is connected to a
Dirac cone at $K$ and connected to each other at $K^{\prime}$. In other
regions of the same figure, the narrow bands are decoupled from the cones with
a finite gap. However, as shown in the dark regions in Fig. S3b if any
$\delta_{1}$ or $\delta_{3}$ is zero, the narrow bands are decoupled from the
Dirac cones but with no gap between them. In Fig. S3d we display the
topological phases as a function of the staggered sublattice potential. In the
region where both $\delta_{1}$ and $\delta_{3}$ are negatives both narrow
bands are isolated with Chern numbers
$\left(\mathcal{C}_{1},\mathcal{C}_{2}\right)=\left(-1,1\right)$, if both are
positives the phase is $\left(1,-1\right)$. The trivial phase
$\left(0,0\right)$ is when both mass terms have different sign.

Figure S3: Mass gap: At the point $\bm{K}$ and as a function of the mass terms
we calculate the gap between a) the upper Dirac cone and the top narrow band,
b) top and bottom narrow bands and c) bottom narrow band and lower Dirac cone.
In d) we show the topological phases as a function of the staggered sublattice
potential. Continuous lines are the boundaries in which there is a topological
phase transition. Dashed diagonal line separates the region within the same
topological phase where the signs of the sublattice potential are inverted.
Figure S4: Magic-angle region: The top plot shows the bandwidth, $W$, as
defined in the text, as a function of angle and displacement field. Bottom
plot shows the four energy eigenvalues closest to charge neutrality of the
mirror-odd bands at $\bar{\Gamma}.$ The band gaps between the middle two flat
bands and the higher-energy bands close when the top two or bottom two
eigenvalues coincide. Region in $\theta$ where magic angle occur is shaded in
light gray.

##  .5 Self-consistent Hartree Interaction

In momentum space, the self-consistent Hartree interaction has matrix elements
given by (Guinea and Walet, 2018; Cea  _et al._ , 2019)

$\rho_{H}(\bm{G})=4V_{C}\left(\bm{G}\right)\int_{BZ}\frac{d^{2}\bm{k}}{V_{mBZ}}\sum_{\bm{G}^{\prime},l}\phi_{k,l}^{{\dagger}}\left(\bm{G}^{\prime}\right)\phi_{k,l}\left(\bm{G}+\bm{G}^{\prime}\right),$
(S13)

where $V_{C}\left(\bm{G}\right)=\frac{2\pi
e^{2}}{\epsilon\left|\bm{G}\right|}$ is the Fourier transform of the Coulomb
potential evaluated at $\bm{G}$, $V_{mBZ}$ is the area of the moiré Brillouin
zone, and the factor 4 takes into account spin/valley degeneracy. The
parameter $l$ is a band index resulting from the diagonalization of the full
Hamiltonian in Eq. (1). In the above equation, $\phi_{k,l}\left(\bm{G}\right)$
is the amplitude for an electron to occupy a state with momentum
$\bm{k}+\bm{G}$. The value of $\rho_{H}(\bm{G})$ in Eq. (S13) depends on the
extent of the wavefunctions in momentum space, hence we can write the Fourier
expansion of the Hartree potential in real space as

$V_{\text{H}}(\bm{r})=V_{0}\sum_{n}\rho_{\text{H}}(\bm{G}_{n})e^{i\bm{G}_{n}\cdot\bm{r}}$
(S14)

with
$\rho_{\text{H}}(\bm{G}_{n})=\left|\rho_{\text{H}}(\bm{G}_{n})\right|e^{i\arg\left[\rho_{\text{H}}(\bm{G}_{n})\right]}$
a complex number with $-\pi<\arg\left[\rho_{\text{H}}(\bm{G}_{n})\right]<\pi$
and $V_{0}=e^{2}/\varepsilon L$ the effective Coulomb potential. As shown in
Ref. (Guinea and Walet, 2018), the dominant contribution to the Hartree
potential comes from the first star,
$\bm{G}_{n}=\pm\bm{G}_{1},\pm\bm{G}_{2},\pm(\bm{G}_{1}+\bm{G}_{2})$, of
reciprocal lattice vectors. In addition, by exploiting the $C_{6}$ symmetry,
the Fourier components of the Hartree potential are equally weighted. To solve
the self-consistent Hartree Hamiltonian, the charge distribution is
approximated as $\rho_{H}=\overline{\rho}_{H}+\delta\rho_{H}$ where
$\overline{\rho}_{H}$ is a constant which takes into account the total density
from all bands not included in the calculations (Guinea and Walet, 2018). The
charge distribution is fixed by considering an homogeneous state at the charge
neutrality point, this is $\rho_{H}=0$. Therefore, the integral in Eq. (S13)
is evaluated only over energy levels with $E_{l}\left(\bm{k}\right)$ between
the charge neutrality point and Fermi level. The matrix elements of the
Hartree potential $V_{H}$ in Eq. (S13) depends implicitly on the filling
fraction $\nu$ of the conduction band. For fully-filled valence and conduction
bands we have $\nu=4$ and when they are both empty we have $\nu=-4$. For a
given value of the filling, we calculate the miniband spectrum and
wavefunctions with self-consistent diagonalization of the main Hamiltonian by
considering a coupling up to 5 stars, this corresponds to a total of $N=91$
vectors in the reciprocal space.

##  .6 The Kohn-Luttinger formalism for TTG

Because of the lack of translational invariance, it is convenient to write the
linearized equation for the order parameter (OP), $\Delta$, in real space:

$\displaystyle\Delta^{i_{1}i_{2}}_{\alpha\beta}(\bm{r}_{1},\bm{r}_{2})=-{\cal
V}^{scr}(\bm{r}_{1},\bm{r}_{2})\int_{\Omega}\,d^{2}\bm{r}_{3}d^{2}\bm{r}_{4}\sum_{i_{3}i_{4}}K_{B}T\sum_{\omega}\mathcal{G}^{i_{1}i_{3}}_{\bm{r}_{1}\bm{r}_{3},\alpha}\left(i\hbar\omega\right)\mathcal{G}^{i_{2}i_{4}}_{\bm{r}_{2}\bm{r}_{4},\beta}\left(-i\hbar\omega\right)\Delta^{i_{3}i_{4}}_{\alpha\beta}(\bm{r}_{3},\bm{r}_{4}),$
(S15)

where the $i$’s denote the layer and sublattice, $\alpha\neq\beta$ is the
index of flavor, which encodes both the valley and spin degrees of freedom:
$\alpha=(\nu,\sigma)$, $\Omega$ is the area of the system, ${\cal V}$ is the
screened Coulomb potential, accounting for the RPA corrections in both the
particle-hole and the electron-phonon channels, as detailed in the App. .7
below, $T$ is the temperature, $\omega$ are fermionic Matsubara frequencies
and $\mathcal{G}$ is the electronic Green’s function calculated in the normal
state. The Eq. (S15) is represented diagrammatically in the Fig. S5, where the
wavy line denotes the screened potential and the straight lines are the
Green’s functions.

Figure S5:  Diagrammatic representation of the equation (S15). The wavy line
denotes the screened potential and the straight lines are the Green’s
functions.

The Green’s functions can be generally written as:

$\displaystyle\mathcal{G}^{ij}_{\bm{r}\bm{r}^{\prime},\alpha}(i\hbar\omega)=\sum_{n\bm{k}}\frac{\Phi^{i}_{n\bm{k},\alpha}(\bm{r})\Phi^{j,*}_{n\bm{k},\alpha}(\bm{r}^{\prime})}{i\hbar\omega+\mu-
E_{n\bm{k},\alpha}},$ (S16)

where $\bm{k}$ is the wave vector in the BZ, $\mu$ is the chemical potential
and $\Phi^{i}_{n\bm{k},\alpha}$ is the Bloch eigenfunction corresponding to
the band $E_{n\bm{k},\alpha}$:

$\displaystyle\Phi^{i}_{n\bm{k},\alpha}(\bm{r})=\frac{e^{i\bm{k}\cdot\bm{r}}}{\sqrt{\Omega}}\sum_{\bm{G}}\phi^{i}_{n\bm{k},\alpha}(\bm{G})e^{i\bm{G}\cdot\bm{r}},$
(S17)

where $\bm{G}=n_{1}\bm{G}_{1}+n_{2}\bm{G}_{2}$ are reciprocal lattice vectors,
with $n_{1},n_{2}$ integers, and $\phi^{i}_{n\bm{k},\alpha}$ are numerical
eigenvectors, normalized according to:
$\sum_{\bm{G}i}\phi^{i,*}_{n\bm{k},\alpha}(\bm{G})\phi^{i}_{m\bm{k},\alpha}(\bm{G})=\delta_{nm}$.
Note that the Green’s function does not depend explicitly on the spin index,
as the continuum Hamiltonian, Eq. (1), does not.

By performing the sum over the Matsubara frequencies in the Eq. (S15), one
obtains:

$\displaystyle\Delta^{i_{1}i_{2}}_{\alpha\beta}(\bm{r}_{1},\bm{r}_{2})$
$\displaystyle=$ $\displaystyle-{\cal
V}^{scr}(\bm{r}_{1},\bm{r}_{2})\sum_{n_{1}\bm{k}_{1}n_{2}\bm{k}_{2}}\Phi^{i_{1}}_{n_{1}\bm{k}_{1},\alpha}(\bm{r}_{1})\Phi^{i_{2}}_{n_{2}\bm{k}_{2},\beta}(\bm{r}_{2})\left[\frac{f\left(-E_{n_{2}\bm{k}_{2},\beta}+\mu\right)-f\left(E_{n_{1}\bm{k}_{1},\alpha}-\mu\right)}{E_{n_{1}\bm{k}_{1},\alpha}+E_{n_{2}\bm{k}_{2},\beta}-2\mu}\right]\times$
$\displaystyle\times$
$\displaystyle\int_{\Omega}\,d^{2}\bm{r}_{3}d^{2}\bm{r}_{4}\sum_{i_{3}i_{4}}\Phi^{i_{3},*}_{n_{1}\bm{k}_{1},\alpha}(\bm{r}_{3})\Phi^{i_{4},*}_{n_{2}\bm{k}_{2},\beta}(\bm{r}_{4})\Delta^{i_{3}i_{4}}_{\alpha\beta}(\bm{r}_{3},\bm{r}_{4}),$

where: $f(\xi)\equiv\left(1+e^{\xi/K_{B}T}\right)^{-1}$ is the Fermi
distribution. Although the translational invariance is not preserved in
general, the screened potential, ${\cal V}^{scr}$, is still translationally
invariant at the moiré scale, which means: ${\cal
V}^{scr}\left(\bm{r}_{1},\bm{r}_{2}\right)={\cal
V}^{scr}\left(\bm{r}_{1}+\bm{R},\bm{r}_{2}+\bm{R}\right)$ for any Bravais
vector, $\bm{R}$, of the moiré lattice. As a consequence, ${\cal V}^{scr}$ can
be generally expressed in the Fourier basis as:

$\displaystyle{\cal
V}^{scr}\left(\bm{r}_{1},\bm{r}_{2}\right)=\frac{1}{\Omega}\sum_{\bm{q}\bm{G}_{1}\bm{G}_{2}}{\cal
V}^{scr}_{\bm{G}_{1},\bm{G}_{2}}\left(\bm{q}\right)e^{i\bm{q}\cdot(\bm{r}_{1}-\bm{r}_{2})}e^{i\bm{G}_{1}\cdot\bm{r}_{1}}e^{-i\bm{G}_{2}\cdot\bm{r}_{2}},$
(S19)

where $\bm{q}$ belongs to the moiré BZ and: ${\cal
V}^{scr}_{\bm{G}_{1},\bm{G}_{2}}\left(\bm{q}+\bm{G}\right)={\cal
V}^{scr}_{\bm{G}_{1}+\bm{G},\bm{G}_{2}+\bm{G}}\left(\bm{q}\right)$. Note that
the off-diagonal elements, $\bm{G}_{1}\neq\bm{G}_{2}$, are triggered by the
Umklapp processes, whereas ${\cal
V}^{scr}_{\bm{G}_{1},\bm{G}_{2}}\left(\bm{q}\right)\sim\delta_{\bm{G}_{1},\bm{G}_{2}}$
in the translational invariant limit. Without loss of generality, we can then
assume a similar Fourier expansion for the OP:

$\displaystyle\Delta^{i_{3}i_{4}}_{\alpha\beta}\left(\bm{r}_{3},\bm{r}_{4}\right)=\frac{1}{\Omega}\sum_{\bm{q}\bm{G}_{3}\bm{G}_{4}}\Delta^{i_{3}i_{4}}_{\alpha\beta;\bm{G}_{3},\bm{G}_{4}}\left(\bm{q}\right)e^{i\bm{q}\cdot(\bm{r}_{3}-\bm{r}_{4})}e^{i\bm{G}_{3}\cdot\bm{r}_{3}}e^{-i\bm{G}_{4}\cdot\bm{r}_{4}},$
(S20)

which implies that only the terms with $\bm{k}_{2}=-\bm{k}_{1}$ survive in the
rhs of the Eq. (.6). Defining:

$\displaystyle\Delta^{n_{1}n_{2}}_{\alpha\beta}(\bm{q})$ $\displaystyle\equiv$
$\displaystyle\int_{\Omega}\,d^{2}\bm{r}_{3}d^{2}\bm{r}_{4}\sum_{i_{3}i_{4}}\Phi^{i_{3},*}_{n_{1}\bm{q},\alpha}(\bm{r}_{3})\Phi^{i_{4},*}_{n_{2}-\bm{q},\beta}(\bm{r}_{4})\Delta^{i_{3}i_{4}}_{\alpha\beta}(\bm{r}_{3},\bm{r}_{4}),$
(S21a) $\displaystyle\tilde{\Delta}^{n_{1}n_{2}}_{\alpha\beta}(\bm{q})$
$\displaystyle\equiv$
$\displaystyle\Delta^{n_{1}n_{2}}_{\alpha\beta}(\bm{q})\times\sqrt{\frac{f\left(-E_{n_{2}-\bm{q},\beta}+\mu\right)-f\left(E_{n_{1}\bm{q},\alpha}-\mu\right)}{E_{n_{2}-\bm{q},\beta}+E_{n_{1}\bm{q},\alpha}-2\mu}},$
(S21b)

we project the Eq. (.6) on the Bloch’s eigenstates, which gives the following
equation for $\tilde{\Delta}^{n_{1}n_{2}}_{\alpha\beta}(\bm{q})$:

$\displaystyle\tilde{\Delta}^{m_{1}m_{2}}_{\alpha\beta}(\bm{k})=\sum_{n_{1}n_{2}}\sum_{\bm{q}}\Gamma^{m_{1}m_{2}}_{n_{1}n_{2};\alpha\beta}(\bm{k},\bm{q})\tilde{\Delta}^{n_{1}n_{2}}_{\alpha\beta}(\bm{q}),$
(S22)

with:

$\displaystyle\Gamma^{m_{1}m_{2}}_{n_{1}n_{2};\alpha\beta}(\bm{k},\bm{q})$
$\displaystyle=$
$\displaystyle-\frac{1}{\Omega}\sum_{\bm{G_{1}}\bm{G_{1}}^{\prime}}\sum_{\bm{G_{2}}\bm{G_{2}}^{\prime}}\sum_{i_{1}i_{2}}{\cal
V}^{scr}_{\bm{G}_{1}-\bm{G}_{1}^{\prime},\bm{G}_{2}-\bm{G}_{2}^{\prime}}\left(\bm{k}-\bm{q}\right)\phi^{i_{1},*}_{m_{1}\bm{k},\alpha}(\bm{G}_{1})\phi^{i_{2},*}_{m_{2}-\bm{k},\beta}\left(\bm{G}_{2}^{\prime}\right)\phi^{i_{1}}_{n_{1}\bm{q},\alpha}\left(\bm{G}_{1}^{\prime}\right)\phi^{i_{2}}_{n_{2}-\bm{q},\beta}(\bm{G}_{2})\times$
(S23) $\displaystyle\times$
$\displaystyle\sqrt{\frac{f\left(-E_{m_{2}-\bm{k},\beta}+\mu\right)-f\left(E_{m_{1}\bm{k},\alpha}-\mu\right)}{E_{m_{2}-\bm{k},\beta}+E_{m_{1}\bm{k},\alpha}-2\mu}}\times\sqrt{\frac{f\left(-E_{n_{2}-\bm{q},\beta}+\mu\right)-f\left(E_{n_{1}\bm{q},\alpha}-\mu\right)}{E_{n_{2}-\bm{q},\beta}+E_{n_{1}\bm{q},\alpha}-2\mu}}.$

The condition for the onset of superconductivity is that the kernel $\Gamma$
has the eigenvalue 1. This defines the critical temperature, $T_{c}$, as the
one at which the largest eigenvalue of $\Gamma$ is equal to 1.

Because the eigenfunctions do not depend explicitly on the spin index, we can
identify two different kinds of OP, depending if $\alpha$ and $\beta$ share
the same or opposite valley indices. These two OP’s describe intra-valley or
inter-valley superconductivity, respectively. In the present work we focus on
the inter-valley case, as we checked that the intra-valley superconductivity
is less robust.

In order to diagonalize the kernel $\Gamma$, we project the Eq. (S22) onto the
two bands in the middle of the spectrum, that give the main contribution to
superconductivity. For each filling, we include into $\Gamma$ the Hartree
corrections on top of the non-interacting band structure, as detailed in the
Refs. (Guinea and Walet, 2018; Cea  _et al._ , 2019).

##  .7 Computation of the screened potential

Here we provide the details concerning the calculation of the screened
potential in the reciprocal space: ${\cal
V}^{scr}_{\bm{G}_{1},\bm{G}_{2}}\left(\bm{q}\right)$. First, we introduce the
unscreened Coulomb potential:

$\displaystyle{\cal V}_{C}(\bm{q})=2\pi
e^{2}\frac{\tanh\left(d_{g}|\bm{q}|\right)}{\epsilon|\bm{q}|},$ (S24)

where $e$ is the electron charge, $d_{g}$ is the distance of the sample from a
metallic gate and $\epsilon$ is the relative dielectric constant of the
environment. In this work we use: $\epsilon=10$ and $d_{g}=40$nm. We consider
the effect of the strain in TTG by means of longitudinal acoustic phonons
coupling to the charge density via the deformation potential:
$v_{d}(\bm{r})=-Ds(\bm{r})$, where $s(\bm{r})$ is the local strain and $D$ is
the electron-phonon coupling, for which we use: $D=16$eV. We compute the
screened potential, including the RPA corrections induced by both the
particle-hole excitations and the electron-phonon coupling. In matrix
notation, the expression for the inverse of ${\cal V}^{scr}$ is given by:

$\displaystyle\left[{\cal V}^{scr}(\bm{q})\right]^{-1}=\left[\hat{{\cal
V}}_{C}(\bm{q})\right]^{-1}-\chi^{0}(\bm{q})+g\chi^{0}(\bm{q})\left[\mathbb{1}+g\chi^{0}(\bm{q})\right]^{-1}\chi^{0}(\bm{q}),$
(S25)

where the entries are indexed by the reciprocal lattice vectors,
$\left[\hat{{\cal V}}_{C}(\bm{q})\right]_{\bm{G}_{1},\bm{G}_{2}}\equiv{\cal
V}_{C}\left(\bm{q}+\bm{G}_{1}\right)\delta_{\bm{G}_{1},\bm{G}_{2}}$,
$g=\frac{D^{2}}{\lambda+2\mu}$, with $\lambda=3.25$ eVÅ-2, $\mu=9.44$ eVÅ-2
the Lamé coefficients of monolayer graphene, and $\chi^{0}(\bm{q})$ is the
static density-density response function, which is given by:

$\displaystyle\chi^{0}_{\bm{G}_{1},\bm{G}_{2}}\left(\bm{q}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{\Omega}\sum_{\bm{k}\bm{G}^{\prime}_{1}\bm{G}^{\prime}_{2}}\sum_{n_{1}n_{2}}\sum_{ij\alpha}\phi^{i}_{n_{1}\bm{k}+\bm{q},\alpha}\left(\bm{G}^{\prime}_{1}+\bm{G}_{1}\right)\phi^{i,*}_{n_{2}\bm{k},\alpha}\left(\bm{G}^{\prime}_{1}\right)\phi^{j,*}_{n_{1}\bm{k}+\bm{q},\alpha}\left(\bm{G}^{\prime}_{2}+\bm{G}_{2}\right)\phi^{j}_{n_{2}\bm{k},\alpha}\left(\bm{G}^{\prime}_{2}\right)\times$
$\displaystyle\times$
$\displaystyle\frac{f\left(E_{n_{2}\bm{k},\alpha}-\mu\right)-f\left(E_{n_{1}\bm{k}+\bm{q},\alpha}-\mu\right)}{E_{n_{2}\bm{k},\alpha}-E_{n_{1}\bm{k}+\bm{q},\alpha}}.$

To compute $\chi^{0}$ numerically, we include the 10 bands closest to the
charge neutrality point.


´

# Dissipative processes at the acoustic horizon

Maria Luisa Chiofalo  Dipartimento di Fisica, Università di Pisa, Polo
Fibonacci, Largo B. Pontecorvo 3, 56127 Pisa, Italy  INFN, Sezione di Pisa,
Largo B. Pontecorvo 3, 56127 Pisa, Italy  Dario Grasso  INFN, Sezione di Pisa,
Largo B. Pontecorvo 3, 56127 Pisa, Italy  Dipartimento di Fisica, Università
di Pisa, Polo Fibonacci, Largo B. Pontecorvo 3, 56127 Pisa, Italy  Massimo
Mannarelli  INFN, Laboratori Nazionali del Gran Sasso, Via G. Acitelli, 22,
I-67100 Assergi (AQ), Italy  Silvia Trabucco  INFN, Laboratori Nazionali del
Gran Sasso, Via G. Acitelli, 22, I-67100 Assergi (AQ), Italy  Gran Sasso
Science Institute, Viale Francesco Crispi, 7, 67100 L’Aquila

###### Abstract

A transonic fluid flow creates an acoustic hole that is the hydrodynamic
analogue of a gravitational black hole. Acoustic holes emit a detectable
thermal radiation of phonons at a characteristic Hawking temperature. We show
that viscosity arises from the thermal phonon emission associated to the
horizon entropy, shedding light on its geometrical origin and on how viscosity
acts as an effective coupling between background (metric) perturbations and
the acoustic horizon. The celebrated Kovtun, Son and Starinets (KSS) universal
lower bound $\eta/s=1/4\pi$ of the shear viscosity to entropy density ratio,
readily follows, and is extended to the longitudinal bulk viscosity.

Black holes (BHs) are perfect playgrounds for quantum mechanics and general
relativity given their understanding requires both [1]. Present quantum
technologies provide accurately controllable platforms where analogue (sonic)
BHs [2] can be realized [3, 4, 5, 6] and investigated in one-to-one
correspondence with precise theoretical predictions. The latter have the
advantage that can be based either on quantum fluid-dynamics, just hinging on
very general conservation laws and symmetries [7, 8, 9, 10], or on microscopic
theories of interacting matter based on accurately known microscopic
Hamiltonians [11], or both. The gravity analogue approach can best provide
hints on those characteristics that can be considered as universal, i.e.
independent of the underlying microscopic structure, with the advantage of
probing the essential ideas in combined, on-purpose designed, theoretical and
experimental efforts, and the disadvantage of hardly having clues on the
microscopic-dependent details.

A number of BHs concepts have been investigated within an analogue-gravity
approach, especially regarding the observation of Hawking radiation and
temperature [3, 6, 2, 12], all supporting the universal nature of these
concepts. One of the most striking predictions of universal behavior concerns
the shear viscosity-to-entropy density ratio $\eta/s$, bounded by the KSS
value $\eta/s\geq 1/4\pi$. Called after Kovtun, Son and Starinets who first
derived it within the AdS/CFT correspondence [13], the KSS bound has been
independently worked out in Rindler causal horizon in flat spacetime [14, 15].
The speed of light not appearing in the bound, the KSS conjecture seems to be
valid for all real fluids, relativistic and not, and can be extended to
classical fluids. The highly debated question here arises, on which are the
necessary conditions for the fluid to provide the minimum $\eta/s$ value, just
hinging on quantum mechanics [11].

The answer to this question is relevant when bringing general relativity and
quantum mechanics together at the black hole horizon, and also as a guidance
to account for discrepancies between microscopic theoretical predictions [11]
and experimental realizations [16]. One general answer to the question comes
from the KSS work itself, i.e. that the fluid must be strongly interacting
with no well-defined quasi-particles [11], something that is found in
relativistic heavy-ion collisions close to the deconfinement transition
temperature, see for instance [17], in pure gauge numerical simulations [18],
in $O(N)$ models [19] as well as in hadronic models [20]. So far, hydrodynamic
transport coefficients such as viscosities were mainly calculated from a
microscopic theory using Kubo-like formulas, which involve finite temperature
Green’s functions of conserved currents [21, 11, 15].

Here, we propose an entirely different, simpler and more transparent
perspective, where the KSS bound is derived within the viewpoint of kinetic
theory [22, 23, 24, 8], enhancing purely geometric considerations. Besides
widening the understanding of the essential idea and its universal nature, our
derivation may be especially useful to design analogue gravity experiments in
current platforms, where the concept can be probed. In essence, we start from
the idea that the low-energy excitations of a transonic fluid in Minkowski
spacetime can be described in terms of scalar particles embedded in an
effective gravitational background [2] with emerging acoustic metric tensor
[25, 26] determined by the fluid’s properties. In analogy with BHs, a
transonic fluid produces a sonic black hole, with excitations confined in the
supersonic region up to the horizon, where the fluid velocity equals the speed
of sound [2, 25, 26, 27, 7, 24]. The acoustic horizon emits a thermal
radiation of phonons [27, 28, 29, 30, 31, 7, 32] which are long-wavelength
sonic vibrations, as it has been confirmed both numerically [9] and in the
laboratory [3, 5] with atomic Bose-Einstein condensates.

The spontaneous phonon emission process at the horizon produces an
irreversible heat increase at the expenses of the bulk fluid kinetic energy,
which can be described in terms of effective shear and bulk viscosities. We
determine the viscosity coefficients, eventually finding that their ratio to
the phonon-entropy density equals $1/4\pi$ [13]. This result is obtained
irrespective of the detailed phonon emission mechanism. To provide a
microscopic description of viscosity and to corroborate our findings, we
derive the KSS bound employing the analogy between hydrodynamics and gravity
to describe the spontaneous emission process by phonon tunneling.

Theoretical setup \- We consider an inviscid, barotropic, irrotational fluid
at vanishing temperature, in a Minkowski spacetime. Adopting from now on units
$\hbar=c=k_{B}=1$, the acoustic metric relevant for the phonon propagation is

$g_{\mu\nu}=\eta_{\mu\nu}+\left(c_{s}^{2}-1\right)v_{\mu}v_{\nu}\,,$ (1)

where $\eta_{\mu\nu}=$(1,-1,-1,-1) is the standard Minkowski metric, $c_{s}$
is the sound speed and $v_{\mu}$ is the fluid four velocity. For definiteness,
we assume that the velocity is oriented along the negative $x$-axis and is
given by

$v=c_{s}-Cx+ky\,,$ (2)

where $C$ and $k$ are two parameters, with $k\ll C$ to ensure a negligible
flow along the $y$ direction. The last term on the r.h.s. of (2) represents a
small shear perturbation to the fluid flow. The horizon position corresponds
to the place where the fluid velocity equals the speed of sound and it is
schematically depicted in Fig. 1, with $\theta\simeq-k/C$ the tilting angle of
the horizon with respect to the $y$-axis.

Figure 1: Two dimensional representation of the horizon determined by the
fluid velocity profile given in Eq. (2). The fluid velocity increases from
right to left and from bottom to top. The blue line corresponds to the horizon
and the coordinate axes $w_{1}$ and $w_{2}$ are orthogonal and parallel to the
horizon, respectively.

By the same kinetic approach derived in [10], it can be shown that the
acoustic horizon determined by the velocity profile in (2) emits phonons at
the Hawking temperature

$T=\frac{1}{2\pi}\left.\left(\frac{c_{s}-|v|}{1-c_{s}|v|}\right)^{\prime}\
\right|_{H}\simeq\frac{C}{2\pi}\,,$ (3)

where the prime indicates the derivative with respect to the $w_{1}$
coordinate, and the last expression is obtained in the non-relativistic limit
and at leading order in $k/C$. Expression (3) reproduces the standard form of
$T$ in the non-relativistic limit, see for instance [7].

In the following we limit our analysis to a small neighbourhood $\ell$ of the
horizon position, satisfying the condition

$L_{c}\ll\ell\ll\lambda_{T}\,,$ (4)

with $L_{c}$ the theory cutoff accounting for the underlying microscopic
structure invisible to our treatment; in acoustic holes realized with Bose-
Einstein condensates it is of the same order of the healing length of the
condensate. The length scale

$\lambda_{T}\equiv\frac{c_{s}}{2\pi T(1-c_{s}^{2})}\,,$ (5)

is a thermal wavelength associated to the Hawking temperature. In turn, taking
into account the above expression of the Hawking temperature, the condition
$\ell\ll\lambda_{T}$ ensures that we consider a small neighborhood of the
horizon, where the fluid velocity is close to $c_{s}$. The considered space
region is the relevant one for studying the spontaneously emitted phonons,
which can be intuitively traced back to thermal effects arising from
dissipative processes occurring close to the horizon.

Bulk and shear viscosities\- We determine the ratio of the bulk and shear
viscosities by equating the viscous stress tensor [33] of the background to
the energy-momentum tensor of the phonons, under the working hypothesis that
both the shear, $\eta$, and bulk, $\zeta$, viscosities are only due to the
spontaneous emission of phonons at the horizon and are thus defined only in a
thin layer around its position.

The viscous stress tensor pertinent to the velocity (2) is

$\sigma_{ik}^{\prime}=\eta\left(\partial_{i}v_{k}+\partial_{k}v_{i}\right)+\zeta_{\text{eff}}\delta_{ix}\delta_{kx}{\bm{\nabla}}\cdot{\bm{v}}\,,$
(6)

where $\zeta_{\text{eff}}=\frac{4}{3}\eta+\zeta$ is the effective bulk
viscosity: in fluid-dynamics it accounts for longitudinal dissipative
processes. At the leading order in $k/C$, we expect that only its $xx$, $yx$
and $xy$ components are nonzero. This is due to the fact that the rotational
symmetry is broken by the acoustic horizon and that the horizon can only
produce a pressure in the $w_{1}$ direction, see Fig. 1. This means, for
instance, that the pressure along the $y$ direction is of order $(k/C)^{2}$.
We shall explicitly see it after using kinetic theory.

The general expression of the phonon energy-momentum tensor is

$T^{\mu}_{\nu}=\int p^{\mu}p_{\nu}f_{\text{ph}}(x,p)d{\cal P}\,,$ (7)

where $d{\cal P}$ is the appropriate covariant momentum measure [22, 23] and
$f_{\text{ph}}$ is the phonon distributon function. Close to the horizon, it
can be cast in the form [10]

$f_{\text{ph}}(x,p)=\frac{1}{e^{E_{+}/T}-1}\delta(p_{w_{2}})\delta(p_{z})\left(\frac{2\pi}{L_{c}}\right)^{2}\,,$
(8)

where

$E_{+}=-p_{x}v+c_{s}\sqrt{p_{x}^{2}+p_{y}^{2}+p_{z}^{2}}\,,$ (9)

is the phonon dispersion law in the non-relativistic limit and
$p_{w_{2}}=p_{y}\cos\theta+p_{x}\sin{\theta}$ is the momentum parallel to the
horizon. The distribution function takes into account that phonons can only be
emitted orthogonally to the horizon. Thus, we have that

$p_{z}=0\,,\qquad p_{y}\simeq-\frac{k}{C}p_{x}\,,$ (10)

and therefore the phonon dispersion law reads

$E_{+}=p_{x}Cw_{1}(1+{\cal O}(k^{2}/C^{2}))\,.$ (11)

Upon substituting this expression in (8), we can readily evaluate the leading
order components of the energy-momentum tensor. The energy density

$\epsilon_{\text{ph}}=\sqrt{-g}T^{0}_{0}\simeq\frac{T}{24L_{c}^{2}w_{1}}\,,$
(12)

takes an expression similar to the one derived in [10], with $L_{c}^{2}w_{1}$
the small volume orthogonal to the horizon. The pressure along the
$x$-direction is

$P_{x}=\sqrt{-g}T^{x}_{x}\simeq\frac{T}{24L_{c}^{2}w_{1}}\,,$ (13)

and it is therefore equal to the energy density. The $T^{y}_{x}$ component
leads instead to a form of transverse pressure:

$P^{y}_{x}=\sqrt{-g}T^{y}_{x}\simeq-\frac{k}{C}P_{x}\,.$ (14)

Finally, the pressure along the $y$-axis turns out to be:

$P_{y}=\sqrt{-g}T^{y}_{y}\propto(k/C)^{2}\,,$ (15)

and thus, as anticipated, negligible to leading order in $k/C$. The same
expressions can be obtained by rotating the reference frame by $\theta$ angle.
We now implement our key idea that the spontaneous phonon emission process at
the horizon produces the irreversible heat increase at the expenses of the
bulk fluid kinetic energy described by $\eta$ and $\zeta_{\text{eff}}$. We
thus equate the viscous stress fluid contributions with the corresponding
phonon terms

$\sigma_{xx}^{\prime}=P_{x}\,,\qquad\sigma_{xy}^{\prime}=P^{y}_{x}\,.$ (16)

Using the thermodynamic relation $\epsilon_{\text{ph}}=-P_{x}+Ts_{\text{ph}}$,
that is valid at the leading order in $k/C$, we finally have that:

$\frac{\zeta_{\text{eff}}}{s_{\text{ph}}}=\frac{\eta}{s_{\text{ph}}}=\frac{1}{4\pi}\,.$
(17)

That is, the KSS bound is saturated for both the shear and the effective bulk
viscosities [13]. Note that we have obtained this result without any reference
to the microscopic emission mechanism of phonons and neglecting any phonon
self-interaction. The viscosities only arise as a consequence of the
spontaneous phonon emission at the horizon.

Phonon tunneling \- We have seen how the viscous coefficients rule the
macroscopic behaviour of the background fluid, i.e. its response to a
perturbation of the velocity field (or of the sound velocity). We now enter
the nature of this process and show that it can be interpreted as the
consequence of phonons tunneling the acoustic horizon. World-lines of phonons
are null geodesics of the acoustic metric, see (1), and their dispersion law
depends on the position with respect to the horizon, see Eq. (11) and the
discussion in [10].

In order to grasp the underlying physics, we use a semiclassical description
of the Hawking emission process [34, 35]. The WKB tunneling amplitude is
determined by the exponential of the imaginary part of the action
$\text{Im\,S}=\text{Im}\int_{x_{\text{in}}}^{x_{\text{out}}}p(x)dx$, where
$x_{\text{in/out}}$ are points inside/outside the horizon and $p(x)$ is the
momentum conjugate to the $x$-coordinate; it corresponds to the phonon
momentum during the emission process. As in [34], we can use the Hamilton-
Jacobi equation $\displaystyle\dot{x}={d\omega}/{dp}$, with $\omega$ the
energy of the outgoing phonon with momentum $p(x)$ to change the integration
variable. The integral can then be evaluated by using the expression of the
null geodesic condition for phonons moving upstream ($x>0$) (see [10])

$\dot{x}=\frac{c_{s}-v}{1-c_{s}v}\,.$ (18)

After a change of variables, we obtain

$\text{Im\,S}=\text{Im}\int_{x_{\text{in}}}^{x_{\text{out}}}dx\int_{0}^{\omega}d\omega\frac{1-c_{s}v}{c_{s}-v}\,,$
(19)

where we took into account that the emitted phonon has initially vanishing
energy. Upon changing the order of the integrals and using the velocity
profile (2), we obtain that the probability of quanta emission with energy
$\omega$ is given by

$P(\omega)\,\propto\,e^{-2\,\text{Im\,S}}=e^{-\omega/{T}}\,,$ (20)

which is the Boltzmann distribution with the Hawking temperature given in Eq.
(3). This result is in agreement with the fluctuation-dissipation theorem
linking viscous dissipation to fluctuations of a thermal equilibrium state.
The thermal phonon emission at the horizon results in an irreversible decrease
of the kinetic energy of the fluid [24, 36]. In turn, we find that this
produces a correction to the Boltzmann distribution of the same form obtained
in [34], with the black hole mass replaced by the fluid mass times the speed
of sound squared. Notice that introducing the backreaction of the acoustic
horizon would not alter the main idea we are focusing on; we leave the
discussion of the backreaction to a future publication.

We instead proceed now to describe the dissipative processes associated to the
spontaneous phonon emission. In particular, we determine the viscosity
coefficients by equating the kinetic energy loss of the fluid with the phonon
energy gain, corresponding to the heat emitted by the horizon. The total
energy of the phonons in a small volume close to the horizon having as basis a
patch of area $L_{c}^{2}$ is

$E_{\text{ph}}=\int\epsilon_{\text{ph}}d^{3}x\simeq
L_{c}^{2}\int_{L_{c}}^{w_{1}}\epsilon_{\text{ph}}dw_{1}=\frac{T}{24}\log\left(\frac{w_{1}}{L_{c}}\right)\,,$
(21)

where the phonon energy density is given in (12). Thus, the emission power of
the horizon is

$\dot{E}_{\text{ph}}=\frac{T}{24}\frac{\dot{w}_{1}}{w_{1}}=\epsilon_{\text{ph}}L_{c}^{2}\dot{w}_{1}\,,$
(22)

where

$\dot{w}_{1}=w_{1}\frac{C^{2}+k^{2}}{C}-\frac{\dot{k}}{C}y\,,$ (23)

which is obtained after taking the non-relativistic limit of (18), considering
the leading order in $k/C$, and taking into account that the small external
perturbation $k$ can be time-dependent.

The rate of kinetic energy loss of the fluid due to dissipative processes is
instead

$\displaystyle\dot{E}_{\text{fluid}}$ $\displaystyle=-\int
d^{3}x\frac{\sigma_{ik}^{\prime}}{2}\left(\partial_{i}v_{k}+\partial_{k}v_{i}\right)$
$\displaystyle\simeq-
L_{c}^{2}w_{1}\frac{\sigma_{ik}^{\prime}}{2}\left(\partial_{i}v_{k}+\partial_{k}v_{i}\right)\,,$
(24)

where in the last expression we only considered the relevant thin layer close
to the horizon. In the stationary case, that is for $\dot{k}=0$, upon equating
the kinetic energy loss of the fluid with the energy gain of the phonon gas we
obtain once again the same bounds as in (17). Including the time-dependence of
the perturbation, the bulk viscosity remains unchanged, while the leading-
order contribution for the shear viscosity is

$\frac{\eta}{s_{\text{ph}}}=\frac{1}{4\pi}\left(1+\frac{\dot{k}}{2\pi
Tk}\right)\,,$ (25)

valid for low-frequency perturbations, i.e. for $\dot{k}/k\ll T$.

Relating $\eta/s$ to quantum entanglement \- Finally, we focus on the ratio
between the shear viscosity coefficient and the entropy surface density of the
horizon. The shear perturbation of the background fluid (hence of effective
spacetime) turns into a deformation of the horizon producing a variation of
the cross-sectional area element of the horizon projected onto the holographic
screen ($yz$ plane). We notice that similarly to BHs, the presence of the
acoustic horizon produces a spatial partition of the vacuum state giving rise
to an entanglement entropy [37]. In our effective framework, this entropy must
be normalized by the cutoff scale $L_{c}$ hence it scales with the horizon
area: $\displaystyle S_{H}=\kappa{A_{H}}/{L_{c}^{2}}$. For $\kappa=1/4$ and
identifying $L_{c}$ with the Planck length, this expression coincides with the
Bekenstein-Hawking entropy of BHs [38], which indeed has been interpreted as
an entanglement entropy [39, 40] (see also [41] and ref.s therein).

We consider the case of a time-dependent horizon tilt with ${\dot{k}}\neq 0$,
which can be originated by an external perturbation of the fluid velocity or
of the sound speed. This gives rise to a time dependent entropy surface
density as seen on a screen perpendicular to the fluid flow. Using
$y=x_{0}-{kx}/{C}$, see Fig. 1, it readily follows that the area change is
given by

$\frac{dA_{H}}{dt}=\frac{k\dot{k}}{C^{2}}\,dydz\,,$ (26)

and the corresponding entropy change is thus
$\displaystyle\dot{S}_{H}=\kappa\dot{A}_{H}/L_{c}^{2}$.

The deformation of the horizon has its counterpart on the phonon world-lines
acquiring a shear. In analogy with the membrane paradigm for BHs, and under
the condition expressed by (4), we identify these phonons as a fluid endowed
with a shear and undergoing quantum tunneling in a “effectively” 2+1
dimensional thin layer about the horizon. The entropy change of the horizon is
transferred to the emitted phonons, giving rise to a shear viscosity of the
fluid. Following [15], we may call it the entanglement viscosity of the
horizon. Now, the phonons are emitted along $w_{1}$, see Fig. 1, with velocity
as in (23) and components

$V_{x}={\dot{w}_{1}}\frac{C}{\sqrt{C^{2}+k^{2}}}\quad
V_{y}={\dot{w}_{1}}\frac{k}{\sqrt{C^{2}+k^{2}}}\quad V_{z}=0\,.$ (27)

Thus, the shear of their flow can be readily computed to be

$\left(\partial_{i}V_{k}+\partial_{k}V_{i}\right)^{2}\simeq
2\frac{k\dot{k}}{C}\,.$ (28)

Hence, using the constitutive equations

$\frac{ds_{H}}{dt}=\frac{2\eta_{H}}{T}\left(\partial_{i}V_{k}+\partial_{k}V_{i}\right)^{2}\,,$
(29)

we obtain

$\eta_{H}=\frac{\kappa}{L_{c}^{2}}\frac{1}{4\pi}\,,$ (30)

where $\eta_{H}$ has dimensions of $1/{\rm length}^{2}$. Since the horizon
entropy density is $s_{H}=\kappa/L_{c}^{2}$, we readily obtain

$\frac{\eta_{H}}{s_{H}}=\frac{1}{4\pi}\,,$ (31)

which is the KSS ratio. This result strengthens the hypothesis that the KSS
ratio could be a fundamental holographic property of space-time (rather than
only of the AdS/BH solutions).

Conclusions - We have shown that the universal Kovtun-Son-Starinets [13] bound
of the shear viscosity to entropy density ratio is saturated in the region
close to the acoustic horizon of a transonic flow. The same happens for the
longitudinal bulk viscosity coefficient. We have first obtained these results
by a general kinetic theory approach where phonons spontaneously emitted at
the horizon irreversibly produce heat. This derivation is quite general, since
it is not based on a microscopic emission process. Then, we have obtained the
same results considering a microscopic description of viscosity where the
emission process is due to phonons tunneling the acoustic horizon. Finally, we
have as well obtained the limiting value of $\eta/s$ considering the variation
of the surface area of the horizon induced by a time dependent perturbation.
The latter derivation bears some resemblance with the one proposed for a
Rindler causal horizon in flat spacetime [15]. Since in all considered cases
there is no known holographic (gauge/gravity) duality, a natural explanation
for this viscosity is in the peculiar properties of quantum entanglement
rather than in the existence of a gravity dual [13].

Our approach makes clear that the viscosity arises from the thermal (Hawking)
phonon emission associated to the horizon entanglement entropy, shedding light
on its geometrical origin and on how viscosity acts as an effective coupling
between background (metric) perturbations and the acoustic horizon. This
perspective helps conceptualizing proper analogue experiments, where the
saturation of the KSS bound can be explored in controllable manners. For
example, the simplest situation can be envisioned to be that of a Bose-
Einstein condensate in effectively 2D circular geometry [6], by impressing a
fluid (or, equivalently, a sound) velocity profile according to (2) and Fig.
1. The entropy density could be measured as in [6]. As to the bulk and shear
viscosities, they could be inferred via their effect on the horizon: in our
understanding, $\zeta_{\text{eff}}$ displaces it by shrinking the supersonic
region, while $\eta$ restores its position by reducing the horizon tilting in
Fig. 1. In such experiments, our derivation can help disentangling
interaction-driven and geometric originated effects.

## References

  * Hawking [1975] S. W. Hawking, Particle Creation by Black Holes, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Unruh [1981] W. Unruh, Experimental black hole evaporation, Phys.Rev.Lett. 46, 1351 (1981).
  * Steinhauer [2016] J. Steinhauer, Observation of quantum Hawking radiation and its entanglement in an analogue black hole, Nature Physics 12, 959–965 (2016).
  * Clark  _et al._ [2017] L. W. Clark, A. Gaj, L. Feng, and C. Chin, Collective emission of matter-wave jets from driven Bose-Einstein condensates, Nature 551, 356–359 (2017).
  * Muñoz de Nova  _et al._ [2019] J. R. Muñoz de Nova, K. Golubkov, V. I. Kolobov, and J. Steinhauer, Observation of thermal Hawking radiation and its temperature in an analogue black hole, Nature 569, 688 (2019), arXiv:1809.00913 [gr-qc] .
  * Hu  _et al._ [2019] J. Hu, L. Feng, Z. Zhang, and C. Chin, Quantum Simulation of Coherent Hawking-Unruh Radiation, Nature Physics 15, 785–789 (2019).
  * Barcelo  _et al._ [2005] C. Barcelo, S. Liberati, and M. Visser, Analogue gravity, Living Rev. Rel. 8, 12 (2005), arXiv:gr-qc/0505065 .
  * Mannarelli and Manuel [2008] M. Mannarelli and C. Manuel, Transport theory for cold relativistic superfluids from an analogue model of gravity, Phys. Rev. D 77, 103014 (2008), arXiv:0802.0321 [hep-ph] .
  * Carusotto  _et al._ [2008] I. Carusotto, S. Fagnocchi, A. Recati, R. Balbinot, and A. Fabbri, Numerical observation of Hawking radiation from acoustic black holes in atomic Bose–Einstein condensates, New J. Phys. 10, 103001 (2008), 0803.0507 .
  * Mannarelli  _et al._ [2021] M. Mannarelli, D. Grasso, S. Trabucco, and M. L. Chiofalo, Hawking temperature and phonon emission in acoustic holes, Phys. Rev. D 103, 076001 (2021), arXiv:2011.00019 [gr-qc] .
  * Zwerger [2011] W. Zwerger,  _The BCS-BEC Crossover and the Unitary Fermi Gas_, Lecture Notes in Physics (Springer Berlin Heidelberg, 2011).
  * Visser [1998] M. Visser, Acoustic black holes: Horizons, ergospheres, and Hawking radiation, Class. Quant. Grav. 15, 1767 (1998), arXiv:gr-qc/9712010 .
  * Kovtun  _et al._ [2005] P. Kovtun, D. T. Son, and A. O. Starinets, Viscosity in strongly interacting quantum field theories from black hole physics, Phys. Rev. Lett. 94, 111601 (2005), arXiv:hep-th/0405231 .
  * Chirco and Liberati [2010] G. Chirco and S. Liberati, Non-equilibrium thermodynamics of spacetime: The role of gravitational dissipation, Physical Review D 81, 024016 (2010), arXiv:0909.4194 [gr-qc] .
  * Chirco  _et al._ [2010] G. Chirco, C. Eling, and S. Liberati, The universal viscosity to entropy density ratio from entanglement, Phys. Rev. D 82, 024010 (2010), arXiv:1005.0475 [hep-th] .
  * Cao  _et al._ [2011] C. Cao, E. Elliott, J. Joseph, H. Wu, J. Petricka, T. Schäfer, and J. E. Thomas, Universal Quantum Viscosity in a Unitary Fermi Gas, Science 331, 58 (2011), arXiv:1007.2625 [cond-mat.quant-gas] .
  * Shuryak [2017] E. Shuryak, Strongly coupled quark-gluon plasma in heavy ion collisions, Rev. Mod. Phys. 89, 035001 (2017), arXiv:1412.8393 [hep-ph] .
  * Meyer [2008] H. B. Meyer, A Calculation of the bulk viscosity in SU(3) gluodynamics, Phys. Rev. Lett. 100, 162001 (2008), arXiv:0710.3717 [hep-lat] .
  * Romatschke [2021] P. Romatschke, Shear Viscosity at Infinite Coupling: A Field Theory Calculation, Phys. Rev. Lett. 127, 111603 (2021), arXiv:2104.06435 [hep-th] .
  * Rais  _et al._ [2020] J. Rais, K. Gallmeister, and C. Greiner, Shear viscosity to entropy density ratio of Hagedorn states, Phys. Rev. D 102, 036009 (2020), arXiv:1909.04522 [hep-ph] .
  * Enss  _et al._ [2011] T. Enss, R. Haussmann, and W. Zwerger, Viscosity and scale invariance in the unitary fermi gas, Annals of Physics 326, 770 (2011).
  * Lindquist [1966] R. W. Lindquist, Relativistic transport theory, Annals of Physics 37, 487 (1966).
  * Stewart [1969] J. Stewart,  _Lecture Notes in Physics_, Lecture Notes in Physics No. v. 10 (Springer-Verlag, 1969).
  * Volovik [2003] G. E. Volovik, What can the quantum liquid say on the brane black hole, the entropy of extremal black hole and the vacuum energy?, Found. Phys. 33, 349 (2003), arXiv:gr-qc/0301043 .
  * Bilic [1999] N. Bilic, Relativistic acoustic geometry, Class. Quant. Grav. 16, 3953 (1999), arXiv:gr-qc/9908002 .
  * Visser and Molina-Paris [2010] M. Visser and C. Molina-Paris, Acoustic geometry for general relativistic barotropic irrotational fluid flow, New J. Phys. 12, 095014 (2010), arXiv:1001.1310 [gr-qc] .
  * Brout  _et al._ [1995] R. Brout, S. Massar, R. Parentani, and P. Spindel, A Primer for black hole quantum physics, Phys. Rept. 260, 329 (1995), arXiv:0710.4345 [gr-qc] .
  * Corley [1998] S. Corley, Computing the spectrum of black hole radiation in the presence of high frequency dispersion: An Analytical approach, Phys. Rev. D 57, 6280 (1998), arXiv:hep-th/9710075 .
  * Saida and Sakagami [2000] H. Saida and M.-a. Sakagami, Black hole radiation with high frequency dispersion, Phys. Rev. D 61, 084023 (2000), arXiv:gr-qc/9905034 .
  * Himemoto and Tanaka [2000] Y. Himemoto and T. Tanaka, Generalization of the model of Hawking radiation with modified high frequency dispersion relation, in  _9th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Gravitation and Relativistic Field Theories (MG 9)_ (2000) pp. 1384–1394.
  * Unruh and Schutzhold [2005] W. G. Unruh and R. Schutzhold, On the universality of the Hawking effect, Phys. Rev. D 71, 024028 (2005), arXiv:gr-qc/0408009 .
  * Balbinot  _et al._ [2005a] R. Balbinot, A. Fabbri, S. Fagnocchi, and R. Parentani, Hawking radiation from acoustic black holes, short distance and back-reaction effects, Riv. Nuovo Cim. 28, 1 (2005a), arXiv:gr-qc/0601079 .
  * Landau and Lifshitz [1980] L. D. Landau and E. M. Lifshitz,  _Course of theoretical physics, Pergamon International Library of Science, Technology, Engineering and Social Studies_ (Pergamon Press, 1980).
  * Parikh and Wilczek [2000] M. K. Parikh and F. Wilczek, Hawking radiation as tunneling, Phys. Rev. Lett. 85, 5042 (2000), arXiv:hep-th/9907001 .
  * Volovik [1999] G. E. Volovik, Simulation of Painleve-Gullstrand black hole in thin He-3 - A film, JETP Lett. 69, 705 (1999), arXiv:gr-qc/9901077 .
  * Balbinot  _et al._ [2005b] R. Balbinot, S. Fagnocchi, A. Fabbri, and G. P. Procopio, Backreaction in acoustic black holes, Phys. Rev. Lett. 94, 161302 (2005b), arXiv:gr-qc/0405096 .
  * Finazzi and Carusotto [2014] S. Finazzi and I. Carusotto, Entangled phonons in atomic Bose-Einstein condensates, Phys. Rev. A 90, 033607 (2014), arXiv:1309.3414 [cond-mat.quant-gas] .
  * Bekenstein [1973] J. D. Bekenstein, Black holes and entropy, Phys. Rev. D 7, 2333 (1973).
  * Sorkin [1984] R. D. Sorkin, ”On the Entropy of the Vacuum outside a Horizon”, in  _10th International Conference on General Relativity and Gravitation_ , Vol. 2 (1984) pp. 734–736, arXiv:1402.3589 [gr-qc] .
  * Srednicki [1993] M. Srednicki, Entropy and area, Phys. Rev. Lett. 71, 666 (1993), arXiv:hep-th/9303048 .
  * Jacobson [2016] T. Jacobson, Entanglement Equilibrium and the Einstein Equation, Phys. Rev. Lett. 116, 201101 (2016), arXiv:1505.04753 [gr-qc] .