# Robust MRI Reconstruction by Smoothed Unrolling (SMUG)

Shijun Liang, , Van Hoang Minh Nguyen, Jinghan Jia, , Ismail Alkhouri, , Sijia
Liu, , Saiprasad Ravishankar S. Liang (corresponding author<EMAIL_ADDRESS>is with the Biomedical Engineering (BME) Department at Michigan State
University (MSU), East Lansing, MI, 48824, USA. J.Jia<EMAIL_ADDRESS>is
with the Computer Science and Engineering (CSE) Department at MSU. M. Nguyen
<EMAIL_ADDRESS>is with the Mathematics Department at MSU. I. Alkhouri
<EMAIL_ADDRESS>&<EMAIL_ADDRESS>is with the Computational Mathematics,
Science & Engineering (CMSE) Department at MSU and a visiting scholar in the
Electrical Engineering & Computer Science Department at the University of
Michigan, Ann Arbor, MI, 48109, USA. S. Liu<EMAIL_ADDRESS>is with the CSE
Department at MSU. S. Ravishankar<EMAIL_ADDRESS>is with the CMSE & BME
Departments at MSU.

###### Abstract

As the popularity of deep learning (DL) in the field of magnetic resonance
imaging (MRI) continues to rise, recent research has indicated that DL-based
MRI reconstruction models might be excessively sensitive to minor input
disturbances, including worst-case additive perturbations. This sensitivity
often leads to unstable, aliased images. This raises the question of how to
devise DL techniques for MRI reconstruction that can be robust to train-test
variations. To address this problem, we propose a novel image reconstruction
framework, termed Smoothed Unrolling (SMUG), which advances a deep unrolling-
based MRI reconstruction model using a randomized smoothing (RS)-based robust
learning approach. RS, which improves the tolerance of a model against input
noises, has been widely used in the design of adversarial defense approaches
for image classification tasks. Yet, we find that the conventional design that
applies RS to the entire DL-based MRI model is ineffective. In this paper, we
show that SMUG and its variants address the above issue by customizing the RS
process based on the unrolling architecture of a DL-based MRI reconstruction
model. Compared to the vanilla RS approach, we show that SMUG improves the
robustness of MRI reconstruction with respect to a diverse set of instability
sources, including worst-case and random noise perturbations to input
measurements, varying measurement sampling rates, and different numbers of
unrolling steps. Furthermore, we theoretically analyze the robustness of our
method in the presence of perturbations. Our code is available at
https://github.com/sjames40/SMUG_journal.

###### Index Terms:

Magnetic resonance imaging, machine learning, deep unrolling, robustness,
randomized smoothing, compressed sensing.

##  I Introduction

Magnetic resonance imaging (MRI) is a popular noninvasive imaging modality,
which involves sequential and slow data collection. As such, MRI scans can be
accelerated by collecting limited data. In this case, the process of image
reconstruction requires tackling an ill-posed inverse problem. To deliver
accurate image reconstructions from such limited information, compressed
sensing (CS) [1] has been extensively used. Conventional CS-MRI assumes the
underlying image’s sparsity (in practice, enforced in the wavelet domain [2]
or via total variation [3]). As further improvement to conventional CS,
various learned sparse signal models have been well-studied, such as involving
patch-based synthesis dictionaries [4, 5], or sparsifying transforms [6, 7].
Learned transforms have been shown to offer an efficient and effective
framework for sparse modeling in MRI [8].

Due to the outstanding representation power of convolutional neural networks
(CNNs), they have been applied in single-modality medical imaging synthesis in
recent years [9, 10, 11, 12]. The U-net neural network, presented in [13] and
referenced in numerous studies, is a popular deep CNN for many tasks involving
image processing. They exhibit two key features: the use of a diminishing path
for gathering contextual information, and a symmetric expansion path in the
U-net for precise localization.

Hybrid-domain DL-based image reconstruction methods, such as Model-based
reconstruction using Deep Learned priors (MoDL) [11], are used to enhance
stability and performance by ensuring data consistency in the training and
reconstruction phases. In MR imaging, data consistency layers are often
essential in reconstruction networks to ensure the image agrees with the
measurements model [14, 15]. Various methods such as [16, 17] maintain this
consistency by deep unrolling-based architectures, which mimic a traditional
iterative algorithm and learn the associated regularization parameters. Other
approaches ensure data consistency by applying methods such as denoising
regularization [18] and plug-and-play techniques [19]. Despite their recent
advancements, DL-based MRI reconstruction models are shown to be vulnerable to
tiny changes in the input, shifts in the measurement sampling rate [20, 21],
and varying iteration numbers in unrolling schemes [22]. In such cases, the
resulting images from DL models are of inferior quality which could possibly
lead to inaccurate diagnoses and, consequently, undesirable clinical
consequences.

It is of much importance in medical imaging applications to learn
reconstruction models that are robust to various measurement artifacts and
noise, or scan or data variations at test time. Although there exist numerous
robustification techniques [23, 24, 25, 26] to tackle the instability of DL
models in image classification tasks, methods to enhance the robustness of DL-
based MRI reconstruction models are less explored due to their regression-
based learning targets. Methods such as randomized smoothing (RS) and its
variations [25, 26, 27], are often used in image classification. They diverge
from traditional defense methods [23, 24] such as adversarial training, which
provide some empirical robustness but are computationally expensive and could
fail under more diverse perturbations. RS ensures the model’s stability within
a radius surrounding the input image [25], which could be critical for medical
use cases such as MRI. Recent early-stage research has begun to apply RS to
DL-based MRI reconstruction in an end-to-end manner [28]. However, the end-to-
end RS approach might not always be an appropriate fit for all image
reconstructors, such as physics-based and hybrid methods.

In our recent conference work [29], we proposed integrating the RS approach
within the MoDL framework for the problem of MR image reconstruction. This is
accomplished by using RS in each unrolling step and at the intermediate
unrolled denoisers in MoDL. This strategy is underpinned by the ‘pre-training
+ fine-tuning’ technique [30, 26]. In [29], we empirically showed that this
approach is effective. In this paper, we provide an analysis and conditions
under which the proposed smoothed unrolling (SMUG) technique is robust against
perturbations. Furthermore, we introduce a novel weighted smoothed unrolling
scheme that learns image-wise weights during smoothing unlike conventional RS.
This approach improves the reconstruction performance. In this work, we
evaluate worst-case additive perturbations in k-space, in contrast to [29],
where image-space perturbations were considered. Moreover, this paper includes
additional comprehensive experimental results.

###  I-A Contributions

The main contributions of this work are as follows:

  * •

We propose SMUG that systematically integrates RS within a deep unrolled
procedure (MoDL).

  * •

We provide a theoretical analysis to demonstrate the robustness of SMUG.

  * •

We enhance the performance of SMUG by introducing weighted smoothing as an
improvement over conventional RS and showcase the resulting gain.

  * •

We compare our methods with the vanilla MoDL [11] and end-to-end RS [28].
Extensive experiments demonstrate the significant improvement of our approach
against various types of reconstruction instabilities.

###  I-B Paper Organization

The remainder of the paper is organized as follows. In Section II, we present
preliminaries and the problem statement. Our proposed method is described in
Section III. Section IV presents experimental results and comparisons, and we
conclude in Section V.

##  II Preliminaries and Problem Statement

Setup of MRI reconstruction. Many medical imaging approaches involve ill-posed
inverse problems such as the work in [31], where the aim is to reconstruct the
original signal $\mathbf{x}\in\mathbb{C}^{q}$ (vectorized image) from
undersampled k-space (Fourier domain) measurements
$\mathbf{y}\in\mathbb{C}^{p}$ with $p<q$. The imaging system in MRI can be
modeled as a linear system $\mathbf{y}\approx\mathbf{A}\mathbf{x}$, where
$\mathbf{A}$ may take on different forms for single-coil or parallel (multi-
coil) MRI, etc. For example, in the single coil Cartesian MRI acquisition
setting, $\mathbf{A}=\mathbf{M}\mathbf{F}$, where $\mathbf{F}$ is the 2D
discrete Fourier transform and $\mathbf{M}$ is a masking operator that
implements undersampling. With the linear observation model, MRI
reconstruction is often formulated as

$\displaystyle\hat{\mathbf{x}}=\underset{\mathbf{x}}{\arg\min}~{}\|\mathbf{A}\mathbf{x}-\mathbf{y}\|^{2}_{2}+\lambda\mathcal{R}(\mathbf{x}),$
(1)

where $\mathcal{R}(\cdot)$ is a regularization function (e.g., $\ell_{1}$ norm
in the wavelet domain to impose a sparsity prior), and $\lambda>0$ is the
regularization parameter.

MoDL [32] is a recent popular supervised deep learning approach inspired by
the MR image reconstruction optimization problem in (1). MoDL combines a
denoising network with a data-consistency (DC) module in each iteration of an
unrolled architecture. In MoDL, the hand-crafted regularizer, $\mathcal{R}$,
is replaced by a learned network-based prior
$\left\|\mathbf{x}-\mathcal{D}_{\bm{\theta}}(\mathbf{x})\right\|_{2}^{2}$
involving a network $\mathcal{D}_{\bm{\theta}}$. MoDL attempts to optimize
this loss by initializing $\mathbf{x}^{0}=\mathbf{A}^{H}\mathbf{y}$, and then
iterating the following process for a number of unrolling steps indexed by
$n\in\\{0,\dots,N-1\\}$. Specifically, MoDL iterations are given by

$\mathbf{x}^{n+1}=\operatorname*{arg\,min}_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{y}\|_{2}^{2}+\lambda\|\mathbf{x}-\mathcal{D}_{\bm{\theta}}(\mathbf{x}^{n})\|_{2}^{2}.$
(2)

After $N$ iterations, we denote the final output of MoDL as
$\mathbf{x}^{N}=\bm{F}_{\text{MoDL}}(\mathbf{x}^{0})$. The weights of the
denoiser are shared across the $N$ blocks and are learned in an end-to-end
supervised manner [11].

Motivation: Lack of robustness in MoDL. In [20], it was demonstrated that deep
learning-based MRI reconstruction can exhibit instability, when confronted
with subtle, nearly imperceptible input perturbations. These perturbations are
commonly referred to as ‘adversarial perturbations’ and have been extensively
investigated in the context of DL-based image classification tasks, as
outlined in [33]. In the context of MRI, these perturbations represent the
worst-case additive perturbations, which may be used to evaluate method
sensitivity and robustness [20, 34]. Let $\bm{\delta}$ denote a small
perturbation of the measurements that falls in an $\ell_{\infty}$ ball of
radius $\epsilon$, i.e., $\|\bm{\delta}\|_{\infty}\leq\epsilon$. Adversarial
disturbances then correspond to the worst-case input perturbation vector
$\bm{\delta}$ that maximizes the reconstruction error, i.e.,

$\begin{array}[]{ll}\max_{\|\bm{\delta}\|_{\infty}\leq\epsilon}\|\bm{F}_{\text{{{MoDL}}}}(\mathbf{A}^{H}(\mathbf{y}+\bm{\delta}))-\mathbf{t}\|_{2}^{2},\end{array}$
(3)

where $\mathbf{t}$ is a ground truth target image from the training set (i.e.,
label). The operator $\mathbf{A}^{H}$ transforms the measurements $\mathbf{y}$
to the image domain, and $\mathbf{A}^{H}\mathbf{y}$ is the input (aliased)
image to the reconstruction model. The optimization problem in (3) can be
effectively solved using the iterative projected gradient descent (PGD) method
[23].

In Fig. 1-(a) and (b), we show reconstructed images using MoDL originating
from a benign (i.e., undisturbed) input and a PGD scheme-perturbed input,
respectively. It is evident that the worst-case input disturbance
significantly deteriorates the quality of the reconstructed image. While one
focus of this work is to enhance robustness against input perturbations,
Fig.1-(c) and (d) highlight two additional potential sources of instability
that the reconstructor (MoDL) can encounter during testing: variations in the
measurement sampling rate (resulting in “perturbations” to the sparsity of the
sampling mask in $\mathbf{A}$) [20] and changes in the number of unrolling
steps [22]. In scenarios where the sampling mask (Fig.1-(c)) or number of
unrolling steps (Fig.1-(d)) deviate from the settings used during MoDL
training, we observe a significant degradation in performance compared to the
original setup (Fig.1-(a)), even in the absence of additive measurement
perturbations. In Section IV, we demonstrate how our method improves the
reconstruction robustness in the presence of different types of perturbations,
including those in Fig.1.

Figure 1: MoDL’s instabilities resulting from perturbations to input data, the
measurement sampling rate, and the number of unrolling steps used at testing
phase shown on an image from the fastMRI dataset [35]. We refer readers to
Section IV for further details about the experimental settings. (a) MoDL
reconstruction from benign (i.e., without additional noise/perturbation)
measurements with $4\times$ acceleration (i.e., 25% sampling rate) and 8
unrolling steps. (b) MoDL reconstruction from disturbed input with
perturbation strength $\epsilon=0.02$ (see Section IV-A) with other settings
identical to (a). (c) MoDL reconstruction from clean measurements with
$2\times$ acceleration (i.e., 50% sampling), and using 8 unrolling steps. (d)
MoDL reconstruction from clean or unperturbed measurements with $4\times$
acceleration and _16_ unrolling steps.

Randomized smoothing (RS). Randomized Smoothing, introduced in [25], enhances
the robustness of DL models against noisy inputs. It is implemented by
generating multiple randomly modified versions of the input data and
subsequently calculating an averaged output from this diverse set of inputs.

Given some function $f(\mathbf{x})$, RS formally replaces $f$ with a smoothed
version

$g(\mathbf{x})\mathrel{\mathop{:}}=\mathbb{E}_{\bm{\eta}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})}[f(\mathbf{x}+\bm{\eta})]\>,$
(4)

