# Quasiprobability distribution of work in the quantum Ising model

Gianluca Francica, Luca Dell’Anna  Dipartimento di Fisica e Astronomia e
Sezione INFN, Università di Padova, via Marzolo 8, 35131 Padova, Italy

###### Abstract

A complete understanding of the statistics of the work done by quenching a
parameter of a quantum many-body system is still lacking in the presence of an
initial quantum coherence in the energy basis. In this case, the work can be
represented by a class of quasiprobability distributions. Here, we try to
clarify the genuinely quantum features of the process by studying the work
quasiprobability for an Ising model in a transverse field. We consider both a
global and a local quench, by focusing mainly on the thermodynamic limit. We
find that, while for a global quench there is a symmetric non-contextual
representation with a Gaussian probability distribution of work, for a local
quench we can get quantum contextuality as signaled by a negative fourth
moment of the work. Furthermore, we examine the critical features related to a
quantum phase transition and the role of the initial quantum coherence as
useful resource.

##  I Introduction

In the last years out-of-equilibrium processes generated by quenching a
parameter of a closed quantum system have been extensively investigated:
Outstanding experiments of this kind have been realized with ultra-cold atoms
[1, 2, 3], and theoretical problems concerning many-body systems have been
examined, such as thermalization and integrability [2, 4], the universality of
the dynamics across a critical point [5] and the statistics of the work done
[6]. In particular, the work statistics can be described in terms of the two-
projective measurement scheme [7] if the initial state is incoherent, i.e.,
there is no initial quantum coherence in the energy basis. In contrast, when
the initial state is not incoherent there may not be a probability
distribution for the work done as proven by a no-go theorem [8]. This is
related to the quantum contextuality as discussed in Ref. [9]. In simple
terms, the problem is similar to looking for a probability distribution in
phase space for a quantum particle in a certain quantum state. Since position
and momentum are not compatible observables, in general we get a
quasiprobability, e.g., the well-known Wigner quasiprobability [10].
Concerning the work, which in a thermally isolated quantum system is equal to
the energy change of the system, the role of position and momentum is played
by the initial and final Hamiltonian of the system. Several attempts have been
made to describe the work statistics, among these, quasiprobabilities have
been defined in terms of full-counting statistics [11] and weak values [12],
which can be viewed as particular cases of a more general quasiprobability
introduced in Ref. [13]. In general, if some fundamental conditions need to be
satisfied, the work will be represented by a class of quasiprobability
distributions [14]. Determining the possible representations of the work has a
fundamental importance: If there is some quasiprobability that is a non-
negative probability, there can be a non-contextual classical representation
of the protocol, i.e., the process can be not genuinely quantum.

Here, we focus on the statistics of the work done by quenching a parameter of
a many-body system starting from a nonequilibrium state having coherence in
the energy basis. Although some investigations on the coherence effects have
already been carried out, e.g., in Refs. [15] and [16] the full-counting
statistics and weak values quasiprobabilities have been examined, the work
statistics still remains rather uninvestigated especially in many-body
systems. Thus, after discussing the statistics of work and the quantum
contextuality in general in Sec. II, we focus on an Ising model, which we
introduce in Sec. III. Our aim is to derive some general features of global
and local quenches present in the thermodynamic limit thanks to the initial
coherence. Furthermore, we are interested to clarify what are the critical
features of the work related to a quantum phase transition: Although several
studies have been performed for initial incoherent states (e.g., on the large-
deviation of work [17, 18, 19] and the Ising model [20, 21, 22, 23, 24, 25,
26, 27, 28], just to name a few), also the initial coherence plays a role, as
found in Ref. [15], which is not entirely clear. Thus, we focus on a global
quench starting from a coherent Gibbs state in Sec. IV, where we show that,
unlike a system of finite size, in the thermodynamic limit the symmetric
quasiprobability representation of the work tends to be non-contextual, in
particular we get a Gaussian probability distribution, even if there are also
other quasiprobabilities that take negative values. In contrast, for a local
quench, since the work is not extensive, there are initial states such that
all the quasiprobabilities of the class can take negative values as signaled
by a negative fourth moment of the work (see Sec. V). Then, these processes
remain genuinely quantum also in the thermodynamic limit. Furthermore, we also
try to clarify the role of initial quantum coherence as useful resource for
the work extraction in Sec. VI, showing that, even when the protocol tends to
be non-contextual, the initial coherence still plays an active role. In the
end, we summarize and discuss further our results in Sec. VII.

##  II Work statistics

We consider a quantum quench, so that the system is initially in the state
$\rho_{0}$ and the time evolution is described by the unitary operator
$U_{t,0}$ which is generated by the time-dependent Hamiltonian
$H(\lambda_{t})$ where the control parameter $\lambda_{t}$ is changed in the
time interval $[0,\tau]$. In detail,
$U_{t,0}=\mathcal{T}e^{-i\int_{0}^{t}H(\lambda_{s})ds}$, where $\mathcal{T}$
is the time order operator and the Hamiltonian can be expressed as
$H(\lambda_{t})=\sum_{k}E_{k}(\lambda_{t})|E_{k}(\lambda_{t})\rangle\langle
E_{k}(\lambda_{t})|$ where $|E_{k}(\lambda_{t})\rangle$ is the eigenstate with
eigenvalue $E_{k}(\lambda_{t})$ at the time $t$. For brevity, we define
$E_{i}=E_{i}(\lambda_{0})$ and $E^{\prime}_{k}=E_{k}(\lambda_{\tau})$. The
average work $\langle w\rangle$ done on the system in the time interval
$[0,\tau]$ can be identified with the average energy change

$\langle
w\rangle=\text{Tr}\left\\{(H^{(H)}(\lambda_{\tau})-H(\lambda_{0}))\rho_{0}\right\\}\,,$
(1)

where, given an operator $A(t)$ we define the Heisenberg time evolved operator
$A^{(H)}(t)=U_{t,0}^{\dagger}A(t)U_{t,0}$. In general, the work performed in
the quench can be represented by a quasiprobability distribution of work. We
recall that if some Gleason-like axioms are satisfied (see Ref. [14] for
details), for the events $E\land F$ we get the quasiprobability
$v(E,F)=\text{Re}\text{Tr}\left\\{EF\rho_{0}\right\\}$, but for more than two
events, i.e., for $E\land F\land G\land\cdots$, the quasiprobability is not
fixed by the axioms. However, we can associate a quasiprobability of the form
$\text{Re}\text{Tr}\left\\{EFG\cdots\rho_{0}\right\\}$ to each of all the
possible decompositions of the form $E\land F$, $F\land G$, $G\land\cdots$. By
considering this notion of quasiprobability, if we require that (W1) the
quasiprobability distribution of work reproduces the two-projective
measurement scheme in the case of initial incoherent states (i.e., for states
$\rho_{0}$ such that $\rho_{0}=\Delta(\rho_{0})$, where we have defined the
dephasing map $\Delta(\rho_{0})=\sum_{i}|E_{i}\rangle\langle
E_{i}|\rho_{0}|E_{i}\rangle\langle E_{i}|$), (W2) the average calculated with
respect to the quasiprobability is equal to Eq. (1), and (W3) the second
moment is equal to

$\langle
w^{2}\rangle=\text{Tr}\left\\{(H^{(H)}(\lambda_{\tau})-H(\lambda_{0}))^{2}\rho_{0}\right\\}\,,$
(2)

the quasiprobability distribution of work belongs to a defined class [13, 14],
i.e., it takes the form

$\displaystyle p_{q}(w)$ $\displaystyle=$
$\displaystyle\sum_{k,j,i}\text{Re}\\{\langle
E_{i}|\rho_{0}|E_{j}\rangle\langle
E_{j}|U^{\dagger}_{\tau,0}|E^{\prime}_{k}\rangle\langle
E^{\prime}_{k}|U_{\tau,0}|E_{i}\rangle\\}$ (3)
$\displaystyle\times\delta(w-E^{\prime}_{k}+qE_{i}+(1-q)E_{j})\,,$

where $q$ is a real parameter. Our aim is to investigate this quasiprobability
for a many-body system. We can focus on the characteristic function which is
defined as $\chi_{q}(u)=\langle e^{iuw}\rangle$ and reads

$\displaystyle\chi_{q}(u)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(X_{q}(u)+X_{1-q}(u)\right)\,,$

where we have defined

$X_{q}(u)=\text{Tr}\left\\{e^{-iuqH(\lambda_{0})}\rho_{0}e^{-iu(1-q)H(\lambda_{0})}e^{iuH^{(H)}(\lambda_{\tau})}\right\\}.$
(4)

The moments of work are $\langle
w^{n}\rangle=(-i)^{n}\partial^{n}_{u}\chi_{q}(0)$, and the higher moments for
$n>2$ depend on the particular representation. In particular we get

$\langle
w^{n}\rangle=(-i)^{n}\partial_{u}^{n}\chi_{q}(0)=\frac{(-i)^{n}\partial_{u}^{n}X_{q}(0)}{2}+\frac{(-i)^{n}\partial_{u}^{n}X_{1-q}(0)}{2}\,,$
(5)

where (see Appendix A)

$\displaystyle(-i)^{n}\partial_{u}^{n}X_{q}(0)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{n}(-1)^{n-k}\binom{n}{k}\sum_{l=0}^{n-k}\binom{n-k}{l}q^{n-k-l}(1-q)^{l}$
(6)
$\displaystyle\times\text{Tr}\left\\{\rho_{0}H(\lambda_{0})^{l}(H^{(H)}(\lambda_{\tau}))^{k}H(\lambda_{0})^{n-k-l}\right\\}\,.$

We can consider the problem if there is a classical representation, i.e., if
there is a non-contextual hidden variables model which satisfies the
conditions about the reproduction of the two-projective-measurement scheme,
the average and the second moment. To introduce the concept of contextuality
at an operational level (see, e.g., Refs. [9, 29]), we consider a set of
preparations procedures $P$ and measurements procedures $M$ with outcomes $k$,
so that we will observe $k$ with probability $p(k|P,M)$. We aim to reproduce
the statistics by using a set of states $\zeta$ that are random distributed in
the set $\mathcal{Z}$ with probability $p(\zeta|P)$ every time the preparation
$P$ is performed. If, for a given $\zeta$, we get the outcome $k$ with the
probability $p(k|\zeta,M)$, we are able to reproduce the statistics if

$p(k|P,M)=\int_{\mathcal{Z}}p(\zeta|P)p(k|\zeta,M)d\zeta\,,$ (7)

and the protocol is called universally non-contextual if $p(\zeta|P)$ is a
function of the quantum state alone, i.e., $p(\zeta|P)=p(\zeta|\rho_{0})$, and
$p(k|\zeta,M)$ depends only on the positive operator-valued measurement
element $M_{k}$ associated to the corresponding outcome of the measurement
$M$, i.e., $p(k|\zeta,M)=p(k|\zeta,M_{k})$. In our case, the outcome $k$
corresponds to the work $w_{k}$, and if the protocol is non-contextual the
work distribution can be expressed as

$p(w)=\sum_{k}p(k|P,M)\delta(w-w_{k})\,,$ (8)

where $p(k|P,M)$ is given by Eq. (7) with $p(\zeta|P)=p(\zeta|\rho_{0})$ and
$p(k|\zeta,M)=p(k|\zeta,M_{k})$, so that for a negative quasiprobability of
work we cannot have a non-contextual protocol. Thus, a process that cannot be
reproduced within any non-contextual protocol will exhibit genuinely non-
classical features. If all the quasiprobabilities in the class take negative
values, the protocol is contextual, whereas if there is a quasiprobability
which is non-negative, there can be a non-contextual representation. We recall
that for an initial incoherent state $\rho_{0}=\Delta(\rho_{0})$, we get the
two-projective measurement scheme that is non-contextual [9]. In contrast, the
presence of initial quantum coherence in the energy basis can lead to a
contextual protocol. Let us investigate the effects of the initial quantum
coherence by considering a Ising model in a transverse field.

##  III Model

We consider a chain of $L$ spin 1/2 described by the Ising model in a
transverse field with Hamiltonian

$H(\lambda)=-\lambda\sum_{i=1}^{L}\sigma^{z}_{i}-\sum_{i=1}^{L}\sigma^{x}_{i}\sigma^{x}_{i+1}\,,$
(9)

where we have imposed periodic boundary conditions
$\sigma^{\alpha}_{L+1}=\sigma^{\alpha}_{1}$, and $\sigma^{\alpha}_{i}$ with
$\alpha=x,y,z$ are the Pauli matrices on the site $i$. We note that the parity
$P=\prod_{i=1}^{L}\sigma^{z}_{i}$ is a symmetry of the model, i.e., it
commutates with the Hamiltonian. The Hamiltonian can be diagonalized by
performing the Jordan-Wigner transformation

$a_{i}=\left(\prod_{j<i}\sigma^{z}_{j}\right)\sigma^{-}_{i}\,,$ (10)

where the fermionic operators $a_{i}$ satisfy the anti-commutation relations
$\\{a_{i},a_{j}^{\dagger}\\}=\delta_{i,j}$, $\\{a_{i},a_{j}\\}=0$. We get the
Hamiltonian of fermions

$\displaystyle H(\lambda)$ $\displaystyle=$
$\displaystyle-\lambda\sum_{i=1}^{L}(2a^{\dagger}_{i}a_{i}-1)-\sum_{i=1}^{L-1}(a_{i}^{\dagger}-a_{i})(a_{i+1}+a_{i+1}^{\dagger})$
(11) $\displaystyle+P(a_{L}^{\dagger}-a_{L})(a_{1}+a_{1}^{\dagger})\,,$

where the parity reads $P=e^{i\pi N}$ and
$N=\sum_{i=1}^{L}a^{\dagger}_{i}a_{i}$ is the number operator. We consider the
projector $P_{\pm}$ on the sector with parity $P=\pm 1$, then the Hamiltonian
reads

$H(\lambda)=P_{+}H_{+}(\lambda)P_{+}+P_{-}H_{-}(\lambda)P_{-}\,.$ (12)

For the sector with odd parity $P=-1$, we get the Kitaev chain

$H_{-}(\lambda)=-\lambda\sum_{i=1}^{L}(2a^{\dagger}_{i}a_{i}-1)-\sum_{i=1}^{L}(a_{i}^{\dagger}-a_{i})(a_{i+1}+a_{i+1}^{\dagger})$
(13)

with periodic boundary conditions $a_{L+1}=a_{1}$. We perform a Fourier
transform $a_{j}=1/\sqrt{L}\sum_{k}e^{-ikj}a_{k}$, where $k=2\pi n/L$ with
$n=-(L-1)/2,\ldots,(L-1)/2$ for $L$ odd and $n=-L/2+1,\ldots,L/2$ for $L$
even. Thus, the Hamiltonian reads

$H_{-}(\lambda)=\sum_{k}\Psi_{k}^{\dagger}\left[-(\lambda+\cos
k)\sigma^{z}+\sin k\sigma^{y}\right]\Psi_{k}\,,$ (14)

where $\sigma^{\alpha}$ with $\alpha=x,y,z$ are the Pauli matrices and we have
defined the Nambu spinor $\Psi_{k}=(a_{k},a_{-k}^{\dagger})^{T}$. In
particular the Hamiltonian can be written as
$H_{-}(\lambda)=\sum_{k}\Psi_{k}^{\dagger}\vec{d}_{k}\cdot\vec{\sigma}\Psi_{k}$,
which, in the diagonal form, reads

