\documentclass[aps,prl,twocolumn,epsfig]{revtex4}

\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{epsfig,times}
\usepackage{dcolumn}
\usepackage{bm}
\usepackage{color}
\usepackage{epstopdf}
\usepackage{amssymb}
\usepackage{amstext}
\usepackage{latexsym}
\usepackage{hyperref}

\usepackage{amsfonts}



\newcommand{\bracket}[2]{\langle{#1}|{#2}\rangle}
\newcommand{\ketbra}[2]{|{#1}\rangle\langle{#2}|}
\newcommand{\proj}[1]{|{#1}\rangle\langle{#1}|}
\newcommand{\comment}[1]{*\textit{[{#1}]}*}
\newcommand{\trace}{{\rm Tr}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\grad}{\vec\nabla}

%%+ User-defined commands +%% Ross commands

	\newcommand{\tr}[1]{\textrm{Tr} \left[ {#1} \right]} % Faster traces

	\newcommand{\e}[1]{\textrm{e}^{ {#1}}} % Faster exponentials

	\newcommand{\s}[2]{\sigma^{#1}_{#2}} % Faster Pauli matrices

	\newcommand{\integral}[3]{{\int^{#2}_{#3} \mathrm{d}{#1} \;}} % Use as \integral{x}{1}{0} f(x) for better version of \int^1_0 dx f(x)



%+ Faster align environments
	\newcommand{\be}{\begin{equation}}
	\newcommand{\ee}{\end{equation}}
%- Faster align environments

%+ Define bra and ket vectors
	\newcommand{\ket}[1]{{\left\vert {#1} \right\rangle}}	
	\newcommand{\bra}[1]{{\left\langle {#1} \right\vert}}
%-----------------------------------------------------------------------------------------------------------------------------

\begin{document}

\title{Supporting material for\\ ``Emergent thermodynamics in a quenched quantum many-body system"}


\author{R.~Dorner$^{1,2}$, J.~Goold$^{2,3}$, C.~Cormick$^{4}$, M.~Paternostro$^{5,6}$ V.~Vedral$^{2,7,8}$}

\affiliation{$^1$Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom}
\affiliation{$^2$Clarendon Laboratory, University of Oxford, Oxford OX1 3PU, United Kingdom}
\affiliation{$^3$ Physics Department, University College Cork, Cork, Ireland}
\affiliation{$^4$ Theoretische Physik, Universit\"{a}t des Saarlandes, D-66123 Saarbr\"{u}cken, Germany}
\affiliation{$^5$Centre for Theoretical Atomic, Molecular and Optical Physics, School of Mathematics and Physics, Queen's University Belfast, BT7 1NN Belfast, United Kingdom}
\affiliation{$^6$Institut f\"ur Theoretische Physik, Albert-Einstein-Allee 11, Universit\"{a}t Ulm, D-89069 Ulm, Germany}
\affiliation{$^7$ Center for Quantum Technologies, National University of Singapore, 3 Science Drive 2, 117543 Singapore, Singapore}
\affiliation{$^8$ Department of Physics, National University of Singapore, 2 Science Drive 3, 117542 Singapore, Singapore}

\date{\today}

\maketitle

\renewcommand{\theequation}{S-\arabic{equation}}
\setcounter{equation}{0}  % reset counter
\section*{Appendix A: Diagonalization of the transverse Ising model}  % use *-form to suppress numbering
\label{dia}

The quantum Ising model in a transverse field describes a lattice of spin-1/2 particles that interact
with their nearest-neighbours via ferromagnetic coupling along the $z$-axis and with an external
field applied along the $x$-axis. For a spatially homogeneous one-dimensional lattice of $N$ spins in a uniform field, the Hamiltonian is
\begin{equation}
H = - \sum_{j=1}^N \lambda \sigma^x_j + \sigma^z_j \sigma^z_{j+1},
\label{eq:Ising1}
\end{equation}
where $\lambda$ is a dimensionless parameter that measures the strength of the external field and the Pauli spin-$1/2$ operators are defined with periodic boundary conditions $\sigma^\alpha_{N+1}=\sigma^\alpha_1$ ($\alpha=x,y,z$). Under the canonical transformation $\sigma^x_j\to\sigma^z_j, \sigma^z_j \to -\sigma^x_j$ $\forall j$ the Hamiltonian Eq.~\eqref{eq:Ising1} becomes
\begin{equation}
H = - \sum_{j=1}^N \lambda \sigma^z_j + \sigma^x_j \sigma^x_{j+1},
\label{eq:Ising2}
\end{equation}
The spin operators are mapped to a spinless fermionic operators by a Jordan-Wigner transformation, thus
\begin{equation}
\begin{aligned}
c_j={}&\frac{1}{2}\prod_{l=1}^{j-1} \s{z}{l}(\s{x}{j}+i\s{y}{j}),
\nonumber \\
c_j^\dagger={}&\frac{1}{2}\prod_{l=1}^{j-1} \s{z}{l}(\s{x}{j}-i\s{y}{j}).
\end{aligned}
\nonumber
\end{equation}
Here, the operators $c_j$ ($c_j^\dagger$)  annihilate (create) a Jordan-Wigner fermion at the $j^\textrm{th}$ lattice site and obey the usual fermionic  anti-commutation relations. This in turn allows the definition of the parity operator
\begin{equation}
\Pi := \prod_{j=1}^N \Big[1-2c_j^\dagger c_j \Big],
\nonumber
%=e^{i\pi \sum_l c^\dagger_j c_j },
\end{equation}
which measures whether the number of fermions in the chain is even ($\Pi=1$) or odd ($\Pi=-1$).
Following the Jordan-Wigner transformation the Hamiltonian Eq.~\eqref{eq:Ising2} factorises into two orthogonal parity subspaces,
\beq
H = P^+ H^+ P^+ + P^- H^- P^-.
\nonumber
\eeq
Here, $P^\pm = (1/2)\left(1\pm\Pi \right)$ are the projectors onto the even (+) and odd (-) parity subspaces and
\begin{equation}
H^{\pm}=-\lambda-\sum_{j}^{N} \left(2\lambda c_j^\dagger c_j-(c_j^\dagger c_{j+1}+c_{j+1} c_{j}+\textrm{h.c.})\right),
\label{eq:ising3}
\end{equation}
are the even and odd parity subspace contributions to the Hamiltonian. The even and odd parity subspace Hamiltonians are identical with the exception that in $H^+$ we impose the boundary condition $c_{N+1}=-c_1$ and in $H^{-}$ we impose $c_{N+1}=c_1$. Note that the parity of the chain is conserved ($\left[H^\pm,\Pi \right]=0$); the Hamiltonian in Eq.~\eqref{eq:ising3} does not mix the parity subspaces. Initialising the system in a state with zero projection onto, say, the odd subspace then restricts the dynamics to the even subspace only. For an initial Gibbs state, the system is a mixture of positive and negative parity states and both subspaces must be accounted for.
Despite this we restrict our attention to the even parity subspace only. This treatment becomes exact in the thermodynamic limit where boundary effects become negligible. This is the correct limit in which to discuss phase transitions, however a discussion of the fluctuation theorems is more suited to a finite chain. The analysis of a chain of arbitrary length in which both parity subspaces are accounted follows from a straightforward extension of what is presented here at the expense of more cumbersome expressions and provides little extra insight. Note that, in the main text, the fact that we consider the even subspace only is denoted by the summation over the set of positive parity pseudomomenta $k \in K^+$ in all relevant expressions. In the main text the `$+$' superscript that is used here to explicitly distinguish between the positive and negative parity subspaces is dropped for the sake of brevity.

The diagonalisation of the Hamiltonian is completed by application of a Fourier transformation followed by a Bogolyubov transformation, which factorizes over the spaces with different pseudomomentum $k$. For the positive parity contribution to the Hamiltonian the Fourier transformation is defined as
\beq
c_j = \frac{e^{-i\pi/4}}{\sqrt{N}} \sum_{k \in K^+} c_k e^{ikj},
\nonumber
\eeq
where the values of pseudomomentum are
\beq
K^+=\left\{k = \pm \frac{\pi}{N}(2n-1), \quad n=1, \ldots, \frac{N}{2}\right\}.
\nonumber
\eeq
After this transformation the Hamiltonian Eq.~\eqref{eq:ising3} takes the form
\beq
H^+ = \sum_{k \in K^+} 2\left(\lambda-\cos(k)\right) c^\dagger_kc_k + \sin(k) (c^\dagger_kc^\dagger_{-k} +c_{-k}c_k) -\lambda.
\nonumber
\eeq
Note that all the terms preserve pseudomomentum so that the remaining step of the diagonalization can be performed within each subspace with assigned value of $\pm k$. The last step is the Bogolyubov transformation
\begin{equation}
c_{\pm k}= \gamma_{\pm k} \cos\left(\frac{\phi_k}{2}\right) \mp \gamma^\dagger_{\mp k} \sin\left(\frac{\phi_k}{2}\right), 
%c_k &= \gamma_k \cos(\phi_k/2) - \gamma^\dagger_{-k} \sin(\phi_k/2), \\
%c_{-k} &= \gamma_{-k} \cos(\phi_k/2) + \gamma^\dagger_k \sin(\phi_k/2),
\label{Aeq:bogrelzz}
\end{equation}
where
\begin{equation}
\begin{aligned}
\cos(\phi_k) &= \frac{\lambda-\cos(k)}{\sqrt{\sin^2(k) + [\lambda-\cos(k)]^2}},\\
\sin(\phi_k) &= \frac{\sin(k)}{\sqrt{\sin^2(k) + [\lambda-\cos(k)]^2}}.
\end{aligned}
\label{Aeq:bogangles}
\end{equation}
With this, the Hamiltonian can be written in the form
\beq
H^+ = \sum_{k \in K^+} \epsilon_k \left(\gamma^\dagger_k\gamma_k-\frac{1}{2}\right),
\nonumber
\eeq
with the dispersion relation
\beq
\epsilon_k = 2\sqrt{\sin^2(k) + [\lambda-\cos(k)]^2}.
\nonumber
\eeq
Note that $\epsilon_k = \epsilon_{-k}>0$ and that the total spectrum is symmetric with respect to the zero of energy.

\section*{Appendix B: Connecting the initial and final Hamiltonians}
To evaluate the characteristic function explicitly, the eigenstates of the initial Hamiltonian $H^+(\lambda_0)$ must be written in terms of the eigenstates of the final Hamiltonian $H^+(\lambda_\tau)$. Inverting Eq.~\eqref{Aeq:bogrelzz} and its hermitian conjugate it is possible to relate the sets of pre- and post-quench Bogolyubov operators. Hence,
\begin{equation}
\begin{aligned}
\tilde\gamma_k &= \gamma_k \cos\left(\frac{\Delta_k}{2}\right) + \gamma^\dagger_{-k} \sin\left(\frac{\Delta_k}{2}\right),
\nonumber \\
\tilde\gamma_{-k} &= \gamma_{-k} \cos\left(\frac{\Delta_k}{2}\right) - \gamma^\dagger_k \sin\left(\frac{\Delta_k}{2}\right).
\end{aligned}
\end{equation}
Here $\Delta_k = \tilde\phi_k - \phi_k$ and the expressions for the pre- and post-quench Bogolyubov angles, $\phi_k$ and $\tilde\phi_k$, have the form given in Eq.\eqref{Aeq:bogangles} with $\lambda=\lambda_0$ and $\lambda_\tau$ respectively.
Using this, the vacuum states in the two representations are related by
\begin{equation}
\ket{0_k,0_{-k}} = \left( \cos\left(\frac{\Delta_k}{2}\right) + \sin\left(\frac{\Delta_k}{2}\right) \tilde{\gamma}^\dagger_k \tilde{\gamma}^\dagger_{-k}\right) \ket{\tilde0_k,\tilde0_{-k}}.
\label{eq:gsrel}
\end{equation}
The expressions for higher energy eigenstates $\ket{n_k,n_{-k}}$ are then obtained by  applying the appropriate creation operators to Eq.~\eqref{eq:gsrel}. 

\section*{Appendix C: Calculation of the average work}
The calculation for the average work done on a quenched transverse Ising model takes advantage of the factorization of $H^+(\lambda_0)$ and $H^+(\lambda_0)$ into blocks of paired pseudomomenta with labels $\pm k$. We start by writing the density matrix of the system as $\rho=\bigotimes_{k>0}\rho_{\pm k}$, so that
\beq
\label{Aeq:aws1}
\langle W \rangle=\tr{\sum_{k>0}\left( H^+_{\pm k}(\lambda_\tau)-H^+_{\pm k}(\lambda_0)\right)\bigotimes_{k'>0} \rho_{\pm k'}},
\eeq  
where $ H^+_{\pm k}(\lambda_\tau) = \epsilon_k(\lambda_\tau)({\tilde\gamma_k}^\dagger\tilde\gamma_k + {\tilde\gamma_{-k}}^\dagger\tilde\gamma_{-k} -1)$ and
$ H^+_{\pm k}(\lambda_0) =  \epsilon_k(\lambda_0)({\gamma_k}^\dagger\gamma_k + {\gamma_{-k}}^\dagger\gamma_{-k} -1)$. With a little effort, Eq.~\eqref{Aeq:aws1} can be rewritten as
\beq
\label{Aeq:aws2}
\langle W \rangle=\sum_{k>0}\trace{\left[\left( H^+_{\pm k}(\lambda_\tau)-H^+_{\pm k}(\lambda_0)\right)\prod_{k'>0}\sigma_{\pm k'}\right]},
\eeq  
where,
\begin{equation}
\sigma_{\pm k}=\sum_{n_{\pm k}=0,1}|n_{k},n_{-k}\rangle\langle n_{k},n_{-k}|\frac{e^{-\beta\epsilon_{k}(\lambda_0)(n_{k}+n_{-k}-1)}}{4\cosh^{2}(\beta \epsilon_{k}(\lambda_0)/2)}.
\nonumber
\end{equation}
Noting that $\tr{\prod_{k'>0}\sigma_{\pm k'}}=1$, Eq.\eqref{Aeq:aws2} reduces to the form $\langle W\rangle=\sum_{k>0}\langle W_k\rangle$, with 
\beq
\langle W_{k}\rangle=\trace{\left[( H^+_{\pm k}(\lambda_\tau)-H^+_{\pm k}(\lambda_0)) \sigma_{\pm k}\right]}.
\nonumber
\eeq
In order to calculate the trace we note that we need only  keep the terms of $H^+_{\pm k}(\lambda_\tau)$ that are diagonal in the 
basis of $H^+_{\pm k}(\lambda_0)$;
\beq
\left[ H^+_{\pm k}(\lambda_\tau) \right]_\textrm{diag}=\epsilon_k(\lambda_\tau) \cos(\Delta_k)({\gamma_k}^\dagger\gamma_k + {\gamma_{-k}}^\dagger\gamma_{-k} -1).
\nonumber
\eeq
With this, $\langle W_{k} \rangle$ takes form
\beq
\langle W_k \rangle=\left(\cos(\Delta_k)\epsilon_{k}(\lambda_\tau)-\epsilon_{k}(\lambda_0)\right)\tr{(n_{k}+n_{-k}-1) \sigma_{\pm k}},
\nonumber
\eeq
which leads straightforwardly to the expression for the average work in the main text.

\begin {thebibliography}{99}

\end{thebibliography}
\end{document}
