\section{I. INTRODUCTION}
Recent advances in fabricating epitaxial thin films have opened a door to utilizing and manipulating the lattice strain by controlling the lattice mismatch with substrate. This strain engineering has a great capability for realizing various physical properties. For example, the epitaxial strain enhances the superconducting transition temperature ($T_c$) in cuprate La$_{2-x}$Sr$_x$CuO$_4$, controls the magnetism and electrical conductivity in manganite La$_{1-x}$Sr$_x$MnO$_3$, induces the room-temperature ferroelectricity in titanate SrTiO$_3$, and realizes a quantum Hall phase in bulk mercury telluride HgTe \cite{Locquet98, Konishi99, Haeni04, Brune11}. Recently a sizable strain effect was reported in iron-based superconductor FeSe \cite {Nabeshima172602, Nie09, Tan13}.

Strain-free bulk FeSe ($T_c$ = 8 K) is a superconductor \cite{Hsu08} which lies in the Bardeen-Cooper-Schrieffer (BCS) and Bose-Einstein condensation (BEC) crossover regime \cite{Kasahara16309}, and exhibits nematic order \cite{Nakayama237001, Shimojima121111, Watson155106, Zhang214503, Huynh144516, Baek15, Bohmer027001} where the C$_4$ rotational symmetry in the electronic structure is broken. It was reported that compressive strain enhances the $T_c$ by 50$\%$ \cite{Nabeshima172602}, while tensile strain kills the superconductivity \cite{Nie09} and simultaneously enhances the nematicity \cite{Tan13}. A question naturally arises as to what triggers such characteristic strain effects. Elucidating this issue would provide a clue to understanding the origins of superconductivity and nematicity as well as their interplay, which are at the heart of intensive investigations of iron-based superconductors. Since the low-energy electronic structure, which governs the key electronic properties, is sensitive to lattice strain \cite{Dhaka14, Gofryk14}, the experimental determination of the electronic structure in strained FeSe is of crucial importance. However, such an experiment has been limited to the nematic phase of the tensile-strained nonsuperconducting regime \cite{Tan13}.

In this article, we traced the evolution of electronic structure with epitaxial strain by angle-resolved photoemission spectroscopy (ARPES) on FeSe. We fabricated a compressive-strained FeSe film on CaF$_2$ substrate (FeSe/CaF$_2$; $T_c$ = 12 K) and a tensile-strained film on SrTiO$_3$ (FeSe/SrTiO$_3$; nonsuperconducting), and compared their electronic structure with that of strain-free bulk FeSe ($T_c$ = 8 K).

 \begin{figure*}
 \includegraphics[width=6.4in]{fig1.eps}
\caption{(color online) Comparison of electronic structure among bulk FeSe, FeSe/CaF$_2$, and FeSe/SrTiO$_3$. (a) ARPES intensity map at $E_F$ as a function of two-dimensional wave vector measured in the nematic phase ({\it T }= 30 K) for strain-free bulk FeSe ($T_c$ $\sim$ 8 K). The intensity was obtained by integrating the spectral intensity within $\pm $10 meV with respect to $E_F$. (b) Near-$E_F$ ARPES intensity plot for bulk FeSe as a function of binding energy and wave vector measured in the normal state ({\it T} = 120 K) along cuts 1 and 2 shown in (a). The data were divided by the Fermi-Dirac function convoluted with the instrumental resolution. Dashed curves are a guide for the eyes to trace the band dispersions. (c) Near-$E_F$ ARPES intensity measured at {\it T} = 30 K along cuts 2 and 3 in (a). (d)-(f) and (g)-(i) Same as (a)-(c), but for compressive-strained FeSe/CaF$_2$ ($T_c$ $\sim$ 12 K) and tensile-strained FeSe/SrTiO$_3$ ($T_c$ $\sim$ 0 K), respectively. The ARPES intensities in (e) and (h) were measured at {\it T} = 120 K and 180 K, respectively. (j) Second-derivative intensity plot measured along cut 4 in (g) at {\it T} = 100 K, obtained after dividing by the Fermi-Dirac function convoluted with the instrumental resolution. Blue dashed lines are a guide for the eyes to trace the Dirac-cone-like dispersion. (k) Schematic illustration of Dirac-cone dispersion around the {\it M} point. Blue, light blue, pink, and red arrows highlight the Lifshitz transition of FS. The former two highlight the case of FeSe/SrTiO$_3$, while the latter two the case of bulk FeSe and FeSe/CaF$_2$.}
\end{figure*}

\section{II. METHODS}
High-quality FeSe thin films were grown on CaF$_2$ and SrTiO$_3$ substrates by pulsed-laser-deposition and molecular-beam-epitaxy methods, respectively. The film thicknesses of FeSe/CaF$_2$ ($T_c$ = 12 K) and FeSe/SrTiO$_3$ (non-superconducting) samples are approximately 300 monolayers (MLs) and 20 MLs, respectively. These values are much larger than the critical thickness of 8 MLs, above which the finite size effect is negligible \cite{Song020503}. It is also noted that the tensile strain is not relaxed at the top surface of 20-ML FeSe/SrTiO$_3$ as demonstrated by the thickness dependence study \cite{Tan13}. For FeSe/CaF$_2$, it has been reported that the compressive strain is applied to relatively thick films but not to ultrathin films because of the Volmer-Weber-type growth; the reported in-plane lattice constant at $\sim$100 MLs or less is close to that of bulk FeSe, while it becomes smaller by $\sim$1\% at $\sim$150 MLs and this compressive strain does not relax at least until 400 MLs \cite{Nabeshima172602}. By taking into account the thickness independence of the compressive strain over a wide range between 150 and 400 MLs, we have chosen the thickness of 300 MLs to apply compressive strain in the present study. The absence of lattice relaxation at the top surface of our 20-ML FeSe/SrTiO$_3$ and 300-ML FeSe/CaF$_2$ is confirmed by the in-plane lattice constant estimated from the Brillouin-zone size determined by the present ARPES experiments (see Supplemental Material for details of the estimation of the lattice constant by ARPES \cite{Suppl}). These ensure the observed differences between thin films and bulk crystals to be essentially of strain-effect origin. It is also noted that all the samples have been kept non-doped to avoid possible complications from the carrier doping; in this regard, they are distinct from heavily electron-doped FeSe/SrTiO$_3$ \cite{Wang037402, Miyata15, Wen16}. High-quality bulk single crystals ($T_c$ = 8 K) were grown by the KCl and AlCl$_3$ flux method. Details of the sample preparation are described elsewhere \cite{Nabeshima172602, Huynh144516, Miyata15}.

ARPES measurements were performed with a Scienta-Omicron SES2002 electron analyzer with a He discharge lamp ($h\nu$ = 21.218 eV) at Tohoku University. The energy and angular resolutions were set at 12-30 meV and 0.2$^{\circ}$, respectively. The Fermi level ($E_F$) of the samples was referenced to that of a gold film evaporated onto the sample holder.

 \section{III. RESULTS and DISCUSSION}
