# Clustering of timed sequences – Application to the analysis of care pathways

Thomas Guyet†, Pierre Pinson‡ and Enoal Gesny†,‡

( $\dagger$. AIstroSight, Inria, Hospices Civils de Lyon, Université Claude
Bernard Lyon 1, 56 Bld Niels Bohr, 69603 Villeurbanne, France
<EMAIL_ADDRESS>
$\ddagger$. Clinical Research Unit, AP-HP, Cochin Academic Hospital, France. )

###### Abstract

Improving the future of healthcare starts by better understanding the current
actual practices in hospitals. This motivates the objective of discovering
typical care pathways from patient data. Revealing homogeneous groups of care
pathways can be achieved through clustering. The difficulty in clustering care
pathways, represented by sequences of timestamped events, lies in defining a
semantically appropriate metric and clustering algorithms.

In this article, we adapt two methods developed for time series to time
sequences: the drop-DTW metric and the DBA approach for the construction of
averaged time sequences. These methods are then applied in clustering
algorithms to propose original and sound clustering algorithms for timed
sequences. This approach is experimented with and evaluated on synthetic and
real use cases.

##  1 Introduction

Timed sequences, composed of dated events, are types of data frequently
encountered in various application domains such as computer log analysis, life
course analysis in social science [6], tourism trajectories [1] or healthcare
pathway analysis [16]. A peculiarity of this type of data is that they do not
have a natural vector representation. In particular, they are formed by two
very different dimensions: symbolic and temporal; and, on the other hand,
their size generally varies from one example to another (e.g., in a healthcare
pathway, some patients have more care events than others). This difficulty in
modeling data in a vector space makes it challenging to adapt standard data
analysis methods (classification, clustering, or forecasting).

The adaptation of data analysis methods raises the question of formalizing a
topology of the objects to be analyzed. Most of the data analytic methods make
the assumption that objects lies in a metric space. This topology assumption
provides intuitive ways to generalize sparse and noisy observations.
Formalizing such a space for timed sequences involves defining a metric
between timed sequences. This problem is very similar to that encountered in
time series analysis (temporal data where the values are numerical, not event-
based, and generally with a regular and dense temporal sampling). Thus, for
time series analysis, a large collection of metrics has been proposed and
allows them to be analyzed [11, 21]. Then, the challenges becomes to identify
the metric that has the properties appropriate for the analysis task one
wishes to carry out (e.g., the properties of a distance) and whose semantics
will give the expected results. Intuitively, a similarity metric (resp.
dissimilarity metric) should produce a high similarity value (resp. low
dissimilarity value) when two examples are considered similar by domain
experts.

For sequential data, few studies have focused on adapting data analysis
methods. [8] and [9] are based on tensor representation of sequential data,
while [22] or [12] use embedding representations. But only a few methods
represent metric time [15, 6]. Beside these end-to-end methods, some works
were dedicated to metrics between temporal sequences [7, 19, 4] or averaging
[3, 2, 20]. Despite these proposals, metrics whose semantics can better take
into account the temporal dimension are still lacking (and in particular the
delays between the occurrence of events). Indeed, most current solutions
generally only take into account the sequentiality in the data.

In this work, we extend the drop-DTW metric [5] which was proposed to evaluate
distances between sequences by taking into account both sequential aspects, by
adapting Dynamic Time Warping [18], and by adding support for meaningless
events, which can be dropped. This metric has the semantic characteristics
that seem appropriate for our application case – i.e., the clustering of care
pathways, but as this metric is not a distance, it does not allow for the
construction of an average sequence nor to be used in classical clustering
methods such as the K-Means algorithm or hierarchical clustering.

The main contribution of this article is the proposal of an algorithm to
compute a mean timed sequences based on drop-DTW. This algorithm is inspired
by the DBA method [14] for time series. The construction of an average
sequence meets both the needs of clustering algorithms (hierarchical
clustering or K-means) and the needs of interpretation of sequence clusters.

This method is applied to discover typical care-pathways from Electronic
Health Record (EHR) datasets. In the context of the OPTISOINS
project111https://www.bernoulli-lab.fr/project/optisoins/, we were
particularly interested in patients who had undergone pulmonary resection
surgery for oncological reasons. Care pathways were reconstructed from data
collected in the EHR data warehouse of the Greater Paris Hospital (AP-HP). The
care pathway is the timed sequence of the key stages of managing the patient
surgery from diagnosis to post-surgery treatments.

The reminder of this article is organized as follow. We first introduce the
notations related to timed sequences. Then, Section 3 introduces the
adaptation of drop-DTW to timed sequences. Section 4 presents the adaptation
of the DBA algorithm for averaging timed sequences based on the drop-DTW. This
averaging method is required in classical clustering algorithms that are
briefly recalled in Section 5. Finally, Sections 6 and 7 present experiments
and results on synthetic data and our use case (OPTISOINS project)
respectively.

##  2 Times Sequences and Probabilistic Timed Sequences

In this section, we introduce the notations for timed sequences and
probabilistic timed sequences. The probabilistic version of timed sequences
enables to adapt the metrics defined for time series to timed sequences (see
next section).

Let $\Sigma$ be a finite set of $n$ event types. A timed sequence is a
sequence of pairs $S=\langle(s_{i},t_{i})\rangle$ where $s_{i}\in\Sigma$ is an
event type and $t_{i}\in\mathbb{R}$ is the event’s timestamp.
$S[i]=(s_{i},t_{i})$ denotes the $i$-th event in the sequence. $|S|$ is the
length of the timed sequences, i.e. its number of events. By convention,
events in a sequence are ordered by their timestamps. Several events may occur
at the same time if their event types are different, i.e.
$(s_{i},t)=(s_{j},t)\Rightarrow s_{i}=s_{j}$.

We propose embedding the space of timed sequences, introduced in the previous
section, into a larger space where events are represented by their one-hot
encoding, which is a vector representation. Instead of representing a sequence
of events, we represent a sequence of event types distributions:

$\langle([d^{1}_{1},\dots d^{n}_{1}],t_{1}),\dots,([d^{1}_{K},\dots
d^{n}_{K}],t_{K})\rangle$

where $\bm{d}\cdot=[d^{1}_{\cdot},\dots,d^{n}_{\cdot}]\in[0,1]^{n}$, such that
$\sum_{i}d^{i}_{\cdot}=1$. $\bm{d}_{\cdot}$ can be seen as a probability
distribution of event types. Thus, events can be represented in a vector space
and the representation of a timed sequence becomes close to the representation
of a (multivariate, irregular) time series.

Figure 1: Illustration of the timed sequences of Example 1 (on the left) and
the cost matrix for the alignment of $s_{1}$ and $s_{2}$ according to the
drop-DTW (on the right). The alignment is illustrated by the grey bars on the
left and by the colored cells in the matrix (on the right). The red-crosses
illustrate that $D$ has been dropped in this alignment, the cells at the
corner are dark because they to no satisfy the constraint $\tau=3.5$.

The example 1 illustrates the representation of a timed sequence as a
probabilistic timed sequence. Conversely, most of the probabilistic timed
sequences can not be represented by a time sequence.

######  Example 1 (Probabilistic representation of a timed sequence).

Let $\Sigma={A,B,C,D}$ be a set of event types, and
$s_{1}=\langle(A,1),(C,2),(B,4.5)\rangle$,
$s_{2}=\langle(C,0),(C,2),(D,3),(A,4),(B,5)\rangle$, and
$s_{3}=\langle(A,0),(C,1),(C,2),(B,4)\rangle$ be three timed sequences
illustrated in Figure 1. These sequences can be represented as probabilistic
sequences:

$\begin{array}[]{l}s_{1}=\langle([\mathsf{1,0,0,0}],1),([\mathsf{0,0,1,0}],2),([\mathsf{0,1,0,0}],4.5)\rangle\\\
s_{2}=\langle([\mathsf{0,0,1,0}],0),([\mathsf{0,0,1,0}],2),([\mathsf{0,0,0,1}],3),([\mathsf{1,0,0,0}],4),([\mathsf{0,1,0,0}],5)\rangle\\\
s_{3}=\langle([\mathsf{1,0,0,0}],0),[\mathsf{0,0,1,0}],1),(B,2),([\mathsf{0,0,1,0}],2),([\mathsf{0,1,0,0}],4)\rangle\end{array}$

##  3 Comparing Timed Sequences

This section introduces the drop-DTW for timed sequences. We first remind the
definition of DTW and drop-DTW and then we adapt the latter for timed
sequences.

###  3.1 DTW and drop-DTW

Dynamic Time Warping (DTW) [18] computes the optimal alignment between two
series subject to certain constraints. Here, we adapt these definitions to the
case of timed sequences. Let
$X=\langle(x_{1},t^{x}_{1}),\dots,(x_{N},t^{x}_{N})\rangle$ and
$Z=\langle(z_{1},t^{z}_{1})\dots,(z_{K},t^{z}_{K})\rangle$ be the input
sequences, where $N,K$ are the respective lengths of the sequences. An
alignment between the two sequences is defined as a binary matrix,
$M\in\\{0,1\\}^{K\times N}$ where $M_{i,j}=1$ if $z_{i}$ is matched to
$x_{j}$, and $0$ otherwise. Matching an element $z_{i}$ to an element $x_{j}$
incurs a cost $C_{i,j}$ between the $i$-th and $j$-th elements of the
sequences. DTW finds the alignment $M^{*}$ between sequences $Z$ and $X$ that
minimizes the global cost:

$M^{*}=\operatorname*{argmin}_{M\in\mathcal{M}}\langle
M,C\rangle=\operatorname*{argmin}_{M\in\mathcal{M}}\sum_{i,j}M_{i,j}C_{i,j},$

where $\langle M,C\rangle$ is the Frobenius inner product and $\mathcal{M}$ is
the set of all feasible alignments satisfying the following constraints:
monotonicity, continuity, and end-point correspondence
($M_{1,1}=M_{K,N}=1,\forall M\in\mathcal{M}$).

drop-DTW [5] extends the set of feasible alignments, $\mathcal{M}$, to those
adhering only to the monotonicity constraint. Consequently, unlike DTW, events
can be dropped from the alignment process. drop-DTW finds the optimal temporal
alignment while allowing for the removal of certain elements. To account for
unmatched elements in the alignment cost calculation, Dvornik et al. [5]
introduce a new drop cost: $\delta\in\mathbb{R}^{+}$.222In the original
approach of Dvornik et al., the cost value depends on the event type. We use a
constant value for the sake of simplification. The optimal matching can then
be defined as follows:

$M^{*}=\operatorname*{argmin}_{M\in\overline{\mathcal{M}}}\langle
M,C\rangle+\delta.\left(P_{z}(M)+P_{x}(M)\right)$

where $\overline{\mathcal{M}}$ is the set of binary matrices satisfying only
the monotonicity constraint, and $P_{x}(M)$ is a vector whose $j$th element is
$1$ if $M_{:,j}=\bm{0}$ and $0$ otherwise; $P_{z}(M)$ is defined similarly,
but along the rows. Similar to DTW, drop-DTW can be efficiently evaluated
using dynamic programming techniques.

###  3.2 Adaptation of drop-DTW for timed sequences

The proposed adaptation of drop-DTW focuses on proposing a metric between
events and adding temporal constraints to the drop-DTW algorithm. The Dynamic
Time Warping algorithm requires defining a cost between two events to compare
sequences. To take into account both time and events, we decided to weigh: a
time distance cost $c_{t}$, and an event difference cost $c_{e}$. Having
introduced the probabilistic representation of timed sequences, the distance
between two events $e_{a}=(\bm{d}_{a}=[d^{1}_{a},\dots d^{n}_{a}],t_{a})$ and
$e_{b}=(\bm{d}_{b}=[d^{1}_{b},\dots d^{n}_{b}],t_{b})$ is given by:

