Code for the numerical calculations in the publication 
"Quantum Control of Excitons for Reversible Heat Transfer", C. N. Murphy and P. R. Eastham, Communications Physics  doi:10.1038/s42005-019-0215-8 (2019).

Authors: C. N. Murphy and P. R. Eastham.

This notebook defines a function heat,which calculates the exciton-phonon heat transfer and related quantitites for a quantum-dot exciton, modelled as a two-level system,driven by a chirped Gaussian laser pulse. It uses units such that ℏ=1 throughout. Times are in picoseconds, and temperatures in Kelvin.

The input parameters for the function heat[] are  : 
	({{a, Spectral chirp of driving pulse ((ps^2))}, {theta, Pulse area of the driving pulse (π, before chirping, denoted  (Subscript[Θ, 0]) in the paper)}, {tau0, Pulse duration before chirping (ps, denoted (Subscript[τ, 0])in the paper)}, {T, Temperature of the phonon bath (K)}, {A,  Phonon coupling parameter as defined in paper (ps (K^(-1)), defaults to 0.0112 ps/K as used in paper)}, {omegac, Phonon cut-off frequency ((ps^(-1)), defaults to 2meV as used in paper)}, {Tc, Cold-temperature reservoir for the heat-engine calculations (K, defaults to 2.7K as used in paper)}, {δ0, Detuning of the pulse centre-frequency below the exciton ((ps^(-1)), defaults to 0)}})

The output is a list, with entries as follows :
({{Total heat absorbed from the phonon reservoir after one pulse((ps^(-1)))}, {Total heat emitted into the cold reservoir in one stroke of the heat engine ((ps^(-1)))}, {Work out in one stroke of the heat engine ((ps^(-1)))}, {Efficiency of the heat engine}, {Dot entropy after the driving pulse ((ps^(-1)) (K^(-1)))}, {Exciton-phonon heat current as a function of time ((ps^(-2)))}, {Effective dot temperature as a function of time (K)}, {Dressed-state splitting as a function of time ((ps^(-1)))}, {Dot entropy as a function of time ((ps^(-1)) (K^(-1)))}, {Rabi frequency as a function of time ((ps^(-1)))}, {Occupation of the upper dressed-state  as a function of time}})

sz={{0.5,0},{0,-0.5}};
sx={{0,0.5},{0.5,0}};
sy={{0,-0.5*I},{0.5*I,0}};
kB=0.131;
heat[a_,theta_,tau0_,T_:20,A_:0.0112,omegac_:2.*1.52,Tc_:2.7,δ0_:0.0]:=Module[{dmat,delta,omega,lambda,tau,rplus,rminus,B1,B2,n,specden,gammaabs,gammaem,absterm,emsterm,hdot,commhdot,dsols,rho11,rho12,rho21,rho22},
dmat[t_]:={{rho11[t], rho12[t]},{rho21[t], rho22[t]}};
delta[t_] = δ0-1.(a*t)/(a^2.+tau0^4.); (* delta is ex-omega(t) *)omega[t_]=((theta*Pi)/(2.*Pi(a^2.+tau0^4.)^0.5)^0.5)*Exp[-((t^2.*tau0^2)/(2.(a^2+tau0^4.)))];lambda[t_] = Sqrt[omega[t]^2.+delta[t]^2.];
tau=(((a^2.)+(tau0^4.))/(tau0^2.))^(0.5);
rplus[t_] = (delta[t]/lambda[t])sx +(omega[t]/lambda[t])sz+I*sy;
rminus[t_]=(delta[t]/lambda[t])sx +(omega[t]/lambda[t])sz-I*sy;
B1[t_]=rminus[t].rplus[t].dmat[t]+dmat[t].rminus[t].rplus[t]-2*rplus[t].dmat[t].rminus[t];
B2[t_]=rplus[t].rminus[t].dmat[t]+dmat[t].rplus[t].rminus[t]-2*rminus[t].dmat[t].rplus[t];
n[t_]= (Exp[t/(kB*T)]-1.)^(-1.);
specden[t_]=(A/(Pi*kB))(t^3.)(Exp[-((t^2.)/(omegac)^2.)]);gammaabs[t_]=2*((omega[t]/(2.*lambda[t]))^2.)*Pi*specden[lambda[t]]*n[lambda[t]];gammaem[t_]=2*((omega[t]/(2.*lambda[t]))^2)*Pi*specden[lambda[t]]*(n[lambda[t]]+1.);absterm[t_]=-(gammaabs[t]/2.)B1[t];emsterm[t_]=-(gammaem[t]/2.)B2[t];
hdot[t_]=delta[t]*sz-omega[t]*sx;
commhdot[t_]=dmat[t].hdot[t]-hdot[t].dmat[t];
dsols=NDSolve[{dmat'[t]== absterm[t] +emsterm[t]+I*commhdot[t],dmat[-4.*tau]=={{0,0},{0,1.}}},{rho11,rho12,rho21,rho22},{t,-4.*tau,4.*tau}];
pup=Tr[rplus[t].rminus[t].dmat[t] /. dsols[[1]]];
pdown=Tr[rminus[t].rplus[t].dmat[t] /. dsols[[1]]];
qoft=lambda[t]*(gammaem[t]*pup-gammaabs[t]*pdown);(* Conor's code uses omega[t]! *)
soft=With[{dmoft=dmat[t] /. dsols[[1]]},kB*Apply[Plus,Map[-(#*Log[#])&,Eigenvalues[dmoft]]]];
q=NIntegrate[qoft,{t,-4*tau,4*tau},Method->{Automatic,"SymbolicProcessing"->0}];
dmfinal=dmat[t] /. dsols[[1]] /. t->4*tau;
sfinal=kB*Apply[Plus,Map[-(#*Log[#])&,Eigenvalues[dmfinal]]];
sirr=sfinal+q/T;
(* Change conventions so qh, qc, w positive for normal operation of heat engine *)
qh=-q;
qc=Tc*sfinal;
wo=Re[qh-qc];
(* Returns qh, qc, work, efficiency, final entropy, heat as a function of time, effective dot temperature as function of time *)
{qh,qc,wo,Max[0,wo]/qh,sfinal,-qoft,-lambda[t]/(kB*Log[pup/pdown]),lambda[t],soft,omega[t],pup}]

Example usage : calculate the efficiency of the heat engine using the parameters of the ' Carnot pulse' and plot the  phonon-exciton heat current as a function of time during this pulse.

carnotpulseresults=heat[10,9,0.5,20,0.0112,2*1.52,2.7,2.5];
Print["Efficiency as a fraction of Carnot efficiency=",Chop[carnotpulseresults[[4]]]/(1-2.7/20)]
Print[]
psm2topw=1.054*^-34*1*^36;  (* Conversion factor from 1/(ps^2) to pW *)
Plot[carnotpulseresults[[6]]*psm2topw,{t,-40,40},AxesLabel->{"t (ps)","Heat current (pW)"}]