# On Exact and Efficient Inference for Many Normal Means

Chuanhai Liu

Department of Statistics, Purdue University

###### Abstract

Inference about the unknown means
$\theta=(\theta_{1},...,\theta_{n})^{\prime}\in\mathbb{R}^{n}$ in the sampling
model $X\sim N_{n}(\theta,I)$ from the observed $X$, known as the many-normal-
means problem, has proven to be fundamental and yet challenging inasmuch as a
satisfactory understanding remains elusive. To tackle exact and efficient
inference about $\theta$, in this paper we propose innovative formulations of
Inferential Models for two kinds of this problem: the classic kind given as is
and the empirical-Bayes kind where $\theta$’s are further assumed to be an
unobservable sample from an unknown non-parametric distribution $\mathbb{G}$.
The formulation for the empirical-Bayes kind via numerical deconvolution
allows for prior-free probabilistic inference with over-parameterization for
the non-parametric model $\mathbb{G}$, whereas the formation for the first
kind utilizes a latent random permutation and as a result provides a sound
reasoning with uncertainty toward a deeper understanding. For uncertainty
quantification with the more familiar frequentist inference framework, the
method of maximum plausibility estimation is used for point estimation. Exact
but conservative interval estimation is obtained based on plausibility, with a
Monte Carlo based adaptive-adjustment approach to constructing shorter
confidence intervals with targeted coverage. These methods are illustrated via
simulation studies and a real-data example. The numerical results show that
for interval estimation, adaptive intervals are satisfactory in both coverage
and efficiency and that for point estimation, the proposed methods outperform
the traditional James-Stein and Efron’s g-modeling in terms of mean square
error. The paper concludes with a few remarks, including future developments
and extensions of the proposed methods.

Key Words: Empirical Bayes, Inferential Models, Latent permutation, Numerical
deconvolution, Stein’s paradox.

##  1 Introduction

The many-normal-means problem concerns simultaneous inference about a
collection of unknown normal means. More precisely, one has the observed data
$X_{1},...,X_{n}$ that are assumed to have the sampling model

$X_{i}=\theta_{i}+Z_{i}\qquad(Z_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,1);\;i=1,...,n)$ (1.1)

This problem has proven to be fundamental in statistics. The most noticeable
contribution, perhaps, dates back to Stein (1956), on the remarkable discovery
that the usual maximum likelihood estimates $\hat{\theta}_{i}=X_{i}$ are
inadmissible under quadratic loss for $n\geq 3$. Understanding of and
developing methods for this problem have been resulting in successful research
in multi-parameter or high-dimensional problems, including, for example, the
James-Stein shrinkage method (James and Stein, 1961), the empirical Bayes
method of Efron and Morris (1972), the wavelet shrinkage method of Donoho and
Johnstone (1995) for non-parametric regression, and LASSO of Tibshirani (1996)
for linear regression. The literature on the many normal means problem itself
and methods built up on shrinkage is massive and is, thereby, impossible to
provide a complete review in this note.

We refer to the problem setting given above as is as the first or classic
kind, reflecting Stein’s discovery. It is perhaps somewhat surprising that the
understanding of the Stein phenomenon has not been very satisfactory, at least
not to a level that would help develop efficient methods based on a sound
reasoning with uncertainty, as Stigler (1990) wrote:

> The best heuristic explanation that has been offered is a Bayesian argument
> … Another explanation that has been offered is that [the James-Stein
> estimator] can be viewed as a relative of a “pre-test” estimator … But
> neither of these explanations is fully satisfactory (although both help
> render the result more plausible); the first because it requires special a
> priori assumptions where Stein did not, the second because it corresponds to
> the result only in the loosest qualitative way. The difficulty of
> understanding the Stein paradox is compounded by the fact that its proof
> usually depends on explicit computation of the risk function or the theory
> of complete sufficient statistics, by a process that convinces us of its
> truth without really illuminating the reasons that it works.

Attributing back to Francis Galton’s knowledge in the 1880’s and viewing Stein
estimation as an inverse regression problem $E(\theta|X)$, Stigler (1990)
provided a nice interpretation of the James-Stein (shrinking toward zero) and
Efron-Morris (shrinking toward the sample mean $\bar{X}$) estimators of
$\theta$. The fundamental problem of this classic kind has however remained
unsolved, due to the difficulty as recognized by Stigler (1990):

> With no distributional assumptions about the $\theta$’s, we are of course
> prevented from looking for an optimal estimate of “$E(\theta|X)$”.

Our innovative formulation in Section 4 via introducing a latent permutation
in addition to unobserved $Z$’s in (1.1) leads to an inferentially valid
$E(\theta|X)$ and thus offers a new understanding, shedding light on this
challenging problem from an Inferential Models perspective.

A popular variant of the classic kind is obtained by modeling the unknown
$\theta_{i}$’s as an unobserved sample from a non-parametric distribution
$\mathbb{G}$ with the cumulative distribution function (cdf) $G(.)$. In this
case, integrating out the unobservable $\theta_{i}$’s leads to inference about
the mixture distribution

$F_{G}(x)=\int\phi(x-\theta)dG(\theta)$ (1.2)

from the observed $X$’s. Note that in this case, $X$’s is necessarily a sample
from the convoluted distribution $F_{G}(x)$. Inference about $\theta$’s is
then carried out with their conditional distribution, given $X_{i}$’s and
estimates of $G$. The central problem in this case is thus to estimate $G(.)$.
There is a considerable literature on this deconvolution problem, as in Laird
(1978), Carroll and Hall (1988), Stefanski and Carroll (1990), Fan (1991),
Zhang (1997), Hall and Meister (2007), Butucea and Comte (2009), and Efron
(2014, 2016). For convenience, we refer to inference via this non-parametric
partial Bayes model as the second or (non-parametric) empirical-Bayes kind, in
contrast to the classic kind.

Here, we consider exact and efficient inference for the two kinds of the many-
normal-means problems by taking the Inferential Models (IMs) approach, which
is to be briefly reviewed in Section 2. The empirical-Bayes kind and classic
kind are investigated in Sections 3 and 4, respectively. The development of
IMs for the empirical-Bayes kind appears relatively straightforward because of
the nice independence structure imposed. Nevertheless, following the
discretization idea in Efron (2014, 2016) but not going further for
$g$-modeling, our approach is to make inference via over-parameterization or
numerical deconvolution without effectively having to impose more restrictions
on $\mathbb{G}$. It should be noted that unlike traditional Frequentist and
Bayesian methods, IMs can provides honest inference even if the unknown
parameters for $\mathbb{G}$ are unidentifiable.

For the classic kind, we introduce a sorted standard normal sample with a
random permutation as the so-called auxiliary variable. This novel approach,
with a full IMs argument, leads to a prior-free probabilistic inference, for
example, by making use of the quantity $E(\theta|X)$ for point estimation. As
a result, it sheds a new light on a satisfactory understanding Stein
estimation.

Predictive random sets for IMs inference are constructed to produce
probabilistic inference in terms of belief and plausibility. It is well known
that valid probabilistic inference says more than but can be used for point
and interval estimations. As a method for point estimation, maximum
plausibility estimation is discussed in detail. For interval estimation in the
more familiar frequentist inference framework, exact but conservative interval
estimation is obtained based on plausibility. A Monte Carlo approach is
proposed as an adaptive adjustment method to create shorter confidence
intervals with targeted coverage. This adaptive method is so straightforward
in the context of frequency calibration that it is simply presented as a
technical remark; see Remark 3.4. This is illustrated in Section 5.2 for a
real-data example.

To evaluate the performance of the proposed methods, a simulation study is
conducted in Section 5 with the computational methods provided in Appendices.
The numerical results show that the proposed methods outperform the
traditional James-Stein and Efron’s g-modeling (Narasimhan and Efron, 2020),
and are comparable to the auto-modeling method proposed recently by Jiang and
Liu (2022). Section 6 concludes the paper with a few remarks, including a
brief discussion of other kinds of the many-normal-means problem, future
potential developments, and extensions.

##  2 Inferential models: a brief review

IMs was proposed in Martin and Liu (2013) as an inferential framework for
prior-free probabilistic inference. It makes use of the so-called unobserved
but predictable auxiliary random variables, or simply, auxiliary variables $U$
that has a known distribution ${\mathbf{P}}_{U}$ on ${\mathbb{U}}$. The setup
for IMs is a system of equations, called association, representing a
statistical model with unknown parameter
$\theta\in\Theta\subseteq{\mathbb{R}}^{p}$ for the observed data
$X\in{\mathbb{X}}$ via an auxiliary variables $U$. For the sake of clarity, we
write this association as

$X=a(U,\theta)\qquad(X\in{\mathbb{X}},\theta\in\Theta,U\sim{\mathbf{P}}_{U})$
(2.1)

IMs relies on so-called valid or frequency-calibrated predictive random sets
(PRS) $\mathcal{S}$ for predicting the unobserved $U$ to produce valid
uncertainty quantification on unknown quantities of interest. Formally, a PRS
$\mathcal{S}$ is said to be valid for $U$ if it is independent of $U$ and

$\mbox{Prob}(\mbox{Prob}(\mathcal{S}\not\ni U)\geq 1-\alpha)\leq\alpha$ (2.2)

for all $\alpha\in(0,1)$. Typically, $\mathcal{S}$ is chosen in such a way
that $\mbox{Prob}(\mathcal{S}\not\ni U)\sim\mbox{Uniform}(0,1)$, the standard
uniform distribution. For example, the $\mathcal{S}(U)=\\{u:|u|\leq|U|\\}$,
$U\sim\mbox{N}(0,1)$, is valid for predicting an unobserved realization from
$\mbox{N}(0,1)$.

Given a PRS $\mathcal{S}$ for predicting $U$ in (2.1), the inverse set-valued
mapping

$\Theta_{X}(\mathcal{S})=\\{\theta:X=a(u,\theta)\,\mbox{ for some
}\,u\in\mathcal{S}\\}$ (2.3)

defines a random set on $\Theta$. For simplicity, here we assume
$\Theta_{X}(\mathcal{S})\neq\emptyset$ for all values of $\mathcal{S}$ (c.f.
Leaf and Liu (2012) for the case where
$\mbox{Prob}(\mathcal{S}=\emptyset)>0$). IMs produces uncertainty assessments
on assertions of interest, $A\subset\Theta$, by computing the following lower
and upper probabilities

$\mbox{bel}_{X}(A)=\mbox{Prob}(\Theta_{X}(\mathcal{S})\subseteq A)$ (2.4)

and

$\mbox{pl}_{X}(A)=1-\mbox{bel}_{X}(A^{c})=\mbox{Prob}(\Theta_{X}(\mathcal{S})\not\subseteq
A^{c})$ (2.5)

The probabilities (2.4) and (2.5) are also known as the belief function at $A$
and plausibility of $A$, respectively. A recent discussion of these IMs lower
and upper probabilities and other existing ones is given in Liu and Martin
(2021); see also references therein.

As described above, IMs is rooted in R. A. Fisher’s fiducial argument in the
sense to make prior-free probabilistic inference about parameters by using
auxiliary variable to represent the sample model. Avoiding the difficulties of
R. A. Fisher’s fiducial argument and the related theories (Zhang and Liu,
2011; Martin et al., 2010; Liu and Martin, 2015), however, IMs emphasizes
frequency-calibrated/valid and efficient probabilistic inference. Together
with Conditional IMs and marginal IMs (Martin and Liu (2015a) and Martin and
Liu (2015b)), subject to further developments, this theory, if so desired,
also provides as a new mathematical tool for developing exact and efficient
frequentist methods. For example, for the familiar Neyman-Pearson school of
thought, it offers its plausibility regions for unknown parameters as exact
Neyman-Pearson confidence regions (Martin and Liu, 2013), and for the
Fisherian school of thought, its plausibility for assertions of interest
provides a plausibility interpretation of or a probabilistic alternative to
Fisher’s P-value (Martin and Liu, 2014).

Nevertheless, research on foundations of statistical inference has been
challenging. Like other existing approaches other than the standard Bayes and
frequentist schools of thought, including generalized fiducial (Hannig, 2009;
Cui and Hannig, 2022), the method of confidence distribution (Xie et al.,
2011; Schweder and Hjort, 2016), and the Dempster-Shafer theory of belief
functions (Shafer, 1976, 1992; Dempster, 2008; Martin et al., 2010), IMs is
subject to scrutinization and acceptance. This is expected to take time.

##  3 IMs for MNMs of the empirical Bayes kind

###  3.1 The sampling model and IMs association