$d(e_{a},e_{b})=\sqrt{p_{t}.||\bm{d}_{a}-\bm{d}_{b}||_{2}^{2}+p_{e}.(t_{b}-t_{a})^{2}}$

where $p_{t},p_{e}\in\mathbb{R}_{+}$ are weights of two terms: a distance
between probabilistic events and a time gap. The choice of this metric turns
out to consider that time events are in an Euclidean space of dimension $n+1$.

We can now adapt the metric between two timed sequences. We have chosen to use
drop-DTW for its ability to handle aberrant or missing events. Algorithm 1
presents the method for calculating the proposed distance between two timed
sequences, where the cost matrix is computed from the event distance, $d$. The
algorithm finds the best alignment of sequence events to minimize the overall
metric. The use of dynamic programming technique ensures the efficiency of the
method. For the sake of computing efficiency, we also add constraints on the
possible event pairings by allowing matches only between events that are close
in time, _i.e._ with a delay between their occurrence lower than a predefined
threshold $\tau$.

Data : $C\in\mathbb{R}^{M\times N}$: cost matrix; $\delta\in\mathbb{R^{+}}$:
drop cost; $\sigma\in\mathbb{N^{+}}$, $\tau\in\mathbb{R^{+}}$:
sequential/temporal proximity constraints

1 $D^{+}_{0,0}\leftarrow 0$ ; $D^{+}_{i,0}\leftarrow+\infty$ ;
$D^{+}_{0,j}\leftarrow+\infty$ ;

2 $D^{-}_{0,0}\leftarrow 0$ ; $D^{-}_{i,0}\leftarrow i\times\delta$ ;
$D^{-}_{0,j}\leftarrow j\times\delta$ ;

3 for  _$i=1$ to $N+1$_ do

4 for  _$j=max(1,i-\sigma+1)$ to $min(m+1,i+\sigma)$_ do // Sakoe-Chiba
constraints

5 if  _$distdate(i-1,j-1)\leq\tau$_ then // Temporal constraint

$D^{+}_{i,j}\leftarrow C_{i-1,j-1}+min(D_{i-1,j},D_{i,j-1},D^{+}_{i-1,j-1})$ ;

// Alignment cost

$D^{-}_{i,j}\leftarrow\delta+D_{i,j-1}$ ;

// Drop cost

$D_{i,j}\leftarrow min(D^{+}_{i,j},D^{-}_{i,j})$ ;

// Final cost

6

7 else

8 $D_{i,j}\leftarrow+\infty$;

9

return _$D_{M,N}$_

Algorithm 1 drop-DTW with proximity constraints

######  Example 2 (Illustration of drop-cost usefulness).

Considering the sequences from Example 1, sequences $s_{1}$ and $s_{2}$ seem
close and may be of interest to be grouped together, while sequence 3 seems
very different. Without drop-cost, sequences 2 and 3 would be grouped before
sequences $s_{1}$ and $s_{2}$. By using the drop-cost, the last event of
sequence 3 would not be considered, and the cost would simply be replaced by
$\delta$. Thus, they would be judged more similar than sequences 2 and 3.

##  4 Construction of an Average Timed Sequence

Computing an average time d sequence consistent with the notion of similarity
is essential for clustering algorithms. In the following, we propose a method
inspired by the DBA algorithm [13] to create average sequences in terms of
drop-DTW, taking into account the possibility of removal. DBA addresses the
problem of generating an average time series in terms of DTW – which does not
hold distance properties. In this case, the generalized average is defined as
a time series that minimizes the global inertia (see Equation 1).

DBA iteratively refines a (randomly initialized) series to minimize its
quadratic distance (DTW) with respect to the series to be averaged. At each
iteration, the algorithm performs “vertical” centroids of the aligned elements
of the time series. The definition of a “vertical” centroid corresponds to
averaging series real values. Our approach enables the adaptation of timed
sequences with event type distributions.

###  4.1 Average Timed Sequences

Let $\mathcal{S}={s_{1},\dots,s_{m}}$ be a set of timed sequences, Algorithm 2
details the steps for generating an average timed sequence inspired by the
principle of vertical averages of DBA [14]. At line 1, the current average
time sequence, $\mathcal{C}^{(0)}$ is initialized with a timed sequence chosen
from $\mathcal{S}$. Then, the distance from $\mathcal{C}^{(0)}$ to each of the
other sequences is calculated to retrieve the optimal alignment $M^{*}$
represented by the indices $i_{r}$ and $i_{s}$ of matching elements in each
timed sequence. Lines 8 to 10 collect, on the one hand, the list of paired
events in the sequences of $\mathcal{S}$ per event in $\mathcal{C}^{(0)}$ and,
on the other hand, their dates per event in $\mathcal{C}^{(0)}$. These lists
are then averaged lines 11-12. The average of events is constructed as the
histogram of event distributions. Finally, line 13 creates a new average
sequence which can be refined in a subsequent iteration until convergence. In
our implementation, we assume convergence after a fixed number of iterations.

Similar to DBA, events from probabilistic sequences are represented in a
vector space. Thus, the new average sequence is always closer (according to
drop-DTW) to the sequences it averages, and the algorithm converges well (see
Section 4.2). Additionally, this algorithm preserves the drop-DTW’s ability to
allow temporal shifts and handle missing or aberrant events.

The choice of the timed sequence in line 1 is crucial for the result. We
choose to randomly select a sequence of the longest length. We will see that
our algorithm decreases the length of the current average timed sequence, but
it can not augment it.

Data : $\mathcal{S}=\\{s_{1},s_{2},...,s_{m}\\}$: list of timed sequences

1 $\mathcal{C}^{(0)}\leftarrow alea(\mathcal{S})$;

2 $k\leftarrow 0$;

3 while  __ do

$\mathcal{E}\leftarrow[\emptyset,\dots,\emptyset]$ ;

// sets of length $|r|$

$\mathcal{T}\leftarrow[\emptyset,\dots,\emptyset]$ ;

// idem

4 for  _$i=1,\dots,m$_ do

$\mu,\pi^{k}_{i}\leftarrow DropDTW(\mathcal{C}^{(k)},s_{i})$ ;

// calculate the drop-DTW and the best alignment between $\mathcal{C}^{(k)}$
and $s_{i}$

5 for  _$(i_{r},i_{s})\in\pi^{k}_{i}$_ do

$\mathcal{E}_{i_{r}}\leftarrow\mathcal{E}_{i_{r}}\cup\\{s_{i}[i_{s}].\bm{d}\\}$
;

// add event from $s_{i}$ aligned at $i_{r}$

$\mathcal{T}_{i_{r}}\leftarrow\mathcal{T}_{i_{r}}\cup\\{s_{i}[i_{s}].t\\}$ ;

// add date from $s_{i}$ aligned at $i_{r}$

6

$\tau\leftarrow\left\\{\tau_{j}\,\mid\,\forall
j\in[|\mathcal{C}^{(k-1)}|],\;\tau_{j}=\frac{1}{|\mathcal{T}_{j}|}\sum_{t\in\mathcal{T}_{j}}t,\,|\mathcal{T}_{j}|\neq
0\right\\}$ ;

// calculate vertical averages

7 $\mathcal{H}\leftarrow\left\\{\bm{h}_{j}\,\mid\,\forall
j\in[|\mathcal{C}^{(k-1)}|],\;\bm{h}_{j}=\frac{1}{|\mathcal{E}_{j}|}\sum_{\bm{d}\in\mathcal{E}_{j}}\bm{d},\,|\mathcal{E}_{j}|\neq
0\right\\}$ ;

$\mathcal{C}^{(k+1)}\leftarrow(\mathcal{H},\tau)$;

// refined average timed sequence

8 $k\leftarrow k+1$;

9

10return _$\mathcal{C}^{(k)}$_ ;

Algorithm 2 TSR: Generation of an Average Timed Sequence

Figure 2 illustrates one iteration of average calculation for the three
sequences of Example 1. It assumes that we already have a current average
sequence $s_{r}$. The figure depicts the alignment by drop-DTW of $s_{r}$ with
each sequence. $s^{\prime}_{r}$ is the resulting updated average sequence. As
the $A$ event has been paired with both $A$ and $C$, the first probabilistic
event of $s^{\prime}_{r}$ combines the two types of events (proportionally).
As the first $B$ event of $s_{r}$ is paired with no event, it is removed from
$s^{\prime}_{r}$. The other events are purely paired with similar event types,
then the probabilistic event is a pure $C$ or $B$ in $s^{\prime}_{r}$ but
their time positions have been updated.

Figure 2: Illustration of one averaging step. At the top, the alignment of the
current average $s_{r}$ is computed with each sequence. Then, the “vertical”
average of the set of events paired to one $s_{r}$ event yields an average
probabilistic event in $s^{\prime}_{r}$ (at the bottom).

###  4.2 Convergence of the Algorithm

The algorithm 2 refines iteratively an average sequence for a set of
sequences, $\left(s_{i}\right)_{i\in[m]}$. In this algorithm,
$\mathcal{C}^{(k)}$ denotes the current average sequence at the $k$-th
iteration. At step $k$, the inertia of the average sequence is defined as:

$I^{(k)}=\frac{1}{m}\sum_{i\in[m]}DropDTW\left(\mathcal{C}^{(k)},s_{i}\right)$
(1)

The inertia represents the average error incurred by using $\mathcal{C}^{(k)}$
as the representative of the set of sequences $\left(s_{i}\right)_{i\in[m]}$.
In the following, we prove that $I^{(k+1)}$ is a decreasing suite, i.e.,
$I^{(k+1)}\leq I^{(k)}$, for all $k\in\mathbb{N}$.

The $i$th term of the inertia is

$DropDTW\left(\mathcal{C}^{(k)},s_{i}\right)=\min_{\pi}\sum_{(x,y)\in\pi}d\left(\mathcal{C}^{(k)}[x],s_{i}[y]\right)+\delta.\left(P_{r}(\pi)+P_{c}(\pi)\right)$

where $d$ is the metric between timed events (see Section 3.2) and $\pi$ plays
the same role as the matrix $M$ introduced in section 3.1. It represents the
alignment between two sequences. In the remainder, we denote by
$\pi^{(k)}_{i}$ the optimal alignment of sequence $s_{i}$ with sequence
$\mathcal{C}^{(k)}$, i.e., the alignment that minimizes the sum of costs.

Let us now introduce
$q_{i}^{(k+1)}=\sum_{(x,y)\in\pi^{k}_{i}}d\left(\mathcal{C}^{(k+1)}[x],s_{i}[y]\right)+\delta.\left(P_{r}(\pi^{k}_{i})+P_{c}(\pi^{k}_{i})\right)$,
for all $i\in[m]$ and $k\in\mathbb{N}^{*}$. Providing $\pi^{k+1}_{i}$ is the
optimal alignment, we have:

$\sum_{(x,y)\in\pi^{k+1}_{i}}d\left(\mathcal{C}^{(k+1)}[x],s_{i}[y]\right)+\delta.\left(P_{r}(\pi^{k+1}_{i})+P_{c}(\pi^{k+1}_{i})\right)\leq
q_{i}^{(k+1)}$

It is worth noticing that we retrieve the inertia of the $k+1$-th average
sequence $I^{(k+1)}$ by summing the left-hand terms over $i$, and then:

$I^{(k+1)}\leq\frac{1}{m}\sum_{j\in[m]}q_{i}^{(k+1)}$ (2)

Let us not have a look at the sum over $i$ of the $q_{i}^{(k+1)}$. This sum
can be rewritten as follow:

