eurm10 msam10

# Thin films spreading under the action of surface acoustic waves

Y. Li1 J.A. Diez2 J. D’Addesa3 O. Manor1 L.J. Cummings3 L. Kondic3 1
Department of Chemical Engineering, Technion - Israel Institute of Technology,
Haifa 32000, Israel
2 Instituto de Física Arroyo Seco, Universidad Nacional del Centro de la
Provincia de Buenos Aires, 7000 Tandil, Argentina
3 Department of Mathematical Sciences and Center for Applied Mathematics and
Statistics, New Jersey Institute of Technology, Newark, New Jersey 07102, USA

###### Abstract

Experimentally and theoretically, we consider the problem of a spreading drop
of millimetric thickness driven by a surface acoustic wave (SAW). The
experimental results show that, after initial transients, the drop spreads
with approximately constant speed while leaving a thin layer of fluid behind.
To better understand the system, we develop a theoretical model that couples
acoustic driving with the fluid mechanics of a spreading drop. An important
feature of the problem is that the strength of the acoustic driving depends on
the fluid thickness, with strong attenuation of the forcing in regions of the
domain covered by fluid. The presented model provides results consistent with
the experiments, particularly regarding drop height and spreading speed. The
constant spreading speed observed in the experiments is also captured by an
approximate traveling-wave solution, which provides predictions for the
dependence of drop thickness and spreading speed on the acoustic amplitude.

##  1 Introduction

Dynamic wetting of a solid substrate by a film of liquid is a common
occurrence. In natural systems it appears in surfactant films that wet the
lungs to promote breathing and the eyes to keep them moist (Veldhuizen &
Haagsman, 2000; Holly & Lemp, 1971), or in the spreading of water drops on
solid surfaces (such as raindrops on windows) (Dussan V. & Davis, 1974; Dussan
V. & Chow, 1983). In technological systems it is harnessed for actuating
microfluidic platforms (Atencia & Beebe, 2005; Stone et al., 2004; Whitesides,
2006), for cooling electronic circuits (Amon et al., 2001; Bar-Cohen et al.,
2006), for desalination (Fletcher et al., 1974) and for a variety of
manufacturing applications (Fendler, 1996; Nagayama, 1996; Wang et al., 2004).
Many different mechanisms govern dynamic wetting, where capillary forces at
the free surface of a liquid film or a meniscus contribute appreciably to the
dynamics of the film. In addition to capillary forces, dynamic wetting may be
powered by gravity, thermal or solutal Marangoni effects, or electrokinetics,
among others, see Oron et al. (1997); Craster & Matar (2009) for comprehensive
reviews. In the present work we focus on the contributions of mechanical
vibrations in the solid substrate, particularly MHz-frequency Rayleigh surface
acoustic waves (SAWs), to dynamic wetting (Morozov & Manor, 2018). Here, we
concentrate on the contribution of Eckart streaming – mass transport that
results from the attenuation of sound or ultrasound waves in the bulk of the
liquid film. While extensively observed and described in the scientific
literature (Shiokawa et al., 1989; Brunet et al., 2007, 2009; Yeo & Friend,
2014), dynamic wetting of a solid substrate by a liquid film, induced by
Eckart streaming, was not studied in detail to date. In this work we use thin
film (long wave) theory alongside physical experiments to elucidate the
phenomenon.

There are three types of acoustic streaming mechanisms by which acoustic waves
or waves at the free surface of a liquid or a solid boundary (surface acoustic
waves) generate streaming – steady flow along the path of the wave. Firstly, a
wave – acoustic wave or surface wave – in the solid or in liquid, invokes
periodic viscous flow within the viscous penetration length
$\delta\equiv\sqrt{2\mu/\rho\omega}$ away from the interface, where
$\mu,~{}\rho$ and $\omega$ are the viscosity and density of the liquid and the
wave angular frequency, respectively. Examples in which this problem was
considered include the work by Rayleigh on a standing acoustic wave grazing a
solid bed; by Schlichting (1932) on standing vibration of waves of infinite
and finite (respectively) wavelength in a solid along its surface; by Longuet-
Higgins (1953) on shallow ocean waves; and by Manor et al. (2012) on
propagating Rayleigh (surface acoustic) waves in a solid substrate. Convective
contributions due to surface waves invoke drift of liquid mass, which does not
attenuate away from the solid surface and is known variously as the Stokes
drift (Stokes, 1847), the Rayleigh law of streaming (Lighthill, 1978), and
Schlichting streaming (Yeo & Friend, 2014); although the physics and
mathematical modeling of the cases studied by Stokes, Rayleigh, and
Schlichting are similar, Rayleigh is the only one among them that explicitly
considered acoustic waves.

Secondly, in the presence of an acoustic wave in the liquid, whether due to an
actuator or an acoustic wave in the solid, one should expect acoustic
radiation pressure (Campbell, 1970; Shiokawa et al., 1989). Convective effects
– Reynolds stresses – within the wave increase the local steady pressure in
the bulk of the liquid; upon reflection of the wave at an interface between
phases of mismatched acoustic impedance (density multiplied by acoustic phase
velocity), the wave imposes convective stress, known as the acoustic radiation
pressure (King, 1934; Borgnis, 2010; Hasegawa et al., 2000), which is capable
of deforming and displacing the interface (Alzuaga et al., 2005; Biwersi et
al., 2000; Rajendran et al., 2022).

Finally, a third (and the most well-known) mechanism that may appear due to
the attenuation of sound or ultrasound waves in fluid is Eckart streaming
(Eckart, 1948; Lighthill, 1978; Nyborg, 1952), where wave intensity is
converted to fluid momentum, characterized by an intense vortical flow field
that appears in the bulk fluid. This is the main acoustic mechanism for the
actuation of fluid in micro-channels (Yeo & Friend, 2014; Wixforth et al.,
2004; Wixforth, 2003) and in drop microfluidics (Guttenberg et al., 2004;
Brunet et al., 2007, 2009, 2010).

These different acoustofluidic mechanisms govern mass transport at different
length scales (Rezk et al., 2014). In the presence of sound or ultrasound
waves in the fluid, the surface effect, i.e., Stokes drift, governs the
transport of mass at separations from the solid that are smaller than the
acoustic wavelength in the fluid. Further away, mass transport in the fluid is
governed by Eckart streaming. Acoustic radiation pressure effects may be
relevant in the former case, and also in the latter case when Eckart streaming
is weak. Moreover, SAW diffraction adds an additional spatial dimension to
flow actuation by SAW. While the SAW wavelength usually determines the lower
spatial resolution of the resulting acoustic flow, SAW diffraction allows for
breaking the SAW wavelength limit as a measure for spatial flow resolution
(Fakhfouri et al., 2018).

There are abundant studies on the dynamic wetting of thin liquid films that
are governed by Stokes drift and acoustic radiation pressure. Rayleigh surface
acoustic waves (SAWs) were recently shown to power the continuous spreading of
sub-micrometer to hundreds of micrometer thick films of fully wetting silicon
oils (Rezk et al., 2012, 2014) and partially wetting water/surfactant mixtures
(Altshuler & Manor, 2015, 2016; Mhatre et al., 2016) along or opposite the
path of the propagating SAW; water mixtures evaporate leaving behind SAW-
induced deposit morphology. However, we are not aware of studies that consider
in detail the acoustic/capillary balance of dynamic wetting in liquid films
that are sufficiently thick that mass transport therein, and hence film
dynamics, are actuated by Eckart streaming.

In this paper, we use physical experiments, theoretical modeling and
simulations to elucidate the physics of dynamic wetting powered by Eckart
streaming under the action of MHz-frequency Rayleigh-type SAWs. We describe
our experimental findings in Sec. 2 and then proceed with describing the
theoretical model in Sec. 3. Theoretical results are presented in Sec. 4; we
first consider computational results in Secs. 4.1 and 4.2 and then proceed
with an exploration of translationally-invariant solutions in Sec. 4.3.
Section 5 is devoted to the summary and conclusions.

##  2 Experiment

Figure 1 shows the experimental setup. We generate a propagating 20 MHz
frequency surface acoustic wave (SAW) by applying a same-frequency sinusoidal
voltage signal to a piezoelectric actuator – a SAW device. The actuator
comprises a 5 nm titanium/1 $\mu$m aluminum interdigitated transducer (IDT)
fabricated atop lithium niobate (LiNbO3, Roditi International, UK) by standard
lift-off photolithography. The substrate used for the SAW device is 11 mm
$\times$ 24 mm in size, 0.5 mm thick, 128° Y-cut, X-propagating, single-
crystal piezoelectric lithium niobate, where X and Y are crystal axes
(Campbell & Jones, 1968). The actuator is integrated into the external
electrical signal using pogo pins (BC201403AD, Interconnect Devices, Inc.)
assembled in a 3D-printed elastomeric stage, which holds the actuator and is
connected to a signal generator (R&S SMB100A microwave signal generator) and
amplifier (model A10160, Tabor Electronics Ltd.). We place a paper cylinder
soaked with glycerol at the far end of the actuator from the IDTs from which
the SAW emanates. The glycerol-soaked paper absorbs the SAW and avoids SAW
reflections. We measure the SAW normal surface vibration at the actuator
surface and verify that it is a propagating wave using a scanning laser
Doppler vibrometer (MSA-500, Polytech).

Figure 1: (A) Upper schematic view of the experimental setup: Surface acoustic
wave (SAW) propagates from the interdigitated transducer (IDT) until it
reaches the acoustic absorber (comprised of glycerol-soaked paper placed on
the actuator under and to the right of the needle). The needle is of known
diameter 510 $\mu$m, placed for identifying spatial resolution in the images.
(a-e) Successive snapshots taken from an experiment monitoring the flow of a
silicone oil film. During the experiment, a drop of silicone oil is placed on
the horizontal surface (a), deforms due to the application of SAW (b), and
moves in the direction of SAW propagation (c-e).

We place an 8 mm3 ($8~{}\mu$l) puddle of silicone oil (50 cSt, 378356, Sigma-
Aldrich) atop the actuator using a pipette, approximately 3.5 mm away from the
IDTs, and introduce electrical signal at different voltage levels to induce
motion of the oil film, which dynamically wets the solid substrate along the
path of the SAW. Figure 2(a-c) shows a time-lapse (side view) of the silicone
oil dynamics, captured using the framework and camera of a goniometer (Data
Physics; OCA 15Pro). Figure 2(d) shows a top view obtained using a camera (EOS
R5, Canon) with a macro lens (RF 100 mm F2.8L MACRO IS USM, Canon). Prior to
the experiment we clean the SAW device using four different solvents: acetone
(AR-b, 99.8%, 67-641, Bio-Lab Ltd.), 2-propanol (AR-b, 99.8%, 67-63-0, Bio-Lab
Ltd.), ethanol (CP-p, 96%, 64-17-5, Bio-Lab Ltd.), and water (HPLC plus,
7732-18-5, Sigma-Aldrich). Finally, we dry the actuator using compressed air.

Figure 2: Two sets of images from experiments monitoring silicone oil powered
by SAWs of three different amplitudes (1.2 nm-1.4 nm), where (A-C) silicone
oil starts from the same location, moves in the direction of the SAW, and
reaches different distances (a-c) after 4.09 seconds. The corresponding side-
view videos at different acoustic power levels are used to analyze the speed
and profile of silicone oil during the movement. (D, d) show top view
snapshots (corresponding to (C,c)) of the oil. A thin film of oil, hardly
visible in the side view, can be seen more clearly behind the main body of
silicone oil in the top views. Figure 3: Successive snapshots taken from the
baseline experiments monitoring silicone oil spreading without a SAW, where
the outline of the oil drop is marked as a red dashed line. The time/radius
stamps are given at the top right corner of each figure part.

Figure 4: (a) Time evolution of the front position for several SAW amplitudes
for $A=0.52$ nm (black), $0.65$ nm (gray), $0.78$ nm (purple), $0.91$ nm
(cyan), $1.04$ nm (green), $1.17$ nm (magenta), $1.30$ nm (red), $1.43$ nm
(brown), $1.56$ nm (orange), $1.69$ nm (blue). The arrow points the direction
of increasing $A$. (b) Time variations of the maximum drop height for SAW
amplitudes $A=0$ nm (green, no SAW applied), $0.5$ nm (blue), $0.7$ nm
(orange), $1.3$ nm (yellow) and $1.6$ nm (purple). The symbols correspond to
the measurements and the lines are simple fits to guide the eye. We expect
errors of $\pm 50~{}\mu$m in the values of the film height due to the limited
resolution of the side view image and tracker software.

In our experiments, we vary the SAW displacement amplitude at the surface of
the actuator by varying the applied voltage level input to the SAW device
(there is a linear relationship between the two; see Appendix A for the
corresponding calibration graph). To confirm that the observed dynamics are
indeed induced by SAW, we also carried out control experiments, shown in Fig.
3, in which SAW is absent and all other parameters are kept the same. The oil
puddle in this experiment, which was repeated three times, slowly spreads
concentrically, as expected in the absence of a directional forcing mechanism.

