Our predicted lifetimes are 10% and 30% longer than the time-domain measurements of Elks and Western, 1999 and Yamasaki et al. [Yamasaki et al., 2003, Yamasaki et al., 2005], respectively (ignoring the outlying $B(3)$ lifetime of Elks and Western, 1999).

§.§ Predissociation broadening

Deperturbed predissociation linewidths of the observed $B(v)$ and $C(v)$ states.
The rotation-dependent widths of $B(8)$ are shown extrapolated to $J=0$.  Open symbols: Fitted parameters. Closed symbols: Fixed to assumed values.
Upper figures: Widths of $\ensuremath{\Omega}=0$ sublevels for all isotopologues. Lower figures: Widths of all $\ensuremath{\Omega}$ levels for ^32S ^16O only.

The deperturbed line broadening fitted to observed v and v bands is plotted in Fig. <ref> and shows significant vibrational, isotopologue and $\ensuremath{\Omega}$ dependence.
The deperturbed widths for $\ensuremath{\Omega}=0$ and 1 $e$-parity levels become rapidly mixed with increasing rotation so that the observed widths of high-$J$ levels converge.

The upper limits of ^32S ^16O $B(v)$ widths estimated by Liu et al., 2006 are in agreement with the values determined here, apart from their limits $\ensuremath{\Gamma}_{B(6)}<1.2$ and $\ensuremath{\Gamma}_{B(11)}<\np[cm^{-1}]{1}$ that fall below our measured widths of $(3.0\pm{}0.3)$ and $3.1-\wnfwhm{3.9}$, respectively.

The modelled predissociation width pattern is likely caused by spin-orbit interaction between and various crossing unbound states known from quantum-chemical calculations.
Yu and Bian, 2011 compute all singlet, triplet, and quintet states dissociating to ground state  and  atoms which may provide such predissociation channels.
They find significant interaction energies mixing with their $(1)\,{}^5\ensuremath{\Pi}$ and $(2)\,{}^5\ensuremath{\Pi}$ states, and the outer limb of an adiabatic $C\,{}^3\ensuremath{\Pi}$ state.
The bound $(1)\,{}^5\ensuremath{\Pi}$ state crosses below the $\ensuremath{\ce{S}({}^3{\rm P})}+\ensuremath{\ce{O}({}^3{\rm P})}$ dissociation limit and provides a threshold dissociation channel beginning with $B(v=4)$.
This crossing likely also explains the predissociation of rotationally-excited $B(v=0-3)$  with their $J$-thresholds known experimentally [Martin, 1932, Clerbaux and Colin, 1994].
The $(2)\,{}^5\ensuremath{\Pi}$ state is predicted to cross near $v=13$ and could provide an explanation for the rapidly increasing predissociation widths of higher levels.
An outer-limb crossing with the adiabatic $C\,{}^3\ensuremath{\Pi}$ state may further enhance predissociation near $v=10$.
The rapid variation of width with $v$ between threshold and these crossings likely results from the varying overlap of bound- and unbound-state radial wavefunctions, as is typical of predissociation by outer-limb crossing repulsive states [Lefebvre-Brion and Field, 2004].
Yu and Bian, 2011 predict additional crossings and spin-orbit interactions between near $v=7$ and 11, and the repulsive outer limbs of two  states.
This purely $\ensuremath{\Omega}=1$ interaction could explain the enhanced $\Omega=1$ dissociation widths we find for $B(9)$ and $B(10)$.
Finally, we note that the $J$- and $\Omega$-dependencies of the $B(v=8)$ linewidths for ^32S ^16O, displayed in Fig. <ref>, are typically characteristic of predissociation by a ${}^3\Pi$ state [Julienne and Krauss, 1975], also observed in some levels of the $B$ state of O2 [Lewis et al., 1994].
Thus, the state is likely to be partially involved in the, possibly complex, $B(v=8)$ predissociation mechanism.

Our analysis of $C(v)$ levels reveals a broad width minimum around $v=4$ and up to factor-of-five differences between $\ensuremath{\Omega}$-substates of the same vibrational level, with either $\ensuremath{\Omega}=0$ or 1 being the most broadened.
The large widths of $C(v=6,\ensuremath{\Omega}=0)$ and $C(v=7)$ are consistent with being nearest to the  avoided crossing shown in Fig. <ref>.

§.§ Rotational line lists and cross sections

Simulations of absorption into the mixed $B(7)\sim C(2)$ levels computed from constants fitted band-by-band and a global electronic-state potential-energy-curve model.

Synthetic photodissociation spectra computed at 360 K.

The band-by-band modelling detailed in Sec. <ref> results in a rotational line list for all observed bands, while the global model described in Secs. <ref> and <ref> also permits the calculation of line frequencies and strengths for bands that are not experimentally observed in some isotopologues.
We combine band-by-band and global-model derived line lists in order to generate a complete line list that is as accurate as possible.

A global-model line list of rotational transitions terminating on $B(v=0-30)$ and $C(v=0-7)$ levels up to $J=50
$ was computed using the potential-energy curves, electronic state interactions, and transition moments presented in Secs. <ref> and <ref> for all measured isotopologues and, for good measure, ^18O-substituted species.
Band-by-band-fitted deperturbed linewidths are added to the global-model line list.
The unmeasured linewidths of S-substituted isotopologues for $B(v>16)$ levels, and all levels of O-substituted isotopologues are set to their ^32S ^16O values.

A cross section for absorption into the strongly-coupled $B(7)\sim C(2)$ levels is computed from the band-by-band constants fitted to the experimental spectra, as well as from the global-model line list.
These are compared at 360 K in Fig. <ref>, and the detailed rotational structure and mass-dependent spin-orbit mixing of levels is well reproduced by the global model.
The cross sections of other SO bands show similar agreement when making this comparison, with the principal differences being due to band-strength differences and small frequency shifts.
The globally-modelled strengths smooth over the background intensity and contaminant absorption uncertainties affecting band-by-band analysis, and are therefore more accurate.
For example, the band-by-band ^36S ^16O cross section in Fig. <ref> is clearly an overestimate when compared with the isotopically self-consistent global model.
Conversely, the band-by-band modelled line frequencies of observed bands are more accurate than the global model.
We construct a final hybrid line list by substituting band-by-band measured line frequencies into the electronic-state model.
Photoabsorption cross sections computed from this hybrid list for all S-isotopes are shown in Fig. <ref>.

Absorption into $B(v>3)$ and all $C(v)$ levels is completely dissociative given the short lifetimes implied by the large measured transition widths.
Additionally, levels of ^32S ^16O $B(v=3)$ with quantum-number $N\geq 10 $ are known to dissociate, with only lower-$N$ levels contributing to the emission spectra of ^32S ^16O [Martin, 1932, Clerbaux and Colin, 1994].
Emission spectra of S-substituted isotopologues have not been measured and will likely have higher-$N$ predissociation thresholds, as found for ^12S ^18O by Clerbaux and Colin, 1994.
We neglect this difference given the rather small contribution of $B(3)$ to the overall photodissociation cross section and adopt the $N=10$ threshold for all species.
Then we compute a photodissociation cross section from the hybrid line list by including only dissociative upper levels.

Band-by-band and global model rotational line lists, and the recommended hybrid of the two, are permanently available in an online data archive [Heays et al., 2022].

§.§ Comparison with other published cross sections

The present ^32S ^16O cross section computed assuming a rotational temperature of 400 K and degraded by convolution with a unit Gaussian of width 50, compared with the low-resolution photoabsorption cross section of Phillips, 1981.

Phillips, 1981 measured the v cross section for ^32S ^16O at low spectral resolution and a digitised version is plotted in Fig. <ref>.
We compare this with a ^32S ^16O cross section computed from our final hybrid line list assuming a rotational temperature of 400 K and broadened by convolution with a Gaussian function of width 50, in order to approximately match the Phillips, 1981 instrumentally-broadened band profiles.

The relative band strengths are similar and integrated cross sections over the region plotted agree within 15%, with the new cross section being larger.

There are absorption features appearing in the cross section of Phillips, 1981 which are not found in our analysis of SO spectra, for example at 48400 and 52400.
These are broadly aligned with SO2 absorption bandheads and likely result from an imperfect subtraction of contaminant SO2.
Danielache et al., 2014 make a similar comparison of their cross section computed with Phillips, 1981, finding reasonable relative agreement but with some differences in relative v band strengths and frequencies, but, more importantly, the cross section was found to have a three-times smaller magnitude overall.

§ SUMMARY

Detailed spectroscopic data on the $\Bstate(v=4-30)$ and $\Cstate(v=0-7)$ states of SO and  its S-substituted isotopologues are determined from high-resolution absorption spectra covering the 43000 to 51000 (190 to 233 nm) spectral region.
This is the first observation of levels above $v=2$ or a or level in any S-substituted isotopologue.

Most of the observed bands are severely blended due to a high line density and predissociation broadening, and their profiles were fitted to minimally-specified effective Hamiltonians including spin-orbit interactions between neighbouring $B(v_B)$ and $C(v_C)$ levels.
The fitted linewidths are quite variable with respect to vibrational level and quantum number $\ensuremath{\Omega}$ and likely arise from further interactions with unbound electronic states.
An empirical model of and electronic states was constructed to ensure globally-reliable band strengths and to enable extrapolation to unmeasured vibrational levels and isotopologues.

The fitted cross section has a band-by-band relative uncertainty within 5% and an additional 10% absolute uncertainty based on published lifetime and radiative data, and their uncertainties, for the $\Astate(v=1)$ level [Lo et al., 1988, Elks and Western, 1999].
The estimated total uncertainty encompasses agreement with a previously measured ^32S ^16O cross section [Phillips, 1981] and $\Bstate(v=0-2)$ lifetimes [Elks and Western, 1999], but somewhat disagrees with a further set of $\Bstate(v=0-3)$ lifetimes [Yamasaki et al., 2005].
However, the new $\Bstate-\Xstate$ transition moment is in significant disagreement with previous calculations [Sarka and Nanbu, 2019, da Silva and Ballester, 2020] that suggest a 50% smaller cross section, while $\Cstate-\Xstate$ transition moments are found to agree well.

The final result is a comprehensive and spectroscopically-accurate line list of photodissociating far-ultraviolet rovibronic transitions for all ^xS^16O isotopologues.
These data, available in an online archive [Heays et al., 2022], are ideally suited for computing SO lifetimes against photodissociation in atmospheres and the interstellar medium, and the resulting likelihood of photolytic S-isotope fractionation.

§ ACKNOWLEDGEMENTS