$\displaystyle\frac{1}{m}\sum_{j\in[m]}q_{j}^{(k+1)}$ $\displaystyle=$
$\displaystyle\frac{1}{m}\sum_{j\in[m]}\sum_{(x,y)\in\pi^{k+1}_{j}}d\left(\mathcal{C}^{(k)}[x],s_{j}[y]\right)+\delta.\left(P_{r}(\pi^{k+1}_{j})+P_{c}(\pi^{k+1}_{j})\right)$
$\displaystyle=$
$\displaystyle\sum_{x}\frac{1}{m}\sum_{j\in[m]}\sum_{\left\\{y\mid(x,y)\in\pi^{k+1}_{j}\right\\}}d\left(\mathcal{C}^{(k)}[x],s_{j}[y]\right)+\delta.\left(P_{r}(\pi^{k+1}_{j})+P_{c}(\pi^{k+1}_{j})\right)$

Let us now have a look at the construction of $\mathcal{C}^{(k+1)}[x]$ in the
algorithm (see lines 11-13). It is the barycenter of all elements $s_{j}[y]$
that have been associated with $\mathcal{C}^{(k)}[x]$ by the alignment
$\pi^{k+1}_{j}$. Thanks to the property of barycenters in an Euclidean space,
we have for any $\bm{v}$:

$\frac{1}{m}\sum_{j\in[m]}\sum_{\left\\{y\mid(x,y)\in\pi^{k}_{j}\right\\}}d(\bm{v},s_{j}[y])\geq\frac{1}{m}\sum_{j\in[m]}\sum_{\left\\{y\mid(x,y)\in\pi^{k}_{j}\right\\}}d(\mathcal{C}^{(k+1)}[x],s_{j}[y])$
(3)

Considering that:

$I^{(k)}=\sum_{x}\frac{1}{m}\sum_{j\in[m]}\sum_{\left\\{y\mid(x,y)\in\pi^{k}_{j}\right\\}}d\left(\mathcal{C}^{(k)}[x],s_{j}[y]\right)+\delta.\left(P_{r}(\pi^{k}_{j})+P_{c}(\pi^{k}_{j})\right)$

and applying inequation 3 with $\bm{v}=\mathcal{C}^{(k)}[x]$, it ensues:

$\displaystyle I^{(k)}$ $\displaystyle\geq$
$\displaystyle\sum_{x}\frac{1}{m}\sum_{j\in[m]}\sum_{\left\\{y\mid(x,y)\in\pi^{k}_{j}\right\\}}d(\mathcal{C}^{(k+1)}[x],s_{j}[y])+\delta.\left(P_{r}(\pi^{k}_{j})+P_{c}(\pi^{k}_{j})\right)$
(4) $\displaystyle\geq$ $\displaystyle\frac{1}{m}\sum_{j\in[m]}q_{j}^{(k+1)}$

In the case $\left\\{y\mid(x,y)\in\pi^{k}_{j}\right\\}=\emptyset$, this means
that for all $s_{i}$, there is an empty row at line $x$ of the alignment with
$\mathcal{C}^{(k+1)}$. For instance, the first $B$ event in Figure 2 is
assigned to none event. Then, this event contribute to the inertia only via a
drop-cost. Removing such an event from $\mathcal{C}^{(k+1)}$ turns out to
remove an empty row for all sequences $s_{j}$, and thus substract one to all
$Pr(\pi^{k+1}_{j})$. Then, it also decreases the true values of
$q_{j}^{(k+1)}$.

Finally, with equations 2 and 4, we conclude that $I^{(k+1)}\leq I^{(k)}$.
This proves that Algorithm 2 actually refines the average timed sequence and
converges. Nevertheless, it does not guarantee to find the global minimum of
inertia. It only ensures to find a local minimum.

###  4.3 Parameters

The proposed method for averaging time sequences has two types of parameters:
the parameters of the drop-DTW metric ($\tau\in\mathbb{R}$, $p_{e}\in[0,1]$,
$p_{t}\in[0,1]$, $\delta\in\mathbb{R}$) and the parameter of the Algorithm 2
(maximum iterations for drop-DTW and iterative averaging).

The proposed metric is mainly parameterized by the ratio $\frac{p_{t}}{p_{e}}$
and the drop-cost $\delta$. We explicitly explain here how to determine these
values relative to the unique threshold value $\tau$. Thus, the expert only
has one parameter to set up: $\tau$, which can be seen as a threshold of
indifference delay. Above this delay, two events cannot be considered as
similar, regardless of their types. An expert can set up $\tau$ based on this
intuition.

The ratio $\frac{p_{t}}{p_{e}}$ can be inferred from the delay $\tau$ by the
following reasoning. For events located within the time interval
$[t-\tau,t+\tau]$, the importance of the event type should be greater than the
importance of time. We must then choose $p_{e}$ to be at least $\tau$ times
greater than $p_{t}$. With lost of generality, the value of $p_{t}$ can be set
at $1$ by default. The drop-cost is also derived from $\tau$, i.e., the delay
from which there is no longer any interest in associating two events. We thus
propose to define the drop-cost as:

$\delta=\frac{1}{\tau^{2}}(\tau+t_{max})^{2}+1$

where $t_{max}$ is the maximum time distance to consider.

##  5 Clustering

The proposed metric between timed sequences and the averaging techniques are
used to adapt two classical clustering techniques: a hierarchical ascendant
clustering (HAC) and a K-means algorithm [10]. These two algorithms require to
predefine the expected number of clusters, $K$.333Model selection based on the
inertia of the clustering may be applied to select the optimal number of
clusters.

HierAsTiSeq denotes the HAC algorithm for timed sequences. The principle of
the algorithm is to iteratively aggregate the two closest timed sequences
according to the drop-DTW. A group of two timed sequences is then averaged
using algorithm 2, and the average timed sequence replaces the two former
ones. The algorithm terminates when it reaches the expected number of
clusters.

The K-means also proceeds iteratively to refine $K$ clusters. K-means
algorithm has two main steps: first the sequences of each cluster are averaged
(algorithm 2). The computed averages are the representatives of their
clusters. Then, new clusters are made by assigning optimally each sequence to
one of the $K$ representative.

##  6 Experiments and Results on Synthetic Data

In this section, we present experiments on synthetic. We conducted several
tests on synthetic data to evaluate the effectiveness of the method under
various scenarios and to study the effect of the parameters.

For this first test, we created two model timed sequences that differ in one
or more of the dimensions we want to study: sequence length, event nature, and
event dates. Therefore, there are $8=2^{3}$ possible variations. A dataset
consists of 15 sequences for each of the 8 possible variations and the base
sequence (thus $135=9\times 15$ sequences in total). The dates of these timed
sequences are obtained and made random by adding uniform noise to the date.
The objective of this experiment is to evaluate the accuracy of identifying
the 9 categories of sequences.

We tested two sets of parameters: first, we chose the ratio
$\frac{p_{t}}{p_{e}}=\frac{1}{9}$, with $\delta=+\infty$. In this
configuration, the cost of an event is more important than the time cost for 3
days around the event. The second set of parameters used was
$\frac{p_{t}}{p_{e}}=\frac{1}{400}$, still with $\delta=+\infty$. This time,
the difference cost between events is more important than the time cost over
the first 20 days.

In the case of the first configuration, we observe that all clusters were
correctly identified (Kappa score of $1$). However, in the case of the second
configuration, we observed a degradation of results with a Kappa score of
$0.62$. This degradation was expected because the maximum time distances
observed between events are around 20. Consequently, the time notion is
practically never taken into account by the metric.

This experiment thus shows that the choice of parameters is crucial but can be
made by taking into account the distance between events of the sequences.

###  6.1 Effect of Drop-Cost

We can then test the effect of the drop-cost on a synthetic dataset, with and
without drop-cost. We take 3 model sequences, 2 of which are similar, except
that one of them has an aberrant value.

  1. 1.

$\langle(D,0.0)(E,3.0)(F,5.0)(D,7.0)(D,10.0)(E,12.0)(F,15.0)\rangle$

  2. 2.

$\langle(E,2.0)(A,4.0)(D,8.0)(C,12.0)\rangle$

  3. 3.

$\langle(D,0.0)(E,3.0)(F,5.0)(D,7.0)(D,10.0)(E,12.0)(F,15.0)(D,1500.0)\rangle$

| Prédit
---|---
Réel | 0 | 15
15 | 15

| Prédit
---|---
Réel | 15 | 0
0 | 30

Table 1: Confusion matrices with supplementary events. On the left:
$p_{t}=\frac{1}{9}$, $p_{e}=1$,$\delta=+\infty$; On the right:
$p_{t}=\frac{1}{9}$, $p_{e}=1$, $\delta=4$

Intuitively, we want to cluster sequences 1 and 3 together. We apply
HierAsTiSeq without drop-cost and obtain the confusion matrix in Table 1 on
the left, with a Kappa score of $-0.5$ (disagreement). We observe that the
algorithm associated the set from sequence 1 with the set from sequence 2.

Adding drop-cost coinciding with the distances between the sequences allows us
to obtain the confusion matrix in Table 1 on the right, with a Kappa score of
$1$ (total agreement). The drop-cost allows neglecting the last event
considered aberrant, and the results match expectations.

| Predicted
---|---
Real | 0 | 15
15 | 15

| Predicted
---|---
Real | 15 | 0
0 | 30

Table 2: Confusion matrices with missing events. On the left:
$p_{t}=\frac{1}{9}$, $p_{e}=1$,$\delta=+\infty$; On the right:
$p_{t}=\frac{1}{9}$, $p_{e}=1$, $\delta=4$

We also conducted a similar experiment with missing values. This time, the
model sequences are as follows:

  1. 1.

$\langle(D,0.0),(E,3.0),(F,5.0),(D,7.0)\rangle$

  2. 2.

$\langle(E,2.0),(A,4.0),(E,7.0),(D,8.0),(D,10.0),(E,12.0)\rangle$

  3. 3.

$\langle(D,0.0),(E,2.0),(E,3.0),(F,5.0),(D,7.0),(F,9.0),(A,13.0)\rangle$

We then notice that all events dated from sequence 1 are present in sequence
3, but they are all different from sequence 2. We would therefore like to
cluster sequences from 1 and 3 together. Keeping the same parameters as
before. Without drop-cost, we obtain the confusion matrix in Table 2 on the
left. This result corresponds to expectations because the sequences from 2 and
3 are then associated as they are closer in time. However, by adding a
$\delta=4$, we obtain the confusion matrix in Table 2 on the right, which is
logical because the events furthest in time are no longer associated.

The tests on synthetic values present the expected results, thus validating
the effectiveness of clustering in the situations described. We now continue
the study of the method with an experiment on real data to evaluate its
ability to identify relevant clusters.

##  7 Application to Real-World Care Pathways

In this section, we apply the proposed method to real-world care pathways and
compare our results with those of Sequence State Analysis (SSA) [6]. SSA,
particularly the TraMineR tool, is commonly used for categorizing care
pathways [17].

This experiment is conducted on the data warehouse of the Greater Paris
Hospital (AP-HP). We identified a population of 3311 patients having a stay
with a lung resection surgery. We collected events that occurred within 90
days around a lung resection surgery as well as the date of death if
applicable. Events further away are clinically less interesting to study. For
clinical reasons, we also decided to keep only the first occurrence of
chemotherapy for each patient. We made the same choice for radiotherapy and
immunotherapy. It is worth noting that deaths are excluded from the the care
pathway. However, we keep this information for the survival analysis of the
clusters.

The parameters chosen are $K=5$ clusters, $\tau=7$, $t_{\text{max}}=16$, and
$\delta=1$. For choosing the number of clusters, different values were tested.
By increasing the number of clusters, it is mainly small-sized clusters that
split. Therefore, there may not necessarily be an interest in increasing this
number. $\tau$ is set to 7 days because, in our dataset, it may be interesting
to associate the same types of events that occurred within a week of each
other. Finally, based on the dataset, we consider that there is no longer any
interest in associating two events more than 16 days apart.

