# Boundary-Bound Reactions: Pattern Formation
with and without Hydrodynamics

Aiden Huffman<EMAIL_ADDRESS>Henry Shum<EMAIL_ADDRESS>University of Waterloo,
Waterloo, Canada

###### Abstract

We study chemical pattern formation in a fluid between two flat plates and the
effect of such patterns on the formation of convective cells. This patterning
is made possible by assuming the plates are chemically reactive or release
reagents into the fluid, both of which we model as chemical fluxes. We
consider this as a specific example of boundary-bound reactions. When the
chemical reagents do not produce significant solutal buoyancy forces, chemical
instabilities such as the diffusion-driven instabilities described by Turing
are still possible with nonlinear reactions. In the other extreme, constant
fluxes at the two plates do not produce diffusion-driven instabilities but
could give rise to hydrodynamic instabilities analogous to Rayleigh–Bénard
instabilities. If we are able to influence the chemical fluxes along the
domain and select suitable reaction systems, then this presents a mechanism
for the control of chemical and hydrodynamic instabilities and pattern
formation. We study a generic class of models and find necessary conditions
for a bifurcation to pattern formation. Afterwards, we present two examples
derived from the Schnakenberg–Selkov reaction. Our results demonstrate that
buoyancy-driven instabilities can occur even in the absence of density
variations in the steady state. We also find parameter combinations that
result in the formation of convective cells whether gravity acts upwards or
downwards relative to the reactive plate. The wavenumber of the cells and the
direction of the flow at regions of high/low concentration depend on the
orientation, hence, different patterns can be elicited by simply inverting the
device. More generally, our results suggest methods for controlling pattern
formation and convection by tuning reaction parameters. As a consequence, we
can drive and alter fluid flow in a chamber without mechanical pumps by
influencing the chemical instabilities.

Chemo-hydrodynamic instabilities; Turing patterns;
Advection–reaction–diffusion systems; Solutal convection; Reactive boundaries

††preprint: APS/123-QED

##  I Introduction

Chemical pattern formation has been of interest since Turing’s original paper
on chemical morphogenesis [1]. Turing demonstrated that nonlinear reactions
with sufficiently asymmetric diffusion rates between chemical reagents can
have unstable homogeneous steady states; concentration fields tend to approach
steady, inhomogeneous distributions. His paper is a classical example of how
breaking symmetry can cause a dissipative process such as diffusion to result
in more order. The characteristics of a diffusion-driven instability are that
the reaction–diffusion system possesses a homogeneous steady state that is
linearly stable to small homogeneous perturbations of the concentration fields
but unstable to spatially varying perturbations. For a two-reagent
reaction–diffusion system given by

$\displaystyle\diffp{Q}{t}$ $\displaystyle=f(Q,R)+\nabla^{2}Q,$ (1a)
$\displaystyle\diffp{R}{t}$ $\displaystyle=g(Q,R)+D\nabla^{2}R,$ (1b)

for which there exist constants $Q^{*},R^{*}$ satisfying
$f(Q^{*},R^{*})=0=g(Q^{*},R^{*})$, the Turing conditions for pattern formation
can be summarized as rClCrCl f_Q+g_R &¡ 0, f_Qg_R-g_Qf_R ¿ 0,
Df_Q+g_R ¿ 0, f_Qg_R-f_Qg_R ¡ (DfQ+gR)24D, where subscripts denote partial
derivatives evaluated at the homogeneous steady state $Q=Q^{*},R=R^{*}$.
Together, these conditions require $D\neq 1$ in order for pattern formation to
occur. Conditions for Turing instabilities in systems of more than two
reacting and diffusing components have also been determined [2].

Recent extensions of Turing-style instabilities include active membrane models
or bulk–surface reaction–diffusion systems [3, 4, 5]. Among many interesting
properties, these works generally require smaller critical diffusion ratios
for pattern formation. In particular, Levine et al. [3] presented a model with
a critical diffusion ratio of one. This result differs from the classical
Turing instability, which generally requires large ratios of diffusion
coefficients, suggesting that patterning may be easier to realize than the
classical model predicts.

Shklyaev et al. [6] studied a system of two boundary-bound reactions taking
place on opposing boundaries. The term ‘boundary-bound’ refers to a chemical
reaction that takes place at the boundary of a domain but does not occur in
the bulk. This type of reaction could be realized by immobilizing enzymes on
the walls of a chamber so that the corresponding reactions, which produce and
remove various chemical species, occur only at the walls. Therefore, we can
model the chemical reaction as a boundary flux. The boundary-bound reaction
modelled by Shklyaev et al. created an oscillatory instability which was
enhanced or suppressed by the buoyancy-driven flow field the reagents
produced.

In early experiments that demonstrated chemical Turing patterns, reactions
took place in a gel reactor [7, 8]. One of the reasons for this was to
suppress convection, which could alter the stability of Turing patterns [9,
10] or give rise to other chemo-hydrodynamic instabilities [11, 12, 13, 14,
15].

More recent experimental work has shown that boundary-bound reactions can
produce meaningful flows through buoyancy effects in microscale chambers [16,
17, 18] and have been used to perform various tasks such as self-assembly,
mixing, and transport [19, 20, 21]. Control over the flow field has typically
been achieved by patterning the enzyme placement along the boundary.

In contrast, the present work assumes a homogeneous distribution of enzyme and
investigates various means of influencing the flow field by changing the
reaction parameters. Parallels in our design can be drawn from the classical
Rayleigh–Bénard instability, the Rayleigh–Taylor instability, and the double
diffusive instability [22, 23]. Here, the chemical reagents play the same role
as the temperature and salt gradients, which produce various convective
instabilities. A unique component to our problem is the nonlinear boundary
conditions of the reactive surfaces. In our design, Turing patterns can be
viewed as a sub-case of the buoyancy driven instability. As a result, the two
instabilities function cooperatively and produce new phenomena.

After presenting the framework and general results for a class of two-reagent
systems, we examine two specific cases based on the Schnakenberg–Selkov
reactions. It should be noted that one paradigm for realizing arbitrary
reaction systems is to use DNA strand displacement reactions, which allow the
flexibility of constructing large collections of modular subsystems [24, 25,
26].

##  II General Model

We consider a domain consisting of a fluid between two flat plates. We will
assume that these plates are separated by a vertical distance $H$ that is
small relative to their horizontal length. As a consequence, we treat the
domain as unbounded in the horizontal directions. The model formulation is
identical in two and three spatial dimensions apart from the physical units of
some parameters. A simple two-dimensional schematic is provided in Figure 1.

Figure 1: Schematic for the boundary-bound reaction examples considered in
this work. The reactions occur on the surface of the top plates. The model
changes depending on which boundary we add reagents into the system through.

If we assume that density variations in the fluid are produced by the presence
of two chemical reagents $(Q,R)$ and that these density variations are small,
then the resulting fluid flows can be modelled using the Boussinesq
approximation for the incompressible Navier–Stokes equations. We will assume
that the reagents undergo exponential decay in the bulk with rate constants
$(\eta_{Q},\eta_{R})$, diffuse with rates $(d_{Q},d_{R})$, and are advected by
the fluid flow. Therefore, the bulk behaviour can be modelled via

$\displaystyle\frac{\mathrm{D}Q}{\mathrm{D}t}$
$\displaystyle=-\eta_{Q}Q+d_{Q}\nabla^{2}Q,$ (2a)
$\displaystyle\frac{\mathrm{D}R}{\mathrm{D}t}$
$\displaystyle=-\eta_{R}R+d_{R}\nabla^{2}R,$ (2b)
$\displaystyle\rho_{0}\frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t}$
$\displaystyle=-\nabla
p+\mu\nabla^{2}\mathbf{u}+\rho_{0}g(\beta_{Q}Q+\beta_{R}R)\mathbf{e}_{z},$
(2c) $\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,$ (2d)

where $\frac{\mathrm{D}}{\mathrm{D}t}=\frac{\partial}{\partial
t}+\mathbf{u}\cdot\nabla$ denotes the convective derivative operator,
$\rho_{0}$ is the fluid density in the absence of reagents, $p$ is the fluid
pressure, $\mu$ is the dynamic viscosity of the fluid, $g$ is the signed
acceleration due to gravity ($g<0$ for gravity acting downwards), and
$\beta_{Q}$ and $\beta_{R}$ are the solutal expansion coefficients for $Q$ and
$R$, respectively. At the boundary, the enzyme-driven reactions (or other
processes that add or remove regents) are described by chemical fluxes. In
combination with the no-slip condition for the fluid flow, the full set of
boundary conditions is,

$\IEEEeqnarraybox[][c]{rClCrCl}-d_{Q}\diffp[1]Qz\bigg{|}_{z=0}&=f_{0}(Q,R),\qquad
d_{Q}\diffp[1]Qz\bigg{|}_{z=H}=f_{1}(Q,R),\\\
-d_{R}\diffp[1]Rz\bigg{|}_{z=0}=g_{0}(Q,R),\qquad
d_{R}\diffp[1]Rz\bigg{|}_{z=H}=g_{1}(Q,R),\IEEEyesnumber\\\
\mathbf{u}\bigg{|}_{z=0}=\mathbf{0},\qquad\mathbf{u}\bigg{|}_{z=H}=\mathbf{0}.$
(3)

We non-dimensionalize the system using the scalings

$\IEEEeqnarraybox[][c]{rClCrClCrCl}\mathbf{x}&=H\hat{\mathbf{x}},\qquad
t=\frac{H^{2}}{d_{Q}}\hat{t},\qquad Q=Q_{s}\hat{Q},\\\
R=R_{s}\hat{R},\quad\mathbf{u}=\frac{d_{Q}}{H}\hat{\mathbf{u}},\quad
p=\frac{d_{Q}^{2}\rho_{0}}{H^{2}}\hat{p},\\\ $ (4)

where $Q_{s}$ and $R_{s}$ are arbitrary concentration scales that will be
determined by reactions at the boundaries. Letting
$\beta^{\max}=\max\\{|\beta^{Q}Q_{s}|,|\beta^{R}R_{s}|\\}$, the dimensionless
parameters are:

$\IEEEeqnarraybox[][c]{rClCrClCrClCrCl}D&=\frac{d_{R}}{d_{Q}},\qquad\hat{\eta}_{Q}=\frac{H^{2}}{d_{Q}}{\eta}_{Q},\qquad\hat{\eta}_{R}=\frac{H^{2}}{d_{Q}}{\eta}_{R},\\\
\mathrm{Sc}=\frac{\mu}{\rho_{0}d_{Q}},\qquad\mathrm{Gr}=\frac{\rho_{0}^{3}g\beta^{\max}H^{3}}{\mu^{2}},\qquad\delta^{Q}=\frac{\beta^{Q}Q_{s}}{\beta^{\max}},\qquad\delta^{R}=\frac{\beta^{R}R_{s}}{\beta^{\max}}.$
(5)

Here, $\mathrm{Sc}$ is the Schmidt coefficient, $\mathrm{Gr}$ is the Grashof
number and $D$ is the diffusion ratio between the two chemical reagents. By
defining $\delta^{Q}$ and $\delta^{R}$ relative to the maximum modulus of the
reagent expansion coefficients, we ensure that one of $\delta^{Q}$ or
$\delta^{R}$ is equal to $\pm 1$ and the other is bounded in absolute value by
$1$. In Section IV we will refer to the case where $|\delta^{Q}|=1$ or
$|\delta^{R}|=1$ as $Q$\- or $R$-driven flow, respectively.

Applying our non-dimensionalization and dropping the hats on dimensionless
variables, (2) becomes

$\displaystyle\frac{\mathrm{D}Q}{\mathrm{D}t}$
$\displaystyle=-\eta_{Q}Q+\nabla^{2}Q,$ (6a)
$\displaystyle\frac{\mathrm{D}R}{\mathrm{D}t}$
$\displaystyle=-\eta_{R}R+D\nabla^{2}R,$ (6b)
$\displaystyle\frac{\mathrm{D}\mathbf{u}}{\mathrm{D}t}$ $\displaystyle=-\nabla
p+\mathrm{Sc}\nabla^{2}\mathbf{u}+\mathrm{Sc}^{2}\mathrm{Gr}(\delta_{Q}Q+\delta_{R}R)\mathbf{e}_{z},$
(6c) $\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0.$ (6d)

The boundary conditions for the non-dimensionalized system are:

$\IEEEeqnarraybox[][c]{rClCrCl}-\diffp[1]Qz\bigg{|}_{z=0}&=\gamma
f_{0}(Q,R),\qquad\diffp[1]Qz\bigg{|}_{z=1}=\gamma f_{1}(Q,R),\\\
-D\diffp[1]Rz\bigg{|}_{z=0}=\gamma g_{0}(Q,R),\qquad
D\diffp[1]Rz\bigg{|}_{z=1}=\gamma g_{1}(Q,R),\IEEEyesnumber\\\
\mathbf{u}\bigg{|}_{z=0}=\mathbf{0},\qquad\mathbf{u}\bigg{|}_{z=1}=\mathbf{0}.$
(7)

Here, we view $\gamma$ as a general positive scaling for the chemical fluxes
at the boundaries. We note that in Figure 1, we presented a schematic for the
examples considered in Section III

Suppose that the system has an unpatterned (i.e., the concentration and
pressure fields depend only on the $z$ coordinate and not on horizontal
coordinates) steady state without flow. We perform a linear stability analysis
around this steady state considering a perturbation of the form,

$\displaystyle Q(\mathbf{x},t)=Q_{0}(z)+\epsilon Q_{1}(\mathbf{x},t),$
$\displaystyle R=R_{0}(z)+\epsilon R_{1}(\mathbf{x},t),$ (8)
$\displaystyle\mathbf{u}(\mathbf{x},t)=\epsilon\mathbf{u}_{1}(\mathbf{x},t),$
$\displaystyle p(\mathbf{x},t)=p_{0}(z)+\epsilon p_{1}(\mathbf{x},t).$

The resulting system can be simplified by eliminating the pressure term and
dropping the horizontal dynamics. Details of this technique can be found in
Appendix A. Considering the ansatz,

