# Dynamics of packed swarms: time-displaced correlators of two dimensional
incompressible flocks

Leiming Chen<EMAIL_ADDRESS>School of Material Science and Physics,
China University of Mining and Technology, Xuzhou Jiangsu, 221116, P. R. China
Chiu Fan Lee<EMAIL_ADDRESS>Department of Bioengineering, Imperial
College London, South Kensington Campus, London SW7 2AZ, U.K.  Ananyo Maitra
<EMAIL_ADDRESS>Laboratoire de Physique Théorique et Modélisation, CNRS
UMR 8089, CY Cergy Paris Université, F-95032 Cergy-Pontoise Cedex, France
Sorbonne Université and CNRS, Laboratoire Jean Perrin, F-75005, Paris, France
John Toner<EMAIL_ADDRESS>Department of Physics and Institute of Theoretical
Science, University of Oregon, Eugene, OR $97403^{1}$ Max Planck Institute for
the Physics of Complex Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany
School of Material Science and Physics, China University of Mining and
Technology, Xuzhou Jiangsu, 221116, P. R. China  Department of Bioengineering,
Imperial College London, South Kensington Campus, London SW7 2AZ, U.K.
Laboratoire de Physique Théorique et Modélisation, CNRS UMR 8089, CY Cergy
Paris Université, F-95032 Cergy-Pontoise Cedex, France  Sorbonne Université
and CNRS, Laboratoire Jean Perrin, F-75005, Paris, France  Department of
Physics and Institute of Theoretical Science, University of Oregon, Eugene, OR
$97403^{1}$ Max Planck Institute for the Physics of Complex Systems,
Nöthnitzer Str. 38, 01187 Dresden, Germany

###### Abstract

We analytically calculate the scaling exponents of a two-dimensional KPZ-like
system: coherently moving incompressible polar active fluids. Using three
different renormalization group approximation schemes, we obtain values for
the “roughness” exponent $\chi$ and anisotropy exponent $\zeta$ that are
extremely near the known exact results. This implies our prediction for the
previously completely unknown dynamic exponent $z$ is quantitatively accurate.

The Mermin-Wagner theorem MW  states that in equilibrium, two-dimensional
systems with finite-range interactions cannot spontaneously break a continuous
symmetry at any finite temperature. One of the earliest and most surprising
results in the field of active matter was that there is no analogue of the
Mermin-Wagner theorem in nonequilibrium matter. In particular, it was
_explicitly_ demonstrated that active polar units moving on a two-dimensional
frictional substrate and with _purely_ short-range interactions _can_
spontaneously break continuous rotation symmetry in two dimensions and form
long-range ordered flocks Vicsek  _et al._ (1995); Toner and Tu (1995, 1998);
toner_pre12 ; Chate_DADAM_rev . Such flocks have a _non-zero_ mean speed
$\langle{\bf v}\rangle$ even at finite noise strengths.

Despite the fact that the _existence_ of long-range-ordered two-dimensional
flocks has been demonstrated analytically, determining their scaling behavior
analytically has proved much more challenging. One class of systems that _has_
proved amenable to analytical treatment is _incompressible_ flocks inc_trans ;
Chen  _et al._ (2016); Chen ; Ano_pol ; CLMT-PRL1 ; CLMT-PRE1 .

Incompressibility can arise in many ways. One way is to make the density
extremely high. In this limit the effective compressibility of the flockers
vanishes, with any departure from the mean density being severely penalised
Chen  _et al._ (2016); Ano_pol . More accurately, as the compressibility goes
to $0$, the lengthscale up to which the system is effectively incompressible
diverges. A flock formed by a suspension of active swimmers in an
incompressible fluid in a narrow channel also inherits the incompressibility
of the fluid and is therefore incompressible. This is technically distinct
from the model we will consider here due to the presence of an extra conserved
quantity – the number of swimmers. However, if the number of active units is
_not_ conserved, for instance, due to birth and death Toner_Malthus ; Chen
_et al._ (2020a, b), the dynamics of such a system does belong to the
“incompressible flock” universality class we consider here.

Our understanding of the equal-time behavior of incompressible two dimensional
flocks is quite complete. Although they are non-equilibrium systems, their
hydrodynamic equations prove to be equivalent, once terms that are
“irrelevant” in the renormalization group sense are dropped, to those of an
equilibrium magnetic system Chen  _et al._ (2016); Kashuba  with long-ranged
interactions. Because of those long-ranged interactions, the Mermin-Wagner
theorem does not apply to these magnetic systems. The partition function for
this equilibrium system can then be further mapped Chen  _et al._ (2016) onto
that for a two-dimensional smectic, whose equal-time scaling laws are known
exactly via a further mapping Golub  onto the $1+1$-dimensional KPZ equation
KPZ . This analysis Chen  _et al._ (2016) gives the scaling law for the equal
time fluctuations $\mathbf{u}(\mathbf{r},t)$ of the local active fluid
velocity $\mathbf{v}(\mathbf{r},t)$ about its mean value
$\langle\mathbf{v}\rangle\equiv v_{0}\hat{\bf x}$, where we’ve defined our
coordinate system so that ${\bf\hat{x}}$ is along the mean velocity
spontaneously chosen by the system. Specifically

$\langle\mathbf{u}(\mathbf{r},t)\cdot\mathbf{u}(\mathbf{0},t)\rangle=|x|^{2\chi}{\cal
G}_{{}_{ET}}\left({|y|\over|x|^{\zeta}}\right)\,,$ (1)

where “ET” stands for “equal-time”. The exponents $\chi$ and $\zeta$ were
determined exactly by the aforementioned mappingsChen  _et al._ (2016) to be
$\chi=-1/2$ and $\zeta=3/2$.

Unfortunately, since this analysis was based entirely on the partition
function for the equivalent equilibrium model, no information about the
dynamics of the system can be obtained by this exact mapping. Most of this
missing information can be encoded in a single additional universal exponent,
namely, the “dynamical” exponent $z$. This can be defined by considering the
unequal\- time correlation function of the velocity fluctuations, which obeys
the scaling law

$\displaystyle\langle\mathbf{u}(\mathbf{r},t)\cdot\mathbf{u}(\mathbf{0},0)\rangle$
(2) $\displaystyle=$ $\displaystyle|x-v_{b}t|^{2\chi}{\cal
G}_{{}_{UET}}\left({|y|\over|x-v_{b}t|^{\zeta}},{|t|\over|x-v_{b}t|^{z}}\right)\,,$

and the subscript “UET” now means “unequal-time”. The ”boost velocity $v_{b}$
is a phenomenological parameter of our model.

In this Letter, we obtain the dynamical exponent $z$ for the polar phase of
two-dimensional, incompressible flocks using three different dynamic
renormalization group (DRG) schemes: two different uncontrolled calculations
exactly in two dimensions, and a $d=(d_{c}-\epsilon)$-expansion. We get

$z=1.65\pm 0.06\,,$ (3)

where the error bars correspond to the averages over the three DRG schemes. We
also calculate the the roughness exponent $\chi$ and the anisotropic exponent
$\xi$ in all three schemes.

As mentioned earlier, we already know the static exponents $\chi$ and $\zeta$
exactly, so our purpose in performing this DRG is only to calculate $z$.
However, our knowledge of the exact values of $\chi$ and $\zeta$ provides us
with a useful check on the quantitative accuracy of our DRG calculations,
because we can compare the approximate values of $\chi$ and $\zeta$ that we
get from those schemes with the known exact values. And they prove to be very
close; indeed, one of our uncontrolled approximations reproduces the known
exact values $\chi=-1/2$ and $\zeta=3/2$. We therefore believe our prediction
(3) is very accurate quantitatively.

Inserting (3) and the exact value for $\chi$ into (2), and considering the
limit $\mathbf{r}\to\mathbf{0}$, we obtain the temporal part of the velocity
correlation:

$\langle\mathbf{u}(\mathbf{0},t)\cdot\mathbf{u}(\mathbf{0},0)\rangle\sim
At^{2\chi\over z}=At^{-0.61\pm 0.02}\,,$ (4)

where $A$ is some non-universal constant.

We’ll now present the drivation of the above results. This begins with the
hydrodynamic equation of motion (EOM) of incompressible polar active fluids
can be constructed purely based on symmetry considerations, as has been
discussed in Chen  _et al._ (2016), and is:

$\partial_{t}\mathbf{v}+\lambda(\mathbf{v}\cdot{\bf\nabla})\mathbf{v}=-{\bf\nabla}\Pi-(\mathbf{v}\cdot{\bf\nabla}\Pi_{1})\mathbf{v}+U(|\mathbf{v}|)\mathbf{v}+\mu_{1}\nabla^{2}\mathbf{v}+\mu_{2}(\mathbf{v}\cdot{\bf\nabla})^{2}\mathbf{v}+\mathbf{f}(\mathbf{r},t)\
,$ (5)

where the “pressure” $\Pi$ is a Lagrange multiplier that enforces the
incompressibility constraint: ${\bf\nabla}\cdot\mathbf{v}=0$, the term
involving $U(|\mathbf{v}|)$—a smooth, analytic function of
$|\mathbf{v}|$—ensures that there is an “ordered” phase in which $\mathbf{v}$
has a nonzero mean magnitude $v_{0}$. That is, we assume that there is a
regime in parameter space of $U(|\mathbf{v}|)$ where $U(|\mathbf{v}|)>0$ for
$|\mathbf{v}|<v_{0}$, $U(|\mathbf{v}|)=0$ for $|\mathbf{v}|=v_{0}$, and
$U(|\mathbf{v}|)<0$ for $|\mathbf{v}|>v_{0}$. The “anisotropic pressure”
$\Pi_{1}$ is another generic function of $|\mathbf{v}|$. Finally,
$\mathbf{f}(\mathbf{r})$ is a zero-mean, Gaussian white noise with the
correlation

$\langle
f_{i}(\mathbf{r},t)f_{j}(\mathbf{r}^{\prime},t^{\prime})\rangle=2D\delta_{ij}\delta^{2}(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})\,,$
(6)

where the indices $i,j$ enumerate the spatial coordinates.

Expanding (5) about an ordered state using ${\bf v}=v_{0}\hat{{\bf x}}+{\bf
u}$ and again retaining only relevant terms (i.e., the terms that are
important in the limit of large time and length scales), we obtain the EOM
governing $\mathbf{u}$

$\displaystyle\partial_{t}u_{i}$ $\displaystyle=$
$\displaystyle-\partial_{i}\Pi+\mu_{\perp}\partial_{y}^{2}u_{i}+\mu_{x}\partial_{x}^{2}u_{i}$
(7) $\displaystyle-\alpha\left(u_{x}+{u_{y}^{2}\over
2v_{0}}\right)(\delta_{ix}+{u_{y}\over v_{0}}\delta_{iy})+f_{i}\,,$

where $\alpha\equiv-v_{0}\left({\rm d}U\over{\rm
d}|\mathbf{v}|\right)_{|\mathbf{v}|=v_{0}}$, $\mu_{{}_{\perp}}\equiv\mu_{1}$,
and $\mu_{x}\equiv\mu_{1}+\mu_{2}v_{0}^{2}$. As discussed in Chen  _et al._
(2016); Toner_Malthus , we have performed a Galilean transformation to a
reference frame ${\bf r}^{\prime}$ which moves with respect to our original
reference frame in the direction of mean flock motion at a speed
$v_{b}=\lambda v_{0}$; that is, ${\bf r}^{\prime}\equiv{\bf r}-\lambda
v_{0}t\hat{{\bf x}}$, or, equivalently, $x^{\prime}=x-v_{b}t$. For simplicity,
in (7) we have dropped the prime in $x$. Note that $u$ is also subject to the
incompressibility constraint ${\bf\nabla}\cdot\mathbf{u}=0$ inherited from
${\bf\nabla}\cdot\mathbf{v}=0$. The somewhat non-trivial power-counting that
leads to (7) is detailed in Chen  _et al._ (2016).

We now perform a DRG analysis to obtain the dynamic exponent for the
incompressible flock. Fourier transforming (7) and acting on the projection
operator $P_{yi}(\mathbf{q})$ to eliminate the pressure term foot1 , we obtain

$-{\rm i}\omega
u_{y}=P_{yx}\left(\mathbf{q}\right)\mathcal{F}_{\tilde{\mathbf{q}}}\left[-\alpha\left(u_{x}+{u_{y}^{2}\over
2}\right)\right]+P_{yy}(\mathbf{q})\mathcal{F}_{\tilde{\mathbf{q}}}\left[-\alpha\left(u_{x}+{u_{y}^{2}\over
2}\right)u_{y}+\mu_{x}\partial_{x}^{2}u_{y}+f_{y}\right]\,,$ (8)

where we have eliminated $v_{0}$ by absorbing it into $\mathbf{u}$ (i.e.,
$\mathbf{u}\to v_{0}\mathbf{u}$), and the symbol
$\mathcal{F}_{\tilde{\mathbf{q}}}$ represents the $\tilde{\mathbf{q}}$th
Fourier component.

Now we implement the standard DRG procedure forster_pra77  on EOM (8). First
we decompose the field $\mathbf{u}$ into “slow” and “fast” parts
$u^{<}(\mathbf{q})$ and $u^{>}(\mathbf{q})$, which are supported at small
$\mathbf{q}$’s and large $\mathbf{q}$’s, respectively. (We will eventually
take the large $\mathbf{q}$’s to lie in an infinitesimally thin shell along
the outer edge of the Brillouin zone).

Next we average out $u^{>}(\mathbf{q})$ to get the effective EOM for
$u^{<}(\mathbf{q})$. In this step the various coefficients in (8) are
renormalized and their corrections can be represented by graphs (i.e., Feynman
diagrams); these corrections are therefore called graphical corrections.
Finally we rescale time, lengths, and fields as