Without real loss of generality, our jumping-off point is to take a numerical
approach to approximate $G(.)$ in (1.2) via discretization; see also Efron
(2014, 2016). Denote by $\vartheta_{1}<...<\vartheta_{K}$ the known space of
$\theta$, and denote by $\gamma_{k}$ the unknown probability mass at
$\vartheta_{k}$, where

$\gamma\equiv(\gamma_{1},...,\gamma_{K})\in\left\\{(\gamma_{1},...,\gamma_{K}):\gamma_{k}\geq
0\mbox{ for all $k=1,...,K$; $\sum_{k=1}^{K}\gamma_{k}=1$}\right\\}$ (3.1)

######  Remark 3.1.

Take $\vartheta_{k}-\vartheta_{k-1}=\delta$ for all $k=1,...,K$. This discrete
approximation can also be view as numerical approach to mixed distributions,
including both discrete and continuous distributions.

In this scenario, $G(.)$ is represented by the $K$-dimensional $\gamma$ vector
in (3.1), that is, in terms of cdf

$G(\theta)=\sum_{\vartheta_{k}\leq\theta}\gamma_{k}$

and $X_{1},...,X_{i}$ is a sample from the convolution of $G(.)$ and $N(0,1)$:

$F_{\vartheta,\gamma}(x)=\sum_{k=1}^{K}\gamma_{k}\Phi(x-\vartheta_{k})$ (3.2)

where $\Phi(.)$ denotes the cdf of the standard norm distribution $N(0,1)$.
Furthermore, given $\gamma$ and $X$, $\theta_{i}$’s are independent of each
other with

$\mbox{Pr}(\theta_{i}=\vartheta_{j}|\gamma,X)=\frac{\gamma_{k}\phi(X_{i}-\theta_{j})}{\sum_{k=1}^{K}\gamma_{k}\phi(X_{i}-\vartheta_{k})}\qquad(j=1,...,K)$
(3.3)

where $\phi(.)$ stands for the probability density function of the standard
norm distribution $N(0,1)$.

From (3.2) and (3.3), we see that inference about $\theta_{i}$s is reduced to
inference about $\gamma_{k}$’s. It is known that
$F_{\vartheta,\gamma}(X_{i})$’s form a sample of size $n$ from
$F_{\vartheta,\gamma}(x)$. For an IMs approach, we take a sorted uniforms of
size of $n$, denoted by $U_{(1)}<...,<U_{(n)}$, as auxiliary variables and
specify the association as

$U_{(i)}=F_{\vartheta,\gamma}(X_{(i)})=\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(i)}-\vartheta_{k})$
(3.4)

###  3.2 A PRS for sorted uniforms

Let $\mathbb{S}_{n}$ denote the sample space of $U_{(i)}$’s, the simplex space
of size $n$:

$\mathbb{S}_{n}=\\{(u_{1},...,u_{n}):0<u_{1}<...<u_{n}<1\\}$

Intuitively, we aim to constructing a PRS for $U_{(i)}$’s that should work
well for the well-known non-parametric one-sample test, where it is known that
the cdf values at the observations serve as the sufficient statistics and form
a sample of uniforms. For simplicity, we use a type of predictive random set
similar to that of Liu and Xie (2014) and Zhang (2010) for
$(\Phi(Z_{(1)}),...,\Phi(Z_{(n)}))$,

$\mathcal{U}=\\{u:u\in\mathbb{S}_{n},B(u)\leq B(U)\\}$ (3.5)

where $U=(U_{(1)},...,U_{(n)})$, the vector of sorted sample of size $n$ from
$\mbox{Uniform}(0,1)$ and

$B(u)=-\sum_{i=1}^{n}\left[a_{i}\ln(u_{i})+b_{i}\ln(1-u_{i})\right]$ (3.6)

with $a_{i}=\frac{1}{n-i+c_{n}}$, $b_{i}=\frac{1}{i-1+c_{n}}$, and
$c_{n}\approx\frac{2}{3}$. Zhang (2010) considered this PRS with $c_{n}=0.7$
as an IMs approach to the non-parametric one-sample problem and showed that
the resulting method is either comparable to or better than traditional
methods such as the Kolmogorov-Smirnov test.

Recall the familiar results on sorted uniforms regarding the marginal
distribution of $U_{(i)}$. We have $U_{(i)}\sim\mbox{Beta}(i,n-i+1)$ which has
its mean and variance given by $\frac{i}{n+1}$ and
$\frac{i(n-i+1)}{(n+1)^{2}(n+2)}$. The PRS defined by (3.5) and (3.6) is
constructed as a weighted logarithms of the marginal densities with a minor
modification so that the medians are approximately most likely (c.f.  Kerman,
2011, for the choice of $c_{n}=\frac{2}{3}$). The weights are set to be about
inversely proportional to their variances and, thereby, pay more attentions to
the tail areas of the underlying sampling distribution. This leads to a
slightly more general class of PRSs is to use alternative weight schemes:

$a_{i}=w_{i}(i-1+c_{n}),\;b_{i}=w_{i}(n-i+c_{n})$ (3.7)

with

$w_{i}=\frac{{[i(n-i+1)]^{-\nu/2}}}{\sum_{j=1}^{n}{[j(n-j+1)]^{-\nu/2}}}$
(3.8)

for $i=1,...,n$, where $\nu\in[0,2]$. The PRS (3.7) with $\nu=2$ is used in
the numeric examples in Section 5.

###  3.3 Induced random sets in the $\gamma$ space

A realization of random set (3.5) induced the random sets in the $\gamma$
space through the convolution (3.9):

$\Gamma(U)=\left\\{(\gamma_{1},...,\gamma_{K}):\;\mbox{$B\left(\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(1)}-\vartheta_{k}),...,\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(n)}-\vartheta_{k})\right)\leq
B(U)$}\right\\}$ (3.9)

This random set for $\gamma$ defines the set of numerical approximations to
the unknown distribution $G(.)$. The usual IMs uncertainty quantification
about $G(.)$ can be carried out as prescribed in Section 2. In the context of
multiple testing, for example, inference may be made about fractional number
of outliers. For MNMs problem, inference about $\theta_{i}$’s is often of
interest and obtained by making use of (3.9) and (3.3). This is discussed in
detail below in next two subsections.

###  3.4 Maximum plausibility estimation

A straightforward application of the above IMs is to produce point estimation.
Given an estimate of $G$, the usual Bayes posterior mean

$\hat{\theta}_{i}(\gamma)=\frac{\sum_{k=1}^{K}\gamma_{k}\phi(X_{i}-\vartheta_{k})\vartheta_{k}}{\sum_{k=1}^{K}\gamma_{k}\phi(X_{i}-\vartheta_{k})}$
(3.10)

provides the most efficient estimator (Brown, 1971), in terms of mean square
error (MSE):

$\mbox{MSE}(\hat{\theta})=\frac{1}{n}\sum_{i=1}^{n}(\hat{\theta}_{i}-\theta_{i})^{2}$
(3.11)

Here we consider MPE of $\gamma$. The smallest PRS with the nested focal
elements, obtained with the method of (elastic) for validity, is given by

$\mathcal{U}_{\min}=\left\\{u:B(u)\leq B_{\min}\right\\}$ (3.12)

where

$B_{\min}=\min_{\gamma}\mbox{$B\left(\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(1)}-\vartheta_{k}),...,\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(n)}-\vartheta_{k})\right)$}$
(3.13)

This induces the maximum plausibility set for $\gamma$

$\Gamma_{\min}=\left\\{(\gamma_{1},...,\gamma_{K}):\;\mbox{$B\left(\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(1)}-\vartheta_{k}),...,\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(n)}-\vartheta_{k})\right)\leq
B_{\min}$}\right\\}$ (3.14)

It follows that the maximum plausibility is necessarily one and is achieved
for all $\gamma\in\Gamma_{\min}$. Since $\Gamma_{\min}$ given in (3.14) is a
set, we have a set of MPEs of $\theta_{i}$, one for each
$\gamma\in\Gamma_{\min}$.

For simple practical applications, we can summarize the above set of MPEs by
two extremes, the lower and upper posterior expectations

$\hat{\theta}_{i}^{(\mbox{lower})}=\min_{\gamma\in\Gamma_{\min}}\hat{\theta}_{i}(\gamma)\qquad\mbox{and}\qquad\hat{\theta}_{i}^{(\mbox{upper})}=\max_{\gamma\in\Gamma_{\min}}\hat{\theta}_{i}(\gamma)$
(3.15)

for $i=1,...,n.$ When a single value estimate is desired, the mean of the two
extremes may be used, that is,

$\hat{\theta}_{i}=\frac{\hat{\theta}_{i}^{(\mbox{lower})}+\hat{\theta}_{i}^{(\mbox{upper})}}{2}$
(3.16)

This one is used to evaluate PRS performance in the simulation study in
Section 5.

###  3.5 Plausibility intervals

Interval estimation can be made to have desirable frequency-calibration
property. Suppose that one is interested in the simultaneous interval
estimation, one for each $\theta_{i}$. Denote by $C_{i}^{(1-\alpha)}(X)$ the
interval estimate of $\theta_{i}$ with coefficient $1-\alpha$. A frequency-
calibration property commonly used in practice is that

$E\left(\frac{1}{n}\sum_{i=1}^{n}I\left\\{\theta_{i}\in
C_{i}^{(1-\alpha)}(X)\right\\}\right)\geq 1-\alpha.$ (3.17)

In words, it requires on average at least $1-\alpha$ percent of the true
values of $\theta_{i}$’s will be covered, in a frequentist sense within and
cross experiments. This, more precisely, the case with only the equality
constraint in (3.18) has been a popular standard in practice (Wahba, 1983;
Brown et al., 2011).

The validity of IMs that guarantees frequency-calibration has been established
in Martin and Liu (2013). Formally, we have the following result.

######  Theorem 1 (IMs Validity).

Let $B_{1-\pi}$ be the $(1-\alpha)$-th quantile of $B(U)$ and denote by
$\Gamma_{1-\pi}$ the $\alpha$-plausibility region

$\Gamma_{1-\pi}=\left\\{(\gamma_{1},...,\gamma_{K}):\;\mbox{$B\left(\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(1)}-\vartheta_{k}),...,\sum_{k=1}^{K}\gamma_{k}\Phi(X_{(n)}-\vartheta_{k})\right)\leq
B_{1-\pi}$}\right\\}.$ (3.18)

Then, $E\left((\gamma_{1},...,\gamma_{n})\in\Gamma_{1-\gamma}\right)\geq
1-\pi$.

Suppose that for each $\gamma\in\Gamma_{1-\pi}$, we construct the two-sided
“confidence interval” with the left-end point

$\theta_{\alpha/2}^{(\mbox{lower})}(\gamma)=\max_{\ell}\left\\{\theta_{\ell}:\hat{\theta}_{i}(\gamma)=\frac{\sum_{k=1}^{\ell}\gamma_{k}\phi(X_{i}-\vartheta_{k})}{\sum_{k=1}^{K}\gamma_{k}\phi(X_{i}-\vartheta_{k})}\leq\frac{\alpha}{2}\right\\}$
(3.19)

and the right-end point

$\theta_{\alpha/2}^{(\mbox{upper})}(\gamma)\min_{r}\left\\{\theta_{r}:\hat{\theta}_{i}(\gamma)=\frac{\sum_{k=r}^{n}\gamma_{k}\phi(X_{i}-\vartheta_{k})}{\sum_{k=1}^{K}\gamma_{k}\phi(X_{i}-\vartheta_{k})}\leq\frac{\alpha}{2}\right\\}.$
(3.20)

This leads to the following plausibility interval

$C_{i}(1-\pi,1-\alpha)\equiv\left[\min_{\gamma\in\Gamma_{1-\pi}}\theta_{\alpha/2}^{(\mbox{lower})}(\gamma),\;\max_{\gamma\in\Gamma_{1-\pi}}\theta_{\alpha/2}^{(\mbox{upper})}(\gamma)\right]$
(3.21)

The following theorem establishes a lower bound for the “confidence”
coefficient of (3.21).

######  Theorem 2 (Coverage Property).

The coverage probability of the confidence interval (3.21) is at least
$(1-\gamma)(1-\pi)$, i.e.,

$\mbox{\rm
Pr}\left(C_{i}(1-\pi,1-\alpha)\ni\theta_{i}\right)\geq(1-\pi)(1-\alpha)$

for all $\pi\in(0,1)$ and $\alpha\in(0,1)$.

###### Proof.

The IMs validity says that $\gamma$ is covered by $\Gamma_{1-\pi}$ $1-\pi$
percent. In such a case, the interval defined by the two end points (3.19) and
(3.20) covers the true value, generated from $G$, with probability $1-\alpha$.
The claimed result follows from the multiplication rule of probability. ∎