$\begin{bmatrix}Q_{1}(\mathbf{x},t)\\\ R_{1}(\mathbf{x},t)\\\
u_{1,z}(\mathbf{x},t)\\\
p_{1}(\mathbf{x},t)\end{bmatrix}=\begin{bmatrix}q(z)\\\ r(z)\\\ w(z)\\\
p(z)\end{bmatrix}e^{\lambda t+i(k_{x}x+k_{y}y)},$ (9)

where $u_{1,z}$ is the vertical component of $\mathbf{u}_{1}$, we obtain the
perturbation problem:

$\displaystyle\left(\diff[2]{}{z}-(\lambda+\eta_{Q}+k^{2})\right)q$
$\displaystyle=w\diff{Q_{0}}{z},$ (10a)
$\displaystyle\left(\diff[2]{}{z}-\left(\frac{\lambda+\eta_{R}}{D}+k^{2}\right)\right)r$
$\displaystyle=w\diff{R_{0}}{z},$ (10b)
$\displaystyle\left[\left(\diff[2]{}{z}-k^{2}\right)^{2}-\frac{\lambda}{\mathrm{Sc}}\left(\diff[2]{}{z}-k^{2}\right)\right]w$
$\displaystyle=k^{2}\mathrm{Ra}(\delta^{Q}q+\delta^{R}r),$ (10c)

where $k^{2}=k_{x}^{2}+k_{y}^{2}$ and
$\mathrm{Ra}=\mathrm{Gr}\cdot\mathrm{Sc}$ is the solutal Rayleigh number,
subject to the boundary conditions,

$\IEEEeqnarraybox[][c]{rClCrCl}-\diffp[1]qz\bigg{|}_{z=0}&=\gamma(f_{0,Q}\,q+f_{0,R}\,r),\qquad\diffp[1]qz\bigg{|}_{z=1}=\gamma(f_{1,Q}\,q+f_{1,R}\,r),\\\
-D\diffp[1]rz\bigg{|}_{z=0}=\gamma(g_{0,Q}\,q+g_{0,R}\,r),\qquad
D\diffp[1]rz\bigg{|}_{z=1}=\gamma(g_{1,Q}\,q+g_{1,R}\,r),\IEEEyesnumber\\\
w\bigg{|}_{z=0}=0,\qquad w\bigg{|}_{z=1}=0.$ (11)

In (11) and below, the notation $h_{j,C}$ represents the partial derivative of
the function $h_{i}$ with respect to $C$ evaluated at the unpatterned steady
state concentrations $Q=Q_{0}(j),R=R_{0}(j)$, for all combinations of
functions $h=f,g$, boundaries $j=0,1$, and concentration labels $C=Q,R$. For
two spatial dimensions, we remove the $y$ variable and set $k_{y}=0$.

To begin, we assume that the Rayleigh number is small and that the density
variations are negligible. This assumption reduces the problem to studying
pattern formation without convection.

##  III Pattern Formation

When it exists, the general solution to (10) without flow is of the form

$\displaystyle q(z)$
$\displaystyle=a_{1}\cosh(\Gamma_{1,\lambda}z)+a_{2}\sinh(\Gamma_{1,\lambda}z),$
(12a) $\displaystyle r(z)$
$\displaystyle=b_{1}\cosh(\Gamma_{2,\lambda}z)+b_{2}\sinh(\Gamma_{2,\lambda}z),$
(12b) $\displaystyle w(z)$ $\displaystyle=0,$ (12c)

where

$\displaystyle\Gamma_{1,\lambda}$
$\displaystyle=\sqrt{\lambda+\eta_{Q}+k^{2}},$ (13a)
$\displaystyle\Gamma_{2,\lambda}$
$\displaystyle=\sqrt{\frac{\lambda+\eta_{R}}{D}+k^{2}}.$ (13b)

Below, we use the notation $\Gamma_{i}$ to indicate the values of
$\Gamma_{i,\lambda}$ with $\lambda=0$, for $i=1,2$. The boundary conditions
(11) can be expressed as a homogeneous linear system of four equations in
$a_{i},\ b_{i}$. This system must be degenerate for non-trivial solutions to
exist.

Consider the case where reactions for both $Q$ and $R$ occur at a shared
(upper) boundary so that $f_{1}$ and $g_{1}$ have non-zero partial
derivatives. At the bottom boundary, we will assume that $f_{0}$ and $g_{0}$
are constant (possibly zero), implying that their partial derivatives vanish.
Given $\lambda$ and $k$, the resulting system of equations has a non-trivial
solution if and only if rCl 0 &= γ^2det(J_1)cosh(Γ_1,λ)cosh(Γ_2,λ) -
γ(DΓ_2,λcosh(Γ_1,λ)sinh(Γ_2,λ)f_1,Q
\+ Γ_1,λsinh(Γ_1,λ)cosh(Γ_2,λ)g_1,R) + DΓ_1,λΓ_2,λsinh(Γ_1,λ)sinh(Γ_2,λ),
where $J_{1}=\begin{bmatrix}f_{1,Q}&f_{1,R}\\\ g_{1,Q}&g_{1,R}\end{bmatrix}$.

Expressed in this form, we can show that $D\neq 1$ is necessary to make a
finite wavenumber the most unstable mode. Suppose that there exist a
(possibly) complex $\lambda$ and real $k>0$ that satisfy (III) with
$\mathrm{Re}(\lambda)\geq 0$. Then, any other pair
$(\lambda^{\prime},k^{\prime})$ such that
$\Gamma_{i,\lambda}(k)=\Gamma_{i,\lambda^{\prime}}(k^{\prime})$ for $i=1,2$
would also satisfy (III). From (13), we note that defining
$\lambda^{\prime}=\lambda+k^{2}$ and $k^{\prime}=0$ satisfies this
requirement. Since $\mathrm{Re}(\lambda^{\prime})>\mathrm{Re}(\lambda)$, the
long wavelength instability $k^{\prime}=0$ is more linearly unstable.

If we are concerned with transitions from a steady unpatterned state to a
steady pattern with finite wavenumber, then we restrict ourselves to the case
where $\lambda$ is real. Close to such a transition, we may assume that
$\Gamma_{1,\lambda}$ and $\Gamma_{2,\lambda}$ are real and positive. Treating
the entries of $J_{1}$ as independent parameters and assuming that
$\det(J_{1})\neq 0$, then (III) reduces to a quadratic equation for $\gamma$.
Requiring that the global flux $\gamma$ be real, we obtain the condition that
either $f_{1,R}$ has the same sign as $g_{1,Q}$ or

$\left|f_{1,Q}-\frac{\Gamma_{1,\lambda}\tanh(\Gamma_{1,\lambda})}{D\Gamma_{2,\lambda}\tanh(\Gamma_{2,\lambda})}g_{1,R}\right|\geq
2\sqrt{-\frac{\Gamma_{1,\lambda}\tanh(\Gamma_{1,\lambda})}{D\Gamma_{2,\lambda}\tanh(\Gamma_{2,\lambda})}f_{1,R}\,g_{1,Q}}.$
(14)

We seek marginally stable solutions, i.e., those with $\lambda=0$. For the
special case $\eta_{R}=0$, (III) can be rearranged to obtain an explicit
formula for the diffusion ratio when a given wavenumber $k>0$ is marginally
stable,

$D_{c}(k)=\frac{\gamma\left(\Gamma_{1}g_{1,R}\tanh{\left(\Gamma_{1}\right)}-f_{1,Q}g_{1,R}\gamma+f_{1,R}g_{1,Q}\gamma\right)}{k\tanh{\left(k\right)}\left(\Gamma_{1}\tanh{\left(\Gamma_{1}\right)}-f_{1,Q}\gamma\right)},$
(15)

which we refer to as the critical diffusion ratio for the wavenumber $k$.

Assuming positive diffusion coefficients, we seek a positive value for
$D_{c}(k)$, which requires that

$(\Gamma_{1}\tanh(\Gamma_{1})-f_{1,Q}\gamma)(\Gamma_{1}g_{1,R}\tanh(\Gamma_{1})-\gamma\text{det}(J_{1}))>0.$
(16)

A critical diffusion cannot be found if this inequality is not satisfied for
some $k$.

In Section III.1, we will present two examples of the case above with a shared
boundary for reactions. If we consider an alternative setup where $g_{0}$ and
$f_{1}$ have non-trivial derivatives while $g_{1}$ and $f_{0}$ are constant
(non-constant reactions on opposite boundaries), then we would obtain the
necessary condition for solutions with $\lambda\approx 0$ that either
$f_{1,R}$ has the same sign as $g_{0,Q}$ or

$\bigg{|}f_{1,Q}-\frac{\Gamma_{1,\lambda}\tanh(\Gamma_{1,\lambda})}{D\Gamma_{2,\lambda}\tanh(\Gamma_{2,\lambda})}g_{0,R}\bigg{|}>2\sqrt{-\frac{\Gamma_{1,\lambda}\tanh(\Gamma_{1,\lambda})}{D\Gamma_{2,\lambda}\cosh(\Gamma_{1,\lambda})\sinh(\Gamma_{2,\lambda})}f_{1,R}g_{0,Q}}\,,$
(17)

and under the same conditions we used in the derivation of $D_{c}(k)$ above,
the critical diffusion ratio for this setup is positive only if

$(\Gamma_{1}\tanh(\Gamma_{1})-f_{1,Q}\gamma)(\Gamma_{1}g_{0,R}\tanh(\Gamma_{1})-\gamma\text{det}(J_{2}))>0,$
(18)

where $J_{2}=\begin{bmatrix}f_{1,Q}&f_{1,R}\operatorname{sech}(\Gamma_{1})\\\
g_{1,Q}\operatorname{sech}(k)&g_{1,R}\end{bmatrix}.$

###  III.1 Schnakenberg–Selkov Examples

Figure 2: (a) Unpatterned steady state concentration profiles for the full
reaction on the upper boundary. (b) Unpatterned steady states for the
separated reaction. In both cases, $a=0.2$, $b=2$, and $\eta_{Q}=\gamma=D=1$.

We present two explicit examples in a two-dimensional domain that correspond
to the setup in Figure 1 with nonlinear reactions for both species on the
upper boundary. In both cases, we assume $\eta_{R}=0$ and thus $D_{c}(k)$ is
determined by (15). The reaction occurring on the upper plate is based on the
Schnakenberg–Selkov [27, 28] reaction system so that the rates of flux (with
physical dimensions of amount of substance per unit length per unit time) are
modeled by:

$\displaystyle f(Q,R)$ $\displaystyle=\alpha-k_{-1}Q+k_{2}Q^{2}R,$ (19a)
$\displaystyle g(Q,R)$ $\displaystyle=\beta-k_{2}Q^{2}R.$ (19b)

We assume that $\alpha,\,\beta,\,k_{-1}$ and $k_{2}$ are all positive so there
is no confusion about which terms are influxes or out-fluxes. This is non-
dimensionalized by taking

$\IEEEeqnarraybox[][c]{rClCrClCrCl}Q&=\sqrt{\frac{k_{-1}}{k_{2}}}\hat{Q},\quad
R=\sqrt{\frac{k_{-1}}{k_{2}}}\hat{R},\quad\gamma=\frac{Hk_{-1}}{d_{Q}},\\\
\hat{\eta}_{Q}=\frac{H^{2}}{d_{Q}}\eta_{Q},\quad
a=\sqrt{\frac{k_{2}}{k^{3}_{-1}}}\alpha,\quad
b=\sqrt{\frac{k_{2}}{k_{-1}^{3}}}\beta.$ (20)

Dropping the hats on the non-dimensionalized variables, (19) becomes,

$\displaystyle f(Q,R)$ $\displaystyle=\gamma(a-Q+Q^{2}R),$ (21a)
$\displaystyle g(Q,R)$ $\displaystyle=\gamma(b-Q^{2}R).$ (21b)

Placing the entire reaction on the upper boundary at $z=1$ yields the first
example, which we will call the _full reaction_ :

$\IEEEeqnarraybox[][c]{rClCrCl}f_{1}(Q,R)&=a-Q+Q^{2}R,\qquad
g_{1}(Q,R)=b-Q^{2}R,\\\ f_{0}(Q,R)=0,\qquad g_{0}(Q,R)=0.$ (22)

Note that these functions are scaled by the parameter $\gamma$ for the rates
of reaction, as defined in (7). The steady state is given by:

$\displaystyle Q_{0}(z)$
$\displaystyle=\gamma\frac{(a+b)}{\Gamma_{Q}\sinh(\Gamma_{Q})+\gamma\cosh(\Gamma_{Q})}\cosh(\Gamma_{Q}z),$
(23a) $\displaystyle R_{0}(z)$ $\displaystyle=\frac{b}{Q_{0}(1)^{2}},$ (23b)

where $\Gamma_{Q}=\sqrt{\eta_{Q}}$.

In the second example, we will break the reactions into two components and
place the constant $(a,b)$ terms on the bottom boundary. We refer to this as
the _separated reaction_. Physically, this can be viewed as changing the
boundary that provides material into the domain while fixing the enzymes that
catalyze the other reactions on the original boundary. The boundary reaction
functions are

$\IEEEeqnarraybox[][c]{rClCrCl}f_{1}(Q,R)&=-Q+Q^{2}R,\qquad
g_{1}(Q,R)=-Q^{2}R,\\\ f_{0}(Q,R)=a,\qquad g_{0}(Q,R)=b,$ (24)

The steady state is

$\displaystyle Q_{0}(z)$
$\displaystyle=\gamma\frac{\Gamma_{Q}a\cosh{\left(\Gamma_{Q}\left(z-1\right)\right)}+\Gamma_{Q}b\cosh{\left(\Gamma_{Q}z\right)}-a\gamma\sinh{\left(\Gamma_{Q}\left(z-1\right)\right)}}{\Gamma_{Q}\left(\Gamma_{Q}\sinh{\left(\Gamma_{Q}\right)}+\gamma\cosh{\left(\Gamma_{Q}\right)}\right)},$
(25a) $\displaystyle R_{0}(z)$
$\displaystyle=\frac{b\left(\Gamma_{Q}\sinh{\left(\Gamma_{Q}\right)}+\gamma\cosh{\left(\Gamma_{Q}\right)}\right)^{2}}{\gamma^{2}\left(a+b\cosh{\left(\Gamma_{Q}\right)}\right)^{2}}+\frac{\gamma
b}{D}(1-z).$ (25b)

