# A characterization of differential operators in the ring of complex
polynomials

Włodzimierz Fechner and Eszter Gselmann

###### Abstract.

The paper aims to provide a full characterization of all operators
$T\colon\mathscr{P}(\mathbb{C})\to\mathscr{P}(\mathbb{C})$ acting on the space
of all complex polynomials that satisfy the Leibniz rule

$T(f\cdot g)=T(f)\cdot g+f\cdot T(g)$

for all $f,g\in\mathscr{P}(\mathbb{C})$. We do not assume the linearity of
$T$. As we will see, contrary to the well-known theorems for function spaces
there are many other solutions here, not only differential operators. From our
main result, we also derive two corollaries, showing that in some special
cases operators that satisfy the Leibniz rule have some particular form.

###### Key words and phrases:

differential operators, operator equations, Leibniz rule, ring of complex
polynomials, derivations

###### 1991 Mathematics Subject Classification:

Primary: 47B92; Secondary: 13B25, 39B22, 47J05, 47A62, 47B38

##  1\. Introduction

Characterizations of basic operators, such as the first-order differential
operator or the Fourier transform by one of its fundamental properties and
without assuming linearity, are one of the most interesting problems in modern
analysis. Let us briefly recall some most important from our point of view
results in this direction.

For instance, if $\mathscr{S}(\mathbb{R}^{N})$ denotes the Schwartz space of
‘rapidly’ decreasing functions $f\colon\mathbb{R}^{n}\to\mathbb{C}$, then in
Artstein-Avidan–Faifman–Milman [1], the authors showed that any _bijective_
transformation
$T\colon\mathscr{S}(\mathbb{R}^{n})\to\mathscr{S}(\mathbb{R}^{n})$ that also
fulfills

$T(f\cdot g)=T(f)\ast
T(g)\qquad\left(f,g\in\mathscr{S}(\mathbb{R}^{N})\right)$

is just a slight modification of the Fourier transform.

In addition, we can mention a characterization theorem of higher order
differential operators due to Peetre [6], see also [2, 3]. Let $n$ be a
positive integer and let $\mathscr{C}^{\infty}_{0}(\mathbb{R}^{n})$ denote the
set of infinitely many times differentiable functions defined on
$\mathbb{R}^{n}$ that vanish at infinity, further let
$\mathscr{C}_{b}(\mathbb{R}^{n})$ stand for the set of all those continuous
functions defined on $\mathbb{R}^{n}$ that are bounded. If the _linear_
operator
$T\colon\mathscr{C}^{\infty}_{0}(\mathbb{R}^{n})\to\mathscr{C}_{b}(\mathbb{R}^{n})$
fulfills

$\mathrm{supp}(Tf)\subset\mathrm{supp}(f)$

for all $f\in\mathscr{C}^{\infty}_{0}(\mathbb{R}^{n})$, then there exists a
positive integer $N$ and for all $\alpha\in\mathbb{N}^{n}$, $|\alpha|\leq N$,
there exists continuous functions $c_{\alpha}$ on $\mathbb{R}^{N}$ such that

$(Tf)(x)=\sum_{|\alpha|\leq N}c_{\alpha}(x)D^{\alpha}f(x)$

for all $f\in\mathscr{C}^{\infty}_{0}(\mathbb{R}^{n})$ and
$x\in\mathbb{R}^{n}$, where

$D^{\alpha}f(x)=\frac{\partial^{|\alpha|}}{\partial
x_{1}^{\alpha_{1}}\cdots\partial x_{n}^{\alpha_{n}}}f(x_{1},\ldots,x_{n}),$

and $x=(x_{1},\ldots,x_{n})\in\mathbb{R}^{n}$.

The main objective of our paper is to prove similar characterization theorems
in polynomial rings. More precisely, we will study here those operators
$T\colon\mathscr{P}(\mathbb{K})\to\mathscr{P}(\mathbb{K})$ that fulfill the
Leibniz rule, i.e.

$T(f\cdot g)=fT(g)+T(f)g\qquad\left(f,g\in\mathscr{P}(\mathbb{K})\right).$

We would like to specifically emphasize that in this work we do not impose any
additional conditions (such as linearity or bijectivity) on the operator $T$.

As we will see, unlike the well-known theorems for function spaces (see
Theorem 1), there are many other solutions here, not only differential
operators. Moreover, if $D^{\alpha}$ is a differential operator acting on
$\mathbb{R}$ or on $\mathbb{C}$, then

$\mathrm{deg}\left(D^{\alpha}p\right)\leq\mathrm{deg}(p)$

holds for any polynomial $p$. Thus after presenting our main result, we focus
on determining the solutions of the Leibniz rule that also have this property,
see Corollaries 1 and 2.

###  1.1. Leibniz rule in $\mathscr{C}^{k}(I)$

Let $I\subset\mathbb{R}$ be an open set, $k$ be a nonnegative integer let

$\mathscr{C}^{k}(I)=\left\\{f\colon I\to\mathbb{R}\,|\,\text{$f$ is $k$-times
continuously differentiable on $I$}\right\\}.$

We denote the space of continuous functions also by $\mathscr{C}(I)$, instead
of writing $\mathscr{C}^{0}(I)$, further we put

$\mathscr{C}^{\infty}(I)=\bigcap_{k}\mathscr{C}^{k}(I).$

On spaces $\mathscr{C}^{k}(I)$ the solutions of the Leibniz rule were
determined by König and Milman [4]. Let us quote their result.

######  Theorem 1 (Leibniz rule on $\mathscr{C}^{k}(I)$).

Let $I\subset\mathbb{R}$ be an open set and $k$ be a nonnegative integer.
Suppose that the operator $T\colon\mathscr{C}^{k}(I)\to\mathscr{C}(I)$
satisfies the Leibniz rule, i.e.,

(1) $T(f\cdot g)=T(f)\cdot g+f\cdot T(g)$

holds for all $f,g\in\mathscr{C}^{k}(I)$. Then there exists continuous
functions $c,d\in\mathscr{C}(I)$ such that

$T(f)=\begin{cases}cf\ln(|f|)+df^{\prime},&\text{ if }k\geq 1\\\
cf\ln(|f|),&\text{ if
}k=0.\end{cases}\qquad\left(f\in\mathscr{C}^{k}(I)\right).$

Conversely, any map $T$ defined by the above formula fulfils the Leibniz rule
on $\mathscr{C}^{k}(I)$.

The proof of this theorem consists of two steps. At first, one has to show
that the operator $T$ is ‘localized on intervals’. In other words, the
operator $T$ is defined pointwise in that sense that there is a function
$F\colon I\times\mathbb{R}^{n+1}\to\mathbb{R}$ such that for all
$f\in\mathscr{C}^{k}(I)$ and $x\in I$ we have

$T(f(x))=F(x,f(x),\ldots,f^{(k)}(x)).$

At this point, no regularity of $F$ is known. Observe that then the Leibniz
rule is equivalent to a functional equation for the representing function $F$.
The second step is to analyze the structure of $F$ and to prove the continuity
of the coefficient functions occurring there, by using the fact that the image
space of the operator $T$ is $\mathscr{C}(I)$.

###  1.2. Derivations

The Leibniz rule, that is, equation (1) plays a key role not only in operator
theory but also in algebra. Here we recall some basic facts from the monograph
Kuczma [5].

######  Definition 1.

Let $(Q;+,\cdot)$ be a commutative ring, and let $(P;+,\cdot)$ be a subring of
$(Q;+,\cdot)$. A function $f:P\to Q$ is called a _derivation_ if it satisfies
both the equations

$f(x+y)=f(x)+f(y),$

and also

$f(xy)=xf(y)+yf(x)$

for all $x,y\in P$.

If $f$ satisfies the second equation, then we will call it the Leibniz
mapping.

######  Example 1.

Let $(F;+,\cdot)$ be a field, and let $P=Q=F\left[x\right]$ be the ring of
polynomials with coefficients from $F$. Let the function $f:F\left[x\right]\to
F\left[x\right]$ be defined as $f(p)=p^{\prime}$, where $p^{\prime}$ is the
derivative of $p$. We have clearly

$f(p+q)=(p+q)^{\prime}=p^{\prime}+q^{\prime}=f(p)+f(q)\,,$
$f(pq)=(pq)^{\prime}=pq^{\prime}+qp^{\prime}=pf(q)+qf(p)\,.$

Consequently, $f$ is a derivation.

######  Example 2.

Let $(F;+,\cdot)$ be a field, and suppose that we are given a derivation
$f:F\to F$. We define a function $f_{0}:F\left[x\right]\to F\left[x\right]$ as
follows. If $p\in F\left[x\right]$, $p(x)=\sum\limits^{n}_{k=0}a_{k}x^{k}$,
consider the polynomial $f_{0}(p)$ defined by (in sequel denoted by $p^{f}$)

$f_{0}(p)=p^{f}(x)=\sum\limits^{n}_{k=0}f(a_{k})x^{k}.$

Then $f_{0}\colon F[x]\to F[x]$ is a derivation.

The derivations described in the above two examples have a fundamental
importance for we have the following lemma, which is a modification of Kuczma
[5, Lemma 14.2.2].

######  Lemma 1.

Let $(K;+,\cdot)$ be a field, and let $(F;+,\cdot)$ be a subfield of
$(K;+,\cdot)$, and let $f:F\to K$ be a derivation. Then we have, for every
$a\in F$ and every polynomial $p\in F\left[x\right]$,

$f\big{(}p(a)\big{)}=p^{f}(a)+f(a)p^{\prime}(a).$

The question of whether there exists a nontrivial, i.e. nonzero derivation on
fields with characteristic zero is far from being obvious but can be answered
affirmatively. The following theorem (see Kuczma [5, Lemma 14.2.4]) plays a
key role in addressing this problem.

######  Theorem 2.

Let $(K;+,\cdot)$ be a field of characteristic zero, let $(F;+,\cdot)$ be a
subfield of $(K;+,\cdot)$, let $S$ be an algebraic base of $K$ over $F$, if it
exists, and let $S=\varnothing$ otherwise. Let $f:F\to K$ be a derivation.
Then, for every function $u:S\to K$, there exists a unique derivation $g:K\to
K$ such that $g\mid F=f$ and $g\mid S=u$.

Let now $\mathbb{K}\in\left\\{\mathbb{R},\mathbb{C}\right\\}$. Since
$\mathrm{algcl}(\mathbb{Q})\neq\mathbb{K}$, there exists an algebraic base of
$\mathbb{K}$ over $\mathbb{Q}$. Therefore we have the following.

######  Theorem 3.

Let $\mathbb{K}\in\left\\{\mathbb{R},\mathbb{C}\right\\}$. There exists a non-
identically zero derivation $d\colon\mathbb{K}\to\mathbb{K}$.

##  2\. Differential operators in polynomial rings

Let $\mathbb{K}\in\left\\{\mathbb{R},\mathbb{C}\right\\}$ and consider

$\mathscr{P}(\mathbb{K})=\left\\{p\colon\mathbb{K}\to\mathbb{K}\,|\,\text{$p$
is a polynomial}\right\\}.$

Below we present some examples that show that on $\mathscr{P}(\mathbb{K})$
there exists a large variety of mappings that satisfy the Leibniz rule.
Contrary to the case of the spaces of smooth functions that was discussed by
König and Milman, see Theorem 1 below.

###  2.1. Examples

######  Definition 2.

Let $p\in\mathscr{P}(\mathbb{K})$ be a polynomial, $x_{0}\in\mathbb{K}$ and
$k$ be a nonnegative integer. We say that the _order of zero_ of the
polynomial $p$ at $x_{0}$ is $k$, if there exists a polynomial
$q\in\mathscr{P}(\mathbb{K})$ such that $q(x_{0})\neq 0$ and

$p(x)=(x-x_{0})^{k}\cdot q(x)\qquad\left(x\in\mathbb{K}\right).$

The order of zero of a polynomial $p$ at the point $x_{0}$ will be denoted by
$n_{x_{0}}(p)$.

######  Example 3.

Let $x_{0}\in\mathbb{K}$ be arbitrarily fixed and define the mapping
$N\colon\mathscr{P}(\mathbb{K})\to\mathscr{P}(\mathbb{K})$ by

$N(p)=n_{x_{0}}(p)\cdot p\qquad\left(p\in\mathscr{P}(\mathbb{K})\right).$

Then $N\colon\mathscr{P}(\mathbb{K})\to\mathscr{P}(\mathbb{K})$ fulfills the
Leibniz rule on $\mathscr{P}(\mathbb{K})$.

######  Example 4.

The mapping $F\colon\mathscr{P}(\mathbb{K})\to\mathscr{P}(\mathbb{K})$ defined
through

$F(p)=\mathrm{deg}(p)\cdot p\qquad\left(p\in\mathscr{P}(\mathbb{K})\right),$

where $\mathrm{deg}(p)$ denotes the degree of the polynomial $p$, fulfills the
Leibniz rule on $\mathscr{P}(\mathbb{K})$.

######  Example 5.

Let us fix a polynomial $p_{0}\in\mathscr{P}(\mathbb{K})$ and consider the
mapping $P\colon\mathscr{P}(\mathbb{K})\to\mathscr{P}(\mathbb{K})$ defined by

$P(p)=p^{\prime}\cdot p_{0}\qquad(p\in\mathscr{P}(\mathbb{K})).$

Then $P$ satisfies the Leibniz rule on $\mathscr{P}(\mathbb{K})$. From this
example, it follows that an action of an operator satisfying the Leibniz rule
on a polynomial can increase its degree.

######  Example 6.

The mapping $E\colon\mathscr{P}(\mathbb{K})\to\mathscr{C}(\mathbb{K})$ defined
through

$E(p)=p\cdot\ln(|p|)\qquad\left(p\in\mathscr{P}(\mathbb{K})\right)$

satisfies the Leibniz rule. Here we adopt the convention $0\cdot\ln(0)=0$.

######  Example 7.

Let $d\colon\mathbb{K}\to\mathbb{K}$ be a derivation, then the mapping

$K(p)=p^{d}$

also fulfills the Leibniz rule on $\mathscr{P}(\mathbb{K})$.

######  Example 8.

Let $f\colon\mathbb{C}\to\mathbb{N}$ be an arbitrary function and let
$q_{0}\in\mathscr{P}(\mathbb{C})$ be a fixed polynomial. Define map
$Q\colon\mathscr{P}(\mathbb{C})\to\mathscr{P}(\mathbb{C})$ as follows. If $p$
is a constant polynomial, then put $Q(p)=0$. If $p$ is nonconstant, then it
decomposes uniquely into linear terms, i.e.,

