\section{Introduction}\label{sec:intro}

Precise measurement of a sufficiently broad range of low temperatures is one of the primary challenges in the practical implementation of quantum technologies~\cite{mehboudi2019thermometry,giovannetti2011advances,PhysRevLett.96.010401,toth2014quantum,celi2016quantum,RevModPhys.80.885,bloch2012quantum}. Quantum thermometry is an emerging field of quantum metrology, which can offer promising solutions to the low temperature measurement problem in terms of quantum thermometers~\cite{de2016local,PhysRevLett.122.030403,PhysRevB.98.045101,PhysRevX.10.041054,olf2015thermometry,PhysRevLett.96.130404,PhysRevA.103.023317,PhysRevLett.119.090603,PhysRevLett.125.080402,latune2020collective,PhysRevLett.127.190402,potts2019fundamental,mehboudi2015thermometry,PhysRevA.82.011611}. An ideal quantum thermometer should be much smaller relative to the sample, and its coupling to the sample should not significantly alter the temperature of the sample. The smallest thermometer can be a quantum two-level system (qubit), which has attracted much attention recently~\cite{PhysRevA.86.012125,PhysRevA.84.032105,bloch2012quantum,PhysRevA.91.012331,PhysRevLett.118.130502,PhysRevA.98.042124,PhysRevA.99.062114,PhysRevResearch.2.033497, PhysRevResearch.2.033394}. A qubit thermometer for ultracold gases was recently demonstrated~\cite{PhysRevX.10.011018}. The increase in the dimensionality of the probe does not significantly improve the thermal sensitivity of the thermometer~\cite{campbell2018precision}. An $N$-level optimal probe at thermal equilibrium must have $N-1$-fold degeneracy in the excited state, making it effectively a two-level system ~\cite{PhysRevLett.114.220405}. Such an optimal probe, however, exhibits a high precision estimation for only a specific temperature, which is indicated by a single peak in the quantum Fisher information.

Finding a specific probe that can precisely measure a broad range of low temperatures is among the central challenges in thermometry~\cite{mok2021optimal}. A probe's range of thermal sensitivity can also be enhanced by considering strong system-bath coupling~\cite{PhysRevA.96.062103, PhysRevB.98.045101,glatthard2022bending}. However, the quantum Fisher information still shows a single peak, indicating the probe remains an optimal thermal sensor for a single temperature. Broader range thermal precision measurements require extremely large degeneracy in the higher excited states of the probe~\cite{campbell2018precision}, which is challenging for physical implementations~\cite{mok2021optimal}. Alternative proposals suggest applying external periodic control on a probe to get multiple peaks in the quantum Fisher information at the cost of making the thermometer non-autonomous~\cite{Mukherjee2019,glatthard2022bending}.

Studies in quantum thermodynamics~\cite{e15062100,binder2018thermodynamics,campbellbook,AsliReview,kurizki_kofman_2022} revealed deeper insight into the interplay between heat, work, quantum information, and in particular, quantum coherence~\cite{science.1078955,pnas.1110234108,PhysRevLett.112.030602,PhysRevLett.113.150402,Lostaglio2015,PhysRevX.5.021001,PhysRevX.5.031044,Hardal2015,Kammerlander2016,Latune2019,Kim2022}. In parallel with the quantum thermodynamical perspectives of quantum coherence, several schemes have been proposed in quantum thermometry to exploit nonclassical features to enhance thermal sensitivity~\cite{PhysRevA.101.032112,PhysRevLett.125.080402,PhysRevLett.123.180602,PhysRevA.96.062103,potts2019fundamental,Campbell2017,PhysRevA.96.012316,PhysRevA.99.062114,PhysRevA.91.012331,PhysRevLett.118.130502,PhysRevA.98.042124,PhysRevResearch.2.033497,PhysRevA.98.050101}. Most of these studies require a dynamical approach in which temperature is measured in transient regime~\cite{PhysRevA.98.042124,PhysRevResearch.2.033497}. This restriction fundamentally exists because, the probe eventually thermalizes in the steady-state, which is a diagonal (Gibbs) state without any coherence, consequently losing all nonclassical features and quantum advantages. For a probe in thermal equilibrium, the temperature information is encoded only in populations, and quantum Fisher information becomes equal to classical Fisher information~\cite{PhysRevLett.114.220405,Paris_2015}. Moreover, previous studies do not consider coherence or entanglement to enhance the thermal sensitivity range. The quantum Fisher information shows only a single peak in these proposals~\cite{PhysRevA.91.012331,PhysRevLett.118.130502,PhysRevA.98.042124,PhysRevResearch.2.033497}.

We propose a precise and wide-range low-temperature measurement scheme based on a simple thermalization process of a multi-qubit system with a sample in thermal equilibrium. The measurement is made on a single probe qubit of the multi-qubit system to benefit from the quantum coherent character of the probe qubit. The sample has no direct access to the probe qubit and only interacts with the rest of the qubits (ancilla qubits). Ancilla qubits interface the probe and the sample, and their task is to generate coherence in the probe qubit in the steady state depending on the temperature of the sample. We analyze the open system dynamics of the ancilla-probe system attached to the sample (thermal environment) by deriving a global Markovian master equation~\cite{Levy_2014}. As opposed to the phenomenological local master equation, the global approach is consistent with the laws of thermodynamics within the validity of Secular approximation~\cite{Levy_2014, PhysRevA.98.052123, Cattaneo_2019}.  In the global master equation, the jump operators which describe the dissipation processes are non-local. Consequently, the sample can indirectly access the probe qubit. As such, the sample can encode its temperature information into the diagonal (population) elements of the probe's density matrix via the non-local dissipation channels. In addition, temperature information is further translated into the off-diagonal elements (coherences) of the probe state due to the direct coherent interaction between the ancilla and probe qubits. By judicious choice of system parameters, temperature-dependent coherence leads to the emergence of multiple peaks in the quantum Fisher information, suggesting the high precision estimation beyond a single temperature. Furthermore, the presence of ancilla qubits enhances the probe's thermal sensitivity at very low temperatures consistent with the dynamical ancilla assisted quantum thermometry for a single-peak enhanced quantum Fisher information~\cite{PhysRevA.98.042124}.

Compared to previous proposals, we summarize the distinct features of our scheme: (i) Our ancilla interfaced qubit thermometer has a high sensitivity for more than a single temperature, which is indicated by multiple peaks in the quantum Fisher information, which typically require highly degenerate multi-level probes~\cite{campbell2018precision}. (ii) Our scheme is based upon the complete thermalization of a multi-qubit system with a thermal sample; hence, it is autonomous in contrast to other proposals for wide-range quantum thermometers~\cite{Mukherjee2019,glatthard2022bending}. (iii) Unlike previous proposals, initial coherence in the probe~\cite{PhysRevA.98.042124}, or initially entangled state~\cite{PhysRevA.91.012331, Tham2016} are not needed for the enhancement of thermal sensitivity. (iv) Our results are robust against the change in bath spectra such as Ohmic and super (sub)-Ohmic, which was not the case in the previous proposal for a multi-peaked quantum Fisher information probe~\cite{Mukherjee2019}.

