toltxlabel

# Efficient and flexible estimation of natural mediation effects under
intermediate confounding and monotonicity constraints

Kara E. Rudolph corresponding author: Department of Epidemiology, Mailman
School of Public Health, Columbia University. Iván Díaz  Division of
Biostatistics, Weill Cornell Medicine.

###### Abstract

Natural direct and indirect effects are mediational estimands that decompose
the average treatment effect and describe how outcomes would be affected by
contrasting levels of a treatment through changes induced in mediator values
(in the case of the indirect effect) or not through induced changes in the
mediator values (in the case of the direct effect). Natural direct and
indirect effects are not generally point-identifiable in the presence of a
treatment-induced confounder, however they may still be identified if one is
willing to assume monotonicity between a treatment and the treatment-induced
confounder. We argue that this assumption may be reasonable in the relatively
common encouragement-design trial setting where intervention is randomized
treatment assignment and the treatment-induced confounder is whether or not
treatment was actually taken/adhered to. We develop efficiency theory for the
natural direct and indirect effects under this monotonicity assumption, and
use it to propose a nonparametric, multiply robust estimator. We demonstrate
the finite sample properties of this estimator using a simulation study, and
apply it to data from the Moving to Opportunity Study to estimate the natural
direct and indirect effects of being randomly assigned to receive a Section 8
housing voucher—the most common form of federal housing assistance—on risk
developing any mood or externalizing disorder among adolescent boys—-possibly
operating through various school and community characteristics.

##  1 Introduction

Researchers are frequently interested in summarizing mediational effects
across individuals (or other units). For example, in trying to understand
reasons underlying an unexpected treatment effect, it may be helpful to
decompose the overall effect into the portion operating through some
intermediate variables (i.e., mediators)—called the indirect
effect—hypothesized to be responsible for the unexpected overall effect, and
the portion not operating through those intermediate variables—called the
direct effect. Natural direct and indirect effects (NDE, NIE, which we also
refer to as natural (in)direct effects) are a type of mediation estimand that
allows for such a decomposition. They are descriptive in nature (Pearl, 2013);
the NIE describes how individuals’ outcomes would be affected by changes in
their mediator values induced by contrasting levels of a treatment/exposure.
(We formally define the NIE and NDE in terms of counterfactual notation in
Section 2.)

Generally, the NDE/NIE are not point-identifiable in the presence of
treatment-induced confounder of the mediator-outcome relationship (Avin et
al., 2005). To give intuition for this, we first define a counterfactual
outcome as an outcome variable had treatment been set to some value $a$ and
the mediator been set to some value, $m$, possibly contrary to fact, denoted
$Y_{A=a,M=m}$ (we provide a more detailed description of notation in Section
2). Similarly, we define a counterfactual mediator as a mediator variable had
treatment been set to some value $a$, possibly contrary to fact, denoted
$M_{A=a}$. NDE/NIE are then generally defined in terms of nested
counterfactual outcomes, such as $Y_{A=a,M_{A=a^{\prime}}}$ (Pearl, 2001).
Identification of the counterfactual outcome $Y_{A=a,M_{A=a^{\prime}}}$ from
observed data becomes a problem, because setting treatment to a certain value
($A=a$) induces a counterfactual treatment-induced confounder under the same
treatment, denoted $Z_{A=a}$. Setting treatment to the other value in the
counterfactual mediator ($A=a^{\prime}$) induces a counterfactual treatment-
induced confounder under the other value of treatment, denoted
$Z_{A=a^{\prime}}$. The two treatment-induced counterfactual confounders,
$Z_{A=a}$ and $Z_{A=a^{\prime}}$ are correlated with each other, because they
share unmeasured common causes, $U_{Z}$. $Z_{A=a}$ is also correlated with
$Y_{A=a,M=m}$, and $Z_{A=a^{\prime}}$ is also correlated with
$M_{A=a^{\prime}}$, which precludes identification.

This general lack of identification in the presence of a treatment-induced
confounder presents a problem for applied mediational research, because such
variables are near-ubiquitous. For example, in trials where an individual or
community is randomized to a particular treatment but cannot be forced to take
the treatment, treatment take-up/ adherence represents a treatment-induced
confounder. In observational data, there may be multiple intermediate
variables linking treatment/exposure to the outcome, but only a subset are of
interest to examine as mediators (for example, in the case discussed above
where one is interested in examining which mediators contribute to the
unexpected overall effect). The remaining intermediate variables would be
considered as treatment-induced confouders.

However, particular cases exist when the NDE and NIE are identifiable even in
the presence of one or more treatment-induced confounders (Miles, 2022;
Tchetgen and VanderWeele, 2014). One of these cases is when one can assume
monotonicity between a treatment and one or more binary treatment-induced
confounders (Tchetgen and VanderWeele, 2014). Assuming monotonicity is also
sometimes referred to as assuming “no defiers” (Angrist et al., 1996). In the
context of an encouragement design intervention, this means that: if treatment
assignment acts to increase the likelihoood treatment take-up, that there is
no one for whom receiving treatment would result in them not taking the
treatment, but not receiving the treatment would result in them taking the
treatment.

Monotonicity may seem like a restrictive assumption, but, in fact, it may be
reasonable in some common scenarios. For example, consider the relatively
common trial setting where the intervention is randomized treatment assignment
and the treatment-induced confounder is whether or not treatment was actually
taken. In this setting, we may feel comfortable assuming that treatment
assignment would not make anyone more likely not to take the treatment.
Indeed, monotonicity in this context is frequently assumed when using
instrumental variables (IV) to identify causal effects (Angrist et al., 1996).
Using IV as an identification strategy is common in the context of randomized
trials where randomized assignment is the IV and treatment uptake or adherence
is the exposure of interest. We consider this same context as our motivating
example, but consider so-called intent-to-treat effects, where randomized
assignment is our treatment of interest (not an IV), and treatment
uptake/adherence is the treatment-induced confounder.

Although monotonicity may be a reasonable assumption in trial settings, we
know of no nonparametric estimators of natural (in)direct effects in this
context (i.e., in the presence of treatment-induced confounders, assuming
monotonicity). Consequently, we develop efficiency theory for the NDE/NIE
under monotonicity, and use it to propose a nonparametric, multiply robust
estimator of NDE/NIE based on solving the efficient influence function (EIF)
estimating equation. Of relevance for real-world analyses, our estimator
allows for continuous and/or multiple mediators and is cross-fitted. This
estimator should be applicable whenever natural direct or indirect effects of
treatment assignment are of interest in the presence potentially imperfect
compliance acting as a single treatment-induced confounder.

This paper is organized as follows. We give notation, the structural causal
model, and review the definition and identification of NDE/NIE in Section 2.
We propose a nonparametric, robust estimator using flexible, data-adaptive
regression methods in Section 3. Section 4 details a limited simulation study
evaluating the finite sample performance of the proposed estimator. In Section
5, we apply the estimator to our motivating example estimating the NDE and NIE
of being randomly assigned to receive a Section 8 housing voucher—the most
common form of federal housing assistance—on risk developing any mood or
externalizing disorder among adolescent boys—-possibly operating through
various school and community characteristics, using longitudinal data from the
Moving to Opportunity Study (MTO) (Sanbonmatsu et al., 2011). Section 6
concludes.

##  2 Notation, structural causal model, estimand definition

Let $O=(W,A,Z,M,Y)$ denote the observed data, and let $O_{1},\ldots,O_{n}$
denote a sample of $n$ i.i.d. observations of $O$. Let $W$ denote a vector of
observed baseline covariates, $W=f(U_{W})$, where $U_{W}$ is unobserved
exogenous error on $W$ and where function $f$ is assumed deterministic but
unknown (Pearl, 2009). Let $A$ denote a binary treatment/exposure variable,
$A=f(W,U_{A})$. Let $Z$ denote single, binary treatment-induced confounder of
the mediator-outcome relationship, e.g., initiation of treatment/ adherence to
treatment assignment, $Z=f(W,A,U_{Z})$, though we note that our results could
be extended to accomodate multiple, binary $Z$. Let $M$ denote a set of
mediating variables, $M=f(W,A,Z,U_{M})$, that may be multiple, multi-valued,
and/or continuous. Finally, let $Y$ denote a continuous or binary outcome,
$Y=f(W,A,Z,M,U_{Y}).$

We use $\mathsf{P}$ to denote the distribution of $O$, and $\mathsf{P}^{c}$ to
denote the distribution of $(O,U)$. We let $\mathsf{P}$ be an element of the
nonparametric statistical model defined as all continuous densities on $O$
with respect to some dominating measure $\nu$. Let $\mathsf{p}$ denote the
corresponding probability density function, and $\mathcal{W,A,Z,M}$ denote the
range of the respective random variables. We let $\mathsf{E}$ and
$\mathsf{E}^{c}$ denote corresponding expectation operators, and define
$\mathsf{P}f=\int f(o)\mathrm{d}\mathsf{P}(o)$ for a given function $f(o)$. We
use $\mathsf{g}(a\mid w)$ to denote the probability mass function of $A=a$
conditional on $W=w$, and $\mathsf{e}(a\mid m,z,w)$ to denote the probability
mass function of $A=a$ conditional on $(M,Z,W)=(m,z,w)$. We use $\mu(m,z,a,w)$
to denote the outcome regression function $\mathsf{E}(Y\mid M=m,Z=z,A=a,W=w)$.
We use $\mathsf{q}(z\mid a,w)$ and $\mathsf{r}(z\mid a,m,w)$ to denote the
corresponding conditional densities of $Z$.

We define counterfactual variables in terms of interventions on the
nonparametric structural equation model (NPSEM). For simplicity of notation,
we will use the random variable with its corresponding intervention in the
index to denote the counterfactual. For example, $Y_{A=1}$ denotes the random
variable $f_{Y}(W,1,Z,M,U_{Y})$. In what follows, we drop the random variable
from the index of the counterfactual whenever the variable that is being
intervened on is clear from context. For example, we simply use $Y_{a,m}$ to
denote $Y_{A=a,M=m}$, and $M_{a}$ to denote $M_{A=a}$.

We are interested in the NDE and NIE, which decompose the overall average
treatment effect (ATE) of the binary treatment/exposure as follows:

$\mathsf{E}^{c}(Y_{1,M_{1}}-Y_{0,M_{0}})=\underbrace{\mathsf{E}^{c}(Y_{1,M_{1}}-Y_{1,M_{0}})}_{\text{natural
indirect effect (through
$M$)}}+\underbrace{\mathsf{E}^{c}(Y_{1,M_{0}}-Y_{0,M_{0}})}_{\text{natural
direct effect (not through $M$)}}.$

Tchetgen and VanderWeele (2014) identified the NDE and NIE from observed data
$O$ in the presence of treatment-induced confounders, $Z$. For a single,
binary $Z$, the parameter $\mathsf{E}^{c}[Y_{a,M_{a^{\prime}}}]$, is
identified as

$\theta(a,a^{\prime})=\theta_{1,1}(a,a^{\prime})+\theta_{1,0}(a,a^{\prime})+\theta_{0,0}(a,a^{\prime}),$

where

$\displaystyle\theta_{1,1}(a,a^{\prime})$ $\displaystyle=\int\mathsf{E}(Y\mid
a,m,1,w)\mathrm{d}\mathsf{P}(m\mid a^{\prime},1,w)\mathsf{P}(Z=1\mid
a^{\prime},w)\mathrm{d}\mathsf{P}(w),$
$\displaystyle\theta_{1,0}(a,a^{\prime})$ $\displaystyle=\int\mathsf{E}(Y\mid
a,m,1,w)\mathrm{d}\mathsf{P}(m\mid a^{\prime},0,w)[\mathsf{P}(Z=1\mid
a,w)-\mathsf{P}(Z=1\mid a^{\prime},w)]\mathrm{d}\mathsf{P}(w),$
$\displaystyle\theta_{0,0}(a,a^{\prime})$ $\displaystyle=\int\mathsf{E}(Y\mid
a,m,0,w)\mathrm{d}\mathsf{P}(m\mid a^{\prime},0,w)\mathsf{P}(Z=0\mid
a,w)\mathrm{d}\mathsf{P}(w),$

under the following assumptions:

######  A1Monotonicity.

If $a^{\prime}<a,$ then $Z_{a^{\prime}}\leq Z_{a}$ for all $i$;

######  A2Sequential Exchangeability.

$Y_{A=a,M=m}\mbox{$\perp\\!\\!\\!\perp$}A\mid W$,
$Y_{A=a,M=m}\mbox{$\perp\\!\\!\\!\perp$}M\mid W,A,Z$, and
$M_{A=a}\mbox{$\perp\\!\\!\\!\perp$}A\mid W$;

######  A3Positivity of treatment/exposure, treatment-induced confounder, and
mediator mechanisms.

Assume:

  * •

$\mathsf{p}(w)>0$ implies $\mathsf{p}(a\mid w)>0$ for $a\in\\{0,1\\}$;

  * •

