text:
\section{Introduction}

Simulation of strongly correlated electrons in the context of quantum chemistry and condensed matter physics is one of the potential areas in which quantum computers will have a significant advantage over their classical counterparts \cite{Cao_2019, Wecker_2015, Ma_2020}. Strong interaction between the ostensibly simple electrons can give rise to novel phases including high-temperature superconductivity \cite{highTc_1, highTc_2}, strange metallic behavior\cite{mcgreevy_2010}, fractional excitations \cite{Laughlin_1999}, and quantum spin liquids\cite{Savary_2016}. Condensed matter physics aims to understand these novel behaviors by studying their phase diagrams. However, failure of perturbation theory and exponential scaling of the Hilbert space for strongly correlated electrons presents a formidable challenge to classical simulation methods such as exact diagonalization, density matrix renormalization group(DMRG) \cite{dmrg1, dmrg2, dmrg3}, quantum Monte-Carlo \cite{qmc} and dynamical mean-field theory\cite{dmft1}.
Whereas, this same challenge provides an exciting opportunity for near-term quantum computers.

\comment{2. Opportunities and challenges for quantum computers. Need to change the organization of this paragraph?
Feynman\cite{Feynman1982} conjectured in 1982 that a quantum computer will be able to efficiently simulate quantum mechanical systems more efficiently than a classical computer. The power of a quantum computer lies in its control over the exponentially scaling Hilbert space of qubits, thus turning the curse of dimensionality into a blessing. }

Harnessing the power of a quantum computer to simulate quantum systems \cite{Feynman1982} requires (i) algorithms that can be executed in a reasonable time and (ii) the ability to learn from quantum experiments without exponentially many measurements. Studying the phases via ground state preparation is a QMA-complete problem \cite{Gorman_2021, Gharibian2019complexity, Bookatz2014QMA, kemp2004complexity}, which cannot be carried out in a reasonable time, even with quantum resources. However, performing dynamics on a quantum state is known to be a BQP-hard problem \cite{Baez2020dynamical, Rudi2020approximating}, possible within polynomial time. Likewise, it has been shown that shadow tomography\cite{Aaronson2018} methods such as classical shadows\cite{huang2020predicting, huang2021provably, Huang2022Learning} are effective at predicting properties using very few measurements.

In classical statistical mechanics, ergodicity provides a link between time averages and statistical averages \cite{Sethna_2012}. In quantum mechanics, unitary time evolution retains memory of the initial state, so the link is different yet still enabled by equilibration of macroscopic observables. This is made possible by certain general conditions on an initial state\cite{von2010proof,  Gogolin_2016,Linden_2009,Reimann_2012,Pintos_2017}, which allow these observables to equilibrate quickly.
Thus, if we can find interesting problems and initial states for which the order parameter (a macroscopic observable that characterizes the existence of a phase) equilibrates before the qubits decohere, we can use quantum computers to perform microcanonical dynamics and study the regions of the phase diagram inaccessible via classical computing resources.

In this manuscript, we present an approach to identifying phase diagrams of strongly correlated systems via dynamics simulations on a quantum computer. It consists of i) identifying an initial state, ii) evolving this state in time using a suitable algorithm, iii) using classical shadow tomography technique to convert the quantum state to classical data, and iv) applying unsupervised machine learning methods to discover phases of matter. We illustrate it with the 1-dimensional Transverse Field Ising Model(1DTFIM).

\begin{figure*}[t]
\includegraphics[scale=.14]{Schematic.pdf}
\caption{Schematic overview of our study. (a) Classical shadows(CS) of ground states  and time averaged classical shadows(TACS) from dynamics of a time-reversal invariant GHZ state are generated using quantum simulation. (b) The classical data from quantum simulation is then fed into diffusion map, an unsupervised machine learning algorithm to learn the phases.}
\label{Fig:Schematic_diagram}
\end{figure*}

A schematic overview of our approach is shown in Fig. \ref{Fig:Schematic_diagram}. It consists of using diffusion maps~\cite{coifman2006diffusion, de2008introduction}, an unsupervised machine learning(UL) algorithm,  to machine learn phase features from unlabeled phase diagram data generated by a quantum computer. Specifically, we apply diffusion maps to CS of ground states from a 100 qubit simulation and show they successfully learn the magnetic phase transition of 1DTFIM (see section II). It does so by \emph{discovering}  the relevant parameters such as magnetization and transverse magnetic field. These results follow similar machine-learning-from-CS success stories in other contexts\cite{ huang2020predicting,  huang2021provably, randomized_toolbox_2022, Huang2022Learning}. We then show (section III) diffusion maps can also identify the phase transition and quantum critical region \emph{when trained on dynamics data from just 20 qubits}. This data is in the form of a representation we introduce: time-averaged classical shadows (TACS)---a shadow tomographic\cite{Aaronson2018} representation of von Neumann's time-averaged density matrix\cite{von2010proof}. TACS are a straightforward extension of CS to capture quantum microcanonical dynamics of finite-size finite-time quantum simulations. Our results thus validate the exciting possibility of studying the phase diagrams of strongly correlated electrons using quantum computers.

\section{Von Neumann's Microcanonical Ensemble}

A central goal of quantum computing is to build qubits that are completely isolated from their environment. While this is not the case today, the current development of quantum error correction techniques\cite{Knill_2005,nielsen_chuang_2010} suggests it is in our future. Simulating quantum systems on a quantum computer will therefore take place \emph{within the microcanonical ensemble}. But quantum microcanonical dynamics, the evolution of a closed quantum system under Schr\"odinger's equation, does not directly produce the microcanonical ensemble.

Following von Neumann's 1929 paper\cite{von2010proof,goldstein2010normal} on the quantum ergodic theorem, it is straightforward to derive the equilibrium distribution of a quantum system using density matrices. Assuming we start from an initial state $|\psi(0)\rangle$ and evolve under a Hamiltonian $H$ via a quantum circuit algorithm to $|\psi(t)\rangle$, the equilibrium distribution is captured by the time-averaged density matrix
\begin{equation}\label{eq:rhoT}
    \rho_T\! =\! \frac{1}{T}\int_0^T\!\! dt |\psi(t)\rangle\langle\psi(t)| \!\underset{T\to\infty}{\longrightarrow}\! \omega\! =\! \sum_k P_k|\psi(0)\rangle\langle\psi(0)|P_k,
\end{equation}
where $P_k$ is a projector onto the $k^\text{th}$ degenerate subspace of the energy eigenvalues i.e. $P_k = \sum_{n\in k}|E_n\rangle\langle E_n|$.
\comment{
We can understand this expression from the snapshot picture of statistical mechanics. An ensemble of states/systems can be thought of as equivalently a history of one system where we record the state of the system in a series of snapshots over time. If the frequency of snapshots is sufficiently low, then each snapshot will be statistically independent of the previous one leading to a probability of $dt/T$ for each snapshot and the above time averaged expression for the density matrix of this ensemble. }
Existence of this equilibrium density matrix $\omega$ results in the ergodic principle that long-time averages of observables should  be captured by the statistical averaging with respect to $\omega$. Specifically, in the Schr\"odinger picture,
\begin{align}
  \langle O\rangle_T &=  \frac{1}{T}\int_0^T dt
    \langle\psi(t)|O|\psi(t)\rangle \\
    &= \frac{1}{T}\int_0^T dt \text{Tr}\left( |\psi(t)\rangle\langle\psi(t)|O\right) \\
    &= \text{Tr} \left(\rho_T O\right) \underset{T \to \infty}{\longrightarrow} \Tr(\omega O)
\end{align}
Thus, the time-averaged density matrix is the link between time averages and statistical averages.

The connection to the microcanonical ensemble, obtained by quantizing the classical microcanonical ensemble, is achieved by taking the thermodynamic limit, measuring only coarse-grained observables, and requiring non-degenerate energy level spacings. Under these circumstances, von Neumann obtained
\begin{equation}
    \omega\to\rho_{mc} = \frac{1}{\Omega}\sum_{E < E_n < E + \Delta} |E_n\rangle\langle E_n|
\end{equation}
where all states in an energy window $[E,E+\Delta]$ are included in the sum for some window size $\Delta$, and $\Omega$ is the phase-space volume of this constant energy shell. So $\rho_T$ approaches $\rho_{mc}$ in the limit of long times and large systems sizes, but for any finite $T$ and system size, $\rho_T$ is the actual microcanonical ensemble defined as the ensemble produced by the Schrödinger equation.

In numerics, two things are of utmost importance---reaching the thermodynamic limit and understanding the errors that may result along the way. Using the above logic, an equilibrium phase diagram with finite size errors associated with a reasonable choice of initial conditions is achievable. For each initial condition, it's possible to check whether the time series was run long enough for relevant observables to reach equilibrium (not all observables equilibrate until an exponential amount of time has passed\cite{Gogolin_2016, Short_2012}). If so, the only errors remaining are of finite size and will be different for each choice of initial conditions. A family of initial conditions that all yield the same microcanonical ensemble, when scaled up, is, therefore, a key resource to study equilibrium statistical mechanics on a quantum computer. Although the equilibrium properties are independent of the initial state, the time required to reach equilibration depends on the initial state via the effective dimension given by $d_\text{eff} = 1/\sum p_k^2$, where $p_k =|\bra{\psi} \ket{E_k}|^2$. Here, $p_k$ measures how many energy eigenstates have a significant weight in $|\psi\rangle$ and subsequently, a large effective dimension results in rapid equilibration with a bound on the equilibration time given by the second largest occupancy \cite{Gogolin_2016, Short_2012}. The distribution of the $p_k$'s also affects equilibration times\cite{margolus2021counting}. These arguments suggest we want to look for an initial state that exhibits a small overlap with most energy levels or a macroscopic occupation of a single energy level. In practice, we are concerned with the distinguishability of the time-averaged state vs. the equilibrium state in terms of experimentally relevant observables. Therefore, we need not worry about the equilibration of the full-density matrix as long as the equilibration of the relevant order parameters can be achieved.

In order to study von Neumann's microcanonical ensemble, we need to compare it to what we already know. As such, in the next section, section \ref{sec:groundstatedata}, we will turn to a ground state study using the methods we will employ for our thermodynamics study to validate these methods first, before applying them to the study of TACS.

\section{Ground state data}\label{sec:groundstatedata}
To verify our approach to phase classification, we first apply it to ground states of the ferromagnetic 1DTFIM defined by the Hamiltonian:
\begin{equation}
    H = -\sum_{\langle i, j \rangle} Z_i Z_j + h_x\sum_{i} X_i,
\end{equation}
where $\langle.\rangle$ denotes the nearest neighbors, $Z_i$ is the Pauli operator in the z-direction, and $h_x$ is a parameter proportional to the transverse magnetic field. At zero temperature, this model has a ferromagnetic phase for $\abs{h_x}<1$, and a paramagnetic phase for $\abs{h_x}>1$. The ground state study aims to detect the second-order phase transition at $h_x = 1$ for a $100$ site  1DTFIM using diffusion maps. The ground states were generated using density matrix renormalization group (DMRG) \cite{dmrg1, dmrg2, dmrg3} with ITensor package\cite{itensor}.
Training datasets were generated using two kinds of measurements on these ground states-  (i) computational basis measurements to obtain the Z dataset and (ii) measurements on a random Pauli basis to obtain the CS dataset.

\subsection{Classical Shadows}

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.9]{CS_diagramatic.pdf}
\caption{Diagrammatic description of classical shadows showing a linear relationship between $S_N[\rho]$ and estimator $\sigma_N[\rho]\approx\rho$. (a) The full-density matrix $\rho$ can be approximated by summing over reduced classical shadows with a coefficient that grows exponentially in the number of remaining qubits. (b) A classical shadow $S_N[\rho]$ is obtained by summing over $N$ measurement outcomes on random Pauli bases. For a given $N$, a reduced density matrix, that involves smaller coefficients in the expansion, can be approximated more accurately compared to the full density matrix.}
\label{Fig:CS_diagramatic}
\end{center}
\end{figure*}

\comment{\michael{In principle, the choice of shadow tomography impacts what observables are visible in the data. Similarly, in finite-time microcanonical dynamics simulations, coarse grained observables equilibrate while fine grained ones do not. So, to study quantum dynamics on a quantum computer, \sout{we should} \eric{one must} adjust the choice of shadow tomography to best observe the behavior of equilibrating observables.}}

Obtaining any useful information from a quantum computer requires performing measurements on a quantum state, which is destructive to the quantum information by nature. Since the dimension of the Hilbert space increases exponentially in the number of qubits, a naive strategy to learn the state requires an exponentially large number of copies. Aaronson\cite{Aaronson2018} introduced an alternative method using the notion of shadow tomography, an approximate classical description of the quantum state, in which $M$ properties of a quantum state can be estimated with error $\epsilon$ by only $O(\frac{M}{\epsilon^2})$ copies of the state. We can think of a shadow as an approximation of a quantum state $\rho$ by summing over measurement outcomes $x$, obtained by performing measurements on bases $b$ for a quantum state $x$ i.e.
\begin{equation}\label{eq:cs}
    S[\rho] = \sum_{b,x}P(b)
        \mathcal{P}_{b,x}
        \rho
        \mathcal{P}_{b,x},
\end{equation}
where $ \mathcal{P}_{b,x}$ is a projector onto the measurement outcome $x$ on basis $b$, and $P(b)$ is the probability of choosing $b$.