Figure 3: Clusters obtained by the clustering methods, from left to right:
HierAsTiSeq, K-means and TraMineR. Clusters are row-wised. Each cluster is
represented by the histogram of event types over time. Each color corresponds
to a specific type of event (see Figure 4 for color legend). The higher the
bar, the more events there are at the given time across all patients. The
vertical red line indicates the index-date of the resection (0 delay).

###  7.1 Comparison with TraMineR

In this section, we selected randomly 500 patients among the cohort of
patients to compare the three clustering approaches: HierAsTiSeq, K-means with
drop-DTW and TraMineR. Regarding the computing time of HierAsTiSeq and
TraMineR, this comparison is conducted only on a small subset of patients.
TraMineR requires converting timed sequences into time series. As TraMineR is
sensitive to sparse series, we increased the time sampling from 1 to 3 days
for calculations, and we added the empty event corresponding to the absence of
data over a period of 3 days. TraMiner uses the hierarchical clustering (ward
criteria) and on optimal matching metric between state sequences.

Figure 3 corresponds to the 5 clusters for the three methods. Each cluster is
represented as an histogram of events over time. One can first notice that
TraMineR provides clusters that look very similar to each others. It is
difficult to identify specificities for each cluster. Contrary to the other
two methods, based on our the drop-DTW metric, the sizes of the clusters are
more homogeneous. For HierAsTiSeq and K-means, there is one large cluster and
several small clusters that gather patients specific timed events. The two
large clusters (cluster 5 for HierAsTiSeq and cluster 1 for K-means) looks
very similar in distribution of events. HierAsTiSeq identifies two clusters
with specific behaviors: the cluster 1 has no post-operative events. This may
correspond to patients that had the surgery in the hospital but their follow-
up in another medical structure; the cluster 2 are patients with short delays
between pre- or post-operative care events around the surgery. The clusters 2
and 3 identified by K-means also gather patients with a late pre-operative
consultation (wrt the resection). According to the distributions, these
patients have two pre-operative consultations: in cluster 2, the first pre-
operative consultation occurs earlier than in cluster 3; and cluster 3 holds
patients having a late post-operative consultation.

We conclude from this experiment that the use of drop-DTW yields clusters of
patients with more specific timed sequences than TraMineR. For this dataset,
the latter method does not seems to be able to reveal interesting clusters of
patients. We performed a quantitative (confusion matrices) and qualitative
analysis. The analysis of confusions between HierAsTiSeq and TraMineR shows
that there is no correspondence between the clusters. This suggests that the
metrics are semantically different. In this case, the semantics of drop-DTW on
timed sequence seems more meaningful, but it does not mean that it would be
the case for all datasets.

It is worth noting that results of K-means with drop-DTW may significantly
differ across different runs due to the randomness of the cluster seeds.
HierAsTiSeq and TraMineR are more deterministic, but the iterative computation
of drop-DTW or optimal matching may introduce some randomness. Nonetheless,
the hierarchical clustering strategy of HierAsTiSeq and TraMineR is resource
demanding. This prevents to apply them on the entire dataset. This was
possible only with K-means and takes only a few minutes (on 1 core of the AP-
HP servers).

Figure 4: Clusters of care-pathways identified by K-means with drop-DTW on
3211 patients. Each row is a cluster with the representation of the barycenter
on the left and the representation of the events distribution over time on the
right. Each color corresponds to a specific type of event.

###  7.2 Cluster Analysis of K-means with drop-DTW

In this section, we present the results of the K-means clustering with drop-
DTW metric applied on the entire cohort. Figure 4 illustrates the 5 clusters
yielded. Each cluster is depicted with its barycenters and the distribution of
events.

Let us first comment on the clusters based on their event distributions. Three
main clusters have been identified (clusters 1, 2 and 5) and share almost the
same size (about 1000 patients). The two others clusters are very small
(number 3 and 4). These small clusters of patients seem to be outliers.
Similarly to the previous results, the K-means identified cluster of patients
with different delays between the pre-operative consultations and resection.
The cluster 2 corresponds to patients having a small delay between this events
and possibly two pre-operative consultations. These patients possibly have a
late stage of the disease and need a quick surgery. We made this hypothesis
because their survival at 6 months is significantly lower than the two others
clusters. The cluster 5 and cluster 1 differs mainly in the delay between the
surgery and the post-operative consultations. The Cluster 1 also holds
patients without exams the day of the resection.

The same observations can be made based on the analysis of the barycenters. We
remind that the barycenter is a probabilistic timed sequence. Each event is
represented by a distribution of event’s types. On the figure, each bar is an
event located at its time of occurrence. A bar stacks the probabilities per
event’s type. We remark that most of the events have one main type. The event
at time 10 for cluster 1 is a mixture of mainly two events types: post-
operative consultation or post-operative multi-disciplinary consultation
meeting (MCM). This possible confusion may be explained by the imprecision of
the records. MCM are identified by consultation reports. If the report is not
correctly tagged as an MCM, then it is simply a consultation. Then, the
flexibility of the metric makes it somehow robust to possible labeling errors
in the data.
We can also notice that there is a clear alignement between the barycenters
events and the event distribution. This illustrates that the averaging method,
based on drop-DTW provides a meaningful average of a collection of complex
timed sequences. More specifically, comparing the clusters by the decay
between consultations and surgery leads to the same conclusions as before.
Finally, one can remark that there are a lot of repeated events in the
barycenter. Furthermore, the presence of repeated events in the barycenter can
be misleading, as it could either indicate a repetition of events or
alternative delays for a single event due to possible drops.

##  8 Conclusion

We proposed a method for clustering timed sequences inspired by existing
techniques for time series. Thus, the use of timed sequences allows co-
occurring events, irregular time sampling, and sparse events. All these
advantages are very interesting from a clinical point of view for describing
care pathways. We showed that by using drop-DTW, the algorithm became less
sensitive to missing and aberrant values, which is also useful in application
to real-world care pathways. Then, we tested our method on real data from the
OPTISOINS study for our approach compared to TraMineR. Based on the initial
results, drop-DTW seems to be clinically promising for identifying average
pathways.

However, the large number of parameters and the required computation time are
limitations of the method. Regarding the number of parameters, our effort
focused on a guide for choosing parameter values, which the experiments tend
to validate. In perspective, a more in-depth study of the effect of parameters
still needs to be carried out. Additionally, we have to validate the clinical
interest of average pathways by clinicians.

#### Acknowledgments

The authors thanks gratefully Dr. Yen-Lân NGuyen from Department of
Anaesthesia and Intensive Care Medicine, AP-HP, Cochin Academic Hospital, for
her time and expertise she provided during this project. Part of the research
presented in this article is funded by the AP-HP Foundation, as part of the
AI-RACLES Chair, and has received approval from the Scientific and Ethical
Committee of the AP-HP CDW (CSE-20-01-OPTISOINS).

## References

  * [1] Guillaume Bernard, Cyril Faucher, and Karell Bertet.  Towards reconstruction of human trajectories in indoor environments.  In EKAW 2018, volume 2262, pages 37–40. CEUR-WS. org, 2018.
  * [2] Philippe Besnard and Thomas Guyet.  Chronicles: Formalization of a Temporal Model.  Springer Briefs in Computer Science, 2023.
  * [3] Yuanzhe Chen, Panpan Xu, and Liu Ren.  Sequence synopsis: Optimize visual summary of temporal event data.  IEEE transactions on visualization and computer graphics, 24(1):45–55, 2017.
  * [4] Paul Compagnon, Grégoire Lefebvre, Stefan Duffner, and Christophe Garcia.  Sequence metric learning as synchronization of recurrent neural networks.  In International Joint Conference on Neural Networks (IJCNN), pages 1–8, Glasgow, United Kingdom, 2021.
  * [5] Mikita Dvornik, Isma Hadji, Konstantinos G Derpanis, Animesh Garg, and Allan Jepson.  Drop-DTW: Aligning common signal between sequences while dropping outliers.  Advances in Neural Information Processing Systems (NIPS), 34:13782–13793, 2021.
  * [6] Alexis Gabadinho, Gilbert Ritschard, Nicolas S. Müller, and Matthias Studer.  Analyzing and visualizing state sequences in R with TraMineR.  Journal of Statistical Software, 40(4):10, 2011.
  * [7] Damien Garreau, Rémi Lajugie, Sylvain Arlot, and Francis Bach.  Metric learning for temporal sequence alignment.  Advances in neural information processing systems, 27, 2014.
  * [8] Dominique Gay, Romain Guigourès, Marc Boullé, and Fabrice Clérot.  TESS: temporal event sequence summarization.  In Proceedings of the International Conference on Data Science and Advanced Analytics (DSAA), pages 1–10, 2015.
  * [9] E. Audureau H. Sebia, T. Guyet.  SWoTTeD: An extension of tensor decomposition to temporal phenotyping.  Machine Learning, page to appear, 2024.
  * [10] John A Hartigan.  Clustering algorithms.  John Wiley & Sons, Inc., 1975.
  * [11] Christopher Holder, Matthew Middlehurst, and Anthony Bagnall.  A review and evaluation of elastic distance functions for time series clustering.  Knowledge and Information Systems, 66(2):765–809, 2024.
  * [12] Dang Nguyen, Wei Luo, Tu Dinh Nguyen, Svetha Venkatesh, and Dinh Phung.  Sqn2vec: Learning sequence representation via sequential patterns with a gap constraint.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 569–584. Springer, 2018.
  * [13] François Petitjean and Pierre Gançarski.  Summarizing a set of time series by averaging: From steiner sequence to compact multiple alignment.  Theoretical Computer Science, 414(1):76–91, 2012.
  * [14] François Petitjean, Alain Ketterlin, and Pierre Gançarski.  A global averaging method for dynamic time warping, with applications to clustering.  Pattern Recognition, 44(3):678–693, 2011.
  * [15] Kishan Rama, Helena Canhão, Alexandra M Carvalho, and Susana Vinga.  Aliclu-temporal sequence alignment for clustering longitudinal clinical data.  BMC medical informatics and decision making, 19(1):1–11, 2019.
  * [16] Kishan Rama, Helena Canhão, Alexandra M. Carvalho, and Susana Vinga.  AliClu - temporal sequence alignment for clustering longitudinal clinical data.  BMC Medical Informatics and Decision Making, 19(1):289, 2019.
  * [17] Jonathan Roux, Olivier Grimaud, and Emmanuelle Leray.  Use of state sequence analysis for care pathway analysis: The example of multiple sclerosis.  Statistical Methods in Medical Research, 28(6):1651–1663, 2019\.
  * [18] Hiroaki Sakoe and Seibi Chiba.  Dynamic programming algorithm optimization for spoken word recognition.  IEEE transactions on acoustics, speech, and signal processing, 26(1):43–49, 1978.
  * [19] Sven Schneider, Maria Maximova, and Holger Giese.  Probabilistic metric temporal graph logic.  In International Conference on Graph Transformation, pages 58–76. Springer, 2022.
  * [20] Nikolaj Tatti and Jilles Vreeken.  The long and the short of it: summarising event sequences with serial episodes.  In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 462–470, 2012.
  * [21] Romain Tavenard, Johann Faouzi, Gilles Vandewiele, Felix Divo, Guillaume Androz, Chester Holtz, Marie Payne, Roman Yurchak, Marc Rußwurm, Kushal Kolar, et al.  Tslearn, a machine learning toolkit for time series data.  The Journal of Machine Learning Research, 21(1):4686–4691, 2020\.
  * [22] Xiufan Yu, Karthikeyan Shanmugam, Debarun Bhattacharjya, Tian Gao, Dharmashankar Subramanian, and Lingzhou Xue.  Hawkesian graphical event models.  In Manfred Jaeger and Thomas Dyhre Nielsen, editors, Proceedings of the 10th International Conference on Probabilistic Graphical Models, volume 138 of Proceedings of Machine Learning Research, pages 569–580. PMLR, 2020.