$[\mathsf{P}(Z=1\mid a^{\prime},w)\times
I(z=z^{\prime}=1)+\\{\mathsf{P}(Z=1\mid a,w)-\mathsf{P}(Z=0\mid
a^{\prime},w)\\}\times I(z=1,z^{\prime}=0)+\mathsf{P}(Z=0\mid a,w)\times
I(z=z^{\prime}=0)]>0$ implies $\mathsf{p}(m\mid z^{\prime},a^{\prime},w)>0$
and $\mathsf{p}(z\mid a,w)>0$ for $(z,z^{\prime})\in\\{(0,0),(1,0),(1,1)\\}$;
and

  * •

$\mathsf{p}(m\mid z^{\prime},a^{\prime},w)>0$ implies $\mathsf{p}(m\mid
z,a,w)>0$ for all observed $m,z,z^{\prime}$ and for
$(a,a^{\prime})\in\\{(0,0),(1,0),(1,1)\\}$.

##  3 Estimation

We propose a one-step estimator for the statistical parameter
$\theta(a,a^{\prime})$ that is the sample average of the estimated uncentered
efficient influence function (EIF). An R package to implement this estimator
is included linkblindedforreview. The EIF for this parameter is given by the
following expressions. Define

$\displaystyle H_{Y,1,1}$
$\displaystyle=\frac{\mathds{1}\\{Z=1,A=a\\}}{\mathsf{P}(a^{\prime}\mid
W)}\frac{\mathsf{P}(a^{\prime}\mid M,1,W)}{\mathsf{P}(a\mid M,1,W)}$
$\displaystyle H_{Y,1,0}$
$\displaystyle=\frac{\mathds{1}\\{Z=1,A=a\\}}{\mathsf{P}(a^{\prime}\mid
W)\mathsf{P}(Z=0\mid a^{\prime},W)}\frac{\mathsf{P}(a^{\prime}\mid
M,1,W)}{\mathsf{P}(a\mid M,1,W)}\frac{\mathsf{P}(Z=0\mid
M,a^{\prime},W)}{\mathsf{P}(Z=1\mid M,a^{\prime},W)}\\{\mathsf{P}(Z=1\mid
a,w)-\mathsf{P}(Z=1\mid a^{\prime},w)\\}$ $\displaystyle H_{Y,0,0}$
$\displaystyle=\frac{\mathds{1}\\{Z=0,A=a\\}}{\mathsf{P}(a\mid
W)}\frac{\mathsf{P}(a^{\prime}\mid M,0,W)}{\mathsf{P}(a\mid
M,0,W)}\frac{\mathsf{P}(Z=0\mid a,W)}{\mathsf{P}(Z=0\mid a^{\prime},W)}$
$\displaystyle H_{M,1,1}$
$\displaystyle=\frac{\mathds{1}\\{A=a^{\prime},Z=1\\}}{\mathsf{P}(a^{\prime}\mid
W)}$ $\displaystyle H_{M,1,0}$
$\displaystyle=\frac{\mathds{1}\\{A=a^{\prime},Z=0\\}}{\mathsf{P}(Z=0\mid
a^{\prime},W)\mathsf{P}(a^{\prime}\mid W)}\\{\mathsf{P}(Z=1\mid
a,w)-\mathsf{P}(Z=1\mid a^{\prime},w)\\}$ $\displaystyle H_{M,0,0}$
$\displaystyle=\frac{\mathds{1}\\{A=a^{\prime},Z=0\\}}{\mathsf{P}(Z=0\mid
a^{\prime},W)\mathsf{P}(a^{\prime}\mid W)}\mathsf{P}(Z=0\mid a,w)$
$\displaystyle H_{Z,1,1}$
$\displaystyle=\frac{\mathds{1}\\{A=a^{\prime}\\}}{\mathsf{P}(a^{\prime}\mid
W)}\int\mathsf{E}(Y\mid a,m,1,W)\mathrm{d}\mathsf{P}(m\mid a^{\prime},1,W)$
$\displaystyle H_{Z,1,0}$
$\displaystyle=\left(\frac{\mathds{1}\\{A=a\\}}{\mathsf{P}(a\mid
W)}-\frac{\mathds{1}\\{A=a^{\prime}\\}}{\mathsf{P}(a^{\prime}\mid
W)}\right)\int\mathsf{E}(Y\mid a,m,1,W)\mathrm{d}\mathsf{P}(m\mid
a^{\prime},0,W)$ $\displaystyle H_{Z,0,0}$
$\displaystyle=-\frac{\mathds{1}\\{A=a\\}}{\mathsf{P}(a\mid
W)}\int\mathsf{E}(Y\mid a,m,0,W)\mathrm{d}\mathsf{P}(m\mid a^{\prime},0,W)$
$\displaystyle H_{W,1,1}$ $\displaystyle=\int\mathsf{E}(Y\mid
a,m,1,w)\mathrm{d}\mathsf{P}(m\mid a^{\prime},1,w)\mathsf{P}(Z=1\mid
a^{\prime},w)$ $\displaystyle H_{W,1,0}$ $\displaystyle=\int\mathsf{E}(Y\mid
a,m,1,W)\mathrm{d}\mathsf{P}(m\mid a^{\prime},0,W)[\mathsf{P}(Z=1\mid
a,W)-\mathsf{P}(Z=1\mid a^{\prime},W)]$ $\displaystyle H_{W,0,0}$
$\displaystyle=\int\mathsf{E}(Y\mid a,m,0,W)\mathrm{d}\mathsf{P}(m\mid
a^{\prime},0,W)\mathsf{P}(Z=0\mid a,W)$

Then, for $(z,z^{\prime})\in\\{(1,1),(1,0),(0,0)\\}$, the EIF of
$\theta_{z,z^{\prime}}(a,a^{\prime})$ is equal to
$\bar{D}_{z,z^{\prime}}(O;\eta)=D_{z,z^{\prime}}(O;\eta)-\theta_{z,z^{\prime}}(a,a^{\prime})$,
where

$\displaystyle D_{z,z^{\prime}}(O)$
$\displaystyle=H_{Y,z,z^{\prime}}\left[Y-\mathsf{E}(Y\mid A,M,Z,W)\right]$
$\displaystyle+H_{Z,z,z^{\prime}}[Z-\mathsf{P}(Z=1\mid A,W)]$
$\displaystyle+H_{M,z,z^{\prime}}\left[\mathsf{E}(Y\mid
a,M,z,W)-\int\mathsf{E}(Y,a,m,z,W)\mathrm{d}\mathsf{P}(m\mid
a^{\prime},z^{\prime},W)\right]$ $\displaystyle+H_{W,z,z^{\prime}}.$

The EIF of $\theta(a,a^{\prime})$ is then
$\bar{D}(O)=D(O)-\theta(a,a^{\prime})$, where
$D(O)=D_{1,1}(O)+D_{1,0}(O)+D_{0,0}(O).$

Let $\rho$ denote $\int\mathsf{E}(Y\mid a,m,z,w)\mathrm{d}\mathsf{P}(m\mid
a^{\prime},z,w)$. Let
$\eta=(\mathsf{g},\mathsf{q},\mathsf{e},\mathsf{r},\mu,\rho)$, where each is
defined as in Section 2. Let $\hat{\eta}$ denote an estimator of $\eta$. Let
$D(O,\mathsf{P})=D(O,\eta),$ as $\eta$ contains all the relevant features of
$\mathsf{P}$. Let $D(O,\hat{\eta})$ denote an estimator of $D(O,\eta)$. The
one-step estimator we propose is the sample average of $D(O,\hat{\eta}).$

We use a cross-fitted version of this estimator. Cross-fitting is a data-
splitting technique that weakens some of the technical assumptions (i.e.,
Donsker-type assumptions) required for asymptotic normality (Klaassen, 1987;
Zheng and van der Laan, 2011; Chernozhukov et al., 2016). We perform
crossfitting for estimation of all the components of $\eta$ as follows. Let
${\cal V}_{1},\ldots,{\cal V}_{J}$ denote a random partition of data with
indices $i\in\\{1,\ldots,n\\}$ into $J$ prediction sets of approximately the
same size such that $\bigcup_{j=1}^{J}{\cal V}_{j}=\\{1,\ldots,n\\}$. For each
$j$, the training sample is given by ${\cal
T}_{j}=\\{1,\ldots,n\\}\setminus{\cal V}_{j}$. $\hat{\eta}_{j}$ denotes the
estimator of $\eta$, obtained by training the corresponding prediction
algorithm using only data in the sample ${\cal T}_{j}$, and $j(i)$ denotes the
index of the validation set which contains observation $i$. We then use these
fits, $\hat{\eta}_{j(i)}(O_{i})$ in computing each efficient influence
function, i.e., we compute $D(O_{i},\hat{\eta}_{j(i)})$.

Thus, our one-step estimator of $\theta(a,a^{\prime})$ can be calculated in
the following steps:

  1. 1.

Let the components of $\eta$ be defined as above. With the exception of
$\rho$, each can be estimated by cross-fitting a regression of the dependent
variable on the independent variables and generating predicted probabilities
(if the dependent variable is binary) or predicted values otherwise, setting
the values of independent variables where indicated. For example,
$\hat{\mathsf{g}}(a^{\prime}\mid w)$ can be estimated by cross-fitting a
logistic regression model of $A$ on $W$ and generating predicted probabilities
that $A=a^{\prime}$ for all observed $w$. One could also use machine learning
in model fitting, which is what we do in the simulations and data analysis.

  2. 2.

To estimate $\rho$,we use the predicted values from $\mu(a,m,z,W)$ as a
pseudo-outcome in a new regression on variables $A,Z,W$, and then generate
predicted values setting $A=a^{\prime}$ and $Z=z^{\prime}.$

  3. 3.

The estimator of $\theta(a,a^{\prime})$ is
$\hat{\theta}(a,a^{\prime})=\frac{1}{n}\sum_{i=1}^{n}D(O_{i},\hat{\eta}_{j(i)}).$

  4. 4.

The variance can be estimated as the sample variance of $D(O,\hat{\eta}).$

Let $\tilde{\eta}$ denote the probability limit of the estimator $\hat{\eta}$
in $L_{2}(\mathsf{P})$ norm. The above estimator is multiply robust such that
$\theta(a,a^{\prime})$ is expected to be estimated consistently under one of
the following:

  * •

$(\tilde{\mu},\tilde{\rho},\tilde{\mathsf{q}})=(\mu,\rho,\mathsf{q})$ or

  * •

$(\tilde{\mu},\tilde{\rho},\tilde{\mathsf{g}})=(\mu,\rho,\mathsf{g})$ or

  * •

$(\tilde{\mu},\tilde{\mathsf{g}},\tilde{\mathsf{q}})=(\mu,\mathsf{g},\mathsf{q})$
or

  * •

$(\tilde{\mathsf{g}},\tilde{\mathsf{e}},\tilde{\mathsf{q}},\tilde{\mathsf{r}})=(\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r}).$

The monotonicity assumption upon which this identification rests may be most
likely tenable in encouragement-design randomized trials like the MTO study
used in the illustrative application. In that case, estimation of $\mathsf{g}$
may not be necessary, so one would expect consistent estimation if either
$(\mu,\mathsf{q})$ or $(\mu,\rho)$ or $(\mathsf{e},\mathsf{q},\mathsf{r})$
were consistently estimated. These robustness conditions are a consequence of
the following Lemma, the proof of which is given in the appendix.

######  Lemma 1.

Let $D_{z,z^{\prime}}(O;\tilde{\eta})$ denote the EIF evaluated at a given
value $\tilde{\eta}$ of the nuisance parameters. Then we have

$\theta_{z,z^{\prime}}(a,a^{\prime})(\tilde{\eta})-\theta_{z,z^{\prime}}(a,a^{\prime})(\eta)=-\mathsf{E}[\bar{D}_{z,z^{\prime}}(O;\tilde{\eta})]+R_{z,z^{\prime}}(\eta,\tilde{\eta}),$

where $R_{z,z^{\prime}}(\eta,\tilde{\eta})$ is a second order term given by:

$\displaystyle R_{z,z^{\prime}}(\eta,\tilde{\eta})$
$\displaystyle=\mathsf{E}\big{\\{}[\tilde{H}_{Y,z,z^{\prime}}-H_{Y,z,z^{\prime}}][\mathsf{E}(Y\mid
A,M,Z,W)-\tilde{\mathsf{E}}(Y\mid A,M,Z,W)]\big{\\}}$
$\displaystyle+\mathsf{E}\left\\{[\tilde{H}_{M,z,z^{\prime}}-H_{M,z,z^{\prime}}]\int\tilde{\mathsf{E}}(Y\mid
a,M,z,W)[\mathrm{d}\mathsf{P}(m\mid
a^{\prime},z^{\prime},W)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},z^{\prime},W)]\right\\}$
$\displaystyle+\mathsf{E}\big{\\{}[\tilde{H}_{Z,z,z^{\prime}}-H_{Z,z,z^{\prime}}][\mathsf{P}(Z=1\mid
A,W)-\tilde{\mathsf{P}}(Z=1\mid A,W)]\big{\\}}$
$\displaystyle+R^{\star}_{z,z^{\prime}}(\eta,\tilde{\eta}),$

