Optical simulation of a quantum thermal machine
M. H. M. Passos Alan C. Santos Marcelo S. Sarandy J. A. O. Huguenin





mhmpassos@id.uff.br
Instituto de Ciências Exatas, Universidade Federal Fluminense, Volta Redonda, Rio de Janeiro, Brazil
Instituto de Física, Universidade Federal Fluminense, Avenida General Milton Tavares de Souza s/n, Gragoatá, 24210-346 Niterói, Rio de Janeiro, Brazil


ac_santos@id.uff.br
Instituto de Física, Universidade Federal Fluminense, Avenida General Milton Tavares de Souza s/n, Gragoatá, 24210-346 Niterói, Rio de Janeiro, Brazil


msarandy@id.uff.br
Instituto de Física, Universidade Federal Fluminense, Avenida General Milton Tavares de Souza s/n, Gragoatá, 24210-346 Niterói, Rio de Janeiro, Brazil


jose_huguenin@id.uff.br
Instituto de Ciências Exatas, Universidade Federal Fluminense, Volta Redonda, Rio de Janeiro, Brazil
Instituto de Física, Universidade Federal Fluminense, Avenida General Milton Tavares de Souza s/n, Gragoatá, 24210-346 Niterói, Rio de Janeiro, Brazil


We introduce both a theoretical and an experimental scheme for simulating a quantum thermal engine through an all-optical approach, 
with the behavior of the working substance and the thermal reservoirs implemented via internal degrees of freedom of a single-photon.
By using the polarization and propagation path, we encode two quantum bits and then implement the thermodynamical steps of an 
Otto cycle. To illustrate the feasibility of our proposal, we experimentally realize such simulation through
an intense laser beam, evaluating heat and work at each individual step of the  thermodynamical cycle. 
In addition, from the analysis of the entropy production during the entire cycle, we can study the amount of quantum friction 
produced in the Otto cycle as a function of the difference of temperature between hot and cold reservoirs. 
Our investigation constitutes, therefore, an all-optical-based thermal machine simulation and opens perspectives for 
other optical simulations in quantum thermodynamics.


Optical simulation of a quantum thermal machine
    J. A. O. Huguenin
    December 30, 2023
===============================================





§ INTRODUCTION


The idealization of microscopic quantum systems allowing for extraction of work and heat is at the heart of quantum thermal engines (QTEs), quantum refrigerators <cit.>, and quantum batteries <cit.>. In analogy with its classical counterpart, a QTE has a quantum system as its working substance, which interacts with thermal reservoirs at different temperatures β_c and β_h. Indeed, a number of works have explored the quantum nature of the working substance in order to investigate whether we can get optimal performance of  QTEs in comparison with their classical counterparts. Concrete proposals of QTEs have been studied in recent years from different approaches <cit.>.  Moreover, experimental verifications of the performance of a QTE 
have recently been achieved, e.g., nuclear spin systems manipulated through nuclear magnetic resonance (NMR) <cit.> and nitrogen-vacancy centers in 
diamond <cit.>. In general, a major difficulty for implementing QTEs in real physical systems is the high controllability required so that robustness against decoherence is achieved. Therefore, there is great interest in designing QTEs from architectures that offer efficient control of reservoirs.

In order to simulate controllable reservoirs, we have to consider the effect of quantum channels acting on quantum information <cit.>. In this context, it is 
fundamental in our approach to take into account the optical implementation of relevant quantum channels, such as amplitude damping, phase-damping (PD), and bit flip channels, among others performed by using single photons <cit.>. On the other hand, degrees of freedom of an intense laser beam have been widely used to simulate single-photon experiments and the results show that such procedure consists in a relevant test-bed for several quantum properties in a rather simple way <cit.>. 
Indeed, it can be shown that such systems can be used to observe violations of Bell's like inequality <cit.> and Mernin's inequality for tripartite systems <cit.>.  Moreover, many other quantum protocols can be investigated, such as quantum key distribution <cit.>, teleportation <cit.> and quantum logical gates <cit.>. As a further implementation of interest here, it is important to highlight the experimental simulation of open quantum systems to investigate environment-induced entanglement <cit.>.
In this paper, we propose an all-optical-based scheme, which allows us to simulate the performance of a thermal machine in quantum mechanics and perform an experimental simulation by using degrees of freedom of an intense laser beam. 
We theoretically show how we can construct a quantum machine by using the dephasing channel to mimic a thermal reservoir. 
We then  implement the Otto cycle for polarization of a single-photon through a simulation via linear optical circuit.

 

§  QUANTUM THERMAL MACHINE





 §.§ Quantum Otto cycle


Let us begin by the definition of heat and work in quantum mechanics. In general, heat and work are not 
quantum observables <cit.>. However, 
for the processes of interest here, either heat or work will be vanishing. In this situation, convenient expressions can be derived from the first law of thermodynamics. 
Indeed, by considering the internal energy  U(t) of a quantum system described by a density operator ρ(t) at instant t as U(t) = H(t)ρ(t), with 
H(t) denoting the Hamiltonian of the system, it is possible to define work δ W(t) and heat δ Q(t) for infinitesimal processes as <cit.>

    δ W(t)     =    Ḣ(t)ρ(t)dt  , 
    δ Q(t)     =    H(t)ρ̇(t)dt  .

Notice that δ W(t) > 0 (δ W(t) < 0) implies that work is being performed on (by) the system, so that its internal energy is increasing (decreasing). 
Similarly, when δ Q(t) > 0 (δ Q(t) < 0) we say that heat is being injected in (extracted from) the system. More details are presented in Appendix <ref>.

We will consider a quantum thermal machine realizing an Otto cycle, whose steps are shown in Fig. <ref> and can be described as follows.

