#  Representation Sparsification with Hybrid Thresholding
for Fast SPLADE-based Document Retrieval

Yifan Qiao  Department of Computer Science, University of CaliforniaSanta
BarbaraCalifornia93106USA ,  Yingrui Yang  Department of Computer Science,
University of CaliforniaSanta BarbaraCalifornia93106USA ,  Shanxiu He
Department of Computer Science, University of CaliforniaSanta
BarbaraCalifornia93106USA and  Tao Yang  Department of Computer Science,
University of CaliforniaSanta BarbaraCalifornia93106USA

(2023)

###### Abstract.

Learned sparse document representations using a transformer-based neural model
has been found to be attractive in both relevance effectiveness and time
efficiency. This paper describes a representation sparsification scheme based
on hard and soft thresholding with an inverted index approximation for faster
SPLADE-based document retrieval. It provides analytical and experimental
results on the impact of this learnable hybrid thresholding scheme.

Learned sparse representations, top-k retrieval, index pruning.

††ccs: Information systems Retrieval efficiency††journalyear: 2023††copyright:
rightsretained††conference: Proceedings of the 46th International ACM SIGIR
Conference on Research and Development in Information Retrieval; July 23–27,
2023; Taipei, Taiwan††booktitle: Proceedings of the 46th International ACM
SIGIR Conference on Research and Development in Information Retrieval (SIGIR
’23), July 23–27, 2023, Taipei, Taiwan††doi: 10.1145/3539618.3592051††isbn:
978-1-4503-9408-6/23/07

##  1\. Introduction

Recently learned sparse retrieval techniques (Zamani et al., 2018; Dai and
Callan, 2020; Mallia et al., 2021; Lin and Ma, 2021; Gao et al., 2021; Formal
et al., 2021b, a, 2022) have become attractive because such a representation
can deliver a strong relevance by leveraging transformer-based models to
expand document tokens with learned weights and can take advantage of
traditional inverted index based retrieval techniques (Mallia et al., 2017,
2022). Its query processing is cheaper than a dense representation which
requires GPU support (e.g. (Xiong et al., 2021; Ren et al., 2021; Santhanam et
al., 2021)) even with efficiency optimization through approximate nearest
neighbor search (Johnson et al., 2019; Zhan et al., 2022; Xiao et al., 2022).

This paper focuses on the SPLADE family of sparse representations (Formal et
al., 2021b, a, 2022) because it can deliver a high MRR@10 score for MS MARCO
passage ranking (Craswell et al., 2020) and a strong zero-shot performance for
the BEIR datasets (Thakur et al., 2021), which are well-recognized IR
benchmarks. The sparsification optimization in SPLADE has used L1 and FLOPS
regularization to minimize non-zero weights during model learning, and our
objective is to exploit additional opportunities to further increase the
sparsity of inverted indices produced by SPLADE. Earlier static inverted index
pruning research (Carmel et al., 2001; Blanco and Barreiro, 2007; Büttcher and
Clarke, 2006) for a lexical model has shown the usefulness of trimming a term
posting list or a document by a limit. Yang et al. (Yang et al., 2021) conduct
top token masking by limiting the top activated weight count uniformly per
document and gradually reducing this weight count limit to a targeted constant
during training. Motivated by these studies (Carmel et al., 2001; Blanco and
Barreiro, 2007; Büttcher and Clarke, 2006; Yang et al., 2021) and since they
have not addressed the learnability of a pruning limit through relevance-
driven training, this paper exploits a learnable thresholding architecture to
filter out unimportant neural weights produced by the SPLADE model through
joint training.

The contribution of this paper is a learnable hybrid hard and soft
thresholding scheme with an inverted index approximation to increase the
sparsity of SPLADE-based document and query feature vectors for faster
retrieval. In addition to experimental validation with MS MARCO and BEIR
datasets, we provide an analysis of the impact of hybrid thresholding with
joint training on index approximation errors and training update
effectiveness.

##  2\. Background

For a query $q$ and a document $d$, after expansion and encoding, they can be
represented by vector $\vec{w}(q)$ and $\vec{w}(d)$ with length $|V|$, where
$V$ is the vocabulary set. The rank score of $q$ and $d$ is computed as
$R(q,d)=\vec{w}(q)\cdot\vec{w}(d)=\sum_{i=1}^{|V|}w_{i}^{q}\times w_{i}^{d}$.
For sparse vectors with many zeros, retrieval can utilize a data structure
called inverted index during online inference for fast score computation
(Mallia et al., 2017, 2022). The SPLADE model uses the BERT token space to
predict the feature vector $\vec{w}$. In its latest SPLADE++ model, it first
calculates the importance of $i$-th input token in $d$ for each $j$ in $V$:
$w_{ij}(\Theta)={\rm Transform}(\vec{h_{i}})^{T}\vec{E_{j}}+b_{j},$ where
$\vec{h_{i}}$ is the BERT embedding of $i$-th token in $d$, $\vec{E_{j}}$ is
the BERT input embedding for $j$-th token. Transform() is a linear layer with
GeLU activation and LayerNorm. The weights in this linear layer,
$\vec{E_{j}}$, and $b_{j}$ are the SPLADE parameters updated during training
and we call them set $\Theta$. Then the $j$-th entry $w_{j}$ of document $d$
(or a query) is max-pooled as $w_{j}(\Theta)=\max_{i\in d}\\{\log(1+{\rm
ReLU}(w_{ij}(\Theta)))\\}.$ Notice that $w_{j}\geq 0$.

The loss function of SPLADE models (Formal et al., 2021b, a, 2022) contains a
per-query ranking loss $L_{R}$ and sparsity regularization. The ranking loss
has evolved from a log likelihood based function for maximizing positive
document probability to margin MSE for knowledge distillation. This paper uses
the loss of SPLADE with a combination that delivers the best result in our
training process. $L_{R}$ is the ranking loss with margin MSE for knowledge
distillation (Hofstätter et al., 2020). The document token regularization
$L_{D}$ is computed on the training documents in each batch based on FLOPS
regularization. The query token regularization $L_{Q}$ is based on L1 norm.
Let $B$ be a set of training queries with $N$ documents involved in a batch.
$L_{Q}=\sum_{j\in V}\frac{1}{|B|}\sum_{q\in B}w_{j}^{q};$ $L_{D}=\sum_{j\in
V}(\frac{1}{N}\sum_{d=1}^{N}w_{j}^{d})^{2}.$

Related work. Other than SPLADE, sparse retrieval studies include SNRM (Zamani
et al., 2018), DeepCT (Dai and Callan, 2020), DeepImpact (Mallia et al.,
2021), and uniCOIL (Lin and Ma, 2021; Gao et al., 2021). The sparsity of a
neural network is studied in the deep learning community. Soft thresholding in
(Kusupati et al., 2020) adopts a learnable threshold with function
$S(x,t)=ReLU(x-t)$ to make parameter $x$ zero under threshold $t$. A hard
thresholding function $H(x,t)=x$ when $x\geq t$ otherwise 0. Approximate hard
thresholding (Park and Lee, 2020) uses a Gauss error function to approximate
$H(x,t)$ with smooth gradients. Dynamic sparse training (LIU et al., 2020)
finds a dynamic threshold with marked layers. These works including the recent
ones (Frantar and Alistarh, 2022) are targeted for sparsification of parameter
edges in a deep neural network. In our context, a token weight $w_{j}$ is an
output node in a network. The sparsification of output nodes is addressed in
activation map compression (Georgiadis, 2019) using ReLU as soft thresholding
together with L1 regularization. The work of (Kurtz et al., 2020) further
boosts sparsity with the Hoyer regularization and a variant of ReLU. The above
techniques have not been investigated in the context of sparse retrieval, and
the impact of thresholding on relevance and query processing time with
inverted indices, requires new design considerations and model structuring for
document retrieval, even the previous work can be leveraged.

##  3\. Hybrid Thresholding (HT)

Figure 1. Hybrid thresholding with an index approximation

Design considerations. To zero out a token weight below a learnable threshold,
there are two options: soft thresholding (Kusupati et al., 2020), and
approximate hard thresholding (Park and Lee, 2020). For query token weights,
we find that soft thresholding does not affect relevance significantly. For
document token weights, our study finds that compared to soft thresholding,
hard thresholding can retain relevance better since it does not change token
weights when exceeding a threshold. Since the subgradient for hard
thresholding with respect to a threshold is always 0, an approximation needs
to be carried out for training. For search index generation, an inverted index
produced with the same approximate hard thresholding as training keeps many
unnecessary non-zero document token weights, slowing down retrieval
significantly. Thus we directly apply hard thresholding with a threshold
learned from training, as shown in Figure 1. There is a gap between trained
document token weights and actual weights used in our inverted index
generation and online inference, and we intend to minimize this gap (called an
index approximation error).

Thus our design takes a hybrid approach that applies soft thresholding to
query token weights during training and inference and applies approximate hard
thresholding to document token weights during training while using hard
thresholding for documents during index generation. For approximate hard
thresholding, we propose to use a logistic sigmoid-based function instead of a
Gauss error function (Park and Lee, 2020). This sigmoid thresholding
simplifies our analysis of the impact of its hyperparameter choice to index
approximation errors, and to training stability.

###  3.1. Trainable and approximate thresholding

Training computes threshold parameters $t_{D}$, and $t_{Q}$ for documents and
queries, respectively. From the output of the SPLADE model, every token weight
of a query is replaced with $S(w_{j}^{q},t_{Q})$, which is
$ReLU(w_{j}^{q}-t_{Q})$, and every document token weight is replaced with
$\widehat{H}(w_{j}^{q},t_{D})$ before their dot product is computed during
training as shown in Figure 1(a). Sigmoid thresholding $\widehat{H}$ is
defined as:

(1)
${\rm\widehat{H}}(w_{j}^{d},t_{D})=w_{j}^{d}\sigma(K(w_{j}^{d}-t_{D}))\mbox{
where }\sigma(x)=\frac{1}{1+e^{-x}}.$

Here $K$ is a hyperparameter to control the slope steepness of step
approximation that jumps from 0 to 1 when exceeding a threshold.

The indexing process uses hard thresholding to replace all document weights
that are below threshold $t_{D}$ as 0 as depicted in Figure 1(b). The above
post processing introduces an index approximation error
$E=|\widehat{H}(w_{j}^{d},t_{D})-H(w_{j}^{d},t_{D})|$. We derive its upper
bound as follows. Notice that $w_{j}\geq 0$, and for any $x\geq 0$, $1+x\leq
e^{x}$.

$E=w_{j}^{d}\sigma(K(w_{j}^{d}-t_{D}))=\frac{w_{j}^{d}}{1+e^{K(t_{D}-w_{j}^{d})}}\leq\frac{w_{j}^{d}}{2+K(t_{D}-w_{j}^{d})}.$