$H_{-}(\lambda)=\sum_{k}\epsilon_{k}\left(\alpha_{k}^{\dagger}\alpha_{k}-\frac{1}{2}\right)=\sum_{k}\epsilon_{k}\alpha_{k}^{\dagger}\alpha_{k}+E_{-}\,,$
(15)

where $E_{-}=-\sum_{k}\epsilon_{k}/2$. In detail we have performed a rotation
with respect to the $x$-axis with an angle $\theta_{k}$ between $\vec{d}_{k}$
and the $z$-axis, corresponding to the Bogoliubov transformation
$\alpha_{k}=\cos(\theta_{k}/2)a_{k}-i\sin(\theta_{k}/2)a^{\dagger}_{-k}$,
where $\epsilon_{k}=2||\vec{d}_{k}||$, or more explicitly,

$\epsilon_{k}=2\sqrt{(\lambda+\cos k)^{2}+\sin^{2}k}\,.$ (16)

For the sector with even parity $P=1$, we get the Hamiltonian $H_{+}(\lambda)$
which is equal to the one in Eq. (13) with antiperiodic boundary conditions
$a_{L+1}=-a_{1}$, thus the only difference is in the momenta $k$ which are
$k=2\pi(n-1/2)/L$. Of course not all eigenstates of the Hamiltonians $H_{\pm}$
are eigenstates of the Hamiltonian $H$, and their parity needs to be
discussed. Let us consider $L$ even. Thus, in the even parity sector, $k\in
K_{+}$, for each $k$ there is $-k$, and the eigenstates of the Hamiltonian are
the states

$\alpha^{\dagger}_{k_{1}}\cdots\alpha^{\dagger}_{k_{2m}}|\tilde{0}_{+}\rangle$
(17)

where $m$ is an integer, $k_{i}\in K_{+}$ and $|\tilde{0}_{+}\rangle$ is the
vacuum state of $\alpha_{k}$ with $k\in K_{+}$. In contrast, in the odd parity
sector, $k\in K_{-}$, for each $k$ there is $-k$ except for $k=0$ and $\pi$.
For $\lambda<-1$ we get $\alpha_{0}=a_{0}$ and $\alpha_{\pi}=a_{\pi}$, for
$\lambda>1$ we get $\alpha_{0}=a_{0}^{\dagger}$ and
$\alpha_{\pi}=a_{\pi}^{\dagger}$, and for $|\lambda|<1$ we get
$\alpha_{0}=a^{\dagger}_{0}$ and $\alpha_{\pi}=a_{\pi}$. Then, for
$|\lambda|>1$ the vacuum state $|\tilde{0}_{-}\rangle$ of $\alpha_{k}$ with
$k\in K_{-}$ has even parity, and the eigenstates of the Hamiltonian are the
states

$\alpha^{\dagger}_{k_{1}}\cdots\alpha^{\dagger}_{k_{2m+1}}|\tilde{0}_{-}\rangle$
(18)

with $k_{i}\in K_{-}$. Conversely, for $|\lambda|<1$ the vacuum state
$|\tilde{0}_{-}\rangle$ of $\alpha_{k}$ with $k\in K_{-}$ has odd parity since
has the fermion $a_{0}$ but not $a_{\pi}$, and the eigenstates of the
Hamiltonian are the states

$\alpha^{\dagger}_{k_{1}}\cdots\alpha^{\dagger}_{k_{2m}}|\tilde{0}_{-}\rangle$
(19)

with $k_{i}\in K_{-}$. Then, for $|\lambda|<1$ both the states
$|\tilde{0}_{+}\rangle$ and $|\tilde{0}_{-}\rangle$ are eigenstates of the
Hamiltonian with energies $E_{+}$ and $E_{-}$, so that the ground-state is
two-fold degenerate in the thermodynamic limit. Thus, at the points
$\lambda=\pm 1$ we get a second-order quantum phase transition.

##  IV Global quench

We start to focus on a sudden global quench of the transverse field $\lambda$,
i.e., $\lambda$ is suddenly changed from the value $\lambda_{0}$ to
$\lambda_{\tau}$, so that $\tau\to 0$ and $U_{\tau,0}=I$. To investigate the
role of initial quantum coherence, we focus on a coherent Gibbs state

$|\Psi_{G}(\beta)\rangle=\frac{1}{\sqrt{Z}}\sum_{j}e^{-\beta
E_{j}/2+i\varphi_{j}}|E_{j}\rangle\,,$ (20)

where $E_{j}$ are the eigenenergies of $H(\lambda_{0})$, $\varphi_{j}$ is a
phase, $Z=Z(\lambda_{0})$ and $Z(\lambda)$ is the partition function defined
as $Z(\lambda)=\text{Tr}\left\\{e^{-\beta H(\lambda)}\right\\}$. Of course,
the incoherent part of the state $|\Psi_{G}(\beta)\rangle$ is
$\Delta(|\Psi_{G}(\beta)\rangle\langle\Psi_{G}(\beta)|)=\rho_{G}(\beta)$,
where $\rho_{G}(\beta)$ is the Gibbs state $\rho_{G}(\beta)=e^{-\beta
H(\lambda_{0})}/Z$. With the aim to calculate the characteristic function for
an arbitrary size $L$, from Eq. (4) by using the relations $\sum_{s}P_{s}=I$,
$P_{s}^{2}=P_{s}$, $[P_{s},H(\lambda)]=0$ and $[P_{s},H_{\pm}(\lambda)]=0$,
where $s=\pm$, it is easy to see that

$X_{q}(u)=\sum_{s}\text{Tr}\big{\\{}e^{-iuqH_{s}(\lambda_{0})}P_{s}\rho_{0}P_{s}e^{-iu(1-q)H_{s}(\lambda_{0})}e^{iuH_{s}^{(H)}(\lambda_{\tau})}\big{\\}}\,.$
(21)

We get $P_{s}\rho_{0}P_{s}=P_{s}\rho^{s}_{0}$, where for the Gibbs state
$\rho^{s}_{0}=e^{-\beta H_{s}(\lambda_{0})}/Z$ and for the coherent Gibbs
state $\rho^{s}_{0}=|\Psi^{s}_{G}\rangle\langle\Psi^{s}_{G}|$. In particular,
we get

$|\Psi^{s}_{G}\rangle=\frac{1}{\sqrt{Z}}\otimes_{k\in
K_{s}}\left(e^{\frac{\beta\epsilon_{k}}{4}}|\tilde{0}_{k}\rangle+e^{-\frac{\beta\epsilon_{k}}{4}+i\phi_{k}}|\tilde{1}_{k}\rangle\right)\,,$
(22)

where we consider a phase such that $\phi_{-k}=\phi_{k}$, with
$|\tilde{n}_{k}\rangle=(\alpha^{\dagger}_{k})^{n_{k}}|\tilde{0}_{k}\rangle$,
where $\epsilon_{k}=\epsilon_{k}(\lambda_{0})$,
$\alpha_{k}=\alpha_{k}(\lambda_{0})$ and $|\tilde{0}_{k}\rangle$ is the vacuum
state for the fermion $\alpha_{k}$. As shown in Appendix B, we get

$X_{q}(u)=\frac{1}{2}\sum_{s}X^{s}_{q}(u)+\eta_{s}X^{\prime s}_{q}(u)\,,$ (23)

where we have defined $\eta_{s}=s\langle\tilde{0}_{s}|e^{i\pi
N}|\tilde{0}_{s}\rangle$ which is $\eta_{+}=1$ and $\eta_{-}=-1$ for
$|\lambda_{0}|>1$ and $\eta_{-}=1$ for $|\lambda_{0}|<1$, and

$X^{s}_{q}(u)=\frac{1}{Z}\prod_{k\in K_{s}:k\geq 0}X^{(k)}_{q}(u)\,.$ (24)

In detail, for $k>0$ and $k\neq\pi$, we get

$X^{(k)}_{q}(u)=X^{(k),th}_{q}(u)+X^{(k),coh}_{q}(u)\,,$ (25)

where $X^{(k),th}_{q}(u)$ is the incoherent contribution, which reads

$\displaystyle X^{(k),th}_{q}(u)$ $\displaystyle=$ $\displaystyle
2\bigg{(}\cos((u-i\beta)\epsilon_{k})\cos(u\epsilon^{\prime}_{k})+\sin((u-i\beta)\epsilon_{k})$
(26)
$\displaystyle\times\sin(u\epsilon^{\prime}_{k})\hat{d}_{k}\cdot\hat{d}^{\prime}_{k}+1\bigg{)}$

and $X^{(k),coh}_{q}(u)$ is the coherent contribution, which reads

$X^{(k),coh}_{q}(u)=-2i\sin(u\epsilon^{\prime}_{k})\sin(u(2q-1)\epsilon_{k}-2\phi_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}\,,$
(27)

where, for brevity we have defined
$\epsilon^{\prime}_{k}=\epsilon_{k}(\lambda_{\tau})$,
$\vec{d}_{k}=\vec{d}_{k}(\lambda_{0})$ and
$\vec{d}^{\prime}_{k}=\vec{d}_{k}(\lambda_{\tau})$. Furthermore, we have

$X^{\prime s}_{q}(u)=\frac{1}{Z}\prod_{k\in K_{s}:k\geq
0}X^{\prime(k)}_{q}(u)$ (28)

with

$X^{\prime(k)}_{q}(u)=X^{(k)}_{q}(u)-4\,.$ (29)

In contrast, for $k=0$ and $k=\pi$, we get

$\displaystyle X^{(0,\pi)}_{q}(u)$ $\displaystyle=$ $\displaystyle
2\cosh\left(\frac{\beta\epsilon_{0,\pi}-iu(s_{0,\pi}\epsilon^{\prime}_{0,\pi}-\epsilon_{0,\pi})}{2}\right)\,,$
(30) $\displaystyle X^{\prime(0,\pi)}_{q}(u)$ $\displaystyle=$ $\displaystyle
2\sinh\left(\frac{\beta\epsilon_{0,\pi}-iu(s_{0,\pi}\epsilon^{\prime}_{0,\pi}-\epsilon_{0,\pi})}{2}\right)\,,$
(31)

where $s_{\pi}=-1$ if $|\lambda_{0}|<1$ and $\lambda_{\tau}>1$ or
$|\lambda_{\tau}|<1$ and $\lambda_{0}>1$, otherwise $s_{\pi}=1$, and
$s_{0}=-1$ if $|\lambda_{0}|<1$ and $\lambda_{\tau}<-1$ or
$|\lambda_{\tau}|<1$ and $\lambda_{0}<-1$, otherwise $s_{0}=1$, while the
partition function is

$Z=\frac{1}{2}\sum_{s}\prod_{k\in
K_{s}}2\cosh(\beta\epsilon_{k}/2)+\eta_{s}\prod_{k\in
K_{s}}2\sinh(\beta\epsilon_{k}/2)\,.$ (32)

If the initial quantum coherence does not contribute, i.e.,
$X^{(k),coh}_{q}(u)=0$, we get $X^{(k)}_{q}(u)=X^{(k),th}_{q}(u)$ and the
characteristic function is the one of the initial Gibbs state
$\rho_{G}(\beta)$. We get $X^{(k),coh}_{q}(u)=0$ for $q=1/2$ and
$\phi_{k}=n\pi/2$, and in this case the quasiprobability is non-negative, in
particular it is equivalent to the two-projective-measurement scheme which is
non-contextual. For $q=1/2$ the initial quantum coherence contributes only for
$\phi_{k}\neq n\pi/2$ with $n$ integer. In this case the quasiprobability can
take negative values. However, in the thermodynamic limit the negativity of
the quasiprobability is always subdominant for $q=1/2$, and we get a Gaussian
probability distribution of work. To prove it, we note that in the
thermodynamic limit we get $Z=\prod_{k\in K_{+}}Z_{k}$ with
$Z_{k}=2\cosh(\beta\epsilon_{k}/2)$, then

$X_{q}(u)=\prod_{k\in K_{+}:k>0}\frac{X^{(k)}_{q}(u)}{Z_{k}^{2}}\,.$ (33)

Basically, in the thermodynamic limit the model is equivalent to the system of
fermions with Hamiltonian $H_{+}$. Thus, we can write

$X_{q}(u)=e^{Lg_{q}(u)}\,,$ (34)

where $g_{q}(u)$ is intensive, so that the work is extensive, i.e., $\langle
w^{n}\rangle\sim L^{n}$. In particular, for the initial coherent Gibbs state
under consideration, $g_{q}(u)$ explicitly reads

$g_{q}(u)=\frac{1}{2\pi}\int_{0}^{\pi}\ln\left(\frac{X^{(k)}_{q}(u)}{Z_{k}^{2}}\right)dk\,.$
(35)

Then, if Eq. (34) holds, regardless of the explicit form of the intensive
function $g_{q}(u)$, as $L\to\infty$ we can consider

$X_{q}(u)\sim
e^{L\left(\partial_{u}g_{q}(0)u+\frac{1}{2}\partial_{u}^{2}g_{q}(0)u^{2}\right)}$
(36)

since in the calculation of the Fourier transform of $X_{q}(u)$ the dominant
contribution of the integral is near $u=0$, so that we can expand $g_{q}(u)$
in Taylor series about $u=0$, and thus the neglected terms in Eq. (36) do not
contribute in the asymptotic formula of the quasiprobability $p_{q}(w)$. We
note that, although the characteristic function $\chi_{q}(u)$ depends on $q$,
the first two moments do not depend on $q$. In particular, we note that the
relative fluctuations of work scale as $\sigma_{w}/\langle w\rangle\sim
1/\sqrt{L}$, where we have defined the variance $\sigma_{w}^{2}=\langle
w^{2}\rangle-\langle w\rangle^{2}$. By noting that $\partial_{u}g_{q}(0)$ does
not depend on $q$ and
$\partial_{u}^{2}g_{1-q}(0)=\partial_{u}^{2}g^{*}_{q}(0)$, we get the
quasiprobability of work

$p_{q}(w)\sim\frac{1}{\sqrt{2\pi}}\text{Re}\left(\frac{e^{-\frac{(w-\bar{w})^{2}}{2v_{q}}}}{\sqrt{v_{q}}}\right)\,,$
(37)

where $\bar{w}=-i\partial_{u}g_{q}(0)L$ and
$v_{q}=-\partial_{u}^{2}g_{q}(0)L$. In particular the average work is $\langle
w\rangle=\bar{w}$ and the variance $\sigma_{w}^{2}$ is the real part of
$v_{q}$, i.e., $v_{q}=\sigma_{w}^{2}+ir_{q}$. As shown in Fig. 1, for $q\neq
1/2$ the asymptotic formula of the quasiprobability can take negative values
due to the presence of the imaginary part $r_{q}$.

Figure 1:  The quasiprobability of work in Eq. (37) for different values of
$q$. We put $L=50$, $\beta=1$, $\lambda_{0}=0.9$, $\lambda_{\tau}=1.1$ and
$\phi_{k}=0$.

In contrast, for $q=1/2$, we get $\chi_{1/2}(u)=X_{1/2}(u)$, from which
$\sigma_{w}^{2}=-\partial_{u}^{2}g_{1/2}(0)L$, i.e., $r_{1/2}=0$ and thus we
get the Gaussian probability distribution

$p_{1/2}(w)\sim\frac{e^{-\frac{(w-\bar{w})^{2}}{2\sigma^{2}_{w}}}}{\sqrt{2\pi}\sigma_{w}}\,.$
(38)