$p(z)=a\cdot\prod_{j=1}^{N}(z-z_{j}).$

In this case, define

$Q(p)=a\cdot\sum_{k=1}^{N}q_{0}^{f(z_{k})}\cdot\prod_{j=1,j\neq
k}^{N}(z-z_{j}).$

We show that $Q$ satisfies the Leibniz rule. Fix
$p,q\in\mathscr{P}(\mathbb{C})$. If $p$ or $q$ is constant, then the assertion
is true by the definition of $Q$ (note that $Q$ is a homogeneous mapping).
Next, assume that $p$ and $q$ are nonconstant. Therefore, they decompose as

$p(z)=a_{1}\cdot\prod_{j=1}^{N}(z-z_{j}),\qquad
q(z)=a_{2}\cdot\prod_{j=N+1}^{M}(z-z_{j}).$

Therefore

$\displaystyle Q(p)\cdot q+p\cdot Q(q)$
$\displaystyle=a_{1}\cdot\left(\sum_{k=1}^{n}q_{0}^{f(z_{k})}\cdot\prod_{j=1,j\neq
k}^{N}(z-z_{j})\right)\cdot a_{2}\cdot\prod_{j=N+1}^{M}(z-z_{j})$
$\displaystyle+a_{2}\cdot\left(\sum_{k=N+1}^{M}q_{0}^{f(z_{k})}\cdot\prod_{j=N+1,j\neq
k}^{M}(z-z_{j})\right)\cdot a_{1}\cdot\prod_{j=1}^{N}(z-z_{j})$
$\displaystyle=a_{1}a_{2}\cdot\left(\sum_{k=1}^{N}q_{0}^{f(z_{k})}\cdot\prod_{j=1,j\neq
k}^{M}(z-z_{j})\right)$
$\displaystyle+a_{1}a_{2}\cdot\left(\sum_{k=N+1}^{M}q_{0}^{f(z_{k})}\cdot\prod_{j=1,j\neq
k}^{M}(z-z_{j})\right)=Q(p\cdot q).$

Since map $f$ can be arbitrary, then the degree of a polynomial $Q(p)$ can be
arbitrarily large, even if $p$ is a polynomial of degree $1$.

######  Remark.

An analogous example can be given for real polynomials, as above. Indeed, if
$p\in\mathscr{P}(\mathbb{R})$, then

$p(z)=a\prod_{i=1}^{n}(z-z_{i})\cdot\prod_{j=1}^{l}(z^{2}+\alpha_{j}z+\beta_{j})\qquad\left(z\in\mathbb{R}\right),$

where for all $j=1,\ldots,l$, the polynomials $z^{2}+\alpha_{j}z+\beta_{j}$
are irreducible over $\mathbb{R}$.

If $n=0$ and $k=0$, so if $p$ is a constant, then we simply put $Q(p)=0$.
Otherwise, let

$Q(p)=a\left(\sum_{k=1}^{n}q_{0}^{f(z_{k})}\cdot\prod_{j=1,j\neq
k}^{n}(z-z_{j})\right)\cdot\prod_{j=1}^{l}(z^{2}+\alpha_{j}z+\beta_{j}).$

A similar computation as was carried out above shows that the mapping
$Q\colon\mathscr{P}(\mathbb{R})\to\mathscr{P}(\mathbb{R})$ fulfils the Leibniz
rule. Again, since $f$ can be any mapping, the polynomial $Q(p)$ can admit
arbitrarily large degree.

######  Remark.

Observe that the Leibniz rule (1) is a homogeneous and linear equation for the
unknown operator $T\colon\mathscr{P}(\mathbb{K})\to\mathscr{C}(\mathbb{K})$.
Thus its solution space forms a linear space. Therefore, any linear
combination of the above operators $N,F,P,E,K$ and $Q$, also satisfies the
same identity.

###  2.2. Complex polynomials

In this subsection, we present our main result, which is a characterization of
all operators on $\mathscr{P}(\mathbb{C})$ that satisfy the Leibniz rule. Thus
here we consider the case $\mathbb{K}=\mathbb{C}$.

Recall, that due to the Theorem of Algebra, if $p\in\mathscr{P}(\mathbb{C})$,
then we have

$p(z)=a\cdot\prod_{j=1}^{N}(z-z_{j})$

with an appropriate nonnegative integer $N$ and with some complex numbers
$a,z_{1},\ldots,z_{N}$.

Let now $T\colon\mathscr{P}(\mathbb{C})\to\mathscr{C}(\mathbb{C})$ be an
operator that fulfills the Leibniz rule. Then, by induction, we obtain that

$T(p_{1}\cdots p_{n})=\sum_{j=1}^{n}\left(\prod_{i\neq j}p_{i}\right)\cdot
T(p_{j})\qquad\left(p_{1},\ldots,p_{n}\in\mathscr{P}(\mathbb{C})\right)$

for all positive integer $n$ and for all
$p_{1},\ldots,p_{n}\in\mathscr{P}(\mathbb{C})$.

Thus we immediately get that

$T(p)(z)=a\cdot\sum_{j=1}^{N}\left(\prod_{i\neq j}(z-z_{i})\right)\cdot
T(z-z_{j})+T(a)\prod_{j=1}^{N}(z-z_{j}).$

holds for all $p\in\mathscr{P}(\mathbb{C})$. From this, we obtain that it is
enough to determine the action of the operator $T$ on complex polynomials that
are of the form

$p(z)=az+b\qquad\left(z\in\mathbb{C}\right)$

with some appropriate complex constants $a,b$.

In the case when the domain of operator $T$ is equal to the space
$\mathscr{C}^{k}(I)$, then a localization lemma [4, Lemma 3.2] plays a key
role. It says that the value $T(p)(z)$ of operator $T$ acting on a polynomial
$p$ at a point $z$ depends only on $z,p(z),p^{\prime}(z),\dots,p^{(k)}(z)$. In
the case of the space $\mathscr{P}(\mathbb{C})$ the localization property is
not so powerful. Thus while proving our main result, which is stated below, we
need to apply a new approach. We begin with the formulation of a localization
result for the space $\mathscr{P}(\mathbb{C})$, which is a consequence of our
Theorem 4 below.

######  Lemma 2 (Localization).

Suppose that the mapping
$T\colon\mathscr{P}(\mathbb{K})\to\mathscr{C}(\mathbb{K})$ fulfills the
Leibniz rule (1) for all $p,q\in\mathscr{P}(\mathbb{K})$. Then $T$ is
_localized_ , i.e., there exists a function
$F\colon\mathbb{K}\times\mathbb{K}\to\mathbb{K}$ such that

$T(p)(z)=F(z,p(z))\qquad\left(z\in\mathbb{K}\right).$

Moreover, if the range of the operator $T$ is $\mathscr{P}(\mathbb{K})$, then
$F$ is a two-variable $\mathbb{K}$-valued _polynomial_.

Now, we will state and prove our main result without assuming the localization
property.

######  Theorem 4.

The mapping $T\colon\mathscr{P}(\mathbb{C})\to\mathscr{P}(\mathbb{C})$
fulfills the Leibniz rule (1) for all $p,q\in\mathscr{P}(\mathbb{C})$ if and
only if there exist sequences of functions $(\psi_{k})_{k\in\mathbb{N}_{0}}$
and $(\tilde{\varphi_{k}})_{k\in\mathbb{N}_{0}}$ such that for all
$k\in\mathbb{N}_{0}$ we have
$\tilde{\varphi}_{k}(ab)=a\tilde{\varphi}_{k}(b)+b\tilde{\varphi}_{k}(a)$ for
all $a,b\in\mathbb{C}$ such that

$T(p)(z)=a\cdot\sum_{j=1}^{N}\left(\prod_{i\neq
j}(z-z_{i})\right)\cdot\left(\sum_{k=0}^{n}\psi_{k}\left(-z_{j}\right)z^{k}+\sum_{k=0}^{n}\tilde{\varphi}_{k}(-z_{j})z^{k}\right)+\left(\sum_{k=0}^{n}\tilde{\varphi}_{k}(a)z^{k}\right)\prod_{j=1}^{N}(z-z_{j}),$

provided that $p(z)=a\prod_{j=1}^{N}(z-z_{j})$ and $n\in\mathbb{N}_{0}$
depends upon $p$.

###### Proof.

Firstly, we consider $\mathscr{P}(\mathbb{C})$-valued solutions of the Leibniz
rule (1). Then

(2)
$T(az+b)=\sum_{k=0}^{n}\varphi_{k}(a,b)z^{k}\qquad\left(a,b,z\in\mathbb{C}\right),$

with an appropriate $n\in\mathbb{N}$ that depends upon $a$ and $b$, and
functions $\varphi_{k}\colon\mathbb{C}^{2}\to\mathbb{C}$, $k=0,1,\ldots,n$.

First, let us consider constant polynomials. We have

$T(b)(z)=\sum_{k=0}^{n}\varphi_{k}(0,b)z^{k}$

for some $n\in\mathbb{N}$ and each $b\in\mathbb{C}$ and $z\in\mathbb{C}$. Due
to (1) we have that for all $b,c\in\mathbb{C}$ there exists some
$n\in\mathbb{N}$ that is equal to the greatest of the degree of the
polynomials $T(b),T(c),T(bc)$ and

$\displaystyle T(bc)(z)$ $\displaystyle=bT(c)(z)+cT(b)(z)$
$\displaystyle=\left(\sum_{k=0}^{n}\varphi_{k}(0,b)z^{k}\right)c+\left(\sum_{k=0}^{n}\varphi_{k}(0,c)z^{k}\right)b$
$\displaystyle=\sum_{k=0}^{n}\left(\varphi_{k}(0,b)c+\varphi_{k}(0,c)b\right)z^{k}$

On the other hand,

$T(bc)(z)=\sum_{k=0}^{n}\varphi_{k}(0,bc)z^{k}.$

Comparing the coefficients we arrive at

$\varphi_{k}(0,bc)=\varphi_{k}(0,b)c+\varphi_{k}(0,c)b$

for all $k\in\mathbb{N}$ and $b,c\in\mathbb{C}$.

Next, substitute in (2) $b=0$. We have

$T(az)=\sum_{k=0}^{n}\varphi_{k}(a,0)z^{k}.$

On the other hand, using (1) first, we get

$T(az)=aT(z)+zT(a)=a\sum_{k=0}^{n}\varphi_{k}(1,0)z^{k}+\sum_{k=0}^{n}\varphi_{k}(0,a)z^{k+1}$

for all $a\in\mathbb{C}$ and $z\in\mathbb{C}$. Comparing both polynomials we
obtain

$\varphi_{0}(a,0)=\varphi_{0}(1,0)a,$
$\varphi_{k}(a,0)=\varphi_{k}(1,0)a+\varphi_{k-1}(0,a)$

for all $k\in\mathbb{N}$ and $a\in\mathbb{C}$.

Let now $a,b\in\mathbb{C}$ be arbitrary scalars such that $a\neq 0$. By (1)
and the previous observations we have

$\displaystyle T(az+b)$
$\displaystyle=T\left(a\cdot\left(z+\frac{b}{a}\right)\right)=aT\left(z+\frac{b}{a}\right)+\left(z+\frac{b}{a}\right)T(a)$
$\displaystyle=\sum_{k=0}^{n}a\varphi_{k}\left(1,\frac{b}{a}\right)z^{k}+\sum_{k=0}^{n}\varphi_{k}\left(0,a\right)z^{k+1}+\sum_{k=0}^{n}\frac{b}{a}\varphi_{k}\left(0,a\right)z^{k}.$

Compare this equality with (2) to derive the recurrences:

$\varphi_{0}(a,b)=\varphi_{0}\left(1,\frac{b}{a}\right)a+\varphi_{0}\left(0,a\right)\frac{b}{a},$
$\varphi_{k}(a,b)=\varphi_{k}\left(1,\frac{b}{a}\right)a+\varphi_{k-1}(0,a)+\varphi_{k}\left(0,a\right)\frac{b}{a},$

For all $k=0,1,\ldots,n$, the mapping

$a\longmapsto\varphi_{k}(0,a)$

fulfills the Leibniz rule (1), we have that

$\varphi_{k}(0,b)=\varphi_{k}\left(0,a\frac{b}{a}\right)=a\varphi_{k}\left(0,\frac{b}{a}\right)+\frac{b}{a}\varphi_{k}\left(0,a\right),$

that is,

$\frac{b}{a}\varphi_{k}\left(0,a\right)=\varphi_{k}(0,b)-a\varphi_{k}\left(0,\frac{b}{a}\right)$

holds for all $a,b\in\mathbb{C}$, $a\neq 0$.

Therefore

$\displaystyle T(az+b)$
$\displaystyle=\sum_{k=0}^{n}\varphi_{k}(a,b)z^{k}=\varphi_{0}(a,b)+\sum_{k=1}^{n}\varphi_{k}(a,b)z^{k}$
$\displaystyle=\varphi_{0}\left(1,\frac{b}{a}\right)a+\varphi_{0}\left(0,a\right)\frac{b}{a}+\sum_{k=1}^{n}\left(\varphi_{k}\left(1,\frac{b}{a}\right)a+\varphi_{k-1}(0,a)+\varphi_{k}\left(0,a\right)\frac{b}{a}\right)z^{k}$
$\displaystyle=\varphi_{0}\left(1,\frac{b}{a}\right)a+\varphi_{0}(0,b)-a\varphi_{0}\left(0,\frac{b}{a}\right)$
$\displaystyle+\sum_{k=1}^{n}\left(\varphi_{k}\left(1,\frac{b}{a}\right)a+\varphi_{k-1}(0,a)+\varphi_{k}(0,b)-a\varphi_{k}\left(0,\frac{b}{a}\right)\right)z^{k}$
$\displaystyle=\psi_{0}\left(\frac{b}{a}\right)a+\varphi_{0}(0,b)+\sum_{k=1}^{n}\left(\psi_{k}\left(\frac{b}{a}\right)a+\varphi_{k-1}(0,a)+\varphi_{k}(0,b)\right)z^{k}$
$\displaystyle=\sum_{k=0}^{n}a\psi_{k}\left(\frac{b}{a}\right)z^{k}+\sum_{k=0}^{n}\varphi_{k}(0,b)z^{k}+\sum_{k=1}^{n}\varphi_{k-1}(0,a)z^{k}$

for all $a,b\in\mathbb{C}$ with $a\neq 0$, where for all $k=0,1,\ldots,n$, the
functions $\psi_{k}\colon\mathbb{C}\to\mathbb{C}$ are defined through