Based on this notion, Huang et al.\cite{Huang2019predicting, huang2021provably, Huang2022Learning} developed an algorithm called classical shadows and showed that it is highly successful at learning the properties of a many-body system. Two kinds of measurement protocols were proposed to construct classical shadows-  (i) random Clifford measurements on the entire Hilbert space; (ii) random single-qubit Pauli measurements. Protocol (ii) results in very shallow measurement circuits and thus is more suitable for the NISQ-era\cite{PreskillNISQ} hardware. \comment{With random Pauli measurements, they show that predicting linear functions with error $\epsilon$ requires at least $\text{log}(M) \text{max}_i \text{tr}(O_i^2)/\epsilon^2$ copies of the state.} After measuring each of the qubits in some random Pauli basis $X$, $Y$ or $Z$ with outcomes $\pm1$, the post-measurement wavefunction is given by the product state $\left|s^{(n)}\right\rangle=\bigotimes_{l=1}^{L}\left|s_{l}^{(n)}\right\rangle$. Here, $\left|s_{l}^{(n)}\right\rangle \in\{|0\rangle,|1\rangle,|+\rangle,|-\rangle,|\mathrm{i}+\rangle,|\mathrm{i}-\rangle\}$ is a Pauli basis state to which the $l^{th}$ qubit has collapsed. The classical shadow $S_N[\rho]$ is obtained by summing over $N$ such randomized measurement outcomes as follows (also see Fig. \ref{Fig:CS_diagramatic}(b))
\begin{align}
    S_N[\rho] &= \frac{1}{N} \sum_{n=1}^N \ket{s^{(n)}} \bra{s^{(n)}}\\
    &= \frac{1}{N} \sum_{n=1}^N \ket{s_1^{(n)}} \bra{s_1^{(n)}} \otimes \cdots \otimes  \ket{s_L^{(n)}}\bra{s_L^{(n)}}.
\end{align}

The underlying quantum state $\rho$ can be approximated by adding the reduced classical shadows (see Fig. \ref{Fig:CS_diagramatic}(a)).  This sum simplifies to the following expression from Ref. \cite{Huang2019predicting, huang2021provably, Huang2022Learning}
\begin{equation}\label{eq:csshots}
    \rho \approx \sigma_{N}(\rho)=\frac{1}{N} \sum_{n=1}^{N} \sigma_{1}^{(n)} \otimes \cdots \otimes \sigma_{L}^{(n)},
\end{equation}
where
\begin{equation}
\sigma_{l}^{(n)}=3\left|s_{l}^{(n)}\right\rangle\left\langle s_{l}^{(n)}\right|-\mathbb{I}.
\end{equation}

The definition of $S_N(\rho)$ presented above is different from Refs. \onlinecite{Huang2019predicting, huang2021provably} which defines it to be\comment{;} the dataset of shots itself and not the density matrix obtained from these shots. But both are useful for Fig. \ref{Fig:CS_diagramatic} shows the density matrix $S_N(\rho)$ defined above is linearly related to the estimator $\sigma_N(\rho)$ of the quantum state $\rho$ obtained by Refs. \onlinecite{Huang2019predicting, huang2021provably, Huang2022Learning}. A derivation of Fig. \ref{Fig:CS_diagramatic} using tensor network diagrams is presented in Appendix A.

Although estimating the exact density matrix requires $N \rightarrow \infty$, we still desire to predict various linear as well as nonlinear functions of $\rho$ (e.g., $\Tr(O\rho)$ and $\Tr (\rho \text{log}(\rho))$ respectively). This can be achieved with $N \propto \log(L) 4^k/\epsilon^2$ copies of the state,  where $k$ is the locality of operator $O$ \cite{randomized_toolbox_2022}. It was shown in Ref. \cite{huang2021provably} that classical machine learning algorithms can efficiently predict the ground state properties of gapped Hamiltonians in finite spatial dimensions after learning the classical shadows from a training set. An example of interest is classifying the quantum phases of matter. Classifying the symmetry-breaking phases is conceptually simple because it involves calculating $\text{tr}(\rho O)$ for some k-local observable $O$, such that $\text{tr}(\rho O) \geq 1 \: \forall \: \rho \in \text{ phase A}$ and  $\text{tr}(\rho O) \leq -1 \: \forall \: \rho \in \text{ phase B}$.

In contrast to the symmetry-breaking phases, capturing continuous phase transitions involves nonlinear observables like thermodynamic potentials. These observables involve measuring a nonlinear function of $\rho$, for example, the entropy. Learning such nonlinear functions requires an expressive ML model. Thus, it is beneficial to devise a mapping from classical shadows to a high-dimensional feature space that includes the polynomial expansion of many-body reduced density matrices and encapsulate this feature space in a kernel function. Using such kernels then establishes a rigorous guarantee that a classical ML model can efficiently classify nonlinear functions, phase transitions, and even topological phases of matter\cite{huang2021provably}.

\subsection{Machine Learning Method: Diffusion Maps}

\comment{Classical shadows provide an approximate information about the full state of a quantum computer where the accuracy of the approximation depends just on the number of shots[ref]. But rather than work directly with this approximation, such as by computing observables that are accurately captured by the data, we will apply machine learning (ML) methods whose success depends on this approximation. As we will see, these ML methods can still identify macroscopic parameters of interest, such as the z-magnetization and $h_{x}$ values of the quantum states, as well as identify phases and phase transitions, even though the data we feed it (the classical shadows) contain far less information than the full description of the quantum states.}

In this paper, we utilized unsupervised machine learning methods called diffusion maps \cite{coifman2006diffusion, de2008introduction} to extract features from the shadows data. A diffusion map is a nonlinear dimensionality reduction technique that relies on learning the underlying manifold from which the data points have been generated. Recently, this method was used to identify phases in systems with complex order parameters, which are difficult to learn using linear methods (such as principal component analysis (PCA)\cite{carleo2019machine}). Examples of such phase identification studies include; topological phases and phase transitions \cite{rodriguez2019identifying}, incommensurate phases, and many-body localized phases in quantum systems \cite{lidiak2020unsupervised}.

In this application of diffusion maps, one can imagine a random walk on a dataset  $X\left(x_{1}, x_{2},\ldots, x_{N}\right)$, where the $x_{n}$ are estimators $\sigma_N(\rho)$ of density matrices $\rho$ obtained from different points in the phase diagram. The transition probability $P(j\backslash i)$ of jumping from $x_{i}$ to $x_{j}$ in a single ``timestep'' is proportional to the kernel function $k(x_i,x_j)$, which is a similarity measure between the two data points. Here we use the classical shadow kernel function prescribed in~\cite{huang2021provably}, defined to be for two points $x$ and $\tilde x$ in the dataset
\begin{equation}\label{eq2}
  \begin{aligned}[b]
    & k^{(shadow)}(x,\tilde x) = k^{(shadow)}(\sigma_N(\rho),\sigma_N(\tilde\rho)) =  {}\\
    & \exp{\!\left(\!\sum_{n,n'=1}^N\!\frac{\tau}{N^2}\!     \exp{\!\left(\!\frac{\gamma}{L}\!\sum_{l=1}^L\! \Tr\Big(\sigma_{l}^{(n)} \tilde\sigma_{l}^{(n')}\!\Big)\!\right)}\! \right)}.
  \end{aligned}