The authors are very pleased to contribute to this special issue in honour of Prof. Wim Ubachs.
We have benefited enormously from his pioneering work in high-resolution and time-resolved spectroscopy and his insightful and energised collaboration on numerous projects.
We thank Bérenger Gans of the Institut des Sciences Moléculaires d'Orsay for permitting the use of the radio-frequency discharge source.
AH was funded by the NASA Postdoctoral Program through the NASA Astrobiology Institute, by grant number 19-03314S of the Czech Science Foundation, and the ERDF/ESF “Centre of Advanced Applied Sciences” (grant number CZ.02.1. 01/0.0/0.0/16_019/0000778).
JRL acknowledges support from the NASA Exobiology program (grant #80NSSC19K0475 to Arizona State University).

[Gottlieb and Ball, 1973]
authorC. A. Gottlieb
and authorJ. A.
Ball, journalAstrophys. J. Lett.
volume184, pagesL59 (year1973).

[Mateen et al., 2006]
authorM. Mateen,
authorP. Hofner, and
authorE. Araya,
journalAstrophys. J. Suppl. Ser. volume167,
pages239 (year2006).

[Guilloteau, S. et al., 2013]
authorGuilloteau, S.,
authorDi Folco, E.,
authorDutrey, A.,
authorSimon, M.,
authorGrosso, N., and
authorPiétu, V., journalAstron.
Astrophys. volume549, pagesA92

[Rivière-Marichalar, P. et al., 2019]
authorRivière-Marichalar, P.,
authorFuente, A.,
authorGoicoechea, J. R.,
authorPety, J., authorLe
Gal, R., authorGratier, P.,
authorGuzmán, V.,
authorRoueff, E.,
authorLoison, J. C.,
authorWakelam, V., et al.,
journalAstron. Astrophys. volume628,
pagesA16 (year2019).

[Lellouch et al., 1996]
authorE. Lellouch,
authorD. F. Strobel,
authorM. J. S. Belton,
authorM. E. Summers,
authorG. Paubert,
and authorR. Moreno,
journalAstrophys. J. Lett. volume459,
pagesL107 (year1996).

[Moullet et al., 2013]
authorA. Moullet,
authorE. Lellouch,
authorR. Moreno,
authorM. Gurwell,
authorJ. H. Black, and
authorB. Butler,
journalAstrophys. J. volume776,
pages32 (year2013).

[de Kleer et al., 2019]
authorK. de Kleer,
authorI. de Pater, and
authorM. Ádámkovics,
journalIcarus volume317, pages104

[Bockelée-Morvan et al., 2000]
authorD. Bockelée-Morvan,
authorD. C. Lis,
authorJ. E. Wink,
authorD. Despois,
authorJ. Crovisier,
authorR. Bachiller,
authorD. J. Benford,
authorN. Biver,
authorP. Colom,
authorJ. K. Davies,
et al., journalAstron. Astrophys.
volume353, pages1101 (year2000).

[Boissier, J. et al., 2007]
authorBoissier, J.,
authorBockelée-Morvan, D.,
authorBiver, N.,
authorCrovisier, J.,
authorDespois, D.,
authorMarsden, B. G., and
authorMoreno, R., journalAstron.
Astrophys. volume475, pages1131

[Na et al., 1990]
authorC. Y. Na,
authorL. W. Esposito,
and authorT. E.
Skinner, journalJ. Geophys. Res. – Atmos.
volume95, pages7485 (year1990).

[Belyaev et al., 2012]
authorD. A. Belyaev,
authorF. Montmessin,
authorJ.-L. Bertaux,
authorA. Mahieux,
authorA. A. Fedorova,
authorO. I. Korablev,
authorE. Marcq,
authorY. L. Yung, and
authorX. Zhang,
journalIcarus volume217, pages740
(year2012), noteadvances in Venus Science.

[Pavlov and Kasting, 2002]
authorA. Pavlov and
authorJ. Kasting,
journalAstrobiology volume2,
pages27 (year2002).

[Ono, 2017]
authorS. Ono, journalAnnu.
Rev. Earth Pl. Sc. volume45, pages301

[Farquhar et al., 2000]
authorJ. Farquhar,
authorH. Bao, and
authorM. Thiemens,
journalScience volume289,
pages756 (year2000).

[Lyons, 2009]
authorJ. R. Lyons,
journalChem. Geo. volume267,
pages164  (year2009), noteadvances in
experimental and theoretical isotope geochemistry.

[Danielache et al., 2014]
authorS. O. Danielache,
authorS. Tomoya,
authorA. Kondorsky,
authorI. Tokue, and
authorS. Nanbu, journalJ.
Chem. Phys. volume140, pages044319

[Martin, 1932]
authorE. V. Martin,
journalPhys. Rev. volume41,
pages167 (year1932).

[Clerbaux and Colin, 1994]
authorC. Clerbaux and
authorR. Colin, journalJ.
Mol. Spectrosc. volume165, pages334

[19]
authorR. Colin, journalCan.
J. Phys. volume47, pages979

[Liu et al., 2006]
authorC.-P. Liu,
authorN. L. Elliott,
authorC. M. Western,
authorY.-P. Lee, and
authorR. Colin, journalJ.
Mol. Spectrosc. volume238, pages213

[Elks and Western, 1999]
authorJ. M. F. Elks and
authorC. M. Western,
journalJ. Chem. Phys. volume110,
pages7699 (year1999).

[Archer et al., 2000]
authorC. P. Archer,
authorJ. M. F. Elks,
and authorC. M.
Western, journalJ. Chem. Phys.
volume112, pages6293 (year2000).

[Yamasaki et al., 2003]
authorK. Yamasaki,
authorF. Taketani,
authorS. Tomita,
authorK. Sugiura, and
authorI. Tokue, journalJ.
Phys. Chem. A volume107, pages2442

[Yamasaki et al., 2005]
authorK. Yamasaki,
authorS. Tomita,
authorT. Hatano,
authorF. Taketani, and
authorI. Tokue,
journalChem. Phys. Lett. volume413,
pages231  (year2005).

[Phillips, 1981]
authorL. F. Phillips,
journalJ. Phys. Chem. volume85,
pages3994 (year1981).

[Clyne and Liddy, 1982]
authorM. A. A. Clyne and
authorJ. P. Liddy,
journalJ. Chem. Soc. Farad. Trans. 2
volume78, pages1127 (year1982).

[Clyne and Tennyson, 1986]
authorM. A. A. Clyne and
authorP. H. Tennyson,
journalJ. Chem. Soc. Farad. Trans. 2
volume82, pages1315 (year1986).

[Lo et al., 1988]
authorG. Lo,
authorR. Beaman, and
authorD. Setser,
journalChem. Phys. Lett. volume149,
pages384  (year1988).

[Stuart et al., 1994]
authorB. C. Stuart,
authorS. M. Cameron,
and authorH. T.
Powell, journalJ. Phys. Chem.
volume98, pages11499 (year1994).

[Dixon et al., 1977]
authorR. Dixon,
authorP. Tasker, and
authorG. Balint-Kurti,
journalMol. Phys. volume34,
pages1455 (year1977).

[Swope et al., 1979]
authorW. C. Swope,
authorY.-P. Lee, and
authorH. F. Schaefer,
journalJ. Chem. Phys. volume71,
pages3761 (year1979).

[Ornellas and Borin, 1998]
authorF. R. Ornellas and
authorA. C. Borin,
journalMol. Phys. volume94,
pages139 (year1998).

[Borin and Ornellas, 1999]
authorA. C. Borin and
authorF. R. Ornellas,
journalChem. Phys. volume247,
pages351  (year1999).

[Yu and Bian, 2011]
authorL. Yu and
authorW. Bian, journalJ.
Comput. Chem. volume32, pages1577

[Sarka and Nanbu, 2019]
authorK. Sarka and
authorS. Nanbu, journalJ.
Phys. Chem. A volume123, pages3697

[Feng and Zhu, 2019]
authorY. Feng and
authorZ. Zhu, journalJ.
Quant. Spectrosc. Radiat. Transfer volume234,
pages98  (year2019).

[da Silva and Ballester, 2020]
authorR. S. da Silva
and authorM. Y.
Ballester, journalTheor. Chem. Acc.
volume139 (year2020).

[Nahon et al., 2012]
authorL. Nahon,
authorN. de Oliveira,
authorG. A. Garcia,
authorJ.-F. Gil,
authorB. Pilette,
authorO. Marcouille,
authorB. Lagarde, and
authorF. Polack,
journalJ. Synchrotron Rad. volume19,
pages508 (year2012).

[Heays et al., 2018]
authorA. N. Heays,
authorN. de Oliveira,
authorB. Gans,
authorK. Ito,
authorS. Boyé-Péronne,
authorS. Douin,
authorK. Hickson,
authorL. Nahon, and
authorJ. Loison,
journalJ. Quant. Spectrosc. Radiat. Transfer
volume204, pages12  (year2018).

[Stark et al., 2018]
authorG. Stark,
authorH. Herde,
authorJ. R. Lyons,
authorA. N. Heays,
authorN. de Oliveira,
authorG. Nave,
authorB. R. Lewis, and
authorS. T. Gibson,
journalJ. Chem. Phys. volume148,
pages244302 (year2018).

[de Oliveira et al., 2011]
authorN. de Oliveira,
authorM. Roudjane,
authorD. Joyeux,
authorD. Phalippou,
authorJ.-C. Rodier,
and authorL. Nahon,
journalNat. Photonics volume5,
pages149 (year2011).

[42]
authorN. de Oliveira,
authorD. Joyeux,
authorM. Roudjane,
authorJ.-F. Gil,
authorB. Pilette,
authorL. Archer,
authorK. Ito, and
authorL. Nahon,
journalJ. Synchrotron Rad. volume23,
pages887 (year2016).

[Freeman et al., 1984]
authorD. Freeman,
authorK. Yoshino,
authorJ. Esmond, and
authorW. Parkinson,
journalPlanet. Space Sci. volume32,
pages1125  (year1984).

[Stark et al., 1999]
authorG. Stark,
authorP. L. Smith,
authorJ. Rufus,
authorA. P. Thorne,
authorJ. C. Pickering,
and authorG. Cox,
journalJ. Geophys. Res. – Planet.
volume104, pages16585

[Heays et al., 2022]
authorA. N. Heays,
authorG. Stark,
authorJ. R. Lyons,
authorN. de Oliveira,
authorB. R. Lewis, and
authorS. T. Gibson,

[Western, 2017]
authorC. M. Western,
journalJ. Quant. Spectrosc. Radiat. Transfer
volume186, pages221  (year2017),
notesatellite Remote Sensing and Spectroscopy: Joint ACE-Odin
Meeting, October 2015.

[Brown et al., 1979]
authorJ. Brown,
authorE. Colbourn,
authorJ. Watson, and
authorF. Wayne, journalJ.
Mol. Spectrosc. volume74, pages294

[Cheung et al., 1986]
authorA.-C. Cheung,
authorK. Yoshino,
authorW. Parkinson,
and authorD. Freeman,
journalJ. Mol. Spectrosc. volume119,
pages1  (year1986).

[Brown and Merer, 1979]
authorJ. M. Brown and
authorA. J. Merer,
journalJ. Mol. Spectrosc. volume74,
pages488 (year1979).

[Lefebvre-Brion and Field, 2004]
authorH. Lefebvre-Brion
and authorR. W. Field,
titleThe spectra and dynamics of diatomic molecules
(publisherElsevier, Amsterdam, year2004).

[Hougen, 1970]
authorJ. T. Hougen,
titleThe calculation of rotational energy levels and
rotational line intensities in diatomic molecules
(publisherPhysics Laboratory Publications, NIST,

[Hansson and Watson, 2005]
authorA. Hansson and
authorJ. K. G. Watson,
journalJ. Mol. Spectrosc. volume233,
pages169  (year2005).

[Heays et al., 2011]
authorA. N. Heays,
authorG. D. Dickenson,
authorE. J. Salumbides,
authorN. de Oliveira,
authorD. Joyeux,
authorL. Nahon,
authorB. R. Lewis, and
authorW. Ubachs,
journalJ. Chem. Phys. volume135,
eid244301 (year2011).

[Johnson, 1977]
authorB. Johnson, journalJ.
Chem. Phys. volume67, pages4086

[Larsson, 1983]
authorM. Larsson,
journalAstron. Astrophys. volume128,
pages291 (year1983).

[Lattanzi et al., 2015]
authorV. Lattanzi,
authorG. Cazzoli, and
authorC. Puzzarini,
journalAstrophys. J. volume813,
pages4 (year2015).

[Le Roy, 1999]
authorR. J. Le Roy,
journalJ. Mol. Spectrosc. volume194,
pages189  (year1999).

[Colin, 1982]
authorR. Colin, journalJ.
Chem. Soc. Farad. Trans. 2 volume78,
pages1139 (year1982).

[Borin and Ornellas, 2000]
authorA. C. Borin and
authorF. R. Ornellas,
journalChem. Phys. Lett. volume322,
pages149  (year2000).

[Morse, 1929]
authorP. M. Morse,
journalPhys. Rev. volume34,
pages57 (year1929).

[Rees, 1947]
authorA. L. G. Rees,
journalProc. Phys. Soc. London volume59,
pages998 (year1947).

[Lewis et al., 2001]
authorB. R. Lewis,
authorS. T. Gibson,
authorF. T. Hawes, and
authorL. Torop,
journalPhys. Chem. Earth volume26,
pages519 (year2001).

[Lewis et al., 2018]
authorB. R. Lewis,
authorS. T. Gibson,
authorG. Stark, and
authorA. N. Heays,
journalJ. Chem. Phys. volume148,
pages244303 (year2018).

[Julienne and Krauss, 1975]
authorP. S. Julienne and
authorM. Krauss,
journalJ. Mol. Spectrosc. volume56,
pages270 (year1975).

[Lewis et al., 1994]
authorB. R. Lewis,
authorS. T. Gibson,
and authorP. M.
Dooley, journalJ. Chem. Phys.
volume100, pages7012 (year1994).


# A master equation incorporating the system-environment correlations present
in the joint equilibrium state

Ali Raza Mirza  School of Science & Engineering, Lahore University of
Management Sciences (LUMS), Opposite Sector U, D.H.A, Lahore 54792, Pakistan
Muhammad Zia  School of Science & Engineering, Lahore University of Management
Sciences (LUMS), Opposite Sector U, D.H.A, Lahore 54792, Pakistan  Adam Zaman
Chaudhry<EMAIL_ADDRESS>School of Science & Engineering, Lahore
University of Management Sciences (LUMS), Opposite Sector U, D.H.A, Lahore
54792, Pakistan

###### Abstract

We present a general master equation, correct to second order in the system-
environment coupling strength, that takes into account the initial system-
environment correlations. We assume that the system and its environment are in
a joint thermal equilibrium state, and thereafter a unitary operation is
performed to prepare the desired initial system state, with the system
Hamiltonian possibly changing thereafter as well. We show that the effect of
the initial correlations shows up in the second-order master equation as an
additional term, similar in form to the usual second-order term describing
relaxation and decoherence in quantum systems. We apply this master equation
to a generalization of the paradigmatic spin-boson model, namely a collection
of two-level systems interacting with a common environment of harmonic
oscillators, as well as a collection of two-level systems interacting with a
common spin environment. We demonstrate that, in general, the initial system-
environment correlations need to be accounted for in order to accurately
obtain the system dynamics.

###### pacs:

03.65.Yz, 05.30.-d, 03.67.Pp, 42.50.Dv

##  I Introduction

The study of quantum systems interacting with their environment, namely open
quantum systems, has gained immense importance lately. This is both due to
practical reasons as well as the relevance of open quantum systems to the
foundations of quantum mechanics [1]. In particular, precise quantum coherent
control and state preparation necessary for quantum computation and
information is only possibly by taking the interaction of a quantum system
with its environment, while the process of decoherence sheds light on the so-
called measurement problem. Amongst the variety of tools and techniques
developed to tackle the dynamics of open quantum systems, the most popular
approach is the quantum master equation approach whereby a differential
equation for the time evolution of the system state is obtained and solved.
The basic approach is that an initial state of the system and the environment
is assumed, the total system-environment Hamiltonian is written down, and then
using this total Hamiltonian, the time evolution of the total system-
environment state is examined. Since we are typically interested only in the
system dynamics, the environment is traced out to obtain the master equation.
Unfortunately, performing this process to obtain the master equation for most
realistic system-environment models involves making a series of
approximations. For example, the system-environment coupling strength is
assumed to be weak so that the joint time evolution operator can be found
perturbatively [1, 2]. The environment is assumed to have a short ‘memory
time’ (the Markovian approximation) - this means that the environment loses
information about the system very quickly [3, 4]. Finally, the initial system-
environment state is assumed to a simple product state, with the system and
the environment independent of one another [3].

With increasingly sophisticated quantum technologies, each of the assumptions
typically made in the derivation of master equations has come under renewed
reexamination. Master equations that allow one to deal with stronger system-
environment coupling strengths have been formulated (see, for example, Ref.
[5]). Measures of non-Markovianity have also been put forward [4]. Most
pertinent for us, the role of the initial system-environments - the
correlations present in the total system-environment state at the initial
state - has been investigated widely [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36]. Such studies have generally been performed using exactly solvable models
such as the pure dephasing model of a two-level system interacting with a
collection of harmonic oscillators [24, 26]. Some efforts have nevertheless
been made to consider master equations that include the effect of the initial
correlations. For example, in Ref. [27], the system and its environment were
allowed to come to thermal equilibrium and thereafter a projective measurement
was performed on the system to prepare the desired initial system state. It
was shown that the effect of the initial correlations appears as an additional
term in the second-order master equation, similar in form to the first term in
the master equation that describes the effect of the free system-Hamiltonian.
This approach was later generalized to higher order system-environment
coupling strengths [34]. Along similar lines, in this paper, we consider the
quantum system and its environment to reach a joint equilibrium state. A
unitary operator is then performed on the system to prepare (approximately) a
desired initial system state, and a time-local master equation, correct to
second-order in the system-environment coupling strength, that describes the
ensuing dynamics of the system is derived. In fact, the system Hamiltonian
before the application of the system Hamiltonian can be different from the
system Hamiltonian after the unitary operator - the former plays a role in the
initial state preparation, while the latter plays a role in the dynamics
thereafter. We show that the effect of the initial correlations is now again
contained in an additional term in the master equation, but now the form of
this additional term is similar to that of the second term in the usual master
equation that describes relaxation and decoherence of open quantum systems. We
then apply our master equation to a collection of two-level systems
interacting with a common environment of harmonic oscillators. We work out the
additional term in the master equation, and perform numerical simulations to
show that the effect of the initial correlations increases as the number of
two-level systems increases. Along similar lines, we also apply our master
equation to analyze the effect of initial correlations for a collection of
two-level systems interacting with a spin environment.

This paper is organized as follows. In Sec. II, we derive our general time-
local second-order master equation. Sec. III discusses the application of this
master equation to the large spin-boson model, while Sec. IV applies the
master equation to a collection of two-level systems interacting with a spin
environment. We then conclude in Sec. V. The appendices consist of some
technical details regarding the derivation of the usual relaxation term in the
master equation, as well as the details of the exactly solvable limit of the
large spin-boson model.

##  II The formalism

We start by briefly discussing the problem we wish to solve. We are given a
quantum system which is interacting with its environment and has reached a
joint equilibrium state. At the initial time, a unitary operation is performed
on the system alone, and the system Hamiltonian itself may also be changed.
After all, one can, for example, apply a large magnetic field in order to
prepare the initial state, and thereafter alter the value of this magnetic
field to generate some desired system dynamics. Our problem is to derive the
master equation, correct to second-order in the system-environment coupling
strength, that describes the system dynamics. We then write the system-
environment Hamiltonian as

$\displaystyle H_{tot}=\begin{cases}H_{S0}+H_{B}+\alpha V&t\leq 0,\\\
H_{S}+H_{B}+\alpha V&t\geq 0.\\\ \end{cases}$ (1)

Here $H_{S}$ is the system Hamiltonian corresponding to coherent evolution of
the system only after the initial time $t=0$ at which the system state is
prepared. $H_{\text{S0}}$ is similar to $H_{S}$ but with possibly different
parameters in order to prepare the initial system state. $H_{B}$ is the
Hamiltonian of the environment, and $V$ corresponds to the system-environment
coupling. $\alpha$ is simply a dimensionless parameter introduced to keep
track of the perturbation order; at the end of the calculation, we will set
$\alpha=1$. Let us now discuss in detail the initial state preparation.

###  II.1 Initial state preparation

We let our system to come to a joint equilibrium state with the environment.
What we mean by this is that the equilibrium state of the system is not simply
proportional to $e^{-\beta H_{S0}}$ \- there are corrections due to the finite
system-environment coupling strength. We instead consider the system and the
environment together in the thermal equilibrium state proportional to
$e^{-\beta H_{\text{tot}}}$ with $H_{\text{tot}}=H_{S0}+H_{B}+\alpha V$; the
system state can then be obtained by simply tracing out the environment. A
unitary operator $\Omega$ is then applied to the system. The initial system-
environment state is then

$\displaystyle\rho_{\text{tot}}(0)$ $\displaystyle=\frac{\Omega e^{-\beta
H_{\text{tot}}}\Omega^{\dagger}}{Z_{\text{tot}}},$ (2)

with $Z_{\text{tot}}=\text{Tr}_{\text{S,B}}\left[e^{-\beta
H_{\text{tot}}}\right]$ the partition function and $\text{Tr}_{S,B}$ denotes
the trace over the system and the environment. Now, assuming the system-
environment coupling strength to be weak, we can use the Kubo identity to
expand the joint state given by Eq. (2). The Kubo identity tells us that for
any two arbitrary operators $X$ and $Y$,

$\displaystyle e^{\beta(X+Y)}$ $\displaystyle=e^{\beta
X}\left[1+\alpha\int_{0}^{\beta}e^{-\lambda
X}Y^{\lambda(X+Y)}d\lambda\right].$ (3)

By setting $X=-(H_{\text{S0}}+H_{B})$ and $Y=-V$, and using the Kubo identity
twice, we find that to second order in the system-environment coupling
strength

$\displaystyle
e^{-\beta(H_{\text{S0}}+H_{B}+V)}=e^{-\beta(H_{\text{S0}}+H_{B})}$
$\displaystyle-\alpha
e^{-\beta(H_{\text{S0}}+H_{B})}\int_{0}^{\beta}e^{\lambda(H_{\text{S0}}+H_{B})}V\times
e^{-\lambda(H_{\text{S0}}+H_{B})}d\lambda$
$\displaystyle+\alpha^{2}e^{-\beta(H_{\text{S0}}+H_{B})}\int_{0}^{\beta}d\lambda
e^{\lambda(H_{\text{S0}}+H_{B})}V$ $\displaystyle\times
e^{-\lambda(H_{\text{S0}}+H_{B})}\int_{0}^{\lambda}e^{\lambda^{\prime}(H_{\text{S0}}+H_{B})}Ve^{-\lambda^{\prime}(H_{\text{S0}}+H_{B})}d\lambda^{\prime}.$
(4)

We now write the system environment coupling $V$ as $F\otimes B$, where $F$
and $B$ are operators living in the system and environment Hilbert space
respectively. The extension to the more general case where
$V=\sum_{\alpha}F_{\alpha}\otimes B_{\alpha}$ is straightforward. Eq. (4) can
then be simplified as

$\displaystyle
e^{-\beta(H_{\text{S0}}+H_{B}+V)}=e^{-\beta(H_{\text{S0}}+H_{B})}-\alpha
e^{-\beta(H_{\text{S0}}+H_{B})}$
$\displaystyle\times\int_{0}^{\beta}F(\lambda)\otimes
B(\lambda)d\lambda+\alpha^{2}e^{-\beta(H_{\text{S0}}+H_{B})}$
$\displaystyle\times\int_{0}^{\beta}d\lambda F(\lambda)\otimes
B(\lambda)\int_{0}^{\lambda}F(\lambda^{\prime})\otimes
B(\lambda^{\prime})d\lambda^{\prime}.$ (5)

where $F\left(\lambda\right)=e^{\lambda H_{\text{S0}}}Fe^{-\lambda
H_{\text{S0}}}$ and $B\left(\lambda\right)=e^{\lambda H_{B}}Be^{-\lambda
H_{B}}$. We now use this in Eq. (2) and thereafter take the trace over the
environment to find the initial system state correct to second order in the
system-environment coupling strength. This is important because our aim is to
derive a master equation correct to second-order in the system-environment
strength. For consistency, the initial system state used to solve this master
equation should also be accurate to second-order in the system-environment
coupling strength. For ease of notation, we write the initial system state as

$\displaystyle\rho(0)=\rho^{(0)}(0)+\rho^{(1)}(0)+\rho^{(2)}(0),$ (6)

where,

$\displaystyle\rho^{(0)}(0)=\frac{\text{Tr}_{B}\left[\Omega\left(e^{-\beta(H_{\text{S0}}+H_{B})}\right)\Omega^{\dagger}\right]}{Z_{\text{tot}}}$
(7) $\displaystyle\rho^{(1)}(0)$
$\displaystyle=\frac{\text{Tr}_{B}\left[-\alpha\Omega\left(e^{-\beta(H_{\text{S0}}+H_{B})}\int_{0}^{\beta}F(\lambda)\otimes
B(\lambda)d\lambda\right)\Omega^{\dagger}\right]}{Z_{\text{tot}}}$ (8)
$\displaystyle\rho^{(2)}(0)=\frac{1}{{Z_{\text{tot}}}}\times\text{Tr}_{B}\Bigg{[}\alpha^{2}\Omega\Big{(}e^{-\beta(H_{\text{S0}}+H_{B})}$
$\displaystyle\times\int_{0}^{\beta}d\lambda F(\lambda)\otimes
B(\lambda)\int_{0}^{\lambda}F(\lambda^{\prime})\otimes
B(\lambda^{\prime})d\lambda^{\prime}\Big{)}\Omega^{\dagger}\Bigg{]}.$ (9)

Let us simplify these relations one by one. $\rho^{(0)}(0)$ can be simplified
as,

$\displaystyle\rho^{(0)}(0)=\frac{e^{-\beta\widetilde{H}_{\text{S0}}}Z_{B}}{Z_{\text{tot}}},$

where $\widetilde{H}_{\text{S0}}=\Omega H_{\text{S0}}\Omega^{\dagger}$ and
$Z_{B}=\text{Tr}_{B}\left[e^{-\beta H_{B}}\right]$. As for $\rho^{(1)}(0)$, we
can write

$\displaystyle\rho^{(1)}(0)=\frac{-\alpha Z_{B}\int_{o}^{\beta}\Omega
e^{-\beta
H_{\text{S0}}}F(\lambda)\Omega^{\dagger}\big{\langle}B(\lambda)\big{\rangle}_{B}d\lambda}{Z_{\text{tot}}},$

where $\langle\ldots\rangle_{B}=\text{Tr}_{B}[e^{-\beta
H_{B}}(\ldots)/Z_{B}]$. Since $\langle B(\lambda)\rangle_{B}$ is zero for most
system-environment models, we simply get that $\rho^{(1)}(0)=0$. Carrying on,
$\rho^{(2)}(0)$ can be simplified as

$\displaystyle\rho^{(2)}(0)=$ $\displaystyle\frac{\alpha^{2}Z_{B}\Omega
e^{-\beta
H_{\text{S0}}}\int_{0}^{\beta}\int_{0}^{\lambda}F(\lambda)F(\lambda^{\prime})\Omega^{\dagger}\big{\langle}B(\lambda)B(\lambda^{\prime})\big{\rangle}_{B}d\lambda^{\prime}d\lambda}{Z_{\text{tot}}}.$

To proceed further, we evaluate the partition function $Z_{\text{tot}}$. We
note that $Z_{\text{tot}}$ has to be such that the trace of the system state
$\rho(0)$ in Eq. (6) is one. It is then clear that

$\displaystyle Z_{\text{tot}}=$ $\displaystyle
Z_{B}\text{Tr}_{\text{S}}\left[e^{-\beta
H_{\text{S0}}}\right]+\alpha^{2}Z_{B}\text{Tr}_{\text{S}}\Bigg{[}\Omega
e^{-\beta H_{\text{S0}}}$
$\displaystyle\times\int_{0}^{\beta}\int_{0}^{\lambda}F(\lambda)F(\lambda^{\prime})\Omega^{\dagger}\big{\langle}B(\lambda)B(\lambda^{\prime})\big{\rangle}_{B}d\lambda^{\prime}d\lambda\Bigg{]}.$

It is then clear that the initial system density matrix can be written as

$\displaystyle\rho(0)$
$\displaystyle=\frac{e^{-\beta\widetilde{H}_{S0}}}{Z_{S0}Z^{\prime}}\Big{[}\mathds{1}+\int_{0}^{\beta}\int_{0}^{\lambda}F^{R}(\lambda)F^{R}(\lambda^{\prime})$
$\displaystyle\times\big{\langle}B(\lambda)B(\lambda^{\prime})\big{\rangle}_{B}d\lambda^{\prime}d\lambda\Big{]}$
(10)

where $F^{R}(\lambda)=\Omega F(\lambda)\Omega^{\dagger}$,
$Z_{S0}=\text{Tr}_{S}[e^{-\beta H_{S0}}]$, and

$\displaystyle Z^{\prime}=1+\int_{0}^{\beta}\int_{0}^{\lambda}\langle
F(\lambda)F(\lambda^{\prime})\rangle_{S}\langle
B(\lambda)B(\lambda^{\prime})\rangle_{B}d\lambda^{\prime}d\lambda,$ (11)

with $\langle\ldots\rangle_{S}=\text{Tr}_{S}[e^{-\beta
H_{S0}}(\ldots)/Z_{S0}]$. With the initial system state, correct to second
order in the system-environment coupling strength available, we now turn our
attention to deriving the second-order master equation.

###  II.2 Derivation of the master equation

We now derive a master equation that describes the time evolution of the
system for $t>0$. The system-environment Hamiltonian is

$\displaystyle H_{\text{tot}}=H_{\text{S}}+H_{\text{B}}+\alpha V\equiv
H_{0}+\alpha V.$

Note that the system Hamiltonian $H_{S}$ can be different from the previous
system Hamiltonian $H_{S0}$. Using perturbation theory, the unitary time
evolution for such an Hamiltonian can be written as,

$\displaystyle U(t)\approx
U_{0}\left(t\right)\left[1-\alpha\int_{0}^{t}U_{0}^{\dagger}(s)VU_{0}(s)\,ds\right]$
(12)

where, $U_{0}(t)\equiv U_{S}(t)\otimes U_{B}(t)$ is the free unitary time
evolution corresponding to $H_{0}.$ The matrix elements of the system density
matrix can be written as
$\rho_{mn}\left(t\right)=\text{Tr}_{S}\left[\left|n\right\rangle\left\langle
m\right|\rho\left(t\right)\right]$, where $\left|n\right\rangle$ and
$\left|m\right\rangle$ are some basis states of the system. Since
$\rho(t)=\text{Tr}_{B}\rho_{\text{tot}}(t)$, we can alternatively write

$\displaystyle\rho_{mn}(t)=\text{Tr}_{\text{S,B}}\left[X_{nm}^{H}\left(t\right)\rho_{\text{tot}}(0)\right],$

where $X_{nm}^{H}(t)=U^{\dagger}(t)(\left|n\right\rangle\left\langle
m\right|\otimes\mathds{1}_{B})U(t)$. The master equation can then be written
in the general form

$\displaystyle\frac{d}{dt}\rho_{mn}(t)$
$\displaystyle=\text{Tr}_{\text{S,B}}\left[\rho_{\text{tot}}(0)\frac{d}{dt}X_{nm}^{H}(t)\right].$
(13)

To make further progress, we note that $X_{nm}^{H}(t)$ is a Heisenberg picture
operator. Using the Heisenberg equation of motion and Eq. (12), it can be
shown that, correct to second order in the system-environment coupling
strength,

$\displaystyle\frac{d}{dt}X_{nm}^{H}(t)$
$\displaystyle=i\left[H_{0}^{H}(t),X_{nm}^{H}(t)\right]+i\alpha\left[\tilde{V}(t),\tilde{X}_{nm}(t)\right]$
$\displaystyle+\alpha^{2}\int_{0}^{t}ds\left[\left[\tilde{V}(t),\tilde{X}_{nm}(t)\right],\tilde{V}(s)\right],$
(14)

where the tildes denote time evolution under the ‘free’ unitary operator
$U_{0}\left(t\right)$ while the superscript ‘$H$’ denotes time evolution with
the full time evolution operator. Using Eq. (II.2) and given the initial
system-environment state $\rho_{\text{tot}}(0)$ in Eq. (6), we can derive our
master equation that describes the time evolution of the quantum system by
simplifying Eq. (13). The first term is very straightforward. We simply have
that

$\displaystyle\text{Tr}_{\text{S,B}}\left[\rho_{\text{tot}}(0)i\left[H_{0}^{H}(t),X_{nm}^{H}(t)\right]\right]$
$\displaystyle=i\text{Tr}_{\text{S,B}}\left[\rho_{\text{tot}}(t)\left[H_{S}+H_{B},(\left|n\right\rangle\left\langle
m\right|\otimes\mathds{1}_{B})\right]\right]$
$\displaystyle=i\text{Tr}_{\text{S}}\left[\rho(t)\left[H_{\text{S}},\left|n\right\rangle\left\langle
m\right|\right]\right]$ $\displaystyle=i\left\langle
m\right|\left[\rho(t),H_{\text{S}}\right]\left|n\right\rangle.$ (15)

This term simply tells us about free system evolution corresponding to the
system Hamiltonian $H_{\text{S}}$.

To calculate the next term in our master equation, that is,

$\displaystyle
i\alpha\text{Tr}_{\text{S,B}}\left[\rho_{\text{tot}}(0)\left[V(t),X_{nm}(t)\right]\right],$

it is useful to write the initial system-environment state as
$\rho_{\text{tot}}(0)=\rho_{\text{tot}}^{(0)}+\rho_{\text{tot}}^{(1)}$, where,
using Eq. (4), it should be clear that

$\displaystyle\rho_{\text{tot}}^{(0)}(0)$ $\displaystyle=\frac{\Omega
e^{-\beta(H_{\text{S0}}+H_{B})}\Omega^{\dagger}}{Z_{\text{tot}}}=\widetilde{\rho}_{S0}\otimes\rho_{B}$
(16) $\displaystyle\rho_{\text{tot}}^{(1)}(0)$
$\displaystyle=\frac{-\alpha\Omega
e^{-\beta(H_{\text{S0}}+H_{B})}Q_{\text{SB}}(\beta)\Omega^{\dagger}}{Z_{\text{tot}}},$
(17)

with $\widetilde{\rho}_{S0}=e^{-\beta\widetilde{H}_{S0}}/Z_{S0}$,
$\rho_{B}=e^{-\beta H_{B}}/Z_{B}$, the partition function
$Z_{\text{tot}}=Z_{S0}Z_{B}$ and
$Q_{\text{SB}}(\beta)=\int_{0}^{\beta}d\lambda F(\lambda)\otimes B(\lambda)$.
We do not need the higher order terms since there is already a factor of
$\alpha$ in $i\alpha\left[V(t),X_{nm}(t)\right]$. Now, the contribution of
$\rho_{\text{tot}}^{(0)}$ is

$\displaystyle
i\alpha\text{Tr}_{\text{S,B}}\left[\rho_{\text{tot}}^{(0)}(0)\left[U_{0}^{\dagger}(t)VU_{0}(t),U_{0}^{\dagger}(t)X_{nm}U_{0}(t)\right]\right]$
$\displaystyle=i\alpha\text{Tr}_{\text{S,B}}\left[\widetilde{\rho}_{S0}\otimes\rho_{B}U_{0}^{\dagger}(t)\left[F\otimes
B,\left|n\right\rangle\left\langle
m\right|\otimes\mathds{1}_{B}\right]U_{0}(t)\right]$
$\displaystyle=i\alpha\text{Tr}_{\text{S}}\left[\tilde{\rho}_{S0}U^{\dagger}(t)\left[F,Y_{nm}\right]U(t)\right]\times\langle
B(t)\rangle_{B}.$

Again, since $\langle B(t)\rangle_{B}$ is usually zero for most system-
environment models, this contribution turns out to be zero.

The most interesting contribution is due to $\rho_{\text{tot}}^{(1)}(0)$. This
is

$\displaystyle
i\alpha\text{Tr}_{\text{S,B}}\Big{[}\rho_{\text{tot}}^{(1)}(0)\left[U_{0}^{\dagger}(t)VU_{0}(t),U_{0}^{\dagger}(t)X_{nm}U_{0}(t)\right]\Big{]}$
$\displaystyle=\frac{-i\alpha^{2}}{Z_{\text{S0}}}\int_{0}^{\beta}\text{Tr}_{\text{S,B}}\Big{[}\rho_{B}{\Omega
e^{-\beta H_{\text{S0}}}F(\lambda)\Omega^{\dagger}\otimes
B(\lambda)}U_{S}^{\dagger}(t)$
$\displaystyle\times\left[F,\left|n\right\rangle\left\langle
m\right|\right]U_{S}(t)U_{B}^{\dagger}(t)BU_{B}(t)\Big{]}d\lambda$
$\displaystyle=\frac{-i\alpha^{2}}{Z_{\text{S0}}}\int_{0}^{\beta}\text{Tr}_{\text{S}}\Big{[}{\Omega
e^{-\beta
H_{\text{S0}}}F(\lambda)\Omega^{\dagger}}U_{S}^{\dagger}(t)\left[F,\left|n\right\rangle\left\langle
m\right|\right]$ $\displaystyle\times
U_{S}(t)\Big{]}\text{Tr}_{B}\Big{[}\rho_{B}B(\lambda)B(t)\Big{]}d\lambda$
$\displaystyle=\frac{-i\alpha^{2}}{Z_{\text{S0}}}\int_{0}^{\beta}\left\langle
m\right|\left[U_{S}(t)\Omega e^{-\beta
H_{\text{S0}}}F(\lambda)\Omega^{\dagger}U_{S}^{\dagger}(t),F\right]\left|n\right\rangle$
$\displaystyle\times B_{\text{corr}}(\lambda,t)\,d\lambda,.$ (18)

where we have defined $B(t)=U_{B}^{\dagger}(t)BU_{B}(t)$ and
$B_{\text{corr}}(\lambda,t)=\text{Tr}_{B}\Big{[}\rho_{B}B(\lambda)B(t)\Big{]}$.
This is additional term in the master equation that takes into account the
effect of the initial correlations, correct to second order in the system-
environment coupling strength. In basis-independent form, we can write this
term as

$\displaystyle-i\left[\widetilde{\rho}(t)J^{R}_{\text{corr}}(\beta,t),F\right],$
(19)

where we have defined

$\displaystyle J^{R}_{\text{corr}}(\beta,t)$
$\displaystyle=\int_{0}^{\beta}\overleftarrow{F}^{R}(\lambda,t)\text{B}_{\text{corr}}(\lambda,t)d\lambda,$
(20) $\displaystyle\overleftarrow{F}^{R}(\lambda,t)$
$\displaystyle=U_{S}(t)\Omega e^{\lambda H_{\text{S0}}}Fe^{-\lambda
H_{\text{S0}}}\Omega^{\dagger}U_{S}^{\dagger}(t).$ (21)

We now write

$\displaystyle-i\left[\widetilde{\rho}(t)J^{R}_{\text{corr}}(\beta,t),F\right]=-\frac{i}{2}\left(\left[\widetilde{\rho}(t)J^{R}_{\text{corr}}(\beta,t),F\right]-\text{h.c.}\right).$
(22)

This is permitted because
$-i\left[\widetilde{\rho}(t)J^{R}_{\text{corr}}(\beta,t),F\right]$ is
hermitian, so $\left[\widetilde{\rho}(t)J^{R}_{\text{corr}}(\beta,t),F\right]$
is anti-hermitian. We now replace $\widetilde{\rho}(t)$ by $\rho(t)$. This
step is permitted since the corrections are of order higher than the second-
order master equation that we are considering. Consequently, the term in the
master equation that takes into account the initial correlations is
$-\frac{i}{2}\left(\left[\rho(t)J^{R}_{\text{corr}}(\beta,t),F\right]-\text{h.c.}\right)$.

We next simplify the contribution of the third term in Eq. (II.2). It is clear
that now only $\rho_{\text{tot}}^{(0)}(0)$ contributes. Similar manipulations
to those performed above lead to (see the appendix for details)

$\alpha^{2}\int_{0}^{t}\left\langle
m\right|\Big{(}\left[\bar{F}(t,s)\widetilde{\rho}_{S}(t),F\right]C_{ts}+\text{h.~{}c.}\Big{)}\left|n\right\rangle\,ds,$

where the environment correlation function is $C_{ts}=\langle
B(t)B(s)\rangle_{B}$,
$\widetilde{\rho}_{S}(t)=U_{S}(t)\widetilde{\rho}_{S0}U_{S}^{\dagger}(t)$,
$\bar{F}(t,s)=U_{S}^{\dagger}(t,s)FU_{S}(t,s)$, and h. c. denotes the
hermitian conjugate. We can further replace $\widetilde{\rho}_{S}(t)$ by
$\rho(t)$ to get

$\alpha^{2}\int_{0}^{t}\left\langle
m\right|\Big{(}\left[\bar{F}(t,s)\rho(t),F\right]C_{ts}+\text{h.~{}c.}\Big{)}\left|n\right\rangle\,ds.$

Once again, this is permitted since the corrections lead to terms of higher
order in the master equation (compared to the second order master equation
that we are considering). We now put all the terms together to arrive at the
general basis-independent form of the master equation given by

$\displaystyle\frac{d}{dt}\rho(t)=i\left[\rho(t),H_{S}\right]-\frac{i}{2}\left(\left[\rho(t)J^{R}_{\text{corr}}(\beta,t),F\right]-\text{h.c.}\right)$
$\displaystyle+\int_{0}^{t}\left(\left[\bar{F}(t,s)\rho(t),F\right]C_{ts}+\text{h.c.}\right)\,ds.$
(23)

##  III Application to the large spin-boson model

We will now apply our derived master equation to a variant of the paradigmatic
spin-boson model [1] with a large number of two-level systems interacting with
a common environment of harmonic oscillators [37, 26, 27]. Recall that the
total system-environment Hamiltonian is given by
$H_{\text{tot}}=H_{\text{S0}}+H_{\text{B}}+V$ for $t<0$, while
$H_{\text{tot}}=H_{S}+H_{B}+V$ for $t\geq 0$. For the large spin-boson model,
we write

$\displaystyle H_{\text{S0}}$
$\displaystyle=\varepsilon_{0}J_{z}+\Delta_{0}J_{x}$ (24) $\displaystyle
H_{S}$ $\displaystyle=\varepsilon J_{z}+\Delta J_{x}$ (25) $\displaystyle
H_{B}$ $\displaystyle=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k},$ (26)
$\displaystyle V$
$\displaystyle=J_{z}\sum_{k}\big{(}g_{k}^{*}b_{k}+g_{k}b_{k}^{\dagger}\big{)}$
(27)

where $J_{x},J_{y},J_{z}$ are the collective spin operators with
$J^{2}=J_{x}^{2}+J_{y}^{2}+J_{z}^{2},\varepsilon$ is energy bias, $\Delta$ is
the tunneling amplitude, $H_{B}$ is the bath of harmonic oscillators (we are
ignoring zero point energy), while $V$ describes the interaction between the
common harmonic oscillator bath and the spin system. We have set $\hbar=1$
throughout and the values of other parameters will be in dimensionless units.
Note that the system operator $F=J_{z},$ and the bath operator
$B=\sum_{k}\big{(}g_{k}^{*}b_{k}+g_{k}b_{k}^{\dagger}\big{)}.$ One imagine
that the large-spin system has been interacting with the environment for a
long time with a relatively large value of $\varepsilon_{0}$ and a small value
of $\Delta_{0}$. In such a situation, the state of the system will be
approximately corresponding to the state with all spins down in the
$z$-direction. At time $t=0$, we then apply a unitary operator to prepare the
desired initial state. For example, if the desired initial state is one with
all spins in the $x$-direction. then the unitary operator that should be
applied is $U_{R}=e^{i\pi J_{y}/2}$. In other words, a $\frac{\pi}{2}$-pulse
is used to prepare the initial system state, with the assumption that this
pulse takes a very short time to be applied. With the initial state
approximately prepared, we can then change the parameters to the system
Hamiltonian to whatever values we desire to generate any required system
evolution. Let us look at how the initial system-environment correlations
appear in this system evolution using our general master equation.

Our objective is to calculate the operator $J^{R}_{\text{corr}}$. To do so, we
first find

$\displaystyle\overleftarrow{F}^{R}(\lambda,t)=U(t)\left[U_{R}\left(e^{\lambda
H_{S0}}Fe^{-\lambda H_{S0}}\right)U^{\dagger}_{R}\right]U^{\dagger}(t)$
$\displaystyle=J_{x}\left[a_{x}d_{x}+a_{y}c_{x}-a_{z}b_{x}\right]+J_{y}\left[a_{x}d_{y}+a_{y}c_{y}-a_{z}b_{y}\right]$
$\displaystyle+J_{z}\left[a_{x}d_{z}+a_{y}c_{z}-a_{z}b_{z}\right],$

with

$\displaystyle a_{x}$
$\displaystyle=\frac{\varepsilon_{0}\Delta_{0}}{\Delta^{\prime
2}}\left\\{1-\cosh\left({\lambda\Delta^{\prime}}\right)\right\\},$
$\displaystyle a_{y}$
$\displaystyle=\frac{-i\Delta_{0}}{\Delta^{\prime}}\sinh\left(\lambda\Delta^{\prime}\right),$
$\displaystyle a_{z}$
$\displaystyle=\frac{\varepsilon_{0}^{2}+\Delta^{2}_{0}\cosh\left({\lambda\Delta^{\prime}}\right)}{\Delta^{\prime
2}},$ $\displaystyle b_{x}$
$\displaystyle=\frac{\Delta^{2}+\varepsilon^{2}\cos\left({\widetilde{\Delta}t}\right)}{\widetilde{\Delta}^{2}},$
$\displaystyle b_{y}$
$\displaystyle=\frac{\varepsilon}{\widetilde{\Delta}}\sin\left(\widetilde{\Delta}t\right),$
$\displaystyle b_{z}$
$\displaystyle=\frac{\varepsilon\Delta}{\widetilde{\Delta}^{2}}\left\\{1-\cos\left({\widetilde{\Delta}t}\right)\right\\},$
$\displaystyle c_{x}$
$\displaystyle=-\frac{\varepsilon}{\widetilde{\Delta}}\sin{\left(\widetilde{\Delta}t\right)},$
$\displaystyle c_{y}$ $\displaystyle=\cos{\left(\widetilde{\Delta}t\right)},$
$\displaystyle c_{z}$
$\displaystyle=\frac{\Delta}{\widetilde{\Delta}}\sin{\left(\widetilde{\Delta}t\right)},$
$\displaystyle d_{x}$
$\displaystyle=\frac{\varepsilon\Delta}{\tilde{\Delta}^{2}}\left\\{1-\cos{\left(\widetilde{\Delta}t\right)}\right\\},$
$\displaystyle d_{y}$
$\displaystyle=-\frac{\Delta}{\widetilde{\Delta}}\sin{\left(\widetilde{\Delta}t\right)},$
$\displaystyle d_{z}$
$\displaystyle=1+\frac{\Delta^{2}}{\tilde{\Delta}^{2}}\left\\{\cos{\left(\widetilde{\Delta}t\right)}-1\right\\}.$

Here $\Delta^{\prime 2}=\varepsilon_{0}^{2}+\Delta_{0}^{2}$ and
$\widetilde{\Delta}^{2}=\varepsilon^{2}+\Delta^{2}$. In short,

$\displaystyle\overleftarrow{F}^{R}(\lambda,t)=\alpha_{1}(\lambda,t)J_{x}+\alpha_{2}(\lambda,t)J_{y}+\alpha_{3}(\lambda,t)J_{z},$
(28)

where,

$\displaystyle\alpha_{1}(\lambda,t)$
$\displaystyle=a_{x}d_{x}+a_{y}c_{x}-a_{z}b_{x},$
$\displaystyle\alpha_{2}(\lambda,t)$
$\displaystyle=a_{x}d_{y}+a_{y}c_{y}-a_{z}b_{y},$
$\displaystyle\alpha_{3}(\lambda,t)$
$\displaystyle=a_{x}d_{z}+a_{y}c_{z}-a_{z}b_{z}.$

It then follows that

$\displaystyle J_{\text{corr}}^{R}(\beta,t)$
$\displaystyle=P(\beta,t)J_{x}+Q(\beta,t)J_{y}+R(\beta,t)J_{z},$ (29)

with,

$\displaystyle P(\beta,t)$
$\displaystyle=\int_{0}^{\beta}\alpha_{1}(\lambda,t)B_{\text{corr}}(\lambda,t)d\lambda$
$\displaystyle Q(\beta,t)$
$\displaystyle=\int_{0}^{\beta}\alpha_{2}(\lambda,t)B_{\text{corr}}(\lambda,t)d\lambda$
$\displaystyle R(\beta,t)$
$\displaystyle=\int_{0}^{\beta}\alpha_{3}(\lambda,t)B_{\text{corr}}(\lambda,t)d\lambda.$

We now calculate bath correlation term $B_{\text{corr}}(\lambda,t)$. First,

$\displaystyle
B(t)=\sum_{k}\left(g_{k}^{*}e^{-i\omega_{k}t}b_{k}+g_{k}e^{i\omega_{k}t}b_{k}^{\dagger}\right),$
(30)

and,

$\displaystyle
B(\lambda)=\sum_{k}\left(g_{k}^{*}e^{-\lambda\omega_{k}}b_{k}+g_{k}e^{\lambda\omega_{k}}b_{k}^{\dagger}\right).$
(31)

Since $B_{\text{corr}}(\lambda,t)=\text{Tr}[\rho_{B}B(\lambda)B(t)]$, we find

$\displaystyle B_{\text{corr}}(\lambda,t)$
$\displaystyle=\sum_{k}|g_{k}|^{2}\Big{\\{}e^{-\omega_{k}\left(\lambda-
it\right)}$
$\displaystyle+2n_{k}\cosh{\left(\lambda\omega_{k}-i\omega_{k}t\right)}\Big{\\}},$
(32)

with $n_{k}$ given by Bose-Einstein statistics as

$\displaystyle n_{k}$
$\displaystyle=\frac{1}{2}\left\\{\coth{\left(\frac{\beta\omega_{k}}{2}\right)}-1\right\\}.$

To perform the sum over the environment modes, we use the spectral density
$J(\omega)$ via
$\sum_{k}|g_{k}|^{2}(\ldots)\rightarrow\int_{0}^{\infty}d\omega\,J(\omega)(\ldots)$.
We generally use an Ohmic spectral density of the form $J(\omega)=G\omega
e^{-\omega/\omega_{c}}$. The integrals are performed numerically to find
$J_{\text{corr}}(\beta,t)$, and the results are incorporated in the numerical
simulations of the master equation. We start by first looking at the pure
dephasing case ($\Delta=\Delta_{0}=0$), since this case can be solved exactly
and serves as a useful benchmark (see the appendix for details regarding the
exact solution). We illustrate our results in Fig. 1 for $N=1$. Two points
should be noted. First, the role played by initial correlations is very small
in this case. Second, our master equation reproduces the exact results very
well. Since we expect that the role of the initial correlations increases with
increasing $N$, we next look at $N=4$ and $N=10$. Results are shown in Figs. 2
and 3. It is clear that as $N$ increases, the initial correlations play a
larger and larger role. Moreover, the extra term in our master equation is
able to take into account the effect of the initial correlations very well.

Figure 1: Behavior of $j_{x}$ versus $t$ for N = 1 using the exact solution
with (blue circled dot) and without (purple squares) initial correlations, as
well as using the master equation with (solid, black line) and without
(dashed, red line)initial correlations. We have used
$\varepsilon=\varepsilon_{0}=4$, $G=0.05$, $\beta=1$ and $\omega_{c}=5$. Here
and in all other figures, the plotted variables are all in dimensionless units
Figure 2: Same as Fig. 1, except that we now have $N=4$. Figure 3: Same as
Fig. 2, except that we now have $N=10$.

Having shown that our master equation is able to reproduce results for the
pure dephasing model, we are now in a position to go beyond the pure dephasing
model and see the effects of the initial correlations. In Fig. 4, we have
shown the dynamics of $j_{x}$ with a non-zero value of the tunneling amplitude
for $N=2$. It is clear that the initial correlations do have a small influence
on the dynamics. This effect becomes more pronounced as we increase $N$ (see
Figs. 5 and 6). We have also looked at how the role played by the initial
correlations changes as the temperature changes. To this end, we compare Fig.
6, where the inverse temperature is $\beta=1$, with Fig. 7 where $\beta=0.5$
and Fig. 8 where $\beta=1.5$. As expected, at higher temperatures, the effect
of the initial correlations decreases, while at lower temperatures, the effect
of the initial correlations increases.

Figure 4: Behavior of $j_{x}$ against $t$ for $N=2$ with (black, solid) and
without (dashed, red) taking into account initial correlations. Here we have
used $\varepsilon_{0}=4$, $\varepsilon=2.5$ and $\Delta=\Delta_{0}=0.5$, while
the rest of the parameters used are the same as Fig. 1. Figure 5: Same as Fig.
4, except that we now have $N=4$. Figure 6: Same as Fig. 5, except that we now
have $N=10$. Figure 7: Behavior of $j_{x}$ against $t$ for $N=10$ with (black,
solid) and without (dashed, red) taking into account initial correlations.
Here we have used the same parameters as Fig. 6, except that $\beta=0.5$.
Figure 8: Same as Figs. 6 and 7, except that we now have $\beta=1.5$.

We now show that the effect of initial correlations are not manifested in the
dynamics of $j_{x}$ alone, we show in Figs. 9 and 10 the dynamics of
$j_{x}^{(2)}\equiv 4\expectationvalue{J_{x}^{2}}/N^{2}$. Such an observable is
relevant in the study of spin squeezing and entanglement. Once again, the
initial correlations can affect the dynamics significantly.

Figure 9: Behavior of $j_{x}^{(2)}$ against $t$ for $N=4$ with (black, solid)
and without (dashed, red) taking into account initial correlations. The rest
of the parameters used are the same as those in Fig. 4. Figure 10: Same as
Fig. 9, except that we now have $N=10$.

Finally, let us also demonstrate the effect of the initial correlations with
sub-Ohmic environment, that is,
$J(\omega)=G\omega^{s}\omega_{c}^{1-s}e^{-\omega/\omega_{c}}$. Since these
environments have longer correlation times, we expect that the effect of the
initial correlations will be greater as well. This is indeed the case, as can
be seen by comparing Figures 11 and 12 with Figs. 5 and Figs. 6 where an Ohmic
environment had been used.

Figure 11: Behavior of $j_{x}$ against $t$ for $N=4$ with (black, solid) and
without (dashed, red) taking into account initial correlations. Here we have
used a sub-Ohmic environment with $s=0.5$. We also have $\varepsilon_{0}=4$,
$\varepsilon=2.5$ and $\Delta=\Delta_{0}=0.5$, while the rest of the
parameters used are the same as Fig. 1. Figure 12: Same as Fig. 11, except
that now $N=10$.

##  IV Application to the spin-spin environment model

We consider now a collection of identical two-level systems interacting with
an environment consisting of two-level systems [38, 39, 40, 41, 42]. We then
have

$\displaystyle H_{\text{S0}}$
$\displaystyle=\varepsilon_{0}J_{z}+\Delta_{0}J_{x}$ $\displaystyle H_{S}$
$\displaystyle=\varepsilon J_{z}+\Delta J_{x}$ $\displaystyle H_{B}$
$\displaystyle=\sum_{k}\frac{\omega_{k}}{2}\sigma_{x}^{(k)},$ $\displaystyle
V$ $\displaystyle=J_{x}\otimes\sum_{k}g_{k}\sigma_{z}^{(k)}$

where $\sigma_{z}^{(k)}$ and $\sigma_{x}^{(k)}$ are the Pauli $z$-spin and
$x$-spin operators of the $k$th environment spin respectively, $\omega_{k}$
denotes the tunneling matrix element for the $k^{\text{th}}$ environment spin
and $g_{k}$ quantifies the coupling strength.
The different environment leads to a different correlation function $C_{ts}$
as well as a different factor $J^{R}_{\text{corr}}(\beta,t)$ that takes into
account the effect of the initial system-environment correlations. We first
note that

$\displaystyle B_{\text{corr}}(\lambda,t)$
$\displaystyle=\sum_{k}|g_{k}|^{2}\Big{\\{}\tanh{\left(\frac{\beta\omega_{k}}{2}\right)}e^{-\omega_{k}\left(\lambda-
it\right)}$
$\displaystyle+2n_{k}\sinh{\left(\lambda\omega_{k}-i\omega_{k}t\right)}\Big{\\}},$
(33)

with,

$\displaystyle
n_{k}=\frac{1}{2}\left\\{\coth{\left(\frac{\beta\omega_{k}}{2}\right)}-1\right\\}.$

Since the factors $\alpha_{1}(\lambda,t)$, $\alpha_{2}(\lambda,t)$, and
$\alpha_{3}(\lambda,t)$ are the same as before, this allows us to work out the
role of the initial correlations [see Eq. (29)]. In a similar manner, the
environment correlation function can also be worked out. As before, to perform
the sum over the environment modes, we use the same procedure as mentioned
earlier, that is, $\sum_{k}|g_{k}|^{2}(...)\rightarrow\int_{0}^{\infty}d\omega
J(\omega)(...)$. Results are shown in Figs. 13 and 14. Once again, the role of
the initial correlations is relatively small for a smaller value of $N$.
However, as $N$ increases, it is clear that we need to take into account the
role of the initial correlations to obtain an accurate picture of the system
dynamics.

Figure 13: Behavior of $j_{x}$ against $t$ for $N=4$ with (black, solid) and
without (dashed, red) taking into account initial correlations. Here we have
$\varepsilon_{0}=4$, $\varepsilon=2.5$, $\Delta=\Delta_{0}=0.5$, $G=0.05$,
$\beta=1$ and $\omega_{c}=5$. Figure 14: Same as Fig. 13, except that now
$N=10$.

##  V Conclusion

To conclude, we have shown that if we start from the joint thermal equilibrium
state of a quantum system and its environment, and then apply a unitary
operation to the system to prepare the system quantum state, the initial
correlations that exist in the joint thermal equilibrium state influence the
subsequent dynamics of the system. We have derived a time-local master
equation, correct to second-order in the system-environment coupling strength,
that takes into account the effect of these correlations, showing therefore
that one need not necessarily be in the strong system-environment coupling
regime to observe the effects of the initial correlations. The structure of
this master equation is very interesting, as the form of the term that takes
into account the initial correlations is the same as the relaxation and
dephasing term. In this sense, one can say that the initial correlations
affect the decoherence and dephasing rates, a fact which has already been
pointed out in studies of the role of initial correlations in pure dephasing
models [24]. Finally, we actually applied our master equation to the large
spin-boson model as well as to a collection of two-level systems interacting
with a spin environment to quantitatively investigate the role of the initial
correlations. We found that when the number of spins is small, then the
initial correlations do not play a significant role. However, for a larger
number of spins, the initial correlations must be accounted for in order to
explain the dynamics accurately.

## acknowledgements

A. R. M, M. Z. and A. Z. C. are grateful for support from HEC under grant No
5917/Punjab/NRPU/R&D/HEC/2016. Support from the LUMS Faculty Initiative Fund
is also acknowledged.

##  Appendix A The relaxation/dephasing term in the master equation

We look at the contribution of the third term in Eq. (II.2). We need to
consider only $\rho_{\text{tot}}^{(0)}(0)$ will contribute since we are
restricted to second-order terms in the master equation only. We obtain

$\displaystyle\alpha^{2}\text{Tr}_{\text{S,B}}\left\\{\rho_{\text{tot}}^{(0)}(0)\int_{0}^{t}\left[\left[\tilde{V}(t),\tilde{X}_{nm}(t)\right],\tilde{V}(s)\right]ds\right\\}$
$\displaystyle=\alpha^{2}\text{Tr}_{\text{S,B}}\left(\\{\tilde{\rho}_{S0}\otimes\rho_{B})\int_{0}^{t}\left[\left[\tilde{V}(t),\tilde{X}_{nm}(t)\right],\tilde{V}(s)\right]ds\right\\}$
$\displaystyle=\alpha^{2}\text{Tr}_{\text{S,B}}\Bigg{\\{}(\tilde{\rho}_{S0}\otimes\rho_{B})\Bigg{(}\int_{0}^{t}\tilde{V}(t)\tilde{X}_{nm}(t)\tilde{V}(s)ds$
$\displaystyle-\int_{0}^{t}\tilde{V}(s)\tilde{V}(t)\tilde{X}_{nm}(t)ds-\int_{0}^{t}\tilde{X}_{nm}(t)\tilde{V}(t)\tilde{V}(s)ds$
$\displaystyle+\int_{0}^{t}\tilde{V}(s)\tilde{X}_{nm}(t)\tilde{V}(t)ds\Bigg{)}\Bigg{\\}}.$
(34)

This splits up into four terms. The first term is

$\displaystyle\alpha^{2}\text{Tr}_{\text{S,B}}\left\\{(\tilde{\rho}_{S0}\otimes\rho_{B})\int_{0}^{t}\tilde{V}(t)\tilde{X}_{nm}(t)\tilde{V}(s)ds\right\\}$
$\displaystyle=\alpha^{2}\int_{0}^{t}\text{Tr}_{\text{S,B}}\Big{\\{}\tilde{\rho}_{S0}\otimes\rho_{B}U_{0}^{\dagger}(t)VU_{0}(t)U_{0}^{\dagger}(t)X_{nm}U_{0}(t)$
$\displaystyle\times U_{0}^{\dagger}(s)VU_{0}(s)\Big{\\}}\,ds$
$\displaystyle=\alpha^{2}\int_{0}^{t}\text{Tr}_{\text{S}}\left\\{\tilde{\rho}_{S0}U_{S}^{\dagger}(t)FY_{nm}U_{S}(t,s)FU_{S}(s)\right\\}$
$\displaystyle\times\text{Tr}_{\text{B}}\left\\{\rho_{B}B(t)B(s)\right\\}ds$
$\displaystyle=\alpha^{2}\int_{0}^{t}\left\langle
m\right|\bar{F}(t,s)\tilde{\rho}_{S0}(t)F\left|n\right\rangle C_{ts}ds.$

In a similar way, we can simplify the other terms. Putting them all back
together, and shifting to the basis-independent representation, we obtain the
third term in Eq. (23).

##  Appendix B Pure dephasing model

For completeness, we now sketch the derivation of the system dynamics for the
large spin-boson model. We have $N$ identical two-level systems interacting
with a common environment of harmonic oscillators. The dynamics of such a
system can be described by the Hamiltonian

$\displaystyle H_{\text{tot}}=H_{S}+H_{B}+V,$

where $H_{S}=\varepsilon J_{z}$,
$H_{B}=\sum_{k}\omega_{k}b_{k}^{\dagger}b_{k}$, and
$V=J_{z}\sum_{k}\left(g_{k}^{*}b_{k}+g_{k}b_{k}^{\dagger}\right)$. In the
interaction picture, the Hamiltonian becomes

$\displaystyle H_{I}\left(t\right)$
$\displaystyle=e^{i\left(H_{B}+H_{S}\right)t}Ve^{-i\left(H_{B}+H_{S}\right)t}$
$\displaystyle=J_{z}\sum_{k}\left(g_{k}^{*}b_{k}e^{-i\omega_{k}t}+g_{k}b_{k}^{\dagger}e^{i\omega_{k}t}\right).$

The corresponding unitary time-evolution operator can be found exactly via the
Magnus expansion. Consequently, the unitary operator for the complete
Hamiltonian is

$\displaystyle U\left(t\right)=e^{-iH_{S}t}e^{-iH_{B}t}U_{I}\left(t\right)$
$\displaystyle U\left(t\right)=e^{-i\varepsilon
J_{z}t}e^{-iH_{B}t}U_{I}\left(t\right),$

where,

$\displaystyle
U_{I}\left(t\right)=\text{exp}\left\\{J_{z}\sum_{k}\left[\alpha_{k}\left(t\right)b_{k}^{\dagger}-\alpha_{k}^{*}\left(t\right)b_{k}\right]-iJ_{z}^{2}t\Delta\left(t\right)\right\\}.$

Here
$\alpha_{k}\left(t\right)=\frac{g_{k}\left(1-e^{-i\omega_{k}t}\right)}{\omega_{k}}$
and
$\Delta\left(t\right)=\frac{1}{t}\sum_{k}\absolutevalue{g_{k}}^{2}\frac{\left[\sin(\omega_{k}t)-\omega_{k}t\right]}{\omega_{k}^{2}}$.
The reduced density operator of the system is then written in the $J_{z}$
eigenbasis as

$\displaystyle\left[\rho\left(t\right)\right]_{mn}=\text{Tr}_{S,B}\left[U\left(t\right)\rho\left(0\right)U^{\dagger}\left(t\right)P_{nm}\right]$
(35)

where $P_{nm}\equiv\left|n\right\rangle\left\langle m\right|,$ such that
$J_{z}\left|n\right\rangle=n\left|n\right\rangle.$ We can write Eq. (35) in
the Heisenberg picture, where
$P_{nm}\left(t\right)=U^{\dagger}\left(t\right)P_{nm}U\left(t\right)$ is the
Heisenberg picture operator. It follows that

$\displaystyle\left[\rho\left(t\right)\right]_{mn}=\text{Tr}_{S,B}\left[\rho\left(0\right)P_{nm}\left(t\right)\right].$

It is straightforward to find that

$\displaystyle
P_{nm}\left(t\right)=e^{-i\varepsilon\left(m-n\right)t}e^{-i\Delta\left(t\right)\left(m^{2}-n^{2}\right)t}e^{-R_{nm}\left(t\right)}P_{nm},$

where,

$\displaystyle
R_{nm}\left(t\right)=\left(n-m\right)\sum_{k}\left[\alpha_{k}\left(t\right)b_{k}^{\dagger}-\alpha_{k}^{*}\left(t\right)b_{k}\right].$

Hence, it follows that,

$\displaystyle\left[\rho\left(t\right)\right]_{mn}=$ $\displaystyle
e^{-i\varepsilon\left(m-n\right)t}e^{-i\Delta\left(t\right)\left(m^{2}-n^{2}\right)t}\times$
$\displaystyle\text{Tr}_{S,B}\left[\rho\left(0\right)e^{-R_{nm}\left(t\right)}P_{nm}\right]$
(36)

This is a general result because we have not yet defined the joint system-
environment initial state $\rho(0)$. Hence, it can be applied to both
uncorrelated and correlated initial states. We will derive an expression for
$\left[\rho_{S}\left(t\right)\right]_{mn}$ for both cases. First, without
considering the initial system-environment correlations,

$\displaystyle\rho_{\text{tot}}\left(0\right)=\rho\left(0\right)\otimes\rho_{B}$

where, $\rho_{B}=\frac{e^{-\beta H_{B}}}{Z_{B}}\ \text{with,}\
Z_{B}=\text{Tr}_{B}\left[e^{-\beta H_{B}}\right]$. We then get

$\displaystyle\left[\rho\left(t\right)\right]_{mn}=$
$\displaystyle\left[\rho\left(0\right)\right]_{mn}e^{-i\varepsilon\left(m-n\right)t}e^{-i\Delta\left(t\right)\left(m^{2}-n^{2}\right)t}\times$
$\displaystyle\text{Tr}_{B}\left[\rho_{B}e^{-R_{nm}\left(t\right)}\right].$
(37)

Simplifying the trace over the environment, we obtain

$\displaystyle\left[\rho\left(t\right)\right]_{mn}=$
$\displaystyle\left[\rho\left(0\right)\right]_{mn}e^{-i\varepsilon\left(m-n\right)t}e^{-i\Delta\left(t\right)\left(m^{2}-n^{2}\right)t}\times$
$\displaystyle e^{-\gamma\left(t\right)\left(m-n\right)^{2}t},$ (38)

with

$\displaystyle\gamma\left(t\right)=\frac{1}{t}\sum_{k}\frac{4\absolutevalue{g_{k}}^{2}\left[1-\cos\left(\omega_{k}t\right)\right]}{\omega_{k}^{2}}\coth{\left(\frac{\beta\omega_{k}}{2}\right)}.$
(39)

Next, we look at the case where the initial system-environment state is of the
form

$\displaystyle\rho\left(0\right)=\frac{\Omega e^{-\beta
H_{\text{tot}}}\Omega^{\dagger}}{Z},$ (40)

where $Z=\text{Tr}_{S,B}\left[\Omega e^{-\beta
H_{\text{tot}}}\Omega^{\dagger}\right],$ while $\Omega$ is a unitary operator.
We simplify $Z$ first by introducing a completeness relation,

$\displaystyle Z$ $\displaystyle=\sum_{l}\text{Tr}_{S,B}\left[\Omega e^{-\beta
H_{\text{tot}}}\left|l\right\rangle\left\langle
l\right|\Omega^{\dagger}\right],$ $\displaystyle=\sum_{l}e^{-\beta\varepsilon
l}\left\langle
l\right|\Omega^{\dagger}\Omega\left|l\right\rangle\text{Tr}_{B}\left[e^{-\beta
H_{B}^{\left(l\right)}}\right],$

where,

$\displaystyle
H_{B}^{\left(l\right)}=H_{B}+l\sum_{k}\left(g_{k}^{*}b_{k}+g_{k}b_{k}^{\dagger}\right).$

Using the displaced harmonic oscillator modes
$B_{k,l}=b_{k}+\frac{lg_{k}}{\omega_{k}}$ and
$B_{k,l}^{\dagger}=b_{k}^{\dagger}+\frac{lg_{k}^{*}}{\omega_{k}}$, we obtain

$\displaystyle Z=\sum_{l}e^{-\beta\varepsilon l}\left\langle
l\right|\Omega^{\dagger}\Omega\left|l\right\rangle e^{\beta
l^{2}\mathcal{C}}Z_{B},$

where, $\mathcal{C}=\sum_{k}\frac{4\absolutevalue{g_{k}}^{2}}{\omega_{k}}$.
Proceeding further, again using a completeness relation, we note that

$\displaystyle\left[\rho\left(t\right)\right]_{mn}$
$\displaystyle=\frac{1}{Z}\sum_{l}e^{-i\varepsilon\left(m-n\right)t}e^{-i\Delta\left(t\right)\left(m^{2}-n^{2}\right)t}\times$
$\displaystyle\left\langle
l\right|\Omega^{\dagger}P_{nm}\Omega\left|l\right\rangle\text{Tr}_{B}\left[{e^{-\beta
H_{B}^{\left(l\right)}}}e^{-R_{nm}\left(t\right)}\right],$

writing $R_{nm}\left(t\right)$ as

$\displaystyle
R_{nm}\left(t\right)=\left(n-m\right)\sum_{k}\left[\alpha_{k}\left(t\right)b_{k}^{\dagger}-\alpha_{k}^{*}\left(t\right)b_{k}\right]+i\Phi_{nm}^{\left(l\right)}\left(t\right),$

where,

$\displaystyle\Phi_{nm}^{\left(l\right)}\left(t\right)$
$\displaystyle=2l\left(n-m\right)\Phi\left(t\right),$
$\displaystyle\Phi\left(t\right)$
$\displaystyle=\sum_{k}\frac{4\absolutevalue{g_{k}}^{2}}{\omega_{k}^{2}}\sin\left(\omega_{k}t\right).$

It can be shown that

$\displaystyle\text{Tr}_{B}\left[{e^{-\beta
H_{B}^{\left(l\right)}}}e^{-R_{nm}\left(t\right)}\right]=e^{-i\Phi_{nm}^{\left(l\right)}\left(t\right)}e^{\beta
l^{2}\mathcal{C}}Z_{B}e^{-\gamma\left(t\right)\left(m-n\right)^{2}t}.$

Hence, we obtain

$\displaystyle\left[\rho\left(t\right)\right]_{mn}=\left[\rho\left(0\right)\right]_{mn}e^{-i\varepsilon\left(m-n\right)t}e^{-i\Delta\left(t\right)\left(m^{2}-n^{2}\right)t}\times$
$\displaystyle\,\,\,\,\,\,\,\,\,e^{-\gamma\left(t\right)\left(m-n\right)^{2}t}\frac{\sum_{l}\left[\left\langle
l\right|\Omega^{\dagger}P_{nm}\Omega\left|l\right\rangle
e^{-i\Phi_{nm}^{\left(l\right)}\left(t\right)}e^{-\beta\varepsilon l}e^{\beta
l^{2}\mathcal{C}}\right]}{\sum_{l}\left[\left\langle
l\right|\Omega^{\dagger}P_{nm}\Omega\left|l\right\rangle e^{-\beta\varepsilon
l}e^{\beta l^{2}\mathcal{C}}\right]}.$ (41)

## References

  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2007).
  * Weiss [2008] U. Weiss,  _Quantum dissipative systems_ (World Scientific, Singapore, 2008).
  * de Vega and Alonso [2017a] I. de Vega and D. Alonso, Dynamics of non-markovian open quantum systems, Rev. Mod. Phys. 89, 015001 (2017a).
  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Colloquium: Non-markovian dynamics in open quantum systems, Rev. Mod. Phys. 88, 021002 (2016).
  * Batalhão  _et al._ [2014] T. B. Batalhão, G. D. de Moraes Neto, M. A. de Ponte, and M. H. Y. Moussa, Nonperturbative approach to system-reservoir dynamics in the strong-coupling regime and non-markovian dynamics, Phys. Rev. A 90, 032105 (2014).
  * Hakim and Ambegaokar [1985] V. Hakim and V. Ambegaokar, Quantum theory of a free particle interacting with a linearly dissipative environment, Phys. Rev. A 32, 423 (1985).
  * Haake and Reibold [1985] F. Haake and R. Reibold, Strong damping and low-temperature anomalies for the harmonic oscillator, Phys. Rev. A 32, 2462 (1985).
  * Grabert  _et al._ [1988] H. Grabert, P. Schramm, and G.-L. Ingold, Quantum brownian motion: The functional integral approach, Phys. Rep. 168, 115 (1988).
  * Smith and Caldeira [1990] C. M. Smith and A. O. Caldeira, Application of the generalized feynman-vernon approach to a simple system: The damped harmonic oscillator, Phys. Rev. A 41, 3103 (1990).
  * Karrlein and Grabert [1997] R. Karrlein and H. Grabert, Exact time evolution and master equations for the damped harmonic oscillator, Phys. Rev. E 55, 153 (1997).
  * Dávila Romero and Pablo Paz [1997] L. Dávila Romero and J. Pablo Paz, Decoherence and initial correlations in quantum brownian motion, Phys. Rev. A 55, 4070 (1997).
  * Lutz [2003] E. Lutz, Effect of initial correlations on short-time decoherence, Phys. Rev. A 67, 022109 (2003).
  * Banerjee and Ghosh [2003] S. Banerjee and R. Ghosh, General quantum brownian motion with initially correlated and nonlinearly coupled environment, Phys. Rev. E 67, 056120 (2003).
  * van Kampen [2004] N. G. van Kampen, J. Stat. Phys. 115, 1057 (2004).
  * Ban [2009] M. Ban, Quantum master equation for dephasing of a two-level system with an initial correlation, Phys. Rev. A 80, 064103 (2009).
  * Campisi  _et al._ [2009] M. Campisi, P. Talkner, and P. Hänggi, Fluctuation theorem for arbitrary open quantum systems, Phys. Rev. Lett. 102, 210401 (2009).
  * Uchiyama and Aihara [2010] C. Uchiyama and M. Aihara, Role of initial quantum correlation in transient linear response, Phys. Rev. A 82, 044104 (2010).
  * Dijkstra and Tanimura [2010] A. G. Dijkstra and Y. Tanimura, Non-markovian entanglement dynamics in the presence of system-bath coherence, Phys. Rev. Lett. 104, 250401 (2010).
  * Smirne  _et al._ [2010] A. Smirne, H.-P. Breuer, J. Piilo, and B. Vacchini, Initial correlations in open-systems dynamics: The jaynes-cummings model, Phys. Rev. A 82, 062114 (2010).
  * Dajka and Łuczka [2010] J. Dajka and J. Łuczka, Distance growth of quantum states due to initial system-environment correlations, Phys. Rev. A 82, 012341 (2010).
  * Zhang  _et al._ [2010] Y.-J. Zhang, X.-B. Zou, Y.-J. Xia, and G.-C. Guo, Different entanglement dynamical behaviors due to initial system-environment correlations, Phys. Rev. A 82, 022108 (2010).
  * Tan and Zhang [2011] H.-T. Tan and W.-M. Zhang, Non-markovian dynamics of an open quantum system with initial system-reservoir correlations: A nanocavity coupled to a coupled-resonator optical waveguide, Phys. Rev. A 83, 032102 (2011).
  * Lee  _et al._ [2012] C. K. Lee, J. Cao, and J. Gong, Noncanonical statistics of a spin-boson model: Theory and exact monte carlo simulations, Phys. Rev. E 86, 021109 (2012).
  * Morozov  _et al._ [2012] V. G. Morozov, S. Mathey, and G. Röpke, Decoherence in an exactly solvable qubit model with initial qubit-environment correlations, Phys. Rev. A 85, 022101 (2012).
  * Semin  _et al._ [2012] V. Semin, I. Sinayskiy, and F. Petruccione, Initial correlation in a system of a spin coupled to a spin bath through an intermediate spin, Phys. Rev. A 86, 062114 (2012).
  * Chaudhry and Gong [2013a] A. Z. Chaudhry and J. Gong, Amplification and suppression of system-bath-correlation effects in an open many-body system, Phys. Rev. A 87, 012129 (2013a).
  * Chaudhry and Gong [2013b] A. Z. Chaudhry and J. Gong, Role of initial system-environment correlations: A master equation approach, Phys. Rev. A 88, 052107 (2013b).
  * Chaudhry and Gong [2013c] A. Z. Chaudhry and J. Gong, The effect of state preparation in a many-body system, Can. J. Chem. 92, 119 (2013c).
  * Reina  _et al._ [2014] J. Reina, C. Susa, and F. Fanchini, Extracting information from qubit-environment correlations, Sci. Rep. 4, 7443 (2014).
  * Zhang  _et al._ [2015] Y.-J. Zhang, W. Han, Y.-J. Xia, Y.-M. Yu, and H. Fan, Role of initial system-bath correlation on coherence trapping, Sci. Rep. 5, 13359 (2015).
  * Chen and Goan [2016] C.-C. Chen and H.-S. Goan, Effects of initial system-environment correlations on open-quantum-system dynamics and state preparation, Phys. Rev. A 93, 032113 (2016).
  * de Vega and Alonso [2017b] I. de Vega and D. Alonso, Dynamics of non-markovian open quantum systems, Rev. Mod. Phys. 89, 015001 (2017b).
  * Halimeh and de Vega [2017] J. C. Halimeh and I. de Vega, Weak-coupling master equation for arbitrary initial conditions, Phys. Rev. A 95, 052108 (2017).
  * Kitajima  _et al._ [2017] S. Kitajima, M. Ban, and F. Shibata, Expansion formulas for quantum master equations including initial correlation, J. Phys. A: Math. Theor 50, 125303 (2017).
  * Buser  _et al._ [2017] M. Buser, J. Cerrillo, G. Schaller, and J. Cao, Initial system-environment correlations via the transfer-tensor method, Phys. Rev. A 96, 062122 (2017).
  * Majeed and Chaudhry [2019] M. Majeed and A. Z. Chaudhry, Effect of initial systemâ€“environment correlations with spin environments, Eur. Phys. J. D 73, 16 (2019).
  * Bhaktavatsala Rao  _et al._ [2011] D. D. Bhaktavatsala Rao, N. Bar-Gill, and G. Kurizki, Generation of macroscopic superpositions of quantum states by linear coupling to a bath, Physical Review Letters 106, 010404 (2011).
  * Cucchietti  _et al._ [2005] F. M. Cucchietti, J. P. Paz, and W. H. Zurek, Decoherence from spin environments, Phys. Rev. A 72, 052113 (2005).
  * Camalet and Chitra [2007] S. Camalet and R. Chitra, Effect of random interactions in spin baths on decoherence, Phys. Rev. B 75, 094434 (2007).
  * Schlosshauer [2007] M. Schlosshauer,  _Decoherence and the quantum-to-classical transition_ (Springer, Berlin, 2007).
  * Villar [2009] P. I. Villar, Spin bath interaction effects on the geometric phase, Phys. Lett. A 373, 206 (2009).
  * Segal [2014] D. Segal, Two-level system in spin baths: Non-adiabatic dynamics and heat transport, The Journal of Chemical Physics 140, 164110 (2014).