Figure 1 displays a comparison of Fermi surface (FS) and representative band dispersions for three momentum cuts among strain-free bulk FeSe [Figs. 1(a)-1(c)], compressive-strained FeSe/CaF$_2$ [Figs. 1(d)-1(f)], and tensile-strained FeSe/SrTiO$_3$ [Figs. 1(g)-1(i)]. In the normal state of bulk FeSe, one can recognize a holelike FS originating from a highly dispersive band ($\alpha'$) around the $\Gamma$ point in the Brillouin zone [cut 1 of Fig. 1(b)]. Another holelike band ($\alpha$) with the top of dispersion below the Fermi level ($E_F$) is also recognized. These bands persist in the nematic phase, keeping their energy positions nearly stationary to temperature \cite{Nakayama237001, Shimojima121111, Watson155106, Zhang214503}. On the other hand, the band structure around the {\it M} point is significantly reconstructed by nematicity. In the normal state [cut 2 of Fig. 1(b)], there exist a shallow electronlike band ($\gamma$) and a holelike band ($\varepsilon$), which degenerate at {\it M} to produce a van Hove singularity in the band structure \cite{Nakayama237001}. In the nematic phase [Fig. 1(c)], the $\varepsilon$ band splits into two branches ($\varepsilon_1$ and $\varepsilon_2$) due to the lifting of $d_{xz}/d_{yz}$-orbital degeneracy, and simultaneously the $\gamma$ band shifts upward \cite{Nakayama237001, Shimojima121111, Watson155106, Zhang214503}. As shown in Figs. 1(d)-1(i), the overall band structure of compressive-strained FeSe/CaF$_2$ and tensile-strained FeSe/SrTiO$_3$ is qualitatively similar to that of the bulk counterpart; one can recognize the semimetallic nature of the band structure regardless of temperature as well as a finite splitting between the $\varepsilon_1$ and $\varepsilon_2$ bands at low temperature, indicative of the commonality of the nematic phase for strain-free and strained FeSe. A closer look at Fig. 1 further reveals some quantitative differences in the Fermi wave vector of the $\alpha'$ band [cut 1 of Figs. 1(b), 1(e), and 1(h)], the position of the van Hove singularity [cut 2 of Figs. 1(b), 1(e), and 1(h)], and the location of the $\gamma$, $\varepsilon_1$, and $\varepsilon_2$ bands in the nematic phase [Figs. 1(c), 1(f), and 1(i)].

Another important consequence of the epitaxial strain is the change in the FS topology in the nematic phase. In FeSe/SrTiO$_3$, the electron pocket at {\it M} is absent [Fig. 1(g)] unlike FeSe/CaF$_2$ and bulk FeSe [Figs. 1(a) and 1(d)], and small electron pockets emerge away from {\it M}. This Lifshitz transition is associated with the Dirac-cone-like band dispersion \cite{Kuroki09, Pandey12}, as highlighted by the ``X''-shaped dispersion in FeSe/SrTiO$_3$ [Fig. 1(j), see also Supplemental Material Fig. S2]. As shown in Fig. 1(k), there exist two Dirac cones originating from the hybridization of the $\varepsilon_1$ and $\delta$ bands (note that the intensity of the $\delta$ band is too faint to resolve at high temperatures). When $E_F$ is below the $\varepsilon_1$-band top (case of FeSe/SrTiO$_3$), each Dirac cone creates a small FS pocket slightly away from {\it M}. On the other hand, when $E_F$ is above the $\varepsilon_1$-band top (cases of bulk FeSe and FeSe/CaF$_2$), two Dirac cones are connected with each other to form a single FS.

\begin{figure}
 \includegraphics[width=3.4in]{fig2.eps}
\caption{(color online) Strain-induced evolution of band dispersions. (a) Comparison of near-$E_F$ band dispersions around the $\Gamma$ point in the normal state ({\it T} = 120-190 K) for FeSe/SrTiO$_3$, bulk FeSe, and FeSe/CaF$_2$. (b) and (c), Same as (a), but measured around the {\it M} point at high temperature ({\it T} = 120-190 K) and low temperature ({\it T} = 30 K), respectively. Filled circles were extracted from the peak maxima in ARPES spectra shown in Fig. 1. Open circles for the $\delta$ band were extracted from the second-derivative plots in supplementary Fig. S2. Dashed curves are a guide for the eyes to trace the band dispersions. Green shade in (a)-(c) highlights the top of the $\alpha'$ band, the bottom of the $\gamma$ band, and the energy scale of $\varDelta E_{nem}$, respectively.}
\end{figure}

To see the strain-induced evolution quantitatively, we show in Fig. 2(a) a side-by-side comparison of the band dispersion around $\Gamma$. Obviously, the $\alpha'$ band shifts {\it upward} with decreasing the in-plane lattice constant {\it a}. In contrast, at {\it M} [Fig. 2(b)], the $\gamma$ and $\varepsilon$ bands shift {\it downward} with decreasing {\it a}. Such an opposite band shift indicates that the energy overlap between the $\alpha'$ and $\gamma$ bands in the normal state, defined here as semimetallic band overlap $\varDelta E_{h-e}$, strongly depends on the strength of strain. In fact, $\varDelta E_{h-e}$ is $\sim$70 meV for FeSe/CaF$_2$ while it is $\sim$20 meV for FeSe/SrTiO$_3$. The observed change in $\varDelta E_{h-e}$ is not due to the $k_z$ dispersion but due to the relative energy shift between the $\alpha'$ and $\gamma$ bands \cite{note}. This demonstrates that the increase of tensile strain significantly reduces both the hole- and electron-carrier concentrations while keeping the semimetallicity. As recognized from Fig. 2(c), an energy separation of the $\varepsilon_1$ and $\varepsilon_2$ bands in the nematic state ($\varDelta E_{nem}$), which reflects the strength of nematicity \cite{Nakayama237001, Shimojima121111, Watson155106, Zhang214503}, is enhanced from $\sim$45 meV in FeSe/CaF$_2$ to $>$70 meV in FeSe/SrTiO$_3$ upon increasing {\it a}. This establishes the anti-correlation between nematicity and superconductivity.

\begin{figure}
 \includegraphics[width=3.4in]{fig3.eps}
 \caption{(color online) Relationship between semimetallic band gap and some key physical parameters. (a)-(c) Schematics of normal-state band structure for tensile-strained, strain-free, and compressive-strained FeSe, respectively. Pink dots indicate the position of van Hove singularity at the {\it M} point, while red allows highlight the magnitude of $\varDelta E_{h-e}$, which is the energy difference between the $\alpha'$-band top and the $\gamma$-band bottom. (d) Electronic phase diagram of FeSe as a function of in-plane lattice constant {\it a} \cite{Nabeshima172602, Tan13, Huynh144516}. (e)-(g) Schematics of the band structure in the nematic state for tensile-strained phase, at a critical point of the Lifshitz transition, and for compressive-strained phase, respectively. Yellow dots indicate the position of the Dirac point. (h) Strain dependence of $T_c$ and hole/electron carrier concentration at {\it T} = 30 K plotted as a function of $\varDelta E_{h-e}$. (i) Plot of $\varDelta E_{nem}$ (energy difference between the $\varepsilon_1$ and $\varepsilon_2$ bands at the {\it M} point) as a function of $\varDelta E_{h-e}$.}
\end{figure}

Having established the strain-induced band evolution, the next question is how it is related to superconductivity and nematicity. We propose that these electronic phases are sensitive to the semimetallic band overlap $\varDelta E_{h-e}$. As illustrated in the schematic normal-state band dispersion in Figs. 3(a)-3(c), the $\varDelta E_{h-e}$ gradually increases with decreasing {\it a}. When we estimate the electron- and hole-carrier concentrations from the FS volume in the nematic phase, they have a positive correlation to $\varDelta E_{h-e}$ [Fig. 3(h)]. Intriguingly, the $T_c$ also has a positive correlation, indicating that a larger $\varDelta E_{h-e}$ which results in a larger carrier concentration may be essential for achieving higher $T_c$ in compensated FeSe.

Next we discuss how $\varDelta E_{h-e}$ is related to the strength of nematicity. The small $\varDelta E_{h-e}$ under tensile strain means that the van Hove singularity with high density of states \cite{Nakayama237001} is located close to $E_F$ in the normal state [red dot in Fig. 3(a)]. In this case, the band splitting at {\it M} is highly promoted in the nematic phase, because the larger splitting between the $\varepsilon_1$ and $\varepsilon_2$ bands (i.e., larger $\varDelta E_{nem}$) and the resultant $E_F$-crossing of the $\varepsilon_1$ band [Fig. 3(e)] would cause a larger energy gain to stabilize the nematic state. This is inferred from the negative correlation of $\varDelta E_{nem}$ to $\varDelta E_{h-e}$ [Fig. 3(i)], which is also consistent with a recent theoretical study reporting the Pomeranchuk-type nematic instability at $\varDelta E_{h-e}$ $\sim$ 0 \cite{Jiang15}.

The Lifshitz transition is also related to $\varDelta E_{h-e}$ through the energy balance with $\varDelta E_{nem}$. Due to the anti-correlation between $\varDelta E_{h-e}$ and $\varDelta E_{nem}$ [Fig. 3(i)], there exists a Lifshitz point at which the $\varepsilon_1$-band top just touches $E_F$ in the nematic phase [Fig. 3(f)]. When the $\varepsilon_1$-band top is further pushed upward into the unoccupied region [tensile-strained case; Fig. 3(e)], the Dirac cones are pushed upward as well. Consequently, the Dirac points are situated in the vicinity of $E_F$. In this regard, the nematic phase could be categorized into two regimes: one with Dirac carriers in the tensile-strained region, and the other without Dirac carriers in the strain-free and compressive-strained regions [Fig. 3(d)]. Our observation thus suggests that the electronic properties of the nematic phase can be sensitive to $\varDelta E_{h-e}$, which would make strained FeSe a promising platform to explore unconventional physical properties arising from Dirac carriers.

While the change in carrier concentration may be essential to the $T_c$ variation [Fig. 3(h)], the Lifshitz transition may also play a role in the suppression of superconductivity. It has been suggested that spin fluctuations are involved in the superconducting pairing in bulk FeSe due to the existence of gap nodes \cite{Kasahara16309, Song1410}. While spin fluctuations promote the pair scattering between the hole and electron pockets via {\it Q} = ($\pi$, 0) \cite{Mazin057003, Kuroki087004}, such a scattering can be significantly suppressed in tensile-strained FeSe due to the poorly-defined {\it Q} = ($\pi$, 0) scattering originating from the shrinkage of the hole pocket at $\Gamma$ and the emergence of Dirac-cone-induced small pockets away from {\it M} [Fig. 1(g)]. It is noted that a similar Lifshitz transition around the Brillouin-zone corner has been reported to coincide with the onset of superconductivity in lightly-doped BaFe$_{2-x}$Co$_x$As$_2$ \cite{Liu10}. This supports a close correlation between the observed Lifshitz transition and the suppression of superconductivity under tensile strain.

Our observation also has an important implication to the BCS-BEC crossover \cite{Kasahara16309}. A key parameter of this crossover is an effective Fermi energy ($\varepsilon_F$), an energy separation between $E_F$ and the top (bottom) of hole (electron) pocket. The relatively large $k_{B}T_{c}/\varepsilon_F$ ratio of $\sim$0.2 in bulk FeSe indicates that strain-free FeSe lies in this crossover regime \cite{Kasahara16309}. We found that an increase of $\varepsilon_F$ under compressive strain reduces $k_{B}T_{c}/\varepsilon_F$ to 0.05-0.07 in FeSe/CaF$_2$. This suggests that $T_c$ increases on approaching the conventional BCS side, signifying that the proximity of the BCS-BEC crossover, which could be a possible route to high-$T_c$ superconductivity \cite{Randeria14}, is unlikely to be essential for achieving higher $T_c$ in FeSe. Our result demonstrates that $\varepsilon_F$ can be tuned via the epitaxial strain through the control of $\varDelta E_{h-e}$ and $\varDelta E_{nem}$. Such tunability provides a useful route to investigating the BCS-BEC crossover physics.

Finally, we discuss the origin of the strain-induced band evolution and its implications. The conventional band calculations do not correctly reproduce the magnitude of experimental $\varDelta E_{h-e}$ even in strain-free FeSe; namely $\varDelta E_{h-e}$ is significantly overestimated, as discussed in a previous study \cite{Terashima14}. One possible clue to reconcile with the experiment would be the electronic correlation effect, e.g., the repulsive interband interactions which have been proposed to reduce $\varDelta E_{h-e}$ \cite{Ortenzi09, Hirschfeld11}. In the present ARPES, $\varDelta E_{h-e}$ decreases as the effective mass $m^*/m_e$ of the $\alpha'$ band increases ($m_e$ is the free-electron mass and $m^*/m_e$ varies from 3.3 in FeSe/CaF$_2$ to 3.8 in FeSe/SrTiO$_3$). This suggests a certain role of high-energy interactions in reducing $\varDelta E_{h-e}$. The strain-induced change in the electronic correlation may also play an important role in controlling the electronic phase diagram, since it has been suggested that superconductivity and nematicity are sensitive to the high-energy interactions as well as the low-energy electronic structure in iron-based superconductors due to the moderately correlated nature. In fact, the change in $\varDelta E_{h-e}$ alone cannot explain the strong material dependence of $T_c$, and hence the change in the strength of pairing interactions should be taken into account.

\section{IV. SUMMARY}
We reported high-resolution ARPES results on epitaxial films and bulk crystals of FeSe to elucidate the evolution of the band structure from the tensile-strained non-superconducting regime to the compressive-strained high-$T_c$ regime. We found that the compressive strain gives rise to large semimetallic band overlap $\varDelta E_{h-e}$ and the resultant expansion of hole and electron Fermi surfaces, leading to higher $T_c$. Under tensile strain, small $\varDelta E_{h-e}$ results in the shrinkage of Fermi surfaces and the enhancement of nematic instability, leading to the suppression of superconductivity accompanied by the stabilization of nematicity. We also observed the reduction of mass renormalization factor under compressive strain due to the weakening of electronic correlation, which may be also responsible for large $\varDelta E_{h-e}$ and higher $T_c$. Our finding related to the semimetallic band overlap provides important insights into the origins of superconductivity and nematicity as well as their interplay in iron-based superconductors.

\section{ACKNOWLEDGEMENTS}
We thank Y. Miyata, S. Kanayama, N. Inami, H. Kumigashira, and K. Ono for their assistance in ARPES measurements. This work was supported by grants from the Japan Society for the Promotion of Science (JSPS), KAKENHI (grants No. JP15H02105, No. JP25287079, No. JP15H05853, and No. JP25107003), the Program for Key Interdisciplinary Research, the Murata Science Foundation, and High Energy Accelerator Research Organization, Photon Factory (KEK-PF) (Proposals No. 2012S2-001 and No. 2015S2-003).

\bibliographystyle{prsty}


\section{Introduction}
\label{Intro}\vskip-0.2cm

Computing by employing heat rather than electricity would offer several advantages compared to standard electronic devices, not least that the unavoidable heat dissipation, intrinsically produced in any computational scheme, may be used as an advantage, rather than an hindrance. On a practical side, additional operations may be done without adding extra heat dissipation. Of course, realizing such a concept requires components that have always been attributed to electronics~\cite{Li12,Skl15,Ben16}. Indeed, in the past few years the thermal counterpart of conventional diodes, transistors, memories, and logic elements have been proposed and discussed thus offering a path to a type of unconventional computing with heat~\cite{Li04,Ote10,Li12,Ben13,Ben14,Kub14,Ord16,Ben16}.

In order to process and store information by phononic heat currents, both thermal logic and thermal memories~\cite{Wan07,Wan08,Li12,Skl15} were initially conceived. However, this phonon-based thermal technology suffers from the limited speed of the heat carriers, i.e., the acoustic phonons, which is some orders of magnitude smaller than the speed of electrons. The first solid-state thermal memory was practically demonstrated in Ref.~\cite{Xie11}.
Subsequently, optical architectures for processing and managing information via thermal photons were also developed~\cite{Kub14,Ben14,Elz14,Dya15,Ben15,Son15,Ben16,Ord16,Ito16,Elz16}. Kubytskyi \emph{et al.}~\cite{Kub14} designed a thermal memory based on far-field radiative effects, in which the time to write the memory states was predicted to be on the order of a few seconds.
Instead, a writing time of several orders of magnitude lower, i.e., $\sim5\;\text{ms}$, was expected in near-field thermal memories~\cite{Dya15}.
Alternatively, Elzouka \emph{et al.}~\cite{Elz14,Elz16} proposed a nano thermo-mechanical memory, by varying the coupling between thermal expansion and near-field radiative heat transfer. In this case, the estimated writing time was on the order of $\sim0.05\;\text{s}$. These operational times are several orders of magnitude too long for practical applications such as thermal logic and/or computation.

\begin{figure}[b!!]
\includegraphics[width=\columnwidth]{Fig00.pdf}
\caption{Thermal fluxes in a magnetically driven SQUID formed by two superconductors, $S_1$ and $S_2$, at temperatures $T_1$ and $T_2$, respectively, tunnel coupled through the junctions $J_a$ and $J_b$. The
applied magnetic flux threading the SQUID loop, $\Phi_{ext}$, drives the temperature $T_2$. The heat current, $P_t$, flowing through the junctions depends on the temperatures and the total flux through the SQUID ring. $P_{e-ph}$ represents the coupling between quasiparticles in $S_2$ and the lattice phonons residing at $T_{bath}$, whereas $P_{heat}$ denotes the power injected into $S_1$ through heating probes in order to impose a fixed quasiparticle temperature $T_1$. The arrows indicate the direction of heat currents for $T_1>T_2>T_{bath}$. The temperature $T_2$ is the observable used to define the logic states 0 and 1 of the thermal memory.}
\label{Fig00}
\end{figure}

In this paper, we make a step forward in the panorama of thermal memories by introducing the concept of superconducting (Josephson) thermal memory, which uses the electronic temperature of an inductive superconducting quantum-interference device (SQUID) to define distinct thermal states. We stress that in a superconductor the electronic temperature follows dynamics appreciably faster than the phononics ones~\cite{Gia06}. Recently, the phase-coherent thermal transport~\cite{Gia12,MarGia13,Mar14,ForGia17,For17}, the negative differential thermal conductance~\cite{ForTim16}, and the hysteretical thermal behavior~\cite{Gua17} in temperature-biased Josephson devices were investigated. Here, we discuss the dynamics of a flux-controlled SQUID, delving into the hysteresis of the steady temperatures~\cite{Gua17} to define the logic states 0 and 1 of a memory device. The hysteretic behavior of the SQUID is a straight consequence of the inductive nature of the device~\cite{Cla04}, while the thermal bistability~\cite{Gua17} results from the coherent thermal transport through a temperature-biased SQUID~\cite{Mak65,Gol13}.

The paper is organized as follows. In Sec.~\ref{SQUIDModel}, the theoretical background used to describe the phase evolution of a thermally biased, magnetically driven SQUID with a non-vanishing ring inductance is discussed. In Sec.~\ref{ThermalModel}, the thermal balance equation and the heat currents are introduced. In Sec.~\ref{Results}, the thermal memory states are defined and the behaviour of the thermal memory is explored. The characteristic switching time of the memory and a fast readout scheme are also discussed. In Sec.~\ref{Conclusions}, conclusions are drawn.

\section{SQUID phase dynamics}
\label{SQUIDModel}\vskip-0.2cm

We first explore the dynamics of a magnetically driven inductive SQUID formed by two Josephson junctions (JJs), see Fig.~\ref{Fig00}. To do so, we rely on the resistively and capacitively shunted junction (RCSJ) model, describing the phase evolution of the JJs. According to the RCSJ model, the current flowing through the $i$-th JJ is given by~\cite{Cla04,Bar82,Gra16}
\begin{equation}
I_i=\frac{\Phi_0}{2\pi}C_i\ddot{\varphi_{i}}+\frac{\Phi_0}{2\pi}\frac{1}{R_i}\dot{\varphi_{i}}+I_{c_i}\sin\varphi_i,
\label{RCSJ}
\end{equation}
where $\Phi _0=h/(2e)\simeq 2\times 10^{-15}\;\textup{Wb}$ is the flux quantum ($e$ and $h$ being the electron charge and the Planck constant, respectively), and $C_i$, $R_i$, $I_{c_i}$, and $\varphi_i$ are the capacitance, the normal resistance, the critical current, and the superconducting phase difference across the $i$-th junction, respectively.

The flux quantization in a superconducting ring interspersed with two JJs imposes the constraint
\begin{equation}
\frac{\varphi_1-\varphi_2}{2}=\pi \frac{\Phi}{\Phi_0}+\pi k,
\label{Fluxquantization}
\end{equation}
where $k$ is an integer representing the amount of enclosed flux quanta in the ring. Here,
\begin{equation}
\Phi =\Phi_{ext}-LI_{circ}
\label{TotalFlux}
\end{equation}
 is the total magnetic flux threading the SQUID loop, where $I_{circ}$ is the circulating current and the superconducting ring inductance $L$ has a geometric contribution as well as a kinetic contribution~\cite{Cla04,Ann10}. The system is driven by the externally applied magnetic flux through the ring, $\Phi_{ext}$.

The dynamics of a SQUID formed by two JJs and biased by a current $I_{bias}$ is determined by the following system of equations~\cite{Gra16}
\begin{eqnarray}\label{RCSJa}
\frac{I_{bias}}{2}+I_{circ}=\frac{\hbar }{2e}C_1\ddot{\varphi_{1}}+\frac{\hbar}{2e}\frac{1}{R_1}\dot{\varphi_{1}}+I_{c_1}\sin\varphi_1\\
\frac{I_{bias}}{2}-I_{circ}=\frac{\hbar }{2e}C_2\ddot{\varphi_{2}}+\frac{\hbar}{2e}\frac{1}{R_2}\dot{\varphi_{2}}+I_{c_2}\sin\varphi_2,
\label{RCSJb}
\end{eqnarray}
with the constraints given by Eqs.~\eqref{Fluxquantization} and~\eqref{TotalFlux}.

For the sake of generality, we suppose $C_1\neq C_2$ and $R_1\neq R_2$ and we introduce the quantities $C_\pm=C_1\pm C_2$, $R_\pm=R_1\pm R_2$, and $I_\pm=I_{c_1}\pm I_{c_2}$. The degree of asymmetry of the SQUID is defined as $r=I_{_{-}}/I_{_{+}}=-R_{_{-}}/R_{_{+}}$ (since $I_{c_i}\propto R_i^{-1}$~\cite{Bar82}). Once the values of $r$ and $R_1$ are chosen, the $S_2$ normal resistance can be estimated according to $R_2=R_1(1+r)/(1-r)$.

We recast Eqs.~\eqref{RCSJa} and~\eqref{RCSJb} in terms of the variables $\varphi=(\varphi_1+\varphi_2)/2$ and $\phi=(\varphi_1-\varphi_2)/2$. By the sum of Eqs.~\eqref{RCSJa} and~\eqref{RCSJb}, the following equation results
\begin{equation}
I_{bias}=\frac{\Phi_0}{2\pi}C_{_{+}}\ddot{\varphi}+\frac{\Phi_0}{2\pi}C_{_{-}}\ddot{\phi}+\frac{\Phi_0}{2\pi}\frac{1}{\mathcal{R}}\dot{\varphi}+\frac{\Phi_0}{2\pi}\frac{r}{\mathcal{R}}\dot{\phi}+I_{_{+}}f_r(\varphi,\phi),
\label{EqSQUIDSum}
\end{equation}
where $\mathcal{R}=R_1R_2/R_+$ and $f_r(\varphi,\phi)=\sin\varphi\cos\phi+r \cos\varphi\sin\phi$.

By proper normalization of the time, i.e., $\tau =2\pi \nu t$ with $\nu$ being the driving frequency, Eq.~\eqref{EqSQUIDSum} becomes
\begin{equation}
\mathcal{K}\left (\frac{\partial^2 \varphi }{\partial \tau^2}+\mathcal{C}\frac{\partial^2 \phi }{\partial \tau^2} \right )+\frac{\partial \varphi }{\partial \tau}+r\frac{\partial \phi }{\partial \tau}+\alpha\left [ f_r\left ( \varphi, \phi \right )-\delta \right ]=0,
\label{EqSQUIDSumNorm}
\end{equation}
where $\mathcal{C}=C_{_{-}}/C_{_{+}}$, $\mathcal{K}=2\pi\nu \mathcal{R}C_{_{+}}$, $\alpha =\frac{\mathcal{R} I_{_{+}}}{\nu\Phi_0}$, and $\delta =\frac{I_{bias}}{I_{_{+}}}$.
By subtracting Eqs.~\eqref{RCSJa} and \eqref{RCSJb}, one obtains
\begin{equation}
2I_{circ}=\frac{\Phi_0}{2\pi}C_{_{+}}\left (\mathcal{C}\ddot{\varphi}+\ddot{\phi} \right )+\frac{\Phi_0}{2\pi}\frac{1}{\mathcal{R}}\left (r\dot{\varphi}+\dot{\phi} \right )+I_{_{+}}g_r(\varphi,\phi),
\label{EqSQUIDSub}
\end{equation}
where $g_r(\varphi,\phi)=r\sin\varphi\cos\phi+\cos\varphi\sin\phi$.

From Eq.~\eqref{TotalFlux}, $\phi(t)-\phi_e(t)=-\frac{\beta}{2}\frac{I_{circ}}{I_{c_1}}$, where $\phi_e=\pi \Phi_{ext}/\Phi_0$ and $\beta=\frac{2\pi}{\Phi_0} L I_{c_1}$
is the SQUID hysteresis parameter~\cite{Cla04}. The higher the value of $\beta$, the greater the hysteretical response of the SQUID~\cite{Gua17}. By proper normalization of Eq.~\eqref{EqSQUIDSub} one obtains
\begin{eqnarray}\nonumber
&&\beta \mathcal{K}\left (\mathcal{C}\frac{\partial^2 \varphi }{\partial \tau^2}+\frac{\partial^2 \phi }{\partial \tau^2} \right )+\beta\left ( r\frac{\partial \varphi }{\partial \tau}+\frac{\partial \phi }{\partial \tau} \right )+\\
&&+\alpha\left \{ \beta g_r\left ( \varphi, \phi \right )+2(r+1)\left [ \phi(t)-\phi_e(t) \right ] \right \}=0.
\label{EqSQUIDSubNorm}
\end{eqnarray}

Eqs.~\eqref{EqSQUIDSumNorm} and~\eqref{EqSQUIDSubNorm} have to be solved numerically to study the behavior of the SQUID, when a non vanishing inductance, i.e., $\beta>0$, is taken into account.
Then, the JJ's phases are calculated as $\varphi_1=\varphi+\phi$ and $\varphi_2=\varphi-\phi$.

\emph{The hysteretic parameter}. $-$ The hysteretic parameter $\beta$ is proportional to both the inductance of the superconducting ring, $L$, and the critical current $I_{c_1}$. $L$ is the sum of both a geometric and a kinetic contribution, $L_G$ and $L_{K_{ring}}$, respectively. Moreover, $\beta$ depends on the temperatures through the kinetic inductance and the critical current~\cite{Gia05,Tir08,Bos16}. Specifically,
\begin{eqnarray}
\beta(T_1,T_2)&=&\frac{2\pi}{\Phi_0}L(T_1,T_2)I_{c_1}(T_1,T_2)=\\
&=&\frac{2\pi}{\Phi_0}\left [ L_G+L_{K_{ring}}(T_1,T_2) \right ]I_{c_1}(T_1,T_2).\nonumber
\label{beta_T}
\end{eqnarray}

The BCS expression of the kinetic inductance of a superconducting strip at temperature $T$ is~\cite{Ann10}
\begin{equation}
L_K(T)=R_{sq}\frac{l}{w}\frac{\hbar}{\pi\Delta(T)\tanh\left [ \Delta(T)/(2k_BT) \right ]},
\label{L_Kwire}
\end{equation}
where $R_{sq}$ is the sheet resistance in the non-superconducting state, $l$ and $w$ are the length and the width of the strip, respectively, so that its normal resistance is $R_{strip}=\left ( \frac{l}{w} \right )R_{sq}$.

Therefore, the kinetic inductance of the superconducting ring of the SQUID with arms residing at temperatures $T_1$ and $T_2$ is
\begin{equation}
L_{K_{ring}}(T_1,T_2)= \sum_{j=1,2}R_{L_j}\frac{\hbar}{\pi\Delta(T_j)\tanh\left [ \frac{\Delta(T_j)}{2k_BT_j} \right ]},
\label{L_K}
\end{equation}
where $R_{L_j}=R_{sq}\left ( \frac{l_j}{w_j} \right )$ is the normal resistance of the $j$-th SQUID arm, with $l_j$ and $w_j$ being its length and width, respectively.

We assume a vanishing geometric inductance (i.e., $L(T_1,T_2)\simeq L_{K_{ring}}(T_1,T_2)$), since we can show that it does not crucially affects the overall behaviour. For $R_{L_1}=R_{L_2}=R_{1}$, if the temperatures are $T_1=6.5\;\text{K}$ and $T_2\in[4.2-4.5]\;\text{K}$, one obtains $\beta(T_1,T_2)\simeq1.91$. This is the value set in all the numerical calculations. Anyway, in the presence of an hysteretical behaviour of the SQUID the precise tuning of the parameters is not required to define a reliable thermal memory.

\emph{The $\pi$-swap dynamics}. $-$ We investigate the behaviour of $\varphi$ and $\Phi$, shown in Fig.~\ref{Fig01}(a), of a slightly asymmetric, not-biased SQUID, i.e., $r=0.01$ and $I_{bias}=0$, driven by the magnetic flux $\Phi_{ext}=\Phi_0|\sin(2\pi\nu t)|$, with $\nu=1\;\text{GHz}$.
The exact shape of the driving flux is not essential for the thermal bistability we will discuss.
Here, we assume $C_1=C_2=10\;\text{fF}$ and $R_1=10\;\Omega$, namely, typical values for a Nb/AlOx/Nb junction~\cite{Sol15,Pat99}, so that $R_2\simeq10.2\;\Omega$ and $I_{_{+}}\simeq0.3\;\text{mA}$. In the calculations, temperature-dependent critical currents are taken into account.

The current circulating through the SQUID tends to compensate and screen the applied flux, see Fig.~\ref{Fig01}. In fact, according to the Faraday-Lenz law, as the external flux $\Phi_{ext}$ increases, the total flux $\Phi$ grows less rapidly than $\Phi_{ext}$, since the flux induced by the circulating current opposes $\Phi_{ext}$, see Fig.~\ref{Fig01}(a). However, as $\Phi$ approaches the critical value $\Phi_0/2$, we observe that the total flux $\Phi$ abruptly changes, i.e., a transition $k\to k\pm 1$~\cite{Gua17} takes place, and the phase difference $\varphi$ jumps from $0$ to $\pi$ (or vice versa), namely, a $\pi$-swap of $\varphi$ occurs, see Fig.~\ref{Fig01}(a).
Then, the resulting path of alternate $\varphi$ jumps is related to the asymmetry of the device, i.e., $r>0$, which is necessary since we fix the external bias to $I_{bias}=0$~\cite{Sol15}. Indeed, in order to guarantee no effects of the electric current on the temperature difference, it is appropriate for this device to consider only non-galvanic scheme, where the SQUID is electrically isolated (floating) from external circuit.

\begin{figure}[t!!]
\includegraphics[width=0.5\textwidth]{Fig01.pdf}
\caption{(a) Behaviour of $\Phi/\Phi_0$ and $\varphi/\pi$ as a function of $t$, for three periods of the drive $\Phi_{ext}$. (b) and (c) Time evolution of $P(t)$ and $T_2$, respectively. The values of the other parameters are: $r=0.01$, $I_{bias}=0$, $\beta=1.91$, $R_1=10\Omega$, $C_1=C_2=10 \;\text{fF}$, $\nu=1\;\text{GHz}$, $T_1=6.5\;\text{K}$, $T_{bath}=4.2\;\text{K}$, $\mathcal{V}=10^{-1}{\mu\textup{m}}^3$, $\Sigma=3\times10^9\textup{W}\textup{m}^{-3}\textup{ K}^{-5}$, and $N_F=10^{47}\textup{ J}^{-1}\textup{ m}^{-3}$. }
\label{Fig01}
\end{figure}

\begin{figure}[t!!]
\includegraphics[width=\columnwidth]{Fig02.pdf}
\caption{Time evolutions of the temperature $T_2(t)$ as the external flux is kept fixed to the value $\Phi_{ref}=0.7\Phi_0$, when this value is reached within the first period during the forward (a) and backward (b) sweep of the drive, for $\nu=1\;\text{GHz}$. Dot-dashed lines represent the asymptotic values approached by $T_2$ in the two cases. In both (a) and (b), the top panel shows the driving flux $\Phi_{ext}/\Phi_0$, the dotted curve tracing the function $|\sin(2\pi\nu t)|$. A red (blue) arrow is used to indicate the increasing (decreasing) sweep direction of the drive, when its value is frozen. }
\label{Fig02}
\end{figure}

\section{Thermal dynamics}
\label{ThermalModel}\vskip-0.2cm

We now suppose a thermal bias across the SQUID, see Fig.~\ref{Fig00}. Specifically, the branch $S_1$ resides at a fixed temperature $T_1$, which is maintained by the good thermal contact with heating probes. The electronic temperature $T_2(t)$ of the branch $S_2$ is the key quantity to define our thermal memory, since it floats and can be driven by the external flux $\Phi_{ext}$. It depends on all the energy exchanges in $S_2$ which is in thermal contact with a phonon bath at a temperature $T_{bath}<T_1$. The thermal balance equation for the incoming, i.e., $P_t\left ( T_1,T_2,t \right )$, and outgoing, i.e., $P_{e-ph,2}\left ( T_2,T_{bath}\right )$, thermal powers in $S_2$ can be written as~\cite{GiaMar12}
\begin{equation}
P_t\left ( T_1,T_2,t \right )-P_{e-ph,2}\left ( T_2,T_{bath}\right )=C_v(T_2)\frac{\mathrm{d} T_2}{\mathrm{d} t}.
\label{ThermalBalanceEq}
\end{equation}
The rhs represents the variations of the internal energy of the system, with $C_v(T)=T \,\partial \mathcal{S}/\partial T$ being the heat capacity. $\mathcal{S}(T)$ is the electronic entropy of the superconductor $S_2$ and is given by~\cite{Rab08,Sol16}
\begin{equation}
\mathcal{S}(T)=-4k_BN_F\mathcal{V}\int_{-\infty}^{\infty} f(\varepsilon,T) \log f(\varepsilon,T) \mathcal{N}(\varepsilon,T) d\varepsilon.
\label{Entropy}
\end{equation}
Here, $k_B$ is the Boltzmann constant, $f ( E ,T )=1/\left (1+e^{E/k_BT} \right )$ is the Fermi distribution function, $N_F$ is the density of states at the Fermi energy, $\mathcal{V}$ is the volume of $S_2$, and $\mathcal{N}_j\left ( \varepsilon ,T_j \right )=\left | \text{Re}\left [ \frac{ \varepsilon +i\gamma_j}{\sqrt{(\varepsilon +i\gamma_j) ^2-\Delta _j\left ( T_j \right )^2}} \right ] \right |$ is the smeared BCS density of states of the $j$-th superconductor, where $\Delta_j\left ( T_j \right )$ and $\gamma_j$ are the BCS energy gap and the Dynes broadening parameter~\cite{Dyn78}, respectively.

The heat current $P_t( T_1,T_2,t)$ flowing from $S_1$ to $S_2$ reads
\begin{eqnarray}\label{Pt}\nonumber
P_t=\sum_{i=1,2}&&P_{qp,i}( T_1,T_2,V_i)+\cos\varphi_i P_{\cos,i}( T_1,T_2,V_i)+\\
&&+\sin\varphi_i P_{\sin,i}( T_1,T_2,V_i),
\end{eqnarray}
where $V_i(t)=\Phi_0/(2\pi)\dot{\varphi_i}$ is the voltage drop across the $i$-th JJ. $P_t( T_1,T_2,t)$ depends, through $\varphi_1(t)$ and $\varphi_2(t)$, on the evolution of the driving flux $\Phi_{ext}$.
In the adiabatic regime~\cite{Gol13}, the quasi-particle and the anomalous heat currents, $P_{qp,i}$, $P_{\cos,i}$, and $P_{\sin,i}$ read, respectively,~\cite{Mak65,Gol13,VirSol17}
\begin{eqnarray}\label{Pqp}\nonumber
&&P_{qp,i}(T_1,T_2,V_i )=\frac{1}{e^2R_i}\int_{-\infty}^{\infty} d\varepsilon \mathcal{N}_1 ( \varepsilon-eV_i ,T_1 )\mathcal{N}_2 ( \varepsilon ,T_2 )\\
&&\times(\varepsilon-eV_i) [ f ( \varepsilon-eV_i ,T_1 ) -f ( \varepsilon ,T_2 ) ],\\ \nonumber
&&P_{\cos,i}( T_1,T_2,V_i )=-\frac{1}{e^2R_i}\int_{-\infty}^{\infty} d\varepsilon \mathcal{N}_1 ( \varepsilon-eV_i ,T_1 )\mathcal{N}_2 ( \varepsilon ,T_2 ) \\
&&\times\frac{\Delta_1(T_1)\Delta_2(T_2)}{\varepsilon}[ f ( \varepsilon-eV_i ,T_1 ) -f ( \varepsilon ,T_2 ) ],\label{Pcos}\\\nonumber
&&P_{\sin,i}(T_1,T_2,V_i)=\frac{eV_i}{2\pi e^2R_i}\iint_{-\infty}^{\infty} d\epsilon_1d\epsilon_2 \frac{\Delta_1(T_1)\Delta_2(T_2)}{E_2}\\
&&\times\left [\frac{1-f(E_1,T_1)-f(E_2,T_2)}{\left ( E_1+E_2 \right )^2-e^2V_i^2}\text{+}\frac{f(E_1,T_1)-f(E_2,T_2)}{\left ( E_1-E_2 \right )^2-e^2V_i^2}\right ]
\end{eqnarray}
with $E_j=\sqrt{\epsilon_j^2+\Delta_j(T_j)^2}$ the Bogoliubov energies.
We observe that the anomalous terms in $P_t$ depend crucially on the junctions phases and are the terms which connect the thermal conduction of the JJs with the phase dynamics.

In Eq.~\eqref{ThermalBalanceEq}, $P_{e-ph,2}$ represents the power loss by the quasiparticles in $S_2$ into the lattice phonons residing at $T_{bath}$~\cite{Pek09}
\begin{eqnarray}\label{Qe-ph}\nonumber
P_{e-ph,2}&=&\frac{-\Sigma \mathcal{V} }{96\zeta(5)k_B^5}\int_{-\infty }^{\infty}dEE\int_{-\infty }^{\infty}d\varepsilon \varepsilon^2\textup{sign}(\varepsilon)M_{_{E,E+\varepsilon}}\\\nonumber
&\times& \Bigg\{ \coth\left ( \frac{\varepsilon }{2k_BT_{bath}}\right ) \left [ \mathcal{F}(E,T_2)-\mathcal{F}(E+\varepsilon,T_2) \right ]\\
&-&\mathcal{F}(E,T_2)\mathcal{F}(E+\varepsilon,T_2)+1 \Bigg\},
\end{eqnarray}
where $M_{E,{E}'}=\mathcal{N}_i(E,T_2)\mathcal{N}_i({E}',T_2)\left [ 1-\Delta ^2(T_2)/(E{E}') \right ]$, $\mathcal{F}\left ( \varepsilon ,T_2 \right )=\tanh\left [ \varepsilon/(2 k_B T_2) \right ]$, $\Sigma$ is the electron-phonon coupling constant, and $\zeta$ is the Riemann zeta function. Hereafter, we impose $\mathcal{V}=10^{-1}{\mu\textup{m}}^3$, $N_F=10^{47}\textup{ J}^{-1}\textup{ m}^{-3}$, $\gamma_1=\gamma_2=10^{-4}\Delta_2(0)$, and $\Delta_1(0)=\Delta_2(0)=1.764k_BT_c$, with $T_c=9.2\;\text{K}$, namely, typical values for an Nb-based SQUID.

The temperature $T_2(t)$ is obtained by solving Eqs.~\eqref{EqSQUIDSumNorm},~\eqref{EqSQUIDSubNorm}, and~\eqref{ThermalBalanceEq} for fixed values of $T_1$ and $T_{bath}$, and represents the observable we use to encode the 0 and 1 logic states of the thermal memory.

\section{Results}
\label{Results}\vskip-0.2cm

We impose that the bath resides at $T_{bath}=4.2\;\text{K}$, and $S_1$ is at a temperature $T_1=6.5\;\text{K}$ kept fixed throughout the computation. However, the thermal logic we are going to discuss is robust to moderate fluctuations of $T_1$. The time evolution of both the power $P_t( T_1,T_2,t)$ injected into $S_2$ and $T_2(t)$ within three periods of the drive $\Phi_{ext}=\Phi_0|\sin(2\pi\nu t)|$ is shown in Figs.~\ref{Fig01}(b) and~\ref{Fig01}(c), respectively. The phase dependency of $P_t$ is clearly visible since it increases and abruptly falls in correspondence of the $\pi$-swaps of $\varphi$.
The temperature $T_2(t)$ increases from $T_2(0)=T_{bath}$ around a quasi-equilibrium value~\cite{Gua17} determined by Eq.~\eqref{ThermalBalanceEq}, with peaks in correspondence of the $P_t$'s jumps.

\begin{figure}[t!!]
\includegraphics[width=\columnwidth]{Fig03.pdf}
\caption{(a) Driving flux giving an initial memory state 1 and four following switches, the latter indicated in all panels by blue dashed lines. A red (blue) arrow is used to indicate the increasing (decreasing) sweep direction of the drive when its value is frozen to the value $\Phi_{ref}=0.7\Phi_0$. (b) Time evolution of the temperature $T_2(t)$. Dot-dashed lines represent the asymptotic values approached by $T_2$ in the two cases. (c) Time evolution of the effective inductance, $\widetilde{L}_{T}$, see Eq.~\eqref{TankeffectiveInductance}, of the tank circuit coupled to the SQUID used for the memory state readout, as $L_T=100\;\text{pH}$.}
\label{Fig03}
\end{figure}

We discuss now how to drive the SQUID to control the thermal memory.
We assume a drive $\Phi_{ext}(t\leq t_{ref})=\Phi_0|\sin(2\pi\nu t)|$ and $\Phi_{ext}(t> t_{ref})=\Phi_{ref}$, with $t_{ref}$ being the time at which $\Phi_{ext}(t_{ref})=\Phi_{ref}$ (i.e., the driving flux is frozen to the value $\Phi_{ref}$ for $t\geq t_{ref}$). Time-dependent drives behaving in this manner are shown as solid lines in top panels of Figs.~\ref{Fig02}(a) and~\ref{Fig02}(b). Here, the dotted line represents the function $|\sin(2\pi\nu t)|$ and the red (blue) arrow indicates the increasing (decreasing) sweep direction of the drive when its value is fixed. Let us set, for instance, $\Phi_{ref}=0.7\Phi_0$. As is clearly shown in Figs.~\ref{Fig01} and~\ref{Fig02}, within each drive period, the condition $\Phi_{ext}=\Phi_{ref}=0.7\Phi_0$ occurs twice, once preceding and once following a phase swap.

The temperature $T_2$ evolves differently according to the increasing, i.e., forward, or decreasing, i.e., backward, sweep direction of the driving flux when its value is kept fixed to $\Phi_{ref}$, as shown in the two driving protocols in Figs.~\ref{Fig02}(a) and~\ref{Fig02}(b), respectively. Although the driving flux finally assumes the same value, in the two cases $T_2$ tends to rapidly converge to different steady temperatures, $T_{2,\bit{1}}\simeq4.5\;\text{K}$ and $T_{2,\bit{0}}\simeq4.4\;\text{K}$ (indicated by two horizontal dot-dashed lines in Fig.~\ref{Fig02}). Hereafter, a subscript enclosed in a rectangle indicates the logical state associated to a specific value of an observable.
Then, in an inductive SQUID the phase bistability, i.e., the $\pi$-swap of $\varphi$, reflects on a temperature bistability. Accordingly, we suggest a thermal memory in which the logic states 0 and 1 are defined by these two distinct values of the electronic temperature $T_2$. The writing operation of these states is performed through the driving flux.

The switch between the logical thermal memory states can be done over a short timescale, as shown in Fig~\ref{Fig03}, by driving the system through a $\pi$-swap and fixing it again as the $\Phi_{ref}$ value is reached anew. Accordingly, a memory state switch corresponds to a change in the sweep direction of the driving flux. Fig.~\ref{Fig03}(a) shows a magnetic drive giving an initial state 1 and four subsequent switches, the latter indicated by blue dashed lines.
The temperature $T_2(t)$ increases from the value $T_2(0)=T_{bath}$, and the logic state 1 is rapidly reached, i.e., $T_{2}\simeq4.5\;\text{K}$, see Fig.~\ref{Fig03}(b). After each switch, the temperature $T_2(t)$ follows a transient regime and then, as the driving flux is kept fixed again, it exponentially approaches the steady temperature which is distinctive of the succeeding logic state, see Fig.~\ref{Fig03}(b).

\emph{Memory switching time}. $-$ The time of the memory writing operation can be evaluated as the characteristic time of this exponential process, namely, $\tau_{wr}\sim0.2\;\text{ns}$. Markedly, a quite good estimate of this time results from first-order expanding the heat current terms in Eq.~\eqref{ThermalBalanceEq} (see Appendix~\ref{AppA}). In doing so, it derives from Eq.~\eqref{ThermalBalanceEq} that $C_v\frac{\mathrm{d} \Delta T_2}{\mathrm{d} t}=(K-G)\Delta T_2$ (see Appendix~\ref{AppA}), where $G$ and $K$ are the electron-phonon~\cite{Vir17} and electron~\cite{Mar14} thermal conductances of the JJs, respectively, and $\Delta T_2(t)$ represent the distance between $T_2(t)$ and its steady value, $T_{2_s}$. Therefore, $T_2$ exponentially approaches $T_{2_s}$, and the characteristic time of this process, i.e., the memory switching time, is $\tau_{sw} =C_v/(G-K)$. For $T_{2_s}=4.5\;\text{K}$, we obtain $\tau_{sw}\sim0.1\;\text{ns}$.

\emph{The readout}. $-$ Recently, it has been shown that that real time temperature measurements in superconducting circuits at nanosecond scale are possible~\cite{Zgi17,Wan17}.
Due to these technical achievements, the readout of the memory state can be done directly with time-dependent calorimetric measurements~\cite{Pek13,Gas15,Gia15,Sai16,Zgi17,Wan17}.
Hereafter, we do not specifically address those proposals and their applicability to the proposed thermal memory, but we suggest an alternative non-invasive indirect procedure. This readout scheme is based on the measurement of the effective inductance $\widetilde{L}_T$ of a tank circuit inductively coupled to our electrically-open double junction SQUID, as depicted in the equivalent electrical circuit schematically shown in Fig.~\ref{Fig07}. The state of the system can be measured via the effective resonance frequency of the tank circuit $\widetilde{f}_T=1\Big /\left ( 2\pi\sqrt{\widetilde{L}_T C_T} \right )$, where $\widetilde{L}_T$ is the effective inductance of the tank circuit (see Appendix~\ref{AppC}) and reads
 %
\begin{equation}
\widetilde{L}_T=L_T\left [1-\textsc{k}^2\frac{L}{L_J(\varphi_1,\varphi_2)+L} \right ],
\label{TankeffectiveInductance}
\end{equation}
where $\textsc{k}$ is the coupling coefficient defined by $M^2=\textsc{k}^2L_TL$, $M$ is the mutual inductance of the system, $C_T$ and $L_T$ are the capacitance and the inductance of the tank circuit, respectively. In Eq.~\eqref{TankeffectiveInductance}, $L$ is the ring inductance [see Eq.~\eqref{beta_T}], $L_J(\varphi_1,\varphi_2)=4L_{\varphi_1}L_{\varphi_2}/(L_{\varphi_1}+L_{\varphi_2})$ is the Josephson contribution, with
\begin{equation}
L_{\varphi_i}(t)=\frac{\Phi_0}{2\pi}\frac{1}{I_{c,i } (T_1,T_2(t))\cos\left [\varphi_i (t) \right ]}
\label{JJKinInductance}
\end{equation}
being the Josephson inductance of the $i$-th JJ.
Assuming $C_T=1\;\text{nF}$, $L_T=100\;\text{pH}$, and $\textsc{k}=0.9$, we obtain for the logic state 1 the steady value $\widetilde{L}_{T,\bit{1}}\simeq83\;\text{pH}$ and a corresponding resonant frequency $\widetilde{f}_{T,\bit{1}}\simeq0.55\;\text{GHz}$, while for the state 0 we obtain $\widetilde{L}_{T,\bit{0}}\simeq63\;\text{pH}$ and $\widetilde{f}_{T,\bit{0}}\simeq0.63\;\text{GHz}$, see Fig.~\ref{Fig03}(c). A cavity with a modest quality factor, $Q\sim10$, should be able to resolve these two memory states, even if it is worthwhile to note that this $Q$ value needs to be increased considerably when the coupling term $\textsc{k}$ is reduced.
In this scheme, the total flux $\Phi_{tot}$ through the SQUID is the sum of the flux $\Phi$, see Eq.\eqref{TotalFlux}, and the probing flux $\Phi_{T,ac}$ due to the tank circuit. For a proper operating point $\Phi_{ref}$ of the external flux, the oscillating component $\Phi_{T,ac}$ should be kept sufficiently small to avoid unwanted memory switches.
We stress that the main contribution to the difference in the effective inductance between the two thermal states is given by the Josephson terms, see Eq.~\eqref{JJKinInductance}, and increases with increasing the temperature difference between the memory states.
Markedly, this detection scheme permits the tuning of the tank circuit effective resonant frequency by tuning the tank circuit parameters.
Finally, the $L_{\varphi}$'s values can be increased, and therefore the visibility in the dispersive mode enhanced, by changing the critical current, namely, by increasing the JJ's resistance $R_j$. Nevertheless, the higher the $R_j$ values, the lower the temperature difference between the two logic states, so an optimal point has to be found.

Notably, the thermal memory shows remarkable robustness against environmental disturbance (see Appendix~\ref{AppB}). In fact, our numerical calculations show that thermal Johnson–Nyquist noise currents~\cite{Bar82} only slightly affects the steady values of the temperatures, while the overall thermal dynamics remains essentially unchanged also in a stochastic Langevin approach.

\begin{figure}[t!!]
\centering
\includegraphics[width=\columnwidth]{Fig07.pdf}
\caption{Thermally biased SQUID (on the left) inductively coupled to a resonant tank circuit (on the right) characterized by an inductance, a resistance, and a capacitance $L_T$, $R_T$, and $C_T$, respectively. The total flux through the SQUID is $\Phi_{tot}=\Phi+\Phi_{T,ac}$ and $M$ is the mutual inductance of the system. }
\label{Fig07}
\end{figure}

\section{Conclusions}
\label{Conclusions}\vskip-0.2cm

In conclusion, we have discussed the bistable thermal behavior of a flux-controlled inductive SQUID and designed a superconducting fast thermal memory based on this effect. We propose a feasible fast readout scheme based on the reading of the effective resonance frequency of a tank circuit inductively coupled to the SQUID. The writing operation of the thermal memory states is performed through the driving flux. Interestingly, a writing time of the memory states on the order of $\sim0.2\;\text{ns}$ is in principle achievable. This is at least seven orders of magnitude faster than any other thermal memory proposed up to now.

The proposed memory is well-placed in the context of superconducting memory elements~\cite{Mat80,Lik12,Rya12,Peo14,Sal17,GuaSol17,Mur17}, significantly pushes forward the thermal memory concept and has a strong relevance in the context of thermal devices.

In addition to the memory encoding, this device will use the heat produced by superconducting electronic circuits as a consequence of the computation to encode memory functionalities. Then, the proposed memory effectively paves the way for a new generation of fast thermal technology, including fundamental devices~\cite{ForGia17} and logic gates~\cite{Pao17}.

We acknowledge P. Virtanen, F. Paolucci, and G. Marchegiani for fruitful discussions.
C.G. and F.G. acknowledges the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013)/ERC Grant agreement No.~615187-COMANCHE for partial financial support.
C.G. and P.S. have received funding from the European Union FP7/2007-2013 under REA Grant agreement No. 630925 -- COHEAT and from MIUR-FIRB2013 -- Project Coca (Grant No.~RBFR1379UX).
A.B. acknowledges the Italian’s MIUR-FIRB 2012 via the HybridNanoDev project under Grant no. RBFR1236VV and CNR-CONICET cooperation programme ``Energy conversion in quantum nanoscale hybrid devices''.
M.D.V. acknowledges support from Department of Energy under Grant No.~DE-FG02-05ER46204.


\section{Introduction}
	In the history of mathematics and computation, the mathematics has often been (far) ahead of the engineering.  In some sense, this is necessary as one won't always know what sort of machine will fit one's needs for specific computations.  Currently (as of 2021) this is certainly the case in quantum computation.  There are, at present, at least three well-known paradigms of computation using quantum mechanical effects.
	\begin{itemize}
		\item The gate model which in some sense attempts to analogize classical computations from a quantum mechanical framework.  This was suggested as early as 1980 \cite{Manin}
		\item Quantum annealing which mimics classical simulated annealing by leveraging quantum tunneling instead of ``temperature" to jump out of local minima.  Classical simulated annealing has been around for quite some time as it is in the suite of Markov chain Monte Carlo methods, while quantum annealing has been around since at least the late 1980s \cite{FACB}.  While quantum annealing is not yet known to be a universal method of computation, it is well suited for optimization problems and problems which rely on random sampling.  Not universal, but still a wide berth of problems and techniques to tackle.  Additionally, since quantum annealing is more specific its mathematics and engineering are further along than the gate model.
		\item Random walks on graphs have found their applications in most fields of modern science and economics.  In the early 2000s, the ``quantization" of random walks lead us to several different paradigms of quantum random walks on graphs, Continuous, Discrete, Coined, and Real \cite{HM, DB, PRR,VBB}.  However, in 2009 Quantum Walks on Graphs were shown to be a universal method of computation \cite{C}.  Despite this, the economic investment in quantum computation by walks on graphs has remained minimal in comparison to the other two paradigms.
	\end{itemize}

	In addition to these well-known paradigms, there are quantum inspired algorithms, noisy intermediate scale computers, hybridized algorithms, and heuristic algorithms.  Among the effects quantum computing has on the world in a real sense is that it has sparked a bit of an arms race to see who can produce a faster probabilistic program, and who can leverage quantum effects without needing the sophisticated hardware that is a ``quantum computer."  Because of this, the true natures of the complexity classes \textbf{BQP} and \textbf{BPP} are begin slowly revealed.

	In this work we seek to add a technique into the canon of quantum walks on graphs.  In this work we consider coined discrete time quantum walks on regular graphs.  In this scenario we simultaneously consider two Hilbert spaces, $\mathcal{H}_\mathcal{P}$ (the position space) and $\mathcal{H}_C$ (the coin space).  One can consider a coined walk on a graph which is not regular, but the computation is far less efficient.  In our case the term ``coin" is an analogy to the classical random walk on a line wherein the walker flips a coin and moves accoding to which face of the coin is showing.  In the situation of a regular graph, each vertex has $d$ edges from which to choose, and thus the quantum ``walker" flips a $d$ sided quantum ``coin."  The quantum walker does not move, but instead allows the wave function to propagate across the graph.  The computational efficiency of this approach comes from the fact that we can encode the dynamics of the system in a (complex) matrix of size $Nd \times Nd$.  More important is that this matrix is also unitary.  So our goal the becomes finding a discrete path in $U(Nd)$.

	This work is split in the following way.  In \S1 we discuss some of the technicalities of a discrete time quantum walk on a regular graph as well as some graph properties necessary to ensure our walk will converge to the desired probability distribution.  We also discuss the parameters we can adjust in a DTQW, which are plentiful, and which
	In \S2 we approach the problem of learning by using the Maximal Torus in $U(d)$ to parameterize our walking space.  We modify the classical algorithm of gradient descent to find our appropriate coins.  It should be noted, that gradient descent is still gradient descent, but we give a very simple implementation by hand which should be understandable to any student who has taken the first semester of calculus and paid attention through the first two weeks of vector calculus.  We approach the problem of trying to fit our target distribution with a single coin, which leaves us with an extremely under determined system and allows us to understand how much we can leverage wave interference to our advantage. In\S3 we use the basis of the Lie algebra $\mathfrak{su}(d)$ which is $d^2-1$ dimensional, thus we essentially square our number of generators.  While for large graphs, and especially ones with good spectral expansion (i.e. the bucky ball, 60 vertices, 3 regular) this is still drastically under determined, we have a better shot and landing on our target distribution.  Consider, for example, a 100 vertex, 10 regular graph.  The quantum walk will have 99 basis elements, for fitting a 100 point vector, so chances are quite good that we can properly fit our distribution. In \S4 and \S5 we give the results and motive future work.  In the end, this particular technique is not extremely useful on its own.  Its utility comes from the fact that we know that there is a plethora of applications for learning a distribution.

	\section{Some Technicalities of Discrete Time Quantum Walks}
	For the sake of consistency, we will use 1-indexing throughout so that our sums appear as $\sum_{i=1}^n$ rather than $\sum_{i=0}^{n-1}$.  The necessary adjustments for 0-indexed languages have been made by the author in the coded implementation of this algorithm and one just needs to pay careful attention throughout to avoid unnecessary headaches with indexing.

	A DTQW on a regular graph is an orbit of a unit vector in $\mathbb{C}^{N}$ via a unitary matrix.  The dimension $N$ is the product of the dimension of two smaller Hilbert spaces.
	\begin{defn}
	Let $G=(V,E)$ be a regular graph with regularity $d$.  Then the two Hilbert spaces one must consider for a DTQW on $G$ are the \emph{position space} $\mathcal{H}_\mathcal{P}$ and the \emph{coin space} $\mathcal{H}_C$.
	\begin{eqnarray}
	\mathcal{H}_\mathcal{P} &=& \text{span} \{ |1 \rangle, \dots ||V|\rangle\}\\
	\nonumber \mathcal{H}_C &=& \text{span} \{|1\rangle,\dots, |d\rangle \}
	\end{eqnarray}
	\end{defn}

	Thus our space $\mathbb{C}^N = \mathbb{C}^{d|V|}$.  We construct a unitary
	\[
	\hat{U} : \mathcal{H}_C \otimes \mathcal{H}_\mathcal{P} \rightarrow \mathcal{H}_C \otimes \mathcal{H}_\mathcal{P}
	\]
	as the product of two other unitaries $\hat{S}$ and $\hat{C}$.

	\begin{defn}
		The two unitary operator $\hat{S}$ and $\hat{C}$ are defined as follows.  For the ``coin operator" $\hat{C}$ we first pick a unitary matrix $C_0$ from $U(d)$ which could be at random or through a parametrization
		and we define
		\begin{equation}
		\hat{C} = C_0 \otimes I_{|V|}
		\end{equation}
		The operator $\hat{S}$ called the ``shifting operator," is built by analogy to a classical random walk.  In a classical random walk one flips a coin, and then moves accordingly.  In the quantum walk, one projects onto a ``con flip" and moves according the the graph structure.
		\begin{equation}
		\hat{S} = \sum_{j=1}^d \sum_{v\in V} |c_j\rangle \langle c_j| \otimes |w\rangle \langle v|
		\end{equation}
		Where $w$ is the $j^{th}$ adjacency to $v$.
	\end{defn}
	\begin{rem}
		One must be careful to respect the graph structure so that $\hat{S}$ remains unitary.  This is done by reading the adjacencies from a consistently labeled rotation map cf\cite{A}
	\end{rem}

	Our discrete time quantum walk can now succinctly written as
	\[
	|\psi_t\rangle = \hat{U}^t |\psi_0\rangle
	\]

	The to compute the probability that our quantum walker will collapse into a classical walker at vertex $i$ we have
	\begin{equation}
	Prob(\psi_t \text{ is at vertex } i) = \sum_{j=1}^d |\langle c_j|\otimes \langle i|\psi)t\rangle|^2
	\end{equation}

	Now if one is observing closely, there are a lot of parameters in play here.  Let's list some of them
	\begin{itemize}
		\item In $\hat{C}$ we have a choice of any unitary matrix.  That gives us $d^2-1$ dimensions of freedom
		\item In $\hat{S}$ we have a  choice of a consistent edge labeling of our graph.  These are not unique, in fact for graphs with a high degree of regularity the number of consistent edge labeling grows algebraically (tight bounds on this growth are not currently known to the author)
		\item We have a choice of starting position.  Generally this is taken as $|1\rangle_{\mathcal{P}}$. This is not necessary, the starting position can be any unit length vector in $\mathbb{C}^{|V|}$.
		\item We have a choice of a starting coin.  This is generally taken as the \emph{Fourier coin}
		\begin{equation}
		|c\rangle = \sum_j e^{2\pi i(j-1)/d}|j\rangle_C
		\end{equation}
		Again, this is not required, it is just a convention, giving us another $d-1$ dimensions of freedom.
		\item Finally, we have the evolution time of $\hat{U}$.  This can be any integer, positive or negative.
	\end{itemize}

    With this many degrees of freedom, it's almost a wonder the DTQW have any utility.  But this is where the difficult mathematics comes in.  For our purposes we will set the number of steps in our evolution ahead of time.  This will be set somewhere near $2*$ diameter$(G)$.  We will also keep the conventions of starting at the first vertex in the graph with a Fourier coin.  While the shifting operator does have some degree of freedom, we will choose a rather simple edge labeling.  In fact, we will choose the first consistent edge labeling we can solve.  As long as we continue with this edge labeling, the operators will remain unitary.  This leaves us only with coin.  Therefore for this work we will stick to ``training" the coin to produce the quantum walk we want.

    One other thing of which we must be aware is that in the situation of a bipartite graph, one cannot guarantee certain distributions.  Therefore we have two potential fixes for this.  The first is to make sure our graph has an odd length cycle.  If this is not available, then we must start with our quantum walker in a superposition of each part of the bipartite graph.  Consider the square ($C_4$) for example.  If we start at $|1\rangle$ and evolve for an odd number of steps, we will not land on $|1\rangle$.  Thus we may start with the distribution $(|1\rangle + |2\rangle)/\sqrt{2}$ or with the uniform distribution.

	\section{Method One: Learning the Angles in a Torus}
	\subsection{A Modified Gradient Descent}
	For the coin operator in DTQW we need to choose some parametrization.  Since the coin operator is of the form $C\otimes I_{|V|}$ we only need to choose a parametrization of the space $U(d)$.  This parametrization need not be topologically dense, but it must simply allow us to move around ``enough" in the space. For the case of a single coin, we have chosen to parametrize the maximal torus of $U(d)$.  This alone does not, however, make a good quantum walking matrix as it is a diagonal matrix and will simply produce a phase shift at each vertex without affecting the probabilities.  On the other hand, the Fourier matrix tends to disperse waves quite well, but is itself only an order 4 operator ($\mathcal{F}^4 = I$) thus we have chosen the conjugate an element of the maximal torus by a Fourier Matrix
	\begin{equation}
	C = \mathcal{F}^* \begin{bmatrix}
	e^{i\theta_1} & & &\\
	& e^{i\theta_2} & &\\
	& & \ddots & \\
	& & & e^{i\theta_d}
	\end{bmatrix} \mathcal{F}
	\end{equation}

	This leaves us with a very tiny amount of parameters to optimize. We shall consider our parameter set as a vector $\vec{\theta}$ and we define the gradient descent algorithm as normal
	\begin{equation}
	\vec{\theta}_{n+1} = \vec{\theta}_n - \gamma \nabla_{\theta} E(\vec{\theta}_n)
	\end{equation}

	where $\gamma$ is the learning rate and $E$ is some error dependent on $\vec{\theta}$.

	None of this is terribly surprising.  Where things get difficult is in actually computing the gradient.  Since we are trying to learn a probability distribution $\vec{\pi}$ we define the error as
	\begin{equation}
	E(\vec{\theta}) = \left(\sum (|\langle j | \hat{U}^t |\psi_0\rangle|^2 - \pi_j)^2\right)^{1/2}
	\end{equation}

	Again nothing terribly surprising.  The difficulty is in trying to take the derivative of this expression with respect to $\theta_j$

	\begin{equation}
	\frac{\partial E}{\partial \theta_j} = \text{ a big mess }
	\end{equation}

	This computation is feasible by hand (as the author has performed it several times), but computationally it's not worth the effort for the mathematical rigor and precision.  One must remember that we are trying to approximate the angles with floats to within some tolerance.  Computers cannot handle mathematical exactness in memory, however, computers can approximate extremely fast, and so we have chosen to take advantage of speed in both simplicitiy of code and ease of reading. We set some tolerance $\varepsilon = 0.01$ and we compute
	\begin{equation}
	\frac{\partial E}{\partial \theta_j} \approx \frac{E(\vec{\theta} + \varepsilon \hat{e}_j) - E(\vec{\theta} - \varepsilon \hat{e}_j)}{ 2\varepsilon}
	\end{equation}

	It is, however, an approximation that we seek, and this gradient can be computed in a few lines of code rather than the several lines it takes to write out the true partial derivatives.
	\begin{rem}
		One of the several difficulties in computing the true partial derivatives is that when considering a matrix
		\[
		\hat{U}^t = \hat{S}\hat{C}\cdots \hat{S}\hat{C}
		\]
		The standard Leibniz Rule from Calculus I fails since $[\hat{S},\hat{C}] \ne 0$.  There are higher order noncommutative effects that must be handled carefully.  If we choose evenly a moderately small $t$ (20 for example) then the formaula for the true partial derivative grows in length to several pages handwritten and grows to unreadable in computer code.
	\end{rem}

	\subsection{A Single Coin}
	In the method of using a single coin in the maximal torus, we have restricted ourselves to exactly $d+1$ parameters.  The $d$ angles along the torus, plus the number of steps we wish to take.  Classical feed forward neural networks have far more parameters than the desired probability distribution.  So our goal is to control the waves emanating from the first vertex in such a way that they constructively and destructively interfere with each other to land on precisely the probability distribution we wish.  It is in this sense that we are leveraging wave interference to take over some of the parameter tuning.

	In \S5 we will look at some results from using a single coin and how preselecting the number of steps in the walk affects the training.  Suffice it to say, using a small number of steps is not sufficient to properly train a DTQW to learn an arbitrary distribution.  One instead needs to solve the equation of picking angles so that $\hat{U}$ has a topologically dense orbit in $\mathbb{S}^{2d|V|-1}$ and then solving the number of steps to get withing a desired tolerance when projecting down to $\mathbb{S}^{2|V|-1}$.

	\subsection{A Single Coin per Step}

	The approach that tends to work empirically is allowing oneself to pick a different coin at every step.  Intuitively this makes much more sense.  One moves a little closer to the desired distribution and then needs to change direction slightly so that in the end one form a discrete path in $\mathbb{S}^{2d|V|-1}$

	We will see again in \S5 how this approach differs from that of allowing only a single coin.  The technical difficult here is training so many more parameters.  It mertis mention that for graphs with good spectral expansion properties, $d$ is relatively small versus $|V|$ and thus using several steps each with $d$ parameters still goes far beneath a general set up in a feed forward neural network.

	Depending on one's computational power, choosing a higher number of steps is preferable as in some experiments the author has seen the last coin or two coins have parameters of 0.  This means, the walk can converge in fewer steps, and the coin matrix, given the current parametrization is simply $I_d\otimes I_{|V|}$.

	\section{Method Two: Learning the Parameters in a Lie Algebra}
	\subsection{A Single Set of Parameters}

	As we mentioned in the first method, we have a severely under determined system, but we used only the Maximal Torus which is dimension $d$.  In the second method we parameterize the whole of $\mathfrak{su}(d)$ by building a version of the Gell-Mann matrices.  While properly the Gell-Mann matrices are the basis of $\mathfrak{su}(3)$ one can extend the basis to include $d$ dimensions.  The basis of $\mathfrak{su}(d)$ is $d^2-1$ dimensions and so we have a significant advantage in using a single coin.  In the computions the author did not use the exact Gell-Mann matrices, but a similar set which is easy to build.  Let's take a small detour to give the construction.

	\subsubsection[basis matrices]{The Basis Matrices of $\mathfrak{su}(d)$}
	We know that the $\mathfrak{su}(d)$ consists of skew-hermitian matrices.  Defining the matrices
	\[E_{ij} = \hat{e}_i\hat{e}_j^t\]
	That is, a matrix of zeros with a 1 in row $i$ column $j$.  Then for $i < j$ we have skew-hermitian matrices
	\[
	E_{ij} + E_{ji} \text{ and } -iE_{ij} + iE_{ji}
	\]

	This gives us $\binom{d}{2} + \binom{d}{2}$ matrices.  We also have the diagonal matrices $E_{ii} - E_{i+1,i+1}$ where $i$ runs from 1 to $d-1$.  This gives us $2\binom{d}{2} + d-1 = d^2-1$ matrices.  It should be noted that these matrices do not satisfy the Lie algebra structure constants in any canonical way, but they do form a basis.  They are also easier to compute and program.  Labeling these matrices as $\lambda_1, \dots, \lambda_n$  we see that our coin can be written
	\begin{equation}
	C = \exp^{-i\vec{\alpha}\cdot \vec{\lambda}}
	\end{equation}

	Thus we seek to learn the $d^2-1$ parameters $\vec{\alpha}$.  This general construction in $U(d)$ allows us to get at some of the results of \cite{VBB} without needing to entangle coins in lower dimensional spaces.  When we consider zig-zag products or even Cartesian Products we automatically increase the dimensionality of our coin Hilbert space.

	As one can imagine this generally gives us more success that training a single coin on the Maximal Torus.

	\subsection{A Single Set of Parameters per Step}

    The difference in the method from \S4.1 to \S4.2 is akin to the differences in \S3.2 and \S3.3, however, given the high number of parameters to tune, this method becomes computationally expensive after even a small graph.  For this reason, we show examples only from \S3.3 and omit training a general walk of $n$ steps required $n(d^2-1)$ parameters to train at every iteration.

	\section{Results and Examples}
	\subsection{Petersen Graph; a Single Coin}
	The first graph we will explore is the Petersen graph.  In this case we have a 10-vertex, 3 regular graph, which means our quantum walking matrix is $30\times 30$.  We were unable to fit a randomly drawn distribution with a single coin, nor were we able to with a single set of parameters.  Just to get a sense of this, consider the following graph:
	\[
	\includegraphics{Petersen_one_coin}
	\]
	This is a graph of the best fit we could get on the Petersen graph using a single coin and 20 steps.  Notice that we have certainly done better than having all the mass accumulated at vertex 1, but still not really teaching the graph the correct distribution.  We can see, however, that it has taken the parameter of 20 steps seriously.  Consider the following graph of errors per step.

	\[
	\includegraphics{Petersen_20_stepsone_coin}
	\]

	We see a clear downward spike at exactly 20.  This means the graph was training to learn how to produce the desired probability using exactly 20 steps.  Part of how we should like to approach quantum image processing is moving one step further or one step less using the ``correct" coinage.

	\subsection{Petersen; One Coin per Step}
	Now let's consider using a single coin per step.  In each case we used 8 steps with 3 parameters per step.  These were fit using a completely randomized initial set of coins, and run for 300 iterations or gradient descent.

	The first graph is a randomly drawn distribution.  The orange line is the probability distribution, and the blue dots are the learned distribution.
	\[
	\includegraphics{DTQW_Learns}
	\]

	The second graph is the most difficult.  We tried to get a quantum walk to converge to a single point. The difficulty is in getting all the wave function to destructively interfere everywhere except a single vertex, but constructively interfere at exactly one vertex.  After 300 iterations we arrived at the following.
	\[
	\includegraphics{vertex7}
	\]

	The third distribution we have chosen to train is the Boltzmann distribution.  This was the best fit of the three on the Petersen graph.

	\[
	\includegraphics{Boltzmann}
	\]

	\subsection{Fullerene $C_{60}$; One Coin per Step}

	Now let's give ourselves a big challenge.  We're going to try to have the fullerene graph $C_{60}$ learn a Gaussian distribution.  To make things even more complicated, let's set a few parameters.  The diameter of $C_{60}$is 9.  So we can't simply use 9 steps.  That's asking too much.  We'll use 13 steps.  There is very little room for error in 13 steps on $C_{60}$.  Additionally, the vertex labeling given in a somewhat random fashion rather than labeling vertices in a Hamiltonian path.  We use the 5-fold symmetric two dimensional projection and begin in the center and label vertices in a counterclockwise direction moving outward.  The Gaussian distribution is centered between vertices 29 and 30, in fact we're using
	\[
	\pi = A\exp(-(x-29.5)^2/(180))
	\]
	where $A$ is a normalizing constant.  That is, we're using $\sigma^2 = 3|V|/2$.  This gives a nice symmetric look in the graph and doesn't concentrate the probabilities too much in the center (making the distribution harder to learn).  So without further ado, let's take a look at some partial results.
	We start at vertex 1 with all the probability giving us this graph
	\[
	\includegraphics{c60_0iterations}
	\]
	We can check in on the progress at several iterations.
	\[
	\includegraphics{c60_120iterations}
	\]
	\[
	\includegraphics{c60_400iterations}
	\]

	The challenge we've given ourselves is immense.  This graph suffers from vanishing gradients quickly.  So after about 1000 iterations, we substituted a mildly stochastic gradient descent.

	\[
	\includegraphics{c60_940iterations}
	\]
	\[
	\includegraphics{c60_1500iterations}
	\]
	\[
	\includegraphics{c60_1800iterations}
	\]

	We now see the probabilities starting to settle onto the distribution.  Given the rate at which errors decrease, we expect this to settle in around 6000 or so iterations.

	\[
	\includegraphics{c60_errors}
	\]

	We can see the small bump in errors starting around 1000 iterations.  This was a modified stochastic gradient descent.

	\subsection{A Picture of Oak Ridge}
	Let's look at one more fun aspect of quantum walking.  Here's a picture of Oak Ridge, TN (not the lab, just some road in autumn)
	\[
	\includegraphics[width = 8cm]{oakridge}
	\]

	Now consider the following related image
	\[
	\includegraphics[width = 8cm]{oakridge_quantum}
	\]

	This may be the first image ever generated by quantum image processing via discrete time quantum walks on graphs.

	The process to create this second image is fairly simple.
	\begin{enumerate}
		\item Using an image library, separate the channels into red, green, blue.
		\item For this image we used only the green channel, and separated it into columns.  This image had channels of size $318\times 474$.  That's too large to run an entire quantum walk on a laptop, so we ran this column by column.
		\item For each column reduce the numerical input to a probability mass function $C \rightarrow C/\sum C$
		\item Initialize a discrete time quantum walk on the cyclic graph $C_{318}$, with initial position vector as the probability of the column
		\item Build a walking matrix of two steps with small, but random coin angles.
		\item Measure the probability of the 2 step quantum walk, then rescale.
		\item Replace the original column with the rounded, scaled quantum walk probabilities.
	\end{enumerate}

	Looking closely at the second image, the wave patterns appear vertically.  This reveals the quantum mechanical nature of the walk, rather than choosing randomly or simply shifting or running a few steps in a Markov chain.

	When the hardware is scaled up, the author would like to train a quantum walk on a large probability mass function and run a quantum walk with initial state as the image itself, and advanced it forward.

	\section{Applications}

	This work is at the intersection of several exciting topics as of early 2021.  Namely, quantum computing and machine learning.  These, however, are only the broadest brushes with which we can paint this method.  Learning a probability distribution, however, is the goal of several famous Monte Carlo methods, for example, Metropolis and Metropolis-Hastings, and Metropolis Coupled Markov Chain Monte Carlo (MCMCMC)\cite{mc3}, Metropolis-adjusted Langevin Algorithm \cite{RT}, etc.  Recently there has been some movement on learning distributions for time series forecasting\cite{BBO}, hurricane forecasting \cite{WT}, image processing, music generation, and natural language processing. \cite{WDL}

	In the actual implementation of this work, the goal is to teach a quantum walk how to control the angles of its interfering waves.  While this is not directly control theory, it is tangentially related.
	An additional area of application which is receiving some attention in the quantum world is ``post-quantum" cryptography.  In \cite{EHMLP} the authors develop cryptography schemes on super singular isogeny graphs.  These are regular graphs whose vertices represent super-singular elliptic curves over finite fields.  The edges represent morphisms between them.  This gives a family of graphs with $q$ vertices and $q(q-1)/4$ edges for $q\equiv1 \mod 4$.  The graph walking methods introduced here are able to handle many of these graphs on digital computers.  One should note, however, that as of 2021 there is no known quantum algorithm which can break this cryptography scheme.  Perhaps quantum walks on super-singular isogeny graphs could be of some use.

	\section{Open Questions and Future Work}

	As this work is relatively early in the intersection of DTQW and machine learning there is a bevy of open questions around it.  Among the most obvious to the author are, how can this be made more efficient.  One practicality in this work is that quantum walks on graphs can be honestly computed  on digital computers if one allows for tiny rounding errors.  However, while attempting to compute a distribution fitting scheme on the cross product of two graphs (a 14-vertex 3-regular with a 7-vertex 4-regular graph) the relatively small size of 98 vertices on a 7-regular graph took an inordinate amount of time.  In using a single coin per step in the maximal torus scheme, only 8 coins causes  one to train 56 parameters per iteration and use a product of fairly sparse matrices of size $686\times 686$.  This is memory intensive and very computationally expensive.  It is not fully clear the role that good spectral expander play in the DTQW universe is yet, however, having a low degree of regularity and small diameter reduces the number of parameters one needs to train, but increases the sparsity of the matrices involved.  Take the buckyball for example.  It's 60-vertex 3-regular structure means that the adjacency matrix has density 180/3600 = 0.05.  Exactly 95\% of the matrix is empty.  It's diameter is 9.  So we need at least 9 coins to even traverse the entire graph.

	A second open question that has come up several times in the course of this research is the existence and implementation of a polynomial time constructive algorithm for consistent roation maps for regular graphs.  \cite{A} shows the necessity of having such maps for ``efficient" construction of DTQW.

	A third question we'd like to explore is the roll of directed regular graphs, and their slight generalization of weighted graphs,for example \cite{wong}.  Weighted graphs have found incredibly many uses in machine learning an optimizations.  With the aid of \cite{S} one expects quantum speed ups of Markov Chain algorithms, but the question becomes what is the necessary hardware to implement an algorithm with any utility outside of pure academic interest.  In particular, the utility of a consistent rotation map on regular graphs reduces the computational space of DTQW from $|V|^2 \times |V|^2$ to $d|V|\times d|V|$ which reduces computational cost by density$^{-2}$.  How far can we take this.

	A fourth question we are considering and working in actively, is the successful implementation of real quantum walks on regular graphs.  This will improve the size of problems that we can tackle since some modern numerical packages (in particular numpy) has a terrible difficulty with complex numbers.  By reducing our coin space from $U(d)$ to $O(d)$ or even $SO(d)$ we cut our dimensions in half.  That is $mathbb{S}^{2d|V|-1} \rightarrow \mathbb{S}^{d|V|}$.

	Another area the author is actively exploring is hybridizing feed forward neural networks by quantum hidden layers.  The quantum hidden layers are multiple random quantum walkers on a graph in between classical Dense layers in a feed forward network.

	A final question and area of future work is the rigorous calculation of material properties using DTQW or CTQW on molecular graphs.

	\subsection{Proposed Additional Methods for Finding Angles}

	In order to increase the efficiency of solving the coins for a quantum walk we propose using several additional methods.  The first method we propose is for the maximal torus angle finding.  Instead of using a diagonal matrix in form of
	\[
	\begin{bmatrix}
	e^{i\theta_1} & & &\\
	& e^{i\theta_2} & &\\
	& & \ddots & \\
	& & & e^{i\theta_d}
	\end{bmatrix}
	\]

	which is $2\pi$ periodic to a matrix with unit periodicities

	\[
	\begin{bmatrix}
	e^{2\pi i\theta_1} & & &\\
	& e^{2\pi i\theta_2} & &\\
	& & \ddots & \\
	& & & e^{2\pi i\theta_d}
	\end{bmatrix}
	\]

	Then we can use either a classic Monte Carlo random sampling, or in lower dimensions, the quasi-Monte Carlo method of sampling by low discrepancy sequences.  In the case of cubic graphs, we have bounds for optimal low discrepancy sequences in $[0,1]^3$. \cite{WLDS}

	We need not use low discrepancy sequences to find the angles exactly, but there is essentially no computation aside from the guess and check method.  The lowest error derived by a low discrepancy sequence can be the initial angles in our coin, thus significantly reducing the number of iterations required in gradient descent to land on our desired probability distribution.

	A second proposed method is to use an evolutionary algorithm, such as a genetic algorithm, simulated annealing, or other MCMC algorithm modified for optimization.  From the memory perspective we will almost always prefer simulated annealing to a genetic algorithm, but for ease of implementation a genetic algorithm will win.

	The basic idea of implementation is to choose our angles at random, measure the error and then guess and measure again.  Whichever is better, we keep going in its direction.  We allow ourselves to step backward from time to time, but that is determined by a parameter (often called temperature).  To change our guess, we draw from a proposed probability distribution (often the normal or uniform distributions, although we could use a DTQW here to draw random samples).  The efficacy of this method comes when we draw small samples and add them to the current samples.  That is
	\[
	f(\vec{x} + \vec{\varepsilon}) \approx f(\vec{x}) + \nabla_x f(\vec{x})\cdot \vec{\varepsilon}
	\]
	which is, of course, a simulated gradient calculation.
	\section*{Appendix: Where to Find the Code}
	The author will periodically release modules of this code to GitHub, however, a better bet is to send a direct communication to the author and a small zip can be sent via email.  Languages are primarily Python and Julia, although Octave can still be found in some places where numerical tests were run.


\section{Introduction}

\IEEEPARstart{T}{he} Dzyaloshinskii-Moriya interaction (DMI) is a critical property of magnetic thin films for the stabilization of chiral magnetic objects, including skyrmions, for use in future spintronic applications.\cite{Dzyaloshinsky1958, Moriya1960,muhlbauer2009,Yu2010,Woo2016} For the interfacial DMI, chiral N\'eel-type domain walls (DWs) are preferred, with the sign of DMI determining the direction of the domain wall internal magnetization.\cite{Thiaville2012} This chirality is critical in governing the efficiency with which both domain walls and skyrmions can be manipulated by spin currents.\cite{Emori2013, ryu2013}. However, direct measurement of the strength of this interaction remains an experimental challenge. Commonly used techniques include Brillioun Light Scattering (BLS) and asymmetric domain expansion, which have offered agreement only in some materials systems.\cite{arora2020variation,Hrabec2014,Lau2016,Lavrijsen2015,Lau2018}

Here, we propose the use of a thin film patterned into the shape of a "nanohorn" that would enable the interfacial DMI to be read directly from the equilibrium position of a magnetic domain wall.  The concept is based on a balance between two energetic driving forces acting on the domain wall.  First, the nanohorn tapering drives the DW towards the narrow end due simply to the favorable reduction in total length.  Second, an in-plane field coupled with the interfacial DMI will drive the DW toward the opposite end of the nanohorn depending on the sign of DMI.  This is due to the Zeeman energy asssociated with the in-plane field and the internal magnetization of the DW itself.  The result is an intermediate position that corresponds to the strength of the interfacial DMI, which we elaborate on below.

\section{Magnetic Nanohorn Design and Modeling}

\begin{figure}[h!]
\includegraphics[width = 3.5 in]{Figure1_v3.pdf}
\caption{\label{Figure1} Schematic of an example magnetic nanohorn with $c/b = 0.1$.  Tick marks indicate the equilibrium position of a DW for different values of DMI.  Lower left inset: coordinate system used for equation 2.}
\end{figure}

The proposed nanohorn design is shown in Fig \ref{Figure1}, which is defined as the area between two curves given by the following polar equations:

\begin{equation}
r(\theta)=a\pm b\mp c(\frac{3\pi}{2}-\theta)
\end{equation}

for $\theta$ ranging from $0$ to $\frac{3\pi}{2}$ (as defined in Fig \ref{Figure1}). This produces two spirals centered at a circle of radius $a$. The wider end of the horn has length $2b$ and the radial cross section of the horn decreases with slope $2c$.

A horn with shape parameters $(a,b,c)=(400 nm,100 nm,10 nm/rad)$ was used to perform micromagnetic energy calculations on a 1nm thin film.  We use the steepest conjugate gradient method built into the software package Mumax$^3$.\cite{vansteenkiste2014design} The cell size was fixed at 1x1x1 nm with the following material properties $M_s=800 kA/m$, $A=15 pJ/m$, and $K_{u}=1\times10^7 J/m^3$. We apply an in-plane field $\mu_{0}H_x=150mT$, which as noted previously, creates a Zeeman energy with the DW internal magnetization. A domain wall was positioned at a range of angles from $\theta=0$ to $\frac{3\pi}{2}$ and at each one, magnetization was relaxed to minimize total energy. The equilibrium position was determined from the minimum in a plot of relaxed energy vs $\theta$ (e.g. Fig \ref{Figure2}). The resulting equilibrium positions as a function of $D_{int}$ are shown in Fig \ref{Figure3} for two possible horn dimensions including both micromagnetic and analytical results.  It is notable that there is a finite range of $D_{int}$ that results in an equilibrium position that does not fall to one of the limits, which is discussed later.

\begin{figure}[h!]
\includegraphics[width = 3.5 in]{Figure2_v2.pdf}
\caption{\label{Figure2} The relaxed energy vs $\theta$ in both a) micromagnetic simulation and b) semi-analytical computation shows a minimum which changes position based on DMI.}
\end{figure}

We now turn towards a semi-analytical treatment of the DW position in the nanohorn. It is well-established that the energy of a magnetic domain wall in the thin film limit can be approximated by the following:\cite{Pellegren2017, Thiaville2012}

\begin{equation}
\begin{split}
    \sigma(\Theta,\phi) = \sigma_0 -\pi D_{int}\cos(\phi-\Theta) - \pi\lambda\mu_0H_xM_s\cos(\phi)\\ +\frac{ln(2)}{\pi}t_f\mu_0M_s^2\cos^2(\phi-\Theta)
\end{split}
\end{equation}

where $\Theta$ and $\phi$ represent the angles of the DW normal and internal magnetization relative to the applied field $H_x$, respectively. $D_{int}$ is the interfacial DMI and $M_s$ is the saturation magnetization. $\sigma_0 = 4\sqrt{AK_{eff}}$ is the Bloch wall energy and $\lambda = \sqrt{A/K_{eff}}$ is the Bloch width, where $A$ is the exchange stiffness and $K_{eff} = K_{u} - \frac{\mu_0}{2}M_s^2$ is the effective anisotropy constant. The length of the domain wall as a function of $\Theta$ in the nanohorn can be derived:

\begin{equation}
\begin{split}
l(\Theta)=2b-2c(\Theta+\frac{\pi}{2})
\end{split}
\end{equation}

Here we ignore possible DW tilting, which is generally consistent with the results of the micromagnetic modeling. In this case, $\theta=\pi-\Theta$.

To compare with the Mumax3 simulations, we used the same horn shape and material parameters. Note that $b$ and $c$ carry the same scale factor and $\frac{c}{b}$ is the parameter which defines tapering. Given equations for DW energy per unit volume and length, (2) and (3), we consider their product the DW energy per unit length. At each particular value of $\Theta$, the DW energy per unit area was minimized with respect to DW internal magnetization $\phi$ to determine the relaxed energy as a function of $\Theta$. As with the micromagnetic calculations, the equilibrium position of the DW is determined from the minimum in energy vs $\Theta$. We repeat this for a range of DMI and compare with the date from Mumax3 simulation, which is superimposed in Fig \ref{Figure3}.

\begin{figure}[h!]
\includegraphics[width = 3.5 in]{Figure3_v4.pdf}
\caption{\label{Figure3} In both micromagnetic simulation (dots) and semi-analytical computation (curves) for multiple horn shapes, equilibrium $\theta$ has a clear relationship with DMI and the $c/b$ tapering parameter.}
\end{figure}

\section{Discussion}

In both the micromagnetic simulations and semi-analytical calculations we show a nonlinear relationship between DMI and equilibrium DW position. Low DMI below a particular threshold is too weak to result in an equilibrium position and instead the domain wall vanishes at the shorter end of the nanohorn. It is noteworthy that this critical value can be shifted to higher values of DMI by increasing the degree of tapering. Above that threshold DMI, there is a steep change in $\theta_{eq}$, which levels off and asymptotically approaches $\pi$. This is a function of the slope of the nanohorn and can be tuned to either reveal equilibria for lower DMI by decreasing the degree of tapering or be more sensitive at higher DMI by increasing the tapering.

For larger nanohorns, there is also the possibility of the domain wall no longer remaining radial with the center of the horn. In other words, the domain wall can tilt, introducing further complexity, which may be the origin of some of the disagreement between the micromagnetic and analytical results. In an experimental context, these nanohorns could be patterned optically and imaged by Kerr microscopy allowing for the simultaneous examination of a multitude of devices.  To account for uncertainty in the initial guess of DMI, it would be reasonable to create an array of devices with varying dimensions to improve the precision of the measured DMI.  We also acknowledge that there may be limitations in driving the wall to its equilibrium. This could be addressed through increased statistics (simultaneously examining many nanohorns of different dimensions) or exciting the DW with AC magnetic fields to overcome edge pinning.

\section{Conclusion}

We present computational and semi-analytical evidence for the use of magnetic nanohorns to measure interfacial DMI. The combination of a tapered shape and applied in-plane field creates an equilibrium position for a domain wall. We show that this equilibrium position varies with DMI, which may provide a technique for determining DMI of a thin film.

\section*{Acknowledgments}

The authors acknowledge valuable discussions with Michael Kitcher, Maxwell Li, and Nisrit Pandey. The authors acknowledge use of computational facilities of the Materials Characterization Facility at CMU under grant \# MCF-677785.

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\bibliographystyle{IEEEtran}