\end{equation}
Then, we construct a transition probability matrix $P$ such that
\begin{equation}
\label{eqn:P-matrix}
   {P(j\backslash i
   }
   = \frac{k^{(shadow)}(i,j)}{\sum_l k^{(shadow)}(i,l)}.
\end{equation}

After $t$ timesteps of the random walk, the transition probabilities are given by the matrix $P^t$, where $P_{ij}^t$ gives the probability of going from $x_{i}$ to $x_{j}$ in $t$ timesteps, it's a sum of the probabilities associated with all of the possible paths to go from $x_{i}$ to $x_{j}$ in $t$ timesteps. As $t$ increases, the diffusion process unfolds, where data points situated along the overall geometric structure of the dataset become more strongly connected because of the abundance of strongly connected intermediate points along the way. The parameter 'diffusion distance' is defined to quantify this idea of connectivity between two data points given by

\begin{equation}
    D_{t}^2(x_{i},x_{j})^2 = \sum_{m=1}^N|P_{im}^t-P_{mj}^t|^2
\end{equation}

where the bigger the diffusion distance, the weaker the connection between them. The data points are then mapped onto a new `diffusion space' so that the diffusion distance in data space is equal to the Euclidean distance in this new space. It can be shown \cite{coifman2006diffusion} that such a map can be made in the following way:
\begin{equation}\label{eq5}
    x_{i}\rightarrow y_{i} = \begin{bmatrix}
       \lambda_{1}^t\psi_{1}(i),  &
       \lambda_{2}^t\psi_{2}(i), &
       \cdots, &
       \lambda_{N-1}^t\psi_{N-1}(i)
    \end{bmatrix}
\end{equation}
where $\lambda_{k}$ and $\psi_{k}$ are eigenvalues and right eigenvectors of the matrix $P$, $\psi_{k}(i)$ is the $i$-the element of the $k$-th eigenvector. $\psi_{k}$ form an orthogonal set of basis vectors for the diffusion space.
Then the diffusion distance is,
\begin{equation}\label{eq6}
    D_{t}^2(x_{i},x_{j}) = |y_{i} - y_{j}|^2 \\
        = \sum_{k=1}^{N-1} \lambda_{k}^{2t}\left[\psi_{k}(i) - \psi_{k}(j)\right]^2
\end{equation}
In equations (\ref{eq5}) and (\ref{eq6}) the $k=0$ component is ignored ($\psi_{0}(i) = \frac{1}{\sqrt{N}},\; \lambda_{0}=1$), as it is constant for all $i$, hence trivial. Equation (\ref{eq6}) suggests a dimensionality reduction, as the terms with bigger $\lambda_{k}$ will dominate the sum increasingly as $t\rightarrow\infty$. So  only the $m$ diffusion space dimensions associated with the bigger eigenvalues ($m << N$) are kept. The data points are then projected onto this m-dimensional space to reveal features in the dataset.

\subsection{Phase Classification of Ground States}
For the ground state study, we used two datasets generated via two different shadow tomographies. The first (Z dataset) contains qubit measurements only along the Z-axis [in Eq. \ref{eq:cs}, $P(b)=1$ for $b=Z$, $P(b)=0$ for all other b]. While the other (CS-dataset) has randomized Pauli measurements using the Classical Shadows method [$P(b)=\frac{1}{3}$ for $b \in \{X, Y, Z\}$]. Both of them contain 200 100-spin 1DTFIM ground state shadows, each state for a different $h_{x}$ value ($h_{x}$ ranging from 0.1 to 100).  Since the $Z$-magnetization is the order parameter for TFIM, the UL algorithm should be able to learn the phases of the model from the Z dataset. Using this knowledge, we compare the Z dataset and the CS data to see if the algorithm can successfully identify phases in it.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth, height=0.3\textheight]{gr7.pdf}
\caption{Learning phases from ground state (CS and Z) data. (a) and (c) show the 10 largest eigenvalues of the $P$-matrix (excluding the trivial $k=0$) for (a) Z dataset and (c) CS dataset. (b) Z datapoints in 5D diffusion space visualized in 2D, using metric MDS. Clustering clearly emerges based on the two phases of the model. (d) CS datapoints in 2D diffusion space, the figure reveals the symmetry-breaking phase transition. In this case, there is a direct correlation of relevant parameters, the z-magnetization and the $h_{x}$-values, with machine-learned diffusion, coordinates dc1 and dc2 respectively. Notice the CS-data shows three clusters in (d), while Z-data shows only two in (b), which is a consequence of the shadow kernel function for the former is not spin-flip symmetric while that of the latter is.}
\label{ground}
\end{figure}

By deploying diffusion maps armed with the shadows kernel function (Eq. \ref{eq2}) as our UL model, we are able to identify the phases from both the Z data and the CS data. In both cases, we set $\tau=1$, $\gamma=1$ in Eq. \ref{eq2}. For the Z data, we chose the first five non-trivial eigenvectors as the diffusion space basis vectors because the $P$-matrix (Eq. \ref{eqn:P-matrix}) eigenvalue spectrum shows the first five eigenvalues to be larger than others (Fig. \ref{ground}(a)). Mapping the states from Z data onto this five-dimensional diffusion space, we found that clear clustering emerges based on the phases of the states. We used multidimensional scaling (MDS), a dimensionality reduction method \cite{borg2005modern} that seeks to preserve point-to-point distances, to project these states onto a 2D plane. We see a clear separation of the two phases even on this 2D reduced space (Fig. \ref{ground}(b)), indicating the machine's success in identifying the two phases.

From the CS data, the unsupervised learning algorithm was also able to learn about the phases and the underlying parameters of the model. We can see in Fig. \ref{ground}(c) that the $P$ matrix has two non-trivial eigenvalues larger than the rest. The eigenvectors corresponding to these two eigenvalues are the basis vectors of the reduced diffusion space. Figure (\ref{ground}(d)) shows all of the ground state classical shadows projected onto this two-dimensional plane. It shows three groups- the top left and the bottom left are the all-up and all-down states, whereas the group on the center-right are states in the disordered phase. This closely resembles the spontaneous symmetry-breaking phase transition of 1DTFIM\cite{sachdev_2011}. The learned diffusion coordinates (dc1 and dc2) have direct correlations with the magnetization $<M_{z}>$ and the field values $h_{x}$ respectively, as shown in appendix C (Fig. \ref{fig:correlation}).

The above-mentioned clustering is dependent on the number of snapshots $N$ for a given state. However, the dimensionality reduction and the subsequent clustering will settle down after a minimum value of $N$ has been reached ($N_{c}$). Increasing the value of $N$ beyond that point does not change the results in any significant way, and in our case ($N_{c}$ (Z data)$< N_{c}$ (CS data)), it should be easier for the algorithm to learn the phase space structure from the Z-data (hence fewer snapshots are required).

Despite the successful identification of the ground state phases with the CS+ML model, the fact remains that the problem of calculating ground-states is a QMA-hard problem \cite{Bookatz2014QMA}. While there is no shortage of algorithms designed to generate ground-states, the necessity of a scalable model can be better shown with a known BQP problem; revealing the microcanonical phase diagram from dynamics data.

\section{Microcanonical Dynamics of the 1DTFIM}
Consider now Fig. \ref{fig:phasediagram}, the thermodynamic phase diagram of the 1DTFIM relating internal energy $E=\langle\psi(0)|\hat H|\psi(0)\rangle$ to the transverse magnetic field $h_x$. This phase diagram is relevant for a microcanonical dynamics study governed by the entropy $S(E,h_x)$. This phase diagram exists even for a simulation over a finite time $T$ and with a finite number of spins $N$ but with finite $T$, finite $N$ errors that round the phase transitions and render the results dependent on initial conditions. We present in this figure our expectations for the phase diagram in this context, pointing out where the phase diagram will sharpen in the thermodynamic limit. We further highlight the path through the phase diagram carried out by our simulations below, showing that we expect it to cross the quantum critical region and so be sensitive to the phase diagram at a rounded level.

\begin{figure}
    \centering
    \includegraphics{PhaseDiagramTikz.pdf}
    \caption{A sketch of the expected 1DTFIM phase diagram at finite $T$ and finite $N$ as a function of internal energy $E$ and transverse magnetic field $h_x$. This diagram is a modification of the canonical ensemble representation of the phase diagram in Ref. \onlinecite{wu2018crossovers}, adapted to the microcanonical ensemble. The spectrum is mirror symmetric about $E=0$ due to the chiral symmetry $\mathcal{C}=ZYZYZY...$ \comment{that anticommutes with the Hamiltonian.}}
    \label{fig:phasediagram}
\end{figure}

\comment{\gaurav{\emph{Where is a good place to add this?} A similar algorithm called thermal shadows has been recently proposed in literature\cite{Predicting_Gibbs}. However, this requires a large number of controlled unitaries and optimization for the quantum signal processing\cite{QSP} subroutine. There is also quantum imaginary time evolution(qite) \cite{QITE}, which approximates a k-local imaginary time evolution operator (non-unitary) by C-local unitary where $C>k$. This involves approximating all the non-unitary operators corresponding to each term in the Hamiltonian by using a larger unitary, whose size grows with the correlations(C)  in the system. Our algorithm seems very straightforward compared to this, so will potentially replace qite. There is also another algorithm based on the zeros of partition function \cite{Partition_zeros}.}}

\subsection{Results}
With the knowledge of expected microcanonical dynamics (see Fig. \ref{fig:phasediagram}), we now present an effective classical representation of quantum data obtained from a microcanonical dynamics on a quantum computer. Noticing that the time averaged integral amounts to an expectation value of the pure state density matrix $|\psi(t)\rangle\langle\psi(t)|$ over the probability distribution $P_T(t) = (1/T)\left(\Theta(T-t\rangle)-\Theta(-t)\right)$, where $\Theta(x)$ is the Heaviside step function, we define \emph{time averaged classical shadows}(TACS) by the map
\begin{equation}\label{eq:tacs}
    TACS\left[\rho \right] = \lim_{T\to\infty}\int dt \sum_{b,\sigma}P_T(t)P(b)
        \mathcal{P}_{b,\sigma}
        |\psi(t)\rangle\langle\psi(t)|
        \mathcal{P}_{b,\sigma}
\end{equation}
Hence by sampling the joint probability distribution $P_T(t)P(b)$ to obtain $(t_i,b_i)$, $i=1\ldots N$, and then measuring one shot $\sigma_i$ from $|\psi(t_i)\rangle$ in basis $b_i$ we obtain a finite-shot TACS via
\begin{equation}\label{eq:tacsshots}
    TACS_N\left[\rho \right] = \comment{\lim_{T\to\infty}}
    \comment{\lim_{N\to\infty}
    \frac{1}{N}}
    \sum_{i=1}^{N}
    |t_ib_i\sigma_i\rangle\langle t_ib_i\sigma_i|
\end{equation}
We ran quantum dynamics of 1DTFIM using the TDVP algorithm\cite{tdvp,timeevomps} starting from the GHZ state $\ket{\psi(0)} = \frac{\ket{000\cdots} + \ket{111\cdots}}{\sqrt{2}}$to generate TACS data from 500 randomly sampled dimensionless time values between $t=10.0$ to $t=20.0$ and 187 randomly sampled $h_x$ field values between $h_x = 0.1$ to $h_x = 10.0$. For each $h_{x}$, a TACS was obtained from the $N = 500$ shots at random times over the fixed $T$ time interval following Eq. \ref{eq:tacsshots}. Then these 187 classical shadows were the data points with which we performed unsupervised learning by constructing the $187\times187$ kernel matrix using the shadow kernel in equation (\ref{eq2}) and then using diffusion maps for dimensionality reduction. Fig. \ref{fig:dynamics} shows our results. The first two eigenvectors were chosen as our diffusion space basis vectors because as Fig. \ref{fig:dynamics}(inset) shows, those are the two dominant non-trivial eigenvalues in the spectrum and the rest go to zero as the number of data points are increased. Projecting the states onto this two-dimensional diffusion space, we see that the states all fall on a curve in this hyperplane along which the value of $h_{x}$ increases monotonically, and the inflection point of the curve coincides with the quantum critical region. Therefore, it is apparent that the unsupervised learning algorithm was able to infer two phases from the data. In contrast to our ground state study, we were not able to map the diffusion coordinates to the phase-defining observables.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth,height=0.25\textheight]{dnmks.pdf}
\caption{The figure above shows phase identification from our dynamics data. The eigenvalue spectrum (inset) obtained from diffusion maps shows the two largest eigenvalues corresponding to the two dominant diffusion coordinates dc1 and dc2 (the trivial point $k=0$, is not shown). The TACS data points largely fall on a curve parameterized by $h_{x}$ in this 2D reduced diffusion space. The region near the critical point ($h_{x}=1$) coincides with the inflection point region (green) of the curve, with points on either side belonging to the two different phases.}
\label{fig:dynamics}
\end{figure}

In a microcanonical ensemble, entropy functions as the thermodynamic potential and thus characterizes the phase transitions. The quantum critical region features the interplay of thermal and quantum fluctuations leading to the entropy being maximized \cite{ Wu_2011_entropyaccumulation, wu2018crossovers}.  To calculate the $2^{nd}$ Renyi entropy, we first compute the purity using eqn \ref{eq:csshots} as follows
\begin{align}
    \gamma\left[\rho \right] &= \text{Tr}\left[ \rho^2 \right] \\
    &\approx \frac{1}{N^2} \sum_{n\neq n' = 1}^{N} \text{Tr}\left[ \sigma_1^{(n)} \sigma_1^{(n')}\right]\times \cdots \times \text{Tr}\left[ \sigma_L^{(n)} \sigma_L^{(n')}\right],
\end{align}
then $S_2= -\log_2 \gamma$.The Renyi entropy has the same limits as the von Neumann entropy- it vanishes for pure states and reaches $N$ for the maximally mixed state.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth,height=0.25\textheight]{Renyi_10_sites.pdf}
\caption{The $2^{nd}$ Renyi entropies for a 10-site 1DTFIM computed from TACS consisting of 100k shots.$n=\infty$ shows the Bayesian inference estimation for the thermodynamic limit and infinite TACS shots. }
\label{fig:renyi}
\end{figure}

Figure \ref{fig:renyi} presents the entropy as calculated from TACS shots on a 10 qubit system using exact diagonalization. We use Bayesian inference, as detailed in appendix \ref{app:bi} (see also Ref. \onlinecite{lukens2021bayesian}). The results show clear evidence that entropy is maximized in the quantum critical region around $h_x=1$ as shown in Fig. 2 of \cite{wu2018crossovers}. Error estimates for these values were obtained in Ref. \onlinecite{Shack}, and given by:
\begin{equation}
    N \geq \frac{4^{k+1}\gamma}{\epsilon^2\delta}
\end{equation}
Where $1-\delta$ is the probability of obtaining a good TACS dataset and $\epsilon$ is the additive error. For a $\delta = 0.33$ (67 percentile), and $N=100,000$ shots, we find an additive error for the $n=5$ entropy curve plotted in Fig. \ref{fig:renyi} at $h_x=1.0$ of $S_2/n \approx 0.5 \pm 0.24$. This is larger than the observed error shown via error bars in Fig. \ref{fig:renyi} but only within a factor of order 1.
\comment{need to say something about n=1 line}

\section{Outlook}
In this paper, we have identified an approach to studying quantum thermodynamics on a quantum computer in a way that is suitable for studying quantum materials and their phase diagrams. This approach consists of
\begin{enumerate}[label=(\roman*)]
    \item preparing an initial state which equilibrates within the coherence time of the quantum computer,
    \item time evolving this state using a quantum algorithm to a randomly chosen time point $t$ within some interval,
    \item extending shadow tomography methods to obtain a physically useful representation of the time-averaged density matrix, such as TACS used in this paper, and
    \item employing an unsupervised machine learning method to discover the phase diagram, with kernel methods such as diffusion maps showing promise.
\end{enumerate}
Our approach parallels statistical mechanics calculations on classical Hamiltonians, where a random initial state is prepared, a Metropolis Monte-Carlo algorithm is run beginning from this state, and data is collected and analyzed using traditional observables and modern machine learning methods \comment{\gaurav{refs? Is this statement important?}}. Our results, demonstrating the existence of a quantum phase transition and the ability to map out regions of the phase diagram by a careful choice of initial conditions, show promise.  Considering that we are just performing time evolution and measurements, our approach is simple compared to the existing methods \cite{Predicting_Gibbs, QITE, Partition_zeros} for studying thermodynamics on a quantum computer. The simplicity not only avoids the issues of optimizing a large number of parameters experienced by these methods but also allows us to exploit the advances in Hamiltonian dynamics. This leads to the possibility of studying thermodynamics on NISQ devices by using variational time evolution algorithms\cite{Yao_2021, Yuan_2019, VFF_2020} such as the variational fast forwarding. These algorithms are relatively robust to errors and can simulate dynamics beyond the coherence time of quantum computers. Since TACS involve only single qubit rotations and computational basis measurements, an extension of our methods to quantum simulators also seems straightforward\cite{Rydberg_atoms}. Lastly, it would be interesting to assess how this algorithm compares with finite temperature Quantum Monte-Carlo and at what system size the expected quantum advantage is reached.

\comment{
There are still a number of improvements that could be made to the approach. We have not scaled both time $T$, and system size $L$ and attempt to extrapolate to the thermodynamic limit where, according to von Neumann's quantum ergodic theorm, Eq. \ref{eq:rhoT} converges to the microcanonical ensemble.  Alternatively, we could phrase this improvement in two steps: 1) the problem of how to obtain the best estimate for $\hat\omega\equiv\lim_{T\to\infty}\hat \rho_T$, given a database of $\rho_T$ at different $T$; 2) the problem of extrapolating to infinite system size.

Another improvement would be to identify a dynamics algorithm that includes error correction to compute $\rho_T$ on a NISQ computer. One approach to this is to use variational time evolution algorithms \cite{Yao_2021, Yuan_2019, VFF_2020} such as the variational fast forwarding. These algorithms are relatively robust to errors, and are able to simulate dynamics beyond the coherence time of quantum computers. It has been shown that the error in variational fast forwarding algorithm scales at worst linearly in the fast forwarding time\cite{VFF_2020}. Post selection may further increase the accuracy of these variational algorithm by rejecting the states which violate the conservation laws. \gaurav{error corrected dynamics algorithms?}

\gaurav{Worth spectulating the use of our results for floquet dynamics. Instead of continuously sampling time, sample within each time interval. Read one of their papers and flesh out the conclusions. Express the time evolution operator in one basis $1/N U^N$ and try to get a similar formula for the time averaged density matrix. If that reduces to something like $\omega$, victory!}

\eric{Would our results be useful for studying quantum chaos? Perhaps only when we have error correction so that errors don't obscure the results.}\michael{Cool idea! Our central premise of simulating microcanonical dynamics presupposes we have error correction and the qubits are isolated from their environment.

Yes, one can study quantum chaos on a quantum computer by just time evolving a quantum system which is a quantized version of a classical system such as the kicked quantum rotor. Then how does equilibrium, studied by $\rho_T$ relate to the chaotic dynamics? Is equilibrium actually the result of quantum chaos? What about systems that don't equilibrate? To write this paragraph in the outlook, we need a reference that discussed equilibrium and quantum chaos and then add our two-cents based on what we know of the subject.}
\gaurav{Thermalization in quantum systems is a direct consequence of non-integrability (and quantum chaos in most cases). Thus, in some sense, thermalization can detect presence/absence of quantum chaos which arises from scrambling. A central tool in study of quantum scrambling is OTOC, which is a complicated quantity to calculate and our TACS does not necessarily lead to OTOC. How else does our method support the chaos community?}}

\section{Example}
	\begin{tikzpicture}[inner sep=1mm]
	\foreach \i in {1,...,5} {
		\node[tensor] (\i) at (\i, 0) {$A_\i$};
		\node (\i spin) at (\i, -0.7) {$s_\i$};
		\draw[-,thick] (\i) -- (\i spin);
	};
	\foreach \i in {1,...,4} {
		\pgfmathtruncatemacro{\iplusone}{\i + 1};
		\draw[-,thick] (\i) -- (\iplusone);
	};

	\end{tikzpicture}

	\section{}
	\begin{align}
	\frac{1}{3} \sum_{b,\sigma}
	\diagram{
		\node[tensor] at (1,0) (A) {$P_{b \sigma}$};
		\node[tensor] at (2,0) (B) {$\rho$};
		\node[tensor] at (3,0) (C) {$P_{b \sigma}$};
		\draw[-,thick] (A)--(B);
		\draw[-,thick] (B)--(C);
		\draw[-,thick] (0,0)--node[above] {$x$} ++(A) ;
		\draw[-,thick] (C)--node[above] {$x$} ++(1,0);}
	\end{align}

	\section{}
	\begin{align}
	\diagram{
		\def1/4{1/4}
		\def1/2{1/2}
		\draw[thick](-1/2,0) -- (1/2,0);
		\draw[thick](-1/4,3/4) -- (1/4,3/4);
		\draw[thick](-1/2,0) .. controls (-1/2-1/4,0) and (-1/2-2*1/4,3/4) .. (-1/4,3/4);
		\draw[thick](1/2,0) .. controls (1/2+1/4,0) and (1/2+2*1/4,3/4) .. (1/4,3/4);
		\node[tensor] at (0,0) (A) {$\rho$};} = 1
	\end{align}

	\section{}
	\begin{align}
	\frac{1}{3}\sum
	\diagram{
		\node[tensor] at (1,0) (A) {$P_{b \sigma}$};
		\node[tensor] at (1,1) (B) {$P^T_{b \sigma}$};
		\draw[-,thick] (0,0)--++(A) ;
		\draw[-,thick] (0,1)--++(B) ;
		\draw[-,thick] (A)--++(1,0);
		\draw[-,thick] (B)--++(1,0);
	} =
	\frac{1}{3}
	\left[ \;
	\diagram{
		\draw[-,thick] (0,0)--(0,1);
		\draw[-,thick] (0.5,0)--(0.5,1);}
	+
	\diagram{
		\draw[-,thick] (0,0)--(1,0);
		\draw[-,thick] (0,0.5)--(1,0.5);}
		\;
	\right]
	\end{align}

	\section{}
	\begin{align}
	\frac{1}{3}
	\left[
	\diagram{
		\node[tensor] at (1,0) (A) {$\rho$};
		\draw[-,thick] (0,0)--(A);
		\draw[-,thick] (A)--(2,0);
	}
	+
	\diagram{\draw[-,thick] (0,0)--(0,0);\draw[-,thick] (0,0.20)--(1,0.20);}
	\right]
	\end{align}

	\section{}
	\begin{align}
	\diagram{
		\node[tensor] at (1,0) (A) {$\rho$};
		\draw[-,thick] (0,0)--(A);
		\draw[-,thick] (A)--(2,0);
	}
	=
	3 \diagram{
		\node[tensor] at (1,0) (A) {$CS(\rho)$};
		\draw[-,thick] (0,0)--(A);
		\draw[-,thick] (A)--(2,0);
	} +
	\diagram{\draw[-,thick] (0,0)--(0,0);\draw[-,thick] (0,0.20)--(1,0.20);}

	\end{align}
\section{}
	\begin{align}
	\diagram{
		\def1/2{3}
		\def1/4{1/4}
		\node[tensor] at (1,0) (A) {$P_{b \sigma}$};
		\node[tensor] at (1,1) (B) {$P_{b \sigma}$};
		\node[tensor] at (2.5,0) (C) {$\rho$};
		\draw[-,thick] (0,2)--(4,2);
		\draw[-,thick] (0,0)--++(A) ;
		\draw[-,thick] (0,1)--++(B) ;
		\draw[-,thick] (A)--(C);
		\draw[-,thick] (B)--++(2.5,0);
		\draw[-,thick] (C)--++(1,0);
		\draw[thick] (3.5,0) .. controls (3.5+1/4,0) and (3.5+2/4,3/4) .. (3.5,1);
		\draw[thick] (0,1) .. controls (0-1/4,1) and (0-2/4,1+3/4) .. (0,2);
		\draw[thick][red,thick,dashed] (0.1,-0.5) rectangle (2,1.5);
	}
	\end{align}
\section{}
\begin{align}
\frac{1}{9}\sum _{\sigma, b} \sum_{\sigma', b'}
\diagram{
	\node[tensor] at (1,0) (A) {$P_{b' \sigma'}$};
	\node[tensor] at (1,1) (B) {$P_{b\sigma}$};
	\draw[-,thick] (0,0)--++(A) ;
	\draw[-,thick] (0,1)--++(B) ;
	\draw[-,thick] (A)--++(1,0);
	\draw[-,thick] (B)--++(1,0);
	\draw[fill=tensorblue, thick] (2,-0.5) rectangle (3,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (3,0)--++(1,0);
	\draw[-,thick]  (3,1)--++(1,0);
	\node[tensor] at (4,0) (C) {$P_{b' \sigma'}$};
	\node[tensor] at (4,1) (D) {$P_{b \sigma}$};
	\draw[-,thick] (C)--++(1,0);
	\draw[-,thick] (D)--++(1,0);
}
\end{align}

\section{}
\begin{align}
\diagram{
	\node[tensor] at (1,0) (A) {$P_{b' \sigma'}$};
	\node[tensor] at (1,1) (B) {$P_{b\sigma}$};
	\node[tensor] at (1,-1) (C) {$P_{b' \sigma'}$};
	\node[tensor] at (1,2) (D) {$P_{b\sigma}$};
	\draw[-,thick] (0,3)--++(4,0);
	\draw[-,thick] (0,-2)--++(4,0);
	\draw[-,thick] (0,0)--++(A);
	\draw[-,thick] (0,1)--++(B);
	\draw[-,thick] (0,-1)--++(C);
	\draw[-,thick] (0,2)--++(D);
	\draw[-,thick] (A)--++(1,0);
	\draw[-,thick] (B)--++(1,0);
	\draw[-,thick] (C)--++(3,0);
	\draw[-,thick] (D)--++(3,0);
	\draw[fill=tensorblue, thick] (2,-0.5) rectangle (3,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (3,0)--++(1,0);
	\draw[-,thick]  (3,1)--++(1,0);
	\draw[thick](4,1) .. controls (4+1/4,1) and (4+2/4,3/4+1) .. (4,2);
	\draw[thick](4,0) .. controls (4+1/4,0) and (4+2/4,-3/4) .. (4,-1);
	\draw[thick](0,2) .. controls (-1/4,2) and (-2/4,3/4+2) .. (0,3);
	\draw[thick](0,-1) .. controls (-1/4,-1) and (-2/4,-3/4-1) .. (0,-2);
	\draw[thick][red,thick,dashed] (0.2,0.5) rectangle (1.8,2.5);
	\draw[thick][red,thick,dashed] (0.2,0.4) rectangle (1.8,-1.5);
}
\end{align}

\section{}
\begin{align}
\frac{1}{9} \left[
\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
}  +
\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,1.7)--++(2,0);
	\draw[-,thick]  (-0.5,2)--++(2,0);
	\draw[thick](-0.5,1) .. controls (-0.5-1/8,1) and (-0.5-2/4,3/6+1) .. (-0.5,1.7);
	\draw[thick](1.5,1) .. controls (1.5+1/8,1) and (1.5+2/8,3/6+1) .. (1.5,1.7);
} \\
&+\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,-0.7)--++(2,0);
	\draw[-,thick]  (-0.5,-1)--++(2,0);
	\draw[thick](-0.5,0) .. controls (-0.5-1/8,0) and (-0.5-2/8,-3/6) .. (-0.5,-0.7);
	\draw[thick](1.5,0) .. controls (1.5+1/8,0) and (1.5+2/8,-3/6) .. (1.5,-0.7);
}
+
\diagram{
	\draw[-,thick]  (-1,0)--++(2,0);
	\draw[-,thick]  (-1,1)--++(2,0);
}
 \right]
\end{align}

\section{}
\begin{align}
	\diagram{
		\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
		\draw[-,thick]  (1,0)--++(0.5,0);
		\draw[-,thick]  (1,1)--++(0.5,0);
		\draw[-,thick]  (-0.5,0)--++(0.5,0);
		\draw[-,thick]  (-0.5,1)--++(0.5,0);
		\draw[-,thick]  (-0.5,2)--++(2,0);
		\draw[thick](1.5,1) .. controls (1.5+1/4,1) and (1.5+2/4,3/4+1) .. (1.5,2);
		\draw[thick](-0.5,1) .. controls (-0.5-1/4,1) and (-0.5-2/4,3/4+1) .. (-0.5,2);
	}
=
\frac{1}{3} \left[
\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,2)--++(2,0);
	\draw[thick](1.5,1) .. controls (1.5+1/4,1) and (1.5+2/4,3/4+1) .. (1.5,2);
	\draw[thick](-0.5,1) .. controls (-0.5-1/4,1) and (-0.5-2/4,3/4+1) .. (-0.5,2);
} +
\diagram{\draw[-,thick]  (0,0)--++(2,0);}
\right]
\end{align}

\begin{align}
\diagram{
		\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
		\draw[-,thick]  (1,0)--++(0.5,0);
		\draw[-,thick]  (1,1)--++(0.5,0);
		\draw[-,thick]  (-0.5,0)--++(0.5,0);
		\draw[-,thick]  (-0.5,1)--++(0.5,0);
	} =
9 \diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
} -
3\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,1.7)--++(2,0);
	\draw[-,thick]  (-0.5,2)--++(2,0);
	\draw[thick] (-0.5,1) .. controls (-0.5-1/8,1) and (-0.5-2/8,3/6+1) .. (-0.5,1.7);
	\draw[thick] (1.5,1) .. controls (1.5+1/8,1) and (1.5+2/8,3/6+1) .. (1.5,1.7);
} -
3\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,-0.7)--++(2,0);
	\draw[-,thick]  (-0.5,-1)--++(2,0);
	\draw[thick] (-0.5,0) .. controls (-0.5-1/8,0) and (-0.5-2/8,-3/6) .. (-0.5,-0.7);
	\draw[thick] (1.5,0) .. controls (1.5+1/8,0) and (1.5+2/8,-3/6) .. (1.5,-0.7);
}
+
 \diagram{
 	\draw[-,thick]  (-1,0)--++(2,0);
 	\draw[-,thick]  (-1,1)--++(2,0);
 }
\end{align}

\begin{figure}[h]
    \center
\begin{tikzpicture}
\tikzset{swapaxes/.style = {rotate=90,yscale=-1}}

\begin{axis}[
    axis line style = very thick,
    xtick = {1,3},
    xticklabels = {$1$,$10$},
    ytick = {-1,1},
    axis y line=left,
    axis x line=middle,
    xlabel=$h_x$,ylabel=$Energy$,
    ymax=2,
    ymin=-2,
    xmax=3.5,
    xmin=0,
]
\addplot[line width = 2pt,name path= ferro,dashed,red,mark=none,
domain=0:1.01,samples=50]
{-x/1.2-1/2};
\addplot[line width = 2pt,name path=GS, smooth,blue,mark=none, domain=0:14,samples=50]{-x/3-1};
\addplot[line width = 2pt,name path= classical, smooth,black,mark=none, domain=0:14,samples=50]{x/3+1};
\addplot[line width = 2pt,name path= para1,dashed,red,mark=none, domain=1:2.52,samples=50]{x/1.2-2.1};
\addplot[line width = 1pt,smooth,black,mark=none, domain=0:12,samples=50]{-1};
\addplot[line width = 1pt,smooth,black,mark=none, domain=0:12,samples=50]{1};
\addplot[name path=zero,smooth,mark=none, domain=0:12,samples=50]{0};
\addplot [
        thick,
        color=red,
        fill=red,
        fill opacity=0.1
    ]
fill between[
    of=ferro and GS,

    ];
\addplot [
    thick,
    color=blue,
    fill=blue,
    fill opacity=0.1
    ]
fill between[
    of=para1 and GS,

    ];
\addplot [
        thick,
        color=green,
        fill=green,
        fill opacity=0.1
    ]
fill between[
    of=classical and zero,

    ];
\end{axis}
\draw[black,-stealth](.5,.1)--(.5,1.6)node[yshift=-50]{Ferromagnetic Phase};
\draw[black,-stealth](2.2,1.9)--(2.2,1.55)node[xshift=-5,yshift=16]{Ferro Initial State};
\draw[black,-stealth](5.5,4.7)--(5.5,4.35)node[xshift=-5,yshift=18]{Anti-ferro Initial State};
\draw[](4.8,1)node[]{Quantum Paramagnetic Phase};
\draw[](4,3.5)node[]{\large Classical Paramagnetic Phase};
\end{tikzpicture}
\end{figure}
\end{document}
\section{Example}
	\begin{tikzpicture}[inner sep=1mm]
	\foreach \i in {1,...,5} {
		\node[tensor] (\i) at (\i, 0) {$A_\i$};
		\node (\i spin) at (\i, -0.7) {$s_\i$};
		\draw[-,thick] (\i) -- (\i spin);
	};
	\foreach \i in {1,...,4} {
		\pgfmathtruncatemacro{\iplusone}{\i + 1};
		\draw[-,thick] (\i) -- (\iplusone);
	};

	\end{tikzpicture}

	\section{}
	\begin{align}
	\frac{1}{3} \sum_{b,\sigma}
	\diagram{
		\node[tensor] at (1,0) (A) {$P_{b \sigma}$};
		\node[tensor] at (2,0) (B) {$\rho$};
		\node[tensor] at (3,0) (C) {$P_{b \sigma}$};
		\draw[-,thick] (A)--(B);
		\draw[-,thick] (B)--(C);
		\draw[-,thick] (0,0)--node[above] {$x$} ++(A) ;
		\draw[-,thick] (C)--node[above] {$x$} ++(1,0);}
	\end{align}

	\section{}
	\begin{align}
	\diagram{
		\def1/4{1/4}
		\def1/2{1/2}
		\draw[thick](-1/2,0) -- (1/2,0);
		\draw[thick](-1/4,3/4) -- (1/4,3/4);
		\draw[thick](-1/2,0) .. controls (-1/2-1/4,0) and (-1/2-2*1/4,3/4) .. (-1/4,3/4);
		\draw[thick](1/2,0) .. controls (1/2+1/4,0) and (1/2+2*1/4,3/4) .. (1/4,3/4);
		\node[tensor] at (0,0) (A) {$\rho$};} = 1
	\end{align}

	\section{}
	\begin{align}
	\frac{1}{3}\sum
	\diagram{
		\node[tensor] at (1,0) (A) {$P_{b \sigma}$};
		\node[tensor] at (1,1) (B) {$P^T_{b \sigma}$};
		\draw[-,thick] (0,0)--++(A) ;
		\draw[-,thick] (0,1)--++(B) ;
		\draw[-,thick] (A)--++(1,0);
		\draw[-,thick] (B)--++(1,0);
	} =
	\frac{1}{3}
	\left[ \;
	\diagram{
		\draw[-,thick] (0,0)--(0,1);
		\draw[-,thick] (0.5,0)--(0.5,1);}
	+
	\diagram{
		\draw[-,thick] (0,0)--(1,0);
		\draw[-,thick] (0,0.5)--(1,0.5);}
		\;
	\right]
	\end{align}

	\section{}
	\begin{align}
	\frac{1}{3}
	\left[
	\diagram{
		\node[tensor] at (1,0) (A) {$\rho$};
		\draw[-,thick] (0,0)--(A);
		\draw[-,thick] (A)--(2,0);
	}
	+
	\diagram{\draw[-,thick] (0,0)--(0,0);\draw[-,thick] (0,0.20)--(1,0.20);}
	\right]
	\end{align}

	\section{}
	\begin{align}
	\diagram{
		\node[tensor] at (1,0) (A) {$\rho$};
		\draw[-,thick] (0,0)--(A);
		\draw[-,thick] (A)--(2,0);
	}
	=
	3 \diagram{
		\node[tensor] at (1,0) (A) {$CS(\rho)$};
		\draw[-,thick] (0,0)--(A);
		\draw[-,thick] (A)--(2,0);
	} +
	\diagram{\draw[-,thick] (0,0)--(0,0);\draw[-,thick] (0,0.20)--(1,0.20);}

	\end{align}
\section{}
	\begin{align}
	\diagram{
		\def1/2{3}
		\def1/4{1/4}
		\node[tensor] at (1,0) (A) {$P_{b \sigma}$};
		\node[tensor] at (1,1) (B) {$P_{b \sigma}$};
		\node[tensor] at (2.5,0) (C) {$\rho$};
		\draw[-,thick] (0,2)--(4,2);
		\draw[-,thick] (0,0)--++(A) ;
		\draw[-,thick] (0,1)--++(B) ;
		\draw[-,thick] (A)--(C);
		\draw[-,thick] (B)--++(2.5,0);
		\draw[-,thick] (C)--++(1,0);
		\draw[thick] (3.5,0) .. controls (3.5+1/4,0) and (3.5+2/4,3/4) .. (3.5,1);
		\draw[thick] (0,1) .. controls (0-1/4,1) and (0-2/4,1+3/4) .. (0,2);
		\draw[thick][red,thick,dashed] (0.1,-0.5) rectangle (2,1.5);
	}
	\end{align}
\section{}
\begin{align}
\frac{1}{9}\sum _{\sigma, b} \sum_{\sigma', b'}
\diagram{
	\node[tensor] at (1,0) (A) {$P_{b' \sigma'}$};
	\node[tensor] at (1,1) (B) {$P_{b\sigma}$};
	\draw[-,thick] (0,0)--++(A) ;
	\draw[-,thick] (0,1)--++(B) ;
	\draw[-,thick] (A)--++(1,0);
	\draw[-,thick] (B)--++(1,0);
	\draw[fill=tensorblue, thick] (2,-0.5) rectangle (3,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (3,0)--++(1,0);
	\draw[-,thick]  (3,1)--++(1,0);
	\node[tensor] at (4,0) (C) {$P_{b' \sigma'}$};
	\node[tensor] at (4,1) (D) {$P_{b \sigma}$};
	\draw[-,thick] (C)--++(1,0);
	\draw[-,thick] (D)--++(1,0);
}
\end{align}

\section{}
\begin{align}
\diagram{
	\node[tensor] at (1,0) (A) {$P_{b' \sigma'}$};
	\node[tensor] at (1,1) (B) {$P_{b\sigma}$};
	\node[tensor] at (1,-1) (C) {$P_{b' \sigma'}$};
	\node[tensor] at (1,2) (D) {$P_{b\sigma}$};
	\draw[-,thick] (0,3)--++(4,0);
	\draw[-,thick] (0,-2)--++(4,0);
	\draw[-,thick] (0,0)--++(A);
	\draw[-,thick] (0,1)--++(B);
	\draw[-,thick] (0,-1)--++(C);
	\draw[-,thick] (0,2)--++(D);
	\draw[-,thick] (A)--++(1,0);
	\draw[-,thick] (B)--++(1,0);
	\draw[-,thick] (C)--++(3,0);
	\draw[-,thick] (D)--++(3,0);
	\draw[fill=tensorblue, thick] (2,-0.5) rectangle (3,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (3,0)--++(1,0);
	\draw[-,thick]  (3,1)--++(1,0);
	\draw[thick](4,1) .. controls (4+1/4,1) and (4+2/4,3/4+1) .. (4,2);
	\draw[thick](4,0) .. controls (4+1/4,0) and (4+2/4,-3/4) .. (4,-1);
	\draw[thick](0,2) .. controls (-1/4,2) and (-2/4,3/4+2) .. (0,3);
	\draw[thick](0,-1) .. controls (-1/4,-1) and (-2/4,-3/4-1) .. (0,-2);
	\draw[thick][red,thick,dashed] (0.2,0.5) rectangle (1.8,2.5);
	\draw[thick][red,thick,dashed] (0.2,0.4) rectangle (1.8,-1.5);
}
\end{align}

\section{}
\begin{align}
\frac{1}{9} \left[
\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
}  +
\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,1.7)--++(2,0);
	\draw[-,thick]  (-0.5,2)--++(2,0);
	\draw[thick](-0.5,1) .. controls (-0.5-1/8,1) and (-0.5-2/4,3/6+1) .. (-0.5,1.7);
	\draw[thick](1.5,1) .. controls (1.5+1/8,1) and (1.5+2/8,3/6+1) .. (1.5,1.7);
} \\
&+\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,-0.7)--++(2,0);
	\draw[-,thick]  (-0.5,-1)--++(2,0);
	\draw[thick](-0.5,0) .. controls (-0.5-1/8,0) and (-0.5-2/8,-3/6) .. (-0.5,-0.7);
	\draw[thick](1.5,0) .. controls (1.5+1/8,0) and (1.5+2/8,-3/6) .. (1.5,-0.7);
}
+
\diagram{
	\draw[-,thick]  (-1,0)--++(2,0);
	\draw[-,thick]  (-1,1)--++(2,0);
}
 \right]
\end{align}

\section{}
\begin{align}
	\diagram{
		\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
		\draw[-,thick]  (1,0)--++(0.5,0);
		\draw[-,thick]  (1,1)--++(0.5,0);
		\draw[-,thick]  (-0.5,0)--++(0.5,0);
		\draw[-,thick]  (-0.5,1)--++(0.5,0);
		\draw[-,thick]  (-0.5,2)--++(2,0);
		\draw[thick](1.5,1) .. controls (1.5+1/4,1) and (1.5+2/4,3/4+1) .. (1.5,2);
		\draw[thick](-0.5,1) .. controls (-0.5-1/4,1) and (-0.5-2/4,3/4+1) .. (-0.5,2);
	}
=
\frac{1}{3} \left[
\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,2)--++(2,0);
	\draw[thick](1.5,1) .. controls (1.5+1/4,1) and (1.5+2/4,3/4+1) .. (1.5,2);
	\draw[thick](-0.5,1) .. controls (-0.5-1/4,1) and (-0.5-2/4,3/4+1) .. (-0.5,2);
} +
\diagram{\draw[-,thick]  (0,0)--++(2,0);}
\right]
\end{align}

\begin{align}
\diagram{
		\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$\rho$};
		\draw[-,thick]  (1,0)--++(0.5,0);
		\draw[-,thick]  (1,1)--++(0.5,0);
		\draw[-,thick]  (-0.5,0)--++(0.5,0);
		\draw[-,thick]  (-0.5,1)--++(0.5,0);
	} =
9 \diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
} -
3\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,1.7)--++(2,0);
	\draw[-,thick]  (-0.5,2)--++(2,0);
	\draw[thick] (-0.5,1) .. controls (-0.5-1/8,1) and (-0.5-2/8,3/6+1) .. (-0.5,1.7);
	\draw[thick] (1.5,1) .. controls (1.5+1/8,1) and (1.5+2/8,3/6+1) .. (1.5,1.7);
} -
3\diagram{
	\draw[fill=tensorblue, thick] (0,-0.5) rectangle (1,1.5) (C)  node[pos=0.5] {$CS(\rho)$};
	\draw[-,thick]  (1,0)--++(0.5,0);
	\draw[-,thick]  (1,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,0)--++(0.5,0);
	\draw[-,thick]  (-0.5,1)--++(0.5,0);
	\draw[-,thick]  (-0.5,-0.7)--++(2,0);
	\draw[-,thick]  (-0.5,-1)--++(2,0);
	\draw[thick] (-0.5,0) .. controls (-0.5-1/8,0) and (-0.5-2/8,-3/6) .. (-0.5,-0.7);
	\draw[thick] (1.5,0) .. controls (1.5+1/8,0) and (1.5+2/8,-3/6) .. (1.5,-0.7);
}
+
 \diagram{
 	\draw[-,thick]  (-1,0)--++(2,0);
 	\draw[-,thick]  (-1,1)--++(2,0);
 }
\end{align}

\begin{figure}[h]
    \center
\begin{tikzpicture}
\tikzset{swapaxes/.style = {rotate=90,yscale=-1}}

\begin{axis}[
    axis line style = very thick,
    xtick = {1,3},
    xticklabels = {$1$,$10$},
    ytick = {-1,1},
    axis y line=left,
    axis x line=middle,
    xlabel=$h_x$,ylabel=$Energy$,
    ymax=2,
    ymin=-2,
    xmax=3.5,
    xmin=0,
]
\addplot[line width = 2pt,name path= ferro,dashed,red,mark=none,
domain=0:1.01,samples=50]
{-x/1.2-1/2};
\addplot[line width = 2pt,name path=GS, smooth,blue,mark=none, domain=0:14,samples=50]{-x/3-1};
\addplot[line width = 2pt,name path= classical, smooth,black,mark=none, domain=0:14,samples=50]{x/3+1};
\addplot[line width = 2pt,name path= para1,dashed,red,mark=none, domain=1:2.52,samples=50]{x/1.2-2.1};
\addplot[line width = 1pt,smooth,black,mark=none, domain=0:12,samples=50]{-1};
\addplot[line width = 1pt,smooth,black,mark=none, domain=0:12,samples=50]{1};
\addplot[name path=zero,smooth,mark=none, domain=0:12,samples=50]{0};
\addplot [
        thick,
        color=red,
        fill=red,
        fill opacity=0.1
    ]
fill between[
    of=ferro and GS,

    ];
\addplot [
    thick,
    color=blue,
    fill=blue,
    fill opacity=0.1
    ]
fill between[
    of=para1 and GS,

    ];
\addplot [
        thick,
        color=green,
        fill=green,
        fill opacity=0.1
    ]
fill between[
    of=classical and zero,

    ];
\end{axis}
\draw[black,-stealth](.5,.1)--(.5,1.6)node[yshift=-50]{Ferromagnetic Phase};
\draw[black,-stealth](2.2,1.9)--(2.2,1.55)node[xshift=-5,yshift=16]{Ferro Initial State};
\draw[black,-stealth](5.5,4.7)--(5.5,4.35)node[xshift=-5,yshift=18]{Anti-ferro Initial State};
\draw[](4.8,1)node[]{Quantum Paramagnetic Phase};
\draw[](4,3.5)node[]{\large Classical Paramagnetic Phase};
\end{tikzpicture}
\end{figure}
\end{document}
\section{Introduction}

Simulation of strongly correlated electrons in the context of quantum chemistry and condensed matter physics is one of the potential areas in which quantum computers will have a significant advantage over their classical counterparts \cite{Cao_2019, Wecker_2015, Ma_2020}. Strong interaction between the ostensibly simple electrons can give rise to novel phases including high-temperature superconductivity \cite{highTc_1, highTc_2}, strange metallic behavior\cite{mcgreevy_2010}, fractional excitations \cite{Laughlin_1999}, and quantum spin liquids\cite{Savary_2016}. Condensed matter physics aims to understand these novel behaviors by studying their phase diagrams. However, failure of perturbation theory and exponential scaling of the Hilbert space for strongly correlated electrons presents a formidable challenge to classical simulation methods such as exact diagonalization, density matrix renormalization group(DMRG) \cite{dmrg1, dmrg2, dmrg3}, quantum Monte-Carlo \cite{qmc} and dynamical mean-field theory\cite{dmft1}.
Whereas, this same challenge provides an exciting opportunity for near-term quantum computers.

\comment{2. Opportunities and challenges for quantum computers. Need to change the organization of this paragraph?
Feynman\cite{Feynman1982} conjectured in 1982 that a quantum computer will be able to efficiently simulate quantum mechanical systems more efficiently than a classical computer. The power of a quantum computer lies in its control over the exponentially scaling Hilbert space of qubits, thus turning the curse of dimensionality into a blessing. }

Harnessing the power of a quantum computer to simulate quantum systems \cite{Feynman1982} requires (i) algorithms that can be executed in a reasonable time and (ii) the ability to learn from quantum experiments without exponentially many measurements. Studying the phases via ground state preparation is a QMA-complete problem \cite{Gorman_2021, Gharibian2019complexity, Bookatz2014QMA, kemp2004complexity}, which cannot be carried out in a reasonable time, even with quantum resources. However, performing dynamics on a quantum state is known to be a BQP-hard problem \cite{Baez2020dynamical, Rudi2020approximating}, possible within polynomial time. Likewise, it has been shown that shadow tomography\cite{Aaronson2018} methods such as classical shadows\cite{huang2020predicting, huang2021provably, Huang2022Learning} are effective at predicting properties using very few measurements.

In classical statistical mechanics, ergodicity provides a link between time averages and statistical averages \cite{Sethna_2012}. In quantum mechanics, unitary time evolution retains memory of the initial state, so the link is different yet still enabled by equilibration of macroscopic observables. This is made possible by certain general conditions on an initial state\cite{von2010proof,  Gogolin_2016,Linden_2009,Reimann_2012,Pintos_2017}, which allow these observables to equilibrate quickly.
Thus, if we can find interesting problems and initial states for which the order parameter (a macroscopic observable that characterizes the existence of a phase) equilibrates before the qubits decohere, we can use quantum computers to perform microcanonical dynamics and study the regions of the phase diagram inaccessible via classical computing resources.

In this manuscript, we present an approach to identifying phase diagrams of strongly correlated systems via dynamics simulations on a quantum computer. It consists of i) identifying an initial state, ii) evolving this state in time using a suitable algorithm, iii) using classical shadow tomography technique to convert the quantum state to classical data, and iv) applying unsupervised machine learning methods to discover phases of matter. We illustrate it with the 1-dimensional Transverse Field Ising Model(1DTFIM).

\begin{figure*}[t]
\includegraphics[scale=.14]{Schematic.pdf}
\caption{Schematic overview of our study. (a) Classical shadows(CS) of ground states  and time averaged classical shadows(TACS) from dynamics of a time-reversal invariant GHZ state are generated using quantum simulation. (b) The classical data from quantum simulation is then fed into diffusion map, an unsupervised machine learning algorithm to learn the phases.}
\label{Fig:Schematic_diagram}
\end{figure*}

A schematic overview of our approach is shown in Fig. \ref{Fig:Schematic_diagram}. It consists of using diffusion maps~\cite{coifman2006diffusion, de2008introduction}, an unsupervised machine learning(UL) algorithm,  to machine learn phase features from unlabeled phase diagram data generated by a quantum computer. Specifically, we apply diffusion maps to CS of ground states from a 100 qubit simulation and show they successfully learn the magnetic phase transition of 1DTFIM (see section II). It does so by \emph{discovering}  the relevant parameters such as magnetization and transverse magnetic field. These results follow similar machine-learning-from-CS success stories in other contexts\cite{ huang2020predicting,  huang2021provably, randomized_toolbox_2022, Huang2022Learning}. We then show (section III) diffusion maps can also identify the phase transition and quantum critical region \emph{when trained on dynamics data from just 20 qubits}. This data is in the form of a representation we introduce: time-averaged classical shadows (TACS)---a shadow tomographic\cite{Aaronson2018} representation of von Neumann's time-averaged density matrix\cite{von2010proof}. TACS are a straightforward extension of CS to capture quantum microcanonical dynamics of finite-size finite-time quantum simulations. Our results thus validate the exciting possibility of studying the phase diagrams of strongly correlated electrons using quantum computers.

\section{Von Neumann's Microcanonical Ensemble}

A central goal of quantum computing is to build qubits that are completely isolated from their environment. While this is not the case today, the current development of quantum error correction techniques\cite{Knill_2005,nielsen_chuang_2010} suggests it is in our future. Simulating quantum systems on a quantum computer will therefore take place \emph{within the microcanonical ensemble}. But quantum microcanonical dynamics, the evolution of a closed quantum system under Schr\"odinger's equation, does not directly produce the microcanonical ensemble.

Following von Neumann's 1929 paper\cite{von2010proof,goldstein2010normal} on the quantum ergodic theorem, it is straightforward to derive the equilibrium distribution of a quantum system using density matrices. Assuming we start from an initial state $|\psi(0)\rangle$ and evolve under a Hamiltonian $H$ via a quantum circuit algorithm to $|\psi(t)\rangle$, the equilibrium distribution is captured by the time-averaged density matrix
\begin{equation}\label{eq:rhoT}
    \rho_T\! =\! \frac{1}{T}\int_0^T\!\! dt |\psi(t)\rangle\langle\psi(t)| \!\underset{T\to\infty}{\longrightarrow}\! \omega\! =\! \sum_k P_k|\psi(0)\rangle\langle\psi(0)|P_k,
\end{equation}
where $P_k$ is a projector onto the $k^\text{th}$ degenerate subspace of the energy eigenvalues i.e. $P_k = \sum_{n\in k}|E_n\rangle\langle E_n|$.
\comment{
We can understand this expression from the snapshot picture of statistical mechanics. An ensemble of states/systems can be thought of as equivalently a history of one system where we record the state of the system in a series of snapshots over time. If the frequency of snapshots is sufficiently low, then each snapshot will be statistically independent of the previous one leading to a probability of $dt/T$ for each snapshot and the above time averaged expression for the density matrix of this ensemble. }
Existence of this equilibrium density matrix $\omega$ results in the ergodic principle that long-time averages of observables should  be captured by the statistical averaging with respect to $\omega$. Specifically, in the Schr\"odinger picture,
\begin{align}
  \langle O\rangle_T &=  \frac{1}{T}\int_0^T dt
    \langle\psi(t)|O|\psi(t)\rangle \\
    &= \frac{1}{T}\int_0^T dt \text{Tr}\left( |\psi(t)\rangle\langle\psi(t)|O\right) \\
    &= \text{Tr} \left(\rho_T O\right) \underset{T \to \infty}{\longrightarrow} \Tr(\omega O)
\end{align}
Thus, the time-averaged density matrix is the link between time averages and statistical averages.

The connection to the microcanonical ensemble, obtained by quantizing the classical microcanonical ensemble, is achieved by taking the thermodynamic limit, measuring only coarse-grained observables, and requiring non-degenerate energy level spacings. Under these circumstances, von Neumann obtained
\begin{equation}
    \omega\to\rho_{mc} = \frac{1}{\Omega}\sum_{E < E_n < E + \Delta} |E_n\rangle\langle E_n|
\end{equation}
where all states in an energy window $[E,E+\Delta]$ are included in the sum for some window size $\Delta$, and $\Omega$ is the phase-space volume of this constant energy shell. So $\rho_T$ approaches $\rho_{mc}$ in the limit of long times and large systems sizes, but for any finite $T$ and system size, $\rho_T$ is the actual microcanonical ensemble defined as the ensemble produced by the Schrödinger equation.

In numerics, two things are of utmost importance---reaching the thermodynamic limit and understanding the errors that may result along the way. Using the above logic, an equilibrium phase diagram with finite size errors associated with a reasonable choice of initial conditions is achievable. For each initial condition, it's possible to check whether the time series was run long enough for relevant observables to reach equilibrium (not all observables equilibrate until an exponential amount of time has passed\cite{Gogolin_2016, Short_2012}). If so, the only errors remaining are of finite size and will be different for each choice of initial conditions. A family of initial conditions that all yield the same microcanonical ensemble, when scaled up, is, therefore, a key resource to study equilibrium statistical mechanics on a quantum computer. Although the equilibrium properties are independent of the initial state, the time required to reach equilibration depends on the initial state via the effective dimension given by $d_\text{eff} = 1/\sum p_k^2$, where $p_k =|\bra{\psi} \ket{E_k}|^2$. Here, $p_k$ measures how many energy eigenstates have a significant weight in $|\psi\rangle$ and subsequently, a large effective dimension results in rapid equilibration with a bound on the equilibration time given by the second largest occupancy \cite{Gogolin_2016, Short_2012}. The distribution of the $p_k$'s also affects equilibration times\cite{margolus2021counting}. These arguments suggest we want to look for an initial state that exhibits a small overlap with most energy levels or a macroscopic occupation of a single energy level. In practice, we are concerned with the distinguishability of the time-averaged state vs. the equilibrium state in terms of experimentally relevant observables. Therefore, we need not worry about the equilibration of the full-density matrix as long as the equilibration of the relevant order parameters can be achieved.

In order to study von Neumann's microcanonical ensemble, we need to compare it to what we already know. As such, in the next section, section \ref{sec:groundstatedata}, we will turn to a ground state study using the methods we will employ for our thermodynamics study to validate these methods first, before applying them to the study of TACS.

\section{Ground state data}\label{sec:groundstatedata}
To verify our approach to phase classification, we first apply it to ground states of the ferromagnetic 1DTFIM defined by the Hamiltonian:
\begin{equation}
    H = -\sum_{\langle i, j \rangle} Z_i Z_j + h_x\sum_{i} X_i,
\end{equation}
where $\langle.\rangle$ denotes the nearest neighbors, $Z_i$ is the Pauli operator in the z-direction, and $h_x$ is a parameter proportional to the transverse magnetic field. At zero temperature, this model has a ferromagnetic phase for $\abs{h_x}<1$, and a paramagnetic phase for $\abs{h_x}>1$. The ground state study aims to detect the second-order phase transition at $h_x = 1$ for a $100$ site  1DTFIM using diffusion maps. The ground states were generated using density matrix renormalization group (DMRG) \cite{dmrg1, dmrg2, dmrg3} with ITensor package\cite{itensor}.
Training datasets were generated using two kinds of measurements on these ground states-  (i) computational basis measurements to obtain the Z dataset and (ii) measurements on a random Pauli basis to obtain the CS dataset.

\subsection{Classical Shadows}

\begin{figure*}[ht]
\begin{center}
\includegraphics[scale=0.9]{CS_diagramatic.pdf}
\caption{Diagrammatic description of classical shadows showing a linear relationship between $S_N[\rho]$ and estimator $\sigma_N[\rho]\approx\rho$. (a) The full-density matrix $\rho$ can be approximated by summing over reduced classical shadows with a coefficient that grows exponentially in the number of remaining qubits. (b) A classical shadow $S_N[\rho]$ is obtained by summing over $N$ measurement outcomes on random Pauli bases. For a given $N$, a reduced density matrix, that involves smaller coefficients in the expansion, can be approximated more accurately compared to the full density matrix.}
\label{Fig:CS_diagramatic}
\end{center}
\end{figure*}

\comment{\michael{In principle, the choice of shadow tomography impacts what observables are visible in the data. Similarly, in finite-time microcanonical dynamics simulations, coarse grained observables equilibrate while fine grained ones do not. So, to study quantum dynamics on a quantum computer, \sout{we should} \eric{one must} adjust the choice of shadow tomography to best observe the behavior of equilibrating observables.}}

Obtaining any useful information from a quantum computer requires performing measurements on a quantum state, which is destructive to the quantum information by nature. Since the dimension of the Hilbert space increases exponentially in the number of qubits, a naive strategy to learn the state requires an exponentially large number of copies. Aaronson\cite{Aaronson2018} introduced an alternative method using the notion of shadow tomography, an approximate classical description of the quantum state, in which $M$ properties of a quantum state can be estimated with error $\epsilon$ by only $O(\frac{M}{\epsilon^2})$ copies of the state. We can think of a shadow as an approximation of a quantum state $\rho$ by summing over measurement outcomes $x$, obtained by performing measurements on bases $b$ for a quantum state $x$ i.e.
\begin{equation}\label{eq:cs}
    S[\rho] = \sum_{b,x}P(b)
        \mathcal{P}_{b,x}
        \rho
        \mathcal{P}_{b,x},
\end{equation}
where $ \mathcal{P}_{b,x}$ is a projector onto the measurement outcome $x$ on basis $b$, and $P(b)$ is the probability of choosing $b$.

Based on this notion, Huang et al.\cite{Huang2019predicting, huang2021provably, Huang2022Learning} developed an algorithm called classical shadows and showed that it is highly successful at learning the properties of a many-body system. Two kinds of measurement protocols were proposed to construct classical shadows-  (i) random Clifford measurements on the entire Hilbert space; (ii) random single-qubit Pauli measurements. Protocol (ii) results in very shallow measurement circuits and thus is more suitable for the NISQ-era\cite{PreskillNISQ} hardware. \comment{With random Pauli measurements, they show that predicting linear functions with error $\epsilon$ requires at least $\text{log}(M) \text{max}_i \text{tr}(O_i^2)/\epsilon^2$ copies of the state.} After measuring each of the qubits in some random Pauli basis $X$, $Y$ or $Z$ with outcomes $\pm1$, the post-measurement wavefunction is given by the product state $\left|s^{(n)}\right\rangle=\bigotimes_{l=1}^{L}\left|s_{l}^{(n)}\right\rangle$. Here, $\left|s_{l}^{(n)}\right\rangle \in\{|0\rangle,|1\rangle,|+\rangle,|-\rangle,|\mathrm{i}+\rangle,|\mathrm{i}-\rangle\}$ is a Pauli basis state to which the $l^{th}$ qubit has collapsed. The classical shadow $S_N[\rho]$ is obtained by summing over $N$ such randomized measurement outcomes as follows (also see Fig. \ref{Fig:CS_diagramatic}(b))
\begin{align}
    S_N[\rho] &= \frac{1}{N} \sum_{n=1}^N \ket{s^{(n)}} \bra{s^{(n)}}\\
    &= \frac{1}{N} \sum_{n=1}^N \ket{s_1^{(n)}} \bra{s_1^{(n)}} \otimes \cdots \otimes  \ket{s_L^{(n)}}\bra{s_L^{(n)}}.
\end{align}

The underlying quantum state $\rho$ can be approximated by adding the reduced classical shadows (see Fig. \ref{Fig:CS_diagramatic}(a)).  This sum simplifies to the following expression from Ref. \cite{Huang2019predicting, huang2021provably, Huang2022Learning}
\begin{equation}\label{eq:csshots}
    \rho \approx \sigma_{N}(\rho)=\frac{1}{N} \sum_{n=1}^{N} \sigma_{1}^{(n)} \otimes \cdots \otimes \sigma_{L}^{(n)},
\end{equation}
where
\begin{equation}
\sigma_{l}^{(n)}=3\left|s_{l}^{(n)}\right\rangle\left\langle s_{l}^{(n)}\right|-\mathbb{I}.
\end{equation}

The definition of $S_N(\rho)$ presented above is different from Refs. \onlinecite{Huang2019predicting, huang2021provably} which defines it to be\comment{;} the dataset of shots itself and not the density matrix obtained from these shots. But both are useful for Fig. \ref{Fig:CS_diagramatic} shows the density matrix $S_N(\rho)$ defined above is linearly related to the estimator $\sigma_N(\rho)$ of the quantum state $\rho$ obtained by Refs. \onlinecite{Huang2019predicting, huang2021provably, Huang2022Learning}. A derivation of Fig. \ref{Fig:CS_diagramatic} using tensor network diagrams is presented in Appendix A.

Although estimating the exact density matrix requires $N \rightarrow \infty$, we still desire to predict various linear as well as nonlinear functions of $\rho$ (e.g., $\Tr(O\rho)$ and $\Tr (\rho \text{log}(\rho))$ respectively). This can be achieved with $N \propto \log(L) 4^k/\epsilon^2$ copies of the state,  where $k$ is the locality of operator $O$ \cite{randomized_toolbox_2022}. It was shown in Ref. \cite{huang2021provably} that classical machine learning algorithms can efficiently predict the ground state properties of gapped Hamiltonians in finite spatial dimensions after learning the classical shadows from a training set. An example of interest is classifying the quantum phases of matter. Classifying the symmetry-breaking phases is conceptually simple because it involves calculating $\text{tr}(\rho O)$ for some k-local observable $O$, such that $\text{tr}(\rho O) \geq 1 \: \forall \: \rho \in \text{ phase A}$ and  $\text{tr}(\rho O) \leq -1 \: \forall \: \rho \in \text{ phase B}$.

In contrast to the symmetry-breaking phases, capturing continuous phase transitions involves nonlinear observables like thermodynamic potentials. These observables involve measuring a nonlinear function of $\rho$, for example, the entropy. Learning such nonlinear functions requires an expressive ML model. Thus, it is beneficial to devise a mapping from classical shadows to a high-dimensional feature space that includes the polynomial expansion of many-body reduced density matrices and encapsulate this feature space in a kernel function. Using such kernels then establishes a rigorous guarantee that a classical ML model can efficiently classify nonlinear functions, phase transitions, and even topological phases of matter\cite{huang2021provably}.

\subsection{Machine Learning Method: Diffusion Maps}

\comment{Classical shadows provide an approximate information about the full state of a quantum computer where the accuracy of the approximation depends just on the number of shots[ref]. But rather than work directly with this approximation, such as by computing observables that are accurately captured by the data, we will apply machine learning (ML) methods whose success depends on this approximation. As we will see, these ML methods can still identify macroscopic parameters of interest, such as the z-magnetization and $h_{x}$ values of the quantum states, as well as identify phases and phase transitions, even though the data we feed it (the classical shadows) contain far less information than the full description of the quantum states.}

In this paper, we utilized unsupervised machine learning methods called diffusion maps \cite{coifman2006diffusion, de2008introduction} to extract features from the shadows data. A diffusion map is a nonlinear dimensionality reduction technique that relies on learning the underlying manifold from which the data points have been generated. Recently, this method was used to identify phases in systems with complex order parameters, which are difficult to learn using linear methods (such as principal component analysis (PCA)\cite{carleo2019machine}). Examples of such phase identification studies include; topological phases and phase transitions \cite{rodriguez2019identifying}, incommensurate phases, and many-body localized phases in quantum systems \cite{lidiak2020unsupervised}.

In this application of diffusion maps, one can imagine a random walk on a dataset  $X\left(x_{1}, x_{2},\ldots, x_{N}\right)$, where the $x_{n}$ are estimators $\sigma_N(\rho)$ of density matrices $\rho$ obtained from different points in the phase diagram. The transition probability $P(j\backslash i)$ of jumping from $x_{i}$ to $x_{j}$ in a single ``timestep'' is proportional to the kernel function $k(x_i,x_j)$, which is a similarity measure between the two data points. Here we use the classical shadow kernel function prescribed in~\cite{huang2021provably}, defined to be for two points $x$ and $\tilde x$ in the dataset
\begin{equation}\label{eq2}
  \begin{aligned}[b]
    & k^{(shadow)}(x,\tilde x) = k^{(shadow)}(\sigma_N(\rho),\sigma_N(\tilde\rho)) =  {}\\
    & \exp{\!\left(\!\sum_{n,n'=1}^N\!\frac{\tau}{N^2}\!     \exp{\!\left(\!\frac{\gamma}{L}\!\sum_{l=1}^L\! \Tr\Big(\sigma_{l}^{(n)} \tilde\sigma_{l}^{(n')}\!\Big)\!\right)}\! \right)}.
  \end{aligned}
\end{equation}
Then, we construct a transition probability matrix $P$ such that
\begin{equation}
\label{eqn:P-matrix}
   {P(j\backslash i
   }
   = \frac{k^{(shadow)}(i,j)}{\sum_l k^{(shadow)}(i,l)}.
\end{equation}

After $t$ timesteps of the random walk, the transition probabilities are given by the matrix $P^t$, where $P_{ij}^t$ gives the probability of going from $x_{i}$ to $x_{j}$ in $t$ timesteps, it's a sum of the probabilities associated with all of the possible paths to go from $x_{i}$ to $x_{j}$ in $t$ timesteps. As $t$ increases, the diffusion process unfolds, where data points situated along the overall geometric structure of the dataset become more strongly connected because of the abundance of strongly connected intermediate points along the way. The parameter 'diffusion distance' is defined to quantify this idea of connectivity between two data points given by

\begin{equation}
    D_{t}^2(x_{i},x_{j})^2 = \sum_{m=1}^N|P_{im}^t-P_{mj}^t|^2
\end{equation}

where the bigger the diffusion distance, the weaker the connection between them. The data points are then mapped onto a new `diffusion space' so that the diffusion distance in data space is equal to the Euclidean distance in this new space. It can be shown \cite{coifman2006diffusion} that such a map can be made in the following way:
\begin{equation}\label{eq5}
    x_{i}\rightarrow y_{i} = \begin{bmatrix}
       \lambda_{1}^t\psi_{1}(i),  &
       \lambda_{2}^t\psi_{2}(i), &
       \cdots, &
       \lambda_{N-1}^t\psi_{N-1}(i)
    \end{bmatrix}
\end{equation}
where $\lambda_{k}$ and $\psi_{k}$ are eigenvalues and right eigenvectors of the matrix $P$, $\psi_{k}(i)$ is the $i$-the element of the $k$-th eigenvector. $\psi_{k}$ form an orthogonal set of basis vectors for the diffusion space.
Then the diffusion distance is,
\begin{equation}\label{eq6}
    D_{t}^2(x_{i},x_{j}) = |y_{i} - y_{j}|^2 \\
        = \sum_{k=1}^{N-1} \lambda_{k}^{2t}\left[\psi_{k}(i) - \psi_{k}(j)\right]^2
\end{equation}
In equations (\ref{eq5}) and (\ref{eq6}) the $k=0$ component is ignored ($\psi_{0}(i) = \frac{1}{\sqrt{N}},\; \lambda_{0}=1$), as it is constant for all $i$, hence trivial. Equation (\ref{eq6}) suggests a dimensionality reduction, as the terms with bigger $\lambda_{k}$ will dominate the sum increasingly as $t\rightarrow\infty$. So  only the $m$ diffusion space dimensions associated with the bigger eigenvalues ($m << N$) are kept. The data points are then projected onto this m-dimensional space to reveal features in the dataset.

\subsection{Phase Classification of Ground States}
For the ground state study, we used two datasets generated via two different shadow tomographies. The first (Z dataset) contains qubit measurements only along the Z-axis [in Eq. \ref{eq:cs}, $P(b)=1$ for $b=Z$, $P(b)=0$ for all other b]. While the other (CS-dataset) has randomized Pauli measurements using the Classical Shadows method [$P(b)=\frac{1}{3}$ for $b \in \{X, Y, Z\}$]. Both of them contain 200 100-spin 1DTFIM ground state shadows, each state for a different $h_{x}$ value ($h_{x}$ ranging from 0.1 to 100).  Since the $Z$-magnetization is the order parameter for TFIM, the UL algorithm should be able to learn the phases of the model from the Z dataset. Using this knowledge, we compare the Z dataset and the CS data to see if the algorithm can successfully identify phases in it.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth, height=0.3\textheight]{gr7.pdf}
\caption{Learning phases from ground state (CS and Z) data. (a) and (c) show the 10 largest eigenvalues of the $P$-matrix (excluding the trivial $k=0$) for (a) Z dataset and (c) CS dataset. (b) Z datapoints in 5D diffusion space visualized in 2D, using metric MDS. Clustering clearly emerges based on the two phases of the model. (d) CS datapoints in 2D diffusion space, the figure reveals the symmetry-breaking phase transition. In this case, there is a direct correlation of relevant parameters, the z-magnetization and the $h_{x}$-values, with machine-learned diffusion, coordinates dc1 and dc2 respectively. Notice the CS-data shows three clusters in (d), while Z-data shows only two in (b), which is a consequence of the shadow kernel function for the former is not spin-flip symmetric while that of the latter is.}
\label{ground}
\end{figure}