where $\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$ denotes a Gaussian
distribution with zero mean and element-wise variance $\sigma^{2}$, and
$\mathbf{I}$ denotes the identity matrix of appropriate size. Prior research
has shown that RS has been effective as an adversarial defense approach in DL-
based image classification tasks [25, 26, 36]. However, the question of
whether RS can significantly improve the robustness of MoDL and other image
reconstructors has not been thoroughly explored. A preliminary investigation
in this area was conducted by [28], which demonstrated the integration of RS
into MR image reconstruction in an end-to-end (E2E) setting. We can formulate
image reconstruction using RS-E2E as

$\displaystyle\mathbf{x}_{\textsc{RS-E2E}}=\mathbb{E}_{\bm{\eta}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})}[\bm{F}_{\text{{MoDL}}}(\mathbf{A}^{H}\mathbf{y}+\bm{\eta})].$
(RS-E2E)

Note that this formulation is a slight modification to the one used in [28],
where the random noise vector $\bm{\eta}$ is directly added to $\mathbf{y}$.
Here, it is added in the image domain (complex-valued) to
$\mathbf{A}^{H}\mathbf{y}$, which is used both as input to MoDL and within
each data-consistency step in MoDL. We found that this approach is slightly
better in practice.

Figure 2: A schematic overview of RS-E2E.

Fig. 2 shows a block diagram of RS-E2E-backed MoDL. This RS-integrated MoDL is
trained with supervision in the standard manner. Although RS-E2E represents a
straightforward application of RS to MoDL, it remains unclear if the
formulation in (RS-E2E) is the most effective method to incorporate RS into
unrolled algorithms such as MoDL, considering the latter’s specialties, e.g.,
the involved denoising and the data-consistency (DC) steps.

As such, for the rest of the paper, we focus on studying the following
questions (Q1)–(Q4).

(Q1): How should RS be integrated into an unrolled algorithm such as MoDL?
(Q2): How do we learn the network $\mathcal{D}_{\bm{\theta}}(\cdot)$ in the
presence of RS operations? (Q3): Can we prove the robustness of SMUG in the
presence of data perturbations? (Q4): Can we further improve the RS operation
in SMUG for enhanced image quality or sharpness?

##  III Methodology

In this section, we address questions (Q1)–(Q4) by taking the unrolling
characteristics of MoDL into the design of an RS-based MRI reconstruction. The
proposed novel integration of RS with MoDL is termed Smoothed Unrolling
(SMUG). We also explore an extension to SMUG.

###  III-A Solution to (Q1): RS at intermediate unrolled denoisers

As illustrated in Fig.,2, the RS operation in RS-E2E is typically applied to
MoDL in an end-to-end manner. This does not shed light on which component of
MoDL needs to be made more robust. Here, we explore integrating RS at each
intermediate unrolling step of MoDL. In this subsection, we present SMUG,
which applies RS to the denoising network. This seemingly simple modification
is related to a robustness certification technique known as “denoised
smoothing” [26]. In this technique, a smoothed denoiser is used, proving to be
sufficient for establishing robustness in the model. We use
$\mathbf{x}_{\textrm{S}}^{n}$ to capture the $n$-th iteration of SMUG.
Starting from $\mathbf{x}_{\textrm{S}}^{0}=\mathbf{A}^{H}\mathbf{y}$, the
procedure is given by

$\mathbf{x}_{\textrm{S}}^{n+1}=\operatorname*{arg\,min}_{\mathbf{x}}\|\mathbf{A}\mathbf{x}-\mathbf{y}\|_{2}^{2}+\lambda\|\mathbf{x}-\mathbb{E}_{\bm{\eta}}\big{[}\mathcal{D}_{\bm{\theta}}(\mathbf{x}_{\textrm{S}}^{n}+\bm{\eta})\big{]}\|_{2}^{2}\>,$
(5)

where $\bm{\eta}$ is drawn from
$\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$. After $N-1$ iterations, the
final output of SMUG is denoted by
$\mathbf{x}^{N}_{\text{S}}=\bm{F}_{\text{SMUG}}(\mathbf{x}^{0})$. Fig. 3
presents the architecture of SMUG.

Figure 3: Architecture of SMUG.

###  III-B Solution to (Q2): SMUG’s pre-training & fine-tuning

In this subsection, we develop the training scheme of SMUG. Inspired by the
currently celebrated “pre-training + fine-tuning” technique [30, 26], we
propose to train SMUG following this learning paradigm. Our rationale is that
pre-training can provide a robustness-aware initialization of the DL-based
denoising network for fine-tuning. To pre-train the denoising network
$\mathcal{D}_{\bm{\theta}}$, we consider a mean squared error (MSE) loss that
measures the Euclidean distance between images denoised by
$\mathcal{D}_{\bm{\theta}}$ and the target images, denoted by $\mathbf{t}$.
This leads to the pre-training step:

$\bm{\theta}_{\mathrm{pre}}=\displaystyle\operatorname*{arg\,min}_{\bm{\theta}}\mathbb{E}_{\mathbf{t}\in\mathcal{T}}[\mathbb{E}_{\bm{\eta}}||\mathcal{D}_{\bm{\theta}}(\mathbf{t}+\bm{\eta})-\mathbf{t}||_{2}^{2}]\>,$
(6)

where $\mathcal{T}$ is the set of ground truth images in the training dataset.
Next, we develop the fine-tuning scheme to improve
$\bm{\theta}_{\mathrm{pre}}$ based on the labeled/paired MRI dataset. Since RS
in SMUG (Fig. 3) is applied to every unrolling step, we propose an unrolled
stability (UStab) loss for fine-tuning $\mathcal{D}_{\bm{\theta}}$:

$\displaystyle\ell_{\mathrm{UStab}}(\bm{\theta};\mathbf{y},\mathbf{t})=\sum_{n=0}^{N-1}\mathbb{E}_{\bm{\eta}}||\mathcal{D}_{\bm{\theta}}(\mathbf{x}^{n}+\bm{\eta})-\mathcal{D}_{\bm{\theta}}(\mathbf{t})||^{2}_{2}\>.$
(7)

The UStab loss in (7) relies on the target images, bringing in a key benefit:
the denoising stability is guided by the reconstruction accuracy of the
ground-truth image, yielding a graceful trade-off between robustness and
accuracy.

Integrating the UStab loss, defined in (7), with the standard reconstruction
loss, we obtain the fine-tuned $\bm{\theta}$ by minimizing
$\mathbb{E}_{(\mathbf{y},\mathbf{t})}[\ell(\bm{\theta};\mathbf{y},\mathbf{t})]$,
where

$\displaystyle\ell(\bm{\theta};\mathbf{y},\mathbf{t})=\ell_{\mathrm{UStab}}(\bm{\theta};\mathbf{y},\mathbf{t})+\lambda_{\ell}\|\bm{F}_{\text{SMUG}}(\mathbf{A}^{H}\mathbf{y})-\mathbf{t}\|_{2}^{2},$
(8)

with $\lambda_{\ell}>0$ is a regularization parameter to strike a balance
between the reconstruction error (for accuracy) and the denoising stability
(for robustness) terms. We initialize $\bm{\theta}$ as
$\bm{\theta}_{\mathrm{pre}}$ when optimizing (8) using standard optimizers
such as Adam [37].

###  III-C Answer to (Q3): Analyzing the robustness of SMUG in the presence of
data perturbations

The following theorem discusses the robustness (i.e., sensitivity to input
perturbations) achieved with SMUG. Note that all norms on vectors (resp.
matrices) denote the $\ell_{2}$ norm (resp. spectral norm) unless indicated
otherwise.

######  Theorem 1.

Assume the denoiser network’s output is bounded in norm. Given the initial
input image $\mathbf{A}^{H}\mathbf{y}$ obtained from measurements
$\mathbf{y}$, let the SMUG reconstructed image at the $n$-th unrolling step be
$\mathbf{x}^{n}_{\text{S}}(\mathbf{A}^{H}\mathbf{y})$ with RS variance of
$\sigma^{2}$. Let $\bm{\delta}$ denote an additive perturbation to the
measurements $\mathbf{y}$. Then,

$\|\mathbf{x}^{n}_{\text{S}}(\mathbf{A}^{H}\mathbf{y})-\mathbf{x}^{n}_{\text{S}}(\mathbf{A}^{H}(\mathbf{y}+\bm{\delta}))\|\leq
C_{n}\|\bm{\delta}\|,$ (9)

where
$C_{n}=\alpha\|\mathbf{A}\|_{2}\begin{pmatrix}\frac{1-\begin{pmatrix}\frac{M\alpha}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}}{1-\frac{M\alpha}{\sqrt{2\pi}\sigma}}\end{pmatrix}+\|\mathbf{A}\|_{2}\begin{pmatrix}\frac{M\alpha}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}$,
with $\alpha=\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}$ and
$M=2\max_{\mathbf{x}}(\|\mathcal{D}_{\bm{\theta}}(\mathbf{x})\|)$.

The proof is provided in the Appendix. Note that the output of SMUG
$\mathbf{x}^{n}_{\text{S}}(\cdot)$ depends on both the initial input (here
$\mathbf{A}^{H}\mathbf{y}$) and the measurements $\mathbf{y}$. We abbreviated
it to $\mathbf{x}^{n}_{\text{S}}(\mathbf{A}^{H}\mathbf{y})$ in the theorem and
proof for notational simplicity. The constant $C_{n}$ depends on the number of
iterations or unrolling steps $n$ as well as the RS standard deviation
parameter $\sigma$. For large $\sigma$, the robustness error bound for SMUG
clearly decreases as the number of iterations $n$ increases. In particular, if
$\sigma>M\alpha/\sqrt{2\pi}$, then as $n\rightarrow\infty$,
$C_{n}\rightarrow\alpha\left\|\mathbf{A}\right\|_{2}/\begin{pmatrix}1-\frac{M\alpha}{\sqrt{2\pi\sigma}}\end{pmatrix}$.
Furthermore, as $\sigma\to\infty$, $C_{n}\to
C\triangleq\alpha\left\|\mathbf{A}\right\|_{2}$. Clearly, if $\alpha\leq 1$
and $\left\|\mathbf{A}\right\|_{2}\leq 1$ (normalized), then $C\leq 1$. Thus,
for sufficient smoothing, the error introduced in the SMUG output due to input
perturbation never gets worse than the size of the input perturbation. Thus,
the output is stable with respect to perturbations. These results corroborate
experimental results in Section IV on how SMUG is robust (whereas other
methods, such as vanilla MoDL, breakdown) when increasing the number of
unrolling steps at test time, and is also more robust for larger $\sigma$
(with good accuracy-robustness trade-off).

The only assumption in our analysis is that the denoiser network output is
bounded in norm. This consideration is handled readily when the denoiser
network incorporates bounded activation functions such as the sigmoid or
hyperbolic tangent. Alternatively, if we expect image intensities to lie
within a certain range, a simple clipping operation in the network output
would ensure boundedness for the analysis.

A key distinction between SMUG and prior works, such as RS-E2E [28], is that
smoothing is performed in every iteration. Moreover, while [28] assumes the
end-to-end mapping is bounded, in MoDL or SMUG, it clearly isn’t because the
data-consistency step’s output is unbounded as $\mathbf{y}$ grows.

###  III-D Solution to (Q4): Weighted Smoothing

In this subsection, we present a modified formulation of randomized smoothing
to improve its performance in SMUG. Randomized smoothing in practice involves
uniformly averaging images denoised with random perturbations. This can be
viewed as a type of mean filter, which leads to oversmoothing of structural
information in practice. As such, we propose weighted randomized smoothing,
which employs an encoder to assess a weighting (scalar) for each denoised
image and subsequently applies the optimal weightings while aggregating images
to enhance the reconstruction performance. Our method not only surpasses the
SMUG technique but also excels in enhancing image sharpness across various
types of perturbation sources. This allows for a more versatile or flexible
and effective approach for improving image quality under different conditions.

The weighted randomized smoothing operation applied on a function $f(\cdot)$
is as follows:

$g_{\textrm{w}}(\mathbf{x}):=\frac{\mathbb{E}_{\bm{\eta}}[w(\mathbf{x}+\bm{\eta})f(\mathbf{x}+\bm{\eta})]}{\mathbb{E}_{\bm{\eta}}[w(\mathbf{x}+\bm{\eta})]}\>,$
(10)

where $w(\cdot)$ is an input-dependent weighting function.