Theorem 2 provides a way of constructing exact confidence intervals. Because
the set for $\mathbb{G}$ or $\gamma$ can be large, rendering conservative
intervals, we provide some relevant comments on this issue, aiming to a method
of adaptive construction of shorter intervals with approximate coverage via
Monte Carlo simulations.

######  Remark 3.2 (Within-experiment Validation).

Within-experiment frequency evaluation of confidence interval is possible when
$n$ is large. The fractional number of $\theta_{i}$’s falling outside of their
$1-\alpha$ confidence intervals is expected to be $n\alpha$ and in such cases,
the corresponding observed $X_{i}$’s have probability 0.5 to fall outside of
the intervals as well. On the basis of this argument, we expect to see that
the number of $X_{i}$ falling outside of corresponding individual $1-\alpha$
confidence intervals for $\theta_{i}$ is about $n\alpha/2$.

######  Remark 3.3 (Selection of $\pi$ as a Way of Modeling).

When confidence intervals are excessively large, modeling can be useful.
Efron’s $g$-modeling, which assumes $\gamma\propto e^{Q\zeta}$ for some fixed
$(K\times p)$ matrix $Q$ with unknown $\zeta\in\mathbb{R}^{p}$, serves a good
example. Alternatively, one can also choose $\pi$ and attribute the ignored
relevant error to the modeling error.

######  Remark 3.4 (Adaptive Adjustment).

Efficient confidence intervals can be constructed via Monte Carlo as follows.
Without loss of generality, suppose that $95\%$ CIs are of interest. So, take
$1-\alpha=\sqrt{0.95}$, evaluate coverage via Monte Carlo simulation with
$\gamma\in...$ $1-\pi=\sqrt{0.95}$ for a sequence of $m$, say, $m=10$,
quantiles $B(U)$ with probabilities $(s/m)\sqrt{0.95}$, $s=1,...,m$. The MPE
estimate of $\mathbb{G}$, for example, can be used to generate data for this
Monte Carlo approach. The approximate intervals for the real data are the
$(s/m)\sqrt{0.95}$ intervals, where $s$ is the index of the simulation
intervals with their actual coverage in the simulated data matching the
targeted coverage. This approach is used in the real-data example in Section
5.2.

##  4  IMs for MNMs of the classic kind

###  4.1 Constructing an association via augmenting latent auxiliary variables

The Stein phenomenon implies that when performance assessment of inferential
methods involves simultaneous inference about $\theta_{i}$’s, such as mean
square errors for point estimation, efficient methods can be developed by
allowing the individual estimators of $\theta_{i}$’s to share risk with each
other. To see this intuitively from an IMs perspective, we write
$Z_{1},...,Z_{n}$ in (1.1) in terms of their sorted values
$Z_{(1)},...,Z_{(n)}$ with an auxiliary random permutation of $(1,...,n)$,
denoted by $(\tau_{1},...,\tau_{n})$:

$X_{i}=\theta_{i}+Z_{(\tau_{i})}$ (4.1)

Since everyone of the $n!$ permutations is equally likely, the sampling models
(1.1) and (4.1) are equivalent to each other. Clearly, efficient simultaneous
inference can be obtained, conditional on $(\tau_{1},...,\tau_{n})$. Here we
develop efficient methods by estimating the unknown permutation $\tau$ as a
way of controlling the level of risk sharing. Before proceeding we record a
result that is useful in what follows, where optimization over
$(Z_{(1)},...,Z_{(n)})$ and $\tau=(\tau_{1},...,\tau_{n})$ can be carried out
over $\theta\in\mathbb{R}^{n}$.

######  Proposition 1.

The observed data $X_{1}$,…, $X_{n}$, the auxiliary variables
$(Z_{(1)},...,Z_{(n)})$ and $\tau=(\tau_{1},...,\tau_{n})$ determine and are
determined by $\theta_{i}$’s.

What is interesting here for IMs is that the number of auxiliary variables is
larger than that of the observed data. For this reason, we call the auxiliary
permutation $\tau$ a latent auxiliary variable. For technical clarity, we
denote by RP the abbreviation of “random permutation” and summarize the above
discussion into the following lemma concerning the equivalence of three
sampling associations, including (1.1) and (4.1).

######  Lemma 1 (Equivalence of Three Sampling Associations).

The sampling associations (1.1), (4.1), and

$X_{\tau_{i}}=\theta_{\tau_{i}}+Z_{i}\qquad(\tau\sim\mbox{RP}(1,...,n);\;Z_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,1))$ (4.2)

are all equivalent in the sense that the induced sampling distributions of
$X$’s, given $\theta$’s, are the same, i.e., $X_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(\theta_{i},1)$ for all $i=1,...,n$.

Lemma 1 can be easily established with standard technical arguments in terms
of sampling distributions. This result allows for flexibility in formulating
IMs for prior-free probabilistic inference. For example, Equation (4.2)
implies that the sampling scheme in the classic kind can be viewed as follows.
There is a collection of $\theta$’s, denoted by
$\\{\theta_{(k)}:k=1,...,n\\}$. One takes $\theta$’s from
$\\{\theta_{(k)}:k=1,...,n\\}$ without replacement, draws a sample of size of
$n$ from $N(0,1)$, and creates $X$’s via either (4.2) or
$X_{\tau_{i}}=\theta_{\tau_{i}}+Z_{(i)}$ for $i=1,...,n$.

Enforcing risk sharing clearly prevents $\tau$ from being integrated out. Note
that different permutations produce observed data samples as a collection of
data points looking like each other or, more precisely, having the same
distribution. Given the observed data, however, the permutation that was
chosen to generate the data would be the best to be conditioned on for
inference via predicting $Z_{(\tau_{1})},...,Z_{(\tau_{n})}$. This suggests
that the collection of data points $X_{i}$’s be used to establish an
association, which we call an augmented association. In other words, we can
make use of the sampling model of the order and sufficient statistics
$X_{(1)}\leq...\leq X_{(n)}$ to create such an augmented association.

Note that the marginal distribution for a random draw from
$\\{X_{(1)},...,X_{(n)}\\}$ is given by

$F_{\theta}(x)=\sum_{k=1}^{n}\frac{1}{n}\Phi(x-\theta_{k})$ (4.3)

In terms of the true realizations of the auxiliary variables we can also write
(4.3) as

$F_{\theta}(x)=\sum_{k=1}^{n}\frac{1}{n}\Phi(x-[X_{k}-Z_{(\tau_{k})}])$ (4.4)

With Equation (4.4) we have the following association

$F_{\theta}(X_{(i)})=\sum_{k=1}^{n}\frac{1}{n}\Phi(X_{(i)}-[X_{k}-Z_{(\tau_{k})}])$
(4.5)

for $i=1,...,n$. Before proceeding we make the following two remarks, which
motivates the PRS construction in the next section.

######  Remark 4.1 (Dependence of Latent Auxiliary Variables).

Comparison of (4.1) and (4.5) reveals that (4.5) unsurprisingly depends on
$\Phi(Z_{i})$’s. This can be seen clearly, for example, by the extreme case
where all the $\theta_{i}$’s are the same and, thereby, in this extreme case
(4.1) and (4.5) are identical.

######  Remark 4.2 (Sampling W/O Replacement).

Although $F_{\theta}(.)$ defined in (4.3) is the distribution function for all
$X_{i}$’s, $\\{X_{(1)},...,X_{(n)}\\}$ is not a sample from $F_{\theta}(.)$.
This is because that $\\{X_{(1)},...,X_{(n)}\\}$ is equivalent to the ordered
sample of observations $X_{i}$’s taken from the $n$ subpopulations
$N(\theta_{i},1)$ without replacement (of subpopulations), whereas a sample
from $F_{\theta}(.)$, i.e., the $n$ subpopulations $N(\theta_{i},1)$, is
necessarily drawn with replacement.

###  4.2 Point and interval estimations

The PRS defined in Section 3 for the sorted uniform samples can be used for
sorted sample of $\Phi(Z_{i})$’s. The remaining challenges concern a
specification of a PRS for the latent permutation $\tau$, which is to be
utilized for inference through (4.5). Intuitively, for largest possible risk
sharing, it is important to control the accuracy of predicting $\tau$ using
PRS. For the extreme case in this way of reasoning, we would ignore predicting
$Z_{(i)}$’s, at least not predicting directly, i.e., formally predicting
$Z_{(i)}$’s with their entire probability sample space. In other words,
predicting $\tau$ alone effectively estimates the collection of $\theta$’s.
What is needed is to create a random bound $\mathcal{B}$ on
$B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))$, working in the same way as
the basic IMs to induce a random set for $\theta$’s. Theoretically, valid
inference can be made as long as we use a random bound $\mathcal{B}$ that is
stochastically not smaller than
$B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))$. This is supported rigorously
by the following theorem.

######  Theorem 3 (Validity of conservative PRSs).

Suppose that $B$ is a random variable that is stochastically larger than
$B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))$ with $X$ and $\theta$ pairs
given in (4.1). Then inference on $\theta$ in the Inferential Models framework
through the induced random set on the parameter space

$\Theta(\mathcal{B})=\left\\{\theta:\,\theta\in\mathbb{R}^{n};\;B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))\leq\mathcal{B}\right\\}$
(4.6)

is valid.

###### Proof.

This result can be proved by using arguments similar to those in the proof of
Theorem 2 of Martin and Liu (2013); see also Martin (2015). ∎

In the simulation study of Section 5, we use
$\mathcal{B}=B(U_{(1)},...,U_{(n)})$, where $U_{(1)},...,U_{(n)}$ is the
sorted $n$ uniforms. A heuristic argument for this is that the sorted sample
of $X$’s without replacement (of $\theta$’s) spreads more evenly than that
with replacement from the same marginal distribution (4.3). Together with
(4.3), this implies that the sorted $X$s have about the same component-wise
locations but smaller variability. As a result,
$\mathcal{B}=B(U_{(1)},...,U_{(n)})$ serves as a random bound in Theorem 3 for
producing valid inference on $\theta$, at least for values of
$B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))$ that are not large; see
Figure 2 for a pictorial illustration that supports this argument. Although it
guarantees validity, the use of $\mathcal{B}=B(U_{(1)},...,U_{(n)})$ can lead
to conservative inference. Since the distribution of
$B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))$ depends on $\theta$’s,
efficient inference can be made by using the so-called assertion-specific PRS
(Martin and Liu, 2015a). More precisely, we summarize this approach into the
following theorem, followed by remark on its potential applications and
computational challenges.

######  Theorem 4 (Validity of Assertion-specific PRSs for MPE and
Constructing CIs).

Suppose that the assertion
$\\{\vartheta^{\prime}_{(k)}:\;k=1,...,n\\}\in\mathbb{R}^{n}$ is of interest.
Define the induced random set on $\theta$

$\Theta(X^{\prime})=\left\\{\theta:\,\theta\in\mathbb{R}^{n};\;B(F_{\theta}(X_{(1)}),...,F_{\theta}(X_{(n)}))\leq
B(F_{\vartheta^{\prime}}(X^{\prime}_{(1)}),...,F_{\vartheta^{\prime}}(X^{\prime}_{(n)}))\right\\}$
(4.7)

where $X^{\prime}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}N(\vartheta^{\prime}_{i},1)$ for $i=1,...,n$. Then inference
based on $\Theta(X^{\prime})$ in the IMs framework is valid.

###### Proof.

This result can also be proved by using arguments similar to those in the
proof of Theorem 2 of Martin and Liu (2013). ∎

######  Remark 4.3.

According to Theorems 3 and 4, a way of constructing efficient confidence
regions is to find the minimum plausibility for any given
$\theta\in\Theta(\mathcal{B})$. This minimum plausibility is the plausibility
of $\Theta(\mathcal{B})$. Evaluation of this plausibility appears intensive.
It would be interesting to develop efficient computational methods for this
approach.

The theory established so far provides inference about the $\theta$’s as a
collection. Recall that given the $\theta$ collection, $X$’s were generated by
drawing $\theta$’s without replacement. Standard algebraic operations of
applying the Bayes theorem lead to the following results.

######  Theorem 5 (Conditional Distribution of $\theta$).

Suppose $n\geq 2$ and assume that the value of $\theta$,
$\\{\vartheta_{k}:\vartheta_{k}=\theta_{(k)},k=1,...,n\\}$, are known. Then
the partial conditional distribution of $\theta_{i}$ given $X_{i}$ has the
probability mass function (pmf)

$\mbox{Pr}(\theta_{i}=\vartheta_{j}|X_{i})=\frac{\phi(X_{i}-\vartheta_{j})}{\sum_{k=1}^{n}\phi(X_{i}-\vartheta_{k})}$
(4.8)