By deploying diffusion maps armed with the shadows kernel function (Eq. \ref{eq2}) as our UL model, we are able to identify the phases from both the Z data and the CS data. In both cases, we set $\tau=1$, $\gamma=1$ in Eq. \ref{eq2}. For the Z data, we chose the first five non-trivial eigenvectors as the diffusion space basis vectors because the $P$-matrix (Eq. \ref{eqn:P-matrix}) eigenvalue spectrum shows the first five eigenvalues to be larger than others (Fig. \ref{ground}(a)). Mapping the states from Z data onto this five-dimensional diffusion space, we found that clear clustering emerges based on the phases of the states. We used multidimensional scaling (MDS), a dimensionality reduction method \cite{borg2005modern} that seeks to preserve point-to-point distances, to project these states onto a 2D plane. We see a clear separation of the two phases even on this 2D reduced space (Fig. \ref{ground}(b)), indicating the machine's success in identifying the two phases.

From the CS data, the unsupervised learning algorithm was also able to learn about the phases and the underlying parameters of the model. We can see in Fig. \ref{ground}(c) that the $P$ matrix has two non-trivial eigenvalues larger than the rest. The eigenvectors corresponding to these two eigenvalues are the basis vectors of the reduced diffusion space. Figure (\ref{ground}(d)) shows all of the ground state classical shadows projected onto this two-dimensional plane. It shows three groups- the top left and the bottom left are the all-up and all-down states, whereas the group on the center-right are states in the disordered phase. This closely resembles the spontaneous symmetry-breaking phase transition of 1DTFIM\cite{sachdev_2011}. The learned diffusion coordinates (dc1 and dc2) have direct correlations with the magnetization $<M_{z}>$ and the field values $h_{x}$ respectively, as shown in appendix C (Fig. \ref{fig:correlation}).

