# Unsupervised and supervised learning of interacting topological phases
from single-particle correlation functions

Simone Tibaldi<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia
“Augusto Righi” dell’Università di Bologna, I-40127 Bologna, Italy  INFN,
Sezione di Bologna, I-40127 Bologna, Italy  Giuseppe Magnifico  Dipartimento
di Fisica e Astronomia “G. Galilei”, Università di Padova, I-35131 Padova,
Italy.  Padua Quantum Technologies Research Center, Università degli Studi di
Padova.  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Padova,
I-35131 Padova, Italy.  Davide Vodola<EMAIL_ADDRESS>Dipartimento di
Fisica e Astronomia “Augusto Righi” dell’Università di Bologna, I-40127
Bologna, Italy  INFN, Sezione di Bologna, I-40127 Bologna, Italy  Elisa
Ercolessi  Dipartimento di Fisica e Astronomia “Augusto Righi” dell’Università
di Bologna, I-40127 Bologna, Italy  INFN, Sezione di Bologna, I-40127 Bologna,
Italy

###### Abstract

The recent advances in machine learning algorithms have boosted the
application of these techniques to the field of condensed matter physics, in
order e.g. to classify the phases of matter at equilibrium or to predict the
real-time dynamics of a large class of physical models. Typically in these
works, a machine learning algorithm is trained and tested on data coming from
the same physical model. Here we demonstrate that unsupervised and supervised
machine learning techniques are able to predict phases of a non-exactly
solvable model when trained on data of a solvable model. In particular, we
employ a training set made by single-particle correlation functions of a non-
interacting quantum wire and by using principal component analysis, k-means
clustering, t-distributed stochastic neighbor embedding and convolutional
neural networks we reconstruct the phase diagram of an interacting
superconductor. We show that both the principal component analysis and the
convolutional neural networks trained on the data of the non-interacting model
can identify the topological phases of the interacting model. Our findings
indicate that non-trivial phases of matter emerging from the presence of
interactions can be identified by means of unsupervised and supervised
techniques applied to data of non-interacting systems.

Figure 1: Phase diagrams and correlation functions. (a) Phase diagram of the
non-interacting Kitaev model obtained from the winding number. The model
presents three phases: A trivial phase with winding number $\nu=0$ (green),
that we distinguish for convenience in TRI- when $\mu<-2$ and TRI+ when
$\mu>+2$, two superconducting topological phases with winding number $\nu=+1$
(TOP+1, orange) and $\nu=-1$ (TOP$-1$, red). The black lines represent phase
transition lines. (b) Correlation functions of the non interacting model
obtained from the analytical solution of the model. Each panel corresponds to
a different region of the phase diagram. (c) Phase diagram of the interacting
model. It presents four phases, of which three are trivial: Schrödinger-cat-
like phase (CAT, dark green dashed line), charge-density wave (CDW, blue),
trivial superconducting phase (TRI, green), and topological superconducting
phase (TOP, yellow). The phase transition points (black squares) were obtained
numerically by looking at the number of edge modes. (d) Correlation functions
of the interacting Kitaev model obtained through the DMRG algorithm. The four
graphs correspond to different regions of the phase diagram.

##  I Introduction

In the last few years research in physics has seen a new series of methods and
practices inspired by and exploiting machine learning Mehta et al. (2019).
These new instruments have proved to be useful in applications in many-body
quantum physics Carleo et al. (2019), in particular for finding a
representation of a wavefunction Carleo and Troyer (2017) and describing its
dynamics Czischek et al. (2018), for reconstructing the wavefunction from
experimental data Torlai et al. (2018), for speeding-up numerical simulations
Czischek et al. (2018), and for classifying quantum phases of matter of both
synthetic Wang (2016); van Nieuwenburg et al. (2017); Broecker et al. (2017);
Hu et al. (2017); Scheurer and Slager (2020); Rodriguez-Nieva and Scheurer
(2019) and experimental data Käming et al. (2021). A particular type of the
latter is given by symmetry protected and topologically ordered phases, whose
classification escapes from the standard Landau theory of spontaneous symmetry
breaking Sachdev (2011). Indeed the combination of symmetries and topology can
lead to new kinds of quantum phases that are characterized by a set of unusual
features, such as non-local order parameters, the appearance of zero energy
states at the boundary of the system, topological invariants, and long-range
entanglement (for reviews see, for example: Bernevig (2013); Zeng et al.
(2019)). From an experimental point of view, topological materials have
attracted much attention also because they represent a promising solution for
physical implementation of qubits, more resilient to decoherence processes
that affect devices based on superconducting or atomic physics technologies.

Due to their intrinsically non-local features and the lack of local order
parameters, the classification of topological phases is considered a very
challenging task to tackle Carleo et al. (2019). However, machine learning
methods have been successfully applied to both non-interacting models where
the topological invariant (winding number) representing each phase was known a
priori Zhang et al. (2018); Sun et al. (2018), and to interacting models where
the topological invariant cannot be obtained easily Huembeli et al. (2018).

Despite their success, in order to perform well, machine learning models
require data sets with a very large number of training data, in the order of
the thousands or millions. However, especially when handling interacting
systems, it is not always easy to build such big data sets from both numerical
simulations or experimental measurements. This led us to study a machine
learning model trained on a data set obtained from a solvable system to be
then applied to an interacting model which is obtained as an interacting
generalization of the former, as in Ref. Molignini et al. (2021). In this way,
insights about the features of a simple dataset can be exploited to
characterize the phases of a more complicated one, saving simulation
resources. Differently from what has been done in previous works Che et al.
(2020); Zhang et al. (2018), our dataset will be constructed out of two-point
single-particle correlation functions. These encode the properties of the
different phases of the model and can be obtained numerically and measured
experimentally Lobo and Gensemer (2008); Parsons et al. (2016); Cheuk et al.
(2016), e.g. by using atom gas microscope with quenches to non-interacting
models Gluza and Eisert (2021) or with randomized measurements Naldesi et al.
(2022).

In this work we show that supervised and unsupervised machine learning models
can classify topological phases of an interacting system by being trained on
data computed from simpler topological models. More specifically, our goal is
to use machine learning techniques to predict the topological and non-
topological quantum phases of two paradigmatic models: the one dimensional
Kitaev chain in its non-interacting Kitaev (2001) and interacting scenario
Stoudenmire et al. (2011); Hassler and Schuricht (2012); Thomale et al.
(2013); Katsura et al. (2015); Miao et al. (2017); Fromholz et al. (2020),
with the idea that one can exploit the knowledge of the easily solvable non-
interacting model in order to extract information also on the interacting
case.
We construct the correlator dataset from the analytical solution of the non-
interacting model and by means of numerical simulations implemented with the
Density Matrix Renormalization Group (DMRG) algorithm for the interacting
system Schollwöck (2005). Firstly, we probe the data with three unsupervised
methods: the principal components analysis (PCA), k-means clustering and
t-distributed stochastic neighbor embedding (t-SNE) Goodfellow et al. (2016);
MacKay (2003); van der Maaten and Hinton (2008). With PCA we investigate to
what extent the knowledge about principal component vectors of the data of the
non-interacting case can be used to learn also the underline pattern that
distinguishes the different phases in the interacting case. We use k-means’
ability to find clusters and centroids to correctly predict the number of
phases as well the location of the phase transition lines for both the non-
interacting and interacting case. Then, we check that also t-SNE is able to
form clusters of data.

Finally we devise an ensemble of convolutional neural networks (CNN) which is
trained on the non-interacting data and then tested to predict the phases of
the interacting model.

The paper is organized as follows: in Sec. II, we introduce the models and the
datasets that will be used for the training and testing of the machine
learning methods. In Sec. III, we apply unsupervised methods (PCA, k-means
clustering and t-SNE) to analyze the internal structure of the datasets. In
Section IV, we apply a supervised model trained on the data of the non
interacting model and test it on the interacting data. Finally, we draw our
conclusions in Sec. V.

##  II Models

In this section we describe the two models, the non-interacting and the
interacting Kitaev chain, whose quantum phases we want to classify, and we
define the standard and anomalous correlation functions that will be used as
indicator of the topological phases, thus providing the training and test
sets.

###  II.1 Non-interacting topological superconductor

We consider the one-dimensional Kitaev model Kitaev (2001) defined on a
lattice with $L$ sites described by the following non-interacting (NI)
Hamiltonian

$H^{\text{NI}}=\sum_{i}(Ja^{\dagger}_{i}a_{i+1}+\Delta\
a_{i}a_{i+1}+\text{h.c.})+\mu\sum_{i}a^{\dagger}_{i}a_{i}.$ (1)

Here the operators $a_{i}$ ($a^{\dagger}_{i}$) annihilate (create) a spinless
fermion on the lattice site $i$. The Hamiltonian $H^{\text{NI}}$ describes a
topological superconductor with nearest neighbour hopping of strength $J$,
p-wave superconducting pairing of strength $\Delta$ and chemical potential
$\mu$. When considering periodic boundary conditions, we can diagonalize
$H^{\text{NI}}$ by going to momentum space by means of standard Fourier
transform, so that it is reduced to a sum over the Brillouin zone (BZ),
$H^{\text{NI}}=\sum_{k\in BZ}H(k)$, of a two-band Hamiltonian
$H(k)=h_{z}(k)\sigma^{z}+h_{y}(k)\sigma^{y}$, where $k$ is the lattice quasi
momentum and

$\displaystyle h_{z}(k)=J\cos k+\mu/2\quad h_{y}(k)=\Delta\sin k$ (2)

and $\sigma^{x},\sigma^{y}$ are Pauli matrices. A Bogoliubov transformation
casts $H^{\text{NI}}$ in diagonal form
$H^{\text{NI}}=\sum_{k}E(k)\eta^{\dagger}_{k}\eta_{k},$ where $\eta_{k}$ are
Bogoliubov operators and the single-particle energy $E(k)$ is given by

$E(k)=\sqrt{h_{z}(k)^{2}+h_{y}(k)^{2}}.$ (3)