It is worth noting that the protocol tends to be non-contextual. To prove it,
we consider the operator $\Delta H=H^{(H)}(\lambda_{\tau})-H(\lambda_{0})$,
and the probability distribution

$p(\Delta E)=\sum_{\mu}\langle\Delta E_{\mu}|\rho_{0}|\Delta
E_{\mu}\rangle\delta(\Delta E-\Delta E_{\mu})\,,$ (39)

where $|\Delta E_{\mu}\rangle$ is the eigenstate of $\Delta H$ with eigenvalue
$\Delta E_{\mu}$. Of course $p(\Delta E)$ is non-contextual, and it is easy to
see that $p_{1/2}(w)\sim p(w)$ as $L\to\infty$. Thus, the work tends to be an
observable with respect to the non-contextual symmetric representation,
differently from finite sizes, where it is not, as was originally noted for
incoherent initial states in Ref. [7]. In particular for the quench
considered, we have $\Delta H=(\lambda_{\tau}-\lambda_{0})S_{z}$, where
$S_{z}=\sum_{j=1}^{L}\sigma^{z}_{j}$, so that the symmetric representation for
$q=1/2$ tends to be equivalent to the distribution probability of the
transverse magnetization $S_{z}$. We emphasize that for small sizes $L$ the
quasiprobability at $q=1/2$ can take negative values, but for large $L$ it is
well described by the Gaussian probability distribution in Eq. (38) (see Fig.
2).

Figure 2:  The histogram of the work distribution. We put $L=10$ in the top
panel, $L=50$ in the bottom panel, $q=1/2$, $\beta=1$, $\lambda_{\tau}=1.5$,
$\lambda_{0}=0.5$ and $\phi_{k}=\pi/4$. The red line corresponds to the
Gaussian distribution probability in Eq. (38). We note that for $L=50$ there
is still some skewness. The histograms are calculated by using the
characteristic function of Eq. (23).

We note that for an arbitrary initial state the representation for $q=1/2$ is
still non-contextual (see Appendix C). In general, the negativity of the
quasiprobability $p_{q}(w)$ can be characterized by the integral

$\mathcal{N}\equiv\int|p_{q}(w)|dw\,,$ (40)

which is equal to one if $p_{q}(w)\geq 0$. In our case,
$\mathcal{N}\sim\frac{(\sigma_{w}^{4}+r^{2}_{q})^{\frac{1}{4}}}{\sigma_{w}}$,
so that $\mathcal{N}=1$ implies that $r_{q}=0$ and thus $p_{q}(w)\geq 0$. We
note that $\mathcal{N}=1$ implies in general that $p_{q}(w)\geq 0$ (see
Appendix D). In the end we note that the effects related to the negativity of
the quasiprobability start to affect the statistics from the fourth moment,
which reads $\langle
w^{4}\rangle\sim\bar{w}^{4}+6\bar{w}^{2}\sigma^{2}_{w}+3\sigma^{4}_{w}-3r_{q}^{2}$.
In contrast the first three moments do not depend on $r_{q}$, explicitly they
read $\langle w\rangle=\bar{w}$, $\langle
w^{2}\rangle=\bar{w}^{2}+\sigma_{w}^{2}$ and $\langle
w^{3}\rangle\sim\bar{w}^{3}+3\bar{w}\sigma_{w}^{2}$. In particular, the
kurtosis is $\text{Kurt}\equiv\langle(w-\langle
w\rangle)^{4}\rangle/\sigma^{4}_{w}\sim 3-3r^{2}_{q}/\sigma_{w}^{4}$ which is
always smaller than 3 if $r_{q}\neq 0$, i.e., the distribution is more ‘flat’
than the normal one. We note that if $\bar{w}\neq 0$, since $\bar{w}\sim L$
and $\sigma^{2}_{w}\sim L$, the fourth moment is always positive. On the other
hand, for $\bar{w}=0$, the fourth moment reads $\langle w^{4}\rangle\sim
3\sigma^{4}_{w}-3r_{q}^{2}$ and becomes negative for $r_{q}>\sigma^{2}_{w}$,
so that in this regime the negativity for $q\neq 1/2$ will be strong. To
conclude our investigation concerning the global quench, we note that the
average work reads

$\bar{w}=\frac{(\lambda_{0}-\lambda_{\tau})L}{\pi}\int_{0}^{\pi}\frac{(\lambda_{0}+\cos
k)\sinh(\beta\epsilon_{k})+\sin
k\sin(2\phi_{k})}{\epsilon_{k}\cosh^{2}(\beta\epsilon_{k}/2)}dk$ (41)

and the variance reads

$\displaystyle\sigma_{w}^{2}$ $\displaystyle=$
$\displaystyle\frac{(\lambda_{0}-\lambda_{\tau})^{2}L}{\pi}\int_{0}^{\pi}\frac{1}{\cosh^{4}(\beta\epsilon_{k}/2)}\bigg{(}\cosh^{2}(\beta\epsilon_{k}/2)$
(42)
$\displaystyle\times\cosh(\beta\epsilon_{k})-\frac{2}{\epsilon_{k}^{2}}\big{(}\sin
k\sin(2\phi_{k})+(\lambda_{0}+\cos k)$
$\displaystyle\times\sinh(\beta\epsilon_{k})\big{)}^{2}\bigg{)}dk\,.$

Both $\bar{w}$ and $\sigma_{w}^{2}$ are not regular at $|\lambda_{0}|=1$ for
$\phi_{k}=\phi\neq n\pi/2$ due to the presence of a quantum phase transition
(see Fig. 3).

Figure 3:  The average work $\bar{w}$ and the variance $\sigma^{2}_{w}$ in the
function of $\lambda_{0}$ for different values of $L$. We put $\beta=1$,
$\delta\lambda=\lambda_{\tau}-\lambda_{0}=0.1$, and $\phi_{k}=\pi/4$. The
values for finite sizes $L$ are calculated by using the characteristic
function of Eq. (23).

Furthermore, concerning the negativity of the quasiprobability of work, we
have

$r_{q}=\frac{2(1-2q)(\lambda_{\tau}-\lambda_{0})L}{\pi}\int_{0}^{\pi}\frac{\sin
k\cos(2\phi_{k})}{\cosh^{2}(\beta\epsilon_{k}/2)}dk\,,$ (43)

which is regular. We deduce that the protocol admits a non-contextual
description, i.e., $r_{q}=0$, for any $q$ and $\phi_{k}=(2n+1)\pi/4$ or for
$q=1/2$. In the end, to investigate the critical features of the work which
can be related to the presence of the quantum phase transition, we introduce
the energy scale $J$ such the Hamiltonian reads

$H_{J}(\lambda)=-J\lambda\sum_{i=1}^{L}\sigma^{z}_{i}-J\sum_{i=1}^{L}\sigma^{x}_{i}\sigma^{x}_{i+1}\,.$
(44)

We focus on $\lambda_{0}\approx 1$ and we start to consider the average work
given by Eq. (41) multiplied by $J$. Then, we change variable
$k^{\prime}=\pi-k$ in the integral and we define $\kappa=k^{\prime}/a$, and
the renormalized couplings $J=c/(2a)$ and $\lambda_{0}=1-mca$. In the scaling
limit $a\to 0$, we get

$\bar{w}\sim\frac{J(\lambda_{0}-\lambda_{\tau})aL}{2\pi}\int_{0}^{\frac{\pi}{a}}\frac{\kappa\sin(2\phi_{\pi})-cm\sinh(\beta
c\omega_{\kappa})}{\omega_{\kappa}\cosh^{2}(\beta
c\omega_{\kappa}/2)}d\kappa\,,$ (45)

where $\omega_{\kappa}=\sqrt{\kappa^{2}+c^{2}m^{2}}$. We note that the
integral extended to the interval $[0,\infty)$ does not converge. Thus the
integral is not determined only by small $\kappa$, and the behavior is not
universal. Similarly, concerning the variance $\sigma_{w}^{2}$, the integral
extended to the interval $[0,\infty)$ does not converge, so that it is not
universal. The coherent contribution to the average work is defined as

$\bar{w}_{coh}=\bar{w}-\bar{w}_{th}\,,$ (46)

where $\bar{w}_{th}$ is the average work corresponding to the initial state
$\rho_{0}=\rho_{G}(\beta)$. Then, the coherent contribution is given by the
term proportional to $\sin(2\phi_{\pi})$ in Eq. (45), i.e.,

$\bar{w}_{coh}\sim\frac{J(\lambda_{0}-\lambda_{\tau})aL}{2\pi}\int_{0}^{\frac{\pi}{a}}\frac{\kappa\sin(2\phi_{\pi})}{\omega_{\kappa}\cosh^{2}(\beta
c\omega_{\kappa}/2)}d\kappa\,.$ (47)

In this case we can extend the integral to the interval $[0,\infty)$, so that
the coherent contribution $\bar{w}_{coh}$ is described by the continuum model,
in this sense it is a universal feature. From Eq. (47), by noting that

$\int_{0}^{\infty}\frac{y}{\sqrt{1+y^{2}}\cosh^{2}(x\sqrt{1+y^{2}}/2)}dy=\frac{4}{(1+e^{|x|})|x|}\,,$
(48)

the coherent contribution to the average work can be expressed as

$\bar{w}_{coh}\sim\frac{(\lambda_{0}-\lambda_{\tau})\sin(2\phi_{\pi})L}{\pi\beta}g_{FD}(\beta
mc^{2})\,,$ (49)

where we have defined the Fermi-Dirac distribution $g_{FD}(x)=1/(1+e^{|x|})$
and $mc^{2}=2J(1-\lambda_{0})$. In the end, let us consider the limit of high
temperatures $\beta\to 0$, so that we get

$\displaystyle g_{q}(u)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{\pi}\ln\frac{1}{2}\bigg{(}\cos(u\epsilon_{k})\cos(u\epsilon^{\prime}_{k})+\sin(u\epsilon_{k})\sin(u\epsilon^{\prime}_{k})$
(50)
$\displaystyle\times\hat{d}_{k}\cdot\hat{d}^{\prime}_{k}+1-i\sin(u\epsilon^{\prime}_{k})\sin(u(2q-1)\epsilon_{k}-2\phi_{k})$
$\displaystyle\times(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}\bigg{)}dk\,.$

For $\phi_{k}=\phi$, we get the closed form of the derivatives

$\displaystyle\partial_{u}g_{q}(0)$ $\displaystyle=$
$\displaystyle-\frac{i(\lambda_{\tau}-\lambda_{0})}{2\pi|\lambda_{0}|}\sin(2\phi)(1+|\lambda_{0}|-|1-|\lambda_{0}||)\,,$
(51) $\displaystyle\partial_{u}^{2}g_{q}(0)$ $\displaystyle=$
$\displaystyle-(\lambda_{\tau}-\lambda_{0})^{2}\bigg{(}1-\frac{1}{8\lambda_{0}^{2}}\left(1+\lambda_{0}^{2}-(1+|\lambda_{0}|)|1-|\lambda_{0}||\right)$
(52)
$\displaystyle\times\sin^{2}(2\phi)\bigg{)}-\frac{4i}{\pi}(\lambda_{\tau}-\lambda_{0})(1-2q)\cos(2\phi)\,,$

from which it is evident that the work statistics is not regular at
$|\lambda_{0}|=1$ for $\phi\neq n\pi/2$. Of course in this limit we can
extract the work $W_{ex}=-\langle w\rangle$, equals to

$W_{ex}=\frac{(\lambda_{\tau}-\lambda_{0})L}{2\pi|\lambda_{0}|}\sin(2\phi)(1+|\lambda_{0}|-|1-|\lambda_{0}||)\,,$
(53)

only because of the presence of the initial coherence, otherwise for an
initial Gibbs state we will get $\langle w\rangle=0$.

##  V Local quench

Things change drastically when the work is non-extensive, e.g., for a local
quench. We focus on the case of a sudden quench in the transverse field, i.e.,
the initial Hamiltonian is $H=H(\lambda_{0})$ and we perform a sudden quench
of the transverse field in a site $l$, so that the final Hamiltonian is
$H^{\prime}=H(\lambda_{0})-\epsilon\sigma^{z}_{l}$. Since we are interested
only to large sizes $L$, we describe the model with the corresponding
fermionic Hamiltonian $H_{+}$. Here we are interested to investigate how
contextuality can emerge in a local quench, thus we focus on the states
$|\Psi_{1}(\beta)\rangle$ and $|\Psi_{2}(\beta)\rangle$, which are defined as

$|\Psi_{1}(\beta)\rangle=\frac{e^{\frac{\beta}{4}\sum_{k}\epsilon_{k}}}{\sqrt{Z_{1}}}\exp\left(\sum_{k}e^{-\frac{\beta\epsilon_{k}}{2}+i\phi_{k}}\alpha_{k}^{\dagger}\right)|\tilde{0}_{+}\rangle$
(54)

and

$\displaystyle|\Psi_{2}(\beta)\rangle$ $\displaystyle=$
$\displaystyle\frac{e^{\frac{\beta}{4}\sum_{k}\epsilon_{k}}}{\sqrt{Z_{2}}}\bigg{(}1+\sum_{k}e^{-\frac{\beta\epsilon_{k}}{2}+i\phi_{k}}\alpha_{k}^{\dagger}$
(55)
$\displaystyle+\frac{1}{2}\sum_{k,k^{\prime}}s_{k,k^{\prime}}e^{-\frac{\beta(\epsilon_{k}+\epsilon_{k^{\prime}})}{2}+i(\phi_{k}+\phi_{k^{\prime}})}\alpha_{k}^{\dagger}\alpha_{k^{\prime}}^{\dagger}\bigg{)}|\tilde{0}_{+}\rangle\,,$

where $s_{k,k^{\prime}}=1$ if $k>k^{\prime}$, $s_{k,k^{\prime}}=-1$ if
$k<k^{\prime}$ and $s_{k,k}=0$, and $Z_{1}$ and $Z_{2}$ are normalization
factors such that $Z\sim Z_{2}\sim Z_{1}$ as $\beta\to\infty$. Indeed,
$|\Psi_{G}(\beta)\rangle\sim|\Psi_{2}(\beta)\rangle\sim|\Psi_{1}(\beta)\rangle$
as $\beta\to\infty$. In general, for these initial states, the function
$X_{q}(u)$ can be calculated with the help of Grassmann variables (see
Appendix E). While for the initial state $|\Psi_{1}(\beta)\rangle$, we find
that the fourth moment of work is positive, for the initial state
$|\Psi_{2}(\beta)\rangle$, we find that the fourth moment of work can be
negative for $\beta$ small enough (see Fig. 4).

Figure 4:  The fourth moment of work $\langle w^{4}\rangle$ in the function of
the local field $\epsilon$ for $q=1/2$ for the states
$|\Psi_{1}(\beta)\rangle$ (top panel) and $|\Psi_{2}(\beta)\rangle$ (bottom
panel). The curves for other values of $q\in[0,1/2]$ are not distinguishable
by eye from the one for $q=1/2$. We put $L=50$, $\lambda_{0}=1$,
$\phi_{k}=\phi_{-k}=\pi$ and $\phi_{k}=\phi_{-k}=0$ for $n$ odd and even,
where $k=2\pi(n-1/2)/L$.

This suggests that to get a contextual protocol with a negative fourth moment
we need to start from an initial state which involves at least couples of
quasiparticles, e.g., such as $|\Psi_{2}(\beta)\rangle$. This result is
corroborated by considering states like $|\Psi_{1}(\beta)\rangle$ but with
random coefficients instead of $e^{-\frac{\beta\epsilon_{k}}{2}+i\phi_{k}}$,
for which we get a non-negative fourth moment for the local quench.