We then employ the side view video, recorded by a goniometer camera, to
capture the position of the advancing three-phase contact line of the SAW-
actuated oil film using the open-source software _Tracker_ (Brown & Cox, 2009)
(see Fig. 4(a) for typical data). We also measure the maximum height of the
moving oil film using the public domain image processing software _ImageJ_
(typical results in Fig. 4(b)). Figure 4(a) shows that, as expected, the
larger the SAW amplitude $A$, the faster the oil film moves due to the
stronger induced leakage wave (Shiokawa et al., 1989). Another notable aspect
of SAW-induced dynamics is the contact line speed for early times ($<1$ s):
initially, it is rather small, and later it increases and then remains
constant. This change of contact line speed appears to be a consequence of the
initial change of shape of the drop (see also Fig. 1); while this is taking
place, contact line speed is small. After the initial deformation, the drop
shape remains unchanged, translating approximately uniformly (this will be
seen in Fig. 9 later in the text). Figure 4(b) shows consistent results for
the maximum drop height, with an increase of the height for early times, and
constant values for longer times. Note that this uniform translation and
constant film height differ significantly from the observed behavior for
spreading under a body force such as gravity (for gravity-driven spreading
such as in Fig. 3 the contact line speed decreases, and the drop becomes
thinner as it spreads). These experimental observations will be explained and
discussed within the framework of our theoretical model, which is described
next.

##  3 Model

To model the experiments described in Sec. 2, we consider a fluid traversed by
a sound wave. The fluid velocity, $\mathbf{V}$, can be written as

$\mathbf{V}=\mathbf{v}+\mathbf{v}_{s},$ (1)

where $\mathbf{v}$ is the “bulk” velocity of the fluid, associated with a
characteristic time of the flow, $t_{c}$, which is on the order of seconds (so
$|\mathbf{v}|\sim\ell/t_{c}$, if $\ell$ is a typical lengthscale of the flow),
and $\mathbf{v}_{s}$ is associated with the much faster timescale of the sound
wave, $T=f^{-1}\ll t_{c}$ (typically $T\sim 10^{-6}$ s). The idea behind Eq.
(1) is that $\mathbf{v}_{s}$ corresponds to an oscillatory flow generated by
the propagating sound wave, while $\mathbf{v}$ is the velocity due to external
forces (or to the wave propagation itself, as we will show below). Therefore,
the fast-time averages, $\langle\cdot\rangle=(1/T)\int_{0}^{T}\cdot\,dt$, of
$\mathbf{v}$ and $\mathbf{v}_{s}$, are given by

$\langle\mathbf{v}\rangle\approx\mathbf{v},\qquad\langle\mathbf{v}_{s}\rangle=0.$
(2)

We write now the Navier–Stokes (NS) equations for the fluid in the compact
form

$\frac{\partial(\rho\mathbf{V})}{\partial
t}+\rho\left(\mathbf{V}\cdot\nabla\right)\mathbf{V}+\mathbf{V}\,\nabla\cdot\left(\rho\mathbf{V}\right)=-\nabla
p+\mu\nabla^{2}\mathbf{V}+\left(\zeta+\frac{\mu}{3}\right)\nabla\left(\nabla\cdot\mathbf{V}\right),$
(3)

where $\rho$ is the density of a compressible fluid and $\mu$, $\zeta$ are the
first and second dynamic viscosities, respectively. Since we seek the
resulting flow field for large times (of the order of $t_{c}$), we take the
fast-time average of Eq. (3) as defined above. Thus, using Eqs. (1) and (2),
we find

$\frac{\partial(\rho\mathbf{v})}{\partial t}=-\nabla
p+\mathbf{F}_{s}+\mu\nabla^{2}\mathbf{v}+\left(\zeta+\frac{\mu}{3}\right)\nabla\left(\nabla\cdot\mathbf{v}\right),$
(4)

where

$\mathbf{F}_{s}=-\langle\rho_{0}\left(\mathbf{v}_{s}\cdot\nabla\right)\mathbf{v}_{s}+\mathbf{v}_{s}\nabla\cdot\left(\rho_{0}\mathbf{v}_{s}\right)\rangle,$
(5)

and $\rho_{0}=\langle\rho\rangle$ is the average (constant) density of the
fluid (note that $\langle\partial(\rho_{0}\mathbf{v}_{s})/\partial
t\rangle=0$). Equation (4) shows that the flow field $\mathbf{v}$, which
develops on the long time scale $t_{c}$, satisfies NS equations similar to
those for $\mathbf{V}$, but with an additional volume force, $\mathbf{F}_{s}$.
This flow field, $\mathbf{v}$, is known as _acoustic streaming flow_ ,
resulting from the high-frequency sound wave traversing the fluid (see, e.g.
Shiokawa et al. (1989)).

When the sound wave is a leaky wave generated at a solid substrate that
propagates in the $(x,z)$–plane (see Fig. 5), the fluid particle displacement,
$\mathbf{u}$, can be written as (Shiokawa et al., 1990)

$\displaystyle u_{x}$ $\displaystyle=$
$\displaystyle\Re\left[A\,e^{\imath(k_{L}x^{\prime}-\omega
t)+\imath\alpha_{1}k_{L}z}\right]=A\,e^{-k_{i}x^{\prime}-\alpha_{1}k_{i}z}\cos\left[k_{r}x^{\prime}+\alpha_{1}k_{r}z-\omega
t\right],$ (6) $\displaystyle u_{y}$ $\displaystyle=$ $\displaystyle 0,$ (7)
$\displaystyle u_{z}$ $\displaystyle=$ $\displaystyle\alpha_{1}u_{x},$ (8)

where $A$ is the wave amplitude, $k_{L}=k_{r}+\imath\,k_{i}$ is the complex
wavenumber of the leaky wave and $x^{\prime}=x-x_{\rm ref}(t)$, where $x_{\rm
ref}(t)$ is a reference $x$–coordinate determining the onset of attenuation,
whose origin will be discussed in Sec. 4. Also, we define

$\alpha_{1}=\sqrt{(c_{s}/c_{\rm oil})^{2}-1},$ (9)

where $c_{s}$ and $c_{\rm oil}$ are the SAW phase velocities in the substrate
and silicon oil, respectively. By choosing (Roditi, 2023) $c_{s}=3880$ m/s and
(PrecisionFlow, 2023) $c_{\rm oil}=1350$ m/s, we obtain that
$\alpha_{1}=2.69$. Thus, the Rayleigh angle, $\theta$, is given by
$\theta=\arctan\alpha_{1}=69.65^{\circ}$. The SAW attenuation factor is given
by (Arzt et al. (1967))

$k_{i}=\frac{1}{l_{i}}=\frac{\rho_{\rm oil}}{\rho_{s}}\frac{c_{\rm
oil}}{c_{s}^{2}}f,$ (10)

where $l_{i}$ is the attenuation length, and $\rho_{s}$ is the density of the
substrate. For $f=20$ MHz, $\rho_{s}=4.65$ g/cm3 and $\rho_{0}=0.96$ g/cm3, we
obtain $l_{i}=2.7$ mm. Note that the SAW wavelengths in the substrate and the
oil are $\lambda_{s}=c_{s}/f=194$ $\mu$m and $\lambda_{\rm oil}=c_{\rm
oil}/f=67.5$ $\mu$m, respectively.

Figure 5: Schematic of a liquid drop (blue curve) forced by a SAW traveling
from left to right, with _rear_ and _front_ contact lines at $x_{r}(t)$ and
$x_{f}(t)$, respectively. The thin red line represents the amplitude $A$ of
the leaky SAW at the substrate and shows how it is attenuated by liquid (the
black dashed line is the envelope of the decaying amplitude in the liquid).
The Rayleigh angle is denoted by $\theta$.

The fluid particle displacement in Eqs. (6)–(8) is related to the Eulerian
velocity field, $\mathbf{v}_{s}=(v_{s,x},v_{s,y},v_{s,z})$, by

$v_{s,x}=\frac{du_{x}}{dt},\qquad v_{s,y}=\frac{du_{y}}{dt},\qquad
v_{s,z}=\frac{du_{z}}{dt},$ (11)

which, on substitution into Eq. (5), leads to explicit expressions for the
force components,

$\displaystyle{\rm F}_{s,x}$ $\displaystyle=$
$\displaystyle\left(1+\alpha_{1}^{2}\right)k_{i}\rho_{0}A^{2}\omega^{2}e^{-2k_{i}\left(x^{\prime}+\alpha_{1}z\right)},$
(12) $\displaystyle{\rm F}_{s,y}$ $\displaystyle=$ $\displaystyle 0,$ (13)
$\displaystyle{\rm F}_{s,z}$ $\displaystyle=$
$\displaystyle\alpha_{1}F_{s,x}.$ (14)

The force $\mathbf{F}_{s}$ is conservative, i.e., it can be expressed as
$\mathbf{F}_{s}=-\nabla\phi_{s}$ where the potential, $\phi_{s}$, is given by

$\phi_{s}(x,z)=\phi_{0}\,e^{-2k_{i}\left(x^{\prime}+\alpha_{1}z\right)},\quad\mbox{with}\quad\phi_{0}=\frac{1+\alpha_{1}^{2}}{2}\rho_{0}A^{2}\omega^{2}.$
(15)

Therefore, the NS equations for the velocity field $\mathbf{v}$, under the
influence of both SAW and gravitational forces, can be written as

$\rho_{0}\frac{\partial\mathbf{v}}{\partial t}=-\nabla
P+\mu\nabla^{2}\mathbf{v},$ (16)

where we assume that the $\mathbf{v}$ field is incompressible
($\nabla\cdot\mathbf{v}=0$), and we have defined the _effective pressure_ ,
$P$, as

$P=p+\phi_{s}+\phi_{g},$ (17)

with $\phi_{s}$ given by Eq. (15) and $\phi_{g}(z)=\rho_{0}gz$, corresponding
to a gravitational body force,
$\mathbf{F}_{g}=-\nabla\phi_{g}=-\rho_{0}g\,\hat{\mathbf{k}}$, where
$\hat{\mathbf{k}}$ is a unit vector in the out-of-plane direction.

In this work, we consider incompressible, low Reynolds number free surface
flows under the long-wave approximation so that all inertial terms on the
left-hand side of Eqs. (16) are negligible with respect to the viscous and
pressure terms. Carrying out a standard long-wave expansion of (16) and
imposing the usual stress boundary conditions at the free surface (vanishing
tangential stress and normal stress balanced by surface tension, see e.g.
Kondic (2003) for details), we obtain the following fourth-order differential
equation for the film height, $z=h(x,y,t)$,

$3\mu\frac{\partial h}{\partial t}-\nabla\cdot\left(h^{3}\nabla P\right)=0,$
(18)

where

$P(x,y,t)=-\gamma\nabla^{2}h+\rho_{0}gh+\phi_{s}(x,h).$ (19)

Here, $\gamma$ is the coefficient of surface tension at the oil-air interface,
assumed constant.

We nondimensionalize the problem using an arbitrary length scale, $\ell$, a
timescale $t_{c}$ such that leading order terms balance in (18), and a
pressure scale $p_{c}$ based on the capillary contribution in (19), so that we
have

$(x,y)=\ell(\tilde{x},\tilde{y}),\quad(x_{\rm ref},h)=\ell(\tilde{x}_{\rm
ref},\tilde{h}),\quad t=t_{c}\tilde{t}=\frac{3\mu\ell}{\gamma}\tilde{t},\quad
P=p_{c}\tilde{P}=\frac{\gamma}{\ell}\tilde{P},$ (20)

where the tilded quantities are dimensionless. Thus, the dimensionless
governing equation is

$\frac{\partial\tilde{h}}{\partial\tilde{t}}-\tilde{\nabla}\cdot\left(\tilde{h}^{3}\tilde{\nabla}\tilde{P}\right)=0,$
(21)

where

$\tilde{P}=-\tilde{\nabla}^{2}\tilde{h}+{\rm Bo}\,\tilde{h}+{\cal
S}\,\tilde{\phi}_{s}(\tilde{x},\tilde{h}),$ (22)

with

$\tilde{\phi}_{s}(\tilde{x},\tilde{h})=e^{-2\tilde{k}_{i}[(\tilde{x}-\tilde{x}_{\rm
ref})+\alpha_{1}\,\tilde{h}]},$ (23)

and we have defined $\tilde{k}_{i}=\ell k_{i}$. The dimensionless parameters
introduced here are given by

${\rm Bo}=\frac{\rho_{0}\ell^{2}}{\gamma}g=\frac{\ell^{2}}{a^{2}},\qquad{\cal
S}=\frac{\ell}{\gamma}\phi_{0},$ (24)