The steady state concentration profiles for the two examples are presented in
Figure 2 for fixed parameter values. There are some key differences between
the two examples. For the separated reaction, we find that $Q_{0}(z)$ is
minimized within the bulk, while for the full reaction, it is minimized at the
non-reactive boundary. Additionally, $R_{0}(z)$ is constant in the full
reaction and linear in the separated reaction.

The partial derivatives appearing in the linearized boundary conditions are
the same for both examples, namely,

$\IEEEeqnarraybox[][c]{rClCrCl}f_{1,Q}&=2Q_{0}(1)R_{0}(1)-1,\qquad
f_{1,R}=Q_{0}(1)^{2},\\\ g_{1,Q}=-2Q_{0}(1)R_{0}(1),\qquad
g_{1,R}=-Q_{0}^{2}(1).$ (26)

Note that, from (23) and (25), $Q_{0}(1)$ and $R_{0}(1)$ are guaranteed to be
positive for $a,b,\gamma>0$. It follows that $f_{1,R}$ and $g_{1,Q}$ have
opposite signs, so real values for $\gamma$ require condition (14) to be
satisfied.

Substituting (26) into inequality (16), we find that the second term on the
left hand side is negative and hence, the first term must also be negative for
the inequality to hold. For the full reaction, substituting in the boundary
concentrations from (23) leads to the condition for positive values of
$D_{c}(k)$ that $b>a$ and

$\gamma>\gamma_{1}(k):=\frac{\Gamma_{1}\tanh(\Gamma_{1})(a+b)-2\Gamma_{Q}b\tanh(\Gamma_{Q})}{b-a}.$
(27)

For the separated reaction, we use (25) to obtain the conditions for positive
values of $D_{c}(k)$ that $b\cosh(\Gamma_{Q})>a$ and

$\gamma>\gamma_{2}(k):=\frac{\Gamma_{1}\tanh(\Gamma_{1})(a+b\cosh(\Gamma_{Q}))-2\Gamma_{Q}b\sinh(\Gamma_{Q})}{b\cosh(\Gamma_{Q})-a}.$
(28)

Since $\Gamma_{1}=\Gamma_{Q}$ when $k=0$, we have that $\gamma_{1}(0)<0$ and
hence, $\gamma_{1}(k)<0$ for all sufficiently small positive wavenumbers.
Similarly, $\gamma_{2}(k)<0$ for all sufficiently small positive $k$.
Additionally, both $\gamma_{1}$ and $\gamma_{2}$ increase monotonically to
$+\infty$ with $k$ so that for any fixed $\gamma>0$, there is a critical
wavenumber $k_{c,i}>0$ at which $\gamma=\gamma_{i}(k_{c,i})$, for $i=1,\,2$.
Hence, positive values of $D_{c}(k)$ exist for all $k$ between 0 and $k_{c,i}$
but not for any wavenumber greater than $k_{c,i}$. In Figure 3 we compare
$D_{c}(k)$ between the two examples. In both cases, after substituting (26)
into (15), $D_{c}(k)$ can be expressed as,

$D_{c}(k)=\frac{\gamma
Q_{0}(1)^{2}(\Gamma_{1}\tanh(\Gamma_{1})+\gamma)}{k\tanh(k)(2\gamma
Q_{0}(1)R_{0}(1)-\gamma-\Gamma_{1}\tanh(\Gamma_{1}))}.$ (29)

We note that $Q_{0}(1)$ and $R_{0}(1)$ depend on $a,b,\gamma$, and $\eta_{Q}$,
but do not depend on $D$ or $k$. By considering the ratio of $D_{c}(k)$ given
by (29) for full and separated reactions over the range of reaction parameter
values and wavenumbers where $D_{c}(k)$ is positive for both, it is possible
to show that the critical diffusion ratio for the full reaction is strictly
larger than for the separated reaction, as shown in Figure 3(a) for one
parameter set. As a consequence, instabilities can arise with smaller
diffusion ratios under the separated reaction compared with the full reaction.
Moreover, we show that $k_{c,1}<k_{c,2}$ for fixed model parameters. Suppose
that $\gamma=\gamma_{1}(k_{c,1})=\gamma_{2}(k_{c,2})$ and let
$\Gamma_{1}^{(i)}=\sqrt{\eta_{Q}+k_{c,i}^{2}}$ denote the value of
$\Gamma_{1}$ corresponding to the threshold wavenumber $k_{c,i}$ for $i=1,2$.
Rearranging (27) and (28), we have that rCl
Γ_1^(2)tanh(Γ_1^(2))-Γ_1^(1)tanh(Γ_1^(1))&=2bΓ_Q sinh(Γ_Q)
(1a+bcosh(ΓQ)-1(a+b)cosh(ΓQ))
\+ 2aγ(1a+b-1a+bcosh(ΓQ)), which is positive if $a,b,\gamma,\Gamma_{Q}>0$.
Since $\Gamma_{1}\tanh(\Gamma_{1})$ increases monotonically with $k$, we
conclude that $k_{c,2}>k_{c,1}$ and hence, a larger range of wavenumbers can
be made unstable using the separated reaction system compared with the full
reaction.

Figure 3: Comparison of critical parameter values for the full and separated
Schnakenberg–Selkov reaction systems. (a) Critical diffusion ratio as a
function of wavenumber $k$. (b) Minimum critical diffusion ratio,
$D_{\min,c}$, as a function of $\eta_{Q}$. We also mark the critical diffusion
ratio for the classical Turing instability with the Schnakenberg–Selkov
reaction system under the same conditions. (c) Global flux thresholds
$\gamma_{1}(k)$ and $\gamma_{2}(k)$ for the full and separated reactions. In
all cases, $(a,b)=(0.2,2.0)$. When unspecified, $\eta_{Q}$ and $\gamma$ are
set to $1$.

####  III.1.1 Pattern Forming Regions and Oscillations

Figure 4:  (a) Concentrations relative to unpatterned steady state, taken
after the system reaches a steady patterned state for the full
Schnakenberg–Selkov reaction with $D=5.1$ and $(a,b)=(0.2,2.0)$. (b) Pattern
forming regions for the full reaction against various $(a,b)$. Contours
indicate the minimum critical diffusion ratio for pattern formation. No
patterns are expected below the line $a=b$. (c) Pattern forming regions for
the separated reaction, with contours for the minimum critical diffusion
ratio. No patterns are expected below the line $a=b\cosh(\Gamma_{Q})$. (d)
Concentrations relative to the unpatterned steady state, taken after the
system reaches a steady patterned state for the separated reaction with
$D=4.2$ and $(a,b)=(0.2,2.0)$. We fix the reaction parameters
$(\eta_{Q},\gamma)=(1,1)$ in all cases.

An example of the patterned state emerging when a finite wavenumber becomes
unstable in the full reaction system without flow is shown in Figure 4(a).
Since we are interested in using the chemical instabilities to influence
convective ones, we focus our attention on parameters $(a,b)$. Varying these
parameters can be interpreted as changing the rate of influx of reagents
through a selectively permeable boundary, or the rate of production of the
reagents from precursors through boundary-bound reactions.

We perform several steps to classify each $(a,b)$ pair in Figure 4(b,c).
First, we calculate the critical diffusion and use it to solve for the
critical decay for chemical oscillations, see Appendix B for a discussion on
chemical oscillations in the model. If the critical decay rate for
oscillations is greater than the corresponding reaction parameter,
$\eta_{Q,c}>\eta_{Q}$, it follows that the oscillatory mode is stable at the
critical diffusion ratio. However, if $\eta_{Q,c}\leq\eta_{Q}$ then it is
unstable or marginally stable at the critical diffusion. Hence, we classify
the regions where $\eta_{Q,c}<\eta_{Q}$ as oscillatory. We also provide a
snapshot of the patterns formed for supercritical diffusion ratios in each
example in Figure 4(a,d).

From Figure 4, we note that there are three regions of the $a$–$b$ parameter
space where we respectively expect stable unpatterned states, oscillatory
behaviour, and steady patterned states to emerge. Thus, we have control over
the chemical instability by varying the influxes at the boundaries, $(a,b)$.
Interestingly, oscillations are only expected to occur for small values of $a$
and for $b$ within a bounded interval that depends on $a$. Steady patterns can
emerge for larger values of $a$ and $b$, provided that the diffusion ratio $D$
is large enough. Below a certain value of $D$, the oscillatory mode is
unstable whenever the patterned state is marginally stable so we do not expect
steady patterns. Comparing the behaviour of the full and separated reaction
systems [Figure 4(b) and (c), respectively], we observe that for any value of
$D$, a larger region of $a$–$b$ parameter space supports oscillations or
patterns. Apart from scaling, the regions of parameter space are qualitatively
the same for the full and separated systems.

##  IV Buoyancy-Driven Instabilities

If the Rayleigh number is not negligible, then we cannot omit the momentum
equation in (6) and must consider the full expression in (10). We again search
for marginally stable states ($\lambda=0$) with $\eta_{R}=0$ so that (10)
reduces to:

$\displaystyle\left(\diff[2]{}{z}-(\eta_{Q}+k^{2})\right)q$
$\displaystyle=w\diff{Q_{0}}{z},$ (30a)
$\displaystyle\left(\diff[2]{}{z}-k^{2}\right)r$
$\displaystyle=w\diff{R_{0}}{z},$ (30b)
$\displaystyle\left(\diff[2]{}{z}-k^{2}\right)^{2}w$
$\displaystyle=k^{2}\mathrm{Ra}(\delta^{Q}q+\delta^{R}r).$ (30c)

A value of $\mathrm{Ra}$ that satisfies (30) is referred to as a critical
Rayleigh numnber for the wavenumber $k$, denoted $\mathrm{Ra}_{c}(k)$. We
solve for $\mathrm{Ra}_{c}(k)$ using a shooting method [29] applied to (30).
The boundary conditions are obtained by calculating (11) for both the full and
separated reactions. From our numerical results, we find that the function
$\mathrm{Ra}_{c}(k)$ has at most one local maximum [see Appendix C, Figure
8(a)], which we refer to as the maximum critical Rayleigh number and denote by
$\mathrm{Ra}_{\max,c}$, and at most one local minimum, denoted
$\mathrm{Ra}_{\min,c}$ and referred to as the minimum critical Rayleigh
number. We expect the unpatterned state to become unstable and transition to
steady patterns with finite wavenumber as the Rayleigh number increases above
the minimum critical Rayleigh number or decreases below the maximum critical
Rayleigh number. The existence and values of these minima and maxima depend on
parameters such as $\delta^{Q}$, $\delta^{R}$, $D$ and $\gamma$. We discuss
the performance of the estimates obtained by this shooting method procedure in
Appendix C.

In Section II we introduced the notion of $Q$-driven flow with $\delta^{Q}=1$.
This corresponds to when $Q$ has a larger effect on the buoyancy force at
similar concentrations compared to $R$ at similar concentrations. Letting
$\delta^{R}=0$, we find that there is a maximum but no minimum critical
Rayleigh number. The maximum, $\mathrm{Ra}_{\max,c}^{Q}$, increases with $D$.
It is negative, then zero, and positive when the diffusion ratio is
respectively below, equal to, and above the minimum critical diffusion ratio,
$D_{\min,c}:=\min_{k}D_{c}(k)$ [Figure 5(a)]. Note that $D_{\min,c}$ is the
diffusion ratio at which the system is marginally stable without flow. When
$D<D_{\min,c}$, the unpatterned state is stable without flow but a finite
wavenumber can be brought to marginal stability with flow if gravity acts
downward (i.e., the Rayleigh number is negative). In this case, coupling with
fluid flow serves to destabilize the unpatterned state. Conversely, when
$D>D_{\min,c}$, the system without flow is unstable but can be stabilized by
flow if gravity acts upwards.

The effect of fluid flow is reversed when $\delta^{Q}=0,\delta^{R}=1$; we find
that there is a minimum critical Rayleigh number, $\mathrm{Ra}_{\min,c}^{R}$,
that is positive when $D<D_{\min,c}$ and negative when $D>D_{\min,c}$ so flow
is destabilizing below the minimum critical diffusion ratio if gravity acts
upwards and flow is stabilizing above the minimum critical diffusion ratio if
gravity acts downwards. The fact that instabilities can be induced by buoyancy
effects when $\delta^{Q}=0$ and $\delta^{R}=1$ is interesting because the
concentration field $R$, and hence the fluid density, is spatially uniform in
the unpatterned steady state. This contrasts with Rayleigh–Bénard
instabilities, in which the instability occurs when the vertical gradient in
density at the steady state exceeds a threshold.

The stabilizing and destabilizing effects of $Q$ and $R$ combine when
$\delta^{Q}$ and $\delta^{R}$ are both nonzero, as shown in Figure 5(b). For
instance, given a fixed diffusion ratio $D<D_{\min,c}$, the maximum critical
Rayleigh number for $\delta^{Q}=1$ becomes more negative as $\delta^{R}$
increases because the stabilizing effect of $R$ (when gravity acts downwards)
offsets the destabilizing effect of $Q$.