When $w_{j}^{d}\geq t_{D}$, we can derive that

$E=w_{j}^{d}(1-\sigma(K(w_{j}^{d}-t_{D})))=w_{j}^{d}\sigma(K(t_{D}-w_{j}^{d})))\leq\frac{w_{j}^{d}}{2+K(w_{j}^{d}-t_{D})}.$

Let $\sigma^{-}$ denote $\sigma(K(w_{i}^{d}-t_{D}))$. $0<\sigma^{-}<1$. In
both of the above cases, the error upper bound is minimized when $K$ is large.
This is consistent with the fact that error $E$ is monotonically decreasing as
$K$ increases because $\frac{\partial E}{\partial K}=-w_{j}^{d}$
$\sigma^{-}(1-\sigma^{-})$ $|w_{j}^{d}-t_{D}|\leq 0$. When $|w_{j}^{d}-t_{D}|$
is big, the error is negligible and when $|w_{j}^{d}-t_{D}|$ is small, the
error could become big with a small $K$ value. But as shown later, an
excessively large $K$ value could cause a big parameter update during a
training step, affecting joint training stability.

Let $Dlen$ and $Qlen$ be the non-zero token weight count of document $d$ and
query $q$, respectively. For our hybrid thresholding,
$Dlen=\sum_{j}\textbf{1}_{w_{j}^{d}\geq
t_{D}},Qlen=\sum_{j}\textbf{1}_{w_{j}^{q}\geq t_{Q}}.$ Here $\textbf{1}_{x\geq
y}$ is an indicator function as 1 if $x\geq y$ otherwise 0. When increasing
$t_{D}$ and $t_{Q}$, $Dlen$ and $Qlen$ decrease. Thus for a batch of training
queries $B$, the original SPLADE loss is extended as:
$L=(\frac{1}{|B|}\sum_{q\in
B}L_{R})+\lambda_{Q}L_{Q}+\lambda_{D}L_{D}+\lambda_{T}L_{T}.$ The extra item
added is $L_{T}=\log(1+e^{-t_{D}})+\log(1+e^{-t_{Q}})$. We retain the original
$L_{Q}$ and $L_{D}$ expressions because as $w_{j}^{q}$ or $w_{j}^{d}$
decreases, more weights can quickly be zeroed out.

###  3.2. Threshold and token weight updating

We study the change of $t_{D}$, $t_{Q}$, $w_{j}^{d}$, and $w_{j}^{q}$ after
each training step with a mini-batch gradient descent update. The analysis
below uses the first-order Taylor polynomial approximation and follows the
fact that sigmoid thresholding $\widehat{H}$ and soft thresholding function
$S$ are used independently for a query and a document in the loss function.
Symbol $\alpha$ is the learning rate. Let “$d\triangleleft q$” mean $d$ is a
positive or negative document of query $q$.

$\Delta t_{D}=t^{new}_{D}-t^{old}_{D}=-\alpha\frac{\partial L}{\partial
t_{D}}=-\alpha\left(\frac{1}{|B|}\sum_{q\in B}\frac{\partial
L_{R}}{\partial\widehat{H}}\frac{\partial\widehat{H}}{\partial
t_{D}}+\lambda_{T}\frac{\partial L_{T}}{\partial t_{D}}\right)\\\
=\alpha\left(\frac{1}{|B|}\sum_{q\in B}\left(K\frac{\partial
L_{R}}{\partial\widehat{H}}\sum_{d\triangleleft
q}\sum_{i}w_{i}^{d}(1-\sigma)^{-}\sigma^{-}\right)+\lambda_{T}\frac{e^{-t_{D}}}{1+e^{-t_{D}}}\right).$

$\Delta t_{Q}=t_{Q}^{new}-t_{Q}^{old}=-\alpha\frac{\partial L}{\partial
t_{Q}}=-\alpha\left(\frac{1}{|B|}\sum_{q\in B}\frac{\partial L_{R}}{\partial
S}\frac{\partial S}{\partial t_{Q}}+\lambda_{T}\frac{\partial L_{T}}{\partial
t_{Q}}\right)=\alpha\left(\frac{1}{|B|}\sum_{q\in B}\left(\frac{\partial
L_{R}}{\partial S}\sum_{i}\textbf{1}_{w_{i}^{q}\geq
t_{Q}}\right)+\lambda_{T}\frac{e^{-t_{Q}}}{1+e^{-t_{Q}}}\right).$

$\Delta
w_{j}^{d}=w_{j}^{d,new}-w_{j}^{d,old}\approx\sum_{\theta\in\Theta}\frac{\partial
w_{j}^{d}}{\partial\theta}\Delta\theta=-\sum_{\theta\in\Theta}\frac{\partial
w_{j}^{d}}{\partial\theta}\alpha\frac{\partial
L}{\partial\theta}=-\alpha\sum_{\theta\in\Theta}\frac{\partial
w_{j}^{d}}{\partial\theta}\left(\frac{1}{|B|}\sum_{q\in B}\left(\frac{\partial
L_{R}}{\partial\widehat{H}}\left(\sum_{d\triangleleft
q}\sum_{i}\frac{\partial\widehat{H}}{\partial w_{i}^{d}}\frac{\partial
w_{i}^{d}}{\partial\theta}\right)+\\\ \frac{\partial L_{R}}{\partial
S}\left(\sum_{i}\frac{\partial S}{\partial w_{i}^{q}}\frac{\partial
w_{i}^{q}}{\partial\theta}\right)\right)+\lambda_{D}\frac{\partial
L_{D}}{\partial\theta}+\lambda_{Q}\frac{\partial
L_{Q}}{\partial\theta}\right).$

Notice that $\frac{\partial\widehat{H}}{\partial
w_{i}^{d}}=\sigma^{-}+Kw_{i}^{d}\sigma^{-}(1-\sigma^{-})$. The above results
indicate:

  * •

A significant number of terms in $\Delta t_{D}$ and $\Delta w_{j}^{d}$ involve
linear coefficient $K$. This is verifiably true also for $\Delta w_{j}^{q}$.
Although a large $K$ value can minimize the index approximation error
$|\widehat{H}(w_{j}^{d},t_{D})-H(w_{j}^{d},t_{D})|$, it can cause an
aggressive change of token weights and thresholds at a training iteration,
making training overshoot and miss the global optimum. Thus $K$ cannot be too
large, and our evaluation further studies this.

  * •

If $\frac{\partial L_{R}}{\partial\widehat{H}}\geq 0$, $\Delta t_{D}\geq 0$,
and the document threshold increases, decreasing $Dlen$. Otherwise document
token threshold may decrease after a parameter update step during training,
and the degree of decreasing is reduced by a positive value
$\frac{e^{-t_{D}}}{1+e^{-t_{D}}}$. Based on the sign of $\frac{\partial
L_{R}}{\partial S}$, we can draw a similar conclusion on $\Delta t_{Q}$.

##  4\. Evaluation

Our evaluation uses MS MARCO passages (Craswell et al., 2020) and BEIR
datasets (Thakur et al., 2021). MS MARCO has 8.8M passages while BEIR has 13
different datasets of varying sizes up-to 5.4M. As a common practice, we
report the relevance in terms of mean reciprocal rank MRR@10 for the MS MARCO
passage Dev query set with 6980 queries, and the normalized discounted
cumulative gain nDCG@10 (Järvelin and Kekäläinen, 2002) for its DL’19 and
DL’20 sets, and also for BEIR. For retrieval with a SPLADE inverted index, we
report the mean response time (MRT) and 99th percentile time ($P_{99}$) in
milliseconds. The query encoding time is not included. For the SPLADE model,
we warm up it following (Formal et al., 2022; Lassance and Clinchant, 2022),
and train it with $\lambda_{Q}=0.01$ and $\lambda_{D}=0.008$, and hybrid
thresholding. We use the PISA (Mallia et al., 2019b) search system to index
documents and search queries using SIMD-BP128 compression (Lemire and Boytsov,
2015) and MaxScore retrieval (Mallia et al., 2019a, 2022). Our evaluation runs
as a single thread on a Linux CPU-only server with Intel i5-8259U 2.3GHz and
32GB memory. Similar retrieval latency results are observed on a 2.3GHz AMD
EPYC 7742 processor. The checkpoints and related code will be released in
https://github.com/Qiaoyf96/HT.

Table 1. Overall results on MS MARCO passages

Methods | MRR | MRT($P_{99}$) | MRT($P_{99}$) | nDCG | nDCG | Dlen
---|---|---|---|---|---|---
| Dev | $k=10$ | $k=1000$ | DL’19 | DL’20 |
SPLADE | 0.3966 | 48.3(228) | 127(408) | 0.7398 | 0.7340 | 351
/DT (Park and Lee, 2020) | 0.3922 | 102(457) | 262(786) | 0.7392 | 0.7319 | 444
/Top305 (Yang et al., 2021) | 0.3962 | 42.4(202) | 114(369) | 0.7353 | 0.7288 | 277
/Top100 (Yang et al., 2021) | 0.3908 | 21.8(106) | 62.5(196) | 0.7192 | 0.7119 | 99
/DCP50% (Büttcher and Clarke, 2006) | 0.3958 | 30.0(145) | 83.9(271) | 0.7385 | 0.7321 | 175
/DCP40% (Büttcher and Clarke, 2006) | 0.3933 | 25.9(124) | 73.3(235) | 0.7335 | 0.7280 | 140
/DCP30% (Büttcher and Clarke, 2006) | 0.3912 | 21.6(101) | 61.8(193) | 0.7287 | 0.7217 | 105
/Cut0.5 | 0.3924 | 21.9(104) | 62.6(195) | 0.7296 | 0.7212 | 144
/Cut0.8 | 0.3885 | 15.6(70.4) | 43.8(128) | 0.7207 | 0.7118 | 112
/HT1 | 0.3955 | 22.8(108) | 62.3(195) | 0.7322 | 0.7210 | 140
/HT3 | 0.3942 | 14.2(67.2) | 40.6(123) | 0.7327 | 0.7228 | 106
/HT1-2GTI (Qiao et al., 2023) | 0.3959 | 10.0(49.1) | 27.6(92.2) | 0.7330 | 0.7210 | 140
/HT3-2GTI (Qiao et al., 2023) | 0.3942 |  6.9(33.9) | 19.3(62.1) | 0.7320 | 0.7228 | 106