where $a=\sqrt{\gamma/(\rho_{0}g)}$ is the capillary length and $\phi_{0}$ was
defined in (15). In our simulations, motivated by a typical drop size in the
experiments, we set $\ell=1$ mm. For simplicity, we omit the tilde from now on
so that all variables will be dimensionless unless specific units are
explicitly written, mainly when comparison with experimental data is required.

##  4 Theoretical results

In this section, we discuss the predictions of the model developed in Sec. 3,
and compare them with the experimental results. In Sec. 4.1 we present results
for a chosen reference case, and then in Sec. 4.2, we discuss the influence of
model parameters, in particular the acoustic amplitude, on the results. In
Sec. 4.3 we discuss (approximate) traveling-wave solutions. All simulations
are carried out in two physical dimensions, assuming translational invariance
in the transverse $y$-direction. The discussion of the influence of the three-
dimensional configuration will be presented elsewhere.

###  4.1 Main features of the results

In our simulations, we solve Eq. (21) numerically using COMSOL™; see Appendix
B for the description of the implementation. As the initial condition, we
consider the two–dimensional parabolic drop of cross-sectional area

$A_{d}=2\int_{0}^{r_{d}}h_{d}\left[1-\left(\frac{x}{r_{d}}\right)^{2}\right]\,dx=\frac{4}{3}h_{d}r_{d},$
(25)

which remains constant during the drop evolution. The drop height parameter
$h_{d}$ is obtained from the known experimental (3D) drop volume, $V_{d}$, as
$h_{d}=2V_{d}/(\pi r_{d}^{2})$, using the formula for the volume $V_{d}$ of a
3D parabolic cap of base radius $r_{d}$ and height $h_{d}$. The 2D parabolic
initial condition is centered at $x=x_{d}$, and is written as