Based on the weighted smoothing operation in (10), we introduce Weighted SMUG.
This approach involves the application of weighted RS at each denoising step,
and the weighting encoder is trained in conjunction with the denoiser during
the fine-tuning stage. Specifically, in the $n$-th unrolling step, we use a
weighting encoder $\mathcal{E}_{\bm{\phi}}$, parameterized by $\bm{\phi}$, to
learn the weight of each image used for (weighted) averaging. Here, we use
$\mathbf{x}^{n}_{\textrm{W}}$ to denote the output of the $n$-th block.
Initializing $\mathbf{x}^{0}_{\textrm{W}}=\mathbf{A}^{H}\mathbf{y}$, the
output of Weighted SMUG w.r.t. $n$ is

$\begin{split}&\mathbf{x}_{\textrm{W}}^{n+1}=\operatorname*{arg\,min}_{\mathbf{x}}~{}\|\mathbf{A}\mathbf{x}-\mathbf{y}\|_{2}^{2}~{}+\\\
&\lambda\begin{Vmatrix}\mathbf{x}-\frac{\mathbb{E}_{\bm{\eta}}[\mathcal{E}_{\bm{\phi}}(\mathbf{x}_{\textrm{W}}^{n}+\bm{\eta})\mathcal{D}_{\bm{\theta}}(\mathbf{x}_{\textrm{W}}^{n}+\bm{\eta})]}{\mathbb{E}_{\bm{\eta}}[\mathcal{E}_{\bm{\phi}}(\mathbf{x}_{\textrm{W}}^{n}+\bm{\eta})]}\end{Vmatrix}_{2}^{2}\>.\end{split}$
(11)

After $N$ iterations, the final output of Weighted SMUG is
$\mathbf{x}^{N}_{\text{W}}=\bm{F}_{\text{wSMUG}}(\mathbf{x}^{0})$. Figure 4
illustrates the architecture of weighted SMUG.

Figure 4: Architecture of weighted SMUG.

Furthermore, we extend the “pre-training + fine-tuning” approach proposed in
Section III-B to the Weighted SMUG method. In this case, we obtain the fine-
tuned $\bm{\theta}$ and $\bm{\phi}$ by using

$\displaystyle\min_{\bm{\theta},\bm{\phi}}\,\mathbb{E}_{(\mathbf{y},\mathbf{t})}[\lambda_{l}\|{F}_{\text{wSMUG}}(\mathbf{A}^{H}\mathbf{y})-\mathbf{t}\|_{2}^{2}+\ell_{\mathrm{UStab}}(\bm{\theta};\mathbf{y},\mathbf{t})].$
(12)

##  IV Experiments

###  IV-A Experimental Setup

Models & datasets: We evaluate the performance of the proposed methods, SMUG
and Weighted SMUG, and compare the results with other schemes such as the
RS-E2E method. We use the fastMRI dataset including brain and knee scans [35].
For the MoDL architecture, we use the recent state-of-the-art denoising
network Deep iterative Down Network, which consists of $3$ down-up blocks
(DUBs) and 64 channels [38]. Additionally, for MoDL, we use $N=8$ unrolling
steps with denoising regularization parameter $\lambda=1$. The conjugate
gradient method [32], with a tolerance level of $10^{-6}$, is utilized to
execute the DC block. We use binary Cartesian k-space undersampling masks
which correspond to 4x (25.0% sampling) and 8x (12.5% sampling) acceleration
of the k-space acquisition. The coil sensitivity maps for all scenarios were
generated with the BART toolbox [39].

Training: For our study, we execute two experimental cases. For the first
case, we utilize the fastMRI knee dataset, with 32 images for validation and
64 unseen images for testing. In the second case, we employ our method for the
fastMRI brain dataset. We used $3000$ training scans in both cases. We use a
batch size of $2$ and $60$ training epochs. The experiments are run using two
A5000 GPUs. The ADAM optimizer [40] is utilized for training the network
weights with momentum parameters of $(0.5,0.999)$ and learning rate of
$10^{-4}$. The stability parameter $\lambda_{\ell}$ in (8) (and (12)) is tuned
so that the standard accuracy of the learned model is comparable to vanilla
MoDL. For RS-E2E, we set the standard deviation of Gaussian noise to
$\sigma=0.01$, and use $10$ Monte Carlo samplings to implement the smoothing
operation. Note that in our experiments, Gaussian noise and corruptions are
added to real and imaginary parts of the data with the indicated $\sigma$.

Testing: We evaluate our methods on clean data (without additional
perturbations), data with randomly injected noise, and data contaminated with
worst-case additive perturbations. The worst-case disturbances allow us to see
worst-case method sensitivity and are generated by the $\ell_{\infty}$-norm
based PGD scheme with 10 steps [20] corresponding to
$\left\|\bm{\delta}\right\|_{\infty}\leq\epsilon$, where $\epsilon$ is set
nominally as the maximum underlying k-space real and imaginary part magnitude
scaled by $0.05$. We will only indicate the scaling factor for $\epsilon$
(e.g., $0.05$) in the results and plots that follow. The quality of
reconstructed images is measured using peak signal-to-noise ratio (PSNR) and
structure similarity index measure (SSIM)[41]. In addition to the worst-case
perturbations and random noise, we evaluate the performance of our methods in
the presence of additional instability sources such as (i) different
undersampling rates, and (ii) different numbers of unrolling steps.

We have also observed the performance of adversarial training. While it
provides consistent performance when handling worst-case additive
perturbations at test time, we observed that its effectiveness degraded for
other perturbations including modified sampling rates and varied unrolling
steps. Moreover, adversarial training is much slower than SMUG. Importantly,
the proposed SMUG and Weighted SMUG are not trained to be robust to any
specific perturbations or instabilities, but are nevertheless effective for
many such scenarios.

###  IV-B Results

TABLE I:  Performance comparison of different smoothing architectures (RS-E2E, SMUG, Weighted SMUG), together with vanilla MoDL at 4x undersampling. Here, ‘Clean Accuracy’, ‘Robust Accuracy by random noise’, and ‘Robust Accuracy by PGD’ refer to PSNR/SSIM evaluated on benign data, random noise-injected (with Gaussian noise) data, and PGD-enabled worst-case perturbed data, respectively. $\uparrow$ signifies that a higher number indicates a better reconstruction accuracy. The result $a$$\pm b$ represents mean $a$ and standard deviation $b$ over 64 testing brain Fast MRI images. The relative (better or worse) performance of methods w.r.t. metric $a$ is reported with respect to that of vanilla MoDL. Models | Clean Accuracy | Robust Accuracy by random noise | Robust Accuracy by PGD
---|---|---|---
Metrics | PSNR $\uparrow$ | SSIM $\uparrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | PSNR $\uparrow$ | SSIM $\uparrow$
Vanilla MoDL | 32.7$\pm$2.87 | 0.921$\pm$0.07 | 31.77$\pm$2.77 | 0.908$\pm$0.07 | 24.45$\pm$2.42 | 0.767$\pm$0.07
RS-E2E |  -0.03$\pm$3.24 |  -0.001$\pm$0.07 | +0.43$\pm$2.90 | +0.010$\pm$0.07 | +0.86$\pm$2.70 | +0.008$\pm$0.08
SMUG (Ours) | -0.05$\pm$3.01 | -0.001$\pm$0.08 | +0.75$\pm$2.98 | +0.015$\pm$0.08 | +3.98$\pm$2.28 | +0.014$\pm$0.11
Weighted SMUG (Ours) | -0.06$\pm$2.25 | -0.001$\pm$0.07 | +0.87$\pm$2.97 | +0.02$\pm$0.08 | +4.82$\pm$2.2 | +0.034$\pm$0.09
TABLE II: Similar setup as Table I, but the training is with knee dataset images, which include 4x and 8x k-space undersampling and testing on knee dataset images with the same undersampling rate for the comparisons. Ax | Models | Clean Accuracy | Robust Accuracy by random noise | Robust Accuracy by PGD
---|---|---|---|---
Factor | Metrics | PSNR $\uparrow$ | SSIM $\uparrow$ | PSNR $\uparrow$ | SSIM $\uparrow$ | PSNR $\uparrow$ | SSIM $\uparrow$
$\mathrm{4x}$ | Vanilla MoDL | 31.82 $\pm$2.44 | 0.918$\pm$0.08 | 30.67$\pm$2.06 | 0.91$\pm$0.09 | 24.23$\pm$1.82 | 0.764$\pm$0.08
RS-E2E | -0.04 $\pm$2.37 | -0.002 $\pm$0.08 | +0.48 $\pm$2.31 | +0.01$\pm$0.07 | +2.07$\pm$2.17 | +0.012$\pm$0.08
SMUG (Ours) | -0.06 $\pm$2.23 | -0.0039$\pm$0.09 | +0.72$\pm$2.32 | +0.018$\pm$0.08 | +4.31 $\pm$2.22 | +0.027$\pm$0.07
Weighted SMUG (Ours) | -0.07 $\pm$2.26 | -0.0041$\pm$0.11 | +0.82 $\pm$2.34 | +0.03$\pm$0.10 | +4.67$\pm$2.25 | +0.029$\pm$0.10
$\mathrm{8x}$ | Vanilla MoDL | 30.83$\pm$3.14 | 0.905$\pm$0.10 | 28.76$\pm$2.87 | 0.872$\pm$0.09 | 21.28$\pm$2.54 | 0.712$\pm$0.10
RS-E2E | -0.06$\pm$3.08 | -0.002$\pm$0.11 | +0.54$\pm$2.91 | +0.017$\pm$0.09 | +1.97$\pm$2.74 | +0.011$\pm$0.09
SMUG (Ours) | -0.09 $\pm$3.04 | -0.019$\pm$0.12 | +0.79$\pm$2.92 | +0.021$\pm$0.09 | +4.84$\pm$2.44 | +0.03$\pm$0.08
Weighted SMUG (Ours) | -0.08$\pm$3.02 | -0.021$\pm$0.10 | +0.94$\pm$2.93 | +0.024$\pm$0.08 | +5.32$\pm$2.48 | +0.033$\pm$0.09
Ground Truth | Vanilla MoDL | RS-E2E | SMUG | Weighted-SMUG
---|---|---|---|---
|  |  |  |
Figure 5: Visualization of ground truth and reconstructed images using different methods for 4x k-space undersampling, evaluated on PGD-generated worst-case inputs of perturbation strength $\epsilon=0.02$. Ground Truth | Vanilla MoDL | RS-E2E | SMUG | Weighted-SMUG
---|---|---|---|---
|  |  |  |
Figure 6: Visualization of ground-truth and reconstructed images using different methods for 4x k-space undersampling, evaluated on PGD-generated worst-case inputs of perturbation strength $\epsilon=0.02$. Ground Truth | Vanilla MoDL | RS-E2E | SMUG | Weighted-SMUG
---|---|---|---|---
|  |  |  |
Figure 7: Visualization of ground truth and reconstructed images using different methods for 8x k-space undersampling, evaluated on PGD-generated worst-case inputs of perturbation scaling $\epsilon=0.02$. Figure 8: PSNR of baseline methods and the proposed method versus perturbation strength (i.e., scaling) $\epsilon$ used in PGD-generated worst-case examples at testing time with 4x k-space undersampling. $\epsilon=0$ corresponds to clean accuracy.  Figure 9:  PSNR results for different MRI reconstruction methods versus different measurement sampling rates (models trained at $4\times$ acceleration).  Figure 10:  PSNR results for different MRI reconstruction methods at 4x k-space undersampling versus number of unrolling steps ($8$ steps used in training). “Clean” and ”Robust” denote the cases without and with added worst-case (for each method) measurement perturbations.  |
---|---
Figure 11: Left: Norm of difference between SMUG and RS-E2E reconstructions
and the ground truth for different choices of $\sigma$ in the smoothing
process. A worst-case PGD perturbation $\bm{\delta}$ computed at $\sigma=0.01$
was added to the measurements in all cases. Right: Robustness error for SMUG
and RS-E2E at various $\sigma$, i.e., norm of difference between output with
the perturbation $\bm{\delta}$ and without it.  Figure 12: Weights predicted
by the weight encoder network in Weighted SMUG (from final layer of unrolling)
plotted against root mean squared error (RMSE) of the corresponding denoised
images for $5$ randomly selected scans (with 4x undersampling). Figure 13:
PSNR vs. worst-case perturbation strength ($\epsilon)$ for SMUG for different
configurations of UStab loss (7).

Table I presents PSNR and SSIM values for different smoothing architectures
with different training schemes, along with vanilla MoDL as a baseline using
the brain dataset. The evaluation is conducted on the clean, noisy (with added
Gaussian noise), and worst-case perturbed (using PGD for each method)
measurements. As indicated by the PSNR and SSIM values, we observe that our
method, SMUG, outperforms RS-E2E and vanilla MoDL in robust accuracy, i.e.,
image quality with perturbations. This observation is consistent with the
visualization of reconstructed images for the brain dataset in Fig. 5.
Furthermore, we observe that Weighted SMUG outperforms SMUG in robustness.
However, weighted SMUG requires longer time for training, which represents a
trade-off. The clean accuracies of Weighted SMUG and SMUG are similar to
vanilla MoDL indicating a good clean accuracy vs. robustness trade-off.