Figure 5:  (a) Comparison of $\mathrm{Ra}_{\max,c}^{Q}$ for $Q$-driven flow
$(\delta^{Q},\ \delta^{R})=(1,0)$ and $\mathrm{Ra}_{\min,c}^{R}$ for
$R$-driven flow $(\delta^{Q},\ \delta^{R})=(0,1)$ for various diffusion
ratios. (b) The variation of $\mathrm{Ra}_{\min,c}^{Q}$ with $D$ and
$\delta^{R}$ for the full reaction system. Reaction parameters are held
constant at $a=0.2$, $b=2$, $\eta_{Q}=1$ and $\gamma=1$. (c)
$\mathrm{Ra}_{\min,c}^{Q}$ as a function of $D_{\min,c}(a,b)$ for various
$\delta^{R}$. We take $(a,b)$ along the line between $(0.268,2.085)$ and
$(0.352,2.19)$ and $(\gamma,\eta_{Q},D)=(1,1,7)$. The dashed line indicates
where the system transitions from pattern forming to non-pattern forming
without buoyancy forces. (d) Critical wavenumber for convective cells as a
function of $\gamma$. We fix $(a,b,\eta_{Q})=(0.2,2,1)$.

In Section III, we discussed how travelling along paths in $(a,b)$ allow us to
control the various chemical instabilities. The relationship between
$\mathrm{Ra}_{c}$ and $D_{\min,c}$ makes it possible to control the convective
instability not only by changing gradients, but also by varying $D_{\min,c}$
through the parameters $(a,b)$. In Figure 5(c), we do precisely this along a
line of $(a,b)$ values and plot it for various $\delta^{R}$ with a fixed $D$.
This plot illustrates the relationship between $D_{\min,c}$ and $\mathrm{Ra}$
and demonstrates how $\delta^{R}$ makes the system more or less sensitive to
changes in $D_{\min,c}(a,b)$. This effect could be predicted from Figure
5(a,b) since the corresponding $R$-driven flow is stable or unstable depending
on whether we are left or right of $D_{\min,c}$.

Whereas for the classical Rayleigh-Bénard instability the critical wavenumber
is fixed, for Turing pattern formation, it can be tuned by the reaction rate.
This parameter corresponds to our $\gamma$, which can be interpreted as the
enzyme density on the reactive plate. Our model shares this characteristic
that the critical wavenumber depends on $\gamma$, as shown in Figure 5(d).
Similar to classical Turing pattern formation, we find a positive correlation
between $\gamma$ and the critical wavenumber.

When we consider the separated reaction, we see other interesting behaviours.
Focusing on $Q$-driven flows, in Figure 6(a) we plot the critical Rayleigh
number, $\mathrm{Ra}_{c}^{Q}(k)$, for various wavenumbers and diffusion
ratios. We find that it is possible for both minimum and maximum critical
Rayleigh numbers to exist, with different associated critical wavenumbers. For
example, when $D=D_{c}/2$, there is a positive minimum critical Rayleigh
number and a negative maximum critical Rayleigh number. Between these Rayleigh
numbers, the unpatterned state is stable and outside this interval, patterns
emerge, illustrated in Figure 6(b)–(c). As demonstrated, different modes can
be made unstable depending on the orientation of the domain relative to the
direction of gravity. We suspect that this behaviour occurs because there is a
local minimum of $Q$ in the bulk at the steady state [see Figure 2(b)]. As a
result, the concentration gradient changes sign within the domain, allowing it
to become unstable regardless of the direction of gravity. We noted above,
however, that concentration gradients at the steady state are not a
requirement for instabilities.

Figure 6: (a) Critical Rayleigh number for the separated reaction with various
diffusion ratios, reaction parameters are $a=0.2$, $b=2.0$, $\eta_{Q}=1$,
$\gamma=1$ subject to $Q$-driven flow $(\delta^{Q},\delta^{R})=(1,0)$. The
critical diffusion ratio is $D_{c}\approx 4.02$. The wavenumbers and Rayleigh
numbers in (b) and (c) are marked with open circles. (b) Simulation of
unstable system for separated reaction with negative Rayleigh number,
$\mathrm{Ra}=-4600$, marked in (a) and $D=D_{c}/2\approx 2.01$. (c) Simulation
of unstable system for separated reaction with positive Rayleigh number
$\mathrm{Ra}=5700$ marked in (a) and $D=D_{c}/2\approx 2.01$.

##  V Conclusions

We presented boundary-bound reactions and their capacity for chemical and
convective pattern formation. We showed that the two types of instabilities
are coupled such that buoyancy effects can induce chemical patterns and
chemical patterns induce convective flows. For specific cases, we determined
necessary conditions for pattern formation in terms of the reaction
parameters. We illustrated these results by considering two examples, both
based on the Schnakenberg–Selkov reaction system. The critical diffusion ratio
was reduced by increasing the bulk-decay $\eta_{Q}$, similar to results for
active-membrane models. We further discussed how the instabilities could be
controlled through the constant boundary flux terms $(a,b)$.

Since the chemical instabilities were a sub-case of the buoyancy-driven
instabilities, this provided a natural transition for discussing mechanisms
for controlling the convective instabilities that form. We also discussed how
the parameter $\gamma$, which in practice could be easily tuned through enzyme
concentrations or temperature, can be used to change the unstable wavenumber
for the convective pattern. Moreover, we observed how the sensitivity to
changes in $D_{\min,c}(a,b)$ is influenced by the ratio
$\delta^{Q}/\delta^{R}$.

Unique to boundary-bound reactions are the non-linear boundary conditions that
influence the buoyancy of the fluid. These non-linearities were essential for
producing behaviour like the buoyancy-driven instability without steady state
density gradients. The use of chemical reagents also made it possible to
present a system where the formation of convective cells took place
independent of the direction of gravity. We presented this as a system which
can be used for detecting the orientation of the device. For instance, in
Figure 6(b,c) the direction of flow near the peaks and troughs changes
depending on the direction of gravity. By measuring the flow field near the
boundary as well as the concentration, we can determine which instability is
present and as a consequence its orientation. The unstable wavenumbers are
also different depending on the direction of gravity. This responsiveness to
orientation means it may be possible to design microfluidic devices which
change their behaviour by a simple act of turning them upside down, rather
than stopping altogether.

We opted to consider reaction parameters that can be thought to physically
represent the amount of material pumped into the domain, or the rate of the
reaction at the boundary. The latter is often tuned by changing the enzyme
density on the surface. In both cases, we demonstrate how one might design
systems where the flow fields are controlled via the reaction parameters.
Recent work in the design of reaction–diffusion programs provides some freedom
in reactions that occur along the boundary. By designing reactions with
properties of interest, such as traveling waves or other pattern forming
instabilities many other phenomena should be possible. We believe that
boundary-bound reactions can facilitate the design of microfluidic devices
with novel complex and self-regulating functionality.

###### Acknowledgements.

We acknowledge the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC), [funding reference number RGPIN-2018-04418]. Cette
recherche a été financée par le Conseil de recherches en sciences naturelles
et en génie du Canada (CRSNG), [numéro de référence RGPIN-2018-04418].

##  Appendix A Reduced Equations for Boundary-Bound Reactions

Substituting Equation (8) into Equation (6) we collect the first order terms:

$\displaystyle\diffp{Q_{1}}{t}$
$\displaystyle=-\eta_{Q}Q_{1}+\nabla^{2}Q_{1}-w\diff{Q_{0}}{z},$ (31a)
$\displaystyle\diffp{R_{1}}{t}$
$\displaystyle=-\eta_{R}R_{1}+\nabla^{2}R_{1}-w\diff{Q_{0}}{z},$ (31b)
$\displaystyle\diffp{\mathbf{u}_{1}}{t}$ $\displaystyle=-\nabla
p_{1}+\mathrm{Sc}\nabla^{2}\mathbf{u}_{1}+\mathrm{Sc}^{2}\mathrm{Gr}\left(\delta^{Q}Q_{1}+\delta^{R}R_{1}\right)\mathbf{e}_{z}.$
(31c)

The only term driving fluid flow is the buoyancy force. Therefore, the
stability of (31) is dependent on the stability of the reduced system
containing only the vertical velocity component. Omitting the horizontal flow
components, if this reduced system is unstable, then the global system will be
as well. To eliminate the pressure term, we consider the divergence of both
sides of Equation (31c). By applying the incompressibility condition, we
obtain the equations

$\nabla^{2}p_{1}=\mathrm{Sc}^{2}\mathrm{Gr}\left(\delta^{Q}\diffp{Q_{1}}{z}+\delta^{R}\diffp{R_{1}}{z}\right).$
(32)

Differentiating with respect to $z$, we obtain

$\nabla^{2}\diffp{p_{1}}{z}=\mathrm{Sc}^{2}\mathrm{Gr}\left(\delta^{Q}\diffp[2]{Q_{1}}{z}+\delta^{R}\diffp[2]{R_{1}}{z}\right).$
(33)

We now have an equation for the Laplacian of $\diffp{p_{1}}{z}$ which we can
use to eliminate the pressure term. In particular, by taking the Laplacian of
the vertical ($z$) component of (31c), we obtain

$\diffp{}{t}(\nabla^{2}u_{1,z})=-\nabla^{2}\diffp{p_{1}}{z}+\mathrm{Sc}\nabla^{4}u_{1,z}+\mathrm{Sc}^{2}\mathrm{Gr}\nabla^{2}\left(\delta^{Q}Q_{1}+\delta^{R}R_{1}\right).\\\
$ (34)

Substituting (33) into this gives

$\diffp{}{t}(\nabla^{2}u_{1,z})=\mathrm{Sc}\nabla^{4}u_{1,z}+\mathrm{Sc}^{2}\mathrm{Gr}\nabla_{H}^{2}\left(\delta^{Q}Q_{1}+\delta^{R}R_{1}\right),$
(35)

where $\nabla^{2}_{H}=\diffp[2]{}{x}+\diffp[2]{}{y}$ is the horizontal Laplace
operator. After substituting the ansatz for the perturbation we can recover
(10).

##  Appendix B On Oscillatory Instabilities

Figure 7:  Critical values of the decay rate $\eta_{Q}$ for the full reaction,
where long wavelength chemical oscillations become marginally stable as a
function of $(\gamma,D)$ with $(a,b)=(0.2,2)$ fixed.

In Turing’s analysis [1], he was able to determine conditions where chemical
oscillations are not possible. Our analysis does not address the stability of
the horizontally uniform mode $k=0$. However, long-wavelength chemical
oscillations are generally possible and should be expected since they are
present in the classical Schnakenberg–Selkov dynamics. To show their
existence, we suppose that the height of the fluid chamber is small compared
with other length scales and consider the vertically averaged version of (22)
and (24), obtaining for both:

$\displaystyle\diffp{\overline{Q}}{t}$
$\displaystyle=\gamma(a-\overline{Q}+\overline{Q}^{2}\overline{R})-\eta_{Q}\overline{Q}+\nabla^{2}_{H}\overline{Q},$
(36a) $\displaystyle\diffp{\overline{R}}{t}$
$\displaystyle=\gamma(b-\overline{Q}^{2}\overline{R})+D\nabla^{2}_{H}\overline{R}.$
(36b)

The homogeneous steady states are given by

$\displaystyle\overline{Q}_{s}$
$\displaystyle=\frac{\gamma(a+b)}{\gamma+\eta_{Q}},$ (37a)
$\displaystyle\overline{R}_{s}$
$\displaystyle=\frac{b(\gamma+\eta_{Q})^{2}}{\gamma^{2}(a+b)^{2}}.$ (37b)

Linearizing this equation at the steady state we consider a long-wavelength
oscillatory instability. The Jacobian for the system is,

$J=\begin{bmatrix}-(\gamma+\eta_{Q})+2\gamma\overline{Q}_{s}\overline{R}_{s}&\gamma\overline{Q}_{s}^{2}\\\
-2\gamma\overline{Q}_{s}\overline{R}_{s}&-\gamma\overline{Q}_{s}^{2}\end{bmatrix}.$
(38)

A quick calculation shows that $\det(J)>0$ when $\gamma>0$. Hence, a Hopf
bifurcaton can occur when $\text{Tr}(J)$ passes through zero [30]. The trace
is

$\text{Tr}(J)=\frac{2b(\eta_{Q}+\gamma)}{(a+b)}-\frac{\gamma^{3}(a+b)^{2}}{(\eta_{Q}+\gamma)^{2}}-(\eta_{Q}+\gamma),$
(39)

which is zero when $b>a$ and the decay rate $\eta_{Q}$ takes the critical
value

$\eta_{Q,c}=\gamma\left(\frac{a+b}{\sqrt[3]{b-a}}-1\right).$ (40)

When $\eta_{Q}>\eta_{Q,c}$, we find $\text{Tr}(J)>0$ yielding an unstable
focus. This implies the existence of chemical oscillations, which we would
also expect to occur in the unaveraged system.

In Figure 7 we consider the critical decay for chemical oscillations in the
full reaction. The estimates are obtained numerically by setting $k=0$ and
assuming the form $\lambda=i\omega$ for some real $\omega$ in (10), without
the fluid momentum equation. We optimize over the $(\eta_{Q},\omega)$-plane to
find the critical values of $\eta_{Q}$ as a function of $(\gamma,D)$. The plot
shows that $\gamma$ increases the critical value of $\eta_{Q}$, which we could
infer from (40), and simultaneously shows that $D$ decreases the critical
value, something that was not captured by the above vertically averaged
analysis.

##  Appendix C Comparison of Numerical and Analytical Estimates for Critical
Values

We implement a bisection algorithm comparing our estimates to simulations in
order to test the critical diffusion and critical Rayleigh number estimates
for the full reaction. Each simulation is done using Dedalus [31], a pseudo-
spectral method to solve (2), with a third order semi-implicit backwards
difference method [32]. To define initial bounds, we initialize the bisection
algorithm with a $1\%$ perturbation above and below the critical values, or
$\pm 1$ when the critical value is zero; we verify that the unpatterned steady
states are unstable and stable, respectively, at the two bracketing points and
widen the window as needed. To determine whether patterns have formed, we
calculate the amplitude of the pattern at the boundary and use it to implement
stopping conditions. The simulation is stopped if the rate of change in the
amplitude is sufficiently small relative to its size or the amplitude is
smaller than a cutoff, both of which were on the order of $10^{-8}$ or
smaller. Additionally, a maximum time was set to end simulations if neither of
these conditions were reached. Once the simulation was finished, we compared
the final amplitude of the pattern to our cutoff; we classify the system as
unstable if it is larger and use this to update the bounds. This means the
algorithm generally overestimates the critical parameters slightly, as it
cannot distinguish between stable states and patterns with amplitudes below
the cutoff. We bisect the interval five times, to obtain an accuracy of $\sim
0.156\%$ assuming that the upper and lower bounds were accepted. We report the
final midpoint values and corresponding deviations from the estimates obtained
analytically or via the shooting method in Tables 1 and 2 . The similarity
between the two critical estimates suggests that our calculations are
accurate.