$h(x,0)=\left\\{\begin{array}[]{ll}h_{p},&x<x_{r}(0),\\\
\left(h_{d}-h_{p}\right)\left[1-\left(\frac{x-x_{d}}{r_{d}}\right)^{2}\right]+h_{p},&x_{r}(0)\leq
x\leq x_{f}(0),\\\ h_{p},&x>x_{f}(0),\end{array}\right.$ (26)

where $x_{r}(0)=x_{d}-r_{d}$ and $x_{f}(0)=x_{d}+r_{d}$ are the initial
positions of the rear and front contact lines of the drop, respectively (see
Fig. 5). Here, $h_{p}$ ($\ll h_{d}$) is the thickness of a precursor film,
introduced for numerical reasons in order to avoid the stress singularity
associated with a moving contact line. It is well known (see, e.g. Diez &
Kondic (2002)), that the step size $\Delta x$ in the spatial discretization
should be similar to $h_{p}$ to ensure numerical convergence, and therefore we
use $\Delta x=h_{p}$ in our simulations. The specific value assigned to
$h_{p}$ has only a minor influence on the results; e.g., if $h_{p}$ is doubled
or halved, the change of spreading speed or drop height is just 1-2%. All
reported results are obtained using $h_{p}=0.01$.

Recalling the introduction of the reference $x$-coordinate origin, $x_{\rm
ref}(t)$, in Sec. 3, we note that the SAW experiences minimal attenuation in
regions of the actuator where fluid is not present, with significant
attenuation only under the bulk droplet. To model this, we assume that the
acoustic potential, $\phi_{s}$, is constant (no attenuation) in regions where
the film height is smaller than some critical thickness $h^{*}$, and
attenuates where $h>h^{*}$,

$\phi_{s}=\left\\{\begin{array}[]{ll}1,&x<x_{1}^{\ast}(t)\,,\\\
e^{-2k_{i}\left[x-x_{1}^{\ast}(t)+\alpha_{1}\left(h-h^{\ast}\right)\right]},&x_{1}^{\ast}(t)\leq
x\leq x_{2}^{\ast}(t)\,,\\\
e^{-2k_{i}\left[x_{2}^{\ast}(t)-x_{1}^{\ast}(t)\right]},&x>x_{2}^{\ast}(t)\,,\end{array}\right.$
(27)

where $x_{1}^{\ast}$, $x_{2}^{\ast}$ are defined by
$h(x_{1}^{\ast},t)=h(x_{2}^{\ast},t)=h^{\ast}$ (this form of $\phi_{s}$
corresponds to the choice $x_{\rm ref}(t)=x_{1}^{*}(t)$). This model is
further motivated by the known fact that the Eckart streaming—the acoustic
streaming that pushes the oil film along the path of SAW in this work—governs
the flow when the oil film is sufficiently thick. Close to the solid surface
the flow is governed by a boundary layer flow (known as the Schlichting
streaming to some and Rayleigh streaming to others), generated by the SAW in
the solid. These mechanisms are observed to support the spreading of
$\lambda_{\rm oil}/4$ thick oil films opposite the path of the SAW. However,
thicker oil films are governed by Eckart streaming and spread along the path
of the SAW. Hence, we account for Eckart streaming-induced acoustic forcing
where the film thickness is greater than $h^{\ast}=\lambda_{\rm oil}/(4\ell)$
(Rezk et al., 2012, 2014; Altshuler & Manor, 2015, 2016).

We start by considering the experimental situation of an oil drop of volume
$8$ mm3 and a SAW pressure of amplitude $A=1.3$ nm, with experimental results
presented in Fig. 4. Table 1 gives the corresponding dimensional and
dimensionless parameters (see also Eqs. (9) and (10)).

Experiment | Simulation
---|---
$\rho_{0}=0.96$ g/cm3 | ${\cal S}=5.087$
$\mu=0.5$ Poise | ${\rm Bo}=0.452$
$\gamma=20.8$ dyn/cm | $A_{d}=2.12$
$\omega=2\pi\times 20$ MHz | $r_{d}=3.2$
$A=1.3$ nm | $h_{d}=0.497$
$l_{i}=2.7$ mm | $k_{i}=0.37$
Table 1: Dimensional parameters of the experiments and the corresponding
dimensionless constants used in the simulations for a $8$ mm3 volume drop. The
length scale is $\ell=1$ mm.

Figure 6(a) shows the evolution of the thickness profile obtained by solving
Eq. (21) for the parameter values in Table 1, subject to the initial condition
given by Eq. (26), and boundary conditions fixing the film thickness to
$h_{p}$ and requiring vanishing derivatives at the domain boundaries. We
observe that the simulated drop settles to a steady shape, which translates
uniformly. The steady shape of the bulk droplet has a pronounced “foot” at its
rear; though we do not yet fully understand the details of why this forms, it
is consistent with the overall drop asymmetry that develops, seen in Fig. 2.
We discuss these features of the results further below.

Figure 6(b) shows a single thickness profile at $t=2$ s (blue curve), which is
compared with the experimental data (black symbols). We note that capturing
the spreading experimental drop geometry (and velocity) requires image
analysis of the captured image data; the analysis adds uncertainty to the
measured data estimated to be $\pm 50$ $\mu$m. Moreover, repeating experiments
show some sensitivity to the initial condition, i.e., the thickness of the oil
drop placed upon the SAW device at the time SAW is initiated, albeit in all
cases the oil reaches a constant speed at long times.

Figure 6: (a) Thickness profiles from $t=0$ (red line) to $t=6$ s every
$\Delta t=0.5$ s for SAW amplitude $A=1.30\,$nm. Here, we use
$h^{\ast}=17\mu{\rm m}/h_{c}\approx(\lambda_{\rm oil}/4)/h_{c}$. (b)
Comparison between simulated (solid blue line) and experimental (black
symbols) thickness profiles at $t=2$ s.

###  4.2 Parametric dependence

In this section, we separate the discussion into two parts: first, we consider
early-time behavior, and then long-time evolution. Figure 7 shows the
computational results for the maximum drop height $h_{\rm max}$ for different
values of the acoustic amplitude, $A$. We show the results for two different
initial conditions, with the purpose of illustrating the relevance of the
initial drop shape. Focusing first on the initial condition used for Fig. 6
(solid curves in Fig. 7) with initial drop height of $h_{d}\approx 0.5$ mm
obtained from the experiments (see also Eq. (25)), we note different behavior
for small and large values of $A$, as could be anticipated from the
experimental results of Fig. 4(b). For the larger $A$-values, $h_{\rm max}$
initially increases rapidly but then decreases. For the two smallest
$A$-values, this initial rapid increase is not observed. In all cases, for
longer times, $h_{\rm max}$ slowly decreases towards a constant
($A$-dependent) value, although it does not completely saturate for the times
shown. One interpretation of these results is that for large $A$, the acoustic
forcing dominates capillary stresses and gravity, leading to an increase of
the drop height, while for smaller $A$, the acoustic forcing is not strong
enough to lead to initial drop height increase. For the other initial
condition (smaller initial drop height, shown by dashed lines in Fig. 7), the
details are different, but we again observe an initial increase of $h_{\rm
max}$ for sufficiently large values of $A$, followed by (in this case) a
slight decrease towards constant $A$-dependent values at large times.

Figure 7: Maximum fluid thickness, $h_{\rm max}$, as a function of time at
different SAW amplitudes, $A$, for two initial conditions of the drop:
$r_{d}=3.2$ mm (solid lines) and $r_{d}=3.9$ mm (dashed lines). The
cross–sectional fluid area, $A_{d}$, is the same, but the latter has a lower
initial maximum thickness than the former (see Eq. (25)). The values of are
$A$ are: $2.00$ nm (black), $1.69$ nm (blue), $1.30$ nm (red), $1.04$ nm
(magenta), $0.65$ nm (brown), $0.20$ nm (green) (same colors correspond to
same values of $A$).

Next, we shift our attention to long-time behavior, focusing first on the
dependence of drop profiles on SAW amplitude $A$. Figure 8 shows the long-time
drop profiles together with the pressure specified by Eq. (22) for four values
of $A$. While the maximum drop height increases with $A$, as expected, we note
that the thickness of the ‘foot’ directly behind the front increases with $A$
as well (varying between approximately 50 and 100 $\mu$m for the considered
values of $A$). In addition, we note a very thin (on the order of the
precursor scale, but a little larger) trailing film left behind the
translating drop. While both the ‘foot’ and the thin trailing film deposited
behind the drop still need to be fully analyzed and understood, we suspect
that these features of the results are related to the pressure (plotted by the
red lines with values shown on the right-hand axis in Fig. 8); note that the
‘foot’ forms in the region where the pressure is O(1), shown by the dashed
lines. On a different note, we point out the strong decrease of the pressure
due to the SAW attenuation, as the film height increases from the precursor
scale, as well as the pressure dip at the front contact line, caused by the
rapidly changing curvature there.

(a) $\,A=1.04$ nm

(b) $\,A=1.30$ nm

(c) $\,A=1.69$ nm

(d) $\,A=2.00$ nm

Figure 8: Simulated drop profiles (blue curves, left axis) and dimensionless
pressure $\tilde{P}=P/p_{c}$ (red curves, right axis) at $t=6$ s for four
different values of SAW amplitude $A$. The horizontal dashed line is at
$\tilde{P}=1$.

Focusing next on the spreading speed, Fig. 9 illustrates that for long times
both experiments and simulations show approximately constant spreading speed
(drop front position $x_{f}$ increasing linearly in time), which seems to be a
characteristic feature of the considered acoustic driving (note that we shift
the experimental results in the figure to illustrate better the long time
behavior; initial transients may account for differences at early times due to
the particulars of the initial conditions). These results suggest that the
spreading drop essentially translates without much of shape change for late
times. Figure 10(a) confirms this expectation: this figure shows that the
width of the spreading drop, ${\rm w}=x_{2}^{*}-x_{1}^{*}$, reaches a plateau
and remains constant in the case of larger acoustic amplitudes $A$; on
sufficiently long time scales (not shown in the figure), the $\rm w$-values
obtained for smaller $A$-values reach a plateau as well.

Figure 9: Comparison between numerical (solid lines) and experimental
(symbols) results for the SAW amplitude $A=1.69\,$nm (blue), $1.30\,$nm (red)
and $1.04\,$nm (magenta). We have added arbitrary $x$–shifts to the simulation
results to facilitate the comparison of the slope of $x_{f}$ versus $t$
between experiments and simulations for long times.

Figure 10: (a) Numerically calculated drop width, $\rm w$, as a function of
time for several values of $A$: $=0.20$ nm (green), $=0.65$ nm (brown),
$=1.04$ nm (magenta), $=1.30$ nm (red), $=1.69$ nm (blue), and $=2.00$ nm
(black). Note that the two smallest values of $A$ have not reached an
asymptotic value of $\rm w$ for $t=6$ s, but will do for later times. (b) Long
time values of the front speed, ${\rm v}_{f}$ (in mm/s, blue filled circles),
maximum thickness, $h_{\rm max}$ (mm, red hollow squares), and drop width,
$\rm w$ (cm, black filled squares), as functions of $A$. The points correspond
to the raw data, and the solid lines connecting them guide the eye.

To summarize our numerical results, Fig. 10(b) plots the numerically-predicted
long-time values (reached for $t\gtrsim 6$ s) of the front velocity, ${\rm
v}_{f}$, the maximum drop height, $h_{\rm max}$, and the drop width, ${\rm
w}=x_{2}^{\ast}-x_{1}^{\ast}$, as a function of $A$. We note that, while ${\rm
v}_{f}$ strongly increases with $A$, the dependence of $h_{\rm max}$ and $\rm
w$ on $A$ is much weaker. These results suggest that a good qualitative
description of the problem might be obtained by considering a traveling wave
type of solution. Such an approach is discussed in the following section.

###  4.3 Traveling wave solution

In view of the results presented so far, we now seek a traveling wave solution
of Eq. (21). Within this framework, we assume that the whole drop translates
by constant speed $U$, so that

$x_{r}(t)=Ut,\qquad x_{f}(t)={\rm w}+Ut,$ (28)

where (constant) $\rm w$ is the width of the moving drop. Here, we will
consider that the rear (front) contact lines of the traveling drop are given
by $x_{1}^{\ast}$ ($x_{2}^{\ast}$) (see Fig. 5), since the SAW force that
drives it corresponds to $h>h^{\ast}$. Then, we define
$\xi=x-x_{1}^{\ast}(t)=x-Ut$ and assume that $h(x,t)=h(\xi)$, so that Eq. (21)
becomes

$U\frac{dh}{d\xi}+\frac{d}{d\xi}\left[h^{3}\frac{d}{d\xi}\left(-\frac{d^{2}h}{d{\xi}^{2}}+{\rm
Bo}\,h+S\phi_{s}(\xi,h)\right)\right]=0$ (29)

where

$\phi_{s}(\xi,h)=e^{-2k_{i}[\xi+\alpha_{1}(h-h^{\ast})]}.$ (30)

This equation can be integrated once to yield

$Uh+h^{3}\frac{d}{d\xi}\left[-\frac{d^{2}h}{d{\xi}^{2}}+{\rm
Bo}\,h+S\phi_{s}(\xi,h)\right]=J,$ (31)

where $J$ represents the flux. The traveling wave solution must be calculated
for $0\leq\xi\leq\xi_{f}$, along with the following boundary conditions at
$\xi=0$

$h(0)=h^{\ast},\qquad h^{\prime}(0)=h^{\prime\prime\prime}(0)=0,$ (32)

where the prime denotes $d/d\xi$. The values of $h^{\prime\prime}(0)$, $U$ and
$\xi_{f}$ are determined by

$h(\xi_{f})=h^{\ast},\quad J(0)=J(\xi_{f}),\quad
A_{d}=\int_{0}^{\xi_{f}}h(\xi)\,d\xi.$ (33)

Note that $J$ corresponds to the flux within a thin film of thickness
$h^{\ast}$ that enters the drop at $\xi_{f}$ and comes out of it at $\xi=0$
(with respect to the reference frame fixed at the drop). Also, note that
$\xi_{f}$ is exactly the width of the traveling drop, $\rm w$.

We develop an iterative scheme to perform the numerical integration of Eq.
(31), under the conditions given by Eq. (33). We start by guessing values of
$h^{\prime\prime}(0)$, $U$ and $\xi_{f}$ to perform the integration of Eq.
(31) and modify them accordingly until these conditions are satisfied within a
small relative error (typically, $10^{-7}$). After obtaining convergence for a
given $A$, we use these converged values of $h^{\prime\prime}(0)$, $U$ and
$\xi_{f}$ as new guess values for $A+\Delta A$ (continuation). $\Delta A$ has
to be very small to get a new convergence (typically, $\Delta A=10^{-5}\,$nm).
Figure 11(a) shows the results of drop profiles for the values of $A$
considered so far. Clearly, drops became taller and narrower as $A$ increases.
To compare the results with the experimental ones, we note that for $A=1.3$
nm, the dimensional drop speed is ${\rm v}_{f}=(\ell/t_{c})U=2.75$ mm/s, which
is of the order of the experimental value. As expected, the value of $J$ is
small compared to unity, $J=3.26\times 10^{-4}$. Figure 11(b) shows the
dependence of drop speed, ${\rm v}_{f}$, its width, $\rm w$, and the maximum
thickness, $h_{\rm max}$, on the SAW amplitude, $A$. Note that the model
predicts that the speed depends approximately quadratically on the amplitude,
$A$, as expected since the drop thickness increases linearly with $A$. We note
that the trends shown in this figure for the travelling wave are similar to
those in Fig. 10(b) obtained using full numerical solution.

This approximate formulation shows that the simple traveling-wave solution can
capture, at least qualitatively, features of the experimental drops. Clearly,
however, this description cannot account for the rear flat film left behind
the drop and the fact that the actual rear contact line remains practically at
rest.

Figure 11: (a) Thickness profile of the traveling wave solution for: $A=0.50$
nm (green), $0.65$ nm (black), $1.04$ nm (magenta), $1.30$ nm (red), $1.69$ nm
(blue). (b) Front speed, ${\rm v}_{f}$ (mm/s, black line), drop width, $\rm w$
(mm, blue line), and maximum height, $h_{\rm max}$ (mm, red line), of the
traveling drop as a function of the SAW amplitude, $A$.

We note that we can use these results also to obtain an approximate value of
$U$ by applying a simple formula to estimate the scale of the front speed,
${\rm v}_{f}$. By setting $J\approx 0$ and neglecting both capillary and
gravity effects in Eq. (29), we have

$U\approx-
Sh^{2}\frac{d\phi_{s}}{d\xi}=2k_{i}Sh^{2}e^{-2k_{i}[\xi+\alpha_{1}(h-h^{\ast})]}.$
(34)

By setting $h=h_{d}/\ell$ and $\xi=2r_{d}/\ell$, this expression yields
$v_{f}=U\ell/t_{c}=0.4$ mm/s for $A=1.3$ nm, which is the right order of
magnitude.

##  5 Conclusions

This paper presents experimental results and a self-consistent model
describing the spreading of a silicone oil drop under the influence of
acoustic driving produced by a surface acoustic wave (SAW). Both simulations
and the theoretical model show consistent features, and in particular constant
spreading speed for long times. This feature of the results is specific to
this type of forcing, and is different from spreading under a body force (such
as gravity) where fluid thickness slowly decreases as a drop spreads, leading
to decreased spreading speed (which typically scales as $h^{2}$, where $h$ is
the drop thickness). Our understanding is that this feature of the results is
a consequence of the particular forcing imposed by SAW, which decreases as a
function of film height; therefore, in contrast to the gravity spreading, here
we have a situation where the drop is essentially pushed from behind, leaving
just a very thin (on the order of the precursor scale) film behind the
spreading drop. The drop profile itself includes a ‘foot’ behind the main
drop; this feature of the results still needs further investigation.

While not all details of the computational and experimental results are fully
understood, in general, we find very reasonable agreement between the
experiments and theoretical results, in particular regarding drop shapes,
spreading speeds, and the dependence of the results on the SAW amplitude. For
long times, both experiments and simulations show that a drop essentially
translates with constant speed and without change of shape. This finding
motivated us to formulate a traveling wave description of the system. The
resulting solution predicts a linear increase in the drop thickness and a
quadratic increase of the spreading speed with the acoustic amplitude $A$.

While many features of the experimental results have been rationalized using
the formulated model, we note that there is still much to be done. In
particular, further investigation is required to understand all the features
of the calculated drop profiles, and also three-dimensional aspects of the
results remain to be simulated and understood. Consideration of such effects
will be the subject of future work.

##  Appendix A Calibration of applied voltage versus SAW displacement

We use the normal displacement amplitude (particle displacement of the SAW),
$A$, to represent the SAW strength. In the experiments, we control this
amplitude by changing the applied voltage of the signal generator. Therefore,
we conducted an experiment to obtain a calibration curve between the applied
voltage and amplitude, as shown in Fig 12. We measure $A$ over a surface of
$1\times 1~{}\text{mm}^{2}$ about 2 mm away from the SAW actuator using laser
doppler vibrometer. We find an approximately linear dependence of $A$ on
voltage, consistently with the literature (Ballantine Jr et al., 1996).

Figure 12: A calibration plot of corresponding SAW displacement when different
voltages are applied to the actuator.

##  Appendix B Computational implementation

Equation (21) is written in a form convenient for the use of
COMSOL™Multiphysics PDE Coefficients Form. This package solves, by finite
elements, a vectorial equation for the unknown vector
$\vec{u}=(u_{1},u_{2},\dots,u_{N})^{\top}$. The equation is of the form

${\bf e}\frac{\partial^{2}\vec{u}}{\partial t^{2}}+{\bf
d}\frac{\partial\vec{u}}{\partial t}+\nabla\cdot\left(-{\bf
c}\nabla\vec{u}-{\bf\alpha}\,\vec{u}+{\bf\gamma}\right)+{\bf\beta}\,\nabla\vec{u}+{\bf
a}\,\vec{u}=\vec{f},$ (35)

where the coefficients of the $N$ scalar equations are in the matrices ${\bf
e}$, ${\bf d}$, ${\bf\gamma}$, ${\bf a}$ (of dimensions $N\times N$),
${\bf\alpha}$, ${\bf\beta}$ (of dimensions $N\times N\times n$), ${\bf c}$ (of
dimensions $N\times N\times n\times n$) and the vector $\vec{f}$ (of dimension
$N$), where $n$ is the spatial dimension of the problem ($n=1,2,3$). In index
notation, this equation reads as

$e_{ij}\frac{\partial^{2}u_{j}}{\partial t^{2}}+d_{ij}\frac{\partial
u_{j}}{\partial t}+\frac{\partial}{\partial
x_{l}}\left(-c_{ijkl}\frac{\partial u_{j}}{\partial
x_{k}}-\alpha_{ijl}u_{j}+\gamma_{il}\right)+\beta_{ijl}\frac{\partial
u_{j}}{\partial x_{l}}+a_{ij}u_{j}=f_{i},$ (36)

where $i,j=1.\ldots,N$ and $k,l=1,\ldots,n$.

The considered system, given by Eq. (21), contains two components
$u=(\tilde{h},\tilde{P})$ ($N=2$) and we use two equations, namely, Eqs. (21)
and (22) for ($n=1)$, corresponding to the solution depending on a single
spatial variable, $x$.

In the following, we list the non-vanishing coefficients (we omit the indexes
$k$ and $l$ for brevity and consider $x_{1}\equiv x$ since $k=l=1$):

  * •

Row $1$ ($i=1$) for Eq. (21)

$d_{11}=1,\quad c_{12}=\tilde{h}^{3}.$ (37)

  * •

Row $2$ ($i=2$) for Eq. (22)

$c_{21}=-1,\quad a_{21}=-{\rm Bo},\qquad a_{22}=1,\qquad f_{2}={\cal
S}{\tilde{\phi}}_{s}$ (38)

At the domain ends, we apply Dirichlet boundary conditions
$\tilde{h}=\tilde{h}_{p}$, where $\tilde{h}_{p}$ is the precursor film
thickness, and $\partial h/\partial x=0$. In order to achieve convergent
numerical simulations, we discretize the domain with steps $\Delta x=h_{p}$.

Declaration of Interests

The author report no conflict of interest.

## References

  * Altshuler & Manor (2015) Altshuler, G. & Manor, O. 2015 Spreading dynamics of a partially wetting water film atop a $\mathrm{MH}$z substrate vibration. Phys. Fluids 27, 102103.
  * Altshuler & Manor (2016) Altshuler, G. & Manor, O. 2016 Free films of a partially wetting liquid under the influence of a propagating mhz surface acoustic wave. Phys. Fluids 28, 072102.
  * Alzuaga et al. (2005) Alzuaga, S., Manceau, J.F. & Bastien, F. 2005 Motion of droplets on solid surface using acoustic radiation pressure. J. Sound Vib. 282, 151–162.
  * Amon et al. (2001) Amon, C.H., Murthy, J., Yao, S.C., Narumanchi, S., Wu, C.F. & Hsieh, C.C. 2001 MEMS-enabled thermal management of high-heat-flux devices EDIFICE: embedded droplet impingement for integrated cooling of electronics. Exp. Thermal Fluid Sci. 25, 231–242.
  * Arzt et al. (1967) Arzt, R.M., Salzmann, E. & Dransfeld, K. 1967 Elastic surface waves in quartz at 316 MHz. Appl. Phys. Lett. 10, 165–167.
  * Atencia & Beebe (2005) Atencia, J. & Beebe, D.J. 2005 Controlled microfluidic interfaces. Nature 437, 648–655.
  * Ballantine Jr et al. (1996) Ballantine Jr, D.S., White, R.M., Martin, S.J., Ricco, A.J., Zellers, E.T., Frye, G.C. & Wohltjen, H. 1996 Acoustic wave sensors: theory, design and physico-chemical applications. Elsevier.
  * Bar-Cohen et al. (2006) Bar-Cohen, A., Arik, M. & Ohadi, M. 2006 Direct liquid cooling of high flux micro and nano electronic components. Proc. IEEE 94, 1549–1570.
  * Biwersi et al. (2000) Biwersi, S., Manceau, J.F. & Bastien, F. 2000 Displacement of droplets and deformation of thin liquid layers using flexural vibrations of structures. Influence of acoustic radiation pressure. J. Acoust. Soc. Am. 107, 661–664.
  * Borgnis (2010) Borgnis, F.E. 2010 Acoustic radiation pressure of plane compressional waves. Rev. Mod. Phys. 25, 653–664.
  * Brown & Cox (2009) Brown, D. & Cox, A.J. 2009 Innovative uses of video analysis. The Physics Teacher 47, 145–150.
  * Brunet et al. (2010) Brunet, P., Baudoin, M., Matar, O.B. & Zoueshtiagh, F. 2010 Droplet displacements and oscillations induced by ultrasonic surface acoustic waves: A quantitative study. Phys. Rev. E 81, 036315.
  * Brunet et al. (2007) Brunet, P., Eggers, J. & Deegan, R.D. 2007 Vibration-induced climbing of drops. Phys. Rev. Lett. 99, 144501.
  * Brunet et al. (2009) Brunet, P., Eggers, J. & Deegan, R.D. 2009 Motion of a drop driven by substrate vibrations. Eur. Phys. J. Special Topics 166, 11.
  * Campbell (1970) Campbell, J.J. 1970 Propagation of surface waves at the boundary between a piezoelectric crystal and a fluid medium. IEEE Trans. Sonics Ultrason. 17, 71–76.
  * Campbell & Jones (1968) Campbell, J.J. & Jones, W.R. 1968 A method for estimating optimal crystal cuts and propagation directions for excitation of piezoelectric surface waves. IEEE Trans. Sonics Ultrason. 15, 209–217.
  * Craster & Matar (2009) Craster, R.V. & Matar, O.K. 2009 Dynamics and stability of thin liquid films. Rev. Mod. Phys. 81, 1131.
  * Diez & Kondic (2002) Diez, J.A. & Kondic, L. 2002 Computing three-dimensional thin film flows including contact lines. J. Comput. Phys. 183, 274.
  * Dussan V. & Chow (1983) Dussan V., E.B. & Chow, R.T.P. 1983 On the ability of drops or bubbles to stick to non-horizontal surfaces of solids. J. Fluid Mech. 137, 1–29.
  * Dussan V. & Davis (1974) Dussan V., Elizabeth B. & Davis, Stephen H. 1974 On the motion of a fluid-fluid interface along a solid surface. J. Fluid Mech. 65, 71–95.
  * Eckart (1948) Eckart, C. 1948 Vortices and streams caused by sound waves. Phys. Rev. 73, 68–76.
  * Fakhfouri et al. (2018) Fakhfouri, A., Devendran, C., Albrecht, T., Collins, D.J., Winkler, A., Schmidt, H. & Neild, A. 2018 Surface acoustic wave diffraction driven mechanisms in microfluidic systems. Lab Chip 18, 2214–2224.
  * Fendler (1996) Fendler, J.H. 1996 Self assembled nanostructured materials. Chem. Materials 8, 1616–1624.
  * Fletcher et al. (1974) Fletcher, L.S., Sernas, V. & Galowin, L.S. 1974 Evaporation from thin water films on horizontal tubes. Ind. Eng. Chem. Process Design Development 13, 265–269.
  * Guttenberg et al. (2004) Guttenberg, Z.V., Rathgeber, A, Keller, S, Rädler, J.O., Wixforth, A, Kostur, M, Schindler, M & Talkner, P 2004 Flow profiling of a surface-acoustic-wave nanopump. Phys. Rev. E 70, 056311.
  * Hasegawa et al. (2000) Hasegawa, T., Kido, T., Iizuka, T. & Matsuoka, C. 2000 A general theory of Rayleigh and Langevin radiation pressure. J. Acoust. Soc. Jap. E 21, 146–152.
  * Holly & Lemp (1971) Holly, F.J. & Lemp, M.A. 1971 Wettability and wetting of corneal epithelium. Experimental Eye Research 11, 239–250.
  * King (1934) King, L.V. 1934 On the acoustic radiation pressure on spheres. Proc. R. Soc. Lond. A 147, 212–240.
  * Kondic (2003) Kondic, L. 2003 Instability in the gravity driven flow of thin liquid films. SIAM Review 45, 95.
  * Lighthill (1978) Lighthill, J. 1978 Acoustic streaming. J. Sound Vib. 61, 391–418.
  * Longuet-Higgins (1953) Longuet-Higgins, M.S. 1953 Mass Transport in Water Waves. Philosophical Transactions of the Royal Society of London A 245, 535–581.
  * Manor et al. (2012) Manor, O., Yeo, L.Y. & Friend, J.R. 2012 The appearance of boundary layers and drift flows due to high-frequency surface waves. J. Fluid Mech. 707, 482–495.
  * Mhatre et al. (2016) Mhatre, S., Zigelman, A., Abezgauz, L. & Manor, O. 2016 Influence of a Propagating Megahertz Surface Acoustic Wave on the Pattern Deposition of Solute Mass off an Evaporating Solution. Langmuir 32, 9611–9618.
  * Morozov & Manor (2018) Morozov, M. & Manor, O. 2018 Vibration-driven mass transfer and dynamic wetting. Current Opinion in Colloid & Interface Science 36, 37–45.
  * Nagayama (1996) Nagayama, K. 1996 Two-dimensional self-assembly of colloids in thin liquid films. Colloids Surf. A .
  * Nyborg (1952) Nyborg, W.L. 1952 Acoustic streaming due to attenuated plane waves. J. Acoust. Soc. Am. 25, 1–8.
  * Oron et al. (1997) Oron, A., Davis, S. H. & Bankoff, S. G. 1997 Long-scale evolution of thin liquid films. Rev. Mod. Phys. 69, 931.
  * PrecisionFlow (2023) PrecisionFlow 2023 Speed of sound in common materials for use with ultrasonic flow meters. http://www.precisionflow.co.uk.
  * Rajendran et al. (2022) Rajendran, V.K., Jayakumar, S., Azharudeen, M. & Subramani, K. 2022 Theory of nonlinear acoustic forces acting on inhomogeneous fluids. J. Fluid Mech. 940, A32.
  * Rezk et al. (2012) Rezk, A.R., Manor, O., Friend, J.R. & Yeo, L.Y. 2012 Unique fingering instabilities and soliton-like wave propagation in thin acoustowetting films. Nat. Commun. 3, 1167.
  * Rezk et al. (2014) Rezk, A.R., Manor, O., Yeo, L.Y & Friend, J.R. 2014 Double flow reversal in thin liquid films driven by MHz order surface vibration. Proc. R. Soc. Lond. A. 470, 20130765.
  * Roditi (2023) Roditi 2023 Lithium niobate wafers - applications and specifications. https://www.roditi.com/SingleCrystal/LiNbO3/LiNbO3-Wafers.html.
  * Schlichting (1932) Schlichting, H. 1932 Calculation of even periodic barrier currents. Phys. Z. 33, 327–335.
  * Shiokawa et al. (1989) Shiokawa, S., Matsui, Y. & Ueda, T. 1989 Liquid streaming and droplet formation caused by leaky Rayleigh waves. Proc. IEEE Ultras. Symp. 1, 643–646.
  * Shiokawa et al. (1990) Shiokawa, S., Matsui, Y. & Ueda, T. 1990 Study on saw streaming and its application to fluid devices. Jap. J. Applied Physics 29, 137–139.
  * Stokes (1847) Stokes, G.G. 1847 On the theory of oscillatory waves. Camb. Trans. 8, 441–473.
  * Stone et al. (2004) Stone, H.A., Stroock, A.D. & Ajdari, A. 2004 Engineering flows in small devices. Ann. Rev. Fluid Mech. 36, 381–411.
  * Veldhuizen & Haagsman (2000) Veldhuizen, E.J.A. & Haagsman, H.P. 2000 Role of pulmonary surfactant components in surface film formation and dynamics. Biochimica et Biophysica Acta (BBA) - Biomembranes 1467, 255–270.
  * Wang et al. (2004) Wang, J.Z., Zheng, Z.H., Li, H.W., Huck, W.T.S. & Sirringhaus, H. 2004 Dewetting of conducting polymer inkjet droplets on patterned surfaces. Nature Materials 3, 171–176.
  * Whitesides (2006) Whitesides, G.M. 2006 The origins and the future of microfluidics. Nature 442, 368–373.
  * Wixforth (2003) Wixforth, A. 2003 Acoustically driven planar microfluidics. Superlattices and Microstructures 33, 389–396.
  * Wixforth et al. (2004) Wixforth, A., Strobl, C., Gauer, C., Toegl, A., Scriba, J. & von Guttenberg, Z. 2004 Acoustic manipulation of small droplets. Anal. Bioanal. Chem. 379, 982–991.
  * Yeo & Friend (2014) Yeo, L.Y. & Friend, J.R. 2014 Surface acoustic wave microfluidics. Ann. Rev. Fluid Mech. 46, 379–406.


# Engineering flat bands in twisted-bilayer graphene
away from the magic angle with chiral optical cavities

Cunyuan Jiang1,2,3 Matteo Baggioli1,2,3<EMAIL_ADDRESS>Qing-Dong
Jiang1,4,5<EMAIL_ADDRESS>1 School of Physics and Astronomy,
Shanghai Jiao Tong University, Shanghai 200240, China  2 Wilczek Quantum
Center, School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China  3 Shanghai Research Center for Quantum Sciences,
Shanghai 201315,China  4 Tsung-Dao Lee Institute, Shanghai Jiao Tong
University, Shanghai 200240, China  5 Shanghai Branch, Hefei National
Laboratory, Shanghai 201315, China

###### Abstract

Twisted bilayer graphene (TBG) is a recently discovered two-dimensional
superlattice structure which exhibits strongly-correlated quantum many-body
physics, including strange metallic behavior and unconventional
superconductivity. Most of TBG exotic properties are connected to the
emergence of a pair of isolated and topological flat electronic bands at the
so-called magic angle, $\theta\approx 1.05^{\circ}$, which are nevertheless
very fragile. In this work, we show that, by employing chiral optical
cavities, the topological flat bands can be stabilized away from the magic
angle in an interval of approximately $0.8^{\circ}<\theta<1.3^{\circ}$. As
highlighted by a simplified theoretical model, time reversal symmetry
breaking, induced by the chiral nature of the cavity, plays a fundamental role
in flattening the isolated bands and gapping out the rest of the spectrum. The
efficiency of the cavity is discussed as a function of the twisting angle, the
light-matter coupling and the optical cavity characteristic frequency. Our
results demonstrate the possibility of engineering flat bands in TBG using
optical devices, extending the onset of strongly-correlated topological
electronic phases in Moiré superlattices to a wider range in the twisting
angle.

Introduction – Controlling and engineering quantum phases of matter is a
central task in condensed matter physics. Inspired by the original discovery
of single-layer graphene Novoselov  _et al._ (2004), two-dimensional (2D)
materials have emerged as a versatile platform to realize strongly-correlated
physics in quantum many body systems Novoselov  _et al._ (2016). Recently,
unconventional superconductivity was discovered in twisted bilayer graphene
(TBG), a two dimensional superlattice where one layer of graphene is stacked
on top of another at a special magic twisting angle, i.e., $\theta\approx
1.05^{\circ}$ Cao  _et al._ (2018a); Yankowitz  _et al._ (2019); Lu  _et al._
(2019); Saito  _et al._ (2020). Galvanized by this breakthrough, several other
stacked two-dimensional systems that host exotic superconductivity, such as
twisted multilayer graphene, have been revealed Chen  _et al._ (2019a, b);
Park  _et al._ (2021); Zhang  _et al._ (2021); Park  _et al._ (2022a). While
the underneath physical mechanism of superconductivity in twisted 2D systems
is still under debate Po  _et al._ (2018); Isobe  _et al._ (2018); Wu  _et
al._ (2018); Xu and Balents (2018); Xie  _et al._ (2020); Song and Bernevig
(2022); Liu  _et al._ (2018); Roy and Juričić (2019), it is clear that the
isolated electronic flat band appearing at the magical angle plays an
essential role. Besides superconductivity, flat bands are also indispensable
for the emergence of strongly-correlated insulating states and the strange-
metal phase near the superconducting dome in the phase diagram of TBG, which
closely mimics that of cuprate superconductors Cao  _et al._ (2018b, 2020);
Liu  _et al._ (2021); Xie and MacDonald (2020); Koshino  _et al._ (2018);
Dodaro  _et al._ (2018); Yuan and Fu (2018); Thomson  _et al._ (2018); Pizarro
_et al._ (2019).

However, despite being a promising platform for studying strongly correlated
physics, the unavoidable and uncontrollable non-uniformity of the twist angle
across the sample, and the consequent difficulty in keeping the twist angle at
its magic value, prevented a wide realization of these phenomena Park  _et
al._ (2022b); Wilson  _et al._ (2020). More precisely, because the magical-
angle configuration is unstable, a little offset (around $\pm 0.1^{\circ}$) of
the twisting angle easily destroys all the emergent exotic properties of TBG.
In this regard, one of the most important challenges in the field is therefore
to achieve superconductivity at non-magic values of the twisting angle. To
achieve this final goal, it is desirable to realize a primary step, namely to
create and stabilize electronic flat bands in a wider range of the twisting
angle Bi  _et al._ (2019); Choi and Choi (2019); Liu  _et al._ (2020);
Sboychakov  _et al._ (2018); Wolf  _et al._ (2019).

In this Letter, we propose a new method to engineer stable flat bands at non-
magic angles by embedding twisted-bilayer graphene in a vacuum chiral cavity
(see top panel in Figure 1 for a cartoon of the setup). Using vacuum cavities
to control materials and molecules has emerged as a fruitful playground
connecting quantum optics to condensed matter and chemistry Hübener  _et al._
(2021); Schlawin  _et al._ (2022); Bloch  _et al._ (2022); Schäfer  _et al._
(2018). Vacuum cavity engineering is superior to the Floquet method, where
external electromagnetic radiation drives the system out of equilibrium, since
this second route inevitably heats up the system, destroying quantum coherence
and inducing transient phenomena away from thermal equilibrium states. In the
past, the usage of vacuum cavities has been proposed to design material
conductivity Rokaj  _et al._ (2022); Moddel  _et al._ (2021), unconventional
superconductivity Sentef  _et al._ (2018); Schlawin  _et al._ (2019); Curtis
_et al._ (2019); Thomas  _et al._ (2019), topological properties Espinosa-
Ortega  _et al._ (2014); Wang  _et al._ (2019); Appugliese  _et al._ (2022),
and even chemical reactivity Galego  _et al._ (2019, 2015). Some of these
proposals have been already successfully realized experimentally.

Figure 1: Top panel. A cartoon of the optical setup considered in this work.
Two skewed sheets of graphene are stacked on top of each other with a twisting
angle $\theta$ creating a characteristic Moiré pattern. They are then put
inside a chiral optical cavity with a light-matter coupling $g$ and a
characteristic frequency $\omega_{c}$. Bottom panel. The dominant Feynman
diagram describing light-matter interaction in the optical chiral cavity and
giving rise to the self-energy $\Sigma_{0}(\epsilon,\mathbf{q})$.

A fundamental property of vacuum chiral cavities is that time-reversal
symmetry is broken without the need of an external driving. Time-reversal
symmetry breaking is essential since, as we will see, quantum fluctuations
alone can not significantly influence the electronic bands in TBG. In single-
layer graphene, a band gap can be induced by quantum fluctuations in a chiral
cavity as well Wang  _et al._ (2019). However, the effect is too small to be
directly observed due to the large bandwidth. As we will demonstrate, the
situation is different in TBG near the magic angle, where the small bandwidth
enables time-reversal symmetry-broken quantum fluctuations to play a
significant role.

In recent years, a number of works have realized the vital impact of symmetry
breaking on quantum-fluctuations-related phenomena, such as anomalous Casimir
effects Butcher  _et al._ (2012); Jiang and Wilczek (2019a), topological gap
generation Espinosa-Ortega  _et al._ (2014); Wang  _et al._ (2019), and
selection of chiral molecules in chemical reactions Ke  _et al._ (2022); Riso
_et al._ (2022). A recent work by one of us Jiang and Wilczek (2019b)
highlighted the combined power of symmetry breaking and quantum fluctuations,
proving that symmetry breaking effects can be transmitted from a material to
its vicinity by vacuum quantum fluctuations. In this scenario, the vacuum in
proximity of a material with broken symmetries is referred to as its Quantum
Atmosphere.

In this Letter, we investigate the band renormalization of TBG in the time-
reversal symmetry broken quantum atmosphere provided by a chiral cavity. We
start from a faithful tight-binding model of TBG and calculate the one-loop
self-energy induced by the light-matter coupling. The bottom panel of Figure 1
displays the specific Feynman diagram considered. We find that, for
experimentally realizable values of the light-matter coupling and cavity
frequency, the topological flat bands in TBG can be stabilized away from the
magic angle in an interval of approximately $0.8^{\circ}<\theta<1.3^{\circ}$.
Our derivation and calculations can be directly generalized to other twisted
2D systems.

Setup and methods – To set the stage, we model the Hamiltonian of the combined
system, TBG and cavity, as follows:

$\displaystyle\hat{H}=\hat{H}_{\text{TBG}}(\mathbf{q}-e\boldsymbol{\hat{A}})+\hbar\omega_{c}\hat{a}^{\dagger}\hat{a},$
(1)

where $H_{\text{TBG}}(\mathbf{q})$ represents the TBG Hamiltonian in
reciprocal space, and $\omega_{c}$ is the cavity photonic mode frequency. TBG
and cavity photonic modes are coupled through Peierls substitution
$\mathbf{q}\mapsto\mathbf{q}-e\boldsymbol{\hat{A}}$, where
$\boldsymbol{\hat{A}}$ can be expressed in terms of photonic creation and
annihilation operators, i.e.,
$\hat{\mathbb{A}}=A_{0}\left(\boldsymbol{\varepsilon}^{*}\hat{a}^{\dagger}+\boldsymbol{\varepsilon}\hat{a}\right)$.
Here, $\boldsymbol{\varepsilon}$ is the polarization tensor of the cavity
photonic modes and $A_{0}=\sqrt{\frac{\hbar}{2\epsilon_{0}V\omega_{c}}}$ is
the mode amplitude in terms of the cavity volume $V$. We focus on chiral
cavities where the photonic polarization is given by
$\boldsymbol{\varepsilon}=\frac{1}{\sqrt{2}}\left(\mathbf{e_{x}}+i\mathbf{e_{y}}\right)$,
with $\mathbf{e_{x(y)}}$ the unit vector in the x(y)-direction. Our setup can
be straightforwardly generalized to the multi-mode case.

To be concrete, let us consider the effective tight-binding Hamiltonian
$H_{\text{TBG}}(\mathbf{q})$ Catarina  _et al._ (2019):

$\left(\begin{matrix}H_{1}(\mathbf{q})&T_{\mathbf{q}_{b}}&T_{\mathbf{q}_{tr}}&T_{\mathbf{q}_{tl}}&\cdots\\\
T_{\mathbf{q}_{b}}^{\dagger}&H_{2}(\mathbf{q}-\mathbf{q}_{b})&0&0&\cdots\\\
T_{\mathbf{q}_{tr}}^{\dagger}&0&H_{2}(\mathbf{q}-\mathbf{q}_{tr})&0&\cdots\\\
T_{\mathbf{q}_{tl}}^{\dagger}&0&0&H_{2}(\mathbf{q}-\mathbf{q}_{tl})\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{matrix}\right)$ (2)

where $\mathbf{q}$ is the wave-vector, and $H_{1,2}(\mathbf{q})$ indicate the
Hamiltonian of the top/bottom layer respectively. Moreover, we have defined:

$\displaystyle\mathbf{q}_{b}=\dfrac{1}{3}(\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m}),\quad\mathbf{q}_{tr}=\dfrac{1}{3}(\mathbf{b}_{1}^{m}+2\mathbf{b}_{2}^{m}),$
$\displaystyle\mathbf{q}_{tl}=\dfrac{1}{3}(-2\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m}),$
(3)