and the full conditional distribution of $\theta$ given $X$ has the pmf

$\mbox{Pr}(\theta_{1}=\vartheta_{\tau_{1}},...,\theta_{n}=\vartheta_{\tau_{n}}|X_{1},...,X_{n})=\frac{\frac{1}{n!}\prod_{j=1}^{n}\phi(X_{j}-\vartheta_{\tau_{j}})}{\sum_{\tau^{\prime}}\frac{1}{n!}\prod_{j=1}^{n}\phi(X_{j}-\vartheta_{\tau^{\prime}_{j}})}$
(4.9)

Theorem 5 makes it possible to estimate $\theta$’s via $E(\theta|X)$,
fulfilling the required operation to overcome the fundamental difficulty of
making inference for the classic kind (c.f., Section 1 and Stigler (1990)).
While it is appealing to use the full conditional distribution (4.9) for MPE,
for computational simplicity here we use the following simple estimator,
obtained from (4.8),

$\hat{\theta}_{i}=\frac{\sum_{k=1}^{n}\vartheta_{k}\phi(X_{i}-\vartheta_{k})}{\sum_{k=1}^{n}\phi(X_{i}-\vartheta_{k})}$
(4.10)

for $i=1,...,n$, with $\\{\vartheta_{(k)}:k=1,...,n\\}$ obtained from the
approach in either Theorem 3 or Theorem 4. The computational strategies for
the empirical-Bayes kind provided in Appendices can be applied here for the
classic kind. More discussion is given in Section 6 on partial conditioning
versus full conditioning in Theorem 5 for estimation $\theta$.

##  5 Numerical illustration

###  5.1 A simulation study

In this simulation study, we consider the three scenarios used in the
simulation study of Jiang and Liu (2022) for their many-normal-means example:

> Single-Mode.
>
>
> $\theta_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N(0,0.01)$, the case in
> favor of James-Stein estimator;
>
> Two-Mode.
>
>
> $\theta_{i}\stackrel{{\scriptstyle
> iid}}{{\sim}}\frac{1}{2}N(-2,0.01)+\frac{1}{2}N(2,0.01)$; and
>
> Single-Mode with Outliers.
>
>
> $\theta_{i}\stackrel{{\scriptstyle iid}}{{\sim}}0.9\delta_{0}+0.1N(-3,1)$,
> where $\delta_{0}$ denotes the degenerated distribution at $0$.

For all simulation studies, three different sample sizes, $n=10,20,$ and $50$
are considered. The numerical deconvolution uses $K=1000$ points equally
spaced in the range from $X_{(1)}+\Phi^{-1}(0.0001)$ to from
$X_{(n)}+\Phi^{-1}(0.9999)$, where $\Phi^{-1}(.)$ stands for the inverse
mapping of the standard normal cdf $\Phi(.)$. For each case, $M=200$ data sets
are generated to compute the Monte Carlo approximations to, for example, mean
square error (MSE).

Table 1: MSE results in three simulation studies with different methods.  Problem Kind /Method |  | $N(0,0.01)$ |  |  | $\frac{1}{2}N(-2,0.01)+\frac{1}{2}N(2,0.01)$
---
|  | $0.9\delta_{0}+0.1N(-3,1)$
---
| $n=10$ | $n=20$ | $n=50$ |  | $n=10$ | $n=20$ | $n=50$ |  | $n=10$ | $n=20$ | $n=50$
The EB Kind |  |  |  |  |  |  |  |  |  |  |  |
$g$-modeling |  | 0.419 | 0.395 | 0.168 |  | 0.748 | 0.724 | 0.737 |  | 0.554 | 0.552 | 0.364
Auto-modeling |  | 0.199 | 0.110 | 0.054 |  | 0.600 | 0.437 | 0.356 |  | 0.420 | 0.418 | 0.312
IMs |  | 0.219 | 0.134 | 0.062 |  | 0.678 | 0.581 | 0.500 |  | 0.462 | 0.370 | 0.361
Classic Kind |  |  |  |  |  |  |  |  |  |  |  |
MLE |  | 1.022 | 0.972 | 0.985 |  | 0.990 | 1.009 | 1.003 |  | 0.971 | 1.021 | 0.983
James-Stein |  | 0.300 | 0.167 | 0.066 |  | 0.876 | 0.850 | 0.826 |  | 0.521 | 0.516 | 0.482
IMs |  | 0.312 | 0.166 | 0.080 |  | 0.719 | 0.557 | 0.422 |  | 0.484 | 0.376 | 0.317
Figure 1: A typical example of IMs for the empirical-Bayes kind in the two-modal case with $n=20$, where X — $X$’s, curves — lower and upper $E(\theta_{i}|X)$, circles — true values of $\theta$’s, and vertical segments — confidence intervals computed at the MPE (set). The coefficients of the confidence intervals is about 97%, according to the simulation results given in Table 2.  Table 2: Actual coverage and average interval length (in parentheses) obtained by the IMs method for the empirical-Bayes kind in three simulation studies.  Targeted Coverage |  | $N(0,0.01)$ |  |  | $\frac{1}{2}N(-2,0.01)+\frac{1}{2}N(2,0.01)$
---
|  | $0.9\delta_{0}+0.1N(-3,1)$
---
|  $n=10$ |  $n=20$ |  |  $n=10$ |  $n=20$ |  |  $n=10$ |  $n=20$
MPE |  | 83.77 (1.57) | 94.14 (1.63) |  | 91.58 (2.95) | 96.94 (2.96) |  | 86.09 (2.00) | 93.84 (1.93)
50% |  | 93.43 (2.14) | 96.79 (1.99) |  | 97.91 (4.34) | 99.08 (3.42) |  | 94.72 (2.67) | 95.73 (2.36)
75% |  | 97.56 (2.47) | 97.93 (2.18) |  | 98.56 (4.68) | 99.23 (3.97) |  | 96.06 (3.03) | 97.01 (2.56)
90% |  | 99.18 (2.84) | 98.93 (2.42) |  | 98.42 (4.89) | 99.44 (4.50) |  | 97.97 (3.43) | 98.41 (2.79)
95% |  | 99.50 (3.19) | 99.14 (2.59) |  | 98.71 (5.09) | 99.29 (4.64) |  | 98.28 (3.76) | 98.54 (2.97)
Figure 2: Illustration of conservative PRSs in four different collections of
$n=20$ $\theta$’s: (1) one simulated values from the single mode model, (2)
one simulated values from the two-mode model, and (3) one simulated values
from the single mode model with outliers.

The MSE performance results of different methods for point estimation are
summarized into Table 1. The results for Efron’s $g$-modeling, the auto-
modeling method of Jiang and Liu (2022), the method of maximum likelihood
estimation (MLE), and the James-Stein shrinkage estimator are taken from Jiang
and Liu (2022). From these simulation-based numerical results, we see that IMs
for the classic kind outperforms the James-Stein estimator. This is not
surprising because of the difficulty recognized, for example, by Stigler
(1990) and the sound reasoning with uncertainty introduced in this paper for
the classic kind. Figure 1 displays the IMs results for a typical data set
generated from the two-mode model. It is interesting to notice that IMs
“shrinks” toward local modes, whereas the James-Stein estimates are of course
linear function of $X$’s.

The results also show that our IMs for the empirical-Bayes kind performs
better than Efron’s $g$-modeling method. Possible explanations may be that
Efron’s $g$-modeling method is subject to a larger model error due to the
$g$-modeling. The use of the method of (penalized) MLE for estimation of the
hyper-parameters introduced in the $g$-modeling could also be a contributing
factor.

In addition, it is encouraging to observe that the IMs are comparable with the
auto-modeling method of Jiang and Liu (2022). Good performance is expected for
the IMs, although they may be further improved a bit when better PRSs are
used. An intuitive explanation for the similar good performance for both the
IMs and the auto-modeling method of Jiang and Liu (2022) is that both methods
take it into consideration that generated data from the estimated model should
look like the observed. The IMs formulated this paper does it through a PRS on
$F_{\theta}(X_{(i)})$’s, whereas auto-modeling explicitly makes this clear
from a modeling perspective.

The results for interval estimation produced by the IMs for the empirical-
Bayes kind are summarized into Table 2 for $n=10$ and $n=20$. The “Targeted
Coverage” is the nominal coefficients, except for the row of $MPE$. The values
in the row of $MPE$ are computed with respect to all $\gamma$’s in the maximum
plausibility region. The confidence coefficient for $\theta$’s conditional on
$\gamma$’s is set to be $1-\alpha=0.95^{1/2}=0.9747.$ The plausibility
coefficient $1-\pi$ is taken to fulfill the nominal level, by making use of
Theorem 2. We see that these intervals have more coverage than their nominal
interval coefficients. This is somewhat expected because of the implicit use
of “confidence” region for the high 1,000-dimensional $\gamma$. More efficient
intervals can be obtained via adaptive construction of shorter intervals; see
Remark 3.4. This can be in fact seen clearly from able 2. This is also used in
Section 5.2 for the real-data example.

###  5.2 A real-data example

Table 3: The data of Rubin (1981) on SAT coaching experiments  |  | Estimated |  | Associated
---|---|---|---|---
School $(i)$ |  |  Treatment-Effects $(y_{i})$ |  |  Standard-Error $(s_{i})$
1 |  | 28.39 |  | 14.9
2 |  | 7.94 |  | 10.2
3 |  | -2.75 |  | 16.3
4 |  | 6.82 |  | 11.0
5 |  | -0.64 |  | 9.4
6 |  | 0.63 |  | 11.4
7 |  | 18.01 |  | 10.4
8 |  | 12.16 |  | 17.6

Rubin (1981) considered assessing effects of SAT coaching programs based on
$n=8$ parallel randomized experiments in 8 schools. The estimated individual
effects and the associated standard deviations are tabulated in Table 3. It is
arguable to use the following sampling model

$y_{i}\stackrel{{\scriptstyle
ind}}{{\sim}}N(\mu_{i},s_{i}^{2})\qquad(i=1,...,n),$ (5.1)

because each $y_{i}$ was obtained from a data set that is large enough for the
acceptance of the normality assumption and the assumption that $s_{i}^{2}$’s
are known.

Let $X_{i}=y_{i}/s_{i}$ and let $\theta_{i}=\mu_{i}/s_{i}$ for $i=1,...,8$.
Inference about $\mu_{i}$ from $y$’s via (5.1) can be obtained from inference
about $\theta$’s from $X$’s based on the sampling model (1.1). In this case,
the Efron-Morris shrinkage estimator over-shrinks $\theta$’s due to the fact
that the sample variance of $X$’s is small, resulting a strong negative
correlation between $X$’s and the shrinkage estimates of their corresponding
$\theta$s. As suggested in the literature, this can be of course easily
modified in the same way of creating the James-Stein positive-part estimator
(James and Stein, 1961).

Figure 3: The results for the real data example by the IMs for the empirical-
Bayes kind, where X — $X$s, curves — lower and upper $E(\theta_{i}|X)$,
circles — Efron-Morris shrinkage estimates, and vertical segments — confidence
intervals computed at the MPE (set). The coefficients of the confidence
intervals is about 97%, suggested by an application of the method of adaptive
adjustment; see Remark 3.4.

IMs results based on the empirical-Bayes kind are shown in Figure 3. Similar
results can be extended to $\mu$’s because of the simple scaling
transformation $\mu_{i}=s_{i}\theta_{i}$. The coefficients of the confidence
intervals is about 97%, suggested by an application of the method of adaptive
adjustment; see Remark 3.4. It is interesting to see that the MPEs of
$\theta$’s are about the same. To verify this directly, consider the mean-
removed statistic

$S^{2}=X^{\prime}(I-\frac{1}{n}J)X=4.2899$

where $I$ is the $n\times n$ identity matrix and $J$ is the $n\times n$ matrix
of ones. Note that the sampling distribution of $S^{2}$ is the non-central
chi-square with 7 degrees of freedom and non-centrality parameter
$\delta^{2}=\theta^{\prime}(I-\frac{1}{n}J)\theta$. The cdf value of $S^{2}$
under the central chi-square distribution is 0.2542. This doesn’t provide much
strong evidence against the assertion that $\delta^{2}=0$ or, equivalently,
all the $\theta$s are the same. So, it is not surprising to see that the MPEs
shown in Figure 3 are about the same; see Leaf et al. (2015) for more
discussion from an IMs perspective.

##  6 Discussion

Taking an IMs approach, we have investigated on making exact and efficient
inference for the two kinds of many-normal-means. For the empirical Bayes
kind, with IMs one can aim to producing efficient inference by balancing
validity and flexibility. Since the IMs framework offers prior-free
probabilistic inference, which says more than confidence, IMs plausibility can
be used for point and interval estimations. In addition to the desired
exactness, efficiency of the resulting methods is satisfactory, as shown in
the simulation study. Furthermore, in the case further $g$-modeling is
preferred, the IMs results could be useful for model building inasmuch as
valid non-parametric methods are useful for model checking; see Tukey (1962)
for more discussion along this line.