$\psi_{k}(c)=\varphi_{k}(1,c)-\varphi_{k}(0,c)\qquad\left(c\in\mathbb{C}\right).$

Define the functions $\tilde{\varphi_{k}}\colon\mathbb{C}\to\mathbb{C}$
through

$\tilde{\varphi_{k}}(a)=\varphi_{k}(0,a)\qquad\left(x\in\mathbb{C}\right).$

Observe that for all $k=0,\ldots,n$, the mapping $\tilde{\varphi_{k}}$
satisfies the Leibniz rule (1). So especially $\tilde{\varphi_{k}}(1)=0$. With
these functions, we can deduce that

$T(az+b)=\sum_{k=0}^{n}a\psi_{k}\left(\frac{b}{a}\right)z^{k}+\sum_{k=0}^{n}\tilde{\varphi}_{k}(b)z^{k}+\sum_{k=1}^{n}\tilde{\varphi_{k-1}}(a)z^{k}$

for all $a,b\in\mathbb{C}$ and $z\in\mathbb{C}$.

Therefore for any fixed $z_{j}\in\mathbb{C}$, we have

$\displaystyle T(z-z_{j})$
$\displaystyle=\sum_{k=0}^{n}\psi_{k}\left(-z_{j}\right)z^{k}+\sum_{k=0}^{n}\tilde{\varphi}_{k}(-z_{j})z^{k}+\sum_{k=1}^{n}\tilde{\varphi}_{k-1}(1)z^{k}$
$\displaystyle=\sum_{k=0}^{n}\psi_{k}\left(-z_{j}\right)z^{k}+\sum_{k=0}^{n}\tilde{\varphi}_{k}(-z_{j})z^{k}$

and also for any fixed $a\in\mathbb{C}$,

$T(a)=\sum_{k=0}^{n}\tilde{\varphi}_{k}(a)z^{k}.$

Together with (1), this implies that if $p\in\mathscr{P}(\mathbb{\mathbb{C}})$
with

$p(z)=a\prod_{j=1}^{N}(z-z_{j}),$

then

$\displaystyle T(p(z))$ $\displaystyle=a\cdot\sum_{j=1}^{N}\left(\prod_{i\neq
j}(z-z_{i})\right)\cdot T(z-z_{j})+T(a)\prod_{j=1}^{N}(z-z_{j})$
$\displaystyle=a\cdot\sum_{j=1}^{N}\left(\prod_{i\neq
j}(z-z_{i})\right)\cdot\left(\sum_{k=0}^{n}\psi_{k}\left(-z_{j}\right)z^{k}+\sum_{k=0}^{n}\tilde{\varphi}_{k}(-z_{j})z^{k}\right)$
$\displaystyle+\left(\sum_{k=0}^{n}\tilde{\varphi}_{k}(a)z^{k}\right)\prod_{j=1}^{N}(z-z_{j})$

Conversely, let $p,q\in\mathscr{P}(\mathbb{C})$ with

$p(z)=a_{1}\prod_{j=1}^{N}(z-z_{j})\quad\text{and}\quad
q(z)=a_{2}\prod_{j=N+1}^{M}(z-z_{j}).$

Then

$\displaystyle T(p)q$
$\displaystyle+T(q)p=\left(a_{1}\cdot\sum_{j=1}^{N}\left(\prod_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{N}(z-z_{i})\right)\cdot
T(z-z_{j})+T(a_{1})\prod_{j=1}^{N}(z-z_{j})\right)\cdot
a_{2}\prod_{j=N+1}^{M}(z-z_{j})$
$\displaystyle+\left(a_{2}\cdot\sum_{j=N+1}^{M}\left(\prod_{\begin{subarray}{c}i=N+1\\\
i\neq j\end{subarray}}^{M}(z-z_{i})\right)\cdot
T(z-z_{j})+T(a_{2})\prod_{j=N+1}^{M}(z-z_{j})\right)\cdot
a_{1}\prod_{j=1}^{N}(z-z_{j})$
$\displaystyle=a_{1}a_{2}\cdot\sum_{j=1}^{M}\left(\prod_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{M}(z-z_{i})\right)\cdot
T(z-z_{j})+(T(a_{1})a_{2}+T(a_{2})a_{1})\prod_{j=1}^{M}(z-z_{j})$
$\displaystyle=a_{1}a_{2}\cdot\sum_{j=1}^{M}\left(\prod_{\begin{subarray}{c}i=1\\\
i\neq j\end{subarray}}^{M}(z-z_{i})\right)\cdot
T(z-z_{j})+T(a_{1}a_{2})\prod_{j=1}^{M}(z-z_{j})$ $\displaystyle=T(pq).$

∎

From the above theorem, we can derive a few corollaries for some special types
of operators. We say that an operator
$T\colon\mathscr{P}(\mathbb{C})\to\mathscr{P}(\mathbb{C})$ _decreases the
degree_ if for any polynomial $p$ of degree $N$, the degree of $T(p)$ is
smaller than $N$. Moreover, operator $T$ does not increase the degree if, for
any polynomial $p$ of degree $N$, the degree of $T(p)$ is not greater than
$N$.

######  Corollary 1.

The mapping $T\colon\mathscr{P}(\mathbb{C})\to\mathscr{P}(\mathbb{C})$
fulfills the Leibniz rule (1) for all $p,q\in\mathscr{P}(\mathbb{C})$ and
decreases the degree if and only if there exist two functions $\psi_{0}$ and
$\tilde{\varphi_{0}}$ such that we have
$\tilde{\varphi}_{0}(ab)=a\tilde{\varphi}_{0}(b)+b\tilde{\varphi}_{0}(a)$ for
all $a,b\in\mathbb{C}$ such that

$T(p)(z)=a\cdot\sum_{j=1}^{N}\left(\prod_{i\neq
j}(z-z_{i})\right)\cdot\left(\psi_{0}\left(-z_{j}\right)+\tilde{\varphi}_{0}(-z_{j})\right).$

provided that $p(z)=a\prod_{j=1}^{N}(z-z_{j})$.

Consequently, for each monomial $p(z)=z^{N}$ there exists a constant
$c_{p}\in\mathbb{C}$ such that

$T(p)=c_{p}\cdot p^{\prime}.$

###### Proof.

Since $T$ decreases the degree, then necessarily $T(b)=0$ for any constant
polynomial $b\in\mathscr{P}(\mathbb{C})$. Therefore,

$\sum_{k=0}^{n}\varphi_{k}(0,b)z^{k}=0.$

Moreover, since the image of each polynomial of the degree $1$ is a constant
map, we have that $n=0$ in the statement of Theorem 4. ∎

######  Corollary 2.

The mapping $T\colon\mathscr{P}(\mathbb{C})\to\mathscr{P}(\mathbb{C})$
fulfills the Leibniz rule (1) for all $p,q\in\mathscr{P}(\mathbb{C})$ and not
increases the degree if and only if there exist four functions $\psi_{0}$,
$\psi_{1}$, $\tilde{\varphi_{0}}$ and $\tilde{\varphi_{1}}$ such that we have
$\tilde{\varphi}_{k}(ab)=a\tilde{\varphi}_{k}(b)+b\tilde{\varphi}_{k}(a)$ for
all $a,b\in\mathbb{C}$ and $k\in\\{0,1\\}$ such that

$T(p)(z)=a\cdot\sum_{j=1}^{N}\left(\prod_{i\neq
j}(z-z_{i})\right)\cdot\left(\sum_{k=0}^{1}\psi_{k}\left(-z_{j}\right)z^{k}+\sum_{k=0}^{1}\tilde{\varphi}_{k}(-z_{j})z^{k}\right)+\tilde{\varphi}_{0}(a)\prod_{j=1}^{N}(z-z_{j}),$

provided that $p(z)=a\prod_{j=1}^{N}(z-z_{j})$.

Consequently, for each monomial $p(z)=z^{N}$ there exist constants
$c_{p},d_{p}\in\mathbb{C}$ such that

$T(p)=c_{p}\cdot p^{\prime}+d_{p}\cdot p.$

###### Proof.

Since $T$ does not increase the degree, then we have that $n=1$ in the
statement of Theorem 4. Using the fact that Leibniz mappings vanish at $0$ and
$1$ we derive the latter part. ∎

######  Acknowledgment.

The research of Eszter Gselmann has been supported by project no. K134191 that
has been implemented with the support provided by the National Research,
Development and Innovation Fund of Hungary, financed under the K_20 funding
scheme.

The research has been accomplished during the stay of Włodzimierz Fechner at
the University of Debrecen, covered by the Visegrad Fellowship #62320104.

## References

  * [1] Shiri Artstein-Avidan, Dmitry Faifman, and Vitali Milman.  On multiplicative maps of continuous and smooth functions.  In Geometric aspects of functional analysis, volume 2050 of Lecture Notes in Math., pages 35–59. Springer, Heidelberg, 2012.
  * [2] Ola Bratteli, George A. Elliott, and Derek W. Robinson.  The characterization of differential operators by locality: classical flows.  Compositio Math., 58(3):279–319, 1986.
  * [3] Martin Engert.  A characterization of differential operators.  Proc. Amer. Math. Soc., 18:87–93, 1967.
  * [4] Hermann König and Vitali Milman.  Operator relations characterizing derivatives.  Birkhäuser / Springer, Cham, 2018.
  * [5] Marek Kuczma.  An introduction to the theory of functional equations and inequalities.  Birkhäuser Verlag, Basel, second edition, 2009.  Cauchy’s equation and Jensen’s inequality, Edited and with a preface by Attila Gilányi.
  * [6] Jaak Peetre.  Une caractérisation abstraite des opérateurs différentiels.  Math. Scand., 7:211–218, 1959.


# Interaction-induced velocity renormalization in magic angle twisted trilayer
graphene

Laura Classen  Condensed Matter Physics and Materials Science Division,
Brookhaven National Laboratory, Upton, New York 11973, USA  Max Planck
Institute for Solid State Research, D-70569 Stuttgart, Germany  J. H. Pixley
Department of Physics and Astronomy, Center for Materials Theory, Rutgers
University, Piscataway, New Jersey 08854, USA  Elio J. König  Max Planck
Institute for Solid State Research, D-70569 Stuttgart, Germany

###### Abstract

Twistronics heterostructures provide a novel route to control the electronic
single particle velocity and thereby to engineer strong effective
interactions. Here we show that the reverse may also hold, i.e. that these
interactions strongly renormalize the band structure. We demonstrate this
mechanism for mirror-symmetric magic angle twisted trilayer graphene at charge
neutrality and in the vicinity of a phase transition which can be described by
an Ising Gross-Neveu critical point corresponding, e.g., to the onset of
valley Hall or Hall order. While the non-interacting model displays massless
Dirac excitations with strongly different velocities, we show that interaction
corrections make them equal in the infrared. However, the RG flow of the
velocities and of the coupling to the critical bosonic mode is strongly non-
monotonic and dominated by the vicinity of a repulsive fixed point. We predict
experimental consequences of this theory for tunneling and transport
experiments and discuss the expected behavior at other quantum critical
points, including those corresponding to intervalley coherent ordering.

Figure 1: a) Schematics of symmetric TTG. b) Mean field phase diagram of the
GN model, Eq. (1), as a function of twist angle $\theta$ for coupling strength
$g_{+}=g_{-}=\bar{g}$. Insets: Schematic dispersion relations in the
semimetallic phase (blue), and in the gapped phase (correlated insulator,
green). At criticality (red line), the dispersion is scale dependent and slow
and fast fermions acquire the same, slow velocity in the infrared.

Introduction: Manipulating the nature of low-energy electronic excitations in
solid-state compounds has remained a long standing goal in condensed matter
physics. Thanks to the recent advances in the isolation and manipulation of
two-dimensional materials, it is now possible to dramatically modify the low
energy band structure by twisting (i.e. rotating) and stacking. As first
demonstrated in twisted bilayer graphene (TBG), the moiré pattern formed from
the interference between the twisted and stacked layers greatly renormalizes
the excitation velocitySuárez Morell  _et al._ (2010); Bistritzer and
MacDonald (2011) allowing a small bare interaction scale to become dominant,
which drives the formation of correlated insulating states and
superconductivityCao  _et al._ (2018a, b). This ‘twistronic’ paradigm has now
been extended well beyond graphene to, e.g., emulating Hubbard models using
transition metal dichalgonidesWang  _et al._ (2020); Tang  _et al._ (2020); Wu
_et al._ (2018); Zhou  _et al._ (2021); Scherer  _et al._ (2021); Schrade and
Fu (2019) or manipulating the superconducting state using cuprates Zhao  _et
al._ (2021); Volkov  _et al._ (2021); Tummuru  _et al._ (2021).

Beyond stacking and twisting two sheets of atoms, recent experiments have
shown it is possible to twist various arrangements of a number of graphene
layers Cao  _et al._ (2020); Shen  _et al._ (2020); Liu  _et al._ (2020); Park
_et al._ (2021a); Zhang  _et al._ (2021) and also produce non-trivial phase
diagrams with insulating and superconducting states. In particular,
symmetrically twisted multi-layer graphene systems, where the layers are
rotated by an alternating angle $\pm\theta$, can be considered as a sum of
twisted bilayers.Khalaf  _et al._ (2019) Thus, they offer an alternative
setting where strongly renormalized narrow bands appear at a magic angle
(where the Dirac velocity vanishes) as in TBG. For an odd number of layers,
these TBG-like bands are supplemented by fast Dirac dispersing bands due to
the symmetric nature of the system.Khalaf  _et al._ (2019); Călugăru  _et al._
(2021) For example, in symmetrically twisted trilayer graphene (TTG) the
mirror eigenvalues in the setup of Fig. 1 a) protect the fast excitations from
being renormalized by the twist. As a result, fast and slow Dirac excitations
coexist in TTG near charge neutrality.

Experiments on TTG have discovered a readily tunable superconducting ground
state, phases of broken symmetry at fillings of $\pm 2$ electrons per moiré
unit cell, and a sizeable resistance at charge neutrality.Park  _et al._
(2021b); Hao  _et al._ (2021) The spontaneous appearance of a gap at charge
neutrality was not reported in TTG systems to date. But it was observed in
some TBG devices without substrate alignmentLu  _et al._ (2019); Wu  _et al._
(2021), where it is believed to arise from interaction effects and expected
based on the exact solution of effective theoriesLian  _et al._ (2021). The
difference between insulating and semi-metallic behavior at charge neutrality
may come from a strain-induced quantum phase transition. Parker  _et al._
(2021) Theoretically, different candidate states for a gap opening at charge
neutrality in TBG have similar energies. Bultinck  _et al._ (2020); Brillaux
_et al._ (2020); Da Liao  _et al._ (2021); Kang and Vafek (2020); Ochi  _et
al._ (2018); Zhang  _et al._ (2020); Liao  _et al._ (2021) Among them,
(sublattice-polarized) quantum valley Hall or intervalley coherent order have
also been discussed to gap out the slow bands in TTG, based on extensions of
the effective, exactly soluble models from TBG to TTGChristos  _et al._
(2021).

