\section{Phase measurements}

\subsection{Absolute phase measurements with independent laser}

In our experiments, the phases of laser pulses are measured via the interference with
the independent CW laser as a local oscillator.
We describe how to estimate the absolute phase of each pulse from the interference with the independent CW laser.

Let $\omega\sub{s}$ denote one of the cavity mode frequencies of the mode-locked laser.
The amplitude $\tilde{A}_k(t)$ of the $k$-th pulse is denoted by
\begin{align}
  \tilde{A}_k(t) = A_k(t){\mathrm{e}}^{{\mathrm{i}} \omega\sub{s}(t + \tau\sub{R}k/N)},  \label{eq:1}
\end{align}
where $A_k(t)$ is the slowly varying amplitude, $\tau\sub{R}$ is the cavity round-trip time, and $N$ is the number of pulses inside the cavity.

Let $\omega\sub{LO}$ denote the angular frequency of the CW laser.
The wavelength of the CW laser is variable 
and we adjusted the wavelength so as to be overlapped with the spectrum of the mode-locked fiber laser.
Since the frequency comb of the mode-locked fiber laser has the angular frequencies with the spacing of $\omega\sub{cav} = 2 \pi \tau\sub{R}^{-1}$,
there exists the angular frequency $\omega\sub{s}$ that satisfies the following condition:
\begin{align}
  |\omega\sub{s} - \omega\sub{LO}| < \frac{\omega\sub{cav}}{2}.  \label{eq:2}
\end{align}

The measured amplitude $A_k\sur{(meas)}(t)$ with reference to the CW laser is written as
\begin{align}
  A_k\sur{(meas)}(t) = \tilde{A}_k(t){\mathrm{e}}^{-{\mathrm{i}} \omega\sub{LO}(t + \tau\sub{R}k/N)} = A_k(t){\mathrm{e}}^{{\mathrm{i}} (\omega\sub{s} - \omega\sub{LO})(t + \tau\sub{R}k/N)}.  \label{eq:3}
\end{align}
To obtain the slowly varying amplitude $A_k(t)$, we need to compensate the phase factor
coming from the frequency difference of $\omega\sub{s}$ and $\omega\sub{LO}$.

We estimated the frequency difference of $\omega\sub{s}$ and $\omega\sub{LO}$ from the results of two-round-trip data of in-phase/quadrature-phase measurements.
In our experiments, the timescale of the phase dynamics is about the order of $10\,\U{us}$ because the fastest phase dynamics is determined by $\gamma\sub{inj} = 72\,\U{kHz}$.
Since the round-trip time $\tau\sub{R} = 101\,\U{ns}$ is much shorter than the timescale of the phase dynamics, we may assume that the slowly varying amplitude $A_k(t)$ is not changed significantly after one round trip of the fiber cavity:
\begin{align}
  A_k(t + \tau\sub{R}) \simeq A_k(t)  \label{eq:4}
\end{align}
Thus, the measured amplitude after one round trip can be written as
\begin{align}
  A_k\sur{(meas)}(t + \tau\sub{R}) \simeq A_k(t){\mathrm{e}}^{{\mathrm{i}}(\omega\sub{s} - \omega\sub{LO})t}{\mathrm{e}}^{{\mathrm{i}} (1 + k/N) \Delta \phi}  \label{eq:5}
\end{align}
where $\Delta\phi := (\omega\sub{s} - \omega\sub{LO})\tau\sub{R}$.
From Eq.~\eqref{eq:2}, the phase difference after one-round trip satisfies
$|\Delta \phi| < \pi$.

We can estimate the phase difference $\Delta \phi$ by comparing the measured amplitude over two round trips:
\begin{align}
  A_k\sur{(meas)}(t + \tau\sub{R})A_k\sur{(meas)}(t)^* = |A_k(t)|^2 {\mathrm{e}}^{{\mathrm{i}} \Delta\phi}  \label{eq:6}