Figure 8:  Critical Rayleigh numbers and diffusion ratios for the full reaction system with $Q$-driven flow $(\delta^{Q},\delta^{R})=(1,0)$. Except where explicitly varied, reaction parameters are $(a,b,\eta_{Q},\gamma)=(0.2,2,1,1)$. (a) Critical Rayleigh number as functions of wavenumber for various diffusion ratios between $1$ and $D_{c}\approx 4.82$, obtained by a shooting method. (b) Comparison of the maximum critical Rayleigh number obtained by determining the local maxima of the functions in (a) with estimates from the bisection method. (c) Critical diffusion estimates vs. $\eta_{Q}$, where $\mathrm{Ra}=0$, compared with estimates from the bisection method. Decay |  Diffusion Estimate |  Bisection Algorithm |  Relative Error
---|---|---|---
0.01 | 38.472 | 38.653 | 0.469%
0.1 | 28.514 | 28.558 | 0.156%
0.2 | 21.286 | 21.319 | 0.156%
0.3 | 16.445 | 16.419 | 0.156%
0.4 | 13.065 | 13.086 | 0.156%
0.5 | 10.624 | 10.641 | 0.156%
0.6 | 8.809 | 8.823 | 0.156%
0.7 | 7.428 | 7.439 | 0.156%
0.8 | 6.353 | 6.363 | 0.156%
0.9 | 5.022 | 5.511 | 0.156%
1.0 | 4.818 | 4.825 | 0.156%
Table 1: Comparison of an analytical estimate and the bisection algorithm for the minimum critical diffusion for the full reaction. Reaction parameters used are $(a,b,\gamma)=(0.2,2,1)$. Results are rounded to three decimal places. Diffusion |  Estimate |  Bisection Algorithm |  Absolute Error
---|---|---|---
$1$ | $-1999.659$ | $-1996.534$ | $3.124$
$1.424$ | $-1767.331$ | $-1765.570$ | $2.761$
$1.848$ | $-1546.015$ | $-1543.599$ | $2.415$
$2.273$ | $-1328.441$ | $-1326.365$ | $2.075$
$2.697$ | $-1111.122$ | $-1109.386$ | $1.736$
$3.121$ | $-892.542$ | $-891.148$ | $1.394$
$3.545$ | $-672.156$ | $-671.106$ | $1.050$
$3.969$ | $-449.864$ | $-449.162$ | $0.702$
$4.393$ | $-225.758$ | $-225.406$ | $0.352$
$4.818$ | $0$ | $0.044$ | $0.044$
Table 2:  Comparison of the shooting method estimate and the bisection
algorithm for $\mathrm{Ra}_{\max,c}^{Q}$ for the full reaction with $Q$-driven
flow with various diffusion ratios $D$ and fixed reaction parameters
$(a,b,\eta_{Q},\gamma)=(0.2,2.0,1,1)$. Results are rounded to three decimal
places.

## References

  * Turing [1952] A. M. Turing, The chemical basis of morphogenesis, Phil. Trans. R. Soc. B 237, 37 (1952).
  * Satnoianu  _et al._ [2000] R. A. Satnoianu, M. Menzinger, and P. K. Maini, Turing instabilities in general systems, J. Math. Biol. 41, 493 (2000).
  * Levine and Rappel [2005] H. Levine and W.-J. Rappel, Membrane-bound Turing patterns, Phys. Rev. E 72, 061912 (2005).
  * Gomez  _et al._ [2021] D. Gomez, S. Iyaniwura, F. Paquin-Lefebvre, and M. J. Ward, Pattern forming systems coupling linear bulk diffusion to dynamically active membranes or cells, Phil. Trans. R. Soc. A. 379, 20200276 (2021).
  * Madzvamuse  _et al._ [2015] A. Madzvamuse, A. H. W. Chung, and C. Venkataraman, Stability analysis and simulations of coupled bulk-surface reaction–diffusion systems, Proc. R. Soc. A 471, 20140546 (2015).
  * Shklyaev  _et al._ [2020] O. E. Shklyaev, V. V. Yashin, S. I. Stupp, and A. C. Balazs, Enhancement of chemical oscillations by self-generated convective flows, Commun. Phys. 3, 70 (2020).
  * Castets  _et al._ [1990] V. Castets, E. Dulos, J. Boissonade, and P. De Kepper, Experimental evidence of a sustained standing Turing-type nonequilibrium chemical pattern, Phys. Rev. Lett. 64, 2953 (1990).
  * Lengyel and Epstein [1993] I. Lengyel and I. R. Epstein, Turing structures in simple chemical reactions, Acc. Chem. Res. 26, 235 (1993).
  * Vasquez  _et al._ [1993] D. A. Vasquez, J. W. Wilder, and B. F. Edwards, Convective Turing patterns, Phys. Rev. Lett. 71, 1538 (1993).
  * Diewald and Brand [1995] M. Diewald and H. R. Brand, Chemically driven convection can stabilize Turing patterns, Phys. Rev. E 51, R5200 (1995).
  * McTaggart and Straughan [1986] C. L. McTaggart and B. Straughan, Chemical surface reactions and nonlinear stability by the method of energy, SIAM J. Math. Anal. 17, 342 (1986).
  * De Wit  _et al._ [2012] A. De Wit, K. Eckert, and S. Kalliadasis, Introduction to the focus issue: Chemo-hydrodynamic patterns and instabilities, Chaos 22, 037101 (2012).
  * Bdzil and Frisch [2008] J. B. Bdzil and H. L. Frisch, Chemically driven convection, J. Chem. Phys. 72, 1875 (2008).
  * Miike  _et al._ [1989] H. Miike, S. C. Müller, and B. Hess, Hydrodynamic flows traveling with chemical waves, Phys. Lett. A 141, 25 (1989).
  * Bees  _et al._ [2001] M. A. Bees, A. J. Pons, P. G. Sørensen, and F. Sagués, Chemoconvection: A chemically driven hydrodynamic instability, J. Chem. Phys. 114, 1932 (2001).
  * Sengupta  _et al._ [2014] S. Sengupta, D. Patra, I. Ortiz-rivera, A. Agrawal, S. Shklyaev, K. K. Dey, U. Córdova-figueroa, T. E. Mallouk, and A. Sen, Self-powered enzyme micropumps, Nature Chem. 6, 415 (2014).
  * Ortiz-Rivera  _et al._ [2016] I. Ortiz-Rivera, H. Shum, A. Agrawal, A. Sen, and A. C. Balazs, Convective flow reversal in self-powered enzyme micropumps, Proc. Natl. Acad. Sci. U.S.A. 113, 2585 (2016).
  * Valdez  _et al._ [2017] L. Valdez, H. Shum, I. Ortiz-Rivera, A. C. Balazs, and A. Sen, Solutal and thermal buoyancy effects in self-powered phosphatase micropumps, Soft Matter 13, 2800 (2017).
  * Shklyaev  _et al._ [2018] O. E. Shklyaev, H. Shum, and A. C. Balazs, Using chemical pumps and motors to design flows for directed particle assembly, Acc. Chem. Res. 51, 2672 (2018).
  * Mathew  _et al._ [2007] G. Mathew, I. Mezić, S. Grivopoulos, U. Vaidya, and L. Petzold, Optimal control of mixing in Stokes fluid flows, J. Fluid Mech. 580, 261–281 (2007).
  * Das  _et al._ [2017] S. Das, O. E. Shklyaev, A. Altemose, H. Shum, I. Ortiz-Rivera, L. Valdez, T. E. Mallouk, A. C. Balazs, and A. Sen, Harnessing catalytic pumps for directional delivery of microparticles in microchambers, Nat. Commun. 8, 14384 (2017).
  * Radko [2013] T. Radko,  _Double-diffusive convection_ (Cambridge University Press, 2013).
  * Kundu and Cohen [2002] P. Kundu and I. Cohen,  _Fluid Mechanics_ , 2nd ed. (Academic Press, 2002).
  * Scalise and Schulman [2014] D. Scalise and R. Schulman, Designing modular reaction-diffusion programs for complex pattern formation, Technology 02, 55 (2014).
  * Soloveichik  _et al._ [2010] D. Soloveichik, G. Seelig, and E. Winfree, DNA as a universal substrate for chemical kinetics, Proc. Natl. Acad. Sci. U.S.A. 107, 5393 (2010).
  * Chen and Seelig [2020] S. Chen and G. Seelig, Programmable patterns in a DNA-based reaction–diffusion system, Soft Matter 16, 3555 (2020).
  * Schnakenberg [1979] J. Schnakenberg, Simple chemical reaction systems with limit cycle behaviour, J. Theor. Biol. 81, 389 (1979).
  * Sel’kov [1968] E. Sel’kov, Self-oscillations in glycolysis 1. A simple kinetic model, Eur. J. Biochem. 4, 79 (1968).
  * Bulirsch and Stoer [2002] R. Bulirsch and J. Stoer,  _Introduction to numerical analysis_, 2nd ed., Vol. 3 (Springer New York, 2002).
  * Ricard and Mischler [2009] R. Ricard and S. Mischler, Turing instabilities at Hopf bifurcation, J. Nonlinear Sci. 19, 467 (2009).
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Research 2, 023068 (2020).
  * Ruuth [1995] S. Ruuth, Implicit-explicit methods for reaction-diffusion problems in pattern formation, J. Math. Biol. 34, 148 (1995).


\frac{e^{-i\mathbf{q}\cdot\boldsymbol{\tau}_{\kappa}}}{\sqrt{\Omega}}\Big{[}\sum_{\beta}iq_{\beta}Z_{\kappa\alpha\beta}+\sum_{\gamma}\frac{q_{\beta}q_{\gamma}}{2}Q_{\kappa\alpha\beta\gamma}+\cdots\Big{]},$
(88)

while the screened kernel is:

$\tilde{W}_{\mathbf{q}}^{\mathcal{L}}(\mathbf{0,0})=\frac{4\pi}{|\mathbf{q}|^{2}}f(\mathbf{q})\left[1-\frac{4\pi
e^{2}}{|\mathbf{q}|^{2}}f(\mathbf{q})\tilde{\chi}^{\mathcal{S}}_{\mathbf{q}}\right]^{-1}$
(89)

where the short range polarizability is analytic in $\mathbf{q}$ and can be
related to the electronic dielectric tensor $\varepsilon_{\alpha\beta}$ in the
long-wavelength limit through:

$\displaystyle\lim_{\mathbf{q}\to\mathbf{0}}4\pi
e^{2}\tilde{\chi}_{\mathbf{q}}^{\mathcal{S}}$
$\displaystyle=-\sum_{\alpha,\beta}{q}_{\alpha}\left(\varepsilon_{\alpha\beta}-\delta_{\alpha\beta}\right){q}_{\beta}+\mathcal{O}(\mathbf{q}^{4})$
(90)

Following similar steps as in the 2D case, we then get the final expression:

$g_{mn\nu}^{\mathcal{L}}(\mathbf{k},\mathbf{q})=\sum_{\kappa\alpha}\bigg{[}\frac{\hbar}{2M_{\kappa}\omega_{\nu}(\mathbf{q})\Omega^{2}}\bigg{]}^{\frac{1}{2}}\frac{4\pi
ef(|\mathbf{q}|)}{|\mathbf{q}|^{2}\tilde{\epsilon}(\mathbf{q})}\\\ \times
e^{-i\mathbf{q}\cdot\boldsymbol{\tau}_{\kappa}}\Big{[}\sum_{\beta}iq_{\beta}Z_{\kappa\alpha\beta}+\sum_{\gamma}\frac{q_{\beta}q_{\gamma}}{2}Q_{\kappa\alpha\beta\gamma}\Big{]}e_{\kappa\alpha\nu}(\mathbf{q})\\\
\bra{\Psi_{m\mathbf{k+q}}}e^{i\mathbf{q}\cdot\mathbf{r}}[1+iq_{\alpha}V^{{\rm
Hxc},{\cal E}_{\alpha}}(\mathbf{r})]\ket{\Psi_{n\mathbf{k}}},$ (91)

where the effective macroscopic dielectric function reads:

$\tilde{\epsilon}(\mathbf{q})=\frac{\mathbf{q}\cdot\boldsymbol{\varepsilon}\cdot\mathbf{q}}{\phantom{{}^{2}}|\mathbf{q}|^{2}}f(|\mathbf{q}|)+1-f(|\mathbf{q}|)~{}.$
(92)

Figure 16:  Phonon dispersion of phosphorene (a-c) and graphene (d-f)
calculated with direct density functional perturbation theory calculations
(black empty circles) compared with Fourier interpolation where the long-range
part of the dynamical matrix (D) includes dipole-dipole (DD), dipole-
quadrupole (DQ) and quadrupole-quadrupole (QQ) and where the yellow line
refers to the scheme of Ref. Sohier et al., 2017b.

##  Appendix D Wannier interpolation

We now report for each material the choice of initial projections for the
Wannier functions as well as interpolated electronic band structures compared
with direct DFT calculations.

For SnS2, we choose to focus on the conduction band only since this band is
separated in energy. We perform the interpolation using maximally localized
Wannier functions (MLWFs) Marzari et al. (2012) as implemented in the
wannier90 software Pizzi et al. (2020), starting from an initial projection
made on the Sn atom with $s$ character and yielding an atom-centred Wannier
function with spread 9.32 Å2. The resulting electronic band structure agrees
with the DFT one as seen in Fig. 4(a).