The above-mentioned clustering is dependent on the number of snapshots $N$ for a given state. However, the dimensionality reduction and the subsequent clustering will settle down after a minimum value of $N$ has been reached ($N_{c}$). Increasing the value of $N$ beyond that point does not change the results in any significant way, and in our case ($N_{c}$ (Z data)$< N_{c}$ (CS data)), it should be easier for the algorithm to learn the phase space structure from the Z-data (hence fewer snapshots are required).

Despite the successful identification of the ground state phases with the CS+ML model, the fact remains that the problem of calculating ground-states is a QMA-hard problem \cite{Bookatz2014QMA}. While there is no shortage of algorithms designed to generate ground-states, the necessity of a scalable model can be better shown with a known BQP problem; revealing the microcanonical phase diagram from dynamics data.

\section{Microcanonical Dynamics of the 1DTFIM}
Consider now Fig. \ref{fig:phasediagram}, the thermodynamic phase diagram of the 1DTFIM relating internal energy $E=\langle\psi(0)|\hat H|\psi(0)\rangle$ to the transverse magnetic field $h_x$. This phase diagram is relevant for a microcanonical dynamics study governed by the entropy $S(E,h_x)$. This phase diagram exists even for a simulation over a finite time $T$ and with a finite number of spins $N$ but with finite $T$, finite $N$ errors that round the phase transitions and render the results dependent on initial conditions. We present in this figure our expectations for the phase diagram in this context, pointing out where the phase diagram will sharpen in the thermodynamic limit. We further highlight the path through the phase diagram carried out by our simulations below, showing that we expect it to cross the quantum critical region and so be sensitive to the phase diagram at a rounded level.