It should be noted that IMs works well for the over-parameterized
$\mathbb{G}$. This can be viewed as a good property of IMs, unlike frequentist
and Bayesian schools of thought. Frequentist statisticians often do inference
by first constructing point estimators and then develop confidence methods
based on point estimators. For over-parameterized models, it is often that
some kinds of regularization via penalty methods are used. Bayesians have to
fully specify a probability model for everything, making it identifiable in a
probabilistic sense with the usual theory of probability. With IMs, our way of
reasoning is done with respect to the auxiliary variables. Non-identifiability
is no longer an issue, because we work with set-valued functions/mappings.

There are also other kinds of the many-normal-means problem. A popular kind is
obviously the so-called normal-normal empirical Bayes, where the unobserved
$\theta$’s are assumed to have followed a normal distribution with unknown
hyper-parameters. This model has been widely used in applied statistics; see,
for example, Rubin (1980), Morris (1983), Gelman et al. (1995), Qiu et al.
(2018), and references therein.

Our understanding of the classic kind from the IMs perspective motivated the
new idea of introducing latent permutation auxiliary variable. While the
proposed methods are satisfactory in the context of our numerical results for
illustration and comparison, there are remaining theoretical challenges as far
as optimal solutions are of interest. This is mainly due to the technical
difficulties in dealing with data collected via sampling without replacement,
whereas investigations of the empirical Bayes kind enjoys the relative
simplicity offered by data collection through sampling with replacement. In
particular, it would be interesting to develop efficient computational
algorithms for evaluation of the full conditional expectation in Theorem 5.
Nevertheless, it is expected that there is a large difference between partial
and full conditional expectations for small $n$ and not so for large $n$.

The importance of the many-normal-means problem suggests immediately that our
new methods can be developed for the many-normal-means type of problems, such
as linear and non-parametric regression. These new methods can also be
potentially extended to situations beyond normal models, including Poisson and
Binomial models. Incidentally, the idea of considering latent associations
would inspire new research on developments of the IMs theory itself as well as
methods of other schools of thought. For example, an alternative to popular
penalty methods for maximum likelihood estimation is to make use of the idea
of using latent associations as a data-adaptive way of regularization, which
has been more or less a standard method in contemporary high-dimensional
statistics. We plan to work on these topics and report results elsewhere.

## Appendices

##  Appendix A Maximum plausibility estimation

The MPE $\hat{\theta}_{i}$ in (3.16) can be computed in the following three
steps:

> Step 1.
>
>
> Find $B_{\min}$ by solving the constrained optimization problem (3.13) with
> the constraints in (3.1) on $\gamma$;
>
> Step 2.
>
>
> Shift $\gamma$ to the left-most and right-most positions within
> $\Gamma_{1-\pi}$ as the starting values for the next step;
>
> Step 3.
>
>
> Find $\hat{\theta}_{i}^{(\mbox{lower})}$ and
> $\hat{\theta}_{i}^{(\mbox{upper})}$ by solving the two corresponding
> constrained optimization problems with the objective function (3.10) and the
> constraints in both (3.1) and (3.14).

These constrained optimization problems are known as nonlinear programming.
Standard numerical methods for such problems are widely available, including
the R package nloptr by (Johnson, 2022, and coauthors).

The algorithms of Svanberg (2002), Conn et al. (1991), and Birgin and Martinez
(2008), made available as options NLOPT_LD_MMA and NLOPT_LD_AUGLAG of the
nloptr function in Johnson (2022), are used for the simulation study in
Section 5. These algorithms performed well with starting values for Step 1 set
as follows:

Step 1.

Set $\gamma_{k}=1/K$ for all $k=1,...,K$;

Step 2.

For each $i=1,...,n$, increase $\gamma_{k}$ by $1/n$, where
$k=\min_{j\in\\{1,...,K\\}}|\vartheta_{j}-X_{i}|$,

followed by the renormalization
$\gamma_{k}=\gamma_{k}/\sum_{j=1}^{K}\gamma_{j}$. To be conservative, we also
used different starting values set as above but for resampled $X_{i}$’s with
replacement. No significant improvements could be found in the simulation
study.

##  Appendix B Plausibility intervals

The confidence intervals can be computed as plausibility intervals in the
following three steps:

> Step 1.
>
>
> Find $B_{1-\pi}$ in (3.18);
>
> Step 2.
>
>
> Apply the methods in Appendix A to find the extreme $\gamma$ values for
> computing extreme means.
>
> Step 3.
>
>
> Take the extreme $\gamma$ values computed in Step 2 as starting values, and
> find the end points (3.19) and (3.20) of the confidence interval.

These constrained optimization problems are known as nonlinear programming.
Standard numerical methods for such problems are widely available, including
the R package nloptr by (Johnson, 2022, and coauthors). The algorithms of
Svanberg (2002) and Runarsson and Yao (2005), made available as options
NLOPT_LD_MMA and NLOPT_NG_ISRES of the nloptr function in Johnson (2022), are
used for the simulation study in Section 5.