This model describes a one-dimensional topological superconductor belonging to
the BDI symmetry class Schnyder et al. (2008); Chiu et al. (2016); Kruthoff et
al. (2017); Slager et al. (2013). Its different phases are classified by the
winding number $\nu$ of the normalized Hamiltonian vector
$\hat{h}(k)=\vec{h}(k)/\|\vec{h}(k)\|$ with $\vec{h}(k)=(h_{y}(k),h_{z}(k))$,
which is a continuous map from the 1D BZ to the circle $S^{1}$. The winding
number $\nu$ is an integer that counts how many times the Hamiltonian vector
$\hat{h}(k)$ turns around the origin when the quasi-momentum $k$ moves from
$0$ to $2\pi$ in the 1D BZ. Panel (a) of Fig. 1 shows the phase diagram of the
model in Eq. (1) in the $(\mu,\Delta)$ plane, having set the energy scale
$J=1$. Notice that the phase diagram is symmetric for
$\mu\leftrightarrow-\mu$. Three different phases appear: a trivial phase with
winding number $\nu=0$ (green) and two non trivial phases with winding number
$\nu=\pm 1$ (orange/red). The winding number corresponds to the number of zero
energy states that the model hosts at the boundaries of the chain, when
considering open boundary conditions, a fact that is known in literature
Bernevig (2013); Asbóth et al. (2015) as bulk-edge correspondence. For $\nu=0$
no edge states will be present, while for $\nu=\pm 1$ an edge state on each
boundary appears.

The information on the different phases can also be extracted from the Fourier
transform of the single-particle standard ($c(k)$) and anomalous ($f(k)$)
correlation functions:

$\displaystyle
c(k)=\sum_{i,j}e^{\mathrm{i}k(i-j)}\braket{a^{\dagger}_{i}a_{j}}$ (4)
$\displaystyle f(k)=\sum_{i,j}e^{\mathrm{i}k(i-j)}\braket{a_{i}a_{j}}$ (5)

where the expectation values are taken on the ground state. We note that
$c(k)$ is real, while $f(k)$ is purely imaginary, due to the antisymmetry of
the expectation value $\braket{a_{i}a_{j}}$ for the exchange $i\leftrightarrow
j$. So we will redefine the latter by taking its imaginary part. For the
Kitaev model of Eq. (1), $c(k)$ and $f(k)$ can be computed analytically and
they take the form

$\displaystyle c(k)=\frac{1}{2}+\frac{\mu/2+J\cos k}{2{E(k)}},$ (6)
$\displaystyle f(k)=\frac{\Delta\sin k}{2{E(k)}}.$ (7)

Notice that they have a similar form to the components of the Hamiltonian
vector $\hat{h}(k)$ from which the winding number is calculated. Their
behaviour in the different phases is shown in panel (b) of Fig. 1.

The correlation functions $c(k)$ and $f(k)$ will be used for building a non-
interacting training set $S$ where each data point is labelled with the
winding number of its corresponding phase.

###  II.2 Interacting topological superconductor

We now add a nearest neighbor interaction term to the Hamiltonian (1) to
obtain:

$H^{\text{I}}=\sum_{i}(Ja^{\dagger}_{i}a_{i+1}+\Delta\
a_{i}a_{i+1}+\text{h.c.})+\mu\sum_{i}n_{i}+V\sum_{i}n_{i}n_{i+1}$ (8)

where $n_{i}=a^{\dagger}_{i}a_{i}$ is the occupation number at site $i$. This
model cannot be solved exactly due to the interacting potential. By means of
the DMRG algorithm Fishman et al. (2020), we have reproduced the phase
diagram, after setting $J=\Delta=1$, which is shown in panel (c) of Fig. 1,
for $\mu>0$ only since the model is symmetric for $\mu\leftrightarrow-\mu$.
The model in Eq. (8) is characterized by only one topological superconducting
phase (TOP, yellow) and three trivial phases: Topological trivial (TRI,
green), Charge Density Wave (CDW, light-blue) and a Schrödinger-cat-like phase
(CAT, light green dashed line) which shows up at the symmetric point $\mu=0$
as a superposition of two trivial superconducting states with different
occupation numbers Hassler and Schuricht (2012); Thomale et al. (2013);
Katsura et al. (2015); Miao et al. (2017). At $V=0$ we recover the non-
interacting case with a critical point at $\mu=2$. The different phases in
Fig. 1(c) are detected from the number of edge states that appear in the
chain: in the TRI, CDW and CAT phases the number of edge states is zero, while
it is one in the TOP phase.

In the interacting case, it is not possible to evaluate exactly the
correlation functions $c(k)$ (Eq. (4)) and $f(k)$ (Eq. (5)) on the ground
state of the Hamiltonian of Eq. (8). Therefore we calculate them by means of
the DMRG algorithm for a lattice of size $L=100$. Some examples of the
correlation functions for the different regions of the phase diagram are shown
in panel (d) of Fig. 1.

##  III Unsupervised training

Having obtained the datasets of the correlation functions for both the non-
interacting and and interacting model, we use three unsupervised methods,
namely PCA, k-means clustering and t-SNE, to extract the relevant information
in both datasets and predict the phases of both models.

###  III.1 Principal Components Analysis

Principal components analysis is a standard technique used in statistics and
machine learning for dimensionality reduction.

In order to apply PCA we start by creating a design matrix

$X=\begin{pmatrix}c_{1}(k_{0})&\dots&c_{1}(k_{L-1})&f_{1}(k_{0})&\dots&f_{1}(k_{L-1})\\\
\lx@intercol\hfil$\vdots$\hfil\lx@intercol\\\
c_{N}(k_{0})&\dots&c_{N}(k_{L-1})&f_{N}(k_{0})&\dots&f_{N}(k_{L-1})\end{pmatrix}$
(9)

where each of its $N$ rows is given by the correlation functions $c(k)$ and
$f(k)$ from Eqs. (4) and (5) for one point ($\mu_{\alpha},\Delta_{\alpha}$),
$\alpha=1,\dots,N$, of the phase diagram. Each column represents the Fourier
components of the correlation functions with quasi-momentum $k_{n}=2\pi n/L$
($n=0,\dots,L-1$) that are interpreted as the features of the data from which
the PCA extracts the principal components. We rescale the columns of $X$ such
that they have zero mean and unit standard deviation and we compute the
eigenvalues $\\{\lambda_{i}\\}$, the explained variances
$\epsilon_{i}=\lambda_{i}/\sum_{j}\lambda_{j}$, and the eigenvectors
$\\{\mathbf{w}_{i}\\}$ of the correlation matrix $\mathcal{S}=X^{T}X$. The
principal components of the data $X$, which are the eigenvectors of
$\mathcal{S}$ corresponding to the largest eigenvalues, are the directions in
a $2L$ dimensional space along which the original data show the largest
variance. A common measure of the projection along the principal components is
the quantified leading component (QLC) Wang (2016); Hu et al. (2017). For both
the non-interacting and the interacting case, we compute the QLC by dividing
the set of the two parameters in $40$ sections, creating a grid of $40\times
40$ subsets, each made of $M$ datapoints, as better specified below. Then, for
each of the subsets we calculate the quantity

$p_{i}=\sum_{s}\frac{\lvert X_{s}\cdot\mathbf{w}_{i}\rvert}{M}$ (10)

where $s$ runs on the elements of each subset, $X_{s}$ is the $s$-th row of
the matrix $X$ corresponding to the $s$-th point of the subset, and
$\mathbf{w}_{i}$ is the $i$-th eigenvector of $\mathcal{S}$. The value of the
QLC for a datapoint shows how much that point is represented by that specific
principal component.

Figure 2: PCA of the non-interacting Hamiltonian. The quantified leading
components (QLC) $p_{i}$ measure the projections of the data onto the space of
each of the first four principal components, corresponding to the eigenvalues
with largest explained variance $\epsilon_{i}$. (a) The first QLC with
$\epsilon_{1}=0.451$ reveals the points of the phase diagram with trivial
winding number. (b) The second QLC, with $\epsilon_{2}=0.421$, highlights the
points of the phase diagram with non-trivial winding number, (c) The third
QLC, with $\epsilon_{3}=0.052$ shows the phase transition lines, (d) The
fourth QLC with $\epsilon_{4}=0.042$ has a high projection on the phase with
$\nu=-1$. The insets show the first (red circles) and the second (orange
triangles) 100 elements of the corresponding principal components
$\mathbf{w}_{i}$ ($i=1,2,3,4$). Figure 3: PCA of the interacting Hamiltonian
Projection of the interacting correlation functions along the principal
components obtained applying PCA to the non interacting dataset. Darker colors
correspond to areas of the phase diagram which are more similar to the
principal component $\mathbf{w}_{i}$ ($i=1,2,3,4$) used for the projection.
QLCs drawn in (a) and (d) highlight the trivial (superconduting and charge
density wave) and topological phases respectively. In both cases, the CAT
phase is highlighted by higher values of $p_{i}$ compared to the other phases.
Component (c) pins down a phase transition line whereas component (b) does not
seem to be informative on the interacting dataset. The projections of all four
plots are rescaled independently. The dots are the same as panel (c) of Fig.
1, added to help recognize the phase transition points. Figure 4: Non-
interacting Hamiltonian: k-means clustering. (a) Values of the silhouette
score $\widetilde{S}$ as a function of the number of clusters $K$. The maximum
value is reached at $K=4$ and suggests the best number of centroids for the
dataset points. (b) Reconstructed phase diagram of the non-interacting model
(with $K=4$ centroids) through the average silhouette $S$ of every data point.
The critical lines are clearly visible because the correlation functions of
critical points have similar distances between two centroids. (c) Centroids
obtained by applying k-means algorithm with $K=4$. Every centroid resembles
the correlation functions of the four different regions: topological with
$\nu=-1$, the two trivial phases with $\nu=0$ and topological with $\nu=+1$,
from left to right. (d) Distances $d_{k}$ of the correlation functions for
each point of the phase diagram from the corresponding $k$-th centroid
(plotted in the corresponding panel above). The points belonging to the same
phase show minimal distance to the same centroid.