# Junction-free microwave two-mode radiation from a kinetic inductance
nanowire

Yufeng Wu  Mingrui Xu  Hong X. Tang<EMAIL_ADDRESS>Department of
Electrical Engineering, Yale University, New Haven, Connecticut 06520, USA

###### Abstract

Parametric down-conversion is a widely exploited technique in optics to
produce entangled states of photons for quantum information processing and
quantum sensing. In the microwave domain, devices based on Josephson
junctions, such as Josephson parametric amplifiers (JPAs) and voltage-biased
Josephson junctions, have been successfully utilized to generate such states.
However, their high susceptibility to magnetic fields has posed challenges in
many applications. Here we demonstrate the generation of two-mode squeezed
states via four-wave-mixing in a superconducting nanowire resonator patterned
from NbN. The NbN nanowire exhibits a strong Kerr nonlinearity, resulting in
the emission of a signal-idler pair with a cross-correlation of
$g^{(2)}(0)=11.9$. Owing to the magnetic resilience and high critical
temperature ($T_{c}$) of NbN, our microwave parametric sources based on
kinetic inductance promise an expanded range of potential applications.

††preprint: APS/123-QED

##  I Introduction

Nonclassical states are quantum states characterized by photon statistics that
do not follow classical descriptions. Among these states, two-mode non-
classical states are particularly intriguing due to their strong quantum
correlation. Two-mode squeezed states [1], in particular, play a vital role in
a wide range of applications, such as quantum teleportation [2, 3], quantum
sensing [4, 5, 6, 7, 8, 9, 10] and continuous-variable quantum computing [11,
12]. In the region where the average photon number is very small, such states
can be approximated as photon-pair states, which is an indispensable component
in Bell inequality testing [13, 14, 15], photonic quantum computing [16] and
global quantum networks [17, 18, 19].

The generation of such non-classical two-mode radiation in the microwave
domain has predominantly relied on Josephson Parametric Amplifiers (JPAs) [1,
20, 21, 22]. These devices utilize parametric down-conversion to transform
pump photons into pairs of signal and idler photons exhibiting strong
correlation. Recently, this process has also been demonstrated in DC-biased
Josephson junctions through inelastic Cooper pair tunneling [23, 24],
traveling wave parametric amplifiers (TWPAs) [25, 26], and access higher
energy level of a superconducting transmon qubit [27]. However, without
exception, all these devices are based on Josephson junctions, which are
sensitive to external flux. Therefore, they are not compatible with quantum
systems that require strong magnetic fields, for example, two-mode squeezed
state enhanced dark-matter axion search [7, 28]. The traditional approach to
mitigating the influence of the magnetic field is to move the parametric
amplifier away from the magnet and protect it with layers of shield [6]. In
some scenarios, flux feedback is required to cancel the excessive magnetic
field [29]. On the other hand, higher-order nonlinearities in JJ-based devices
have been identified as a limiting factor for the amplification and squeezing
[30]. These limitations have prompted researchers to explore alternative
approaches to overcome these challenges and expand the potential applications
of these devices in various quantum systems.

Niobium nitride (NbN) and niobium titanium nitride (NbTiN) have been gathering
attention due to their high kinetic inductance nonlinearity properties.
Utilizing a nanowire structure, these materials can function as parametric
amplifiers with quantum-limited added noise [31, 32]. Furthermore, NbN
exhibits a high critical temperature (around 10 K), facilitating kinetic
inductance traveling-wave parametric amplifiers at 4 Kelvin with an added
noise temperature of 6.5 K [33]. Additionally, it demonstrates robust in-plane
magnetic field resilience, allowing NbN based resonators to operate at
magnetic field up to 6 T [34], and quantum amplifiers to operate up to 2 T
[32, 35].

In this work, we observed non-classical two-mode radiation from a kinetic
inductance nanowire based on NbN. We send two symmetric detuned pumps to the
device to drive the parametric down-conversion, and demodulate the two
itinerant modes that are evenly detuned from the resonance frequency. The non-
classicality is characterized by a violation of the classical bound of the
second-order cross-correlation. Our investigation reveals an second-order
cross-correlation reaching $11.91\pm 0.68$ with MHz photon-pair generation
rate. These results solidify the benefits of kinetic inductance devices as
potential replacements for JJ-based parametric amplifiers.

##  II Hamiltonian and Photon Correlations

Figure 1: Device and illustration of the two-mode generation. (a) Optical
image of the NKPA. The brighter area represents NbN, while the darker area
signifies the silicon substrate. The inset SEM image displays a zoomed-in view
of the nanowire. (b) Pumping schematic of the NKPA device. The two arrows
illustrate the symmetric pump, and the two shaded regions represent the signal
and idler modes to be demodulated. (c) Artistic illustration of the two-mode
generation from an NbN nanowire. The inserted is the equivalent lumped element
circuit diagram of the NKPA.

The nanowire kinetic-inductance parametric amplifier (NKPA) device is
fabricated from NbN with a thickness of 4 nm deposited on silicon substrate
via atomic-layer deposition (ALD) [32]. The device is designed as an LC
resonator that incorporates highly nonlinear kinetic inductance from the
nanowire, shown in Fig.1 (a). Under a symmetric two-tone drive, as shown in
Fig. 1 (b), the system Hamiltonian in the interacting picture becomes
(ignoring the self-Kerr and cross-Kerr component) [32]

$\displaystyle H_{\text{int}}=\hbar\left(\epsilon
a^{\dagger}b^{\dagger}+\epsilon^{*}ab\right),$ (1)

where $\epsilon=KA_{1}A_{2}e^{i(\phi_{1}+\phi_{2})}$ is the parametric down-
conversion coefficient, and $A_{1},A_{2}$, $\phi_{1},\phi_{2}$ are the two-
pump amplitudes and phases, $K=2\pi\times 110$ kHz is the single-photon Kerr
nonlinearity. Detailed derivation can be found in Appendix A. This represents
a two-mode squeezing Hamiltonian. Non-classical correlations emerge for two
modes symmetrically positioned around the resonance frequency. The relative
position and total momentum between the two modes are squeezed below the
standard vacuum limit [1]. However, due to the added noise from the amplifier
in the subsequent circuit, this can only be observed when operating the
amplifier in the high-gain region, and may sometimes necessitate an additional
quantum-limited preamplifier [21]. In this work, we operate the amplifier in
the low-gain region, the system serves as a photon-pair source, as illustrated
in Fig. 1 (c). The non-classical nature of the two-mode squeezing state is
characterized by the second-order correlation

$g^{(2)}_{\alpha\beta}(\tau)=\frac{\langle\beta^{\dagger}(0)\alpha^{\dagger}(\tau)\alpha(\tau)\beta(0)\rangle)}{\langle\alpha^{\dagger}(\tau)\alpha(\tau)\rangle\langle\beta^{\dagger}(0)\beta(0)\rangle},$
(2)

with $\alpha,\beta\in\\{a,b\\}$. In analogy with optics, the numerator
represents the coincidence count rate and the denominator represents the
product of the count rate of the two respective modes. It is important to note
that while these values depend on the system loss, their quotient
$g^{(2)}(\tau)$ remains independent of system loss [36]. For a classical
state, zero-delay cross-second order correlation is smaller than the average
of the two auto-second order correlation

$g^{(2)}_{ab}(0)\leq\frac{g^{(2)}_{aa}(0)+g^{(2)}_{bb}(0)}{2}.$ (3)

A two-mode squeezed state can be expressed as [37]

$|\Phi_{TM}\rangle=\frac{1}{\cosh{r}}\sum_{n}(\tanh{r})^{n}|n,n\rangle,$ (4)

where $r=\epsilon\tau$ is the squeezing parameter, $\tau$ is the cavity decay
time. Under a weak drive, where $r$ is small, the two-mode squeezed state is
reduced to a mix of vacuum and photon-pair states. Therefore, a higher-than-2
$g^{2}_{ab}(0)$ can be observed and thus violate Eq. (3), inferring the non-
classical characteristics.

##  III Experimental Setup

Figure 2: Schematic of the experimental setup. (a) Cryogenic electrical
circuit connection schematic. The two pumps are combined and directed toward
the NKPA device located in the MXC chamber. (b) The room-temperature
measurement setup. The output signal is split into two demodulation channels
with a $90^{\circ}$ hybrid. The intermediate-frequency filter is centered at
$74.6$ MHz with a bandwidth of $3$ MHz.

The cryogenic circuit connection schematic is shown in Fig. 2 (a). The NKPA
device was placed in the mixing chamber of a dilution refrigerator and cooled
to a temperature of 10 mK. To send the two pump tones to the device, they were
symmetrically detuned from the resonance frequency of $7.359$ GHz by $57$ MHz
and subsequently combined. In order to prevent the impact of room-temperature
thermal noise and pump phase noise, the pump tone is attenuated 60 dB to the
device. The output signal is amplified by a high electron mobility transistor
(HEMT) amplifier and room-temperature (RT) amplifier.

Figure 2 (b) depicts the room-temperature measurement setup. The output signal
is directed towards a $90^{\circ}$ hybrid coupler, where one of the input
ports is terminated by a 50 $\Omega$ resistor, dividing the signal into two
channels. The microwaves are down-converted to $75$ MHz using a mixer, with
local oscillator frequencies
$\omega_{\mathrm{LO1,2}}=f_{c}\pm(75\,\mathrm{MHz}+\delta$), where $\delta=2$
MHz represents the detuning of the signal/idler tone from the resonance
center. The down-converted signal and idler tone are processed through
bandpass filters with a center frequency of $74.60$ MHz and a bandwidth of
$3.88$ MHz. Following intermediate frequency amplification, the signal and
idler tones are transmitted to two detection channels of the digitizer
operating at a sampling rate of 100 MSp/s. This aliases both the signal and
idler frequency to an effective frequency of $25$ MHz, while the signal-to-
noise ratio (SNR) is not degraded due to the narrow-bandwidth bandpass filter.
The digitizer converts the analog signal to digital voltage with an 8-bit
precision and sends it to a graphics processing unit (GPU) for real-time data
streaming and processing. To further reduce the spectral overlap between the
signal and idler and increase the mode selectivity, a digital filter with a
bandwidth of 0.86 MHz is applied to the data in both channels, after which the
second-order auto- and cross-correlation is calculated and the result is
transmitted to the CPU.

Notably, the high noise temperature of the HEMT amplifier causes the
correlation to be obscured by the amplified noise, making it challenging to
extract directly. To overcome this issue, we employed the post-processing
method of “ON/OFF” measurement, which involves taking measurements under two
different conditions - the “ON” and “OFF” states [38, 23, 24]. Specifically,
the “ON” measurement, represents the statistics of both the signal and noise,
while the “OFF” measurement, with one of the pump shift 80 MHz away to destroy
the frequency-matching condition, only characterizes the statistics of the
noise. By combining both sets of measurements, we could extract the statistics
corresponding to the signal alone. The detailed discussion can be found in
Appendix C.

##  IV Result

Figure 3: (color online) (a) Parametric amplifier gain (top panel) and zero-
delayed second-order correlations $g^{(2)}_{aa}(0)$ (red circles),
$g^{(2)}_{bb}(0)$ (blue circles), and $g^{(2)}_{ab}(0)$ (black triangles) with
different pump power (bottom panel). The inserted figure in the bottom panel
is the time-delayed correlation corresponding to the pump power of $-102.36$
dBm. The grey dashed line represents the theoretical value of the auto-
correlation, and the black dashed line is the predicted cross-correlation from
the parametric amplifier gain. (b) Time-delayed second-order cross-
correlations at different pump powers. The dashed lines are the theoretical
fitting of the cross-correlation profile.

