# Activity Detection for Massive Random Access using Covariance-based Matching
Pursuit

Leatile Marata, Member, IEEE, Esa Ollila, Senior Member, IEEE, and Hirley
Alves, Member, IEEE Leatile Marata, and Hirley Alves are with Centre for
Wireless Communications – Radio Technologies, FI-90014, University of Oulu,
Finland. e-mail: {leatile.marata<EMAIL_ADDRESS>Esa Ollila is with
the Department of Information and Communications Engineering, Aalto
University, FI-00076 Aalto, Finland (e-mail: esa.ollila@aalto.fi)This work is
supported by the Research Council of Finland (former Academy of Finland)
(Grants n.319485, n.340171, n.346208 (6G Flagship), n.338408, n.353093,
n.359848). Leatile Marata’s work was partly supported by the Riitta ja Jorma
J. Takanen Foundation, the Nokia Foundation, the Finnish Foundation for
Technology Promotion, and the Botswana International University of Science and
Technology.

###### Abstract

The Internet of Things paradigm heavily relies on a network of a massive
number of machine-type devices (MTDs) that monitor changes in various
phenomena. Consequently, MTDs are randomly activated at different times
whenever a change occurs. This essentially results in relatively few MTDs
being active simultaneously compared to the entire network, resembling
targeted sampling in compressed sensing. Therefore, signal recovery in
machine-type communications is addressed through joint user activity detection
and channel estimation algorithms built using compressed sensing theory.
However, most of these algorithms follow a two-stage procedure in which a
channel is first estimated and later mapped to find active users. This
approach is inefficient because the estimated channel information is
subsequently discarded. To overcome this limitation, we introduce a novel
covariance-learning matching pursuit algorithm that bypasses explicit channel
estimation. Instead, it focuses on estimating the indices of the active users
greedily. Simulation results presented in terms of probability of miss
detection, exact recovery rate, and computational complexity validate the
proposed technique’s superior performance and efficiency.

###### Index Terms:

Activity detection, compressed sensing, covariance-learning, grant-free, NOMA,
matching pursuit, sporadic activation.

##  I Introduction

Machine-type communications (MTC) is a framework that enables nonhuman devices
to exchange information with each other without requiring human intervention
[1, 2, 3]. This kind of communication is predominantly uplink-driven and is
characterized by devices transmitting their data to a central point such as a
base station (BS) in the cellular Internet of Things (IoT). Indeed, this is
particularly crucial for the IoT, where many machine-type devices (MTDs),
e.g., sensors and actuators work hand in hand to collect measurements of
different phenomena and perform automatic tasks [1, 4]. For instance, the
sixth generation (6G) wireless communications standard is expected to host a
minimum ratio of $10^{8}$ MTDs per square kilometer. This makes it necessary
to develop novel efficient signal processing techniques as well as adapt old
ones to work in MTC [5]. Given the massive communication scale, next-
generation multiple access techniques are fundamental in achieving ubiquitous
connectivity [6]. Such techniques include the grant-free (GF) non-orthogonal
multiple access (NOMA) techniques, which have been widely adopted as a key
enabler for MTC [1, 2]. These schemes bypass the conventional handshaking
procedure by permitting uncoordinated transmissions and, as such, leading to
reduced network access latency [7]. Despite this, such latency reduction comes
at the expense of an increased number of transmitted data collisions, which
are ultimately detrimental to the performance of BS receivers [8]. In the end,
the efficiency of MTC (and consequently IoT) networks is subject to how well
the collisions in media access can be resolved. To address this, Liu et al. in
[9] proposed the compressed sensing (CS)- multi-user detection (MUD) framework
for MTC receivers. For instance, the CS-MUD frameworks in [10, 11, 12] follow
a two-stage procedure where the sparse channel matrix is first estimated,
followed by a thresholding step. The thresholding step facilitates mapping
each row of the channel matrix to either a one or zero, thus revealing whether
an MTD is active or inactive, respectively. Among the most widely adopted
solutions are the variants of approximate message passing (AMP), sparse
Bayesian learning (SBL), and simultaneous orthogonal matching pursuit (SOMP)
[13].

Though the two-stage signal recovery approaches mentioned in the previous
paragraph are applicable in many cases, they are inefficient when the main
interest is to identify the active devices and not to estimate the channels of
the MTDs. For example, applications like industrial process monitoring
typically require collective decision-making based on averages of some
measurements. As such, estimating indices of hardwired codebook entries for
each MTD is typically sufficient, thereby eliminating the need for channel
estimation [14, 15, 16]. Equally important is to note that the counterpart of
the two-stage signal recovery called non-coherent detection makes it necessary
to estimate the indices explicitly [17]. Moreover, as further deployments are
being undertaken, an overarching challenge in MTC is not just algorithms that
solve collisions but also those that are efficient [18]. As such, bypassing
explicit channel estimation can facilitate the development of faster
algorithms that use less computational resources. To this end, our work is
focused on devising a novel receiver algorithm that facilitates optimal active
user identification by estimating the active user indices in MTC. To bring
context, we summarise the state-of-the-art works related to our work in the
next sub-section.

###  I-A State of the Art

Research on the problem of active user detection (AUD) and channel estimation
for MTC has been gaining popularity due to its importance in enabling IoT.
There are many approaches to solving this problem. For example, the Bayesian
formulation that relies on the definition of probability distributions of the
variables of interest such as sparsity-inducing parameters has been presented
in [13, 9]. Notably, Zhang et al. proposed an AMP algorithm combined with
variational Bayesian inference for joint user detection and data recovery in
[19]. It is worth mentioning that in contrast with SBL, algorithms developed
following the AMP framework are highly sensitive to non-Gaussian sensing
matrices, thus limiting their practical implementation. Meanwhile, Zhang et
al. proposed a fast inversion SBL for GF-NOMA for MTC in [20], aimed at
reducing the computational complexity of the classical SBL. On the other hand,
Wang et al. proposed a variational Bayesian inference algorithm in [21] for
joint user activity detection in mMTC that is supported by a cloud-based radio
access network.

Another school of thought to solving the AUD problem follows the convex
optimization framework, which uses the $\ell_{1}$-norm minimization. Among
some of the works in this direction include [22], where Djelouat et al.
proposed the alternating direction method of multipliers (ADMM) for device
activity detection in MTC that has correlated channels. Sun et al. also
proposed a joint user activity detection optimization algorithm in [23] by
exploiting the joint sparsity of the sensing matrix and the effective channels
of the MTDs. Similarly, Zhu et al. proposed a joint user activity detection in
[24], using convex optimization relying on atomic norm minimization. In
addition to the aforementioned approaches, the covariance-based activity
detection approach has also gained popularity due to its appealing scalability
[25, 15, 26]. For instance, Liu et al. proposed a maximum likelihood
estimation (MLE) problem in [27] for device activity detection under channels
with the line of sight (LoS). The proposed solution is shown to reduce the
probability of error by $50\%$ with reduced computation times. On the other
hand, Fengler et al. proposed a low complexity MLE algorithm in [28], together
with the scaling law under which the problem can be solved. Meanwhile, You et
al. proposed a covariance-based algorithm for unsourced random access in [29].
On the other hand, Yu et al. in [30] proposed sufficient conditions for
solving the maximum likelihood problem via a covariance-based approach.

Despite all the solutions discussed above, most of them first estimate the
effective channels induced by sporadic activation and not the explicit
estimation of the active users, except for [27, 15]. Unfortunately, works [27,
15] have inherent limitations of the coordinate descent minimization as they
rely on iterative solutions and thus can suffer from slow convergence or
convegence to poor local solution. Moreover, most existing algorithms, save
for SBL, require some form of parameter tuning and are not robust to some
pilot sequences, thus limiting their applicability. The downside of SBL is
that it has a prohibitively high computational complexity that is unsuitable
for mMTC. We therefore endeavor to propose a novel solution for receivers in
MTC using a greedy search approach. Contrary to coordinate descent approaches,
greedy search is not subject to convergence issues and thus yields faster
solutions. In the end, the proposed solution yields results under different
sensing matrices and does not require parameter tuning, just like SBL, yet has
very low complexity, as we present next.

###  I-B Contributions

We consider a mMTC scenario where sporadically activated MTDs transmit their
data in the uplink and devise a covariance-learning (CL) matching pursuit (MP)
algorithm, termed CL-MP, to recover the indices of active devices. The work
departs from classical covariance-based techniques that rely on coordinate-
wise optimization (CWO) algorithms such as [25, 28]. Our main contributions
are summarized as follows:

  * •

We formulate the device activity detection problem using the Gaussian negative
log-likelihood function and solve it using the CL-MP algorithm. It is worth
noting that SBL, CWO, and CL-MP have similar cost functions. Consequently, the
CL-MP algorithm offers the same performance in terms of user identification as
SBL and CWO. However, CL-MP is more efficient compared to SBL and CWO.

  * •

We evaluate the CL-MP algorithm using various pilot sequences, including
structured pilot sequences proposed in [31, 4], as well as Bernoulli and
Gaussian pilot sequences. While the performance may vary across different
pilot sequences, CL-MP consistently yields reasonable results under these
different pilot sequences.

  * •

Finally, we compare our proposed solution with existing techniques such as
SOMP, SBL, and VAMP. Results reveal the superior performance of CL-MP compared
to SOMP and VAMP while having comparable performance to SBL. However, the
runtime results show that CL-MP is faster than SBL.

###  I-C Organization and Notation

The remainder of this paper is organized as follows: Section II introduces the
system model and the AUD problem. In Section III, we present the solution of
the device activity detection using CL-MP. Lastly, in Section IV, we present
the results and discussions, while in Section V, we summarize the paper and
provide some future research directions.

Notation: boldface lowercase and boldface uppercase letters denote column
vectors and matrices, respectively. Moreover, for a matrix $\bm{A}$,
${\bm{a}}_{i}$ is the $i^{\text{th}}$ column, while $a_{ij}$ denotes the
element in the $i^{\text{th}}$ row, $j^{\text{th}}$ column, respectively,
whereas $a_{i}$ denotes the $i^{\text{th}}$ element of a vector ${\bm{a}}$. We
denote the transpose and conjugate transpose operations by superscripts
$(\cdot)^{\mbox{\scriptsize T}}$ and $(\cdot)^{\mbox{\scriptsize H}}$,
respectively, while $\mathbb{C}$ and $\mathbb{R}$ refer to complex and real
domains, respectively. The circularly symmetric $M$-variate complex Gaussian
distribution with mean ${\bm{a}}$ and covariance matrix $\bm{\Sigma}$ is
denoted by $\mathcal{CN}_{\\!M}({\bm{a}},\bm{\Sigma})$. Let
$\operatorname{diag}(\bm{a})$ denote a diagonal matrix whose diagonal elements
are vector ${\bm{a}}$. Meanwhile $\bm{I}$ is the identity matrix, while
$(\cdot)_{+}=\max(\cdot,0)$ is a non-negativity function, and
$[\\![M]\\!]=\\{1,\ldots,M\\}$.

##  II SYSTEM MODEL

Figure 1: A massive machine-type communications network scenario, where $N$
MTDs transmit in the uplink to a BS equipped with $M$ antennas.

We consider a communication network scenario depicted by Figure 1, where a BS
equipped with $M$ antennas serves a set of $N$ MTDs that transmit data in the
uplink. To align with practical MTC, we assume that in a given coherence
interval of $T$ symbols, an unknown set $\mathcal{K}\subset[\\![N]\\!]$ of
cardinality $K=|\mathcal{K}|$ is active. Note that from these $T$ symbols of
the coherence interval, $L$ symbols are used for metadata processing to
identify the $K$ active MTDs. Meanwhile, the remaining $T-L$ symbols are used
to transmit their payload data. Furthermore, we assume that each
$n^{\text{th}}$ MTD is connected to the BS through a quasi-static fading
uplink channel $\bm{h}_{n}\in\mathbb{C}^{M\times 1}$. For simplicity, we
assume an independent Rayleigh small-scale fading model so that the channel
vectors are statistically independent of each other and spatially white, i.e.,
$\bm{h}_{n}\sim\mathcal{CN}_{\\!M}(\bm{0},\bm{I})$.

To account for the random activation of the MTDs, we define the activity
indicator function as

$\alpha_{n}\\!=\\!\begin{cases}1,\leavevmode\nobreak\ \text{if device $n$ is
active}\\\ 0,\leavevmode\nobreak\ \text{otherwise}\\\ \end{cases}.$ (1)

As a result of (1), $\alpha_{n}$ induces an effective channel
$\bm{x}_{n}=\sqrt{\gamma_{n}}\bm{h}_{n}$ for the $n^{\text{th}}$ MTD. Here,
$\gamma_{n}=\alpha_{n}\rho_{n}\beta_{n}$ is the effective power of the channel
of the $n$-th MTD, where $\beta_{n}>0$ is the large-scale fading component
(LSFC) of the $n^{\text{th}}$ device due to path-loss and shadowing, while
$\rho_{n}$ is the device’s uplink transmission power. Hence, the concatenation
of all the effective power levels of the channels of the MTDs in the network
is defined by
$\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{N})^{\top}\in\mathbb{R}^{N}$.
Similarly, let $\bm{X}=(\bm{x}_{1}\,\cdots\,\bm{x}_{N})^{\mbox{\scriptsize
T}}\in\mathbb{C}^{N\times M}$ denote the concatenation of all the effective
channel vectors. Since only a small set of MTDs are active, $\bm{X}$ is row-
sparse. This leads to the uplink signal being a “compressively sensed” signal
coming exclusively from active devices, as will be seen later. As such, the
number of active MTDs is equivalent to the row-support of
$\bm{X}\in\mathbb{C}^{N\times M}$. The index set of rows containing the non-
zero elements is defined as