When Dirac electrons reside close to a phase transition and interact with
critical bosonic modes, Lorentz symmetry is often emergentRoy  _et al._
(2016). This implies that different bare velocities must renormalize to become
equal, and the average velocity may vanish Sitte  _et al._ (2009), diverge
König and Pixley (2018), or approach a finite value, as it occurs, for
example, to the two Dirac velocities in spin 3/2 systems. Roy  _et al._
(2018); Boettcher (2020); Mandal (2021) While in TTG devices, the vicinity to
external gates implies effective short-range interactions, we highlight that
velocity renormalization due to Coulomb interactions in suspended graphene was
predicted theoreticallySheehy and Schmalian (2007); Barlas  _et al._ (2007);
Hwang  _et al._ (2007) and subsequently confirmed experimentallyElias  _et
al._ (2011); Yu  _et al._ (2013); Ulybyshev  _et al._ (2021) more than a
decade ago. In contrast to the single-particle magic-angle phenomenon, many-
body effects are crucial to the physics in the examples above. This raises the
question how interactions affect the velocity renormalization in TTG. Do the
velocities of fast and slow Dirac modes approach each other and does this work
for or against the correlation effects in the TBG-like bands which are rooted
in the vanishing Dirac velocity?

In this paper we address the velocity renormalization combining the single-
particle effects in a moiré system with the many-body physics of electron-
electron interactions. To this end, we study an Ising Gross-Neveu (GN) theory
with a generic ratio of fast and slow Dirac fermions employing a
renormalization group (RG) procedure in $D=2+1$ dimensions. This describes the
low energy properties of TTG near charge neutrality in the vicinity of a
quantum critical point, which separates a Dirac semimetal phase from a gapped
(valley) quantum Hall state. The critical interaction strength can be smoothly
varied by changing the twist, which allows us to connect the untwisted
critical point to the ground state at the magic-angleBrillaux  _et al._
(2020), Fig. 1 b). We demonstrate that slow and fast Dirac excitations cannot
be considered separately as often done in TTG. As a result of their strong
influence on each other the slow excitations become faster and the fast ones
slower ensuring emergent Lorentz symmetry at the asymptotic critical point.
Overall, the system becomes more correlated in the sense that the mean
velocity decreases. More generally, the asymptotic mean velocity is determined
by the ratio of the number of fast and slow modes. We estimate the velocity
renormalization to be a measureable effect on length scales accessible in
experiment and predict the simultaneous suppression of the quasiparticle
weight.

Model. As an effective model for TTG near charge neutrality, we consider the
generalized Gross-Neveu theory ($\hbar=1$ throughout and $\sigma_{x,y,z}$ are
Pauli matrices in sublattice space)

$\mathcal{L}=\sum_{\pm}\sum_{\alpha=1}^{N_{\pm}}\bar{\psi}_{\alpha,\pm}[\partial_{\tau}+v_{\pm}\vec{p}\cdot\vec{\sigma}+g_{\pm}\phi\sigma_{z}]\psi_{\alpha,\pm}+\frac{N}{2}\phi^{2},$
(1)

which describes fast (+) and slow (-) Dirac fermions with $2N=(N_{+}+N_{-})$
flavors in total and different velocities $v_{\pm}$. To leading order, they
are related to the twist angle via

$\frac{v_{-}}{v_{+}}=\frac{1-3\alpha^{2}}{1+6\alpha^{2}},$ (2)

where $\alpha=\theta_{c}/\theta{\sqrt{3}}$.Bistritzer and MacDonald (2011);
Khalaf  _et al._ (2019) The magic angle defined by vanishing velocity occurs
at $\theta=\theta_{c}\approx 1.6^{\circ}$. Interactions are mediated by the
Ising order parameter field $\phi$, which is introduced via a Hubbard-
Stratonovich transformation and corresponds to the interaction channel for a
(sublattice polarized) quantum Hall or valley quantum Hall instability. By
symmetry, $\phi$ couples simultaneously to both slow and fast fields, see
Supplement Sup . The case $N_{+}=4$ (spin + valley) and $N_{-}=8$ (spin +
valley + mini-valley) applies to TTG in the vicinity of the magic angle. More
generally, the model captures the dynamical mass generation in Dirac materials
with an arbitrary number of slow and fast flavors near the onset of
spontaneous order that breaks an Ising symmetry.

Mean field phase diagram. In the large $N$ limit, we can integrate out the
fermions and calculate the effective bosonic mass in the saddle-point
approximation

$m=\frac{1}{g_{+}g_{-}}-\sum_{\pm}\frac{N_{\pm}}{2N\pi}\frac{g_{\pm}}{v_{\pm}g_{\mp}}\Lambda,$
(3)

where we use $\Lambda$ as a UV momentum cut-off scale. The second term in Eq.
(3) stems from the fermionic susceptibility, Fig. 2 a). A sign change in $m$
signals the dynamical mass generation and defines the critical coupling where
the mean field phase transition occurs. Exploiting the relation in Eq. (2)
allows us to perturbatively determine the interaction driven transition as a
function of twist angle and leads to the phase diagram plotted in Fig. 1 b).
Importantly, this shows that at the magic-angle an infinitesimal interaction
strength is sufficient to drive the mean field transition, a hallmark of the
renormalized flat band that occurs in the slow Dirac excitations. Right at the
magic angle of TTG, the kinetic part of the slow sector is dominated by a
quadratic band touching Hejazi  _et al._ (2019) where the critical coupling
vanishes. Therefore, logarithmic self-energy corrections introduce an
additional sharpeningRay  _et al._ (2018) of the phase transition near
$\theta=\theta_{c}$ (not shown).

Quantum fluctuations. We here summarize the strategy to incorporate quantum
fluctuations near $m=0$, in order to correctly capture the nature of the
quantum critical point (leaving technical details for the supplementSup ).
This allows us to unveil an interaction driven renormalization of the
velocities making them equal in the infared. The contributions beyond mean
field theory are organized in a systematic way via a large-$N$ expansion,
Rosenstein  _et al._ (1989, 1991); Khveshchenko and Paaske (2001); Gracey
(2018); Boyack  _et al._ (2019) which has the advantage of performing
calculations directly in $D=2+1$ and complements commonly used expansions near
the upper critical dimension.Roy  _et al._ (2018); Boettcher (2020); Mandal
(2021) As a starting point, it is exploited that the one-loop bosonic self
energy (i.e. the susceptibility, Fig. 2 a)) overpowers a quadratic kinetic
term leading to an effective RPA resummed bosonic propagator to leading order
in $N$

$D(\omega,\vec{q})=\frac{1}{N}\left(m+\sum_{\pm}\alpha_{\pm}\sqrt{\omega^{2}+v_{\pm}^{2}\vec{q}^{2}}\right)^{-1},$
(4)

where $\alpha_{\pm}=N_{\pm}\rho_{g}^{\pm 1}\rho_{v}^{\mp 1}/N$ and we defined
$\rho_{g}=g_{+}/g_{-}$ and $\rho_{v}=v_{+}/v_{-}$. The linear scaling of
$D^{-1}$ leads to logarithmic divergencies in the next-to-leading-order
diagrams in Fig. 2 b) - e), which are formally cured by the inclusion of
counter terms in the bare action. These counter terms renormalize all coupling
constants as well as the fields in Eq. (1) and determine the corresponding
renormalization group (RG) equations. In addition to this field-theoretic RG
procedure, we also performed Wilsonian momentum-shell-RG, yielding the same
result.Sup

Figure 2: a) The polarization operator calculated with the non-interacting
fermionic Green’s function. It is $\mathcal{O}(N)$, while other diagrams are
sub-leading at $N\rightarrow\infty$. This controls an RPA resummation, which
is denoted by a double dashed line. b) Fermionic self-energy and vertex
correction [$\mathcal{O}(1/N)$]. d) Leading two-loop corrections to the
bosonic self-energy [$\mathcal{O}(1)$].

Renormalization group equations. There are four running coupling constants,
$v_{\pm},g_{\pm}$, in addition to field renormalizations
$Z_{\psi}^{\pm},Z_{\phi}$, which introduce non-trivial scaling dimensions for
three fields. It is convenient to express the RG equations in terms of mean
and relative coupling constants
$\rho_{g}=g_{+}/g_{-},\rho_{v}=v_{+}/v_{-},\bar{v}=\sqrt{v_{+}v_{-}}$ and the
bosonic gap $m$. We find to leading order in the $1/N$ expansion Sup

$\displaystyle\frac{d\rho_{g}}{d\ln b}$
$\displaystyle=\rho_{g}\sum_{\pm}\frac{\pm
4}{3\pi^{2}N}\left[f_{2}^{\pm}-3f_{1}^{\pm}\right],$ (5a)
$\displaystyle\frac{d\rho_{v}}{d\ln b}$
$\displaystyle=\rho_{v}\sum_{\pm}\frac{\pm
2}{\pi^{2}N}\left[f_{2}^{\pm}-f_{1}^{\pm}\right],$ (5b)
$\displaystyle\frac{d\bar{v}}{d\ln b}$
$\displaystyle=\bar{v}\sum_{\pm}\frac{1}{\pi^{2}N}\left[f_{2}^{\pm}-f_{1}^{\pm}\right],$
(5c) $\displaystyle\frac{dm}{d\ln b}$
$\displaystyle=\left\\{1-\frac{8}{3\pi^{2}N}\left[3f_{3}-\sum_{\pm}\frac{3f_{1}^{\pm}-f_{2}^{\pm}}{2}\right]\right\\}m,$
(5d) with $b=\Lambda/\Lambda^{\prime}$ where $\Lambda^{\prime}$ is the running
scale. We have also absorbed a shift into the mass in Eq. (5d) as compared to
Eq. (3). In addition, bosonic and fermionic fields are rescaled at one-loop
level $\displaystyle\frac{d\ln Z_{\psi}^{\pm}}{d\ln b}$
$\displaystyle=-\frac{2}{3\pi^{2}N}[3f_{1}^{\pm}-2f_{2}^{\pm}],$ (5e)
$\displaystyle\frac{d\ln Z_{\phi}}{d\ln b}$
$\displaystyle=\frac{4}{3\pi^{2}N}\sum_{\pm}[3f_{1}^{\pm}-f_{2}^{\pm}].$ (5f)

In these equations, we have introduced the five functions
$f_{1,2}^{\pm}=f_{1,2}^{\pm}(\rho_{g},\rho_{v},N_{+}/N_{-}),f_{3}=f_{3}(\rho_{g},\rho_{v},N_{+}/N_{-})$,
which are implicitly defined through the following integrals

$\displaystyle f_{1}^{\pm}$ $\displaystyle=\rho_{g}^{\pm
1}\int_{-\infty}^{\infty}dx\frac{1}{({x^{2}+\rho_{v}^{\pm
1}})[\sum_{\xi}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}]},$ (6a) $\displaystyle
f_{2}^{\pm}$ $\displaystyle=3\rho_{g}^{\pm
1}\int_{-\infty}^{\infty}dx\frac{\sum_{\xi=\pm}\alpha_{\xi}\frac{\rho_{v}^{\xi}/2}{\sqrt{x^{2}+\rho_{v}^{\xi}}}}{({x^{2}+\rho_{v}^{\pm
1}})[\sum_{\xi=\pm}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}]^{2}},$ (6b)
$\displaystyle f_{3}$
$\displaystyle=\sum_{\pm}\int_{-\infty}^{\infty}dx\frac{\alpha_{\pm}}{\sqrt{x^{2}+\rho_{v}^{\pm
1}}[\sum_{\xi=\pm}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}]^{2}}.$ (6c)

The non-trivial coefficients of the RG equations are due to the presence of an
unequal number of Dirac excitations that have distinct velocities and coupling
constants (i.e. $\rho_{v}\neq 1$ and $\rho_{g}\neq 1$). Note that in the
isotropic limit $\rho_{v}=\rho_{g}=1$ we find that all these functions reduce
to $f_{1,2}^{\pm}=f_{3}=1$ and we reproduce results from the
literature.Rosenstein  _et al._ (1991); Khveshchenko and Paaske (2001); Moshe
and Zinn-Justin (2003)

Figure 3: a) RG Flow, Eqs. (5) in the plane of ratios of velocities
$\rho_{v}=v_{+}/v_{-}$ and couplings $\rho_{g}=g_{+}/g_{-}$, The color coding
represents the beta function of the average velocity
$\bar{v}=\sqrt{v_{+}v_{-}}$. The apparent asymmetry under
$\rho_{v,g}\leftrightarrow 1/\rho_{v,g}$ vanishesSup  in the limit
$N_{-}=N_{+}$. b) Integration of RG equations Eq. (5a)-(5c) for starting
values $\rho_{v}=50,\rho_{g}=5,\bar{v}=1$. c) Field renormalizations, Eq.
(5e)-(5f). In panels b), c) we indicate the first (second) stage of RG by
yellow (blue) background color. In all panels we set $N_{-}=2N_{+}=8$, as in
TTG.

Analysis of RG flow. The renormalization group equations of the ratio of
velocities and coupling constant are closed, i.e. the right hand side of Eqs.
(5a),(5b) only depends on the two running coupling constants
$\rho_{g},\rho_{v}$ themselves. We solve these equations numerically and
present the solution in Fig. 3 a). We find an attractive fixed point
$\rho_{v}=\rho_{g}=1$, where Lorentz symmetry is emergent and two lines of
repulsive fixed points where only one of the two fermion species couples to
the critical boson, $\rho_{g}\in\\{0,\infty\\}$. In particular, the two fixed
points at $\rho_{v}=0=\rho_{g}$ and $\rho_{v}=\infty=\rho_{g}$ are connected
to the attractive fixed point by a separatrix. Generally, during the evolution
the two coupling constants, $\rho_{v,g}$ approach this separatrix in a first,
quick RG stage, and then slowly approach isotropy. This shows that,
independent on the difference between bare velocities and coupling constants,
the fast and slow Dirac fermions cannot be considered separately once
correlation effects are taken into account.

To describe the RG flow for TTG, we choose $N_{-}=2N_{+}$ and starting values
$\rho_{v}\gg 1$, $\rho_{g}\sim 1$ on the right in Fig. 3 a). We Taylor expand
the RG equations about this limit using numerically evaluated series
coefficients, and obtain