##  VI Initial quantum coherence

To conclude we investigate further the role of initial coherence by focusing
on an initial state $\rho_{0}$ with a thermal incoherent part, i.e.,
$\Delta(\rho_{0})=\rho_{G}(\beta)$. In general, we have the equality (see Ref.
[13])

$\langle e^{-\beta w-C}\rangle=e^{-\beta\Delta F}\,,$ (56)

where $\Delta F=F(\lambda_{\tau})-F(\lambda_{0})$ is the change in the
equilibrium free energy, where $F(\lambda)=-\beta^{-1}\ln Z(\lambda)$, and $C$
is the random quantum coherence that has the probability distribution

$p_{c}(C)=\sum_{i,n}R_{n}\left|\langle
E_{i}|R_{n}\rangle\right|^{2}\delta(C+\ln\langle
E_{i}|\rho_{0}|E_{i}\rangle-\ln R_{n})\,,$ (57)

where we have considered the decomposition
$\rho_{0}=\sum_{n}R_{n}|R_{n}\rangle\langle R_{n}|$. In detail, the average of
$C$ is the relative entropy of coherence $\langle
C\rangle=S(\Delta(\rho_{0}))-S(\rho_{0})$, where $S(\rho)$ is the von Neumann
entropy defined as $S(\rho)=-\text{Tr}\left\\{\rho\ln\rho\right\\}$, and we
have the equality $\langle e^{-C}\rangle=1$. In particular, from Eq. (56), we
get the inequality $\langle w\rangle\geq\Delta F-\beta^{-1}\langle C\rangle$,
and we note that Eq. (56) reduces to the Jarzynski equality [30] $\langle
e^{-\beta w}\rangle=e^{-\beta\Delta F}$ when $\rho_{0}=\rho_{G}(\beta)$. From
Eq. (56) we get

$\Delta F=\beta^{-1}\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n!}\kappa_{n}(s)\,,$
(58)

where $\kappa_{n}(s)$ is the nth cumulant of $s=\beta w+C$ which of course it
can be expressed in terms of expectation values of work and coherence: The
cumulants $\kappa_{n}(C)$ of $C$ cancel in the sum due to the equality
$\langle e^{-C}\rangle=1$, and only work cumulants $\kappa_{n}(w)$ (e.g., the
variance $\sigma_{w}^{2}$) and correlation terms (e.g., the covariance
$\sigma_{w,C}=\langle wC\rangle-\langle w\rangle\langle C\rangle$) are
present. For instance, if work and coherence are uncorrelated, we get
$\kappa_{n}(s)=\beta^{n}\kappa_{n}(w)+\kappa_{n}(C)$ and so $\Delta
F=\sum_{n=1}^{\infty}(-1)^{n+1}\beta^{n-1}\kappa_{n}(w)/n!$ and the coherence
does not appear. If we consider a Gaussian probability distribution for the
random variable $s$, we get

$\Delta F=\langle w\rangle-\frac{\beta\sigma^{2}_{w}}{2}-\sigma_{w,C}\,.$ (59)

For a given free energy change $\Delta F$, from Eq. (59) we see that the
average work extracted $W_{ex}=-\langle w\rangle$ in the process increases as
the fluctuation of work becomes weak, i.e., the variance $\sigma^{2}_{w}$
decreases, and the work and coherence become strongly negative correlated,
i.e., $\sigma_{w,C}<0$, which clarifies the role of initial quantum coherence
as useful resource. However, we note that Eq. (59) cannot be exactly satisfied
for a global quench because we have to take into account also higher work
cumulants and correlations which will contribute to the series in Eq. (58) due
to large deviations. In particular, if we focus on the high temperature limit
$\beta\to 0$, Eq. (58) reduces to

$\langle w\rangle=\Delta
F+\sum_{k=1}^{\infty}\frac{i^{k+1}}{k!}\partial^{k}_{t}\partial_{u}G(0,0)\,,$
(60)

where we have defined the function $G(u,t)=\ln\langle e^{iuw+itC}\rangle$. The
derivatives are correlation terms, e.g.,
$\partial_{t}\partial_{u}G(0,0)=-\sigma_{w,C}$,
$\partial^{2}_{t}\partial_{u}G(0,0)=2i\langle
C\rangle\sigma_{w,C}-i\sigma_{w,C^{2}}$ and
$\partial^{3}_{t}\partial_{u}G(0,0)=3(2\langle C\rangle^{2}-\langle
C^{2}\rangle)\sigma_{w,C}-3\langle C\rangle\sigma_{w,C^{2}}+\sigma_{w,C^{3}}$.
For the initial state
$\rho_{0}=\eta|\Psi_{G}(0)\rangle\langle\Psi_{G}(0)|+(1-\eta)\rho_{G}(0)$, we
get the characteristic function of the coherence (see Appendix F)

$\langle
e^{itC}\rangle=D^{it}\left(\left(\eta+\frac{1-\eta}{D}\right)^{it+1}+(D-1)\left(\frac{1-\eta}{D}\right)^{it+1}\right)\,,$
(61)

where $D$ is the dimension of the Hilbert space. Furthermore, by considering

$\displaystyle\langle e^{iuw+itC}\rangle$ $\displaystyle=$
$\displaystyle\text{Tr}\\{\rho_{0}e^{it\ln\rho_{0}}e^{-iuH/2-it\ln\Delta(\rho_{0})/2}e^{iuH^{\prime}}$
(62) $\displaystyle\times e^{-iuH/2-it\ln\Delta(\rho_{0})/2}\\}\,,$

where for brevity we have defined $H=H(\lambda_{0})$ and
$H^{\prime}=H^{(H)}(\lambda_{\tau})$, we get

$-i\partial_{u}G(0,t)=\frac{\left(\eta+\frac{1-\eta}{D}\right)^{it+1}w_{1}+(D-1)\left(\frac{1-\eta}{D}\right)^{it+1}w_{2}}{\left(\eta+\frac{1-\eta}{D}\right)^{it+1}+(D-1)\left(\frac{1-\eta}{D}\right)^{it+1}}\,,$
(63)

where $w_{1}=\langle\Psi_{G}(0)|(H^{\prime}-H)|\Psi_{G}(0)\rangle$ is the
average work done starting from the coherent Gibbs state, which can be
expressed as $w_{1}=\left(\langle w\rangle-(1-\eta)\Delta F\right)/\eta$, and
$w_{2}=D\Delta F-w_{1}$, where $\Delta
F=\text{Tr}\left\\{H(\lambda_{\tau})-H(\lambda_{0})\right\\}/D$. Thus, the
terms in Eq. (60) can be obtained by calculating the derivatives of Eq. (63)
with respect to $t$. We note that for the Ising model we get $\Delta F=0$, so
that in this limit the work extracted, i.e., Eq. (53) multiplied by $\eta$,
completely comes from the correlations between work and coherence. Of course,
the same situation occurs for a cyclic change of any Hamiltonian, i.e., such
that $H(\lambda_{\tau})=H(\lambda_{0})$.

##  VII Conclusions

We investigated the effects of the initial quantum coherence in the energy
basis to the work done by quenching a transverse field of a one-dimensional
Ising model. The work can be represented by considering a class of
quasiprobability distributions. To study how the work statistics changes with
the increasing of the system size, we calculated the exact formula of the
characteristic function of work for an arbitrary size by imposing periodic
boundary conditions. Then, we focused on the thermodynamic limit, and we
showed that, for an initial coherent Gibbs state, by neglecting subdominant
terms for the symmetric value $q=1/2$ we get a Gaussian probability
distribution of work, and so a non-contextual protocol. However, for $q\neq
1/2$, the quasiprobability of work can take negative values depending on the
initial state. In contrast, for a local quench there are initial states such
that any quasiprobability representation in the class is contextual as
signaled by a negative fourth moment. We note that the quasiprobability
distribution can be measured experimentally in different ways [13, 14], also
by using a qubit (see Appendix G). In the end, beyond the fundamental purposes
of the paper, it is interesting to understand if the contextuality can be
related to some advantages from a thermodynamic point of view, however further
investigations are needed to going in this direction. In particular, although
the protocol tends to be non-contextual in the thermodynamic limit for a
global quench, the initial quantum coherence can be still a useful resource
for the work extraction in the protocol when it is correlated with the work.

### Acknowledgements

The authors acknowledge financial support from the project BIRD 2021
”Correlations, dynamics and topology in long-range quantum systems” of the
Department of Physics and Astronomy, University of Padova and from the
EuropeanUnion-NextGenerationEU within the National Center for HPC, Big Data
and Quantum Computing (Project No. CN00000013, CN1 Spoke 10 Quantum
Computing).

##  Appendix A Work moments

Let us derive a closed formula for the work moments. We define
$H=H(\lambda_{0})$ and $H^{\prime}=H^{(H)}(\lambda_{\tau})$. The nth work
moment can be calculated as

$\langle
w^{n}\rangle=(-i)^{n}\partial_{u}^{n}\chi_{q}(0)=\frac{(-i)^{n}\partial_{u}^{n}X_{q}(0)}{2}+\frac{(-i)^{n}\partial_{u}^{n}X_{1-q}(0)}{2}$
(64)

To calculate $(-i)^{n}\partial_{u}^{n}X_{q}(0)$, we note that

$X_{q}(u)=\text{Tr}\left\\{\rho_{0}(u)e^{iuH^{\prime}}\right\\}$ (65)

where we have defined

$\rho_{0}(u)=e^{-iuqH}\rho_{0}e^{-iu(1-q)H}$ (66)

Then

$(-i)^{n}\partial_{u}^{n}X_{q}(u)=\sum_{k=0}^{n}\binom{n}{k}\text{Tr}\left\\{((-i)^{n-k}\partial_{u}^{n-k}\rho_{0}(u))H^{\prime
k}e^{iuH^{\prime}}\right\\}$ (67)

where we have noted that $(-i)^{k}\partial_{u}^{k}e^{iuH^{\prime}}=H^{\prime
k}e^{iuH^{\prime}}$. It is easy to see that

$(-i)^{n}\partial_{u}^{n}\rho_{0}(u)=(-1)^{n}\sum_{k=0}^{n}\binom{n}{k}(qH)^{n-k}\rho_{0}(u)((1-q)H)^{k}$
(68)

from which

$\displaystyle(-i)^{n}\partial_{u}^{n}X_{q}(0)$ $\displaystyle=$
$\displaystyle\sum_{k=0}^{n}(-1)^{n-k}\binom{n}{k}\sum_{l=0}^{n-k}\binom{n-k}{l}q^{n-k-l}(1-q)^{l}$
(69) $\displaystyle\times\text{Tr}\left\\{H^{n-k-l}\rho_{0}H^{l}H^{\prime
k}\right\\}$

##  Appendix B Quasiprobability of work

We consider two different initial states, a Gibbs state $\rho_{G}=e^{-\beta
H(\lambda_{0})}/Z$, and a coherent Gibbs state $|\Psi_{G}\rangle$. In
particular, for $\phi_{j}=0$, the state $|\Psi^{s}_{G}\rangle$ in Eq. (22)
reads

$|\Psi^{s}_{G}\rangle=\frac{1}{\sqrt{Z}}\otimes_{k\in
K_{s}}\left(e^{\frac{\beta\epsilon_{k}}{4}}|\tilde{0}_{k}\rangle+e^{-\frac{\beta\epsilon_{k}}{4}}|\tilde{1}_{k}\rangle\right)$
(70)

It can be expressed as

$\displaystyle|\Psi^{+}_{G}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{Z}}\left(\otimes_{k>0}|\Psi_{k}\rangle\right)$
(71) $\displaystyle|\Psi^{-}_{G}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{Z}}\left(\otimes_{k>0}|\Psi_{k}\rangle\right)\otimes|\Psi_{0}\rangle\otimes|\Psi_{\pi}\rangle$
(72)

where
$|\Psi_{k}\rangle=\left(|\tilde{0}_{k}\rangle+e^{-\frac{\beta\epsilon_{k}}{2}}|\tilde{1}_{k}\rangle\right)\otimes\left(e^{\frac{\beta\epsilon_{k}}{2}}|\tilde{0}_{-k}\rangle+|\tilde{1}_{-k}\rangle\right)$.
Thus, by noting that $P_{s}=(I+se^{i\pi N})/2$, and $e^{i\pi
N}=\langle\tilde{0}_{s}|e^{i\pi N}|\tilde{0}_{s}\rangle e^{i\pi\sum_{k\in
K_{s}}\alpha^{\dagger}_{k}\alpha_{k}}$, we get

$\displaystyle X_{q}(u)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s}\text{Tr}\big{\\{}e^{-iuqH_{s}(\lambda_{0})}\rho^{s}_{0}e^{-iu(1-q)H_{s}(\lambda_{0})}e^{iuH_{s}^{(H)}(\lambda_{\tau})}\big{\\}}$
(73)
$\displaystyle+\eta_{s}\text{Tr}\big{\\{}e^{-iuqH_{s}(\lambda_{0})}e^{i\pi\sum_{k\in
K_{s}}\alpha^{\dagger}_{k}\alpha_{k}}\rho^{s}_{0}e^{-iu(1-q)H_{s}(\lambda_{0})}$
$\displaystyle\times e^{iuH_{s}^{(H)}(\lambda_{\tau})}\big{\\}}$

where we have defined $\eta_{s}=s\langle\tilde{0}_{s}|e^{i\pi
N}|\tilde{0}_{s}\rangle$. Let us focus on the first term in the sum over $s$,
which is

$X^{s}_{q}(u)=\text{Tr}\big{\\{}e^{-iuqH_{s}(\lambda_{0})}\rho^{s}_{0}e^{-iu(1-q)H_{s}(\lambda_{0})}e^{iuH_{s}^{(H)}(\lambda_{\tau})}\big{\\}}$
(74)

Then, e.g., for $s=-$, to evaluate the trace we can consider the basis formed
by the vectors
$|\\{n_{k}\\}\rangle=(\otimes_{k>0}|n_{k}n_{-k}\rangle)\otimes|n_{0}\rangle\otimes|n_{\pi}\rangle$,
with $n_{k}=0,1$, where
$|n_{k}n_{-k}\rangle=(a^{\dagger}_{k})^{n_{k}}(a^{\dagger}_{-k})^{n_{-k}}|0_{k}0_{-k}\rangle$,
where $|0_{k}\rangle$ is the vacuum state for the fermion $a_{k}$. Of course
$\\{|n_{k}n_{-k}\rangle\\}$ generates an invariant dynamically subspace, and
in this subspace the Hamiltonian $H_{s}(\lambda)$ is the matrix
$H_{k}(\lambda)$ such that

$\displaystyle H_{k}(\lambda)|0_{k}0_{-k}\rangle$ $\displaystyle=$
$\displaystyle 2(\lambda+\cos k)|0_{k}0_{-k}\rangle-2i\sin
k|1_{k}1_{-k}\rangle$ (75) $\displaystyle H_{k}(\lambda)|1_{k}1_{-k}\rangle$
$\displaystyle=$ $\displaystyle-2(\lambda+\cos k)|1_{k}1_{-k}\rangle+2i\sin
k|0_{k}0_{-k}\rangle$ (76) $\displaystyle H_{k}(\lambda)|0_{k}1_{-k}\rangle$
$\displaystyle=$ $\displaystyle 0$ (77) $\displaystyle
H_{k}(\lambda)|1_{k}0_{-k}\rangle$ $\displaystyle=$ $\displaystyle 0$ (78)