\begin{figure}
    \centering
    \includegraphics{PhaseDiagramTikz.pdf}
    \caption{A sketch of the expected 1DTFIM phase diagram at finite $T$ and finite $N$ as a function of internal energy $E$ and transverse magnetic field $h_x$. This diagram is a modification of the canonical ensemble representation of the phase diagram in Ref. \onlinecite{wu2018crossovers}, adapted to the microcanonical ensemble. The spectrum is mirror symmetric about $E=0$ due to the chiral symmetry $\mathcal{C}=ZYZYZY...$ \comment{that anticommutes with the Hamiltonian.}}
    \label{fig:phasediagram}
\end{figure}

\comment{\gaurav{\emph{Where is a good place to add this?} A similar algorithm called thermal shadows has been recently proposed in literature\cite{Predicting_Gibbs}. However, this requires a large number of controlled unitaries and optimization for the quantum signal processing\cite{QSP} subroutine. There is also quantum imaginary time evolution(qite) \cite{QITE}, which approximates a k-local imaginary time evolution operator (non-unitary) by C-local unitary where $C>k$. This involves approximating all the non-unitary operators corresponding to each term in the Hamiltonian by using a larger unitary, whose size grows with the correlations(C)  in the system. Our algorithm seems very straightforward compared to this, so will potentially replace qite. There is also another algorithm based on the zeros of partition function \cite{Partition_zeros}.}}