The rest of the paper is organized as follows: In Sec.~\ref{sec:II}, we introduce our model of quantum thermometer and derive the global master equation. In addition, we briefly discuss quantum Fisher information (QFI), a figure of merit in quantum metrology.
In Sec.~\ref{sec III}, we present and discuss the QFI for the case of single and two ancilla qubits. The concluding remarks are given in Sec.~\ref{sec IV}.

	\section{Model and Preliminaries}\label{sec:II}
 In this section, we introduce our theoretical model, which can boost the thermal sensitivity of a two-level system probe at low temperatures via coherence generation. Then, we briefly discuss QFI, a key figure of merit for the thermal sensitivity of a quantum thermometer.

	\subsection{Model description}\label{sec:IIA}

We consider a qubit-probe of transition frequency $\omega_{p}$ coupled to $N$ ancilla qubits of frequency $\omega_{k}$ ($k= 1, 2,..., N$). The ancilla qubits are coupled directly to a thermal bath of temperature $T$ whose temperature is to be estimated. The qubit probe is outside the thermal bath and does not have direct access to it, as shown in Fig.~\ref{fig1}. The Hamiltonian of the ancilla-probe system is given by~\cite{PhysRevA.98.042124, PhysRevE.99.042121, Naseem_2020, PhysRevResearch.2.033285} (for simplicity, we take $\hbar=1$)
\begin{equation}
  \hat{H}=\frac{1}{2}\omega_p\hat{\sigma}^z_p+\sum\limits_{k=1}^{N}\frac{1}{2}\omega_k\hat{\sigma}^z_k+\sum\limits_{k=1}^{N}g_k\hat{\sigma}^z_k\hat{\sigma}^x_p,\label{eq1}
\end{equation}
here $g_k$ being the coupling strengths between the ancilla and probe, and
$\hat{\sigma}^\zeta_{j}$ ($j=p, k$) with $\zeta_=x,y,z$ are the Pauli matrices. In Eq.~(\ref{eq1}), the first, second, and third terms denote the Hamiltonians of the probe, ancilla, and their interaction, respectively. We note that the asymmetric interaction between the qubits given in Eq.~(\ref{eq1}) has previously been employed for the proposals of ancilla-assisted thermometry~\cite{PhysRevA.98.042124}, quantum thermal diode~\cite{PhysRevE.99.042121}, quantum absorption refrigerator~\cite{Naseem_2020}, and quantum heat manager~\cite{PhysRevResearch.2.033285}.

We consider asymmetric interaction between the ancilla and probe qubits, which leads to temperature-dependent coherences in the probe (see Eq.~(\ref{eq:coherences})). Such temperature-dependent coherences lead to the enhancement of the thermal sensitivity of the probe (Fig.~\ref{fig2}). We note that initial coherences or coherent driving on the probe to induce coherence do not enhance the thermal sensitivity of the probe because the coherences, in this case, are temperature independent~\cite{PhysRevLett.114.220405}. In addition, other symmetric or anti-symmetric interactions between the ancilla and probe qubits may not generate coherence in the probe. For instance, we have shown in Appendices~\ref{app:A} and \ref{app:B} that dipole-dipole interaction and Dzyaloshinskii-Moriya interaction do not generate probe coherence, respectively.

  \begin{figure}[t]
  	\centering
  	\subfloat{
  		\includegraphics[scale=0.3]{model.pdf}}
  	\caption{A schematic illustration of our quantum thermometry model. The ancilla qubits of transition frequencies $\omega_k$ are coupled to a thermal bath of temperature $T$. The probe-qubit of transition frequency $\omega_p$ has no direct access to the thermal bath; instead, it is coupled with the ancilla qubits and indirectly probes the temperature of the thermal bath.
 It is important to note that the qubit outside the bath is our quantum probe for the estimation of $T$, and the energy measurements are only performed over it.}\label{fig1}
  \end{figure}

To investigate the dynamics of the ancilla-probe system, we employ the standard Born-Markov master equation~\cite{Breuer}. The derivation of the master equation for coupled systems requires diagonalization of the total system Hamiltonian. The Hamiltonian of the ancilla-probe system can be diagonalized using the unitary transformation~\cite{PhysRevE.99.042121}
\begin{equation}\label{eq:transformation}
U :=\text{exp}[-\frac{i}{2}\sum\limits_{k=1}^{N}\theta_k\hat{\sigma}^z_k\hat{\sigma}^y_p],
\end{equation}
this gives the diagonalized Hamiltonian in the form
\begin{equation}
\tilde{H}=\sum\limits_{k=1}^{N}\frac{\omega_k}{2}\tilde{\sigma}^z_k+\frac{\Omega}{2}\tilde{\sigma}^z_p.\label{eq:DiagHamil}
\end{equation}
The transformed frequency $\Omega$ of the probe-qubit depends on the number of ancilla qubits, and for a single and two ancilla qubits, the explicit expressions are given in Eqs.~(\ref{eq:singletrans}), and ~(\ref{eq:omega2}), respectively.
After the transformation, the Pauli matrices of the ancilla and probe qubits take the form
\begin{eqnarray}\label{eq:Paulitrans}
\tilde{\sigma}^x_k &=& \cos\theta_k\hat{\sigma}^x_k+\sin\theta_k\hat{\sigma}^y_k\hat{\sigma}^y_p, \nonumber \\
\tilde{\sigma}^y_k &=& \cos\theta_k\hat{\sigma}^y_p-\sin\theta_k\hat{\sigma}^x_k\hat{\sigma}^y_p, \quad\quad \tilde{\sigma}^z_k = \hat{\sigma}^z_k \nonumber \\
\tilde{\sigma}^x_p &=& \cos\theta_k\hat{\sigma}^x_p-\sin\theta_k\hat{\sigma}^z_k\hat{\sigma}^z_p, \quad\quad \tilde{\sigma}^y_p=\hat{\sigma}^y_p, \nonumber\\
\tilde{\sigma}^z_p &=&\cos\theta_k\hat{\sigma}^z_p+\sin\theta_k\hat{\sigma}^z_k\hat{\sigma}^x_p.
\end{eqnarray}

As an example, we assume a bosonic thermal bath---whose temperature is to be probed---coupled to ancilla qubits. The interaction Hamiltonian of ancilla and bath has the form
\begin{equation}\label{eq:int}
	\hat{H}_\text{int} = \sum_{l,k} g_{l,k}\hat{\sigma}^{x}_{k}(\hat{a}_{k}+\hat{a}^{\dagger}_{k}),
\end{equation}
where, $\hat{a}_{l}$ ($\hat{a}^{\dagger}_{l}$) is the annihilation (creation) bosonic operator of the  $l$-th bath mode, and $g_{k,l}$ indicates the coupling strength between the $l$-th bath mode with $k$-th ancilla qubit. The master equation for the coupled ancilla-probe system after the Born-Markov and Secular approximations is given by~\cite{PhysRevE.99.042121, Naseem_2020, PhysRevResearch.2.033285} (here for notational simplicity, we use the summation convention to drop the summation sign such that summation over repeated index $k$ is implied.)
	\begin{equation}
\begin{aligned}
\dot{\tilde{\rho}} =&-i[\tilde{H},\tilde{\rho}]+\cos^2\theta_k \mathcal{G}(\omega_k)\Big(\mathcal{D}[\tilde{\sigma}^-_k]+e^{-\omega_k/T}\mathcal{D}[\tilde{\sigma}^+_k]\Big)\\
& +\sin^2\theta_k \mathcal{G}(\omega_{-k})\Big(\mathcal{D}[\tilde{\sigma}^-_k\tilde{\sigma}^+_p]+e^{-\omega_{-k}/T}\mathcal{D}[\tilde{\sigma}^+_k\tilde{\sigma}^-_p]\Big)\\
&
+\sin^2\theta_k \mathcal{G}(\omega_{+k})\Big(\mathcal{D}[\tilde{\sigma}^-_k\tilde{\sigma}^-_p]+e^{-\omega_{+k}/T}\mathcal{D}[\tilde{\sigma}^+_k\tilde{\sigma}^+_p]\Big)\label{eq:master}.
\end{aligned}
	\end{equation}