Overall results with MS MARCO. Table 1 is a comparison with the baselines on
MS MARCO passage Dev set, DL’19, and DL’20. It lists the average $Dlen$ value,
and top-$k$ retrieval time with depth $k=10$ and 1000. Row 3 is for original
SPLADE trained by ourselves with an MRR number higher than 0.38 reported in
(Formal et al., 2022; Lassance and Clinchant, 2022). Rows 12 and 13 list the
result of our hybrid thresholding marked as HT${}_{\lambda_{T}}$ and $K=25$.
With $\lambda_{T}=1$, SPLADE/HT1 converges to a point where $t_{Q}=0.4$ and
$t_{D}=0.5$, which is about 2x faster in retrieval. HT3 with $\lambda_{T}=3$
converges at $t_{Q}=0.7$ and $t_{D}=0.8$, resulting 3.1x speedup than SPLADE
while having a slightly lower MRR@10 0.3942. No statistically significant
degradation in relevance has been observed at the 95% confidence level for
both HT1 and HT3. The inverted index size reduces from 6.4GB for original
SPLADE to 2.8GB and 2.2GB for HT1 and HT3 respectively. When applying two-
level guided traversal 2GTI (Qiao et al., 2023) with its fast configuration,
Rows 14 and 15 show a further latency reduction to 6.9ms or 19.3ms.

We discuss other baselines listed in this table. Row 4 named DT uses the
thresholding scheme from (Park and Lee, 2020). Its training does not converge
with its loss function, and its retrieval is much slower. Rows 5 and 6 follow
joint training of top-$k$ masking (Yang et al., 2021) with the top 305 tokens
as suggested in (Yang et al., 2021) and with the top 100 tokens. Rows 7, 8 and
9 marked with DCP$x$ follow document centric pruning (Büttcher and Clarke,
2006) that keeps $x$ of top tokens per document where $x$=50%, 40%, and 30%.
We did not list term centric pruning (Carmel et al., 2001; Blanco and
Barreiro, 2007) because (Büttcher and Clarke, 2006) shows DCP is slightly
better in relevance under the same latency constraint. Rows 10 and 11 with
“/Cut0.5” and “/Cut0.8” apply a hard threshold with 0.5 and 0.8 in the output
of original SPLADE without joint training. The index pruning options without
learning from Rows 5 to 11 can either reduce the latency to the same level as
HT, but their relevance score is visibly lower; or have a relevance similar to
HT but with much slower latency. This illustrates the advantage of learned
hybrid thresholding with joint training.

.

Table 2. Zero-shot performance on BEIR datasets

| SPLADE |  SPLADE/HT1 |  SPLADE/HT3
---|---|---|---
Dataset | nDCG | MRT | nDCG | MRT | nDCG | MRT
DBPedia | 0.430 | 135 | 0.435 | 64.2 | 0.426 | 32.3
FiQA | 0.354 | 6.5 | 0.345 | 4.0 | 0.336 | 3.2
NQ | 0.547 | 81.8 | 0.545 | 45.9 | 0.539 | 28.6
HotpotQA | 0.678 | 481 | 0.680 | 265 | 0.678 | 140
NFCorpus | 0.351 | 0.5 | 0.352 | 0.3 | 0.346 | 0.2
T-COVID | 0.719 | 16.0 | 0.730 | 10.1 | 0.695 | 7.5
Touche-2020 | 0.307 | 15.0 | 0.306 | 9.3 | 0.313 | 4.5
ArguAna | 0.440 | 20.8 | 0.463 | 7.8 | 0.500 | 4.1
C-FEVER | 0.234 | 1375 | 0.219 | 681 | 0.213 | 332
FEVER | 0.781 | 1584 | 0.778 | 559 | 0.764 | 264
Quora | 0.806 | 17.5 | 0.776 | 9.2 | 0.792 | 4.5
SCIDOCS | 0.151 | 6.9 | 0.155 | 3.0 | 0.151 | 2.0
SciFact | 0.676 | 5.7 | 0.681 | 2.4 | 0.672 | 1.4
Average | 0.498 | - | 0.497 | 2.0x | 0.494 | 3.6x

Table 2 lists the zero-shot performance of HT when $k=1000$ by applying the
SPLADE/HT model learned from MS MARCO to the BEIR datasets without any
additional training. HT1 has a similar nDCG@10 score as SPLADE without HT,
while having a 2x MRT speedup on average. HT3 is even faster with 3.6x
speedup, and its nDCG@10 drops in some degree to 0.494.

(a)

(b)

Figure 2.  Weight/threshold/sparsity changes during training

Figure 2 depicts the average values of $w_{j}^{d}$, $t_{D}$, and $Dlen$ on the
left and $w_{j}^{q}$, $t_{Q}$, and $Qlen$ on the right during MS MARCO
training under $HT_{1}$. $x$-axis is the training epoch number. It shows that
$Dlen$ and $Qlen$ decrease while $t_{D}$ and $t_{Q}$ increase as training
makes progress and SPLADE/HT1 converges after about 20 epochs.

Design options. Table 3 lists performance under 4 thresholding combinations
from Row 3 to Row 7. $S[x]$ means soft thresholding function $S()$ is applied
to $x$ for both training and indexing where $x$ can be documents (D) or
queries (Q). $\widehat{H}[$x$]$ means sigmoid thresholding $\widehat{H}$ is
applied in both training and indexing. $\widehat{H}H[x]$ means $\widehat{H}$
is applied in training and $H$ is applied in indexing. $\phi[x]$ means no
thresholding is applied to $x$ during training and indexing. When thresholding
is not applied to queries, $\widehat{H}H[D]$ is 1.3x faster than $S[D]$ when
$k=10$ and $k=1000$ while their relevance scores are similar. Shifting of
document weight distribution by soft thresholding significantly affects
retrieval time. Rows 6 and 7 fix $\widehat{H}H[D]$ setting, and show that soft
thresholding is more effective in relevance than hard thresholding for query
tokens. Shifting of query weight distribution has less effect on latency while
gaining more relevance through model consistency between training and
indexing.

Table 3. Impact of design options. MS MARCO passages.

HT Config. | MRR | MRT($P_{99}$) | MRT($P_{99}$) | Qlen | Dlen
---|---|---|---|---|---
$\lambda_{T}=1$ |  | $k=10$ | $k=1000$ |  |
Soft vs. hard thresholding in 4 combinations. Fix $K=25$.
$\phi[Q],S[D]$ | 0.3941 | 31.7(157) | 91.5(315) | 14.3 | 145
$\phi[Q],\widehat{H}H[D]$ | 0.3942 | 24.1(111) | 70.7(219) | 13.5 | 142
$S[Q],\widehat{H}H[D]$ | 0.3955 | 22.8(108) |  62.3(195) | 11.3 | 140
$\widehat{H}H[Q],\widehat{H}H[D]$ | 0.3904 | 24.9(106) |  62.6(182) | 9.0 | 142
Vary $K$. $\widehat{H}[D]$ vs. $\widehat{H}H[D]$. Fix $S[Q]$.
$\widehat{H}H[D]$, $K=2.5$ | 0.3947 | 22.8(110) | 62.6(199) | 11.5 | 149
$\widehat{H}[D]$, $K=2.5$ | 0.3963 | 41.4(198) | 112(358) | 11.5 | 421
$\widehat{H}H[D]$, $K=25$ | 0.3955 | 22.8(108) | 62.3(195) | 11.3 | 140
$\widehat{H}[D]$, $K=25$ | 0.3961 | 28.7(136) | 76.9(239) | 11.3 | 208
$\widehat{H}H[D]$, $K=250$ | 0.3946 | 21.9(102) | 60.5(189) | 11.2 | 135
$\widehat{H}[D]$, $K=250$ | 0.3947 | 23.1(112) | 63.9(203) | 11.2 | 159
Usefulness of $L_{Q}$ and $L_{D}$. Fix $S[Q],\widehat{H}H[D]$, and $K=25$.
w/o $L_{Q}$ | 0.3956 | 56.2(245) | 166(502) | 20.1 | 138
w/o $L_{Q}$, $L_{D}$ | 0.3954 | 99.4(434) | 254(772) | 25.9 | 421

Hyperparameter $K$ in sigmoid thresholding $\widehat{H}$. Table 3 compares
$\widehat{H}H[D]$ with $\widehat{H}[D]$ when varying $K$ from Row 8 to Row 14.
In these cases, training always uses $\widehat{H}$ while indexing uses
$\widehat{H}$ or $H$. When $K$ is small as 2.5, applying $\widehat{H}$ to both
training and indexing yields good relevance, but retrieval is about 1.8x
slower because much more non-zero weights are kept in the index. When $K$
becomes large as 250, training does not converge to the global optimum due to
large update sizes, resulting in an MRR score lower than $K$=25 even with no
index approximation. $K=25$ has a reasonable MRR while $\widehat{H}H[D]$ is
up-to 26% faster than $\widehat{H}[D]$.

Retaining $L_{Q}$ and $L_{D}$. Last three rows of Table 3 shows that the query
length is higher when $L_{Q}$ is removed from the loss function, and documents
get longer when $L_{D}$ is removed further. The result means $L_{Q}$ and
$L_{D}$ are useful in sparsity control together with $L_{T}$.

##  5\. Concluding Remarks

Our evaluation shows that learnable hybrid thresholding with index
approximation can effectively increase the sparsity of inverted indices with
2-3x faster retrieval and competitive or slightly degraded relevance (0.28% -
0.6% MRR@10 drop). Its trainability allows relevance and sparsity guided
threshold learning and it can outperform index pruning without such an
optimization. Our scheme retains a non-uniform number of non-zero token
weights per vector based on a trainable weight and threshold difference for
flexibility in relevance optimization. Our analysis shows that hyperparameter
$K$ in sigmoid thresholding needs to be chosen judiciously for a small index
approximation error without hurting training stability.

If a small relevance tradeoff is allowed, more retrieval time reduction is
possible when applying other related orthogonal efficiency optimization
techniques (Lassance and Clinchant, 2022; Mallia et al., 2022; Qiao et al.,
2022, 2023; Lin and Trotman, 2015; Mackenzie et al., 2021). Applying hybrid
thresholding HT3 to a checkpoint of a recent efficiency-driven SPLADE model
(Lassance and Clinchant, 2022) with 0.3799 MRR@10 on the MS MARCO passage Dev
set, decreases the response time from 36.6ms to 21.7ms (1.7x faster) when
$k$=1000 while having 0.3868 MRR@10. This latency can be further reduced to
14.2ms with the same MRR@10 number (0.3868) when 2GTI (Qiao et al., 2023) is
applied to the above index.

A future study is to investigate the use of the proposed hybrid thresholding
scheme for other learned sparse models (Mallia et al., 2021; Lin and Ma, 2021;
Gao et al., 2021).

Acknowledgments. We thank Wentai Xie and anonymous referees for their valuable
comments and/or help. This work is supported in part by NSF IIS-2225942 and
has used computing resource of NSF’s ACCESS program. Any opinions, findings,
conclusions or recommendations expressed in this material are those of the
authors and do not necessarily reflect the views of the NSF.