The second-order auto- and cross-correlation are calculated simultaneously in
real-time. The results for each drive power are averaged over 500
computational buffers, with $2^{21}$ samples per buffer. To estimate the error
bars, the 500 buffers are divided into 10 segments of equal numbers, and the
standard deviation is computed from the results of these segments. The second-
order cross-correlation result is shown in Fig. 3 (a). By sweeping the pump
power from $-100.82$ dBm to $-102.79$ dBm, we observe an increase on
$g^{(2)}_{ab}(0)$ from $2.21\pm 0.01$ to $11.91\pm 0.68$. This shows a
violation of the classical relation shown in Eq. (3) and thus confirms the
generation and non-classicality of the photon pairs. The increasing
$g^{(2)}_{ab}(0)$ shows higher photon-pair purity with reduced squeezing
parameter in Eq. (4). The inserted figure is the time-delayed second-order
correlation corresponding to $P_{\mathrm{drive}}=-102.36$ dBm. The auto-
correlation remains consistent at around 2 over the power sweep. With a lower
drive ($P_{\mathrm{drive}}<-102.79$ dBm), it is hard to resolve auto-
correlation with a low error bar due to strong thermal-noise correlation and
limitations of the measurement setup (See Appendix. C for detail). The cross
correlation $g^{(2)}_{ab}$, on the other hand, is less susceptible to the
background noise as the noise from the two modes is not correlated.

The parametric amplifier gain $G_{\mathrm{NKPA}}$ is connected to the mean
photon number in each itinerant mode as $\langle a^{\dagger}a\rangle=\langle
b^{\dagger}b\rangle=G_{\mathrm{NKPA}}-1$. The zero-delayed second-order cross-
correlating can be expressed as

$g^{(2)}_{ab}(0)=2+\frac{1}{2\langle
a^{\dagger}(t)a(t)\rangle}=2+\frac{\eta}{2(G_{\mathrm{NKPA}}-1)}.$ (5)

In this context, we have introduced a unitless parameter $\eta$, which ranges
from 0 to 1, to compensate for the imperfections observed in the resulting
two-mode states. These imperfections include non-zero intracavity photon
occupancy and imperfect symmetry between the two modes. In an ideal situation,
$\eta=1$. The definition of $\eta$ is similar to the noise-reduction-factor
(NRF) defined in the reference [23]. We fit $g^{(2)}_{ab}(0)$ with
corresponding gain, as shown in Fig. 3 (a), and get $\eta=0.72$.

In Fig. 3 (b), we show the time-delayed cross-correlation $g^{(2)}_{ab}(\tau)$
with pump power at $-102.79$ dBm, $-102.36$ dBm, and $-101.95$ dBm. The dashed
line shows a fitting of the decay equation
$g^{(2)}_{ab}(\tau)=1+\gamma_{c}/(2R)\cdot e^{-\gamma_{c}|\tau|}$, where
$\gamma_{c}$ is the characteristic decay rate and $R$ is the detected photon-
pair rate [39]. At a pump power of $-102.79$ dBm, the fitted $\gamma_{c}$ is
$5.90$ MHz, and $R$ is $0.21$ MHz. Employing a similar approach, we observe
the photon-pair rates at pump powers of $-102.36$ dBm and $-101.95$ dBm to be
$0.52$ MHz and $0.99$ MHz, respectively.

Figure 4: The illustration compares $g^{(2)}_{\mathrm{ab}}(\tau)$ under
different measurement bandwidths as determined by the physical filter
(represented in brown) and the digital filter (depicted in cyan). The inserts
showcase the filter response function within the frequency domain. The left
section visualizes the physical bandpass filter, which has a bandwidth of 3.88
MHz, while the right section portrays the digital bandpass filter, with a
bandwidth of 0.86 MHz.

The measurement bandwidth, which is set by the physical and digital filters
used to improve the SNR, constrains the detectable rates of photon-pair
occurrences. To demonstrate this, we adjust the measurement bandwidth by
either incorporating or excluding the narrower digital bandpass filter, which
is executed on the GPU. In both scenarios, the physical bandpass filter is
kept in place, and therefore, the measurement bandwidth is set to be 3.88 MHz
without the digital filter and 0.88 MHz with the digital filter. The results
of the second-order cross-correlation are depicted in Fig. 4. The fitting for
both cases yields a very consistent ratio of $\gamma_{c}/(2R)$, with
respective values of $3.94$ and $3.99$. However, the characteristic decay
rate, $\gamma_{c}$, varies with $20.94$ MHz and $7.026$ MHz, respectively,
indicating photon-pair rates of $2.66$ MHz and $0.88$ MHz. These ratios
closely correspond to their relative measurement bandwidths.

The two modes can be demodulated at different frequencies, contingent upon the
device’s bandwidth. The results presented in Fig. 3 are based on a $2$ MHz
detuning. As depicted in Fig. 5, we increment the detuning from 2 MHz to 14
MHz with a fixed pump power of -101.97 dBm and -100.82 dBm, and subsequently
measure the cross-correlation between the two modes. The zero-delayed
correlation augments due to the reduction in the squeezing parameter.
Remarkably, non-classicality is preserved up to 14 MHz, signifying that two
modes separated by 28 MHz still retain non-classical behavior. The bandwidth
of the source is constrained by the gain-bandwidth product, which in our case
is 59 MHz.

In conclusion, we have demonstrated a novel kinetic-inductance device capable
of generating non-classical microwave photon pairs. In comparison to Josephson
junction (JJ)-based devices, this innovative design offers easier fabrication
and eliminates the need for magnetic field shielding when installed within a
cryogenic environment. Furthermore, the device’s high transition temperature
facilitates radiative cooling [40, 41]. As a result, this technology promises
significant reductions in both space requirements and heat budget when
installed in a dilution refrigerator.

Figure 5: Zero-time delayed second-order cross-correlation corresponding to
different demodulation detuning for pump power at $-101.87$ dBm and $-100.82$
dBm. The inserted illustrate the change of demodulation detuning for the two
modes.

###### Acknowledgements.

The authors would like to thank Dr. Yong Sun, Dr. Lauren McCabe, Mr. Kelly
Woods, and Dr. Michael Rooks for assistance in device fabrication. We
acknowledge funding support from DOE Office of Science, National Quantum
Information Science Research Centers, Co-design Center for Quantum Advantage
(C2QA), Contract No. DE-SC0012704. We acknowledge partial funding support from
the Office of Naval Research on the development of nitride-based
superconductors (under Grant No. N00014-20-1-2126).

##  Appendix A Hamiltonian derivation

The Hamiltonian of the NKPA can be expressed in the interaction picture as

$\begin{split}\frac{H_{\mathrm{int}}}{\hbar}&=-\frac{1}{4\hbar}\frac{L_{k0}}{I^{*2}}I_{ZPF}^{4}(a+a^{\dagger})^{4}\\\
&\approx-\frac{6}{4\hbar}\frac{L_{k0}}{I^{*2}}I_{ZPF}^{4}(2a^{\dagger}a+a^{\dagger
2}a^{2}),\end{split}$ (6)

where $a$ is the bosonic operator, $L_{k0}$ is the kinetic inductance of the
nanowire, $I_{ZPF}$ is the zero-point current of the resonator, and $I^{*}$ is
the characteristic current of the nanowire. The approximation is applied by
dropping the fast-rotating term and constant [31]. We consider two symmetrical
detuned pumps as classical drive, and assume the two correlated modes as
$\delta a$ and $\delta b$. The NKPA intracavity field $a$ can be expressed in
the rotating frame of the resonance frequency as

$a=A_{1}e^{i(\Delta t+\phi_{1})}+A_{2}e^{i(-\Delta t+\phi_{2})}+\delta
ae^{i\delta t}+\delta be^{-i\delta t},$ (7)

where $A_{1,2}$, $\phi_{1,2}$ are the pump amplitudes and phases. Plugging Eq.
(7) into Eq. (6), we derive the following Hamiltonian after rotating-wave
approximation

$\begin{split}\frac{H_{\mathrm{int}}}{\hbar}=&-\frac{K}{4}\left(\delta
a^{2}\delta a^{\dagger 2}+\delta b^{2}\delta b^{\dagger 2}\right)\\\ &-K\delta
a\delta a^{\dagger}\delta b\delta b^{\dagger}\\\
&+K\left(A_{1}^{2}+A_{2}^{2}\right)\left(\delta a^{\dagger}\delta a+\delta
b^{\dagger}\delta b\right)\\\ &+\left(\epsilon\delta a^{\dagger}\delta
b^{\dagger}+\epsilon^{*}\delta a\delta b\right),\end{split}$ (8)

where the vacuum Kerr nonlinearity $K$ and interaction strength $\epsilon$ are
defined as

$\begin{split}K&=\frac{6}{\hbar}\frac{L_{k0}}{I^{*2}}I_{ZPF}^{4},\\\
\epsilon&=KA_{1}A_{2}e^{i(\phi_{1}+\phi_{2})}.\end{split}$ (9)

The first term in Eq. (8) is due to the self-Kerr effect and is negligible
when $\delta a$, $\delta b$ are weak signal. The second and third term refere
to the cross-Kerr effect, and the last term represents non-degenerate
parametric down-conversion.

To relate the pump power $P$ and the amplifier gain $G_{\mathrm{NKPA}}$, we
consider the pump detuning $\Delta$, we can write the interaction strength as
$\epsilon=KP\kappa_{\mathrm{ex}}^{2}/(\hbar\omega_{s}(\Delta^{2}+\kappa^{2}))$,
where $\kappa_{\mathrm{ex}}$ is the external loss rate, and
$\kappa=\kappa_{\mathrm{ex}}+\kappa_{\mathrm{in}}$ is the total loss rate. The
squeezing coefficient $r=\epsilon\tau$, where $\tau=1/\kappa$ is the cavity
decay time. Therefore, we can derive

$G_{\mathrm{NKPA}}=(\cosh\left[K\frac{P}{\hbar\omega_{s}}\frac{\kappa_{\mathrm{ex}}^{2}}{\kappa(\Delta^{2}+\kappa^{2})}\right])^{2}.$
(10)

##  Appendix B Cryogenic setup and control instrument

Figure 6: Schematic representation of the cryogenic wiring along with the
room-temperature control and measurement arrangement. A Quinstar double
junction circulator (QCY-G0400802) is employed to offer an isolation of 40 dB
from the output line. The elements labeled as “att1” and “att2” are
attenuators, utilized to prevent power saturation in the subsequent mixers.

The detailed schematics of the cryogenic wiring along with the room-
temperature measurement setup are depicted in Fig. 6. The two pump signals are
generated by a Zurich Instrument SHFSG signal generator, combined, and
attenuated by a total of 60 dB before being fed into the NKPA device. The
emitted two-mode radiation, in conjunction with the residual pump, is
reflected and channeled to the output line through a single junction
circulator (QCY-G0400802). A double junction circulator (QCY-G0400802) is
employed to further reject any reflected tone from the output line. The
microwave radiation output is then amplified by a HEMT (Low Noise Factory,
LNF-LNC$\\_$16B) situated at the 4K plate and transferred to the room-
temperature setup.

