Quantum work distributions associated with the dynamical Casimir effect

TL;DR

This paper analyzes the quantum work distributions in a cavity with an oscillating boundary, demonstrating the validity of fluctuation theorems and the transition from quantum to classical behavior at different temperatures.

Contribution

It provides analytical expressions for work distributions in a dynamical Casimir effect setup, extending fluctuation theorem validation to particle-creating processes.

Findings

Work distributions obey fluctuation theorems in nonequilibrium conditions.

High temperature limit yields classical work behavior.

Low temperature limit shows significant quantum effects similar to Casimir phenomena.

Abstract

We study the joint probability distribution function of the work and the change of photon number of the nonequilibrium process of driving the electromagnetic (EM) field in a three-dimensional cavity with an oscillating boundary. The system is initially prepared in a grand canonical equilibrium state and we obtain the analytical expressions of the characteristic functions of work distributions in the single-resonance and multiple-resonance conditions. Our study demonstrates the validity of the fluctuation theorems of the grand canonical ensemble in nonequilibrium processes with particle creation and annihilation. In addition, our work illustrates that in the high temperature limit, the work done on the quantized EM field approaches its classical counterpart; while in the low temperature limit, similar to Casimir effect, it differs significantly from its classical counterpart.

Figures

Tables

Table 1. Table 1: Moments of work in two limits of the initial temperature. In the low temperature limit, the nonzero works ⟨ w ⟩ 1 subscript delimited-⟨⟩ 𝑤 1 \left\langle w\right\rangle_{1} and ⟨ w ⟩ 2 subscript delimited-⟨⟩ 𝑤 2 \left\langle w\right\rangle_{2} are manifestations of the quantum nature of the EM field (zero point energy). When ℏ → 0 → Planck-constant-over-2-pi 0 \hbar\to 0 , the works ⟨ w ⟩ 1 subscript delimited-⟨⟩ 𝑤 1 \left\langle w\right\rangle_{1} and ⟨ w ⟩ 2 subscript delimited-⟨⟩ 𝑤 2 \left\langle w\right\rangle_{2} become zero, which agrees with our intuition that no work is done on the classical vacuum.

	Low-temperature limit ($\beta \to \infty$)	High-temperature limit ($\beta \to 0$)
${⟨w⟩}_1$	$\hslash \omega {\sinh }^2{g}_1\tau$	$2{\beta }^{-1}{\sinh }^2{g}_1\tau$
${\sigma }_{w,1}^2$	$\frac{1}{2}{\left(\hslash \omega \right)}^2{\sinh }^{2}2g_1\tau$	$4{\beta }^{-2}\cosh 2g_1{\tau }_1{\sinh }^2{g}_1\tau$
${⟨w⟩}_2$	$\left(1+r\right)\hslash \omega {\sinh }^2{g}_2\tau$	${\beta }^{-1}\frac{{\left(1+r\right)}^2}{r}{\sinh }^2{g}_2\tau$
${\sigma }_{w,2}^2$	$\frac{{\left(1+r\right)}^2}{4}{\left(\hslash \omega \right)}^2{\sinh }^{2}2g_2\tau$	${\beta }^{-2}\frac{{\left(1+r\right)}^4}{r^2}\left({\sinh }^2{g}_2{\tau }_2+\frac{2r}{{\left(1+r\right)}^2}\right){\sinh }^2{g}_2\tau$
${⟨w⟩}_3$	0	${\beta }^{-1}\frac{{\left(r-1\right)}^2}{r}{\sin }^2{g}_3\tau$
${\sigma }_{w,3}^2$	0	${\beta }^{-2}\frac{{\left(1-r\right)}^4}{r^2}\left({\sin }^2{g}_3{\tau }_3+\frac{2r}{{\left(1-r\right)}^2}\right){\sin }^2{g}_3\tau$

Equations

$$⟨e^{-\beta w}⟩=e^{-\beta \Delta F}.$$

$$⟨e^{-\beta (w-\mu \Delta N)}⟩=e^{-\beta \Delta \Phi },$$

$$\varphi =-\frac{1}{\epsilon }\nabla \cdot {\mathbf{\Pi }}_e,\mathbf{A}=\mu \frac{\partial {\mathbf{\Pi }}_{\mathrm{e}}}{\partial t}+\nabla \times {\mathbf{\Pi }}_{\mathrm{m}},$$

$$\left({\nabla }^2-{\epsilon }_0{\mu }_0{\partial }_t^2\right){\psi }^{\mathrm{TE},\mathrm{TM}}=0,$$

$$\mathcal{L}(\mathbf{r},t)=\frac{1}{2}\left(-\dot{\psi }{\nabla }_{\perp }^2\dot{\psi }-{\nabla }^2\psi {\nabla }_{\perp }^2\psi \right).$$

$${\psi |}_{z=0,L_z}=0,{\frac{\partial \psi }{\partial x}|}_{x=0,L_x}={\frac{\partial \psi }{\partial y}|}_{y=0,L_y}=0,$$

$$\left[\hat{\psi }(\mathbf{r},t),\hat{\psi }({\mathbf{r}}^{\prime },t)\right]$$

$$\left[\hat{\psi }(\mathbf{r},t),\hat{\pi }({\mathbf{r}}^{\prime },t)\right]$$

$$\hat{\pi }(\mathbf{r},t)=-{\nabla }_{\perp }^2\frac{\partial \hat{\psi }(\mathbf{r},t)}{\partial t}.$$

$${\psi }_{\mathbf{k}}(\mathbf{r})$$

$$\hat{\psi }\left(\mathbf{r},t\right)$$

$$\hat{\pi }(\mathbf{r},t)$$

$$\left[{\hat{a}}_{\mathbf{k}},{\hat{a}}_{{\mathbf{k}}^{\prime }}\right]=0,\left[{\hat{a}}_{\mathbf{k}},{\hat{a}}_{{\mathbf{k}}^{\prime }}^{†}\right]={\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }}.$$

$$\hat{H}=\int d\mathbf{r}\left[\frac{\partial \hat{\psi }}{\partial t}\hat{\pi }-\hat{L}\right]=\sum _{\mathbf{k}}\frac{{\omega }_{\mathbf{k}}}{2}\left({\hat{a}}_{\mathbf{k}}^{†}{\hat{a}}_{\mathbf{k}}+{\hat{a}}_{\mathbf{k}}{\hat{a}}_{\mathbf{k}}^{†}\right).$$

$${\nabla }^2{\psi }_{\mathbf{k},\lambda }\left(\mathbf{r}\right)+{\omega }_{\mathbf{k},\lambda }^2{\psi }_{\mathbf{k},\lambda }(\mathbf{r})=0,$$