$\displaystyle\mathcal{K}$
$\displaystyle=\mathrm{supp}(\bm{X})\triangleq\\{i\in[\\![N]\\!]:x_{nm}\neq
0\,\mbox{ for some $m\in[\\![M]\\!]$}\\},$

thus collecting the indices of the active devices, i.e.,
$\mathcal{K}=\\{n\in[\\![N]\\!]:\alpha_{n}=1\\}$. Note that due to the low
likelihood of the MTDs being active simultaneously,
$\lvert\mathcal{K}\rvert<N$ is a practical network property exploited in the
CS-MUD framework, as will be shown next.

As mentioned earlier, $L$ symbols are used as pilot sequences, i.e., used
mainly for metadata processing, such as to estimate the active MTDs. To avoid
signaling overheads, these pilot sequences are normally preallocated by the BS
and stored in the local memory of the MTDs. Let
${\bm{a}}_{n}\in\mathbb{C}^{L\times 1}$ denote the normalized pilot sequence
of the $n^{\text{th}}$ MTD. These can be collected into a matrix
$\bm{A}=({\bm{a}}_{1}\,\cdots\,{\bm{a}}_{N})\in\mathbb{C}^{L\times N}$. In the
end, the first $L$ symbols of the noisy signal received at the BS, exclusively
coming from the $K$ active MTDs, can now be modeled as

$\displaystyle\bm{Y}$
$\displaystyle=\sum_{k\in\mathcal{K}}\sqrt{\rho_{k}\beta_{k}}\alpha_{k}\bm{a}_{k}\bm{h}_{k}^{\mbox{\scriptsize
T}}+\bm{W}$ $\displaystyle=\bm{A}\bm{X}+\bm{W},$ (2)

where $\bm{W}\in\mathbb{C}^{L\times M}$ is an additive noise matrix whose
elements are independent and identically distributed (i.i.d.) Gaussian random
variables with known variance $\sigma^{2}>0$, i.e.,
$w_{lm}\sim\mathcal{CN}(0,\sigma^{2})$. Owing to the massiveness of the MTDs,
the length of the pilot sequences is normally very short compared to the
number of devices, i.e., $L\ll N$, essentially leading to non-orthogonal pilot
sequences that can have noticeable mutual coherence [3, 4]. Therefore, the
task of the BS is to identify the active users by estimating the row-support
$\mathcal{K}$ of $\bm{X}$. However, the condition $L\ll N$ indicates (2) to
follow an underdetermined system111An undetermined system of linear equations
is characterized by having more unknowns than the number of equations. of
equations, thus leading to an infinite number of solutions. To overcome this
challenge, one needs to resort to CS-MUD methodology, such as those that were
mentioned under Section I-A. In the sequel, we propose greedy pursuit approach
for solving the Gaussian log-likelihood function, as we discuss next.

Since there is no correlation in the receive antennas, the columns
$\\{\bm{y}_{m}\in\mathbb{C}^{L}\\}_{m=1}^{M}$ of $\bm{Y}$ are independent of
each other. This means that
$\bm{y}_{m}\sim\mathcal{CN}_{L}(\bm{0},\bm{\Sigma})$ has a positive definite
Hermitian (PDH) $L\times L$ covariance matrix

$\displaystyle\bm{\Sigma}$
$\displaystyle=\bm{A}\bm{\Gamma}\bm{A}^{\mathsf{H}}+\sigma^{2}\bm{I}$
$\displaystyle=\sum_{n=1}^{N}\gamma_{n}{\bm{a}}_{n}{\bm{a}}_{n}^{\mathsf{H}}+\sigma^{2}\bm{I},$
(3)

where $\bm{\Gamma}=\operatorname{diag}(\bm{\gamma})$ is a diagonal matrix
containing the effective powers of the MTDs. Note that
$\mathcal{K}=\mathrm{supp}(\bm{X})=\mathsf{supp}(\bm{\gamma})$, which makes it
possible to estimate the support using (II). As such, _covariance learning_
(CL-) based support recovery algorithms can be constructed by minimizing the
negative log-likelihood function (LLF) of the data $\bm{Y}$. The negative LLF
is defined by (after scaling by $1/M$ and ignoring additive constants)

$\displaystyle\ell(\bm{\gamma})=\operatorname{tr}((\bm{A}\bm{\Gamma}\bm{A}^{\mathsf{H}}+\sigma^{2}\bm{I})^{-1}\hat{\bm{\Sigma}})+\log|\bm{A}\bm{\Gamma}\bm{A}^{\mathsf{H}}+\sigma^{2}\bm{I}|,$
(4)

where $\hat{\bm{\Sigma}}$ is the sample covariance matrix (SCM),

$\hat{\bm{\Sigma}}=\frac{1}{M}\sum_{m=1}^{M}\bm{y}_{m}\bm{y}_{m}^{\mathsf{H}}={M}^{-1}\bm{Y}\bm{Y}^{\mathsf{H}},$

while $\operatorname{tr}(\cdot)$ and $|\cdot|$ denote the matrix trace and
determinant, respectively. Based on the LLF in (4), we next construct a
covariance-based greedy matching algorithm that estimates $\mathcal{K}$ by
solving for $\bm{\gamma}$.

##  III Covariance based matching pursuit

Consider the conditional likelihood of (4) where device powers $\gamma_{n}$
for $n\neq i$ are known, and define

$\bm{\Sigma}_{\backslash i}=\sum_{n\neq
i}\gamma_{n}{\bm{a}}_{n}{\bm{a}}_{n}^{\mathsf{H}}+\sigma^{2}\bm{I}=\bm{\Sigma}-\gamma_{i}{\bm{a}}_{i}{\bm{a}}_{i}^{\mathsf{H}}$
(5)

as the covariance matrix of $\bm{y}_{m}$-s when the contribution from the
$i^{\text{th}}$ device is removed. As $\gamma_{n}$ for $n\neq i$ are assumed
to be given, $\bm{\Sigma}_{\backslash i}$ is a known fixed matrix. It then
follows that the conditional negative LLF for the single unknown
$i^{\text{th}}$ device power is defined as

$\displaystyle\ell_{i}(\gamma)$
$\displaystyle=\operatorname{tr}((\bm{\Sigma}_{\setminus
i}+\gamma{\bm{a}}_{i}{\bm{a}}_{i}^{\mathsf{H}})^{-1}\hat{\bm{\Sigma}})+\log|\bm{\Sigma}_{\setminus
i}+\gamma{\bm{a}}_{i}{\bm{a}}_{i}^{\mathsf{H}}|.$ (6)

This function has a unique optimal value that is given by [32]:

$\gamma_{i}=\max\left(\frac{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}_{\setminus
i}^{-1}\hat{\bm{\Sigma}}\bm{\Sigma}_{\setminus
i}^{-1}{\bm{a}}_{i}-{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}_{\setminus
i}^{-1}{\bm{a}}_{i}}{({\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}_{\setminus
i}^{-1}{\bm{a}}_{i})^{2}},0\right),$ (7)

where $\bm{\Sigma}_{\setminus i}$ is already defined in (5).

This result will be the foundation for developing the CL-MP algorithm. We
follow the generic matching pursuit strategy (see, e.g., [33]), and our
proposed CL-based method is an adaptation of the CL orthogonal matching
pursuit (CL-OMP) algorithm developed in [34]. There are two important
distinctions to [34]. First, we assume that the noise power level $\sigma^{2}$
is known as is common in mMTC applications. Second, unlike in [34], we follow
the MP strategy instead of the OMP as in [34]. This allows us to skip the
“update provisional solution” step of the OMP algorithm and to use the
computationally light update rule for the inverse of the covariance matrix,
thus, speeding up the computations. We denote the $k^{\text{th}}$ update of
$\bm{\gamma}$ and $\mathcal{K}$ by $\bm{\gamma}^{(k)}$ and
$\mathcal{K}^{(k)}$, respectively. Therefore, starting with $k=0$, our CL-MP
algorithm follows the steps below.

Initialization phase. Set $k=0$, and $\bm{\gamma}^{(0)}=\bm{0}_{N\times 1}$,
and $\mathcal{K}^{(0)}=\mathsf{supp}(\bm{\gamma}^{(0)})=\emptyset$ as initial
solutions of device powers and the signal support, respectively. Then
$(\bm{\Sigma}^{(0)})^{-1}=(1/\sigma^{2})\bm{I}$ is the initial inverse
covariance matrix at the start of iterations.

Main Iteration phase consists of the following steps:

1) Sweep: Compute the errors

$\displaystyle\epsilon_{i}$ $\displaystyle=\min_{\gamma\geq
0}\ell_{i}(\gamma)$ (8)

for each $i\in\mathcal{K}^{\complement}=[\\![N]\\!]\setminus\mathcal{K}^{(k)}$
using its unique optimal value

$\gamma_{i}=\max\Big{(}\frac{{\bm{a}}_{i}^{\mathsf{H}}(\bm{\Sigma}^{(k)})^{-1}\hat{\bm{\Sigma}}(\bm{\Sigma}^{(k)})^{-1}{\bm{a}}_{i}-{\bm{a}}_{i}^{\mathsf{H}}(\bm{\Sigma}^{(k)})^{-1}{\bm{a}}_{i}}{({\bm{a}}_{i}^{\mathsf{H}}(\bm{\Sigma}^{(k)})^{-1}{\bm{a}}_{i})^{2}},0\Big{)}.$

2) Update support: Find a minimizer, $i_{k}$ of $\epsilon_{i}$: $\forall
i\in\mathcal{K}^{\complement}$, $\epsilon_{i_{k}}\leq\epsilon_{i}$, and update
the support $\mathcal{K}^{(k+1)}=\mathcal{K}^{(k)}\cup\\{i_{k}\\}$.

3) Update the inverse covariance matrix: Since the new covariance matrix after
support update is

$\bm{\Sigma}^{(k+1)}=\bm{\Sigma}^{(k)}+\gamma_{i_{k}}{\bm{a}}_{i_{k}}{\bm{a}}_{i_{k}}^{\mathsf{H}},$

its inverse covariance matrix $(\bm{\Sigma}^{(k+1)})^{-1}$ can be obtained
using Sherman-Morrison formula222The Sherman-Morrison formula states that
$(\bm{A}+\bm{u}\bm{v}^{\mathsf{H}})^{-1}=\bm{A}^{-1}-\bm{A}^{-1}\bm{u}\bm{v}^{\mathsf{H}}\bm{A}^{-1}/(1+\bm{v}^{\mathsf{H}}\bm{A}^{-1}\bm{u})$
as shown in line 5 of algorithm 1.

4) Stopping rule: If stopping rule is not met, then increment $k$ by 1 and
repeat steps 1)-3). Several stopping rules can be used. For example, one can
stop iterations after the fixed number (say $K_{\max}$) or when the added
index has power level $\gamma_{i_{k}}$ that falls below some threshold level.
In the former case, $K_{\max}$ serves as the maximum number of active devices
expected to be active.

The error term $\epsilon_{i}$ in (8) has the following closed-form expression:

$\epsilon_{i}=c+\log(1+\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i})-\gamma_{i}\frac{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}\hat{\bm{\Sigma}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}{1+\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}},$
(9)

where we have for simplicity of notation written
$\bm{\Sigma}=\bm{\Sigma}^{(k)}$ and where $c$ denotes an irrelevant constant
that is not dependent on $\gamma_{i}$. Thus, w.lo.g., we set $c=0$. Then using
that $\gamma_{i}$ is given by (7) we can write

$\displaystyle
1+\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}$
$\displaystyle=1+\frac{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}\hat{\bm{\Sigma}}\bm{\Sigma}^{-1}{\bm{a}}_{i}-{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}{({\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i})^{2}}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}$
$\displaystyle=1+\frac{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}\hat{\bm{\Sigma}}\bm{\Sigma}^{-1}{\bm{a}}_{i}-{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}$
$\displaystyle=\frac{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}\hat{\bm{\Sigma}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}\mbox{
when $\gamma_{i}>0$}.$

Substituting this into the denominator of the last term in (9), we obtain

$\epsilon_{i}=\log(1+\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i})-\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}.$
(10)

If $\gamma_{i}=0$, then $\epsilon_{i}=0$. This explains line 3 in algorithm 1.

Input : $\hat{\bm{\Sigma}}=M^{-1}\bm{Y}\bm{Y}^{\mathsf{H}}$, $\bm{A}$,
$\sigma^{2}$