where $\mathbf{b}_{1}^{m}$ and $\mathbf{b}_{2}^{m}$ are Moiré reciprocal
vectors. The twisting angle $\theta$ is hidden in these vectors; see Fig.6.11
in Ref.Catarina  _et al._ (2019) for more details. The hopping matrix elements
are given by,

$\displaystyle T_{\mathbf{q}_{b}}=t\left(\begin{matrix}1&1\\\
1&1\end{matrix}\right),\quad
T_{\mathbf{q}_{tr}}=t\left(\begin{matrix}e^{i\phi}&1\\\
e^{-i\phi}&e^{i\phi}\end{matrix}\right),$ $\displaystyle
T_{\mathbf{q}_{tl}}=t\left(\begin{matrix}e^{-i\phi}&1\\\
e^{i\phi}&e^{-i\phi}\end{matrix}\right),$ (4)

where $\phi=2\pi/3$ and $t=0.11$ is the hopping parameter. For more details
about the TBG Hamiltonian we refer to the Supplementary Information (SI). Once
the effective Hamiltonian is known, the bare electron propagator can be
obtained using Wang  _et al._ (2019)

$G_{0}^{-1}(\epsilon,\textbf{q})=\left[\left(\epsilon+i0^{+}\right)\mathbf{I}-H_{\text{TBG}}(\textbf{q})\right]^{-1}.$
(5)