where, for each $(z,z^{\prime})$,
$R^{\star}_{z,z^{\prime}}(\eta,\tilde{\eta})$ is equal to:

$\displaystyle R^{\star}_{1,1}(\eta,\tilde{\eta})$
$\displaystyle=\int[\mathsf{E}(Y\mid a,m,1,w)-\tilde{\mathsf{E}}(Y\mid
a,m,1,w)][\mathsf{P}(Z=1\mid a^{\prime},w)-\tilde{\mathsf{P}}(Z=1\mid
a^{\prime},w)]\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},1,w)\mathrm{d}\mathsf{P}(w)$ $\displaystyle+\int\mathsf{E}(Y\mid
a,m,1,w)[\mathsf{P}(Z=1\mid a^{\prime},w)-\tilde{\mathsf{P}}(Z=1\mid
a^{\prime},w)][\mathrm{d}\mathsf{P}(m\mid
a^{\prime},1,w)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},1,w)]\mathrm{d}\mathsf{P}(w)$ $\displaystyle
R^{\star}_{1,0}(\eta,\tilde{\eta})$ $\displaystyle=\int[\mathsf{E}(Y\mid
a,m,1,w)-\tilde{\mathsf{E}}(Y\mid
a,m,1,w)][\mathsf{q}_{Z}(w)-\tilde{\mathsf{q}}_{Z}(w)]\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},0,w)\mathrm{d}\mathsf{P}(w)$ $\displaystyle+\int\mathsf{E}(Y\mid
a,m,1,w)[\mathsf{q}_{Z}(w)-\tilde{\mathsf{q}}_{Z}(w)][\mathrm{d}\mathsf{P}(m\mid
a^{\prime},0,w)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},0,w)]\mathrm{d}\mathsf{P}(w)$ $\displaystyle
R^{\star}_{0,0}(\eta,\tilde{\eta})$ $\displaystyle=\int[\mathsf{E}(Y\mid
a,m,0,w)-\tilde{\mathsf{E}}(Y\mid a,m,0,w)][\mathsf{P}(Z=1\mid
a^{\prime},w)-\tilde{\mathsf{P}}(Z=0\mid
a^{\prime},w)]\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},0,w)\mathrm{d}\mathsf{P}(w)$ $\displaystyle+\int\mathsf{E}(Y\mid
a,m,0,w)[\mathsf{P}(Z=0\mid a^{\prime},w)-\tilde{\mathsf{P}}(Z=0\mid
a^{\prime},w)][\mathrm{d}\mathsf{P}(m\mid
a^{\prime},0,w)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},0,w)]\mathrm{d}\mathsf{P}(w),$

for $\mathsf{q}_{Z}(w)=\mathsf{P}(Z=1\mid a,w)-\mathsf{P}(Z=1\mid
a^{\prime},w)$.

Inspection of the remainder terms reveals that the one-step estimator should
be consistent in the configurations outlined above, as in such cases all the
remainder terms become null. In addition, this lemma will allow us to prove an
asymptotic normality result stating that the estimator converges in
distribution to a normal random variable with variance equal to the non-
parametric efficiency bound. Importantly, this result holds even when the
nuisance parameters $\eta$ are estimated using using flexible regression
techniques such as those available in the machine learning literature, as long
as the regressions are cross-fitted as detailed above. The only requirement is
that the second order terms of the previous lemma converge to zero in
probability at rate $n^{-1/2}$. This can occur, for example, if all the
nuisance regressions are estimated consistently at rate $n^{-1/4}$. This rate
is much slower than the convergence rate of parametric models and is
achievable under a flexible regression framework.

######  Theorem 1 (Asymptotic normality of estimators).

Let $R=R_{1,1}+R_{1,0}+R_{0,0}$. Assume $R(\eta,\hat{\eta})=o_{P}(n^{-1/2})$.
Let $\hat{H}_{Y,z,z^{\prime}}$, $\hat{H}_{M,z,z^{\prime}}$,
$\hat{H}_{Z,z,z^{\prime}}$ denote the estimates of $H_{Y,z,z^{\prime}}$,
$H_{M,z,z^{\prime}}$, $H_{Z,z,z^{\prime}}$ constructed by plugging in
estimates $\hat{\eta}$. Assume that there exists a constant $c$ such that
$\mathsf{P}(\hat{H}_{Y,z,z^{\prime}}<c)=\mathsf{P}(\hat{H}_{M,z,z^{\prime}}<c)=\mathsf{P}(\hat{H}_{Z,z,z^{\prime}}<c)=1$.
Then we have

$\hat{\theta}(a,a^{\prime})-\theta(a,a^{\prime})=\frac{1}{n}\sum_{i=1}^{n}\bar{D}(O_{i},\eta)+o_{P}(n^{-1/2}).$

The proof of this result for estimators that allow a second order expansion as
in Lemma 1 follows standard empirical process theory and is given, for
example, in Díaz et al. (2021).

This result together with the central limit theorem shows that
$n^{1/2}(\hat{\theta}(a,a^{\prime})-\theta(a,a^{\prime}))$ converges to a
normal distribution with variance
$\operatorname{\mathsf{Var}}[\bar{D}(O,\eta)]$ equal to the non-parametric
efficiency bound. Furthermore, this result is useful to compute standard
errors and Wald-type confidence intervals for contrasts of the parameter
$\theta(a,a^{\prime})$ for varying values of $(a,a^{\prime})$, by simply
applying the Delta method to the corresponding contrast.

##  4 Simulation

We conducted a limited simulation study to: 1) verify that our programmed
estimator had the theoretical properties of: a) consistency under the
robustness conditions established in Lemma 1, b) efficiency, and c) confidence
interval coverage; and 2) illustrate the estimator’s finite sample
performance. We caution that these simulations should not be taken as
representative of the estimator’s general performance.

We considered the following data-generating mechanism:

$\displaystyle P(W_{1}=1)$ $\displaystyle=0.6$ $\displaystyle P(W_{2}=1)$
$\displaystyle=0.3$ $\displaystyle P(W_{3}=1\mid W_{1},W_{2})$
$\displaystyle=\min(0.2+0.33(W_{1}+W_{2}),1)$ $\displaystyle P(A=1)$
$\displaystyle=0.5$ $\displaystyle P(Z=1\mid A,W)$
$\displaystyle=\operatorname{expit}(-\log(1.3)\times(W_{1}+W_{2}+W_{3})/3+2A-1)$
$\displaystyle P(M=1\mid Z,W)$
$\displaystyle=\operatorname{expit}(-\log(1.1)W_{3}+2Z-0.9)$ $\displaystyle
P(Y=1\mid M,Z,W)$
$\displaystyle=\operatorname{expit}(-\log(1.3)\times(W_{1}+W_{2}+W_{3})/3+Z+M)$

This DGM was constructed to align with the illustrative application in which:
$A$ is randomly assigned and adheres to the exclusion restriction assumption
by only having an effect on $M$ and $Y$ through $Z$ (Angrist et al., 1996),
and $Z$ is monotonic with respect to $A$.

We conducted 1,000 simulations for sample sizes $n\in\\{1,000,10,000\\}$ and
under correct specification of nuisance parameters in $\eta$ and under several
misspecifications of parameters in $\eta$: 1) where
$(\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r})$ are correct but $(\mu,\rho)$
are misspecified, 2) where $(\mu,\rho,\mathsf{g})$ are correct but
$(\mathsf{e},\mathsf{q},\mathsf{r})$ are misspecified, and 3) where
$(\mu,\rho,\mathsf{q})$ are correct but $(\mathsf{g},\mathsf{e},\mathsf{r})$
are misspecified. We would still expect consistent estimation of
$\theta(a,a^{\prime})$ under each of these misspecifications given the
robustness results that are a consequence of Lemma 1. Also, note that $A$ is
exogenous in this simulation, so although we estimate $\mathsf{g}$ in this
simulation, its estimation is not necessary.

We fit parameters with the highly adaptive lasso (Hejazi et al., 2020;
Benkeser and van der Laan, 2016; van der Laan, 2017) for the correctly
specified scenarios, and fit parameters using an intercept-only model for the
incorrectly specified scenarios.

We considered estimator performance in terms of absolute bias, absolute bias
scaled by $\sqrt{n}$, influence curve-based standard error relative to the
Monte Carlo-based standard error, standard deviation of the estimator relative
to the efficiency bound scaled by $\sqrt{n}$, mean squared error relative to
the efficiency bound scaled by $n$, and 95% confidence interval (CI) coverage.

Table 1 shows simulation results for the natural direct effect, and Table 2
shows simulation results for the natural indirect effect under correct
specification of all nuisance parameters and various misspecifications,
comparing sample sizes of 10,000 and 1,000. We expect consistent estimation in
all scenarios. Focusing on the sample size of 10,000 and correct specification
of all nuisance parameters, we see very little bias; relative standard error,
relative standard deviation, and relative mean squared error close to 1; and
95% confidence interval coverage close to 95% for both the natural direct
effect and natural indirect effect. For this same sample size of 10,000 but
misspecification of some parameters, we maintain relatively consistent
estimation, with bias increasing only slightly. However, the relative standard
error, relative standard deviation, and relative mean squared error decrease
markedly in both the $(\mu,\rho,g)$ and $(\mu,\rho,q)$ correct specifications,
particularly in the case of the natural indirect effect. The 95% confidence
interval coverage also decreases under these misspecified scenarios,
particularly again, for the natural indirect effect. In fact, the 95% CI
covers only 46% of the time in the $(\mu,\rho,g)$ correct scenario for the
natural indirect effect. The results under sample size of 1,000 are similar,
albeit slightly worse on average, compared to the larger sample size.

Table 1: Simulation results for the natural direct effect, $\theta(1,0)-\theta(0,0)=0.1036$, efficiency bound= 1.7858. Abbreviations: relse = influence curve-based standard error relative to the Monte-Carlo-based standard error, relsd = standard deviation of the estimator relative to the efficiency bound scaled by $\sqrt{n}$, relrmse = mean squared error relative to the efficiency bound scaled by $n$, 95%CI Cov = 95% confidence interval coverage. Correctly specified nuisance parameters |  $|\text{bias}|$ |  $\sqrt{n}|\text{bias}|$ |  relse |  relsd |  relrmse |  95%CI Cov
---|---|---|---|---|---|---
N=10,000
All |  0.0004 |  0.0406 |  0.972 |  0.982 |  0.982 |  0.945
$\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r}$ |  0.0012 |  0.124 |  1.01 |  0.964 |  0.968 |  0.959
$\mu,\rho,\mathsf{g}$ |  0.0011 |  0.112 |  0.806 |  0.857 |  0.861 |  0.885
$\mu,\rho,\mathsf{q}$ |  0.0010 |  0.100 |  0.966 |  0.879 |  0.882 |  0.938
N=1,000
All |  0.0019 |  0.0592 |  0.9302 |  0.9737 |  0.9743 |  0.9273
$\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r}$ |  0.0047 |  0.1477 |  0.9941 |  0.9339 |  0.9400 |  0.9530
$\mu,\rho,\mathsf{g}$ |  0.0018 |  0.0567 |  0.7962 |  0.8743 |  0.8749 |  0.8768
$\mu,\rho,\mathsf{q}$ |  0.0070 |  0.2199 |  0.9449 |  0.8957 |  0.9103 |  0.9303
Table 2: Simulation results for the natural indirect effect, $\theta(1,1)-\theta(1,0)=0.0827$, efficiency bound= 0.9293. Abbreviations: relse = influence curve-based standard error relative to the Monte-Carlo-based standard error, relsd = standard deviation of the estimator relative to the efficiency bound scaled by $\sqrt{n}$, relrmse = mean squared error relative to the efficiency bound scaled by $n$, 95%CI Cov = 95% confidence interval coverage. Nuisance Parameters Correctly Specified |  $|\text{bias}|$ |  $\sqrt{n}|\text{bias}|$ |  relse |  relsd |  relrmse |  95%CI Cov
---|---|---|---|---|---|---
N=10,000
Correct |  0.0004 |  0.0412 |  0.951 |  0.974 |  0.975 |  0.932
$\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r}$ |  0.0012 |  0.116 |  0.988 |  0.916 |  0.924 |  0.95
$\mu,\rho,\mathsf{g}$ |  0.0012 |  0.120 |  0.314 |  0.727 |  0.738 |  0.464
$\mu,\rho,\mathsf{q}$ |  0.0010 |  0.0993 |  0.634 |  0.776 |  0.782 |  0.791
N=1,000
Correct |  0.0021 |  0.0672 |  0.8628 |  0.9743 |  0.9763 |  0.9010
$\mathsf{g},\mathsf{e},\mathsf{q},\mathsf{r}$ |  0.0038 |  0.1200 |  0.9154 |  0.8786 |  0.8870 |  0.9270
$\mu,\rho,\mathsf{g}$ |  0.0027 |  0.0867 |  0.3220 |  0.7412 |  0.7463 |  0.4859
$\mu,\rho,\mathsf{q}$ |  0.0067 |  0.2131 |  0.6157 |  0.7922 |  0.8222 |  0.7444