Here, $\omega_{\pm k} = \omega_{k}\pm\Omega$, and mixing angle $\theta_{k}$ depends on the number of ancilla qubits, for a single and two ancilla qubits $\theta_{k}$ is given in Eqs.~(\ref{eq:singletrans}), and ~(\ref{eq:angle2}), respectively. In Eq.~(\ref{eq:master}), $\mathcal{D}[\hat{c}]$ is the Lindblad super-operator and is defined as
\begin{equation}
\mathcal{D}[\tilde{c}]=\tilde{c}\tilde{\rho}\tilde{c}^\dagger-\frac{1}{2}(\tilde{c}^\dagger \tilde{c}\tilde{\rho}+\tilde{\rho}\tilde{c}^\dagger \tilde{c}).
\end{equation}
The raising (lowering) spin operators $\tilde{\sigma}^+(\tilde{\sigma}^-)$ are defined as
\begin{equation}
\tilde{\sigma}^\pm_{k,p}=\frac{1}{2}(\tilde{\sigma}^x_{k,p}\pm i\tilde{\sigma}^y_{k,p}).
\end{equation}
The bath spectral response function is denoted by $\mathcal{G}(\omega_k)$ in Eq.~(\ref{eq:master}).
An appealing feature of our scheme is that the thermal sensitivity of the probe does not depend on the bath spectrum. The ancilla-probe system evolves to a steady state independent of the bath spectrum (see Sec.~\ref{sec III} for details). In the master equation (\ref{eq:master}), we have neglected cross dissipators of the form $\mathcal{D}[\tilde{\sigma}_{k},\tilde{\sigma}_{k^{\prime}}]$; this is well-justified if the distance between the ancilla qubits is much larger compared to the wavelength of the bath quanta (for example photons)~\cite{Damanet_2016,PhysRevA.99.052105,PhysRevE.102.042111}.

\subsection{Quantum Fisher information} \label{sec IIB}
In quantum thermometry, one requires a quantum system for the estimation of the unknown temperature of a sample. A small thermometer is desired compared to the sample so that it should not significantly affect the temperature of the sample. In this respect, the smallest possible thermometer can be a two-level quantum system (such as a qubit) which is small yet efficient as a probe for estimation purposes \cite{bloch2012quantum, PhysRevA.91.012331,PhysRevLett.118.130502,PhysRevA.98.042124,PhysRevA.99.062114,PhysRevResearch.2.033497, PhysRevResearch.2.033394}. The qubit probe initially prepared in a state $\hat{\rho}_p$ allows interaction with a bath whose temperature is to be estimated. The bath information is imprinted onto the probe state to be mapped via some estimator $\hat{T}$ \cite{paris2009quantum}. Fundamental bounds limit the precision of parameter estimation. The optimal estimators in quantum estimation theory are those saturating the quantum Cram\'er-Rao bound inequality such that \cite{cramer1999mathematical,helstrom1969quantum,PhysRevLett.72.3439}
\begin{equation}
\Delta T^2\geq\frac{1}{\mathcal{N}\mathcal{F}_Q(T)}.
\end{equation}
Here, $\Delta T^2$ is the temperature variance, and $\mathcal{N}$ denotes the number of measurements. $\mathcal{F}_Q(T)$ is the quantum Fisher information which quantifies the information about temperature $T$ encoded onto the state of the probe, and it is defined as
\begin{equation}
\mathcal{F}_Q(T)=\text{Tr}[\rho(T) L_\text{T}].
\end{equation}
$L_\text{T}$ being symmetric logarithmic derivative satisfying the equation $2\partial_\text{T}\rho(T)=\{L_\text{T},\rho(T)\}$. As depicted in Fig.~(\ref{fig1}), our quantum thermometer is a single qubit probe interacting indirectly with the thermal bath via an interface of ancilla qubits. Accordingly, we will only focus on the QFI of a single qubit, which is given by~\cite{dittmann1999explicit,dittmann1993riemannian}
\begin{equation}\label{eq:TLSfisher}
\mathcal{F}_Q(\hat{\rho}_p)=\text{Tr}\Big[\Big(\frac{\partial \hat{\rho}_p}{\partial T}\Big)^2\Big]+\frac{1}{\text{Det}(\hat{\rho}_p)}\text{Tr}\Big[\hat{\rho}_p\Big(\frac{\partial\hat{\rho}_p}{\partial T}\Big)^2\Big].
\end{equation}
Mathematical basis of our scheme can be seen in this expression, which can be used for both thermal and non-thermal qubit states. If there are any bath temperature dependent coherences in the probe qubit state, they can be used to induce multiple peaks and enhancements in QFI.

\section{results}\label{sec III}
In this section, we will present results for the low-temperature quantum thermometry for the model described in Fig.~\ref{fig1} by considering two (Sec.~\ref{sec IIIA}) and three (Sec.~\ref{sec IIIB}) ancilla qubits. We show that multiple peaks in the quantum Fisher information can be obtained by harvesting quantum coherence in the probe induced by ancilla qubits.

\subsection{Enhancement of low-T sensitivity with a single ancilla qubit ($N=1$)}\label{sec IIIA}

In the case of a single ancilla qubit, there is no summation and $k=1$ in the master equation (\ref{eq:master}). The diagonalized Hamiltonian of the ancilla-qubit system $\tilde{H}$ is given in Eq.~(\ref{eq:DiagHamil}). The transformed frequency of the probe qubit and the mixing angle $\theta_{1}$ are given by
	\begin{equation}\label{eq:singletrans}
     \Omega = \sqrt{\omega^{2}_{1}+4g^2_{1}},\quad  \text{and}\quad\theta_{1} = \text{arctan}\bigg(\frac{2g_1}{\omega_p}\bigg),
	\end{equation}
respectively. The second term in the master equation~(\ref{eq:master}) describes the energy exchange between the ancilla and the thermal bath via the local dissipation channel. The last two terms, on the contrary, correspond to the non-local dissipation channels between the bath and ancilla-probe system. These two terms show that the probe qubit can exchange energy with the bath via non-local dissipation processes, though it is not directly coupled with the bath. The evolution of the ancilla-qubit system by the master equation~(\ref{eq:master}) leads to a thermal equilibrium state of the ancilla-qubit system, and at steady-state, which is given by
\begin{equation}
\tilde{\rho}_\text{ss}=\tilde{\rho}_1\otimes\tilde{\rho}_p.
\end{equation}
Here $\tilde{\rho}_{1}$ and $\tilde{\rho}_{p}$ are the reduced states of the ancilla and probe qubits, which are given by
\begin{equation}
	\hat{\rho}_1=\frac{e^{-\frac{1}{2}\beta\omega_{1}\tilde{\sigma}^{z}_{1}}}{\mathcal{Z}},\quad\text{and}\quad \hat{\rho}_p=\frac{e^{-\frac{1}{2}\beta\Omega\tilde{\sigma}^{z}_{p}}}{\mathcal{Z}}