$${{\psi }_{\mathbf{k},\lambda }(\mathbf{r})|}_{z=0,\lambda }=0,$$

$${\frac{\partial {\psi }_{\mathbf{k},\lambda }(\mathbf{r})}{\partial x}|}_{x=0,L_x}={\frac{\partial {\psi }_{\mathbf{k},\lambda }(\mathbf{r})}{\partial y}|}_{y=0,L_y}=0,$$

$$\hat{\psi }(\mathbf{r},t)$$

$$\hat{\pi }(\mathbf{r},t)$$

$$\frac{d{\hat{Q}}_{\mathbf{k}}(t)}{dt}$$

$$\frac{d{\hat{\mathcal{P}}}_{\mathbf{k}}(t)}{dt}$$

$${\tilde{g}}_{\mathbf{kp}}(t)$$

$$\displaystyle g_{\mathbf{kp}}=\left\{\begin{array}[]{ll}\frac{(-1)^{k_{z}+p_{z}}2k_{z}p_{z}}{k_{z}^{2}-p_{z}^{2}}\delta_{{\mathbf{k}}_{\perp},\mathbf{p}_{\perp}},&k_{z}\neq p_{z}\\ 0,&k_{z}=p_{z}\end{array}\right..$$

$${\hat{H}}_{\mathrm{eff}}(t)$$

$${\hat{a}}_{\mathbf{k}}(t)$$

$${\hat{a}}_{\mathbf{k}}^{†}(t)$$

$$\frac{d}{dt}{\hat{a}}_{\mathbf{k}}(t)$$

$$h_{\mathbf{kp}}(t)$$

$$d_{\mathbf{kp}}(t)$$

$${\hat{H}}_{\mathrm{eff}}(t)$$

$$\lambda (t)={\lambda }_0\left[1+ϵ\sin (\Omega t)\right],$$

$${\hat{V}}_1^{\mathrm{I}}$$

$$g_1$$

$${\hat{V}}_2^{\mathrm{I}}$$

Full text

I Introduction

Fluctuation theorems have attracted considerable attentions in the field of nonequilibrium thermodynamics in the last two decades. One of the most important results is the Jarzynski equality no1997

[TABLE]

It connects the equilibrium free energy difference $\Delta F$ with the fluctuating nonequilibrium work $w$ done on a system initially prepared in a canonical equilibrium state at the inverse temperature $\beta=1/k_{B}T$. This equality is first derived in the classical regime and later generalized to the quantum regime aq2000 ; ja2000 ; wo2007 . The validity has been tested experimentally in various systems  eq2002 ; ve2005 ; ex2015 .

For the grand canonical equilibrium initial state, besides the trajectory work $w$, the particle number along every “trajectory” is also a fluctuating quantity and the Jarzynski equality takes a similar form st2007 ; th2009 ; wo2012

[TABLE]

where $\mu$ is the chemical potential of the initial state, $\Delta N$ and $\Delta\Phi$ are the change of the particle numbers and the difference of the grand potentials respectively before and after the force protocol. This equality has been discussed in chemical reaction networks st2007 , exchange fluctuation systems qu2011 ; cl2004 ; th2009 ; sy2008 and isolated systems without particle number change wo2011 ; no2011 ; wo2012 . However fluctuation theorems of the grand canonical ensemble in isolated systems in nonequilibrium processes with particles creation and annihilation have not been studied so far (but see Ref. ja2018 ). The reason may be that usually the mass gap of massive particles is so large that the energy input during the driving of the boundary is too small to create massive particles.

Nevertheless, photons are massless particles, which makes it possible to create photons with relatively low energy input. One example is the dynamical Casimir effect (DCE). The DCE is a quantum effect which describes the generation of photons due to the EM field in the presence of time-dependent boundaries. The study was initiated by Moore qu1970 and followed by many researchers  fulling1976 ; ef1994 ; cole1995 ; ge1996 ; tr1998 ; no2001 ; re2001 ; qu2002 ; he2005 ; th2006 ; ph2006 ; dodonov2010 ; dalvit2011 ; nation2012 . The first experimental verification was carried out in 2011 johansson2010 ; ob2011 . The researchers used a modified SQUID to mimic a mirror moving at the required relativistic velocity and observed the DCE in a superconducting circuit. The main concern in the studies of DCE is about the average number of photon creation. The studies of the DCE inspire us to investigate the distributions of work and photon number change in the nonequilibrium processes with the photon creation and annihilation. It is worth mentioning that classical ideal gas insider a cylinder lua2005 has been a prototype model for the study of thermodynamics. EM field in a cavity is an analogue quan2009 ; ja2009 of gas insider a cylinder, but with more complicated dynamics, because it incorporates effects of quantum mechanics, quantum statistics, and special relativity. Nevertheless, for certain specific protocols, we are able to obtain the analytical results of the joint distribution of work and the change of photon number. The analytical results enable us to study the quantum-classical correspondence of the work done on the EM field in the high temperature limit and the quantum nature (Casimir effect) of the work in the low temperature limit.

In this paper, we focus on the DCE in a three-dimensional cavity with one oscillating boundary. Because the energy levels of the system are not equidistant, only a limited number of modes of the EM field are coupled for a specific oscillation frequency of the boundary. In addition, relevant to Fermi’s golden rule, there are only three kinds of resonance conditions, which significantly simplify the calculation of the characteristic function of the joint probability distribution. We further obtain the analytical expressions in the DCE with several different time-dependent geometries (the rectangular, cylindrical and spherical cavities with one moving boundary). Analytical expressions of work distributions for an arbitrary nonequilibrium process is extremely rare. This example, which incorporates effects of quantum mechanics, quantum statistics and special relativity has pedagogical value, and deepens our understanding of the quantum trajectory work and the validity of the fluctuation theorems in nonequilibrium processes with photon creation and annihilation.

We notice that Ref. ja2009 also discussed the Jarzynski equality for the photon gas. However, we believe that the particles considered in Ref. ja2009 are not real photons but some relativistic massless particles, since the wave character of light is ignored and no particle creation and annihilation occur during the force protocol.

This paper is organized as follows: In Sec. II, we introduce the effective Hamiltonian of a quantized EM field in a three-dimensional cavity with a moving boundary. In Sec. III, we clarifies the concept of trajectory work and the validity of the Jarzynski equality in the system. In Sec. IV, we obtain the analytical expressions of the characteristic functions by utilizing the matrix representation technique. And then we analyze the characteristic functions in the single-resonance and multiple-resonance conditions. We summarize this work and make conclusions in Sec. V.

II Effective Hamiltonian and resonance conditions for the EM field in a trembling cavity

We begin with the quantization of the EM field confined in a 3-D rectangular cavity, expressed in the scalar Hertz potentials he2005 . Then we derive the effective Hamiltonian, which describes the dynamics of the EM field when one of the boundaries is moving with time. Finally, we introduce some thermodynamic concepts relevant to the non-equilibrium driving process, where the photons are created or annihilated, including the two-point measurement and the characteristic function ${\mathcal{G}}(u,v)$. The study can be straightforwardly extended to the EM field in other geometries, such as the cylindrical and the spherical cavity (see Appendix A).

II.1 Hertz Potential Formalism

It is known that the EM field and the Maxwell’s equations can be formulated in terms of the scalar potential $\phi$ and the vector potential ${\mathbf{A}}$. Alternatively, the Hertz potentials, ${\mathbf{\Pi}}_{\rm e}$ and ${\mathbf{\Pi}}_{\rm m}$, offer an equivalent, more convenient for EM field in a cavity, formalism. In the Lorentz gauge, the relations between the two formulations can be written as follows,

[TABLE]

with $\varepsilon$ and $\mu$ being the permittivity and the permeability of the medium.

In source-free vacuum, the vector Hertz potentials become two scalar fields, i.e. ${\mathbf{\Pi}}_{\rm e}=\psi^{\rm TM}{\mathbf{e}}_{z}$ and ${\mathbf{\Pi}}_{m}=\psi^{\rm TE}{\mathbf{e}}_{z}$, where $\psi^{\rm TM}$ and $\psi^{\rm TE}$ are the transverse magnetic (TM) and the transverse electric (TE) field with respect to the longitudinal $z$ axis, with ${\mathbf{e}}_{z}$ being the unit vector along the $z$ axis. The Maxwell’s equations in the form of the scalar Hertz potentials can be written in the following form,

[TABLE]

with $\varepsilon_{0}$ and $\mu_{0}$ being the permittivity and the permeability of the vacuum.

Without loss of generality, we only consider the TE field in the following. The same procedure can be straightforwardly applied to the study of the TM field (see Appendix B). Also, we hereafter drop the superscript (TE) and set $\varepsilon_{0}=\mu_{0}=\hbar=1$ for simplicity, if not explicitly stated otherwise.

The Lagrangian density of the TE field is

[TABLE]

Meanwhile, the TE field in a 3-D cavity satisfies the following boundary conditions,

[TABLE]

where the boundaries of the cavity locate at 0 and $L_{\alpha}$, with $\alpha=x,y,z$.

II.2 Quantization of the scalar Hertz Potential

We quantize the scalar Hertz potential $\psi$ by promoting it from ordinary numbers to an operator $\hat{\psi}$ and imposing the following canonical commutation relations,

[TABLE]

where the conjugate momentum operator $\hat{\pi}({\mathbf{r}},t)$ is obtained from the Lagrangian density in Eq. (5),

[TABLE]

Solving the wave equation (4) under the boundary conditions (6), we obtain an orthonormal basis $\left\{\psi_{\mathbf{k}}({\mathbf{r}})\right\}$, where the basis state $\psi_{\mathbf{k}}({\mathbf{r}})$ is

[TABLE]

with the subscript ${\mathbf{k}}\equiv(k_{x},k_{y},k_{z})\in{\mathbb{N}}^{3}$, and the eigen-frequency of the ${\mathbf{k}}$-th mode is $\omega_{\mathbf{k}}=\sqrt{\left(\frac{k_{x}\pi}{L_{x}}\right)^{2}+\left(\frac{k_{y}\pi}{L_{y}}\right)^{2}+\left(\frac{k_{z}\pi}{L_{z}}\right)^{2}}$.

Then we expand the filed operator $\hat{\psi}({\mathbf{r}},t)$ and its conjugate momentum operator $\hat{\pi}({\mathbf{r}},t)$ as follows,

[TABLE]

where ${\mathbf{k}}_{\perp}=(k_{x},k_{y})$ and $\omega_{\mathbf{k}_{\perp}}=\sqrt{\left(\frac{k_{x}\pi}{L_{x}}\right)^{2}+\left(\frac{k_{y}\pi}{L_{y}}\right)^{2}}$, and $C_{\mathbf{k}}$ are the normalization constants. Note that $\hat{a}_{\mathbf{k}}$ ($\hat{a}_{\mathbf{k}}^{\dagger}$) is the annihilation (creation) operator for the ${\mathbf{k}}$-th mode, which satisfies the canonical commutation relations,

[TABLE]

Using Eqs. (5) and (10) and setting $C_{\mathbf{k}}\equiv\left(\sqrt{2\omega_{\mathbf{k}}}\omega_{{\mathbf{k}}_{\perp}}\right)^{-1}$, we obtain the Hamiltonian of the TE field as follows he2005 ,

[TABLE]

II.3 Effective Hamiltonian for Driving Processes

Now we consider that one of the boundaries of the cavity along the longitudinal $z$ direction is moving according to a prefixed time-dependent function $\lambda(t)$, i.e. $L_{z}(t)\equiv\lambda(t)$.

We first define an instantaneous orthonormal basis $\{\psi_{{\mathbf{k}},\lambda}({\mathbf{r}})\}$. The basis functions $\psi_{{\mathbf{k}},\lambda}({\mathbf{r}})$ satisfy the following Helmholtz equation,

[TABLE]

with the time-dependent boundary conditions qu1970

[TABLE]

where $\psi_{{\mathbf{k}},\lambda}({\mathbf{r}})$ have the same form as $\psi_{\mathbf{k}}({\mathbf{r}})$ in Eq. (9) except that $L_{z}$ is replaced by $\lambda$ and $\omega_{{\mathbf{k}},\lambda}=\sqrt{\left(\frac{k_{x}\pi}{L_{x}}\right)^{2}+\left(\frac{k_{y}\pi}{L_{y}}\right)^{2}+\left(\frac{k_{z}\pi}{\lambda}\right)^{2}}$. Note that both $\psi_{{\mathbf{k}},\lambda}({\mathbf{r}}$ and $\omega_{{\mathbf{k}},\lambda}$ depend on time when the boundary is moving with time. We then expand the field operators $\hat{\psi}({\mathbf{r}},t)$ and $\hat{\pi}({\mathbf{r}},t)$ with the instantaneous basis $\left\{\psi_{{\mathbf{k}},\lambda(t)}({\mathbf{r}})\right\}$,

[TABLE]

Taking time derivatives of Eqs. (15) and (16), we obtain the following equations of $\hat{Q}_{\mathbf{k}}$ and $\hat{\mathcal{P}}_{\mathbf{k}}$,

[TABLE]

where the coupling coefficients $\tilde{g}_{\mathbf{kp}}$ can be expressed as follows,

[TABLE]

with

[TABLE]

treating Eqs. (17) and (18) as the Heisenberg equations of motion, we obtain the following effective Hamiltonian,

[TABLE]

In order to move into the Fock representation, we introduce the following ladder operators,

[TABLE]

The first-order derivative of $\hat{a}_{\mathbf{k}}(t)$ can be written as follows,

[TABLE]

with $\gamma_{\mathbf{k}}(t)=-\frac{\dot{\omega}_{\mathbf{k},\lambda(t)}}{4\omega_{\mathbf{k},\lambda(t)}}$ and

[TABLE]

Treating Eq. (25) as the equation of motion, we obtain the effective Hamiltonian in terms of $\hat{a}_{\mathbf{k}}\equiv\hat{a}_{\mathbf{k}}(0)$ and $\hat{a}_{\mathbf{k}}^{\dagger}\equiv\hat{a}_{\mathbf{k}}^{\dagger}(0)$ as follows,

[TABLE]

This effective Hamiltonian (26) will be essential for our later analysis.

II.4 Perturbative Driving Processes

For simplicity, we consider the perturbative periodic driving protocol, in which the work parameter $\lambda(t)$ takes the following form,

[TABLE]

where $\Omega$ is the oscillation frequency, and the oscillation amplitude $\epsilon$ is assumed to be small, i.e. $\epsilon\ll 1$. We then expand all relevant physical quantities to the first order of $\epsilon$.

In the following, we turn to the interaction picture defined by the free Hamiltonian $\hat{H}_{0}=\sum_{\mathbf{k}}\omega_{\mathbf{k}}\left(\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\frac{1}{2}\right)$, with $\omega_{\mathbf{k}}\equiv\omega_{\mathbf{k},\lambda(0)}$, and treat $\hat{V}(t)=\hat{H}_{\rm eff}(t)-\hat{H}_{0}$ as the perturbation. Similar to Fermi’s golden rule, if $\hat{V}(t)$ is a periodic function of time with the angular frequency $\Omega$, the transition is into states with energies that differ by $\hbar\Omega$ from the energy of the initial state.

Using the rotating-wave approximation (RWA), we obtain time-independent Hamiltonians $\hat{V}_{j}^{\rm I}$ ($j=1,2,3$) in the interaction picture under the following three kinds of resonance conditions:

The Double-frequency (DoF) resonance: $\Omega=2\omega_{\mathbf{k}}$

[TABLE] 2. 2.

The Sum-frequency (SuF) resonance: $\Omega=\omega_{\mathbf{k}}+\omega_{\mathbf{p}}$

[TABLE] 3. 3.

The Difference-frequency (DiF) resonance: $\Omega=\left|\omega_{\mathbf{k}}-\omega_{\mathbf{p}}\right|$

[TABLE]

In addition, we can draw intuitions from Eqs. (28)-(30). The DoF resonance condition (28) corresponds to the process of simultaneously creating two photons with the same frequency $\omega_{\mathbf{k}}$, while the SuF resonance condition (29) corresponds to the process of creating one photon with the frequency $\omega_{\mathbf{k}}$ and another photon with the frequency $\omega_{\mathbf{p}}$. The DiF resonance condition (30) corresponds to the process of creating one photon with the frequency $\omega_{\mathbf{k}}$ meanwhile annihilate one photon with the frequency $\omega_{\mathbf{p}}$. It is worth mentioning that the evolution of other modes are quantum adiabatic with no photons being created or annihilated in these modes. Also, if the driving frequency $\Omega$ does not satisfy any of these resonance conditions (28)-(30), the evolution of the EM field is quantum adiabatic, i.e., no photons will be created or annihilated during the driving process.

III Trajectory work and the Jarzynski equality associated with the DCE

We are interested in the statistics of the work and the photon number difference of the EM field confined in a 3-D cavity driven by a perturbatively oscillating boundary. Before $t=0$ and after $t=\tau$, with $\tau$ being the driving duration, the boundary stops at $\lambda_{0}\equiv\lambda(0)$ and $\lambda_{\tau}\equiv\lambda(\tau)$ respectively, while it oscillates according to Eq. (27) when $t\in[0,\tau]$. During the driving process, both the work and the photon number fluctuate, which requires the grand-canonical description wo2011 .

When the boundary is fixed at $\lambda$, the Hamiltonian $\hat{H}_{\lambda}=\sum_{\mathbf{k}}\omega_{{\mathbf{k}},\lambda}\hat{a}_{{\mathbf{k}},\lambda}^{\dagger}\hat{a}_{{\mathbf{k}},\lambda}$ and the photon-number operator $\hat{N}_{\lambda}=\sum_{\mathbf{k}}\hat{a}_{{\mathbf{k}},\lambda}^{\dagger}\hat{a}_{{\mathbf{k}},\lambda}$ commute with each other, i.e. $\left[\hat{H}_{\lambda},\hat{N}_{\lambda}\right]=0$. Thus there exists a set of common eigenstates $\left|\left\{n_{\mathbf{k}}\right\}\right\rangle_{\lambda}\equiv\otimes_{\mathbf{k}}\left|n_{\mathbf{k}}\right\rangle_{\lambda}$, which satisfy

[TABLE]

where the total energy $E_{\lambda}\left(\left\{n_{\mathbf{k}}\right\}\right)=\sum_{\mathbf{k}}\omega_{{\mathbf{k}},\lambda}n_{\mathbf{k}}$ and the total photon number $N\left(\left\{n_{\mathbf{k}}\right\}\right)=\sum_{\mathbf{k}}n_{\mathbf{k}}$, with $n_{\mathbf{k}}\in{\mathbb{N}}$ being the number of photons in the ${\mathbf{k}}$-th mode. The density matrix and the partition function of the grand canonical ensemble of the photon gas is written as follows,

[TABLE]

with $\beta$ being the inverse temperature of the initial equilibrium state. Here we would like to emphasize that although the chemical potential of the photon gas is equal to zero, i.e. $\mu=0$, the total photon number is indefinite in thermal equilibrium.

To obtain the joint probability distribution of work and photon-number difference, ${\mathcal{P}}(w,\Delta N)$, the conceptual procedure of the two-point measurement protocol is prescribed as follows:

Prepare the system in the thermal equilibrium by connecting it with a heat bath at the temperature $\beta^{-1}$. Then remove the heat bath so that the system is isolated. 2. 2.

Perform the first projective measurements of $\hat{H}_{\lambda_{0}}$ and $\hat{N}_{\lambda_{0}}$. Project the system to one of the common eigenstates, i.e. $\left|\left\{n_{\mathbf{k}}\right\}\right\rangle_{\lambda_{0}}$, and record the eigenvalues $E_{\lambda_{0}}\left(\left\{n_{\mathbf{k}}\right\}\right)$ and $N\left(\left\{n_{\mathbf{k}}\right\}\right)$. 3. 3.

Control the boundary according to Eq. (27) for a prefixed duration $\tau$. The frequency $\Omega$ is chosen such that one of the resonance conditions (28–30) is satisified, 4. 4.

Perform the second projective measurements of $\hat{H}_{\lambda_{\tau}}$ and $\hat{N}_{\lambda_{\tau}}$ and record the eigenvalues $E_{\lambda_{\tau}}\left(\left\{n^{\prime}_{\mathbf{k}}\right\}\right)$ and $N\left(\left\{n^{\prime}_{\mathbf{k}}\right\}\right)$.

Ideally, the above procedure is repeated infinitely many times to obtain a good statistics of the joint probability distribution, which is defined as

[TABLE]

with $\delta(\cdot)$ and $\delta_{\cdot,\cdot}$ being the Dirac and the Kronecker delta functions, respectively.

Alternatively, we can calculate the characteristic function of the work and the photon-number difference, which is defined as the Fourier transform of ${\mathcal{P}}(w,\Delta N)$,

[TABLE]

where the Heisenberg-picture operators are defined as $\hat{O}^{\rm H}(t)=\hat{U}^{\dagger}(t)\hat{O}\hat{U}(t)$ and the ensemble average is taken over $\hat{\rho}_{\beta}$ in Eq. (32), i.e. $\left\langle\cdot\right\rangle_{\beta}={\rm Tr}\left[\cdot\hat{\rho}_{\beta}\right]$. Then the joint probability distribution ${\mathcal{P}}(w,\Delta N)$ can be obtained via the inverse Fourier transform.

Although there exists no Schördinger picture description in our system during the driving process qu1970 , the Hamiltonian and the total photon number operators in the Heisenberg picture before and after the driving process are well-defined. Especially at the end of the driving process, they are defined as follows,

[TABLE]

where the time-dependent ladder operators $\hat{a}_{\mathbf{k}}(t)$ at the beginning ($t=0$) and the end ($t=\tau$) of the driving process are related via the unitary transformation defined by $\hat{H}_{\rm eff}$ in Eq. (26),

[TABLE]

with $V^{\rm I}$ being one of $V_{j}^{\rm I}$ in Eqs. (28–30), depending on the type of the resonance condition satisfied in the driving process.

The Jarzynski equality for a grand canonical ensemble,

[TABLE]

is obtained by setting $u=i\beta$ and $v=-i\beta\mu$, where we have defined the grand potential difference $\Delta\Phi=\Phi_{\lambda_{\tau}}-\Phi_{\lambda_{0}}$ with the grand potential $\Phi_{\lambda_{\tau,0}}=-\beta^{-1}\mathrm{ln}{\rm Tr}\left[e^{-\beta\left(\hat{H}_{\lambda_{\tau}}-\mu\hat{N}_{\lambda_{\tau}}\right)}\right]$.

The equality in Eq. (37) reproduces the original Jarzynski equality, $\langle e^{-\beta w}\rangle=e^{-\beta\Delta\Phi}$, in the following two cases

[TABLE]

We would like to emphasize that Refs. no2011 ; wo2011 ; wo2012 have discussed systems that belong to case (i), while our system belongs to case (ii).

IV analytical results of work distributions under various resonance conditions

Generally, it is impossible to obtain closed-form expressions of ${\mathcal{G}}(u,v)$ for an arbitrary non-equilibrium driving processes of a quantum many-body system, let alone the case with particle creation and annihilation. However, we find that for the photon gas confined in a cavity subject to a perturbative resonant driving (27), the characteristic function ${\mathcal{G}}(u,v)$ can be written into the expectation value of a product of a series of exponentials of a quadratic form in the Boson operators. By substituting Eqs. (35) and (36) into Eq. (III) and then by applying the matrix representation technique en1978 , we obtain the analytical results of the characteristic function ${\mathcal{G}}(u,v)$ (for more technical details see Appendix C).

IV.1 Characteristic Functions in Single-Resonance Conditions

If $\Omega$ is chosen such that one of the three resonance conditions (28)-(30) is satisfied, the perturbation Hamiltonian becomes independent of time (see Eqs. (28–30)). In these cases, we can analytically calculate the characteristic function ${\mathcal{G}}(u,v)$ in Eq. (III). The results are given below:

The DoF resonance: $\Omega=2\omega_{\mathbf{k}}$

[TABLE] 2. 2.

The SuF resonance: $\Omega=\omega_{\mathbf{k}}+\omega_{\mathbf{p}}$

[TABLE] 3. 3.

The DiF resonance: $\Omega=\left|\omega_{\mathbf{k}}-\omega_{\mathbf{p}}\right|$

[TABLE]

Note that we restrict ourselves to the cases with $\lambda_{0}=\lambda_{\tau}$ for simplicity. The analytic solutions for more general cases with $\lambda_{0}\neq\lambda_{\tau}$ are more cumbersome and can be found in Appendix D.

Note that the characteristic function ${\mathcal{G}}(u,v)$ for an externally-driven closed quantum system should satisfy the following requirements:

The normalization of the joint probability ${\mathcal{P}}(w,\Delta N)$, i.e. $\sum_{\Delta N}\int dw{\mathcal{P}}(w,\Delta N)=1$, requires that ${\mathcal{G}}(0,0)=1$. 2. 2.

The grand-canonical quantum Jarzynski equality, i.e. $\left\langle e^{-\beta\left(w-\mu\Delta N\right)}\right\rangle=e^{-\beta\Delta\Phi}$, requires that ${\mathcal{G}}(i\beta,-i\beta\mu)=e^{-\beta\Delta\Phi}$. 3. 3.

The grand-canonical Crooks’ fluctuation theorem, i.e. $\frac{{\mathcal{P}}_{\rm F}\left(w,\Delta N\right)}{{\mathcal{P}}_{\rm R}\left(-w,-\Delta N\right)}=e^{\beta\left(w-\nu\Delta N-\Delta\Phi\right)}$, where the forward and the reverse processes are respectively described by $\lambda(t)$ and $\lambda(\tau-t)$, requires that ${\mathcal{G}}_{\rm R}(-u,-v)={\mathcal{G}}_{\rm F}(u+i\beta,v-i\beta\mu)e^{\beta\Delta\Phi}$. 4. 4.

The discreteness of the distribution function of work for a closed quantum system requires that ${\mathcal{G}}(u,v)$ is periodic in $u$. More specifically, a delta-function peak $\delta(w-w_{0})$ in ${\mathcal{P}}(w)$ implies that ${\mathcal{G}}\left(u+\frac{2\pi}{\omega_{0}},v\right)={\mathcal{G}}(u,v)$. 5. 5.

The indivisibility of the constituent photons, $\Delta N\in{\mathbb{Z}}$, requires that ${\mathcal{G}}(u,v+2\pi)={\mathcal{G}}(u,v)$. Furthermore, if only $n>0$ photons can be created or annihilated simultaneously, ${\mathcal{G}}\left(u,v+\frac{2\pi}{n}\right)={\mathcal{G}}(u,v)$. Finally, if the photon creation and annihilation processes are prohibited, ${\mathcal{G}}(u,v)$ is independent of $v$, i.e. ${\mathcal{G}}(u,v)={\mathcal{G}}(u)$.

For the particular situation we are interested in, the EM field confined in a 3-D cavity resonantly driven by a trembling boundary, the chemical potential vanishes and the positions of the trembling boundary at the initial and final time are exactly the same. It can be checked that Eqs. (39–41) satisfy all the above requirements.

IV.2 Marginal Distributions of the work and the photon-number difference

The joint distribution function ${\mathcal{P}}(w,\Delta N)$ can be obtained from the inverse Fourier transform of the characteristic function ${\mathcal{G}}(u,v)$,

[TABLE]

from which we obtain the following marginal distributions ${\mathcal{P}}(w)$ and ${\mathcal{P}}(\Delta N)$,

[TABLE]

In Fig. 1, we show the marginal distributions for all of the three kinds of resonance conditions. Since we are considering closed systems with discrete energy levels, all of the distributions are discrete, consisting of a series of the Dirac delta-functions with different weights. The work performed on the system is always an integer multiple of $2\omega_{\mathbf{k}}$, $\omega_{\mathbf{k}}+\omega_{\mathbf{p}}$ or $\omega_{\mathbf{k}}-\omega_{\mathbf{p}}$, when the DoF (28), the SuF (29) or the DiF (30) resonance condition is satisfied, respectively. As to the photon number difference $\Delta N$, the photons are created and annihilated in pairs in the first two cases. And the numbers of the photons created and annihilated are equal to each other in the DiF resonance condition footnote1 . These features are consistent with our understandings of the three resonance conditions in Eqs. (28–30).

IV.3 Classical Limit of the Work Distributions

The classical characteristic functions can be obtained from the closed-form expressions in Eqs. (39–41) by taking the classical limit, i.e. $\hbar\rightarrow 0$ footnote3 . Similar to Ref. wo2011 , we introduce $\tilde{u}=u/\beta$.

The classical characteristic functions for the three resonance conditions can be written as follows footnote5 ,

[TABLE]

with the ratio being defined as $r\equiv\omega_{\mathbf{p}}/\omega_{\mathbf{k}}$ in the latter two cases. Note that the classical work distributions for the latter two cases can be obtained analytically as follows,

[TABLE]

with $j=2,3$, where $\Theta(\cdot)$ is the Heaviside step function and

[TABLE]

Fig. 2 shows the consistency between the semiclassical work distributions and the exact ones when the temperature of the initial state is high. Also, we show the discrepancy between the exact marginal distributions of work and the Gaussian fitting with the same average value and the same standard deviation of the work distribution. It can be seen that the work distributions obviously deviate from the Gaussian distribution.

IV.4 Quantum-to-Classical Crossover

The analytic expressions of the average value and the standard deviation of the work distributions are defined as

[TABLE]

[TABLE]

which can be obtained from Eqs. (39–41). For the convenience of analysis, we explicitly write down the Planck constant $\hbar$. It can be anticipated that when $\hbar\rightarrow 0$, the system loses all of its quantum features footnote2 .

In Fig. 3, we demonstrate the average work performed on the system in these three resonance conditions. Inspecting the figures, we find that the average work decreases to zero as $\hbar$ vanishes in the low-temperature limit for the DoF and the SuF cases, which is consistent with the fact that no work can be done on the classical vacuum. Also, the dependence of $\left\langle w\right\rangle_{j}$ on $\hbar$ becomes weaker for all three resonance conditions as the temperature increases, which implies a crossover from the quantum to the classical regime.

Inspired by the observations in Fig. 3, we list in TABLE 1 the average value and the standard deviation of the work distributions for the three resonance conditions in two limits, i.e. the low-temperature ($\beta\rightarrow\infty$) and the high-temperature ($\beta\rightarrow 0$) limits. It can be seen clearly that the energy scale is $\hbar\omega$ and $\beta^{-1}$ in the low-temperature and the high-temperature limits, respectively. In the high temperature limit, the work is equal to the work done on a classical EM field, which can be seen from the fact that the work value does not depend on $\hbar$.

IV.5 Multiple-Resonance cases

If $\Omega$ is chosen such that for all modes in the cavity field, at least two resonance conditions are satisfied simultaneously, we are dealing with a “multiple-resonance” condition. The multiple-resonance condition can be double-resonant, triple-resonant, and so on.

For the multiple-resonance condition, we consider the following two cases: the uncoupled and the coupled. For the uncoupled case, the involved modes are different, e.g. $\Omega=2\omega_{\mathbf{k}}=\omega_{\mathbf{s}}-\omega_{\mathbf{p}}$ with $\mathbf{k}\neq\mathbf{s}\neq\mathbf{p}$. Because the resonance conditions are mutually independent, the characteristic function can be decomposed into a product of characteristic functions, each of which can be obtained analytically in the corresponding resonance condition. Let us denote the characteristic function for a multiple-resonance case as ${\mathcal{G}}_{\mathbf{k}_{1},\mathbf{k}_{2},\cdots,\mathbf{s}_{1}\pm\mathbf{p}_{1},\mathbf{s}_{2}\pm\mathbf{p}_{2},\cdots}(u,v)$ for resonance conditions $\Omega=2\omega_{\mathbf{k}_{1}}=2\omega_{\mathbf{k}_{2}}=\cdots=\omega_{\mathbf{s}_{1}}\pm\omega_{\mathbf{p}_{1}}=\omega_{\mathbf{s}_{2}}\pm\omega_{\mathbf{p}_{2}}=\cdots$. The characteristic function for this case reads

[TABLE]

and the joint probability distribution function ${\mathcal{P}}(w,\Delta N)$ can be obtained by the inverse Fourier transform.

For the coupled case, more than one resonance conditions involve the same modes, for example $\Omega=2\omega_{\mathbf{k}}=\left|\omega_{\mathbf{k}}-\omega_{\mathbf{p}}\right|$. Since the mode $\mathbf{k}$ appears in both resonance conditions, these two resonance conditions cannot be considered separately. Instead, the perturbation Hamiltonian in the interaction picture becomes $\hat{V}^{\mathrm{I}}=-ig_{1}(\hat{a}_{\mathbf{k}}^{{\dagger}}\hat{a}_{\mathbf{k}}^{\dagger}-\hat{a}_{\mathbf{k}}\hat{a}_{\mathbf{k}})-ig_{3}(\hat{a}_{\mathbf{k}}^{{\dagger}}\hat{a}_{\mathbf{p}}-\hat{a}_{\mathbf{p}}^{{\dagger}}\hat{a}_{\mathbf{k}})]$, with which the characteristic function ${\mathcal{G}}(u,v)$ can also be obtained by the matrix representation technique. The discussions can be straightforwardly extended to the multiple-resonance case with more than two resonance conditions are satisfied simultaneously.

V Conclusions

To calculate the work distribution in an arbitrary nonequilibrium process in a quantum many-body system is usually very cumbersome due to the interplay of effects of quantum mechanics and quantum statistics gong2014 . In very rare situations, one is able to obtain the analytical solution to the distribution of work. These analytical results deepens our understanding of quantunm trajectory work and fluctuation theorems in the nonequilibrium processes.

In our current study, we investigated the work distribution of a quantized EM field in a three-dimensional cavity with an oscillating boundary. This system incorporates the effects of not only quantum mechanics and quantum statistics, but also special relativity. For the periodic perturbative driving protocol (27), under the RWA, we obtained the effective Hamiltonian in the interaction picture. Also, we analytically evaluate the characteristic function ${\mathcal{G}}(u,v)$ in the single-resonance (39–41) and multiple-resonance conditions using the matrix representation technique. If $\Omega$ is chosen such that none of the single resonance conditions is satisified, the evolution of the EM field is quantum adiabatic. We discussed the general properties of ${\mathcal{G}}(u,v)$ and verified various fluctuation theorems in the nonequilibrium processes with photon creation and annihilation. From the analytical result of the work distribution, we can clearly see that nonzero work is done at zero temperature, which is a manifestation of the quantum nature (Casimir effect) of the EM field. However, the work vanishes when $\hbar\to 0$, which agrees with our intuition that no work is done on the classical vacuum when the boundary is driven. We also obtained the approximate expression of the work distribution ${\mathcal{P}}(w)$ and the moments of work at high temperature, which is consistent with the work done on a classical EM field. Our study has pedagogical value because analytical solutions to the work distribution in a quantum many-body system is very rare. Last but not least, the dynamical Casimir effect has been experimentally tested in a superconducting circuit johansson2010 ; ob2011 . Hopefully, our theoretical predictions about the work distributions and the validity of the Jarzynski equality can be experimentally tested in the superconducting circuit ob2011 in the future.

Acknowledgements.

H. T. Quan gratefully acknowledges support from the National Science Foundation of China under grants 11775001, 11534002, and The Recruitment Program of Global Youth Experts of China. Jing-Ning Zhang gratefully acknowledges the support from the National Natural Science Foundation of China (Grants No. 11504197).

Appendix A: EM field in various geometries

Our results of the EM field in a 3-D trembling rectangular cavity can be generalized to the EM field in a trembling cylindrical or spherical cavity. We skip the derivation which is similar to that of the rectangular cavity and list the main results here. The work parameter $\lambda(t)$ takes the form $\lambda(t)=\lambda_{0}[1+\epsilon\mathrm{sin}(\Omega t)]$ in the following, and we take the first-order approximation when doing expansion in $\epsilon$.

In a cylindrical cavity, we introduce the cylindrical coordinates $(\rho,\phi,z)$. If the boundaries are at the radial coordinate $\rho=R$ and the longitudinal coordinate $z=(0,\lambda(t))$, we have

Hertz potentials

[TABLE] 2. 2.

Boundary conditions

[TABLE]

[TABLE] 3. 3.

Instantaneous orthonormal bases

[TABLE]

[TABLE]

where $J_{n}$ denotes the Bessel function of the first kind of the $n$th order, $y_{nm}$ is the $m$th positive root of the equation $J^{\prime}_{n}(y)=0$, and $x_{nm}$ is the $m$th root of the equation $J_{n}(x)=0$. 4. 4.

Eigenvalues

[TABLE]

[TABLE] 5. 5.

Coupled strengths

[TABLE]

[TABLE]

Similarly if the boundaries are at the radial coordinate $\rho=\lambda(t)$ and the longitudinal coordinate $z=(0,L_{z})$, we have

Boundary conditions

[TABLE]

[TABLE] 2. 2.

Instantaneous orthonormal bases

[TABLE]

[TABLE] 3. 3.

Eigenvalues

[TABLE]

[TABLE] 4. 4.

Coupled strengths

[TABLE]

[TABLE]

In a spherical cavity, we introduce the spherical coordinates $(r,\theta,\phi)$. If the boundary is at the radial coordinate $r=\lambda(t)$, we have

Hertz potentials

[TABLE]

where we have used the Debye potentials ph2006 . 2. 2.

Boundary conditions

[TABLE]

[TABLE] 3. 3.

Instantaneous orthonormal bases

[TABLE]

[TABLE]

where $j_{l}$ denotes the spherical Bessel function of the $l$th order, $Y_{lm}(\theta,\phi)$ denotes the normalized spherical harmonics of degree $l$ and order $m$, $j_{ln}$ denotes the $n$th zero for $j_{l}(x)=0$, and $\kappa_{ln}$ denotes the $n$th zero of $\partial_{x}[xj_{l}(x)]=0$. 4. 4.

Eigenvalues

[TABLE]

[TABLE] 5. 5.

Coupled strengths

[TABLE]

[TABLE]

Appendix B: TM field in a trembling cavity

The TM field in the trembling cavity satisfies different boundary conditions from the TE field. According to Refs. qu2002 ; he2005 , the boundary conditions are

[TABLE]

Similar to the TE field, we also define the following instantaneous orthonormal basis $\{\psi^{TM}_{\mathbf{k},\lambda}({\mathbf{r}})\}$

[TABLE]

which satisfies the Helmholtz equation (Eq. 13) and the boundary conditions (Eq. (LABEL:ea1)) to the first order of $\epsilon$ (notice $\left.\dot{\lambda}(t)\partial_{t}\psi^{\mathrm{TM}}_{\mathbf{k}}\right|_{z=\lambda(t)}\sim\epsilon^{2}$). Following the same procedure, we obtain the expression of the effective Hamiltonian, which is the same as that of the TE field (Eq. (26)) except

[TABLE]

Appendix C: the matrix representation technique

The matrix representation technique is a mathematical technique that can be used to calculate the trace of the products of several exponentials of a quadratic form in Boson operators. For simplicity, we first consider the product of two exponentials, i.e. $\hat{J}_{1}\hat{J}_{2}$, where $\hat{J}_{i},i=1,2$ is an exponential of a quadratic form in Boson operator

[TABLE]

Here, $\hat{\bm{\alpha}}=(\hat{a}_{1},\hat{a}_{2},\ldots,\hat{a}_{n},\hat{a}^{{\dagger}}_{1},\hat{a}^{{\dagger}}_{2},\ldots,\hat{a}^{\dagger}_{n})$ and $\hat{a}_{j},j=1,\ldots,n$ satisfies the bosonic commutation relations and the ($2n\times 2n$) matrix $S_{i}$ is a complex symmetric matrix. For later convenience, we introduce a characteristic matrix $[\hat{J}_{i}]$ corresponding to $\hat{J}_{i}$.

[TABLE]

where

[TABLE]

and $I$ is the $n\times n$ identity matrix. Now, let us define a set $\mathcal{\hat{J}}$ which includes all operators of the type (Eq. (C.1)) and the arbitrary product of these operators. Also, we define a set $\mathcal{J}$ which includes all matrices of the type (Eq. (C.2)) and the arbitrary product of these matrices. Then, it can be proved that both sets are the representations of the $2n$-dimensional complex symplectic group with operator and matrix multiplications respectively no1969 . Also, the product of operators is related to the matrix multiplication, i.e. if $\hat{J}_{3}=\hat{J}_{1}\hat{J}_{2}$, we have $[\hat{J}_{3}]=[\hat{J}_{1}][\hat{J}_{2}]$, where $\hat{J}_{i}\in\mathcal{\hat{J}}$, and $[J_{i}]\in\mathcal{J},i=1,2,3.$ This result can be straightforwardly generalized to the product of more than two exponentials. So, the product of several exponentials of a quadratic form in Boson operators can be related to one exponential of a quadratic form in Boson operator of which the characteristic matrix is known. Finally, the trace of an exponential of a quadratic form in Boson operator can be calculated by its characteristic matrix en1978 by

[TABLE]

where

[TABLE]

Thus, we transform the trace of the product of several exponentials of a quadratic form in Boson operators to the calculation of the characteristic matrices which is easy to deal with.

Appendix D: analytical solutions for the general case with $\lambda_{0}\neq\lambda_{\tau}$

The matrix representation technique can also be applied to more general cases in which the trembling boundary starts and stops at different positions $\lambda_{0}\neq\lambda_{\tau}$. In the single-resonance conditions, the modified characteristic function ${\overline{\mathcal{G}}}(u,v)={\mathcal{G}}(u,v)e^{-iu\Delta\Phi}$ reads

The DoF resonance: $\Omega=2\omega_{\mathbf{k}}$

[TABLE] 2. 2.

The SuF resonance: $\Omega=\omega_{\mathbf{k}}+\omega_{\mathbf{p}}$

[TABLE] 3. 3.

The DiF resonance: $\Omega=\left|\omega_{\mathbf{k}}-\omega_{\mathbf{p}}\right|$

[TABLE]

where $\Delta\omega_{\bm{s}}=\omega_{\bm{s},\lambda_{\tau}}-\omega_{\bm{s},\lambda_{0}},\bm{s}={\mathbf{k}},{\mathbf{p}}.$

It is worth mentioning that the Crook’s fluctuation theorem in these cases is non-trivial. For later convenience, we consider the moving boundary takes a more general form

[TABLE]

where $\varphi$ denotes the initial phase of the boundary. In this case, the effective Hamiltonian after the RWA becomes $e^{-\frac{i\varphi}{2}\hat{N}_{\lambda_{0}}}\hat{V}^{\mathrm{I}}e^{\frac{i\varphi}{2}\hat{N}_{\lambda_{0}}}$. The characteristic function will not change after the unitary transformation because of the cyclic property of trace.

In the nonzero initial phase $\varphi$ cases, the modified characteristic functions associated with the forward process $\lambda(t)$ and the reverse process $\lambda(\tau-t)$ read

[TABLE]

[TABLE]

where $f|_{a\leftrightarrow b}$ means the value of the parameters $a$ and $b$ in function $f$ is interchanged. Then, we have ${\mathcal{G}}_{\rm R}(-u,-v)={\mathcal{G}}_{\rm F}(u+i\beta,v-i\beta\mu)e^{\beta\Delta\Phi}$, which is equivalent to the Crook’s fluctuation theorem $\frac{{\mathcal{P}}_{\rm F}\left(w,\Delta N\right)}{{\mathcal{P}}_{\rm R}\left(-w,-\Delta N\right)}=e^{\beta\left(w-\nu\Delta N-\Delta\Phi\right)}$ and the time reversal symmetry of the effective Hamiltonian $\hat{H}_{\mathrm{eff}}(t)$ th2009

[TABLE]

where the time-reversal operator $\mathcal{\hat{T}}$ is an anti-linear and anti-unitary operator and which implies changing the sign of all odd operators. We would like to emphasize that $\mathcal{\hat{T}}^{-1}\hat{\psi}^{\rm TE}(\bm{r},t)\mathcal{\hat{T}}=-\hat{\psi}^{\rm TE}(\bm{r},-t)$ and $\mathcal{\hat{T}}^{-1}\hat{\psi}^{\rm TM}(\bm{r},t)\mathcal{\hat{T}}=\hat{\psi}^{\rm TM}(\bm{r},-t)$, which indicate that the electric field $\bm{E}$ (the magnetic field $\bm{B}$) is an even (odd) operator. It is interesting to notice that the time derivative of the work parameter $\dot{\lambda}(t)$ also appears in the effective Hamiltonian in Ref. re2017 .