However, it is convenient to consider the initial eigenstates
$|\tilde{n}_{k}\tilde{n}_{-k}\rangle$ such that

$H_{k}(\lambda_{0})|\tilde{n}_{k}\tilde{n}_{-k}\rangle=(\epsilon_{k}n_{k}+\epsilon_{k}(n_{-k}-1))|\tilde{n}_{k}\tilde{n}_{-k}\rangle$
(79)

For our two initial states it is equal to

$X^{s}_{q}(u)=\frac{1}{Z}\prod_{k\in K_{s}:k\geq 0}X^{(k)}_{q}(u)$ (80)

For the Gibbs state, for $k>0$ and $k\neq\pi$, we have

$\displaystyle X^{(k)}_{q}(u)$ $\displaystyle=$
$\displaystyle\sum_{n_{k},n_{-k}}e^{(-iu-\beta)(\epsilon_{k}n_{k}+\epsilon_{k}(n_{-k}-1))}$
(81)
$\displaystyle\times\langle\tilde{n}_{k}\tilde{n}_{-k}|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\tilde{n}_{k}\tilde{n}_{-k}\rangle$

To evaluate $X^{(k)}_{q}(u)$, we note that

$e^{iuH_{k}(\lambda_{\tau})}=e^{-iu\epsilon^{\prime}_{k}\hat{d}^{\prime}_{k}\cdot\vec{\tau}}=(\cos(u\epsilon^{\prime}_{k})I-i\sin(u\epsilon^{\prime}_{k})\hat{d}^{\prime}_{k}\cdot\vec{\tau})\oplus
I$ (82)

where $\hat{d}^{\prime}_{k}=\hat{d}_{k}(\lambda_{\tau})$,
$\epsilon^{\prime}_{k}=\epsilon_{k}(\lambda_{\tau})$,
$\vec{\tau}=(\tau_{1},\tau_{2},\tau_{3})^{T}$, where $\tau_{i}$ are the Pauli
matrices, i.e., $\tau_{3}=|0_{k}0_{-k}\rangle\langle
0_{k}0_{-k}|-|1_{k}1_{-k}\rangle\langle 1_{k}1_{-k}|$, and so on. We have to
calculate

$\displaystyle\langle\tilde{n}_{k}\tilde{n}_{-k}|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\tilde{n}_{k}\tilde{n}_{-k}\rangle=\cos(u\epsilon^{\prime}_{k})$
$\displaystyle-i\sin(u\epsilon^{\prime}_{k})\langle\tilde{n}_{k}\tilde{n}_{-k}|U_{\tau,0}^{\dagger}\hat{d}^{\prime}_{k}\cdot\vec{\tau}U_{\tau,0}|\tilde{n}_{k}\tilde{n}_{-k}\rangle$
(83)

with $(n_{k},n_{-k})=(0,0)$ and $(n_{k},n_{-k})=(1,1)$, while
$\langle\tilde{0}_{k}\tilde{1}_{-k}|U^{\dagger}_{\tau,0}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\tilde{0}_{k}\tilde{1}_{-k}\rangle=\langle\tilde{1}_{k}\tilde{0}_{-k}|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\tilde{1}_{k}\tilde{0}_{-k}\rangle=1$.
In particular, since $\hat{d}^{\prime}_{k}\cdot\vec{\tau}$ is traceless, we
get
$\langle\tilde{0}_{k}\tilde{0}_{-k}|U^{\dagger}_{\tau,0}\hat{d}^{\prime}_{k}\cdot\vec{\tau}U_{\tau,0}|\tilde{0}_{k}\tilde{0}_{-k}\rangle+\langle\tilde{1}_{k}\tilde{1}_{-k}|U_{\tau,0}^{\dagger}\hat{d}^{\prime}_{k}\cdot\vec{\tau}U_{\tau,0}|\tilde{1}_{k}\tilde{1}_{-k}\rangle=0$,
from which we get $X^{(k)}_{q}(u)=X^{(k),th}_{q}(u)$ with

$\displaystyle
X^{(k),th}_{q}(u)=2\bigg{(}\cos((u-i\beta)\epsilon_{k})\cos(u\epsilon^{\prime}_{k})+\sin((u-i\beta)\epsilon_{k})$
$\displaystyle\times\sin(u\epsilon^{\prime}_{k})\langle\tilde{0}_{k}\tilde{0}_{-k}|U^{\dagger}_{\tau,0}\hat{d}^{\prime}_{k}\cdot\vec{\tau}U_{\tau,0}|\tilde{0}_{k}\tilde{0}_{-k}\rangle)+1\bigg{)}$
(84)

In contrast, for the coherent Gibbs state, for $k>0$ and $k\neq\pi$ we get

$X^{(k)}_{q}(u)=\langle\Psi_{k}(q-1)|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\Psi_{k}(q)\rangle$
(85)

where

$|\Psi_{k}(q)\rangle=\left(|\tilde{0}_{k}\rangle+e^{-iuq\epsilon_{k}-\frac{\beta\epsilon_{k}}{2}}|\tilde{1}_{k}\rangle\right)\otimes\left(e^{iuq\epsilon_{k}+\frac{\beta\epsilon_{k}}{2}}|\tilde{0}_{-k}\rangle+|\tilde{1}_{-k}\rangle\right)$
(86)

Thus, we get

$\displaystyle X^{(k)}_{q}(u)$ $\displaystyle=$ $\displaystyle
2\bigg{(}\cos((u-i\beta)\epsilon_{k})\cos(u\epsilon^{\prime}_{k})-\frac{i}{2}\sin(u\epsilon^{\prime}_{k})$
(87)
$\displaystyle\times\langle\tilde{\Psi}_{k}(q-1)|U_{\tau,0}^{\dagger}\hat{d}^{\prime}_{k}\cdot\vec{\tau}U_{\tau,0}|\tilde{\Psi}_{k}(q)\rangle+1\bigg{)}$

where
$|\tilde{\Psi}_{k}(q)\rangle=e^{iuq\epsilon_{k}+\beta\epsilon_{k}/2}|\tilde{0}_{k}\tilde{0}_{-k}\rangle+e^{-iuq\epsilon_{k}-\beta\epsilon_{k}/2}|\tilde{1}_{k}\tilde{1}_{-k}\rangle$.
We get

$X^{(k)}_{q}(u)=X^{(k),th}_{q}(u)+X^{(k),coh}_{q}(u)$ (88)

where the coherent contribution is

$\displaystyle X^{(k),coh}_{q}(u)$ $\displaystyle=$
$\displaystyle-2i\sin(u\epsilon^{\prime}_{k})\text{Re}\big{(}e^{-iu(2q-1)\epsilon_{k}}$
(89)
$\displaystyle\times\langle\tilde{0}_{k}\tilde{0}_{-k}|U_{\tau,0}^{\dagger}\hat{d}^{\prime}_{k}\cdot\vec{\tau}U_{\tau,0}|\tilde{1}_{k}\tilde{1}_{-k}\rangle\big{)}$

To calculate the second term in the sum over $s$ in Eq. (73), we note that

$\displaystyle e^{i\pi\sum_{k\in K_{s}}\alpha^{\dagger}_{k}\alpha_{k}}$
$\displaystyle=$ $\displaystyle(-1)^{\frac{L}{2}}e^{i\pi\sum_{k\in
K_{s}}\left(\alpha^{\dagger}_{k}\alpha_{k}-\frac{1}{2}\right)}$ (91)
$\displaystyle=$ $\displaystyle(-1)^{\frac{L}{2}}e^{iq\pi\sum_{k\in
K_{s}}\left(\alpha^{\dagger}_{k}\alpha_{k}-\frac{1}{2}\right)}$
$\displaystyle\times e^{i(1-q)\pi\sum_{k\in
K_{s}}\left(\alpha^{\dagger}_{k}\alpha_{k}-\frac{1}{2}\right)}$

then the second term is $\eta_{s}X^{\prime s}_{q}(u)$, where $X^{\prime
s}_{q}(u)$ is obtained by multiplying $X^{s}_{q}(u)$ by $(-1)^{\frac{L}{2}}$
and by performing the substitution $u\epsilon_{k}\mapsto u\epsilon_{k}-\pi$ so
that

$X^{\prime s}_{q}(u)=\frac{1}{Z}\prod_{k\in K_{s}:k\geq
0}X^{\prime(k)}_{q}(u)$ (92)

with

$X^{\prime(k)}_{q}(u)=X^{(k)}_{q}(u)-4$ (93)

for $k>0$ and $k\neq\pi$. Then, we get

$X_{q}(u)=\frac{1}{2}\sum_{s}X^{s}_{q}(u)+\eta_{s}X^{\prime s}_{q}(u)$ (94)

The partition function can be calculated as

$\displaystyle Z$ $\displaystyle=$ $\displaystyle\text{Tr}\left\\{e^{-\beta
H(\lambda_{0})}\right\\}=\sum_{s}\text{Tr}\left\\{P_{s}e^{-\beta
H_{s}(\lambda_{0})}\right\\}$ (97) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s}\text{Tr}\left\\{e^{-\beta
H_{s}(\lambda_{0})}\right\\}$
$\displaystyle+\eta_{s}\text{Tr}\left\\{e^{-\beta
H_{s}(\lambda_{0})+i\pi\sum_{k\in
K_{s}}\alpha_{k}^{\dagger}\alpha_{k}}\right\\}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s}\prod_{k\in
K_{s}}2\cosh(\beta\epsilon_{k}/2)$ $\displaystyle+\eta_{s}\prod_{k\in
K_{s}}2\sinh(\beta\epsilon_{k}/2)$

Concerning the quasiprobability distribution of work $p_{q}(w)$, it can be
calculated from the characteristic function as

$\displaystyle p_{q}(w)$ $\displaystyle=$
$\displaystyle\int\frac{e^{-iuw}}{2\pi}\chi_{q}(u)du$ (98) $\displaystyle=$
$\displaystyle\frac{1}{2}\int\frac{e^{-iuw}}{2\pi}\left(X_{q}(u)+X_{1-q}(u)\right)du$
(99)

Let us focus on the thermodynamic limit. For $|\lambda_{0}|>1$, we get
$Z\sim\prod_{k\in K_{+}}Z_{k}$ with $Z_{k}=2\cosh(\beta\epsilon_{k}/2)$, and
$X_{q}(u)\sim X^{+}_{q}(u)$. Then $X_{q}(u)$ is the product of the
characteristic functions having quasiprobability distributions

$p^{(k)}_{q}(w)=\frac{1}{Z_{k}^{2}}\int\frac{e^{-iuw}}{2\pi}X^{(k)}_{q}(u)du$
(100)

thus the quasiprobability distribution of work reads

$\displaystyle p_{q}(w)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int\left(\prod_{k>0}p^{(k)}_{q}(w_{k})+\prod_{k>0}p^{(k)}_{1-q}(w_{k})\right)$
(101)
$\displaystyle\times\delta\left(w-\sum_{k>0}w_{k}\right)\prod_{k>0}dw_{k}$

We note that the average work can be calculated as

$\langle
w\rangle=-i\partial_{u}\chi_{q}(0)=-i\sum_{k>0}\frac{1}{Z_{k}^{2}}\partial_{u}X^{(k)}_{q}(0)$
(102)

On the other hand, for $|\lambda_{0}|<1$, we get $Z\sim\prod_{k\in
K_{+}}Z_{k}+\prod_{k\in K_{+}}Z^{\prime}_{k}$ with
$Z^{\prime}_{k}=2\sinh(\beta\epsilon_{k}/2)$, and $X_{q}(u)\sim
X^{+}_{q}(u)+X^{\prime+}_{q}(u)$ from which

$X_{q}(u)=\gamma\prod_{k>0}\frac{X^{(k)}_{q}(u)}{Z_{k}^{2}}+(1-\gamma)\prod_{k>0}\frac{X^{\prime(k)}_{q}(u)}{Z^{\prime
2}_{k}}$ (103)

where $\gamma=(\prod_{k>0}Z_{k}^{2})/Z$. In the thermodynamic limit, we get

$\gamma=\frac{e^{\frac{2L}{\pi}\int_{0}^{\pi}\cosh^{2}(\beta\epsilon_{k}/2)dk}}{e^{\frac{2L}{\pi}\int_{0}^{\pi}\cosh^{2}(\beta\epsilon_{k}/2)dk}+e^{\frac{2L}{\pi}\int_{0}^{\pi}\sinh^{2}(\beta\epsilon_{k}/2)dk}}\to
1$ (104)

for a non-zero temperature, since
$\int_{0}^{\pi}\cosh^{2}(\beta\epsilon_{k}/2)dk>\int_{0}^{\pi}\sinh^{2}(\beta\epsilon_{k}/2)dk$.
In contrast, for $\beta\to\infty$, we get $Z^{\prime}_{k}\sim Z_{k}$ so that
$\gamma=1/2$, and $X^{\prime(k)}_{q}(u)\sim X^{(k)}_{q}(u)$. Then we get the
same expression of the quasiprobability of work of Eq. (101).

###  B.1 Sudden quench

Let us consider a sudden quench, i.e., the limit $\tau\to 0$, so that
$U_{\tau,0}=I$. For the Gibbs state we get $X^{(k)}_{q}(u)=X^{(k),th}_{q}(u)$
with

$\displaystyle X^{(k),th}_{q}(u)$ $\displaystyle=$ $\displaystyle
2\bigg{(}\cos((u-i\beta)\epsilon_{k})\cos(u\epsilon^{\prime}_{k})+\sin((u-i\beta)\epsilon_{k})$
(105)
$\displaystyle\times\sin(u\epsilon^{\prime}_{k})\hat{d}_{k}\cdot\hat{d}^{\prime}_{k}+1\bigg{)}$

by noting that

$\displaystyle\langle\tilde{0}_{k}\tilde{0}_{-k}|\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{0}_{k}\tilde{0}_{-k}\rangle=\text{Tr}\left\\{\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{0}_{k}\tilde{0}_{-k}\rangle\langle\tilde{0}_{k}\tilde{0}_{-k}|\right\\}$
$\displaystyle=\frac{1}{2}\text{Tr}\left\\{\hat{d}^{\prime}_{k}\cdot\vec{\tau}\hat{d}_{k}\cdot\vec{\tau}\right\\}=\hat{d}_{k}\cdot\hat{d}^{\prime}_{k}$
(106)

since
$\hat{d}_{k}\cdot\vec{\tau}=|\tilde{0}_{k}\tilde{0}_{-k}\rangle\langle\tilde{0}_{k}\tilde{0}_{-k}|-|\tilde{1}_{k}\tilde{1}_{-k}\rangle\langle\tilde{1}_{k}\tilde{1}_{-k}|$.
On the other hand, for the coherent Gibbs state we get

$X^{(k)}_{q}(u)=X^{(k),th}_{q}(u)+X^{(k),coh}_{q}(u)$ (107)

To evaluate the coherent contribution $X^{(k),coh}_{q}(u)$, we note that

$\displaystyle\langle\tilde{0}_{k}\tilde{0}_{-k}|\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{1}_{k}\tilde{1}_{-k}\rangle=\text{Tr}\left\\{\hat{d}_{k}\cdot\vec{\tau}\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{1}_{k}\tilde{1}_{-k}\rangle\langle\tilde{0}_{k}\tilde{0}_{-k}|\right\\}$
$\displaystyle=i\text{Tr}\left\\{(\hat{d}_{k}\times\hat{d}^{\prime}_{k})\cdot\vec{\tau}|\tilde{1}_{k}\tilde{1}_{-k}\rangle\langle\tilde{0}_{k}\tilde{0}_{-k}|\right\\}$
(108)