Initialize : $\bm{\Sigma}^{-1}=(1/\sigma^{2})\bm{I}$, $\mathcal{K}=\emptyset$

1

2for  _$k=1,\ldots,N$_ do

3

4
$\gamma_{i}\leftarrow\max\Big{(}\frac{{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}\hat{\bm{\Sigma}}\bm{\Sigma}^{-1}{\bm{a}}_{i}-{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}}{({\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i})^{2}},0\Big{)}$,
$\forall i\in\mathcal{K}^{\complement}$

5
$\epsilon_{i}\leftarrow\log(1+\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i})-\gamma_{i}{\bm{a}}_{i}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i}$,
$\forall i\in\mathcal{K}^{\complement}$

6 $\mathcal{K}\leftarrow\mathcal{K}\cup\\{i_{k}\\}$ with
$i_{k}\leftarrow\arg\min_{i\in\mathcal{K}^{\complement}}\epsilon_{i}$

7
$\bm{\Sigma}^{-1}\leftarrow\bm{\Sigma}^{-1}-\dfrac{\gamma_{i_{k}}\bm{\Sigma}^{-1}{\bm{a}}_{i_{k}}{\bm{a}}_{i_{k}}^{\mathsf{H}}\bm{\Sigma}^{-1}}{1+\gamma_{i_{k}}{\bm{a}}_{i_{k}}^{\mathsf{H}}\bm{\Sigma}^{-1}{\bm{a}}_{i_{k}}}$

8 if  _stopping rule is met (see text)_ then

break

end if

end for

Output : A set $\mathcal{K}$, containing $I$ indices, where $I$ is the number
of iterations completed.

Algorithm 1 Covariance Learning Matching Pursuit (CL-MP) algorithm

###  III-A Computational complexity comparisons