Non-interacting Hamiltonian – For the Hamiltonian of Eq. (1), we create the
desing matrix $X^{(\text{NI})}$ from data points generated in the range
$\mu\in[-8,8[$ and $\Delta\in[-2,2[$ with a step of $0.1$ for $\mu$ and 0.05
for $\Delta$, for a system with $L=100$ sites. This subdivision of the ranges
of $\Delta$ and $\mu$ corresponds to have a total of $N=12800$ data points.
The sum $\sum_{i=1}^{4}\epsilon_{i}$ of the explained variances of the first
four eigenvalues results in $96.6\%$ meaning that most of the information of
the data $X$ is contained in the first four eigenvectors. For this reason, in
Fig. 2, we show the first four QLCs. These are calculated as in Eq. (10) by
grouping the $N=12800$ datapoints in $40\times 40$ subsets corresponding to
$M=8$ datapoints per subset. The explained variances $\epsilon_{i}$ of each
eigenvector are reported for each of the first four components. In panel (a)
we see that $p_{1}$ is large for the points with $\lvert\mu\rvert>2$, this
means that the first principal component allows us to find the points of the
phase diagram with winding number $\nu=0$ (that belong to the trivial phases
TRI$-$ and TRI$+$). This is due to the shape of the first principal component
$\mathbf{w}_{1}$ (depicted in the inset of Fig. 2(a)) that resembles the shape
of the correlation functions $c(k)$ and $f(k)$ shown in Fig. 1(b) (TOP$+1$
phase). In Fig. 2(b), we see that, instead, $p_{2}$ allows us to extract the
points of the phase diagram with winding number $\nu=\pm 1$ as the shape of
the second principal component $\mathbf{w}_{2}$ (depicted in the inset of Fig.
2(b)) is similar to the correlation functions $c(k)$ and $f(k)$ shown in Fig.
1(b) (TOP$\pm 1$ phases). This analysis shows that the first two principal
components are sensitive to the trivial and non-trivial phases. All the other
$198$ QLCs have a total explained variance of the order of 13%, in particular
the third QLC (Fig. 2(c)) has explained variance $\epsilon_{3}=5.2\%$ and
recognizes the phase transition lines, while the fourth (Fig. 2(d)) has
explained variance $\epsilon_{4}=4.2\%$ and has a larger projection on the
phase with winding $\nu=-1$.

Interacting Hamiltonian – We are interested in understanding if the principal
components of the non-interacting model computed before can be used to
distinguish among the phases of the interacting Hamiltonian. To this end, we
construct a design matrix ${X}^{(\text{I})}$ with data obtained from the
interacting Hamiltonian in the range $\mu\in[0,5[$, $V\in[-4,4[$ with a step
of 0.1. Then, we calculate the QLC by projecting these data along the
principal components $\mathbf{w}_{i}$ of the non-interacting Hamiltonian. The
first four QLC are shown in Fig. 3. The first (panel (a)) and the fourth
(panel (d)) principal components $\mathbf{w}_{1}$ and $\mathbf{w}_{4}$ of the
non-interacting data show higher projections on the interacting data and are
thus able to discriminate between the trivial and topological phases,
respectively. We will comment in Sec. IV about the appearance of a ligther
colored region in the CDW phase, just on the right of the phase transition
line for large values of $\mu$. The third principal components
$\mathbf{w}_{3}$ (panel (c)) recovers only the transition line between the
trivial and the non-trivial superconducting phases while the second one seems
to highlight only the region for low $\mu$ and $V\in[-1,0]$. We note that,
even though the CAT phase is present only on the single line $\mu=0$, three
out of four principal components are able to recognize it (Fig. 3(a), (c),
(d)).

We can say with fairness that the PCA is able to learn the underlying pattern
which distinguishes a trivial phase from a topological one. This is a good
indication that even a supervised method can exploit the representation it
learns of the non-interacting data to classify the interacting ones.

Figure 5: Interacting Hamiltonian: k-means clustering (a) Values of the
silhouette score $\widetilde{S}$ as a function of the number of clusters $K$.
The peak is at $K=4$ in accordance to the 4 phases of the original phase
diagram. (b) Reconstructed phase diagram of the interacting model (with $K=4$
centroids) through the average silhouette value of every data point. Critical
lines separating the phases are clearly visible, CAT phase is also
highlighted. (c) Centroids obtained by applying k-means algorithm with $K=4$.
Every centroid resembles the correlation functions of the four different
phases, as indicated by the label in the bottom right corner. (d) Distances of
the correlation functions for each point of the phase diagram from the
corresponding centroid (plotted in the panel above). The points belonging to
the same phase show minimal distance to the same centroid.

###  III.2 K-means clustering

K-means clustering is a machine learning method for finding clusters and
cluster centers in a set of unlabelled data MacKay (2003). The algorithm
starts by choosing a number of cluster centers called centroids and then it
iteratively moves the centers to minimize the total variance within cluster.
The centroids are the central point of every cluster that are calculated by
the algorithm autonomously, so we interpret them as the most representative
point of the cluster.

In order to settle on a value of $K$ without any a priori knowledge of our
data we choose to calculate the silhouette score $\widetilde{S}$ which is a
value representing the average quality of the clusterization for each $K$.
Calling $a$ the distance of a sample point to its centroid and $b$ its
distance to the second closest centroid we can calculate the silhouette of a
data point as $S=(b-a)/\max(a,b)$. The values for $S$ lie in the range
$[-1,1]$ where negative numbers correspond to a wrongly assigned point.
Positive values indicate the right classification and the quality increases
reaching 1.

We run k-means algorithm with different $K$ values and select the $K$ with the
largest silhouette score. Due to the sensitivity of the algorithm to initial
conditions we run k-means 10 times for each $K$ and then collect the average
silhouette value $\widetilde{S}$ of every point over the 10 runs, whereas the
centroids and projections correspond to the largest silhouette only.

Figure 6: Phase diagrams with k-means. Reconstruction of the phase diagrams of
the non-interacting (a) and the interacting (b) models obtained with the
labels assigned by k-means. For each data point we plot its label predicted by
k-means running with $k=4$, the value chosen according to the best silhouette
as explained in the text.

Non-interacting Hamiltonian –  The analysis of the silhouette for the non-
interacting data turns out to be very informative. Firstly, after multiple
iterations as explained above, we find the maximum value of the silhouette
score at $K=4$ as shown in Fig. 4(a). This suggests that the most reasonable
way to divide the data points is in 4 classes, which might correspond to 4
different phases. Secondly, we see that the silhouette of a point can be used
for identifying the phase transition lines. In fact, the points lying near a
phase transition might be reasonably associable to the two different phases of
the transition, that is they might show properties of both phases. So we
expect their silhouette value to be close to 0, indicating that their
classification might be non ideal. This is indeed the case, as it is seen from
Fig. 4(b) which shows the average silhouette of every point calculated over 10
iterations of k-means algorithm applied with $K=4$, as suggested by the
previous analysis. We can see that points that are inside the phases have a
larger silhouette than points lying at the phase transition.

Let us now move to the analysis of the centroids obtained by applying k-means
algorithm with $K=4$. In Fig. 4(c) and (d) we plot the 4 centroids and the
distance of the points to each centroid, respectively. Each centroid seems to
exactly represent the features of the datapoints, i.e. the correlation
functions, of the four different regions of the phase diagram. Also, all the
points of a phase have a distance very close to zero to their centroid,
showing that we can separate the phases of the diagram with ease.

Interacting Hamiltonian – The same analysis can be repeated on the interacting
dataset, obtaining similar results which are collected in Fig. 5. In
particular the silhouette reaches its peak at $K=4$ corresponding to the four
phases of the interacting model shown in Fig. 1(c).

In Fig. 6, for completeness, we show the scatter plots of the labels assigned
by k-means ran with $k=4$ for both datasets. The reconstruction of the phase
diagrams is neat and in accordance with the results of Figs. 4 and 5.

We can summarize the results of this section saying that, in both the non-
interacting and the interacting models, the k-means clustering algorithm, by
trying to separate points, is able to learn the characteristic shape of the
correlation functions of each phase and to identify the boundaries where the
phase transitions occur.

Figure 7: Clustering with t-SNE. The plot shows the visualization of the data
projected in 2D applying t-SNE (components $y_{1},y_{2}$). In both panels we
assign a color to each point corresponding to its phase in the model while the
brightness of each point corresponds to its silhouette (see Section III.B).
(a) Projection of the non-interacting model data points. Four clusters appear
evidently and they correspond to the 3 phases of the non-interacting model
with the two trivial phases (TRI+, TRI-) separated. The silhouette of the
points (brightness of the color) highlights the points close to a phase
transition, which happen to be at the borders of each cluster. (b) Projection
of the interacting model data points. The separation of the clusters is less
neat compared to (a) but thanks to the silhouette we can see the borders of
each cluster. Figure 8: Structure of the CNNs forming the ensemble. The input
for the CNNs is a $2\times L$ dataset with the correlation functions
$c(k),f(k)$. First, the input is filtered by a layer of $m_{1}$ 2D
convolutional weights (in yellow) trained independently on each network. Each
weight produces one array of length $N_{1}=L-1$ (the first set of green
columns denoted by $\ell_{1}$) for a total of $m_{1}$ arrays. Then each
network has a different number of additional hidden layers (columns in red):
$\ell_{2},\ell_{3},\ell_{4}$ from the top. Each of these layers performs a 1D
convolution using the weights in yellow. The last convolution always has one
single weight in order to produce a single array. The last layer is a feed-
forward layer (FF$\ell$) with one final output neuron that gives a real number
which is the predicted topological indicator $\omega$ by the network. All the
resulting topological indicators from the different CNNs that constitute the
ensemble are then averaged to give the predicted topological indicator
$\widetilde{\omega}$ (Eq. (15)). Figure 9: Convolutional neural networks. (a)
Predictions $\widetilde{\omega}$ of the phase diagram of the interacting
Hamiltonian carried out by means of the CNN ensemble. The topological phase is
correctly predicted (dark red region). (b) Two cuts of the phase diagram
showing the predictions $\widetilde{\omega}$ and their relative standard
errors as a function of $V$ for different $\mu$. For the points with
$\mu\gtrsim 1$ in the CDW phase it is possible to see the appearance of
discontinuities (indicated by dashed vertical lines) for the values
$\widetilde{\omega}$ that correspond to the onset of the incommensurate CDW
phase.

###  III.3 t-distributed stochastic neighbor embedding

t-SNE van der Maaten and Hinton (2008) is a popular dimensionality reduction
technique that creates a low-dimensional distribution of data which is
faithful to the original one and thus helps visualizing clusters of data in 2
or 3 space. It starts by creating a probability distribution from the
Euclidean distances between data points in their original space. In particular
given two points $x_{i},x_{j}\in\mathcal{D}$ where $\mathcal{D}$ is the
dataset, we interpret the conditional probability $p_{i|j}$ as a measure of
similarity between the two points under a Gaussian centered at $x_{i}$:

$p_{j|i}=\frac{\text{exp}(||x_{i}-x_{j}||^{2}/2\sigma^{2})}{\sum_{k\neq
i}\text{exp}(||x_{i}-x_{k}||^{2}/2\sigma^{2})}$ (11)

where $\sigma$ is a parameter that we will discuss later. t-SNE creates a
symmetric joint probability $P$ from 11 by taking
$p_{ij}=(p_{i|j}+p_{j|i})/2$. Then, it gives to each point $x_{i}$ a new set
of coordinates in a lower-dimensional space $y_{i}\in\mathbb{R}^{d}$ for $d=2$
or 3, where the similarity between points is given by a $t$-Student
distribution $Q$ in the form:

$q_{ij}=\frac{||y_{i}-y_{j}||^{2}}{\sum_{k\neq l}||y_{k}-y_{l}||^{2}}$ (12)

In order to make $q_{ij}$ a faithful low-dimensional representation of the
distribution $p_{ij}$, t-SNE minimizes the Kullback-Leibler divergence

$C=KL(P||Q)=\sum_{i}\sum_{j}p_{ij}\log\frac{p_{ij}}{q_{ij}}$ (13)

Typically this is done updating the position of points $y_{i}$ by gradient
descent:

$y_{i}=y_{i}-\eta\frac{\partial C}{\partial y_{i}}$ (14)

with $\eta$ the learning rate. The parameter $\sigma$ appearing in 11 is
related to another hyperparameter called perplexity which is the one that
needs to be set when running the t-SNE. Perplexity can be seen as the average
number of neighbors we expect the points to have in the original space and
typical values range from 5 to 50 van der Maaten and Hinton (2008).

We ran t-SNE on our datasets varying both perplexity and learning rate $\eta$
without seeing major changes to the final result. Also, no relevant
distinction was found in either reducing the dimensions to 3 or 2, so we stick
to 2D.

In Figure 7 we plot the results of t-SNE on the non-interacting/interacting
dataset. For both models we plot the dimensional reduction to 2 dimensions. In
the non-interacting case we see four well-separated clusters and each one
corresponds to one of the phases of the model. In particular t-SNE separates
the two trivial phases (TRI+, TRI-) of the non interacting model in the same
way as K-means does. For the interacting case the clusters are not as well-
separated except for the CAT phase. Although the shape of the of the low-
dimensional representation is not meaningful in t-SNE cluster, we are
interested in finding the position of the points close to the phase transition
lines. For this reason we adjust the brightness of each data point according
to its silhouette value that was calculated in section III.B (compare with
Figs. 4 and 5). In this way we are able to see that the points closer to a
phase transition are close to the margin of the cluster.

##  IV Supervised Training

In this section we exploit the information gained from the unsupervised
analysis of the data to use a supervised model for predicting the different
(trivial or topological) phases via the analysis of correlation functions. Our
aim is to train the network on the data of the non-interacting Hamiltonian and
test it on the interacting data. This is an approach that has been recently
exploited in the context of machine learning applied to systems without
analytical solution, e.g. Valenti et al. (2019).

For this reasons we decided to employ an ensemble method which is favourable
even when tackling difficult classification problems Carleo et al. (2019). The
base element of the ensemble we use is a convolutional neural network (CNN), a
type of machine learning model that has found many successful applications in
image recognition problems Krizhevsky et al. (2012) and vastly employed in
applications to quantum many body problems Carrasquilla and Melko (2017); Ming
et al. (2019); Sun et al. (2018); Zhang et al. (2018). This network uses small
matrices of weights, called features, which perform a series of discrete
convolutions and are able to extract local properties of the data, for details
see Goodfellow et al. (2016).

CNN ensemble – The ensemble of classifiers is made of $C=180$ CNNs, as
schematically shown in Fig. 8. Each CNN has an initial hidden layer $\ell_{1}$
(in green) of $m_{1}=100$ filters (weights) of size $2\times 2$ (yellow in the
picture) and stride 1, which produce $m_{1}=100$ arrays of size
$N_{1}=L-1=99$. Each network can have three different internal structures
whose schemes are reported in Fig. 8. The number of additional hidden layers
$\ell_{\alpha}$ (in red) can vary from 1 to 3 while each of the hidden layers
in a network has the same number of weights per layer which can be
$m=25,50,100$. Finally, the networks have been trained with different batch
sizes (512, 1024) and different random initial training configuration of the
weights centered at zero (10 different starting seeds for each network). Each
network accepts as input the $2\times L$ matrix made by stacking $c(k)$ and
$f(k)$ from Eqs. (6)-(7) and produces a single real number $\omega_{\alpha}$
in output which is the estimate for the topological indicator of the input.
All layers neurons have ReLU activation ($f(x)=\max\\{0,x\\}$) except for the
last one which is linear ($f(x)=x$). The loss function is given by the mean
squared error. The output of each network is then averaged to produce:

$\widetilde{\omega}=\frac{1}{C}\sum_{\alpha=1}^{C}\omega_{\alpha}.$ (15)

In the following we will show that this quantity can be used as a topological
indicator in order to reconstruct the phase diagram of the interacting model.

Training – The networks are trained separately. The dataset of the non-
interacting model is made of $2\times 10^{5}$ points covering the whole phase
diagram of the non interacting Hamiltonian of Fig. 1(a). Each network is
trained with Adam gradient descent method by using an early stopping technique
until convergence is achieved. This typically requires around 100 steps.

Testing – For the test dataset we considered $4\times 10^{3}$ evenly spaced
samples in the grid $[-4,4]\times[0,5]$ for $V$ and $\mu$ in order to fully
reconstruct the interacting phase diagram of Fig. 1(c). The ensemble, trained
on the non-interacting data, efficiently evaluates the topological
classification of the test set resulting in the predicted phase diagram shown
in Fig. 9(a). We note that the predictions $\omega_{\alpha}$ of the 180 CNNs
are very consistent: they are constrained in the range $[0,1]$ and their
standard deviation is $\sim 9\%$.

All the points of the SC and CAT phases are correctly classified as non
topological and the predicted values of $\tilde{\omega}$ show a sharp
transition between the TRI and TOP phases, meaning that the ensemble was able
to learn how to identify the two different phases of the superconductor. On
the other hand, there is a less marked transition between the TOP and CDW
phases. In Fig. 9(b) we show two horizontal cuts of the phase diagram taken
for $\mu=0.5$ and $\mu=4$. Interestingly, in the CDW phase for $\mu\gtrsim 1$,
the values of $\widetilde{\omega}$ present discontinuities as a function of
the potential $V$ (indicated by the dashed vertical lines for $\mu=4$ in Fig.
9(b)) that correspond to the boundaries of an incommensurate CDW that has been
found e.g. in Thomale et al. (2013). The appearence of this additional phases
is also caught by the PCA (Fig. 3(a, d)) and k-means (Fig. 5(b)) results.

##  V Conclusions

In this work we have shown how one can use machine learning techniques to
predict the topological and non-topological quantum phases of a paradigmatic
model, namely the interacting Kitaev chain, by exploiting the knowledge of the
easily solvable non-interacting model. More specifically, we have constructed
the non-interacting dataset from the analytical solutions of the standard and
anomalous correlation functions of the model, while we have used the DMRG
algorithm to calculate the correlation functions in the interacting case.

At first, we have probed our data with PCA, k-means clustering and t-SNE. With
the former, we have confirmed that the non-interacting data contain enough
information to predict the main features of the interacting dataset as well.
With k-means, we were able to identify the right number of phases as well as
correctly locate the phase transition lines. Also, we have seen that the
correlation functions of the four centroids reproduce the pattern of the
correlation functions of the different phases. With t-SNE, we were able to see
clusters of data corresponding to the phases of the models and, thanks to the
silhouette values obtained from k-means, the localization of the phase
transition points at the borders of the various clusters.

Then we have used an ensemble of CNN, which was trained on the non-interacting
dataset and then tested to calculate and predict the phases of the interacting
model. The ensemble does efficiently evaluate the topological class of this
model, with only the data points of the topological superconducting phase
producing a topological indicator equal to one on average. This means that the
ensemble network is able to reconstruct the phase diagram indifferently of the
shapes and phases of the input data.

These findings clearly indicate that non-trivial phases of matter that emerge
in presence of interactions can be identified by means of unsupervised and
supervised techniques applied to data of non-interacting systems. These
results offer a number of advantages since the data of non-interacting or
solvable quantum many-body systems can be easily generated by means of
analytical computation, numerical simulations, or directly measured in state-
of-the-art experiments. Furthermore, all the protocols developed in this work
could be easily generalized to higher-dimensional systems or finite-
temperature regimes.

###### Acknowledgements.

We thank E. Tignone, F. Dell’Anna, S. Pradhan for useful discussions. This
research is funded by Istituto Nazionale di Fisica Nucleare (INFN) through the
project “QUANTUM”, the International Foundation Big Data and Artificial
Intelligence for Human Development (IFAB) through the project “Quantum
Computing for Applications” and the QuantERA 2020 Project “QuantHEP”. G. M. is
partially supported by the Italian PRIN2017 and the Horizon 2020 research and
innovation programme under grant agreement No 817482 (Quantum Flagship -
PASQuanS). We acknowledge the use of computational resources from the parallel
computing cluster of the Open Physics Hub
(https://site.unibo.it/openphysicshub/en) at the Physics and Astronomy
Department in Bologna.

## References

  * Mehta et al. (2019) P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to Machine Learning for physicists, Physics Reports 810, 1 (2019).
  * Carleo et al. (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019).
  * Carleo and Troyer (2017) G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Czischek et al. (2018) S. Czischek, M. Gärttner, and T. Gasenzer, Quenches near Ising quantum criticality as a challenge for artificial neural networks, Phys. Rev. B 98, 024311 (2018).
  * Torlai et al. (2018) G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-network quantum state tomography, Nature Physics 14, 447 (2018).
  * Wang (2016) L. Wang, Discovering phase transitions with unsupervised learning, Phys. Rev. B 94, 195105 (2016).
  * van Nieuwenburg et al. (2017) E. P. L. van Nieuwenburg, Y.-H. Liu, and S. D. Huber, Learning phase transitions by confusion, Nature Physics 13, 435 (2017).
  * Broecker et al. (2017) P. Broecker, J. Carrasquilla, R. G. Melko, and S. Trebst, Machine learning quantum phases of matter beyond the fermion sign problem, Scientific Reports 7, 8823 (2017).
  * Hu et al. (2017) W. Hu, R. R. P. Singh, and R. T. Scalettar, Discovering phases, phase transitions, and crossovers through unsupervised machine learning: A critical examination, Phys. Rev. E 95, 062122 (2017).
  * Scheurer and Slager (2020) M. S. Scheurer and R.-J. Slager, Unsupervised Machine Learning and Band Topology, Phys. Rev. Lett. 124, 226401 (2020).
  * Rodriguez-Nieva and Scheurer (2019) J. F. Rodriguez-Nieva and M. S. Scheurer, Identifying topological order through unsupervised machine learning, Nature Physics 15, 790 (2019).
  * Käming et al. (2021) N. Käming, A. Dawid, K. Kottmann, M. Lewenstein, K. Sengstock, A. Dauphin, and C. Weitenberg, Unsupervised machine learning of topological phase transitions from experimental data, Machine Learning: Science and Technology 2, 035037 (2021).
  * Sachdev (2011) S. Sachdev, Quantum Phase Transitions, 2nd ed. (Cambridge University Press, Cambridge, England, 2011).
  * Bernevig (2013) B. A. Bernevig, Topological Insulators and Topological Superconductors (Princeton University Press, 2013).
  * Zeng et al. (2019) B. Zeng, X. Chen, D.-L. Zhou, and X.-G. Wen, Quantum Information Meets Quantum Matter (Springer, 2019).
  * Zhang et al. (2018) P. Zhang, H. Shen, and H. Zhai, Machine Learning Topological Invariants with Neural Networks, Phys. Rev. Lett. 120, 066401 (2018).
  * Sun et al. (2018) N. Sun, J. Yi, P. Zhang, H. Shen, and H. Zhai, Deep learning topological invariants of band insulators, Phys. Rev. B 98, 085402 (2018).
  * Huembeli et al. (2018) P. Huembeli, A. Dauphin, and P. Wittek, Identifying quantum phase transitions with adversarial neural networks, Phys. Rev. B 97, 134109 (2018).
  * Molignini et al. (2021) P. Molignini, A. Zegarra, E. van Nieuwenburg, R. Chitra, , and W. Chen, A supervised learning algorithm for interacting topological insulators based on local curvature, SciPost Phys. 11, 73 (2021).
  * Che et al. (2020) Y. Che, C. Gneiting, T. Liu, and F. Nori, Topological quantum phase transitions retrieved through unsupervised machine learning, Phys. Rev. B 102, 134213 (2020).
  * Lobo and Gensemer (2008) C. Lobo and S. D. Gensemer, Technique for measuring correlation functions in interacting gases, Phys. Rev. A 78, 023618 (2008).
  * Parsons et al. (2016) M. F. Parsons, A. Mazurenko, C. S. Chiu, G. Ji, D. Greif, and M. Greiner, Site-resolved measurement of the spin-correlation function in the Fermi-Hubbard model, Science 353, 1253 (2016).
  * Cheuk et al. (2016) L. W. Cheuk, M. A. Nichols, K. R. Lawrence, M. Okan, H. Zhang, E. Khatami, N. Trivedi, T. Paiva, M. Rigol, and M. W. Zwierlein, Observation of spatial charge and spin correlations in the 2D Fermi-Hubbard model, Science 353, 1260 (2016).
  * Gluza and Eisert (2021) M. Gluza and J. Eisert, Recovering Quantum Correlations in Optical Lattices from Interaction Quenches, Phys. Rev. Lett. 127, 090503 (2021).
  * Naldesi et al. (2022) P. Naldesi, A. Elben, A. Minguzzi, D. Clément, P. Zoller, and B. Vermersch, Fermionic correlation functions from randomized measurements in programmable atomic quantum devices, (2022), arXiv:2205.00981 .
  * Kitaev (2001) A. Y. Kitaev, Unpaired Majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001).
  * Stoudenmire et al. (2011) E. M. Stoudenmire, J. Alicea, O. A. Starykh, and M. P. Fisher, Interaction effects in topological superconducting wires supporting Majorana fermions, Phys. Rev. B 84, 014503 (2011).
  * Hassler and Schuricht (2012) F. Hassler and D. Schuricht, Strongly interacting Majorana modes in an array of Josephson junctions, New Journal of Physics 14, 125018 (2012).
  * Thomale et al. (2013) R. Thomale, S. Rachel, and P. Schmitteckert, Tunneling spectra simulation of interacting Majorana wires, Phys. Rev. B 88, 161103 (2013).
  * Katsura et al. (2015) H. Katsura, D. Schuricht, and M. Takahashi, Exact ground states and topological order in interacting Kitaev/Majorana chains, Phys. Rev. B 92, 115137 (2015).
  * Miao et al. (2017) J.-J. Miao, H.-K. Jin, F.-C. Zhang, and Y. Zhou, Exact Solution for the Interacting Kitaev Chain at the Symmetric Point, Phys. Rev. Lett. 118, 267701 (2017).
  * Fromholz et al. (2020) P. Fromholz, G. Magnifico, V. Vitale, T. Mendes-Santos, and M. Dalmonte, Entanglement topological invariants for one-dimensional topological superconductors, Phys. Rev. B 101, 085136 (2020).
  * Schollwöck (2005) U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, and A. Courville, Deep Learning (MIT Press, 2016).
  * MacKay (2003) D. J. C. MacKay, Information Theory, Inference, and Learning Algorithms (Cambridge University Press, 2003).
  * van der Maaten and Hinton (2008) L. van der Maaten and G. Hinton, Viualizing data using t-SNE, Journal of Machine Learning Research 9, 2579 (2008).
  * Schnyder et al. (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * Chiu et al. (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Kruthoff et al. (2017) J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Topological Classification of Crystalline Insulators through Band Structure Combinatorics, Phys. Rev. X 7, 041069 (2017).
  * Slager et al. (2013) R.-J. Slager, A. Mesaros, V. Juričić, and J. Zaanen, The space group classification of topological band-insulators, Nature Physics 9, 98 (2013).
  * Asbóth et al. (2015) J. Asbóth, L. Oroszlány, and A. Pályi, A Short Course on Topological Insulators: Band-structure topology and edge states in one and two dimensions, Lecture Notes in Physics, Vol. 919 (Springer Nature, 2015).
  * Fishman et al. (2020) M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, (2020), arXiv:2007.14822 .
  * Valenti et al. (2019) A. Valenti, E. van Nieuwenburg, S. Huber, and E. Greplova, Hamiltonian learning for quantum error correction, Phys. Rev. Research 1, 033092 (2019).
  * Krizhevsky et al. (2012) A. Krizhevsky, I. Sutskever, and G. E. Hinton, ImageNet Classification with Deep Convolutional Neural Networks, in Advances in Neural Information Processing Systems, Vol. 25, edited by F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger (Curran Associates, Inc., 2012).
  * Carrasquilla and Melko (2017) J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nature Physics 13, 431 (2017).
  * Ming et al. (2019) Y. Ming, C.-T. Lin, S. D. Bartlett, and W.-W. Zhang, Quantum topology identification with deep neural networks and quantum walks, npj Computational Materials 5, 88 (2019).


# Spectral and temporal characterization of a fused-quartz microresonator
optical frequency comb

Scott B. Papp and Scott A. Diddams  National Institute of Standards and
Technology, Boulder, Colorado 80305, USA

###### Abstract

We report on the fabrication of high-$Q$, fused-quartz microresonators and the
parametric generation of a frequency comb with 36 GHz line spacing using them.
We have characterized the intrinsic stability of the comb in both the time and
frequency domains to assess its suitability for future precision metrology
applications. Intensity autocorrelation measurements and line-by-line comb
control reveal near-transform-limited picosecond pulse trains that are
associated with good relative phase and amplitude stability of the comb lines.
The comb’s 36 GHz line spacing can be readily photodetected, which enables
measurements of its intrinsic and absolute phase fluctuations.

Femtosecond-laser optical frequency combs have revolutionized frequency
metrology and precision timekeeping by providing a dense set of absolute
reference lines spanning more than an octave. These sources exhibit sub-
femtosecond timing jitter corresponding to an ultralow phase noise of $<100$
$\mu$rad on high harmonics of their typically 100’s of MHz repetition
frequency (line spacing) Diddams  _et al._ (2004). This remarkable level of
performance has enabled measurements of atomic transition frequencies at the
17th digit Rosenband  _et al._ (2008) and the generation of ultralow noise
signals Fortier  _et al._ (2011). Beyond their natural application as an
optical clockwork, frequency combs are used in diverse applications including
precision measurements of fundamental physics Diddams  _et al._ (2004);
Rosenband  _et al._ (2008), direct spectroscopy and real-time trace gas
detection Thorpe  _et al._ (2006), molecular fingerprinting Diddams  _et al._
(2007), astronomical spectrograph calibration Murphy  _et al._ (2007), and
optical arbitrary waveform generationHuang  _et al._ (2008). An even broader
range of applications may be possible if the bulk and complexity of a
femtosecond comb could be reduced without significantly increasing
fluctuations in the comb spectrum.

A new class of frequency combs has recently emerged based on monolithic
microresonators, henceforth denoted microcombs Del’Haye  _et al._ (2007);
Savchenkov  _et al._ (2008a); Levy  _et al._ (2010); Razzari  _et al._ (2010);
Kippenberg  _et al._ (2011). Here the comb generation relies on parametric
conversion provided simply by third-order nonlinear optical effects and is
enabled by recent advances in the quality factor $Q$ and the small volume of
microresonators. These devices require only a single continuous-wave laser
source, but the usable frequency span of the comb depends on low dispersion,
making material properties critical. Microcombs also present a unique platform
for creating large line spacings (10’s of GHz to THz). To date microcombs
possessing 100’s of lines each spaced by 100’s of GHz have been created with
silica microtoroids Del’Haye  _et al._ (2007) and silicon-nitride microrings
Levy  _et al._ (2010); Razzari  _et al._ (2010). Uniformity Del’Haye  _et al._
(2007) and control Del’Haye  _et al._ (2008) of silica microcombs have also
been demonstrated. Microcombs with mode spacings of 10 to 25 GHz have also
been realized with machined crystalline resonators Savchenkov  _et al._
(2008a) and aspects of their microwave-domain spectral purities have been
analyzed Savchenkov  _et al._ (2008a, b).

Figure 1: High-$Q$ fused-quartz microresonators. (a) A fused quartz rod is
rotated in a lathe (not shown) and rubbed against a fixed metal form charged
with diamond abrasive. The ground form is then flame polished with household-
grade propane and oxygen (not shown). (b) Image showing an array of fused
quartz microresonators after grinding and flame polishing. The disks of 1.9 mm
diameter are separated by 0.45 mm. (c) Frequency scan of a mode for disk four
(indicated by the arrow) demonstrating an optical $Q_{0}=5.2\times 10^{8}$.
(d) Cavity ringdown at $K=1$ to determine the optical $Q$ of disk four. The
dashed lines indicate a decay time of 260 ns, corresponding to
$Q_{0}=6.2\times 10^{8}$.

Future frequency metrology applications of microcombs will require a line
spacing (repetition rate, $f_{\rm{rep}}$) in the 10’s of GHz (millimeter-scale
resonator), comb span approaching an octave, and low absolute phase and
frequency fluctuations that can be further reduced by wideband comb control
Kippenberg  _et al._ (2011). Additionally, a threshold power for comb
generation in the mW range, and the capability for integration with chip-based
photonic circuits would enable novel portable applications. Here we present a
platform towards attaining these requirements. Our fabrication approach takes
advantage of precision mechanical shaping and polishing possible with fused
quartz, a versatile amorphous material with low optical losses throughout the
visible and into the infrared. We created resonators that feature high $Q$ and
low effective optical mode area, which allowed comb generation for a $<$10 mW
pump laser. By pumping above threshold, a comb with $\sim 150$ lines spaced by
36 GHz was generated. We harnessed full line-by-line control (amplitude and
phase) of the comb to generate near-transform limited optical waveforms and to
understand the comb’s relative phase stability. In the frequency domain, we
characterized the comb’s mode-spacing fluctuations to assess prospects for
stabilizing it to optical atomic standards.

Figure 2: Microresonator lineshape and corresponding comb spectra for laser
pump power of 62 mW. (a) Thermally broadened lineshape with a mostly sawtooth
shape. Power levels L1, L2, L3 respectively correspond to the parametric
threshold, the initiation of a 36 GHz comb but with unstable behavior, and the
region of stable 36 GHz comb generation. (b) Comb spectrum at L1. (c) Comb
spectrum at L2. (d) Comb spectrum at L3.

We fabricated an array of disk-like microresonators on a common substrate by
use of a combination of diamond grinding and flame polishing. A schematic of
our procedure is shown in Fig. 1a. Fused quartz rods with $\approx 2$ mm
diameter were rotated in a ball-bearing spindle. To create the basic shape of
a disk, the glass was rubbed against a metal form with triangular cross
section. Diamond abrasive between them removed material from both. At a
diamond abrasive size of 6 $\mu$m, $\sim$30 minutes of grinding was required
to create the underlying shape shown in Fig. 1b. The resulting surface was far
too rough to support high $Q$ whispering-gallery modes. By polishing the disks
with a propane-oxygen flame we achieved an extremely smooth surface.
Substantial melting of the quartz was required to alleviate subsurface damage
caused by grinding. Figure 1b shows an image of a completed disk array. Free-
spectral range measurements of several disks revealed a diameter uniformity of
0.1 $\%$.

After processing, the fused quartz rod was placed in a temperature-stabilized
mount ($<$0.01 K stability) and whispering-gallery modes of the disks were
accessed using standard fiber tapers Knight  _et al._ (1997); Cai  _et al._
(2000); Spillane  _et al._ (2003). The coupling between disk modes and the
taper can be characterized by the parameter $K=\tau_{0}/\tau_{e}$, where
$\tau_{0}$ is the intrinsic photon lifetime and $\tau_{e}$ is the photon loss
rate due to coupling. The frequency linewidth of the resonator-taper system is
$\gamma=\tau_{0}^{-1}+\tau_{e}^{-1}$. Figure 1c shows a frequency scan over a
disk mode using a tunable diode laser at 1560 nm. A low optical power of
$\sim$10 $\mu$W was used to minimize thermal broadening of the resonance
Carmon  _et al._ (2004) and taper coupling was small ($K\approx 0.03$). The
linewidth of the Lorentzian transmission feature is 0.37 MHz, which
corresponds to an optical $Q_{0}=\omega\,\tau_{0}$ of $5.2\times 10^{8}$,
where $\omega$ is the optical carrier frequency 111We have not observed any
scattering-induced mode splitting in these devices Gorodetsky  _et al._
(2000); Kippenberg  _et al._ (2002).. We confirmed this $Q_{0}$ measurement by
recording the resonator energy ringdown time following a rapid frequency scan
across resonance with $K=1$; see Fig. 1d. A fit to the cavity energy decay
yields $Q_{0}=6.2\times 10^{8}$, in good agreement with our linewidth
measurement. Material absorption in the fused-quartz samples we use imposes a
limit to $Q_{0}\sim 2\times 10^{10}$. We speculate that some contamination
occurs during the flame-polishing step Kuzuu  _et al._ (2003); *Kuzuu2004.

To achieve the long-term stability in $K$ needed for studies of microcomb
generation, we place the taper directly in contact with the resonator surface.
Translating the taper along the axis of the quartz rod allows us to tune the
resonator-taper coupling from under- to over-coupling. Moreover, measurements
of $K$ at various locations along the disk edge provide a qualitative picture
of the mode’s spatial profile. Specifically, the resonator mode studied in
Fig. 1 c and d has two field nodes in the axial direction and is presumably of
lowest radial order. Numerical calculations based on an ideal disk model
Borselli  _et al._ (2005) suggest the effective area of this mode
($A_{\rm{eff}}$) is approximately 500 $\mu$m2. The threshold for parametric
oscillation threshold is
$P_{\rm{th}}=\frac{1}{4}\,\frac{1+K}{K}\,\frac{n}{n_{2}}\,\frac{\omega}{\Delta\nu}\,\frac{A_{\rm{eff}}}{Q^{2}}$,
where $n$ ($n_{2}$) is the refractive index (nonlinear index) of fused quartz
and $\Delta\nu$ is the free spectral range. For our disk resonators we expect
$P_{\rm{th}}=$ 4 mW for $K=0.2$.

We have generated a microcomb at, and significantly above, parametric
threshold. Figure 2a-d shows a disk resonator transmission lineshape and
typical comb optical spectra. With constant laser power, we adjusted and
stabilized the intracavity power by tuning the pump laser frequency for a
fixed $K=0.2$. Here the lineshape is thermally broadened and assumes a
characteristic “sawtooth” shape Carmon  _et al._ (2004). Power level L1
indicates the observed 3.7 mW threshold for parametric oscillation, which is
in excellent agreement with our prediction. Notably, the signal and idler
modes are spaced by 828 GHz ($23\times 36$ GHz); indicating the point at which
resonator-induced phase mismatch of pump, signal, and idler balances the
mismatch created by nonlinear effects Agrawal (2007); Kippenberg  _et al._
(2004); Savchenkov  _et al._ (2008a). Increasing the intracavity power further
resulted in a widened comb span up to $\pm$25 nm about 1560 nm and a reduction
of the comb spacing toward the fundamental 36 GHz of the disk resonator, which
is consistent with the simulations of Chembo et al. Chembo  _et al._ (2010);
Chembo and Yu (2010). However, our observations indicate that 36 GHz comb
generation is complex with competition between different modes of operation.
Specifically, we present two instances of 36 GHz combs with dramatically
different properties: (1) At power level L2 (Fig. 2c) an $\approx$36 GHz comb
spacing is initiated, but a finer analysis (see Fig. 6c) reveals numerous comb
line spacing frequencies and (2) in a range about power level L3 (Fig. 2d) a
robust comb pattern exists with a single line spacing of 36.0012 GHz.
Interestingly, when we tuned to power level L2 a deviation in the sawtooth
shape of the transmission lineshape was observed. This behavior depends on
resonator-taper coupling and laser polarization, but not the laser frequency
scan rate, which we varied from 0.005 to 1 s. Moreover, in this range stable
time-domain waveforms and low noise frequency-domain performance were not
evident under the conditions of Fig. 2 222The numerical simulation of
parametric comb generation by Chembo appears to explain the spectrum of Fig.
2d with 36 GHz lines modulated at an approximately 828 GHz period, but not the
sawtooth deviation at power level L2 or the associated instability of comb
generation here Chembo  _et al._ (2010); Chembo and Yu (2010). About L3 we
observed the resumption of a sawtooth lineshape.

Figure 3: Time-domain measurements of a three line comb. (a) Schematic of our
setup to control and analyze the microcomb spectrum. The amplitude and phase
of light exciting the resonator via tapered fiber is manipulated with a
programmable filter. An intensity autocorrelator is used to record the fields’
waveforms. (b) Optical spectrum following the filter and amplifier. (c) Time-
averaged autocorrelation signal for $\Delta\phi=0$ (black circles) and
$\Delta\phi=\pi$ (blue circles), which respectively demonstrate constructive
and destructive interference of the three comb fields. The solid lines show
the results of an analytic model. (d) Visibility of constructive interference
as a function of $\Delta\phi$.

The spectral structure similar to Fig. 2 has been observed in a variety of
parametric microcombs, but little is known about the temporal structure of the
output as determined by the phase relationship between the comb elements.
Measurements and control of optical waveforms can reveal crucial information
about the internal fields of the microcomb Foster  _et al._ (2011); Ferdous
_et al._ (2011). Using time-domain diagnostics, we studied the relative phase
stability of the first three comb lines (pump, signal, and idler) at
parametric threshold, and the phase stability of many lines far above
threshold. For a comb composed of only three equal-amplitude lines its time
waveform has two distinct behaviors, which depend primarily on the relative
phases of the lines: constructive interference, resulting in a rudimentary
pulse at $f_{\rm{rep}}$, or partial destructive interference, resulting in a
beat signal at $2f_{\rm{rep}}$. A comb with many lines exhibits more
complicated behavior and offers high peak power waveforms.

Figure 4: Generation of picosecond pulses utilizing line-by-line comb control.
(a) Measured time-domain waveform (black line) following amplitude and phase
optimization, and simulation (red dashed line) based on the measured
amplitudes of comb lines used and assuming phase alignment of all the fields.
The inset shows the optical spectrum after optimization. (b) Optical spectrum
following broadening in highly nonlinear fiber. The comb was amplified to 1.6
W.

Our experimental starting point was a comb spectrum similar to that shown in
Fig. 2b, except with spacing of 864 GHz 333The mode spacing at threshold can
be tuned over a small range near 864 GHz with $K$.. Using a programable
amplitude and phase filter featuring 10 GHz resolution (Fig. 3a), we selected
the central three lines of this comb, equalized their amplitudes to better
than 7 %, and adjusted their phases. The optical spectrum we obtained
following amplification with a low-dispersion erbium fiber amplifier is shown
in Fig. 3b. It was delivered via a short section of optical fiber to an
intensity autocorrelator for characterization. By varying the relative phase
($\Delta\phi$) of the lowest wavelength line, we were able to explore the
transition between constructive and deconstructive interference of the comb
lines. The relative phases between the center and the highest wavelength line
were set to zero. Figure 3c shows the time-waveform autocorrelation signal
($P(\tau)$), where $\tau$ is the differential path delay, for $\Delta\phi=0$
(black open circles) and $\Delta\phi=\pi$ (blue open circles). The solid lines
indicate the results of an analytic model based on three phase-coherent
fields, which are in excellent agreement with the data. To characterize the
interference between the three fields we introduce a visibility $V$ based on
$P(\tau)$, which is defined as
$V=|\frac{P(\pi/\Delta\omega)-P(0)}{P(\pi/\Delta\omega)+P(0)}|$. The filled
points in Fig. 3d are measurements of $V$ as a function of $\Delta\phi$, and
the red line shows our analytic model. We measured a maximum (minimum)
visibility of 0.7 (0.007), and the model predicts the visibility is 0.768
(0.002). This indicates good phase coherence of the comb lines. Fitting the
model to the data with free parameter $\Delta\phi$ yields an uncertainty in
$\Delta\phi$ of 20 mrad, which we associate with an upper limit on the
fluctuations of the relative phase of the comb elements over the measurement
time of 1700 seconds 444We note this value also includes un-characterized
fluctuations intrinsic to the programable filter and the measurement
apparatus..

We also carried out time-domain waveform generation with a microcomb far above
parametric threshold. The amplitude and phase filter was used to adjust the
central 15 lines (0.54 THz) of a 36 GHz comb similar to that of Fig. 2d. Here
a manual line-by-line procedure was used to optimize the peak waveform
intensity. This resulted in stable near-transform-limited 2.5 ps optical
pulses, as demonstrated by the time-averaged autocorrelation signal in Fig.
4a. Hence, the relative phases of the 15 comb lines were presumed to be
constant across the spectrum and were correspondingly stationary in time.
Moreover, the autocorrelation signal remained unchanged despite large
perturbations to the microcomb system, such as re-locking the pump laser
frequency or re-placing the tapered fiber in contact with the disk. These
observations indicate that the internal fields of the resonator and the pulse
generation mechanism are deterministic and repeatable. Significantly, the
creation of time-domain waveforms, such as shown in Fig. 4a, implies that such
a compact and simple device could be competitive with more conventional mode-
locked lasers a source of high repetition rate ultrashort optical pulses that
would be useful for a variety of time- and frequency-domain applications
Kippenberg  _et al._ (2011). For example, nonlinear broadening in fibers could
be a route to ultra-broadband spectra for comb self-referencing and comb
spectroscopy. Using a 101 m length of highly nonlinear fiber, we have
broadened our 15 line comb by more than a factor of 10, as shown in Fig. 4b.
This demonstrates sufficient peak power to drive nonlinear processes external
to the microresonator.

Figure 5: Residual phase noise of microcomb line spacing. (a) Schematic of our
experimental setup. The microcomb spectrum is split (BS) into two paths, which
are separately filtered at $\lambda_{1}=1560$ nm and $\lambda_{2}=1553.6$ nm.
Photodectors PD1,2 convert the comb spacing of the independent paths into
electrical signals near 36 GHz. PD1 (PD2) generates 2.15 mA (0.04 mA) of
photocurrent. Double-balanced mixers M1 and M2 convert these signals to
baseband, and mixer M3 offsets one path by 11 MHz. This offset frequency was
referenced to a hydrogen maser. At mixer $M_{4}$ the two paths are interfered
yielding residual noise. (b) Single-sideband phase noise of the microcomb 36
GHz tone (black line) and the contribution from the 11 MHz offset oscillator
(red line). (c) Histogram of zero-dead-time counter measurements of the signal
exiting M4. (d) Allan deviation of the microcomb line spacing (filled points).
The solid line shows the $1/\tau$ scaling exhibited at most averaging times.

Frequency-domain techniques offer significant advantages in characterizing
microcomb fluctuations, including access to the power spectral density, wide
measurement bandwidth, and extremely high precision. A key feature of our
system is a comb line spacing of 36 GHz, which enables direct photodetection
and analysis of the resulting microwave signal in the frequency domain. Here
we present measurements of the microcomb’s residual and absolute phase
fluctuations. Residual (i.e. intrinsic) noise indicates the degree that the
frequencies of different comb lines are correlated, and absolute noise
indicates the stability of the lines with respect to a fixed frequency
reference. Moreover, we have characterized the relationship between line
spacing fluctuations and pump laser intensity noise. This information provides
a benchmark for what level of control will be required for future applications
of stabilized microcombs.

To understand the intrinsic stability of the microcomb, we directly compared
the line spacing at two independent portions of its spectrum. We created two
copies of the microcomb spectrum shown in Fig. 2d, using a fiber beamsplitter,
and separately bandpass filtered them at 1560 nm and 1553.6 nm; see Fig. 5a.
Independent high-speed photodetectors created electronic signals at 36 GHz.
The frequencies of these signals were reduced to baseband with a common 35.1
GHz oscillator at mixers M1 and M2, and the path corresponding to filter
$\lambda_{1}$ was further shifted by 11 MHz with a low-phase-noise oscillator.
A final mixer (M4) interfered the two baseband signals to reject common-mode
fluctuations. Hence, only residual phase noise of the microcomb (and noise
associated with the independent paths) appeared at the output of M4. Using a
commercial phase noise analyzer we recorded the single-sideband phase noise
spectral density $L_{\phi}$, which is shown by the black line in Fig. 5b. For
reference, the red line in Fig. 5b shows the phase noise contribution of the
11 MHz offset signal. The microcomb’s residual phase noise is extremely low.
At carrier offsets greater than 10 kHz, $L_{\phi}$ is dominated by
photodetection noise and would be reduced with greater microcomb optical
power. Close to carrier the residual phase noise is mostly below -100 dBc/Hz,
a value comparable to the absolute phase noise of the best optical and
microwave oscillators Fortier  _et al._ (2011). This indicates that with
appropriate frequency control such a microcomb could be harnessed for portable
low-noise synthesis applications.

Figure 6: Microcomb fluctuations measured with respect to fixed references.
For all of these experiments the entire microcomb spectrum was delivered to
PD2. (a) Free-running microwave spectrum of the microcomb operating at power
level L3 obtained via photodetection. (b) Power spectral density (PSD) of
$L_{\phi}$ (black), line-spacing signal amplitude (red), and pump laser RIN
(blue). The RIN attains its shot noise limited value of -151 dBc/Hz beyond 100
kHz. The cyan shaded region shows a scaling of the RIN to phase noise using
the measured conversion factor $\gamma_{P}$ for our system. The top range of
this area corresponds to a constant $\gamma_{P}$ with frequency, and the
bottom range accounts for the measured frequency response of pump power to
line spacing. The green line shows the phase noise contribution of the 35.5
GHz oscillator. (c) Microwave spectrum at L2 with a 30 kHz resolution
bandwidth. Under these conditions the comb exhibits multiple line spacings
near 36 GHz. (d) $L_{\phi}$ (black) and amplitude (red) noise of the microcomb
at L2. The blue trace shows pump laser RIN.

We have assessed whether the microcomb line spacing is the same at 1560 nm and
1553.6 nm by electronically counting the 11 MHz signal exiting M4. Here an
offset from 11 MHz would indicate non-equidistance of the microcomb spectrum.
Previous work by Del’Haye et. al. studied microcomb equidistance by way of an
auxiliary fiber-laser frequency comb Del’Haye  _et al._ (2007). With access to
the 36 GHz comb spacing, our measurements directly characterize microcomb
equidistance for the first time. We acquired 546 consecutive (zero dead time)
one-second long measurements; Fig. 5c shows a histogram of the frequency
difference ($\Delta f$) between the microcomb and reference oscillator. The
data exhibit a Gaussian distribution with mean of -0.65 $\mu$Hz with a 25
$\mu$Hz width. Hence, the mode spacing of the comb does not change
fractionally by more than 3$\times 10^{-17}$ across 6.4 nm. The Allan
deviation $\sigma_{A}$, shown in Fig. 5d, demonstrates the intrinsic stability
of the microcomb line spacing with $\sigma_{A}$=$1\times 10^{-17}$ after only
400 s of averaging.

Many future applications of a microcomb will require a stable spectrum with
respect to fixed frequency references, such as atomic clocks. Here we show the
microcomb line spacing’s microwave spectrum, absolute phase noise, amplitude
noise, and its dependence on pump laser power. Importantly, different behavior
of these was observed at power levels L2 and L3 that was not revealed in
coarse measurements of the optical spectrum (Fig. 2c and d). The experimental
setup for these measurements was similar to Fig. 5a, except that the
beamsplitter, optical filters, and mixers M1,3,4 were removed. The
photocurrent generated at PD2 was converted to baseband with a 35.5 GHz
oscillator for analysis. For phase noise measurements we used a digital
prescaler (divide-by-265) following M1 to reject spurious amplitude
fluctuations.

First, we present our observations at L3. Figure 6a shows a single microwave
tone 55 dB above a broad asymmetric pedestal. The width of this feature is
limited by the resolution bandwidth of 10 kHz. The amplitude and phase
fluctuations of this tone are shown by the red and black traces in Fig. 6b,
respectively. Notably, the absolute phase noise is between 20 and 120 dB above
the residual noise, depending on the carrier offset. For comparison the green
line shows the phase noise of our commercial 35.5 GHz oscillator. The pump
laser power, with the relative intensity noise (RIN) spectrum shown by the
blue trace in Fig. 6b, significantly influences the comb line spacing. We
measured a line spacing-power dependence of $|\gamma_{P}|$ = 0.3 MHz/mW with a
3 dB bandwidth of $\sim$20 Hz. The correlation between pump laser RIN and
microcomb amplitude and phase noise is evident in Fig. 6b. Scaling the RIN by
$\gamma_{P}$ results in the phase noise prediction shown by the cyan-shaded
area. Pump laser power fluctuations explain a significant fraction of the line
spacing phase noise, and either passive or active reduction of the RIN will be
crucial for future experiments. In particular, obtaining a $f$-2$f$ heterodyne
beatnote Diddams  _et al._ (2000); Jones  _et al._ (2000); Del’Haye  _et al._
(2009), a key outstanding milestone in microcomb systems, will be challenging
in the face of significant comb noise.

Figure 6c and d shows a similar analysis of the microcomb line spacing’s
spectrum (c), amplitude (red trace in (d)) and phase fluctuations (black trace
in (d)) at power level L2. Under this operating mode, the comb exhibits
multiple line spacings in a few MHz range around 36 GHz. We observe
significantly larger (10’s of dB) amplitude and phase noise, even though the
pump laser RIN is the same as at L3. The correspondence between pump RIN and
amplitude, phase noise remains evident. We have not been able to develop a
detailed understanding of, nor any strategy to mitigate, the poor performance
at L2. However, these observations emphasize that measurements beyond those of
the optical spectra will be critical to understanding the properties of the
microcomb in all its operating regimes.

In summary, we have fabricated mm-scale, high-$Q$ optical microresonators with
fused quartz. By using them, a frequency comb with 36 GHz line spacing, 50 nm
span, and $<$10 mW threshold was created near 1560 nm. We studied the
microcomb spectrum in both the time and frequency domains. Picosecond optical
pulses were generated by way of full line-by-line control of up to 15 comb
lines. Precise frequency-domain techniques enabled a direct test of comb line
fractional equidistance at the 3$\times 10^{-17}$ level, and a
characterization of the microcomb’s absolute line spacing phase noise.
Quantifying the noise of microresonator frequency combs is crucial in
assessing whether they can distribute modern microwave and optical atomic
standards. In the future we will focus generating a wider comb span, and
stabilizing the comb to Rb transitions at 780 nm by way of higher bandwidth
comb control.

We acknowledge useful conversations with Kerry Vahala and Tobias Kippenberg,
and we thank Gabe Ycas and Nathan Lemke for a thoughtful reading of the
manuscript. This work is a contribution of NIST and is not subject to
copyright in the United States.

## References

  * Diddams  _et al._ (2004) S. A. Diddams, J. C. Bergquist, S. R. Jefferts, and C. W. Oates, Science 306, 1318 (2004).
  * Rosenband  _et al._ (2008) T. Rosenband, D. B. Hume, P. O. Schmidt, C. W. Chou, A. Brusch, L. Lorini, W. H. Oskay, R. E. Drullinger, T. M. Fortier, J. E. Stalnaker, S. A. Diddams, W. C. Swann, N. R. Newbury, W. M. Itano, D. J. Wineland, and J. C. Bergquist, Science 319, 1808 (2008).
  * Fortier  _et al._ (2011) T. M. Fortier, M. S. Kirchner, F. Quinlan, J. Taylor, J. C. Bergquist, T. Rosenband, N. Lemke, A. Ludlow, Y. Jiang, C. W. Oates, and S. A. Diddams, arXiv:1101.3616v3 (2011).
  * Thorpe  _et al._ (2006) M. J. Thorpe, K. D. Moll, R. J. Jones, B. Safdi, and J. Ye, Science 311, 1595 (2006).
  * Diddams  _et al._ (2007) S. A. Diddams, L. Hollberg, and V. Mbele, Nature 445, 627 (2007).
  * Murphy  _et al._ (2007) M. T. Murphy, T. Udem, R. Holzwarth, A. Sizmann, L. P. C. Araujo-Hauck, H. Dekker, S. D’Odorico, M. Fischer, T. Hansch, and A. Manescau, Mon. Not. R. Astron. Soc. 380, 839 (2007).
  * Huang  _et al._ (2008) C.-B. Huang, Z. Jiang, D. Leaird, J. Caraquitena, and A. Weiner, Laser & Photonics Reviews 2, 227 (2008).
  * Del’Haye  _et al._ (2007) P. Del’Haye, A. Schliesser, O. Arcizet, T. Wilken, R. Holzwarth, and T. J. Kippenberg, Nature 450, 1214 (2007).
  * Savchenkov  _et al._ (2008a) A. A. Savchenkov, A. B. Matsko, V. S. Ilchenko, I. Solomatine, D. Seidel, and L. Maleki, Phys. Rev. Lett. 101, 093902 (2008a).
  * Levy  _et al._ (2010) J. S. Levy, A. Gondarenko, M. A. Foster, A. C. Turner-Foster, A. L. Gaeta, and M. Lipson, Nat Photon 4, 37 (2010).
  * Razzari  _et al._ (2010) L. Razzari, D. Duchesne, M. Ferrera, R. Morandotti, S. Chu, B. E. Little, and D. J. Moss, Nat Photon 4, 41 (2010).
  * Kippenberg  _et al._ (2011) T. J. Kippenberg, R. Holzwarth, and S. A. Diddams, Science 332, 555 (2011).
  * Del’Haye  _et al._ (2008) P. Del’Haye, O. Arcizet, A. Schliesser, R. Holzwarth, and T. J. Kippenberg, Phys. Rev. Lett. 101, 053903 (2008).
  * Savchenkov  _et al._ (2008b) A. A. Savchenkov, E. Rubiola, A. B. Matsko, V. S. Ilchenko, and L. Maleki, Optics Express 16, 4130 (2008b).
  * Knight  _et al._ (1997) J. C. Knight, G. Cheung, F. Jacques, and T. A. Birks, Opt. Lett. 22, 1129 (1997).
  * Cai  _et al._ (2000) M. Cai, O. Painter, and K. J. Vahala, Phys. Rev. Lett. 85, 74 (2000).
  * Spillane  _et al._ (2003) S. M. Spillane, T. J. Kippenberg, O. J. Painter, and K. J. Vahala, Phys. Rev. Lett. 91, 043902 (2003).
  * Carmon  _et al._ (2004) T. Carmon, L. Yang, and K. Vahala, Opt. Express 12, 4742 (2004).
  * Note (1) We have not observed any scattering-induced mode splitting in these devices Gorodetsky  _et al._ (2000); Kippenberg  _et al._ (2002).
  * Kuzuu  _et al._ (2003) N. Kuzuu, H. Horikoshi, T. Nishimura, and Y. Kokubo, J. Appl. Phys. 93, 9062 (2003).
  * Kuzuu  _et al._ (2004) N. Kuzuu, Y. Kokubo, T. Nishimura, I. Serizawa, L.-H. Zeng, K. Fujii, M. Yamaguchi, K. Saito, and A. J. Ikushima, J. Non-Cryst. Solids 333, 115 (2004).
  * Borselli  _et al._ (2005) M. Borselli, T. Johnson, and O. Painter, Opt. Express 13, 1515 (2005).
  * Agrawal (2007) G. P. Agrawal,  _Nonlinear Fiber Optics_ (Academic Press, 2007).
  * Kippenberg  _et al._ (2004) T. J. Kippenberg, S. M. Spillane, and K. J. Vahala, Phys. Rev. Lett. 93, 083904 (2004).
  * Chembo  _et al._ (2010) Y. K. Chembo, D. V. Strekalov, and N. Yu, Phys. Rev. Lett. 104, 103902 (2010).
  * Chembo and Yu (2010) Y. K. Chembo and N. Yu, Phys. Rev. A 82, 033801 (2010).
  * Note (2) The numerical simulation of parametric comb generation by Chembo appears to explain the spectrum of Fig. 2d with 36 GHz lines modulated at an approximately 864 GHz period, but not the sawtooth deviation at power level L2 or the associated instability of comb generation here Chembo  _et al._ (2010); Chembo and Yu (2010).
  * Foster  _et al._ (2011) M. A. Foster, J. S. Levy, O. Kuzucu, K. Saha, M. Lipson, and A. L. Gaeta, arXiv:1102.0326v1 (2011).
  * Ferdous  _et al._ (2011) F. Ferdous, H. Miao, D. E. Leaird, K. Srinivasan, J. Wang, L. Chen, L. T. Varghese, and A. M. Weiner, arXiv:1103.2330v1 (2011).
  * Note (3) The mode spacing at threshold can be tuned over a small range near 864 GHz with $K$.
  * Note (4) We note this value also includes un-characterized fluctuations intrinsic to the programable filter and the measurement apparatus.
  * Diddams  _et al._ (2000) S. A. Diddams, D. J. Jones, J. Ye, S. T. Cundiff, J. L. Hall, J. K. Ranka, R. S. Windeler, R. Holzwarth, T. Udem, and T. W. H nsch, Phys. Rev. Lett. 84, 5102 (2000).
  * Jones  _et al._ (2000) D. J. Jones, S. A. Diddams, J. K. Ranka, A. Stentz, R. S. Windeler, J. L. Hall, and S. T. Cundiff, Science 288, 635 (2000).
  * Del’Haye  _et al._ (2009) P. Del’Haye, T. Herr, E. Gavartin, R. Holzwarth, and T. Kippenberg, arXiv:0912.4890v1 (2009).
  * Gorodetsky  _et al._ (2000) M. L. Gorodetsky, A. D. Pryamikov, and V. S. Ilchenko, J. Opt. Soc. Am. B 17, 1051 (2000).
  * Kippenberg  _et al._ (2002) T. J. Kippenberg, S. M. Spillane, and K. J. Vahala, Opt. Lett. 27, 1669 (2002).