## References

  * (1)
  * Blanco and Barreiro (2007) Roi Blanco and Alvaro Barreiro. 2007.  Boosting Static Pruning of Inverted Files. In _Proc. of SIGIR_ (Amsterdam, The Netherlands) _(SIGIR ’07)_. Association for Computing Machinery, New York, NY, USA, 777–778.  https://doi.org/10.1145/1277741.1277904
  * Büttcher and Clarke (2006) Stefan Büttcher and Charles L. A. Clarke. 2006.  A Document-Centric Approach to Static Index Pruning in Text Retrieval Systems. In _Proceedings of the 15th ACM International Conference on Information and Knowledge Management_ (Arlington, Virginia, USA) _(CIKM ’06)_. Association for Computing Machinery, New York, NY, USA, 182–189.  https://doi.org/10.1145/1183614.1183644
  * Carmel et al. (2001) David Carmel, Doron Cohen, Ronald Fagin, Eitan Farchi, Michael Herscovici, Yoelle S. Maarek, and Aya Soffer. 2001.  Static Index Pruning for Information Retrieval Systems. In _Proc. of SIGIR_ (New Orleans, Louisiana, USA) _(SIGIR ’01)_. Association for Computing Machinery, New York, NY, USA, 43–50.  https://doi.org/10.1145/383952.383958
  * Craswell et al. (2020) Nick Craswell, Bhaskar Mitra, Emine Yilmaz, Daniel Fernando Campos, and Ellen M. Voorhees. 2020\.  Overview of the TREC 2020 Deep Learning Track.  _ArXiv_ abs/2102.07662 (2020).
  * Dai and Callan (2020) Zhuyun Dai and Jamie Callan. 2020.  Context-Aware Term Weighting For First Stage Passage Retrieval.  _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_ (2020).
  * Formal et al. (2021a) Thibault Formal, C. Lassance, Benjamin Piwowarski, and Stéphane Clinchant. 2021a.  SPLADE v2: Sparse Lexical and Expansion Model for Information Retrieval.  _ArXiv_ abs/2109.10086 (2021).
  * Formal et al. (2022) Thibault Formal, Carlos Lassance, Benjamin Piwowarski, and Stéphane Clinchant. 2022.  From Distillation to Hard Negative Sampling: Making Sparse Neural IR Models More Effective.  _Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval_ (2022).
  * Formal et al. (2021b) Thibault Formal, Benjamin Piwowarski, and Stéphane Clinchant. 2021b.  SPLADE: Sparse Lexical and Expansion Model for First Stage Ranking.  _Proceedings of the 44th International ACM SIGIR Conference on Research and Development in Information Retrieval_ (2021).
  * Frantar and Alistarh (2022) Elias Frantar and Dan Alistarh. 2022.  SPDY: Accurate Pruning with Speedup Guarantees. In _Proceedings of the 39th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 162)_ , Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato (Eds.). PMLR, 6726–6743.  https://proceedings.mlr.press/v162/frantar22a.html
  * Gao et al. (2021) Luyu Gao, Zhuyun Dai, and Jamie Callan. 2021.  COIL: Revisit Exact Lexical Match in Information Retrieval with Contextualized Inverted List.  _NAACL_ (2021).
  * Georgiadis (2019) Georgios Georgiadis. 2019\.  Accelerating convolutional neural networks via activation map compression. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 7085–7095.
  * Hofstätter et al. (2020) Sebastian Hofstätter, Sophia Althammer, Michael Schröder, Mete Sertkan, and Allan Hanbury. 2020\.  Improving Efficient Neural Ranking Models with Cross-Architecture Knowledge Distillation.  _ArXiv_ abs/2010.02666 (2020).
  * Järvelin and Kekäläinen (2002) Kalervo Järvelin and Jaana Kekäläinen. 2002.  Cumulated gain-based evaluation of IR techniques.  _ACM Transactions on Information Systems (TOIS)_ 20, 4 (2002), 422–446.
  * Johnson et al. (2019) Jeff Johnson, Matthijs Douze, and Hervé Jégou. 2019\.  Billion-scale similarity search with GPUs.  _IEEE Trans. on Big Data_ 7, 3 (2019), 535–547.
  * Kurtz et al. (2020) Mark Kurtz, Justin Kopinsky, Rati Gelashvili, Alexander Matveev, John Carr, Michael Goin, William Leiserson, Sage Moore, Nir Shavit, and Dan Alistarh. 2020\.  Inducing and Exploiting Activation Sparsity for Fast Inference on Deep Neural Networks. In _Proceedings of the 37th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 119)_ , Hal Daumé III and Aarti Singh (Eds.). PMLR, 5533–5543.  https://proceedings.mlr.press/v119/kurtz20a.html
  * Kusupati et al. (2020) Aditya Kusupati, Vivek Ramanujan, Raghav Somani, Mitchell Wortsman, Prateek Jain, Sham Kakade, and Ali Farhadi. 2020.  Soft threshold weight reparameterization for learnable sparsity. In _International Conference on Machine Learning_. PMLR, 5544–5555.
  * Lassance and Clinchant (2022) Carlos Lassance and Stéphane Clinchant. 2022.  An efficiency study for SPLADE models. In _Proceedings of the 45th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 2220–2226.
  * Lemire and Boytsov (2015) Daniel Lemire and Leonid Boytsov. 2015.  Decoding billions of integers per second through vectorization.  _Softw. Pract. Exp._ 45, 1 (2015), 1–29.
  * Lin and Trotman (2015) Jimmy Lin and Andrew Trotman. 2015.  Anytime Ranking for Impact-Ordered Indexes. In _Proceedings of the 2015 International Conference on The Theory of Information Retrieval_ (Northampton, Massachusetts, USA) _(ICTIR ’15)_. Association for Computing Machinery, New York, NY, USA, 301–304.  https://doi.org/10.1145/2808194.2809477
  * Lin and Ma (2021) Jimmy J. Lin and Xueguang Ma. 2021.  A Few Brief Notes on DeepImpact, COIL, and a Conceptual Framework for Information Retrieval Techniques.  _ArXiv_ abs/2106.14807 (2021).
  * LIU et al. (2020) Junjie LIU, Zhe XU, Runbin SHI, Ray C. C. Cheung, and Hayden K.H. So. 2020.  Dynamic Sparse Training: Find Efficient Sparse Network From Scratch With Trainable Masked Layers. In _International Conference on Learning Representations_.
  * Mackenzie et al. (2021) Joel Mackenzie, Matthias Petri, and Alistair Moffat. 2021\.  Anytime Ranking on Document-Ordered Indexes.  _ACM Trans. Inf. Syst._ 40, 1, Article 13 (sep 2021), 32 pages.
  * Mallia et al. (2021) Antonio Mallia, O. Khattab, Nicola Tonellotto, and Torsten Suel. 2021.  Learning Passage Impacts for Inverted Indexes.  _SIGIR_ (2021).
  * Mallia et al. (2022) Antonio Mallia, Joel Mackenzie, Torsten Suel, and Nicola Tonellotto. 2022.  Faster Learned Sparse Retrieval with Guided Traversal. In _Proceedings of the 45th Annual International ACM SIGIR Conference on Research and Development in Information Retrieval (SIGIR 2022)_. 1901–1905.
  * Mallia et al. (2017) Antonio Mallia, Giuseppe Ottaviano, Elia Porciani, Nicola Tonellotto, and Rossano Venturini. 2017\.  Faster BlockMax WAND with Variable-sized Blocks. In _Proc. of the 40th International ACM SIGIR Conference on Research and Development in Information Retrieval_. 625–634.
  * Mallia et al. (2019b) Antonio Mallia, Michal Siedlaczek, Joel Mackenzie, and Torsten Suel. 2019b.  PISA: Performant indexes and search for academia.  _Proceedings of the Open-Source IR Replicability Challenge_ (2019).
  * Mallia et al. (2019a) Antonio Mallia, Michal Siedlaczek, and Torsten Suel. 2019a.  An Experimental Study of Index Compression and DAAT Query Processing Methods. In _Proc. of 41st European Conference on IR Research, ECIR’ 2019_. 353–368.
  * Park and Lee (2020) Jin-Woo Park and Jong-Seok Lee. 2020.  Dynamic Thresholding for Learning Sparse Neural Networks.  In _ECAI 2020_. IOS Press, 1403–1410.
  * Qiao et al. (2022) Yifan Qiao, Yingrui Yang, Haixin Lin, Tianbo Xiong, Xiyue Wang, and Tao Yang. 2022\.  Dual Skipping Guidance for Document Retrieval with Learned Sparse Representations.  _ArXiv_ abs/2204.11154 (April 2022).
  * Qiao et al. (2023) Yifan Qiao, Yingrui Yang, Haixin Lin, and Tao Yang. 2023\.  Optimizing Guided Traversal for Fast Learned Sparse Retrieval. In _Proceedings of the ACM Web Conference 2023_ _(WWW ’23)_. ACM, Austin, TX, USA.
  * Ren et al. (2021) Ruiyang Ren, Yingqi Qu, Jing Liu, Wayne Xin Zhao, QiaoQiao She, Hua Wu, Haifeng Wang, and Ji-Rong Wen. 2021\.  RocketQAv2: A Joint Training Method for Dense Passage Retrieval and Passage Re-ranking. In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_. Association for Computational Linguistics, Punta Cana, Dominican Republic, 2825–2835.
  * Santhanam et al. (2021) Keshav Santhanam, O. Khattab, Jon Saad-Falcon, Christopher Potts, and Matei A. Zaharia. 2021\.  ColBERTv2: Effective and Efficient Retrieval via Lightweight Late Interaction.  _ArXiv_ abs/2112.01488 (16 Dec. 2021).
  * Thakur et al. (2021) Nandan Thakur, Nils Reimers, Andreas Rücklé, Abhishek Srivastava, and Iryna Gurevych. 2021\.  BEIR: A Heterogeneous Benchmark for Zero-shot Evaluation of Information Retrieval Models. In _Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2)_.  https://openreview.net/forum?id=wCu6T5xFjeJ
  * Xiao et al. (2022) Shitao Xiao, Zheng Liu, Weihao Han, Jianjin Zhang, Defu Lian, Yeyun Gong, Qi Chen, Fan Yang, Hao Sun, Yingxia Shao, Denvy Deng, Qi Zhang, and Xing Xie. 2022\.  Distill-VQ: Learning Retrieval Oriented Vector Quantization By Distilling Knowledge from Dense Embeddings.  _Proc. of SIGIR_ (2022).
  * Xiong et al. (2021) Lee Xiong, Chenyan Xiong, Ye Li, Kwok-Fung Tang, Jialin Liu, Paul N. Bennett, Junaid Ahmed, and Arnold Overwijk. 2021\.  Approximate Nearest Neighbor Negative Contrastive Learning for Dense Text Retrieval. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=zeFrfgyZln
  * Yang et al. (2021) Jheng-Hong Yang, Xueguang Ma, and Jimmy Lin. 2021\.  Sparsifying Sparse Representations for Passage Retrieval by Top-k Masking.  _CoRR_ abs/2112.09628 (2021).  arXiv:2112.09628 https://arxiv.org/abs/2112.09628
  * Zamani et al. (2018) Hamed Zamani, Mostafa Dehghani, William Bruce Croft, Erik G. Learned-Miller, and J. Kamps. 2018\.  From Neural Re-Ranking to Neural Ranking: Learning a Sparse Representation for Inverted Indexing.  _Proceedings of the 27th ACM International Conference on Information and Knowledge Management_ (2018).
  * Zhan et al. (2022) Jingtao Zhan, Jiaxin Mao, Yiqun Liu, Jiafeng Guo, Min Zhang, and Shaoping Ma. 2022\.  Learning Discrete Representations via Constrained Clustering for Effective and Efficient Dense Retrieval. In _Proc. of Fifteenth ACM International Conference on Web Search and Data Mining_ _(WSDM ’22)_. ACM, New York, NY, USA, 1328–1336.