$\displaystyle\frac{d\ln\rho_{g}}{d\ln b}$
$\displaystyle\approx\frac{4}{N\pi^{2}},\qquad\frac{d\ln m}{d\ln b}\approx
1-\frac{8}{N\pi^{2}}-12\frac{\rho_{g}-1}{N\pi^{2}},$ $\displaystyle\frac{d\ln
Z_{\psi}^{-}}{d\ln b}$
$\displaystyle\approx-\frac{1}{N\pi^{2}},\quad\frac{d\ln Z_{\phi}}{d\ln
b}\approx-\frac{4}{N\pi^{2}},$ (7)

while

$\displaystyle\frac{d\ln\rho_{v}}{d\ln b}=\frac{d\ln\bar{v}}{d\ln
b}=\frac{d\ln Z_{\psi}^{+}}{d\ln b}=\mathcal{O}(1/\rho_{v}).$ (8)

Thus $\rho_{v}$ and $\bar{v}$ are only weakly affected in comparison to
$\rho_{g}$, because the RG leads to a massive enhancement of $\rho_{g}$.
However, we can still observe that $v_{-}$ gets enhanced and $v_{+}$ reduced
(Fig. 3 b). We also see that all field renormalizations start to decrease so
that bosonic and fermionic excitations become less coherent, Fig. 3 c). This
first stage of RG, illustrated by a yellow background in Fig. 3 b) and c),
lasts until $\rho_{g}(b)\sim\rho_{v}(0)$ corresponding to length scales below
$L_{\rm crossover}\sim[\rho_{v}(0)/\rho_{g}(0)]^{N\pi^{2}/4}/\Lambda$, which
can be very large.

In the second stage of RG, depicted by a blue background in Fig. 3 b) and c),
both $\rho_{v,g}$ very slowly approach unity and there is a moderate
renormalization of $\bar{v}$. Of course, if starting values of
$\rho_{g}(0)\simeq\rho_{v}(0)$, the system may only explore the second regime
in which the coupling constants slowly approach $\rho_{v}=\rho_{g}=1$. In the
vicinity of the attractive fixed point, the flow along the separatrix, which
is asymptotically defined by $\rho_{g}=(3\rho_{v}-1)/2$, is characterized by a
scaling dimension $-{8}/({15\pi^{2}N})\stackrel{{\scriptstyle
N=6}}{{\approx}}-0.009$. corresponding to the least irrelevant operator. At
the attractive fixed point, the average velocity $\bar{v}$ approaches a finite
(non-universal) value, which is determined by the relative difference between
fast and slow flavors $(N_{+}-N_{-})/(N_{+}+N_{-})$Sup . Thus for TTG, where
$N_{-}>N_{+}$, $\bar{v}$ decreases.

Observables and consequences for experiments. Observables which allow to
experimentally access the velocity renormalization include quantum
oscillations, as exploited to uncover the velocity renormalization in
suspended grapheneElias  _et al._ (2011). Another probe is the tunneling
density of states that within our RG calculation is given by

$\nu(E)=\sum_{\pm}N_{\pm}\frac{Z_{\psi}^{\pm}(b_{E})}{2\pi
v_{\pm}^{2}(b_{E})}|E|.$ (9)

It not only measures the velocity renormalization, but also the suppression of
the quasi-particle weight at the corresponding energy scale
$b_{E}=\bar{v}\Lambda/E$, see Fig. 3 c.

Based on our estimate for the extent of the first stage of the RG flow, we
expect that experiments mainly explore this first regime. For example, for
transport experiments, the maximum RG time can be estimated as $\ln(b_{\rm
final})-\ln(b_{\rm initial})\sim\ln(500K/10mK)\sim 10$, where $50meV\sim 500K$
is the mini bandwidth of the fast bands and $10mK$ the base temperature of
dilution refrigerators.

At the same time, the bare coupling constants are a priori unknown and
variations in twist angle affect the bare $\rho_{v}$ in TTG. In addition,
other Dirac materials also have different starting values. In these situations
other segments of the RG flow can be probed. This is the rationale behind the
above complete theoretical analysis of the RG flow. Assuming starting values
as in Fig. 3 b), we highlight that in the initial, in TTG experimentally
accessible regime $\ln(b)\in(0,10)$, $\rho_{g}$ doubles and $\rho_{v}$
increases by 20 % - both are certainly measurable effects.

Another implication of our analysis regards the quantum-critical regime as the
system is tuned across the quantum phase transition, where experiments can
probe quantum critical fluctuations akin to what has been observed in strongly
correlated materials. To this end, due to the parametrically slow RG flow, we
define ’local’ scaling exponents for the correlation length and anomalous
dimensions via $1/\tilde{\nu}=d\ln m/d\ln b$, $\tilde{\eta}_{\psi\pm}=-d\ln
Z_{\psi}^{\pm}/2d\ln b$ and $\tilde{\eta}_{\phi}=1-d\ln Z_{\phi}/d\ln b$ using
Eqs. (5d)-(5f). These are not the standard universal critical exponents but
scale-dependent quantities, unless the RG flow reaches one of the fixed
points. However, they establish an analytical relationship between the
velocity renormalization and experimentally observable exponents. At the
attractive Lorentz-symmetric fixed point, they reduce to the universal scaling
exponents of the Ising GN modelRosenstein  _et al._ (1991); Khveshchenko and
Paaske (2001): $\nu=1+8/(3\pi^{2}N)$,
$\eta_{\psi+}=\eta_{\psi-}=1/(3\pi^{2}N)$, $\eta_{\phi}=1-16/(3\pi^{2}N)$.

Effect of displacement field and Dirac node offset. A major experimental
tuning knob in twisted trilayer graphene is the displacement field $D_{0}$.
The displacement field generates a hybridization of order $D_{0}$ of fast and
slow modes located at the mini-$K$ point. This hybridization vertically splits
the two Dirac points to energies $\pm D_{0}$ and generates a Fermi surface of
size $p_{F}=D_{0}/\sqrt{v_{+}v_{-}}$ at charge neutrality. The renormalization
group flow calculated above describes the system at energies larger than
$D_{0}$. The evaluation of running coupling constants at the scale $D_{0}$
serves as bare values for the physics at lower energies, which is not of
Gross-Neveu type (instead, e.g., logarithmic Fermi surface instabilities may
occur in Cooper and density wave channels).

We also comment on a possible offset in energy between the Dirac node of fast
and slow fermions. Such an offset would also provide an infrared cut-off for
the integral of RG equations, beyond which the impact of the Fermi surface in
one of the two carrier types becomes relevant. Ab initioPark  _et al._ (2021b)
calculations for TTG predict an offset of the order of meV which is thus small
as compared to the mini-band width.

Conclusions: Motivated by the recent discoveries in TTG, we have presented the
leading renormalization group equations for an Ising Gross-Neveu theory in
$D=2+1$ dimensions where the order parameter field couples to $N_{+}$ fast and
$N_{-}$ slow fermion fields. We found a two-stage RG flow, with an ultimate
Lorentz-symmetric fixed point at which all fermion velocities are equal, but
the mean velocity is reduced, thereby enhancing correlation physics. We
predict that this velocity renormalization leads to an increase of the slow
and decrease of the fast Dirac velocity in TTG. We remark that a U(1)
generalization of our Gross-Neveu theory would correspond to intervalley
coherent order in TTG Christos  _et al._ (2021), which is expected to display
qualitatively similar behavior to the velocity renormalization and critical
exponents we have found based on studies of the XY-Gross-Neveu modelRosenstein
_et al._ (1993); Roy  _et al._ (2016); Zerf  _et al._ (2017, 2016); Li  _et
al._ (2017); Classen  _et al._ (2017); Iliesiu  _et al._ (2018). The
generalization of our theory to experimentally relevantPark  _et al._ (2021a);
Zhang  _et al._ (2021) twisted pentalayar graphene and other van-der-Waals
materials with more than two Dirac velocities is left for future work. Beyond
our application to TTG at charge neutrality, where $N_{+}=4$ and $N_{-}=8$, we
highlight that our theory may also apply to TTG in the vicinity of correlated
insulating states at integer filling per moiré unit cell. Following the
paradigm of a cascadeZondiner  _et al._ (2020); Wong  _et al._ (2020) of
resetting Dirac node energies, the effective theory may again be described by
Eq. (1), yet with a reduced number of slow modes $N_{-}<8$.

Acknowledgments: It is a pleasure to thank Lukas Janssen, Pablo Jarillo-
Herrero, Eslam Khalaf, Walter Metzner, Pavel Ostrovsky for useful discussions.
JHP is partially supported by the Air Force Office of Scientific Research
under Grant No. FA9550-20-1-0136 and the Alfred P. Sloan Foundation through a
Sloan Research Fellowship. JHP and EJK acknowledge hospitality by the Aspen
Center for Physics, where part of this work was completed and which is
supported by National Science Foundation grant PHY1607611. Work at BNL is
supported by the U.S. Department of Energy (DOE), Office of Basic Energy
Sciences, under Contract No. DE- SC0012704.

## References

  * Suárez Morell  _et al._ (2010) E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Phys. Rev. B 82, 121407 (2010).
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * 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, _et al._ , Nature 556, 80 (2018a).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b).
  * 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, _et al._ , Nature materials 19, 861 (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).
  * Wu  _et al._ (2018) F. Wu, T. Lovorn, E. Tutuc, and A. H. MacDonald, Phys. Rev. Lett. 121, 026402 (2018).
  * Zhou  _et al._ (2021) Y. Zhou, D. N. Sheng, and E.-A. Kim, “Quantum phases of transition metal dichalcogenide moiré systems,” (2021), arXiv:2105.07008 .
  * Scherer  _et al._ (2021) M. M. Scherer, D. M. Kennes, and L. Classen, “$\mathcal{N}=4$ chiral superconductivity in moiré transition metal dichalcogenides,” (2021), arXiv:2108.11406 .
  * Schrade and Fu (2019) C. Schrade and L. Fu, Phys. Rev. B 100, 035413 (2019).
  * Zhao  _et al._ (2021) S. Zhao, N. Poccia, X. Cui, P. A. Volkov, H. Yoo, R. Engelke, Y. Ronen, R. Zhong, G. Gu, S. Plugge, _et al._ , (2021), arXiv:2108.13455 .
  * Volkov  _et al._ (2021) P. A. Volkov, S. Y. F. Zhao, N. Poccia, X. Cui, P. Kim, and J. Pixley, “Josephson effects in twisted nodal superconductors,” (2021), arXiv:2108.13456 .
  * Tummuru  _et al._ (2021) T. Tummuru, S. Plugge, and M. Franz, “Josephson effects in twisted cuprate bilayers,” (2021), arXiv:2108.13457 .
  * 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).
  * Shen  _et al._ (2020) C. Shen, Y. Chu, Q. Wu, N. Li, S. Wang, Y. Zhao, J. Tang, J. Liu, J. Tian, K. Watanabe, _et al._ , Nature Physics 16, 520 (2020).
  * Liu  _et al._ (2020) X. Liu, Z. Hao, E. Khalaf, J. Y. Lee, Y. Ronen, H. Yoo, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, _et al._ , Nature 583, 221 (2020).
  * Park  _et al._ (2021a) J. M. Park, Y. Cao, L. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, (2021a), arXiv:2112.10760 .
  * Zhang  _et al._ (2021) Y. Zhang, R. Polski, C. Lewandowski, A. Thomson, Y. Peng, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, J. Alicea, _et al._ , “Ascendance of superconductivity in magic-angle graphene multilayers,” (2021), arXiv:2112.09270 .
  * Khalaf  _et al._ (2019) E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019).
  * 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).
  * Park  _et al._ (2021b) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021b).
  * Hao  _et al._ (2021) Z. Hao, A. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * 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, _et al._ , Nature 574, 653 (2019).
  * Wu  _et al._ (2021) S. Wu, Z. Zhang, K. Watanabe, T. Taniguchi, and E. Y. Andrei, Nature Materials 20, 488 (2021).
  * Lian  _et al._ (2021) B. Lian, Z.-D. Song, N. Regnault, D. K. Efetov, A. Yazdani, and B. A. Bernevig, Phys. Rev. B 103, 205414 (2021).
  * Parker  _et al._ (2021) D. E. Parker, T. Soejima, J. Hauschild, M. P. Zaletel, and N. Bultinck, Phys. Rev. Lett. 127, 027601 (2021).
  * Bultinck  _et al._ (2020) N. Bultinck, E. Khalaf, S. Liu, S. Chatterjee, A. Vishwanath, and M. P. Zaletel, Phys. Rev. X 10, 031034 (2020).
  * Brillaux  _et al._ (2020) E. Brillaux, D. Carpentier, A. A. Fedorenko, and L. Savary, “Nematic insulator at charge neutrality in twisted bilayer graphene,” (2020), arXiv:2008.05041 .
  * Da Liao  _et al._ (2021) Y. Da Liao, J. Kang, C. N. Breiø, X. Y. Xu, H.-Q. Wu, B. M. Andersen, R. M. Fernandes, and Z. Y. Meng, Phys. Rev. X 11, 011014 (2021).
  * Kang and Vafek (2020) J. Kang and O. Vafek, Phys. Rev. B 102, 035161 (2020).
  * Ochi  _et al._ (2018) M. Ochi, M. Koshino, and K. Kuroki, Phys. Rev. B 98, 081102 (2018).
  * Zhang  _et al._ (2020) Y. Zhang, K. Jiang, Z. Wang, and F. Zhang, Phys. Rev. B 102, 035136 (2020).
  * Liao  _et al._ (2021) Y.-D. Liao, X.-Y. Xu, Z.-Y. Meng, and J. Kang, Chinese Physics B 30, 017305 (2021).
  * Christos  _et al._ (2021) M. Christos, S. Sachdev, and M. S. Scheurer, “Correlated insulators, semimetals, and superconductivity in twisted trilayer graphene,” (2021), arXiv:2106.02063 .
  * Roy  _et al._ (2016) B. Roy, V. Juricić, and I. F. Herbut, Journal of High Energy Physics 2016, 18 (2016).
  * Sitte  _et al._ (2009) M. Sitte, A. Rosch, J. S. Meyer, K. A. Matveev, and M. Garst, Phys. Rev. Lett. 102, 176404 (2009).
  * König and Pixley (2018) E. J. König and J. H. Pixley, Phys. Rev. Lett. 121, 083402 (2018).
  * Roy  _et al._ (2018) B. Roy, M. P. Kennett, K. Yang, and V. Juricić, Phys. Rev. Lett. 121, 157602 (2018).
  * Boettcher (2020) I. Boettcher, Phys. Rev. Lett. 124, 127602 (2020).
  * Mandal (2021) I. Mandal, Physics Letters A 418, 127707 (2021).
  * Sheehy and Schmalian (2007) D. E. Sheehy and J. Schmalian, Phys. Rev. Lett. 99, 226803 (2007).
  * Barlas  _et al._ (2007) Y. Barlas, T. Pereg-Barnea, M. Polini, R. Asgari, and A. H. MacDonald, Phys. Rev. Lett. 98, 236601 (2007).
  * Hwang  _et al._ (2007) E. H. Hwang, B. Y.-K. Hu, and S. Das Sarma, Phys. Rev. Lett. 99, 226801 (2007).
  * Elias  _et al._ (2011) D. Elias, R. V. Gorbachev, A. Mayorov, S. Morozov, A. Zhukov, P. Blake, L. Ponomarenko, I. V. Grigorieva, K. S. Novoselov, F. Guinea, _et al._ , Nature Physics 7, 701 (2011).
  * Yu  _et al._ (2013) G. L. Yu, R. Jalil, B. Belle, A. S. Mayorov, P. Blake, F. Schedin, S. V. Morozov, L. A. Ponomarenko, F. Chiappini, S. Wiedmann, U. Zeitler, M. I. Katsnelson, A. K. Geim, K. S. Novoselov, and D. C. Elias, Proceedings of the National Academy of Sciences 110, 3282 (2013).
  * Ulybyshev  _et al._ (2021) M. Ulybyshev, S. Zafeiropoulos, C. Winterowd, and F. Assaad, “Bridging the gap between numerics and experiment in free standing graphene,” (2021), arXiv:2104.09655 .
  * (47) See supplementary materials for a motivation of our model in connection to twisted trilayer graphene, the explicit calculations of all diagrams, as well as an analysis of the RG flow.
  * Hejazi  _et al._ (2019) K. Hejazi, C. Liu, H. Shapourian, X. Chen, and L. Balents, Phys. Rev. B 99, 035111 (2019).
  * Ray  _et al._ (2018) S. Ray, M. Vojta, and L. Janssen, Phys. Rev. B 98, 245128 (2018).
  * Rosenstein  _et al._ (1989) B. Rosenstein, B. J. Warr, and S. H. Park, Phys. Rev. Lett. 62, 1433 (1989).
  * Rosenstein  _et al._ (1991) B. Rosenstein, B. J. Warr, and S. H. Park, Physics Reports 205, 59 (1991).
  * Khveshchenko and Paaske (2001) D. V. Khveshchenko and J. Paaske, Phys. Rev. Lett. 86, 4672 (2001).
  * Gracey (2018) J. A. Gracey, International Journal of Modern Physics A 33, 1830032 (2018).
  * Boyack  _et al._ (2019) R. Boyack, A. Rayyan, and J. Maciejko, Phys. Rev. B 99, 195135 (2019).
  * Moshe and Zinn-Justin (2003) M. Moshe and J. Zinn-Justin, Physics Reports 385, 69 (2003).
  * Rosenstein  _et al._ (1993) B. Rosenstein, Hoi-Lai Yu, and A. Kovner, Physics Letters B 314, 381 (1993).
  * Zerf  _et al._ (2017) N. Zerf, L. N. Mihaila, P. Marquard, I. F. Herbut, and M. M. Scherer, Phys. Rev. D 96, 096010 (2017).
  * Zerf  _et al._ (2016) N. Zerf, C.-H. Lin, and J. Maciejko, Phys. Rev. B 94, 205106 (2016).
  * Li  _et al._ (2017) Z.-X. Li, Y.-F. Jiang, S.-K. Jian, and H. Yao, Nature Communications 8, 314 (2017).
  * Classen  _et al._ (2017) L. Classen, I. F. Herbut, and M. M. Scherer, Phys. Rev. B 96, 115132 (2017).
  * Iliesiu  _et al._ (2018) L. Iliesiu, F. Kos, D. Poland, S. S. Pufu, and D. Simmons-Duffin, Journal of High Energy Physics 2018, 36 (2018).
  * Zondiner  _et al._ (2020) U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, _et al._ , Nature 582, 203 (2020).
  * Wong  _et al._ (2020) D. Wong, K. P. Nuckolls, M. Oh, B. Lian, Y. Xie, S. Jeon, K. Watanabe, T. Taniguchi, B. A. Bernevig, and A. Yazdani, Nature 582, 198 (2020).