then

$\langle\tilde{0}_{k}\tilde{0}_{-k}|\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{1}_{k}\tilde{1}_{-k}\rangle=i(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}\langle\tilde{0}_{k}\tilde{0}_{-k}|\tau_{1}|\tilde{1}_{k}\tilde{1}_{-k}\rangle$
(109)

since $\hat{d}_{k}\times\hat{d}^{\prime}_{k}$ has only x-component. Since
$\langle\tilde{0}_{k}\tilde{0}_{-k}|\tau_{1}|\tilde{1}_{k}\tilde{1}_{-k}\rangle=1$,
we get

$\langle\tilde{0}_{k}\tilde{0}_{-k}|\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{1}_{k}\tilde{1}_{-k}\rangle=i(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}$
(110)

Thus we get

$X^{(k),coh}_{q}(u)=-2i\sin(u\epsilon^{\prime}_{k})\sin(u(2q-1)\epsilon_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}$
(111)

For $\phi_{j}\neq 0$, we have
$|\Psi_{k}\rangle=\left(|\tilde{0}_{k}\rangle+e^{i\phi_{k}-\frac{\beta\epsilon_{k}}{2}}|\tilde{1}_{k}\rangle\right)\otimes\left(e^{\frac{\beta\epsilon_{k}}{2}}|\tilde{0}_{-k}\rangle+e^{i\phi_{-k}}|\tilde{1}_{-k}\rangle\right)$,
with $\phi_{-k}=\phi_{k}$. Thus, by considering the corresponding state
$|\Psi_{k}(q)\rangle$, the only effect of the phase $\phi_{k}$ is the shift
$uq\epsilon_{k}\to uq\epsilon_{k}-\phi_{k}$, then we get

$X^{(k),coh}_{q}(u)=-2i\sin(u\epsilon^{\prime}_{k})\sin(u(2q-1)\epsilon_{k}-2\phi_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}$
(112)

###  B.2 Arbitrary quench

In the end, let us consider an arbitrary quench. The time evolution acts as a
rotation of the vector $\vec{d}_{k}(\lambda_{\tau})$, so that
$U_{\tau,0}^{\dagger}\hat{d}_{k}(\lambda_{\tau})\cdot\vec{\tau}U_{\tau,0}=\hat{d}^{\prime}_{k}(\lambda_{\tau})\cdot\vec{\tau}$,
where for brevity $\hat{d}^{\prime}_{k}=\hat{d}^{\prime}_{k}(\lambda_{\tau})$.
Then, $X^{(k),th}_{q}(u)$ is still given by Eq. (105) with the new vector
$\hat{d}^{\prime}_{k}$, and if $\hat{d}_{k}\times\hat{d}^{\prime}_{k}$ has
also a non-zero y-component, then