$\displaystyle t\to te^{z{\rm d}\ell},~{}~{}x\to xe^{{\rm d}\ell},~{}~{}y\to
ye^{\zeta{\rm d}\ell},$ (9a) $\displaystyle u_{y}\to u_{y}e^{\chi{\rm
d}\ell},~{}~{}u_{x}\to u_{x}e^{\left(\chi+1-\zeta\right){\rm d}\ell}\,,$ (9b)

to restore the Brillouin zone to its original size. Note that the form of the
rescaling in $u_{x}$ is imposed by the incompressibility condition.

This whole process can be repeated iteratively, which leads to recursion
relations for the three coefficients $\alpha$, $\mu_{x}$, and $D$ in (8) [Note
that $D$ is hidden in the correlations of the noise $\mathbf{f}$ (6)].

In deriving these recursion relations, we use three distinct DRG schemes. Two
of these are one-loop order DRG calculations in exactly two dimensions, which
are uncontrolled approximations since there is no limit in which they become
exact. The third scheme is an $\epsilon=5/2-d$ expansion to $O(\epsilon)$, in
which we analytically continue our calculation to $d>2$ by treating the $y$
direction as $d-1$ dimensional. This is a controlled approximation, since it
formally becomes exact in the limit $\epsilon\to 0$. Obviously, it is also
only approximate in $d=2$. The first of the two uncontrolled approximations
uses a Brillouin zone (BZ) $-\infty<q_{y}<\infty$, $-\Lambda<q_{x}<\Lambda$,
where $\Lambda$ is the ultraviolet cut-off, and the second uses the BZ
$-\infty<q_{x}<\infty$, $-\Lambda<q_{y}<\Lambda$. The $\epsilon$-expansion
result is independent of the shape of the Brillouin zone. All three approaches
yield values of the dynamical exponent that are fairly close to each other.
Furthermore, all three obtain values of the already exactly known exponents
$\chi$ and $\zeta$ which prove to be very close to those known exact values;
indeed, the first uncontrolled approximation we use yields the exact values of
the $\chi$ and $\zeta$.

Crucially, in all three DRG calculations, we make use the symmetry properties
of the EOM for ${\bf u}$ to make two important simplifications. First, we
notice that rotation invariance of our hydrodynamic EOM is ensured by choosing
the values of $\chi$ and $\zeta$ to keep all three “$\alpha$”s appearing in
(8) the same under rescaling, i.e.,

$\chi=\zeta-1\,.$ (10)

This simplification reduces total number of the recursion relations to three.
Specifically,

$\displaystyle\frac{{\rm d}\ln\alpha}{{\rm d}\ell}$ $\displaystyle=$
$\displaystyle z-2\zeta+2+\eta_{\alpha}\,,$ (11) $\displaystyle\frac{{\rm
d}\ln\mu_{x}}{{\rm d}\ell}$ $\displaystyle=$ $\displaystyle z-2+\eta_{\mu}\,,$
(12)

and

$\frac{{\rm d}\ln D}{{\rm d}\ell}=z-3+\zeta+\eta_{D}$ (13)

for the uncontrolled approximations, or

$\frac{{\rm d}\ln D}{{\rm d}\ell}=z-3+(3-d)\zeta+\eta_{D}$ (14)

for the $\epsilon$ expansion, where $\eta_{\alpha,\,u,D}$ denote graphical
corrections.

Next, the fact that the dynamics of ${\bf u}$ obeys detailed balance
introduces another simplification. This becomes clear if we formally introduce
a friction coefficient in the dynamics for ${\bf u}$;
$\Gamma\partial_{t}u_{i}=-\delta H/\delta u_{i}+f_{i}$, which is just the
time-dependent-Ginzberg-Landau model, or model $A$ PCMP ; Hohenberg , with
$H({\mathbf{M}})$ the Hamiltonian for a ordered divergence-free two-
dimensional magnet expanded around its minimum at nonzero magnetization, where
${\mathbf{M}}=(v_{0}+u_{x})\hat{x}+u_{y}\hat{y}$. $\Gamma=1$ in (8) but it
doesn’t retain that value under renormalization. This would appear as a
renormalization of the coefficient of $-i\omega u_{y}$; in general, this is an
independent quantity. However, here detailed balance implies that the ratio
$D/\Gamma$ _cannot_ change under renormalization i.e., $\Gamma$ and $D$ must
renormalize in the same way. This implies that the correction to $\ln\Gamma$
are the same as the correction to $\ln D$. If we denote the direct graphical
correction to the annealed noise $\eta^{dir}_{D}$ and the direct graphical
correction to $\Gamma$ as $\eta_{\omega}$, the argument just given implies
$\eta_{\omega}=\eta^{dir}_{D}$. To avoid retaining the friction coefficient as
another extra parameter, we divide both sides of the renormalized EOM by the
coefficient of $-i\omega u_{y}$ to fix the coefficient of $-i\omega u_{y}$ at
$1$. This effectively introduces an additional correction $-\eta_{D}$ to both
$\alpha$ and $\mu$, and $-2\eta_{D}$ to $D$. (The factor of $2$ arises because
$D$ is proportional to the correlation of two noises).

Therefore, the overall graphical corrections to $\alpha$, $\mu$, and $D$ are
given by

$\displaystyle\eta_{\alpha}=\eta^{dir}_{\alpha}-\eta^{dir}_{D}\,,~{}~{}~{}\eta_{\mu}=\eta^{dir}_{\mu}-\eta^{dir}_{D}\,,$
(15) $\displaystyle\eta_{D}=\eta^{dir}_{D}-2\eta^{dir}_{D}=-\eta^{dir}_{D}\,.$
(16)

Using these considerations, and explicitly evaluating the graphical
corrections to one-loop order, we find the following recursion relations for
the first uncontrolled approximation, i.e., using the BZ
$-\infty<q_{y}<\infty$, $-\Lambda<q_{x}<\Lambda$ in exactly $d=2$:

$\displaystyle\frac{{\rm d}\ln\alpha}{{\rm d}\ell}$ $\displaystyle=$
$\displaystyle z-2\zeta+2-\frac{3g^{U1}}{4},$ (17) $\displaystyle\frac{{\rm
d}\ln\mu_{x}}{{\rm d}\ell}$ $\displaystyle=$ $\displaystyle
z-2+\frac{g^{U1}}{4},$ (18) $\displaystyle\frac{{\rm d}\ln D}{{\rm d}\ell}$
$\displaystyle=$ $\displaystyle z-3+\zeta-\frac{g^{U1}}{4},$ (19)

where

$g^{U1}=\frac{\alpha^{1/2}D}{4\mu^{3/2}\pi{\Lambda}}\,.$ (20)

Using the recursion relations (11), (12), (13), and the definition of
$g^{U1}$, implies the following recursion relation for $g^{U1}$:

$\frac{{\rm d}\ln g^{U1}}{{\rm
d}\ell}=\frac{2+2\eta_{D}+\eta_{\alpha}-3\eta_{\mu}}{2}\approx\left({1-g^{U1}}\right),$
(21)

where the first equality is exact, and the second, approximate, equality is
valid only to one-loop order, obtained from combining (20) with (17,18,19) .
This implies that there is a stable fixed point at $(g^{U1})^{*}=1$. While
this fixed point is only valid to one-loop order, since $g^{U1}$ _does_ have a
non-zero fixed point value at _all_ orders, and since the first equality in
(21) is valid to _all_ orders, setting the L.H.S. of (21) to 0 leads to an
_exact_ identity between $\eta_{\alpha,\mu,D}$:

$2+2\eta_{D}+\eta_{\alpha}-3\eta_{\mu}=0\,.$ (22)

An exactly analogous calculation for the second uncontrolled calculation,
using the BZ $-\infty<q_{x}<\infty$, $-\Lambda<q_{y}<\Lambda$, obtains the
following recursion relations to one-loop order in exactly $d=2$:

$\displaystyle\frac{{\rm d}\ln\alpha}{{\rm
d}\ell}=\left(z-2\zeta+2-\frac{7g^{U2}}{8}\right),$ (23)
$\displaystyle\frac{{\rm d}\ln\mu_{x}}{{\rm
d}\ell}=\left(z-2+\frac{5g^{U2}}{8}\right),$ (24) $\displaystyle\frac{{\rm
d}\ln D}{{\rm d}\ell}=\left(z-3+\zeta-\frac{3g^{U2}}{8}\right),$ (25)

where

$g^{U2}=\frac{\zeta\alpha^{1/4}D}{8\sqrt{2}\mu^{5/4}\pi\sqrt{\Lambda}}\,.$
(26)

This definition of $g^{U2}$ implies the recursion relation

$\frac{{\rm d}\ln g^{U2}}{{\rm
d}\ell}=\frac{2\zeta+4\eta_{D}+\eta_{\alpha}-5\eta_{\mu}}{4}\approx\frac{4\zeta-11g^{U2}}{8},$
(27)

where, again, the first equality is formally valid to all orders, but the
second, approximate equality, to only one-loop order. Since $\zeta>0$, this
implies a stable fixed point at $(g^{U2})^{*}={4\zeta^{*}/11}$ to one loop.

Finally, analytically continuing our calculation to $d>2$ by treating the $y$
direction as $d-1$ dimensional, we obtain

$\displaystyle\frac{{\rm d}\ln\alpha}{{\rm
d}\ell}=\left(z-2\zeta+2-\frac{7g^{\epsilon}}{8}\right)\,,$ (28)
$\displaystyle\frac{{\rm d}\ln\mu_{x}}{{\rm
d}\ell}=\left(z-2+\frac{5g^{\epsilon}}{8}\right)\,,$ (29)
$\displaystyle\frac{{\rm d}\ln D}{{\rm
d}\ell}=\left(z-3+(3-d)\zeta-\frac{3g^{\epsilon}}{8}\right)$ (30)

where

$g^{\epsilon}=\frac{\zeta\alpha^{1/4}DS_{d-1}\Lambda^{d-5/2}}{8\sqrt{2}(2\pi)^{d-1}\mu^{5/4}}$
(31)

with $\Lambda$ being the large wavenumber cut-off in $q_{y}$ and the closed
recursion relation for $g^{\epsilon}$ is

$\frac{{\rm d}\ln g^{\epsilon}}{{\rm
d}\ell}=\frac{(10-4d)\zeta+4\eta_{D}+\eta_{\alpha}-5\eta_{\mu}}{4}\approx\frac{5-2d}{2}\zeta-\frac{11}{8}g^{\epsilon}.$
(32)

For $d<{5/2}$, these flow equations imply a fixed point at
$(g^{\epsilon})^{*}=\frac{8\epsilon}{11}\zeta^{*}$ where $\epsilon=5/2-d$.

We now use the trajectory integral-matching formalism Nelson_traj  to
calculate
$C\left(\mathbf{r},t\right)=\langle\mathbf{u}(\mathbf{r},t)\cdot\mathbf{u}(\mathbf{0},0)\rangle$,
by relating the correlators in the original system to those of the rescaled
system, via

$\displaystyle C\left(\alpha_{0},\mu_{x0},D_{{0}},\mathbf{r},t\right)$ (33)
$\displaystyle=$ $\displaystyle
e^{2\chi\ell}C\left[\alpha(\ell),\mu_{x}(\ell),D(\ell),{|y|\over
e^{\zeta\ell}},{|x-v_{b}t|\over e^{\ell}},{|t|\over
e^{z\ell}}\right]\,,~{}~{}~{}~{}$

with $\alpha$, $\mu_{x}$ and $D$ controlling the magnitude of
$C\left(\mathbf{r},t\right)$, and the subscript “0” denotes the bare values of
the parameters. Note the argument $x-v_{b}t$ appears in these expressions
because we have undone the aforementioned Galilean boost to return to the
laboratory coordinate system.

We will illustrate this calculation in detail for the first uncontrolled
approximation; for the other approximations the calculation is very similar.
We choose $\zeta$ and $z$ to fix $\alpha$, $\mu_{x}$, and $D$ (of course, due
to the exact relation (10), fixing $\zeta$ also fixes $\chi$). This gives

$z^{U1}={7\over 4}\,,~{}~{}~{}\zeta^{U1}={3\over
2}\,,~{}~{}~{}\chi^{U1}=-{1\over 2}\,.$ (34)

Then setting $\ell=\ln\left(\Lambda|x{-v_{b}t}|\right)$ casts the R.H.S. of
(33) in the form (2) where

$\displaystyle{\cal
G}_{{}_{UET}}\left({|t|\over|x{-v_{b}t}|^{z}},{|y|\over|x{-v_{b}t}|^{\zeta}}\right)$
$\displaystyle\equiv$
$\displaystyle\Lambda^{2\chi}C\left[\alpha_{0},{\mu_{x}}_{0},D_{0},{|y|\over(|x{-v_{b}t}|\Lambda)^{\zeta}},{1\over\Lambda},{|t|\over(|x{-v_{b}t}|\Lambda)^{z}}\right]\,.$

Since the static limit of the correlator $C\left(\mathbf{r},t=0\right)$ has
been obtained in Chen  _et al._ (2016), we focus here on the dynamic limit
$C\left({\mathbf{0}},t\right)$. In this limit we expect the correlator to be a
power law of $t$ only. As a result, the scaling function ${\cal G}_{{}_{UET}}$
must behave as $\left(t/|x-v_{b}t|^{z}\right)^{\chi\over z}$ to cancel out the
prefactor $|x-v_{b}t|^{2\chi}$ in (2). Thus we obtain (4). We remind the
reader that the static exponents $\zeta^{U1}$ and $\chi^{U1}$ in this
approximation, coincidentally, coincide with the exact static exponents
obtained in Chen  _et al._ (2016). We however note that as these are one-loop
calculations, the static exponents have no reason to coincide with their exact
values and, indeed, they do not in the other two approximations.

Likewise, the second uncontrolled approximation yields

$z^{U2}={23\over 14}\,,~{}~{}~{}\zeta^{U2}={11\over
7}\,,~{}~{}~{}\chi^{U2}=-{4\over 7}\,,$ (36)

and the epsilon expansion