\end{align}
Taking the angle of the summation of the inner products gives us the estimated value as
\begin{align}
  \Delta \phi\sur{(est)} &:= \mathrm{arg}
  \left[
    \sum_{k=1}^N A_k\sur{(meas)}(t + \tau\sub{R})A_k\sur{(meas)}(t)^*
  \right] \nonumber \\
  &= \mathrm{arg} \left[
    {\mathrm{e}}^{{\mathrm{i}}\Delta\phi} \sum_{i=1}^N |A_k(t)|^2 {\mathrm{d}} t
  \right]
  = \Delta \phi  \label{eq:7},
\end{align}
The absolute phases of the pulses were obtained by compensating the frequency difference of $\omega\sub{s}$ and $\omega\sub{LO}$ by using $\Delta \phi\sur{(est)}$ as
\begin{align}
  A_k(t) = A_k\sur{(meas)}(t){\mathrm{e}}^{-{\mathrm{i}} \Delta \phi\sur{(est)} k / N}{\mathrm{e}}^{-{\mathrm{i}} \Delta \phi\sur{(est)} t/\tau\sub{R} }.  \label{eq:8}
\end{align}
Since the second phase factor ${\mathrm{e}}^{-{\mathrm{i}} \Delta \phi\sur{(est)} t/\tau\sub{R}}$ is common for all pulses at the same time $t$, only the first phase factor ${\mathrm{e}}^{-{\mathrm{i}} \Delta \phi\sur{(est)} k / N}$ was compensated in our experiments.

\subsection{Relative phase measurement of adjacent pulses}

\begin{figure}
  \centering
  \includegraphics[width=400pt]{./rel_phase_meas.pdf}
  \caption{One-interval delay IQ interferometer for relative phase measurements.
  The input pulses (shown by the red line) are aligned so as to have the horizontal (H) polarization
  and pass through the first polarizing beam splitter (PBS). Then the pulses are
  split into two beams by the following non-polarizing beam splitter (NBS).
  The two beams travel through two interferometer arms that have different optical path lengths. The beam travelling trough the longer arm is delayed by one pulse interval, and its polarization is rotated into the vertical (V) direction with a half-wave plate.
  The two beams are recombined with the second polarizing beam splitter. Then the relative phase between adjacent pulses is mapped onto the relative phase between horizontal and vertical polarization in the single pulse.
  The converted pulses are introduced into a single-mode fiber and are split by the 50:50 coupler. The polarization of those pulses are tuned by polarization controllers (PCs), and their diagonal and circular polarization components are measured by the polarization beam splitters (PBSs) followed by the balanced photodetectors (BPDs).
  The diagonal and circular polarization components of the converted pulses correspond to the in-phase and quadrature phase components of the relative phases between adjacent pulses, respectively.
  The free-space interferometer is stabilized with the CW laser (shown by the blue lines), which is the same laser as that used for stabilizing the optical path length of $\pm 1$-interval delay lines for the mutual injection. The CW laser with diagonal polarization is inputted from the other output port of the second polarization beam splitter. The polarization of the CW laser is converted into vertical polarization in both arms and the two beams interfere when they are combined with the non-polarization beam splitter. The output beam is reflected at the first polarization beam splitter and is measured by the photodetector. The interference signal is fed back to the piezoelectric actuator driving the mirror placed in the longer arm of the free-space interferometer.}
  \label{fig:rel_phase_meas}
\end{figure}

The relative phases of adjacent pulses were measured by using a one-interval delay IQ interferometer.
The experimental setup is shown in Fig.~\ref{fig:rel_phase_meas}.

The one-interval delay IQ interferometer is composed of two parts. 
The first part consists of the free-space delay-line interferometer.
This part splits the pulses into two beams and delays one of them by one pulse interval. Furthermore, the polarization of the beam is rotated into the orthogonal polarization in the longer arm. Then two beams are recombined.
As a result, the relative phase between the two adjacent pulses is converted into the relative phase between the horizontal and vertical polarization in the single pulse.
The length of the delay line interferometer is stabilized with the CW laser which travels through the interferometer in the opposite direction to the measured pulses.