$\langle\tilde{0}_{k}\tilde{0}_{-k}|\hat{d}^{\prime}_{k}\cdot\vec{\tau}|\tilde{1}_{k}\tilde{1}_{-k}\rangle=i(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}+(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{y}$
(113)

from which the coherence contribution has a further term and reads

$\displaystyle X^{(k),coh}_{q}(u)$ $\displaystyle=$
$\displaystyle-2i\sin(u\epsilon^{\prime}_{k})\big{(}\sin(u(2q-1)\epsilon_{k}-2\phi_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}$
(114)
$\displaystyle+\cos(u(2q-1)\epsilon_{k}-2\phi_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{y}\big{)}$

###  B.3 Histogram

To determinate the quasiprobability distribution of work from the
characteristic function $\chi_{q}(u)$ we consider the intervals
$I_{n}=[w_{n}-\Delta w/2,w_{n}+\Delta w/2]$, where $w_{n}=n\Delta w$ with $n$
integer. Then, we can determinate the histogram by calculating the probability

$p_{n}=\int_{I_{n}}p_{q}(w)dw=\frac{\Delta
w}{2\pi}\int\chi_{q}(u)\text{sinc}\left(\frac{u\Delta
w}{2}\right)e^{-iuw_{n}}du$ (115)

where $\text{sinc}(x)=\sin(x)/x$. To calculate the integral we can focus on
the interval $[-2\pi K/\Delta w,2\pi K/\Delta w]$ with $K$ large enough. Of
course $p_{q}(w_{n})\approx p_{n}/\Delta w$ for $\Delta w$ small enough.

##  Appendix C Superposition of two coherent Gibbs states

For simplicity we consider the fermionic Hamiltonian $H_{+}$. We focus on the
initial state

$|\Psi\rangle=a|\Psi^{+}_{G,1}\rangle+b|\Psi^{+}_{G,2}\rangle$ (116)

where $|\Psi^{+}_{G,i}\rangle$ is the coherent Gibbs state

$|\Psi^{+}_{G,i}\rangle=\otimes_{k>0}\frac{|\Psi_{i,k}\rangle}{Z_{i,k}}$ (117)

where

$|\Psi_{i,k}\rangle=\left(|\tilde{0}_{k}\rangle+e^{i\phi_{i,k}-\frac{\beta_{i}\epsilon_{k}}{2}}|\tilde{1}_{k}\rangle\right)\otimes\left(e^{\frac{\beta_{i}\epsilon_{k}}{2}}|\tilde{0}_{-k}\rangle+e^{i\phi_{i,-k}}|\tilde{1}_{-k}\rangle\right)$
(118)

We will get

$\displaystyle X^{+}_{q}(u)$ $\displaystyle=$
$\displaystyle|a|^{2}\prod_{k>0}\frac{X^{(k)}_{q,1}(u)}{Z_{1,k}^{2}}+|b|^{2}\prod_{k>0}\frac{X^{(k)}_{q,2}(u)}{Z_{2,k}^{2}}$
(119)
$\displaystyle+ab^{*}\prod_{k>0}\frac{Y^{(k)}_{q,1}(u)}{Z_{1,k}Z_{2,k}}+a^{*}b\prod_{k>0}\frac{Y^{(k)}_{q,2}(u)}{Z_{1,k}Z_{2,k}}$

where

$\displaystyle X^{(k)}_{q,i}(u)$ $\displaystyle=$
$\displaystyle\langle\Psi_{i,k}(q-1)|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\Psi_{i,k}(q)\rangle$
(120) $\displaystyle Y^{(k)}_{q,1}(u)$ $\displaystyle=$
$\displaystyle\langle\Psi_{2,k}(q-1)|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\Psi_{1,k}(q)\rangle$
(121) $\displaystyle Y^{(k)}_{q,2}(u)$ $\displaystyle=$
$\displaystyle\langle\Psi_{1,k}(q-1)|U_{\tau,0}^{\dagger}e^{iuH_{k}(\lambda_{\tau})}U_{\tau,0}|\Psi_{2,k}(q)\rangle$
(122)

with

$|\Psi_{i,k}\rangle=\left(|\tilde{0}_{k}\rangle+e^{-iuq\epsilon_{k}+i\phi_{i,k}-\frac{\beta_{i}\epsilon_{k}}{2}}|\tilde{1}_{k}\rangle\right)\otimes\left(e^{iuq\epsilon_{k}+\frac{\beta_{i}\epsilon_{k}}{2}}|\tilde{0}_{-k}\rangle+e^{i\phi_{i,-k}}|\tilde{1}_{-k}\rangle\right)$
(123)

Then we can write

$X^{+}_{q}(u)=|a|^{2}e^{Lg_{1,q}(u)}+|b|^{2}e^{Lg_{2,q}(u)}+ab^{*}e^{Ly_{1,q}(u)}+a^{*}be^{Ly_{2,q}(u)}$
(124)

where

$\displaystyle g_{i,q}(u)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{\pi}\ln\left(\frac{X^{(k)}_{q,i}(u)}{Z_{i,k}^{2}}\right)dk$
(125) $\displaystyle y_{i,q}(u)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{0}^{\pi}\ln\left(\frac{Y^{(k)}_{q,i}(u)}{Z_{1,k}Z_{2,k}}\right)dk$
(126)

As $L\to\infty$, the Fourier transform of $e^{Ly_{i,q}(u)}$ gives

$p^{\prime}_{i,q}(w)\sim\frac{e^{Ly_{i,q}(0)}}{\sqrt{2\pi
x^{(2)}_{i,q}}}e^{-\frac{\left(w-x^{(1)}_{i,q}\right)^{2}}{2x^{(2)}_{i,q}}}$
(127)

where $x^{(1)}_{i,q}=-i\partial_{u}y_{i,q}(0)L$ and
$x^{(2)}_{i,q}=-\partial^{2}_{u}y_{i,q}(0)L$. For $q=1/2$, we get
$x^{(i)}_{2,1/2}=\left(x^{(i)}_{1,1/2}\right)^{*}$, then the quasiprobability
distribution of work reads

$\displaystyle p_{1/2}(w)$ $\displaystyle\sim$
$\displaystyle|a|^{2}\frac{e^{-\frac{(w-\bar{w}_{1})^{2}}{2v_{1,q}}}}{\sqrt{2\pi
v_{1,q}}}+|b|^{2}\frac{e^{-\frac{(w-\bar{w}_{2})^{2}}{2v_{2,q}}}}{\sqrt{2\pi
v_{2,q}}}$ (128)
$\displaystyle+2\text{Re}\left(ab^{*}p^{\prime}_{1,1/2}(w)\right)$

where $\bar{w}_{i}=-i\partial_{u}g_{i,q}(0)L$ and
$v_{i,q}=-\partial^{2}_{u}g_{i,q}(0)L$, so that $p_{1/2}(w)$ can take negative
values. However, since the real part of $y_{i,q}(0)$ is negative,
$p^{\prime}_{i,q}(w)$ tends exponentially to zero in the thermodynamic limit
and $p_{1/2}(w)$ is the convex combination of two Gaussian probability
distributions, which is positive.

###  C.1 Generalized coherent Gibbs state

For an arbitrary quench from the initial coherent Gibbs state, from Eq. (114),
we get

$\displaystyle r_{q}$ $\displaystyle=$
$\displaystyle\frac{(2q-1)L}{2\pi}\int_{0}^{\pi}\frac{\epsilon_{k}\epsilon^{\prime}_{k}}{\cosh^{2}(\beta\epsilon_{k}/2)}(\cos(2\phi_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}$
(129)
$\displaystyle+\sin(2\phi_{k})(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{y})dk$

which is zero for $q=1/2$. Let us focus on an initial state of the form

$|\Psi^{+}\rangle=\otimes_{k>0}|\Psi_{k}\rangle$ (130)

which generalizes the coherent Gibbs state $|\Psi^{+}_{G}\rangle$, where we
have defined the states

$|\Psi_{k}\rangle=\sum_{n_{k},n_{-k}}c_{n_{k}n_{-k}}|\tilde{n}_{k}\tilde{n}_{-k}\rangle$
(131)

This implies that $X_{q}(u)$ has the form in Eq. (33) with $Z_{k}=1$, where
$X^{(k)}_{q}(u)$ can be calculated from Eq. (85) with

$\displaystyle|\Psi_{k}(q)\rangle$ $\displaystyle=$ $\displaystyle
c_{00}e^{iuq\epsilon_{k}}|\tilde{0}_{k}\tilde{0}_{-k}\rangle+c_{11}e^{-iuq\epsilon_{k}}|\tilde{1}_{k}\tilde{1}_{-k}\rangle$
(132)
$\displaystyle+c_{01}|\tilde{0}_{k}\tilde{1}_{-k}\rangle+c_{10}|\tilde{1}_{k}\tilde{0}_{-k}\rangle$

Then, the representation for $q=1/2$ will be non-contextual. Let us show
explicitly that $r_{1/2}=0$. $X^{(k)}_{q}(u)$ reads

$X^{(k)}_{q}(u)=X^{(k)}_{no}(u)+\delta X^{(k)}_{q}(u)$ (133)

where $X^{(k)}_{no}(u)$ does not depend on $q$, and

$\displaystyle\delta X^{(k)}_{q}(u)$ $\displaystyle=$
$\displaystyle-2i\sin(u\epsilon^{\prime}_{k})\text{Re}\bigg{(}c_{00}^{*}c_{11}e^{-iu(2q-1)\epsilon_{k}}\big{(}i(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{x}$
(134)
$\displaystyle+(\hat{d}_{k}\times\hat{d}^{\prime}_{k})_{y}\big{)}\bigg{)}$

Then, $\partial^{2}_{u}\delta X^{(k)}_{q}(0)$ is imaginary and
$\partial^{2}_{u}\delta X^{(k)}_{q}(0)\propto(1-2q)$. Similarly, it is easy to
see that $\partial^{2}_{u}X^{(k)}_{no}(0)$ is real. Furthermore,
$\partial_{u}X^{(k)}_{q}(0)$ is imaginary, then $r_{q}$ is obtained by
calculating an integral with respect to $k$ of $\partial^{2}_{u}\delta
X^{(k)}_{q}(0)$, so that we get $r_{q}\propto(1-2q)$, which is zero for
$q=1/2$. In the end, we note that a linear combination of states of the form
in Eq. (130) will give for $q=1/2$ a convex combination of Gaussian
probability distributions, which is positive.

##  Appendix D Negativity

To prove that $\mathcal{N}=1$ implies in general that $p_{q}(w)\geq 0$, we can
proceed ad absurdum. We write $p_{q}(w)=p(w)+\delta p(w)$ where $p(w)\geq 0$,
$\int p(w)dw=1$ and $\int\delta p(w)dw=0$. If $p_{q}(w)\geq 0$ for $w\in I$
and $p_{q}(w)<0$ for $w\in I^{\prime}$, then $\delta p(w)<0$ for $w\in
I^{\prime}$ and $I=I_{+}\cup I_{-}$ such that $\delta p(w)\geq 0$ for $w\in
I_{+}$ and $\delta p(w)<0$ for $w\in I_{-}$. Then, from $\mathcal{N}=1$, we
get the condition $p(I)-p(I^{\prime})+\delta p(I_{+})+\delta p(I_{-})-\delta
p(I^{\prime})=1$, where $p(I)=\int_{I}p(w)dw$ and so on, thus we get the
system

$\begin{cases}p(I)+p(I^{\prime})=1\\\ p(I)\geq 0\\\ p(I^{\prime})\geq 0\\\
\delta p(I_{+})+\delta p(I-)+\delta p(I^{\prime})=0\\\ \delta p(I_{+})\geq
0\\\ \delta p(I_{-})<0\\\ \delta p(I^{\prime})<0\\\ p(I)-p(I^{\prime})+\delta
p(I_{+})+\delta p(I_{-})-\delta p(I^{\prime})=1\\\ \end{cases}$ (135)

which admits as solution $p(w)$ such that $0\leq p(I)<1$ and
$p(I^{\prime})=1-p(I)$ and $\delta p(w)$ such that $\delta
p(I_{+})>(1-p(I)+p(I^{\prime}))/2$, $\delta p(I_{-})=(1-2\delta
p(I_{+})-p(I)+p(I^{\prime}))/2$ and $\delta p(I^{\prime})=-\delta
p(I_{-})-\delta p(I_{+})$. Then $\delta p(I^{\prime})=-p(I^{\prime})$, so that
$p_{q}(I^{\prime})=0$, which implies that $p_{q}(w)$ is non-negative.

##  Appendix E General quadratic form in Fermi operators

We consider the initial Hamiltonian

$H=\sum_{i,j}\left(a_{i}^{\dagger}A_{ij}a_{j}+\frac{1}{2}\left(a_{i}^{\dagger}B_{ij}a_{j}^{\dagger}+H.c.\right)\right)-\frac{1}{2}\sum_{i}A_{ii}$
(136)

where $A$ and $B$ are real matrices such that $A^{T}=A$ and $B^{T}=-B$. The
Hamiltonian can be diagonalized by performing the transformation

$\alpha_{k}=\sum_{i}g_{ki}a_{i}+h_{ki}a^{\dagger}_{i}$ (137)

so that

$H=\sum_{k}\epsilon_{k}\left(\alpha^{\dagger}_{k}\alpha_{k}-\frac{1}{2}\right)$
(138)

In detail the matrices $g$ and $h$ are such that $\phi=g+h$ and $\psi=g-h$,
where $\phi$ and $\psi$ are orthogonal matrices such that
$\psi^{T}\epsilon\phi=A+B$, where $\epsilon$ is the diagonal matrix with
entries $\epsilon_{k}$. The final time-evolved Hamiltonian is $H^{\prime}$
with matrices $A^{\prime}$ and $B^{\prime}$, and will be diagonalized by
performing the transformation

$\alpha^{\prime}_{k}=\sum_{i}g^{\prime}_{ki}a_{i}+h^{\prime}_{ki}a^{\dagger}_{i}$
(139)

so that

$H^{\prime}=\sum_{k}\epsilon^{\prime}_{k}\left(\alpha^{\prime\dagger}_{k}\alpha^{\prime}_{k}-\frac{1}{2}\right)$
(140)

Let us proceed with our investigation by considering the initial state

$|\Psi_{1}\rangle=\frac{e^{\frac{\beta}{4}\sum_{k}\epsilon_{k}}}{\sqrt{Z_{1}}}\exp\left(\sum_{k}e^{-\frac{\beta\epsilon_{k}}{2}}\alpha_{k}^{\dagger}\right)|\tilde{0}\rangle$
(141)

We note that for $\beta\to\infty$ we get $Z_{1}\sim
Z=\prod_{k}2\cosh(\beta\epsilon_{k}/2)$ and
$|\Psi_{1}\rangle\sim|\Psi_{G}\rangle$. We aim to calculate

$X_{q}(u)=\langle\Psi_{1}|e^{-iu(1-q)H}e^{iuH^{\prime}}e^{-iuqH}|\Psi_{1}\rangle$
(142)

We consider the vacuum state $|\tilde{0}^{\prime}\rangle$ of the fermions
$\alpha^{\prime}_{k}$, we get the relation

$|\tilde{0}\rangle=Ke^{\frac{1}{2}\sum_{k,k^{\prime}}G_{kk^{\prime}}\alpha^{\prime\dagger}_{k}\alpha^{\prime\dagger}_{k^{\prime}}}|\tilde{0}^{\prime}\rangle$
(143)

where $G$ is solution of the equation $\tilde{g}G+\tilde{h}=0$, where
$\tilde{g}=gg^{\prime T}+hh^{\prime T}$ and $\tilde{h}=gh^{\prime
T}+hg^{\prime T}$. In particular,

$\alpha_{k}=\sum_{k^{\prime}}\tilde{g}_{kk^{\prime}}\alpha^{\prime}_{k^{\prime}}+\tilde{h}_{kk^{\prime}}\alpha^{\prime\dagger}_{k^{\prime}}$
(144)

We get

$\displaystyle
X_{q}(u)=|K|^{2}\frac{e^{(\beta+iu)\sum_{k}\epsilon_{k}/2-iu\sum_{k}\epsilon^{\prime}_{k}/2}}{Z_{1}}\langle\tilde{0}^{\prime}|\exp\bigg{(}-\frac{1}{2}\sum_{k,k^{\prime}}G_{kk^{\prime}}$
$\displaystyle\times\alpha^{\prime}_{k}\alpha^{\prime}_{k^{\prime}}\bigg{)}\exp\bigg{(}\sum_{k}u_{k}\alpha^{\prime}_{k}+v_{k}\alpha^{\prime\dagger}_{k}\bigg{)}\exp\bigg{(}\sum_{k}u^{\prime}_{k}\alpha^{\prime\dagger}_{k}+v^{\prime}_{k}\alpha^{\prime}_{k}\bigg{)}$
$\displaystyle\times\exp\bigg{(}\frac{1}{2}\sum_{k,k^{\prime}}\tilde{G}_{kk^{\prime}}\alpha^{\prime\dagger}_{k}\alpha^{\prime\dagger}_{k^{\prime}}\bigg{)}|\tilde{0}^{\prime}\rangle$
(145)

where
$\tilde{G}_{kk^{\prime}}=G_{kk^{\prime}}e^{iu(\epsilon^{\prime}_{k}+\epsilon^{\prime}_{k^{\prime}})}$
and

$\displaystyle u_{k}$ $\displaystyle=$
$\displaystyle\sum_{k^{\prime}}e^{-(\beta/2+iu(1-q))\epsilon_{k^{\prime}}}\tilde{g}_{k^{\prime}k}$
(146) $\displaystyle v_{k}$ $\displaystyle=$
$\displaystyle\sum_{k^{\prime}}e^{-(\beta/2+iu(1-q))\epsilon_{k^{\prime}}}\tilde{h}_{k^{\prime}k}$
(147) $\displaystyle u^{\prime}_{k}$ $\displaystyle=$
$\displaystyle\sum_{k^{\prime}}e^{-(\beta/2+iuq)\epsilon_{k^{\prime}}+iu\epsilon^{\prime}_{k}}\tilde{g}_{k^{\prime}k}$
(148) $\displaystyle v^{\prime}_{k}$ $\displaystyle=$
$\displaystyle\sum_{k^{\prime}}e^{-(\beta/2+iuq)\epsilon_{k^{\prime}}-iu\epsilon^{\prime}_{k}}\tilde{h}_{k^{\prime}k}$
(149)

We note that

$\displaystyle\exp\bigg{(}\sum_{k}u_{k}\alpha^{\prime}_{k}+v_{k}\alpha^{\prime\dagger}_{k}\bigg{)}\exp\bigg{(}\sum_{k}u^{\prime}_{k}\alpha^{\prime\dagger}_{k}+v^{\prime}_{k}\alpha^{\prime}_{k}\bigg{)}=1$
$\displaystyle+\sum_{k,k^{\prime}}u_{k}u^{\prime}_{k^{\prime}}\alpha^{\prime}_{k}\alpha^{\prime\dagger}_{k^{\prime}}+u_{k}v^{\prime}_{k^{\prime}}\alpha^{\prime}_{k}\alpha^{\prime}_{k^{\prime}}+v_{k}u^{\prime}_{k^{\prime}}\alpha^{\prime\dagger}_{k}\alpha^{\prime\dagger}_{k^{\prime}}-v^{\prime}_{k}v_{k^{\prime}}\alpha^{\prime}_{k}\alpha^{\prime\dagger}_{k^{\prime}}$
$\displaystyle+\sum_{k}v_{k}v^{\prime}_{k}+\cdots$ (150)

where we have omitted terms linear in the Fermi operators. Then, the overlap
in Eq. (E) can be easily calculated by using the coherent states $|\xi\rangle$
such that $\alpha^{\prime}_{k}|\xi\rangle=\xi_{k}|\xi\rangle$. By using the
identity $\int d\xi^{*}d\xi
e^{-\sum_{k}\xi^{*}_{k}\xi_{k}}|\xi\rangle\langle\xi|=1$, we get

$\displaystyle
X_{q}(u)\sim|K|^{2}\frac{e^{(\beta+iu)\sum_{k}\epsilon_{k}/2-iu\sum_{k}\epsilon^{\prime}_{k}/2}}{Z_{1}}\bigg{[}\int
d\xi^{*}d\xi$
$\displaystyle\times\exp\bigg{(}-\frac{1}{2}\sum_{k,k^{\prime}}G_{kk^{\prime}}\xi_{k}\xi_{k^{\prime}}+\sum_{k,k^{\prime}}(u_{k}u^{\prime}_{k^{\prime}}\xi_{k}\xi^{*}_{k^{\prime}}+u_{k}v^{\prime}_{k^{\prime}}\xi_{k}\xi_{k^{\prime}}$
$\displaystyle+v_{k}u^{\prime}_{k^{\prime}}\xi^{*}_{k}\xi^{*}_{k^{\prime}}-v^{\prime}_{k}v_{k^{\prime}}\xi_{k}\xi^{*}_{k^{\prime}})-\sum_{k}\xi^{*}_{k}\xi_{k}+\frac{1}{2}\sum_{k,k^{\prime}}\tilde{G}_{kk^{\prime}}\xi^{*}_{k}\xi^{*}_{k^{\prime}}\bigg{)}$
$\displaystyle+\sum_{k}v_{k}v^{\prime}_{k}\int
d\xi^{*}d\xi\exp\bigg{(}-\frac{1}{2}\sum_{k,k^{\prime}}G_{kk^{\prime}}-\sum_{k}\xi^{*}_{k}\xi_{k}$
$\displaystyle+\frac{1}{2}\sum_{k,k^{\prime}}\tilde{G}_{kk^{\prime}}\xi^{*}_{k}\xi^{*}_{k^{\prime}}\bigg{)}\bigg{]}$
(151)

By performing the integral, we get

$X_{q}(u)\sim
Ce^{iu\sum_{k}(\epsilon_{k}-\epsilon^{\prime}_{k})/2}\bigg{(}\sqrt{\det(\Gamma(u))}+\sqrt{\det(\Gamma_{0}(u))}\sum_{k}v_{k}v^{\prime}_{k}\bigg{)}$
(152)

where

$\Gamma_{0}(u)=\left(\begin{array}[]{cc}G&-I\\\ I&-\tilde{G}\\\
\end{array}\right)$ (153)

and

$\Gamma(u)=\left(\begin{array}[]{cc}G-M_{1}&-I-M_{2}\\\
I+M_{2}^{T}&-\tilde{G}-M_{3}\\\ \end{array}\right)=\Gamma_{0}(u)+M(u)$ (154)

where
$M_{1,kk^{\prime}}=u_{k}v^{\prime}_{k^{\prime}}-u_{k^{\prime}}v^{\prime}_{k}$,
$M_{2,kk^{\prime}}=u_{k}u^{\prime}_{k^{\prime}}-v^{\prime}_{k}v_{k^{\prime}}$
and
$M_{3,kk^{\prime}}=v_{k}u^{\prime}_{k^{\prime}}-v_{k^{\prime}}u^{\prime}_{k}$.
The constant $C$ can be determined by requiring that $X_{q}(0)=1$. The exact
expression of $X_{q}(u)$ can be obtained by expanding $\sqrt{\det(\Gamma(u))}$
at the first order in $M(u)$, i.e.,

$\displaystyle X_{q}(u)$ $\displaystyle=$ $\displaystyle
Ce^{iu\sum_{k}(\epsilon_{k}-\epsilon^{\prime}_{k})/2}\sqrt{\det(\Gamma_{0}(u))}\bigg{(}1+\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}(u)M(u)\right\\}$
(155) $\displaystyle+\sum_{k}v_{k}v^{\prime}_{k}\bigg{)}$

Concerning the coherent Gibbs state, for low temperatures $\beta\to\infty$ we
get

$|\Psi_{G}\rangle\sim\frac{e^{\frac{\beta}{4}\sum_{k}\epsilon_{k}}}{\sqrt{Z}}\left(1+\sum_{k}e^{-\frac{\beta\epsilon_{k}}{2}}\alpha_{k}^{\dagger}+\sum_{k>k^{\prime}}e^{-\frac{\beta(\epsilon_{k}+\epsilon_{k^{\prime}})}{2}}\alpha_{k}^{\dagger}\alpha_{k^{\prime}}^{\dagger}\right)|\tilde{0}\rangle$
(156)

We define

$\displaystyle u_{kq}$ $\displaystyle=$ $\displaystyle
e^{-(\beta/2+iu(1-q))\epsilon_{k}}\tilde{g}_{kq}$ (157) $\displaystyle v_{kq}$
$\displaystyle=$ $\displaystyle
e^{-(\beta/2+iu(1-q))\epsilon_{k}}\tilde{h}_{kq}$ (158) $\displaystyle
u^{\prime}_{kq}$ $\displaystyle=$ $\displaystyle
e^{-(\beta/2+iuq)\epsilon_{k}+iu\epsilon^{\prime}_{q}}\tilde{g}_{kq}$ (159)
$\displaystyle v^{\prime}_{kq}$ $\displaystyle=$ $\displaystyle
e^{-(\beta/2+iuq)\epsilon_{k}-iu\epsilon^{\prime}_{q}}\tilde{h}_{kq}$ (160)

so that $u_{k}=\sum_{k^{\prime}}u_{k^{\prime}k}$ and so on, then the matrices
$V_{1}$, $V_{2}$, $V_{3}$, $V^{\prime}_{1}$, $V^{\prime}_{2}$ and
$V^{\prime}_{3}$ with elements
$V_{1,qq^{\prime}}=\sum_{k,k^{\prime}}s_{k,k^{\prime}}u_{kq}u_{k^{\prime}q^{\prime}}$,
$V_{2,qq^{\prime}}=\sum_{k,k^{\prime}}s_{k,k^{\prime}}u_{kq}v_{k^{\prime}q^{\prime}}$,
$V_{3,qq^{\prime}}=\sum_{k,k^{\prime}}s_{k,k^{\prime}}v_{kq}v_{k^{\prime}q^{\prime}}$,
$V^{\prime}_{1,qq^{\prime}}=\sum_{k,k^{\prime}}s_{k,k^{\prime}}v^{\prime}_{kq}v^{\prime}_{k^{\prime}q^{\prime}}$,
$V^{\prime}_{2,qq^{\prime}}=\sum_{k,k^{\prime}}s_{k,k^{\prime}}v^{\prime}_{kq}u^{\prime}_{k^{\prime}q^{\prime}}$,
$V^{\prime}_{3,qq^{\prime}}=\sum_{k,k^{\prime}}s_{k,k^{\prime}}u^{\prime}_{kq}u^{\prime}_{k^{\prime}q^{\prime}}$,
where $s_{k,k^{\prime}}=1$ if $k>k^{\prime}$, $s_{k,k^{\prime}}=-1$ if
$k<k^{\prime}$ and $s_{k,k}=0$. Thus, by proceeding similarly, we get at the
second order

$\displaystyle X_{q}(u)$ $\displaystyle\sim$ $\displaystyle
Ce^{iu\sum_{k}(\epsilon_{k}-\epsilon^{\prime}_{k})/2}\sqrt{\det(\Gamma_{0}(u))}\bigg{(}1+\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}(u)M(u)\right\\}$
(161)
$\displaystyle+\sum_{k}v_{k}v^{\prime}_{k}+\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}(u)(V(u)-V^{\prime}(u))\right\\}$
$\displaystyle+\frac{1}{2}\text{Tr}\left\\{V_{2}-V^{\prime}_{2}\right\\}\bigg{)}$

where we have defined the matrices

$V(u)=\left(\begin{array}[]{cc}V_{1}&V_{2}\\\ -V_{2}^{T}&V_{3}\\\
\end{array}\right)\,,\quad
V^{\prime}(u)=\left(\begin{array}[]{cc}V^{\prime}_{1}&V^{\prime}_{2}\\\
-V^{\prime T}_{2}&V^{\prime}_{3}\\\ \end{array}\right)$ (162)