The full electron propagator, taking into account the interaction with the
vacuum cavity, can then be derived from the Dyson equation

$G^{-1}(\epsilon,\mathbf{q})=G_{0}^{-1}(\epsilon,\mathbf{q})-\Sigma_{0}(\epsilon,\mathbf{q}),$
(6)

where $\Sigma_{0}$ is the self-energy (see bottom panel of Fig.1) given by

$\Sigma_{0}(\epsilon,\mathbf{q})=-\dfrac{g^{2}}{\beta}\sum_{m=1}^{\infty}G_{0}(\epsilon+i\omega_{m},\mathbf{q})D_{0}(\omega_{m}).$
(7)

Here, $\omega_{m}=2\pi mk_{B}T$ is the $m$th Matsubara frequency and
$g=v_{F}eA_{0}$ denotes the electron-photon coupling strength with $v_{F}$ the
Fermi velocity of monolayer graphene and $e$ the electromagnetic coupling.
Finally, $D_{0}(\omega_{m})$ is the photon propagator given by

$D_{0}(\omega_{m})=\left(\begin{matrix}\dfrac{-1}{i\omega_{m}+\omega_{c}}&0&\cdots\\\
0&\dfrac{1}{i\omega_{m}-\omega_{c}}&\cdots\\\
\vdots&\vdots&\ddots\end{matrix}\right),$ (8)

with $\omega_{c}$ the cavity frequency. It should be noticed here that all
quantities $G_{0}$, $G$, $D_{0}$, $\Sigma_{0}$ are matrices with the same
dimension of the effective TBG Hamiltonian. With the full propagator in
Eq.(6), the spectral function, giving the renormalized electronic band
structure, can be calculated from

$A(\epsilon,\textbf{q})=-\frac{1}{\pi}\mathrm{Im}\mathrm{Tr}G(\epsilon,\textbf{q}).$
(9)

More details about the TBG Hamiltonian, the structure in reciprocal space and
the numerical methods employed can be found in the SI. In the rest of this
manuscript, we will be mainly interested in studying the electronic spectral
function $A(\epsilon,\mathbf{q})$ as a function of the twisting angle
$\theta$, the light-matter coupling $g$ and the cavity frequency $\omega_{c}$.
For convenience, we define the dimensionless coupling $\tilde{g}\equiv
g/k_{B}T$.

Figure 2: The spectral function $A(\epsilon,q)$ showing the band structure of
TBG with different light-matter coupling strength $g$ and twisting angle
$\theta$. (a): $\theta=0.8^{\circ},\,\tilde{g}=0$; (b):
$\theta=0.8^{\circ},\,\tilde{g}=2$; (c): $\theta=1.5^{\circ},\,\tilde{g}=0$;
(d): $\theta=1.5^{\circ},\,\tilde{g}=2$. The optical cavity characteristic
frequency is fixed to $\omega_{c}=0.3$eV. The $q$-path in reciprocal space is
chosen along the wave vector $\mathbf{q}_{b}$.

Electronic spectrum – It is well known that TBG exhibits a pair of topological
flat bands at the magic angle, $\theta\approx 1.05^{\circ}$ Bistritzer and
MacDonald (2011); Suárez Morell  _et al._ (2010); Andrei and MacDonald (2020),
which play a key role for the underlying strongly-correlated physics. Away
from that magic value, the flat bands, and most of the exotic and interesting
related properties of TBG, disappear. In panels (a) and (c) of Fig.2, we show
two examples of the electronic spectrum of TBG at respectively
$\theta=0.8^{\circ}$ and $\theta=1.5^{\circ}$. As already mentioned, no
isolated flat bands are present anymore in the electronic spectrum just moving
of $\pm 0.3^{\circ}$ from the magic angle. In other words, the flat bands are
very fragile and sensitive to the twisting angle. In panels (b) and (d) of
Fig.2, we show the same electronic spectra in presence of the chiral optical
cavity, with a coupling $\tilde{g}=2$ and a characteristic frequency
$\omega_{c}=0.3$eV. A pair of nearly-flat bands re-appear away from the magic
angle thanks to the coupling to the chiral cavity. Importantly, the two bands
are not anymore degenerate as their energy is shifted from the Fermi energy
and grows with the light-matter coupling $\tilde{g}$. This is a direct
consequence of the breaking of time reversal symmetry induced by the chiral
cavity. At the same time, the other bands, similarly to the case of the Dirac
cone in monolayer graphene (see SI), are also gapped away as a consequence of
the same symmetry breaking pattern.