Supplementary materials on

“Interaction induced velocity renormalization in magic angle twisted trilayer
graphene”
Laura Classen1,2, J.H. Pixley3, Elio J. König2

1 Condensed Matter Physics and Materials Science Division, Brookhaven National
Laboratory, Upton, New York 11973, USA
2 Max Planck Institute for Solid State Research, D-70569 Stuttgart, Germany
3Department of Physics and Astronomy, Center for Materials Theory, Rutgers
University, Piscataway, NJ 08854

These supplementary materials contain a motivation of our model, Sec. S1, in
connection to twisted trilayer graphene, and the explicit calculations for
mean field solution, Sec. S2, and RG equations, Sec. S3, stated in the main
text, as well as an analysis of the RG flow, Sec. S4.

##  S1 Derivation of effective model

We consider symmetrically twisted trilayer graphene in zero displacement field
near charge neutrality, which can be described by (uncoupled) fast and slow
Dirac fermions on the single-particle level. We assume that interactions lead
to an instability in a particle-hole channel with Ising symmetry. In that
case, the most general attraction between different fermion species in the
interaction channel under consideration is

$\displaystyle\mathcal{L}_{\rm int}$
$\displaystyle=-\sum_{\xi,\xi^{\prime}=\pm
1}\frac{g_{\xi\xi^{\prime}}^{2}}{2N}(\bar{\psi}_{\xi}\sigma_{z}\psi_{\xi})(\bar{\psi}_{\xi^{\prime}}\sigma_{z}\psi_{\xi^{\prime}})$
$\displaystyle\doteq\frac{N}{2}\vec{\phi}^{T}\left(\begin{array}[]{cc}g_{++}^{2}&g_{+-}^{2}\\\
g_{+-}^{2}&g_{--}^{2}\end{array}\right)^{-1}\vec{\phi}$ (S3)
$\displaystyle+\left(\bar{\psi}_{+}\sigma_{z}\psi_{+},\bar{\psi}_{-}\sigma_{z}\psi_{-}\right)\vec{\phi}.$
(S4)

Here, we have introduced the Hubbard-Stratonovich field
$\vec{\phi}=(\phi_{+},\phi_{-})^{T}$, which is composed of the order parameter
field in fast and slow sectors and we suppress the flavor index $\alpha$.
Importantly, the order parameter fields in fast and slow sectors have the same
symmetry properties, such that the microscopic interaction generally couples
them, i.e. $g_{+-}\neq 0$.

In the large N limit, the instability occurs when the inverse RPA propagator

$D^{-1}(\omega=0,\vec{q}=0)=\left(\begin{array}[]{cc}g_{++}^{2}&g_{+-}^{2}\\\
g_{+-}^{2}&g_{--}^{2}\end{array}\right)^{-1}+\frac{1}{N}\left(\begin{array}[]{cc}\Pi_{+}&0\\\
0&\Pi_{-}\end{array}\right)$ (S5)

acquires a zero mode (here $\Pi_{\pm}$ is the fermionic susceptibility at zero
external frequency and momentum). One may expand
$\vec{\phi}=\phi\hat{e}_{0}+\phi_{\perp}\hat{e}_{\perp}$, where
$D^{-1}(\omega=0,\vec{q}=0)\hat{e}_{0}=0$ at criticality, while
$D^{-1}(\omega=0,\vec{q}=0)\hat{e}_{\perp}=C\hat{e}_{\perp}$, $C>0$. Thus
$\hat{e}_{0}\hat{e}_{0}^{T}$ is the projector on the interaction channel which
first acquires a zero mode, as determined by the condition
$\det[D^{-1}(\omega=0,\vec{q}=0)]=0$ or

$1+(g_{++}^{2}\Pi_{+}+g_{--}^{2}\Pi_{-})/N=(g_{+-}^{4}-g_{++}^{2}g_{--}^{2})\Pi_{+}\Pi_{-}/N^{2}.$
(S6)

Clearly, this condition reduces the number of independent coupling constants
$g_{++},g_{--},g_{+-}$ from three to two.

In the main text, we tacitly project the interaction onto the critical channel
$\hat{e}_{0}$ and drop the massive modes $\hat{e}_{\perp}$ altogether (the
impact of virtual fluctuations can be assumed to be weak without changing
qualitative features of the theory). Technically, the projection amounts to
keeping only fluctuations $\vec{\phi}=\phi\hat{e}_{0}$ in Eq. (S4), i.e.

$\displaystyle\mathcal{L}_{\rm int}$
$\displaystyle=\frac{N\phi^{2}}{2}\hat{e}_{0}^{T}\left(\begin{array}[]{cc}g_{++}^{2}&g_{+-}^{2}\\\
g_{+-}^{2}&g_{--}^{2}\end{array}\right)^{-1}\hat{e}_{0}$ (S9)
$\displaystyle+\phi\left(\bar{\psi}_{+}\sigma_{z}\psi_{+},\bar{\psi}_{-}\sigma_{z}\psi_{-}\right)\hat{e}_{0}.$
(S10)

We define
$\bar{g}^{-2}=\hat{e}_{0}^{T}\left(\begin{array}[]{cc}g_{++}^{2}&g_{+-}^{2}\\\
g_{+-}^{2}&g_{--}^{2}\end{array}\right)^{-1}\hat{e}_{0}$ and absorb
$\phi\rightarrow\bar{g}\phi$. Then, by comparison with Eq. (1) of the main
text, $g_{+}=\bar{g}(1,0)^{T}\hat{e}_{0}$,
$g_{-}=\bar{g}(0,1)^{T}\hat{e}_{0}$.

We conclude this section with an estimate of $g_{+}$, $g_{-}$ in terms of
$g_{++},g_{+-},g_{--}$ in the important limit $|\Pi_{++}|\ll|\Pi_{--}|\sim
g_{++}\sim g_{+-}\sim g_{--}$. Criticality is driven by slow modes, leading to
$\hat{e}_{0}=(g_{-+}^{2},g_{--}^{2})/\sqrt{g_{-+}^{4}+g_{--}^{4}}$,
$\hat{e}_{\perp}=(g_{--}^{2},-g_{-+}^{2})/\sqrt{g_{-+}^{4}+g_{--}^{4}}$, with
the approximate mean field transition at $1+g_{--}\Pi_{-}/N=0$. By
consequence, we can estimate the parameters of Eq. (1) of the main text as
$g_{+}=g_{-+}^{2}/g_{--},g_{-}=g_{--}$.

##  S2 Model and mean field solution

As outlined in the main text and in the previous section of this supplement,
we consider the model

$\mathcal{L}=\sum_{\pm}\sum_{\alpha=1}^{N_{\pm}}\bar{\psi}_{\alpha,\pm}[\partial_{\tau}+v_{\pm}\vec{p}\cdot\vec{\sigma}+g_{\pm}\phi\sigma_{z}]\psi_{\alpha,\pm}+\frac{N}{2}\phi^{2}.$
(S11)

We assume $N_{+}+N_{-}=2N\gg 1$. It is convenient to parametrize

$g_{\pm}=\bar{g}\rho_{g}^{\pm 1/2},\quad v_{\pm}=\bar{v}\rho_{v}^{\pm 1/2},$
(S12)

where $\bar{g}=\sqrt{g_{+}g_{-}}$ and $\rho_{g}=g_{+}/g_{-}$ and analogously
for velocity. We absorb $\bar{g}$ into $\phi$, then

$\mathcal{L}=\sum_{\pm}\sum_{\alpha=1}^{N_{\pm}}\bar{\psi}_{\alpha,\pm}[\partial_{\tau}+v_{\pm}\vec{p}\cdot\vec{\sigma}+\rho_{g}^{\pm
1/2}\phi\sigma_{z}]\psi_{\alpha,\pm}+\frac{N}{2\bar{g}^{2}}\phi^{2}.$ (S13)

In the case, $N_{+}=N_{-}=N$ and $v_{+}=v_{-}=\rho_{g}=1$, this model is the
same as the Gross-Neveu(-Yukawa) model, see Ref. Rosenstein  _et al._ , 1991
of the main text. Note that the energy dimension of $\phi$ and $\psi$ is one,
and of $\bar{g}^{2}$ is -1.

###  S2.1 Fermionic Green’s function

When the $\phi$ field condenses, the inverse fermionic Green’s functions are

$G_{\pm}^{M_{\pm}}(\textbf{p})^{-1}=(i\epsilon-
v_{\pm}\vec{p}\cdot\vec{\sigma}-M_{\pm}\sigma_{z}),$ (S14)

where $M_{\pm}=\rho_{g}^{\pm 1/2}\langle\phi\rangle\equiv\rho_{g}^{\pm
1/2}\bar{M}$. We use the notation of a bold symbol
$\mathbf{p}=(\epsilon,\vec{p})$ for energy-momentum vectors. Inverting the
matrix yields

$G_{\pm}^{M_{\pm}}(\textbf{p})=-\frac{i\epsilon+v_{\pm}\vec{p}\cdot\vec{\sigma}+M_{\pm}\sigma_{z}}{\epsilon^{2}+v_{\pm}^{2}\vec{p}^{2}+M_{\pm}^{2}}.$
(S15)

We use non-zero $M_{\pm}$ to simplify some calculations, but generally we are
interested in the symmetric phase and consider $M_{\pm}\rightarrow 0$.

###  S2.2 Polarization operator and effective boson propagator

The polarization operator (or susceptibility), Fig. 2 a) of the main text, of
each fermionic species is

$\displaystyle\Pi_{\pm}(\textbf{q})$ $\displaystyle=$ $\displaystyle
N_{\pm}\rho_{g}^{\pm
1}\int_{\textbf{p}}\text{tr}^{\sigma}[G_{\pm}^{M_{\pm}}(\textbf{p}+\textbf{q}/2)G_{\pm}^{M_{\pm}}(\textbf{p}-\textbf{q}/2)]$
(S16) $\displaystyle=$ $\displaystyle\frac{N_{\pm}}{\pi}\frac{\rho_{g}^{\pm
1}}{v_{\pm}^{2}}\Big{[}-({v_{\pm}\Lambda/2}-M_{\pm}/2)$
$\displaystyle+\frac{4M_{\pm}^{2}+q_{\pm}^{2}}{4q_{\pm}}\arctan\left(\frac{q_{\pm}}{2M_{\pm}}\right)\Big{]}.$