$\displaystyle z^{\epsilon}$ $\displaystyle=$ $\displaystyle 2-{10\over
11}\epsilon={17\over 11}\,,$ (37) $\displaystyle\zeta^{\epsilon}$
$\displaystyle=$ $\displaystyle 2-{12\over 11}\epsilon={16\over 11}\,,$ (38)
$\displaystyle\chi^{\epsilon}$ $\displaystyle=$ $\displaystyle-1+{12\over
11}\epsilon=-{5\over 11}\,,$ (39)

where $\epsilon=5/2-d$ and in the second equality $d=2$ is taken.

Taking an average over the three different sets of values of the exponents
leads to the value $z=1.65\pm 0.06$ quoted in the introduction, $\zeta=1.51\pm
0.03$, and $\chi=-0.51\pm 0.03$. Note that the predictions for $\zeta$ and
$\chi$ are extremely close to their exact values $\zeta=1.5$ and $\chi=-0.50$,
which implies our prediction for $z$ is almost certainly also very accurate.
This implies that though the linear theory suggested that velocity
fluctuations in incompressible flocks are diffusive, nonlinearities _lower_
the dynamical exponent, which makes the dynamics _super_ diffusive: that is,
distance scales like $t^{1/z}$ with $z<2$, which means distance grows faster
with time than the $t^{1/2}$ diffusive law.

In summary, we have calculated the dynamical exponent $z$ characterizing the
scaling of velocity fluctuations with time in incompressible two-dimensional
flocks. We did so using _three_ different approximation schemes, and
demonstrated that the dynamics is significantly modified; specifically, the
value of $z$ is substantially changed by nonlinearities, which turn the
fluctuations superdiffusive. Interestingly, due to the mapping Chen  _et al._
(2016) between this system and the model-A dynamics of a divergence-free
magnet Kashuba , this also implies that magnetisation fluctuations of a
divergence-free magnet or a two-dimensional magnet with two-dimensional
dipolar interactions are also superdiffusive.

###### Acknowledgements.

L.C. acknowledges support by the National Science Foundation of China (under
Grant No. 11874420), and thanks the MPI-PKS, where the early stage of this
work was performed, for their support. J. T. thanks the Max Planck Institute
for the Physics of Complex Systems, Dresden, Germany, for their support
through a Martin Gutzwiller Fellowship during this period. AM was supported by
a TALENT fellowship from CY Cergy Paris Université.

## References

  * (1) N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg Models,” Physical Review Letters 17, 1133 (1966); P.C. Hohenberg, “Existence of long-range order in one and two dimensions,” Physical Review 158, 383 (1967).
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, “Novel Type of Phase Transition in a System of Self-Driven Particles,” Physical Review Letters 75, 1226–1229 (1995).
  * Toner and Tu (1995) J. Toner and Y. Tu, “Long-Range Order in a Two-Dimensional Dynamical $XY$ Model: How Birds Fly Together,” Physical Review Letters 75, 4326–4329 (1995).
  * Toner and Tu (1998) J. Toner and Y. Tu, “Flocks, herds, and schools: A quantitative theory of flocking,” Physical Review E 58, 4828–4858 (1998).
  * (5) J. Toner, “Reanalysis of the hydrodynamic theory of fluid, polar-ordered flocks,” Physical Review E 86, 031918 (1012).
  * (6) H. Chaté, “Dry Aligning Dilute Active Matter,” Annual Review of Condensed Matter Physics 11, 189 (2020).
  * (7) L. Chen, J. Toner, and C. F. Lee, “Critical phenomenon of the order-disorder transition in incompressible flocks,” New Journal of Physics 17, 042002 (2015).
  * Chen  _et al._ (2016) L. Chen, C. F. Lee, and J. Toner, “Mapping two-dimensional polar active fluids to two-dimensional soap and one-dimensional sandblasting,” Nature Communications 7, 12215 (2016.)
  * (9) L. Chen, C. F. Lee, and J. Toner, “Incompressible polar active fluids in the moving phase in dimensions $d>2$,” New Journal of Physics 20, 113035 (2018).
  * (10) A. Maitra, P. Srivastava, M. C. Marchetti, S. Ramaswamy, M. Lenz, “Swimmer suspensions on substrates: anomalous stability and long-range order,” Physical Review Letters 124, 028002 (2020).
  * (11) L. Chen, C. F. Lee, A. Maitra, and J. Toner, “Packed Swarms on Dirt: Two-Dimensional Incompressible Flocks with Quenched and Annealed Disorder,” Physical Review Letters 129, 188004 (2022).
  * (12) L. Chen, C. F. Lee, A. Maitra, and J. Toner, “Hydrodynamic theory of two-dimensional incompressible polar active fluids with quenched and annealed disorder, ” Physical Review E 106, 044608 (2022).
  * (13) J. Toner, “Birth, Death, and Flight: A Theory of Malthusian Flocks,” Physical Review Letters 108, 088102 (2012).
  * Chen  _et al._ (2020a) L. Chen, C. F. Lee, and J. Toner, “Moving, Reproducing, and Dying Beyond Flatland: Malthusian Flocks in Dimensions $d>2$,” Physical Review Letters 125, 098003 (2020a), arXiv:2001.01300 .
  * Chen  _et al._ (2020b) L. Chen, C. F. Lee, and J. Toner, “Universality class for a nonequilibrium state of matter: A $d=4-\epsilon$ expansion study of Malthusian flocks,” Physical Review E 102, 022610 (2020b).
  * (16) A. Kashuba, “Exact scaling of spin-wave correlations in the 2D XY ferromagnet with dipolar forces,” Physical Review Letters 73, 2264 (1994).
  * (17) L. Golubović and Z.-G. Wang, “Anharmonic elasticity of smectics A and the Kardar-Parisi-Zhang model,” Physical Review Letters 69, 2535 (1992); “Kardar-Parisi-Zhang model and anomalous elasticity of two- and three-dimensional smectic-$A$ liquid crystals,” Physical Review E 49, 2567 (1994).
  * (18) M. Kardar, G. Parisi and Y.-C. Zhang, “Dynamic scaling of growing interfaces,” Physical Review Letters 56, 889-892 (1986);
  * (19) Alternatatively, we could act on eqn. (7) with $P_{xi}(\mathbf{q})$ to get the EOM for $u_{x}(\tilde{\mathbf{q}})$. However, the resultant equation is equivalent to (8), since $u_{x}(\tilde{\mathbf{q}})$ and $u_{y}(\tilde{\mathbf{q}})$ are not independent, but rather related by the incompressibility constraint $q_{x}u_{x}(\tilde{\mathbf{q}})+q_{y}u_{y}(\tilde{\mathbf{q}})=0$.
  * Forster  _et al._ (1977) D. Forster, D. R. Nelson, and M. J. Stephen, “Large-distance and long-time properties of a randomly stirred fluid,” Physical Review A 16, 732–749 (1977).
  * (21) P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, (Cambridge University Press, Cambridge, UK, 1995).
  * (22) P. C. Hohenberg and B. I. Halperin, “Theory of dynamic ctitical phenomena,” Review of Modern Physics 49, 435 (1977).
  * (23) D. R. Nelson, “Crossover scaling functions and renormalization-group trajectory integrals,” Physical Review B11, 3504 (1975).

Dynamics of packed swarms: time-displaced correlators of two dimensional
incompressible flocks
Supplementary material

##  Appendix A Review of the DRG procedure

Our approach is that developed by forster_pra77 , to which the reader is
referred for a more detailed discussion.

We start with the EOM in Fourier space:

$-{\rm i}\omega
u_{y}(\tilde{\mathbf{q}})=P_{yx}\left(\mathbf{q}\right)\mathcal{F}_{\tilde{\mathbf{q}}}\left[-\alpha\left(u_{x}+{u_{y}^{2}\over
2}\right)\right]+P_{yy}\left(\mathbf{q}\right)\mathcal{F}_{\tilde{\mathbf{q}}}\left[-\alpha\left(u_{x}+{u_{y}^{2}\over
2}\right)u_{y}+\mu_{x}\partial_{x}^{2}u_{y}+f_{y}\right]\,,$ (1)

where

$\tilde{\mathbf{q}}\equiv(\omega,\mathbf{q})\,,~{}~{}~{}P_{yy}(\mathbf{q})=q_{x}^{2}/q^{2}\,,~{}~{}~{}P_{xy}(\mathbf{q})=-q_{x}q_{y}/q^{2}\,.$
(2)

We can solve (1) for $u_{y}(\tilde{\mathbf{q}})$ to get

$\displaystyle u_{y}(\tilde{\mathbf{q}})$ $\displaystyle=$ $\displaystyle
G(\tilde{\mathbf{q}})\Bigg{[}P_{yy}(\mathbf{q})f_{y}(\tilde{\mathbf{q}})+\left(\frac{\alpha}{2}\right)P_{yx}(\mathbf{q})\int_{\tilde{\mathbf{k}}}u_{y}(\tilde{\mathbf{k}})u_{y}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})-\alpha
P_{yy}(\mathbf{q})\int_{\tilde{\mathbf{k}}}u_{y}(\tilde{\mathbf{k}})u_{x}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})$
(3)
$\displaystyle-\left(\frac{\alpha}{2}\right)P_{yy}(\mathbf{q})\int_{\tilde{\mathbf{k}},\tilde{\mathbf{k}}^{\prime}}u_{y}(\tilde{\mathbf{k}})u_{y}(\tilde{\mathbf{k}}^{\prime})u_{y}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}}-\tilde{\mathbf{k}}^{\prime})\Bigg{]}\,,$

where

$\displaystyle\tilde{\mathbf{k}}\equiv(\Omega,\mathbf{k})\,,~{}~{}~{}\tilde{\mathbf{k}}^{\prime}\equiv(\Omega^{\prime},\mathbf{k}^{\prime})\,,$
(4) $\displaystyle\int_{\tilde{\mathbf{k}}}\equiv\int_{-\infty}^{\infty}{{\rm
d}\Omega\over\sqrt{2\pi}}\int{{\rm
d}^{d}k\over(2\pi)^{d/2}}\,,~{}~{}~{}\int_{\tilde{\mathbf{k}}^{\prime}}\equiv\int_{-\infty}^{\infty}{{\rm
d}\Omega^{\prime}\over\sqrt{2\pi}}\int{{\rm
d}^{d}k^{\prime}\over(2\pi)^{d/2}}\,,$ (5)

where the integrals of $\mathbf{k}$ and $\mathbf{k}^{\prime}$ are over the
entire Brillouin zone. The propagator $G(\tilde{\mathbf{q}})$ is defined as:

$G(\tilde{\mathbf{q}})=\frac{1}{-{\rm
i}\omega+\mu_{x}q_{x}^{2}+\mu_{\perp}q_{y}^{2}+\alpha\frac{q_{y}^{2}}{q^{2}}}\approx\frac{1}{-{\rm
i}\omega+\mu_{x}q_{x}^{2}+\alpha\frac{q_{y}^{2}}{q_{x}^{2}}}\ ,$ (7)

where the second, approximate, equality is asymptotically exact in the long
wave length limit $q\to 0$.

The statistics of the random force $f_{y}$ in Fourier space are given by

$\langle
f_{y}(\tilde{\mathbf{q}})f_{y}(\tilde{\mathbf{q}}^{\prime})\rangle=2D\delta^{d}(\mathbf{q}+\mathbf{q}^{\prime})\delta(\omega+\omega^{\prime})\
.$ (8)

In the linear approximation, Eq. (3) reduces to

$u_{y}(\tilde{\mathbf{q}})=G(\tilde{\mathbf{q}})P_{yy}(\mathbf{q})f_{y}(\tilde{\mathbf{q}})\,.$
(9)

Autocorrelating both sides of (11) we get

$\langle
u_{y}(\tilde{\mathbf{q}})u_{y}(\tilde{\mathbf{q}}^{\prime})\rangle=2D\left(q_{x}^{2}\over
q^{2}\right){C(\tilde{\mathbf{q}})\delta(\omega+\omega^{\prime})\delta(\mathbf{q}+\mathbf{q}^{\prime})}\,,$
(10)

where

$C(\tilde{\mathbf{q}})=\frac{2D}{\omega^{2}+\left(\mu_{x}q_{x}^{2}+\alpha{q_{y}^{2}\over
q_{x}^{2}}\right)^{2}}\ .$ (11)

Combining this with the incompressibility constraint
$u_{x}(\tilde{\mathbf{q}})=-q_{y}u_{y}(\tilde{\mathbf{q}})/q_{x}$ leads to the
correlators for other components of $\mathbf{u}(\tilde{\mathbf{q}})$:

$\displaystyle\langle
u_{x}(\tilde{\mathbf{q}})u_{x}(\tilde{\mathbf{q}}^{\prime})\rangle=2D\left(q_{y}^{2}\over
q^{2}\right){C(\tilde{\mathbf{q}})\delta(\omega+\omega^{\prime})\delta(\mathbf{q}+\mathbf{q}^{\prime})},$
(12a) $\displaystyle\langle
u_{x}(\tilde{\mathbf{q}})u_{y}(\tilde{\mathbf{q}}^{\prime})\rangle=-2D\left(q_{x}q_{y}\over
q^{2}\right)C(\tilde{\mathbf{q}}){\delta(\omega+\omega^{\prime})\delta(\mathbf{q}+\mathbf{q}^{\prime})}\,.$
(12b)

Adding (10) and (12a) together, we get the full auto-correlator of
$\mathbf{u}$:

$\langle\mathbf{u}(\tilde{\mathbf{q}})\cdot\mathbf{u}(\tilde{\mathbf{q}}^{\prime})\rangle=2D{C(\tilde{\mathbf{q}})\delta(\omega+\omega^{\prime})\delta(\mathbf{q}+\mathbf{q}^{\prime})}\,.$
(13)

We see that in the low frequency and small wavenumber limit (i.e., $\omega\to
0$, $q\to 0$)
$\langle\mathbf{u}(\tilde{\mathbf{q}})\cdot\mathbf{u}(\tilde{\mathbf{q}}^{\prime})\rangle$
diverges most strongly in the regime $\omega\sim q_{y}\sim q_{x}^{2}$. In this
dominant regime, $P_{yy}(\mathbf{q})$ can be approximated as 1. Therefore, the
EOM (3) can be simplified as