Gap expansion step – Initially a quantum bit (qubit) is prepared in a thermal state of the reference Hamiltonian H_e(0) = ω_iniσ_y, at inverse temperature β_c. Thus, the system undergoes a unitary dynamics driven by the Hamiltonian H_e(t) = ħ[ ω_ini f(t) + ω_fin g(t) ]σ_y, with functions {g,f}:t∈→ g,f∈ satisfying g(0)=f(τ) = 1 and g(τ)=f(0) = 0 and |ω_ini|<|ω_fin|.  This constitutes the branch A→ B of the thermal engine cycle  in Fig. <ref>, where an amount of work W_A→B is performed by the engine.

Thermalization with hot reservoir – At this stage, the system is coupled to a thermal reservoir at inverse temperature β_h, thermalizing with it. Therefore, the final state at this step is a thermal state of the Hamiltonian H_e(τ) at inverse temperature β_h. In this step, the system exchange heat with the reservoir, but no work is performed. In this branch B→ C  (see Fig. <ref>), heat Q_B→C is introduced into engine.

Gap compression step – Now, we switch off the interaction between the system and the reservoir. Thus, we drive the system by a time-dependent Hamiltonian H_c(t) = ħ[ ω_fin f(t) + ω_ini g(t) ]σ_y. In this branch C→ D  (see Fig. <ref>), an amount of work  W_C→D is being performed on the system.

Thermalization with cold reservoir – To end, the system is coupled to the cold reservoir, in which the final state is the thermal state of the Hamiltonian H_c(τ) at inverse temperature β_c. No work is performed, but heat is exchanged between the system and the cold reservoir. This is the last branch D→ A 
of the cycle  (see Fig. <ref>), with an amount of heat Q_D→A extracted from the engine.
From Eqs. (<ref>) and (<ref>) we can compute heat and work for each step of this Otto cycle (shown in Fig. <ref>) as

    W_A→B   =    - ħ (ω_fin -ω_ini)tanh (ħω_iniβ_c)  , 
    
    Q_B→C   =   ħω_fin[ tanh (ħω_iniβ_c) - tanh (ħω_finβ_h) ]  , 
    
    W_C→D   =   ħ (ω_fin -ω_ini)tanh (ħω_finβ_h)  , 
    
    Q_D→A   =    - ħω_ini[ tanh (ħω_iniβ_c) - tanh (ħω_finβ_h) ]  ,

where we can derive the condition between β_c and β_h in order to get Q_D→A < 0 as ω_iniβ_c > ω_finβ_h. Such condition establishes the relation between the parameters of the reservoir and Hamiltonian as T_h/T_c > ω_fin/ω_ini.





 §.§ Simulation of a thermal reservoir with the phase damping channel


We can simulate the required reservoirs for implementing a quantum  thermal machine by using a PD quantum channel for a single-photon. To see this, we first need to realize that our system is initially prepared in a thermal state of H_ini at inverse temperature β_c, which reads 

    ρ_ini^th = 1/2[ - tanh (ħω_iniβ_c)σ_y]  .


Thus, due to the contact of our system with a thermal reservoir at inverse temperature β_h, under action of the Hamiltonian H_fin, the state after thermalization will be 

    ρ_fin^th = 1/2[  - tanh (ħω_finβ_h)σ_y]  .

Thus, the reservoir just changes the off-diagonal elements of the initial state, from tanh (ħω_iniβ_c) to tanh (ħω_finβ_h). On the other hand, given a density matrix ρ with elements ρ_nm, we know that a  PD channel acts over the elements ρ_01 and ρ_10 as ρ_01→ρ_01 e^-γτ_d and ρ_10→ρ_10 e^-γτ_d, respectively, where γ is the dephasing rate and τ_d is the total time interval of interaction of our system with the decohering reservoir <cit.>. To conclude, by applying this map to the state ρ_ini^th we have tanh (ħω_iniβ_c) → e^-γτ_dtanh (ħω_iniβ_c), so that we can adjust the parameter γτ_d to get the parameter β_h from

    ħω_finβ_h = arctanh[e^-γτ_dtanh (ħω_iniβ_c) ]  .

Thus, one can use the phase-damping channel to simulate the thermal reservoir in a heat engine, where we set the parameter γτ_d 
to encode the hot reservoir temperature.  It is important to mention here that the state ρ_fin^th  in Eq. (<ref>), obtained  after imposing ħω_finβ_h as given by Eq. (<ref>), 
does not represent  an actual thermal (Gibbs) state, since we do not have  actual thermal baths 
 in contact with the quantum system. Its correct meaning should be understood as a simulated thermal state, 
 which is achieved by mapping the desired temperature in terms of the dephasing parameters.

In several schemes of QTEs <cit.>, both steps of compression and expansion are performed by slow (adiabatic) unitary evolution, so that an amount of work is performed on/by the system and no heat is exchanged. However, since any unitary dynamics suppresses the heat exchange  (closed system), we can implement a fast evolution in this step <cit.>. In single-photon experiments we can simulate the dynamics of a quantum system through unitary operators, thus the expansion and compression steps are implemented by unitary U_e(τ) and U_c(τ), respectively, where τ is the total compression/expansion time interval (adopted to be the same in both steps). By writing the expansion/compression Hamiltonian as H_c/e(t) = ħω_c/e(t) σ_y, the unitary evolution operator is given as U_c/e(ω̅τ) = e^-iω̅_c/eτσ_y, where we denote ω̅_c/e = (1/τ)∫_0^τω_c/e(t)dt. 



§ EXPERIMENTAL IMPLEMENTATION


In this section we discuss the optical encoding of the Otto cycle discussed above so that we can simulate it with our particular system, with a general schematic representation shown in Fig. <ref>. The working substance  and an auxiliary (ancilla) system are encoded in the degrees of freedom of a laser beam. The qubit associated with the machine, in which we will extract/introduce heat and work, are the two independent photon polarization states |V⟩ (vertical) and |H⟩ (horizontal).



 §.§  PD channel with linear optical circuits