\subsection{Results}
With the knowledge of expected microcanonical dynamics (see Fig. \ref{fig:phasediagram}), we now present an effective classical representation of quantum data obtained from a microcanonical dynamics on a quantum computer. Noticing that the time averaged integral amounts to an expectation value of the pure state density matrix $|\psi(t)\rangle\langle\psi(t)|$ over the probability distribution $P_T(t) = (1/T)\left(\Theta(T-t\rangle)-\Theta(-t)\right)$, where $\Theta(x)$ is the Heaviside step function, we define \emph{time averaged classical shadows}(TACS) by the map
\begin{equation}\label{eq:tacs}
    TACS\left[\rho \right] = \lim_{T\to\infty}\int dt \sum_{b,\sigma}P_T(t)P(b)
        \mathcal{P}_{b,\sigma}
        |\psi(t)\rangle\langle\psi(t)|
        \mathcal{P}_{b,\sigma}
\end{equation}
Hence by sampling the joint probability distribution $P_T(t)P(b)$ to obtain $(t_i,b_i)$, $i=1\ldots N$, and then measuring one shot $\sigma_i$ from $|\psi(t_i)\rangle$ in basis $b_i$ we obtain a finite-shot TACS via
\begin{equation}\label{eq:tacsshots}
    TACS_N\left[\rho \right] = \comment{\lim_{T\to\infty}}
    \comment{\lim_{N\to\infty}
    \frac{1}{N}}
    \sum_{i=1}^{N}
    |t_ib_i\sigma_i\rangle\langle t_ib_i\sigma_i|