$\displaystyle u_{y}(\tilde{\mathbf{q}})$ $\displaystyle=$ $\displaystyle
G(\tilde{\mathbf{q}})\Bigg{[}f_{y}(\tilde{\mathbf{q}})+\left(\frac{\alpha}{2}\right)P_{yx}(\mathbf{q})\int_{\tilde{\mathbf{k}}}u_{y}(\tilde{\mathbf{k}})u_{y}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})-\alpha\int_{\tilde{\mathbf{k}}}u_{y}(\tilde{\mathbf{k}})u_{x}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})$
(14)
$\displaystyle-\left(\frac{\alpha}{2}\right)\int_{\tilde{\mathbf{k}},\tilde{\mathbf{k}}^{\prime}}u_{y}(\tilde{\mathbf{k}})u_{y}(\tilde{\mathbf{k}}^{\prime})u_{y}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}}-\tilde{\mathbf{k}}^{\prime})\Bigg{]}\,.$

Now we perform the perturbative DRG. First we partition the field
$\mathbf{u}(\tilde{\mathbf{q}})$ into “slow” and “fast” parts
$\mathbf{u}^{<}(\tilde{\mathbf{q}})$ and $\mathbf{u}^{>}(\tilde{\mathbf{q}})$,
which are supported on small $\mathbf{q}$’s and large $\mathbf{q}$’s,
respectively. This can be represented by Feynman diagrams. The solutions for
$u_{y}(\mathbf{q})$ (14) and for $u_{y}^{<}(\mathbf{q})$ after the
decomposition are represented by the graphs in Figs. 1 and 3, respectively,
and the meanings of the various elements of the graphs are depicted in Fig. 2.
Next we iteratively solve for $\mathbf{u}^{>}(\tilde{\mathbf{q}})$. Then we
average out all the short wave length (i.e., “$>$”) fluctuations. For
instance, the third graph on the RHS in Fig. 4 represents a contribution to
the large wave length noise
$G(\tilde{\mathbf{q}})f_{y}^{<}(\tilde{\mathbf{q}})$. The correction to the
noise strength can be calculated by autocorrelating this contribution, which
is represented by the graphs (a) and (b) in Fig. 5. Also, for the fourth graph
on the RHS in Fig. 4 we do a pairwise average over the short wave length
noise. This is represented by connecting the two legs labeled with “$>$”. The
resultant graph is just the one (d) in Fig. 6, representing a contribution to
the terms linear in $\mathbf{u}^{<}(\tilde{\mathbf{q}})$. Following this
procedure we get all the other one-loop graphs listed in Figs. 5 and 6.

Eventually, eliminating $\mathbf{u}^{>}(\tilde{\mathbf{q}})$ leads to an
effective EOM for $\mathbf{u}^{<}(\tilde{\mathbf{q}})$,

$\displaystyle u^{<}_{y}(\tilde{\mathbf{q}})$ $\displaystyle=$ $\displaystyle
G(\tilde{\mathbf{q}})\Bigg{[}f^{<}_{y}(\tilde{\mathbf{q}})+\delta
f^{<}_{y}(\tilde{\mathbf{q}})+{\rm i}\omega\left(\eta_{\omega}{\rm
d}\ell\right)u^{<}_{y}(\tilde{\mathbf{q}})-\alpha\left(\eta_{\alpha}^{dir}{\rm
d}\ell\right)P_{yx}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})-\mu_{x}\left(\eta_{\mu}^{dir}{\rm
d}\ell\right)q_{x}^{2}u^{<}_{y}(\tilde{\mathbf{q}})\Bigg{]}$ (15)
$\displaystyle+G(\tilde{\mathbf{q}})\Bigg{\\{}\alpha\left(1+\eta_{\alpha}^{dir}{\rm
d}\ell\right)\mbox{NL}\left[\mathbf{u}^{<}(\tilde{\mathbf{q}})\right]\Bigg{\\}}\,,$

where ${\rm d}\ell$ is an infinitesimal number quantifying the thickness of
the momentum shell occupied by the large wavevectors, and
$\mbox{NL}\left[\mathbf{u}^{<}(\tilde{\mathbf{q}})\right]$ stands for the
nonlinear part:

$\displaystyle\mbox{NL}\left[\mathbf{u}(\tilde{\mathbf{q}})\right]=\left(\frac{1}{2}\right)P_{yx}(\mathbf{q})\int_{\tilde{\mathbf{k}}}u_{y}(\tilde{\mathbf{k}})u_{y}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})-\int_{\tilde{\mathbf{k}}}u_{y}(\tilde{\mathbf{k}})u_{x}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})-\left(\frac{1}{2}\right)\int_{\tilde{\mathbf{k}},\tilde{\mathbf{k}}^{\prime}}u_{y}(\tilde{\mathbf{k}})u_{y}(\tilde{\mathbf{k}}^{\prime})u_{y}(\tilde{\mathbf{q}}-\tilde{\mathbf{k}}-\tilde{\mathbf{k}}^{\prime})\,.$
(16)

The autocorrelation of the noise contribution $\delta
f^{<}_{y}(\tilde{\mathbf{q}})$ is given by

$\langle 2\delta f^{<}_{y}(\tilde{\mathbf{q}})\delta
f^{<}_{y}(\tilde{\mathbf{q}}^{\prime})\rangle=2\left(1+\eta_{D}^{dir}{\rm
d}\ell\right)D\delta^{d}(\mathbf{q}+\mathbf{q}^{\prime})\delta(\omega+\omega^{\prime})\,.$
(17)

The detailed calculation of $\eta_{\alpha,D,\mu}^{dir}$ to one-loop order is
given in Sec. B.

We multiply both sides of (15) by $\left[G(\tilde{\mathbf{q}})\right]^{-1}$
and rearrange the terms to get

$\displaystyle-{\rm i}\omega(1+\eta_{\omega}{\rm
d}\ell)u^{<}_{y}(\tilde{\mathbf{q}})$ $\displaystyle=$
$\displaystyle\Bigg{[}f^{\prime}_{y}(\tilde{\mathbf{q}})-\alpha\left(1+\eta_{\alpha}^{dir}{\rm
d}\ell\right)P_{yx}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})-\mu_{x}\left(1+\eta_{\mu}^{dir}{\rm
d}\ell\right)q_{x}^{2}u^{<}_{y}(\tilde{\mathbf{q}})\Bigg{]}$ (18)
$\displaystyle+\Bigg{\\{}\alpha\left(1+\eta_{\alpha}^{dir}{\rm
d}\ell\right)\mbox{NL}\left[\mathbf{u}^{<}(\tilde{\mathbf{q}})\right]\Bigg{\\}}\,.$

Then we divide both sides of (18) by $(1+\eta_{\omega}{\rm d}\ell)$. This
gives

$\displaystyle-{\rm i}\omega u^{<}_{y}(\tilde{\mathbf{q}})$ $\displaystyle=$
$\displaystyle\Bigg{[}{f^{\prime}_{y}(\tilde{\mathbf{q}})\over
1+\eta_{\omega}{\rm d}\ell}-{\alpha\left(1+\eta_{\alpha}^{dir}{\rm
d}\ell\right)\over 1+\eta_{\omega}{\rm
d}\ell}P_{yx}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})-{\mu_{x}\left(1+\eta_{\mu}^{dir}{\rm
d}\ell\right)\over 1+\eta_{\omega}{\rm
d}\ell}q_{x}^{2}u^{<}_{y}(\tilde{\mathbf{q}})\Bigg{]}$ (19)
$\displaystyle+\Bigg{\\{}{\alpha\left(1+\eta_{\alpha}^{dir}{\rm
d}\ell\right)\over 1+\eta_{\omega}{\rm
d}\ell}\mbox{NL}\left[\mathbf{u}^{<}(\tilde{\mathbf{q}})\right]\Bigg{\\}}\,,$

which is essentially (14) with “renormalized” noise strength and coefficients:

$\displaystyle D_{r}={\left(1+\eta_{D}^{dir}{\rm
d}\ell\right)D\over(1+\eta_{\omega}{\rm
d}\ell)^{2}}=\left[1+\left(\eta_{D}^{dir}-2\eta_{\omega}\right){\rm
d}\ell\right]D\,,$ (20)
$\displaystyle\alpha_{r}={\left(1+\eta_{\alpha}^{dir}{\rm
d}\ell\right)\alpha\over 1+\eta_{\omega}{\rm
d}\ell}=\left[1+\left(\eta_{\alpha}^{dir}-\eta_{\omega}\right){\rm
d}\ell\right]\alpha\,,$ (21)
$\displaystyle\mu_{x}^{r}={\left(1+\eta_{\mu}^{dir}{\rm
d}\ell\right)\mu_{x}^{r}\over 1+\eta_{\omega}{\rm
d}\ell}=\left[1+\left(\eta_{\mu}^{dir}-\eta_{\omega}\right){\rm
d}\ell\right]\mu_{x}\,.$ (22)

where

$\eta_{D}=\eta^{dir}_{D}-2\eta_{\omega}\,;~{}~{}~{}\eta_{\alpha}=\eta^{dir}_{\alpha}-\eta_{\omega};\,~{}~{}~{}\eta_{\mu}=\eta^{dir}_{\mu}-\eta_{\omega}\,..$
(23)

Note that the graphical correction for each parameter has two parts: the
$\eta$ with superscript “dir” and $\eta_{\omega}$ (or 2$\eta_{\omega}$ in the
case of $D$). The former comes from the “direct” renormalization of the
parameter in the process of eliminating $\mathbf{u}^{>}(\tilde{\mathbf{q}})$,
while the latter comes from rescaling the renormalized coefficient of $-{\rm
i}\omega u_{y}^{<}(\tilde{\mathbf{q}})$ back to 1. Furthermore, the graphical
corrections to the three $\alpha$’s are identical. This is not a coincidence,
but rather due to the symmetry (i.e., rotation invariance) of (1), which must
be preserved under the DRG transformation. This symmetry argument implies that
we do not have to explicitly calculate graphical corrections to the quadratic
and cubic vertices themselves, thus simplifying our calculation considerably.

Also note that EOM (1) can be written as a detailed-balance-obeying dynamics
of the form

$\partial_{t}u_{l}=-\Gamma\frac{\delta\mathcal{H}}{\delta u_{l}}+f_{l}$ (24)

where $\mathcal{H}$ is the Hamiltonian for a divergence-free magnet, and the
friction coefficient $\Gamma$ is taken to be unity. This dynamics obeys
fluctuation-dissipation theorem which _must_ be preserved under
renormalization. This imposes a constraint: the ratio of the remormalized
noise strength to the renormalized friction coefficient is fixed The
renormalized friction coefficient appears as the common factor
$1/(1+\eta_{\omega}{\rm d}\ell)$ on the RHS in (19). The renormalized noise
strength is given by (20). Then it follows from the constraint imposed by
fluctuation-dissipation theorem that

$\eta_{\omega}=\eta_{D}^{dir}\,.$ (25)

Inserting (25) into (23) we get

$\displaystyle\eta_{D}=-\eta^{dir}_{D}\,,$ (26)
$\displaystyle\eta_{\alpha}=\eta^{dir}_{\alpha}-\eta_{D}^{dir}\,,$ (27)
$\displaystyle\eta_{\mu}=\eta^{dir}_{\mu}-\eta_{D}^{dir}\,.$ (28)

That is, because of (25) , we do not need to explicitly calculate the
graphical corrections to $-i\omega u^{<}_{y}(\tilde{\mathbf{q}})$.

Finally we rescale time, lengths, and fields as

$\displaystyle t\to te^{z{\rm d}\ell},~{}~{}x\to xe^{{\rm d}\ell},~{}~{}y\to
ye^{\zeta{\rm d}\ell},$ (29) $\displaystyle u_{y}\to u_{y}e^{\chi{\rm
d}\ell},~{}~{}u_{x}\to u_{x}e^{\left(\chi+1-\zeta\right){\rm d}\ell}\,,$ (30)

to restore the Brillouin zone to its original size.

Figure 1: The graphical representation of the EOM (14). Figure 2: Definitions
of the various elements in the Feynman diagrams: (a)
$={q_{y}\left(1-2\delta_{ix}\right)\over
q_{i}}G\left(\tilde{\mathbf{q}}\right)$, (b) $=u_{i}(\tilde{\mathbf{q}})$, (c)
$={q_{y}\left(1-2\delta_{ix}\right)\over
q_{i}}G\left(\tilde{\mathbf{q}}\right)f_{y}(\tilde{\mathbf{q}})$, (d)
$=P_{ij}(\mathbf{q})C(\tilde{\mathbf{q}})$, (e)
${\rm{the~{}circle}}=-\alpha\left(1-{1\over
2}\delta_{yi}\right)\left[\delta_{xi}+P_{yx}\left(\mathbf{q}\right)\delta_{yi}\right]$,
(f) the square $=-\alpha/2$. Figure 3: The graphical representation of the
solution for $u^{<}(\tilde{\mathbf{q}})$. Figure 4: The graphical
representation of the solution for $u^{<}(\tilde{\mathbf{q}})$ after the
iteration of $u^{>}(\tilde{\mathbf{q}})$. Figure 5: The graphical
representations of the correction to $D$. Figure 6: The graphical
representations of the correction to the propagator.

To take further the advantage of rotational invariance of the EOM, we choose
$\chi=1-\zeta$ so that the three $\alpha$’s remain identical upon
renormalization. This relation is, unsurprisingly, obeyed by their exact
values $\chi=-1/2$ and $\zeta=3/2$.

So far we have completed one cycle of DRG procedure. Repeating this cycle ad
infinitum leads to the recursion relations for $D$, $\alpha$, and $\mu_{x}$
given in the main text.

##  Appendix B Calculation of graphical corrections

###  B.1 Direct noise correction