In Table I, we present the computational complexity analysis of CL-MP and some
benchmark algorithms. In the complexity calculations, it is assumed that the
number of antennas ($M$) is smaller than the number of users ($N$) and the
pilot length ($L$) is smaller than $N$, which are typical scenarios in mMTC.
In SOMP, the index $k$ denotes the iteration index (i.e., the number of atoms
selected). The dominant complexity of CL-MP Algorithm 1 is the matrix-vector
multiplications in lines 2, 3, and 5, whose complexity is
$\mathcal{O}(L^{2})$. In lines 2 and 3, one performs computations for
(essentially) all $N$ devices and hence the total complexity of these steps is
$\mathcal{O}(NL^{2})$. In line 3, one searches the smallest value in a set of
$N$ scalars with linear complexity $\mathcal{O}(N)$. In line 5, the inverse
$\bm{\Sigma}^{-1}$ is computed using a rank-one update, which further improves
the computational efficiency as inverting an $L\times L$ matrix is not needed
(which is of complexity $\mathcal{O}(L^{3})$ instead of $\mathcal{O}(L^{2})$.
Hence, the full complexity of Algorithm 1 per iteration is
$\mathcal{O}(NL^{2})$. Thus the complexity of CL-MP is linear in $N$ and
quadratic in $L$. This is a desired feature for low-latency mMTC scenarios
where the number of MTDs $N$ is large while the pilot length $L$ is often
small. Note that SBL and VAMP grow quadratic in $N$.

It can also be noted that the computational complexity of CL-MP is the same as
that of CWO per iteration. However, CL-MP is a greedy approach, while CWO is
iterative, whose best performance is obtained after the algorithm’s
convergence. Unlike iterative approaches, greedy methods do not suffer from
convergence issues, which often happens when the nominal assumptions do not
hold (e.g., when Gaussianity of data does not hold or in the presence of
outliers).

It should also be stressed that the final complexity of all algorithms is the
complexity per iteration times the number of iterations $I$ needed for
convergence. For example, for SBL the full complexity is
$\mathcal{O}(IN^{2}L)$. It is well known that the Expectation Maximization
(EM) algorithm utilized by SBL suffers from slow convergence, and thus $I$ can
often be very large. Such problems are not faced by greedy pursuit approaches
which aim to find the minimum of the LLF greedily instead of using an
iterative optimization algorithm.

Table I: Computational complexity of CS-MUD algorithms. Algorithm | Complexity (per iteration)
---|---
CL-MP | $\mathcal{O}(NL^{2})$
CWO [28] | $\mathcal{O}(NL^{2})$
VAMP [35] | $\mathcal{O}(N^{2}M)$
SBL [36] | $\mathcal{O}(N^{2}L)$
SOMP [37] | $O(L\cdot\max(NM,k^{2})+k^{3})$

##  IV Empirical comparisons

This section validates the performance of the proposed CL-MP algorithm under a
variety of settings. As performance metrics, we used missed detection and
exact recovery. Missed detection (MD) occurs when a device is active but
detected inactive. Therefore, the (empirical) probability of misdetection is
the Monte Carlo (MC) average ratio of the number of such devices relative to
the total number of active MTDs and is defined as

$\text{P}_{md}=\frac{1}{Q}\sum_{q=1}^{Q}\frac{|\mathcal{K}\setminus\mathcal{K}^{(q)}|}{|\mathcal{K}|},$
(11)

where $Q$ is the number of MC trials, and $\mathcal{K}^{(q)}$ is the estimate
of the support $\mathcal{K}$ of the $q^{\text{th}}$ MC trial. On the other
hand, exact recovery (ER) of devices occurs when active devices are correctly
detected, and it is defined by

$\text{P}_{er}=\frac{1}{Q}\sum_{q=1}^{Q}\mathbbm{1}(\mathcal{K}^{(q)}=\mathcal{K}).$
(12)

We note that the indices of active devices $\mathcal{K}$ (with
$|\mathcal{K}|=K$) for each $q^{\text{th}}$ MC trial are also randomly
selected from set $[\\![N]\\!]$ without replacement.

###  IV-A The effect of $M$, $L$ and $K$

$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$\text{P}_{md}$CWO
$K=10$$K=20$$K=30$$K=40$$K=50$CL-MP $K=10$$K=20$$K=30$$K=40$$K=50$

$10$$15$$20$$25$$30$$35$$40$$45$$50$$55$$60$$0$$0.2$$0.4$$0.6$$0.8$$1$$M$
(number of antennas)$\text{P}_{er}$

Figure 2: Probability of misdetection (top) and exact recovery (bottom) vs $M$
for different levels of active devices ($K$); $L=64$, $N=1000$, and LSFC-s
$\beta_{n}\sim Unif(-15,0)$ on dB scale.

We first investigate the effects of $L$, $M$, and $K$ on the AUD performance.
We let the LSFCs $\beta_{n}$, $\forall n$ to follow uniform distribution in dB
scale, i.e., $10\log_{10}(\beta_{n})$ is uniformly distributed in range
$[10\log_{10}(\beta_{\text{min}}),10\log_{10}(\beta_{\text{max}})]$. We
compare the proposed CL-MP in algorithm 1 against the CWO algorithm [28,
Algorithm 1]. We set the maximum number of iterations of the CWO algorithm to
15. To make the comparisons fair, we assume that the number of active devices
$K$ is known to both algorithms. Hence the CL-MP algorithm is terminated after
$K$ iterations and directly returns the estimated support $\hat{\mathcal{K}}$
of size $|\hat{\mathcal{K}}|=K$. In the case of the CWO algorithm,
thresholding is needed. Here, the active users are identified by picking the
positions of the $K$ largest entries of the final estimate
$\hat{\bm{\gamma}}$. The number of MC runs is $2000$.

$10$$20$$30$$40$$50$$0.00$$2.00$$4.00$$6.00$$8.00$$\cdot 10^{-2}$$K$ (number
of active devices)Running time [s]CL-MPCWO

$30$$40$$50$$60$$70$$80$$0.00$$0.02$$0.04$$0.06$$0.08$$0.10$$L$ (pilot
length)Running time [s]CWO$30$$40$$50$$60$$70$$80$$L$ (pilot length)CL-
MP$K=10$$K=20$$K=30$$K=40$$K=50$

Figure 3: Top:Average running time in seconds for different $K$ when $L=64$
and $N=1000$. Bottom: Average running time for different $L$ and $K$ when
$M=40$ for CWO (left panel) and CL-MP (right panel) methods.

Figure 2 displays the probability of misdetection and exact recovery as a
function of number of antennas ($M$) and for different numbers of active
devices ($K$). Bernoulli pilot sequence of length $L=64$ are used and the
number of devices is $N=1000$. The pilots are normalized to unit energy per
symbol, i.e., $\|{\bm{a}}_{n}\|_{2}^{2}=L$. LSFC coefficients follow a uniform
distribution with a range $[-15,0]$ dB. Both methods show similar performance,
but the CL-MP has a slight advantage in all levels of $K$ and $M$. In addition
to this, the running times (on MacBook Pro laptop 2.3GH Intel i9) are shown in
the top panel of Figure 3 demonstrate that it achieves high-level performance
in significantly lower running time than the CWO algorithm. The differences
between the methods are more pronounced when the number of active devices is
small. Since CWO and CL-MP algorithms do not scale with $M$, the CPU bars are
averages over all instances of $M$.

Figure 4 displays the probability of misdetection and exact recovery as a
function of pilot sequence length ($L$) and for different levels of active
devices ($K$) when the number of antennas is fixed ($M=40$). Again, we can
note that the proposed CL-MP performed better for all levels of $K$ and $L$.
The running times (on MacBook Pro laptop 2.3GH Intel i9) for both methods are
displayed at the bottom panel of Figure 3. As noted, CL-MP scales better with
$L$ than CWO, and the advantage is, as expected, more significant when $K$ is
small.

$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$\text{P}_{md}$CWO
$K=10$$K=20$$K=30$$K=40$$K=50$CL-MP $K=10$$K=20$$K=30$$K=40$$K=50$

$30$$35$$40$$45$$50$$55$$60$$65$$70$$75$$80$$0$$0.2$$0.4$$0.6$$0.8$$1$$L$
(pilot length)$\text{P}_{er}$

Figure 4: Probability of misdetection (top) and exact recovery (bottom) vs $L$
for different levels of active devices ($K$); $N=1000$, $M=40$ and LSFC-s
$\beta_{n}\sim Unif(-15,0)$ on dB scale.

###  IV-B MIMO uplink network simulation

In this subsection, we present the results of the performance of our proposed
CL-MP in the mMTC network. For these simulations, we consider a single-cell
massive MIMO uplink network, where a BS serves $1000$ MTDs randomly placed
within the cell between a radius of $25$ m and $250$ m. We adopt a log-
distance path loss model where LSFC is $\beta_{n}=-130-37.6\log_{10}d_{n}$,
where $d_{n}$ is the distance between the $n^{\text{th}}$ MTD and the BS. For
resource fairness, we assume uplink power control following channel inversion,
such that each MTD transmits with power[8]

$\rho_{n}=\frac{p^{\mathrm{max}}\beta_{\mathrm{min}}}{\beta_{n}},$ (13)

where $\beta_{\mathrm{min}}$ is the large-scale effect at the edge of the cell
and $p^{\mathrm{max}}$ is the maximum allowable power for each MTD.
Consequently, the received signal-to-noise ratio (SNR) for each MTD is defined
by

$\text{SNR}_{n}=\frac{p^{\mathrm{max}}\beta_{\mathrm{min}}}{\sigma^{2}}.$ (14)

Unless otherwise stated, the simulation parameters used in this paper are
summarized in Table II.

For the performance benchmark, we compare CL-MP with Vector AMP (VAMP) [35],
SBL [36], and SOMP [37]. We acknowledge that all these benchmark solutions
first estimate $\bm{X}$, then map it to $\hat{\bm{\alpha}}$. To ensure
fairness in our comparison, we input the number of active MTDs $K$ into each
algorithm. This ensures that $K$ indices corresponding to maximum row norms
are selected from matrix $\bm{X}$ after the convergence of each algorithm.
Therefore, the set of indices of the selected $K$ rows will serve as an
estimate $\hat{\mathcal{K}}$ of $\mathcal{K}$ at the end of each algorithm.
Each algorithm will have two stopping criteria, i.e. if the maximum number of
iterations $T_{max}$ has been reached or the relative change between two
consecutive estimates $\Delta$ becomes insignificant. For random sensing
matrices, we use the Gaussian pilot sequences, which we denote by G-MAT, and
the Bernoulli pilot sequences, which we denote by B-MAT. For structured pilot
sequences, we generate pilot sequences following the restrictively combined
fixed orthogonal block matrix from [4], which we denoted by RC-FOBM as well as
restively combined random orthogonal block matrix, which we denote by RC-ROBM.

Table II: Simulation parameters Parameter |  Value
---|---
Cell radius |  $250$ m
Number of MTDs $(N)$ |  $1000$
Bandwidth |  $20$ MHz
Noise power ($\sigma^{2}$) |  $2\\!\times\\!10^{-13}$ W
Coherence interval $(T)$ |  $300$
Length of the pilot sequences $(L)$ |  $64$
Number of BS antennas $(M)$ |  $32$
Activation probability $(\epsilon)$ |  $0.01$
$\text{Maximum MTD power}(p^{\text{max}})$ |  $0.1$ W
Average SNR |  $10$ dB
Error tolerance $(\Delta)$ |  $10^{-6}$
Monte Carlo $(\text{MC})$ |  $2000$
Maximum Iterations $(T_{max})$ |  $150$

Fig. 5 presents the results in terms of the $\text{P}_{md}$ as a function of
the SNR, for different algorithms. The general trend shows that increasing SNR
improves the $\text{P}_{md}$ of all the algorithms. This is because as the SNR
increases, the desired signal component, which in our case is represented by
$\bm{A}\bm{X}$, has higher power relative to the noise. As a result, it
becomes easier to recover and identify the active users. The results show that
both SBL and CL-MP have comparable performance as the SNR increases and their
performance is superior to that of the SOMP and VAMP. The comparable
performance between SBL and CL-MP is due to the equivalence of their objective
functions as mentioned in Section I-B. However, the superior performance of
CL-MP validates the proposed solutions. Since CL-MP has a similar performance
as SBL, we next evaluate its performance for different sensing matrices.

$-16$$-15$$-14$$-13$$-12$$-11$$-10$$-9$$-8$$-7$$-6$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$SNR$\text{P}_{md}$VAMPSOMPSBLCL-
MP Figure 5: $\text{P}_{md}$ as a function of SNR for different algorithms
$-20$$-15$$-10$$-5$$0$$5$$10$$15$$20$$10^{-2}$$10^{-1}$$10^{0}$SNR$\text{P}_{md}$RC-
FOBMRC-ROBMG-MATB-MAT Figure 6: $\text{P}_{md}$ of CL-MP as a function of the
SNR.

Fig. 6 presents the performance in terms of the probability of miss detection
($\text{P}_{md}$) as a function of the SNR and the structure of the pilot
sequences used. Similar to the results in Fig. 5, the performance generally
improves as the SNR increases. Despite this, the proposed solution CL-MP
performs differently for varying sensing matrices. This shows the non-
triviality of choosing pilot sequences in MTC [4, 3]. For instance, the
performance is superior when using Bernoulli pilot sequences and least when
using the Gaussian pilot sequences. This is because pilot sequences such as
Gaussian tend to increase the collision probability of the pilot sequences,
thus, degrading the performance [4, 3, 8]. On the other hand, Bernoulli, RC-
FOBM, and RC-ROBM have reduced collision probabilities and thus have better
performance than the Gaussian pilot sequences. Despite this, the proposed CL-
MP performs reasonably well under the different set of pilot sequences.
Another interesting result is that the CL-MP performs reasonably well at very
low SNR values, reaching $P_{md}=0.01$ for BMAT at SNR of $0$ dB.

$0$$5$$10$$15$$20$$25$$30$$35$$40$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$M$
(number of antennas)$\text{P}_{md}$VAMPSOMPSBLCL-MP Figure 7: $\text{P}_{md}$
as a function of the number of antennas

Fig. 7 shows the performance in terms of $\text{P}_{md}$ as a function of the
number of antennas for different algorithms. The general trend shows that the
performance of all the algorithms improves with an increase in the number of
antennas ($M$). This is mainly because an increase in $M$ imposes more
structural sparsity on the matrix $\bm{X}$, thus increasing the accuracy of
the estimation. This reveals the advantage of having multiple BS antennas,
which leads to a multiple measurement vector (MMV) problem in compressed
sensing[8]. However, it can be observed that both CL-MP and SBL have
equivalent and superior performance compared to VAMP and SOMP, thus,
confirming the superiority of our proposed solution. Though both SBL and CL-MP
are comparable in terms of performance in $\text{P}_{md}$, it is worth
emphasizing that CL-MP has lower computational complexity than SBL, as shown
in Table I, making it a more practical solution in MTC. To substantiate this,
we compare the run-times for varying $M$.

$0$$5$$10$$15$$20$$25$$30$$35$$40$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$M$
(number of antennas)Running time [s]VAMPSOMPSBLCL-MP Figure 8: Run-time as a
function of the number of antennas

Figure 8 presents the results of the run-times of each algorithm for a
different number of antennas of the BS. Effectively, this reveals how each
algorithm scales up. From the results, it can be observed that the proposed
CL-MP has superior scale-up properties as compared to all the benchmarks. On
the other hand, it is worth noting that, even though SBL has similar
performance as the CL-MP in terms of $\text{P}_{md}$, it takes a longer time
to run and is thus not favorable in mMTC.

$10$$15$$20$$25$$30$$35$$40$$45$$50$$55$$60$$65$$70$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$L$
(pilot length)$\text{P}_{md}$VAMPSOMPSBLCL-MP Figure 9: $\text{P}_{md}$ as a
function of the pilot lengths for different algorithms.

Fig. 9 shows the performance as a function of the pilot lengths ($L$) for
different algorithms. The trend reveals that the $\text{P}_{md}$ reduces as
the pilot length increases. This is mainly due to better pilot structure as
the length increases. Essentially, the signal detection phase is done over a
longer period, thus improving the accuracy of the channel estimation for SOMP,
VAMP and SBL, as well as the index estimation of the CL-MP. Similar to the
results in Fig 7 and Fig. 5, CL-MP and SBL perform better than all the other
benchmarks. Consistent with the results of Fig. 8, CL-MP has the added
advantage of being more scalable.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$10^{-2}$$10^{-1}$$10^{0}$$N$ (number of
MTDs)Running time [s]$L=16$$L=32$$L=48$$L=56$$L=72$$\times 10^{3}$ Figure 10:
Runtime of CL-MP as a function of the number $N$ of the MTDs in the network

To conclude the results, Fig. 10 presents the runtime of CL-MP as a function
of the number of MTDs. In general, the runtime increases with the number of
MTDs, which is expected because the size of $\bm{\gamma}$ increases. Despite
this, the CL-MP has lower runtimes for a large number of users. For instance,
it takes approximately $1$ second for CL-MP to recover and perform AUD from
the signal received from $N=6000$ users when the pilot length is $L=72$. This
translates to better runtime than SBL, which is shown to require more than $5$
seconds to process signals from $N=1000$ users in Fig. 8. To further support
the results of Fig. 10, in Fig. 11, we present the results of $\text{P}_{md}$
for varying sizes of the networks. In general, the proposed CL-MP algorithm is
effective in solving AUD problems, with very low $\text{P}_{md}$ even in very
large MTC networks, rendering our proposal a practical solution for the
rapidly growing MTC networks.

$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$0$$0.2$$0.4$$0.6$$0.8$$1$$N$ (number of
MTDs)$\text{P}_{md}$$L=16$$L=32$$L=48$$L=56$$L=72$$\times 10^{3}$ Figure 11:
$\text{P}_{md}$ as a function of the number $N$ of the MTDs in the network.

##  V Conclusions and Future Directions

This paper introduces a novel framework for estimating the active user indices
in MTC. The proposed technique uses the covariance of the received signal to
formulate an optimization problem for the sparse effective power levels in MTC
and solves it using a greedy approach. The results demonstrate that the
proposed framework offers highly accurate device activity detection with very
low runtimes. However, while efficient, the proposed solution does not
estimate the activation probability. Therefore, the results of this work can
be extended to incorporate activation probability estimation as future
research direction. Moreover, considering that the IoT framework will
eventually be served by Low Earth Orbit (LEO) satellites, one potential
extension of this work could be the integration of CL-MP in satellite
receivers, where the channel has a line of sight.

## References

  * [1] H. Shariatmadari, R. Ratasuk, S. Iraji, A. Laya, T. Taleb, R. Jäntti, and A. Ghosh, “Machine-type communications: Current status and future perspectives toward 5G systems,” _IEEE Communications Magazine_ , vol. 53, no. 9, pp. 10–17, 2015.
  * [2] E. Dutkiewicz, X. Costa-Perez, I. Z. Kovacs, and M. Mueck, “Massive machine-type communications,” _IEEE Network_ , vol. 31, no. 6, pp. 6–7, 2017\.
  * [3] L. Marata, O. L. A. López, H. Djelouat, M. Leinonen, H. Alves, and M. Juntti, “Joint coherent and non-coherent detection and decoding techniques for heterogeneous networks,” _IEEE Transactions on Wireless Communications_ , vol. 22, no. 3, pp. 1730–1744, 2022.
  * [4] L. Marata, O. L. A. López, A. Hauptmann, H. Djelouat, and H. Alves, “Joint activity detection and channel estimation for clustered massive machine type communications,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2023.
  * [5] J. Gao, Y. Wu, T. Li, and W. Zhang, “Energy Efficiency of MIMO Massive Unsourced Random Access with Finite Blocklength,” _IEEE Wireless Communications Letters_ , 2023.
  * [6] Y. Liu, S. Zhang, X. Mu, Z. Ding, R. Schober, N. Al-Dhahir, E. Hossain, and X. Shen, “Evolution of NOMA toward next generation multiple access (NGMA) for 6G,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 4, pp. 1037–1071, 2022.
  * [7] J. Choi, J. Ding, N.-P. Le, and Z. Ding, “Grant-free random access in machine-type communication: Approaches and challenges,” _IEEE Wireless Communications_ , vol. 29, no. 1, pp. 151–158, 2021.
  * [8] K. Senel and E. G. Larsson, “Grant-free massive MTC-enabled massive MIMO: A compressive sensing approach,” _IEEE Transactions on Communications_ , vol. 66, no. 12, pp. 6164–6175, 2018.
  * [9] L. Liu, E. G. Larsson, W. Yu, P. Popovski, C. Stefanovic, and E. De Carvalho, “Sparse signal processing for grant-free massive connectivity: A future paradigm for random access protocols in the Internet of Things,” _IEEE Signal Processing Magazine_ , vol. 35, no. 5, pp. 88–99, 2018.
  * [10] B. Li, J. Zheng, and Y. Gao, “Compressed sensing based multiuser detection of grant-free NOMA with dynamic user activity,” _IEEE Communications Letters_ , vol. 26, no. 1, pp. 143–147, 2021.
  * [11] J. Liu, G. Wu, X. Zhang, S. Fang, and S. Li, “Modeling, analysis, and optimization of grant-free NOMA in massive MTC via stochastic geometry,” _IEEE Internet of Things Journal_ , vol. 8, no. 6, pp. 4389–4402, 2020.
  * [12] X. Zhang, P. Fan, J. Liu, and L. Hao, “Bayesian learning-based multiuser detection for grant-free NOMA systems,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 8, pp. 6317–6328, 2022.
  * [13] R. B. Di Renna, C. Bockelmann, R. C. de Lamare, and A. Dekorsy, “Detection techniques for massive machine-type communications: Challenges and solutions,” _IEEE Access_ , vol. 8, pp. 180 928–180 954, 2020.
  * [14] Y. Polyanskiy, “A perspective on massive random-access,” in _IEEE International Symposium on Information Theory (ISIT)_. IEEE, 2017, pp. 2523–2527.
  * [15] X. Shao, X. Chen, D. W. K. Ng, C. Zhong, and Z. Zhang, “Cooperative activity detection: Sourced and unsourced massive random access paradigms,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 6578–6593, 2020.
  * [16] A. E. Kalør, G. Duris, S. Coleri, S. Parkvall, W. Yu, A. Mueller, and P. Popovski, “Wireless 6G Connectivity for Massive Number of Devices and Critical Services,” _arXiv preprint arXiv:2401.01127_ , 2024.
  * [17] J. Ni and J. Zheng, “Non-coherent grant-free NOMA through pilot-and channel block-index modulation,” _IEEE Wireless Communications Letters_ , vol. 10, no. 4, pp. 705–709, 2020.
  * [18] O. L. López, N. H. Mahmood, M. Shehab, H. Alves, O. M. Rosabal, L. Marata, and M. Latva-Aho, “Statistical tools and methodologies for ultrareliable low-latency communication—a tutorial,” _Proceedings of the IEEE_ , vol. 111, no. 11, pp. 1502–1543, 2023.
  * [19] Z. Zhang, Q. Guo, Y. Li, M. Jin, and C. Huang, “Variational Bayesian inference clustering based joint user activity and data detection for grant-free random access in mMTC,” _IEEE Internet of Things Journal_ , 2023\.
  * [20] X. Zhang, F. Labeau, L. Hao, and J. Liu, “Joint active user detection and channel estimation via Bayesian learning approaches in MTC communications,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 6, pp. 6222–6226, 2021.
  * [21] J. Wang, J. Yi, R. Han, L. Bai, and J. Choi, “Variational Bayesian inference for channel estimation and user activity detection in C-RAN,” _IEEE Wireless Communications Letters_ , vol. 9, no. 7, pp. 953–956, 2020.
  * [22] H. Djelouat, M. Leinonen, L. Ribeiro, and M. Juntti, “Joint user identification and channel estimation via exploiting spatial channel covariance in mMTC,” _IEEE Wireless Communications Letters_ , vol. 10, no. 4, pp. 887–891, 2021.
  * [23] G. Sun, M. Cao, W. Wang, W. Xu, and C. Studer, “Joint active user detection, channel estimation, and data detection for massive grant-free transmission in cell-free systems,” in _IEEE 24th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_. IEEE, 2023, pp. 406–410.
  * [24] L. Zhu, K.-H. Liu, L. Wan, and L. Sun, “Active user detection and channel estimation via fast ADMM,” in _IEEE Wireless Communications and Networking Conference (WCNC)_. IEEE, 2023, pp. 1–6.
  * [25] Z. Chen, F. Sohrabi, Y.-F. Liu, and W. Yu, “Covariance based joint activity and data detection for massive random access with massive MIMO,” in _IEEE International Conference on Communications (ICC)_. IEEE, 2019, pp. 1–6.
  * [26] Y.-F. Liu, W. Yu, Z. Wang, Z. Chen, and F. Sohrabi, “Grant-free random access via covariance-based approach,” _Next Generation Multiple Access_ , pp. 391–414, 2024.
  * [27] W. Liu, Y. Cui, F. Yang, L. Ding, and J. Sun, “MLE-based device activity detection under Rician fading for massive grant-free access with perfect and imperfect synchronization,” _IEEE Transactions on Wireless Communications_ , 2024.
  * [28] A. Fengler, S. Haghighatshoar, P. Jung, and G. Caire, “Non-Bayesian activity detection, large-scale fading coefficient estimation, and unsourced random access with a massive MIMO receiver,” _IEEE Transactions on Information Theory_ , vol. 67, no. 5, pp. 2925–2951, 2021.
  * [29] J. You, W. Wang, S. Liang, W. Han, and B. Bai, “An efficient two-stage SPARC decoder for massive MIMO unsourced random access,” _IEEE Transactions on Wireless Communications_ , 2023.
  * [30] N. Y. Yu and W. Yu, “Joint activity and data detection for massive grant-free access using deterministic non-orthogonal signatures,” _IEEE Transactions on Wireless Communications_ , 2024.
  * [31] L. Marata, “Advanced signal processing techniques for machine type communications,” 2024.
  * [32] T. Yardibi, J. Li, P. Stoica, M. Xue, and A. B. Baggeroer, “Source localization and sensing: A nonparametric iterative adaptive approach based on weighted least squares,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 46, no. 1, pp. 425–443, 2010.
  * [33] M. Elad, _Sparse and redundant representations_. New York: Springer, 2010.
  * [34] E. Ollila, “Sparse signal recovery and source localization via covariance learning,” _arXiv preprint arXiv:2401.13975_ , 2024.
  * [35] L. Liu and W. Yu, “Massive connectivity with massive MIMO—Part I: Device activity detection and channel estimation,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 11, pp. 2933–2946, 2018.
  * [36] D. P. Wipf and B. D. Rao, “Sparse bayesian learning for basis selection,” _IEEE Transactions on Signal processing_ , vol. 52, no. 8, pp. 2153–2164, 2004.
  * [37] J. A. Tropp, “Algorithms for simultaneous sparse approximation. Part II: Convex relaxation,” _Signal Processing_ , vol. 86, no. 3, pp. 589–602, 2006\.


# Fast, low-loss all-optical phase modulation in warm rubidium vapour

Will Davis  Centre for Photonics and Photonic Materials, Department of
Physics, University of Bath, Bath, BA2 7AY, UK  Cameron McGarry  Centre for
Photonics and Photonic Materials, Department of Physics, University of Bath,
Bath, BA2 7AY, UK  Tabijah Wasawo  Centre for Photonics and Photonic
Materials, Department of Physics, University of Bath, Bath, BA2 7AY, UK  Peter
J. Mosley  Centre for Photonics and Photonic Materials, Department of Physics,
University of Bath, Bath, BA2 7AY, UK  Josh Nunn  Centre for Photonics and
Photonic Materials, Department of Physics, University of Bath, Bath, BA2 7AY,
UK  ORCA Computing Ltd. 30 Eastbourne Terrace, London W2 6LA UK

###### Abstract

High-speed switching with low loss would be a versatile tool for photonic
quantum technologies, with applications in state generation, multiplexing, and
the implementation of quantum gates. Phase modulation is one method of
achieving this switching, but existing optical phase modulators either achieve
high bandwidth or low loss, but not both. We demonstrate fast
($100\text{\,}\mathrm{MHz}$ bandwidth), low-loss ($74(2)\text{\,}\mathrm{\char
37\relax}$ transmission) phase shifting ($\Delta\phi=($0.90(5)\text{\,}$)\pi$)
in a signal field, induced by a control field, and mediated by the two-photon
$5S_{1/2}\rightarrow{}5P_{3/2}\rightarrow{}5D_{5/2}$ transition in 87Rb
vapour. We discuss routes to enhance both performance and scalability for
application to a range of quantum and classical technologies.

##  I Introduction

Photonics has revolutionised telecommunications since the development of fibre
optics in the 1980s. Photonic data buses are supplanting electronics in high
performance computing [1], and more recently photonic platforms for machine
learning are emerging [2]. Looking forwards, photonics can provide a platform
for communication with enhanced security by quantum key distribution [3] and
support the transfer of quantum information between nodes [4]. All of these
applications require high speed switching which can be achieved by phase
modulation of an optical signal, and fibre-integrated electro-optical
modulators are commercially mature. Nevertheless the insertion losses of these
devices (typically $3\text{\,}\mathrm{dB}$) add a practical overhead:
mitigating these losses requires increased input power, intermediate
amplifiers, and waste heat management [5]. Further, increasing demand on
switching speeds could lead to the obsolescence of existing semiconductor-
based telecommunications devices [6]. More efficient technologies for optical
modulation are thus desirable across a range of application areas.

The primary motivation for our work in this paper is in the area of photonic
quantum computing [7], where engineered non-classical optical states are used
to solve computational problems that are intractable with classical (i.e. non-
quantum) resources. Photonic quantum computing is appealing for a number of
reasons, including room-temperature operation of all or many components, high
clock-rates, high connectivity, insensitivity to stray fields and modular
construction. But a key technical challenge remains: the requirement to switch
and dynamically re-route photons with high speed and extremely low loss. This
is an essential stage in a variety of processes for photonic quantum
computing, such as implementing: loop memories [8, 9], synchronisation [10] or
multiplexing [11, 12, 13] of single photon sources and demultiplexing for
graph state generation [14]. Amplification destroys quantum coherence and so
cannot be used to mitigate losses in a quantum system. The lifetime of photons
in a waveguide is limited, and so high bandwidth is required for scalability.
For these reasons, quantum systems have extremely stringent tolerances for
speed and loss [15], which motivates an exploration of alternative platforms
that could ultimately deliver better performance than current fibre-integrated
electro-optic modulators.

Figure 1: (a) Ladder scheme used for phase modulation of weak signal field
(red) by presence of strong control field (blue) detuned from
$5S_{1/2}\rightarrow{}5P_{3/2}$ and $5P_{3/2}\rightarrow{}5D_{5/2}$
transitions by $\Delta_{s}$ and $\Delta_{c}$ respectively. (b) Transmission
with control off (black, solid) and on (orange, dotted) and phase shift (red,
dashed) through 87Rb vapour for various $\Delta_{s}$, as determined by the
theoretical model described in the main text. Control detuning is fixed at
$\Delta_{c}=$-1.6\text{\,}\mathrm{GHz}$$. We identify the experimental region
of interest (shaded), where there is low-loss for control on and off, as well
as a high phase shift. Note also the two-photon absorption feature at
$\Delta_{s}=$-6\text{\,}\mathrm{GHz}$$, which is away from
$\Delta_{s}=-\Delta_{c}$ due to the a.c. Stark-shift.

In this letter, we describe and demonstrate efficient, all-optical phase
modulation of light. Ideal implementation of phase modulation (a phase shift
of $\pi$ radians with no loss) would be equivalent to implementation of an
optical switch, by embedding the phase modulator in one arm of a Mach-Zehnder
interferometer. Our scheme, depicted in Fig. 1 (a), takes advantage of the
$5S_{1/2}\rightarrow{}5P_{3/2}\rightarrow{}5D_{5/2}$ two-photon transition in
87Rb. This is the same transition used in electromagnetically-induced
transparency (EIT) optical control [16, 17] and memory [18] schemes. A weak
signal field, detuned from resonance with the $5S_{1/2}\rightarrow{}5P_{3/2}$
transition by frequency $\Delta_{s}$, counter-propagates through a 87Rb vapour
cell with a high-intensity, pulsed control field. This control field is
detuned from resonance with the $5P_{3/2}\rightarrow{}5D_{5/2}$ transition by
frequency $\Delta_{c}$, and its presence induces a change in the
susceptibility of the 87Rb, as experienced by the signal. Hence, the control
field modulates the phase of the signal field as it traverses the vapour cell.
We will first outline a simple theoretical model of this phase modulation,
before describing the experimental results for a continuous wave (CW) and
pulsed signal field.

##  II Theory of optical phase modulation in atoms

In general, light with wavenumber $k$ travelling through a medium of length
$L$ will accumulate a phase

$\phi=kL\real(\sqrt{1+\chi}),$ (1)

where $\chi$ is the susceptibility of the medium. The transmission of the
medium is given by $T=\exp(-\alpha L)$, where $\alpha$ is the extinction
coefficient,

$\alpha=-2k\imaginary(\sqrt{1+\chi}).$ (2)

The difference in phase between two possible paths is the phase shift
($\Delta\phi$). The phase shift and loss experienced by a field, $E$, can be
extracted by interference with the original field. We perform such
interference using a Mach-Zehnder interferometer, so at the two output ports
we have fields $E_{1,2}$, whose amplitudes are

$\begin{split}E_{1}&=\sqrt{T}E\cos(\Delta\phi/2),\\\
E_{2}&=\sqrt{T}E\sin(\Delta\phi/2).\end{split}$ (3)

We can then extract $\Delta\phi$ and $T$ from measurement of the two fields’
intensities.

The phase shift through the 87Rb vapour can be understood as a function of
$\chi$, which in turn depends on the properties of both the control field and
the vapour (specifically its density and temperature). We model $\chi$ by
numerical solution of the Maxwell-Bloch equations [19, 20, 21, 22, 23] (see
supplementary material), allowing us to identify suitable signal and control
detunings where we can achieve both high phase shift and low loss. We consider
a simplified model in which the system is in a steady state and the field
across the medium is uniform – however, both of these assumptions may lead to
disparity with experimental results where we have a pulsed control field and
observe a broadening of the Stark shift across the mode due to the radial
variation of the fields.

The phase shift and transmission for a representative experiment is shown in
Fig. 1 (b), and can be used to identify a region of low loss and high phase
shift (grey shading,
$\Delta_{s}=$-4.8\text{\,}\mathrm{GHz}-2\text{\,}\mathrm{GHz}$$). In this
region there is a flat, non-zero phase shift of $\Delta\phi=-15\pi$, and since
it far detuned from both the single- and two-photon absorption features
($\Delta_{c}=$1.6\text{\,}\mathrm{GHz}$$) it also exhibits low loss. For this
calculation the CW control power is $11\text{\,}\mathrm{W}$, which is
comparable to the peak power of the $47\text{\,}\mathrm{nJ}$,
$4\text{\,}\mathrm{ns}$ FWHM duration control pulses used in Section III. The
signal is also taken to be CW with negligible power. Both beams have a
$125\text{\,}\mathrm{\SIUnitSymbolMicro m}$ waist and travel through a cell of
length $7\text{\,}\mathrm{cm}$ at a temperature of
$95\text{\,}\mathrm{\SIUnitSymbolCelsius}$. It is assumed the atoms are
equally spread across the Zeeman states of both the $F=2$ and $F=1$ ground
states at room temperature [24].

##  III Phase modulation of a CW signal field

First we discuss the experimental implementation in which a CW signal field is
modulated by a pulsed control field. This allows the demonstration of a phase
shift that varies with the signal detuning $\Delta_{s}$ for a fixed control
detuning $\Delta_{c}$. The experiment is shown in Fig. 2, and uses the CW
signal preparation in the upper dashed box. The CW signal is generated by an
external cavity diode laser (MOGLabs CEL) whose frequency is scanned across
the region of interest shown in Fig. 1. Control pulses are generated by pulse
carving, amplifying and frequency doubling a pulsed telecomms laser operating
near $1552\text{\,}\mathrm{nm}$ (details are given in the supplementary
material). The control pulses are square, with $4\text{\,}\mathrm{ns}$
duration and $47\text{\,}\mathrm{nJ}$ pulse energy. The signal field average
power is $1.6\text{\,}\mathrm{mW}$. Since the timescale of the pulses is much
shorter than that of the signal frequency scan, we assume that $\Delta_{s}$ is
effectively constant for each control pulse.

Figure 2: Experimental setup, showing control and signal fields counter-
propagating through 87Rb vapour cell in free space, the generation scheme for
pulses and CW signal field, and detection region. In the CW case, a
beamsplitter (BS) is used to calibrate the external cavity diode laser (ECDL)
by saturation absorption spectroscopy of Rb and in the pulsed case an optical
chopper (OC) is inserted to modulate the control pulses. Signal and control
pulses can be generated by combining telecomms seed lasers with a wavelength-
division multiplexer (WDM) before pulse carving using a LiNbO3 intensity
modulator (IM) before amplification by erbium-doped fibre amplifier (EDFA).
The two wavelengths are separated at a second WDM, before frequency doubling
at periodically poled LiNbO3 crystals (PPLN) and sent to the experiment. A
delay line is used to synchronise the arrival time of signal and control
pulses. Further details of the pulse generation are given in the supplementary
material. Polarisation controllers (PC) and waveplates ($\lambda/2$,
$\lambda/4$) ensure that the polarisation of the fields is $\sigma^{+}$
circularly polarised at the vapour cell, and separates the fields at a
polarising beam splitter (PBS) for phase-sensitive signal detection by a
Franson interferometer. The purpose of fibre delay lines with duration $\tau$
is explained in the main text.

The signal field counter propagates with control pulses through the
$7\text{\,}\mathrm{cm}$ long, isotopically pure 87Rb vapour cell (Precision
Glassblowing) heated to $97.6\text{\,}\mathrm{\SIUnitSymbolCelsius}$. Half-
and quarter-waveplates are set so that the light is circularly polarised,
maximising the strength of the interaction [25]. Following this, the signal
field is separated into the detection region, which consists of a Franson
interferometer [26] (see Fig. 3) with polarisation control used to maximise
the fringes. Note that the interferometer has a free spectral range of
$200\text{\,}\mathrm{MHz}$, however for simplicity we consider only values of
$\Delta_{s}$ where the output through port 1 is maximised. A phase shift is
detected by a change in the intensity at each output.

Figure 3:  Time-bin interferometer detection scheme for CW operation. Regions
of the signal overlapping with a control pulse in the cell are shaded blue.
Blue hashing indicates regions of interference at the output.

We observe fringes described by Eqn. 3 at the interferometer output, from
which we extract $\Delta\phi$ and $T$ for a range of signal detunings and
control pulses detuned by $\Delta_{c}=$1.6\text{\,}\mathrm{GHz}$$. The results
are shown in Fig. 4. We observe that for $\Delta_{s}$ close to zero there is
increased absorption due to broadening of the single-photon resonance, and
that there is significant absorption from the two-photon effect around
$\Delta_{s}=$-3\text{\,}\mathrm{GHz}$$. We identify a range of
$\Delta_{s}=$-1.85\text{\,}\mathrm{GHz}-1.5\text{\,}\mathrm{GHz}$$ where the
transmission is highest and the phase shift is relatively flat. We achieve a
phase shift $\Delta\phi/\pi=$0.90(5)\text{\,}$$ with end-to-end transmission
$T=$74(2)\text{\,}\mathrm{\char 37\relax}$$, which includes the insertion loss
through the vapour cell of $4.5\text{\,}\mathrm{\char 37\relax}$.

Figure 4:  The upper three plots show the signal arriving at detectors 1
(green, dashed) and 2 (blue, solid) for various values of $\Delta_{s}$ (cf.
Fig. 3). In each of these traces, $t=0$ corresponds to the detection time of
the control pulse. Analysis of these peaks allows us to determine the
transmission through the vapour cell ($T$) as well as the phase shift
($\Delta\phi$), both of which are shown for the full scan of $\Delta_{s}$ in
the second row.

##  IV Phase modulating a pulsed signal field

In this section we discuss the phase modulation of a $\sim 100$ MHz-bandwidth
pulsed signal field, as required for clocked photonic systems. Signal pulses
are generated by the same means as for control pulses, as was described in the
previous section and detailed in the supplementary material. The same EDFA is
used simultaneously for both wavelengths, which reduces the total power in the
control. Signal field pulses are split, with a time delay of
$\tau=$10\text{\,}\mathrm{ns}$$, as shown in the lower dashed box of Fig. 2.
The delay line is used to ensure that both signal and control pulses arrive at
the cell together, such that the control overlaps with only the second pulse.
The control beam is chopped at $1\text{\,}\mathrm{kHz}$ to allow extraction of
the control-induced phase shift on a time-scale fast compared to the
interferometer drift. The optical amplifier gain is now divided between signal
and control fields, which reduces the available control power and consequently
the achievable phase shift.

We set $\Delta_{c}=$2\text{\,}\mathrm{GHz}$$, and ramp $\Delta_{s}$ over a
range $-2\text{\,}\mathrm{GHz}0\text{\,}\mathrm{GHz}$ at a frequency of
$1\text{\,}\mathrm{Hz}$, much slower than the chopper frequency. Ramping
$\Delta_{s}$ allows the single-photon absorption line to be resolved when the
control is off, so that the stability of the laser can be monitored, and
$\Delta_{s}$ can be determined precisely. Any signal pulses that experience a
partially chopped control are discarded from the analysis. As the chopper
frequency is three orders of magnitude larger than the signal ramp frequency,
the signal laser frequency is roughly constant across the control pulses
contained within a single window of the chopper.

The results for modulation of the pulsed signal field are shown in Fig. 5. The
cell is again held at $97.6\text{\,}\mathrm{\SIUnitSymbolCelsius}$, with
control pulse energy of $33\text{\,}\mathrm{nJ}$, signal pulse energy of
$200\text{\,}\mathrm{pJ}$, $\Delta_{c}=$2\text{\,}\mathrm{GHz}$$,
$\Delta_{s}=-$1.9\text{\,}\mathrm{GHz}$$, and pulse duration
$4\text{\,}\mathrm{ns}$. The switch of the interferometer output from detector
2 to detector 1 is due to the presence of the control field, and is a clear
indicator of phase modulation. Comparison of the two right hand peaks shows
that the control pulses are changing the shape and magnitude of the signal
pulse, due to a combination of absorption and dispersion induced by the
control.

Figure 5:  Demonstration of phase modulation for pulsed control and signal. We
show the interferometer output for control off (top panel) and on (bottom
panel) for signal pulses as described in the main text. Here we set $t=0$ as
the time that the first pulse arrives. Note the change in the central peak
between the control on and off, corresponding to
$\Delta\phi/\pi=$0.53(6)\text{\,}$$ The inset shows how we interpret the three
peaks arriving in at the signal detection arm shown in Fig. 2, with
polarisation controllers omitted for clarity. A pulse which has encountered a
control pulse in the vapour cell is indicated by blue shading. Interference on
the central peak at the output is indicated by blue-white hashing, and is used
to extract transmission and phase information as in the CW case.

The average transmission and phase shift across these pulses is calculated by
integration. This yields a transmission of $T=$81(6)\text{\,}\mathrm{\char
37\relax}$$, including vapour cell insertion loss, and phase shift of
$\Delta\phi/\pi=$0.53(6)\text{\,}$$. As the control pulse is the same duration
as the signal pulse, then the signal pulse does not experience a uniform
control intensity distribution in time. The incomplete temporal overlap
reduces the phase modulation effect, as can be seen in the edges of the the
detector 2 peak when the control pulse is on in Fig. 5.

##  V Conclusion

We have demonstrated low-loss, all-optical phase modulation operating on the
timescale of nanoseconds. Our scheme utilises the change in susceptibility of
warm 87Rb vapour in the presence of a strong control field, inducing a phase
shift in a weak signal field. We are able to demonstrate phase shifts close to
$\Delta\phi=\pi$, the critical value required for optical switching. Our
theoretical model predicts that reaching $|\Delta\phi|>10\pi$ is possible with
low loss, and that simultaneously increasing the temperature, control power
and control detuning is predicted to further improve the performance [27].

However a complete description of the phase modulation behaviour, accounting
for the time-variation of the control field, as well as the dependence on the
population of the various hyperfine states would allow further improvements to
the scheme [28]. This may come in the form of different detunings, or
optimisation of the hyperfine occupation by optical pumping [18].

We propose that improvement to phase shift and transmission can be made by
implementing the phase modulation scheme inside a hollow core optical fibre
[29] where the core is filled with 87Rb vapour to realise an in-fibre vapour
cell [30, 31, 32]. This would allow for significantly higher intensity, and
hence higher phase shift, of the control field, due to the small
($\sim$1000\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}$$) size of the mode
which can be maintained over lengths of optical fibre much greater than the
Rayleigh range. Such a device would have potential for scalability, due to
ability to easily integrate with other optical fibre systems [33]. An
alternative route to enhancement is by implementation in an optical cavity,
where the interaction of the control field is increased by a factor of the
cavity finesse [34, 35, 36]. This is particularly appealing for
microresonators where high quality factors ($Q>10^{6}$) are readily achievable
[37, 38, 39, 40, 41].