Table II compares image reconstruction quality with different methods at two
sampling acceleration factors for the knee dataset. We observe quite similar
and congruous outcomes to those reported in Table I. Figs. 6 and 7 show
reconstructed images by different methods for knee scans at 4x and 8x
undersampling, respectively. We observe that SMUG and Weighted SMUG show fewer
artifacts, sharper features, and fewer errors when compared to Vanilla MoDL
reconstruction in the presence of the worst-case perturbations.

Fig. 8 presents PSNR results for the considered models under different scales
of worst-case perturbations (i.e., attack strength $\epsilon$). We used the
knee dataset for this and following results. We observe that SMUG and weighted
SMUG outperform other methods across perturbation strengths.

Next, we evaluate the effectiveness of the proposed MRI reconstruction methods
when faced with sampling rates and unrolling steps that are different from the
settings used in the training phase. Figs. 9 and 10 present the evaluation
results of SMUG and Weighted SMUG along with the two baselines, vanilla MoDL,
and RS-E2E. As observed, SMUG and Weighted SMUG achieve remarkable
improvements in robustness with different sampling rates and unrolling steps,
which MoDL and RS-E2E fail to achieve. Although we do not intentionally design
our method to mitigate MoDL’s instabilities against different sampling rates
and unrolling steps, the SMUG approaches nevertheless provide improved PSNRs
over other baselines. This indicates broader value for the robustification
strategies incorporated in our schemes.

We conduct additional studies on the unrolled stability loss in our scheme to
show the importance of integrating target image denoising into SMUG’s training
pipeline in (7). Fig. 13 presents PSNR values versus perturbation
strength/scaling ($\epsilon$) when using different alternatives to
$\mathcal{D}_{\bm{\theta}}(\mathbf{t})$ in (7), including $\mathbf{t}$ (the
original target image), $\mathcal{D}_{\bm{\theta}}(\mathbf{x}_{n})$ (denoised
output of each unrolling step), and variants when using the fixed, vanilla
MoDL’s denoiser $\mathcal{D}_{\bm{\theta}_{\text{{MoDL}}}}$ instead. As we can
see, the performance of SMUG varies when the UStab loss (7) is configured
differently. The proposed $\mathcal{D}_{\bm{\theta}}(\mathbf{t})$ outperforms
other baselines. A possible reason is that it infuses supervision of target
images in an adaptive, denoising-friendly manner, i.e., taking the influence
of $\mathcal{D}_{\bm{\theta}}$ into consideration.

To comprehensively assess the influence of the introduced noise during
smoothing, denoted as $\bm{\eta}$, on the efficacy of the suggested
approaches, we undertake an experiment involving varying noise standard
deviations $\sigma$. The outcomes, documented in terms of RMSE, are showcased
in Fig.11. The accuracy (reconstruction quality w.r.t. ground truth) and
robustness error (error between with and without measurement perturbation
cases) are shown for both SMUG and RS-E2E. We notice a notable trend: as the
noise level $\sigma$ increases, the accuracy for both methods improves before
beginning to degrade. Importantly, SMUG consistently outperforms end-to-end
smoothing. Furthermore, the robustness error continually drops as $\sigma$
increases (corroborating with our analysis/bound in Section III), with more
rapid decrease for SMUG.

In our final study, we analyze the behavior of the Weighted SMUG algorithm. We
delve into the nuances of weighted smoothing, which can assign different
weights to different images during the smoothing process. The aim is to gauge
how the superior performance of Weighted SMUG arises from the variations in
learned weights. Our findings indicate that among the $10$ Monte Carlo
samplings implemented for the smoothing operation, those with lower denoising
RMSE when compared to the ground truth images generally receive higher
weights, as illustrated in Fig. 12.

##  V Discussion and Conclusion

In this work, we proposed a scheme for improving the robustness of DL-based
MRI reconstruction. In particular, we investigated deep unrolled
reconstruction’s weaknesses in robustness against worst-case or noise-like
additive perturbations, sampling rates, and unrolling steps. To improve the
robustness of the unrolled scheme, we proposed SMUG with a novel unrolled
smoothing loss. We also provided a theoretical analysis on the robustness
achieved by our proposed method. Compared to the vanilla MoDL approach and
other schemes, we empirically showed that our approach is effective and can
significantly improve the robustness of a deep unrolled scheme against a
diverse set of external perturbations. We also further improved SMUG’s
robustness by introducing weighted smoothing as an alternative to conventional
RS, which adaptively weights different images when aggregating them. In future
work, we hope to apply the proposed schemes to other imaging modalities and
evaluate robustness against more types of realistic perturbations.

##  Appendix A Proof of Theorem 1

###  A-A Preliminary of Theorem 1

######  Lemma 1.

Let $f:\mathbb{R}^{d}\to\mathbb{R}^{m}$ be any bounded function. Let
$\bm{\eta}\sim\mathcal{N}(0,\sigma^{2}\mathbf{I})$. We define
$g:\mathbb{R}^{d}\to\mathbb{R}^{m}$ as

$g(\mathbf{x})=\mathbb{E}_{\bm{\eta}}[f(\mathbf{x}+\bm{\eta})].$

Then, $g$ is an $\frac{M}{\sqrt{2\pi}\sigma}$-Lipschitz map, where
$M=2\max_{\mathbf{x}\in\mathbb{R}^{d}}(\|f(\mathbf{x})\|_{2})$. In particular,
for any $\mathbf{x},\bm{\delta}\in\mathbb{R}^{d}$:

$\|g(\mathbf{x})-g(\mathbf{x}+\bm{\delta})\|_{2}\leq\frac{M}{\sqrt{2\pi}\sigma}\|\bm{\delta}\|_{2}.$

###### Proof.

The proof of this bound follows recent work [28], with a modification on $M$.
Let $\mu$ be the probability distribution function of random variable
$\bm{\eta}$. By the change of variables $\bm{w}=\mathbf{x}+\bm{\eta}$ and
$\bm{w}=\mathbf{x}+\bm{\eta}+\bm{\delta}$ for the integrals constituting
$g(\mathbf{x})$ and $g(\mathbf{x}+\bm{\delta})$, we have
$\|g(\mathbf{x})-g(\mathbf{x}+\bm{\delta})\|_{2}=\|\int_{\mathbb{R}^{d}}f(\bm{w})[\mu(\bm{w}-\mathbf{x})-\mu(\bm{w}-\mathbf{x}-\bm{\delta})]\
d\bm{w}\|_{2}$. Then, we have
$\|g(\mathbf{x})-g(\mathbf{x}+\bm{\delta})\|_{2}$

$\displaystyle\leq\int_{\mathbb{R}^{d}}\|f(\bm{w})[\mu(\bm{w}-\mathbf{x})-\mu(\bm{w}-\mathbf{x}-\bm{\delta})]\|_{2}\
d\bm{w},$

which is a standard result for the norm of an integral. We further apply
Holder’s inequality to upper bound
$\|g(\mathbf{x})-g(\mathbf{x}+\bm{\delta})\|_{2}$ with

$\displaystyle\max_{\mathbf{x}\in\mathbb{R}^{d}}(\|f(\mathbf{x})\|_{2})\int_{\mathbb{R}^{d}}|\mu(\bm{w}-\mathbf{x})-\mu(\bm{w}-\mathbf{x}-\bm{\delta})\
|d\bm{w}.$ (13)

Observe that $\mu(\bm{w}-\mathbf{x})\geq\mu(\bm{w}-\mathbf{x}-\bm{\delta})$ if
$\|\bm{w}-\mathbf{x}\|_{2}\leq\|\bm{w}-\mathbf{x}-\bm{\delta}\|_{2}$. Let
$D=\\{\bm{w}:\|\bm{w}-\mathbf{x}\|_{2}\leq\|\bm{w}-\mathbf{x}-\bm{\delta}\|_{2}\\}$.
Then, we can rewrite the above bound as

$\displaystyle=\max_{\mathbf{x}\in\mathbb{R}^{d}}(\|f(\mathbf{x})\|_{2})\cdot
2\int_{D}[\mu(\bm{w}-\mathbf{x})-\mu(\bm{w}-\mathbf{x}-\bm{\delta})]\ d\bm{w}$
(14) $\displaystyle=\frac{M}{2}\begin{pmatrix}2\int_{D}\mu(\bm{w}-\mathbf{x})\
d\bm{w}-2\int_{D}\mu(\bm{w}-\mathbf{x}-\bm{\delta})\ d\bm{w}.\end{pmatrix}$
(15)

Following Lemma 3 in [42], we obtain the bound

$\displaystyle 2\int_{D}\mu(\bm{w}-\mathbf{x})\
d\bm{w}-2\int_{D}\mu(\bm{w}-\mathbf{x}-\bm{\delta})\
d\bm{w}\leq\frac{2}{\sqrt{2\pi}\sigma}\|\bm{\delta}\|_{2}\>,$ (16)

which implies that
$\|g(\mathbf{x})-g(\mathbf{x}+\bm{\delta})\|_{2}\leq\frac{2\max_{\mathbf{x}\in\mathbb{R}^{d}}(\|f(\mathbf{x})\|_{2})}{\sqrt{2\pi}\sigma}\|\bm{\delta}\|_{2}=\frac{M}{\sqrt{2\pi}\sigma}\|\bm{\delta}\|_{2}$.
This completes the proof. ∎

###  A-B Proof of Theorem 1

###### Proof.

Assume that the data consistency step in MoDL at iteration $n$ is denoted by
$\mathbf{x}^{n}_{\text{M}}(\mathbf{A}^{H}\mathbf{y})$. We will sometimes drop
the input and $\mathbf{y}$ dependence for notational simplicity. Then

$\displaystyle\mathbf{x}^{1}_{\text{M}}=(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}(\mathbf{A}^{H}\mathbf{y}+\mathcal{D}_{\bm{\theta}}(\mathbf{A}^{H}\mathbf{y}))\>,$
(17)
$\displaystyle\mathbf{x}^{n}_{\text{M}}=(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}(\mathbf{A}^{H}\mathbf{y}+\mathcal{D}_{\bm{\theta}}(\mathbf{x}^{n-1}_{\text{M}})),$
(18)

where $\mathcal{D}_{\bm{\theta}}$ is the denoiser function. For the sake of
simplicity and consistency with the experiments, we use the weighting
parameter $\lambda=1$ (in the data consistency step). We note that the proof
works for arbitrary $\lambda$. SMUG introduces an iteration-wise smoothing
step into MoDL as follows:

$\displaystyle\mathbf{x}_{\text{S}}^{1}=((\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}(\mathbf{A}^{H}\mathbf{y}+\mathbb{E}_{\bm{\eta}_{1}}[\mathcal{D}_{\bm{\theta}}(\mathbf{A}^{H}\mathbf{y}+\bm{\eta}_{1})])$
(19)
$\displaystyle\mathbf{x}_{\text{S}}^{n}=((\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}(\mathbf{A}^{H}\mathbf{y}+\mathbb{E}_{\bm{\eta}_{n}}[\mathcal{D}_{\bm{\theta}}(\mathbf{x}_{\text{S}}^{n-1}+\bm{\eta}_{n})])$
(20)
$\displaystyle=(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}(\mathbf{A}^{H}\mathbf{y})+$
(21)
$\displaystyle(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\mathbb{E}_{\bm{\eta}_{n}}[\mathcal{D}_{\bm{\theta}}(\mathbf{x}_{\text{S}}^{n-1}+\bm{\eta}_{n})],$

where we apply the expectation to the denoiser $\mathcal{D}_{\bm{\theta}}$ at
each iteration. We use $\bm{\eta}_{n}$ to denote the noise during smoothing at
iteration $n$. The robustness error of SMUG after $n$ iterations is
$\|\mathbf{x}_{\text{S}}^{n}(\mathbf{A}^{H}\mathbf{y})-\mathbf{x}_{\text{S}}^{n}(\mathbf{A}^{H}(\mathbf{y}+\bm{\delta}))\|$.
We apply Lemma 1 and properties of the norm (e.g., triangle inequality) to
bound
$\|\mathbf{x}_{\text{S}}^{n}(\mathbf{A}^{H}\mathbf{y})-\mathbf{x}_{\text{S}}^{n}(\mathbf{A}^{H}(\mathbf{y}+\bm{\mathbf{\delta}}))\|$
as

$\displaystyle\leq\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\mathbf{A}^{H}\bm{\delta}\|$
(22)
$\displaystyle+\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\cdot\bigl{(}\mathbb{E}_{\bm{\eta}_{n}}[\mathcal{D}_{\bm{\theta}}\bigl{(}\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}\mathbf{y})+\bm{\eta}_{n}\bigr{)}]-$
$\displaystyle\mathbb{E}_{\bm{\eta}_{n}}[\mathcal{D}_{\bm{\theta}}\bigl{(}\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}(\mathbf{y}+\bm{\delta}))+\bm{\eta}_{n}\bigr{)}]\bigr{)}\|$
$\displaystyle\leq\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}\|\mathbf{A}^{H}\bm{\delta}\|_{2}$
$\displaystyle+\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}\|\mathbb{E}_{\bm{\eta}_{n}}[\mathcal{D}_{\bm{\theta}}\bigl{(}\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}\mathbf{y})+\bm{\eta}_{n}\bigr{)}]-$
$\displaystyle\mathbb{E}_{\bm{\eta}_{n}}[\mathcal{D}_{\bm{\theta}}\bigl{(}\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}(\mathbf{y}+\bm{\delta}))+\bm{\eta}_{n}\bigr{)}]\|$
$\displaystyle\leq\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}\|\mathbf{A}^{H}\bm{\delta}\|_{2}+\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}\times$
(23)
$\displaystyle\begin{pmatrix}\frac{M}{\sqrt{2\pi}\sigma}\end{pmatrix}\|\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}\mathbf{y})-\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}(\mathbf{y}+\bm{\mathbf{\delta}}))\|.$