In the case of h-BN, two sets of MLWFs are considered separately for the
valence (blue) and conduction bands (red), as shown in Fig. 4(d). For the
valance bands, the initial projection is on the N atom with $p$ character and
yields two in-plane Wannier functions with spread 1.02 Å2 and one out-of-plane
with spread 1.35 Å2. For the conduction band, the initial projections are
$p$-orbitals on both B and N atoms, with a resulting average spread of 5.12
Å2.

For MoS2, we wannierize together the top 2 valence bands and the lowest 8
conduction bands by initially projecting on Mo $d$ orbitals for both spin
polarizations, yielding Wannier functions with spread 4.9 Å2, 5.0 Å2, and 5.1
Å2. The corresponding interpolated band structure is shown in Fig. 4(b).

As far as InSe is concerned, we address separately valence and conduction
bands. For the valance bands (blue solid lines in Fig. 4(e)), the initial
projection is on one of the two Se atoms with $p_{z}$ character and yields two
spin-degenerate Wannier functions with spread 15.79 Å2. For the conduction
bands, the initial projection is an $s$-orbital (for each spin) on one of the
two In atoms and gives two Wannier functions with spread of 11.37 Å2.

For graphene, a single set of 10 MLWFs is created to interpolate both valence
and conduction bands, as shown in Fig. 4(c) with blue lines. The initial
Wannier projections consist of hybrid $sp^{2}$ orbitals on one carbon atom and
$p_{z}$ orbitals on both carbon atoms, resulting in 6 (including spin
degeneracy) MLWFs with spread 0.61 Å2 and 4 with spread 0.98 Å2.

Finally, for phosphorene we consider valence and conduction bands together by
constructing 32 (including spin) MLWFs that accurately reproduce DFT results
around the band gap, as shown in Fig. 4(f). As initial projections on each of
the four P atoms we take $p$\- and $s$-orbitals, giving 8 MLWFs with spread
1.92 Å2, 16 MLWFs with spread 2.06 Å2, and 8 MLWFs with spread 2.36 Å2.

##  Appendix E Additional phonon dispersion and mobility results

For completeness, we report here the figures showing the DFPT phonon
dispersions and their interpolations for SnS2, MoS2 and InSe in Fig. 15, and
for phosphorene and graphene in Fig. 16. Finally, we give explicit room
temperature mobility values for various integration fine grids as well as
their extrapolated values for all the materials investigated here in Table
LABEL:table2 as well as a comparison between different long-range treatments
in Table LABEL:table3.

Table 2:  Room temperature fine grid convergence of drift and Hall mobility using adaptive smearing and dipole velocity. The symbol $\infty$ for the fine grids means that the value is obtained by extrapolating the linear dependence on inverse grid size as explained in Ref. Poncé et al., 2021. SERTA mobility refers to the mobility in the self-energy relaxation time approximation. The “-e” or “-h” suffix after the material’s name refers to electron or hole mobility, respectively. In all cases we report the mobility including quadrupoles, except for graphene. For phosphorene, results along the two principal directions ($xx$ and $yy$) are reported separately. Material | Window | Coarse grids | Fine grids | SERTA mobility | BTE mobility | BTE Hall | Hall mobility
---|---|---|---|---|---|---|---
| eV | k/q | k/q | cm2/Vs | cm2/Vs | factor | cm2/Vs
SnS2-e | 0.3 | 162 | 3002 | 16.79 | 23.52 | 0.979 | 23.03
|  |  | 4002 | 16.78 | 23.58 | 0.992 | 23.39
|  |  | 5002 | 16.68 | 23.45 | 0.999 | 23.43
|  |  | 6002 | 16.69 | 23.49 | 1.000 | 23.49
|  |  | 7002 | 16.62 | 23.42 | 1.004 | 23.51
|  |  | $\infty$ | 16.53 | 23.38 | 1.023 | 23.92
MoS2-e | 0.2 | 182 | 3002 | 116.66 | 134.39 | 1.066 | 143.26
|  |  | 4002 | 114.81 | 131.09 | 1.069 | 140.14
|  |  | 5002 | 115.95 | 132.37 | 1.069 | 141.50
|  |  | 6002 | 117.39 | 133.19 | 1.072 | 142.78
|  |  | 8002 | 116.84 | 132.75 | 1.074 | 142.57
|  |  | $\infty$ | 117.20 | 131.78 | 1.078 | 142.06
MoS2-h | 0.25 | 182 | 3002 | 60.04 | 72.70 | 1.576 | 114.58
|  |  | 4002 | 60.35 | 73.18 | 1.588 | 116.21
|  |  | 5002 | 60.55 | 73.33 | 1.594 | 116.89
|  |  | 6002 | 61.16 | 73.41 | 1.604 | 117.75
|  |  | $\infty$ | 61.97 | 74.14 | 1.629 | 120.77
BN-e | 0.4 | 162 | 2002 | 131.87 | 141.18 | 0.912 | 128.76
|  |  | 3002 | 130.64 | 143.28 | 0.903 | 129.38
|  |  | 4002 | 128.06 | 139.26 | 0.912 | 127.01
|  |  | 5002 | 127.27 | 138.87 | 0.909 | 126.23
|  |  | 6002 | 126.42 | 138.10 | 0.911 | 125.81
|  |  | $\infty$ | 124.06 | 137.08 | 0.909 | 124.61
BN-h | 0.4 | 162 | 2002 | 539.45 | 665.66 | 0.990 | 659.00
|  |  | 3002 | 501.15 | 641.54 | 1.016 | 651.80
|  |  | 4002 | 495.52 | 650.36 | 1.016 | 660.77
|  |  | 5002 | 485.26 | 640.82 | 1.015 | 650.43
|  |  | 6002 | 480.64 | 638.16 | 1.017 | 649.01
|  |  | $\infty$ | 447.65 | 626.50 | 1.017 | 637.15
InSe-e | 0.3 | 162 | 2002 | 68.90 | 122.20 | 1.085 | 132.59
|  |  | 3002 | 64.03 | 115.28 | 1.127 | 129.92
|  |  | 4002 | 63.05 | 113.07 | 1.132 | 128.00
|  |  | 5002 | 62.21 | 111.18 | 1.136 | 126.30
|  |  | 6002 | 61.67 | 110.26 | 1.139 | 125.59
|  |  | 8002 | 61.38 | 110.09 | 1.144 | 125.94
|  |  | $\infty$ | 59.70 | 106.34 | 1.153 | 122.61
InSe-h | 0.3 | 162 | 1002 | 0.42 | 0.56 | 1.534 | 0.86
|  |  | 2002 | 0.42 | 0.56 | 1.432 | 0.80
|  |  | 3002 | 0.42 | 0.56 | 1.420 | 0.80
|  |  | 4002 | 0.42 | 0.56 | 1.417 | 0.79
|  |  | $\infty$ | 0.42 | 0.56 | 1.400 | 0.78
P-e ($xx$) | 0.2 | 162 | 2002 | 31.58 | 31.44 | 1.172 | 36.85
|  |  | 3002 | 31.63 | 32.26 | 1.187 | 38.29
|  |  | 4002 | 31.69 | 32.25 | 1.193 | 38.47
|  |  | 5002 | 31.57 | 32.23 | 1.194 | 38.48
|  |  | 6002 | 31.62 | 32.22 | 1.197 | 38.57
|  |  | $\infty$ | 31.64 | 32.18 | 1.206 | 38.81
P-e ($yy$) | 0.2 | 162 | 2002 | 251.93 | 252.68 | 1.172 | 296.14
|  |  | 3002 | 249.18 | 254.97 | 1.187 | 302.65
|  |  | 4002 | 249.58 | 254.28 | 1.193 | 303.36
|  |  | 5002 | 247.70 | 253.73 | 1.194 | 302.95
|  |  | 6002 | 248.59 | 253.80 | 1.197 | 303.80
|  |  | $\infty$ | 247.23 | 252.41 | 1.206 | 304.41
P-h ($xx$) | 0.3 | 162 | 2002 | 21.36 | 26.57 | 1.107 | 29.41
|  |  | 3002 | 21.40 | 26.39 | 1.108 | 29.24
|  |  | 4002 | 21.34 | 26.40 | 1.127 | 29.75
|  |  | 5002 | 21.21 | 26.27 | 1.147 | 30.13
|  |  | 6002 | 21.15 | 26.20 | 1.148 | 30.08
|  |  | $\infty$ | 20.92 | 26.05 | 1.194 | 31.10
P-h ($yy$) | 0.3 | 162 | 2002 | 553.77 | 504.52 | 1.107 | 558.50
|  |  | 3002 | 525.16 | 491.12 | 1.108 | 544.16
|  |  | 4002 | 511.78 | 480.22 | 1.127 | 541.21
|  |  | 5002 | 494.05 | 465.65 | 1.147 | 534.10
|  |  | 6002 | 494.57 | 464.75 | 1.148 | 533.53
|  |  | 7002 | 490.12 | 460.06 | 1.156 | 531.83
|  |  | $\infty$ | 461.40 | 435.64 | 1.193 | 519.72
Graphene-e+h | 0.5 | 122 | 10002 | 1437330 | 1224730 | 3.019 | 3697460
|  |  | 16002 | 1371750 | 1160470 | 3.423 | 3972289
|  |  | 20002 | 1366710 | 1174670 | 4.005 | 4704553
|  |  | 24002 | 1308210 | 1091250 | 3.820 | 4168575
|  |  | $\infty$ | 1249560 | 1051630 | 4.549 | 4783865
Table 3:  Room temperature drift and Hall mobility using adaptive smearing and dipole velocity with and without spin-orbit coupling (SOC) included. SERTA mobility refers to the mobility in the self-energy relaxation time approximation, while BTE refers to the mobility calculated using the iterative solution of the Boltzmann transport equation. The “-e” or “-h” suffix after the material’s name refers to electron or hole mobility, respectively. The k/q-point grids reported in parenthesis refers to the coarse and fine electron and phonon momentum grids, respectively. We compare various long-range treatments where the long-range part of the dynamical matrix (D) and electron-phonon matrix elements (G) includes dipole-dipole (DDS) and monopole-dipole (eDS) with the scheme of Refs. Sohier et al., 2016, 2017b, dipole-dipole (DD), dipole-quadrupole (DQ), quadrupole-quadrupole (QQ), monopole-dipole (eD), and monopole-quadrupole (eQ). We also report if the pseudopotential includes non-linear core corrections (NLCC) or not. |  | No SOC | SOC
---|---|---|---
|  | Drift mobility | Hall mobility | Drift mobility | Hall mobility
| Long-range | SERTA | BTE | SERTA | BTE | SERTA | BTE | SERTA | BTE
Material | treatment | cm2/Vs | cm2/Vs | cm2/Vs | cm2/Vs | cm2/Vs | cm2/Vs | cm2/Vs | cm2/Vs
SnS2-e | No long-range | 31.97 | 29.74 | 33.18 | 31.00 | - | - | - | -
(162 k/q) | D(DDS) + G(eDS) | 16.35 | 23.28 | 17.59 | 23.38 | - | - | - | -
(6002 k/q) | D(DD) + G(eD) | 16.70 | 23.49 | 17.98 | 23.56 | - | - | - | -
without NLCC | D(DD+DQ+QQ) + G(eD) | 16.70 | 23.49 | 17.98 | 23.57 | - | - | - | -
| D(DD) + G(eD+eQ) | 16.68 | 23.49 | 17.95 | 23.51 | - | - | - | -
| D(DD+DQ+QQ) + G(eD+eQ) | 16.69 | 23.49 | 17.95 | 23.51 | - | - | - | -
MoS2-e | No long-range | 178.18 | 159.67 | 200.41 | 175.30 | 183.13 | 163.82 | 212.73 | 185.26
(182 k/q) | D(DDS) + G(eDS) | 157.58 | 152.71 | 178.53 | 166.51 | 160.76 | 163.55 | 187.41 | 189.14
(5002 k/q) | D(DD) + G(eD) | 158.52 | 153.03 | 179.49 | 166.88 | 161.75 | 156.35 | 188.48 | 175.08
with NLCC | D(DD+DQ+QQ) + G(eD) | 158.70 | 153.19 | 179.70 | 167.09 | 160.73 | 155.27 | 186.42 | 173.70
| D(DD) + G(eD+eQ) | 113.78 | 129.90 | 119.90 | 136.51 | 115.98 | 132.40 | 124.33 | 141.58
| D(DD+DQ+QQ) + G(eD+eQ) | 114.24 | 130.22 | 120.44 | 136.90 | 115.95 | 132.37 | 124.37 | 141.57
MoS2-h | No long-range | 14.38 | 14.18 | 34.52 | 32.65 | 119.03 | 108.08 | 431.59 | 345.76
(182 k/q) | D(DDS) + G(eDS) | 13.90 | 13.98 | 31.97 | 30.51 | 97.80 | 96.22 | 229.36 | 203.17
(5002 k/q) | D(DD) + G(eD) | 13.92 | 13.99 | 32.08 | 30.60 | 98.54 | 96.60 | 233.18 | 205.61
with NLCC | D(DD+DQ+QQ) + G(eD) | 13.95 | 14.02 | 32.07 | 30.60 | 98.90 | 96.48 | 231.06 | 203.24
| D(DD) + G(eD+eQ) | 11.63 | 13.07 | 22.10 | 22.54 | 59.84 | 72.45 | 94.22 | 94.437
| D(DD+DQ+QQ) + G(eD+eQ) | 11.62 | 13.06 | 22.19 | 22.64 | 60.55 | 73.33 | 96.49 | 117.42
BN-e |  |  |  |  |  |  |  |  |
(162 k/q) | D(DDS) + G(eDS) | 96.94 | 61.86 | 166.26 | 175.55 | - | - | - | -
(5002 k/q) | D(DD) + G(eD) | 98.32 | 62.88 | 168.55 | 178.98 | - | - | - | -
without NLCC | D(DD+DQ+QQ) + G(eD) | 127.27 | 138.87 | 100.87 | 126.20 | - | - | - | -
BN-h |  |  |  |  |  |  |  |  |
(162 k/q) | D(DDS) + G(eDS) | 303.17 | 861.28 | 299.22 | 969.31 | - | - | - | -
(5002 k/q) | D(DD) + G(eD) | 308.31 | 878.68 | 303.40 | 984.70 | - | - | - | -
without NLCC | D(DD+DQ+QQ) + G(eD+eQ) | 485.26 | 640.82 | 472.67 | 650.22 | - | - | - | -
BN-e | No long-range | 148.07 | 168.58 | 313.96 | 405.16 | 148.06 | 168.57 | 313.94 | 405.14
(162 k/q) | D(DDS) + G(eDS) | 49.48 | 96.73 | 61.10 | 132.94 | 49.48 | 96.73 | 61.10 | 132.94
(5002 k/q) | D(DD) + G(eD) | 50.64 | 98.49 | 62.52 | 137.16 | 50.64 | 98.49 | 62.52 | 137.15
with NLCC | D(DD+DQ+QQ) + G(eD) | 53.25 | 97.08 | 65.66 | 133.33 | 53.25 | 65.66 | 97.08 | 133.32
| D(DD) + G(eD+eQ) | 100.50 | 151.81 | 120.67 | 206.72 | 100.50 | 151.80 | 120.66 | 206.70
| D(DD+DQ+QQ) + G(eD+eQ) | 101.38 | 128.89 | 153.62 | 235.15 | 101.38 | 128.88 | 153.62 | 235.13
BN-h | No long-range | 1619.14 | 1369.60 | 1777.93 | 1469.98 | 1619.14 | 1369.60 | 1777.93 | 1469.98
(162 k/q) | D(DDS) + G(eDS) | 255.05 | 866.39 | 254.10 | 968.47 | 255.05 | 866.39 | 254.10 | 968.47
(5002 k/q) | D(DD) + G(eD) | 262.21 | 893.49 | 260.09 | 997.86 | 262.21 | 893.48 | 260.09 | 997.85
with NLCC | D(DD+DQ+QQ) + G(eD) | 276.68 | 903.90 | 273.30 | 1020.49 | 276.68 | 903.90 | 273.30 | 1020.49
| D(DD) + G(eD+eQ) | 613.67 | 596.56 | 738.28 | 745.53 | 613.67 | 596.56 | 738.28 | 745.53
| D(DD+DQ+QQ) + G(eD+eQ) | 636.38 | 837.63 | 617.34 | 847.13 | 636.38 | 837.62 | 617.34 | 847.12
InSe-e |  |  |  |  |  |  |  |  |
(162 k/q) | D(DDS) + G(eDS) | 56.44 | 100.58 | 71.02 | 112.47 | 57.25 | 101.95 | 72.25 | 114.28
(6002 k/q) | D(DD) + G(eD) | 58.64 | 103.50 | 73.89 | 115.69 | 59.49 | 104.94 | 75.10 | 117.52
with NLCC | D(DD+DQ+QQ) + G(eD+eQ) | 60.88 | 108.95 | 78.85 | 123.86 | 61.67 | 110.28 | 79.96 | 125.53
InSe-h |  |  |  |  |  |  |  |  |
(162 k/q) | D(DDS) + G(eDS) | 0.27 | 0.16 | 0.28 | 0.14 | 0.41 | 0.55 | 0.43 | 0.79
(4002 k/q) | D(DD) + G(eD) | 0.45 | 0.36 | 0.59 | 0.69 | 0.42 | 0.55 | 0.44 | 0.80
with NLCC | D(DD+DQ+QQ) + G(eD+eQ) | 0.46 | 0.36 | 0.59 | 0.69 | 0.42 | 0.56 | 0.44 | 0.80