# Relaxing graphene plasmon excitation constraints through the use of an
epsilon-near-zero substrate

Vinicius Tadin Alvarenga<EMAIL_ADDRESS>D. A. Bahamon
MackGraphe, Mackenzie Presbyterian Institute, Sao Paulo, SP 01302-907, Brazil
School of Engineering, Mackenzie Presbyterian University, Sao Paulo, SP
01302-907, Brazil  Nuno M. R. Peres  Department and Centre of Physics,
University of Minho, Campus of Gualtar, 4710-057, Braga, Portugal
International Iberian Nanotechnology Laboratory (INL), Braga, Portugal
Christiano J. S. de Matos<EMAIL_ADDRESS>MackGraphe, Mackenzie
Presbyterian Institute, Sao Paulo, SP 01302-907, Brazil  School of
Engineering, Mackenzie Presbyterian University, Sao Paulo, SP 01302-907,
Brazil

###### Abstract

Graphene plasmons have attracted significant attention due to their
tunability, potentially long propagation lengths and ultracompact wavelengths.
However, the latter characteristic imposes challenges to light-plasmon
coupling in practical applications, generally requiring sophisticated coupling
setups, extremely high doping levels and/or graphene nanostructuting close to
the resolution limit of current lithography techniques. Here, we propose and
theoretically demonstrate a method for alleviating such a technological strain
through the use of a practical substrate whose low and negative dielectric
function naturally enlarges the graphene polariton wavelength to more
manageable levels. We consider silicon carbide (SiC), as it exhibits a
dielectric function whose real part is between -1 and 0, while its imaginary
part remains lower than 0.05, in the 951 to 970 cm-1 mid-infrared spectral
range. Our calculations show hybridization with the substrate’s phonon
polariton, resulting in a polariton wavelenth that is an order of magnitude
longer than obtained with a silicon dioxide substrate, while the propagation
length increases by the same amount.

###### pacs:

Valid PACS appear here

##  I Introduction

Plasmons, the collective oscillation of free electrons in a conductive
material, have attracted great attention in the field of nanophotonics due to
its high confinement of electromagnetic energy, achieving subdifractional
spatial extensionsabajoch ; maier , which allows for the construction of
photonic devices at the nanometer scale. The applications range from
biosensingbio1 ; bio2 ; bio3  to enhanced spectroscopy techniques such as
SERSsers1 ; sers2 ; sers3 , TERSters1 ; ters2 ; ters3 ; ters4  and SNOMsnom1 ;
snom2 ; snom3 ; snom5 ; snom5 ; snom6 .

Graphene, a two-dimensional material consisting of a single layer of carbon
atoms in a hexagonal latticegeim-nov-2004 , has emerged as a promising
material for plasmonic applications due to its high electronic mobility and
null bandgap, which allows for high Fermi level tunability and significant
changes on the optical conductivitykoppens ; pol-2d . Since the plasmon
dispersion is dependent on the electronic density, graphene plasmons can be
actively controlled more easily than those in conventional conductorsfang .
Also, graphene can be transfered to a range of dielectric surfaces, allowing
for the integration of graphene-based plasmonic devices with photonic devices,
such as silicon-based waveguideslivro-nuno .

The plasmon wavevector in graphene tends to be up to 2 orders of magnitude
larger than that of free space radiation in the mid-infrared range. However,
efficient radiation-plasmon coupling, for the excitation of surface plasmon
polaritons (SPPs), requires circumventing the large wavevector mismatch. For
bridging such a high mismatch, finding a suitable high-index dielectric
material for prism couplingotto ; kre  becomes impractical, while diffraction
grating coupling has required structures with periodicity of tens of
nanometersgrating1 ; grating2 ; grating3 , posing a strain on the device
fabrication step. Alternatively, localized surface plasmons (LSPs) have been
excited in graphene patterned into nanostructures such as nanoribbons or
nanodisksju ; Rodrigo ; fang2 . Again, due to the large wavevector mismatch,
the features in these structures need to be in the tens to hundreds of
nanometers range, for operation in the mid-infraredRodrigo ; fang ; vasic .
Even though such structures have been fabricated and used to excite graphene
plasmons, the fabrication of such devices is challenging due to their size
being close to the spatial resolution limit of current nanofabrication
techniquesnikitin .

The plasmon dispersion in graphene is dependent on the dielectric function of
the adjacent media. Therefore, the choice of the substrate plays a key role in
the wavevector matching for plasmon excitation. In particular a decrease in
the magnitude of the plasmon wavevector is expected for a substrate with a low
permittivityabajoch ; livro-nuno . The plasmon wavevector is proportional to
the sum of the dielectric functions of the substrate and superstrate. If we
consider the latter to be air, meaning a dielectric function equal to 1, the
optimal choice of substrate would be one with dielectric function of -1.
However, a common drawback of materials with a negative dielectric function is
that the imaginary part of the dielectric function cannot be neglected,
leading to plasmon losses.

Epsilon-near-zero (ENZ) materials have received considerable recent attention
for a number of interesting properties, including increased nonlinear optical
susceptibilitiesAlam , which can be used to enhance nonlinear effects such as
second harmonic generationpilar , and wavelength enlargement, which allows
light tunneling through narrow channelssilveirinha . As in plasmonics, ENZ
applications also require a low imaginary dielectric function component at the
spectral range where the real dielectric function is low. As a consequence,
ENZ materials are also good candidates for plasmonic applications. Indeed,
they have been studied in the field of plasmonics for being capable of
directly coupling SPPs to free-space lighttraviss .

Silicon carbide is an ENZ material in the mid-infrared due to its optical
phononstiwald ; engheta1 ; engheta2 . The range between the transverse optical
phonon (TO), at 797 cm-1, and the longitudinal optical phonon (LO), at 970
cm-1, called the Reststrahlen band, exhibits a real negative dielectric
function, with the function crossing zero at both these wavenumbers.caldwell ;
koch . In particular, a low imaginary part of epsilon is obtained near the LO
phonon wavenumber. Within the Reststrahlen band, the material also supports
surface phonon polaritons (SPhP)koch ; huber ; sic1 ; sic2 . It has been shown
that when graphene plasmons are excited within the Reststrahlen band of a SiC
substrate, hybrid plasmon-phonon modes arise, which maintain the
characteristics of both excitations, exhibiting a longer wavelength than a
pure SPP and keeping losses lowhy1 ; hy2 ; hy3 ; sic3 ; sic4 . The use of SiC
is also convenient since graphene can be epitaxially grown directly onto this
substrateepi1 ; epi2 ; epi3 ; epi4 ; epi5 . Other polar dielectrics have also
been studied as substrates for graphene plasmonics due to their ENZ naturehy1
; hy2 ; hy3 , however, ENZ substrates have not been used for increasing the
graphene plasmon wavelength with the purpose of reducing the wavevector
mismatch for optical excitation. In this work, we propose the use of SiC as a
graphene substrate for mitigating this mismatch in the mid-infrared. Our
results show an increase from hundreds of nanometers to a few micrometers in
polariton wavelength on SiC when compared with a SiO2 substrate. This
simplifies the polaritonic excitation while keeping the polariton modal volume
3 orders of magnitude smaller than that of free-space radiation. Although
there is a loss of confinement, these polaritonic waves become easier to
excite, and the fabrication of polariton-based devices becomes easier. The
conditions for this trade-off can be adjusted by tailoring the size of the
structures used for excitation and controlling the graphene doping.

Figure 1: Dielectric functions for silicon carbide (A) and silicon dioxide
(B). The data for SiC was taken from ref. dong  and for SiO2 from ref. karimi

##  II Methods

To obtain the graphene plasmon dispersion, we used the condition for SPP
formationabajoch :

$\frac{\varepsilon_{1}}{\sqrt{\beta^{2}-k_{0}^{2}\varepsilon_{1}}}+\frac{\varepsilon_{2}}{\sqrt{\beta^{2}-k_{0}^{2}\varepsilon_{2}}}+\frac{i\sigma}{\omega\epsilon_{0}}=0,$
(1)

where $\beta$ is the plasmon wavenumber, $k_{0}$ is the free-space radiation
wavenumber, $\sigma$ is the graphene optical conductivity, $\varepsilon_{1}$,
$\varepsilon_{2}$ are the dielectric functions of the superstrate, substrate,
respectively, and $\epsilon_{0}$ is the vacuum permittivity.

We assume graphene to be on the interface between air ($\varepsilon_{1}$ = 1)
and either SiC or SiO2. We use the Drude model for the graphene
conductivityabajoch ; livro-nuno :

$\sigma(\omega)=\frac{e^{2}E_{F}\tau}{\hbar^{2}\sqrt{\pi}(1-i\tau\omega)},$
(2)

where $e$ is the electron charge, $E_{F}$ is the Fermi energy and $\tau$ is
the relaxation rate. This is a good approximation when the temperature is low
and electronic density is high. Also, non-local effects are neglected. These
are considered important only for small nanostructures and when graphene is
near metallic surfacesnature-nuno . The value of the relaxation time is
assumed to be such that $\hbar\tau^{-1}=2.4$ meVfei1  . The dielectric
function for SiO2 was taken to be given bykarimi .

$\varepsilon_{SiO_{2}}(\omega)=\varepsilon_{\infty,SiO_{2}}\sum_{j}\frac{s_{j}^{2}}{\omega_{j}^{2}-\omega^{2}-i\Gamma_{j}\omega_{j}},$
(3)

with $\varepsilon_{\infty,SiO_{2}}$ = 2.3 and using the parameters shown in
Table 1. The dielectric function of SiC was taken to be given bydong

$\varepsilon_{SiC}(\omega)=\varepsilon_{\infty,SiC}\frac{\omega^{2}-\omega^{2}_{l}+i\Gamma\omega}{\omega^{2}-\omega_{t}^{2}+i\Gamma\omega},$
(4)