The second part is composed of fiber optics, and used for the measurements of the in-phase and quadrature phase components of the relative phase difference between horizontal and vertical polarization. 
The input pulses are coupled into two fibers with a 50:50 fiber beam splitter (BS).
One of the fiber outputs is used for the in-phase measurement and the other is used for the quadrature-phase measurement.
In each fiber, the polarization of the pulses is controlled by the polarization controllers (PCs) so as to be measured in a proper basis by the following polarizing beam splitters (PBSs) and the balanced photodetector (BPDs).

The in-phase and quadrature phase components of the one-interval delay IQ measurements are written as
\begin{align}
  I_i &= \mathop{\mathrm{Re}}[A_iA_{i+1}^*] \simeq n\sub{s}\cos(\theta_i - \theta_{i+1})  \label{eq:9} \\
  Q_i &= \mathop{\mathrm{Im}}[A_iA_{i+1}^*] \simeq n\sub{s}\sin(\theta_i - \theta_{i+1})  \label{eq:10}
\end{align}
Thus, the summation over one-round trip of the in-phase components 
is proportional to the Hamiltonian of the one-dimensional XY ring:
\begin{align}
  \sum_{i=1}^N I_i = n\sub{s}\sum_{i=1}^N \cos(\theta_i - \theta_{i+1}).  \label{eq:11}
\end{align}
In our experiments, this value was directly measured by averaging the output over about five round trips with a low-pass filter (cutoff frequency: $1.9\, \U{MHz}$).

\section{Equilibration time of one-dimensional XY model}

In our experiments, the effective temperature estimated from the winding number distribution was much different from the expected temperature determined from the
ratio between injection and diffusion.
This difference comes from the slow equilibration of the Langevin dynamics for the one-dimensional XY model.
In this section, we numerically analyze the equilibration of the winding number distribution for various inverse temperature.

To evaluate the long-term dynamics, we further simplified the Langevin dynamics of the coupled lasers so that only phases of lasers are treated as dynamical variables.
Assuming all lasers have the same steady-state photon number $n\sub{s}$, the dynamics of phases of lasers are described by the following Langevin equations:
\begin{align}
  \frac{{\mathrm{d}} \theta_i}{{\mathrm{d}} t} &= -\frac{\gamma\sub{inj}}{2}\frac{\partial H}{\partial \theta_i} + \sqrt{D_\theta} \eta(t),  \label{eq:12} \\
  H(\vct{\theta}) &= - \sum_{i < j} J_{ij} \cos(\theta_i - \theta_j),  \label{eq:13}
\end{align}
where $D_\theta = D/n\sub{s}$ and $\eta(t)$ is a white noise satisfying $\bracket{\eta(t)\eta(t')} = \delta(t-t')$.

In our numerical simulation, the number of spins is set as $N=100$ and the Hamiltonian~\eqref{eq:13} is chosen to be the one-dimensional ferromagnetic ring with $J_{i, i+1} = 1$ for all $i$.
The phase diffusion constant is set as $D_\theta = 0.5 \,\U{kHz}$.
The injection ratio is set as $\gamma\sub{inj} = \beta D_\theta$ depending on the inverse temperature $\beta$ to be simulated.
We used the Euler-Maruyama method to simulate Eq.~\eqref{eq:12}.
The simulated range of $\beta$ was from $\beta = 1$ to $\beta = 10$.
The numerical simulation was repeated over 1,000 runs for each $\beta$.


\begin{figure}
  \centering
  \includegraphics[width=480pt]{./w_hist_time_series.pdf}
  \caption{Time evolution of the winding number distribution for $\beta = 5$.
  (a) Probability distribution of winding number at each running time.
  The horizontal axis shows the running time $t$ plotted in a log scale.
  The color bar represents the probability of obtaining the corresponding winding number. (b) Probability distribution of winding number at $t=0.01, 1, 100 \,\U{s}$.}
  \label{fig:w_hist_time_series}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=350pt]{./beta_vs_beta_est.pdf}
  \caption{Running time dependence of the estimated inverse temperature.
  Each solid line represents the relationship between the simulated inverse temperature $\beta = \gamma\sub{inj}/D_\theta$ and the estimated inverse temperature for the respective running time. The dashed line represents $\beta\sur{(est)} = \beta$.}
  \label{fig:beta_vs_beta_est}