\end{equation}
respectively. Here, $\mathcal{Z}=\text{Tr}[\exp(-\beta\omega_{j}\tilde{\sigma}^{z}_{j}/2)]$ is the partition function and $\beta=1/k_BT$ is the  inverse temperature. Henceforth, we take the Boltzmann constant $k_B=1$. We emphasize that the steady-state density matrix is written in the
eigenstates of the dressed Hamiltonian $\tilde{H}$, and these are given by the individual eigenstates of the ancilla and probe qubits as,
	\begin{eqnarray}\label{eq:eigenstate1}
	\ket{1} = \cos\frac{\theta_1}{2}\ket{++} - \sin\frac{\theta_1}{2}\ket{+-}, \nonumber \\
	\ket{2} = \sin\frac{\theta_1}{2}\ket{++} + \cos\frac{\theta_1}{2}\ket{+-}, \nonumber \\
	\ket{3} = \cos\frac{\theta_1}{2}\ket{-+} + \sin\frac{\theta_1}{2}\ket{--}, \nonumber \\
	\ket{4} = \cos\frac{\theta_1}{2}\ket{--} - \sin\frac{\theta_1}{2}\ket{-+}\label{eq:eigenstate4},
	\end{eqnarray}
	with their corresponding eigenvalues $\epsilon_1 = (\omega_{1} + \Omega)/2$,  $\epsilon_2 = (\omega_{1} - \Omega)/2$, $\epsilon_3 = (-\omega_{1} + \Omega)/2$, and $\epsilon_4 = (-\omega_{1} - \Omega)/2$, respectively.

	\begin{figure}[t]
		\centering
		\subfloat[]{
			\includegraphics[scale=0.98]{new.pdf}}\\
		\subfloat[]{
			\includegraphics[scale=0.92]{Fig2b.pdf}}
		\caption{ (a) Quantum Fisher information $\mathcal{F}_{Q}$ associated with the reduced state of the probe qubit $\hat{\rho}_{p}$ by performing local measurements on the two-level system probe as a function of the temperature $T$ for different coupling strengths between the ancilla and probe qubits. (b) The absolute values of coherence $\abs{c}$ as a function of thermal bath temperature $T$ for different coupling strengths between the ancilla and probe qubits. The solid, dot-dashed, dashed, and dotted curves are for $g_{1}=0.01, 0.02, 0.03, 0.04$, respectively. In both panels, the other system parameters are the same: $\omega_p=1$ and $\omega_1 =0.04$.}\label{fig2}
	\end{figure}

Our probe is a single qubit placed outside the thermal bath (Fig.~\ref{fig1}), and we wish to perform local measurements on this probe qubit. Therefore, by using the transformation given in Eq.~(\ref{eq:transformation}), the reduced state of the probe qubit in the local basis is given by
	\begin{equation}
	\hat{\rho}_{p}= \text{Tr}_1[U\tilde{\rho}_\text{ss}U^\dagger]=\frac{1}{2}\begin{pmatrix}
	1-\chi & c\\
	 c & 1 +\chi
	\end{pmatrix},
\end{equation}\label{eq:DensityProb}
where the diagonal (populations) and off-diagonal (coherences) terms are determined by
\begin{eqnarray}
\chi &=&\cos(\theta_1)\tanh(\frac{\Omega }{2 T}), \quad\text{and} \\
c &=& \sin(\theta_1)\tanh(\frac{\Omega}{2T})\tanh (\frac{\omega_1}{2 T}),
\label{eq:coherences}
\end{eqnarray}
respectively. In addition to free Hamiltonian and coherent interaction parameters $\omega_p,\omega_1$ and $g_1$, through $\theta_1$ and $\Omega$, both the coherence and populations of the probe qubit
carry information on the temperature of the sample. In contrast to thermalized probes with the sample,
such a non-thermal probe state allows us to manipulate QFI with additional terms associated with the temperature dependent coherences. We shall explore different system parameter regimes to search for
multiple peaks and enhancement thermal sensitivity effects in the low temperature regimes in QFI,
computed by using Eq.~(\ref{eq:TLSfisher}). It is worth noting that the QFI is computed by considering the local measurements on our probe qubit.

The complete expression for the QFI is too cumbersome to report here. We depict QFI $\mathcal{F}_{Q}$, and absolute coherence $\abs{c}$ in the probe as a function of the bath temperature $T$ for different values of coupling strength $g$ in Figs.~\ref{fig2}(a) and~\ref{fig2}(b), respectively. It is immediately evident that the lower peak in the QFI is associated with the coherence developed in the probe qubit. By increasing the ancilla-probe coupling strength, coherence in the probe qubit become larger, thus enhancing the probe's thermal sensitivity at lower temperatures. It is often useful to evaluate the relative error bound in estimating the temperature for the performance of a given probe~\cite{Mukherjee2019,glatthard2022bending}. We discuss the relative error bound in the case of our scheme in Appendix~\ref{app:relative}.

We emphasize that, for a two-level system probe at thermal equilibrium with the bath, it is not possible to have multiple peaks in the QFI. To get better insight into the underlying working mechanism of our probe and find the locations of the peaks in the QFI, we first briefly revisit a two-level system probe at thermal equilibrium. For a two-level system probe of frequency $\omega_{0}$ at thermal equilibrium with the bath of temperature $T$, the density matrix of the probe is given by
	\begin{equation}
	\hat{\rho}_\text{TLS}=\frac{1}{2}\begin{pmatrix}
	1-\text{tanh}\big(\frac{\omega_0}{2T}\big) & 0\\
	 0 & 1+\text{tanh}\big(\frac{\omega_0}{2T}\big)
	\end{pmatrix}.
	\end{equation}\label{eq:DensityTLS}
QFI associated with this two-level system probe can be evaluated using Eq.~(\ref{eq:TLSfisher}), which yields
	\begin{equation}
	\mathcal{F}_\text{TLS} = \frac{\omega^2_0\text{sech}^2\big(\frac{\omega_0}{2T}\big)}{4T^4}.
	\end{equation}
This expression indicates that the QFI associated with a two-level system probe at thermal equilibrium has a single peak. Accordingly, a two-level system with a given energy gap $\omega_{0}$, can only be an optimal probe for a single temperature \cite{campbell2018precision}. At this temperature, $T^{*}$, QFI has a maximum value, and the temperature $T^{*}$ can be found by \cite{campbell2018precision}
	\begin{equation}
	T^{*} = \gamma \omega_{0}, \quad \text{here} \quad 2\gamma = \text{tanh}\bigg(\frac{1}{\gamma}\bigg).
	\end{equation}

The diagonal elements of our probe qubit density matrix appear similar to the two-level system probe at thermal equilibrium, except for an additional factor of $\cos\theta_1$. This should not be surprising, as the nonlocal dissipation channels in the master equation (\ref{eq:master}) provide indirect access of the thermal bath to our probe. Accordingly, the information about the temperature $T$ of the bath is partially imprinted on the diagonal elements of our probe qubit. This is witnessed by the emergence of a peak in the QFI associated with higher bath temperature $T^{*}_p$ (see the right peak in Fig.~\ref{fig2}(a)). The coherent interaction between the probe and ancilla provides an additional channel of indirect energy exchange between the bath and the probe. This information is imprinted on the off-diagonal elements of our probe's density matrix, which is responsible for the emergence of an additional peak at lower temperature $T^*_{1}$ in the QFI.