##  5 Empirical Illustration

Next, we apply our proposed estimators to estimate natural direct and indirect
effects among adolescent boys whose families participated in the Moving to
Opportunity Study (MTO) (Sanbonmatsu et al., 2011). Briefly, MTO randomized
families living in public housing in five cities in the United States (US),
who volunteered to participate, to either receive a Section 8 housing voucher
or not. (Families were actually randomized to one of three groups—two of which
received a housing voucher, but we collapse these for simplicity, as others
have done (Rudolph et al., 2021, 2018b, 2018a; Osypuk et al., 2012).) Section
8 housing vouchers are the primary form of federal housing assistance, and
they serve to subsidize rents on the private rental market for low-income
households, thus making it more financially viable for a family to move out of
public housing. In this example, we are interested in decomposing the average
effect of being randomized to receive a Section 8 housing voucher ($A$) in
early childhood on 10-15-year risk of developing any psychiatric disorder by
adolescence ($Y$) (defined using the Diagnostic and Statistical Manual,
Version IV (DSM-IV) criteria, as measured by the CIDI-SF (Kessler et al.,
1998; Kessler and Üstün, 2004)) into the portion that operates through changes
in the school environment and residential and school instability ($M$)—the
natural indirect effect—and the portion that does not—the natural direct
effect.

MTO is an example of an encouragement intervention, because receipt of a
Section 8 housing voucher (the randomized intervention) encourages families to
move out of public housing and into a rental on the private market by
subsidizing their rent. In the context of an encouragement intervention study
design, the monotonicity assumption is reasonable. Randomized receipt of the
encouragement—in this case, the housing voucher—should only serve to increase
(as opposed to increasing for some and decreasing for others) the likelihood
of complying with the intervention—in this case, moving within the first year
to a lower poverty neighborhood ($Z$).

In this example, we only consider boys whose families participated in MTO, as
there have been marked differences in MTO’s health effects between girls and
boys in which being in the intervention group resulted in higher rates of
post-traumatic stress disorder, mood disorder, any psychiatric disorder,
smoking, and problematic drug use for boys, but not for girls (Clampet-
Lundquist et al., 2011; Sanbonmatsu et al., 2011; Schmidt et al., 2017;
Kessler et al., 2014; Rudolph et al., 2018b; Kling et al., 2007). We also
exclude the Baltimore site, because of evidence that the intervention was
meaningfully different in this city versus the others (Rudolph et al., 2018b),
likely due to concurrent housing related-interventions in that city. These
restrictions resulted in a rounded sample size of N=2,100. We adjust for
numerous baseline covariates, $W$, which we detail in the appendix. We provide
additional detail of mediator variables $M$ in the appendix as well. For
simplicity, we used one imputed dataset, which was imputed using multiple
imputation by chained equations (Buuren and Groothuis-Oudshoorn, 2010). The
mediator of school rank had the most missingness at 12%, other mediators had
8-9% missingness or no missingness. The outcome of any DSM-IV disorder had 8%
missingness. The randomized voucher assignment variable and moving variable
had no missingness. Two baseline covariates had 2% missing data
(race/ethnicity and baseline neighborhood poverty), and the rest had no
missing data.

We estimated the natural direct and indirect effects of being randomized to
the Section 8 voucher group on risk of having a psychiatric disorder disorder
in adolescence, 10-15 years later, mediated through features of the
neighborhood and school environments (in the case of the indirect effect), and
not (in the case of the direct effect), using the cross-fitted one-step
estimator proposed here, with 2 folds. We used the SuperLearner ensemble
method of combining machine learning algorithms in fitting the nuisance
parameters, implemented with the sl3 package (Coyle et al., 2021); this
approach weights the algorithms to minimize the 10-fold cross-validated
prediction error (Van der Laan et al., 2007). We included the following
algorithms: intercept-only regression, generalized linear regression, lasso
(Tibshirani, 1996), and gradient boosted machines (Chen and Guestrin, 2016).
Columbia University determined this analysis of deidentified data to be non-
human subjects research.

Figure 1 shows the point estimates and 95% confidence intervals for the
natural direct effect and natural indirect effect. While the natural direct
effect appears null, the natural indirect effect shows evidence of an
unintended harmful path from being randomized to the housing voucher group on
later risk of developing a psychiatric disorder operating through voucher-
induced changes in the school and neighborhood environments and the
instability of those environments, among boys. This indirect effect was
estimated to contribute to a 3.04 percentage point increased risk of
developing such a disorder (95% CI: 0.36, 5.72 percentage points), which is
in-line with previous estimates using population interventional indirect
effects as the estimand (Rudolph et al., 2021).

Figure 1: Natural direct and indirect effects estimates of being randomized to
the Section 8 voucher group on risk of having a psychiatric disorder disorder
in adolescence, 10-15 years later, mediated through features of the
neighborhood and school environments, among boys. Estimates and 95% CIs. All
results were approved for release by the U.S. Census Bureau, authorization
number CBDRB-FY22-CES018-010.

##  6 Conclusion

Although natural direct and indirect effects are not generally identifiable in
the presence of treatment-induced confounders, assuming a monotonic
relationship between the treatment and treatment-induced confounder achieves
identifiability, as shown by Tchetgen and VanderWeele (2014). However, we are
not aware of epidemiologists or other applied researchers using this
identification strategy in practice. One reason could be that such scenarios
are believed to be rare, niche circumstances. However, the monotonicity
assumption is likely satisfied in any trial or intervention where the
intervention assigns individuals to a treatment condition, but it is up to the
individual whether or not to comply with the treatment they were assigned.
Indeed, such study designs are common, and have been the setting for much
mediation-related research (Vo et al., 2020). In this encouragement design
setting, intervention take-up or adherence becomes the treatment-induced
confounder, and one can then make the assumption that treatment assignment is
monotonic with treatment take-up. If the remaining treatment-induced variables
are treated as mediators, as we do in the applied example, then one can use
our proposed one-step estimator. Our cross-fitted estimator is multiply robust
and efficient, incorporates data-adaptive machine learning into model-fitting,
and is available as an R package linkblindedforreview.

## References

  * Angrist et al. (1996) J. Angrist, G. Imbens, and D.B. Rubin.  Identification of causal effects using instrumental variables.  _Journal of the American Statistical Association_ , 91:444–555, 1996.
  * Avin et al. (2005) Chen Avin, Ilya Shpitser, and Judea Pearl.  Identifiability of path-specific effects.  In _IJCAI International Joint Conference on Artificial Intelligence_ , pages 357–363, 2005.
  * Benkeser and van der Laan (2016) David Benkeser and Mark van der Laan.  The highly adaptive lasso estimator.  In _2016 IEEE International Conference on Data Science and Advanced Analytics (DSAA)_ , pages 689–696. IEEE, 2016.
  * Buuren and Groothuis-Oudshoorn (2010) S van Buuren and Karin Groothuis-Oudshoorn.  mice: Multivariate imputation by chained equations in r.  _Journal of statistical software_ , pages 1–68, 2010.
  * Chen and Guestrin (2016) Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939785.  URL http://doi.acm.org/10.1145/2939672.2939785.
  * Chernozhukov et al. (2016) Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian Hansen, et al.  Double machine learning for treatment and causal parameters.  _arXiv preprint arXiv:1608.00060_ , 2016.
  * Clampet-Lundquist et al. (2011) Susan Clampet-Lundquist, Kathryn Edin, Jeffrey R Kling, and Greg Duncan.  Moving teenagers out of high-risk neighborhoods: How girls fare better than boys.  _American Journal of Sociology_ , 116(4):1154–1189, 2011.
  * Coyle et al. (2021) Jeremy R Coyle, Nima S Hejazi, Ivana Malenica, Rachael V Phillips, and Oleg Sofrygin.  sl3: Modern pipelines for machine learning and Super Learning.  https://github.com/tlverse/sl3, 2021.  URL https://doi.org/10.5281/zenodo.1342293.  R package version 1.4.2.
  * Díaz et al. (2021) Iván Díaz, Nima S Hejazi, Kara E Rudolph, and Mark J van Der Laan.  Nonparametric efficient causal mediation with intermediate confounders.  _Biometrika_ , 108(3):627–641, 2021.
  * Hejazi et al. (2020) Nima S Hejazi, Jeremy R Coyle, and Mark J van der Laan.  hal9001: Scalable highly adaptive lasso regression inr.  _Journal of Open Source Software_ , 5(53):2526, 2020.
  * Kessler and Üstün (2004) Ronald C Kessler and T Bedirhan Üstün.  The world mental health (wmh) survey initiative version of the world health organization (who) composite international diagnostic interview (cidi).  _International Journal of Methods in Psychiatric Research_ , 13(2):93–121, 2004.
  * Kessler et al. (1998) Ronald C Kessler, Gavin Andrews, Daniel Mroczek, Bedirhan Ustun, and Hans-Ulrich Wittchen.  The world health organization composite international diagnostic interview short-form (cidi-sf).  _International journal of methods in psychiatric research_ , 7(4):171–185, 1998.
  * Kessler et al. (2014) Ronald C Kessler, Greg J Duncan, Lisa A Gennetian, Lawrence F Katz, Jeffrey R Kling, Nancy A Sampson, Lisa Sanbonmatsu, Alan M Zaslavsky, and Jens Ludwig.  Associations of housing mobility interventions for children in high-poverty neighborhoods with subsequent mental disorders during adolescence.  _JAMA_ , 311(9):937–948, 2014.
  * Klaassen (1987) Chris AJ Klaassen.  Consistent estimation of the influence function of locally asymptotically linear estimators.  _The Annals of Statistics_ , pages 1548–1562, 1987.
  * Kling et al. (2007) Jeffrey R Kling, Jeffrey B Liebman, and Lawrence F Katz.  Experimental analysis of neighborhood effects.  _Econometrica_ , 75(1):83–119, 2007.
  * Miles (2022) Caleb H Miles.  On the causal interpretation of randomized interventional indirect effects.  _arXiv preprint arXiv:2203.00245_ , 2022.
  * Osypuk et al. (2012) Theresa L Osypuk, Nicole M Schmidt, Lisa M Bates, Eric J Tchetgen-Tchetgen, Felton J Earls, and M Maria Glymour.  Gender and crime victimization modify neighborhood effects on adolescent mental health.  _Pediatrics_ , 130(3):472–481, 2012.
  * Pearl (2001) Judea Pearl.  Direct and indirect effects.  In _Proceedings of the seventeenth conference on uncertainty in artificial intelligence_ , pages 411–420. Morgan Kaufmann, 2001.
  * Pearl (2009) Judea Pearl.  Myth, Confusion, and Science in Causal Analysis.  Technical Report R-348, Cognitive Systems Laboratory, Computer Science Department University of California, Los Angeles, Los Angeles, CA, May 2009.
  * Pearl (2013) Judea Pearl.  Direct and indirect effects.  _arXiv preprint arXiv:1301.2300_ , 2013.
  * Rudolph et al. (2018a) Kara E Rudolph, Nicole M Schmidt, M Maria Glymour, Rebecca E Crowder, Jessica Galin, Jennifer Ahern, and Theresa L Osypuk.  Composition or context: Using transportability to understand drivers of site differences in a large-scale housing experiment.  _Epidemiology_ , 29(2):199–206, 2018a.
  * Rudolph et al. (2018b) Kara E Rudolph, Oleg Sofrygin, Nicole M Schmidt, Rebecca Crowder, M Maria Glymour, Jennifer Ahern, and Theresa L Osypuk.  Mediation of neighborhood effects on adolescent substance use by the school and peer environments.  _Epidemiology_ , 29(4):590–598, 2018b.
  * Rudolph et al. (2021) Kara E Rudolph, Catherine Gimbrone, and Iván Díaz.  Helped into harm: Mediation of a housing voucher intervention on mental health and substance use in boys.  _Epidemiology_ , 32(3):336–346, 2021.
  * Sanbonmatsu et al. (2011) Lisa Sanbonmatsu, Jens Ludwig, Lawrence F Katz, Lisa A Gennetian, Greg J Duncan, Ronald C Kessler, Emma Adam, Thomas W McDade, and Stacy Tessler Lindau.  _Moving to Opportunity for Fair Housing Demonstration Program–Final Impacts Evaluation_.  US Department of Housing and Urban Development, Office of Policy Development and Research, Washington, DC, 2011.
  * Schmidt et al. (2017) Nicole M Schmidt, M Maria Glymour, and Theresa L Osypuk.  Adolescence is a sensitive period for housing mobility to influence risky behaviors: An experimental design.  _Journal of Adolescent Health_ , 60(4):431–437, 2017.
  * Tchetgen and VanderWeele (2014) Eric J Tchetgen Tchetgen and Tyler J VanderWeele.  On identification of natural direct effects when a confounder of the mediator is directly affected by exposure.  _Epidemiology (Cambridge, Mass.)_ , 25(2):282, 2014.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society. Series B (Methodological)_ , pages 267–288, 1996.
  * van der Laan (2017) Mark van der Laan.  A generally efficient targeted minimum loss based estimator based on the highly adaptive lasso.  _The international journal of biostatistics_ , 13(2), 2017\.
  * Van der Laan et al. (2007) Mark J Van der Laan, Eric C Polley, and Alan E Hubbard.  Super learner.  _Statistical Applications in Genetics and Molecular Biology_ , 6(1), 2007.
  * Vo et al. (2020) Tat-Thang Vo, Cecilia Superchi, Isabelle Boutron, and Stijn Vansteelandt.  The conduct and reporting of mediation analysis in recently published randomized controlled trials: results from a methodological systematic review.  _Journal of clinical epidemiology_ , 117:78–88, 2020.
  * Zheng and van der Laan (2011) Wenjing Zheng and Mark J van der Laan.  Cross-validated targeted minimum-loss-based estimation.  In _Targeted Learning_ , pages 459–474. Springer, 2011.