According to the rules illustrated in Fig. 2, the graphs in Fig. 5 give
contributions to the autocorrelation of
$G(\tilde{\mathbf{q}})f^{<}_{y}(\tilde{\mathbf{q}})$. Specifically, graph (a)
gives

$\alpha^{2}G(\tilde{\mathbf{q}})G(-\tilde{\mathbf{q}})\int{{\rm d}\Omega\over
2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{y}^{2}}{k^{2}}\frac{(q_{x}-k_{x})^{2}}{|\mathbf{q}-\mathbf{k}|^{2}}C(\tilde{\mathbf{k}})C(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})\,,$
(31)

where $\int_{>}{\rm d}^{d}k$ integrates over large $\mathbf{k}$ (i.e., the
$\mathbf{k}$ in the momentum shell); graph (b) gives

$\alpha^{2}G(\tilde{\mathbf{q}})G(-\tilde{\mathbf{q}})\int{{\rm d}\Omega\over
2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{x}k_{y}}{k^{2}}\frac{(q_{x}-k_{x})(q_{y}-k_{y})}{|\mathbf{q}-\mathbf{k}|^{2}}C(\tilde{\mathbf{k}})C(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})\,;$
(32)

graph (c) gives

$-2\alpha^{2}P_{xy}(\mathbf{q})G(\tilde{\mathbf{q}})G(-\tilde{\mathbf{q}})\int{{\rm
d}\Omega\over 2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{x}k_{y}}{k^{2}}\frac{(q_{x}-k_{x})^{2}}{|\mathbf{q}-\mathbf{k}|^{2}}C(\tilde{\mathbf{k}})C(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})\,;$
(33)

graph (d) gives

$2\alpha^{2}\left[P_{xy}(\mathbf{q})\right]^{2}G(\tilde{\mathbf{q}})G(-\tilde{\mathbf{q}})\int{{\rm
d}\Omega\over 2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{x}^{2}}{k^{2}}\frac{(q_{x}-k_{x})^{2}}{|\mathbf{q}-\mathbf{k}|^{2}}C(\tilde{\mathbf{k}})C(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})\,;$
(34)

Since we are interested in the contributions to the white noise, we can set
the wavevector $\mathbf{q}$ inside the loop integral to zero. Since the
prefactors $P_{xy}(\mathbf{q})$ and $\left[P_{xy}(\mathbf{q})\right]^{2}$ are
both much less than 1 in the dominant region $q_{y}\sim q_{x}^{2}$, the
contributions (33) and (34) are negligible compared to (31) and (32).

Upon setting set the wavevector $\mathbf{q}$ inside the loop integral to zero,
(31) and (32) become identical, and their sum is

$2\alpha^{2}G(\tilde{\mathbf{q}})G(-\tilde{\mathbf{q}})\int{{\rm d}\Omega\over
2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{y}^{2}}{k_{x}^{2}}C(\tilde{\mathbf{k}})C(-\tilde{\mathbf{k}})=2\alpha^{2}D^{2}G(\tilde{\mathbf{q}})G(-\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{y}^{2}k_{x}^{4}}{\left(\mu k_{x}^{4}+\alpha
k_{y}^{2}\right)^{3}}\,,$ (35)

which implies the following correction to the noise strength:

$\delta D^{dir}=\alpha^{2}D^{2}I_{1}\,,$ (36)

where

$I_{1}=\int{{\rm d}^{d}k\over(2\pi)^{d}}\frac{k_{y}^{2}k_{x}^{4}}{\left(\mu
k_{x}^{4}+\alpha k_{y}^{2}\right)^{3}}\,.$ (37)

Using the results for $I_{1}$ calculated in Sec. C we get

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta
D^{dir}=\frac{3}{8}Dg^{u}_{\parallel}{\rm d}\ell\,,$ (38)
$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta
D^{dir}=\frac{1}{4}Dg^{u}_{\perp}{\rm d}\ell\,,$ (39)
$\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta
D^{dir}=\frac{3}{8}Dg_{h}{\rm d}\ell\,,$ (40)

where

$g^{u}_{\parallel}=\frac{\zeta\alpha^{1/4}D}{8\sqrt{2}\mu^{5/4}\pi\sqrt{\Lambda}}\,,~{}~{}~{}g^{u}_{\perp}=\frac{\alpha^{1/2}D}{4\mu^{3/2}\pi{\Lambda}}\,,~{}~{}~{}g_{h}=\frac{\zeta\alpha^{1/4}DS_{d-1}\Lambda^{d-5/2}}{8\sqrt{2}(2\pi)^{d-1}\mu^{5/4}}\,.$
(41)

This implies

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{D}^{dir}=\frac{3}{8}g^{u}_{\parallel}\,,$
(42) $\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{D}^{dir}=\frac{1}{4}g^{u}_{\perp}\,,$
(43) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{D}^{dir}=\frac{3}{8}g_{h}\,.$
(44)

Inserting (42), (43), and (44) into (26), respectively, we obtain

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{D}=-\frac{3}{8}g^{u}_{\parallel}\,,$
(45) $\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{D}=-\frac{1}{4}g^{u}_{\perp}\,,$
(46) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{D}=-\frac{3}{8}g_{h}\,.$
(47)

###  B.2 Direct $\alpha$ correction

The only graphical contribution to the term
$G(\tilde{\mathbf{q}})P_{yx}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})$ in (15)
comes from graph (i) in Fig. 6, which corresponds to

$\alpha^{2}G(\tilde{\mathbf{q}})P_{xy}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})\int{{\rm
d}\Omega\over 2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{k_{x}^{2}}{k^{2}}G(\tilde{\mathbf{k}})C(\tilde{\mathbf{q}}-\tilde{\mathbf{k}})\,.$
(48)

Setting both $\omega$ and $\mathbf{q}$ to zero inside the integrals and
integrating over $\Omega$, we have

$\frac{\alpha^{2}D}{2}G(\tilde{\mathbf{q}})P_{xy}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{{k_{x}^{6}}}{(\mu k_{x}^{4}+\alpha
k_{y}^{2})^{2}k^{2}}\,.$ (49)

Noting that the integral over the region $k_{x}^{2}/k^{2}\approx 1$ dominates,
we can rewrite (49) as

$\frac{\alpha^{2}D}{2}G(\tilde{\mathbf{q}})P_{xy}(\mathbf{q})u_{x}^{<}(\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{{k_{x}^{4}}}{(\mu k_{x}^{4}+\alpha
k_{y}^{2})^{2}}\,,$ (50)

which implies

$\delta\alpha^{dir}=-\frac{\alpha^{2}D}{2}I_{2}\,,$ (51)

where $I_{2}$ is defined as

$I_{2}=\int{{\rm d}^{d}k\over(2\pi)^{d}}\frac{{k_{x}^{4}}}{(\mu
k_{x}^{4}+\alpha k_{y}^{2})^{2}}\,.$ (52)

Using the results for $I_{2}$ calculated in Sec. C we get

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta\alpha^{dir}=-{1\over
2}\alpha g^{u}_{\parallel}{\rm d}\ell\,,$ (53)
$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta\alpha^{dir}=-{1\over
2}\alpha g^{u}_{\perp}{\rm d}\ell\,,$ (54) $\displaystyle\mbox{$\epsilon$
expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta\alpha^{dir}=-{1\over
2}\alpha g_{h}{\rm d}\ell\,,$ (55)

which implies

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\alpha}^{dir}=-{1\over
2}g^{u}_{\parallel}\,,$ (56) $\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\alpha}^{dir}=-{1\over
2}g^{u}_{\perp}\,,$ (57) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\alpha}^{dir}=-{1\over
2}g_{h}\,.$ (58)

Inserting (56), (57), and (58) into (27), respectively, we obtain

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\alpha}=-\frac{7}{8}g^{u}_{\parallel}\,,$
(59) $\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\alpha}=-\frac{3}{4}g^{u}_{\perp}\,,$
(60) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\alpha}=-\frac{7}{8}g_{h}\,.$
(61)

###  B.3 Direct $\mu$ correction

The graphical contribution to the term
$G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}^{<}(\tilde{\mathbf{q}})$ in (15) comes
from graphs (a), (b), (c), and (d) in Fig. 6. Note that in the previous two
calculations how $\mathbf{q}$ splits between the two legs inside the loop does
not affect the results since we set $\mathbf{q}={\mathbf{0}}$. It does now
since the loop integrals have to be expanded to $O(q_{x}^{2})$ to get the
factor $q_{x}^{2}$. We make the choice that $\mathbf{q}$ splits equally on the
two internal legs. The same choice was made in the seminal paper for “randomly
stirred incompressible fluids”forster_pra77 . For notational convenience we
define the wavevectors on the two internal leg as the follows:

${\bf p}={{\bf q}\over 2}+{\bf k}\,,~{}~{}~{}{\bf h}={{\bf q}\over 2}-{\bf
k}\,.$ (62)

Then graph (a) in Fig. 6 can be expressed as

$\alpha^{2}G(\tilde{\mathbf{q}})u_{y}(\tilde{\mathbf{q}})\int{{\rm
d}\Omega\over 2\pi}\int{{\rm
d}^{d}k\over(2\pi)^{d}}\left(\frac{p_{y}}{p_{x}}\right)^{2}G(\tilde{\mathbf{p}})C(\tilde{\mathbf{h}})\,.$
(63)

We set $\omega$ and $q_{y}$ to zero inside the integral and integrate over
$\Omega$. Next we expand the resultant formula to quadratic order in $q_{x}$
and focus on the $q_{x}^{2}$ piece. This gives the following contribution to
$G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}^{<}(\tilde{\mathbf{q}})$ in (15):

$\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\left(5\alpha^{2}k_{y}^{6}-10\alpha\mu
k_{x}^{4}k_{y}^{4}+\mu^{2}k_{x}^{8}k_{y}^{2}\right)}{\left(\alpha
k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})(5I_{3}-10I_{4}+I_{5})\,,$
(64)

where

$\displaystyle I_{3}=\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\alpha^{2}k_{y}^{6}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}\,,$ (65) $\displaystyle I_{4}=\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\alpha\mu k_{x}^{4}k_{y}^{4}}{\left(\alpha
k_{y}^{2}+\mu k_{x}^{4}\right)^{4}},$ (66) $\displaystyle I_{5}=\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\mu^{2}k_{x}^{8}k_{y}^{2}}{\left(\alpha
k_{y}^{2}+\mu k_{x}^{4}\right)^{4}}\,.$ (67)

Similarly, graph (b) in Fig. 6 represents

$\alpha^{2}G(\tilde{\mathbf{q}})u_{y}(\tilde{\mathbf{q}})\int_{\tilde{\bf
k}}\left(\frac{p_{y}}{p_{x}}\right)\left(\frac{h_{y}}{h_{x}}\right)G(\tilde{\mathbf{p}})C(\tilde{\mathbf{h}})\,,$
(68)

from which we get the following contribution to
$G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}^{<}(\tilde{\mathbf{q}})$ in (15):

$\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\left(-\alpha^{2}k_{y}^{6}-14\alpha\mu
k_{x}^{4}k_{y}^{4}+3\mu^{2}k_{x}^{8}k_{y}^{2}\right)}{\left(\alpha
k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})(-I_{3}-14I_{4}+3I_{5})\,.$
(69)

Graph (c) in Fig. 6 represents

$\alpha^{2}G(\tilde{\mathbf{q}})u_{y}(\tilde{\mathbf{q}})\int_{\tilde{\bf
k}}\left(\frac{h_{y}}{h_{x}}\right)^{2}G(\tilde{\mathbf{p}})C(\tilde{\mathbf{h}})$
(70)

from which we get the following contribution to
$G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}^{<}(\tilde{\mathbf{q}})$ in (15):

$\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\left(-3\alpha^{2}k_{y}^{6}-10\alpha\mu
k_{x}^{4}k_{y}^{4}+9\mu^{2}k_{x}^{8}k_{y}^{2}\right)}{\left(\alpha
k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})(-3I_{3}-10I_{4}+9I_{5})\,.$
(71)

Graph (d) in Fig. 6 represents

$\alpha^{2}G(\tilde{\mathbf{q}})u_{y}(\tilde{\mathbf{q}})\int_{\tilde{\bf
k}}\left(\frac{p_{y}}{p_{x}}\right)\left(\frac{h_{y}}{h_{x}}\right)G(\tilde{\mathbf{p}})C(\tilde{\mathbf{h}})$
(72)

which is exactly the same as (68) for graph (b). Therefore, this graph also
yields

$\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\left(-\alpha^{2}k_{y}^{6}-14\alpha\mu
k_{x}^{4}k_{y}^{4}+3\mu^{2}k_{x}^{8}k_{y}^{2}\right)}{\left(\alpha
k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})(-I_{3}-14I_{4}+3I_{5})\,.$
(73)

Summing up (64), (69), (71), and (73), we get the total contribution to
$G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}^{<}(\tilde{\mathbf{q}})$ in (15):

$\frac{\alpha^{2}D}{8}G(\tilde{\mathbf{q}})q_{x}^{2}u_{y}(\tilde{\mathbf{q}})(-48I_{4}+16I_{5})\,,$
(74)

which implies

$\delta\mu^{dir}=-2{\alpha^{2}D}(-3I_{4}+I_{5})\,.$ (75)

Using the results for $I_{4,5}$ calculated in Sec. C we get

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta\mu^{dir}=\mu
g^{u}_{\parallel}{\rm d}\ell\,,$ (76) $\displaystyle\mbox{Uncontrolled
approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta\mu^{dir}={1\over 2}\mu
g^{u}_{\perp}{\rm d}\ell\,,$ (77) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\delta\mu^{dir}=\mu g_{h}{\rm
d}\ell\,,$ (78)

which implies

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\mu}^{dir}=g^{u}_{\parallel}\,,$
(79) $\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\mu}^{dir}={1\over
2}g^{u}_{\perp}\,,$ (80) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\mu}^{dir}=g_{h}\,.$ (81)

Inserting (79), (80), and (81) into (28), respectively, we obtain