In the weak coupling limit of ancilla and probe qubits $g\ll (\omega_p,\omega_1)$, we can replace $\text{sin}\theta_1\approx\theta_1$, and $\text{cos}\theta_1\approx1$ in the master equation (\ref{eq:master}). In this case, we can write a simple approximate expression for the QFI,
	\begin{eqnarray}\label{eq:QFIapprox}
	\mathcal{F}_{Q} &\approx & \mathcal{F}_{0} + \mathcal{F}_{1}, \quad\text{where}\\
	\mathcal{F}_{0} &=& \frac{\omega^2_p\text{sech}^2\big(\frac{\omega_p}{2T}\big)}{4T^4}, \quad \mathcal{F}_{1} = \frac{\frac{1}{2}\theta^2_{1}\omega^2_1\text{sech}^2\big(\frac{\omega_1}{2T}\big)}{4T^4}.
	\end{eqnarray}

	\begin{figure}[t!]
	\centering
		\includegraphics[scale=0.9]{Fig3.pdf}
	\caption{Quantum Fisher information $\mathcal{F}_{Q}$ associated with the qubit probe---considering weak coupling between the ancilla and probe---as a function of temperature $T$ of the thermal bath. The blue-dashed and red-solid lines correspond to approximate and exact QFI, respectively. The approximate expression for the QFI in the weak coupling limit is given by Eq.~(\ref{eq:QFIapprox}), and the exact QFI is evaluated by considering the full expression of QFI, which is too cumbersome to report here. Parameters: $\omega_{p}=1$, $\omega_1=0.04$, and $g=0.01$.}\label{fig3}
	\end{figure}

This approximate expression of the QFI can be used to find the location of the peaks in the QFI. The temperatures associated with the left and right peaks in Fig.~\ref{fig2}(a) are given by $T^*_{1} = \gamma \omega_{1}$, and $T^*_{p} = \gamma \omega_{p}$, respectively.

		In Fig.~\ref{fig3}, we depict approximate and exact QFI given in Eq.~(\ref{eq:QFIapprox}), and Fig.~\ref{fig2}(a), respectively. It is evident that our approximate results for QFI overlap with the exact QFI within the weak coupling limit. The location of the lower peak can be controlled by ancilla qubit energy gap $\omega_1$, and baths at lower temperatures can be probed with increased accuracy by making the energy gap smaller. In addition, the coupling strength between ancilla and probe qubits can be used to enhance the amplitude of the lower peak, thus enhancing the thermal sensitivity of the probe. The right peak in Fig.~\ref{fig3} is identical to the QFI for a two-level system probe at thermal equilibrium.

	\begin{figure}[b]
	\centering
		\includegraphics[scale=0.9]{Fig4.pdf}
	\caption{ Quantum Fisher information $\mathcal{F}_{Q}$  with respect to the bath temperature $T$ in the case of two ancilla qubits. The solid, dot-dashed, and dashed lines are for $\omega_p=0.26$, $\omega_\text{p}=0.3$, and $\omega_\text{p}=0.4$, respectively. We have added the QFI (dotted line) of a two-level system at thermal equilibrium for a comparison.
Rest of the parameters are: $\omega_{1}=0.09$, $\omega_{2}=0.17$, $g_1=0.003$, and $g_2=0.05$.} \label{fig4}
	\end{figure}

It is important to note that for a single ancilla qubit, the QFI can be approximately given by the sum of two $\text{sech}^2\theta$ functions with different amplitudes and centers (Eq.~(\ref{eq:QFIapprox})). As a result, two peaks emerge in the QFI, each of which can be identified by the frequencies of the ancilla and probe qubits. This hints that including additional ancilla qubits may increase the number of peaks in the QFI. Similar results are reported in a periodically driven probe~\cite{Mukherjee2019}; however, the physical mechanism of multiple peaks in the QFI is entirely different in these non-autonomous schemes. In the next section, we will investigate the possibility of increasing the optimal sensing points of the probe by increasing the number of ancilla qubits.

\subsection{Multiple peaks in QFI with 2 ancilla qubits ($N=2$)}\label{sec IIIB}

	To reveal the full advantange of our scheme for the enhancement of low temperature thermal sensitivity, here we consider two ancilla qubits. In this case, the diagnolized Hamiltonian of ancilla-probe system is given by
	\begin{equation}
	\tilde{H} = \sum^2_{k=1}\omega_{k}\tilde{\sigma}^{z}_{k} + \tilde{\Omega}\tilde{\sigma}^{z}_{p},
	\end{equation}
here the transformed Pauli matrices are given in Eq.~(\ref{eq:Paulitrans}). The transformed frequency of the probe qubit is given by
	\begin{eqnarray}\label{eq:omega2}
	\tilde{\Omega} = \Omega_{+} + \Omega_{-}\quad \text{where}\quad  \Omega_{\pm} &=& \sqrt{\omega^2_p\pm4(g_1+g_2)^2}.\nonumber\\
	\end{eqnarray}
The master equation for this case remains the same as given in Eq.~(\ref{eq:master}), only the mixing angle $\theta_{k}$ is changed to
\begin{eqnarray}\label{eq:angle2}
\theta_1&=&\frac{1}{2}[\arctan(\frac{2(g_1+g_2)}{\omega_p})+\arctan(\frac{2(g_1-g_2)}{\omega_p})],\nonumber\\
\theta_2&=&\frac{1}{2}[\arctan(\frac{2(g_1+g_2)}{\omega_p})-\arctan(\frac{2(g_1-g_2)}{\omega_p})].
\end{eqnarray}

At steady-state, the ancilla-probe system reaches an equilibrium state under the dissipative dynamics governed by the master equation~(\ref{eq:master}). We are only interested in the reduced state of the probe, as we wish to perform local measurements on the probe qubit. The state of the probe in the local basis is given in Appendix~\ref{App:2Ancilla}, which shows the presence of temperature-dependent coherence induced by both ancilla qubits. By the judicious choice of system parameters, QFI exhibits three peaks, as shown in Fig.~\ref{fig4}. The emergence of the third peak is associated with additional ancilla qubit in this case. The addition of further ancilla qubits will enhance the number of peaks in the QFI, hence enhancing the accuracy of temperature estimation and broadening the range of thermal sensitivity.

\section{Conclusion}\label{sec IV}
We have investigated the effects of the probe's genuine quantum features (quantum coherence) on its low temperature thermal sensitivity. We have considered a two-level system (qubit) as a probe coupled indirectly with the sample via ancilla qubits. Typically, a probe is assumed to be in thermal equilibrium with the sample. In this case, temperature information is encoded only on the diagonal elements of the probe's density matrix. Accordingly, the QFI becomes equal to the classical information, and other than the discreteness of energy level, the quantum probe does not exhibit any quantumness. A two-level system probe at thermal equilibrium can only estimate a single temperature with optimal accuracy, which can be seen by the existence of a single peak in the QFI. On the contrary, we have shown that in the presence of ancilla qubits, two different channels imprint temperature information on the probe. (i) The nonlocal dissipation channel revealed by deriving a global master equation imprints information on the diagonal elements. (ii) Importantly, ancilla qubits assist the exchange of energy between the sample and the probe and encode the temperature information on the off-diagonal elements of the probe's density matrix. We have then shown that the coherence in the probe's qubit, at steady-state, can lead to the emergence of multiple peaks in the QFI at low temperatures depending on the number of the ancilla qubits. Thus, enhancing the estimation accuracy and broadening the thermal sensitivity range of the quantum probe. Our analysis revealed that quantum features such as coherence could assist a quantum thermometer in going beyond the single low-temperature temperature optimal estimation.
Contrary to the other studies, our scheme does not require quantum probes with highly degenerate excited states or energetically costly non-autonomous schemes based upon dynamically controlled probes to generate multiple peaks of QFI.
\section*{Acknowledgment }
\"{O}.~E.~M.~acknowledges support from TUBITAK Grant No.~122F371.