Theoretical model – The emergence of the two isolated quasi-flat bands and
their energy splitting, shown in Fig.2, are intimately connected to the chiral
nature of the optical cavity, which plays a fundamental role in this regard.
To explicitly prove this statement we construct a simplified analytical model
which, as we will see, possesses all the minimal ingredients to describe our
setup. In order to model the effects of the chiral cavity on TBG, we consider
the following deformed hamiltonian

$H_{\text{TBG}+\tau}(\mathbf{q})=H_{\text{TBG}}(\mathbf{q})+\tau\left(\begin{matrix}\sigma_{z}&0&0&0&\cdots\\\
0&\sigma_{z}&0&0&\cdots\\\ 0&0&\sigma_{z}&0&\cdots\\\ 0&0&0&\sigma_{z}\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{matrix}\right),$ (10)

where the coupling $\tau$ parameterizes the breaking of time-reversal symmetry
on top of the original TBG Hamiltonian in Eq.(2). This is unlikely the most
general deformation which breaks time-reversal symmetry, but it will be
sufficient to qualitatively reproduce the numerical results displayed in Fig.2
and identify the main underlying physical principle behind them. By
diagonalizing the above Hamiltonian $H_{\text{TBG}+\tau}(\mathbf{q})$, the
band structure with broken time reversal symmetry can be obtained. The results
are shown in Fig.3 for different strength of the time-reversal symmetry
breaking $\tau$. As clearly demonstrated, the effects of $\tau$ is to gap away
the higher energy bands and create a pair of isolated quasi-flat bands wth
non-degenerate energy. At least at a qualitative level, these results are in
perfect agreement with the more realistic scenario of TBG in a chiral cavity
shown in Fig.2, where the light-matter coupling $\tilde{g}$ plays an analogous
role of the phenomenological parameter $\tau$ in Eq.(10). This simplified but
tractable analytical model highlights the fundamental role of time reversal
symmetry breaking, induced by the chiral cavity, in stabilizing the flat band
of TBG.

Figure 3: The band structure of TBG at $\theta=1.5^{\circ}$ obtained from the
deformed Hamiltonian in Eq.(10). Panels (a), (b), (c), (d) respectively
correspond to an increasing value for the time reversal breaking parameter
$\tau=0,0.1,0.2,0.3$.

Phase diagram and quasi-flat bands – In order to explore the effects of the
chiral cavity in more detail, we have performed an extensive analysis of the
band structure for different values of $\theta$ and $\tilde{g}$ covering a
wide range of the phase diagram around the magic-angle value. To give a
quantitative estimate of the flatness of the bands, we define the energy
bandwidth parameter

$\Delta\epsilon(g,\theta)\equiv\epsilon_{b}^{\text{max}}(g,\theta)-\epsilon_{b}^{\text{min}}(g,\theta)$
(11)

which quantifies the variation of the energy along the isolated nearly-flat
bands. As a reference, a completely flat band would correspond to
$\Delta\epsilon(g,\theta)=0$.

Figure 4: The value of the energy bandwidth $\Delta\epsilon(\tilde{g},\theta)$
for different $\tilde{g}$ and $\theta$. Darker color corresponds to a flatter
band. The black solid lines indicates a few constant $\Delta\epsilon$ values
as reference. The optical cavity frequency is set to $\omega_{c}=0.3$eV.

Our results are shown in Fig.4 in the interval
$0.8^{\circ}<\theta<1.5^{\circ}$ and $0<\tilde{g}<2.5$ for a fixed and
reasonable value of the cavity characteristic frequency $\omega_{c}=0.3$eV. By
increasing the value of the dimensionless coupling $\tilde{g}$, the energy
bandwidth becomes smaller and therefore the isolated bands can be flattened
away from the magic angle. As expected from simple intuition, the isolated
bands can be flattened more efficiently for angles which are closer to the
magic value. Interestingly, we find that it is easier to flatten the bands for
angles smaller than the magic one compared to angles larger than the latter.

In general, we observe that “quasi-flat” bands, with a variation of the energy
within $0.01$eV, can be easily tuned using reasonable values of the light-
matter coupling $\tilde{g}\sim\mathcal{O}(1)$ for angles of $\pm 0.3^{\circ}$
away from the magic value $\approx 1.05^{\circ}$.

Conclusions – In this Letter, we have revealed the possibility of extending
the onset of topological flat-bands in twisted bilayer graphene away from the
magic angle by using optical chiral cavities. We have demonstrated that the
effects of light-matter coupling can stabilize and flatten the topological
flat bands for a large range of the twisting angle without the need of fine-
tuning. Using physical values for the optical cavity frequency and the
strength of the light-matter coupling, we have estimated that quasi-flat bands
can be achieved at least in an interval of $0.8^{\circ}<\theta<1.3^{\circ}$.
From a theoretical point of view, taking advantage of a simplified analytical
model, we have identified the breaking of time-reversal symmetry as the
fundamental ingredient behind the achieved flattening.

One immediate future task is to verify whether all the interesting strongly-
correlated physics related to the topological flat bands survive in presence
of the chiral cavity or how that is modified. For example, it would be
interesting to understand further the effects of time-reversal symmetry
breaking on the emergent exotic superconductivity of TBG. A second and more
pressing point is to verify our theoretical predictions within an experimental
setup. Following our preliminary estimates, we conclude that the results shown
in this Letter might already be within experimental reach.

In general, we expect the combination of twistronics and photonics engineering
to become a powerful platform to study strongly-correlated electronic systems
and topological matter beyond the case of twisted bilayer graphene.

Acknowledgments – We would like to thank Dario Rosa for collaboration at an
early stage of this project. CJ and MB acknowledge the support of the Shanghai
Municipal Science and Technology Major Project (Grant No.2019SHZDZX01). MB
acknowledges the support of the sponsorship from the Yangyang Development
Fund. Q.-D. Jiang was sponsored by Pujiang Talent Program 21PJ1405400, Jiaoda
2030 program WH510363001-1, and the Innovation Program for Quajtum Science and
Technology Grant No.2021ZD0301900.

## References

  * Novoselov  _et al._ (2004) K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * Novoselov  _et al._ (2016) K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, Science 353, aac9439 (2016).
  * Cao  _et al._ (2018a) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018a).
  * Yankowitz  _et al._ (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. Young, and C. Dean, Science 363, eaav1910 (2019).
  * 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).
  * Saito  _et al._ (2020) Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. Young, Nature Physics 16 (2020), 10.1038/s41567-020-0928-3.
  * Chen  _et al._ (2019a) G. Chen, L. Jiang, S. Wu, B. Lyu, H. Li, B. L. Chittari, K. Watanabe, T. Taniguchi, Z. Shi, J. Jung, Y. Zhang, and F. Wang, Nature Physics 15 (2019a), 10.1038/s41567-018-0387-2.
  * Chen  _et al._ (2019b) G. Chen, A. Sharpe, P. Gallagher, I. Rosen, E. 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, 1 (2019b).
  * Park  _et al._ (2021) J. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 1 (2021).
  * Zhang  _et al._ (2021) X. Zhang, K.-T. Tsai, Z. Zhu, W. Ren, Y. Luo, S. Carr, M. Luskin, E. Kaxiras, and K. Wang, Phys. Rev. Lett. 127, 166802 (2021).
  * Park  _et al._ (2022a) J. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature Materials (2022a), 10.1038/s41563-022-01287-1.
  * Po  _et al._ (2018) H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Phys. Rev. X 8, 031089 (2018).
  * Isobe  _et al._ (2018) H. Isobe, N. F. Q. Yuan, and L. Fu, Phys. Rev. X 8, 041041 (2018).
  * Wu  _et al._ (2018) F. Wu, A. H. MacDonald, and I. Martin, Phys. Rev. Lett. 121, 257001 (2018).
  * Xu and Balents (2018) C. Xu and L. Balents, Phys. Rev. Lett. 121, 087001 (2018).
  * Xie  _et al._ (2020) F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Phys. Rev. Lett. 124, 167002 (2020).
  * Song and Bernevig (2022) Z.-D. Song and B. A. Bernevig, Phys. Rev. Lett. 129, 047601 (2022).
  * Liu  _et al._ (2018) C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Phys. Rev. Lett. 121, 217001 (2018).
  * Roy and Juričić (2019) B. Roy and V. Juričić, Phys. Rev. B 99, 121407 (2019).
  * Cao  _et al._ (2018b) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Nature 556, 80 (2018b).
  * Cao  _et al._ (2020) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020).
  * Liu  _et al._ (2021) X. Liu, Z. Wang, K. Watanabe, T. Taniguchi, O. Vafek, and J. Li, Science 371, 1261 (2021).
  * Xie and MacDonald (2020) M. Xie and A. H. MacDonald, Phys. Rev. Lett. 124, 097601 (2020).
  * 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).
  * Dodaro  _et al._ (2018) J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun, and C. Wang, Phys. Rev. B 98, 075154 (2018).
  * Yuan and Fu (2018) N. F. Q. Yuan and L. Fu, Phys. Rev. B 98, 045103 (2018).
  * Thomson  _et al._ (2018) A. Thomson, S. Chatterjee, S. Sachdev, and M. S. Scheurer, Phys. Rev. B 98, 075109 (2018).
  * Pizarro  _et al._ (2019) J. M. Pizarro, M. Calderón, and E. Bascones, Journal of Physics Communications 3, 035024 (2019).
  * Park  _et al._ (2022b) J. M. Park, Y. Cao, L.-Q. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature Materials 21, 877 (2022b).
  * Wilson  _et al._ (2020) J. H. Wilson, Y. Fu, S. Das Sarma, and J. H. Pixley, Phys. Rev. Res. 2, 023325 (2020).
  * Bi  _et al._ (2019) Z. Bi, N. F. Q. Yuan, and L. Fu, Phys. Rev. B 100, 035448 (2019).
  * Choi and Choi (2019) Y. W. Choi and H. J. Choi, Physical Review B 100, 201402 (2019).
  * Liu  _et al._ (2020) Y.-W. Liu, Y. Su, X.-F. Zhou, L.-J. Yin, C. Yan, S.-Y. Li, W. Yan, S. Han, Z.-Q. Fu, Y. Zhang, _et al._ , Physical Review Letters 125, 236102 (2020).
  * Sboychakov  _et al._ (2018) A. O. Sboychakov, A. V. Rozhkov, A. L. Rakhmanov, and F. Nori, Phys. Rev. Lett. 120, 266402 (2018).
  * Wolf  _et al._ (2019) T. M. R. Wolf, J. L. Lado, G. Blatter, and O. Zilberberg, Phys. Rev. Lett. 123, 096802 (2019).
  * Hübener  _et al._ (2021) H. Hübener, U. De Giovannini, C. Schäfer, J. Andberger, M. Ruggenthaler, J. Faist, and A. Rubio, Nature materials 20, 438 (2021).
  * Schlawin  _et al._ (2022) F. Schlawin, D. M. Kennes, and M. A. Sentef, Applied Physics Reviews 9, 011312 (2022).
  * Bloch  _et al._ (2022) J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, Nature 606, 41 (2022).
  * Schäfer  _et al._ (2018) C. Schäfer, M. Ruggenthaler, and A. Rubio, Physical Review A 98, 043801 (2018).
  * Rokaj  _et al._ (2022) V. Rokaj, M. Ruggenthaler, F. G. Eich, and A. Rubio, Physical Review Research 4, 013012 (2022).
  * Moddel  _et al._ (2021) G. Moddel, A. Weerakkody, D. Doroski, and D. Bartusiak, Physical Review Research 3, L022007 (2021).
  * Sentef  _et al._ (2018) M. A. Sentef, M. Ruggenthaler, and A. Rubio, Science advances 4, eaau6969 (2018).
  * Schlawin  _et al._ (2019) F. Schlawin, A. Cavalleri, and D. Jaksch, Physical review letters 122, 133602 (2019).
  * Curtis  _et al._ (2019) J. B. Curtis, Z. M. Raines, A. A. Allocca, M. Hafezi, and V. M. Galitski, Physical review letters 122, 167002 (2019).
  * Thomas  _et al._ (2019) A. Thomas, E. Devaux, K. Nagarajan, T. Chervy, M. Seidel, D. Hagenmüller, S. Schütz, J. Schachenmayer, C. Genet, G. Pupillo, _et al._ , arXiv preprint arXiv:1911.01459 (2019).
  * Espinosa-Ortega  _et al._ (2014) T. Espinosa-Ortega, O. Kyriienko, O. Kibis, and I. Shelykh, Physical Review A 89, 062115 (2014).
  * Wang  _et al._ (2019) X. Wang, E. Ronca, and M. A. Sentef, Phys. Rev. B 99, 235156 (2019).
  * Appugliese  _et al._ (2022) F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, and J. Faist, Science 375, 1030 (2022).
  * Galego  _et al._ (2019) J. Galego, C. Climent, F. J. Garcia-Vidal, and J. Feist, Physical Review X 9, 021057 (2019).
  * Galego  _et al._ (2015) J. Galego, F. J. Garcia-Vidal, and J. Feist, Physical Review X 5, 041022 (2015).
  * Butcher  _et al._ (2012) D. T. Butcher, S. Y. Buhmann, and S. Scheel, New Journal of Physics 14, 113013 (2012).
  * Jiang and Wilczek (2019a) Q.-D. Jiang and F. Wilczek, Physical Review B 99, 165402 (2019a).
  * Ke  _et al._ (2022) Y. Ke, Z. Song, and Q.-D. Jiang, arXiv preprint arXiv:2211.11132 (2022).
  * Riso  _et al._ (2022) R. R. Riso, L. Grazioli, E. Ronca, T. Giovannini, and H. Koch, arXiv preprint arXiv:2209.01987 (2022).
  * Jiang and Wilczek (2019b) Q.-D. Jiang and F. Wilczek, Physical Review B 99, 201104 (2019b).
  * Catarina  _et al._ (2019) G. Catarina, B. Amorim, E. V. Castro, J. M. V. P. Lopes, J. M. V. P. Lopes, and N. Peres, “Twisted bilayer graphene: Low-energy physics, electronic and optical properties,” in  _Handbook of Graphene Set_ (John Wiley & Sons, Ltd, 2019) Chap. 6, pp. 177–231.
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011).
  * 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).
  * Andrei and MacDonald (2020) E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020).