with $\varepsilon_{\infty,SiC}=6.7$, $\omega_{l}$ = 120.5 meV, $\omega_{t}$ =
98403 meV and $\Gamma$ = 5.9 meV.

$\omega_{j}$ (meV) |  $s_{j}^{2}$ (meV2) |  $\Gamma_{j}$(meV)
---|---|---
142 | 812 | 7.4
133 | 7832 | 5.4
100 | 537 | 4
57 | 3226 | 6.2
47 | 1069 | 24.5
Table 1:  parameters used for calculating the dielectric function of SiO2
using equation 3

We numerically solved equation 1 for obtaining the in-plane wavevector
$\beta$, from which we can calculate both the polariton wavelength
($\lambda_{sp}=2\pi/Re\\{\beta\\}$) and the effective propagation length
($L_{eff}=1/2Im\\{\beta\\}$). To better understand the plasmon-phonon
interaction, we also calculated the loss function, given by the imaginary part
of the p-polarized reflectivity rplivro-nuno :

$r_{p}=\frac{k_{1}\varepsilon_{2}-k_{2}\varepsilon_{1}+k_{1}k_{2}\sigma/(\varepsilon_{0}\omega)}{k_{1}\varepsilon_{2}+k_{2}\varepsilon_{1}+k_{1}k_{2}\sigma/(\epsilon_{0}\omega)},$
(5)

where $k_{1}$ and $k_{2}$ are the wavevectors of light propagating in air and
the substrate, respectively.

We have also calculated the reflectance and absorbance spectra for an array of
graphene ribbons in order to determine the corresponding LSP resonance
condition. Since the conductivity of a regular array of graphene nanoribbons
is spatially periodic ($\sigma(x+R)=\sigma(x)$, where $R$ is the periodicity),
the reflected and transmitted fields can be expressed by a Fourier-Floquet
series. Matching the boundary conditions at the interface, we arrive at the
linear system of equations: livro-nuno

$\displaystyle\left(\frac{\varepsilon_{1}}{\kappa_{1,n}}+\frac{\varepsilon_{2}}{\kappa_{2,n}}\right)E^{(1)}_{x,n}+$
$\displaystyle\frac{i}{\omega\varepsilon_{0}}\sum_{l}\tilde{\sigma}_{n-l}E^{(1)}_{x,l}=i2\frac{\varepsilon_{2}}{k_{z}}E^{\text{inc}}_{x}\delta_{n,0}$
(6)

where $\beta=|k|\sin{\theta}$, $k_{z}=|k|\cos{\theta}$,
$\kappa_{j,n}=\sqrt{(\beta+nG)^{2}-\varepsilon_{j}\omega^{2}/c^{2}}$ and
$G=2\pi/R$. The Kronecker delta signals that the incident field
$E^{\text{inc}}_{x}$ only appears for the diffraction order $n=0$. The Fourier
coefficients of the conductivity $\tilde{\sigma}_{l}=\int_{0}^{R}\sigma
e^{-iGlx}dx/R$ are calculated assuming that the width of the graphene
nanoribbon in the unit cell is $d_{g}\leq R$, and has the same conductivity of
a graphene sheet ($\sigma$). Once, the fields are obtained the reflectance,
transmittance and absorbance are respectively calculated as:

$\displaystyle\mathcal{T}=\sum_{n}\frac{\varepsilon_{1}k_{z}}{\varepsilon_{2}|\kappa_{1,n}|}\left|\frac{E^{(1)}{x,n}}{E^{\text{inc}}_{x}}\right|^{2}$
(7)
$\displaystyle\mathcal{R}=\sum_{n}\frac{k_{z}}{|\kappa_{2,n}|}\left|\frac{E^{(2)}{x,n}}{E^{\text{inc}}_{x}}\right|^{2}$
(8) $\displaystyle\mathcal{A}=1-\mathcal{T}-\mathcal{R}$ (9)

##  III Results and discussion

Figures 1A and B show the dielectric functions of SiC and SiO2, respectively.
For silicon carbide the real part of the dielectric function is equal to zero
at 970 cm-1. At 951 cm-1 the real dielectric function is -1 and the
corresponding imaginary part is $<0.04$, meaning that polaritonic losses are
expected to be low. In the same spectral region, SiO2 has a real dielectric
function of 4 to 5, with a much higher immaginary part of around 0.4 to 0.8,
meaning that the graphene plasmon wavelength is expected to be shorter and
propagation losses are expected to be higher.

Figures 2A and B respectively show the polariton wavelength and effective
propagation length as functions of free-space wavenumber calculated by solving
eq. 1 for an air/SiC interface with and without graphene with a Fermi energy
of 0.4 eV. We focus this analysis in the spectral region where the real
dielectric function of the substrate is between 0 and -1. Without graphene
this system should not, in principle, exhibit surface phonon polaritons for an
excitation wavenumber beyond 951 cm-1, as
$Re[\varepsilon_{1}+\varepsilon_{2}]>0$, not allowing for polaritons. However,
because of the imaginary part of $\varepsilon_{2}$, a mode arises with very
high losses between 951 and 970 cm-1, where the sum of the real dielectric
functions is close to zero. The addition of graphene to the interface allows
for the excitation of lower loss plasmon-phonon polaritons above 951 cm-1,
which is expressed by a shorter wavelength and a longer propagation length.
For lower wavenumbers, within the Reststrahlen band, graphene imposes only a
minor change to the wavelength, which suggests that in that region the
polariton would exhibit a dominant SPhP character.

Figure 2: Polariton wavelength (A) and effective propagation length (B) on an
air-SiC interface with (red line) and without (black line) graphene. Polariton
wavelength (C) and propagation length (D) for graphene on SiO2 (black) and SiC
(red).Figures A-D consider graphene with a 0.4 eV doping. Polariton wavelength
(E) and propagation loss (F) for graphene on a SiC-air interface for various
doping levels.

Figures 2C and D show a comparison of the polariton wavelength and propagation
length using SiC and SiO2 as the substrate. At the spectral range where the
real dielectric function of SiC lies between -1 and 0, an order of magnitude
increase in polariton wavelength is obtained. At 951 cm-1, for example, the
plasmon wavelength is 0.4 $\mu$m on SiO2 and 5.8 $\mu$m on SiC. As for
propagation length, it is 0.2 $\mu$m for SiO2 and 9 $\mu$m for SiC. Note that
for SiO2, the propagation length is shorter than the polariton wavelentgh,
while for SiC, the propagation length is longer. The combination of longer
wavelength and propagation lengths makes polaritonic devices based on graphene
on SiC simultaneously easier to fabricate and more efficient. Note that
despite the polariton wavelength increase, it remains about half the light
wavelength at the same frequency, allowing for subdiffractional devices to be
designed. The optimal wavelength for a polaritonic device can be tailored by
changing the free-space wavenumber or adjusting graphene’s Fermi level.

Figures 2E and F respectively show the polariton wavelength and propagation
length for various graphene Fermi energies, ranging from 0.05 eV to 0.8 eV. As
expected, $\lambda_{sp}$ increases with EF. At 951 cm-1, for example, Fermi
energies of 0.2, 0.4 and 0.8 eV result in plasmon wavelengths of approximately
4.9, 6 and 7.2 $\mu$m, respectively. Once again we notice that the propagation
length also increases. For the same Fermi energies, we get propagation lengths
of 5.8, 10 and 17.3 $\mu$m. Even for a low graphene doping such as 0.2 eV we
get a much longer polariton wavelength and propagation length, on the order of
micrometers, than with a SiO2 substrate, which results in wavelengths of
hundreds of nanometers, with much higher losses.

In order to easily visualize the wavevector mismatch between free-space
radiation and the polariton wave in graphene on the two different substrates,
we plot, in Figure 3, $\Delta K=|\frac{2\pi}{\lambda_{0}}-Re\\{\beta\\}|$,
where $\lambda_{0}$ is the free-space light wavelength, for both SiC and SiO2
as substrates. We consider a graphene with 0.4 eV Fermi energy. A decrease of
more than an order of magnitude in $\Delta$K is obtained by using SiC. The
lower wavevector mismatch between polariton and free-space radiation makes the
fabrication of polaritonic devices simple. Even if we consider a graphene with
a low doping of 0.2 eV, as our results show, the polariton wavelength is 4.9
$\mu$m, meaning an excitation grating, for example, with a periodicity of 2.45
$\mu$m (considering it to be half the wavelength). Such a feature size is
easily fabricated with simple photolitography, rather than electron beam
lithography, commonly used for conventional graphene plasmonic devices.

Figure 3: Phase mismatch between free-space propagating light and surface
polariton for a graphene sheet at a SiC-air (red line) and a SiO2-air (black
line) interface with a 0.4 eV Fermi energy.

.

Since patterned graphene is one of the most convenient ways to build plasmonic
devices, we now analyze the case of localized surface plasmons in
microstructured graphene. We consider ribbon widths, $d_{g}$, from 0.5 to 2
µm, and fix the ribbon array periodicity at $R=2D$. Figures 4A and B show the
reflectance and absorbance spectra for various ribbon widths with a 0.4 eV
Fermi energy. Since silicon carbide has a negative dielectric function below
970 cm-1, for lower wavenumbers the material has a high reflectance, which is
seen as a step in reflectance spectra. The polaritonic resonance is observed
as a localized dip in reflectance and as peaks in absorbance. The resonance
point can be tuned by adjusting the ribbon width and for excitation at 951
cm-1 2 $\mu$m ribbons are required. With $d_{g}=1.0$ and 0.5 $\mu$m the
polariton resonance is at 951 and 964 cm-1, respectively. Figures 4C and D
show the same calculations for different Fermi energies considering 1 $\mu$m
wide ribbons. Fermi energies of 0.2, 0.4 and 0.8 eV present polaritonic
resonances at 952, 954 and 961 cm-1 respectively. Our calculations, therefore,
show that the resonances in this spectral region can be excited using
structures on the scale of micrometers, which alleviate the strain on
fabrication. The absorbance peak that remains still ($\sim
970~{}\text{cm}^{-1}$) while the doping level and period of the graphene array
are modified is produced by bulk evanescent waves in the substrate (SiC). It
is important to notice that this peak is present even without the excitation
of the surface polaritons.C8NR01706A

Figure 4: Polariton resonance calculations for graphene at an SiC-air
interface. Reflectance (B) and absorbance (C) spectra for different graphene
ribbon widths considering 0.4 eV Fermi energy. Reflectance (D) and absorbance
(E) spectra for different Fermi energies considering 1$\mu$m wide graphene
ribbons.