Here, $M=2\max_{\mathbf{x}}(\|\mathcal{D}_{\bm{\theta}}(\mathbf{x})\|)$. Then
we plug in the expressions for
$\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}\mathbf{y})$ and
$\mathbf{x}_{\text{S}}^{n-1}(\mathbf{A}^{H}(\mathbf{y}+\bm{\mathbf{\delta}}))$
(from (20)) and bound their normed difference with
$\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\mathbf{A}^{H}\bm{\delta}\|+\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\cdot\bigl{(}\mathbb{E}_{\bm{\eta}_{n-1}}[\mathcal{D}_{\bm{\theta}}\bigl{(}\mathbf{x}_{\text{S}}^{n-2}(\mathbf{A}^{H}\mathbf{y})+\bm{\eta}_{n-1}\bigr{)}]-\mathbb{E}_{\bm{\eta}_{n-1}}[\mathcal{D}_{\bm{\theta}}\bigl{(}\mathbf{x}_{\text{S}}^{n-2}(\mathbf{A}^{H}(\mathbf{y}+\bm{\delta}))+\bm{\eta}_{n-1}\bigr{)}]\bigr{)}\|$.
This is bounded above similarly as for (22). We repeat this process until we
reach the initial $\mathbf{x}_{\text{S}}^{0}$ on the right hand side. This
yields the following bound involving a geometric series.

$\displaystyle\|\mathbf{x}_{\text{S}}^{n}(\mathbf{A}^{H}\mathbf{y})-\mathbf{x}_{\text{S}}^{n}(\mathbf{A}^{H}(\mathbf{y}+\bm{\mathbf{\delta}}))\|$
(24)
$\displaystyle\leq\|\mathbf{A}^{H}\bm{\delta}\|_{2}\begin{pmatrix}\sum_{j=1}^{n}\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}^{j}\cdot\begin{pmatrix}\frac{M}{\sqrt{2\pi}\sigma}\end{pmatrix}^{j-1}\end{pmatrix}$
$\displaystyle+\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}^{n}\begin{pmatrix}\frac{M}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}\|\mathbf{A}^{H}\bm{\delta}\|_{2}$
(25)
$\displaystyle\leq\|\mathbf{A}\|_{2}\|\bm{\delta}\|_{2}\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}\begin{pmatrix}\frac{1-\begin{pmatrix}\frac{M}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}^{n}}{1-\frac{M}{\sqrt{2\pi}\sigma}\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}}\end{pmatrix}$
$\displaystyle+\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}^{n}\begin{pmatrix}\frac{M}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}\|\mathbf{A}\|_{2}\|\bm{\delta}\|_{2}\leq
C_{n}\|\bm{\delta}\|_{2},$ (26)

where we used the geometric series formula, and
$C_{n}=\alpha\|\mathbf{A}\|_{2}\begin{pmatrix}\frac{1-\begin{pmatrix}\frac{M\alpha}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}}{1-\frac{M\alpha}{\sqrt{2\pi}\sigma}}\end{pmatrix}+\|\mathbf{A}\|_{2}\begin{pmatrix}\frac{M\alpha}{\sqrt{2\pi}\sigma}\end{pmatrix}^{n}$,
with $\alpha=\|(\mathbf{A}^{H}\mathbf{A}+\mathbf{I})^{-1}\|_{2}$. ∎

## References

  * [1] M. Lustig, D. L. Donohoand, et al.,  “Compressed Sensing MRI,”  IEEE signal processing magazine, 2008.
  * [2] M. Kivanc Mihcak, I. Kozintsev, K. Ramchandran, and P. Moulin,  “Low-complexity image denoising based on statistical modeling of wavelet coefficients,”  IEEE Signal Processing Letters, vol. 6, no. 12, pp. 300–303, 1999.
  * [3] S. Ma, W. Yin, Y. Zhang, and A. Chakraborty,  “An efficient algorithm for compressed MR imaging using total variation and wavelets,”  in 2008 IEEE Conference on Computer Vision and Pattern Recognition, 2008, pp. 1–8.
  * [4] S. Ravishankar and Y. Bresler,  “MR image reconstruction from highly undersampled k-space data by dictionary learning,”  IEEE Transactions on Medical Imaging, vol. 30, no. 5, pp. 1028–1041, 2011.
  * [5] S. G. Lingala and M. Jacob,  “Blind compressive sensing dynamic MRI,”  IEEE Transactions on Medical Imaging, vol. 32, no. 6, pp. 1132–1145, 2013.
  * [6] S. Ravishankar and Y. Bresler,  “Learning sparsifying transforms,”  IEEE Transactions on Signal Processing, vol. 61, no. 5, pp. 1072–1086, 2012.
  * [7] S. Ravishankar, J. C. Ye, and J. A.Fessler,  “Image reconstruction: From sparsity to data-adaptive methods and machine learning,”  Proceedings of the IEEE, vol. 108, no. 1, pp. 86–109, 2020.
  * [8] Ravishankar. B, Wen.and S, Pfister. L, and Bresler Y,  “Transform learning for magnetic resonance image reconstruction: From model-based learning to building neural networks,”  IEEE Signal Processing Magazine, vol. 37, no. 1, pp. 41–53, 2020.
  * [9] J. Schlemper, C. Qin, J. Duan, R. M. Summers, and K. Hammernik,  “Sigma-net: Ensembled iterative deep neural networks for accelerated parallel MR image reconstruction,”  arXiv preprint arXiv:1912.05480, 2019.
  * [10] S. Ravishankar, A. Lahiri, C. Blocker, and J. A. Fessler,  “Deep dictionary-transform learning for image reconstruction,”  in 2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018). IEEE, 2018, pp. 1208–1212.
  * [11] H. K. Aggarwal, M. P. Mani, and M. Jacob,  “Modl: Model-based deep learning architecture for inverse problems,”  IEEE transactions on medical imaging, vol. 38, no. 2, pp. 394–405, 2018.
  * [12] J. Schlemper, J. Caballero, J. V. Hajnal, A. N. Price, and D. Rueckert,  “A deep cascade of convolutional neural networks for dynamic MR image reconstruction,”  IEEE Trans. Med. Imaging, vol. 37, no. 2, pp. 491–503, Feb. 2018.
  * [13] O. Ronneberger, P. Fischer, and T. Brox,  “U-net: Convolutional networks for biomedical image segmentation,”  in Medical Image Computing and Computer-Assisted Intervention – MICCAI 2015, 2015, pp. 234–241.
  * [14] H. Zheng, F. Fang, and G. Zhang,  “Cascaded dilated dense network with two-step data consistency for MRI reconstruction,”  in NeurIPS, 2019.
  * [15] J. Schlemper, J. Caballero, J. V. Hajnal, A. Price, and D. Rueckert,  “A Deep Cascade of Convolutional Neural Networks for Dynamic MR Image Reconstruction,”  IEEE Transactions on Medical Imaging, vol. 37, no. 2, pp. 491–503, 2018.
  * [16] Y. Yang, J. Sun, H. Li, and Z. Xu,  “Deep ADMM-Net for compressive sensing MRI,”  in Advances in Neural Information Processing Systems, 2016, pp. 10–18.
  * [17] K. Hammernik, T. Klatzer, E. Kobler, M. P. Recht, D K. Sodickson, Thomas Pock, and Florian Knoll,  “Learning a variational network for reconstruction of accelerated MRI data,”  Magnetic resonance in medicine, vol. 79, no. 6, pp. 3055–3071, 2018.
  * [18] Y. Romano, M. Elad, and P. Milanfar,  “The Little Engine That Could: Regularization by Denoising (RED),”  SIAM Journal on Imaging Sciences, vol. 10, no. 4, pp. 1804–1844, 2017.
  * [19] G. T. Buzzard, S. H. Chan, S. Sreehari, and C. A. Bouman,  “Plug-and-play unplugged: optimization-free reconstruction using consensus equilibrium,”  SIAM J. Imaging Sci., vol. 11, no. 3, pp. 2001–20, Jan. 2018.
  * [20] V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen,  “On instabilities of deep learning in image reconstruction and the potential costs of AI,”  Proceedings of the National Academy of Sciences, 2020.
  * [21] C. Zhang, J. Jia, et al.,  “On Instabilities of Conventional Multi-Coil MRI Reconstruction to Small Adverserial Perturbations,”  arXiv preprint arXiv:2102.13066, 2021.
  * [22] D. Gilton, G. Ongie, and R. Willett,  “Deep equilibrium architectures for inverse problems in imaging,”  IEEE Transactions on Computational Imaging, vol. 7, pp. 1123–1133, 2021.
  * [23] A. Madry, A. Makelov, L. Schmidt, D. Tsipras, and A. Vladu,  “Towards deep learning models resistant to adversarial attacks,”  arXiv preprint arXiv:1706.06083, 2017.
  * [24] Y. Zhang, H.and Yu, J. Jiao, E. P. Xing, L. E. Ghaoui, and M. I. Jordan,  “Theoretically principled trade-off between robustness and accuracy,”  International Conference on Machine Learning, 2019.
  * [25] J. Cohen, E. Rosenfeld, and Z. Kolter,  “Certified adversarial robustness via randomized smoothing,”  in International Conference on Machine Learning. PMLR, 2019, pp. 1310–1320.
  * [26] H. Salman, M. Sun, G. Yang, A. Kapoor, and J. Z. Kolter,  “Denoised smoothing: A provable defense for pretrained classifiers,”  Advances in Neural Information Processing Systems, vol. 33, 2020.
  * [27] Y. Zhang, Y.and Yao, J. Jia, J. Yi, M. Hong, S. Chang, and S. Liu,  “How to robustify black-box ml models? a zeroth-order optimization perspective,”  arXiv preprint arXiv:2203.14195, 2022.
  * [28] A. Wolf,  “Making medical image reconstruction adversarially robust,”  2019\.
  * [29] H. Liu, J. Jia, S. Liang, Y. Yao, S. Ravishankar, and S. Liu,  “Smug: Towards robust mri reconstruction by smoothed unrolling,” 2023.
  * [30] B. Zoph, G. Ghiasi, T. Lin, Y. Cui, H. Liu, E. D. Cubuk, and Q. Le,  “Rethinking pre-training and self-training,”  Advances in Neural Information Processing Systems, vol. 33, 2020.
  * [31] D.L. Donoho,  “Compressed sensing,”  IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
  * [32] H. K. Aggarwal, M. P. Mani, and M. Jacob,  “MoDL: Model-based deep learning architecture for inverse problems,”  IEEE Trans. Med. Imaging, vol. 38, no. 2, pp. 394–405, Feb. 2019.
  * [33] Goodfellow I, Shlens. J, and Szegedy C,  “Explaining and harnessing adversarial examples,”  2015 ICLR, vol. arXiv preprint arXiv:1412.6572, 2015.
  * [34] J. Jia, M. Hong, Y.Zhang, M.Akcakaya, and S. Liu,  “On the robustness of deep learning-based MRI reconstruction to image transformations,”  in Workshop on Trustworthy and Socially Responsible Machine Learning, NeurIPS 2022, 2022.
  * [35] J. Zbontar, F. Knoll, A. Sriram, T. Murrell, Z. Huang, M. J. Muckley, A. Defazio, R. Stern, P. Johnson, M. Bruno, et al.,  “fastmri: An open dataset and benchmarks for accelerated mri,”  arXiv preprint arXiv:1811.08839, 2018.
  * [36] Y. Zhang., Y. Yao, J. Jia, J. Yi, M. Hong, S. Chang, and S. Liu,  “How to robustify black-box ML models? a zeroth-order optimization perspective,”  in International Conference on Learning Representations, 2022.
  * [37] D.P. Kingma and J.Ba,  “Adam: A method for stochastic optimization,”  2015 ICLR, vol. arXiv preprint arXiv:1412.6980, 2015.
  * [38] S. Yu, B. Park, and J. Jeong,  “Deep iterative down-up cnn for image denoising,”  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops, 2019, pp. 0–0.
  * [39] J. I. Tamir, F. Ong, J. Y. Cheng, M. Uecker, and M. Lustig,  “Generalized magnetic resonance image reconstruction using the berkeley advanced reconstruction toolbox,”  in ISMRM Workshop on Data Sampling & Image Reconstruction, Sedona, AZ, 2016.
  * [40] D. P Kingma and J. Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint arXiv:1412.6980, 2014.
  * [41] Z. Wang, A.C. Bovik, H.R. Sheikh, and P.E. Simoncelli,  “Image quality assessment: From error visibility to structural similarity,”  IEEE Trans. Image Process., vol. 13, no. 4, pp. 600–612, 2004.
  * [42] H. Lakshmanan, F. De, and P. Daniela,  “Decentralized resource allocation in dynamic networks of agents,”  SIAM Journal on Optimization, vol. 19, no. 2, pp. 911–940, 2008.