## Supplementary Information

###  SM1 The simpler case of graphene

In this section, as a warp-up exercise, we report the results related to the
effects of the chiral optical cavity on monolayer graphene (see Wang  _et al._
(2019)). As for the case of TBG presented in the main text, a simple
theoretical model based on the breaking of time reversal symmetry will play an
important role in understanding the outcomes of our computations.

For monolayer graphene, near the Dirac point, the effective Hamiltonian can be
expressed as,

$H(\mathbf{q})=v_{F}(\sigma_{x}q_{x}+\sigma_{y}q_{y}),$ (S1)

where $\mathbf{\sigma}$ are the Pauli matrices. The Hamiltonian gives a linear
dispersion relation $\epsilon(\mathbf{q})=\pm v_{F}\sqrt{q_{x}^{2}+q_{y}^{2}}$
as shown in the left panel of Fig.S1. Time reversal symmetry breaking can be
modelled by allowing for a new term in the Hamiltonian $\propto\tau$, which
becomes

$H(\mathbf{q})=v_{F}(\sigma_{x}q_{x}+\sigma_{y}q_{y})+\tau\sigma_{z}.$ (S2)

Solving the eigenvalue equation for the above Hamiltonian, the dispersion
relation under broken time reversal symmetry reads $\epsilon(\mathbf{q})=\pm
v_{F}\sqrt{q_{x}^{2}+q_{y}^{2}+\tau^{2}}$. The resulting band structure
acquire a band gap $\Delta=2v_{F}\tau$ at the Dirac point, as shown in the
right panel of Fig.S1.

In Fig.S1 we show the band structure of graphene given by the model based on
Eq.(S1) and the one obtained numerically by coupling graphene to the chiral
cavity, as explained in the main text for TBG. The two results match very well
and show that the chiral cavity opens a gap at the Dirac point as well, as a
natural consequence of time reversal symmetry breaking.

Figure S1: The band structure near the Dirac point for monolayer graphene with
time reversal symmetry (left panel) and broken time reversal symmetry (right
panel). Units are arbitrary. The orange solid line is the result given by the
simple theoretical model based on Eq.(S1). The color map represents the result
of coupling monolayer graphene to a chiral optical cavity and computing
numerically the spectral function $A(\epsilon,\mathbf{q})$ as outlined in the
main text for TBG.

###  SM2 The role of the cavity frequency

In order to complete our analysis, we now discuss briefly the role of the
characteristic frequency of the chiral optical cavity $\omega_{c}$. In the
main text, we have fixed $\omega_{c}=0.3$eV, which is a reasonable and
experimentally realizable value. In Fig.S2, we study how the electronic
spectrum is modified for larger values of such a frequency. Panels (a) and (b)
of Fig.S2 have to be compared with panel (d) of Fig.2 in the main text, which
shows the same results for smaller $\omega_{c}$. We observe that increasing
the cavity frequency the electronic bands becomes broader and the quasi-flat
bands get closer to the other gapped bands at higher energy. Eventually, as a
consequence of this strong broadening, the quasi-flat bands are not isolated
anymore. This indicates that a large optical cavity frequency is detrimental
in engineering isolated topological flat bands.

Figure S2: The spectral function $A(\epsilon,\mathbf{q})$ showing the band
structure of TBG in an optical chiral cavity for $g=2$ and
$\theta=1.5^{\circ}$. Panels (a) and (b) correspond respectively to
$\omega_{c}=0.8$eV, $\omega_{c}=1.5$eV. These figures must be compared with
panel (d) in Fig.2 which shows the same spectrum for $\omega_{c}=0.3$eV.

###  SM3 Efficiency of the cavity

In order to discuss the effect of the chiral cavity in more detail, we can
define the “flattening-efficiency” as,

$\eta(g,\theta)=1-\dfrac{\Delta\epsilon(g,\theta)}{\Delta\epsilon(0,\theta)},$
(S3)

where $\Delta\epsilon(g,\theta)$ is the bandwidth parameter discussed in the
main text. From this definition it can be found that $\eta$ takes values in
the range from $0$ to $1$. The minimum value $\eta=0$ means that the chiral
cavity has no effect on the flatness of the isolated bands. On the contrary,
the maximum value $\eta=1$ correspond to a perfect flattening of the band
after having introduce the cavity, $\Delta\epsilon(g,\theta)=0$. More in
general, a large value of $\eta$ implies that the bands are strongly flattened
by the chiral cavity with respect to their structure at $g=0$.

Fig.S3 shows the flattening parameter $\eta$ as a function of the twisting
angle $\theta$ and the dimensionless coupling $\tilde{g}$ across the phase
diagram, around the magic angle value $\theta=1.1^{\circ}$. We observe that,
except for small $\tilde{g}$ or near the magic angle, the chiral cavity can
always greatly improve the flatness of band. The results in Fig.S3 are
consistent with those reported in the main text in Fig.4.

Figure S3: The flattening efficiency parameter $\eta$ as a function of the
dimensionless light-matter coupling $\tilde{g}$ and the twisting angle
$\theta$. From blue to red, the efficiency grows larger. Because the band is
already flat near the magic angle, the efficiency drops to zero in there.

###  SM4 Numerical methods and details of the setup

The Hamiltonian of monolayer graphene can be obtained by considering the
hopping of electrons from one carbon atom to its three nearest neighboring
atoms. The electronic band structure shows a characteristic linear dispersion
relation around Fermi energy. The low energy Hamiltonian reads Catarina  _et
al._ (2019)

$H_{1}(\mathbf{q})=v_{F}\left(\begin{matrix}0&q_{x}-iq_{y}\\\
q_{x}+iq_{y}&0\end{matrix}\right),$ (S4)

where $v_{F}$ is the Fermi velocity which is around
$5.944eV\cdot\textup{\AA}$.

However, the situation in twisted bilayer graphene (TBG) is more complicated
because the unit cell of the Moiré pattern is much larger than that of
graphene. Hence, the Moiré Brillouin zone is quite small in reciprocal space
and therefore one should consider many hopping channels between a group of
Moiré Brillouin zones in the two layers. More details on this point can be
found in Catarina  _et al._ (2019). Below is the effective low energy
Hamiltonian for TBG used in this manuscript:

$H_{TBG}(\mathbf{q})=\left(\begin{matrix}H_{1}(\mathbf{q})&T_{\mathbf{q}_{b}}&T_{\mathbf{q}_{tr}}&T_{\mathbf{q}_{tl}}&\cdots\\\
T_{\mathbf{q}_{b}}^{\dagger}&H_{2}(\mathbf{q}-\mathbf{q}_{b})&0&0&\cdots\\\
T_{\mathbf{q}_{tr}}^{\dagger}&0&H_{2}(\mathbf{q}-\mathbf{q}_{tr})&0&\cdots\\\
T_{\mathbf{q}_{tl}}^{\dagger}&0&0&H_{2}(\mathbf{q}-\mathbf{q}_{tl})\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{matrix}\right).$ (S5)

In this effective Hamiltonian, $H_{1,2}(\mathbf{q})$ indicate the Hamiltonian
of the top/bottom layer. Moreover,
$\mathbf{q}_{b}=\dfrac{1}{3}(\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m})$,
$\mathbf{q}_{tr}=\dfrac{1}{3}(\mathbf{b}_{1}^{m}+2\mathbf{b}_{2}^{m})$,
$\mathbf{q}_{tl}=\dfrac{1}{3}(-2\mathbf{b}_{1}^{m}-\mathbf{b}_{2}^{m})$, where
$\mathbf{b}_{1}^{m}$ and $\mathbf{b}_{2}^{m}$ are Moirè reciprocal vectors.
The twisting angle $\theta$ is hidden in these vectors. The hopping matrix
elements are given by,

$T_{\mathbf{q}_{b}}=t\left(\begin{matrix}1&1\\\ 1&1\end{matrix}\right),\quad
T_{\mathbf{q}_{tr}}=t\left(\begin{matrix}e^{i\phi}&1\\\
e^{-i\phi}&e^{i\phi}\end{matrix}\right),\quad
T_{\mathbf{q}_{tl}}=t\left(\begin{matrix}e^{-i\phi}&1\\\
e^{i\phi}&e^{-i\phi}\end{matrix}\right),$ (S6)

where $\phi=2\pi/3$ and $t=0.11$ is the hopping parameter.

It should be noticed here that the Hamiltonian matrix of TBG shown in Eq.(S5)
is not of dimension $8\times 8$. Indeed, one should extend the matrix Eq.(S5)
by taking into account more bases to obtain the flat band at the magic angle.
For example, in the basis $|l,(m,n)\rangle$, $l=1,2$ indicates the top/bottom
layer, while $m,n=-2,-1,0,1,2$ indicates the combination of two Moire
reciprocal vectors, namely $\langle
l,(m,n)|l,(m,n)\rangle=H_{l}(\mathbf{q}+(l-1)\mathbf{q}_{b}+m\mathbf{b}_{1}^{m}+n\mathbf{b}_{2}^{m})$.
According to this formula, the Hamiltonian matrix should be $100\times 100$
dimesional.

As outlined in the main text, in principle the self-energy is given by an
infinite sum over the Matsubara frequencies as

$\Sigma_{0}(\epsilon,\mathbf{q})=-\dfrac{g^{2}}{\beta}\sum_{m=1}^{\infty}G_{0}(\epsilon+i\omega_{m},\mathbf{q})D_{0}(\omega_{m}),$
(S7)

where $\omega_{m}=2\pi mk_{B}T$ is the Matsubara frequency and $g$ the
coupling parameter between photon and electron. Nevertheless, from a
computational perspective, the Matsubara frequency summation in Eq.(S7) is
performed by introducing a finite UV cutoff $N$

$\Sigma_{0}(\epsilon,\mathbf{q})=-\dfrac{g^{2}}{\beta}\sum_{m=1}^{N}G_{0}(\epsilon+i\omega_{m},\mathbf{q})D_{0}(\omega_{m}).$
(S8)

The computations shown in the main text are for $N=100$. To make sure that the
sum converged for the value of the cutoff used in the main text, we have
computed the same spectral function with different cutoffs as shown in Fig.S4.
The results show only little difference upon changing the cutoff from $N=50$
to $N=150$ confirming that the numerical method is reliable.

Finally, let us emphasize that the propagator used in the main text is just
the lowest order self-energy without any self consistent calculation. To
obtain more accurate results, one should obtain the two solutions, $G[\Sigma]$
and $\Sigma[G]$, consistently by iteration algorithm. However, the lowest
order self energy is enough to grasp the main features.

Figure S4: The results for the self-energy Eq.(S8) using different values of
the cutoff. From left to right, $N=50$, $N=100$, $N=150$ while
$\theta=1.5^{\circ}$, $g=2$. This figure emphasizes the stability of our
numerical routine and its convergence as a function of the cutoff $N$.