Supplementary Materials for Efficient and flexible estimation of natural
mediation effects under intermediate confounding and monotonicity constraints

##  Appendix S1 Proofs

Proof  The proof proceeds as follows. Define
$R_{z,z^{\prime}}(\eta,\tilde{\eta})=\theta_{z,z^{\prime}}(\tilde{\eta})-\theta_{z,z^{\prime}}(\eta)+\mathsf{E}[D_{z,z^{\prime}}(O;\tilde{\eta})]$.
This is equal to:

$\displaystyle R_{z,z^{\prime}}(\eta,\tilde{\eta})$
$\displaystyle=E\big{\\{}[\tilde{H}_{Y,z,z^{\prime}}-H_{Y,z,z^{\prime}}][\mathsf{E}(Y\mid
A,M,Z,W)-\tilde{\mathsf{E}}(Y\mid A,M,Z,W)]\big{\\}}$
$\displaystyle+\mathsf{E}\left\\{[\tilde{H}_{M,z,z^{\prime}}-H_{M,z,z^{\prime}}]\int\tilde{\mathsf{E}}(Y\mid
a,M,z,W)[\mathrm{d}\mathsf{P}(m\mid
a^{\prime},z^{\prime},W)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},z^{\prime},W)]\right\\}$
$\displaystyle+\mathsf{E}\big{\\{}[\tilde{H}_{Z,z,z^{\prime}}-H_{Z,z,z^{\prime}}][\mathsf{P}(Z=1\mid
A,W)-\tilde{\mathsf{P}}(Z=1\mid A,W)]\big{\\}}$
$\displaystyle+\mathsf{E}\big{\\{}H_{Y,z,z^{\prime}}[\mathsf{E}(Y\mid
A,M,Z,W)-\tilde{\mathsf{E}}(Y\mid A,M,Z,W)]\big{\\}}$ (1)
$\displaystyle+\mathsf{E}\left\\{H_{Z,z,z^{\prime}}\int\tilde{\mathsf{E}}(Y\mid
a,M,z,W)[\mathrm{d}\mathsf{P}(m\mid
a^{\prime},z^{\prime},W)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},z^{\prime},W)]\right\\}$ (2)
$\displaystyle+\mathsf{E}\big{\\{}H_{Z,z,z^{\prime}}[\mathsf{P}(Z=1\mid
A,W)-\tilde{\mathsf{P}}(Z=1\mid A,W)]\big{\\}}$ (3)
$\displaystyle+\mathsf{E}\big{\\{}\tilde{H}_{W,z,z^{\prime}}-H_{W,z,z^{\prime}}\big{\\}}.$
(4)

It remains to show that the sum of terms (1)-(4) is equal to
$R^{\star}_{z,z^{\prime}}(\eta,\tilde{\eta})$. We show the proof for
$(z,z^{\prime})=(1,1)$, the proofs for the other terms follow similar
arguments. Term (1) is equal to

$\displaystyle(\ref{eq:term1})$ $\displaystyle=\int
H_{Y}(\tilde{a},m,z,w)[\mathsf{E}(Y\mid a,m,1,w)-\tilde{\mathsf{E}}(Y\mid
a,m,1,w)]\mathrm{d}\mathsf{P}(\tilde{a},m,z,w)$
$\displaystyle=\int\frac{\mathds{1}\\{z=1,\tilde{a}=a\\}}{\mathsf{P}(a^{\prime}\mid
w)}\frac{\mathsf{P}(a^{\prime}\mid m,z,w)}{\mathsf{P}(\tilde{a}\mid
m,z,w)}[\mathsf{E}(Y\mid a,m,1,w)-\tilde{\mathsf{E}}(Y\mid
a,m,1,w)]\mathrm{d}\mathsf{P}(\tilde{a}\mid m,z,w)\mathrm{d}\mathsf{P}(m,z,w)$
$\displaystyle=\int\frac{\mathds{1}\\{z=1\\}}{\mathsf{P}(a^{\prime}\mid
w)}\mathsf{P}(a^{\prime}\mid m,z,w)[\mathsf{E}(Y\mid
a,m,1,w)-\tilde{\mathsf{E}}(Y\mid a,m,1,w)]\mathrm{d}\mathsf{P}(m,z,w)$
$\displaystyle=\int\frac{\mathds{1}\\{z=1\\}}{\mathsf{P}(a^{\prime}\mid
w)}\frac{\mathsf{P}(a^{\prime},m,z,w)}{\mathsf{p}(m,z,w)}[\mathsf{E}(Y\mid
a,m,1,w)-\tilde{\mathsf{E}}(Y\mid a,m,1,w)]\mathrm{d}\mathsf{P}(m,z,w)$
$\displaystyle=\int\frac{\mathds{1}\\{z=1\\}}{\mathsf{P}(a^{\prime}\mid
w)}\mathsf{P}(a^{\prime},m,z,w)[\mathsf{E}(Y\mid
a,m,1,w)-\tilde{\mathsf{E}}(Y\mid a,m,1,w)]\mathrm{d}\nu(m,z,w)$
$\displaystyle=\int\mathsf{P}(Z=1\mid a^{\prime},w)[\mathsf{E}(Y\mid
a,m,1,w)-\tilde{\mathsf{E}}(Y\mid a,m,1,w)]\mathrm{d}\mathsf{P}(m\mid
a^{\prime},1,w)\mathrm{d}\mathsf{P}(w)$

Term (2) equals

$\displaystyle(\ref{eq:term2})$ $\displaystyle=\int
H_{M}(\tilde{a},z,w)\tilde{\mathsf{E}}(Y\mid
a,m,1,w)[\mathrm{d}\mathsf{P}(m\mid
a^{\prime},1,w)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},1,w)]\mathrm{d}\mathsf{P}(\tilde{a},z,w)$
$\displaystyle=\int\frac{\mathds{1}\\{z=1,\tilde{a}=a^{\prime}\\}}{\mathsf{P}(a^{\prime}\mid
w)}\tilde{\mathsf{E}}(Y\mid a,m,1,w)[\mathrm{d}\mathsf{P}(m\mid
a^{\prime},1,w)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},1,w)]\mathrm{d}\mathsf{P}(\tilde{a},z,w)$
$\displaystyle=\int\mathsf{P}(Z=1\mid a^{\prime},w)\tilde{\mathsf{E}}(Y\mid
a,m,1,w)[\mathrm{d}\mathsf{P}(m\mid
a^{\prime},1,w)-\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},1,w)]\mathrm{d}\mathsf{P}(w)$

Term (3) equals

$\displaystyle(\ref{eq:term3})$ $\displaystyle=\int
H_{Z}(a,z,w)[\mathsf{P}(Z=1\mid a^{\prime},w)-\tilde{\mathsf{P}}(Z=1\mid
a^{\prime},w)]\mathrm{d}\mathsf{P}(a,w)$ $\displaystyle=\int\mathsf{E}(Y\mid
a,m,1,w)\mathrm{d}\mathsf{P}(m\mid a^{\prime},1,w)[\mathsf{P}(Z=1\mid
a^{\prime},w)-\tilde{\mathsf{P}}(Z=1\mid
a^{\prime},w)]\mathrm{d}\mathsf{P}(w)$

And term (4) equals

$(\ref{eq:term4})=\int[\tilde{\mathsf{E}}(Y\mid
a,m,1,w)\mathrm{d}\tilde{\mathsf{P}}(m\mid
a^{\prime},1,w)\tilde{\mathsf{P}}(Z=1\mid a^{\prime},w)-\\\ \mathsf{E}(Y\mid
a,m,1,w)\mathrm{d}\mathsf{P}(m\mid a^{\prime},1,w)\mathsf{P}(Z=1\mid
a^{\prime},w)]\mathrm{d}\mathsf{P}(w).$

Factorizing these four terms gives the desired result. ∎

##  Appendix S2 Details of variables used in the empirical illustration

Baseline covariates used:

  * •

Adolescent characteristics: site (Boston, Chicago, LA, NYC), age,
race/ethnicity (categorized as black, latino/Hispanic, white, other), number
of family members (categorized as 2, 3, or 4+), someone from school asked to
discuss problems the child had with schoolwork or behavior during the 2 years
prior to baseline, child enrolled in special class for gifted and talented
students.

  * •

Adult household head characteristics included: high school graduate, marital
status (never vs ever married), whether had been a teen parent, work status,
receipt of AFDC/TANF, whether any family member has a disability.

  * •

Neighborhood characteristics: felt neighborhood streets were unsafe at night;
very dissatisfied with neighborhood; poverty level of neighborhood.

  * •

Reported reasons for participating in MTO: to have access to better schools.

  * •

Moving-related characteristics: moved more then 3 times during the 5 years
prior to baseline, previous application for Section 8 voucher.

Mediator variables were duration-weighted (i.e., calculated as a weighted mean
where weights were proportionate to the length of follow-up time (e.g., length
of time that the youth attended each school, lived in each neighborhood,
etc.)) in the 10-15 years between randomization and the final follow-up
timepoint when the outcome was measured. These variables included:

  * •

School rank

  * •

School student-to-teacher ratio

  * •

School’s percent of students receiving free or reduced lunch

  * •

Proportion of Title-I schools attended

  * •

Number of moves

  * •

Number of schools attended

  * •

Number of school changes within the school year

  * •

Whether or not the most recent school was in the same school district as
baseline

  * •

Neighborhood poverty


# Effective Model for Fractional Topological Corner Modes in Quasicrystals

Citian Wang  School of Physics, Peking University, Beijing 100871, China  Feng
Liu<EMAIL_ADDRESS>Department of Materials Science and Engineering,
University of Utah, Salt Lake City, Utah 84112, USA  Huaqing Huang
<EMAIL_ADDRESS>School of Physics, Peking University, Beijing 100871,
China  Collaborative Innovation Center of Quantum Matter, Beijing 100871,
China  Center for High Energy Physics, Peking University, Beijing 100871,
China

###### Abstract

High-order topological insulators (HOTIs), as generalized from topological
crystalline insulators (TCIs), are characterized with lower-dimensional
metallic boundary states protected by spatial symmetries of a crystal, whose
theoretical framework based on band inversion at special $k$-points cannot be
readily extended to quasicrystals because quasicrystals contain rotational
symmetries that are not compatible with crystals, and momentum is no longer a
good quantum number. Here, we develop a low-energy effective model underlying
HOTI states in 2D quasicrystals for all possible rotational symmetries. By
implementing a novel Fourier transform developed recently for quasicrystals
and approximating the long-wavelength behavior by their large-scale average,
we construct an effective $k\cdot p$ Hamiltonian to capture the band inversion
at the center of a pseudo-Brillouin zone (PBZ). We show that an in-plane
Zeeman field can induce mass-kinks at the intersection of adjacent edges of a
2D quasicrystal TI and generate corner modes (CMs) with fractional charge,
protected by rotational symmetries. Our model predictions are confirmed by
numerical tight-binding calculations. Furthermore, when the quasicrystal is
proximitized by an s-wave superconductor, Majorana CMs can also be created by
tuning the field strength and chemical potential. Our work affords a generic
approach to studying the low-energy physics of quasicrystals, in association
with topological excitations and fractional statistics.