$\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\mu}=\frac{5}{8}g^{u}_{\parallel}\,,$
(82) $\displaystyle\mbox{Uncontrolled approximations
with}~{}k_{y}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\mu}=\frac{1}{4}g^{u}_{\perp}\,,$
(83) $\displaystyle\mbox{$\epsilon$ expansion
with}~{}k_{x}\in\\{-\infty,\infty\\}:~{}~{}~{}\eta_{\mu}=\frac{5}{8}g_{h}\,.$
(84)

##  Appendix C Calculation of integrals

In this section we will calculate the integrals $I_{1,2,3,4,5}$ for the three
different schemes.

###  C.1 $I_{1}$

This integral is

$I_{1}=\int{{\rm d}^{d}k\over(2\pi)^{d}}\frac{k_{y}^{2}k_{x}^{4}}{\left(\mu
k_{x}^{4}+\alpha k_{y}^{2}\right)^{3}}.$ (85)

####  C.1.1 Uncontrolled $d=2$ with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{1}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm d}k_{y}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{k_{y}^{2}k_{x}^{4}}{\left(\mu k_{x}^{4}+\alpha
k_{y}^{2}\right)^{3}}=\frac{3\zeta}{64\sqrt{2}\pi}\Lambda^{-1/2}\alpha^{-7/4}\mu^{-5/4}{\rm
d}\ell=\frac{3}{8}\frac{1}{\alpha^{2}D}g^{u}_{\parallel}{\rm d}\ell$ (86)

####  C.1.2 Uncontrolled $d=2$ with BZ $k_{y}\in\\{-\infty,\infty\\}$

$I_{1}=\int_{\Lambda e^{-{\rm d}\ell}}^{\Lambda}{{\rm d}k_{x}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{y}\over
2\pi}\frac{k_{y}^{2}k_{x}^{4}}{\left(\mu k_{x}^{4}+\alpha
k_{y}^{2}\right)^{3}}=\frac{1}{16\pi}\Lambda^{-1}\alpha^{-3/2}\mu^{-3/2}{\rm
d}\ell=\frac{1}{4}\frac{1}{\alpha^{2}D}g^{u}_{\perp}{\rm d}\ell$ (87)

####  C.1.3 Epsilon expansion with BZ $k_{x}\in\\{-\infty,\infty\\}$

We extend the integral to $d$ dimensions by taking $y$ to be $d-1$ dimensional
and $x$ 1 dimensional. Thus we have

$I_{1}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm
d}^{d-1}k_{y}\over(2\pi)^{d-1}}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{k_{y}^{2}k_{x}^{4}}{\left(\mu k_{x}^{4}+\alpha
k_{y}^{2}\right)^{3}}=\frac{3S_{d-1}\zeta}{64\sqrt{2}(2\pi)^{d-1}}\Lambda^{d-5/2}\alpha^{-7/4}\mu^{-5/4}d\ell=\frac{3}{8}\frac{1}{\alpha^{2}D}g_{h}{\rm
d}\ell$ (88)

###  C.2 $I_{2}$

This integral is

$I_{2}=\int{{\rm d}^{d}k\over(2\pi)^{d}}\frac{{k_{x}^{4}}}{(\mu
k_{x}^{4}+\alpha k_{y}^{2})^{2}}.$ (89)

####  C.2.1 Uncontrolled $d=2$ with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{2}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm d}k_{y}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over 2\pi}\frac{{k_{x}^{4}}}{(\mu
k_{x}^{4}+\alpha
k_{y}^{2})^{2}}=\frac{\zeta}{8\sqrt{2}\pi}\Lambda^{-1/2}\alpha^{-3/4}\mu^{-5/4}{\rm
d}\ell=\frac{1}{\alpha D}g^{u}_{\parallel}{\rm d}\ell$ (90)

####  C.2.2 Uncontrolled $d=2$ with BZ $k_{y}\in\\{-\infty,\infty\\}$

$I_{2}=\int_{\Lambda e^{-{\rm d}\ell}}^{\Lambda}{{\rm d}k_{x}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{y}\over 2\pi}\frac{{k_{x}^{4}}}{(\mu
k_{x}^{4}+\alpha
k_{y}^{2})^{2}}=\frac{1}{4\pi}\Lambda^{-1}\alpha^{-1/2}\mu^{-3/2}{\rm
d}\ell=\frac{1}{\alpha D}g^{u}_{\perp}{\rm d}\ell$ (91)

####  C.2.3 Epsilon expansion with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{2}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm
d}^{d-1}k_{y}\over(2\pi)^{d-1}}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{{k_{x}^{4}}}{(\mu k_{x}^{4}+\alpha
k_{y}^{2})^{2}}=\frac{S_{d-1}\zeta}{8\sqrt{2}(2\pi)^{d-1}}\Lambda^{d-5/2}\alpha^{-3/4}\mu^{-5/4}{\rm
d}\ell=\frac{1}{\alpha D}g_{h}{\rm d}\ell$ (92)

###  C.3 $I_{3}$

This integral is

$I_{3}=\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\alpha^{2}k_{y}^{6}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}.$ (93)

####  C.3.1 Uncontrolled $d=2$ with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{3}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm d}k_{y}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{\alpha^{2}k_{y}^{6}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{77\zeta}{256\sqrt{2}\pi}\Lambda^{-1/2}\alpha^{-7/4}\mu^{-1/4}{\rm
d}\ell=\frac{77}{32}\frac{\mu}{\alpha^{2}D}g^{u}_{\parallel}{\rm d}\ell$ (94)

####  C.3.2 Uncontrolled $d=2$ with BZ $k_{y}\in\\{-\infty,\infty\\}$

$I_{3}=\int_{\Lambda e^{-{\rm d}\ell}}^{\Lambda}{{\rm d}k_{x}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{y}\over
2\pi}\frac{\alpha^{2}k_{y}^{6}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{5}{32\pi}\Lambda^{-1}\alpha^{-3/2}\mu^{-1/2}{\rm
d}\ell=\frac{5}{8}\frac{\mu}{\alpha^{2}D}g^{u}_{\perp}{\rm d}\ell$ (95)

####  C.3.3 Epsilon expansion with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{3}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm
d}^{d-1}k_{y}\over(2\pi)^{d-1}}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{\alpha^{2}k_{y}^{6}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{77S_{d-1}\zeta}{256\sqrt{2}(2\pi)^{d-1}}\Lambda^{d-5/2}\alpha^{-7/4}\mu^{-1/4}{\rm
d}\ell=\frac{77}{32}\frac{\mu}{\alpha^{2}D}g_{h}{\rm d}\ell$ (96)

###  C.4 $I_{4}$

This integral is

$I_{4}=\int{{\rm d}^{d}k\over(2\pi)^{d}}\frac{\alpha\mu
k_{x}^{4}k_{y}^{4}}{\left(\alpha k_{y}^{2}+\mu k_{x}^{4}\right)^{4}}.$ (97)

####  C.4.1 Uncontrolled $d=2$ with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{4}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm d}k_{y}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over 2\pi}\frac{\alpha\mu
k_{x}^{4}k_{y}^{4}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{7\zeta}{256\sqrt{2}\pi}\Lambda^{-1/2}\alpha^{-7/4}\mu^{-1/4}{\rm
d}\ell=\frac{7}{32}\frac{\mu}{\alpha^{2}D}g^{u}_{\parallel}{\rm d}\ell$ (98)

####  C.4.2 Uncontrolled $d=2$ with BZ $k_{y}\in\\{-\infty,\infty\\}$

$I_{4}=\int_{\Lambda e^{-{\rm d}\ell}}^{\Lambda}{{\rm d}k_{x}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{y}\over 2\pi}\frac{\alpha\mu
k_{x}^{4}k_{y}^{4}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{1}{32\pi}\Lambda^{-1}\alpha^{-3/2}\mu^{-1/2}{\rm
d}\ell=\frac{1}{8}\frac{\mu}{\alpha^{2}D}g^{u}_{\perp}{\rm d}\ell$ (99)

####  C.4.3 Epsilon expansion with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{4}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm
d}^{d-1}k_{y}\over(2\pi)^{d-1}}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{\alpha\mu k_{x}^{4}k_{y}^{4}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{7S_{d-1}\zeta}{256\sqrt{2}(2\pi)^{d-1}}\Lambda^{d-5/2}\alpha^{-7/4}\mu^{-1/4}{\rm
d}\ell=\frac{7}{32}\frac{\mu}{\alpha^{2}D}g_{h}{\rm d}\ell$ (100)

###  C.5 $I_{5}$

This integral is

$I_{5}=\int{{\rm
d}^{d}k\over(2\pi)^{d}}\frac{\mu^{2}k_{x}^{8}k_{y}^{2}}{\left(\alpha
k_{y}^{2}+\mu k_{x}^{4}\right)^{4}}.$ (101)

####  C.5.1 Uncontrolled $d=2$ with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{5}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm d}k_{y}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{\mu^{2}k_{x}^{8}k_{y}^{2}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{5\zeta}{256\sqrt{2}\pi}\Lambda^{-1/2}\alpha^{-7/4}\mu^{-1/4}{\rm
d}\ell=\frac{5}{32}\frac{\mu}{\alpha^{2}D}g^{u}_{\parallel}{\rm d}\ell$ (102)

####  C.5.2 Uncontrolled $d=2$ with BZ $k_{y}\in\\{-\infty,\infty\\}$

$I_{5}=\int_{\Lambda e^{-{\rm d}\ell}}^{\Lambda}{{\rm d}k_{x}\over
2\pi}\int_{-\infty}^{\infty}{{\rm d}k_{y}\over
2\pi}\frac{\mu^{2}k_{x}^{8}k_{y}^{2}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{1}{32\pi}\Lambda^{-1}\alpha^{-3/2}\mu^{-1/2}{\rm
d}\ell=\frac{1}{8}\frac{\mu}{\alpha^{2}D}g^{u}_{\perp}{\rm d}\ell$ (103)

####  C.5.3 Epsilon expansion with BZ $k_{x}\in\\{-\infty,\infty\\}$

$I_{5}=\int_{\Lambda e^{-\zeta{\rm d}\ell}}^{\Lambda}{{\rm
d}^{d-1}k_{y}\over(2\pi)^{d-1}}\int_{-\infty}^{\infty}{{\rm d}k_{x}\over
2\pi}\frac{\mu^{2}k_{x}^{8}k_{y}^{d}}{\left(\alpha k_{y}^{2}+\mu
k_{x}^{4}\right)^{4}}=\frac{5S_{d-1}\zeta}{256\sqrt{2}(2\pi)^{d-1}}\Lambda^{d-5/2}\alpha^{-7/4}\mu^{-1/4}{\rm
d}\ell=\frac{5}{32}\frac{\mu}{\alpha^{2}D}g_{h}{\rm d}\ell$ (104)

## References

  * (1) D. Forster, D. R. Nelson, and M. J. Stephen, “Large- distance and long-time properties of a randomly stirred fluid”, Physical Review A 16, 732-749 (1977).


# Pseudo Electric Field and Pumping Valley Current in Graphene Nano-bubbles

Naif Hadadi  Department of Physics, King Fahd University of Petroleum and
Minerals, 31261 Dhahran, Saudi Arabia  Adel Belayadi  Department of Physics,
University of Science And Technology Houari Boumediene, Bab-Ezzouar, Algeria
Ahmed AlRabiah  Department of Physics, King Fahd University of Petroleum and
Minerals, 31261 Dhahran, Saudi Arabia  Department of physics, Northwestern
university, USA  Ousmane Ly  Université de Nouakchott, Faculté des Sciences et
Techniques, Département de Physique, Avenue du Roi Faiçal, 2373, Nouakchott,
Mauritanie  Collins Ashu Akosa  Department of Applied Physics, School of
Advanced Science and Engineering, Waseda University, Japan  Michael Vogl
Department of Physics, King Fahd University of Petroleum and Minerals, 31261
Dhahran, Saudi Arabia  Quantum Computing Group, IRC of Intelligent Secure
systems, King Fahd University of Petroleum and Minerals,31261 Dhahran, Saudi
Arabia  Hocine Bahlouli  Department of Physics, King Fahd University of
Petroleum and Minerals, 31261 Dhahran, Saudi Arabia  Aurelien Manchon  Aix-
Marseille Université, CNRS, CINaM, Marseille, France  Adel Abbout
<EMAIL_ADDRESS>Department of Physics, King Fahd University of
Petroleum and Minerals, 31261 Dhahran, Saudi Arabia  Quantum Computing Group,
IRC of Intelligent Secure systems, King Fahd University of Petroleum and
Minerals,31261 Dhahran, Saudi Arabia

###### Abstract

Abstract: The extremely high pseudo-magnetic field emerging in strained
graphene suggests that an oscillating nano-deformation will induce a very high
current even without electric bias. In this paper, we demonstrate the sub-
terahertz (THz) dynamics of a valley-current and the corresponding charge
pumping with a periodically excited nano-bubble. We discuss the amplitude of
the pseudo-electric field and investigate the dependence of the pumped valley
current on the different parameters of the system. Finally, we report the
signature of extra-harmonics generation in the valley current that might lead
to potential modern devices development operating in the nonlinear regime.

##  1 INTRODUCTION

Fig. 1: A mechanical deformation of a graphene sheet. The bump which can have
different shapes (spherical, Gaussian, …) induces an out-of-plane stretching
of the hoppings between neighboring atoms.

An interesting and powerful property of graphene is its ability to be
stretched elastically up to 25$\%$, and to control the induced strains in
different ways [1, 2]. In addition, the unique coupling between mechanical
deformation and electronic structure along with the possibility of deformation
makes graphene attracts considerable attention [3, 4]. This interplay was
examined by observing the Landau levels that can form in graphene due to an
induced strain [4]. Indeed, deformed graphene (Fig. 1) can generate an
effective gauge field $\vec{A}$ [4], to which one can associate a pseudo
magnetic field (PMF), $\vec{B}=\vec{\nabla}\times\vec{A}$. This stimulated
pseudo magnetic field would allow electrons to behave as if they were
subjected to a strong real magnetic field with strength a few orders of
magnitude higher than the one generated by superconducting magnets. [5].

Currently, there is a great interest in generating and controlling the valley
degree of freedom of electrons in semiconductors. Indeed, the ability to
manipulate valley electrons can potentially enable advanced valley-resolved
electronic devices. In particular, this valley controllability opens up the
possibility of using the momentum state of electrons, holes, or excitons as a
completely new paradigm in information processing [6]. Moreover, pure valley
currents are interestingly non-dissipative currents with no accompanying net
charge flow, akin to pure spin currents [7]. This property is very useful in
seeking ultra-low power devices. Different routes have been proposed for
generating valley currents such as using quantum pumps with a Dirac gap [8],
and using electrically induced Berry curvature in bilayer graphene [7].

In the quantum pumping technique, a cyclic change of the AC voltage of the
gates leads to the variation in the scattering matrix of the device and the
generation of a DC current [8]. This method of generating a charge DC current
without bias voltage between two electrodes [9, 8] can be generalized to
account for the spin or even the valley degree of freedom. In this work, a
time-dependent nano-bubble in graphene, that can be created by a time-
dependent voltage of an Atomic Force Microscopy (AFM) tip, capacitively
coupled to the device [10], will be utilized as a quantum pumping device to
examine the possibilities of generating a nonzero valley current with a zero
net charge current. The deformation itself can be initially created by
different techniques such as corrugated substrate engineering [11, 12] or gas
inflation [13].
The pseudo magnetic field associated with the deformation in graphene is
valley dependent and changes its sign between the two valley points K and
K’.This feature allows the PMF to preserve time-reversal symmetry unlike a
real external magnetic field. Consequently, a time manipulation of the PMF
will induce different behaviour for the K and K’ electrons and thus a valley-
resolved study of the current is necessary [14]. It is worth mentioning that
different studies have looked at the effect of other dynamical strains on the
electronic transport. Ref. [15] for example, considered in-plane strain in a
gapped graphene to induce a topological current, Ref. [16] looked at the
symmetry requirements for pumping valley current, while Ref. [17] investigated
the pseudo quantum Hall regime due to strain. Other shapes or mechanisms to
get the valley filtering were also investigated [18, 19]

##  2 MODEL

Fig. 2: Tight-binding model for a quasi-1D graphene system with zigzag edges.
a) shows the local deformation as a local bump in the center of the graphene
sheet, altering the translational invariance of the quasi 1D waveguide. b)
shows a face side of the Gaussian strain bump which is centered in the middle
of the system. The green links in the zoom of c) represent the interface
through which the pumped current is calculated. The red parts in c) are the
leads (reservoirs) that indicate that the system is infinite in that
direction. An animation sketching the system can be found in [20]