The experimental implementation of PD channels that simulate the thermal reservoirs has been performed by using linear optical circuits. In our experiment, instead of a single-photon source, we used an intense laser beam that can be described by a coherent state with a macroscopic photon number. This approach has been successfully explored in literature in different scenarios <cit.>. For this reason, we will present the experiment by using Dirac notation for polarization states once the discussion for single-photon states is straightforward. We encoded the qubit in the polarization degree of freedom and the environment in the propagation direction (path). For the polarization states, we have a two-level system, where we can associate the horizontal polarization as a ground state (|H⟩_S≡|0⟩_S) and the vertical polarization as an excited state (|V⟩_S≡|1⟩_S). In case of the propagation direction, we encoded the path also as a two-level system, with orthogonal directions, k⃗_0 and k⃗_1, representing the reservoir ground (|0⟩_R) and excited state (|1⟩_R), respectively.

The scheme for the PD channel is shown in Fig. <ref> (PD1, red square in the circuit). To describe the channel action on the polarization states, let us consider, without loss of generality, an incoming laser beam described by a right-circular polarized state 

    |ψ_RC⟩ = 1/√(2)( |H⟩ - i |V⟩)  ,

which will interact with the reservoir. It is worth mentioning that the density matrix associated with above state is written as ρ_RC = |ψ_RC⟩⟨ψ_RC|=(1/2)( - σ_y). Thus, by considering  thermal states such that k_BT_c≈ħω_ini/3, we can use the approximation tanh(ħω_iniβ_c) ≈ 1 in Eq. (<ref>) to see that we (approximately) get the 
same the density matrix as that for the state |ψ_RC⟩. Therefore, as the state |ψ_RC⟩ arrives at the channel, the polarization beam splitter PBS1 transmits (reflects) the horizontal (vertical) polarization state. In this way, the H-polarization component (H) goes to the half-wave plate HWP4@θ_h (θ_h = 0^∘), where no change occurs in the polarization component (H) of transmitted arm. On the other hand, for the reflected arm, the V-polarization component passes through HWP3@θ_v, implementing the transformation

    |V⟩_S→sin(2θ_v)|H⟩_S + cos(2θ_v)|V⟩_S .

Moreover, in this reflected arm, we introduced a piezoelectric ceramic (PZT) placed in the mirror for adjusting the difference of phase (Δϕ) between the two arms.
In this way, by adjusting Δϕ = 0, we have the state of polarization of the transmitted arm going out to PBS2 in the path |0⟩_R. 
For the reflected arm, after the transformation implemented by HWP3, the H-polarization component of (<ref>) leaves PBS2 in the path |1⟩_R 
and the V-polarization component of (<ref>) is reflected to the path |0⟩_R.
The last stage of the channel is implemented by another half-wave plate (HWP5@45^∘) introduced in the path |1⟩_R. This device turns  |H⟩_S→|V⟩_S. Thereby, the transformations implemented by this channel in the initial state |ψ_RC⟩  can be written as the map (by using the notation |x⟩|y⟩ = |x⟩_S|y⟩_R)

    (|H⟩ - i|V⟩) |0⟩→|H⟩|0⟩ - i[cos2θ_v|V⟩|0⟩ + sin2θ_v|V⟩|1⟩]  ,

up to a normalization factor 1/√(2) on both sides. If we consider the definition of the PD channel in terms of its Kraus operators <cit.>, we obtain the map

    (|0⟩ - i|1⟩)|0⟩→   |0⟩|0⟩  - i[ 1-p(τ_d) ]^1/2|1⟩|0⟩
       - i p^1/2(τ_d) |1⟩|1⟩ ,

where p(τ_d) = 1 - e^-γτ_d, being γ the decay rate. By comparing Eqs. (<ref>) and (<ref>), we get 

    cos^2 (2θ_v) = 1-p(τ_d)  .
 
Therefore, HWP3 simulates the time evolution during the PD channel. 
For the initial condition p(τ_d=0)=0, where the system does not interact with the reservoir, we have θ_v = 0^∘. 
In this case, HWP3 does not implement any change in the polarization state and, as expected, nothing happens with the initial state. 
Consequently, coherence does not decrease. On the other hand, for the asymptotic behavior, p(τ_d→∞)=1, HWP3 implements 
the maximum rotation in the polarization state and the state completely loses its coherence.



 §.§ Otto cycle with linear optical circuits


In order to realize the Otto cycle, we start with the state preparation. As shown in Fig <ref>, a vertically polarized DPSS laser (1.5 mW power, λ = 532 nm) passes through a quarter-wave plate QWP1@-45^∘ to produce a right-circular polarization that is the analog of the initial state |ψ_RC⟩. The laser beam passes into a spatial filter SF in order to improve the fundamental transverse mode quality. The initial state is verified by performing state tomography in the polarization of the laser beam at point T_A. Polarization tomography can be performed by following Ref. <cit.>. A detailed procedure to perform optical polarization tomography is presented in Appendix <ref>. 


The adiabatic expansion AB is performed by the unitary evolution operator U_e(ω̅τ), where the gap expansion given by ω_e(t) = ω_0(1 - t/τ) + 2ω_0(t/τ) is realized by two half-wave plates HWP1 and HWP2, with their fast axes performing an angle of θ between them. As we show in Appendix <ref>, by using the Jones matrices S(α) for polarization manipulation <cit.>, 
 it follows that the dimensionless quantity 
ω_0τ can be associated with the angle θ between HWP1 and HWP2 through (3/2)ω_0τ = θ. 
In our experiment, we consider θ = 3π/2, so that ωτ = π. Note that, from the initial state, the circular polarization remains unchanged up to a global phase, which corresponds to the evolution of the eigenstates of the Hamiltonian. At point T_B, we perform the tomography of the evolved state in the same way as performed in T_A (see Appendix <ref>).