\end{figure}
Figure~\ref{fig:w_hist_time_series} shows the time evolution of the winding number distribution for $\beta = 5$, as an example.
The winding number distribution is first spread over broad range and gradually converges to a sharp distribution around $m = 0$.
The equilibration time for this case is around $t = 1\,\U{s}$.
The inverse temperature was estimated to be $\beta\sur{(est)} = 4.64 \pm 0.06$ at $t = 100 \,\U{s}$.

Figure~\ref{fig:beta_vs_beta_est} shows the estimated inverse temperature $\beta\sur{(est)}$ for various simulated inverse temperatures $\beta$ and running time $t$.
The estimated inverse temperature approaches the line $\beta\sur{(est)} = \beta$
as the running time becomes larger.
The slight deviation from the line $\beta\sur{(est)} = \beta$ possibly comes from the inaccuracy of the estimation method.
We assumed the states with winding number $m$ always have the same energy as that of the local minimum $E_m = - N \cos(2\pi m/N)$. We then estimated the inverse temperature by fitting the distribution with $p(m) \propto \exp(-\beta\sur{(est)} E_m)$.
However, the actually sampled configurations deviate slightly from the configuration of the local minimum and
the energies also deviate from the assumed one. This made the estimation slightly inaccurate.

From Fig.~\ref{fig:beta_vs_beta_est}, the distribution for larger inverse temperature $\beta$ takes a longer running time to reach the true equilibrium. The equilibration time seems to depends exponentially on $\beta$.
Even when the running time $t = 100\,\U{s}$, we can obtain the true steady-state distribution for up to around $\beta = 7$.

\section{Theory of one-dimensional XY model}

Various statistical features of the one-dimensional XY ring can be exactly calculated by using transfer matrix approach \cite{Mattis1984}.
We describe the way to calculate the partition function, the correlation function,
and the probability distribution of the relative angle between adjacent phases.

\subsection{Partition function}

The partition function $Z$ of the one-dimensional ferromagnetic XY ring is given by
\begin{align}
  Z &= \mathrm{Tr}[{\mathrm{e}}^{- \beta H}],  \label{eq:14} \\
  H(\theta) &= J \sum_{i = 1}^N \cos(\theta_i - \theta_{i+1}),  \label{eq:15}
\end{align}
with the periodic boundary condition $\theta_{N+1} = \theta_1$.