Figure 5 shows the calculated loss function for graphene at an air-SiC
interface at various doping levels, through which the polariton dispersion
relation can be visualized as bright yellow/white lines. Figures 5A-C show the
polariton dispersion for graphene with a Fermi energy of 0.05, 0.2 and 0.4 eV,
respectively. If there was no graphene at the interface, the only polariton
solution would correspont to a straight horizontal line at the SiC’s SPhP
frequency (represented by the dashed green line at 5D). When graphene is added
to the interface, plasmon-phonon polariton hibridization takes place upon an
anti-crossing of the dispersion relations, yielding 2 polariton branches,
which shift towards higher free-space wavenumbers as the graphene Fermi level
increases. As a general rule, the hybridized modes behave more like an SPP for
low free-space wavenumbers and more like an SPhP for higher free-space
wavenumberslivro-nuno . Figure 5D shows a zoom of Figure 5C at the area
indicated by the yellow rectangle. The dashed cyan lines represent the
wavenumbers of the LO and TO SiC phonons, the dashed green line represents the
point where Re$\\{\varepsilon\\}$ = -1, which is also the point where SPhP
would appear for SiC without graphene and the magenta line represents light
propagation through silicon carbide. The long polariton wavelengths reported
here correspond to the upper hibrid polariton branch.

It should be noted that similar hybridization takes place near the optical
phonon wavenumbers of other polar dielectrics used as substrates. In fact, the
same effect can be observed in SiO2 at 1184 cm-1, where its dielectric
function equals -1. However, we calculate that the polariton wavelength
reaches a more moderate value of 1.2 $\mu$m, with a corresponding sub-
wavelength propagation length of 0.9 $\mu$m. When compared with SiC, the
shorter polariton wavelength and propagation lengths are a consequence of
SiO2’s imaginary dielectric function, which remains high, 0.6 at the 1184
cm-1, at the epsilon near zero point. Therefore we can conclude that, although
polariton wavelength increase can be achieved for other polar dielectrics, SiC
is a more efficient substrate due to its lower losses.

Figure 5: Loss function for graphene on a SiC/air interface with doping levels
of 0.05 eV (A), 0.2 eV (B) and 0.4 eV (C). (D) zoom of (C) in the area
indicated by the yellow square. The dashed cyan lines indicate the wavenumbers
of SiC’s longitudinal and transversal optical phonons. The dashed green line
represents the wavenumber at which the real dielectric function is equal to
-1, which is also the SiC’s SPhP frequency. The magenta line represents the
light line in SiC.

##  IV Conclusions

We have demonstrated a method for increasing the wavelength of graphene
polaritons in order to relax the conditions for matching their wavevector with
that of free-space radiation. For this purpose, we choose a substrate with a
real dielectric function between -1 and 0 in the mid-infrared region: silicon
carbide. In the 951 to 970 cm-1 spectral range, we observe an order of
magnitude polariton wavalength increase relative to SiO2 as a substrate, which
is accompanied by an order of magnitude increase in propagation length. Due to
the Reststrahlen band formed by phonons in the crystal, graphene plasmons
hybridize with SiC phonons, allowing for the wavelength increase, which makes
the excitation of polaritons with launching gratings easier. Along the same
lines graphene structures of up to 2 $\mu$m can be used for localized surface
plasmon excitation. These result confirm SiC as a convenient substrate for
graphene plasmonics, allowing for simple and large scale microfabrication
methods to be employed, which can contribute to the more widespread use of
graphene-based subdifractional photonics in real world applications.

##  V Acknowledgements

This work was funded by FAPESP (grant nos. 2015/11779-4, 2018/07276-5 and
2018/25339-4), the Brazilian Nanocarbon Institute of Science and Technology
(INCT/Nanocarbon), CAPES-PRINT (grant no. 88887.310281/2018-00) and CAPES-
PROSUC. NMRP acknowledges support by the Portuguese Foundation for Science and
Technology (FCT) in the framework of the Strategic Funding UIDB/04650/2020,
support from the European Commission through the project “Graphene-Driven
Revolutions in ICT and Beyond” (Ref. No. 881603, CORE 3), COMPETE 2020,
PORTUGAL 2020, FEDER and the FCT through project POCI-01-0145-FEDER-028114.