# Quantization of mode shifts in nanocavities integrated with atomically thin
sheets

N. Fang<EMAIL_ADDRESS>Nanoscale Quantum Photonics Laboratory, RIKEN
Cluster for Pioneering Research, Saitama 351-0198, Japan  D. Yamashita
<EMAIL_ADDRESS>Quantum Optoelectronics Research Team, RIKEN Center
for Advanced Photonics, Saitama 351-0198, Japan  S. Fujii  Quantum
Optoelectronics Research Team, RIKEN Center for Advanced Photonics, Saitama
351-0198, Japan  K. Otsuka  Nanoscale Quantum Photonics Laboratory, RIKEN
Cluster for Pioneering Research, Saitama 351-0198, Japan  Department of
Mechanical Engineering, The University of Tokyo, Tokyo 113-8656, Japan  T.
Taniguchi  International Center for Materials Nanoarchitectonics, National
Institute for Materials Science, Ibaraki 305-0044, Japan  K. Watanabe
Research Center for Functional Materials, National Institute for Materials
Science, Ibaraki 305-0044, Japan  K. Nagashio  Department of Materials
Engineering, The University of Tokyo, Tokyo 113-8656, Japan  Y. K. Kato
<EMAIL_ADDRESS>Nanoscale Quantum Photonics Laboratory, RIKEN Cluster
for Pioneering Research, Saitama 351-0198, Japan  Quantum Optoelectronics
Research Team, RIKEN Center for Advanced Photonics, Saitama 351-0198, Japan

###### Abstract

The unique optical properties of two-dimensional layered materials are
attractive for achieving increased functionality in integrated photonics.
Owing to the van der Waals nature, these materials are ideal for integrating
with nanoscale photonic structures. Here we report on carefully designed air-
mode silicon photonic crystal nanobeam cavities for efficient control through
two-dimensional materials. By systematically investigating various types and
thickness of two-dimensional materials, we are able to show that enhanced
responsivity allows for giant shifts of the resonant wavelength. With
atomically precise thickness over a macroscopic area, few-layer flakes give
rise to quantization of the mode shifts. We extract the dielectric constant of
the flakes and find that it is independent of the layer number down to a
monolayer. Flexible reconfiguration of a cavity is demonstrated by stacking
and removing ultrathin flakes. With an unconventional cavity design, our
results open up new possibilities for photonic devices integrated with two-
dimensional materials.

##  I Introduction

Two-dimensional (2D) layered materials such as graphene, transition metal
dichalcogenides (TMDs), and hexagonal boron nitride (h-BN) have attracted
considerable attention due to their exotic physical properties and potential
for diverse applications. Monolayer graphene exhibits the quantum Hall effect
even at room temperature [1], while ultra-thin TMD based field-effect
transistors are promising for continuation of the Moore’s Law [2]. Moiré
patterns can be created by layers with twisted angles, which has led to the
observations of superconductivity in bilayer graphene [3] and ferroelectricity
in bilayer h-BN [4]. Many TMDs are found to become direct bandgap
semiconductors at the monolayer limit, making them the thinnest optical
emitters [5]. Furthermore, rich optical phenomena exist in these atomically
thin flakes, such as spin-valley coupling [6], Moiré excitons [7], and layer
number dependent nonlinearity [8].

To further utilize the optical properties of TMDs, microcavities can be used
to enhance the light-matter interaction. By forming TMD/cavity hybrids,
nanolasers have been demonstrated [9, 10, 11] and second-order nonlinear
response has significantly improved [12, 13]. Evidence for strong coupling and
formation of microcavity polaritons have also been reported [14]. Compared to
the successful use of TMDs as an active material within microcavities, it is a
challenge to control the cavity modes with 2D materials. Because of the
atomically thin nature, cavity modes are barely affected and the resonant
wavelength can only be tuned slightly [15, 16, 17, 18, 19]. With atomically
precise thickness over macroscopic areas and extensive compatibility offered
by van der Waals interface, photonic devices controlled through single atomic
layers would present a new direction in nanotechnology. Cavity structures
designed specifically for enhanced sensitivity to ultra-thin materials is
required to unravel the full functionality and capability of the TMD/cavity
hybrids.

Here we substantially boost the responsiveness of microcavities to 2D
materials by using specially designed photonic crystal nanobeam structures.
Various 2D material types and thickness dependence are investigated, and we
find that it is possible to achieve giant wavelength shifts of the cavity
mode. By precisely controlling tungsten diselenide (WSe2) thickness down to a
monolayer, we observe quantized mode shifts of the 2D/cavity hybrids. The
dielectric constant in the atomically thin limit can be extracted by comparing
to finite-difference time-domain (FDTD) calculations. Stacking and removal of
WSe2 flakes are also performed, demonstrating the flexible reconfiguration
capability of the 2D/cavity hybrids.

##  II Results and discussion

### Design of 2D/cavity hybrids.

Figure 1:  Design of PC nanobeam cavities. (a) A schematic of a 2D/cavity
hybrid. (b) A top-view SEM image of a bare nanobeam cavity. Spatial
distribution of $E_{y}$ for the fundamental cavity modes of (c) dielectric-
mode and (d) air-mode structures with top views (top) and cross-section views
(bottom). PL spectra of the fundamental mode for (e) a dielectric-mode
nanobeam and (f) an air-mode nanobeam with (red) and without (black) a 9.0-nm-
thick WSe2 flake on top. The excitation power is 300 $\mu$W and the excitation
wavelength is 780 nm. Spectra are plotted with respect to the resonant
wavelengths before the transfer of WSe2 flakes, which are 1405.6 nm and 1477.4
nm for (e) and (f), respectively. We note that the initial resonant wavelength
in the present range barely affects the shift values in the FDTD simulation.
Dots are data and solid lines are Lorentzian peak fits. Arrows indicate the
shift directions. Scale bars in (b–d) are 1 $\mu$m.

We begin by designing cavities suited for interaction with 2D materials. Most
microcavities, such as microsphere resonators [20, 21], are designed based on
refractive index contrast between two media, and the modes are confined in the
high-index medium using total internal reflection. The evanescent fields are
typically small outside of the cavity, and coupling to 2D materials would be
weak. Reducing the size of the cavity would force the mode to extend outside,
and therefore smaller mode volumes are in general favorable for enhanced
evanescent fields. For example, toroidal microcavities having small diameters
show more of the modes occupying the space external to the cavity [22].

We propose to use photonic crystal (PC) nanobeam cavities with ultra-small
mode volumes [23], consisting of an array of air holes on a silicon waveguide
as shown in Fig. 1a, b. Bragg reflectors are formed by periodic arrays of
holes, and localized modes are formed by modulating the periodicity at the
cavity center to break the Bragg condition. The nanobeam cavities feature
ultra-small mode volumes on the order of $10^{-2}\lambda^{3}$ where $\lambda$
is the cavity wavelength. In addition, the air holes act to reduce the average
refractive index, weakening the confinement and increasing the evanescent
fields.

We further utilize the characteristics of the photonic bands to engineer the
cavity modes for enhancing the interaction with 2D materials. The periodic
modulation of the refractive index gives rise to a photonic bandgap,
separating the lower frequency dielectric band from the higher frequency air
band [24]. The dielectric bands are characterized by field amplitudes
maximized within the dielectric material, while the air bands exhibit large
fields in the air holes. The dielectric-band modes can be confined by locally
reducing the lattice constant, as the frequency of the modes will become
higher and the photons will be surrounded by the photonic bandgap [24].
Similarly, air-band modes can be confined by introducing a larger lattice
constant region [25, 26]. We expect that the air mode would have increased
coupling with 2D materials because of stronger fields in air. Furthermore, the
air band is closer to the light line than the dielectric band, corresponding
to a lower effective refractive index with stronger evanescent fields.

Figure 2:  Resonant wavelength shifts by different types of 2D materials. PL
spectra of an air-mode nanobeam before (black) and after (red) transfer of (a)
a 12.5-nm-thick h-BN flake, (b) a 11.0-nm-thick MoS2 flake, (c) a 10.8 nm-
thick WSe2 flake, and (d) a 11.9 nm-thick MoTe2 flake on top. The excitation
power is 300 $\mu$W and the excitation wavelength is 780 nm. Arrows indicate
the shift directions. Insets are optical images for the samples, and the scale
bars are 20 $\mu$m. Dots are data and solid lines are Lorentzian peak fits.

The mode profiles in the two different cavity designs are mapped out by
performing FDTD calculations. The $y$-component of the electric field $E_{y}$
for the fundamental transverse-electric mode in the dielectric-mode cavity is
shown in Fig. 1c, and higher field amplitudes within silicon are confirmed as
illustrated in the top view. In the air-mode cavity, the fields are mostly
distributed within the holes as shown in Fig. 1d. The mode profiles are
distinctly different not only in the $x$–$y$ plane but also for the evanescent
fields in the $z$ direction as can be seen in the cross sectional views. The
dielectric-mode cavity does show some evanescent fields, but the air-mode
cavity exhibits even larger evanescent fields that extend farther out into
free space as expected (see Supplementary Fig. 1).

We experimentally compare the dielectric mode and the air mode cavities by
transferring WSe2 flakes onto the devices. The cavities are designed to have
the fundamental mode in the telecommunication band based on the FDTD
calculations, and we fabricate the nanobeams from a silicon-on-insulator
substrate using electron beam lithography and inductively-coupled plasma
etching. WSe2 flakes with the same thickness of 9.0 nm are transferred on the
two types of cavities using a conventional polydimethylsiloxane (PDMS) stamp
method [27]. The interaction with the 2D material is evaluated with a home-
built confocal photoluminescence (PL) microscopy system at room temperature by
comparing the fundamental mode before and after the transfer of WSe2. All
measurements are performed in dry nitrogen to avoid environmental effects (see
Supplementary Fig. 2).

PL spectra for the dielectric-mode and the air-mode cavities are shown in Fig.
1e and f, respectively. A single sharp peak is observed in each spectrum and
is identified as the fundamental mode. The dielectric-mode cavity shows a
redshift of $13.5$ nm after the transfer, which is attributed to a change in
the average dielectric constant. In comparison, the air-mode cavity has an
increased shift value by almost twofold. The larger shift indicates the
enhanced responsivity, consistent with the simulations showing stronger fields
for the air mode cavities. We note that higher-order modes show even more
wavelength shifts because of the larger fraction of the mode volumes in free
space (see Supplementary Fig. 3).

We now survey the interaction of various 2D materials with the air-mode
cavities. Four typical 2D materials are studied including insulating h-BN and
semiconducting TMDs ranging from molybdenum disulfide (MoS2), WSe2, to
molybdenum ditelluride (MoTe2). We select flakes with thickness close to $11$
nm, and compare PL spectra before and after the transfer. As shown in Fig. 2,
various 2D materials show distinct shift values because of the different
dielectric constants. The smallest redshift of only $7.6$ nm is observed for
h-BN, which is a widely studied 2D material for photonics [28, 29]. For the
TMDs, the redshifts increase in the order of MoS2, WSe2, and MoTe2. The tail
of the MoTe2 emission spectrum shows some overlap with the cavity mode, which
is expected to lead to considerable cavity loss by absorption. In the
following study, we make use of WSe2 which exhibits large shifts with
negligible emission in the present wavelength range [30].

### Giant wavelength tuning in WSe2/cavity hybrids.

Figure 3:  Thickness dependence of resonant wavelength shifts. (a) Normalized
PL spectra of the fundamental mode for air-mode nanobeams before (red) and
after (black) transfer of WSe2 flakes. The data from top to bottom correspond
to thickness of 40.3, 35.7, 22.0, 15.5 nm. Spectra are plotted with respect to
the resonant wavelengths before the transfer, which are 1378.1, 1420.0,
1460.4, 1439.7 nm, respectively. The excitation power is 300 $\mu$W and the
excitation wavelength is 780 nm. Arrows indicate the shift directions. Dots
are data and solid lines are Lorentzian peak fits. (b) Resonant wavelength
shifts as a function of WSe2 thickness. Dots are data and lines are FDTD
simulations with $\epsilon$ of 18.0 (purple), 19.0 (green), and 20.0 (orange).
We note that several samples with the thin WSe2 flakes deviate from the
simulations, which could be due to the formation of an air space between the
2D material and the nanobeam (see Supplementary Fig. 4).

To explore the capability of the enhanced responsiveness in WSe2/cavity
hybrids, modifications of the cavity mode by thick WSe2 flakes are
investigated. We transfer flakes with different thickness $t$ onto the
cavities and mode shifts are measured by PL spectroscopy. As shown in Fig. 3a,
the redshifts become larger with increasing $t$ owing to a larger mode
overlap. The largest value of $201.8$ nm is achieved for $t=40.3$ nm, where
the mode wavelength is close to the detection limit of our detector. We note
that the tuning range covers a wide spectrum corresponding to the
telecommunication E-band to L-band.