\section{Introduction}

Ocean swell in the wake of a storm or driven by the force of strong wind is seen as a collection of almost periodic unidirectional waves from a beach or an ocean liner.
The waves that do not vary in the direction transverse to propagation are known as Stokes waves.
Such waves were discovered by Stokes~\cite{stokes1847theory} who  conjectured~\cite{Stokes1880_314} that the wave of greatest height, referred to as the extreme Stokes wave, has a crest angle of $120$\degree.
That conjecture was proven
more than 100 years later, see \cite{AmickFraenkelTolandActaMath1982, PlotnikovDinamikaSploshnSredy1982}.
The study of Stokes waves has been the subject of numerous works, see e.g.~\cite{michell1893,Nekrasov1921,MalcolmGrantJFM1973LimitingStokes,SchwartzJFM1974,Longuet-HigginsFoxJFM1977,Longuet-HigginsFoxJFM1978,toland1978existence,PlotnikovDinamikaSploshnSredy1982,AmickFraenkelTolandActaMath1982,Williams1981,WilliamsBook1985,TanveerProcRoySoc1991,gandzha2002steep,Longuet-HigginsWaveMotion2008,DyachenkoLushnikovKorotkevichJETPLett2014,DyachenkoLushnikovKorotkevichPartIStudApplMath2016,LushnikovStokesParIIJFM2016,LushnikovDyachenkoSilantyevProcRoySocA2017}.
Stokes waves are a one-parameter family of solutions often characterized by the steepness $s=H/L,$ where $H$ is the crest-to-trough height and $L$ is the wave length.
The steepness of the extreme Stokes wave is reached at $s_{lim}=0.14106348398\ldots$~\cite{LushnikovDyachenkoSilantyevProcRoySocA2017} (see also~\cite{dhs2022almost} for a more accurate value).
However, the curvature radius at the crest of the Stokes wave vanishes as the extreme wave is approached. In the ocean, once the capillary length scale $\simeq 1.7$cm (for ocean water) is reached, the surface tension effects cannot be disregarded.
This limits the physical relevance of Stokes waves in the immediate vicinity of the extreme wave.
Nevertheless, the wavelength of gravity waves in ocean swell often ranges from a dozen centimeters to a few kilometers, see e.g.~\cite{YoungBookOceanWaves1999,liu2003effects}.
Thus, without oversimplification, a wide range of Stokes waves is relevant for modeling oceanic swell.
The Stokes waves found in~\cite{DyachenkoLushnikovKorotkevichPartIStudApplMath2016,LushnikovDyachenkoSilantyevProcRoySocA2017} were used to determine the radius of curvature at the crest.
Surface tension effects are small for $s \lesssim s_{50} = 0.14100$ and $s\lesssim s_{100} = 0.141034$ for Stokes waves
with wave length $50$ and $100$ meters respectively.
In what follows we provide stability results for $s\le  0.1408627$ that are well below these applicability estimates.
The nonlinear stage of instability, as observed in experiments or computation, reveals that wave breaking leads to the rapid vanishing of the curvature radius to the point that surface tension becomes non-negligible, and whitecaps form through ejection of capillary waves from spilling/plunging breakers, see~\cite{longuet1963generation,duncan2001spilling,dyachenko2016whitecapping}.

\begin{figure}
\includegraphics[width=0.495\textwidth]{nlin.eps}
\vspace*{-0.25in}
\caption{Initial stage of instability development at the crest of a Stokes wave (solid black) from the nonlinear dynamics (red line) compared to a perturbation of the Stokes wave by the eigenfunction corresponding to a localized instability (black dashed line). The later stage of the dynamics reveals the formation of a breaker (orange \& green) and overturning (inset).
}
\vspace*{-0.25in}
\label{fig:breaking}
\end{figure}

Since \cite{KharifRamamonjiarisoaPhysFluids1988} indicates that transverse perturbations of steep waves have smaller growth rates than do longitudinal ones,
we examine the growth rates of longitudinal instabilities of Stokes waves in water of infinite depth, ignoring the effects of surface tension, vorticity, wind and dissipation.
We consider {\em all} quasiperiodic perturbations, {\em i.e.}, we include {\it superharmonic}  (with period $L$) and {\it subharmonic} (with larger period, an integer multiple of $L$) perturbations, using the terminology of~\cite{Longuet-HigginsFoxProcRoySoc1978superharmonics, Longuet-HigginsFoxProcRoySoc1978subharmonics}.
The subharmonic instability of small-amplitude Stokes waves (in water waves and other fields) has been studied since 1965, see Refs.~\cite{benjamin1967disintegration, bridgesmielke, Lighthill1965,whitham1967non,Zakharov1968,ZakharovOstrovskyPhysD2009} and is now well understood \cite{berti2021full, creedon2022ahigh, nguyenstrauss}.
 \begin{figure*}
    \centering
    \includegraphics[width=0.364\textwidth]{hamiltonian.eps}
    \includegraphics[width=0.621\textwidth]{hamiltonian2.eps}
    \vspace*{-0.1in}
    \caption{(Left panel) Hamiltonian $\cal H$ given by~\eqref{HamConf} as a function of steepness $s$. Green denotes the region of dominance of the BF instability, red corresponds to dominance of the localized branch.
    %
    (Right panel) $\cal H$ as a function of $-\ln\left( s_{lim} - s\right)$ oscillates near the extreme wave. The points marked $s_{max,1},s_{min,1},s_{max,2}$ are the first three extrema of $\cal H$. For details about $s_1, s_2, s_3$ and $s_4$ we refer the reader to the main body of the text.
    }
    \label{fig:ham}
    \vspace*{-0.2in}
\end{figure*}
This instability is referred to as the Benjamin-Feir (BF) or modulational instability.
The growth rate of the BF instability increases with $s$, passes through its maximum and decays to zero.
The so-called {\em localized} instability branch
appears for steep waves before the BF instability ceases.
The eigenfunctions associated with this branch change rapidly in the vicinity of the wave-crest.
Increasing $s$ gives rise to additional localized instability branches.
These branches correspond to both the sub-
\cite{Longuet-HigginsFoxProcRoySoc1978subharmonics,McLeanMaMartinSaffmanYuenPRL1981,KharifRamamonjiarisoaPhysFluids1988}
and the superharmonic disturbances (the latter form at the extrema of the Hamiltonian, {see Fig.~\ref{fig:ham}})
\cite{Longuet-HigginsFoxProcRoySoc1978superharmonics,tanaka1983stability,LonguetHigginsTanakaJFM1997,KorotkevichLushnikovSemenovaDyachenkoStudApplMath2022accepted}.
See also \cite{KharifRamamonjiarisoaPhysFluids1988,KorotkevichLushnikovSemenovaDyachenkoStudApplMath2022accepted} for earlier numerical studies of different branches.

It has been unclear as to which of the instabilities dominates the dynamics of almost-extreme Stokes waves.
The present work elucidates this long-standing question and establishes the localized branch as the dominant one for steep waves.
Depending on the value of $s$, either a perturbation co-periodic with the Stokes wave, or the one having twice its period is the most unstable.
The dominance of the latter perturbations interchanges repeatedly as the extreme wave is approached.
Growth rates observed for near-extreme waves achieve large values, resulting in instabilities that would swiftly disintegrate wave trains of large-amplitude Stokes waves.
This is crucial for understanding why long-lived oceanic swell is never observed in the ocean.