## References

  * [1] F. Javier García de Abajo.  Graphene plasmonics: Challenges and opportunities.  ACS Photonics, 1(3):135–152, 2014.
  * [2] SA Maier.  Plasmonics: Fundamentals And Applications, page 245.  01 2007.
  * [3] Chengjun Huang, Jian Ye, Shuo Wang, Tim Stakenborg, and Liesbet Lagae.  Gold nanoring as a sensitive plasmonic biosensor for on-chip dna detection.  Applied Physics Letters, 100(17):173114, 2012.
  * [4] Lin Wu, Hong-Son Chu, Wee Shing Koh, and Er-Ping Li.  Highly sensitive graphene biosensors based on surface plasmon resonance.  Optics express, 18:14395–400, 07 2010.
  * [5] Shuwen Zeng, Kandammathe Valiyaveedu Sreekanth, Jingzhi Shang, Ting Yu, Chih-Kuang Chen, Feng Yin, Dominique Baillargeat, Philippe Coquet, Ho-Pui Ho, Andrei V. Kabashin, and Ken-Tye Yong.  Graphene–gold metasurface architectures for ultrasensitive plasmonic biosensing.  Advanced Materials, 27(40):6163–6169, 2015.
  * [6] Kyle C. Bantz, Audrey F. Meyer, Nathan J. Wittenberg, Hyungsoon Im, Özge Kurtuluş, Si Hoon Lee, Nathan C. Lindquist, Sang-Hyun Oh, and Christy L. Haynes.  Recent progress in sers biosensing.  Phys. Chem. Chem. Phys., 13:11551–11567, 2011.
  * [7] Shan Cong, Zhen Wang, Wenbin Gong, Zhigang Chen, Weibang Lu, John Lombardi, and Zhigang Zhao.  Electrochromic semiconductors as colorimetric sers substrates with high reproducibility and renewability.  Nature Communications, 10, 02 2019.
  * [8] Hiang Kwee Lee, Yih Hong Lee, Charlynn Sher Lin Koh, Gia Chuong Phan-Quang, Xuemei Han, Chee Leng Lay, Howard Yi Fan Sim, Ya-Chuan Kao, Qi An, and Xing Yi Ling.  Designing surface-enhanced raman scattering (sers) platforms beyond hotspot engineering: emerging opportunities in analyte manipulations and hybrid materials.  Chem. Soc. Rev., 48:731–756, 2019.
  * [9] Alyssa B. Zrimsek, Naihao Chiang, Michael Mattei, Stephanie Zaleski, Michael O. McAnally, Craig T. Chapman, Anne-Isabelle Henry, George C. Schatz, and Richard P. Van Duyne.  Single-molecule chemistry with surface- and tip-enhanced raman spectroscopy.  Chemical Reviews, 117(11):7583–7613, 2017.  PMID: 28610424.
  * [10] Prabhat Verma.  Tip-enhanced raman spectroscopy: Technique and recent advances.  Chemical Reviews, 117(9):6447–6466, 2017.  PMID: 28459149.
  * [11] Elena Bailo and Volker Deckert.  Tip-enhanced raman scattering.  Chem. Soc. Rev., 37:921–930, 2008.
  * [12] Weitao Su, Naresh Kumar, Andrey Krayev, and Marc Chaigneau.  In situ topographical chemical and electrical imaging of carboxyl graphene oxide at the nanoscale.  Nature Communications, 9, 12 2018.
  * [13] Peter Hermann, Bernd Kästner, Arne Hoehl, Vyacheslavs Kashcheyevs, Piotr Patoka, Georg Ulrich, Jörg Feikes, Markus Ries, Tobias Tydecks, Burkhard Beckhoff, Eckart Rühl, and Gerhard Ulm.  Enhancing the sensitivity of nano-ftir spectroscopy.  Opt. Express, 25(14):16574–16588, Jul 2017.
  * [14] Sergiu Amarie, Paul Zaslansky, Yusuke Kajihara, Griesshaber Erika, Wolfgang Schmahl, and Fritz Keilmann.  Nano-ftir chemical mapping of minerals in biological materials.  Beilstein journal of nanotechnology, 3:312–23, 04 2012.
  * [15] Alexandre Dazzi and Craig B. Prater.  Afm-ir: Technology and applications in nanoscale infrared spectroscopy and chemical imaging.  Chemical Reviews, 117(7):5146–5173, 2017.  PMID: 27958707.
  * [16] Hans A. Bechtel, Eric A. Muller, Robert L. Olmon, Michael C. Martin, and Markus B. Raschke.  Ultrabroadband infrared nanospectroscopic imaging.  Proceedings of the National Academy of Sciences, 111(20):7191–7196, 2014.
  * [17] Ingrid D. Barcelos, Hans A. Bechtel, Christiano J. S. de Matos, Dario A. Bahamon, Bernd Kaestner, Francisco C. B. Maia, and Raul O. Freitas.  Probing polaritons in 2d materials with synchrotron infrared nanospectroscopy.  Advanced Optical Materials, 8(5):1901091, 2020.
  * [18] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov.  Electric field effect in atomically thin carbon films.  Science, 306(5696):666–669, 2004.
  * [19] Frank H. L. Koppens, Darrick E. Chang, and F. Javier García de Abajo.  Graphene plasmonics: A platform for strong light–matter interactions.  Nano Letters, 11(8):3370–3377, 2011.  PMID: 21766812.
  * [20] D. N. Basov, M. M. Fogler, and F. J. Garcia de Abajo.  Polaritons in van der waals materials.  Science, 354(6309), 2016.
  * [21] Zheyu Fang, Sukosin Thongrattanasiri, Andrea Schlather, Zheng Liu, Lulu Ma, Yumin Wang, Pulickel M. Ajayan, Peter Nordlander, Naomi J. Halas, and F. Javier García de Abajo.  Gated tunability and hybridization of localized plasmons in nanostructured graphene.  ACS Nano, 7(3):2388–2395, 2013.  PMID: 23390960.
  * [22] Paulo André Gonçalves and Nuno Peres.  An Introduction to Graphene Plasmonics.  06 2016.
  * [23] Andreas Otto.  Excitation of nonradiative surface plasma waves in silver by the method of frustrated total reflection.  Zeitschrift für Physik A Hadrons and nuclei, 216(4):398–410, Aug 1968.
  * [24] E. Kretschmann and H. Raether.  Notizen: Radiative Decay of Non Radiative Surface Plasmons Excited by Light.  Zeitschrift Naturforschung Teil A, 23:2135–2136, December 1968\.
  * [25] Xiaolong Zhu, Wei Yan, Peter Uhd Jepsen, Ole Hansen, N. Asger Mortensen, and Sanshui Xiao.  Experimental observation of plasmons in a graphene monolayer resting on a two-dimensional subwavelength silicon grating.  Applied Physics Letters, 102(13):131101, 2013.
  * [26] Tatiana L. Zinenko, Akira Matsushima, and Alexander I. Nosich.  Surface-plasmon, grating-mode, and slab-mode resonances in the h- and e-polarized thz wave scattering by a graphene strip grating embedded into a dielectric slab.  IEEE Journal of Selected Topics in Quantum Electronics, 23(4):1–9, 2017.
  * [27] Weiwei Luo, Wei Cai, Yinxiao Xiang, Lei Wang, Mengxin Ren, Xinzheng Zhang, and Jingjun Xu.  Flexible modulation of plasmon-induced transparency in a strongly coupled graphene grating-sheet system.  Opt. Express, 24(6):5784–5793, Mar 2016.
  * [28] Long Ju, Baisong Geng, Jason Horng, Caglar Girit, Michael Martin, Zhao Hao, Hans Bechtel, Xiaogan Liang, Alex Zettl, Y Shen, and Feng Wang.  Graphene plasmonics for tunable terahertz metamaterials.  Nature nanotechnology, 6:630–4, 09 2011.
  * [29] Daniel Rodrigo, Odeta Limaj, Davide Janner, Dordaneh Etezadi, F. Javier Garcia de Abajo, Valerio Pruneri, and Hatice Altug.  Mid-infrared plasmonic biosensing with graphene.  Science, 349(6244):165–168, 2015.
  * [30] Zheyu Fang, Yumin Wang, Andrea E. Schlather, Zheng Liu, Pulickel M. Ajayan, F. Javier García de Abajo, Peter Nordlander, Xing Zhu, and Naomi J. Halas.  Active tunable absorption enhancement with graphene nanodisk arrays.  Nano Letters, 14(1):299–304, 2014.  PMID: 24320874.
  * [31] Borislav Vasić, Goran Isić, and Radoš Gajić.  Localized surface plasmon resonances in graphene ribbon arrays for sensing of dielectric environment at infrared frequencies.  Journal of Applied Physics, 113(1):013110, 2013.
  * [32] A. Yu. Nikitin, P. Alonso-González, and R. Hillenbrand.  Efficient coupling of light to graphene plasmons by compressing surface polaritons with tapered bulk materials.  Nano Letters, 14(5):2896–2901, 2014.  PMID: 24773123.
  * [33] M. Zahirul Alam, Israel De Leon, and Robert W. Boyd.  Large optical nonlinearity of indium tin oxide in its epsilon-near-zero region.  Science, 352(6287):795–797, 2016.
  * [34] Pilar G. Vianna, Aline dos S. Almeida, Rodrigo M. Gerosa, Dario A. Bahamon, and Christiano J. S. de Matos.  Second-harmonic generation enhancement in monolayer transition-metal dichalcogenides by using an epsilon-near-zero substrate.  Nanoscale Adv., 3:272–278, 2021.
  * [35] Mario Silveirinha and Nader Engheta.  Tunneling of electromagnetic energy through subwavelength channels and bends using $\epsilon$-near-zero materials.  Phys. Rev. Lett., 97:157403, Oct 2006.
  * [36] Daniel Traviss, Roman Bruck, Ben Mills, Martina Abb, and Otto L. Muskens.  Ultrafast plasmonics using transparent conductive oxide hybrids in the epsilon-near-zero regime.  Applied Physics Letters, 102(12):121112, 2013.
  * [37] Thomas E. Tiwald, John A. Woollam, Stefan Zollner, Jim Christiansen, R. B. Gregory, T. Wetteroth, S. R. Wilson, and Adrian R. Powell.  Carrier concentration and lattice absorption in bulk and epitaxial silicon carbide determined using infrared ellipsometry.  Phys. Rev. B, 60:11464–11474, Oct 1999.
  * [38] Jongbum Kim, Aveek Dutta, Gururaj V. Naik, Alexander J. Giles, Francisco J. Bezares, Chase T. Ellis, Joseph G. Tischler, Ahmed M. Mahmoud, Humeyra Caglayan, Orest J. Glembocki, Alexander V. Kildishev, Joshua D. Caldwell, Alexandra Boltasseva, and Nader Engheta.  Role of epsilon-near-zero substrates in the optical response of plasmonic antennas.  Optica, 3(3):339–346, Mar 2016.
  * [39] I. Liberal and Nader Engheta.  Near-zero refractive index photonics.  Nature Photonics, 11:149–158, 03 2017.
  * [40] Joshua D. Caldwell, Orest J. Glembocki, Yan Francescato, Nicholas Sharac, Vincenzo Giannini, Francisco J. Bezares, James P. Long, Jeffrey C. Owrutsky, Igor Vurgaftman, Joseph G. Tischler, Virginia D. Wheeler, Nabil D. Bassim, Loretta M. Shirey, Richard Kasica, and Stefan A. Maier.  Low-loss, extreme subdiffraction photon confinement via silicon carbide localized surface phonon polariton resonators.  Nano Letters, 13(8):3690–3697, 2013.  PMID: 23815389.
  * [41] R. J. Koch, Th. Seyller, and J. A. Schaefer.  Strong phonon-plasmon coupled modes in the graphene/silicon carbide heterosystem.  Phys. Rev. B, 82:201413, Nov 2010.
  * [42] A. Huber, N. Ocelic, D. Kazantsev, and R. Hillenbrand.  Near-field imaging of mid-infrared surface phonon polariton propagation.  Applied Physics Letters, 87(8):081103, 2005.
  * [43] Miles V. Klein, B. N. Ganguly, and Priscilla J. Colwell.  Theoretical and experimental study of raman scattering from coupled lo-phonon-plasmon modes in silicon carbide.  Phys. Rev. B, 6:2380–2388, Sep 1972.
  * [44] M. Chafai, A. Jaouhari, A. Torres, R. Anton, E. Martin, J. Jimenez, and W. C. Mitchel.  Raman scattering from lo phonon-plasmon coupled modes and hall-effect in n-type silicon carbide 4h–sic.  Journal of Applied Physics, 90(10):5211–5215, 2001.
  * [45] Ke Li, Jamie M. Fitzgerald, Xiaofei Xiao, Joshua D. Caldwell, Cheng Zhang, Stefan A. Maier, Xiaofeng Li, and Vincenzo Giannini.  Graphene plasmon cavities made with silicon carbide.  ACS Omega, 2(7):3640–3646, 2017.
  * [46] Patrick Rufangura, Thomas G Folland, Arti Agrawal, Joshua D Caldwell, and Francesca Iacopi.  Towards low- loss on-chip nanophotonics with coupled graphene and silicon carbide: a review.  Journal of Physics: Materials, 3(3):032005, jul 2020.
  * [47] Yue Li, Iñigo Liberal, and Nader Engheta.  Structural dispersion–based reduction of loss in epsilon-near-zero and surface plasmon polariton waves.  Science Advances, 5(10), 2019.
  * [48] R. J. Koch, S. Fryska, M. Ostler, M. Endlich, F. Speck, T. Hänsel, J. A. Schaefer, and Th. Seyller.  Robust phonon-plasmon coupling in quasifreestanding graphene on silicon carbide.  Phys. Rev. Lett., 116:106802, Mar 2016.
  * [49] Yu Liu and R. F. Willis.  Plasmon-phonon strongly coupled mode in epitaxial graphene.  Phys. Rev. B, 81:081406, Feb 2010.
  * [50] Walt A. de Heer, Claire Berger, Xiaosong Wu, Phillip N. First, Edward H. Conrad, Xuebin Li, Tianbo Li, Michael Sprinkle, Joanna Hass, Marcin L. Sadowski, Marek Potemski, and Gérard Martinez.  Epitaxial graphene.  Solid State Communications, 143(1):92 – 100, 2007.  Exploring graphene.
  * [51] Konstantin Emtsev, Aaron Bostwick, Karsten Horn, Johannes Jobst, Gary Kellogg, L. Ley, J. Mcchesney, Taisuke Ohta, Sergey Reshanov, Jonas Röhrl, Eli Rotenberg, Andreas Schmid, Daniel Waldmann, Heiko Weber, and Thomas Seyller.  Towards wafer-size graphene layers by atmospheric pressure graphitization of silicon carbide.  Nature materials, 8:203–7, 03 2009.
  * [52] Zhen-Yu Juang, Chih-Yu Wu, Chien-Wei Lo, Wei-Yu Chen, Chih-Fang Huang, Jenn-Chang Hwang, Fu-Rong Chen, Keh-Chyang Leou, and Chuen-Horng Tsai.  Synthesis of graphene on silicon carbide substrates at low temperature.  Carbon, 47(8):2026 – 2031, 2009.
  * [53] Sergey Kopylov, Alexander Tzalenchuk, Sergey Kubatkin, and Vladimir I. Fal’ko.  Charge transfer between epitaxial graphene and silicon carbide.  Applied Physics Letters, 97(11):112109, 2010.
  * [54] Walt A. de Heer, Claire Berger, Ming Ruan, Mike Sprinkle, Xuebin Li, Yike Hu, Baiqian Zhang, John Hankinson, and Edward Conrad.  Large area and structured epitaxial graphene produced by confinement controlled sublimation of silicon carbide.  Proceedings of the National Academy of Sciences, 108(41):16900–16905, 2011.
  * [55] Jian Dong, J. Zhao, and L. H. Liu.  Near-field radiative heat transfer between clusters of dielectric nanoparticles.  Journal of Quantitative Spectroscopy and Radiative Transfer, 197:114–122, 08 2017.
  * [56] F. Karimi, A. H. Davoody, and I. Knezevic.  Dielectric function and plasmons in graphene: A self-consistent-field calculation within a markovian master equation formalism.  Phys. Rev. B, 93:205421, May 2016.
  * [57] Paulo André Gonçalves, Thomas Christensen, Nuno Peres, Antti-Pekka Jauho, Itai Epstein, Frank Koppens, Marin Soljačić, and N. Mortensen.  Quantum surface-response of metals revealed by acoustic graphene plasmons.  Nature Communications, 12, 06 2021.
  * [58] Zhe Fei, Aleksandr Rodin, G. Andreev, Wenzhong Bao, Alexander McLeod, M Wagner, L Zhang, Z Zhao, M Thiemens, Gerardo Dominguez, Michael Fogler, Antonio Castro Neto, Jeanie Lau, Fritz Keilmann, and D Basov.  Gate-tuning of graphene plasmons revealed by infrared nano-imaging.  Nature, 487:82–5, 06 2012.
  * [59] Dongxue Chen, Jianjie Dong, Jianji Yang, Yilei Hua, Guixin Li, Chuanfei Guo, Changqing Xie, Ming Liu, and Qian Liu.  Realization of near-perfect absorption in the whole reststrahlen band of sic.  Nanoscale, 10:9450–9454, 2018.