We note that for an initial state that is the ground-state of $H$, we get the
characteristic function

$\chi^{(0)}(u)=e^{iu\sum_{k}(\epsilon_{k}-\epsilon^{\prime}_{k})/2}\sqrt{\frac{\det(\Gamma_{0}(u))}{\det(\Gamma_{0}(0))}}$
(163)

which is obtained from $X_{q}(u)$ in the limit $\beta\to\infty$.
Alternatively, by considering $\theta^{T}=(\xi^{T},{\xi^{*}}^{T})$, Eq. (155)
can be derived with the help of the identity

$\int
d\theta\theta_{i}\theta_{j}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}=-\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}X_{ij}\right\\}\sqrt{\det(\Gamma_{0})}$
(164)

where $X_{ij}=|i\rangle\langle j|-|j\rangle\langle i|$, and $|i\rangle$ is the
unit vector with only the i-th component which is nonzero. Actually
$\sqrt{\det(\Gamma_{0})}$ is the Pfaffian of $\Gamma_{0}$. To prove it, we
note that

$\displaystyle\int
d\theta\theta_{i}\theta_{j}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}$
$\displaystyle=$ $\displaystyle\frac{1}{\epsilon}\bigg{(}\int
d\theta(1+\epsilon\theta_{i}\theta_{j})e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}$
(165) $\displaystyle-\int d\theta
e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}\bigg{)}$

The second integral is $\int d\theta
e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}=\sqrt{\det(\Gamma_{0})}$. By
considering the limit $\epsilon\to 0$, we get

$\int
d\theta\theta_{i}\theta_{j}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}\sim\frac{1}{\epsilon}\left(\int
d\theta
e^{\frac{\epsilon}{2}(\theta_{i}\theta_{j}-\theta_{j}\theta_{i})}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}-\sqrt{\det(\Gamma_{0})}\right)$
(166)

then

$\int
d\theta\theta_{i}\theta_{j}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}\sim\frac{1}{\epsilon}\left(\sqrt{\det(\Gamma_{0}-\epsilon
X_{ij})}-\sqrt{\det(\Gamma_{0})}\right)$ (167)

by evaluating the limit $\epsilon\to 0$, we get Eq. (164). Similarly, we have
the identity

$\displaystyle\int
d\theta\theta_{i}\theta_{j}\theta_{k}\theta_{l}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}=-\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}X_{ij}\Gamma_{0}^{-1}X_{kl}\right\\}\sqrt{\det(\Gamma_{0})}$
$\displaystyle+\frac{1}{4}\text{Tr}\left\\{\Gamma_{0}^{-1}X_{ij}\right\\}\text{Tr}\left\\{\Gamma_{0}^{-1}X_{kl}\right\\}\sqrt{\det(\Gamma_{0})}$
(168)

To prove it, we consider that

$\displaystyle\int
d\theta\theta_{i}\theta_{j}\theta_{k}\theta_{l}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}$
$\displaystyle=$ $\displaystyle\frac{1}{\epsilon}\bigg{(}\int
d\theta\theta_{i}\theta_{j}(1+\epsilon\theta_{k}\theta_{l})e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}$
(169) $\displaystyle-\int
d\theta\theta_{i}\theta_{j}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}\bigg{)}$

which, in the limit $\epsilon\to 0$ can be evaluated with the help of the
identity in Eq. (164). We get

$\displaystyle\int
d\theta\theta_{i}\theta_{j}\theta_{k}\theta_{l}e^{-\frac{1}{2}\theta^{T}\Gamma_{0}\theta}\sim\frac{1}{2\epsilon}\bigg{(}\text{Tr}\left\\{\Gamma_{0}^{-1}X_{ij}\right\\}\sqrt{\det(\Gamma_{0})}$
$\displaystyle-\text{Tr}\left\\{(\Gamma_{0}-\epsilon
X_{kl})^{-1}X_{ij}\right\\}\sqrt{\det(\Gamma_{0}-\epsilon X_{kl})}\bigg{)}$
(170)

by evaluating the limit $\epsilon\to 0$, we get Eq. (E). In the end, we
consider the initial state in Eq. (156), which is

$|\Psi_{2}\rangle=\frac{e^{\frac{\beta}{4}\sum_{k}\epsilon_{k}}}{\sqrt{Z_{2}}}\left(1+\sum_{k}e^{-\frac{\beta\epsilon_{k}}{2}}\alpha_{k}^{\dagger}+\frac{1}{2}\sum_{k,k^{\prime}}s_{k,k^{\prime}}e^{-\frac{\beta(\epsilon_{k}+\epsilon_{k^{\prime}})}{2}}\alpha_{k}^{\dagger}\alpha_{k^{\prime}}^{\dagger}\right)|\tilde{0}\rangle$
(171)

By using the identities in Eqs. (164) and (E), we get

$\displaystyle X_{q}(u)$ $\displaystyle=$ $\displaystyle
Ce^{iu\sum_{k}(\epsilon_{k}-\epsilon^{\prime}_{k})/2}\sqrt{\det(\Gamma_{0}(u))}\bigg{(}1+\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}(u)M(u)\right\\}$
(172)
$\displaystyle+\sum_{k}v_{k}v^{\prime}_{k}+\frac{1}{2}\text{Tr}\left\\{\Gamma_{0}^{-1}(u)(V(u)-V^{\prime}(u))\right\\}$
$\displaystyle+\frac{1}{2}\text{Tr}\left\\{V_{2}-V^{\prime}_{2}\right\\}-\frac{1}{4}\text{Tr}\left\\{V_{2}\right\\}\text{Tr}\left\\{V^{\prime}_{2}\right\\}$
$\displaystyle-\frac{1}{4}\text{Tr}\left\\{V_{2}\right\\}\text{Tr}\left\\{\Gamma^{-1}_{0}(u)V^{\prime}(u)\right\\}$
$\displaystyle-\frac{1}{4}\text{Tr}\left\\{V^{\prime}_{2}\right\\}\text{Tr}\left\\{\Gamma^{-1}_{0}(u)V(u)\right\\}-\frac{1}{2}\text{Tr}\left\\{V_{3}V^{\prime}_{1}\right\\}$
$\displaystyle+\frac{1}{2}\text{Tr}\left\\{\Gamma^{-1}_{0}(u)V^{\prime\prime}(u)\right\\}+\frac{1}{2}\text{Tr}\left\\{\Gamma^{-1}_{0}(u)V(u)\Gamma^{-1}_{0}(u)V^{\prime}(u)\right\\}$
$\displaystyle-\frac{1}{4}\text{Tr}\left\\{\Gamma^{-1}_{0}(u)V(u)\right\\}\text{Tr}\left\\{\Gamma^{-1}_{0}(u)V^{\prime}(u)\right\\}\bigg{)}$

where we have defined

$V^{\prime\prime}(u)=\left(\begin{array}[]{cc}V_{2}V^{\prime}_{1}+V^{\prime}_{1}V^{T}_{2}&V_{2}V^{\prime}_{2}-V^{\prime}_{1}V_{3}\\\
V_{3}V^{\prime}_{1}-V^{\prime T}_{2}V_{2}^{T}&V_{3}V^{\prime}_{2}+V^{\prime
T}_{2}V_{3}\\\ \end{array}\right)$ (173)

If we introduce a relative phase $\phi_{k}$, we have to multiply $u_{kq}$ and
$v_{kq}$ by $e^{-i\phi_{k}}$ and $u^{\prime}_{kq}$ and $v^{\prime}_{kq}$ by
$e^{i\phi_{k}}$.

If $A$ and $B$ are complex matrices, we get $g$ and $h$ complex. In this case
we have same formulas, with $\tilde{g}=gg^{\prime\dagger}+hh^{\prime T}$ and
$\tilde{h}=gh^{\prime T}+hg^{\prime\dagger}$, and in $\Gamma_{0}$ in Eq.
(153), we have $G^{*}$ instead of $G$, and in $u_{k}$, $v_{k}$, $u_{kq}$ and
$v_{kq}$ we have $\tilde{g}^{*}$ and $\tilde{h}^{*}$ instead of $\tilde{g}$
and $\tilde{h}$.

##  Appendix F Initial quantum coherence

We consider the initial state
$\rho_{0}=\eta|\Psi_{G}(0)\rangle\langle\Psi_{G}(0)|+(1-\eta)\rho_{G}(0)$, we
get

$\langle
e^{itC}\rangle=\text{Tr}\left\\{\rho_{0}e^{i\ln\rho_{0}}e^{i\ln\rho_{G}(0)t}\right\\}=\text{Tr}\left\\{\rho_{0}^{1+it}\right\\}D^{it}$
(174)

since $\rho_{G}(0)=I/D$ is the completely mixed state. Then, since the
eigenvalues of $\rho_{0}$ are $\eta+(1-\eta)/D$ and $(1-\eta)/D$ which is
$D-1$ fold degenerate, by evaluating the trace we get

$\langle
e^{itC}\rangle=D^{it}\left((\eta+(1-\eta)/D)^{1+it}+(D-1)((1-\eta)/D)^{1+it}\right)$
(175)

which is Eq. (61). Concerning $\langle e^{iuw+itC}\rangle$ can be easily
derived from the joint quasiprobability distribution of the work and coherence
given in Ref. [13]. By doing a symmetric choice of the parameters
$q=q^{\prime}=1/2$ we get Eq. (62), from which

$-i\partial_{u}G(0,t)=-i\partial_{u}\ln\langle
e^{itC+iuw}\rangle|_{u=0}=\frac{\text{Tr}\left\\{\rho_{0}e^{i\ln\rho_{0}t}(H^{\prime}-H)\right\\}}{\text{Tr}\left\\{\rho_{0}e^{i\ln\rho_{0}t}\right\\}}$
(176)

and by proceeding similarly we get Eq. (63).

##  Appendix G Measuring the characteristic function

The characteristic function can be measured as observed in Ref. [13]. Here we
note the detector can be a qubit in the initial state $\rho_{D}(t_{i})$ with
Hamiltonian $H_{D}=\omega|e\rangle\langle e|$. We consider the interactions
with the system described by $H_{I}=-\delta_{e}|e\rangle\langle
e|-\delta_{g}|g\rangle\langle g|$ and
$H^{\prime}_{I}=-\delta^{\prime}_{e}|e\rangle\langle
e|-\delta^{\prime}_{g}|g\rangle\langle g|$, where $|g\rangle$ is the ground-
state of the qubit and $|e\rangle$ is the excited state. The total system is
in the initial state $\rho_{D}(t_{i})\otimes\rho_{0}$ at the initial time
$t_{i}=-t_{D}$, in the time interval $(-t_{D},0)$ the time-evolution is
generated by the total Hamiltonian $H_{tot}=H(\lambda_{0})+H_{D}+H_{I}$. Then,
in the time interval $(0,\tau)$ the qubit and the system do not interact and
the quench is performed. Finally, in the time interval
$(\tau,\tau+t^{\prime}_{D})$ the time-evolution is generated by the total
Hamiltonian $H^{\prime}_{tot}=H(\lambda_{\tau})+H_{D}+H^{\prime}_{I}$. The
coherence of the qubit at the final time $t_{f}=\tau+t^{\prime}_{D}$ reads

$\displaystyle\langle e|\rho_{D}(t_{f})|g\rangle$ $\displaystyle=$
$\displaystyle\langle e|\rho_{D}(t_{i})|g\rangle
e^{-i\omega(t_{f}-t_{i})}\text{Tr}\bigg{\\{}e^{-i(1-\delta_{e})t_{D}H(\lambda_{0})}\rho_{0}$
(177) $\displaystyle\times
e^{i(1-\delta_{g})t_{D}H(\lambda_{0})}U^{\dagger}_{\tau,0}e^{i(\delta^{\prime}_{e}-\delta^{\prime}_{g})t^{\prime}_{D}H(\lambda_{\tau})}U_{\tau,0}\bigg{\\}}$

from which we can determine $X_{q}(u)$.

## References

  * [1] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * [2] A. Polkovnikov, K. Sengupta, A. Silva, and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011).
  * [3] M. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, Rev. Mod. Phys. 83, 1405 (2011).
  * [4] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, Adv. Phys. 65, 239 (2016).
  * [5] J. Dziarmaga, Adv. Phys. 59, 1063 (2010).
  * [6] J. Goold, F. Plastina, A. Gambassi, and A. Silva, ’The role of quantum work statistics in many-body physics,’ in Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions, 317–336 (2018).
  * [7] P. Talkner, E. Lutz, and P. Hänggi, Phys. Rev. E 75, 050102(R) (2007).
  * [8] M. Perarnau-Llobet, E. Bäumer, K. V. Hovhannisyan, M. Huber, and A. Acin, Phys. Rev. Lett. 118, 070601 (2017).
  * [9] M. Lostaglio, Phys. Rev. Lett. 120, 040602 (2018).
  * [10] E. Wigner, Phys. Rev. 40, 749 (1932).
  * [11] P. Solinas and S. Gasparinetti, Phys. Rev. E 92, 042150 (2015).
  * [12] A. E. Allahverdyan, Phys. Rev. E 90, 032137 (2014).
  * [13] G. Francica, Phys. Rev. E 105, 014101 (2022).
  * [14] G. Francica, Phys. Rev. E 106, 054129 (2022).
  * [15] Bao-Ming Xu, Jian Zou, Li-Sha Guo, and Xiang-Mu Kong, Phys. Rev. A 97, 052122 (2018).
  * [16] M. G. Díaz, G. Guarnieri, and M. Paternostro, Entropy 22, 1223 (2020).
  * [17] A. Gambassi, and A. Silva, Phys. Rev. Lett. 109, 250602 (2012).
  * [18] S. Sotiriadis, A. Gambassi, and A. Silva, Phys. Rev. E 87, 052129, (2013).
  * [19] G. Perfetto, L. Piroli, and A. Gambassi, Phys. Rev. E 100, 032114, (2019).
  * [20] A. Silva, Phys. Rev. Lett. 101, 120603 (2008).
  * [21] R. Dorner, J. Goold, C. Cormick, M. Paternostro, and V. Vedral, Phys. Rev. Lett. 109, 160601 (2012).
  * [22] L. Fusco, S. Pigeon, T. J. G. Apollaro, A. Xuereb, L. Mazzola, M. Campisi, A. Ferraro, M. Paternostro, and G. De Chiara, Phys. Rev. X 4, 031029 (2014).
  * [23] A. Russomanno, S. Sharma, A. Dutta, and G. E Santoro, Journal of Statistical Mechanics: Theory and Experiment 2015, P08030 (2015).
  * [24] N. O. Abeling and S. Kehrein, Phys. Rev. B 93, 104302 (2016).
  * [25] K. Zawadzki, A. Kiely, G. T. Landi, and S. Campbell, Phys. Rev. A 107, 012209 (2023).
  * [26] Z. Fei and H. T. Quan, Phys. Rev. Research 1, 033175 (2019).
  * [27] Q. Wang, D. Cao, and H. T. Quan, Phys. Rev. E 98, 022107 (2018).
  * [28] Z. Fei, N. Freitas, V. Cavina, H. T. Quan, and M. Esposito, Phys. Rev. Lett. 124, 170603 (2020).
  * [29] R. W. Spekkens, Phys. Rev. Lett. 101, 020401 (2008).
  * [30] C. Jarzynski, Phys. Rev. Lett. 78, 2690 (1997).