Following the cycle, the step BC corresponds to the hot reservoir. This part is simulated by PD1. Note that the amount of heat Q_B→C is related to the angle θ_V, as described above to the evolution in the PD channel. A new tomography tracing out the environment is performed in T_C in the two outputs of PD1 exactly as depicted in the measurement box at the end of the circuit.  The couples HWP_13/QWP_2 and HWP_14/QWP_3 are responsible by basis choices of tomography measurements while HWP_15 and the incidence of the two arms in PBS_TR correspond to tracing out the environment as discussed in Ref. <cit.>.

The step CD corresponds to the adiabatic compression and is also realized by two half-wave plates. Note that each output of PD1 passes through a couple of wave plates (HWP_6HWP_7 and HWP_8HWP_9) at the same angles of the couple HWP_1HWP_2 of step AB. This set of angle simulates the compression to the same initial volume, since the gap compression is considered as ω_c(t) = ω_0(t/τ) + ω_0(1 - t/τ). A new tomography of the state is performed at T_D. 

In order to complete the cycle, at step DA, it is necessary to return to the initial state. Therefore, the action of the PD channel should be undone. In order to accomplish such an assignment, we designed an optical circuit, namely, the inverted-PD channel (I-PD) to simulate the last step of the thermal engine (See the block I-PD in Fig. <ref>). Note that HWP_10 at 45^∘  undoes the action of HWP_5. In the interferometer, HWP_11 at θ_V2 = θ_V and HWP_12  at θ_H2=0 undo the action of HWP_3  at + θ_V. PBS_4 regroups the arms |0⟩ and |1⟩, with the relative phase controlled by PZT_2 in order to obtain the initial state (circular polarization).
It is worth mentioning that our apparatus is not an actual thermalization process, but it is able to reproduce the desired output
thermal states at the end of the evolution. To this end, we used the standard procedure to simulate
decoherence in quantum circuits by taking the path degree of freedom as an ancilla, where the phase damping is recorded after the PD channel and recovered by the I-PD.  
These simulated thermal baths are fundamental in our implementation since we do not have a natural source
of thermal reservoir in our system. A tomography is performed at point T_A' that corresponds to the initial state at T_A.

The thermodynamic quantities are evaluated by adopting regularly spaced θ_V, namely 0^∘, 8^∘, 16^∘, 22.5^∘, 29^∘, 37^∘, and 45^∘. The minimum angle corresponds to the absence of interaction between the qubit and the reservoir, namely, 
θ_V=0^∘ means that no heat is exchanged, while θ_V=45^∘ (maximum angle) means maximum heat exchange obtained for complete dephasing of the initial state.For each θ_V, the whole cycle is performed so that hot and cold reservoirs with different temperature ratios are simulated.






§ RESULTS AND DISCUSSION





Let us start with the characterization of the initial state (right circular polarization). After preparation, we perform a tomographic measurement according to the procedure detailed in Appendix A. The images associated to each measured basis are presented in Fig. <ref>. Note that for the R/L basis only one port (I_β) is lightened, as expected. With the intensities I_α and I_β we can  reconstruct the density matrix analog to that of the quantum state |ψ_RC⟩. The density matrix for the initial state ρ_ini^exp is


    ρ_ini^exp =  ( [            0.5134 0.0033 + 0.4999 i; 0.0033 - 0.4999 i            0.4865;                   ] )  .

This result is in very good agreement with the theoretical prediction, which is given by

    ρ_ini =  ( [   0.5  0.5i; -0.5i   0.5;       ] )  .

From Eqs. (<ref>–<ref>) we compute heat Q and work W from the internal energy variation Δ U at each step of the Otto cycle.  
By using the definition of internal energy as U = Hρ, for some reference Hamiltonian H, we evaluate U from the experimental density matrix 
after each thermodynamical process, which is obtained by performing quantum state tomography. Then, we have

    W    =   ρ(τ) H(τ) - ρ(0) H(0) , 
      
    Q    =   H(τ)[ρ(0) - ρ(τ)] .

where W is the work performed during a unitary evolution and Q is obtained through a nonunitary process keeping the Hamiltonian constant (See Appendix <ref>). In our experimental implementation, we start from the thermal state of the Hamiltonian H_e(0) = ħω_0σ_y, with tanh (ħω_0β_c) ≈ 1 (corresponding to ħω_0β_c≈ 3). In this case, the initial state is ρ_ini^th≈|ψ_RC⟩⟨ψ_RC|. The results are presented in Fig <ref>. In Fig <ref>, we show the  internal energy variation Δ U/ħω_0 as a function of ħω_0 β_h for each step of the cycle for seven different values of the inverse hot temperature β_h. In this plot, work and heat have experimentally been obtained by Eqs. (<ref>) and (<ref>), respectively. For a closed cycle, Δ U must be zero, which can be observed by summing up W and Q for all the curves at a 
fixed value of β_h. Fig. <ref> shows the extracted work, quantified from difference |Q_B→C|-|Q_D→A|, due to the coupling of the system with thermal baths at different inverse temperatures β_c and β_h. As expected, the extracted work decreases as the inverse hot temperature 
β_h increases. In addition, note that the energy balance Δ U_↻, is kept close to zero, as theoretically predicted.

To study the entropy production in the thermodynamic cycle, we consider the analysis of the expansion/compression step followed by the 
subsequent thermalization process. The irreversible contribution to the entropy variation is then given by <cit.>

    ⟨Σ⟩_e/c = Δ S_e/c - β_h/cQ_e/c ,