\begin{figure*}
    \centering
    \includegraphics[width=0.395\textwidth]{growth2.eps}
    \includegraphics[width=0.597\textwidth]{multiA.eps}

    \vspace*{-0.25in}
    \caption{(Left panel) The maximal instability growth rate $\gamma$ as a function of steepness $s$. BF (green) is the dominant instability for small steepness, and the localized branch (red) dominates for large steepness. The blue curve is the theoretical curve from~\cite{creedon2022ahigh}.
    The BF figure eight detaches at steepness $s_1$ and begins to shrink. At $s=s_2$, the dominant instability transitions to the localized branch with real eigenvalue. The inset shows how the
    transition to the localized branch occurs. The spectral plane
    for waves $(a)$-$(d)$ is shown on the right.
    %
    (Right panel) Spectrum of the problem (real part vs. imaginary part) in the transition region $(a)$-$(d)$ from the inset of the left panel. Red circles correspond to the subharmonic instability ($\mu = 0.5$) of the localized branch and green circles are the eigenvalues with largest real part of the BF remnant.
    }
    \vspace*{-0.1in}
    \label{fig:growth}
\end{figure*}

\noindent {\bf Problem Formulation}. We consider potential motion of an ideal fluid in a 2D domain, ${\bf r} = (x,y)^T\in\mathcal{D}(t) \subset \mathbb{R}^2$.
The fluid velocity ${\bf v}$ is the gradient of the potential $\varphi({\bf r},t)$: ${\bf v}= \nabla \varphi$ with $\nabla\equiv(\partial_x,\partial_y)$ which satisfies $\nabla^2 \varphi = 0$.
The fluid domain is bounded by a $2\pi$-periodic free surface and extends downward to infinity.
The Hamiltonian of the Euler equations with a free surface is recast in terms of surface variables~\cite{Zakharov1968}, and the conformal mapping approach from \cite{TanveerProcRoySoc1991,DyachenkoZakharovKuznetsovPlasmPhysRep1996} is used.
Here $w=u+iv\in \mathbb{C}^-$ parameterizes the fluid domain through the conformal mapping $z(w,t) = x(w,t) + i y(w,t) \in \mathcal{D}(t)$.

The Hamiltonian of the fluid flow is the sum of the kinetic and  potential energy due to gravity.
The Eulerian dynamics with a free surface is obtained from the least action principle~\cite{DyachenkoEtAl1996}.
The Hamiltonian in conformal
variables is
\begin{align}
    \mathcal{H} = -\frac{1}{2}\int\limits_{-\pi}^\pi \psi \hat H \psi_u \,du + \frac{g}{2} \int\limits_{-\pi}^\pi y^2 x_u \,du, \label{HamConf}
\end{align}
where $g$ is the acceleration of gravity, and $\hat H$ is the circular Hilbert transform. The Lagrangian has the form
\begin{align}
\mathcal{L} = \int\limits_{-\pi}^\pi \psi \left(y_t x_u - y_u x_t\right)\,du - \mathcal {H}.
\end{align}
Varying the action $\int \mathcal{L}\,dt$, with respect to $x$, $y$ and $\psi$, holding the Cauchy-Riemann relations for $x$ and $y$ as a constraint, one arrives at the equations of motion~\cite{DyachenkoEtAl1996,DyachenkoLushnikovZakharovJFM2019}:
\begin{align}
\label{impbase}
\begin{array}{l}
y_t x_u - y_u x_t = -\hat H\psi_u,  \\
x_t \psi_u \!-\! x_u \psi_t \!-\!\hat H\left[y_t \psi_u -  y_u\psi_t \right] =
g(x_u y \!-\! \hat H \left[yy_u\right]),
\end{array}
\end{align}
where $x = u - \hat H y$.
 It has been shown~\cite{DLZ1995Fivewave} that $y$
and $\mathcal{P} = x_u\psi - \hat H\left[y_u\psi\right]$ are canonical conformal variables.

The Stokes wave is a solution of \eqref{impbase}, traveling at constant speed $c$, measured in units of $c_0=\sqrt{g/k}$, the speed of linear gravity waves.
The time variable is normalized to the frequency of linear gravity waves and is measured in units $(gk)^{-1/2}$.
After rescaling the spatial and conformal coordinates we set the
period of the Stokes wave to be $2\pi$ in $u$ (and $x$) without loss of generality.
Writing
\eqref{impbase} in the moving frame using $y(u-ct)$ and $\psi(u-ct)$, we linearize
around the Stokes wave in canonical conformal variables $y$ and $\mathcal{P}$ as in~\cite{DyachenkoSemenova2022}.
This results in a linear  integro-differential equation with variable periodic coefficients.
Quasi-periodic solutions are sought in the form~\cite{deconinck2006computing}%
\begin{equation}\label{deltapsideltay}
\begin{pmatrix} \delta \mathcal{P}(u,t) \\
\delta y(u,t)
\end{pmatrix}=e^{i\mu u+\lambda t}\begin{pmatrix}\delta \mathcal{P}_0(u) \\
\delta y_0(u)
\end{pmatrix}
\end{equation}
where $\delta y_0(u)$ and $\delta \mathcal{P}_0(u)$ are $2\pi$-periodic functions, $\mu\in(-1/2,1/2]$ is the Floquet exponent, and $\lambda$ is a complex eigenvalue whose real part is the growth rate $\gamma = \mbox{Re}\, \lambda$ (if positive).
The stability spectrum of a Stokes wave is found as the union of point spectra for all values of $\mu$~\cite{deconinck2006computing, kapitulapromislow}.

\begin{figure*}
    \centering
    \includegraphics[width=0.495\textwidth]{growth2log.eps}
    \includegraphics[width=0.495\textwidth]{efgh_multiplot.eps}
    \vspace*{-0.25in}
    \caption{(Left panel) The maximal growth rate $\gamma$ as a function of the variable $-\ln(s_{lim}-s)$.
    %
    %
    The red curve represents the localized instability with real eigenvalues for $\mu = 0$ and $\mu = 1/2$.
    %
    %
    %
    At $s=s_3$, the dominant instability becomes superharmonic ($\mu=0$). At $s_4$ it switches back to subharmonic.
    (Right panel) Eigenvalues for steepness (e)-(h) in the vicinity of $s_{max,1}$: (from left to right) $s = 0.1341061714$, $s = 0.1364173038$, $s = 0.1366051511$ and $s = 0.1378879221$. As the steepness increases, superharmonic eigenvalues (green) move along the imaginary axis toward the origin and collide. As a result a pair of symmetric instability isole (h) appear and move away from the origin.
    }
    \vspace*{-0.2in}
    \label{fig:all}
\end{figure*}

The fast Fourier transform (FFT) is used to apply the linearization operator to the periodic functions appearing in the eigenvalue problem.
The resulting spectral problem in canonical conformal variables $y$ and $\mathcal{P}$~\cite{DyachenkoSemenova2022} is amenable to the shift-invert method, applied
in $O(m N \ln N)$ flops, where $N$ is the number of Fourier
modes and $m$ is the number of iterations to solve the linear system embedded in the shift-invert operator,
using the MINRES method~\cite{saad1992numerical}. Combined with the Krylov-Schur algorithm~\cite{lehoucq1989arpack}, the numerical solution of the eigenvalue problem is solved with up to $N = 10^{6}$ Fourier modes.