We introduced the notation $q_{\pm}=\sqrt{\omega^{2}+v_{\pm}^{2}\vec{q}^{2}}$,
the ultraviolet cut-off $\Lambda$ ($\vec{p}^{2}<\Lambda^{2}$) and
$\textbf{q}=(\omega,\vec{q})$, $\textbf{p}=(\epsilon,\vec{p})$. Thus, to
leading order in $N\gg 1$ we get the effective bosonic theory

$\mathcal{L}_{\rm
eff}[\phi]=\frac{N}{2}\int_{\textbf{q}}\phi(-\textbf{q})\underbrace{[\frac{1}{\bar{g}^{2}}+\frac{1}{N}\sum_{\pm}\Pi_{\pm}(\textbf{q})]}_{\equiv
D^{-1}(\textbf{q})/N}\phi(\textbf{q}).$ (S17)

We see that, in the symmetric phase where $\bar{M}=0$, the effective mass of
bosons is given by

$\displaystyle m$
$\displaystyle=\frac{1}{\bar{g}^{2}}-\frac{\Lambda}{2\pi\bar{v}}\sum_{\pm}\alpha_{\pm}\rho_{v}^{\pm
1/2},$ (S18)

where $\alpha_{\pm}=\frac{N_{\pm}}{N}\frac{\rho_{g}^{\pm 1}}{\rho_{v}^{\pm
1}}$. We stated this in Eq. (3) of the main text.

$m=0$ defines the mean field transition. We remark that in Ref. Rosenstein
_et al._ , 1991 of the main text, the UV cutoff is defined as
${\frac{\bar{v}\pi\Lambda}{2}}={\Lambda_{\rm Rosenstein}}$. From Eq. (S17), we
obtain the effective bosonic propagator

$D(\textbf{q})=\frac{1}{N}\frac{1}{m+\sum_{\pm}\alpha_{\pm}q_{\pm}/(8\bar{v}^{2})}.$
(S19)

The dimension of the propagator is 1/energy. At the isotropic point we can set
$v_{+}=v_{-}$ and $\alpha_{\pm}=1$.

For later reference, we also remark the following relations here. In the
calculation of the fermionic self-energy we need

$D(\textbf{q}+\textbf{p})\stackrel{{\scriptstyle\textbf{p}\ll\textbf{q}}}{{\simeq}}D(\textbf{q})\left(1-D(\textbf{q})\frac{N}{4\bar{v}^{2}}\sum_{\pm}\frac{\alpha_{\pm}}{2}\frac{\langle\textbf{p},\textbf{q}\rangle_{\pm}}{q_{\pm}}\right),$
(S20)

with the scalar product
$\langle\textbf{p},\textbf{q}\rangle_{\pm}=\epsilon\omega+v_{\pm}^{2}\vec{p}\cdot\vec{q}$.
For the calculation of the two-loop correction to the bosonic mass, Fig. 2
d),e) of the main text, we calculate the second derivative of the polarization
operator with respect to the mass $\bar{M}=\langle\phi\rangle$, see Fig. S1.
We find

$\displaystyle\frac{1}{2}\frac{\partial^{2}}{\partial\bar{M}^{2}}\Pi_{\pm}(\textbf{q})$
$\displaystyle=$ $\displaystyle\frac{\rho_{g}^{\pm
1}}{2}\frac{\partial^{2}}{\partial M_{\pm}^{2}}\Pi_{\pm}(\textbf{q})$ (S21)
$\displaystyle=$ $\displaystyle\frac{N_{\pm}\rho_{g}^{\pm 2}}{\pi
v_{\pm}^{2}q_{\pm}}\left(\arctan\left(\frac{q_{\pm}}{2M_{\pm}}\right)-\frac{2M_{\pm}q_{\pm}}{q_{\pm}^{2}+4M_{\pm}^{2}}\right)$
$\displaystyle\simeq$ $\displaystyle N\frac{\rho_{g}^{\pm
1}}{\bar{v}^{2}}\frac{\alpha_{\pm}}{2q_{\pm}}.$

Figure S1: Diagrammatic representation of Eq. (S21). The differentiation with
respect to $\bar{M}$ introduces two external legs at zero momentum. In a
subsequent step, to calculate the diagrams of Fig. 2 d), e) of the main text,
the legs at finite q are contracted. This leads to the loop correction of the
zero momentum susceptibility.

##  S3 Loop Corrections and renormalization group

We here present calculations of the loop corrections to the fermionic, and
bosonic self energy, as well as details on the RG procedure.

###  S3.1 Fermionic self energy and vertex correction

The fermionic self energy and the vertex correction, Fig. 2 b) and c) of the
main text, are most efficiently calculated by studying the self energy in the
symmetry-broken phase

$\displaystyle\Sigma_{\pm}(\textbf{p})$
$\displaystyle=\int_{\textbf{q}}D(\textbf{q}-\textbf{p})\sigma_{z}G_{\pm}(\textbf{q})\sigma_{z}\rho_{g}^{\pm
1}$
$\displaystyle=-\int_{\textbf{q}}\frac{D(\textbf{q})}{q_{\pm}^{2}+M_{\pm}^{2}}M_{\pm}\sigma_{z}\rho_{g}^{\pm
1}$ $\displaystyle-\int_{\textbf{q}}\frac{D(\textbf{q})^{2}(i\omega-
v_{\pm}\vec{q}\cdot\vec{\sigma})}{q_{\pm}^{2}+M_{\pm}^{2}}\rho_{g}^{\pm
1}\sum_{\xi=\pm
1}{N\alpha_{\xi}/8\bar{v}^{2}}\frac{\langle\textbf{p},\textbf{q}\rangle_{\xi}}{q_{\xi}}$
$\displaystyle=-I_{1}^{\pm}M_{\pm}\sigma_{z}\rho_{g}^{\pm
1}-I_{2}^{\pm}i\epsilon\rho_{g}^{\pm 1}+I_{3}^{\pm}\rho_{g}^{\pm
1}v_{\pm}\vec{p}\cdot\vec{\sigma},$ (S22)

where we repeat from the main text that $\alpha_{\pm}=N_{\pm}\rho_{g}^{\pm
1}/(N\rho_{v}^{\pm 1})$. The leading logarithmic contributions are encoded in
the following dimensionless integrals

$\displaystyle I_{1}^{\pm}$
$\displaystyle=\int^{\Lambda}\frac{d^{2}q}{(2\pi)^{2}}\int\frac{d\omega}{2\pi}\frac{D(\textbf{q})}{q_{\pm}^{2}},$
(S23) $\displaystyle I_{2}^{\pm}$
$\displaystyle=\frac{N}{8\bar{v}^{2}}\int^{\Lambda}\frac{d^{2}q}{(2\pi)^{2}}\int\frac{d\omega}{2\pi}\frac{D(\textbf{q})^{2}}{q_{\pm}^{2}}\sum_{\xi}\alpha_{\xi}\frac{\omega^{2}}{q_{\xi}},$
(S24) $\displaystyle I_{3}^{\pm}$
$\displaystyle=\frac{N}{8\bar{v}^{2}}\int^{\Lambda}\frac{d^{2}q}{(2\pi)^{2}}\int\frac{d\omega}{2\pi}\frac{D(\textbf{q})^{2}}{q_{\pm}^{2}}\sum_{\xi}\alpha_{\xi}\frac{1}{2}\frac{v_{\xi}^{2}\vec{q}^{2}}{q_{\xi}}.$
(S25)

To evaluate these integrals, we change the variables
$\omega=\bar{v}x|\vec{q}|$, so that we obtain

$\displaystyle I_{1}^{\pm}$
$\displaystyle=\frac{2}{N\pi^{2}}\ln\left(\frac{\Lambda}{m\bar{v}}\right)g_{1}^{\pm}(v_{+},v_{-},g_{+},g_{-},N_{+},N_{-}),$
(S26) $\displaystyle I_{2}^{\pm}$
$\displaystyle=\frac{2}{3N\pi^{2}}\ln\left(\frac{\Lambda}{m\bar{v}}\right)g_{2}^{\pm}(v_{+},v_{-},g_{+},g_{-},N_{+},N_{-}),$
(S27) $\displaystyle I_{3}^{\pm}$
$\displaystyle=\frac{2}{3N\pi^{2}}\ln\left(\frac{\Lambda}{m\bar{v}}\right)g_{3}^{\pm}(v_{+},v_{-},g_{+},g_{-},N_{+},N_{-}),$
(S28)

where all of the following terms are unity in the isotropic limit

$\displaystyle g_{1}^{\pm}$
$\displaystyle=\int_{-\infty}^{\infty}dx\frac{1}{({x^{2}+\rho_{v}^{\pm
1}})[\sum_{\xi}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}]},$ (S29)
$\displaystyle g_{2}^{\pm}$
$\displaystyle=3\int_{-\infty}^{\infty}dx\frac{\sum_{\xi=\pm}\alpha_{\xi}\frac{x^{2}}{\sqrt{x^{2}+\rho_{v}^{\xi}}}}{({x^{2}+\rho_{v}^{\pm
1}})[\sum_{\xi=\pm}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}]^{2}},$ (S30)
$\displaystyle g_{3}^{\pm}$
$\displaystyle=3\int_{-\infty}^{\infty}dx\frac{\sum_{\xi=\pm}\alpha_{\xi}\frac{\rho_{v}^{\xi}/2}{\sqrt{x^{2}+\rho_{v}^{\xi}}}}{({x^{2}+\rho_{v}^{\pm
1}})[\sum_{\xi=\pm}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}]^{2}}.$ (S31)

Note that none of $g_{1,2,3}^{\pm}$ depends explicitly on $\bar{g}$ or
$\bar{v}$ and that $g_{2}^{\pm}=3g_{1}^{\pm}-2g_{3}^{\pm}$. In the isotropic
case $\rho_{v}=\rho_{g}=1$, we get $g_{1}^{\pm}=g_{2}^{\pm}=g_{3}^{\pm}=1$. We
also remark that we dropped terms with logarithmically slow dependence on
$\rho_{v},\rho_{g}$.

###  S3.2 Two-loop correction to bosonic self-energy

We calculate the leading correction to the polarization operator at zero
external momentum, Fig. 2 d) and e) of the main text, by exploiting Eq. (S21)
and obtain

$\displaystyle\delta\Pi_{\pm}(0)$
$\displaystyle=\int_{\textbf{q}}D(\textbf{q})N\frac{\rho_{g}^{\pm
1}}{\bar{v}^{2}}\frac{\alpha_{\pm}}{2q_{\pm}}$
$\displaystyle=\frac{1}{8\pi^{2}\bar{v}^{2}}\int dx\frac{\rho_{g}^{\pm
1}\alpha_{\pm}}{\sqrt{x^{2}+\rho_{v}^{\pm 1}}}$
$\displaystyle\times\int^{\Lambda}dq\frac{q}{m+\frac{q}{8\bar{v}}\sum_{\xi}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}}$
$\displaystyle=\frac{\Lambda}{2\pi\bar{v}}\rho_{g}^{\pm
1}\alpha_{\pm}g_{4}^{\pm}-m\frac{4\rho_{g}^{\pm
1}\alpha_{\pm}}{\pi^{2}}g_{5}^{\pm}\ln\left(\frac{\Lambda}{m\bar{v}}\right).$
(S32)

Again, we dropped slow, logarithmic dependence on $\rho_{v},\rho_{g}$ and
introduced two more functions which are unity in the isotropic limit

$\displaystyle g_{4}^{\pm}$
$\displaystyle=\int_{-\infty}^{\infty}dx\frac{2/\pi}{\sqrt{x^{2}+\rho_{v}^{\pm
1}}\sum_{\xi}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}}},$ (S33) $\displaystyle
g_{5}^{\pm}$
$\displaystyle=\int_{-\infty}^{\infty}dx\frac{2}{\sqrt{x^{2}+\rho_{v}^{\pm
1}}(\sum_{\xi}\alpha_{\xi}\sqrt{x^{2}+\rho_{v}^{\xi}})^{2}}.$ (S34)

This concludes the calculation of divergent leading order diagrams.

###  S3.3 Counter-terms and derivation of RG equations

We next present details on the derivation of the RG equations. We introduce
renormalized quantities to obtain the renormalized action

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{\pm,\alpha}\bar{\psi}_{\alpha,\pm}^{R}[Z_{1}^{\pm}\partial_{\tau}+Z_{1}^{\pm}Z_{v}^{\pm}v_{\pm}^{R}\vec{p}\cdot\vec{\sigma}+Z_{2}^{\pm}(\rho_{g}^{R})^{\pm
1/2}\phi^{R}\sigma_{z}]\psi_{\alpha,\pm}^{R}$
$\displaystyle+\frac{Z_{3}Z_{2}^{+}Z_{2}^{-}}{Z_{1}^{+}Z_{1}^{-}}\frac{N}{2(\bar{g}^{R})^{2}}(\phi^{R})^{2},$
(S35)

where the renormalization factors $Z_{i}$ contain counter-terms such that the
effective generating functional in terms of renormalized coupling constants
and source fields conjugate to renormalized fields (indicated by a superscript
$R$) is finite. The counterterms are introduced to render the loop diagrams
calculated above finite in the limit $\Lambda\rightarrow\infty$, but come at
the expense of endowing physical quantities with a non-trivial scaling with
respect to a running scale $\mu$. We can read off the renormalization factors
curing the one-loop diagrams from Eq. (S22)

$\displaystyle Z_{1}^{\pm}$
$\displaystyle=1-\ln(\Lambda/\mu)\frac{2}{3\pi^{2}N}g_{2}^{\pm}\rho_{g}^{\pm
1},$ (S36a) $\displaystyle Z_{v}^{\pm}$
$\displaystyle=1-\ln(\Lambda/\mu)\frac{2}{3\pi^{2}N}[g_{3}^{\pm}-g_{2}^{\pm}]\rho_{g}^{\pm
1},$ (S36b) $\displaystyle Z_{2}^{\pm}$
$\displaystyle=1+\ln(\Lambda/\mu)\frac{2}{\pi^{2}N}g_{1}^{\pm}\rho_{g}^{\pm
1}.$ (S36c)

Next, we determine $Z_{3}$. To this end we consider the bosonic mass