Finally, we suggest that phase modulation at other wavelengths may be possible
by applying our same scheme at other wavelengths, particularly those that
would be compatible with existing telecommunication infrastructure. Of
particular interest may be the
$5S_{1/2}\rightarrow{}5P_{3/2}\rightarrow{}4D_{5/2}$ transition also in 87Rb,
which has recently been used to demonstrate an EIT memory for classical light
[42] and single photons [43]. Here, the signal field is on the
$5P_{3/2}\rightarrow{}4D_{5/2}$ transition at $1529\text{\,}\mathrm{nm}$,
which is close to the telecommunication C-band.

## Acknowledgements

This work was supported by the UK Hub in Quantum Computing and Simulation,
part of the UK National Quantum Technologies Programme with funding from UKRI
EPSRC Grant EP/T001062/1. This work is partially funded by Innovate UK Quantum
Data Centre of the Future grant 10004793.

##  VI Supplementary material

###  VI.1 Equations of motion

Here we derive the equations of motion for the Rb atoms in the presence of
signal and control fields. This follows the similar derivations for the off-
resonant cascaded absorption (ORCA) ladder scheme [23] or the derivations for
a $\Lambda$-type Raman memory [21]. We begin with the Heisenberg equation

$\partial_{t}\rho=i[\mathcal{H},\rho].$ (4)