Introduction.— A first-order $\mathbb{Z}_{2}$ TI is characterized with a
correspondence between a $d$-dimensional gapped bulk state and a
$(d-1)$-dimensional gapless boundary state protected by time-reversal (TR)
symmetry, while a HOTI is instead featured with a $(d-2)$-dimensional gapless
boundary state protected by spatial symmetries of a crystal, such as mirror
and rotation symmetry Benalcazar  _et al._ (2017); Schindler  _et al._
(2018a); Langbehn  _et al._ (2017); Song  _et al._ (2017); Schindler  _et al._
(2018b); Yue  _et al._ (2019); Xu  _et al._ (2019); Zhang  _et al._ (2020);
Sheng  _et al._ (2019); Lee  _et al._ (2020); Liu  _et al._ (2019a); Park  _et
al._ (2019); Liu  _et al._ (2021); Huang and Liu (2021); Pan  _et al._ (2022);
Ezawa (2018a); Liu  _et al._ (2019b); Wang  _et al._ (2019); Hsu  _et al._
(2019); Ren  _et al._ (2020); Chen  _et al._ (2020a); Zhao  _et al._ (2021).
In the context of crystalline symmetry protected topological boundary states,
the HOTI can be viewed as a generalization of TCI Fu (2011). In general,
$\mathbb{Z}_{2}$ invariant can be calculated by the product of band inversion
indices at the TR invariant momenta (TRIM) in the subspace of occupied bands,
the eigenvalues of inversion symmetry operator, in the first BZ Fu and Kane
(2007). One can generalize this scheme to identify a HOTI by calculating the
topological invariant, the eigenvalues of a spatial symmetry at all
$\mathbf{k}$-points linked by this spatial symmetry Benalcazar  _et al._
(2017); Schindler  _et al._ (2018a). Namely, a high-order topology is defined
in the subspace of a crystal operated by spatial symmetries. Apparently, this
approach is not applicable to quasicrystals which do not have a BZ and
momentum ($\mathbf{k}$) is no longer a good quantum number. Moreover,
quasicrystals contain rotation symmetries not compatible with translational
symmetry.

Alternatively, a HOTI can be viewed by gapping the $(d-1)$-dimensional gapless
boundary (surfaces or edges) of a $\mathbb{Z}_{2}$ TI, but the band degeneracy
is locally protected at the $(d-2)$-dimensional boundary (hinges and corners)
by spatial symmetry. For example, in a 2D HOTI, corner modes (CMs) can be
viewed as topological Jackiw-Rebbi domain-wall states Jackiw and Rebbi (1976),
with opposite Dirac masses between two edges enforced by a mirror symmetry
Geier  _et al._ (2018); Trifunovic and Brouwer (2019, 2020), and a variety of
2D HOTI systems have been proposed by implementing such mirror-invoked mass-
inversion mechanism Ezawa (2018a); Liu  _et al._ (2019b); Wang  _et al._
(2019); Hsu  _et al._ (2019); Ren  _et al._ (2020); Chen  _et al._ (2020a);
Zhao  _et al._ (2021), including interestingly CMs in quasicrystals Varjas
_et al._ (2019); Chen  _et al._ (2020b); Hua  _et al._ (2020); Huang  _et al._
(2021); Fan and Huang (2022). More generally, a CM of domain-wall state can be
protected by rotation symmetry, as derived from edge network theory Wang and
Moore (2019). Instead of opposite edge masses encoded by mirror eigenvalues
($\pm 1$), an edge-dependent Dirac masses emerges with a phase difference of
$2\pi/n$, termed as mass kink Goldstone and Wilczek (1981); Jackiw and
Semenoff (1983), defined by eigenvalues of the $C_{n}$ rotation
$(e^{\frac{2m\pi i}{n}},\,m=1,\dots,n)$, giving rise to a fractional charged
CMs of $e/n$. We note that the mass inversion mechanism, as also applied to
quasicrystals Varjas  _et al._ (2019); Chen  _et al._ (2020b); Hua  _et al._
(2020); Huang  _et al._ (2021); Fan and Huang (2022), is a special case of
mass link with $C_{2}$ rotation with a phase shift of $\pi$. Therefore, it is
very interesting to explore if the mass kink approach can be generalized to
quasicrystals, because they have rotation symmetries (such as 5-fold rotation)
that are not compatible with translational symmetry, which may hinder the
realization of topological CMs with fractional charges that do not exist in
crystals.

In this Letter, we develop a low-energy effective model for quasicrystals,
based on a novel Fourier transform developed recently Jiang and Zhang (2014)
by representing a quasicrystal as a projection of a hypercrystal from higher
dimensions, from a 4D hyper-BZ to a 2D PBZ. Then, an effective $k\cdot p$
Hamiltonian is constructed at the center ($\Gamma$) of the PBZ, under long-
wavelength approximation by large-scale average with quasicrystalline
symmetry. As band inversion occurs at $\Gamma$, charged and Majarona CMs
arising from fractional mass kinks can emerge by applying an in-plane Zeeman
field. Taking the Penrose-tiling quasicrystal TI as an example, we show that
the field induces a fractional mass kink with a phase shift of $2\pi/5$ at the
intersection of adjacent edges, generating five in-gap localized CMs with a
fractional charge of $e/5$, as displayed in Fig. 1. Our scheme can be easily
generalized to $q$-fold quasicrystals ($q=5,8,10,12,\cdots$), leading to the
corner charge fractionalization by $e/q$ that is disallowed in crystals. In
addition, when the 2D quasicrystal TI is in proximity with an s-wave SC,
Majorana CMs can be generated by tuning the in-plane Zeeman field and chemical
potential (see Fig. 1).

Figure 1:  Schematic illustration of the Zeeman-field-induced topological
phase transitions in a 2D pentagonal quasicrystal. Starting from a TI phase
with helical edge states, the quasicrystal is driven to a HOTI phase with five
charged corner modes (CMs) by an in-plane Zeeman field $\delta_{M}\bm{m}$.
When the system is in proximity with an s-wave superconductor (SC) with
pairing $\Delta_{0}$, Majorana CMs can be generated by tuning the Zeeman field
and chemical potential.

Model.— We construct quasicrystal lattices based on the rhombic Penrose and
Ammann-Beenker tilings which have 5- and 8-fold rotational symmetry,
respectively Penrose (1974); Beenker (1982). We assume the atoms have three
atomic orbitals ($s,p_{x},p_{y}$) at the vertices of tiling and consider
hoppings between neighboring vertices connected by edges or the shorter
diagonals of the rhombi. The tight-binding (TB) Hamiltonian is given by

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{i\alpha}\epsilon_{\alpha}c_{i\alpha}^{\dagger}c_{i\alpha}+\sum_{\langle
i\alpha,j\beta\rangle}t_{\alpha,\beta}(\mathbf{r}_{ij})c_{i\alpha}^{\dagger}c_{j\beta}$
$\displaystyle+i\lambda\sum_{i}(c_{ip_{y}}^{\dagger}s_{z}c_{ip_{x}}-c_{ip_{x}}^{\dagger}s_{z}c_{ip_{y}})+\sum_{i\alpha}\delta_{\alpha}c_{i\alpha}^{\dagger}(\bm{m}\cdot\bm{s})c_{i\alpha},$

where
$c_{i\alpha}^{\dagger}=(c_{i\alpha\uparrow}^{\dagger},c_{i\alpha\downarrow}^{\dagger})$
is electron creation operators on the $\alpha(=s,p_{x},p_{y})$ orbital at the
$i$-th site, $\epsilon_{\alpha}$ is the on-site energy, and
$t_{\alpha,\beta}(\mathbf{r}_{ij})$ is the Slater-Koster hopping integral
which depends on the orbital type ($\alpha$ and $\beta$) and the vector
$\mathbf{r}_{ij}$ between sites $i$ and $j$. $\lambda$ is the spin-orbit
coupling (SOC) strength and $\bm{s}=(s_{x},s_{y},s_{z})$ are the Pauli
matrices. The last term represents a Zeeman field along the direction of
$\bm{m}$. $\delta_{\alpha}$ depends on field strength and the g-factor of
$\alpha$ orbital. For simplicity, we take a uniform value
$\delta_{\alpha}=\delta_{M}$, which would not change the main physics we
discuss here. Experimentally, the Zeeman term can be introduced by a magnetic
field, coupling to a magnetic substrate, or depositing magnetic adatoms on
quasicrystal substrates. It is well known that by considering a band inversion
between s and ($p_{x}$, $p_{y}$) orbitals, topological states can be realized
in quasicrystal lattices Huang and Liu (2018a, b, 2019). Hence, we will use
the same settings and focus on $n_{occ}=2/3$ filling of electron states
hereafter.

Figure 2:  (a) The set of all reciprocal wave vectors consists of five wave
vectors $\mathbf{g}_{j}$, forming a dense set of points in k space. The
perpendicular bisectors of the 10 principal wave vectors $\pm\mathbf{g}_{j}$
form a decagonal boundary to the pseudo-Brillouin zone (PBZ). (b) The
effective band structure around the center of the PBZ of the Penrose-tiling
quasicrystal. (c) Schematic illustration of an edge along
$\vec{e}_{1}=\cos\zeta\hat{e}_{x}+\sin\zeta\hat{e}_{y}$ direction, where
$\zeta$ is the angle between the edge and positive x direction.

Low-energy effective theory.— Existing approaches based on the analysis of
states at high-symmetry $k$-points in the BZ of a crystal are ruled out
because momentum is no longer a good quantum number for quasicrystals.
Instead, in the following we elucidate the HOTI in quasicrystals by
establishing a low-energy effective theory in a continuum model and performing
analysis in terms of coupled edge modes.

Generally, with the quasiperiodicity, the Fourier transform of any function
$f(\mathbf{r})$ [e.g., the particle density $\rho(\mathbf{r})$ and
quasicrystalline potential $U(\mathbf{r})$] of a 2D pentagonal quasicrystal
can be expressed as

$f(\mathbf{r})=\sum_{\mathbf{G}\in\mathcal{L}}\hat{f}(\mathbf{G})e^{i\mathbf{G}\cdot\mathbf{r}},$
(2)

where $\mathcal{L}$ is a countable set of reciprocal wave vectors that consist
of $\mathbf{G}=\sum_{j=1}^{5}n_{j}\mathbf{g}_{j}$ ($n_{j}\in\mathbb{Z}$)
filling densely the 2D reciprocal space. The five principal reciprocal vectors
in the $k_{x}$-$k_{y}$ plane are

$\mathbf{g}_{j}=2\pi(\cos[2\pi(j-1)/5],\sin[2\pi(j-1)/5]),\quad j=1,\ldots,5,$
(3)

among which only four are independent due to the linear dependence of
$\sum_{j=1}^{5}\mathbf{g}_{j}=0$. Note that all the momentum values are scaled
in the unit of inversed bond length $a^{-1}$, and we set $a=1$ without loss of
generality. The wave vectors in $\mathcal{L}$ can be divided into $n$-th order
with $n=\sum_{j=1}^{5}|n_{j}|$, corresponding to the order of diffraction
peaks [See Fig. 2(a)]. The first order of $\mathcal{L}$, which contains ten
principal wave vectors as $\pm\mathbf{g}_{j}$, defines the PBZ Janssen  _et
al._ (2007), in analogy to the conventional first BZ of crystals. According to
the gap labelling Bellissard (1992); Gambaudo and Vignolo (2014), the leading-
order gap opens at the boundary of the PBZ with the gap size determined by the
first-order Fourier coefficient $\hat{U}(\pm\mathbf{g}_{j})$ of
quasicrystalline potential. Whereas nearby the $\Gamma$ point of the PBZ, the
higher-order gaps appear in a hierarchy, which corresponds to multiple
scattering processes having a rapidly decreasing size. This enables a
continuum description based on a low-energy effective model with a proper
truncation of $\mathcal{L}$.

Specifically, we employ the projection method, for which a quasicrystal can be
viewed as a higher-dimensional crystal in hyperspace Steurer and Deloudi
(2009). Then, an equivalent 4D representation of the 2D quasicrystal can be
made by implementing the novel Fourier expansion proposed by Jiang and Zhang
(JZ) Jiang and Zhang (2014),

$f(\mathbf{r})=\sum_{\mathbf{\Pi}}\hat{f}(\mathbf{\Pi})e^{i\left[(\mathcal{S}\cdot\mathbf{\Pi})^{T}\cdot\mathbf{r}\right]},$
(4)

where $\mathbf{\Pi}=\sum_{i=1}^{4}m_{i}\mathbf{Q}_{i}\in\mathbb{R}^{4}$,
$m_{i}\in\mathbb{Z}$, and $\mathbf{Q}_{i}$ are the 4D primitive reciprocal
vectors of the hypercrystal. $\mathcal{S}$ is the projection matrix connecting
the 2D physical space with the 4D hyperspace. Mathematically, for any 2D
quasiperiodic function $f(\mathbf{r})$, Eq. (4) has the property Jiang and
Zhang (2014)

$\lim_{V\rightarrow\infty}\frac{1}{V}\int
d\mathbf{r}f(\mathbf{r})=\left.\hat{f}(\mathbf{\Pi})\right|_{\mathbf{\Pi}=0}.$
(5)