## References

  * Birgin and Martinez (2008) Birgin, E. G. and J. M. Martinez (2008).  Improving ultimate convergence of an augmented lagrangian method.  Optimization Methods and Software 23(2), 177–195.
  * Brown (1971) Brown, L. D. (1971).  Admissible estimators, recurrent diffusions, and insoluble boundary value problems.  The Annals of Mathematical Statistics 42(3), 855–903.
  * Brown et al. (2011) Brown, L. D., X. Fu, and L. H. Zhao (2011).  Confidence intervals for nonparametric regression.  Journal of Nonparametric Statistics 23(1), 149–163.
  * Butucea and Comte (2009) Butucea, C. and F. Comte (2009).  Adaptive estimation of linear functionals in the convolution model and applications.  Bernoulli 15(1), 69–98.
  * Carroll and Hall (1988) Carroll, R. J. and P. Hall (1988).  Optimal rates of convergence for deconvolving a density.  Journal of the American Statistical Association 83(404), 1184–1186.
  * Conn et al. (1991) Conn, A. R., N. I. Gould, and P. Toint (1991).  A globally convergent augmented lagrangian algorithm for optimization with general constraints and simple bounds.  SIAM Journal on Numerical Analysis 28(2), 545–572.
  * Cui and Hannig (2022) Cui, Y. and J. Hannig (2022).  Demystifying inferential models: A fiducial perspective.  arXiv preprint arXiv:2205.05612.
  * Dempster (2008) Dempster, A. P. (2008).  The dempster–shafer calculus for statisticians.  International Journal of approximate reasoning 48(2), 365–377.
  * Donoho and Johnstone (1995) Donoho, D. L. and I. M. Johnstone (1995).  Adapting to unknown smoothness via wavelet shrinkage.  Journal of the american statistical association 90(432), 1200–1224.
  * Efron (2014) Efron, B. (2014).  Two modeling strategies for empirical bayes estimation.  Statistical science: a review journal of the Institute of Mathematical Statistics 29(2), 285.
  * Efron (2016) Efron, B. (2016).  Empirical bayes deconvolution estimates.  Biometrika 103(1), 1–20.
  * Efron and Morris (1972) Efron, B. and C. Morris (1972).  Empirical bayes on vector observations: An extension of stein’s method.  Biometrika 59(2), 335–347.
  * Fan (1991) Fan, J. (1991).  On the optimal rates of convergence for nonparametric deconvolution problems.  The Annals of Statistics, 1257–1272.
  * Gelman et al. (1995) Gelman, A., J. B. Carlin, H. S. Stern, and D. B. Rubin (1995).  Bayesian data analysis.  Chapman and Hall/CRC.
  * Hall and Meister (2007) Hall, P. and A. Meister (2007).  A ridge-parameter approach to deconvolution.  The Annals of Statistics 35(4), 1535–1558.
  * Hannig (2009) Hannig, J. (2009).  On generalized fiducial inference.  Statistica Sinica 19, 491–544.
  * James and Stein (1961) James, W. and C. Stein (1961).  Estimation with quadratic loss.  Proc. 4th Berkeley Symp. Math. Stat. Probab. 1, 361-379 (1961).
  * Jiang and Liu (2022) Jiang, Y. and C. Liu (2022).  Estimation of over-parameterized models via fitting to future observations.  arXiv preprint arXiv:2206.01824.
  * Johnson (2022) Johnson, S. G. (2022).  The nlopt nonlinear-optimization package, https://nlopt.readthedocs.io/en/latest/.  R Package.
  * Kerman (2011) Kerman, J. (2011).  A closed-form approximation for the median of the beta distribution.  arXiv preprint arXiv:1111.0433.
  * Laird (1978) Laird, N. (1978).  Nonparametric maximum likelihood estimation of a mixing distribution.  Journal of the American Statistical Association 73(364), 805–811.
  * Leaf and Liu (2012) Leaf, D. E. and C. Liu (2012).  Inference about constrained parameters using the elastic belief method.  International Journal of Approximate Reasoning 53(5), 709–727.
  * Leaf et al. (2015) Leaf, D. E., H. Yun, and C. Liu (2015).  Valid prior-free probabilistic inference and its applications in medical statistics.  In Advanced Medical Statistics, pp. 1301–1337. World Scientific.
  * Liu and Martin (2015) Liu, C. and R. Martin (2015).  Frameworks for prior-free posterior probabilistic inference.  Wiley Interdisciplinary Reviews: Computational Statistics 7(1), 77–85.
  * Liu and Martin (2021) Liu, C. and R. Martin (2021).  Comment: Settle the unsettling: An inferential models perspective.  Statistical Science 36(2), 196–200.
  * Liu and Xie (2014) Liu, C. and J. Xie (2014).  Probabilistic inference for multiple testing.  International Journal of Approximate Reasoning 55(2), 654–665.
  * Martin (2015) Martin, R. (2015).  Plausibility functions and exact frequentist inference.  Journal of the American Statistical Association 110(512), 1552–1561.
  * Martin and Liu (2013) Martin, R. and C. Liu (2013).  Inferential models: A framework for prior-free posterior probabilistic inference.  Journal of the American Statistical Association 108(501), 301–313.
  * Martin and Liu (2014) Martin, R. and C. Liu (2014).  A note on p-values interpreted as plausibilities.  Statistica Sinica, 1703–1716.
  * Martin and Liu (2015a) Martin, R. and C. Liu (2015a).  Conditional inferential models: combining information for prior-free probabilistic inference.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 77(1), 195–217.
  * Martin and Liu (2015b) Martin, R. and C. Liu (2015b).  Marginal inferential models: prior-free probabilistic inference on interest parameters.  Journal of the American Statistical Association 110(512), 1621–1631.
  * Martin et al. (2010) Martin, R., J. Zhang, and C. Liu (2010).  Dempster–shafer theory and statistical inference with weak beliefs.  Statist. Sci., 72–87.
  * Morris (1983) Morris, C. N. (1983).  Parametric empirical bayes inference: theory and applications.  Journal of the American statistical Association 78(381), 47–55.
  * Narasimhan and Efron (2020) Narasimhan, B. and B. Efron (2020).  deconvolveR : A $g-$modeling program for deconvolution and empirical bayes estimation.  Journal of Statistical Software 94(11).
  * Qiu et al. (2018) Qiu, Y., L. Zhang, and C. Liu (2018).  Exact and efficient inference for partial bayes problems.  Electronic Journal of Statistics 12(2), 4640–4668.
  * Rubin (1980) Rubin, D. B. (1980).  Using empirical bayes techniques in the law school validity studies.  Journal of the American Statistical Association 75(372), 801–816.
  * Rubin (1981) Rubin, D. B. (1981).  Estimation in parallel randomized experiments.  Journal of Educational Statistics 6(4), 377–401.
  * Runarsson and Yao (2005) Runarsson, T. P. and X. Yao (2005).  Search biases in constrained evolutionary optimization.  IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews) 35(2), 233–243.
  * Schweder and Hjort (2016) Schweder, T. and N. L. Hjort (2016).  Confidence, likelihood, probability, Volume 41.  Cambridge University Press.
  * Shafer (1976) Shafer, G. (1976).  Mathematical theory of evidence [mi. princeton.  Princeton University Press 35(1), 35–37.
  * Shafer (1992) Shafer, G. (1992).  Dempster-shafer theory.  Encyclopedia of artificial intelligence 1, 330–331.
  * Stefanski and Carroll (1990) Stefanski, L. A. and R. J. Carroll (1990).  Deconvolving kernel density estimators.  Statistics 21(2), 169–184.
  * Stein (1956) Stein, C. (1956).  Inadmissibility of the usual estimator for the mean of a multivariate normal distribution.  Proc. 3rd Berkeley Sympos. Math. Statist. Probability 1, 197-206 (1956).
  * Stigler (1990) Stigler, S. M. (1990).  The 1988 neyman memorial lecture: a galtonian perspective on shrinkage estimators.  Statistical Science, 147–155.
  * Svanberg (2002) Svanberg, K. (2002).  A class of globally convergent optimization methods based on conservative convex separable approximations.  SIAM journal on optimization 12(2), 555–573.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society. Series B (Methodological), 267–288.
  * Tukey (1962) Tukey, J. W. (1962).  The future of data analysis.  The Annals of Mathematical Statistics, 1–67.
  * Wahba (1983) Wahba, G. (1983).  Bayesian “confidence intervals” for the cross-validated smoothing spline.  Journal of the Royal Statistical Society: Series B (Methodological) 45(1), 133–150.
  * Xie et al. (2011) Xie, M., K. Singh, and W. E. Strawderman (2011).  Confidence distributions and a unifying framework for meta-analysis.  Journal of the American Statistical Association 106(493), 320–333.
  * Zhang (1997) Zhang, C.-H. (1997).  Empirical bayes and compound estimation of normal means.  Statistica Sinica 7(1), 181–193.
  * Zhang (2010) Zhang, J. (2010).  Statistical inference with weak beliefs.  Ph.D. thesis, Purdue Univer- sity, West Lafayette, IN..
  * Zhang and Liu (2011) Zhang, J. and C. Liu (2011).  Dempster-shafer inference with weak beliefs.  Statistica Sinica, 475–494.


# A mesoscopic device for a realization of the Topological Kondo effect

Saheli Sarkar  Alexei M. Tsvelik  Division of Condensed Matter Physics and
Materials Science, Brookhaven National Laboratory, Upton, NY 11973-5000, USA

###### Abstract

The search of anyons is a field of immense interest owing to its potential
application in the field of quantum information. Quantum critical Kondo
impurities present one possible platform for realization of anyons. In this
paper we discuss practical steps for realization of Topological Kondo effect
which, in contrast to the well known multichannel Kondo one, remains critical
even in the presence of perturbations.

##  I Introduction

Anyons - fractionalized particles with exotic statistics, are considered an
essential part of quantum information systems Nayak  _et al._ (2008). There
are different suggestions of their realization, one of them being the Kondo
effect. Existence of quantum critical Kondo (QCK) effect with a non-Fermi
liquid (NFL) behavior, first described and studied theoretically Nozieres and
Blandin (1980); Tsvelick and Wiegmann (1984); Affleck and Ludwig (1991) has
been well established experimentally Potok  _et al._ (2007); Keller  _et al._
(2015); Iftikhar  _et al._ (2015, 2018); Karki  _et al._ (2023). It is well
understood that the corresponding critical ground state contains non-Abelian
anyons Lopes  _et al._ (2020); Komijani (2020); Ge and Komijani (2022).
However, the experiments implement the so-called multichannel Kondo effect
where the criticality is sensitive to asymmetry between the scattering
channels. Hence a realization of the quantum critical point (QCP) requires a
fine tuning of the coupling parameters which limits its potential
applications.

Béri and Cooper et.al. Béri and Cooper (2012) proposed a new type of non-Fermi
liquid Kondo effect where the criticality is stable against most common
perturbations. The Topological Kondo effect (TKE) is predicted to arise as a
result of coupling between bulk conduction electrons and Majorana Zero modes
[MZM, also called Majorana Bound States (MBSs) ] located inside of the so-
called Majorana-Cooper box (MCB). They also provided a sketch of the MCB setup
for TKE Béri and Cooper (2012) in the form of a device consisting of a
mesoscopic superconducting island with proximitized to it semiconductor
nanowires containing the MBSs. The Coulomb blockade in the island transforms
entangled MBSs into an effective quantum spin, which interacts with the
conduction electrons of the leads via superexchange interaction facilitating
the TKE. In contrast to the conventional QCK effect, the TKE is robust against
various perturbations Altland and Egger (2013); Altland  _et al._ (2014).

Apart from being an anyon platform, the TKE with its characteristic features
may resolve persistent controversies about the existence of MBS. The theory
Béri and Cooper (2012); Altland and Egger (2013) suggests that in the device
with $M$ leads, the linear conductance ($\sigma_{ij}$) between different leads
$i$ and $j$ saturates at small temperatures at the universal value
$\frac{2e^{2}}{Mh}$, following a nontrivial power law temperature dependence,
different from the $T^{2}$ Fermi liquid one. This universal behavior of the
conductance constitutes a manifestation of the anyonic character of the TKE
ground state, which therefore can serve as a smoking gun for the presence of
MBS.

As we have mentioned above, the suggestions of realization of TKE contained in
earlier theoretical studies described heterostructure devices containing
nanowires with strong spin-orbit interaction in the presence of magnetic field
proximatized to a mesoscopic size superconducting island Béri and Cooper
(2012); Altland  _et al._ (2014); Väyrynen  _et al._ (2020); Liu  _et al._
(2021). We suggest another route to TKE: to use intrinsic topological
superconductor (TSC) Beenakker (2013); Sato and Ando (2017); Xu  _et al._
(2016) based devices, which may naturally host Majorana fermions. This may
drastically simplify the fabrication process providing a material realization
for the TKE. Promising candidates include iron-based superconductors (FeSC)
Fe(Tex,Se${}_{1-x})$ ($T_{c}=14.5$ K) Wang  _et al._ (2018); Zhang  _et al._
(2018); Machida  _et al._ (2019), (Li0.84Fe0.16)OHFeSe ($T_{c}=42$ K) Liu  _et
al._ (2018) and CaKFe2As4 ($T_{c}=35$ K) Liu  _et al._ (2020a). All of them
exhibit signatures of intrinsic topological superconductivity and have
relatively high superconducting critical temperatures Tc. Thus, it is
conceivable that MCB based on such materials can provide a viable path towards
observing the TKE. Observation of TKE in turn can serve as a smoking gun
signature of MBS, whose detection in various systems has been so far largely
focused on observing a zero-bias conductance peak (ZBCP)Yin  _et al._ (2015);
Zhu  _et al._ (2020) in local probe techniques. Additionally, finding
signatures of non-local character Alicea (2012) of MBS through TKE will be
important for quantum computation applications.

The aim of this paper is to provide detailed suggestions for realization of
MCBs for TKE focusing on FeSC based systems, with a final goal to use it for
quantum information applications. The rest of the paper is organized as
follows. In section II, we provide different scenarios for realizing MCBs
using FeSC. In section III, we present a FeSC based mesoscopic device and the
corresponding parameters characterizing the device for realizing TKE in more
realistic experimental situations. In section IV, we depict an arrangement for
a TKE-based chiral Kondo lattice which contains multiple non-Abelian anyons.
Section V contains a summary.

##  II Different scenario for realization of Majorana-Cooper boxes

There are multiple experiments indicating a possible existence of MBSs in
various iron-based superconductors. Non of them are decisive, but as we
suggest, in the context of TKE the corresponding systems may provide smoking
gun evidence for MBS. Below we consider different arrangements weighting their
merits as possible candidates for realization of TKE. In all cases the
corresponding MCBs are conceived as mesoscopic devices consisting of
superconductor with charging energy $E_{C}$ containing MBSs.

###  II.1 Vortices

One way of coexistence of MBSs and superconductivity is through vortices.
There are multiple observations of ZBCP trapped in vortex cores of
superconductors Fe(Tex,Se${}_{1-x})$ Wang  _et al._ (2018); Zhang  _et al._
(2018); Machida  _et al._ (2019), (Li0.84Fe0.16)OHFeSe Liu  _et al._ (2018)
and CaKFe2As4 Liu  _et al._ (2020a). These experiments are interpreted as
evidence for MBS. However, there are several difficulties in clearly
distinguishing the MBSs, for e.g. due to the presence of topologically trivial
Caroli-de Gennes- Matricon (CdGM) Caroli  _et al._ (1964) states at the vortex
cores. The CdGM states are observed to produce a broad peak centering at the
zero energy, as the energy separation between the states is of the order of
$\Delta^{2}/E_{F}$, where $E_{F}$ is the Fermi energy. Thus as the energy
separation can be significantly smaller than the instrumental energy
resolution, they appear as a broad peak at the zero energy, being different in
origin from that of the MBS. Moreover, scanning-tunneling microscopy (STM)
experiments Machida  _et al._ (2019); Kong  _et al._ (2019) in FeSe0.45Te0.55
have reported the presence of the ZBCP only in a fraction of the vortex cores
present in the system, depending on the magnetic field, rendering the MBSs
questionable.

###  II.2 Magnetic point defects

Quantum anomalous vortices Jiang  _et al._ (2019) may nucleate at magnetic
impurities at zero magnetic field. In the presence of a surface topological
superconductivity they can support MBSs at the vortex centers, manifesting
themselves as ZBCP in tunneling experiments. These types of sharp zero energy
peaks located at interstitial iron impurities (IFI) have been observed by
(STM) in the superconducting state of Fe1+x(Te,Se) Yin  _et al._ (2015) with
superconducting gap of about 2 meV, below superconducting temperature
$T_{c}=14.5$ K. The zero bias peak intensity decays exponentially from the
center of the peak with correlation length $\xi\approx 3.5$ Å, which is much
smaller than the superconducting correlation length $\sim 25$ Å. The STM have
found IFIs on the exposed $(Te,Se)$ surface located right at the mid position
between four neighboring (Te,Se) atoms. They manifest themselves as sharp
peaks at zero bias; the peaks remain robust in applied magnetic field up to 8
T. Similar ZBCP have also been found later in LiFeAs Zhang  _et al._ (2020)
and in IFI on monolayers of FeSe and FeSe0.5Te0.5 on /SrTiO3(STO) Liu  _et
al._ (2020b). One potential problem is that the peak has an intrinsic width in
all of the systems [for e.g., in Fe1+x(Te,Se) Yin  _et al._ (2015) the width
$\sim 0.6$ meV at $T=1.5$ K] whose origin is unclear. Perhaps, this widening
can be taken into account phenomenologically as a result of a coupling of the
MBS to the bath. In any case, it will likely to have a damaging effect on TKE.
Additionally, in a STM experiment Fan  _et al._ (2021) in FeSe0.45Te0.55 with
IFI, the ZBCP was found only near some iron adatoms, the rest of them having
Yu-Shiva- Rusinov (YSR) bound states at finite energies.

###  II.3 Line defects

The STM study Chen  _et al._ (2020) reports a discovery of zero-energy bound
states simultaneously appearing at both ends of a one-dimensional atomic line
defect in monolayer iron-based high-temperature superconductor FeTe0.5Se0.5
films grown on SrTiO3(001) substrates , which has $T_{c}\approx 62$ K. These
line defects naturally emerge during the growth process and correspond to
lines of missing Te/Se atoms.

Figure 1: A schematic representation of a mesoscopic setup for the formation
of a paired Majorana bound states: The device contains a slab of Fe(Se,Te)
superconductor and a ferromagnetic film (FM) of width W deposited on the
surface of the superconductor, with W being much smaller than the size of the
superconductor. This effectively gives rise to a superconductor-ferromagnet-
superconductor (S-F-S) junction.

Another recent STM Wang  _et al._ (2020) study has found evidence of Majorana
fermions in a certain type of crystalline domain walls associated to the half-
unit cell shift of the Se atom in superconducting FeSe0.45Te0.55 [Fe(Se,Te)].
It was established in Zaki  _et al._ (2021), this material also develops
surface ferromagnetism. The first principle calculations Song  _et al._ (2022)
have suggested that these types of crystalline domain walls can develop an in-
plane ferromagnetism along the domain wall orientation and may support
Majorana modes. Moreover, it has been also suggested, that if the
magnetization can be tuned, MBS can be trapped by defects such as
ferromagnetic domain wall Song  _et al._ (2022).

Motivated by these, we consider a model device consisting of a mesoscopic
Fe(Se,Te) and a ferromagnetic film, which can be used to study TKE with a
higher Topological Kondo temperature scale, thus providing advantage of
previously proposed setup based on complex heterostructures.

Below we describe in detail the setup of a device giving rise to a pair of
MBS’s. It is composed of a slab of Fe(Se,Te) material and a thin ferromagnetic
film (FM) deposited on the surface of the Fe(Se,Te) or any similar material as
shown in the Fig.1. Below the Tc, one part of the device becomes
superconducting, in the thin film area, however, the ferromagnetism persists.
Hence, the system can be simply modeled as a superconductor-ferromagnet-
superconductor (S-F-S) junction [see Fig.1]. In the Fig. 1, the x-axis points
to the normal to the junction between the ferromagnet and the superconductor,
where as the y-axis points parallel to the orientation of the junction. The
semi-infinite superconducting regions occupy intervals $x<-W/2$ and $x>W/2$,
while the FM region occupies interval $-W/2<x<W/2$. We consider that the
system size along y-direction as infinite, the width W remains finite. The
Fe(Se,Te) is known to have non-trivial topological spin-helical Dirac surface
states Zhang  _et al._ (2018) as well as the Rashba type spin-orbit coupling
Christensen  _et al._ (2019); Mascot  _et al._ (2022). The bulk of Fe(Se,Te)
possesses a node-less and almost isotropic s-wave type superconducting gap
Sarkar  _et al._ (2017). The FM produces an effective Zeeman field $h_{y}$
parallel to the junction.

Consequently, the surface S-F-S junction of width W is described by the
following Dirac-Bogolyubov-de-Gennes (BdG) Hamiltonian Song  _et al._ (2022)
in the Nambu basis
$\Psi=(\psi^{\dagger}_{k,\uparrow},\psi^{\dagger}_{k,\downarrow},\psi_{-k,\uparrow},\psi_{-k,\downarrow})$
:

$\displaystyle H_{BdG}$
$\displaystyle=\begin{pmatrix}H_{0}(k)&i\sigma_{y}\Delta(r)\\\
-i\sigma_{y}\Delta^{*}(r)&-H_{0}^{*}(-k).\end{pmatrix}$ (1)

In Eqn.(1), $H_{0}(k)$ is a $2\times 2$ matrix,
$H_{0}(k)=H_{SO}-\mu+H_{Z}(x)$. The $H_{0}(k)$ contains the Rashba spin-orbit
coupling term $H_{SO}=\alpha({\bm{\sigma}}\times{\bm{k}}).\hat{z}$, associated
to the Dirac surface states of the Fe(Se,Te), a chemical potential $\mu$ and a
Zeeman energy term due to the FM given by
$H_{Z}(x)=\sigma_{y}h_{y}\Theta[(W/2)-|x|]$. Here, $\alpha=v_{F}\hbar$, is the
coupling strength, and $v_{F}$ is the Fermi velocity. $\sigma$ are the Pauli
matrices in the spin- space. With $\hbar=1$, the Rashba term $H_{SO}$ becomes
$v_{F}(\sigma_{x}k_{y}+i\sigma_{y}\partial_{x})$, where $k_{y}$ is momentum
along the y-direction that is conserved in the system. The spatial variation
of the superconducting gap function on the surface satisfies
$\Delta(r)=\Delta_{0}e^{isgn(x)(\theta/2)}\Theta[|x|-W/2]$, where $\theta$ is
the phase difference between the left and the right superconducting regime.

It has been shown in Song  _et al._ (2022), that the S-F-S junction described
by the Hamiltonian Eqn.(1) can support one dimensional counter propagating
Majorana modes for a constant phase difference ($\Delta\phi=\pi$) between the
two superconducting regimes on the left and right arising due to the in-plane
magnetic moment in the FM. Moreover, the localized MBSs are formed Fu and Kane
(2008)from the coupling between these dispersive Majorana modes, if the phase
difference $\Delta\phi(y)$ varies continuously along the orientation $y$ of
the junction such that, $\Delta\phi(y)$ becomes an odd integer multiple of
$\pi$. This follows from the analysis of the low energy effective Hamiltonian
derived from the BdG Hamiltonian Eqn.(1), by projecting it onto the basis of
the one-dimensional Majorana modes Song  _et al._ (2022):

$\displaystyle H_{eff}$
$\displaystyle=v_{m}k_{y}\tau_{y}-\Delta_{0}\cos\left[\frac{\Delta\phi(y)}{2}\right]\tau_{z},$
(2)

where, $v_{m}=v_{F}\left[\cos k_{F}W+\frac{\Delta_{0}}{\mu}\sin
k_{F}W\right]\frac{\Delta^{2}_{0}}{(\mu^{2}+\Delta^{2}_{0})}$,
$\Delta\phi(y)=\frac{2h_{y}W}{v_{F}}$ and $\tau$ are the Pauli matrices in the
particle-hole space. Note that the Eqn.(2) has the form of a Su-Schrieffer-
Heeger model. It is also well-known from the Jackiw-Rabbi problem, that a zero
energy localized MBS forms when the mass term
$\Delta_{0}\cos(\frac{\Delta\phi(y)}{2})$ in Eqn.(2) smoothly changes sign or
$\Delta\phi(y)$ becomes an odd integer multiple of $\pi$. The Bogoliubov
quasiparticle operator associated with this state is a Majorana fermion, which
satisfies $\gamma_{0}=\gamma^{\dagger}_{0}$.

Therefore in this setup, to tune the phase difference $\Delta\phi(y)$ for
obtaining a pair of MBSs, either the width $W$ of the FM can be changed along
the $y$ direction or a soft magnetic material for the FM can be used where one
can vary the magnetization strength $h_{y}$.

##  III A $FeSC$ -based device for observation of Topological Kondo effect

Figure 2: A schematic representation of a mesoscopic setup for Topological
Kondo effect with $N=4$ MBSs and $M=3$ metallic leads: The device contains a
Coulomb blockaded superconducting island made of Fe(Se,Te) superconductor and
FM. The island is connected to the ground through a capacitor and hence
possesses a charging energy $E_{c}$. 3 MBSs are tunnel coupled to the metallic
leads with coupling constants $t_{j}$’s (j= 1,2,3).

The MCB device considered in our work for the realization of the Topological
Kondo effect (TKE) comprises of a mesoscopic Fe(Se,Te) (or a similar material)
and two ferromagnetic films deposited on it and resulting in a pair of MBS in
each of the ferromagnetic films as shown in Fig. 2. Realization of the TKE
suggested in the Ref.Béri and Cooper (2012) is based on the MCB which was
described, for example, in Liu  _et al._ (2019), containing $N\geq 4$ MZMs.
Following the suggestions, we consider a floating mesoscopic MCB island,
grounded through a capacitor. To get the critical ground state, one needs to
connect MZMs to just three external metallic leads (this is a minimal number)
via tunneling contacts as schematically depicted in Fig. 2. The charging
energy $E_{C}=\frac{e^{2}}{2C}$, with $C$ being the geometric capacitance of
the box, contributes to the Hamiltonian of the MCB island as

$H_{c}=E_{c}(2N_{c}-Q_{0})^{2};~{}~{}2N_{c}=i\partial/\partial\phi.$ (3)

In Eqn.(3), $N_{c}$ is the number of Cooper pairs in the island, $\phi$ is the
phase of the superconducting order parameter, and $Q_{0}$ is the background
charge determined by the voltage across the capacitor, connected to the MCB.

In general, if there are a total N number of MBS in the MCB, there will be
$N/2$ number of zero energy fermionic modes. As the total charge of the island
is fixed by a large charging energy $E_{c}$, to even or odd number, thus gives
a parity constraint to the MBS operators. This leads to a ground state of the
system with $2^{(N/2)-1}$ fold degeneracy. In this work, we consider the
arrangement where MCB island contains four Majorana zero modes $N=4$, hence
the ground state is two-fold degenerate. The two degenerate quantum ground
states, $|\downarrow\rangle$ ($|\uparrow\rangle$), have $N_{0}$ ($N_{0}-2$)
particles in the condensate and empty (filled) pairs of Majorana modes, thus
encoding a qubit Fu (2010); Béri and Cooper (2012). It is this topological
degeneracy which will lead to the effective spin degeneracy for the effective
Kondo model. Now, we assume, no overlap between the MBS, hence the degeneracy
remains exact, otherwise a “Zeeman coupling” Altland  _et al._ (2014) term can
arise.

The TKE occurs when the MCB is coupled to conduction electrons through
external metallic leads as shown in Fig. 2. We work at the energy scales much
smaller than the superconducting gap, as well as below the energies of any
other subgap excitations of non-Majorana character. In this case the
Hamiltonian of the system is given by $H_{TK}=H_{cond}+H_{c}+H_{tun}$, where
$H_{cond}$ is the Hamiltonian of the conduction electrons of the leads. The
tunneling Hamiltonian $H_{tun}$ describes the low-energy coupling between the
leads and the MCB island. $H_{tun}$ is given by

$\displaystyle
H_{tun}=\exp(\mbox{i}\phi/2)\sum_{j=1}^{3}t_{jj}\gamma_{j}\psi_{j}+\mbox{H.c.},$
(4)

where $\phi$ is the phase of the superconducting order parameter, $t_{jj}$ is
the tunneling amplitude, $\psi_{j}$ is the electron annihilation operator and
$\gamma_{j}$ is the Majorana operator at lead j. Such tunneling process
explicitly excludes the possibility of exciting quasiparticles Plugge  _et
al._ (2016). As was demonstrated in Béri and Cooper (2012); Galpin  _et al._
(2014), if the charging energy of the box is much greater than the
characteristic value of the tunneling matrix elements, i.e. if $E_{C}>>t_{j}$,
one can integrate out the phase fluctuations of the condensate which results
in the exchange interaction between such MCB and the electrons of the leads
(it is supposed that they are spin polarized):

$\displaystyle
H_{ex}=J_{K}^{ij}(\psi^{+}_{i}\psi_{j}-\psi^{+}_{j}\psi_{i})\gamma_{i}\gamma_{j},~{}~{}J_{K}^{ij}\sim
t_{i}t_{j}/E_{C},$ (5)

where indices $i,j$ correspond to the leads. For a single MCB this interaction
gives rise to the TKE, where the leads serve as the bulk. The spin operator
$S^{i}$ realized by the Majorana operators is

$\displaystyle S^{i}=\frac{\rm
i}{2}\epsilon_{ijk}\gamma_{j}\gamma_{k},~{}~{}\\{\gamma_{k},\gamma_{j}\\}=\delta_{jk}.$
(6)

The model comprised of the lead Hamiltonian and the effective exchange
Hamiltonian [Eqn.(5)] defines an antiferromagnetic Kondo problem.

The Topological Kondo temperature $T_{K}$ at which the perturbative approach
breaks down due to the RG flow towards strong coupling is given by Béri and
Cooper (2012); Altland and Egger (2013); Altland  _et al._ (2014),

$\displaystyle T_{K}\sim E_{c}\exp(-1/\rho\bar{J}_{K}),$ (7)

where $\rho$ is the density of states of the lead electrons at the Fermi
energy and $\bar{J}_{K}$ is the average value of the exchange coupling.
Consequently, $T_{K}$ sets an energy scale between the trivial regime and the
Topological Kondo regime, below which a robust non-Fermi liquid (NFL) behavior
can be observed.

One problem is that in all cases apart from tunneling into MBS (4) there may
be a direct tunneling of bulk electrons into the superconductor. There are,
however, reasons why this process is irrelevant. Indeed, the corresponding
Hamiltonian is

$\displaystyle
H_{tun,2}=g\Big{[}A_{lj}\exp(\mbox{i}\phi)\psi_{l}\psi_{j}+H.c.\Big{]},~{}~{}A_{lj}=-A_{jl}.$
(8)

When we integrate over $\phi$, however, the result is a local four fermion
interaction which is highly irrelevant.

Below we give estimates of the parameters pertinent to the realization of TKE
for the device setup depicted in Fig.2.

We start by summarizing the criteria to be satisfied are: (i) four well-
localized MBS, (ii) the MBS must be well-separated, such that the
hybridization between them is insignificant, (iii) a charging energy $E_{c}$
of the order of the superconducting gap $\Delta$. We first evaluate the decay
length ($\lambda$) of the MBS wave function inside the ferromagnetic film
region. The decay length is given by $\lambda=\frac{v_{m}}{\Delta_{0}}$ Song
_et al._ (2022). For FeSe0.45Te0.55, with typical values of $v_{F}\sim 216$
meV ÅWang  _et al._ (2018), $\Delta_{0}\sim 1.8$ meV, $\mu\sim 5\Delta_{0}$
($\mu>>\Delta_{0}$), and considering width of the ferromagnetic region $W\sim
10~{}$nm, the decay length is found to be $\lambda\sim 5$ Å.

With a superconducting island of area $A\sim 60$ nm $\times 60$ nm and
thickness d $\sim 1$ nm, and assuming permittivity of an insulating material
being $\epsilon\sim 2\times 55e^{2}ev^{-1}(\mu m)^{-1}$, the charging energy
$E_{c}=\frac{e^{2}}{2C}$, with $C=\frac{\epsilon A}{d}$ will be of the order
of 1.38 meV $\equiv 16K\sim\Delta_{0}$. Therefore, we see that a mesoscopic
MCB island with well-separated and localized MBS can be achieved with a large
$T_{K}$ [Eqn.(7)]. As the $T_{K}$ achieved in the Fe(Se,Te) based device can
be significantly higher than the heterostructures systems based on
semiconductor nanowire and conventional superconductors Béri and Cooper
(2012), this device will significantly improve the experimental feasibility of
the TKE.

As we have mentioned in the introduction, that one possible way to detect the
TKE is to measure the conductance between different leads, as the TKE leads to
unusual temperature dependencies for the two-terminal linear conductance
$\sigma_{ij}$ ($i\neq j$)Altland and Egger (2013); Altland  _et al._ (2014).
In the temperature-regime larger than $T_{K}$, the $\sigma_{ij}$ grows as
$\frac{1}{ln^{2}(T/T_{K})}$ with lowering the temperature. For temperature-
regime smaller than $T_{K}$, $\sigma_{ij}$ approaches the universal value
$\frac{2e^{2}}{3h}$ with $(T/T_{K})^{2/3}$ temperature dependence. Both the
universal value and the temperature dependence of the conductance are unique
characteristics of the NFL ground state providing an alternative pathway for
detecting the Majorana fermions in the Fe(Se,Te) based device.

##  IV Arrays of TKE

Figure 3: A schematic depiction of the arrangement for the Chiral Kondo
effect. The triangle is the superconducting (Cooper pair) box, the red dots
are Majorana zero modes coupled by tunneling to chiral edges of topological
insulators.

We conclude our paper with a brief description of a device which can
incorporate multiple anyons. This is chiral Kondo lattice first described in
Lopes  _et al._ (2020); Gabay  _et al._ (2022); Lotem  _et al._ (2022).Here
the itinerant electrons are chiral, they move in one direction and hence such
lattice does not have Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction. In
fact, it is equivalent to the bunch of impurities and as a consequence the
ground state remains critical as for the single impurity case.

Chiral Kondo lattice for Topological Kondo effect can serve as a platform for
anyon-based quantum computation Lopes  _et al._ (2020); Gabay  _et al._
(2022); Lotem  _et al._ (2022). An elementary block of such Kondo lattice is
depicted on Fig.(3). The suitable material for chiral edges containing
polarized electrons is MnBi2Te4 Ying  _et al._ (2022), which is distinguished
by a relatively large bulk gap, allowing one to operate at temperatures of
order of several degrees.

##  V Summary

In this work, we have provided a new potential route to achieve the
Topological Kondo effect (TKE) through various iron-based superconductors such
as FeSe0.45Te0.55 [Fe(Se,Te)] based Majorana-Cooper boxes (MCB). The proposed
device set up has two marked differences from the earlier proposals set ups
which include heterostructures containing semiconductor nanowires proximitized
to conventional superconductors. First, there is a significant enhancement of
the Topological Kondo temperature scale and the second, there is a
simplification of the device structure.

##  VI Acknowledgements

This work was supported by U.S. Department of Energy (DOE) the Office of Basic
Energy Sciences, Materials Sciences and Engineering Division under Contract
No. DE-SC0012704.

## References

  * Nayak  _et al._ (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Nozieres and Blandin (1980) P. Nozieres and A. Blandin, Journal de Physique 41, 193 (1980).
  * Tsvelick and Wiegmann (1984) A. Tsvelick and P. Wiegmann, Zeitschrift für Physik B Condensed Matter 54, 201 (1984).
  * Affleck and Ludwig (1991) I. Affleck and A. W. Ludwig, Nuclear Physics B 352, 849 (1991).
  * Potok  _et al._ (2007) R. Potok, I. Rau, H. Shtrikman, Y. Oreg, and D. Goldhaber-Gordon, Nature 446, 167 (2007).
  * Keller  _et al._ (2015) A. Keller, L. Peeters, C. Moca, I. Weymann, D. Mahalu, V. Umansky, G. Zaránd, and D. Goldhaber-Gordon, Nature 526, 237 (2015).
  * Iftikhar  _et al._ (2015) Z. Iftikhar, S. Jezouin, A. Anthore, U. Gennser, F. Parmentier, A. Cavanna, and F. Pierre, Nature 526, 233 (2015).
  * Iftikhar  _et al._ (2018) Z. Iftikhar, A. Anthore, A. Mitchell, F. Parmentier, U. Gennser, A. Ouerghi, A. Cavanna, C. Mora, P. Simon, and F. Pierre, Science 360, 1315 (2018).
  * Karki  _et al._ (2023) D. Karki, E. Boulat, W. Pouse, D. Goldhaber-Gordon, A. K. Mitchell, and C. Mora, Physical Review Letters 130, 146201 (2023).
  * Lopes  _et al._ (2020) P. L. S. Lopes, I. Affleck, and E. Sela, Phys. Rev. B 101, 085141 (2020).
  * Komijani (2020) Y. Komijani, Phys. Rev. B 101, 235131 (2020).
  * Ge and Komijani (2022) Y. Ge and Y. Komijani, Phys. Rev. Lett. 129, 077202 (2022).
  * Béri and Cooper (2012) B. Béri and N. R. Cooper, Phys. Rev. Lett. 109, 156803 (2012).
  * Altland and Egger (2013) A. Altland and R. Egger, Phys. Rev. Lett. 110, 196401 (2013).
  * Altland  _et al._ (2014) A. Altland, B. Béri, R. Egger, and A. M. Tsvelik, Phys. Rev. Lett. 113, 076401 (2014).
  * Väyrynen  _et al._ (2020) J. I. Väyrynen, A. E. Feiguin, and R. M. Lutchyn, Phys. Rev. Res. 2, 043228 (2020).
  * Liu  _et al._ (2021) D. Liu, Z. Cao, X. Liu, H. Zhang, and D. E. Liu, Phys. Rev. B 104, 205125 (2021).
  * Beenakker (2013) C. Beenakker, Annu. Rev. Condens. Matter Phys. 4, 113 (2013).
  * Sato and Ando (2017) M. Sato and Y. Ando, Reports on Progress in Physics 80, 076501 (2017).
  * Xu  _et al._ (2016) G. Xu, B. Lian, P. Tang, X.-L. Qi, and S.-C. Zhang, Phys. Rev. Lett. 117, 047001 (2016).
  * Wang  _et al._ (2018) D. Wang, L. Kong, P. Fan, H. Chen, S. Zhu, W. Liu, L. Cao, Y. Sun, S. Du, J. Schneeloch, R. Zhong, G. Gu, L. Fu, H. Ding, and H.-J. Gao, Science 362, 333 (2018).
  * Zhang  _et al._ (2018) P. Zhang, K. Yaji, T. Hashimoto, Y. Ota, T. Kondo, K. Okazaki, Z. Wang, J. Wen, G. Gu, H. Ding, _et al._ , Science 360, 182 (2018).
  * Machida  _et al._ (2019) T. Machida, Y. Sun, S. Pyon, S. Takeda, Y. Kohsaka, T. Hanaguri, T. Sasagawa, and T. Tamegai, Nature materials 18, 811 (2019).
  * Liu  _et al._ (2018) Q. Liu, C. Chen, T. Zhang, R. Peng, Y.-J. Yan, C.-H.-P. Wen, X. Lou, Y.-L. Huang, J.-P. Tian, X.-L. Dong, G.-W. Wang, W.-C. Bao, Q.-H. Wang, Z.-P. Yin, Z.-X. Zhao, and D.-L. Feng, Phys. Rev. X 8, 041056 (2018).
  * Liu  _et al._ (2020a) W. Liu, L. Cao, S. Zhu, L. Kong, G. Wang, M. Papaj, P. Zhang, Y.-B. Liu, H. Chen, G. Li, _et al._ , Nature Communications 11, 5688 (2020a).
  * Yin  _et al._ (2015) J.-X. Yin, Z. Wu, J. Wang, Z. Ye, J. Gong, X. Hou, L. Shan, A. Li, X. Liang, X. Wu, _et al._ , Nature Physics 11, 543 (2015).
  * Zhu  _et al._ (2020) S. Zhu, L. Kong, L. Cao, H. Chen, M. Papaj, S. Du, Y. Xing, W. Liu, D. Wang, C. Shen, _et al._ , Science 367, 189 (2020).
  * Alicea (2012) J. Alicea, Reports on progress in physics 75, 076501 (2012).
  * Caroli  _et al._ (1964) C. Caroli, P. De Gennes, and J. Matricon, Physics Letters 9, 307 (1964).
  * Kong  _et al._ (2019) L. Kong, S. Zhu, M. Papaj, H. Chen, L. Cao, H. Isobe, Y. Xing, W. Liu, D. Wang, P. Fan, _et al._ , Nature Physics 15, 1181 (2019).
  * Jiang  _et al._ (2019) K. Jiang, X. Dai, and Z. Wang, Phys. Rev. X 9, 011033 (2019).
  * Zhang  _et al._ (2020) S. S. Zhang, J.-X. Yin, G. Dai, L. Zhao, T.-R. Chang, N. Shumiya, K. Jiang, H. Zheng, G. Bian, D. Multer, M. Litskevich, G. Chang, I. Belopolski, T. A. Cochran, X. Wu, D. Wu, J. Luo, G. Chen, H. Lin, F.-C. Chou, X. Wang, C. Jin, R. Sankar, Z. Wang, and M. Z. Hasan, Phys. Rev. B 101, 100507 (2020).
  * Liu  _et al._ (2020b) C. Liu, C. Chen, X. Liu, Z. Wang, Y. Liu, S. Ye, Z. Wang, J. Hu, and J. Wang, Science advances 6, eaax7547 (2020b).
  * Fan  _et al._ (2021) P. Fan, F. Yang, G. Qian, H. Chen, Y.-Y. Zhang, G. Li, Z. Huang, Y. Xing, L. Kong, W. Liu, _et al._ , Nature communications 12, 1348 (2021).
  * Chen  _et al._ (2020) C. Chen, K. Jiang, Y. Zhang, C. Liu, Y. Liu, Z. Wang, and J. Wang, Nature Physics 16, 536 (2020).
  * Wang  _et al._ (2020) Z. Wang, J. O. Rodriguez, L. Jiao, S. Howard, M. Graham, G. Gu, T. L. Hughes, D. K. Morr, and V. Madhavan, Science 367, 104 (2020).
  * Zaki  _et al._ (2021) N. Zaki, G. Gu, A. Tsvelik, C. Wu, and P. D. Johnson, Proceedings of the National Academy of Sciences 118, e2007241118 (2021).
  * Song  _et al._ (2022) R. Song, P. Zhang, and N. Hao, Phys. Rev. Lett. 128, 016402 (2022).
  * Christensen  _et al._ (2019) M. H. Christensen, J. Kang, and R. M. Fernandes, Phys. Rev. B 100, 014512 (2019).
  * Mascot  _et al._ (2022) E. Mascot, S. Cocklin, M. Graham, M. Mashkoori, S. Rachel, and D. K. Morr, Communications Physics 5, 188 (2022).
  * Sarkar  _et al._ (2017) S. Sarkar, J. Van Dyke, P. O. Sprau, F. Massee, U. Welp, W.-K. Kwok, J. C. S. Davis, and D. K. Morr, Phys. Rev. B 96, 060504 (2017).
  * Fu and Kane (2008) L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Liu  _et al._ (2019) C.-X. Liu, D. E. Liu, F.-C. Zhang, and C.-K. Chiu, Phys. Rev. Appl. 12, 054035 (2019).
  * Fu (2010) L. Fu, Phys. Rev. Lett. 104, 056402 (2010).
  * Plugge  _et al._ (2016) S. Plugge, A. Zazunov, E. Eriksson, A. M. Tsvelik, and R. Egger, Phys. Rev. B 93, 104524 (2016).
  * Galpin  _et al._ (2014) M. R. Galpin, A. K. Mitchell, J. Temaismithi, D. E. Logan, B. Béri, and N. R. Cooper, Phys. Rev. B 89, 045143 (2014).
  * Gabay  _et al._ (2022) D. Gabay, C. Han, P. L. S. Lopes, I. Affleck, and E. Sela, Phys. Rev. B 105, 035151 (2022).
  * Lotem  _et al._ (2022) M. Lotem, E. Sela, and M. Goldstein, Phys. Rev. Lett. 129, 227703 (2022).
  * Ying  _et al._ (2022) Z. Ying, S. Zhang, B. Chen, B. Jia, F. Fei, M. Zhang, H. Zhang, X. Wang, and F. Song, Phys. Rev. B 105, 085412 (2022).