where $\rho$ is the density matrix, and $\mathcal{H}$ is the Hamiltonian of
the system. We write the latter as

$\mathcal{H}=\mathcal{H}_{\textnormal{A}}+\mathcal{H}_{\textnormal{AF}}$ (5)

where we have the atomic Hamiltonian in the absence of any external fields,

$\mathcal{H}_{\textnormal{A}}=\sum_{j}\hbar\omega_{j}\ket{j}\bra{j}=\sum_{j}\hbar\omega_{j}\rho_{jj}$
(6)

and the light-atom interaction Hamiltonian

$\mathcal{H}_{\textnormal{AF}}=-\mathbf{d}\cdot\mathbf{E}.$ (7)

The electric field is defined as the sum of the signal and control fields
($\mathbf{E}_{s}$ and $\mathbf{E}_{c}$ respectively),

$\displaystyle\begin{split}\mathbf{E}=&\mathbf{E}_{s}+\mathbf{E}_{c}\\\
=&\bm{\nu}_{\textnormal{c}}E_{\textnormal{c}}e^{i\omega_{\textnormal{c}}(t+\frac{z}{c})}+\bm{\nu}^{*}_{\textnormal{c}}E^{*}_{\textnormal{c}}e^{-i\omega_{\textnormal{c}}(t+\frac{z}{c})}\\\
&+\bm{\nu}_{\textnormal{s}}E_{\textnormal{s}}e^{i\omega_{\textnormal{s}}(t-\frac{z}{c})}+\bm{\nu}^{*}_{\textnormal{s}}E^{*}_{\textnormal{s}}e^{-i\omega_{\textnormal{s}}(t-\frac{z}{c})},\end{split}$
(8)

where $\bm{\nu}_{\textnormal{c}}$ and $\bm{\omega}_{\textnormal{c}}$
($\bm{\nu}_{\textnormal{s}}$ and $\bm{\omega}_{\textnormal{s}}$) are the
polarisation vector and angular frequency of the control (signal) fields
respectively. The dipole operator $\mathbf{d}$ is given by

$\displaystyle\begin{split}\mathbf{d}=&\sum_{j,k}\bra{j}\mathbf{d}\ket{k}\rho_{jk}\
+\bra{j}\mathbf{d}\ket{k}^{*}\rho^{\dagger}_{jk}\\\
=&\sum_{j,k}\bra{j}\mathbf{d}\ket{k}\ket{j}\bra{k}+\bra{j}\mathbf{d}\ket{k}^{*}\ket{k}\bra{j}.\end{split}$
(9)

We write the dipole matrix elements in the usual way,

$d_{jk}=\bra{j}\mathbf{d}\ket{k}.$ (10)

We will treat the Rb atom as a three level system, so that, $j,k=1,2,3$
(corresponding to the $S$, $P$, and $D$ states of interest respectively), and
the operators defined above can be written as $3\times 3$ matrices. We now
calculate each matrix element of the Hamiltonian

$\mathcal{H}_{jk}=\bra{j}\mathcal{H}\ket{k}.$ (11)