Graphene consists of carbon atoms, that are arranged in an hexagonal honeycomb
structure. This arrangement is a triangular lattice with a basis of two atoms
per unit cell where the lattice vector for this structure are
$\overrightarrow{\mathbf{a}}_{1}=\frac{a}{2}(3,\sqrt{3})$ and
$\overrightarrow{\mathbf{a}}_{2}=\frac{a}{2}(3,-\sqrt{3})$ (Fig. 2) and the
constant a is carbon-carbon distance with value $d_{0}$ $\sim$ 1.42 $\AA$
[21]. The electronic structure of strained graphene is studied by utilizing a
first nearest neighbor tight-binding model approximation,
$\text{H}=-\sum_{\langle ij\rangle}t_{ij}c_{i}^{\dagger}c_{j}$ where the sum
runs over the nearest neighboring atoms. Figure 2 shows the semi-infinite
graphene system, where the plane waves are coming along the zigzag direction
from left and right leads. Experimentally, the deformation can be created in a
graphene sheet by using an atomic force microscope’s tip [10, 22, 23]. In the
presence of a nano-bubble the nearest neighbors hopping parameters become all
different. Thus, the strain can be included in the system by altering the
hopping parameter to the strained bond length. The adjustment of the hopping
in the presence of the strain can be described by

$t_{ij}=t_{0}e^{-\beta(\frac{d_{ij}}{d_{0}}-1)},$ (1)

where $t_{0}=2.7\text{ eV}$, $\beta=3.37$ for graphene structure [24].
$d_{ij}$ represents the length of the strained lattice bonds whereas $d_{0}$
is the lattice constant in the absence of strain.

##  3 RESULTS

###  3.1 Pseudo Magnetic Field

In our model, the graphene out-of-plane deformation is local with a
characteristic width $\sigma=5\text{ nm}$, a maximum height $h_{0}=3.5\text{
nm}$ [4, 25, 3] and an overall Gaussian shape modeled by
$z=h_{0}e^{-\frac{x^{2}+y^{2}}{2\sigma^{2}}}$. In the limit of low-energy
carriers, the deformation is equivalent to an induced effective vector
potential given by [1]
$\vec{A}=-\frac{\hbar\beta}{2ea_{0}}((\epsilon_{xx}-\epsilon_{yy})\hat{{\bf\it
x}}-2\epsilon_{xy}\hat{{\bf\it y}})$, where the stress tensors $\epsilon_{ij}$
can be calculated by
$\epsilon_{ij}=\frac{1}{2}[(\partial_{i}z)(\partial_{j}z)]$. A straightforward
calculation leads to the following form of the vector potential,

$\vec{A}_{0}=\frac{-\hbar\beta}{ea_{0}\sigma^{4}}z^{2}(\left(x^{2}-y^{2}\right)\hat{x}-2xy\hat{y}).$
(2)

The corresponding pseudo magnetic field can be obtained by
$\vec{B}_{0}=\vec{\nabla}\times\vec{A}_{0}$:

$\vec{B}_{0}=\frac{4\hbar\beta}{ea_{0}\sigma^{6}}z^{2}(y^{3}-3x^{2}y)\hat{{\bf\it
z}}.$ (3)

The pseudo magnetic field shows alternating positive and negative regions in
the deformation position as shown in Fig. 3. This magnetic field is valley
dependent and changes its sign between the high symmetry points $K$ and
$K^{\prime}$ and is responsible for the valley filtering in graphene nano-
bubbles [3].

Fig. 3: The color map indicates the threefold symmetric pseudo-magnetic field
caused by the circularly symmetric Gaussian deformation. The vector field
illustrates the stimulated electric field by the time-dependent pseudo-
magnetic field. ’a’ is the lattice constant.

Using eq.1, we can demonstrate that the maximum relative change in the hopping
parameter is proportional to $\beta h_{0}^{2}/\sigma^{2}$. This implies that
for a small deformation width $\sigma$, it becomes necessary to reduce the
maximum height to remain within the constraints of stretchable graphene [1,
2]. In our particular case, a numerical examination of max($\Delta
t_{ij}/t_{ij}$) confirms that we are, indeed, within this specified limit.

###  3.2 Induced Pseudo Electric Field

The dynamics of the charge carriers in graphene under a dynamical external
mechanical strain is a challenging problem. The generated frequencies of the
perturbation projects into the sub-Terahertz domain with reported frequencies
up to 200 GHz [26]. We note in passing that such frequencies are much smaller
than the bandwidth of graphene. Therefore, standard techniques like Floquet
theory [27, 28, 29, 30, 31, 32, 33] are hard to apply, which makes a full
numerical treatment necessary and which will be the subject of our work.
Things become more interesting when the strain is a nano-bubble exhibiting a
few hundred Teslas pseudo magnetic field. Indeed, an oscillating bump with a
height $h=h_{0}+\delta h(t)$ will change the vector potential to a new form
$\vec{A}=(1+\frac{\delta h(t)}{h_{0}})^{2}\vec{A}_{0}$. In the Weyl gauge, the
Maxwell equations show an induced pseudo electric field in the limit
$\frac{\delta h(t)}{h_{0}}\ll 1$, to the lowest order:

$\vec{E}(x,y)=-2\frac{\delta\accentset{\mbox{\large.}}{h}(t)}{h_{0}}\vec{A}_{0},$
(4)

where the dot represents the time derivative. This result accounts for the
small perturbations due to a small amplitude in the out-of-plane vibration of
the nano-bubble which allows us to keep our interpretation within the linear
response theory. It is worth noting that although the pseudo magnetic field is
strong (hundreds of Teslas), the induced electric field is only a few kV/cm.
This is mainly explained by the fact that these large values of $B_{0}$ are
rather due to the large gradient of $\vec{A}_{0}$ (variation on very short
scales) than to its magnitude. The pseudo-electric field map is depicted in
Fig. 3 where a threefold symmetry emerges and the system alternates between an
inward and outward pseudo electric field for the $K$ and $K^{\prime}$ valleys.

Fig. 4: The generated currents through the interface shown in Fig. 2. The two
types of currents at $K$ and $K^{\prime}$ are in phase and average to zero.
The charge current is obtained for $h_{0}=3.5\text{ nm}$ and $\delta
h=0.35\text{ nm}$. $E_{F}=0.28t_{0}$. T is the period of oscillation of the
bubble.

###  3.3 Charge and Valley Currents

As, in the case of charge-pumping by the dynamics of a skyrmion deformation
[34, 9] and spin-pumping by the precession of magnetic textures [35, 36], an
oscillating nano-bubble is expected to alter the charge density and induce a
charge/valley current through the interface between the system and the lead
[37]. To investigate this property, we adopt a nano-bubble oscillating at
hundreds GHz with $\delta h(t)=\delta h_{0}\sin(\omega_{0}t+\phi)$ with a
height modulation strength $h_{0}=3.5\text{ nm}$, a characteristic width
$\sigma=5\text{ nm}$ and a perturbation magnitude $\delta h_{0}=0.35\text{
nm}$. The way to generate these oscillations is unimportant at this stage but
can be initiated by a capacitively coupled atomic force microscope’s tip or a
laser on top of the strained part. The first step towards the calculation of
the charge current is to solve the time-dependent Schrodinger equation to
obtain the wave functions at different times and then calculate the different
currents in the system in one-body formalism [38, 39]. It is important to
mention that a more accurate procedure involves an integration over the
different contributions from all energies in the Fermi sea. This was safely
avoided due to the small driving frequencies (compared to the Fermi level
$E_{F}=270\text{ meV}$) and the small excitations considered in this study
($\hbar\omega/E_{F}\ll 1$). This reduces the needed numerical resources and
keeps the problem at the Fermi surface (similar to the case of the spin-
pumping theory). Figure 4 shows the result of the pumped charge current
$I_{\text{charge}}=I_{K}+I_{K^{\prime}}$, which expectedly averages to zero as
explained by the scattering approach to parametric pumping. Indeed, varying
one parameter is not enough to generate a sustainable charge current [9].

Fig. 5: The charge current (Left) and the valley current (Right) pumped
through an interface with the lead (see figure 2). The average of each of them
is zero. The charge oscillates at the same frequency as the nano-bubble,
whereas the pumped valley current shows higher order harmonics.The blue dots
represent the times at which the current maps were plotted in Fig. 7. T is the
period of oscillation of the nanobubble Fig. 6:  The current map of the pumped
charge current (top) and the valley pumped current at different times for 6
modes with $E_{F}=0.07t_{0}$. The nano-bubble size is $\sigma=5\text{ nm}$ and
the height $\delta h_{0}=0.1h_{0}$ . The position is expressed in nm.The times
at which the maps are plotted are inside one period of oscillation and
increase from left to right. (the 4 different phases of an oscillation). The
blue dots in Fig. 5 show the time where the maps are plotted. ( for these
figures, the time increases from left to right.)

The frequency of the signal is the same as that of the oscillating nano-
bubble. So far, these results look like the same as for the spin pumping of
$I_{x}$, $I_{y}$ and $I_{\text{charge}}$ in Ferromagnets (precessing around
the z-axis): zero average and the same frequency as the driving perturbation.
Since the pseudo electric field is valley dependent, the contribution coming
from each K-point to the charge current is investigated and plotted in Fig. 4.
The two currents, $I_{K}$, and $I_{K^{\prime}}$ have amplitudes of the same
order of magnitude and both oscillate at the same angular frequency
$\omega_{0}$. The two currents are not exactly the same, because for a Fermi
energy slightly higher than zero, the modes in the two different valleys are
non-symmetric and have slightly different wavenumbers. The valley current,
which is defined as $I_{\text{valley}}=I_{K}-I_{K^{\prime}}$ is of huge
interest in low energy consumption devices. The reason behind it is that the
net charge transferred can be zero whereas at the same time the valley current
is not vanishing. This eliminates the dissipation of energy as heat while
transporting information through the valleys’ degrees of freedom. The valley
current is well defined in the interface region (away from the bubble),
because the lead is uniform and the conducting modes (incoming wave functions
from the lead) do not mix. One has to compute the time-dependent wave function
for each mode, using T-kwant package [38, 39], and identify the different
modes in each valley. The contribution of each mode to the current between
sites $a$ and $b$ reads:

$I^{ab}=i\left(\boldsymbol{\psi}_{b}^{\dagger}\left(H_{ab}\right)^{\dagger}\boldsymbol{\psi}_{a}-\boldsymbol{\psi}_{a}^{\dagger}H_{ab}\boldsymbol{\psi}_{b}\right)$
(5)

where, $\boldsymbol{\psi}_{a(b)}$ is the wave function of a given mode at site
$a(b)$. $H_{ab}$ is the matrix element between sites a and b (in our case, it
represents the hopping parameter between a and b in graphene). The
contribution of modes of a given valley are summed together to finally express
$I_{\text{valley}}$. Inside the scattering region, this procedure can still be
used, but the current is not exactly a pure valley current. Indeed, the
scattering between the valleys is possible, yet it remains small (the points K
and K’ are far apart in the Brillouin zone) [14]. For this reason, the
approximation of the local current in the scattering region is still valid
[14]. Another approach is to project the wave functions on the valley-resolved
ones of pristine graphene [24].