\end{equation}
We ran quantum dynamics of 1DTFIM using the TDVP algorithm\cite{tdvp,timeevomps} starting from the GHZ state $\ket{\psi(0)} = \frac{\ket{000\cdots} + \ket{111\cdots}}{\sqrt{2}}$to generate TACS data from 500 randomly sampled dimensionless time values between $t=10.0$ to $t=20.0$ and 187 randomly sampled $h_x$ field values between $h_x = 0.1$ to $h_x = 10.0$. For each $h_{x}$, a TACS was obtained from the $N = 500$ shots at random times over the fixed $T$ time interval following Eq. \ref{eq:tacsshots}. Then these 187 classical shadows were the data points with which we performed unsupervised learning by constructing the $187\times187$ kernel matrix using the shadow kernel in equation (\ref{eq2}) and then using diffusion maps for dimensionality reduction. Fig. \ref{fig:dynamics} shows our results. The first two eigenvectors were chosen as our diffusion space basis vectors because as Fig. \ref{fig:dynamics}(inset) shows, those are the two dominant non-trivial eigenvalues in the spectrum and the rest go to zero as the number of data points are increased. Projecting the states onto this two-dimensional diffusion space, we see that the states all fall on a curve in this hyperplane along which the value of $h_{x}$ increases monotonically, and the inflection point of the curve coincides with the quantum critical region. Therefore, it is apparent that the unsupervised learning algorithm was able to infer two phases from the data. In contrast to our ground state study, we were not able to map the diffusion coordinates to the phase-defining observables.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth,height=0.25\textheight]{dnmks.pdf}
\caption{The figure above shows phase identification from our dynamics data. The eigenvalue spectrum (inset) obtained from diffusion maps shows the two largest eigenvalues corresponding to the two dominant diffusion coordinates dc1 and dc2 (the trivial point $k=0$, is not shown). The TACS data points largely fall on a curve parameterized by $h_{x}$ in this 2D reduced diffusion space. The region near the critical point ($h_{x}=1$) coincides with the inflection point region (green) of the curve, with points on either side belonging to the two different phases.}
\label{fig:dynamics}
\end{figure}

In a microcanonical ensemble, entropy functions as the thermodynamic potential and thus characterizes the phase transitions. The quantum critical region features the interplay of thermal and quantum fluctuations leading to the entropy being maximized \cite{ Wu_2011_entropyaccumulation, wu2018crossovers}.  To calculate the $2^{nd}$ Renyi entropy, we first compute the purity using eqn \ref{eq:csshots} as follows
\begin{align}
    \gamma\left[\rho \right] &= \text{Tr}\left[ \rho^2 \right] \\
    &\approx \frac{1}{N^2} \sum_{n\neq n' = 1}^{N} \text{Tr}\left[ \sigma_1^{(n)} \sigma_1^{(n')}\right]\times \cdots \times \text{Tr}\left[ \sigma_L^{(n)} \sigma_L^{(n')}\right],
\end{align}
then $S_2= -\log_2 \gamma$.The Renyi entropy has the same limits as the von Neumann entropy- it vanishes for pure states and reaches $N$ for the maximally mixed state.

\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth,height=0.25\textheight]{Renyi_10_sites.pdf}
\caption{The $2^{nd}$ Renyi entropies for a 10-site 1DTFIM computed from TACS consisting of 100k shots.$n=\infty$ shows the Bayesian inference estimation for the thermodynamic limit and infinite TACS shots. }
\label{fig:renyi}
\end{figure}

Figure \ref{fig:renyi} presents the entropy as calculated from TACS shots on a 10 qubit system using exact diagonalization. We use Bayesian inference, as detailed in appendix \ref{app:bi} (see also Ref. \onlinecite{lukens2021bayesian}). The results show clear evidence that entropy is maximized in the quantum critical region around $h_x=1$ as shown in Fig. 2 of \cite{wu2018crossovers}. Error estimates for these values were obtained in Ref. \onlinecite{Shack}, and given by:
\begin{equation}
    N \geq \frac{4^{k+1}\gamma}{\epsilon^2\delta}
\end{equation}
Where $1-\delta$ is the probability of obtaining a good TACS dataset and $\epsilon$ is the additive error. For a $\delta = 0.33$ (67 percentile), and $N=100,000$ shots, we find an additive error for the $n=5$ entropy curve plotted in Fig. \ref{fig:renyi} at $h_x=1.0$ of $S_2/n \approx 0.5 \pm 0.24$. This is larger than the observed error shown via error bars in Fig. \ref{fig:renyi} but only within a factor of order 1.
\comment{need to say something about n=1 line}

\section{Outlook}
In this paper, we have identified an approach to studying quantum thermodynamics on a quantum computer in a way that is suitable for studying quantum materials and their phase diagrams. This approach consists of
\begin{enumerate}[label=(\roman*)]
    \item preparing an initial state which equilibrates within the coherence time of the quantum computer,
    \item time evolving this state using a quantum algorithm to a randomly chosen time point $t$ within some interval,
    \item extending shadow tomography methods to obtain a physically useful representation of the time-averaged density matrix, such as TACS used in this paper, and
    \item employing an unsupervised machine learning method to discover the phase diagram, with kernel methods such as diffusion maps showing promise.
\end{enumerate}
Our approach parallels statistical mechanics calculations on classical Hamiltonians, where a random initial state is prepared, a Metropolis Monte-Carlo algorithm is run beginning from this state, and data is collected and analyzed using traditional observables and modern machine learning methods \comment{\gaurav{refs? Is this statement important?}}. Our results, demonstrating the existence of a quantum phase transition and the ability to map out regions of the phase diagram by a careful choice of initial conditions, show promise.  Considering that we are just performing time evolution and measurements, our approach is simple compared to the existing methods \cite{Predicting_Gibbs, QITE, Partition_zeros} for studying thermodynamics on a quantum computer. The simplicity not only avoids the issues of optimizing a large number of parameters experienced by these methods but also allows us to exploit the advances in Hamiltonian dynamics. This leads to the possibility of studying thermodynamics on NISQ devices by using variational time evolution algorithms\cite{Yao_2021, Yuan_2019, VFF_2020} such as the variational fast forwarding. These algorithms are relatively robust to errors and can simulate dynamics beyond the coherence time of quantum computers. Since TACS involve only single qubit rotations and computational basis measurements, an extension of our methods to quantum simulators also seems straightforward\cite{Rydberg_atoms}. Lastly, it would be interesting to assess how this algorithm compares with finite temperature Quantum Monte-Carlo and at what system size the expected quantum advantage is reached.

\comment{
There are still a number of improvements that could be made to the approach. We have not scaled both time $T$, and system size $L$ and attempt to extrapolate to the thermodynamic limit where, according to von Neumann's quantum ergodic theorm, Eq. \ref{eq:rhoT} converges to the microcanonical ensemble.  Alternatively, we could phrase this improvement in two steps: 1) the problem of how to obtain the best estimate for $\hat\omega\equiv\lim_{T\to\infty}\hat \rho_T$, given a database of $\rho_T$ at different $T$; 2) the problem of extrapolating to infinite system size.

Another improvement would be to identify a dynamics algorithm that includes error correction to compute $\rho_T$ on a NISQ computer. One approach to this is to use variational time evolution algorithms \cite{Yao_2021, Yuan_2019, VFF_2020} such as the variational fast forwarding. These algorithms are relatively robust to errors, and are able to simulate dynamics beyond the coherence time of quantum computers. It has been shown that the error in variational fast forwarding algorithm scales at worst linearly in the fast forwarding time\cite{VFF_2020}. Post selection may further increase the accuracy of these variational algorithm by rejecting the states which violate the conservation laws. \gaurav{error corrected dynamics algorithms?}

\gaurav{Worth spectulating the use of our results for floquet dynamics. Instead of continuously sampling time, sample within each time interval. Read one of their papers and flesh out the conclusions. Express the time evolution operator in one basis $1/N U^N$ and try to get a similar formula for the time averaged density matrix. If that reduces to something like $\omega$, victory!}

\eric{Would our results be useful for studying quantum chaos? Perhaps only when we have error correction so that errors don't obscure the results.}\michael{Cool idea! Our central premise of simulating microcanonical dynamics presupposes we have error correction and the qubits are isolated from their environment.

Yes, one can study quantum chaos on a quantum computer by just time evolving a quantum system which is a quantized version of a classical system such as the kicked quantum rotor. Then how does equilibrium, studied by $\rho_T$ relate to the chaotic dynamics? Is equilibrium actually the result of quantum chaos? What about systems that don't equilibrate? To write this paragraph in the outlook, we need a reference that discussed equilibrium and quantum chaos and then add our two-cents based on what we know of the subject.}
\gaurav{Thermalization in quantum systems is a direct consequence of non-integrability (and quantum chaos in most cases). Thus, in some sense, thermalization can detect presence/absence of quantum chaos which arises from scrambling. A central tool in study of quantum scrambling is OTOC, which is a complicated quantity to calculate and our TACS does not necessarily lead to OTOC. How else does our method support the chaos community?}}

redpajama_set_name: RedPajamaArXiv