The atomic Hamiltonian has non-zero elements only when $j=k$. The electric
dipole interaction has odd parity, and only off-diagonal elements for the
electric dipole interaction remain are non-zero. The on-diagonal elements in
the density matrix are the populations, and the off-diagonal elements are the
coherences. For a weak signal field in which the number of rubidium atoms in
the mode is much larger than the number of photons. We can assume that any
change in populations of the states are negligible. This means that
$\rho_{11}=1$, as virtually the entire atomic population will remain in the
ground state, and thus $\rho_{22}$ and $\rho_{33}=0$ as there is negligible
population in the upper states. The change in the populations of the states
are also negligible and so $\partial_{t}\rho_{11}\ =~{}\partial_{t}\rho_{22}\
=~{}\partial_{t}\rho_{33}\ =~{}0$.

The ground state energy is defined as $\hbar\omega_{1}$ = 0 and it is assumed
there is no coupling between the excited states, so that $d_{13}=d_{31}$ = 0.
Applying the above assumptions we can write down the matrices

$\mathcal{H}=\begin{pmatrix}0&-\mathbf{d}_{12}.\mathbf{E}&0\\\
-\mathbf{d}_{21}.\mathbf{E}&\hbar\omega_{2}&-\mathbf{d}_{23}.\mathbf{E}\\\
0&-\mathbf{d}_{32}.\mathbf{E}&\hbar\omega_{3}\end{pmatrix},$ (12)
$\rho=\begin{pmatrix}1&\rho_{12}&\rho_{13}\\\ \rho_{21}&0&\rho_{23}\\\
\rho_{31}&\rho_{32}&0\end{pmatrix}.$ (13)

The angular frequencies are defined as $\omega_{2}=-\omega_{12}$,
$\omega_{3}=-\omega_{13}$ and $\omega_{23}=\omega_{2}-\omega_{3}$, where
$\omega_{jk}=\omega_{j}-\omega_{k}$ is the transition frequency of the
$\ket{j}\xrightarrow{}\ket{k}$ transition.

Doppler broadening can be taken into account by applying the chain rule to the
Heisenberg equation of motion, changing the partial differential to a normal
time differential [22, 20]. Considering only propagation and motion in the
z-direction leads to the equation

$\displaystyle\frac{\differential\rho(z,t)}{\differential t}$
$\displaystyle=\partial_{t}\rho|_{z}+\partial_{z}\rho|_{t}\frac{\differential
z}{\differential t}$ (14)
$\displaystyle=-i[\mathcal{H},\rho]-v\partial_{z}\rho$ (15)

Here $v$ is the velocity of an atom in the $z$ direction. The equations of
motion are then [21, 23],

$\displaystyle\begin{split}\differential_{t}\rho_{13}=&\left(-{\frac{\gamma_{3}}{2}}-i\omega_{13}-v\partial_{z}\right)\rho_{13}\\\
&+i\Omega^{*}_{\textnormal{c}}\rho_{12}e^{-i\omega_{\textnormal{c}}(t+\frac{z}{c})},\end{split}$
(16)

and

$\displaystyle\begin{split}\differential_{t}\rho_{12}=&\left(-{\frac{\gamma_{2}}{2}}-i\omega_{12}-v\partial_{z}\right)\rho_{12}\\\
&+i\kappa_{\textnormal{s}}S+i\Omega_{\textnormal{c}}\rho_{13}e^{i\omega_{\textnormal{c}}(t+\frac{z}{c})}.\end{split}$
(17)

Here we have introduced the Rabi frequency of the control field

$\Omega_{\textnormal{c}}={\frac{\mathbf{d}^{*}_{23}\cdot\bm{\nu}_{c}}{\hbar}}E_{\textnormal{c}},$
(18)

a related parameter for the signal field

$\kappa_{\textnormal{s}}={\frac{\mathbf{d}_{12}\cdot\bm{\nu}_{s}}{\hbar}},$
(19)

and we define for convenience

$S=E_{\textnormal{s}}e^{i\omega_{\textnormal{s}}(t-\frac{z}{c})}$ (20)

along with $\gamma_{2}$ as the decay rate of the $5P_{3/2}$ state and
$\gamma_{3}$ as the decay rate of the $5D_{5/2}$ state.

It is notable that the equation for $\rho_{12}$ has terms for the single-
photon transition frequency, the decay rate of the single-photon coherence,
and will depend on the single-photon Doppler shift. The equation for
$\rho_{13}$ has a term for the two-photon transition frequency, the decay rate
of the two-photon coherence, and will depend on the two-photon Doppler shift.

####  VI.1.1 Fourier transform solution

The equations are solved numerically by Fourier transform from time and space
to frequency and reciprocal space. This leads to two equations in frequency
and reciprocal space (using a hat to represent the transformed operator)

$\left({\frac{\gamma_{3}}{2}}+i(\omega+\omega_{13})+ivq\right)\hat{\rho}_{13}(q,\omega)=i\Omega^{*}_{\textnormal{c}}\hat{\rho}_{12}(q+k_{\textnormal{c}},\omega+\omega_{\textnormal{c}}),$
(21)
$\left({\frac{\gamma_{2}}{2}}+i(\omega+\omega_{12})+ivq\right)\hat{\rho}_{12}(q,\omega)=i\kappa_{\textnormal{s}}\hat{S}(q,\omega)+i\Omega_{\textnormal{c}}\hat{\rho}_{13}(q-k_{\textnormal{c}},\omega-\omega_{\textnormal{c}}).$
(22)

We now perform a co-ordinate transform on equation 21 and note that the signal
field, $S(q,\omega)$ is a pulse, with most of its energy concentrated around
the central frequency $\omega_{\textnormal{s}}$ and central wavevector
$k_{\textnormal{s}}$. Therefore, it is a good approximation to define
$q=k_{\textnormal{s}}$ and $\omega=\omega_{\textnormal{s}}$. We therefore have
that

$\left({\frac{\gamma_{3}}{2}}+i(\omega_{\textnormal{s}}+\omega_{\textnormal{c}}+\omega_{13})+iv(k_{\textnormal{s}}-k_{\textnormal{c}})\right)\hat{\rho}_{13}(k_{\textnormal{s}}-k_{\textnormal{c}},\omega_{\textnormal{s}}-\omega_{\textnormal{c}})=i\Omega^{*}_{\textnormal{c}}\hat{\rho}_{12}(k_{\textnormal{s}},\omega_{\textnormal{s}}),$
(23)
$\left({\frac{\gamma_{2}}{2}}+i(\omega_{\textnormal{s}}+\omega_{12})+ivk_{\textnormal{s}}\right)\hat{\rho}_{12}=i\kappa_{\textnormal{s}}\hat{S}+i\Omega_{\textnormal{c}}\hat{\rho}_{13}.$
(24)

We now note that that $\omega_{13}=-\omega_{21}-\omega_{32}$, which allows us
to write the signal detuning as
$\Delta_{\textnormal{s}}=-\omega_{21}+\omega_{\textnormal{s}}$ and control
detuning as $\Delta_{\textnormal{c}}=-\omega_{32}+\omega_{\textnormal{c}}$. We
can now re-write the above expressions as

$\left(i{\frac{\gamma_{3}}{2}}-(\Delta_{\textnormal{s}}+\Delta_{\textnormal{c}})-v(k_{\textnormal{s}}-k_{\textnormal{c}})\right)\hat{\rho}_{13}=-\Omega^{*}_{\textnormal{c}}\hat{\rho}_{12},$
(25)
$\left(i{\frac{\gamma_{2}}{2}}-\Delta_{\textnormal{s}}-vk_{\textnormal{s}}\right)\hat{\rho}_{12}=-\kappa_{\textnormal{s}}\hat{S}-\Omega_{\textnormal{c}}\hat{\rho}_{13}.$
(26)

Rearranging equation 25 allows an expression for $\hat{\rho}_{13}$ in terms of
$\hat{\rho}_{12}$ to be obtained

$\hat{\rho}_{13}=\frac{\Omega^{*}_{\textnormal{c}}\hat{\rho}_{12}}{(-i{\frac{\gamma_{3}}{2}}+(\Delta_{\textnormal{s}}+\Delta_{\textnormal{c}})+v(k_{\textnormal{s}}-k_{\textnormal{c}}))}.$
(27)

Substituting in for $\hat{\rho}_{13}$ in equation 26 and re-arranging provides
an expression in terms of $\hat{\rho}_{12}$

$\hat{\rho}_{12}=\frac{\kappa_{\textnormal{s}}\hat{S}(q,\omega)}{-i{\frac{\gamma_{2}}{2}}+\Delta_{\textnormal{s}}+vk_{\textnormal{s}}-\frac{|\Omega^{2}_{\textnormal{c}}|}{-i{\frac{\gamma_{3}}{2}}+(\Delta_{\textnormal{s}}+\Delta_{\textnormal{c}})+v(k_{\textnormal{s}}-k_{\textnormal{c}})}}.$
(28)

Next, we define the polarisation density as

$\mathbf{P}(z,t)=n_{\textnormal{d}}\mathbf{d}_{12}\rho_{12}(z,t),$ (29)

where $n_{\textnormal{d}}$ is the number density of the vapour. Including
Doppler broadening the expression is given by

$\mathbf{P}(z,t,v)=n_{\textnormal{d}}\mathbf{d}_{12}\int
p(v)\rho_{12}(z,t,v)\differential v,$ (30)

where $p(v)$ is the normal distribution of velocities of particle with mass
$m$, and temperature, $T_{\text{cell}}$

$p(v)=\sqrt{{\frac{m}{2\pi
k_{B}T_{\text{cell}}}}}e^{\frac{-mv^{2}}{2k_{B}T_{\text{cell}}}}.$ (31)

The polarisation density is given by

$\mathbf{P}(z,\omega)=\epsilon_{0}\chi(\omega)\mathbf{E_{\textnormal{s}}}(z,\omega).$
(32)

Fourier transforming the above polarisation density relations to the frequency
and reciprocal space domains leads to the expressions

$\hat{\mathbf{P}}(q,\omega,v)=n_{\textnormal{d}}\mathbf{d}_{12}\int
p(v)\hat{\rho}_{12}(q,\omega,v)\differential v,$ (33)
$\hat{\mathbf{P}}(q,\omega)=\epsilon_{0}\chi(\omega)\mathbf{\hat{E}_{\textnormal{s}}}(q,\omega).$
(34)

Substituting in for $\hat{\rho}_{12}(v)$ and re-arranging for the
susceptibility, $\chi$, leads to the following expression in the steady-state
for the three-level system

$\chi(\omega)={\frac{n_{\textnormal{d}}d_{12}\kappa_{\textnormal{s}}}{\epsilon_{0}}}\int{\frac{p(v)}{-i{\frac{\gamma_{2}}{2}}+\Delta_{\textnormal{s}}+vk_{\textnormal{s}}-\frac{|\Omega^{2}_{\textnormal{c}}|}{-i{\frac{\gamma_{3}}{2}}+(\Delta_{\textnormal{s}}+\Delta_{\textnormal{c}})+v(k_{\textnormal{s}}-k_{\textnormal{c}})}}}\differential
v.$ (35)

This equation can then be solved numerically to find the real and imaginary
parts of the susceptibility, by numerically integrating over all velocities.
The real and imaginary parts of the susceptibility are related to the complex
refractive index by

$\tilde{n}=n-iK\approx\sqrt{1+\chi}.$ (36)

Here we have introduced the real numbers $n$ and $K$ to distinguish the real
and imaginary parts of the complex refractive index. The former provides
information on phase velocity and phase shifts, whereas the latter is the
extinction coefficient.

Calculating the refractive index with and without the presence of the control
field allows the change in the refractive index, and therefore the change in
the phase to be determined using the equation

$\Delta\phi={\frac{2\pi L}{\lambda_{0}}}(n_{\mathrm{c}}-n_{0})$ (37)

where $n_{\mathrm{c}}$ ($n_{0}$) is the refractive index with the control
field switched on (off), $\lambda_{0}$ is the signal field wavelength, and $L$
is the interaction length. For pulses that are longer than the cell this will
be the length of the vapour cell. The absorption coefficient $\alpha$ relates
to the extinction coefficient by the following relation

$\alpha=\frac{4\pi K}{\lambda},$ (38)

for signal frequency, $f$, extinction coefficient, $K$, and speed of light,
$c$. Assuming that the vapour cell is much shorter than the control pulse,
then the absorption coefficient experienced by the signal will be uniform
across the cell. The transmission in this case is then given by the Beer-
Lambert law, giving an exponential decay of the cell transmission

$T=e^{-\alpha L}.$ (39)

###  VI.2 Pulse generation

Figure 6: Setup used to generate signal and control pulses. SMF, single-mode
fibre; PMF, polarisation-maintaining fibre; FPC, fibre polarisation
controller; AWG, arbitary waveform generator; IM, intensity modulator; PD,
photodiode; PPLN, periodically poled lithium niobate; WDM, wavelength division
multiplexer; EDFA, erbium doped fibre amplifier.