Set $K = \beta J$ and define the matrix
\begin{align}
  \hat{V} = \frac{1}{2\pi}\oint \oint {\mathrm{d}} \theta' {\mathrm{d}} \theta \exp[K\cos(\theta' - \theta)] \ket{\theta'}\bra{\theta},  \label{eq:16}
\end{align}
where $\bracket{\theta|\theta'} = \delta(\theta - \theta')$.
The partition function of the one-dimensional ferromagnetic XY ring can be written as
\begin{align}
  Z = \mathrm{Tr}[\hat{V}^N].  \label{eq:17}
\end{align}
We can use the following expansion to diagonalize the matrix $\hat{V}$:
\begin{align}
  {\mathrm{e}}^{K \cos\theta} = \sum_{n = -\infty}^{\infty} I_n(K){\mathrm{e}}^{{\mathrm{i}} n \theta},  \label{eq:18}
\end{align}
where $I_n$ is the modified Bessel function of the first kind.
Define the basis vector $\ket{n}$ as
\begin{align}
   \ket{n} = \frac{1}{(2\pi)^{1/2}} \oint {\mathrm{e}}^{{\mathrm{i}} n \theta} \ket{\theta} {\mathrm{d}} \theta ,  \label{eq:19}
\end{align}
then we have
\begin{align}
  \hat{V} = \sum_{n=-\infty}^{\infty} I_n(K) \ket{n}\bra{n}  \label{eq:20}
\end{align}
Thus, the partition function can be given by
\begin{align}
  Z = \sum_{n=-\infty}^{\infty} I_n(K)^N.  \label{eq:21}
\end{align}

\subsection{Correlation function}

Similar to the calculation of the partition function,
the correlation function can be given by the following form:
\begin{align}
  \bracket{\exp({\mathrm{i}} m (\theta_1 - \theta_{k+1}))}
  = \frac{1}{Z}\mathrm{Tr} [\hat{W}_m^{k} \hat{V}^{N-k}],  \label{eq:22}
\end{align}
where
\begin{align}
  \hat{W}_m = \frac{1}{2\pi}\oint \oint {\mathrm{d}} \theta' {\mathrm{d}} \theta \exp({\mathrm{i}} m (\theta' - \theta)) \exp[K\cos(\theta' - \theta)] \ket{\theta'}\bra{\theta}.  \label{eq:23}
\end{align}
and $m$ is an integer number.
Using the basis ${\ket{n}}$, the matrix $\hat{W}_m$ can be expressed as
\begin{align}
  \hat{W}_m = \sum_{n= -\infty}^{\infty} I_{n-m}(K) \ket{n}\bra{n}.  \label{eq:24}
\end{align}
Thus, we have
\begin{align}
  \bracket{\exp({\mathrm{i}} m (\theta_1 - \theta_{k+1}))}
  = \frac{1}{Z} \sum_{n=-\infty}^{\infty} I_{n-m}(K)^k I_n(K)^{N-k}.  \label{eq:25}
\end{align}
Since the right-hand side of the equation is real-valued, we have
\begin{align}
  \bracket{\cos(m(\theta_1 - \theta_{k+1}))} 
  = \frac{1}{Z} \sum_{n=-\infty}^{\infty} I_{n-m}(K)^k I_n(K)^{N-k}.  \label{eq:26}
\end{align}

\subsection{Probability distribution of relative angle}

The probability distribution of the relative angle between adjacent XY spins can be given as
\begin{align}
  p(\theta) &= \bracket{\delta(\theta_2 - \theta_1 - \theta)} \nonumber \\
  &= \frac{1}{Z}\frac{1}{(2\pi)^N} \oint {\mathrm{d}} \theta_1 \cdots \oint {\mathrm{d}} \theta_N \delta(\theta_2 - \theta_1- \theta) \prod_{i=1}^N \exp [K \cos(\theta_i - \theta_{i+1}) ] \nonumber \\
  &= \frac{\exp(K\cos\theta)}{Z} \frac{1}{(2\pi)^N} \oint {\mathrm{d}} \theta_2 \cdots \oint {\mathrm{d}} \theta_N 
   \exp[K\cos(\theta_N - \theta_2 + \theta)]
  \prod_{i=2}^{N-1} \exp [K \cos(\theta_i - \theta_{i+1}) ].  \label{eq:27}
\end{align}
Define the following matrix
\begin{align}
  \hat{X}(\theta) = \frac{1}{2\pi}\oint \oint {\mathrm{d}} \theta' {\mathrm{d}} \theta'' \exp[K\cos(\theta' - \theta'' + \theta)] \ket{\theta'}\bra{\theta''},  \label{eq:28}  
\end{align}
Then $p(\theta)$ can be written as
\begin{align}
  p(\theta) &= \frac{\exp(K\cos\theta)}{2\pi Z}\mathrm{Tr}[\hat{V}^{N-2}\hat{X}(\theta)]  \label{eq:29}
\end{align}
The matrix $\hat{X}(\theta)$ can be diagonalized as
\begin{align}
  \hat{X}(\theta)
   = \sum_{n=-\infty}^{\infty} I_n(K) {\mathrm{e}}^{{\mathrm{i}} n \theta} \ket{n}\bra{n},  \label{eq:30}
\end{align}
and then we have
\begin{align}
  p(\theta) = \frac{\exp(K\cos\theta)}{2\pi Z}\sum_{n=-\infty}^{\infty} {\mathrm{e}}^{{\mathrm{i}} n \theta} I_n(K)^{N-1}.  \label{eq:31}
\end{align}
Using the relationship $I_n(K) = I_{-n}(K)$, finally we obtain
\begin{align}
  p(\theta) = \frac{\exp(K\cos\theta)}{2\pi Z}\sum_{n=-\infty}^{\infty} \cos(n\theta) I_n(K)^{N-1}.  \label{eq:32}
\end{align}