Consequently, the large-scale average property of quasicrystals is well
captured by the contribution around
$\mathbf{k}=\mathcal{S}\cdot\mathbf{\Pi}=0$, i.e., the $\Gamma$ point of the
PBZ. In fact, it has been proved that the behavior of a quasicrystal for
excitations of any kind with long-wavelength modes can be related to its
average structure Janner and Janssen (1977); Levine and Steinhardt (1986);
Socolar and Steinhardt (1986); Wolny (1998a); Xu and Mai (1998); Wolny
(1998b); Steurer and Haibach (1999); Aragón  _et al._ (2002). Since the band
inversion happens around $\Gamma$ (see Fig. S1 in Supplemental Material 111See
Supplemental Material at http://link.aps.org/supplemental/xxx, for more
details about the derivation of the effective Hamiltonian of the quasicrystal
and numerical results, which include Refs. Jiang and Zhang (2014); Steurer and
Deloudi (2009); Hiller (1985); Rabson  _et al._ (1991); Mermin (1992); Dräger
and Mermin (1996); Slater and Koster (1954); Janner and Janssen (1977); Levine
and Steinhardt (1986); Socolar and Steinhardt (1986); Wolny (1998a); Xu and
Mai (1998); Wolny (1998b); Steurer and Haibach (1999); Aragón  _et al._
(2002); König  _et al._ (2008); Qi and Zhang (2011); Ezawa (2018b); Hou  _et
al._ (2007); Teo and Kane (2010); Wang and Moore (2019); Huang and Liu (2018a,
b)), a low-energy effective model at the long-wavelength limit, which can be
approximated by the average structure of quasicrystals, is sufficient to
describe the relevant topological physics.

To derive the effective Hamiltonian in the pseudo $k$-space, we first apply
the JZ Fourier expansion to Eq. (Effective Model for Fractional Topological
Corner Modes in Quasicrystals) without the Zeeman field, which yields (see
Supplemental Material ††footnotemark: ),

$\displaystyle H(\mathbf{\Pi})$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\epsilon_{\alpha}c_{\mathbf{\Pi},\alpha}^{\dagger}c_{\mathbf{\Pi},\alpha}+\sum_{\alpha,\beta}t_{\alpha,\beta}(\mathbf{\Pi})c_{\mathbf{\Pi},\alpha}^{\dagger}c_{\mathbf{\Pi},\beta}$
(6) $\displaystyle+$ $\displaystyle
i\lambda\left(c_{\mathbf{\Pi},p_{y}}^{\dagger}s_{z}c_{\mathbf{\Pi},p_{x}}-c_{\mathbf{\Pi},p_{x}}^{\dagger}s_{z}c_{\mathbf{\Pi},p_{y}}\right).$

Here we have adopted the long-wavelength approximation to calculate an average
hopping as

$t_{\alpha,\beta}\left(\mathbf{\Pi}\right)\approx\lim_{V\rightarrow\infty}\frac{1}{V}\int
d\mathbf{r}\mathcal{P}(\mathbf{r})t_{\alpha,\beta}\left(\mathbf{r}\right)e^{i\left[(\mathcal{S}\cdot\mathbf{\Pi})^{T}\cdot\mathbf{r}\right]}.$
(7)

Note that around $\Gamma$, electron scattering is not affected by the local
details of the quasicrystalline potential but only feels an average effective
potential Aragón  _et al._ (2002); Wolny and Kozakowski (2003). Here
$\mathcal{P}(\mathbf{r})$ is the statistical average distribution of
interatomic vectors in quasicrystals, which is also known as the Patterson
function and can be extracted from diffraction data Patterson (1934); van
Smaalen (1989). We then downfold the Hamiltonian to the two-orbital subspace
around the Fermi level based on the Löwding perturbation method Winkler
(2003); Huang  _et al._ (2013), followed by projecting the two-orbital
Hamiltonian to 2D space by taking $\mathbf{k}=\mathcal{S}\cdot\mathbf{\Pi}$
and expanding $\mathbf{k}$ around $\Gamma$. Finally, we obtain the effective
$k\cdot p$ Hamiltonian without the Zeeman term as ††footnotemark: ,

$\displaystyle H_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\left(m-bk^{2}\right)\sigma_{z}+[ak_{y}+a^{\prime}(k_{y}^{3}+k_{x}^{2}k_{y})]\sigma_{x}s_{z}$
(8) $\displaystyle+$
$\displaystyle[ak_{x}+a^{\prime}(k_{x}^{3}+k_{x}k_{y}^{2})]\sigma_{y},$

where $\sigma$ and $s$ are Pauli matrices acting on the orbital and spin
degrees of freedom, respectively. It is noteworthy that Eq. (8) still
satisfies the $C_{5v}$ symmetry of the pentagonal quasicrystal. More
importantly, its solution of band structure shows a band inversion between s
and p orbitals around $\Gamma$ [see Fig. 2(b)], indicating a nontrivial
topology.

Now, let us consider a generic open boundary along
$\vec{e}_{1}=\cos\zeta\hat{e}_{x}+\sin\zeta\hat{e}_{y}$ direction with the
normal vector $\vec{e}_{2}=-\sin\zeta\hat{e}_{x}+\cos\zeta\hat{e}_{y}$ [see
Fig. 2(c)]. To further derive the low-energy Hamiltonian for the edge states
along $\vec{e}_{1}$, we perform a rotational transformation to Eq. (8), and
replace $k_{2}\rightarrow-i\partial_{x_{2}},k_{1}\rightarrow 0$. After some
algebra (see Supplemental Material ††footnotemark: ), we obtain a pair of
edge-state solutions for two spin channels, and arrive at the 1D edge model to
the leading order in $k_{1}$: $\mathcal{H}_{\mathrm{edge}}=-ak_{1}s_{z}$. It
indicates the existence of a pair of spin-polarized gapless edge states
protected by time-reversal symmetry.

When an in-plane Zeeman field of
$H_{\mathrm{in}}^{\prime}=\delta_{M}(\cos\theta\sigma_{0}s_{x}+\sin\theta\sigma_{0}s_{y})$
is applied, the edge state $\mathcal{H}_{\mathrm{edge}}$ will generally be
gapped by a mass term

$\mathcal{M}(\zeta_{i})\sim\cos\phi_{i}s_{x}+\sin\phi_{i}s_{y},$ (9)

where $\phi_{i}=\zeta_{i}+\theta-\pi/2$ is the generalized phase of the
effective Dirac mass Hou  _et al._ (2007); Teo and Kane (2010), which depends
on the orientation angle $\zeta_{i}$ of the i-th edge. From the work of Jackiw
and Rebbi Jackiw and Rebbi (1976), a phase shift of
$\Delta\phi=\Delta\zeta=\pi$ between two edges results in a CM with fractional
charge $Q=e/2$ due to mass inversion. Moreover, according to Moore’s theory
Wang and Moore (2019), a kink arising from the effective mass term at the
corner gives rise to a localized CM with fractional charge of
$Q=e|\Delta\phi/2\pi|=e|\Delta\zeta/2\pi|$ Goldstone and Wilczek (1981);
Jackiw and Semenoff (1983). Remarkably, for adjacent edges of a pentagonal
Penrose-tiling quasicrystal, the angle difference $\Delta\zeta=2\pi/5$ leads
to a fractional charge $Q=e/5$. On the contrary, an out-of-plane Zeeman field
$H_{\mathrm{out}}^{\prime}=\delta_{M}\sigma_{0}s_{z}$ only contributes an
energy-shift term of $\delta_{M}s_{z}$ which cannot gap the edge state
$\mathcal{H}_{\mathrm{edge}}$ (see Supplemental Material ††footnotemark: ).

For the Ammann-Beenker-tiling quasicrystal, a similar low-energy theory can be
derived and the phase shift becomes $\Delta\phi=\pi/4$ at corners of the
octagonal sample, giving rise to a fractional charge of $Q=e/8$ at each
corner. Thus, our model is valid for quasicrystals with odd- as well as even-
rotational symmetries. We stress that our approach of realizing CM in
quasicrystals via the fractional mass kink is fundamentally different from
previous works based on the mass-inversion mechanism Chen  _et al._ (2020b);
Varjas  _et al._ (2019); Hua  _et al._ (2020), because their CMs rely on an
alternating sign of the artificial mass term at the boundary, and hence does
not work for odd rotations.

Figure 3:  Energy spectrum of a finite (a) pentagonal Penrose-tiling
quasicrystal sample with 2146 atoms and (b) octagonal Ammann-Beenker-tiling
quasicrystal sample with 2241 atoms in the presence of an in-plane Zeeman
field with $\delta_{M}=0.12$ eV and $\delta_{M}=0.15$ eV, respectively. The
parameters used here are
$\epsilon_{s}=0.7,\epsilon_{p}=-2.3,V_{ss\sigma}=-V_{sp\sigma}=-0.17,V_{pp\sigma}=V_{pp\pi}=0.34$,
and $\lambda=1.0$ eV. Insets show the spatial distributions of CMs.

Topological CMs in quasicrystals.— The above effective model predictions are
confirmed by numerical TB calculations. In the absence of the Zeeman field,
the model (Effective Model for Fractional Topological Corner Modes in
Quasicrystals) describes a TI state in the pentagonal quasicrystal, which is
verified by the calculation results of a nonzero spin Bott index ($B_{s}=1$)
Huang and Liu (2018a, b) and time-reversal symmetry-protected edge states
residing inside the bulk gap (see Fig. S3 and S4 in Supplemental Material
††footnotemark: ). In the presence of an in-plane Zeeman field along the
x-axis with $\delta_{M}=0.12$ eV, the energy spectrum of the finite pentagonal
quasicrystal is gapped and five states appear at the Fermi level and separate
from other states, as shown in Fig. 3(a). We plot the spatial distribution of
these states [see inset of Fig. 3(a)], and found that they are indeed CMs
localized at corners of the pentagonal quasicrystal. This implies that the
system becomes a HOTI, although the CMs are not located at the mid-gap
position due to the lack of chiral or particle-hole symmetry. Remarkably,
since four of the five CMs are occupied for a charge-neutral system, a
fractional charge of $e/5$ localized at each corner can be realized by adding
one electron, resulting in a fractionalized charge distribution due to the
filling anomaly Benalcazar  _et al._ (2019); Takahashi  _et al._ (2021).
Similarly, for the octagonal Ammann-Beenker-tiling quasicrystal, we found
eight CMs in the gap and a charge fractionalization of $e/8$ per corner if one
extra electron is added [see Fig. 3 (b)]. It is worth noting that the corner
charge in crystals are fractionally quantized module of $e/n$ with $n=2,3,4$
and $6$, covering all the allowed rotational symmetries by the
crystallographic restriction theorem. Our results extend the possible values
of the corner charge fractionalization with $e/q$, where $q=5,8$ and other
rotational orders in quasicrystals.

We further investigated numerically the phase evolution with the field
strength $\delta_{M}$ for the Penrose-tiling quasicrystal. By increasing
$\delta_{M}$, the bulk energy gap, which is estimated from the calculations of
quasicrystal approximants with 1364 atoms ††footnotemark: , decreases
gradually and closes eventually when $\delta_{M}>0.2$ eV. Whilst the CMs exist
in the gap only when $\delta_{M}<0.15$ eV, beyond which they merge into the
bulk spectrum. In addition, we studied the effect of arbitrary Zeeman field
orientations in the x-y plane, e.g., $\bm{m}=(\cos\theta,\sin\theta,0)$, and
found that the CMs persist regardless of $\theta$. This can be understood by
simply performing a rotation of spin about the z axis to make the Zeeman field
pointing along the $x$ direction Ren  _et al._ (2020).

Figure 4:  Majorana CMs in the quasicrystal TI/s-wave SC heterostructure with
an in-plane Zeeman field. (a) Quasiparticle energy spectrum vs the Zeeman-
field strength $\delta_{M}$ at fixed $\Delta_{0}=0.04$ eV and $\mu=0.02$ eV.
(b) Same as (a) vs the chemical potential $\mu$ at fixed $\Delta_{0}=0.04$ eV
and $\delta_{M}=0.13$ eV. (c) Energy spectrum of a finite pentagonal sample at
$(\delta_{M},\mu,\Delta_{0})=(0.13,0.02,0.04)$ eV. Other parameters are the
same as Fig. 3. The inset shows the real-space distribution of Majorana CMs.

Majorana CMs in quasicrystal TI and s-wave SC heterostructures.— In addition
to fractional topological corner charge, our model and approach are also
applicable to create Majorana CMs in the TI quasicrystals, under an in-plane
Zeeman field and in proximity with s-wave SCs. The physics of the
heterostructure can be described by an effective Bogoliubov-de Gennes (BdG)
Hamiltonian

$H_{\mathrm{BdG}}=\left(\begin{array}[]{cc}H-\mu&\Delta_{0}\\\
\Delta_{0}^{\dagger}&-H^{*}+\mu\\\ \end{array}\right),$ (10)