where Δ S_e/c = S^fin_e/c - S^ini_e/c is the (von Neumann) entropy variation during the thermalization process after expansion/compression step, with S_e/c = - ρ^hot/cold_thlnρ^hot/cold_th, and Q_e/c is the amount of heat exchanged during such a  process. It is possible to show that the above equation can be written as the quantum relative entropy (Kullback–Leibler divergence) ⟨Σ⟩_e/c = [ρ_e/c(τ)||ρ^hot/cold_th] =  ρ_e/c(τ) lnρ_e/c(τ) - ρ_e/c(τ) lnρ^hot/cold_th, where ρ_e/c(τ) is the state after expansion/compression step and ρ^hot/cold_th is the 
thermal state at inverse temperature β_h/c to be achieved as a subsequent thermalization process <cit.> (see Appendix <ref> for more details). As originally proposed, ⟨Σ⟩_e/c quantifies the lag between the non-equilibrium state after the expansion/compression unitary step and the 
target equilibrium thermal state. For a recent experimental implementation in NMR, see Ref. <cit.>. For thermodynamical cycles, the entropy production accounts for the dissipated energy during the expansion and compression steps, which may quantify 
quantum friction during the quantum evolution <cit.>. 
The results are shown in Fig. <ref>. Observe that the individual amounts of entropy production ⟨Σ⟩_e and ⟨Σ⟩_c associated with the thermalization after the (unitary) expansion and compression steps, respectively, are nonvanishing, which implies a nonvanishing total entropy production ⟨Σ⟩_↻  for thermal 
baths with distinct inverse temperatures β_c and β_h. Notice that, as theoretically predicted, ⟨Σ⟩_↻ vanishes as β_h gets nearer β_c. The experimental images and the reconstructed density matrix for each step  considering θ_V=22.5^∘ are presented in Appendix <ref>.  Note that the first experimental points in Fig.<ref>), <ref>) and <ref>) correspond to θ_V=45^∘ (β_h / β_c < 10^-4) and the last ones corresponds to θ_V=0^∘ (β_h / β_c = 1).



§ CONCLUSIONS


In summary, we introduced a map from a single-qubit thermal machine into a single-photon setup. The feasibility of this proposal has been experimentally tested by an all-optical 
experiment realized through an intense laser beam. By using the polarization degree of freedom of the laser beam, we encoded a qubit as the working substance, while the two thermal baths are simulated by  correlating the polarization with an auxiliary degree of freedom, which was the propagation path in our experiment. We have then shown how different thermal baths can be implemented with optical devices, with the difference of temperatures controllable through a dimensionless parameter associated with a combination of half-wave plates. 
Agreement between experimental and theoretical results is remarkable, with errors within a 5% range.  It is worth emphasizing that we are proposing  
an all-optical simulation of a thermal machine, which aims at reproducing analogues of Gibbs states, heat transfers and 
work extraction.  However, our investigation opens perspectives for 
optical implementations of other protocols in quantum thermodynamics. In this scenario, the experimental discussion of the performance of quantum refrigerators <cit.> with optical devices and optical quantum thermometers <cit.> is left as future research.



§ ACKNOWLEDGMENTS

The authors acknowledge financial support from the Brazilian funding agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), 
Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), 
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) (Finance Code 001), 
and the Brazilian National Institute for Science and Technology of Quantum Information (INCT-IQ).





§ HEAT AND WORK IN THERMODYNAMIC PROCESSES DRIVEN BY TIME-LOCAL MASTER EQUATIONS
 

Let us consider a quantum system evolving driven by a time-dependent Hamiltonian H(t) and coupled to a reservoir governed by a time-local master equation

    ρ̇(t) = 1/iħ [H(t),ρ(t)] + [ρ(t)]  ,

where [∙] describes the coupling between the environment and the system, so that we can recover a free-decohering evolution as [ρ(t)]=0 for all t∈[0,τ_ev], with τ_ev being the total evolution time. The internal energy of the system can be computed from U(t) = ρ(t)H(t), 
yielding

    dU(t)    = d/dt[ρ(t)H(t)]dt 
       = ρ̇(t)H(t)dt + ρ(t)Ḣ(t)dt  .

Now, let us consider two distinct situations. First, consider that the system evolves under the action of a time-independent Hamiltonian. Thus, we get

    dU(t)|_H(t)=H = ρ̇(t)Hdt  .

We can then use the Eq. (<ref>) and write

    dU(t)|_H(t)=H   = [1/iħ [H,ρ(t)] + [ρ(t)]]Hdt 
       = 1/iħ[H,ρ(t)]Hdt + [ρ(t)]Hdt  .

By using the cyclic property of the trace, we obtain [H,ρ(t)]H = 0, so that

    dU(t)|_H(t)=H   = [ρ(t)]Hdt .

This means that any internal energy variation will be due to fluctuations in the energy level populations arising from the coupling of the system with its environment. 
For this reason, we identify this change in the internal energy as heat and we write

    dU(t)|_H(t)=H   = dQ = [ρ(t)]Hdt .

Thus, the total exchanged heat is obtained by integration of the above equation as

    Q = ∫_τ_0^τdQ = ∫_τ_0^τ dU(t)|_H(t)=H .

By using that ∫_t_1^t_2df(t) = f(t_2) - f(t_1) for any analytical function f in interval = [t_1,t_2], we get

    Q = ∫_τ_0^τ dU(t)|_H(t)=H = H[ρ(τ_0) - ρ(τ)] .

On the other hand, consider the case where we have a time-dependent Hamiltonian, but the system is decoupled from its reservoir [[ρ(t)]=0]. 
Thus, from Eq. (<ref>) we get

    dU(t)|_[ρ(t)]=0   = ρ̇(t)H(t)dt + ρ(t)Ḣ(t)dt  ,

where we can use the Eq. (<ref>) to write

    ρ̇(t)H(t)   = 1/iħ [H(t),ρ(t)] H(t)
       = 1/iħ[  H(t)ρ(t)H(t) - ρ(t)H(t)H(t)] 
       = 1/iħ[ ρ(t)H(t)H(t) - ρ(t)H(t)H(t)] = 0 .

Therefore

    dU(t)|_[ρ(t)]=0   = ρ(t)Ḣ(t)dt  .