Within the room-temperature configuration, the signal is initially amplified
using a low noise amplifier (SLNA 180-38-25-SMA), followed by a high power
amplifier (Mini-circuit ZVE-3W-183+). Subsequently, the signal is split via a
$90^{\circ}$ hybrid (RF-Lambda RFHB02G18GVT). It should be noted that up until
this point, the output microwave still retains the components of the two pump
tones. As such, an attenuator is connected to each line to reduce the overall
power, thus preventing the mixer (Mini-Circuit ZX05-153LH-S+, 1 dB compression
power at 2 dBm) from reaching saturation. The local oscillators for the two
mixers are sourced from the SHFSG and are phase coherent with the pump tones.
The signal is down-converted to 75 MHz, which is well-separated from the down-
converted pump tone, enabling efficient filtering by a narrow bandpass filter
(Texscan Microwave Products 6BL75.1/3.0-EE). The filtered signal is further
amplified (using SR400) and dispatched to the digitizer (AlazarTech ATS9870).

The 8-bit digitizer operates at 100 MSp/s with a $\pm 40$ mV input range. To
commence data acquisition, the SHFSG sends TTL triggers to the digitizer at a
frequency of 5 Hz. Each received trigger initiates the acquisition of a record
consisting of $2^{21}$ data points. Two such records constitute a data buffer
sent to the GPU for processing. These two records correspond to two triggers,
where the first trigger has the SHFSG sourcing the two-pump tone, and the
second trigger shifts one of the pump tones by 80 MHz, negating gain on the
NKPA. These records correspond to the “ON” and “OFF” data states. This data
buffer is subsequently transferred to a GPU (NVIDIA GeForce RTX 2080 Ti) for
correlation calculations.

##  Appendix C Data Processing

The signal is demodulated to $75$ MHz before being sent to the digitizer,
which operates at a sampling rate of $100$ MSp/s. As a result, the signal
frequency is aliased to $25$ MHz. It is important to note that this aliasing
process does not decrease the signal-to-noise ratio since the narrow bandpass
filter has already filtered out the noise at frequencies higher than $100$
MHz. Therefore, the signal-to-noise ratio is preserved in the digitized signal
at $25$ MHz.

The second-order cross-correlation is calculated with Eq. (2)

$\begin{split}g^{(2)}_{ab}(\tau)&=\frac{\langle
b^{\dagger}(0)a^{\dagger}(\tau)a(\tau)b(0)\rangle)}{\langle
a^{\dagger}(\tau)a(\tau)\rangle\langle b^{\dagger}(0)b(0)\rangle}\\\
&=\frac{\langle n_{a}^{s}(\tau)n_{b}^{s}(0)\rangle}{\langle
n_{a}^{s}(0)\rangle\langle n_{b}^{s}(0)\rangle}\\\ &=\frac{\langle
n_{a}(\tau)n_{b}(0)\rangle-\langle n_{a}(\tau)\rangle\langle
n_{b}^{n}(0)\rangle}{\langle n_{a}^{s}(0)\rangle\langle
n_{b}^{s}(0)\rangle}\\\ &+\frac{-\langle n_{a}^{n}(\tau)\rangle\langle
n_{b}(0)\rangle+\langle n_{a}^{n}(\tau)\rangle\langle
n_{b}^{n}(0)\rangle}{\langle n_{a}^{s}(0)\rangle\langle
n_{b}^{s}(0)\rangle},\end{split}$ (11)

where $n^{n}(\tau)$ and $n(\tau)$ represent the output field first-order
autocorrelation of the “OFF” and “ON” process, and
$n^{s}(\tau)=n(\tau)-n^{n}(\tau)$ is first-order autocorrelation of the output
signal field. It is important to note that we assume no correlation between
the noise terms in the two modes, enabling us to express the mean of the
product of two variables as the product of their means.

To calculate the second-order auto-correlation, it is no longer valid to
implement the calculations similar to Eq. 11 due to the emergence of non-
commuting terms. To illustrate, we write the field operator as
$S_{a}(t)=\sqrt{G_{\mathrm{line}}}(a(t)+h_{a}^{\dagger}(t))$, where
$G_{\mathrm{line}}$ is the gain of the output line, $\langle
S_{a}^{\dagger}(t)S_{a}(0)\rangle=n_{a}(t)$, $G_{\mathrm{line}}\langle
a^{\dagger}(t)a(0)\rangle=n^{s}_{a}(t)$, and $G_{\mathrm{line}}\langle
h_{a}^{\dagger}(t)h_{a}(0)\rangle=n^{n}_{a}(t)$. The noise term
$h_{a}^{\dagger}(t))$ is determined by the HEMT. Therefore, the first order
correlation

$\begin{split}n_{a}^{s}(t)&=G_{\mathrm{line}}\langle
a^{\dagger}(t)a(0)\rangle\\\ &=\langle S_{a}^{\dagger}(t)S_{a}(0)\rangle-
G_{\mathrm{line}}\langle h^{\dagger}_{a}(t)h_{a}(0)\rangle-\delta(t)\\\
&=n_{a}(t)-n_{a}^{n}(t)-\delta(t),\end{split}$ (12)

where $\delta_{a}$ comes from the commutation relation
$\left[h_{a}(t),h^{\dagger}_{a}(0)\right]=\delta(t)$. The $\delta$ function
can be eliminated by the limited measurement bandwidth [42, 38]. The second-
order correlation can be calculated in a similar fashion by rearranging the
field operators with the commutation relation [38]

$\begin{split}\langle n_{a}(\tau)n_{a}(0)\rangle&=\langle
S_{a}^{\dagger}(0)S_{a}^{\dagger}(\tau)S_{a}(\tau)S_{a}(0)\rangle\\\ &=\langle
n_{a}^{s}(\tau)n_{a}^{s}(0)\rangle+\langle
n_{a}^{s}(\tau)n_{a}^{n}(0)\rangle\\\ &+\langle
n_{a}^{n}(\tau)n_{a}^{s}(0)\rangle+\langle
n_{a}^{n}(\tau)n_{a}^{n}(0)\rangle\\\ &+2\langle
n_{a}^{n}(0)n_{a}^{s}(0)\rangle.\\\ \end{split}$ (13)

This gives

$\begin{split}g^{(2)}_{aa}(\tau)&=\frac{\langle
n_{a}^{s}(\tau)n_{a}^{s}(0)\rangle}{\langle n_{a}^{s}(0)\rangle\langle
n_{a}^{s}(0)\rangle}\\\ &=\frac{\langle n_{a}(\tau)n_{a}(0)\rangle-\langle
n_{a}^{s}(\tau)\rangle\langle n_{a}^{n}(0)\rangle}{\langle
n_{a}^{s}(0)\rangle\langle n_{a}^{s}(0)\rangle}\\\ &+\frac{-\langle
n_{a}^{n}(\tau)\rangle\langle n_{a}^{s}(0)\rangle-\langle
n_{a}^{n}(\tau)n_{a}^{n}(0)\rangle}{\langle n_{a}^{s}(0)\rangle\langle
n_{a}^{s}(0)\rangle}\\\ &+\frac{-2\langle n_{a}^{n}(0)\rangle\langle
n_{a}(0)\rangle}{\langle n_{a}^{s}(0)\rangle\langle
n_{a}^{s}(0)\rangle},\end{split}$ (14)

where we have dropped the terms with the $\delta$ function. The expression for
$g^{(2)}_{bb}(\tau)$ can be derived similarly. In practice, we replace field
$a$ and $b$ with their corresponding complex envelop operator
$S_{a,b}(t)=X_{a,b}(t)+iP_{a,b}(t),$ where $X_{a,b}$ and $P_{a,b}$ are the two
quadratures of the complex envelope. To establish the correspondence between
the digitized voltage and the quadratures, we can rewrite the voltage as

$V(t)=X(t)\cos(2\pi ft)+P(t)\sin(2\pi ft),$ (15)

with $f=75$ MHz is the carrying frequency. With a sampling rate of $100$
MSp/s, the sampled voltage can be rewritten as

$V_{n}=V(n\Delta T)\begin{cases}X_{n}(-1)^{\frac{n}{2}}&\text{if $n$ is
even}\\\ P_{n}(-1)^{\frac{n+1}{2}}&\text{if $n$ is odd}\\\ \end{cases}$ (16)

where $\Delta T=$ 10 ns is the time interval between two samples. Eq. (16)
shows that the sample voltages are the two quadratures with alternating signs.
If we manually correct the sign, we are effectively demodulate the signal to
zero frequency. This process is known as frequency translation in digital
signal processing [43].

In order to proficiently filter $1/f$ noise, a Finite Impulse Response (FIR)
bandpass filter, incorporating a Hann window, has been utilized. In the
frequency domain, the filter’s frequency response, $H(\omega)$, is multiplied
with the signal spectrum. Concurrently, in the time domain, the signal is
convolved with $h(n)=\mathcal{F}^{-1}{H(\omega)}$, where the order, $N$, is
$200$. The frequency response can be viewed in Fig. 4 (inserted on the left)
and the bandwidth is 0.86 MHz.

The ability to measure the correlation at a smaller power (than -102.79 dBm)
is limited by two factors. Firstly, the quantization noise emanating from the
analog-to-digital conversion (ADC) imposes a fundamental limitation. Our
digitizer incorporates an 8-bit ADC, enabling the translation of analog
signals into integers ranging from -127 to 128. This conversion process
inevitably leads to quantization errors. The quantization error can be
quantified as quantization signal-to-noise ratio (SNR) and is given by
$\mathrm{SNR}\approx 6.08\times Q+1.76$ dB, where Q is 8 for 8 bits ADC [44].
In other word, this set a lower bound of the noise power due to the
quantization, defined by the signal power that takes the full dynamic range of
the ADC divided by the quantization SNR. This means for signal power smaller
than the quantization noise power would be challenging to measure. The second
factor is the emergence of the delta terms. As discussed in earlier, due to
the commutation relation, there will be $\delta(t)$ terms in the second-order
auto-correlation. The $\delta(t)$ terms would be small with a finite
bandwidth. However, if the signal component is small, the contribution from
$\delta(t)$ terms cannot be neglected. This will make the measurement very
unstable and not accurate [38]. In conclusion, we are not able to decrease the
pump power as the result becomes very volatile and not accurate, due to the
limitation of the measurement setups.

##  Appendix D Non-classicality analysis

In this section, we give a detailed analysis on the relation between the
parametric amplifier gain and the second-order cross-correlation. The input
and output relation of the interactivity mode can be expressed as

$\begin{split}a&=\sqrt{G_{\mathrm{NKPA}}}a_{\mathrm{in}}+\sqrt{G_{\mathrm{NKPA}}-1}b_{\mathrm{in}}^{\dagger},\\\
b&=\sqrt{G_{\mathrm{NKPA}}}b_{\mathrm{in}}+\sqrt{G_{\mathrm{NKPA}}-1}a_{\mathrm{in}}^{\dagger},\\\
\end{split}$ (17)

where $G_{\mathrm{NKPA}}=(\cosh{r})^{2}$ is the parametric gain of NKPA,
$a_{\mathrm{in}}$ and $b_{\mathrm{in}}$ are the intracavity modes. The mean
photon number of the output two-mode squeezed state is

$\begin{split}\langle a^{\dagger}a\rangle&=G_{\mathrm{NKPA}}\langle
a_{\mathrm{in}}^{\dagger}a_{\mathrm{in}}\rangle+(G_{\mathrm{NKPA}}-1)\langle
b_{\mathrm{in}}b_{\mathrm{in}}^{\dagger}\rangle\\\ \langle
b^{\dagger}b\rangle&=G_{\mathrm{NKPA}}\langle
b_{\mathrm{in}}^{\dagger}b_{\mathrm{in}}\rangle+(G_{\mathrm{NKPA}}-1)\langle
a_{\mathrm{in}}a_{\mathrm{in}}^{\dagger}\rangle\\\ \end{split}$ (18)

Considering intracavity photon occupancy $\langle
a_{\mathrm{in}}^{\dagger}a_{\mathrm{in}}\rangle=\langle
b_{\mathrm{in}}^{\dagger}b_{\mathrm{in}}\rangle=n$, then the above expression
can be simplified as