$\displaystyle\frac{D^{-1}(0)}{N}$
$\displaystyle=\frac{Z_{3}Z_{2}^{+}Z_{2}^{-}}{Z_{1}^{+}Z_{1}^{-}[\bar{g}^{R}]^{2}}+\frac{1}{N}\sum_{\pm}\left[\frac{Z_{2}^{\pm}}{Z_{1}^{\pm}}\right]^{2}\Pi_{\pm}(0)$
$\displaystyle=\frac{Z_{3}}{[\bar{g}^{R}]^{2}}-\sum_{\pm}\frac{\Lambda}{2\pi\bar{v}}\alpha_{\pm}(\rho_{v}^{\pm
1/2}-\frac{1}{N}\rho_{g}^{\pm 1}g_{4}^{\pm})$
$\displaystyle-m\ln\left(\frac{\Lambda}{\mu}\right)\frac{8\tilde{f}_{3}}{3\pi^{2}N}-\frac{2}{3\pi^{2}N}\frac{\Lambda}{2\pi\bar{v}}\ln\left(\frac{\Lambda}{\mu}\right)\left(\sum_{\pm}\pm(6g_{1}^{\pm}-2g_{2}^{\pm})\rho_{g}^{\pm
1}\right)\left(\sum_{\pm}\pm\alpha_{\pm}\rho^{\pm
1/2}\right)+\mathcal{O}(\ln(\mu/m)).$ (S37)

While it will be of no importance for the final RG equation, we remark that
all velocities entering this equation are bare velocities $v^{\rm
bare}=Z_{v}^{\pm}v^{R}_{\pm}$ containing infinities in the form of counter
terms of $\mathcal{O}(1/N)$. We introduced the function

$\tilde{f}_{3}=\sum_{\pm}\frac{3\alpha_{\pm}g_{5}^{\pm}-(3g_{1}^{\pm}-g_{3}^{\pm})}{2}\rho_{g}^{\pm
1}.$ (S38)

For the calculation of the RG equation, it is more convenient to express Eq.
(S37) using

$\displaystyle m_{\rm bare}$
$\displaystyle=\frac{Z_{3}}{[\bar{g}^{R}]^{2}}-\sum_{\pm}\frac{\Lambda}{2\pi\bar{v}}\alpha_{\pm}(\rho_{v}^{\pm
1/2}-\frac{1}{N}\rho_{g}^{\pm 1}g_{4}^{\pm})$
$\displaystyle-\frac{2}{3\pi^{2}N}\frac{\Lambda}{2\pi\bar{v}}\ln\left(\frac{\Lambda}{\mu}\right)\left(\sum_{\pm}\pm(6g_{1}^{\pm}-2g_{2}^{\pm})\rho_{g}^{\pm
1}\right)$ $\displaystyle\times\left(\sum_{\pm}\pm\alpha_{\pm}\rho^{\pm
1/2}\right).$ (S39)

Then, Eq. (S37) becomes

$\displaystyle m_{\rm physical}$ $\displaystyle=m_{\rm bare}-m_{\rm
bare}\frac{8\tilde{f}_{3}}{3\pi^{2}N}\ln\left(\frac{\Lambda}{\mu}\right).$
(S40)

###  S3.4 Renormalization group equations

Next, we express the renormalizations of $v_{\pm},g_{\pm}$ by renormalization
of relative and mean couplings, i.e.

$\displaystyle\frac{d\ln(\bar{v})}{d\ln(\mu)}$
$\displaystyle=-\frac{1}{3\pi^{2}N}[g_{3}^{+}\rho_{g}^{+1}-g_{2}^{+}\rho_{g}^{+1}+g_{3}^{-}\rho_{g}^{-1}-g_{2}^{-}\rho_{g}^{-1}],$
(S41) $\displaystyle\frac{d\rho_{v}}{d\ln(\mu)}$
$\displaystyle=\frac{2}{3\pi^{2}N}[-g_{3}^{+}\rho_{g}^{+1}+g_{3}^{-}\rho_{g}^{-1}+g_{2}^{+}\rho_{g}^{+1}-g_{2}^{-}\rho_{g}^{-1}]\rho_{v},$
(S42) $\displaystyle\frac{d\rho_{g}}{d\ln(\mu)}$
$\displaystyle=\frac{2}{3\pi^{2}N}[3g_{1}^{+}\rho_{g}^{+1}-3g_{1}^{-}\rho_{g}^{-1}+g_{2}^{+}\rho_{g}^{+1}-g_{2}^{-}\rho_{g}^{-1}]\rho_{g},$
(S43) $\displaystyle\frac{dm}{d\ln(\mu)}$
$\displaystyle=-\frac{8\tilde{f}_{3}}{3\pi^{3}N}m.$ (S44)

In the main text, we use the notation $g_{1}^{\pm}\rho_{g}^{\pm
1}=f_{1}^{\pm}$, $g_{3}^{\pm}\rho_{g}^{\pm 1}=f_{2}^{\pm}$,
$\tilde{f}_{3}=3f_{3}-\sum_{\pm}(3f_{1}^{\pm}-f_{2}^{\pm})/2$ in the
presentation of these RG equations, where we also use
$g_{2}^{\pm}\rho_{g}^{\pm 1}=(3g_{1}^{\pm}-2g_{3}^{\pm})\rho_{g}^{\pm
1}=3f_{1}^{\pm}-2f_{2}^{\pm}$.

###  S3.5 Wilsonian RG

The leading order $1/N$ corrections to $\rho_{v},\rho_{g},\bar{v}$ as well as
the scaling dimensions of the fields all stem from one-loop diagrams.
Therefore, the corresponding RG equations can be obtained using a Wilsonian RG
protocol including and iterative integration of fast modes with momenta
$q\in(\mu=\Lambda/b,\Lambda)$. However, the renormalization of the bosonic
mass stems from two-loop diagrams, Fig. 2 c) of the main text. Conventional
wisdom states that the Wilsonian RG may not be applicable to such higher-order
corrections, because of potential miscounting of fast degrees of freedom in
the phase space of the momentum shell approach. Here we explicitly rederive
the RG equations in the Wilsonian formalism and find full agreement with the
field theoretical counter-term approach. At the same time, as we believe that
the Wilsonian approach is more transparent to the condensed matter readership,
we keep it in this supplement despite the redundancy.

The integration of fast modes in the momentum shell readily lead to the
effective action

$\displaystyle S_{\rm eff}$ $\displaystyle=\int^{<}d^{2}x\int
d\tau\Big{[}\frac{Z_{1}^{+}Z_{1}^{-}}{Z_{3}^{W}Z_{2}^{+}Z_{2}^{-}}\frac{N}{2\bar{g}^{2}}(\phi^{<})^{2}$
$\displaystyle+\sum_{\pm,\alpha}\bar{\psi}_{\alpha,\pm}^{<}[(Z_{1}^{\pm})^{-1}\partial_{\tau}+(Z_{1}^{\pm}Z_{v}^{\pm})^{-1}v_{\pm}\vec{p}\cdot\vec{\sigma}$
$\displaystyle+(Z_{2}^{\pm})^{-1}(\rho_{g})^{\pm
1/2}\phi^{<}\sigma_{z}]\psi_{\alpha,\pm}^{<}\Big{]}.$ (S45)

Here, the notation $\int^{<}d^{2}x$ implicitly defines the real space
correspendent to momentum integrals $\int^{<}d^{2}p\dots=\int
d^{2}p\Theta(\Lambda/b-|\vec{p}|)\dots$. Similarly, fields with superscript <
are slow fields. The renormalization of coupling constants $Z_{1}^{\pm}$,
$Z_{v}^{\pm}$ and $Z_{2}^{\pm}$ directly follow from diagrams of the form Fig.
2 b),c) of the main text and are given by Eqs. (S36), keeping in mind that
$\mu=\Lambda/b$. The bosonic term has a prefactor

$\displaystyle\frac{Z_{1}^{+}Z_{1}^{-}}{Z_{3}^{W}Z_{2}^{+}Z_{2}^{-}\bar{g}^{2}}$
$\displaystyle=\frac{1}{\bar{g}^{2}}-\frac{\Lambda/b}{2\pi\bar{v}}\sum_{\pm}\alpha_{\pm}(\rho_{v}^{\pm
1/2}-g_{4}^{\pm}\rho_{g}^{\pm 1}/N)(b-1)$
$\displaystyle-\frac{4m}{\pi^{2}N}\ln(b)\sum_{\pm}\alpha_{\pm}\rho_{g}^{\pm
1}g_{5}^{\pm}.$ (S46)

To obtain this quantum correction, we consider diagrams Fig. 2 d) and e) of
the main text at a fast bosonic energy momentum with a RPA resummation keeping
only bubbles of fast fermions. However, since the momentum is fast, inclusion
of fast and slow fermions yields the same result, so that effectively we can
use the polarization operator introduced above.

Next, we rescale fields and space time in a manner to reproduce the appearance
of the original action

x $\displaystyle=b\textbf{x}_{R},$
$\displaystyle\psi_{\pm}^{R}(\textbf{x}_{R})$
$\displaystyle=b(Z_{1}^{\pm})^{-1/2}\psi_{\pm}^{<}(\textbf{x})$ (S47)
$\displaystyle\phi^{R}(\textbf{x}_{R})$
$\displaystyle=b\sqrt{\frac{Z_{1}^{+}Z_{1}^{-}}{Z_{2}^{+}Z_{2}^{-}}}\phi^{<}(\textbf{x}),$
(S48)

to get

$\displaystyle S$ $\displaystyle=\int
d^{2}x_{R}d\tau_{R}\Big{[}\frac{N}{2\bar{g}^{2}_{R}}[\phi^{R}(\textbf{x}_{R})]^{2}$
$\displaystyle+\sum_{\pm}\sum_{\alpha=1}^{N_{\pm}}\bar{\psi}^{R}_{\alpha,\pm}(\textbf{x}_{R})\\{\partial_{\tau_{R}}-iv_{\pm}^{R}\vec{\nabla}_{x_{R}}\cdot\vec{\sigma}$
$\displaystyle+[\rho_{g}^{R}]^{\pm
1/2}\phi^{R}(\textbf{x}_{R})\sigma_{z}\\}\psi^{R}_{\alpha,\pm}(\textbf{x}_{R})\Big{]}.$
(S49)

Now the momentum integral associated to $\int d^{2}x_{R}$ runs up to the
original cut-off $\Lambda.$ The renormalized coupling constants are

$\displaystyle v_{\pm}^{R}$
$\displaystyle=\left[1+\ln(b)\frac{2[g_{3}^{\pm}-g_{2}^{\pm}]\rho_{g}^{\pm
1}}{3\pi^{2}N}\right]v_{\pm},$ (S50) $\displaystyle\rho_{g}^{R}$
$\displaystyle=\left[1-\ln(b)\sum_{\pm}(\pm
1)\frac{2[g_{2}^{\pm}+3g_{1}^{\pm}]\rho_{g}^{\pm
1}}{3\pi^{2}N}\right]\rho_{g},$ (S51) $\displaystyle m^{R}$
$\displaystyle=bm\left(1-\ln(b)\frac{8\tilde{f}_{3}}{3\pi^{2}N}\right),$ (S52)

where the relationship between coupling constant and bosonic mass is

$m^{R}=\frac{1}{[\bar{g}^{R}]^{2}}-\frac{\Lambda}{2\pi\bar{v}}\sum_{\pm}\alpha_{\pm}(\rho_{v}^{\pm
1/2}-\frac{g_{4}^{\pm}\rho_{g}^{\pm 1}}{N}).$ (S53)

The renormalization group equations presented in the main text follow directly
from these equations by means of differentiation with respect to $\ln(b)$. To
get the renormalization of fields, we define $Z_{\psi}^{\pm}=Z_{1}^{\pm}$ and
$Z_{\phi}=\prod_{\pm}Z_{2}^{\pm}/Z_{1}^{\pm}$. Note that the renormalization
group equations reproduce the results obtained in the previous section on the
basis of the field theoretical counter-term method.

##  S4 Analysis of RG flow

In this section we present additional details on the analysis of the RG flow.

Figure S2: RG flow, Eqs. (5a), (5b) in the plane spanned by
$\rho_{v},\rho_{g}$ in the case $N_{+}=N_{-}$. The color code illustrates the
strength of the renormalization of the average velocity $\bar{v}$.

###  S4.1 RG flow in symmetric case

In the main text we discussed the RG flow for the case $N_{+}\neq N_{-}$ which
is most relevant to twisted trilayer graphene. Here, we briefly discuss the
case $N_{+}=N_{-}$, see Fig. S2. Note that in contrast to Fig. 3 of the main
text, the RG flow is now symmetric under $\rho_{v,g}\leftrightarrow
1/\rho_{v,g}$.

###  S4.2 RG flow near $\rho_{v}=1$

The RG equations can be calculated exactly to leading order in
$\rho_{v}=1+\delta\rho_{v}$, in which case

$\displaystyle\frac{d\delta\rho_{v}}{d\ln(b)}$
$\displaystyle=-\frac{32\rho_{g}^{2}}{15N\pi^{2}\left[n\left(\rho_{g}^{2}-1\right)+\rho_{g}^{2}+1\right]^{2}}\delta\rho_{v},$
(S54) $\displaystyle\frac{d\rho_{g}}{d\ln(b)}$
$\displaystyle=-\frac{16\rho_{g}\left(\rho_{g}^{2}-1\right)}{3N\pi^{2}\left[n\left(\rho_{g}^{2}-1\right)+\rho_{g}^{2}+1\right]}$
$\displaystyle+\frac{64\rho_{g}^{3}}{5N\pi^{2}\left[n\left(\rho_{g}^{2}-1\right)+\rho_{g}^{2}+1\right]^{2}}\delta\rho_{v},$
(S55) $\displaystyle\frac{d\ln(\bar{v})}{d\ln(b)}$
$\displaystyle=\frac{16n\rho_{g}^{2}}{15N\pi^{2}\left[n\left(\rho_{g}^{2}-1\right)+\rho_{g}^{2}+1\right]^{2}}\delta\rho_{v},$
(S56)

where $n=(N_{+}-N_{-})/(N_{+}+N_{-})$. Near the isotropic point, we can expand
these RG equations for $\delta\vec{\rho}=(\delta\rho_{v},\delta\rho_{g})$,
both $\delta\rho_{v,g}=\rho_{v,g}-1$,

$\displaystyle\frac{d\delta\vec{\rho}}{d\ln(b)}$
$\displaystyle=-\frac{8}{15\pi^{2}N}\left(\begin{array}[]{cc}1&0\\\
-6&10\end{array}\right)\delta\vec{\rho}.$ (S59)

This leads to scaling dimensions $-\frac{8}{15\pi^{2}N}$
[$-\frac{80}{15\pi^{2}N}$] with eigenvectors $(3,2)/\sqrt{13}$ [(0,1)],
independently on $n$.