Many-body work distributions

TL;DR

This paper derives exact formulas for the work distribution in non-interacting quantum many-body systems under external classical influences, accounting for time-dependent frequencies and boundary movements.

Contribution

It provides the first exact expressions for the characteristic function of work distributions in such systems, including arbitrary switching and coupling functions.

Findings

Derived exact characteristic functions for work distributions.

Analyzed effects of time-dependent frequencies and boundary movements.

Discussed limiting cases and their physical implications.

Abstract

The work distribution function for a non-relativistic, non-interacting quantum many-body system interacting with classical external sources is investigated. Exact expressions for the characteristic function corresponding to the work distribution function is obtained for arbitrary switching function and coupling functions. The many-body frequencies are assumed to be generally time-dependent in order to take into account the possibility of moving the boundaries of the system in a predefined process linking the characteristic function to the fluctuation-induced energies in confined geometries. Some limiting cases are considered and discussed.

Equations

$$\hat{H}$$

$$G(t)=\left\{\begin{array}[]{ll}g(t),&0\leq t\leq\tau;\\ 0,&t>\tau.\end{array}\right.$$

$$\sum _{\vec{k}}\to \int \frac{d^3\vec{k}}{(2\pi {)}^3}.$$

$${\dot{\hat{a}}}_{\mathbf{k}}=-i{\omega }_{\mathbf{k}}(t){\hat{a}}_{\mathbf{k}}-iG(t)\mathcal{F}(\mathbf{k}),$$

$$\hat{a}_{\mathbf{k}}(t)=e^{-i\zeta_{\mathbf{k}}(t)}\big{(}\hat{a}_{\mathbf{k}}(0)-\xi_{\mathbf{k}}(t)\big{)},$$

$${\zeta }_{\mathbf{k}}(t)$$

$${\xi }_{\mathbf{k}}(t)$$

$${\hat{D}}_{\mathbf{k}}({\xi }_{\mathbf{k}})=e^{{\xi }_{\mathbf{k}}(t){\hat{a}}_{\mathbf{k}}^{†}(0)-{\xi }_{\mathbf{k}}^{\ast }(t){\hat{a}}_{\mathbf{k}}},$$

$${\hat{D}}_{\mathbf{k}}({\xi }_{\mathbf{k}}){\hat{a}}_{\mathbf{k}}(0){\hat{D}}_{\mathbf{k}}^{†}({\xi }_{\mathbf{k}})={\hat{a}}_{\mathbf{k}}(0)-{\xi }_{\mathbf{k}}(t).$$

$${\hat{a}}_{\mathbf{k}}(t)={\hat{U}}_{\mathbf{k}}^{†}(t){\hat{a}}_{\mathbf{k}}(0){\hat{U}}_{\mathbf{k}}(t),$$

$${\hat{U}}_{\mathbf{k}}(t)=e^{i\theta (t)}{e}^{-i{\zeta }_{\mathbf{k}}(t)({\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}+1/2)}{\hat{D}}_{\mathbf{k}}^{†}({\xi }_{\mathbf{k}}).$$

$${\hat{U}}_{0\mathbf{k}}(t)$$

$${\hat{U}}_I(t)$$

$${\hat{D}}_{\mathbf{k}}({\alpha }_{\mathbf{k}})H_0^{\mathbf{k}}(t){\hat{D}}_{\mathbf{k}}^{†}({\alpha }_{\mathbf{k}})$$

$${\alpha }_{\mathbf{k}}(t)=-\frac{G(t)}{{\omega }_{\mathbf{k}}(t)}\mathcal{F}(\mathbf{k}),$$

$${\hat{H}}^{\mathbf{k}}(t)$$

$$E_n^{\prime }=\hslash {\omega }_{\mathbf{k}}(t)(n+\frac{1}{2})-\frac{\hslash G^2(t)|\mathcal{F}(\mathbf{k}){|}^2}{{\omega }_{\mathbf{k}}(t)}.$$

$$P_k(W,t)=\sum _{i,j}{p}_0(E_i)p_t(|E_i⟩\to |E_j^{\prime }⟩)\delta (W-E_j^{\prime }+E_i),$$

$${\tilde{\rho }}_{\mathbf{k}}(0)=\sum _{i}p(E_i)|E_i⟩⟨E_i|,$$

$$G_{\mathbf{k}}(\nu ,t)$$

$$p_t(|E_i⟩\to |E_j^{\prime }⟩)=|⟨E_j^{\prime }|{\hat{U}}_{\mathbf{k}}(t)|E_i⟩{|}^2,$$

$$G_{\mathbf{k}}(\nu ,t)$$

$$G_{\mathbf{k}}^{(n)}(\nu ,t)$$

$${△}_{\mathbf{k}}(t)$$

$${\eta }_{\mathbf{k}}(t)$$

$$G_{\mathbf{k}}^{(n)}(\nu ,t)$$

$$L_n(x)=\sum _{l=0}^n\left(\genfrac{}{}{0ex}{}{n}{l}\right)\frac{(-x{)}^l}{l!}.$$

$$P_k^{(n)}(W,t)$$

$$q_s(n,|{\eta }_{\mathbf{k}}{|}^2)=\sum _{l=0}^n\sum _{p=\text{mbox}max(l-s,0)}^{2l}(\genfrac{}{}{0ex}{}{n}{l}\left)\right(\genfrac{}{}{0ex}{}{2l}{p})(-1{)}^p\frac{e^{-|{\eta }_{\mathbf{k}}{|}^2}|{\eta }_{\mathbf{k}}{|}^{2(s+p)}}{l!(s+p-l)!}.$$

$$z$$

Full text

1 Introduction

While discussing quantum thermodynamics of a small system, quantum fluctuations being of the same order of magnitude as expectation values, play an important role. Fluctuation theorems describe the connection between non-equilibrium fluctuations and thermal equilibrium states of small systems. They also describe the nonlinear response of a system under the influence of a time-dependent external source [1, 2]. Various fluctuation and work theorems have been studied [3, 4]. Fluctuation relations were initially derived for classical systems [1, 5, 6, 7] and experimentally confirmed in some classical systems [8, 9, 10, 11]. The generalization of the fluctuation theorems to the realm of small quantum mechanical systems has led to considerable progress in formulating quantum thermodynamics [12, 13, 14, 15, 16, 17, 18, 19].

An important quantity in this context is the amount of work that can be extracted or done on a quantum system in a nonequilibrium process. But contrary to classical statistical mechanics work is a subtle concept in quantum thermodynamics and it is not an observable [20]. Therefore, indirect methods should be applied to measure work. The most established method is the two-point measurement scheme which is based on projective measurements of energy at the end-points of a process [21].

Projective measurements when applied to relativistic quantum systems cause inconsistencies. Since due to the relativistic structure, the projective measurements can not be localized [22] and lead to the possibility of superluminal signaling [23], so they should be banished from any consistent relativistic quantum theory [23, 24, 25]. An effective way to remove these inconsistencies is by coupling the main system locally to auxiliary systems or ancillae like detectors, atoms or qubits. Using a controllable ancilla interacting with the main system have enabled the first experimental characterization of quantum fluctuation relations [26]. In a recent work [27], the work distributions on a relativistic quantum field has been investigated in the framework of Ramsey interferometric scheme [26, 28, 29].

The Hamiltonian Eq. (1) is a collection of non-interacting quantum driven harmonic oscillators specified by an index $\mathbf{k}$ that can be a collection of particle characteristics like polarization ($\lambda$), spin ($s$) or wavenumber ($\vec{k}$), ($\mathbf{k}=\{\lambda,s,\vec{k}\}$). The distribution of work for a driven quantum harmonic oscillator with a time-independent frequency has been investigated in [18, 30, 31, 32].

In this letter, we investigate the work distribution function for a non-relativistic quantum many-body system interacting with classical external sources. The duration of the interaction of external sources with the many-body system is controlled by inserting a switching function into the interaction part of the Hamiltonian. For generality, the many-body frequencies are assumed to be time-dependent. This can arise as a result of applying external sources to the system known as Lamb-shifts or due to moving the boundaries of the system. Here we ignore from the Lamb-shifts caused by external sources and just consider the effects due to moving boundaries. In the absence of external sources, the change in the energy of the system due to the variations in the boundary of the system is studied in the context of the Casimir physics [33] and we have shown that it is connected to the characteristic function. We have found an exact characteristic function for arbitrary coupling functions and switching function and considered two interesting limiting cases.

2 The model

Consider a non-relativistic, non-interacting quantum many-body system interacting with external sources described by the second quantized Hamiltonian

[TABLE]

where $\mathscr{F}(\mathbf{k})(\mathscr{F}^{*}(\mathbf{k}))$ are coupling functions between many-body system and external sources and $G(t)$ is a switching time-dependent function defined by

[TABLE]

Note that the Hamiltonian Eq. (1) is the sum of independent Hamiltonians $\hat{H}^{\mathbf{k}}(t)$, ($[\hat{H}^{\mathbf{k}}(t),\hat{H}^{\mathbf{k}^{\prime}}(t)]=0$). The possibility of moving the boundaries of the many-body system and its effect on the energy spectrum of the system has been taken into account by considering the frequencies $\omega_{\mathbf{k}}(t)$ to be generally time-dependent. So, one can affect the many-body system by two processes: (i) By applying external classical sources, like the interaction of external electromagnetic fields with the system. (ii) By moving the boundaries of the system, for example in the case of a bosonic gas confined between parallel conducting plates, the distance between the plates can be considered as a varying parameter.

Here we have adopted the lab quantization, so the wave number in $\mathbf{k}$, denoted by $\vec{k}$ takes on discrete values, if the volume of the lab tends to infinity, the discrete sums over $\vec{k}$ should be replaced by the integral

[TABLE]

3 The time-evolution operator

Let us focus on the $\mathbf{k}$th subsystem in the Hamiltonian Eq. (1), from the Heisenberg equations of motion, we find for the annihilation operator $\hat{a}_{\mathbf{k}}$

[TABLE]

with the formal solution

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

From Eq. (2) it is clear that $\xi_{\mathbf{k}}(t)$ is time-independent for $t\geq\tau$, ($\xi_{\mathbf{k}}(t)=\xi_{\mathbf{k}}(\tau)$).

An important operator in quantum optics terminology is the displacement operator

[TABLE]

with displacement parameter $\xi_{\mathbf{k}}(t)$, which acts on the annihilation operator as

[TABLE]

The displacement operator when applied to the vacuum state $|0\rangle$ produces a coherent state $|\xi_{\mathbf{k}}(t)\rangle=\hat{D}_{\mathbf{k}}(\xi_{\mathbf{k}})|0\rangle$ which is an eigenket of the annihilation operator $\hat{a}_{\mathbf{k}}|\xi_{\mathbf{k}}(t)\rangle=\xi_{\mathbf{k}}(t)|\xi_{\mathbf{k}}(t)\rangle$.

In Heisenberg picture, the time-evolution of the annihilation operator is given by

[TABLE]

where $\hat{U}_{\mathbf{k}}(t)$ is the time-evolution operator corresponding to the kth subsystem with the Hamiltonian $\hat{H}^{\mathbf{k}}(t)$. One easily finds (A)

[TABLE]

The time-dependent Hamiltonian $\hat{H}_{0}^{\mathbf{k}}(t)$ satisfies $[\hat{H}_{0}^{\mathbf{k}}(t)$, $\hat{H}_{0}^{\mathbf{k}}(t^{\prime})]=0$, therefore, in the absence of the external sources ($\mathscr{F}(\mathbf{k})=0$), the time-evolution operator for the $\mathbf{k}$th subsystem is given by

[TABLE]

The time-evolution operator for the $\mathbf{k}$th subsystem in the interaction picture is by definition

[TABLE]

that is the time-evolution operator in the interaction picture is the adjoint of the displacement operator.

4 Characteristic function for a single subsystem

A quantum system responds to external sources that can be characterized by finding the variation of the energy of the system including the interaction with external sources. In this context, the work $W$ is defined as the difference between the final and initial energies of the system plus interaction. Let us focus on the $\mathbf{k}$th subsystem with Hamiltonian $\hat{H}^{\mathbf{k}}(t)$, and let $E_{i},\,(i=1,2,\cdots)$ be the energy spectrum of the Hamiltonian $\hat{H}^{\mathbf{k}}(0)$, at $t=0$, and let $E^{\prime}_{i},\,(i=1,2,\cdots)$, be the energy spectrum of the Hamiltonian $\hat{H}^{\mathbf{k}}(t)$ at $t>0$. To find the energy spectrum of the Hamiltonian $\hat{H}^{\mathbf{k}}(t)$, we consider the unitary transformation

[TABLE]

by setting

[TABLE]

we will find

[TABLE]

where $\alpha_{\mathbf{k}}$ is given by Eq. (15). According to Eq. (16), if $|\psi_{n}\rangle$ is an eigenket of the Hamiltonian $H^{\mathbf{k}}_{0}(t)$ with eigenvalue $E_{n}=\hbar\omega_{\mathbf{k}}(t)(n+1/2)$ then $\hat{D}_{\mathbf{k}}(\alpha_{\mathbf{k}})|\psi_{n}\rangle$ is an eigenket of the Hamiltonian $\hat{H}^{\mathbf{k}}(t)$ with eigenvalue

[TABLE]

The characteristic function corresponding to the $\mathbf{k}$th subsystem is denoted by $G_{\mathbf{k}}(\nu,t)$ and provides a complete statistical description of the work performed on or extracted from the subsystem. The characteristic function is defined as the Fourier transform of the probability density of the work $P_{k}(W,t)$ [21]

[TABLE]

where $p_{0}(E_{i})$ is the probability that the kth subsystem be in the eigenstate $|E_{i}\rangle$ with the eigenvalue $E_{i}$ at $t=0$, and $p_{t}(|E_{i}\rangle\rightarrow|E^{\prime}_{j}\rangle)$ is the transition-probability from the initial state $|E_{i}\rangle$ to the final state $|E^{\prime}_{j}\rangle$ at time $t$. Let the diagonal part of the initial density matrix $\hat{\rho}_{\mathbf{k}}(0)$, in the basis of energy eigenkets $\{|E_{i}\rangle\}$, be denoted by

[TABLE]

then from the definition of the characteristic function

[TABLE]

and the transition probability

[TABLE]

we can rewrite Eq. (20) as (B)

[TABLE]

4.1 Example 1.

Let the initial state be the number state $\hat{\rho}_{\mathbf{k}}(0)=|n\rangle\langle n|=\tilde{\rho}_{\mathbf{k}}(0)$, then

[TABLE]

after straightforward calculations, one finds (C)

[TABLE]

where $L_{n}(x)$ is the Laguerre polynomial of degree $n$ and

[TABLE]

In the special case where the boundaries of the system are fixed, the frequencies are time-independent ($\omega_{\mathbf{k}}(t)=\omega_{\mathbf{k}}(0)$), so $\triangle_{\mathbf{k}}(t)=0$ and the Eq. (24) tends to the same result reported in [32].

Now let us find the work distribution function corresponding to the characteristic function Eq.(24). For this purpose, we rewrite Eq. (24) as follows

[TABLE]

where we used the following expansion for the Laguerre polynomial of the order n

[TABLE]

The work distribution function can be obtained from Eq. (26) by taking the inverse Fourier transform. After straightforward calculations, we will find (D)

[TABLE]

From Eq. (4.1) the weight functions are defined by [32]

[TABLE]

In Eq. (29) the dimensionless quantity $|\eta_{\mathbf{k}}|$ corresponds to the rapidity parameter defined in [32] for $\triangle_{\mathbf{k}}(t)=0$. Note that for $t>\tau$ we have $G(t)=\alpha_{\mathbf{k}}=0$ (see Eqs. (2), (15)), therefore, from Eqs. (4.1), (7) we have

[TABLE]

For the case where the frequencies are time-independent $\omega_{\mathbf{k}}(t)=\omega_{\mathbf{k}}(0)$, Eq. (30) can be simplified as

[TABLE]

where $\tilde{G}(\omega_{\mathbf{k}}(0))$ is the Fourier transform of the switching function $G(t)$.

The weight functions $q_{s}(3,z)$ in terms of the rapidity $z$ have been depicted in Fig.1 for $s=-3,...,3$, ($\rho_{k}(0)=|3\rangle\langle 3|$) and in Fig.2 for $s=0,1,2,3$, ($\rho_{k}(0)=|0\rangle\langle 0|$). In both figures the location of the maximum of weight function has been shifted to the larger values of $z$ by increasing $|s|$.

4.2 Example 2.

Let the initial density matrix be a thermal state with inverse temperature $\beta=1/k_{B}T$ ($k_{B}$ is the Boltzmann constant)

[TABLE]

where $z_{k}(0)$ is the partition function

[TABLE]

then we will find the characteristic function $G^{th}_{\mathbf{k}}(\nu,t)$ as (D)

[TABLE]

which is in agreement with the result reported in [32] when $\triangle_{\mathbf{k}}=0$.

An interesting case arises when the external sources are switched off ($\xi_{\mathbf{k}}=0$) but the boundaries of the many-body system change with time in a predefined process. In this case $\alpha_{\mathbf{k}}=\eta_{\mathbf{k}}=\xi_{\mathbf{k}}=0$, therefore,

[TABLE]

To find the work distribution function we rewrite Eq. (35) as

[TABLE]

and by taking the inverse Fourier transform we find

[TABLE]

with the weight functions

[TABLE]

The probabilities $q_{r}$ correspond to the initial thermal state Eq. (32) and the sign of the work $W$ depend on the sign of $\triangle_{\mathbf{k}}(t)$. If work is done on a system the eigenfrequency increases ($\triangle_{\mathbf{k}}(t)>0$) and when work is delivered by the system the eigenfrequency decreases ($\triangle_{\mathbf{k}}(t)<0$).

4.3 Example 3.

Let the initial state be a coherent state $\hat{\rho}_{\mathbf{k}}(0)=|\alpha\rangle\langle\alpha|$,

[TABLE]

then

[TABLE]

and the characteristic function is

[TABLE]

where $G^{(n)}_{\mathbf{k}}(\nu,t)$ is given by Eq. (26). Therefore, the corresponding work distribution function is

[TABLE]

where $P^{(n)}_{k}(W,t)$ is given by Eq. (4.1). From Eq. (42) it is clear that the work distribution function in this case is the average of $P^{(n)}_{\mathbf{k}}(W,t)$ with respect to the Poisson distribution function. By inserting Eq. (24) into Eq. (41) we have

[TABLE]

now by making use of the identity [34]

[TABLE]

we finally find

[TABLE]

where $J_{0}(x)$ is the Bessel function of order zero. The Eq. (43) is in agreement with the result reported in [32] when $\triangle_{\mathbf{k}}=0$. In the absence of external sources we have

[TABLE]

By taking the inverse Fourier transform of Eq. (46), we find the work distribution function as

[TABLE]

where the weight functions have a poisson distribution

[TABLE]

5 The characteristic function of the total system

To find the characteristic function of the total system we can easily proof that if the total Hamiltonian $\hat{H}$ is the sum of noninteracting Hamiltonians $\hat{H}_{\mathbf{k}}$

[TABLE]

then the characteristic function of the total system is the multiplication of the characteristic functions of the subsystems given by (E)

[TABLE]

The Eq. (50), can be rewritten as

[TABLE]

where $G_{\mathbf{k}}(\nu,t)$ is given by Eq. (22).

Now from the Eqs. (22), (50), and the fact that the subsystems defined by Hamiltonians $\hat{H}^{\mathbf{k}}(t)$ are non-interacting, we can generalise the Eq. (22) as

[TABLE]

where

[TABLE]

Also from Eq. (20) and its inverse, we find

[TABLE]

and

[TABLE]

6 Expectation values and cumulant function

From the definition of the characteristic function Eq. (20), the nth moment of work can be obtained as

[TABLE]

To facilitate the calculations, let us consider the characteristic function Eq. (34) in zero temperature ($\beta\rightarrow\infty$), in this limit, we have

[TABLE]

The first and the second moments of the work are respectively given by

[TABLE]

with the variance

[TABLE]

In general, by inserting Eq. (51) into the cumulant function defined by [21]

[TABLE]

we find

[TABLE]

From the cumulant function, we find the nth cumulant as

[TABLE]

where $C_{1}=\langle W\rangle$ and $C_{2}=(\Delta W)^{2}$. If the characteristic function is given by Eq. (57), the skewness, which determines the deviation of $P(W)$ from a Gaussian distribution, is

[TABLE]

6.1 Casimir energy

Consider an arbitrary fluctuating field confined in a region by imposing Dirichlet boundary conditions. The region can be contracted or expanded by varying external parameters. Casimir energy is defined as the difference between zero-point energies of the system in the presence and absence of boundaries. As an illustration, consider a photonic gas confined between parallel, neutral, perfect conductors with unit area. The external parameter in this case is the distance between conductors ($d$). The Hamiltonian of the photonic gas is like Eq. (1), without external sources, where the discrete index $\mathbf{k}$ refers to the wave number and polarization of photons. In zero temperature ($\beta\rightarrow\infty$), from Eq. (35) we have

[TABLE]

the total characteristic function is

[TABLE]

and the work distribution function is obtained as

[TABLE]

with zero variance, that is

[TABLE]

as expected. The sum over $\mathbf{k}$ in Eq. (67) can be achieved in different ways and an explicit dependence on $d$ is obtained as [33]

[TABLE]

In general, by setting $\beta=i\nu$ in Eq. (52), the quantum version of Jarzynski equality is recovered [20]

[TABLE]

where $\triangle F=F(\tau)-F(0)$, represents the difference of free energies at times $t=\tau$ and $t=0$, or equivalently for different configurations of boundaries. Therefore,

[TABLE]

By inserting Eq. (35 into Eq. (70), one can find the temperature correction to the Casimir energy.

7 Summary and conclusions

We have considered a non-relativistic non-interacting quantum many-body system under the influence of external classical sources. The external sources are coupled to the many-body system by the coupling functions $\mathscr{F}(\mathbf{k})\,(\mathscr{F}^{*}(\mathbf{k}))$ in a finite interval of time $\tau$ controlled by the switching function $G(t)$. The switching function and the coupling functions are arbitrary functions and one can assign them desired explicit forms depending on the characteristics of the many-body system and the external sources. We have also taken into account the possibility of moving boundaries of the many-body system in a predefined process by considering the frequencies of the system to be generally time-dependent. We have obtained the closed form expression Eq. (51) for the characteristic function $G(\nu,t)$ corresponding to the work distribution function $P(W)$. Two interesting limiting cases are considered (i) There are external sources but the boundary of the system do not change, in this case the system frequencies are time-independent and the characteristic function of the system can be obtained using Eqs.(34) and Eq. (51). (ii) The external sources are switched off but the boundaries of the system can be moved in a predefined process. The characteristic function in this case is given by Eqs. (35) and (51). In a process where the configuration of the boundaries of the many-body system at the end points are different but initial and final states of the system are thermal equilibrium states, the change in free energy is given by Eq. (70) which is a useful formula in the context of the fluctuation-induced forces.

Appendix A Proof of Eq. (11)

By setting

[TABLE]

in the Baker–Campbell–Hausdorff (BCH) formula

[TABLE]

we find

[TABLE]

and similarly by setting

[TABLE]

we find

[TABLE]

Now by making use of Eqs. (72) and (75) we have

[TABLE]

therefore, $\hat{U}_{\mathbf{k}}(t)$ is the quantum propagator corresponding to the Hamiltonian $\hat{H}^{\mathbf{k}}(t)$. By inserting $\hat{U}_{\mathbf{k}}(t)$ into the Schrödinger equation

[TABLE]

we find the phase factor $\theta(t)$ in Eq. (11) as

[TABLE]

Appendix B Proof of Eq. (22)

From Eqs. (20), (21) we have

[TABLE]

Appendix C Proof of Eq. (24)

From Eq. (23) we have

[TABLE]

also from Eq. (16) we have

[TABLE]

so

[TABLE]

By making use of the identity

[TABLE]

we will find

[TABLE]

therefore,

[TABLE]

finally, using the formulas

[TABLE]

and [35]

[TABLE]

one finds

[TABLE]

where $L_{n}(x)$ is the Laguerre polynomial of degree $n$.

Appendix D Proof of Eqs. (4.1) and (34)

D.1 Eq.(4.1

From Eq. (26 we have

[TABLE]

Now by taking the inverse Fourier transform (see Eq. (20)) and a redefinition of variables $s=r+l-p$, we recover Eq. (4.1).

D.2 Eq.(34)

For a thermal state we have

[TABLE]

where

[TABLE]

therefore,

[TABLE]

Now using the identity [35]

[TABLE]

we finally find Eq. (34).

Appendix E Proof of Eq. (50)

For notational simplicity, and with no lose of generality, let us replace the index $\mathbf{k}$ of mutually-independent Hamiltonians ($\hat{H}^{\mathbf{k}}$) by a single index $j$ as $\hat{H}_{j}$. Let $|E_{k_{j}}\rangle$ be an eigenket of the Hamiltonian $\hat{H}_{j}(0)$ with eigenvalue $E_{k_{j}}$ and $|E^{\prime}_{k_{j}}\rangle$ be an eigenket of the Hamiltonian $\hat{H}_{j}(t)$ with eigenvalue $E^{\prime}_{k_{j}}$. Also, let the probability that the jth subsystem being initially at the eigenket $|E_{i_{j}}\rangle$ be denoted by $P_{j}(E_{i_{j}})$. Then, from the definition of work distribution function we can write

[TABLE]

The corresponding characteristic function $G(\nu,t)$ is the Fourier transform of $P(W)$

[TABLE]

By inserting Eq. (94) into Eq. (95) and doing the integral over $W$, we find

[TABLE]

therefore,

[TABLE]