Fig. 5 shows the result for the valley current pumped through the same
interface. The average is again zero and the amplitude is much smaller than
that of the charge current. The most important thing to notice is that the
frequency of the $I_{\text{valley}}$ is not exactly the same as that of the
nano-bubble. As we can see in Fig. 7, the valley current exhibits a main
frequency of $3\omega_{0}$ with small contributions from $\omega_{0}$,
$2\omega_{0}$ and $5\omega_{0}$. The appearance of higher-order harmonics was
recently reported for spin-pumping in the presence of spin-orbit interaction
[40] or in the presence of non-collinear magnetic structure [41]. In our case,
none of the discussed causes in [40, 41] are present in our system.

Fig. 7:  Left) The Fourier transform of the Valley current.The signal exhibits
higher order frequencies (multiple of $\omega_{0}$). Right) (blue curve) The
amplitude of the pumped electric current vs the amplitude of the bump
oscillation. (green curve) The amplitude of the valley current vs the
amplitude of the bump oscillation. $h_{0}=3.5\text{ nm}$, $\sigma=5\text{ nm}$
and the Fermi energy is $E_{F}=270\text{ meV}$. The width of the system is
$W=50\text{ nm}$.

The charge and valley pumping depend on the parameters of the system. As we
can see it from the expression of the pseudo electric field, the latter is
proportional to the perturbation $\frac{\delta h_{o}}{h_{0}}$. Naturally, this
is expected to be the same for the charge and valley currents. Fig. 7 b) shows
the amplitude of the $I_{\text{charge}}$ and $I_{\text{valley}}$ which are
clearly linearly varying with the change $\frac{\delta h_{o}}{h_{0}}$. For
stronger perturbations, higher non-linear effects are expected to manifest.
Another way to increase the different pumped currents is to play with the
frequency and use a Terahertz spectrum. This is a challenging task, but very
promising results suggest that this projection is reachable [26]. The zero-
average of the valley current makes it very difficult to exploit in devices
since no accumulation can be obtained and only oscillating current
redistribution of the valleys is achieved. With proper engineering, we still
can think of it as the basis for a potentially new Field-effect transistor
[42].

###  3.4 Local currents

Since the pseudo electric field is local and mainly appears in the central
system, investigating the current at different regions than the interface with
the lead might be of a great interest. The plot of the current map will show
how the regions of the maximum pseudo magnetic field influence the current
direction and generation. The wave functions are calculated at all sites of
the systems which allows us to express the local currents. Fig. 6 illustrates
the map obtained at different times during one period. The current lines for
the charge and the valley are symmetric which is just the reflection of the
symmetry of the Gaussian bump considered in this work. This will be different
if a triaxial strain was adopted. In fact, it is known that a Gaussian
deformation filters the valleys whereas the triaxial one splits them [3]. The
Fermi energy is chosen low to reduce the number of modes due to the limited
computing resources in this kind of 2D maps (current is calculated between
each two links, at each time) and this doesn’t change the general conclusions.
The maps for the valley current and the charge current are intrinsically
different: for the charge current, it is hard to see the positions of the
valley-dependent maxima of the pseudo-magnetic field. The situation is
different for the valley current. We clearly notice loops of different
orientations showing the dependence of the pseudo-magnetic field on the valley
of the carriers. The size of the current’s loops can be changed via the
parameter $\sigma$ ie, creating very peaked bumps or smoothly varying ones
over a larger region. The number of modes for the chosen Fermi energy is very
small. In order to increase the current, one needs to increase it either by
changing the width of the system or by increasing the Fermi energy.
It’s worth mentioning that we did investigate deformations with elliptical
symmetry, deliberately tilting them to disrupt left-right symmetry.
Unfortunately, this approach did not lead to a sustained average valley
current.

## CONCLUSION

The pseudo magnetic field generated by a strained graphene sheet was analyzed.
In the presence of a time-dependent oscillating strain, a corresponding pseudo
electric field is stimulated giving rise to a redistribution of the quantum
states in momentum space, leading to charge and valley pumping. We
demonstrated that the valley current exhibits extra-harmonics and that on
average, the pseudo electric field did not sustain a direct current, and
rather a zero average was obtained. This work opens the door to further
interesting studies like the pseudo-electric field generated by the
propagation of mechanical waves in graphene or the distribution of local
electric currents due to phonons and temperature effects on graphene.

## ACKNOWLEDGMENT

The authors acknowledge computing time on the supercomputer SHAHEEN at KAUST
Supercomputing Centre and the team assistance. A.A. and M.V. gratefully
acknowledge the support provided by the Deanship of Research Oversight and
Coordination (DROC) at King Fahd University of Petroleum and Minerals (KFUPM)
for funding this work through exploratory research grant No. ER221002.

## REFERENCES

  * Masir  _et al._ [2013] M. R. Masir, D. Moldovan, and F. Peeters, Pseudo magnetic field in strained graphene: Revisited, Solid state communications 175, 76 (2013).
  * Guinea  _et al._ [2010a] F. Guinea, A. Geim, M. Katsnelson, and K. Novoselov, Generating quantizing pseudomagnetic fields by bending graphene ribbons, Physical Review B 81, 035408 (2010a).
  * Settnes  _et al._ [2016] M. Settnes, S. R. Power, M. Brandbyge, and A.-P. Jauho, Graphene nanobubbles as valley filters and beam splitters, Physical review letters 117, 276801 (2016).
  * De Juan  _et al._ [2011] F. De Juan, A. Cortijo, M. A. Vozmediano, and A. Cano, Aharonov–bohm interferences from local deformations in graphene, Nature Physics 7, 810 (2011).
  * Guinea  _et al._ [2010b] F. Guinea, M. Katsnelson, and A. Geim, Energy gaps and a zero-field quantum hall effect in graphene by strain engineering, Nature Physics 6, 30 (2010b).
  * Vitale  _et al._ [2018] S. A. Vitale, D. Nezich, J. O. Varghese, P. Kim, N. Gedik, P. Jarillo-Herrero, D. Xiao, and M. Rothschild, Valleytronics: opportunities, challenges, and paths forward, Small 14, 1801483 (2018).
  * Shimazaki  _et al._ [2015] Y. Shimazaki, M. Yamamoto, I. V. Borzenets, K. Watanabe, T. Taniguchi, and S. Tarucha, Generation and detection of pure valley current by electrically induced berry curvature in bilayer graphene, Nature Physics 11, 1032 (2015).
  * Wang  _et al._ [2014] J. Wang, Z. Lin, and K. S. Chan, Pure valley current generation in graphene with a dirac gap by quantum pumping, Applied Physics Express 7, 125102 (2014).
  * Brouwer [1998] P. Brouwer, Scattering approach to parametric pumping, Physical Review B 58, R10135 (1998).
  * Klimov  _et al._ [2012] N. N. Klimov, S. Jung, S. Zhu, T. Li, C. A. Wright, S. D. Solares, D. B. Newell, N. B. Zhitenev, and J. A. Stroscio, Electromechanical properties of graphene drumheads, Science 336, 1557 (2012).
  * et al [2014] A. R. P. et al, Strain superlattices and macroscale suspension of graphene induced by corrugated substrates, Nano Lett 14, 5044 (2014).
  * S. Viola Kusminskiy and Guinea [2011] A. H. C. N. S. Viola Kusminskiy, D. K Campbell and F. Guinea, Pinning of a two-dimensional membrane on top of a patterned substrate: The case of graphene, Phys. rev. B 83 (2011).
  * Bunch  _et al._ [2008] J. S. Bunch, S. S. Verbridge, J. S. Alden, A. M. van der Zande, J. M. Parpia, H. G. Craighead, and P. L. McEuen, Impermeable atomic membranes from graphene sheets, Nano lett. 8, 2458 (2008).
  * Nedell  _et al._ [2022] J. G. Nedell, J. Spector, A. Abbout, M. Vogl, and G. A. Fiete, Deep learning of deformation-dependent conductance in thin films: Nanobubbles in graphene, Phys. Rev. B 105, 075425 (2022).
  * Vaezi  _et al._ [2013] A. Vaezi, N. Abedpour, R. Asgari, A. Cortijo, and M. A. H. Vozmediano, Topological electric current from time-dependent elastic deformations in graphene, Phys. Rev. B 88, 125406 (2013).
  * Jiang  _et al._ [2013] Y. Jiang, T. Low, K. Chang, M. I. Katsnelson, and F. Guinea, Generation of pure bulk valley current in graphene, Phys. Rev. Lett. 110, 046601 (2013).
  * Sela  _et al._ [2020] E. Sela, Y. Bloch, F. von Oppen, and M. B. Shalom, Quantum hall response to time-dependent strain gradients in graphene, Phys. Rev. Lett. 124, 026602 (2020).
  * Amasay and Sela [2021] J. Amasay and E. Sela, Transport through dynamic pseudogauge fields and snake states in a corbino geometry, Phys. Rev. B 104, 125428 (2021).
  * Belayadi  _et al._ [2023] A. Belayadi, N. A. Hadadi, P. Vasilopoulos, and A. Abbout, Valley-dependent tunneling through electrostatically created quantum dots in heterostructures of graphene with hexagonal boron nitride, Phys. Rev. B 108, 085419 (2023).
  * et al. [2023] A. A. et al., Supplemental material on website, Phys. Rev. B xxx, xxxx (2023).
  * Carrillo-Bastos  _et al._ [2014] R. Carrillo-Bastos, D. Faria, A. Latgé, F. Mireles, and N. Sandler, Gaussian deformations in graphene ribbons: Flowers and confinement, Physical Review B 90, 041411 (2014).
  * Lee  _et al._ [2008] C. Lee, X. Wei, J. W. Kysar, and J. Hone, Measurement of the elastic properties and intrinsic strength of monolayer graphene, science 321, 385 (2008).
  * Nemes-Incze  _et al._ [2017] P. Nemes-Incze, G. Kukucska, J. Koltai, J. Kürti, C. Hwang, L. Tapasztó, and L. P. Biró, Preparing local strain patterns in graphene by atomic force microscope based indentation, Scientific reports 7, 3035 (2017).
  * Stegmann and Szpak [2018] T. Stegmann and N. Szpak, Current splitting and valley polarization in elastically deformed graphene, 2D Materials 6, 015024 (2018).
  * Pengfei  _et al._ [2019] J. Pengfei, C. Wenjing, Q. Jiabin, Z. Miao, Z. Xiaohu, X. Zhongying, L. Rongda, T. Chuanshan, H. Lin, D. Zengfeng, and W. Xi, Programmable graphene nanobubbles with three-fold symmetric pseudo-magnetic fields, Nature Communications 10, 3127 (2019).
  * et al [2016] Y. W. et al, 200 ghz maximum oscillation frequency in cvd graphene radio frequency transistors, ACS Appl. Mater. Interfaces 8, 25645 (2016).
  * Eckardt and Anisimovas [2015] A. Eckardt and E. Anisimovas, High-frequency approximation for periodically driven quantum systems from a floquet-space perspective, New Journal of Physics 17, 093039 (2015).
  * Bukov  _et al._ [2015] M. Bukov, L. D'Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to floquet engineering, Advances in Physics 64, 139 (2015).
  * Vogl  _et al._ [2022] M. Vogl, S. Chaudhary, and G. A. Fiete, Light driven magnetic transitions in transition metal dichalcogenide heterobilayers, Journal of Physics: Condensed Matter 35, 095801 (2022).
  * Abanin  _et al._ [2017] D. A. Abanin, W. D. Roeck, W. W. Ho, and F. Huveneers, Effective hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems, Physical Review B 95 (2017).
  * Rodriguez-Vega  _et al._ [2021] M. Rodriguez-Vega, M. Vogl, and G. A. Fiete, Low-frequency and moiré–floquet engineering: A review, Annals of Physics 435, 168434 (2021).
  * Vogl  _et al._ [2019] M. Vogl, P. Laurell, A. D. Barr, and G. A. Fiete, Analog of hamilton-jacobi theory for the time-evolution operator, Physical Review A 100 (2019).
  * Vogl  _et al._ [2020] M. Vogl, M. Rodriguez-Vega, and G. A. Fiete, Effective floquet hamiltonian in the low-frequency regime, Phys. Rev. B 101, 024303 (2020).
  * Abbout  _et al._ [2018] A. Abbout, J. Weston, X. Waintal, and A. Manchon, Cooperative charge pumping and enhanced skyrmion mobility, Phys. Rev. Lett. 121, 257203 (2018).
  * Tserkovnyak  _et al._ [2002] Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, Spin pumping and magnetization dynamics in metallic multilayers, Phys. Rev. B 66, 224403 (2002).
  * Tatara and Mizukami [2017] G. Tatara and S. Mizukami, Consistent microscopic analysis of spin pumping effects, Phys. Rev. B 96, 064423 (2017).
  * [37] This current might vanish on average, .
  * Kloss  _et al._ [2021] T. Kloss, J. Weston, B. Gaury, B. Rossignol, C. Groth, and X. Waintal, Tkwant: a software package for time-dependent quantum transport, New Journal of Physics 23, 023025 (2021).
  * Groth  _et al._ [2014] C. W. Groth, M. Wimmer, A. R. Akhmerov, and X. Waintal, Kwant: a software package for quantum transport, New Journal of Physics 16, 063065 (2014).
  * Ly and Manchon [2022] O. Ly and A. Manchon, Spin-orbit coupling induced ultrahigh-harmonic generation from magnetic dynamics, Phys. Rev. B 105, L180415 (2022).
  * Ly [2023] O. Ly, Noncollinear antiferromagnetic textures driven high-harmonic generation from magnetic dynamics in the absence of spin-orbit coupling, Journal of Physics: Condensed Matter 35, 125802 (2023).
  * [42] W. Y. Choi, H. jun Kim, J. Chang, S. H. Han, A. Abbout, H. B. M. Saidaoui, A. Manchon, K. J. Lee, and H. C. Koo, Ferromagnet-free all-electric spin hall transistors, Nano Lett. 2018, 18, 12, 7998–8002 .