$\begin{split}\langle a^{\dagger}a\rangle=\langle
b^{\dagger}b\rangle&=2(G_{\mathrm{NKPA}}-1)n+G_{\mathrm{NKPA}}-1\\\
&=(G_{\mathrm{NKPA}}-1)/\eta,\end{split}$ (19)

where $\eta=1/(1+2n)$, and $\eta=1$ when both modes are at the ground state.
Here we consider that all the imperfection come from the non-zero intracavity
photon occupancy. The zero-delayed second-order correlation is linked to the
intracavity photon number as established by [23]

$g^{(2)}_{ab}(0)=2+\frac{\eta}{2(G_{\mathrm{NKPA}}-1)}.$ (20)

In the main text, the fitting result gives $\eta=0.72$. That means if all the
imperfections are from the non-zero intracavity photon occupancy, then
$n=(1/\eta-1)/2=0.19$. The unitless parameter $\eta$ can also be connected
with noise-reduction factor defined in reference [23] with relation
$\textrm{NRF}=1-\eta/2$. In our case, this yield $\textrm{NRF}=0.64$.

## References

  * Eichler  _et al._ [2011] C. Eichler, D. Bozyigit, C. Lang, M. Baur, L. Steffen, J. Fink, S. Filipp, and A. Wallraff, Observation of two-mode squeezing in the microwave frequency domain, Physical Review Letters 107, 113601 (2011).
  * Furusawa  _et al._ [1998] A. Furusawa, J. L. Sørensen, S. L. Braunstein, C. A. Fuchs, H. J. Kimble, and E. S. Polzik, Unconditional quantum teleportation, science 282, 706 (1998).
  * Milburn and Braunstein [1999] G. Milburn and S. L. Braunstein, Quantum teleportation with squeezed vacuum states, Physical Review A 60, 937 (1999).
  * Zhuang  _et al._ [2018] Q. Zhuang, Z. Zhang, and J. H. Shapiro, Distributed quantum sensing using continuous-variable multipartite entanglement, Physical Review A 97, 032329 (2018).
  * Lawrie  _et al._ [2019] B. J. Lawrie, P. D. Lett, A. M. Marino, and R. C. Pooser, Quantum sensing with squeezed light, Acs Photonics 6, 1307 (2019).
  * Backes  _et al._ [2021] K. M. Backes, D. A. Palken, S. A. Kenany, B. M. Brubaker, S. Cahn, A. Droster, G. C. Hilton, S. Ghosh, H. Jackson, S. K. Lamoreaux, _et al._ , A quantum enhanced search for dark matter axions, Nature 590, 238 (2021).
  * Wurtz  _et al._ [2021] K. Wurtz, B. Brubaker, Y. Jiang, E. Ruddy, D. Palken, and K. Lehnert, Cavity entanglement and state swapping to accelerate the search for axion dark matter, PRX Quantum 2, 040350 (2021).
  * Bienfait  _et al._ [2017] A. Bienfait, P. Campagne-Ibarcq, A. Kiilerich, X. Zhou, S. Probst, J. Pla, T. Schenkel, D. Vion, D. Estève, J. Morton, _et al._ , Magnetic resonance with squeezed microwaves, Physical Review X 7, 041011 (2017).
  * Didier  _et al._ [2015] N. Didier, A. Kamal, W. D. Oliver, A. Blais, and A. A. Clerk, Heisenberg-limited qubit read-out with two-mode squeezed light, Physical review letters 115, 093604 (2015).
  * Liu  _et al._ [2022] G. Liu, X. Cao, T.-C. Chien, C. Zhou, P. Lu, and M. Hatridge, Noise reduction in qubit readout with a two-mode squeezed interferometer, Physical Review Applied 18, 064092 (2022).
  * Lloyd and Braunstein [1999] S. Lloyd and S. L. Braunstein, Quantum computation over continuous variables, Physical Review Letters 82, 1784 (1999).
  * Gu  _et al._ [2009] M. Gu, C. Weedbrook, N. C. Menicucci, T. C. Ralph, and P. van Loock, Quantum computing with continuous-variable clusters, Physical Review A 79, 062318 (2009).
  * Tittel  _et al._ [1998] W. Tittel, J. Brendel, H. Zbinden, and N. Gisin, Violation of bell inequalities by photons more than 10 km apart, Physical review letters 81, 3563 (1998).
  * Weihs  _et al._ [1998] G. Weihs, T. Jennewein, C. Simon, H. Weinfurter, and A. Zeilinger, Violation of bell’s inequality under strict einstein locality conditions, Physical Review Letters 81, 5039 (1998).
  * Storz  _et al._ [2023] S. Storz, J. Schär, A. Kulikov, P. Magnard, P. Kurpiers, J. Lütolf, T. Walter, A. Copetudo, K. Reuer, A. Akin, _et al._ , Loophole-free bell inequality violation with superconducting circuits, Nature 617, 265 (2023).
  * Zhong  _et al._ [2020] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, _et al._ , Quantum computational advantage using photons, Science 370, 1460 (2020).
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Simon [2017] C. Simon, Towards a global quantum network, Nature Photonics 11, 678 (2017).
  * Wehner  _et al._ [2018] S. Wehner, D. Elkouss, and R. Hanson, Quantum internet: A vision for the road ahead, Science 362, eaam9288 (2018).
  * Bergeal  _et al._ [2012] N. Bergeal, F. Schackert, L. Frunzio, and M. Devoret, Two-mode correlation of microwave quantum noise generated by parametric down-conversion, Physical review letters 108, 123902 (2012).
  * Flurin  _et al._ [2012] E. Flurin, N. Roch, F. Mallet, M. H. Devoret, and B. Huard, Generating entangled microwave radiation over two transmission lines, Physical review letters 109, 183901 (2012).
  * Barzanjeh  _et al._ [2020] S. Barzanjeh, S. Pirandola, D. Vitali, and J. M. Fink, Microwave quantum illumination using a digital receiver, Science advances 6, eabb0451 (2020).
  * Westig  _et al._ [2017] M. Westig, B. Kubala, O. Parlavecchio, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, D. Esteve, M. Hofheinz, _et al._ , Emission of nonclassical radiation by inelastic cooper pair tunneling, Physical Review Letters 119, 137001 (2017).
  * Peugeot  _et al._ [2021] A. Peugeot, G. Ménard, S. Dambach, M. Westig, B. Kubala, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, _et al._ , Generating two continuous entangled microwave beams using a dc-biased josephson junction, Physical Review X 11, 031008 (2021).
  * Esposito  _et al._ [2022] M. Esposito, A. Ranadive, L. Planat, S. Leger, D. Fraudet, V. Jouanny, O. Buisson, W. Guichard, C. Naud, J. Aumentado, _et al._ , Observation of two-mode squeezing in a traveling wave parametric amplifier, Physical Review Letters 128, 153603 (2022).
  * Qiu  _et al._ [2023] J. Y. Qiu, A. Grimsmo, K. Peng, B. Kannan, B. Lienhard, Y. Sung, P. Krantz, V. Bolkhovsky, G. Calusine, D. Kim, _et al._ , Broadband squeezed microwaves and amplification with a josephson travelling-wave parametric amplifier, Nature Physics , 1 (2023).
  * Gasparinetti  _et al._ [2017] S. Gasparinetti, M. Pechal, J.-C. Besse, M. Mondal, C. Eichler, and A. Wallraff, Correlations and entanglement of microwave photons emitted in a cascade decay, Physical Review Letters 119, 140504 (2017).
  * Jiang  _et al._ [2023] Y. Jiang, E. Ruddy, K. Quinlan, M. Malnou, N. Frattini, and K. Lehnert, Accelerated weak signal search using mode entanglement and state swapping, PRX Quantum 4, 020302 (2023).
  * Brubaker  _et al._ [2017] B. Brubaker, L. Zhong, Y. Gurevich, S. Cahn, S. Lamoreaux, M. Simanovskaia, J. Root, S. Lewis, S. Al Kenany, K. Backes, _et al._ , First results from a microwave cavity axion search at 24 $\mu$ ev, Physical review letters 118, 061302 (2017).
  * Boutin  _et al._ [2017] S. Boutin, D. M. Toyli, A. V. Venkatramani, A. W. Eddins, I. Siddiqi, and A. Blais, Effect of higher-order nonlinearities on amplification and squeezing in josephson parametric amplifiers, Physical Review Applied 8, 054030 (2017).
  * Parker  _et al._ [2022] D. J. Parker, M. Savytskyi, W. Vine, A. Laucht, T. Duty, A. Morello, A. L. Grimsmo, and J. J. Pla, Degenerate parametric amplification via three-wave mixing using kinetic inductance, Physical Review Applied 17, 034064 (2022).
  * Xu  _et al._ [2023a] M. Xu, R. Cheng, Y. Wu, G. Liu, and H. X. Tang, Magnetic field-resilient quantum-limited parametric amplifier, PRX Quantum 4, 010322 (2023a).
  * Malnou  _et al._ [2022] M. Malnou, J. Aumentado, M. Vissers, J. Wheeler, J. Hubmayr, J. Ullom, and J. Gao, Performance of a kinetic inductance traveling-wave parametric amplifier at 4 kelvin: Toward an alternative to semiconductor amplifiers, Physical Review Applied 17, 044009 (2022).
  * Yu  _et al._ [2021] C. X. Yu, S. Zihlmann, G. Troncoso Fernández-Bada, J.-L. Thomassin, F. Gustavo, É. Dumur, and R. Maurand, Magnetic field resilient high kinetic inductance superconducting niobium nitride coplanar waveguide resonators, Applied Physics Letters 118, 054001 (2021).
  * Khalifa and Salfi [2023] M. Khalifa and J. Salfi, Nonlinearity and parametric amplification of superconducting nanowire resonators in magnetic field, Physical Review Applied 19, 034024 (2023).
  * Guo  _et al._ [2017] X. Guo, C.-l. Zou, C. Schuck, H. Jung, R. Cheng, and H. X. Tang, Parametric down-conversion photon-pair source on a nanophotonic chip, Light: Science & Applications 6, e16249 (2017).
  * Lvovsky [2015] A. I. Lvovsky, Squeezed light, Photonics: Scientific Foundations, Technology and Applications 1, 121 (2015).
  * Grimm [2015] A. Grimm,  _Josephson photonics: Statistics of photons emitted by inelastic Cooper pair tunneling_ , Ph.D. thesis, Grenoble University (2015).
  * Förtsch  _et al._ [2013] M. Förtsch, J. U. Fürst, C. Wittmann, D. Strekalov, A. Aiello, M. V. Chekhova, C. Silberhorn, G. Leuchs, and C. Marquardt, A versatile source of single photons for quantum information processing, Nature communications 4, 1818 (2013).
  * Xu  _et al._ [2020] M. Xu, X. Han, C.-L. Zou, W. Fu, Y. Xu, C. Zhong, L. Jiang, and H. X. Tang, Radiative cooling of a superconducting resonator, Physical Review Letters 124, 033602 (2020).
  * Xu  _et al._ [2023b] M. Xu, Y. Wu, W. Dai, and H. X. Tang, Radiatively-cooled quantum microwave amplifiers, arXiv preprint arXiv:2308.02106 (2023b).
  * Da Silva  _et al._ [2010] M. P. Da Silva, D. Bozyigit, A. Wallraff, and A. Blais, Schemes for the observation of photon correlation functions in circuit qed with linear detectors, Physical Review A 82, 043804 (2010).
  * Lyons [1997] R. G. Lyons,  _Understanding digital signal processing, 3/E_ (Pearson Education India, 1997).
  * Kester [2009] W. Kester, Taking the mystery out of the infamous formula,” snr= 6.02 n+ 1.76 db,” and why you should care, Analog Devices Tutorial, MT-001 Rev. A 10 (2009).