\noindent {\bf Main Results.}
The stability of Stokes waves of small to moderate steepness was studied
using variations of the Fourier-Floquet-Hill method \cite{deconinck2006computing} by many authors, {\em e.g.}, \cite{DeconinckOliverasJFM2011, KharifRamamonjiarisoaPhysFluids1988,mclean1982instabilities}.
Our results confirm that for $s < 0.1381465$, all Stokes waves are found to be unstable, see Fig.~\ref{fig:growth}.
For higher steepness waves, our results are new for subharmonic perturbations ($\mu \neq 0$). They are in agreement with \cite{KorotkevichLushnikovSemenovaDyachenkoStudApplMath2022accepted} for superharmonic perturbations ($\mu = 0$).

For steepness $0<s< s_2=0.12894$, we find that the dominant instability of Stokes waves is the modulational (BF) instability, corresponding to subharmonic perturbations, see Fig.~\ref{fig:growth}. The figure-$8$
detaches from the origin at $s_1 = 0.10423$, shortly after the maximum growth rate, $\gamma_{BF} = 0.02357$ is attained. At $s_2 = 0.12894$, the dominant instability switches (Fig.~\ref{fig:growth}) to the {\it localized} branch, described next.

The localized instability branch first appears at the steepness $s_{loc,1} = 0.12890308$ from a
collision of eigenvalues with $\mu = 1/2$ at the origin, resulting in an oval of eigenvalues, see Fig.~\ref{fig:all}e. As $s$ increases, the oval is pinched vertically (Fig.~\ref{fig:all}f) as more eigenvalues from the imaginary axis flood onto it. This results in formation of a figure ``infinity'' at $s=s_{max,1}$ (Fig.~\ref{fig:all}g), for the Stokes wave corresponding to the first maximum of the Hamiltonian (Fig.~\ref{fig:ham}). With further increase of $s$, the figure infinity is pinched at the origin, its two isole drifting away along the real axis, as their diameter decreases (Fig.~\ref{fig:all}h). We refer to the localized instability branch as such because its eigenvalues for $\mu\in(-\frac{1}{2},\frac{1}{2}]$ have eigenfunctions that are localized near the wave crest as the extreme wave is approached, see Fig.~\ref{fig:eigf0}.
After detachment ($s > s_{max,1}$) the entire range of the Floquet parameter $\mu \in(-\frac{1}{2}, \frac{1}{2}]$ covers the localized instability branch, in contrast to the high-frequency instabilities discussed in \cite{DeconinckOliverasJFM2011, cdo}.

On the unstable localized branch isola, the subharmonic mode with $\mu=1/2$ has the dominant growth rate for $s_{max,1} < s < s_3 = 0.139492$. Through a sequence of topological changes to the isola (Fig.~\ref{fig:bubble1}i-l), the superharmonic mode ($\mu = 0$) becomes dominant for $s_3 < s < s_4 = 0.140613$. For Stokes waves with larger steepness, these switches repeat. One such additional switch is illustrated in Fig.~\ref{fig:bubble1}m-p. We conjecture that the dominance of unstable mode switches between $\mu=0$ and $\mu=1/2$ an infinite number of times as the extreme wave is approached.

As $s$ increases, a secondary modulational instability branch appears and vanishes, repeating the stages of the BF instability.
Similarly, a secondary localized branch appears through the collision of a pair of $\mu=1/2$ eigenvalues at the origin, repeating the stages seen in Fig.~\ref{fig:bubble1}i-l. The figure infinity detaches from the origin at the next extremum of the Hamiltonian at $s = s_{min,1}$. We conjecture that
new modulational and new localized instability branches appear infinitely many times
as the extreme wave is approached.

\begin{figure}
\centering
\includegraphics[width=0.49\textwidth]{ijkl_multiplot.eps}
\includegraphics[width=0.49\textwidth]{mnop_multiplot.eps}
\vspace*{-0.25in}
\caption{Spectrum in transition regions $(i)$-$(l)$ and $(m)$-$(p)$ from the
left panel of Fig.~\ref{fig:all}: around $s_3$ (top panel, from left to right) $s = 0.13946478$,
$0.13948945$, $0.13951484$ and $0.13953804$; and around $s_4$ (bottom panel,
from left to right) $s = 0.1405850778$, $0.1406007221$,
$0.1406080087$ and $0.1406384552$. Red circles correspond to $\mu = 0.5$ (subharmonic) and green circles to $\mu=0$ (superharmonic). This demonstrates the switch between the two types of instabilities. At $s_5 = 0.14077 < s_{min,1}$, the interchange of dominance between $\mu =0$ and $\mu = 1/2$ occurs again.}
\vspace*{-0.1in}
\label{fig:bubble1}
\end{figure}

\begin{figure}
\includegraphics[width=0.485\textwidth]{eigf2mu012.eps}
\vspace*{-0.2in}
\caption{
Eigenfunctions of the localized instability branch for $\mu=0$ (red) and $\mu = 1/2$ (green) of the Stokes wave with steepness $s = 0.1406007221$ (blue). The corresponding wave is sketched in a dashed line.  Insets show zooms near peaks.
}
\label{fig:eigf0}
\vspace*{-0.2in}
\end{figure}

Eigenfunctions on the localized instability branch are shown in Fig.~\ref{fig:eigf0} for $\mu=0$ (two periods) and $\mu=1/2$ (one period).
As $s$ increases, the eigenfunctions are increasingly localized near the crest of the wave.
When used as a perturbation to seed the nonlinear dynamics about a Stokes wave, the localized instability leads to the steepening and overturning of the wave, see Fig.~\ref{fig:breaking}.
Simultaneously, the dominant unstable isola drifts to infinity, as the extreme wave is approached. An infinite growth rate is not physical, indicating the breakdown of our model. Indeed, as increasingly steep
waves approach the extreme 120\degree\, Stokes wave, the effects of
surface tension, among others, play a role and have to be incorporated~\cite{banner2014linking}.

\noindent {\bf Conclusion}\,
While the instability spectrum is quite complicated, the instability consequences for highly nonlinear Stokes waves are straightforward. Unlike for small-amplitude waves \cite{DeconinckOliverasJFM2011}, the dominant eigenvalue is real and is associated with
eigenfunctions ({\em i.e.}, perturbations) that are co-periodic
with the Stokes wave or have double its period. In either case, the instabilities that dominate the dynamics are localized near the wave's crest.

Our numerical observations are particularly significant because they
answer long-standing questions of which instability dominates the
evolution of steep ocean waves. Further, our computations illustrate
the complexity of the spectral stability problem governing the different
instabilities of large-steepness Stokes waves. We have focused on the
dominant instabilities of steep waves, but the stability spectrum of such
waves has many sub-dominant components, which will be discussed elsewhere.

\noindent {\bf Acknowledgements}
The authors thank T. J. Bridges and A. O. Korotkevich for helpful discussions.
The work of P.M.L. was supported by NSF-DMS-1814619.
P.M.L. and S.D. thank the Isaac Newton Institute for Mathematical
Sciences, Cambridge, UK for support and hospitality during the programme ``Dispersive hydrodynamics''
where work on this paper was partially undertaken.
S.D. thanks the Department of Applied Mathematics at the U. of Washington for hospitality. A.S and S.D. acknowledge the FFTW project and its authors~\cite{frigo2005design} as well as the entire GNU Project.
A.S. thanks the Institute for Computational and Experimental Research in Mathematics, Providence, RI, being resident during the ``Hamiltonian Methods in Dispersive and Wave Evolution Equations" program supported by NSF-DMS-1929284.