where $\mu$ is the chemical potential and $\Delta_{0}$ denotes s-wave SC
pairing gap by proximity effect. Although the normal-state part of
$H_{\mathrm{BdG}}$ is topologically nontrivial with helical edge states, the
proximity induced s-wave SC pairing necessarily trivializes the full BdG model
and gaps out the edge states Zhang  _et al._ (2019). As shown in Fig. 4(a),
with the increasing in-plane Zeeman field $\delta_{M}$ at a fixed
$\Delta_{0}$, the quasiparticle energy gap first closes and then reopens
accompanied by the emergence of localized modes at corners, implying the
nontrivial topology is resumed. Moreover, these CMs can be fine-tuned to zero
energy by adjusting the chemical potential, giving rise to Majorana CMs Yan
(2019) [see Fig. 4(b)]. As shown in Fig. 4(c), in the topological phase at
$(\delta_{M},\mu,\Delta_{0})=(0.13,0.02,0.04)$ eV, five pairs of Majorana zero
modes emerge inside the gap. The inset of Fig. 4(c) shows the spatial
distribution of these zero modes, confirming that they are localized around
the corners of a finite pentagonal quasicrystal.

Conclusion.— We have devised a low-energy theory of quasicrystals under the
long-wavelength approximation, and demonstrated that higher-order topological
CMs with fractional charges can be generated by in-plane Zeeman-field-induced
fractional mass kinks. Our model predictions are further confirmed by
numerical TB calculations, which show also emergence of Majorana CMs in TI
quasicrystals in proximity with an s-wave SC. Our work greatly extends the
higher-order topological physics for mass-kink induced domain-wall states to
quasicrystals and establishes a generic theoretical framework to study the
low-energy physics of quasicrystalline systems.

###### Acknowledgements.

This work was supported by the National Key R&D Program of China (No.
2021YFA1401600), the National Natural Science Foundation of China (Grant No.
12074006), and the start-up fund from Peking University. F.L. acknowledges the
support of U.S. DOE-BES (Grant No. DE-FG02-04ER46148). The computational
resources were supported by the high-performance computing platform of Peking
University.

## References

  * Benalcazar  _et al._ (2017) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017).
  * Schindler  _et al._ (2018a) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Sci. Adv. 4, eaat0346 (2018a).
  * Langbehn  _et al._ (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * Song  _et al._ (2017) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * Schindler  _et al._ (2018b) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Nat. Phys. 14, 918 (2018b).
  * Yue  _et al._ (2019) C. Yue, Y. Xu, Z. Song, H. Weng, Y.-M. Lu, C. Fang, and X. Dai, Nat. Phys. 15, 577 (2019).
  * Xu  _et al._ (2019) Y. Xu, Z. Song, Z. Wang, H. Weng, and X. Dai, Phys. Rev. Lett. 122, 256402 (2019).
  * Zhang  _et al._ (2020) R.-X. Zhang, F. Wu, and S. Das Sarma, Phys. Rev. Lett. 124, 136407 (2020).
  * Sheng  _et al._ (2019) X.-L. Sheng, C. Chen, H. Liu, Z. Chen, Z.-M. Yu, Y. X. Zhao, and S. A. Yang, Phys. Rev. Lett. 123, 256402 (2019).
  * Lee  _et al._ (2020) E. Lee, R. Kim, J. Ahn, and B.-J. Yang, npj Quantum Mater. 5, 1 (2020).
  * Liu  _et al._ (2019a) B. Liu, G. Zhao, Z. Liu, and Z. Wang, Nano Lett. 19, 6492 (2019a).
  * Park  _et al._ (2019) M. J. Park, Y. Kim, G. Y. Cho, and S. Lee, Phys. Rev. Lett. 123, 216803 (2019).
  * Liu  _et al._ (2021) B. Liu, L. Xian, H. Mu, G. Zhao, Z. Liu, A. Rubio, and Z. F. Wang, Phys. Rev. Lett. 126, 066401 (2021).
  * Huang and Liu (2021) H. Huang and F. Liu, Natl. Sci. Rev. 2021, nwab170 (2021).
  * Pan  _et al._ (2022) M. Pan, D. Li, J. Fan, and H. Huang, npj Comput. Mater. 8, 1 (2022).
  * Ezawa (2018a) M. Ezawa, Phys. Rev. Lett. 120, 026801 (2018a).
  * Liu  _et al._ (2019b) F. Liu, H.-Y. Deng, and K. Wakabayashi, Phys. Rev. Lett. 122, 086804 (2019b).
  * Wang  _et al._ (2019) Z. Wang, B. J. Wieder, J. Li, B. Yan, and B. A. Bernevig, Phys. Rev. Lett. 123, 186401 (2019).
  * Hsu  _et al._ (2019) C.-H. Hsu, X. Zhou, T.-R. Chang, Q. Ma, N. Gedik, A. Bansil, S.-Y. Xu, H. Lin, and L. Fu, Proc. Natl. Acad. Sci. 116, 13255 (2019).
  * Ren  _et al._ (2020) Y. Ren, Z. Qiao, and Q. Niu, Phys. Rev. Lett. 124, 166804 (2020).
  * Chen  _et al._ (2020a) C. Chen, Z. Song, J.-Z. Zhao, Z. Chen, Z.-M. Yu, X.-L. Sheng, and S. A. Yang, Phys. Rev. Lett. 125, 056402 (2020a).
  * Zhao  _et al._ (2021) G. Zhao, H. Mu, H. Zhang, and Z. Wang, arXiv preprint arXiv:2107.05910 (2021).
  * Fu (2011) L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Fu and Kane (2007) L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007).
  * Jackiw and Rebbi (1976) R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * Geier  _et al._ (2018) M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Phys. Rev. B 97, 205135 (2018).
  * Trifunovic and Brouwer (2019) L. Trifunovic and P. W. Brouwer, Phys. Rev. X 9, 011012 (2019).
  * Trifunovic and Brouwer (2020) L. Trifunovic and P. W. Brouwer, Phys. Status Solidi B 258, 2000090 (2020).
  * Varjas  _et al._ (2019) D. Varjas, A. Lau, K. Pöyhönen, A. R. Akhmerov, D. I. Pikulin, and I. C. Fulga, Phys. Rev. Lett. 123, 196401 (2019).
  * Chen  _et al._ (2020b) R. Chen, C.-Z. Chen, J.-H. Gao, B. Zhou, and D.-H. Xu, Phys. Rev. Lett. 124, 036803 (2020b).
  * Hua  _et al._ (2020) C.-B. Hua, R. Chen, B. Zhou, and D.-H. Xu, Phys. Rev. B 102, 241102 (2020).
  * Huang  _et al._ (2021) H. Huang, J. Fan, D. Li, and F. Liu, Nano Lett. 21, 7056 (2021).
  * Fan and Huang (2022) J. Fan and H. Huang, Front. Phys. 17, 13203 (2022).
  * Wang and Moore (2019) Y.-Q. Wang and J. E. Moore, Phys. Rev. B 99, 155102 (2019).
  * Goldstone and Wilczek (1981) J. Goldstone and F. Wilczek, Phys. Rev. Lett. 47, 986 (1981).
  * Jackiw and Semenoff (1983) R. Jackiw and G. Semenoff, Phys. Rev. Lett. 50, 439 (1983).
  * Jiang and Zhang (2014) K. Jiang and P. Zhang, J. Comput. Phys. 256, 428 (2014).
  * Penrose (1974) R. Penrose, Bull. Inst. Math. Appl. 10, 266 (1974).
  * Beenker (1982) F. P. M. Beenker,  _Algebraic theory of non-periodic tilings of the plane by two simple building blocks: a square and a rhombus_ (Eindhoven University of Technology Eindhoven, The Netherlands, 1982).
  * Huang and Liu (2018a) H. Huang and F. Liu, Phys. Rev. Lett. 121, 126401 (2018a).
  * Huang and Liu (2018b) H. Huang and F. Liu, Phys. Rev. B 98, 125130 (2018b).
  * Huang and Liu (2019) H. Huang and F. Liu, Phys. Rev. B 100, 085119 (2019).
  * Janssen  _et al._ (2007) T. Janssen, D. Janssen, G. Chapuis, and M. de Boissieu,  _Aperiodic Crystals: From Modulated Phases to Quasicrystals_ (Oxford University Press, Oxford, 2007).
  * Bellissard (1992) J. Bellissard, in  _From number theory to physics_ (Springer, Berlin, Heidelberg, 1992) pp. 538–630.
  * Gambaudo and Vignolo (2014) J.-M. Gambaudo and P. Vignolo, New J. Phys. 16, 043013 (2014).
  * Steurer and Deloudi (2009) W. Steurer and S. Deloudi,  _Crystallography of Quasicrystals: Concepts, Methods and Structures_ , Vol. 126 (Springer Series in Materials Science, 2009).
  * Janner and Janssen (1977) A. Janner and T. Janssen, Phys. Rev. B 15, 643 (1977).
  * Levine and Steinhardt (1986) D. Levine and P. J. Steinhardt, Phys. Rev. B 34, 596 (1986).
  * Socolar and Steinhardt (1986) J. E. S. Socolar and P. J. Steinhardt, Phys. Rev. B 34, 617 (1986).
  * Wolny (1998a) J. Wolny, Philos. Mag. A 77, 395 (1998a).
  * Xu and Mai (1998) L. Xu and Z.-H. Mai, Phys. Rev. B 58, 23 (1998).
  * Wolny (1998b) J. Wolny, Acta Crystallogr. A 54, 1014 (1998b).
  * Steurer and Haibach (1999) W. Steurer and T. Haibach, Acta Crystallogr. A 55, 48 (1999).
  * Aragón  _et al._ (2002) J. L. Aragón, G. G. Naumis, and M. Torres, Acta Crystallogr. A 58, 352 (2002).
  * Note (1) See Supplemental Material at http://link.aps.org/supplemental/xxx, for more details about the derivation of the effective Hamiltonian of the quasicrystal and numerical results, which include Refs. Jiang and Zhang (2014); Steurer and Deloudi (2009); Hiller (1985); Rabson  _et al._ (1991); Mermin (1992); Dräger and Mermin (1996); Slater and Koster (1954); Janner and Janssen (1977); Levine and Steinhardt (1986); Socolar and Steinhardt (1986); Wolny (1998a); Xu and Mai (1998); Wolny (1998b); Steurer and Haibach (1999); Aragón  _et al._ (2002); König  _et al._ (2008); Qi and Zhang (2011); Ezawa (2018b); Hou  _et al._ (2007); Teo and Kane (2010); Wang and Moore (2019); Huang and Liu (2018a, b).
  * Wolny and Kozakowski (2003) J. Wolny and B. Kozakowski, Acta Crystallogr. A 59, 54 (2003).
  * Patterson (1934) A. L. Patterson, Phys. Rev. 46, 372 (1934).
  * van Smaalen (1989) S. van Smaalen, Phys. Rev. B 39, 5850 (1989).
  * Winkler (2003) R. Winkler,  _Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems_ , Vol. 191 (Springer Nature, Switzerland AG, 2003).
  * Huang  _et al._ (2013) H. Huang, W. Duan, and Z. Liu, New J. Phys. 15, 023004 (2013).
  * Hou  _et al._ (2007) C.-Y. Hou, C. Chamon, and C. Mudry, Phys. Rev. Lett. 98, 186809 (2007).
  * Teo and Kane (2010) J. C. Y. Teo and C. L. Kane, Phys. Rev. B 82, 115120 (2010).
  * Benalcazar  _et al._ (2019) W. A. Benalcazar, T. Li, and T. L. Hughes, Phys. Rev. B 99, 245151 (2019).
  * Takahashi  _et al._ (2021) R. Takahashi, T. Zhang, and S. Murakami, Phys. Rev. B 103, 205123 (2021).
  * Zhang  _et al._ (2019) R.-X. Zhang, W. S. Cole, X. Wu, and S. Das Sarma, Phys. Rev. Lett. 123, 167001 (2019).
  * Yan (2019) Z. Yan, Phys. Rev. B 100, 205406 (2019).
  * Hiller (1985) H. Hiller, Acta Crystallogr. A 41, 541 (1985).
  * Rabson  _et al._ (1991) D. A. Rabson, N. D. Mermin, D. S. Rokhsar, and D. C. Wright, Rev. Mod. Phys. 63, 699 (1991).
  * Mermin (1992) N. D. Mermin, Rev. Mod. Phys. 64, 3 (1992).
  * Dräger and Mermin (1996) J. Dräger and N. D. Mermin, Phys. Rev. Lett. 76, 1489 (1996).
  * Slater and Koster (1954) J. C. Slater and G. F. Koster, Phys. Rev. 94, 1498 (1954).
  * König  _et al._ (2008) M. König, H. Buhmann, L. W. Molenkamp, T. Hughes, C.-X. Liu, X.-L. Qi, and S.-C. Zhang, J. Phys. Soc. Japan 77, 031007 (2008).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Ezawa (2018b) M. Ezawa, Phys. Rev. Lett. 121, 116801 (2018b).