This means that, when we have a unitary evolution, changes in the internal energy may be achieved by variations in the instantaneous energy eigenspectrum of the time-dependent Hamiltonian. So, by taking the Hamiltonian as controlled by external fields governing the gap expansion/compression, we identify dU(t)|_[ρ(t)]=0 as an amount 
of work performed by/on these fields and we write

    dU(t)|_[ρ(t)]=0   = dW(t) = ρ(t)Ḣ(t)dt  .

Hence, the work performed on/by the system can be obtained by integration of the above equation as

    W = ∫_τ_0^τdW(t) = ∫_τ_0^τdU(t)|_[ρ(t)]=0 ,

and therefore we conclude that

    W = Δ U(τ) = ρ(τ) H(τ) - ρ(τ_0) H(τ_0) .




§ OPTICAL POLARIZATION TOMOGRAPHY

	
The state of one qubit can be described in terms of the density operator 

    ρ = 1/2( + ∑_i=1^3 r_i σ_i)  ,

where the matrices (σ_1,σ_2,σ_3) are the Pauli matrices and (r_1,r_2,r_3) are the components of the Bloch vector. 
Considering optical polarization, it is possible  to rewrite the density operator given by Eq. (<ref>) in terms of the Stokes parameters 
S_i (i=0,1,2,3) <cit.> so that we have

    ρ = 1/2(S_0  + ∑_i=1^3 S_i σ_i)  ,

where S_0 = ρ=1 and each S_i (i=0,1,2,3) is defined by

    S_i = σ_iρ .


As provided in Ref. <cit.>, in optical polarization tomography, 
each component S_i can be obtained by defining a set of projective measurements in the polarization states of light as 

    S_0    =    P_|H ⟩ + P_|V ⟩ , 
     
    S_1    =    P_|D ⟩ - P_|AD ⟩ , 
     
    S_2    =    P_|L ⟩ - P_|R ⟩ , 
     
    S_3    =    P_|H ⟩ - P_|V ⟩ ,

where |H ⟩ (|V ⟩) represents the horizontal (vertical) polarization, |D ⟩ (|AD ⟩) represents the linear diagonal (linear anti-diagonal) polarization,  
|L ⟩ (|R ⟩) represents the left-circular (right-circular) polarization of light, and P_|α⟩ is the probability of obtain the component 
|α⟩, with α= {H,V}, {D, AD}, {R, L}.

Fig. <ref>  shows how we can implement the polarization tomography. Experimentally, to obtain the Stokes parameters S_0 and S_3, we make a projective measurement in the basis {H, V} using a Polarization Beam Splitter (PBS). To obtain the Stokes parameter S_1 the projective measure is performed by using the HWP@22.5^∘ in association with the PBS, that corresponds to measurement in the basis {D, AD}. Finally, to obtain the Stokes parameters S_2 we use the QWP@0^∘, the HWP@22.5^∘, and the PBS in order to measurement in the basis {R,L}. Lastly, a charged-coupled device (CCD) records one single image with the intensity of each component projected on a screen.

Considering our apparatus, we measure simultaneously the intensity of each projected components I_α and I_β, 
being β the complementary (orthogonal) basis component to α, e.g., for α=H, then β=V. 
We can associate the probability P_|α⟩ with the normalized intensity as

    P_|α⟩ = I_|α⟩= I_α/I_α + I_β.

Then, the Stokes parameters S_i can be obtained from the intensities of tomographic measurements and we can reconstruct the density matrix for any polarization state by using Eq. (<ref>).



§ SIMULATING EXPANSION/COMPRESSION STEP BY JONES MATRIX
 

The polarization rotation Jones matrix is defined by

    S(α) = [   cos (α) - sin (α);   sin (α)   cos (α) ] .

 Eq. (<ref>) may be used to encode changes in character of light as it passes through a partial polarizer, for example, with α denoting the rotation angle between the initial and final polarizations. The matrix S(α) is in the SO(2) group. Then, it can be used to provide a transformation of a vector |ψ⟩ to a transformed vector |ψ^'⟩ such that the norm of |ψ⟩ is preserved. We can write S(α) = e^-i J α, with the generator J given by <cit.>

    J = [   0 - i;   i   0 ] ,




On the other hand, the evolution operator of the working substance driven by Hamiltonian H = ħω(t)σ_y is given by

    U(t,t_0) = exp[ -i/ħ∫_t_0^t H(ξ) d ξ]  ,

which can be rewritten as

    U(t,t_0) = exp[ - i ω̅ (t-t_0) σ_y]  ,

where ω̅ = (t-t_0)^-1∫_t_0^tω(ξ) d ξ. In this case, σ_y can be viewed as the dynamical generator associated with 
the evolution operator U(t,t_0), 
being exactly the same as the generator of J. Thus, we conclude that there is a correspondence between evolution operator U(t,t_0) and the Jones matrix S(α) by taking

    α = ω̅ (t-t_0) = ∫_t_0^tω(ξ) d ξ .


Therefore, by setting the expansion/compression duration as t-t_0 = τ and the frequency ω, it is possible to adjust α to simulate the desired dynamics. 
In particular, in the expansion step considered in main text we have ω_e(t) = ω_0(1 - t/τ) + nω_0(t/τ) (with t_0=0), where n>1 in 
order to obtain ω_e(0)<ω_e(τ). Then, we get the expansion Jones parameter α_e as

    α_e = ∫_0^τω_e(ξ) d ξ =  (n+1) ω_0τ/2 .

On the other hand, in the compression step, we need to get ω_c(0)>ω_c(τ), so that we write ω_c(t) = nω_0(1 - t/τ) + ω_0(t/τ). Therefore the compression Jones parameter α_c reads

    α_c = ∫_0^τω_c(ξ) d ξ =  (n+1) ω_0τ/2 .