The experimental setup starts with generating two trains of pulses. Strong
control pulses at 776 $\mathrm{nm}$ and weak signal pulses at 780.2
$\mathrm{nm}$. The pulse generation scheme is shown in Fig. 6. Two tunable
fibre coupled seed lasers in the telecomms C band at 1552 $\mathrm{nm}$
(control) and 1560 $\mathrm{nm}$ (signal) (ID Photonics CoBrite DX1, and
Eblana EP1560-0-DM-B05-FM, respectively) are used as seed lasers. Both
wavelengths are combined using a wavelength division multiplexer (Oz Optics
WDM-12N-111-1552/1560-8/125-PPP-60-3A3A3A-1-1). The seed lasers are pulsed
carved using a waveguided Lithium Niobate intensity modulator (EOSpace - AZ-
DS5-10-PFA-PFA-LV), stabilised by a DC bias controller (MBC-HER-PD-3A-0V) at
the operating point for minimum transmission. Fast electrical pulses are
applied by a radio-frequency (RF) pulse generator (Picoquant PPG512) with the
repetition rate controlled by an arbitrary waveform generator. This allows a
stable train of signal and control pulses to be carved, and the pulse duration
can be varied between 600 ps and 4 $\mathrm{ns}$, or longer.

The signal and control pulses are amplified in an Erbium doped fibre amplifier
(EDFA). The amplified pulses are then split out using a second wavelength
division multiplexer (Oz Optics
WDM-12N-111-1552/1560-8/125-SSS-60-3A3A3A-1-1). The 1552 $\mathrm{nm}$ control
field passes into a waveguided periodically poled lithium niobate crystal
(Covesion WGP-H-1552-40) for frequency doubling to 776 $\mathrm{nm}$. The 1560
$\mathrm{nm}$ signal field is coupled into free-space and is focused into a
bulk periodically poled lithium niobate crystal (Covesion MSHG1550-1.0-40) for
frequency doubling to 780.2 $\mathrm{nm}$. Any remaining 1560 $\mathrm{nm}$
light is removed by a short-pass filter. The 780 $\mathrm{nm}$ signal light
passes through a variable optical delay line before being coupled back into
fibre. The optical delay line allows the relative delay between the signal and
control fields to be adjusted to ensure that pulses meet inside the Rb vapour
cell.

## References

  * Biberman and Bergman [2012] A. Biberman and K. Bergman, Optical interconnection networks for high-performance computing systems, Reports on Progress in Physics 75, 046402 (2012).
  * Miscuglio and Sorger [2020] M. Miscuglio and V. J. Sorger, Photonic tensor cores for machine learning, Applied Physics Reviews 7, 031404 (2020), https://pubs.aip.org/aip/apr/article-pdf/doi/10.1063/5.0001942/14577276/031404_1_online.pdf .
  * Gisin and Thew [2007] N. Gisin and R. Thew, Quantum communication, Nature Photonics 1, 165 (2007).
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Amin  _et al._ [2021] R. Amin, R. Maiti, J. B. Khurgin, and V. J. Sorger, Performance analysis of integrated electro-optic phase modulators based on emerging materials, IEEE Journal of Selected Topics in Quantum Electronics 27, 1 (2021).
  * Wada [2004] O. Wada, Femtosecond all-optical devices for ultrafast communication and signal processing, New Journal of Physics 6, 183 (2004).
  * Slussarenko and Pryde [2019] S. Slussarenko and G. J. Pryde, Photonic quantum information processing: A concise review, Applied Physics Reviews 6, https://doi.org/10.1063/1.5115814 (2019).
  * Pittman and Franson [2002] T. B. Pittman and J. D. Franson, Cyclical quantum memory for photonic qubits, Phys. Rev. A 66, 062302 (2002).
  * Makino  _et al._ [2016] K. Makino, Y. Hashimoto, J.-i. Yoshikawa, H. Ohdan, T. Toyama, P. van Loock, and A. Furusawa, Synchronization of optical photons for quantum information processing, Science Advances 2, 10.1126/sciadv.1501772 (2016), arXiv:1509.04409 .
  * Kaneda  _et al._ [2017] F. Kaneda, F. Xu, J. Chapman, and P. G. Kwiat, Quantum-memory-assisted multi-photon generation for efficient quantum information processing, Optica 4, 1034 (2017), arXiv:1704.00879 .
  * Collins  _et al._ [2013] M. Collins, C. Xiong, I. Rey, T. Vo, J. He, S. Shahnia, C. Reardon, T. Krauss, M. Steel, A. Clark, and B. Eggleton, Integrated spatial multiplexing of heralded single-photon sources, Nature Communications 4, 2582 (2013).
  * Mendoza  _et al._ [2016] G. J. Mendoza, R. Santagati, J. Munns, E. Hemsley, M. Piekarek, E. Martín-López, G. D. Marshall, D. Bonneau, M. G. Thompson, and J. L. O’Brien, Active temporal and spatial multiplexing of photons, Optica 3, 127 (2016).
  * Francis-Jones  _et al._ [2016] R. J. A. Francis-Jones, R. A. Hoggarth, and P. J. Mosley, All-fiber multiplexed source of high-purity single photons, Optica 3, 1270 (2016).
  * Li  _et al._ [2020] J.-P. Li, J. Qin, A. Chen, Z.-C. Duan, Y. Yu, Y. Huo, S. Höfling, C.-Y. Lu, K. Chen, and J.-W. Pan, Multiphoton Graph States from a Solid-State Single-Photon Source, ACS Photonics 7, 1603 (2020).
  * Pankovich  _et al._ [2023] B. Pankovich, A. Kan, K. H. Wan, M. Ostmann, A. Neville, S. Omkar, A. Sohbi, and K. Brádler, High photon-loss threshold quantum computing using ghz-state measurements (2023), arXiv:2308.04192 [quant-ph] .
  * Hendrickson  _et al._ [2010] S. M. Hendrickson, M. M. Lai, T. B. Pittman, and J. D. Franson, Observation of two-photon absorption at low power levels using tapered optical fibers in rubidium vapor, Phys. Rev. Lett. 105, 173602 (2010).
  * Jones  _et al._ [2015] D. E. Jones, J. D. Franson, and T. B. Pittman, Ladder-type electromagnetically induced transparency using nanofiber-guided light in a warm atomic vapor, Phys. Rev. A 92, 043806 (2015).
  * Finkelstein  _et al._ [2018] R. Finkelstein, E. Poem, O. Michel, O. Lahad, and O. Firstenberg, Fast, noise-free memory for photon synchronization at room temperature, Science Advances 4, eaap8598 (2018), https://www.science.org/doi/pdf/10.1126/sciadv.aap8598 .
  * Metcalf and van der Straten [1999] H. J. Metcalf and P. van der Straten,  _Laser cooling and trapping_, 1st ed. (Springer-Verlag New York, 1999).
  * Page [2022] C. Page,  _Atomic Switchable Cavity Quantum Memory Through Induced Two Photon Dispersion_ , Ph.D. thesis, University of Bath (2022).
  * Kaczmarek [2017] K. T. Kaczmarek,  _Towards an integrated noise-free Quantum Memory_ , Ph.D. thesis, University of Oxford (2017).
  * Firstenberg  _et al._ [2008] O. Firstenberg, M. Shuker, R. Pugatch, D. R. Fredkin, N. Davidson, and A. Ron, Theory of thermal motion in electromagnetically induced transparency: Effects of diffusion, Doppler broadening, and Dicke and Ramsey narrowing, Physical Review A 77, 043830 (2008).
  * Nunn [2008] J. Nunn,  _Quantum Memory in Atomic Ensembles_ , Ph.D. thesis, St. John’s College, Oxford (2008).
  * Siddons  _et al._ [2008] P. Siddons, C. S. Adams, C. Ge, and I. G. Hughes, Absolute absorption on rubidium D lines: comparison between theory and experiment, Journal of Physics B: Atomic, Molecular and Optical Physics 41, 155004 (2008), arXiv:0805.1139v1 .
  * Olson  _et al._ [2006] A. J. Olson, E. J. Carlson, and S. K. Mayer, Two-photon spectroscopy of rubidium using a grating-feedback diode laser, American Journal of Physics 74, 218 (2006).
  * Franson [1989] J. D. Franson, Bell inequality for position and time, Phys. Rev. Lett. 62, 2205 (1989).
  * Lahad and Firstenberg [2017] O. Lahad and O. Firstenberg, Induced cavities for photonic quantum gates, Phys. Rev. Lett. 119, 113601 (2017).
  * Burdekin  _et al._ [2020] P. Burdekin, S. Grandi, R. Newbold, R. A. Hoggarth, K. D. Major, and A. S. Clark, Single-photon-level sub-doppler pump-probe spectroscopy of rubidium, Phys. Rev. Appl. 14, 044046 (2020).
  * Yu and Knight [2016] F. Yu and J. C. Knight, Negative Curvature Hollow-Core Optical Fiber, IEEE Journal of Selected Topics in Quantum Electronics 22, 146 (2016).
  * Perrella  _et al._ [2018] C. Perrella, P. S. Light, S. A. Vahid, F. Benabid, and A. N. Luiten, Engineering Photon-Photon Interactions within Rubidium-Filled Waveguides, Physical Review Applied 9, 044001 (2018).
  * Slepkov  _et al._ [2008] A. D. Slepkov, A. R. Bhagwat, V. Venkataraman, P. Londero, and A. L. Gaeta, Generation of large alkali vapor densities inside bare hollow-core photonic band-gap fibers, Optics Express 16, 18976 (2008).
  * Sprague  _et al._ [2014] M. R. Sprague, P. S. Michelberger, T. F. M. Champion, D. G. England, J. Nunn, X.-M. Jin, W. S. Kolthammer, A. Abdolvand, P. S. J. Russell, and I. A. Walmsley, Broadband single-photon-level memory in a hollow-core photonic crystal fibre, Nature Photonics 8, 287 (2014).
  * Suslov  _et al._ [2021] D. Suslov, M. Komanec, E. R. Numkam Fokoua, D. Dousek, A. Zhong, S. Zvánovec, T. D. Bradley, F. Poletti, D. J. Richardson, and R. Slavík, Low loss and high performance interconnection between standard single-mode fiber and antiresonant hollow-core fiber, Scientific Reports 11, 1 (2021).
  * Gorshkov  _et al._ [2007] A. V. Gorshkov, A. André, M. D. Lukin, and A. S. Sørensen, Photon storage in $\Lambda$-type optically dense atomic media. i. cavity model, Phys. Rev. A 76, 033804 (2007).
  * O’Shea  _et al._ [2013] D. O’Shea, C. Junge, J. Volz, and A. Rauschenbeutel, Fiber-optical switch controlled by a single atom, Phys. Rev. Lett. 111, 193601 (2013).
  * Nunn  _et al._ [2017] J. Nunn, J. H. D. Munns, S. Thomas, K. T. Kaczmarek, C. Qiu, A. Feizpour, E. Poem, B. Brecht, D. J. Saunders, P. M. Ledingham, D. V. Reddy, M. G. Raymer, and I. A. Walmsley, Theory of noise suppression in $\mathrm{\Lambda}$-type quantum memories by means of a cavity, Phys. Rev. A 96, 012338 (2017).
  * Del’Haye  _et al._ [2007] P. Del’Haye, A. Schliesser, O. Arcizet, T. Wilken, R. Holzwarth, and T. J. Kippenberg, Optical frequency comb generation from a monolithic microresonator, Nature 450, 1214 (2007).
  * Lee  _et al._ [2013] H. Lee, M.-G. Suh, T. Chen, J. Li, S. A. Diddams, and K. J. Vahala, Spiral resonators for on-chip laser frequency stabilization, Nature Communications 4, 2468 (2013).
  * Sumetsky  _et al._ [2010] M. Sumetsky, Y. Dulashko, and R. S. Windeler, Optical microbubble resonator, Opt. Lett. 35, 898 (2010).
  * Pfeifer  _et al._ [2022] H. Pfeifer, L. Ratschbacher, J. Gallego, C. Saavedra, A. Faßbender, A. von Haaren, W. Alt, S. Hofferberth, M. Köhl, S. Linden, and D. Meschede, Achievements and perspectives of optical fiber fabry–perot cavities, Applied Physics B 128, 29 (2022).
  * Svela  _et al._ [2020] A. Ø. Svela, J. M. Silver, L. Del Bino, S. Zhang, M. T. M. Woodley, M. R. Vanner, and P. Del’Haye, Coherent suppression of backscattering in optical microresonators, Light: Science & Applications 9, 204 (2020).
  * Thomas  _et al._ [2022] S. E. Thomas, S. Sagona-Stophel, Z. Schofield, I. A. Walmsley, and P. M. Ledingham, A single-photon-compatible telecom-c-band quantum memory in a hot atomic gas (2022), arXiv:2211.04415 [quant-ph] .
  * Thomas  _et al._ [2023] S. E. Thomas, L. Wagner, R. Joos, R. Sittig, C. Nawrath, P. Burdekin, T. Huber-Loyola, S. Sagona-Stophel, S. Höfling, M. Jetter, P. Michler, I. A. Walmsley, S. L. Portalupi, and P. M. Ledingham, Deterministic storage and retrieval of telecom quantum dot photons interfaced with an atomic quantum memory (2023), arXiv:2303.04166 [quant-ph] .