###### Acknowledgements.

The authors would like to thank J.-M. Lihm, G. Brunin and X. Gonze for useful
discussions; S.P. acknowledges support from the F.R.S.-FNRS as well as from
the European Unions Horizon 2020 Research and Innovation Program, under the
Marie Skłodowska-Curie Grant Agreement (SELPH2D, No. 839217); N.M.
acknowledges support from the Swiss National Science Foundation and the NCCR
MARVEL; M.S. and M.R. acknowledge support from Ministerio de Ciencia y
Innovación (MICINN-Spain) through Grant No. PID2019-108573GB-C22; from Severo
Ochoa FUNFUTURE center of excellence (CEX2019-000917-S); from Generalitat de
Catalunya (Grant No. 2017 SGR1506); and from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(Grant Agreement No. 724529). M.G. acknowledges support from the Italian
Ministry for University and Research through the Levi-Montalcini program.
Computational resources have been provided by the PRACE-21 resources
MareNostrum at BSC-CNS, the Supercomputing Center of Galicia (CESGA), and by
the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds
de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No.
2.5020.11 and by the Walloon Region as well as computational resources awarded
on the Belgian share of the EuroHPC LUMI supercomputer.

## References

  * Wang et al. (2013) L. Wang, I. Meric, P. Y. Huang, Q. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. A. Muller, J. Guo, et al., Science 342, 614 (2013).
  * Fiori et al. (2014) G. Fiori, F. Bonaccorso, G. Iannaccone, T. Palacios, D. Neumaier, A. Seabaugh, S. K. Banerjee, and L. Colombo, Nature Nanotechnology 9, 768 (2014).
  * Lv et al. (2016) W. Lv, Z. Li, Y. Deng, Q.-H. Yang, and F. Kang, Energy Storage Materials 2, 107 (2016).
  * Chhowalla et al. (2016) M. Chhowalla, D. Jena, and H. Zhang, Nature Reviews Materials 1, 16052 (2016).
  * Kang et al. (2020) S. Kang, D. Lee, J. Kim, A. Capasso, H. S. Kang, J.-W. Park, C.-H. Lee, and G.-H. Lee, 2D Materials 7, 022003 (2020).
  * Radisavljevic et al. (2011) B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti, and A. Kis, Nat. Nanotechnol. 6, 147 (2011).
  * Manzeli et al. (2017) S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, Nature Reviews Materials 2, 17033 (2017).
  * Reyes-Retana and Cervantes-Sodi (2016) J. A. Reyes-Retana and F. Cervantes-Sodi, Sci. Rep. 6, 24093 (2016).
  * Ahn (2020) E. C. Ahn, npj 2D Materials and Applications 4, 17 (2020).
  * Sierra et al. (2021) J. F. Sierra, J. Fabian, R. K. Kawakami, S. Roche, and S. O. Valenzuela, Nature Nanotechnology 16, 856 (2021).
  * Zeng et al. (2012) H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Nat. Nanotechnol. 7, 490 (2012).
  * Schaibley et al. (2016) J. R. Schaibley, H. Yu, G. Clark, P. Rivera, J. S. Ross, K. L. Seyler, W. Yao, and X. Xu, Nature Reviews Materials 1, 16055 (2016).
  * Sohier et al. (2020) T. Sohier, M. Gibertini, and N. Marzari, 2D Materials 8, 015025 (2020).
  * Cheng et al. (2020) L. Cheng, C. Zhang, and Y. Liu, Phys. Rev. Lett. 125, 177701 (2020).
  * Gonze and Lee (1997) X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997).
  * Baroni et al. (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
  * Giustino (2017) F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
  * Sohier et al. (2017a) T. Sohier, M. Calandra, and F. Mauri, Phys. Rev. B 96, 075448 (2017a).
  * Li (2015) W. Li, Phys. Rev. B 92, 075405 (2015).
  * Zhou and Bernardi (2016) J.-J. Zhou and M. Bernardi, Phys. Rev. B 94, 201201 (2016).
  * Ma et al. (2018) J. Ma, A. S. Nissimagoudar, and W. Li, Phys. Rev. B 97, 045201 (2018).
  * Macheda and Bonini (2018) F. Macheda and N. Bonini, Phys. Rev. B 98, 201201(R) (2018).
  * Poncé et al. (2018) S. Poncé, E. R. Margine, and F. Giustino, Phys. Rev. B 97, 121201 (2018).
  * Poncé et al. (2019a) S. Poncé, D. Jena, and F. Giustino, Phys. Rev. B 100, 085204 (2019a).
  * Li et al. (2019) W. Li, S. Poncé, and F. Giustino, Nano Letters 19, 1774 (2019).
  * Poncé et al. (2019b) S. Poncé, D. Jena, and F. Giustino, Phys. Rev. Lett. 123, 096602 (2019b).
  * Macheda et al. (2020) F. Macheda, S. Poncé, F. Giustino, and N. Bonini, Nano Letters 20, 8861 (2020).
  * Protik and Broido (2020) N. H. Protik and D. A. Broido, Phys. Rev. B 101, 075202 (2020).
  * Jalil et al. (2020) O. Jalil, S. Ahmad, X. Liu, K. W. Ang, and U. Younis, EPL (Europhysics Letters) 130, 57001 (2020).
  * D’Souza et al. (2020) R. D’Souza, J. Cao, J. D. Querales-Flores, S. Fahy, and I. Savić, Phys. Rev. B 102, 115204 (2020).
  * Brunin et al. (2020a) G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Phys. Rev. B 102, 094308 (2020a).
  * Brunin et al. (2020b) G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Phys. Rev. Lett. 125, 136601 (2020b).
  * Poncé and Giustino (2020) S. Poncé and F. Giustino, Phys. Rev. Research 2, 033102 (2020).
  * Jhalani et al. (2020) V. A. Jhalani, J.-J. Zhou, J. Park, C. E. Dreyer, and M. Bernardi, Phys. Rev. Lett. 125, 136602 (2020).
  * Park et al. (2020) J. Park, J.-J. Zhou, V. A. Jhalani, C. E. Dreyer, and M. Bernardi, Phys. Rev. B 102, 125203 (2020).
  * Poncé et al. (2021) S. Poncé, F. Macheda, E. R. Margine, N. Marzari, N. Bonini, and F. Giustino, Phys. Rev. Research 3, 043022 (2021).
  * Kaasbjerg et al. (2012) K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Phys. Rev. B 85, 115317 (2012).
  * Li et al. (2013) X. Li, J. T. Mullen, Z. Jin, K. M. Borysenko, M. Buongiorno Nardelli, and K. W. Kim, Phys. Rev. B 87, 115418 (2013).
  * Kaasbjerg et al. (2013) K. Kaasbjerg, K. S. Thygesen, and A.-P. Jauho, Phys. Rev. B 87, 235312 (2013).
  * Zhang et al. (2014) W. Zhang, Z. Huang, W. Zhang, and Y. Li, Nano Research 7, 1731 (2014).
  * Park et al. (2014) C.-H. Park, N. Bonini, T. Sohier, G. Samsonidze, B. Kozinsky, M. Calandra, F. Mauri, and N. Marzari, Nano letters 14, 1113 (2014).
  * Sohier et al. (2014) T. Sohier, M. Calandra, C.-H. Park, N. Bonini, N. Marzari, and F. Mauri, Phys. Rev. B 90, 125414 (2014).
  * Gunst et al. (2016) T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge, Phys. Rev. B 93, 035414 (2016).
  * Sohier et al. (2018) T. Sohier, D. Campi, N. Marzari, and M. Gibertini, Phys. Rev. Materials 2, 114010 (2018).
  * Gaddemane et al. (2021) G. Gaddemane, S. Gopalan, M. V. de Put, and M. V. Fischetti, J. Comput. Elect. 20, 49 (2021).
  * Guo et al. (2019) F. Guo, Z. Liu, M. Zhu, and Y. Zheng, Phys. Chem. Chem. Phys. 21, 22879 (2019).
  * Sio and Giustino (2022) W. H. Sio and F. Giustino, Phys. Rev. B 105, 115414 (2022).
  * Efetov and Kim (2010) D. K. Efetov and P. Kim, Physical Review Letters 105, 256805 (2010).
  * Poncé et al. (2020) S. Poncé, W. Li, S. Reichardt, and F. Giustino, Reports on Progress in Physics 83, 036501 (2020).
  * Lundstrom (2009) M. S. Lundstrom, _Fundamentals of carrier transport_ (Cambridge University Press, 2009).
  * Ganose et al. (2021) A. M. Ganose, J. Park, A. Faghaninia, R. Woods-Robinson, K. A. Persson, and A. Jain, Nature Communications 12 (2021).
  * Claes et al. (2022) R. Claes, G. Brunin, M. Giantomassi, G.-M. Rignanese, and G. Hautier, Phys. Rev. B 106, 094302 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.106.094302.
  * Giustino et al. (2007) F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B 76, 165108 (2007).
  * Calandra et al. (2010) M. Calandra, G. Profeta, and F. Mauri, Phys. Rev. B 82, 165111 (2010).
  * Poncé et al. (2016) S. Poncé, E. Margine, C. Verdi, and F. Giustino, Computer Physics Communications 209, 116 (2016).
  * Gonze et al. (2020) X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, et al., Comput. Phys. Commun. 248, 107042 (2020).
  * Romero et al. (2020) A. H. Romero, D. C. Allan, B. Amadon, G. Antonius, T. Applencourt, L. Baguet, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, et al., The Journal of Chemical Physics 152, 124102 (2020).
  * Zhou et al. (2021) J.-J. Zhou, J. Park, I.-T. Lu, I. Maliyov, X. Tong, and M. Bernardi, Computer Physics Communications 264, 107970 (2021).
  * Marzari et al. (2012) N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Rev. Mod. Phys. 84, 1419 (2012).
  * Vogl (1976) P. Vogl, Phys. Rev. B 13, 694 (1976).
  * Sjakste et al. (2015) J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Phys. Rev. B 92, 054307 (2015).
  * Verdi and Giustino (2015) C. Verdi and F. Giustino, Phys. Rev. Lett. 115, 176401 (2015).
  * Sohier et al. (2016) T. Sohier, M. Calandra, and F. Mauri, Phys. Rev. B 94, 085415 (2016).
  * Sohier et al. (2017b) T. Sohier, M. Gibertini, M. Calandra, F. Mauri, and N. Marzari, Nano Letters 17, 3758 (2017b).
  * Deng et al. (2021) T. Deng, G. Wu, W. Shi, Z. M. Wong, J.-S. Wang, and S.-W. Yang, Phys. Rev. B 103, 075410 (2021).
  * Royo and Stengel (2021) M. Royo and M. Stengel, Phys. Rev. X 11, 041027 (2021).
  * Zein (1984) N. Zein, Fiz. Tverd. Tela 26, 3028 (1984).
  * Baroni et al. (1987) S. Baroni, P. Giannozzi, and A. Testa, Phys. Rev. Lett. 58, 1861 (1987).
  * Gonze (1997) X. Gonze, Phys. Rev. B 55, 10337 (1997).
  * Martin et al. (2016) R. M. Martin, L. Reining, and D. M. Ceperley, _Interacting Electrons_ (Cambridge University Press, 2016).
  * Royo et al. (2020) M. Royo, K. R. Hahn, and M. Stengel, Phys. Rev. Lett. 125, 217602 (2020).
  * Gonze et al. (1994) X. Gonze, J.-C. Charlier, D. Allan, and M. Teter, Phys. Rev. B 50, 13035 (1994).
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Tsirkin (2021) S. S. Tsirkin, npj Computational Materials 7 (2021).
  * Stengel and Spaldin (2006) M. Stengel and N. A. Spaldin, Phys. Rev. B 73, 075121 (2006).
  * Reggiani et al. (1983) L. Reggiani, D. Waechter, and S. Zukotynski, Phys. Rev. B 28, 3550 (1983).
  * Popovic (1991) R. Popovic, _Hall Effect Devices: Magnetic Sensors and Characterization of Semiconductors_ (Taylor & Francis: Boca Raton, 1991).
  * Giannozzi et al. (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi et al. (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, et al., Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Gonze et al. (2016) X. Gonze, F. Jollet, F. A. Araujo, D. Adams, B. Amadon, T. Applencourt, C. Audouze, J.-M. Beuken, J. Bieder, A. Bokhanchuk, et al., Computer Physics Communications 205, 106 (2016).
  * Royo and Stengel (2019) M. Royo and M. Stengel, Phys. Rev. X 9, 021050 (2019).
  * Pizzi et al. (2020) G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, et al., Journal of Physics: Condensed Matter 32, 165902 (2020).
  * Hamann (2013) D. R. Hamann, Phys. Rev. B 88, 085117 (2013).
  * van Setten et al. (2018) M. van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze, and G.-M. Rignanese, Computer Physics Communications 226, 39 (2018).
  * Perdew and Wang (1992) J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992).
  * Liu et al. (2019) D. Liu, X. Chen, Y. Yan, Z. Zhang, Z. Jin, K. Yi, C. Zhang, Y. Zheng, Y. Wang, J. Yang, et al., Nature Commun. 10, 1188 (2019).
  * Zhu et al. (2015) H. Zhu, Y. Wang, J. Xiao, M. Liu, S. Xiong, Z. J. Wong, Z. Ye, Y. Ye, X. Yin, and X. Zhang, Nature Nanotechnology 10, 151 (2015).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Molina-Sánchez et al. (2015) A. Molina-Sánchez, K. Hummer, and L. Wirtz, Surface Science Reports 70, 554 (2015).
  * Li and Li (2015) W. Li and J. Li, Nano Res. 8, 3796 (2015).
  * Sucharitakul et al. (2015) S. Sucharitakul, N. J. Goble, U. R. Kumar, R. Sankar, Z. A. Bogorad, F.-C. Chou, Y.-T. Chen, and X. P. A. Gao, Nano Letters 15, 3815 (2015).
  * Bandurin et al. (2017) D. A. Bandurin, A. V. Tyurnina, G. L. Yu, A. Mishchenko, V. Zólyomi, S. V. Morozov, R. K. Kumar, R. V. Gorbachev, Z. R. Kudrynskyi, S. Pezzini, et al., Nature Nanotechnology 12, 223 (2017).
  * Ho et al. (2017) P.-H. Ho, Y.-R. Chang, Y.-C. Chu, M.-K. Li, C.-A. Tsai, W.-H. Wang, C.-H. Ho, C.-W. Chen, and P.-W. Chiu, ACS Nano 11, 7362 (2017).
  * Li et al. (2018) M. Li, C.-Y. Lin, S.-H. Yang, Y.-M. Chang, J.-K. Chang, F.-S. Yang, C. Zhong, W.-B. Jian, C.-H. Lien, C.-H. Ho, et al., Advanced Materials 30, 1803690 (2018).
  * Sun et al. (2016) C. Sun, H. Xiang, B. Xu, Y. Xia, J. Yin, and Z. Liu, Applied Physics Express 9, 035203 (2016).
  * Hu et al. (2017) T. Hu, J. Zhou, and J. Dong, Phys. Chem. Chem. Phys. 19, 21722 (2017).
  * Peng et al. (2017) Q. Peng, R. Xiong, B. Sa, J. Zhou, C. Wen, B. Wu, M. Anpo, and Z. Sun, Catal. Sci. Technol. 7, 2744 (2017).
  * Wang et al. (2019) Q. Wang, L. Han, W. L., T. Zhang, S. Li, and P. Lu, Nanoscale Research Letters p. 287 (2019).
  * Gopalan et al. (2019) S. Gopalan, G. Gaddemane, M. L. Van de Put, and M. V. Fischetti, Materials p. 4210 (2019).
  * Geim and Novoselov (2007) A. K. Geim and K. S. Novoselov, Nature Materials 6, 183 (2007).
  * Peres (2010) N. M. R. Peres, Rev. Mod. Phys. 82, 2673 (2010).
  * Marzari et al. (1999) N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne, Phys. Rev. Lett. 82, 3296 (1999).
  * Yang et al. (2018) G. Yang, L. Li, W. B. Lee, and M. C. Ng, Science and Technology of Advanced Materials 19, 613 (2018).
  * Xia et al. (2014) F. Xia, H. Wang, and Y. Jia, Nature Communications 5, 4458 (2014).
  * Liu et al. (2014) H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, ACS Nano 8, 4033 (2014).
  * Li et al. (2014) L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nature Nanotechnology 9, 372 (2014).
  * Long et al. (2016a) G. Long, D. Maryenko, J. Shen, S. Xu, J. Hou, Z. Wu, W. K. Wong, T. Han, J. Lin, Y. Cai, et al., Nano Letters 16, 7768 (2016a).
  * Long et al. (2016b) G. Long, S. Xu, J. Shen, J. Hou, Z. Wu, T. Han, J. Lin, W. K. Wong, Y. Cai, R. Lortz, et al., 2D Materials 3, 031001 (2016b).
  * Zollner et al. (1990) S. Zollner, S. Gopalan, and M. Cardona, Journal of Applied Physics 68, 1682 (1990).
  * Shafique et al. (2017) A. Shafique, A. Samad, and Y.-H. Shin, Phys. Chem. Chem. Phys. 19, 20677 (2017).
  * Zschieschang et al. (2014) U. Zschieschang, T. Holzmann, B. V. Lotsch, and H. Klauk, in _72nd Device Research Conference_ (2014), pp. 267–268.
  * Shibata et al. (1991) T. Shibata, Y. Muranushi, T. Miura, and T. Kishi, Journal of Physics and Chemistry of Solids 52, 551 (1991).
  * Song et al. (2013) H. S. Song, S. L. Li, L. Gao, Y. Xu, K. Ueno, J. Tang, Y. B. Cheng, and K. Tsukagoshi, Nanoscale 5, 9666 (2013).
  * Huang et al. (2014) Y. Huang, E. Sutter, J. T. Sadowski, M. Cotlet, O. L. A. Monti, D. A. Racke, M. R. Neupane, D. Wickramaratne, R. K. Lake, B. A. Parkinson, et al., ACS Nano 8, 10743 (2014).
  * Giri et al. (2019) B. Giri, M. Masroor, T. Yan, K. Kushnir, A. D. Carl, C. Doiron, H. Zhang, Y. Zhao, A. McClelland, G. A. Tompsett, et al., Advanced Energy Materials 9, 1901236 (2019).
  * Yu et al. (2014) Z. Yu, Y. Pan, Y. Shen, Z. Wang, Z. Y. Ong, T. Xu, R. Xin, L. Pan, B. Wang, L. Sun, et al., Nature Communications 5, 5290 (2014).
  * Liu et al. (2015) Y. Liu, H. Wu, H.-C. Cheng, S. Yang, E. Zhu, Q. He, M. Ding, J. Li, D. Guo, W. N. O., Y. Huang, et al., Nano Lett. 15, 3030 (2015).
  * Yu et al. (2016) Z. Yu, Z.-Y. Ong, Y. Pan, Y. Cui, R. Xin, Y. Shi, B. Wang, Y. Wu, T. Chen, Y.-W. Zhang, et al., Advanced Materials 28, 547 (2016).
  * Sanne et al. (2015) A. Sanne, R. Ghosh, A. Rai, H. C. P. Movva, A. Sharma, R. Rao, L. Mathew, and S. K. Banerjee, Applied Physics Letters 106, 062101 (2015).
  * Kang et al. (2015) K. Kang, S. Xie, L. Huang, Y. Han, P. Y. Huang, K. F. Mak, C. J. Kim, D. Muller, and J. Park, Nature 520, 656 (2015).
  * Cui et al. (2015) X. Cui, G.-H. Lee, Y. D. Kim, G. Arefe, P. Y. Huang, C.-H. Lee, D. A. Chenet, X. Zhang, L. Wang, F. Ye, et al., Nature Nanotechnology 10, 534 (2015).
  * Huo et al. (2018) N. Huo, Y. Yang, Y.-N. Wu, X.-G. Zhang, S. T. Pantelides, and G. Konstantatos, Nanoscale 10, 15071 (2018).
  * Momose et al. (2018) T. Momose, A. Nakamura, M. Daniel, and M. Shimomura, AIP Advances 8, 025009 (2018).
  * Sohier et al. (2019) T. Sohier, E. Ponomarev, M. Gibertini, H. Berger, N. Marzari, N. Ubrig, and A. F. Morpurgo, Phys. Rev. X 9, 031019 (2019).
  * Doan et al. (2014) T. C. Doan, J. Li, J. Y. Lin, and H. X. Jiang, AIP Advances 4, 107126 (2014).
  * Grenadier et al. (2021) S. Grenadier, A. Maity, J. Li, J. Lin, and H. Jiang, Semiconductors and Semimetals 107, 393 (2021).
  * Grenadier et al. (2019) S. Grenadier, A. Maity, J. Li, J. Y. Lin, and H. X. Jiang, Applied Physics Letters 115, 072108 (2019).
  * Martin (1972) R. M. Martin, Phys. Rev. B 5, 1607 (1972).
  * Feng et al. (2014) W. Feng, W. Zheng, W. Cao, and P. Hu, Advanced Materials 26, 6587 (2014).
  * Yang et al. (2017) Z. Yang, W. Jie, C.-H. Mak, S. Lin, H. Lin, X. Yang, F. Yan, S. P. Lau, and J. Hao, ACS Nano 11, 4225 (2017).
  * Chang et al. (2018) H.-C. Chang, C.-L. Tu, K.-I. Lin, J. Pu, T. Takenobu, C.-N. Hsiao, and C.-H. Chen, Small 14, 1802351 (2018).
  * Jiang et al. (2019) J. Jiang, J. Li, Y. Li, J. Duan, L. Li, Y. Tian, Z. Zong, H. Zheng, X. Feng, Q. Li, et al., npj 2D Materials and Applications 3, 29 (2019).
  * Li and Giustino (2020) W. Li and F. Giustino, Phys. Rev. B 101, 035201 (2020).
  * Shi et al. (2019) L.-B. Shi, S. Cao, M. Yang, Q. You, K.-C. Zhang, Y. Bao, Y.-J. Zhang, Y.-Y. Niu, and P. Qian, Journal of Physics: Condensed Matter 32, 065306 (2019).
  * Fischetti and Vandenberghe (2016) M. V. Fischetti and W. G. Vandenberghe, Phys. Rev. B 93, 155413 (2016).
  * Rybkovskiy et al. (2014) D. V. Rybkovskiy, A. V. Osadchy, and E. D. Obraztsova, Phys. Rev. B 90, 235302 (2014).
  * Lugovskoi et al. (2019) A. V. Lugovskoi, M. I. Katsnelson, and A. N. Rudenko, Phys. Rev. Lett. 123, 176401 (2019).
  * Kibirev et al. (2018) I. A. Kibirev, A. V. Matetskiy, A. V. Zotov, and A. A. Saranin, Applied Physics Letters 112, 191602 (2018).
  * Zólyomi et al. (2014) V. Zólyomi, N. D. Drummond, and V. I. Fal’ko, Phys. Rev. B 89, 205416 (2014).
  * Du et al. (2008) X. Du, I. Skachko, A. Barker, and E. Y. Andrei, Nature Nanotechnology 3, 491 (2008).
  * Bolotin et al. (2008) K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. Stormer, Solid State Communications 146, 351 (2008).
  * Castro et al. (2010) E. V. Castro, H. Ochoa, M. I. Katsnelson, R. V. Gorbachev, D. C. Elias, K. S. Novoselov, A. K. Geim, and F. Guinea, Phys. Rev. Lett. 105, 266601 (2010).
  * Liao et al. (2015) B. Liao, J. Zhou, B. Qiu, M. S. Dresselhaus, and G. Chen, Phys. Rev. B 91, 235419 (2015).
  * Rudenko et al. (2016) A. N. Rudenko, S. Brener, and M. I. Katsnelson, Phys. Rev. Lett. 116, 246401 (2016).
  * Jin et al. (2016) Z. Jin, J. T. Mullen, and K. W. Kim, Applied Physics Letters 109, 053108 (2016).
  * Trushkov and Perebeinos (2017) Y. Trushkov and V. Perebeinos, Phys. Rev. B 95, 075436 (2017).
  * Gaddemane et al. (2018) G. Gaddemane, W. G. Vandenberghe, M. L. Van de Put, S. Chen, S. Tiwari, E. Chen, and M. V. Fischetti, Phys. Rev. B 98, 115416 (2018).