This means that we can simulate the expansion and compression steps with the same Jones parameter. In an experimental approach, the Jones matrix in Eq. (<ref>) can be implemented by using an arrangement of two wave plates HWP_θ, where the parameter θ is associated with the angle between the fast axis and the vertical direction. The matrix representation of the HWP is

    HWP_θ = [  cosθ  sinθ;  sinθ -cosθ ] ,

so that we can obtain the following result

    S(α) = HWP_2α·HWP_α = [   cos (α) - sin (α);   sin (α)   cos (α) ] .

In conclusion, the HWP angle needs to be adjusted as α = 3 τω_0/2, so that the first HWP is 3 τω_0/2 and the second one is 3 τω_0, 
implying a relative angle between the HWPs given by θ_r = 3 τω_0/2.




§ EXPERIMENTAL IMAGES AND DENSITY MATRIX FOR A COMPLETE CYCLE


In this Appendix, we will present the resulting density matrix for each step of the engine for θ_V=22.5^∘, in order to show a complete measurement set of the experiment. Fig. <ref> presents the images obtained from tomographic measurement for each step of the cycle. 


 







 





Regarding the Otto's cycle in Fig. <ref>, the images of step AB are presented in Fig. <ref>-(a). 
The theoretical values of Stokes parameters for this case are S_1=0, S_2=1, and S_3=0. Note that only  R/L basis measurement should have a value different from zero (I_α≠ I_β), which can be inferred from Fig. <ref>-(a). Then, the experimental density matrix for the step AB is given by

    ρ_A→ B^exp =  ( [           0.49663 0.0045 + 0.4966 i; 0.0045 - 0.4966 i            0.5033;                   ] )  .

The resulting images for the end of the step BC are shown in Fig. <ref>-(b). Here, the expected theoretical values for the Stokes parameters are 
S_1 = 0, S_2 =-0.7071, and S_3 = 0, that means I_α≠ 0, as observed in Fig. <ref>-(b). The resulting density matrix is given by

    ρ_B→ C^exp =  ( [            0.5057 0.0174 + 0.3407 i; 0.0174 - 0.3407 i            0.4942;                   ] )  .

Note the change of the state due to the PD channel at this stage. 
Following the cycle, the images for the end of the step  CD are shown in Fig. <ref>-(c). The obtained density matrix associated with the stage is

    ρ_C→ D^exp =  ( [            0.5190 0.0041 + 0.3482 i; 0.0041 - 0.3482 i            0.4809;                   ] )  .

As can be noted ρ_C→ D^exp is very similar to ρ_B→ C^exp, as expected, since the performed unitary operation does not change the polarization. 
The last step of the Otto's cycle (DA) gives us the images presented in Fig. <ref>-(d). The respective density matrix is given by

    ρ_D→ A^exp =  ( [            0.5118 0.0115 + 0.4503 i; 0.0115 - 0.4503 i            0.4881;                   ] )  .

Note that, if we compare the experimental density matrices ρ_ini^exp and ρ_D→ A^exp we can see that they are approximately equal, with the small differences due to the losses after the laser beam passes through all the optical elements of the experiment.



§ ENTROPY PRODUCTION AND QUANTUM RELATIVE ENTROPY
 

In order to study the relation between entropy production and quantum relative entropy, we will follow a similar procedure as preformed in 
Refs. <cit.>. This approach has also been previously discussed in Ref. <cit.>.

Consider the system initially prepared in a thermal state ρ^ini_th at temperature β and internal Hamiltonian H_ini. The system can then be driven from the thermal state ρ^ini_th to another thermal state ρ^fin_th with final Hamiltonian H_fin. This process can be done through a sequence of equilibrium states or through a far-from equilibrium process. In a non-equilibrium evolution, we decouple the system from the thermal bath and implement a unitary evolution with a time-dependent driving Hamiltonian H(t), where the boundary conditions are H(0) = H_ini and H(τ) = H_fin. The final state achieved at time τ is  
denoted by ρ(τ). Then, from the non-equilibrium state ρ(τ), we drive the system to the thermal state associated with a bath at inverse temperature given by β_fin. 
During this process, an amount of irreversible entropy arises, reading

    Σ  = β( w - Δ F )  ,

where Δ F is the Helmholtz free energy variation and w is the work realized on/by the system during the unitary evolution. In Ref. <cit.>, it is shown 
that Eq. (<ref>) can be expressed in terms of the relative entropy between the states ρ(τ) and ρ^fin_th. Namely, we have

    Σ = ρ(τ) lnρ(τ) - ρ(τ)lnρ^fin_th = [ρ(τ)||ρ^fin_th]  .

On the other hand, by computing the total internal energy variation of the system we can write Δ U = w + Q^th, where Q^th is the amount of 
heat exchanged between system and reservoir during the non-unitary process towards the thermal state at β_fin. Now, we can use the relation Δ U = Δ F + (1/β)Δ S, where Δ S is the von Neumann entropy variation during the thermalization, to find

    β( w - Δ F ) = Δ S - βQ^th ,

which is the standard expression for the (irreversible) entropy production adopted in Eq. (<ref>) of the main text. 
Now, by using Eqs. (<ref>) and  (<ref>) in Eq. (<ref>), we obtain 

    Δ S - βQ^th = [ρ(τ)||ρ^fin_th]  .