We are able to extract the dielectric constant $\epsilon$ from the thickness
dependent shifts by comparing to FDTD simulations (Fig. 3b). The calculated
shifts are almost linearly dependent on the thickness when $t<10.0$ nm, and
become superlinear in the thicker region which may be due to changes in the
mode profile. The experimental data is reproduced well when $\epsilon=19\pm
1$, whose value corresonds to the dielectric constant of bulk WSe2. We note
that in-plane $\epsilon$ is important for determining the redshifts because
the fundamental mode is transverse electric (see Supplementary Fig. 5).

We also study the quality factor $Q$ in our WSe2/cavity hybrids since large
wavelength shifts are often accompanied by increased optical losses. As the
spectral resolution of the present PL spectroscopy system is insufficient,
transmittance measurements are performed to determine the $Q$ factor [31, 32].
For $t=22.0$ nm sample, the resonant wavelength is shifted by over $110$ nm
while the $Q$ factor remains almost unaltered at $2.2\times 10^{4}$ (see
Supplementary Fig. 6). The negligible effect on the $Q$ factor could be
attributed to a homogeneous interface with suppressed scattering loss as well
as the absence of absorption from the WSe2 flake at the cavity wavelength.

### Quantized wavelength shifts by atomically thin WSe2.

Figure 4:  Quantized wavelength shifts induced by atomically thin WSe2 flakes.
(a) A top-view SEM image of a nanobeam cavity with a monolayer WSe2 flake
transferred above. Scale bar is 2 $\mu$m. (b) Normalized PL spectra of the
fundamental modes for air-mode nanobeams with a monolayer (blue), bilayer
(green), trilayer (orange), and quadlayer WSe2 flake (red) on top. The
horizontal axis is the shift with respect to the resonant wavelength before
the transfer, which are 1461.6, 1454.9, 1464.2, 1468.4 nm, respectively. The
excitation power is 300 $\mu$W and the excitation wavelength is 780 nm. Dots
are data and solid lines are Lorentzian plus linear fits. (c) Resonant
wavelength shifts for all the measured samples. (d) Resonant wavelength shifts
as a function of layer number. Dots are data and the broken line is a linear
fit. (e) Extracted dielectric constant as a function of the layer number.
Thickness of each layer is taken to be 0.65 nm in the analysis.

With the enhanced interaction, even a minor perturbation by 2D materials would
influence the cavity mode and become detectable. By controlling the thickness
of WSe2 flakes with atomic precision, we demonstrate extreme sensitivity down
to the monolayer limit. Since ultra-thin 2D materials are too fragile for
sustaining the strain induced by PMDS stamps, deformation including rupture
often occurs in the suspended region [19, 33]. Here we utilize the anthracene-
assisted transfer method [34], where the organic molecule single crystal
supports and protects 2D flakes during the process (see Supplementary Fig.
7-9). Figure 4a shows a scanning electron microscopy (SEM) image of the
prepared monolayer WSe2/cavity hybrid, where a uniform morphology over the
nanobeam can be observed.

The typical PL spectra of nanobeam cavities that are integrated with quad-
(4L), tri- (3L), bi- (2L), and mono-layer (1L) WSe2 flakes are shown in Fig.
4b. The wavelength shift decreases by $\sim$3 nm as the thickness is reduced
layer by layer. It is notable that the 1L WSe2, the limit of the scaling,
still tunes the mode of the nanobeam cavity effectively. The observed shifts
are almost two orders of magnitude larger than the linewidth, assuming $Q\sim
10^{4}$.

The roughly constant spacings between the cavity resonances suggest that the
mode shifts are quantized due to the discrete thickness of atomically thin
layers. To verify that the mode shifts are indeed governed by the layer
number, we have prepared three sets of samples for each thickness and repeated
the measurements (Fig. 4c). Clearly resolved steps are observed, confirming
the quantization of the mode shifts. It is remarkable that the reproducible
resonant shifts are obtained, benefitting from atomically precise thickness of
the 2D material over the entire cavity.

The observed discrete wavelength shifts offer an opportunity to study the
dielectric constant even in the ultra-thin limit. We compare experimental
results with FDTD simulation to extract $\epsilon$ for 1L to 4L WSe2. Since
calculated wavelength shifts are linearly dependent on $t$ in the thin limit
(Fig. 3b), we fit the experimental results with a simple linear model (Fig.
4d). We observe a small intercept of 0.7 nm, which may come from anthracene
residues as well as adsorption of other molecules on the surface of the 2D
material. We have also performed FDTD simulations for various $\epsilon$, and
find that wavelength shift $\Delta\lambda$ is also linearly dependent on
$\epsilon$ (see Supplementary Fig. 5). The data can therefore be described by

$\Delta\lambda=\alpha(\epsilon-1)t+\lambda_{0},$ (1)

where $\alpha=0.21$ is the coefficient determined by the FDTD calculations and
$\lambda_{0}$ is the intercept. Considering that anthracene residue or surface
adsorption should occur similarly in all the samples, we take
$\lambda_{0}=0.7$ nm. The values of $\epsilon$ are extracted from data using
eq. (1), and the results are shown in Fig. 4e. The obtained values are
comparable with the bulk value, revealing that the dielectric screening in the
atomically thin limit is as effective.

Since the cavity resonances arise from transverse-electric modes, the
extracted values correspond to the in-plane $\epsilon$. It is reported that
in-plane $\epsilon$ measured by spectroscopic ellipsometry decreases with
increasing the layer number for ultra-thin TMDs [35, 36]. In comparison,
first-principle calculations show that in-plane $\epsilon$ is insensitive to
the layer number since the main contribution comes from in-plane atomic
displacements [37]. Our results show no apparent thickness dependence down to
the monolayer limit, which is consistent with the calculations. The behavior
of in-plane $\epsilon$ is important for 2D materials based electronic devices,
for example, the tunneling length through Schottky barriers in metal/2D
contacts is strongly affected by in-plane $\epsilon$ [38]. We note that
extrinsic effects from processing and adsorption have been taken into account
through $\lambda_{0}$ in our analysis, and the extracted values should
therefore represent the intrisic dielectric properties of atomically thin
WSe2.

### Flexible reconfiguration of 2D/cavity hybrids.

Figure 5:  Reconfiguration of a cavity by stacking and removing of WSe2
flakes. PL spectra of the fundamental mode of the as-fabricated cavity (grey),
after the first transfer of a 2L WSe2 flake (blue), after the second transfer
of a 3L flake (purple), and after the removal of the stacked flakes (green).
The excitation power is 300 $\mu$W and the excitation wavelength is 780 nm.
Dots are data and solid lines are Lorentzian plus linear fits. Insets
illustrate the configuration of the device during the measurements.

One of the advantages of using 2D materials comes from their versatile
manipulation capabilities, where flakes can be placed in the desired order
without the restriction of lattice matching. In addition, the removal of
transferred flakes is possible owing to the weak van der Waals force. In Fig.
5, we demonstrate additive and reversible control with a WSe2/cavity hybrid. A
2L WSe2 flake is first transferred on top of a nanobeam, and we observe a mode
shift by $6.5$ nm (blue curve). Next, another transfer of a 3L flake is
performed, which further shifts the mode by $9.5$ nm (purple curve). The
shifts are consistent with the quantized values shown in Fig. 4c. We note that
the linewidth increases slightly, suggesting additional loss from scattering
induced by 2D-2D interfacial bubbles. By picking up the stack using an
anthracene stamp (see Supplementary Fig. 8), the $Q$ factor is recovered and
the mode is shifted back.

The high compatibility between 2D materials and cavities demonstrated here
opens up various applications in photonics. The integration with 2D materials
could be used as a post-fabrication tuning technology, where the giant
wavelength tunable range of over $200$ nm can cover the telecommunication
E-band through L-band. By using ultra-thin flakes, fine tuning can be done
with high predictability and reproducibility because of the quantized shifts
arising from the atomically precise thickness. This is in contrast to
continuous mode shifts obtained in conventional methods such as chemical
etching [39], coating of photochromic films [40], and condensation of gases
[41]. In addition, the stackable and removable 2D materials offer a new degree
of freedom for reconfigurable cavities. Finally, as demonstrated in the
measurement of $\epsilon$, 2D/cavity hybrids are an ideal platform for
investigating the optical properties of atomically thin structures due to the
high sensitivity of the air-mode cavities. Sensing of a small number of
molecules could become possible by monitoring a high-$Q$ cavity resonance,
which would extend the limits of the 2D materials based sensing technologies
[42].

In conclusion, we have developed an air-mode cavity for integration with 2D
materials. By investigating various 2D material types and thickness
dependence, we show that the fabricated 2D/cavity hybrids have a giant
wavelength tunable range of over $200$ nm. We observe a clear quantization
behavior of the wavelength shift for WSe2 where each step represents an effect
from an additional monolayer. The in-plane dielectric constant of WSe2 can be
extracted because of the high surface sensitivity of the air-mode cavity, and
shows a thickness independent characteristic down to the monolayer. By
stacking and removing the transferred flakes, we demonstrate that the hybrids
benefit from the versatile manipulation capabilities of the 2D materials.
These findings reveal fascinating features of the cavity that is efficiently
controlled by 2D materials, and provide a universal design strategy for
enhancing the light-matter interaction with nanomaterials.

## Methods

### Nanobeam cavity fabrication.

The nanobeam cavities are fabricated on a silicon-on-insulator substrate with
a 260-nm-thick top silicon layer and a 1-$\mu$m-thick buried-oxide (BOX)
layer. After defining the PC pattern on a resist mask by electron beam
lithography, the pattern is transferred to the top silicon slab through an
inductively coupled plasma process using C4F8 and SF6 gases. Following resist
removal, the BOX layer is etched with 20% hydrofluoric acid to form an air-
suspended nanobeam structure.

### PL measurements.

A homebuilt confocal microscopy system is used to perform PL measurements at
room temperature in dry nitrogen gas. We use a wavelength-tunable Ti:sapphire
laser for excitation with its power controlled by neutral density filters.
Polarization angle is adjusted to match the cavity mode by a half-wave plate.
The laser beam is focused on the samples using an objective lens with a
numerical aperture of 0.65 and a working distance of 4.5 mm. The $1/e^{2}$
beam diameter and the collection spot size defined by a confocal pinhole are
$\sim$1.2 and 5.4 $\mu$m, respectively. PL is collected through the same
objective lens and detected using a liquid-nitrogen-cooled 1024 pixel InGaAs
diode array attached to a spectrometer. PL is dispersed by a 150 lines/mm
grating and the dispersion is 0.52 nm/pixel at the center wavelength of 1340
nm.

### Transfer of 2D materials by PDMS stamps.

We transfer thick flakes of 2D materials with $t\textgreater 5$ nm on the
cavities by a conventional PDMS stamp method to form 2D/cavity hybrids [27].
2D material flakes are prepared on a PDMS sheet (Gelfilm by Gelpak) by
mechanical exfoliation and then transferred on the target nanobeam cavity at
120∘C by using a micromanipulator system. The thickness of the transferred 2D
material flakes is determined by atomic force microscopy. h-BN crystals are
supplied by NIMS and other 2D material crystals are purchased from HQ
graphene.

### Transfer of 2D materials by anthracene crystals.

Anthracene crystals are grown with the same method as in Ref. [34]. 2D flakes
are prepared on standard 90-nm-thick SiO2/Si substrates by mechanical
exfoliation and the layer number is determined by optical contrast. An
anthracene single crystal is picked up with a glass-supported PDMS sheet to
form an anthracene/PDMS stamp. 2D flakes are picked up by pressing the
anthracene/PDMS stamp against a substrate with the target 2D flakes, followed
by quick separation ($\textgreater 10$ mm/s) so that the anthracene crystal
remains attached to the PDMS sheet. The stamp is then pressed on a receiving
substrate with the cavity. By slowly peeling off the PDMS ($\textless 0.2$
$\mu$m/s), the anthracene crystal with the target 2D flake is released.
Sublimation of anthracene in air at 110∘C for 10 min leaves behind clean 2D
flakes because contamination from solvents is absent in the all-dry process.

###### Acknowledgements.

Parts of this study are supported by JSPS (KAKENHI JP20H02558, JP20K15199,
JP20J00817, JP20K15137, JP20K15120, JP19K23593, JP18H03864, JP19H00755,
JP21H05237, JP21H05232), MIC (SCOPE 191503001), and MEXT (Nanotechnology
Platform JPMXP09F19UT0075). The growth of hBN crystals is supported by the
Elemental Strategy Initiative conducted by the MEXT, Japan (Grant Number
JPMXP0112101001) and JSPS KAKENHI (Grant Numbers JP19H05790, JP20H00354 and
JP21H05233). D.Y. are supported by JSPS (Research Fellowship for Young
Scientists). N.F. and S.F. are supported by RIKEN Special Postdoctoral
Researcher Program. We thank the Advanced Manufacturing Support Team at RIKEN
for technical assistance. We also thank H. Kumazaki for helpful discussion of
the dimpled fiber fabrication.

## Author Contributions

N.F. performed transfer of 2D materials as well as measurements of 2D
materials and 2D/cavity hybrids. D.Y. performed simulation and fabrication of
nanobeam cavities as well as measurements of 2D/cavity hybrids. S.F. developed
the transmittance measurement system for 2D/cavity hybrids. K.O. helped
develop anthracene-assisted dry method. T.T. and K.W. provided h-BN crystals.
K.N. helped measurements of 2D materials. Y.K.K. supervised the project. N.F.,
D.Y., and Y.K.K. co-wrote the manuscript. All the authors commented on the
manuscript. N.F. and D.Y. contributed equally to this work.

## References

  * Novoselov  _et al._ [2007] K. S. Novoselov, Z. Jiang, Y. Zhang, S. V. Morozov, H. L. Stormer, U. Zeitler, J. C. Maan, G. S. Boebinger, P. Kim, and A. K. Geim, Room-temperature quantum hall effect in graphene, Science 315, 1379 (2007).
  * Desai  _et al._ [2016] S. B. Desai, S. R. Madhvapathy, A. B. Sachid, J. P. Llinas, Q. Wang, G. H. Ahn, G. Pitner, M. J. Kim, J. Bokor, C. Hu, H. S. Wong, and A. Javey, MoS2 transistors with 1-nanometer gate lengths, Science 354, 99 (2016).
  * Cao  _et al._ [2018] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018).
  * Yasuda  _et al._ [2021] K. Yasuda, X. Wang, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Stacking-engineered ferroelectricity in bilayer boron nitride, Science 372, 1458 (2021).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically thin MoS2: a new direct-gap semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * Xiao  _et al._ [2012] D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Coupled spin and valley physics in monolayers of MoS2 and other group-VI dichalcogenides, Phys. Rev. Lett. 108, 196802 (2012).
  * Jin  _et al._ [2019] C. Jin, E. C. Regan, A. Yan, M. I. B. Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng, S. Shi, K. Watanabe, T. Taniguchi, S. Tongay, A. Zettl, and F. Wang, Observation of Moiré excitons in WSe2/WS2 heterostructure superlattices, Nature 567, 76 (2019).
  * Säynätjoki  _et al._ [2017] A. Säynätjoki, L. Karvonen, H. Rostami, A. Autere, S. Mehravar, A. Lombardo, R. A. Norwood, T. Hasan, N. Peyghambarian, H. Lipsanen, _et al._ , Ultra-strong nonlinear optical processes and trigonal warping in MoS2 layers, Nat. Commun. 8, 893 (2017).
  * Li  _et al._ [2017] Y. Li, J. Zhang, D. Huang, H. Sun, F. Fan, J. Feng, Z. Wang, and C. Z. Ning, Room-temperature continuous-wave lasing from monolayer molybdenum ditelluride integrated with a silicon nanobeam cavity, Nat. Nanotechnol. 12, 987 (2017).
  * Wu  _et al._ [2015] S. Wu, S. Buckley, J. R. Schaibley, L. Feng, J. Yan, D. G. Mandrus, F. Hatami, W. Yao, J. Vučković, A. Majumdar, and X. Xu, Monolayer semiconductor nanocavity lasers with ultralow thresholds, Nature 520, 69 (2015).
  * Shang  _et al._ [2017] J. Shang, C. Cong, Z. Wang, N. Peimyoo, L. Wu, C. Zou, Y. Chen, X. Y. Chin, J. Wang, C. Soci, W. Huang, and T. Yu, Room-temperature 2D semiconductor activated vertical-cavity surface-emitting lasers, Nat. Commun. 8, 543 (2017).
  * Busschaert  _et al._ [2020] S. Busschaert, R. Reimann, M. Cavigelli, R. Khelifa, A. Jain, and L. Novotny, Transition metal dichalcogenide resonators for second harmonic signal enhancement, ACS Photonics 7, 2482 (2020).
  * Fryett  _et al._ [2016] T. K. Fryett, K. L. Seyler, J. Zheng, C.-H. Liu, X. Xu, and A. Majumdar, Silicon photonic crystal cavity enhanced second-harmonic generation from monolayer WSe2, 2D Mater. 4, 015031 (2016).
  * Liu  _et al._ [2015] X. Liu, T. Galfsky, Z. Sun, F. Xia, E.-c. Lin, Y.-H. Lee, S. Kéna-Cohen, and V. M. Menon, Strong light–matter coupling in two-dimensional atomic crystals, Nat. Photonics 9, 30 (2015).
  * Zhao  _et al._ [2018] C. Zhao, X. Gan, Q. Yuan, S. Hu, L. Fang, and J. Zhao, High-performance volatile organic compounds microsensor based on few-Layer MoS2-coated photonic crystal cavity, Adv. Opt. Mater. 6, 1 (2018).
  * Majumdar  _et al._ [2013] A. Majumdar, J. Kim, J. Vuckovic, and F. Wang, Electrical control of silicon photonic crystal cavity by graphene, Nano Lett. 13, 515 (2013).
  * He  _et al._ [2021] J. He, I. Paradisanos, T. Liu, A. R. Cadore, J. Liu, M. Churaev, R. N. Wang, A. S. Raja, C. Javerzac-Galy, P. Roelli, D. D. Fazio, B. L. Rosa, S. Tongay, G. Soavi, A. C. Ferrari, and T. J. Kippenberg, Low-loss integrated nanophotonic circuits with layered semiconductor materials, Nano Lett. 21, 2709 (2021).
  * Datta  _et al._ [2020] I. Datta, S. H. Chae, G. R. Bhatt, M. A. Tadayon, B. Li, Y. Yu, C. Park, J. Park, L. Cao, D. Basov, _et al._ , Low-loss composite photonic platform based on 2D semiconductor monolayers, Nat. Photonics 14, 256 (2020).
  * Fryett  _et al._ [2018] T. K. Fryett, Y. Chen, J. Whitehead, Z. M. Peycke, X. Xu, and A. Majumdar, Encapsulated silicon nitride nanobeam cavity for hybrid nanophotonics, ACS Photonics 5, 2176 (2018).
  * Righini  _et al._ [2011] G. Righini, Y. Dumeige, P. Féron, M. Ferrari, G. N. Conti, D. Ristic, and S. Soria, Whispering gallery mode microresonators: fundamentals and applications, La Riv. del Nuovo Cim. 34, 435 (2011).
  * Hall  _et al._ [2015] J. M. Hall, V. S. Afshar, M. R. Henderson, A. François, T. Reynolds, N. Riesen, and T. M. Monro, Method for predicting whispering gallery mode spectra of spherical microresonators, Opt. Express 23, 9924 (2015).
  * Min  _et al._ [2007] B. Min, L. Yang, and K. Vahala, Perturbative analytic theory of an ultrahigh-Q toroidal microcavity, Phys. Rev. A 76, 013823 (2007).
  * Foresi  _et al._ [1997] J. Foresi, P. R. Villeneuve, J. Ferrera, E. Thoen, G. Steinmeyer, S. Fan, J. Joannopoulos, L. Kimerling, H. I. Smith, and E. Ippen, Photonic-bandgap microcavities in optical waveguides, Nature 390, 143 (1997).
  * Notomi  _et al._ [2008] M. Notomi, E. Kuramochi, and H. Taniyama, Ultrahigh-${Q}$ nanocavity with 1D photonic gap, Opt. Express 16, 11095 (2008).
  * Zhang  _et al._ [2011] Y. Zhang, I. Bulu, W.-M. Tam, B. Levitt, J. Shah, T. Botto, and M. Loncar, High-$Q/V$ air-mode photonic crystal cavities at microwave frequencies, Opt. Express 19, 9371 (2011).
  * Miura  _et al._ [2014] R. Miura, S. Imamura, R. Ohta, A. Ishii, X. Liu, T. Shimada, S. Iwamoto, Y. Arakawa, and Y. K. Kato, Ultralow mode-volume photonic crystal nanobeam cavities for high-efficiency coupling to individual carbon nanotube emitters, Nat. Commun. 5, 5580 (2014).
  * Castellanos-Gomez  _et al._ [2014] A. Castellanos-Gomez, M. Buscema, R. Molenaar, V. Singh, L. Janssen, H. S. Van Der Zant, and G. A. Steele, Deterministic transfer of two-dimensional materials by all-dry viscoelastic stamping, 2D Mater. 1, 011002 (2014).
  * Wong  _et al._ [2015] D. Wong, J. Velasco, L. Ju, J. Lee, S. Kahn, H.-Z. Tsai, C. Germany, T. Taniguchi, K. Watanabe, A. Zettl, _et al._ , Characterization and manipulation of individual defects in insulating hexagonal boron nitride using scanning tunnelling microscopy, Nat. Nanotechnol. 10, 949 (2015).
  * Hayee  _et al._ [2020] F. Hayee, L. Yu, J. L. Zhang, C. J. Ciccarino, M. Nguyen, A. F. Marshall, I. Aharonovich, J. Vučković, P. Narang, T. F. Heinz, and J. A. Dionne, Revealing multiple classes of stable quantum emitters in hexagonal boron nitride with correlated optical and electron microscopy, Nat. Mater. 19, 534 (2020).
  * Kozawa  _et al._ [2014] D. Kozawa, R. Kumar, A. Carvalho, K. K. Amara, W. Zhao, S. Wang, M. Toh, R. M. Ribeiro, A. C. Neto, K. Matsuda, _et al._ , Photocarrier relaxation pathway in two-dimensional semiconducting transition metal dichalcogenides, Nat. Commun. 5, 4543 (2014).
  * McCutcheon  _et al._ [2011] M. W. McCutcheon, P. B. Deotare, Y. Zhang, and M. Lončar, High-${Q}$ transverse-electric/transverse-magnetic photonic crystal nanobeam cavities, Appl. Phys. Lett. 98, 111117 (2011).
  * Tetsumoto  _et al._ [2015] T. Tetsumoto, Y. Ooka, and T. Takasumi, High-Q coupled resonances on a PhC waveguide using a tapered nanofiber with high coupling efficiency, Opt. Express 23, 16256 (2015).
  * Onodera  _et al._ [2020] M. Onodera, S. Masubuchi, R. Moriya, and T. Machida, Assembly of van der Waals heterostructures: exfoliation, searching, and stacking of 2D materials, Jpn. J. Appl. Phys. 59, 010101 (2020).
  * Otsuka  _et al._ [2021] K. Otsuka, N. Fang, D. Yamashita, T. Taniguchi, K. Watanabe, and Y. K. Kato, Deterministic transfer of optical-quality carbon nanotubes for atomically defined technology, Nat. Commun. 12, 3138 (2021).
  * Yu  _et al._ [2015] Y. Yu, Y. Yu, Y. Cai, W. Li, A. Gurarslan, H. Peelaers, D. E. Aspnes, C. G. Van De Walle, N. V. Nguyen, Y. W. Zhang, and L. Cao, Exciton-dominated dielectric function of atomically thin MoS2 films, Sci. Rep. 5, 16996 (2015).
  * Song  _et al._ [2019] B. Song, H. Gu, M. Fang, X. Chen, H. Jiang, R. Wang, T. Zhai, Y. T. Ho, and S. Liu, Layer-dependent dielectric function of wafer-scale 2D MoS2, Adv. Opt. Mater. 7, 1801250 (2019).
  * Laturia  _et al._ [2018] A. Laturia, M. L. V. de Put, and W. G. Vandenberghe, Dielectric properties of hexagonal boron nitride and transition metal dichalcogenides: From monolayer to bulk, npj 2D Mater. Appl. 2, 6 (2018).
  * Schulman  _et al._ [2018] D. S. Schulman, A. J. Arnold, and S. Das, Contact engineering for 2D materials and devices, Chem. Soc. Rev. 47, 3037 (2018).
  * Kuwabara  _et al._ [2019] M. Kuwabara, S. Noda, and Y. Takahashi, Ultrahigh-Q photonic nanocavity devices on a dual thickness SOI substrate operating at both 1.31- and 1.55-µm telecommunication wavelength bands, Laser Photonics Rev. 13, 1800258 (2019).
  * Sridharan  _et al._ [2010] D. Sridharan, E. Waks, G. Solomon, and J. T. Fourkas, Reversible tuning of photonic crystal cavities using photochromic thin films, Appl. Phys. Lett. 96 (2010).
  * Mosor  _et al._ [2005] S. Mosor, J. Hendrickson, B. C. Richards, J. Sweet, G. Khitrova, H. M. Gibbs, T. Yoshie, A. Scherer, O. B. Shchekin, and D. G. Deppe, Scanning a photonic crystal slab nanocavity by condensation of xenon, Appl. Phys. Lett. 87, 141105 (2005).
  * Mao  _et al._ [2017] S. Mao, J. Chang, H. Pu, G. Lu, Q. He, H. Zhang, and J. Chen, Two-dimensional nanomaterial-based field-effect transistors for chemical and biological sensing, Chem. Soc. Rev. 46, 6872 (2017).