10
	
	Geva:92
	E. Geva and R. Kosloff,
	J. Chem. Phys. 96, 3054 (1992).
	
	Lutz:16
	O. Abah and E. Lutz,
	EPL (Europhysics Letters) 113, 60002 (2016).
	
	Noah:10
	N. Linden, S. Popescu, and P. Skrzypczyk,
	Phys. Rev. Lett. 105, 130401 (2010).
	
	Alicki:13
	R. Alicki and M. Fannes,
	Phys. Rev. E 87, 042123 (2013).
	
	Binder:15
	F. C. Binder, S. Vinjanampathy, K. Modi, and J. Goold,
	New J. Phys. 17, 075015 (2015).
	
	Kieu:04
	T. D. Kieu,
	Phys. Rev. Lett. 93, 140403 (2004).
	
	Scully:03
	M. O. Scully, M. S. Zubairy, G. S. Agarwal, and H. Walther,
	Science 299, 862 (2003).
	
	Kosloff:84
	R. Kosloff,
	The Journal of chemical physics 80, 1625 (1984).
	
	Rezek:06
	Y. Rezek and R. Kosloff,
	New J. Phys. 8, 83 (2006).
	
	Henrich:07
	M. J. Henrich, G. Mahler, and M. Michel,
	Phys. Rev. E 75, 051118 (2007).
	
	Geusic:59
	J. Geusic, E. Schulz-Du Bois, R. De Grasse, and H. Scovil,
	Journal of Applied Physics 30, 1113 (1959).
	
	Quan:05
	H. T. Quan, P. Zhang, and C. P. Sun,
	Phys. Rev. E 72, 056110 (2005).
	
	Lutz:16-Science
	J. Roßnagel et al.,
	Science 352, 325 (2016).
	
	Scovil:59
	H. Scovil and E. Schulz-DuBois,
	Phys. Rev. Lett. 2, 262 (1959).
	
	Peterson:18
	J. P. Peterson et al.,
	arXiv preprint arXiv:1803.06021  (2018).
	
	Klatzow:17
	J. Klatzow et al.,
	Phys. Rev. Lett. 122, 110601 (2019).
	
	Nielsen:Book
	M. A. Nielsen and I. L. Chuang,
	Quantum Computation and Quantum Information: 10th Anniversary
		Edition, 10th ed. (Cambridge University Press, New York, NY, USA, 2011).
	
	PRA.78.Davidovich
	A. Salles et al.,
	Phys. Rev. A 78, 022322 (2008).
	
	PRL.99.Topo
	C. E. R. Souza, J. A. O. Huguenin, P. Milman, and A. Z. Khoury,
	Phys. Rev. Lett. 99, 160401 (2007).
	
	PRA.82.Borges
	C. V. S. Borges, M. Hor-Meyll, J. A. O. Huguenin, and A. Z. Khoury,
	Phys. Rev. A 82, 033833 (2010).
	
	NaturePhot.Kagalwala-2012
	K. H. Kagalwala, G. Di Giuseppe, A. F. Abouraddy, and B. E. A. Saleh,
	Nature Photonics 7 (2012).
	
	Opt.Lett.Balthazar-2016
	W. F. Balthazar et al.,
	Optics Letters 41, 5797 (2016).
	
	PRA.77.Cadu-Cripto
	C. E. R. Souza et al.,
	Phys. Rev. A 77, 032345 (2008).
	
	PRA.83.Zela-Teleport
	A. Z. Khoury and P. Milman,
	Phys. Rev. A 83, 060301 (2011).
	
	Opt.Exp.Souza-Cnot-2010
	C. E. R. Souza and A. Z. Khoury,
	Optics Express 18, 9207 (2010).
	
	JOPS-B.Cod.Op.Balthazar-2016
	W. F. Balthazar and J. A. O. Huguenin,
	Journal of the Optical Society of America B 33, 1649 (2016).
	
	PRA.97.Enviro-Passos
	M. H. M. Passos et al.,
	Phys. Rev. A 97, 022321 (2018).
	
	Talkner:07
	P. Talkner, E. Lutz, and P. Hänggi,
	Phys. Rev. E 75, 050102 (2007).
	
	Alicki:79
	R. Alicki,
	J. Phys. A: Math. Gen. 12, L103 (1979).
	
	Anders:13
	J. Anders and V. Giovannetti,
	New J. Phys. 15, 033022 (2013).
	
	Lutz:18
	O. Abah and E. Lutz,
	Phys. Rev. E 98, 032121 (2018).
	
	Photonic-State-Tomography
	J. B. Altepeter, E. R. Jeffrey, and P. G. Kwiat,
	Adv. Atom. Mol. Opt. Phy. 52, 105 (2005).
	
	Jones-i:41
	R. C. Jones,
	J. Opt. Soc. Am. 31, 488 (1941).
	
	Procaccia:76
	I. Procaccia and R. Levine,
	J. Chem. Phys. 65, 3357 (1976).
	
	Jarzynski:97
	C. Jarzynski,
	Phys. Rev. Lett. 78, 2690 (1997).
	
	Parrondo:09
	J. M. Parrondo, C. Van den Broeck, and R. Kawai,
	New J. Phys. 11, 073008 (2009).
	
	Lindblad:Book
	C. Lindblad,
	Non-equilibrium entropy and irreversibility (Reidel Publishing
	Company, Dordrecht, Holland, 1983).
	
	Batalhao:15
	T. B. Batalhão et al.,
	Phys. Rev. Lett. 115, 190601 (2015).
	
	Adolfo:14
	A. del Campo, J. Goold, and M. Paternostro,
	Sci. Rep. 4, 6208 (2014).
	
	Plastina:14
	F. Plastina et al.,
	Phys. Rev. Lett. 113, 260601 (2014).
	
	Brask:15
	J. B. Brask and N. Brunner,
	Phys. Rev. E 92, 062101 (2015).
	
	Haupt:14
	F. Haupt, A. Imamoglu, and M. Kroner,
	Phys. Rev. Applied 2, 024001 (2014).
	
	Hofer:17
	P. P. Hofer, J. B. Brask, M. Perarnau-Llobet, and N. Brunner,
	Phys. Rev. Lett. 119, 090603 (2017).
	
	Tung:Book
	W.-K. Tung,
	Group Theory in Physics (World Scientific, 1985).
	
	Deffner:10
	S. Deffner and E. Lutz,
	Phys. Rev. Lett. 105, 170402 (2010).
	
	Apollaro:15
	T. J. Apollaro, G. Francica, M. Paternostro, and M. Campisi,
	Phys. Scr. 2015, 014023 (2015).
	


