# unyt: Handle, manipulate, and convert data with units in Python

Nathan J. Goldbaum  National Center for Supercomputing Applications,
University of Illinois at Urbana-Champaign  John A. ZuHone  Harvard-
Smithsonian Center for Astrophysics  Matthew J. Turk  National Center for
Supercomputing Applications, University of Illinois at Urbana-Champaign
Kacper Kowalik  National Center for Supercomputing Applications, University of
Illinois at Urbana-Champaign  Anna L. Rosen  Harvard-Smithsonian Center for
Astrophysics

††margin:  DOI: xxx Software • Review • Repository • Archive Submitted: 6 June
2018
Published: xxx Licence
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License (CC-BY).

## Summary

Software that processes real-world data or that models a physical system must
have some way of managing units. This might be as simple as the convention
that all floating point numbers are understood to be in the same physical unit
system (for example, the SI MKS units system). While simple approaches like
this do work in practice, they also are fraught with possible error, both by
programmers modifying the code who unintentionally misinterpret the units, and
by users of the software who must take care to supply data in the correct
units or who need to infer the units of data returned by the software.
Famously, NASA lost contact with the Mars Climate Orbiter spacecraft after it
crash-landed on the surface of Mars due to the use of English Imperial units
rather than metric units in the spacecraft control software (Board 1999).

The unyt library is designed both to aid quick calculations at an interactive
python prompt and to be tightly integrated into a larger Python application or
library. The top-level unyt namespace ships with a large number of predefined
units and physical constants to aid setting up quick calculations without
needing to look up unit data or the value of a physical constant. Using the
unyt library as an interactive calculation aid only requires knowledge of
basic Python syntax and awareness of a few of the methods of the unyt_array
class - for example, the unyt_array.to() method to convert data to a different
unit. As the complexity of the usage increases, unyt provides a number of
optional features to aid these cases, including custom unit registries
containing both predefined physical units as well as user-defined units,
built-in output to disk via the pickle protocol and to HDF5 files using the
h5py library (Collette 2013), and round-trip converters for unit objects
defined by other popular Python unit libraries.

Physical units in the unyt class are defined in terms of the dimensions of the
unit, a string representation, and a floating point scaling to the MKS unit
system. Rather than implementing algebra for unit expressions, we rely on the
SymPy symbolic algebra library (Meurer et al. 2017) to handle symbolic
algebraic manipulation. The unyt.Unit object can represent arbitrary units
formed out of base dimensions in the SI unit system: time, length, mass,
temperature, luminance, and electric current. We currently treat units such as
mol with the seventh SI base dimension, amount of substance, as dimensionless,
although we are open to changing this based on feedback from users. In
addition, unyt supports forming quantities defined in other unit systems - in
particular CGS Gaussian units common in astrophysics as well as geometrized
“natural” units common in relativistic calculations. In addition, unyt ships
with a number of other useful predefined unit systems based, including
imperial units, Planck units, a unit system for calculations in the solar
system, and a “galactic” unit system based on the solar mass, kiloparsecs, and
Myr, a convention common in galactic astronomy.

In addition to the unyt.Unit class, unyt also provides a two subclasses of the
NumPy (Oliphant 2006) ndarray (Walt, Colbert, and Varoquaux 2011),
unyt.unyt_array and unyt.unyt_quantity to represent arrays and scalars with
units attached, respectively. The unyt library also provides a
unyt.UnitRegistry class to allow custom systems of units, for example to track
the internal unit system used in a simulation. These subclasses are tightly
integrated with the NumPy ufunc system, which ensures that algebraic
calculations that include data with units automatically check to make sure the
units are consistent, and allow automatic converting of the final answer of a
calculation into a convenient unit.

We direct readers interested in usage examples and a guide for integrating
unyt into an existing Python application or workflow to the unyt documentation
hosted at http://unyt.readthedocs.io/en/latest/.

## Comparison with Pint and astropy.units

The scientific Python ecosystem has a long history of efforts to develop a
library to handle unit conversions and enforce unit consistency. For a
relatively recent review of these efforts, see (Bekolay 2013). While we won’t
exhaustively cover extant Python libraries for handling units in this paper,
we will focus on Pint (Grecco 2018) and astropy.units (The Astropy
Collaboration et al. 2018), which both provide a robust implementation of an
array container with units attached and are commonly used in research software
projects. At time of writing a GitHub search for import astropy.units returns
approximately 10,500 results and a search for import pint returns
approximately 1,500 results.

While unyt provides functionality that overlaps with astropy.units and Pint,
there are important differences which we elaborate on below. In addition, it
is worth noting that all three codebases had origins at roughly the same time
period. Pint initially began development in 2012 according to the git
repository logs. Likewise astropy.units began development in 2012 and was
released as part of astropy 0.2 in 2013, although the initial implementation
was adapted from the pynbody library (Pontzen et al. 2013), which started its
units implementation in 2010 according to the git repository logs. In the case
of unyt, it originated via the dimensionful library (Stark 2012) in 2012.
Later, dimensionful was elaborated on and improved to become yt.units, the
unit system for the yt library (Turk et al. 2011) at a yt developer workshop
in 2013 and was subsequently released as part of yt 3.0 in 2014. One of the
design goals for the yt unit system was the ability to dynamically define
“code” units (e.g. units internal to data loaded by yt) as well as units that
depend on details of the dataset - in particular cosmological comoving units
and the “little $h$” factor, used to parameterize the Hubble constant in
cosmology calculations (Croton 2013). For cosmology simulations in particular,
comparing data with different unit systems can be tricky because one might
want to use data from multiple outputs in a time series, with each output
having a different mapping from internal units to physical units. This despite
the fact that each output in the time series represents the same physical
system and common workflows involve combining data from multiple outputs. This
requirement to manage complex custom units and interoperate between custom
unit systems drove the yt community to independently develop a custom unit
system solution. We have decided to repackage and improve yt.units in the form
of unyt to both make it easier to work on and improve the unit system and
encourage use of the unit system for scientific python users who do not want
to install a heavy-weight dependency like yt.

Below we present a table comparing unyt with astropy.units and Pint. Estimates
for lines of code in the library were generated using the cloc tool (Danial
2018); blank and comment lines are excluded from the estimate. Test coverage
was estimated using the coveralls output for Pint and astropy.units and using
the codecov.io output for unyt.

Library | unyt | astropy.units | Pint
---|---|---|---
Lines of code | 5128 | 10163 | 8908
Lines of code excluding tests | 3195 | 5504 | 4499
Test Coverage | 99.91% | 93.63% | 77.44%

We offer lines of code as a very rough estimate for the “hackability” of the
codebase. In general, smaller codebases with higher test coverage have fewer
defects (Lipow 1982; Koru, Zhang, and Liu 2007; Gopinath, Jensen, and Groce
2014). This comparison is somewhat unfair in favor of unyt in that
astropy.units only depends on NumPy and Pint has no dependencies, while unyt
depends on both SymPy and NumPy. Much of the reduction in the size of the unyt
library can be attributed to offloading the handling of algebra to SymPy
rather than needing to implement the algebra of unit symbols directly in unyt.
For potential users who are wary of adding SymPy as a dependency, that might
argue in favor of using Pint in favor of unyt.

### astropy.units

The astropy.units subpackage provides a PrefixUnit class, a Quantity class
that represents both scalar and array data with attached units, and a large
number of predefined unit symbols. The preferred way to create Quantity
instances is via multiplication with a PrefixUnit instance. Similar to unyt,
the Quantity class is implemented via a subclass of the NumPy ndarray class.
Indeed, in many ways the everyday usage patterns of astropy.units and unyt are
similar, although unyt is not quite a drop-in replacement for astropy.units as
there are some API differences. The main functional difference between
astropy.units and unyt is that astropy.units is a subpackage of the larger
astropy package. This means that depending on astropy.units requires
installing a large collection of astronomically focused software included in
the astropy package, including a substantial amount of compiled C code. This
presents a barrier to usage for potential users of astropy.units who are not
astronomers or do not need the observational astronomy capabilities provided
by astropy.

### Pint

The Pint package provides a different API for accessing units compared with
unyt and astropy.units. Rather than making units immediately importable from
the Pint namespace, Pint instead requires users to instantiate a UnitRegistry
instance (unrelated to the unyt.UnitRegistry class), which in turn has Unit
instances as attributes. Just like with unyt and astropy.units, creating a
Quantity instance requires multiplying an array or scalar by a Unit instance.
Exposing the UnitRegistry directly to all users like this does force users of
the library to think about which system of units they are working with, which
may be beneficial in some cases, however it also means that users have a bit
of extra cognitive overhead they need to deal with every time they use Pint.

Figure 1: A benchmark comparing the ratio of the time to apply units to lists
and NumPy ndarray instances to the time to interpret the same list or ndarray
to an ndarray. This ratio, $T_{\rm package}/T_{\rm numpy}$, corresponds to the
overhead of converting data to work with one of the three packages. Values
close to unity correspond to zero or negligible overhead, while values larger
than unity correspond to measureable overhead. Optimally all values would be
near zero. In practice, applying units to small arrays incurs substantial
overhead. Each test is shown for three different sizes of input data,
including inputs with size 3, 1,000, and 1,000,000. The black lines at the top
of the bars indicate the sample standard deviation. The $T_{\rm numpy}$ time
is calculated by benchmarking the time to perform np.asarray(data) where data
is either a list or an ndarray.

In addition, the Quantity class provided by Pint is not a subclass of NumPy’s
ndarray. Instead, it is a wrapper around an internal ndarray buffer. This
simplifies the implementation of Pint by avoiding the somewhat arcane process
for creating an ndarray subclass, although the Pint Quantity class must also
be careful to emulate the full NumPy ndarray API so that it can be a drop-in
replacement for ndarray.

Finally, in comparing the output of our benchmarks of Pint, astropy.units, and
unyt, we found that in-place operations making use of a NumPy ufunc will
unexpectedly strip units in Pint. For example, if a and b are Pint Quantity
instances, np.add(a, b, out=out)) will operate on a and b as if neither have
units attached. Interestingly, without the out keyword, Pint does get the
correct answer, so it is possible that this is a bug in Pint, and we have
reported it as such upstream (see https://github.com/hgrecco/pint/issues/644).

### Performance Comparison

Checking units will always add some overhead over using hard-coded unit
conversion factors. Thus a library that is entrusted with checking units in an
application should incur the minimum possible overhead to avoid triggering
performance regressions after integrating unit checking into an application.
Optimally, a unit library will add zero overhead regardless of the size of the
array. In practice that is not the case for any of the three libraries under
consideration, and there is a minimum array size above which the overhead of
doing a mathematical operation exceeds the overhead of checking units. It is
thus worth benchmarking unit libraries in a fair manner, comparing with the
same operation implemented using plain NumPy.

Figure 2: A benchmark comparing the time to square an array and to take the
square root of an array. See Figure 1 for a detailed explanation of the plot
style.

Here we present such a benchmark. We made use of the perf (Stinner 2018)
Python benchmarking tool, which not only provides facilities for establishing
the statistical significance of a benchmark run, but also can tune a linux
system to turn off operating system and hardware features like CPU throttling
that might introduce variance in a benchmark. We made use of a Dell Latitude
E7270 laptop equipped with an Intel i5-6300U CPU clocked at 2.4 Ghz. The
testing environment was based on Python 3.6.3 and had NumPy 1.14.2, sympy
1.1.1, fastcache 1.0.2, Astropy 3.0.1, and Pint 0.8.1 installed. fastcache
(Brady 2017) is an optional dependency of SymPy that provides an optimized LRU
cache implemented in C that can substantially speed up SymPy. The system was
instrumented using perf system tune to turn off CPU features that might
interfere with stable benchmarks. We did not make any boot-time Linux kernel
parameter changes.

Figure 3: A benchmark comparing the time to perform various binary arithmetic
operations on input operans that have different but dimensionallty compatible
units. See Figure 1 for a detailed explanation of the plot style.

For each of the benchmarks below, we show the ratio of the time to perform an
operation with one of unyt, Pint, and astopy.units, $T_{\rm package}$, to the
time it takes for NumPy to perform the equivalent operation, $T_{\rm numpy}$.
For example, for the comparison of the performance of np.add(a, b) where a and
b have different units with the same dimension, the corresponding benchmark to
generate $T_{\rm numpy}$ would use the code np.add(a, c*b) where a and b would
be ndarray instances and c would be the floating point conversion factor
between the units of a and b. Much of the time in $T_{\rm package}$ relative
to $T_{\rm numpy}$ is spent in the respective packages calculating the
appropriate conversion factor c. Thus the comparisons below depict very
directly the overhead for using a unit library over an equivalent operation
that uses hard-coded unit-conversion factors.

Figure 4: A benchmark comparing the time to perform various binary arithmetic
operations on input operands that have the same units. See Figure 1 for a
detailed explanation of the plot style.

#### Applying units to data

In Figure 1 we plot the overhead for applying units to data, showing both
Python lists and NumPy ndarray instances as the input to apply data to. Since
all three libraries eventually convert input data to a NumPy ndarray, the
comparison with array inputs more explicitly shows the overhead for _just_
applying units to data. When applying units to a list, all three libraries as
well as NumPy need to first copy the contents of the list into a NumPy array
or a subclass of ndarray. This explains why the overhead is systematically
lower when starting with a list.

Figure 5: A benchmark comparing the overhead for computing various NumPy ufunc
operations. The operands of all binary ufuncs have the same units. See Figure
1 for a detailed explanation of the plot style.

In all cases, unyt either is fastest by a statistically significant margin, or
ties with astropy. Even for large input arrays, Pint still has statistically
significant overhead, while both unyt and astropy.units have negligible
overhead once the input array size reaches $10^{6}$ elements.

#### Unary arithmetic operations

Expressions involving powers of data with units, including integer and
fractional powers, are very common in the physical sciences. It is therefore
very important for a library that handles units to be able to track this case
in a performant way. In Figure 2 we present a benchmark comparing Pint, unyt,
and astropy.units for the squaring and square root operation. In all cases,
unyt has the lowest overhead, with Pint coming in second, and astropy.units
trailing. Note that the y-axis is plotted on a log scale, so astropy is as
much as 4 times slower than unyt for these operations.

#### Binary arithmetic operations

Binary operations form the core of arithmetic. It is vital for a library that
handles unit manipulation to both transparently convert units when necessary
and to ensure that expressions involving quantities with units are
dimensionally consistent. In Figure 3 and 4 we present benchmarks for binary
arithmetic expressions, both with input data that has the same units and with
input data with different units but the same dimensions. In most cases, unyt
has less overhead than both astropy and Pint, although there are a few
anomalies that are worth explaining in more detail. For comparison operations,
Pint exhibits a slowdown even on large input arrays. This is not present for
other binary operations, so it is possible that this overhead might be
eliminated with a code change in Pint. For multiplication on large arrays, all
three libraries have measured overhead of $\sim 0.5$ than that of the
“equivalent” Numpy operation. That is because all three libraries produce
results with units given by the product of the input units. That is, there is
no need to multiply the result of the multiplication operation by an
additional constant, while a pure NumPy implementation would need to multiply
the result by a constant to ensure both operands of the operation are in the
same units. Finally, for division, both Pint and astropy.units exhibit the
same behavior as for multiplication, and for similar reasons: the result of
the division operation is output with units given by the ratio of the input
units. On the other hand, unyt will automatically cancel the dimensionally
compatible units in the ratio and return a result with dimensionless units.

Figure 6: The same as Figure 5, but with in-place ufunc operations. See Figure
1 for a detailed explanation of the plot style.

#### NumPy ufunc performance

Lastly, In Figures 4 and 5, we present benchmarks of NumPy ufunc operations. A
NumPy ufunc is a fast C implementation of a basic mathematical operation. This
includes arithmetic operators as well as trigonometric and special functions.
By using a ufunc directly, one bypasses the Python object protocol and short-
circuits directly to the low-level NumPy math kernels. We show both directly
using the NumPy ufunc operators (Figure 4) as well as using the same operators
with a pre-allocated output array to benchmark in-place operations.

As for the other benchmarks, unyt tends to have the lowest amount of overhead,
although there are some significant exceptions. For np.power, Pint has the
lowest overhead, except for very large input arrays, where the overhead for
all three libraries is negligible. On the other hand, for np.sqrt, np.equal,
np.add, and np.subtract, Pint still has statistically significant overhead for
large input arrays. Finally, for the in-place ufunc comparison, Pint has the
lowest overhead for all operations. However, as discussed above, this is
because of a bug in Pint which causes the library to ignore units when calling
a ufunc with the out keyword set.

## Conclusions

In this paper we present the unyt library, giving background on the reasons
for its existence and some historical context for its origin. We also present
a set of benchmarks for common arithmetic operations, comparing the
performance of unyt with Pint and astropy.units. In general, we find that unyt
either outperforms or matches the performance astropy.units and Pint,
depending on the operation and size of the input data. We also demonstrate
that the unyt library constitutes a smaller codebase with higher test coverage
than both Pint and astropy.units.

## Acknowledgements

NJG would like to thank Brandon Carswell and Alex Farthing of the NCSA IT
staff for providing a laptop with Linux installed for the performance
benchmark. This work was supported by NSF grant OAC-1663914 (NJG, MJT), by the
Gordon and Betty Moore Foundation’s Data-Driven Discovery Initiative through
Grant GBMF4561 (MJT) and by NASA through Einstein Postdoctoral Fellowship
grant number PF7-180166 awarded by the Chandra X-ray Center, which is operated
by the Smithsonian Astrophysical Observatory for NASA under contract
NAS8-03060 (ALR).

## References

Bekolay, Trevor. 2013. “A Comprehensive Look at Representing Physical
Quantities in Python.” https://www.youtube.com/watch?v=N-edLdxiM40.

Board, Mars Climate Orbiter Mishap Investigation. 1999. _Mars Climate Orbiter
Mishap Investigation Board: Phase I Report_. Jet Propulsion Laboratory.
https://books.google.com/books?id=4OMIHQAACAAJ.

Brady, Peter. 2017. “Fastcache.” _GitHub Repository_.
https://github.com/pbrady/fastcache; GitHub.

Collette, Andrew. 2013. _Python and Hdf5_. O’Reilly.

Croton, D. J. 2013. “Damn You, Little h! (Or, Real-World Applications of the
Hubble Constant Using Observed and Simulated Data)” 30 (October): e052.
https://doi.org/10.1017/pasa.2013.31.

Danial, Al. 2018. “Cloc.” _GitHub Repository_.
https://github.com/AlDanial/cloc; GitHub.

Gopinath, Rahul, Carlos Jensen, and Alex Groce. 2014. “Code Coverage for Suite
Evaluation by Developers.” In _Proceedings of the 36th International
Conference on Software Engineering_ , 72–82. ICSE 2014. New York, NY, USA:
ACM. https://doi.org/10.1145/2568225.2568278.

Grecco, Hernan E. 2018. “Pint.” _GitHub Repository_.
https://github.com/hgrecco/pint; GitHub.

Koru, A. G., D. Zhang, and H. Liu. 2007. “Modeling the Effect of Size on
Defect Proneness for Open-Source Software.” In _Predictor Models in Software
Engineering, 2007. PROMISE’07: ICSE Workshops 2007. International Workshop on_
, 10–10. https://doi.org/10.1109/PROMISE.2007.9.

Lipow, M. 1982. “Number of Faults Per Line of Code.” _IEEE Trans. Softw. Eng._
8 (4). Piscataway, NJ, USA: IEEE Press: 437–39.
https://doi.org/10.1109/TSE.1982.235579.

Meurer, Aaron, Christopher P. Smith, Mateusz Paprocki, Ondřej Čertík, Sergey
B. Kirpichev, Matthew Rocklin, AMiT Kumar, et al. 2017. “SymPy: Symbolic
Computing in Python.” _PeerJ Computer Science_ 3 (January): e103.
https://doi.org/10.7717/peerj-cs.103.

Oliphant, Travis E. 2006. _A Guide to Numpy_. Trelgol Publishing.

Pontzen, A., R. Roškar, G. S. Stinson, R. Woods, D. M. Reed, J. Coles, and T.
R. Quinn. 2013. “pynbody: Astrophysics Simulation Analysis for Python.”

Stark, Casey W. 2012. “Dimensionful.” _GitHub Repository_.
https://github.com/caseywstark/dimensionful; GitHub.

Stinner, Victor. 2018. “Perf.” _GitHub Repository_.
https://github.com/vstinner/perf; GitHub.

The Astropy Collaboration, A. M. Price-Whelan, B. M. Sipőcz, H. M. Günther, P.
L. Lim, S. M. Crawford, S. Conseil, et al. 2018. “The Astropy Project:
Building an inclusive, open-science project and status of the v2.0 core
package.” _ArXiv E-Prints_ , January.

Turk, M. J., B. D. Smith, J. S. Oishi, S. Skory, S. W. Skillman, T. Abel, and
M. L. Norman. 2011. “yt: A Multi-code Analysis Toolkit for Astrophysical
Simulation Data.” _The Astrophysical Journal Supplement Series_ 192 (January):
9. https://doi.org/10.1088/0067-0049/192/1/9.

Walt, S. van der, S. C. Colbert, and G. Varoquaux. 2011. “The Numpy Array: A
Structure for Efficient Numerical Computation.” _Computing in Science
Engineering_ 13 (2): 22–30. https://doi.org/10.1109/MCSE.2011.37.


# Using the Environment to Understand non-Markovian Open Quantum Systems

Dominic Gribben  SUPA, School of Physics and Astronomy, University of St
Andrews, St Andrews KY16 9SS, United Kingdom  Aidan Strathearn  School of
Mathematics and Physics, The University of Queensland, St Lucia, Queensland
4072, Australia  Gerald E. Fux  SUPA, School of Physics and Astronomy,
University of St Andrews, St Andrews KY16 9SS, United Kingdom  Peter Kirton
Department of Physics and SUPA, University of Strathclyde, Glasgow G4 0NG,
United Kingdom  Brendon W. Lovett  SUPA, School of Physics and Astronomy,
University of St Andrews, St Andrews KY16 9SS, United Kingdom

###### Abstract

Tracing out the environmental degrees of freedom is a necessary procedure when
simulating open quantum systems. While being an essential step in deriving a
tractable master equation it represents a loss of information. In situations
where there is strong interplay between the system and environmental degrees
of freedom this loss makes understanding the system’s dynamics challenging.
These dynamics, when viewed in isolation, are non-Markovian and memory effects
induce complex features that are difficult to interpret. Here we exploit a
numerically exact approach to simulating the system dynamics of the spin-boson
model, which is based on the construction and contraction of tensor network
that represents the process tensor of this open quantum system. We are then
able to find any system correlation function exactly. We show that we can use
these to infer any correlation function of a Gaussian environment, so long as
the coupling between system and envinroment is linear. This not only allows
reconstruction of the full dynamics of both system and environment, but also
opens avenues into studying the effect of a system on its environment. To
demonstrate the applicability of our method we show how heat moves between
different modes of a bosonic bath when coupled to a two-level system that is
subject to an off-resonant drive.

##  1 Introduction

When modelling a realistic quantum system the effect of its environment must
be captured in some form [1]. The environment consists of many degrees of
freedom and is often approximated as an infinite collection of harmonic
oscillators. Thus, an explicit description of the environment is usually
impossible. Most methods for describing open quantum systems rely on tracing
out the environment to obtain an effective description of the system involving
only a tractable number of degrees of freedom. When the system-environment
coupling is strong, and/or when the spectral density of the environment is
structured, the full dynamics involves a significant interplay between the
system and environment degrees of freedom [2, 3, 4, 5, 6, 7, 8, 9, 10]. When
the environment is traced out in this case, the dynamics of the system are
non-Markovian [11, 12], whose simulation generally requires sophisticated
numerical methods. In addition, such system dynamics are often complex to
interpret.

Analyzing the behaviour of the environment is useful for providing insights
into processes where the system-environment interaction is non-trivial or for
any application where the goal is to manipulate the environment via the
system. An important example of this class of problem occurs in modelling
quantum thermal machines [13, 14] where distinct thermodynamic effects beyond
weak coupling can be observed [15, 16]. Access to the state of the environment
could also be used to explicitly track the formation of bound states such as
polarons [17, 18, 19, 20] or polaritons [17, 21]. There is also potential for
simulating phase transitions in the environment, such as to lasing or
superradiant states [22] in multimode cavity QED [23]. As well as these, the
technique we introduce here allows for calculation of out-of-time-order
correlations of the environment; these are known to characterise the degree of
information scrambling in many-body quantum systems [24, 25, 26].

Some existing techniques allow varying levels of insight into the behaviour of
the environment. The reaction coordinate mapping explicitly calculates the
dynamics of a collective coordinate that gives a sense of the behaviour of the
environment [27, 28]. It has also been shown that the auxiliary density
operators calculated in the hierarchy equations of motion method can be used
to calculate the same collective coordinate dynamics [29] as well as higher
order moments of the total bath coupling operator [30, 31, 32]. Counting field
techniques can also be used to calculate moments of the total bath coupling
operator [33, 34, 35]. Techniques such as TEDOPA [36] map the star topology of
an environment to a chain of modes with nearest neighbour couplings. Dynamics
within this chain can also be used to characterise the behaviour of the
environment [37] and these dynamics can be passed through an inverse chain-
mapping to extract the full environmental dynamics [38, 39]. Here, we give a
prescription that allows calculation of any multi-time bath correlation
function using only system moments of equal order: no auxiliary modes or
additional degrees of freedom are required. This allows their calculation
using even those methods where the entire environment is traced out.

Calculating the required correlation functions for a system undergoing non-
Markovian evolution is numerically challenging for all but the simplest cases.
It is necessary to use a technique which is capable of recording the history
of the system dynamics while allowing for the application of arbitrary
sequences of system operators. Here we employ a version of the time-evolving
matrix product operator (TEMPO) algorithm [40]. This has been shown to
efficiently capture the required time non-local influences by representing the
system’s discretised trajectories as a matrix product state (MPS). This is
very similar to the use of MPS in simulating 1D quantum many-body problems
which can be compressed by truncating the sizes of the matrices, keeping only
the most relevant contributions. However, even with the TEMPO algorithm,
calculating all of the necessary correlation functions would require many
separate calculations. This difficulty can be overcome by building on the
approach developed in [41] where a modification of the network contracted in
the TEMPO algorithm allowed efficient generation of an object called the
process tensor. From a single process tensor it is possible to access not just
the system dynamics but also all system correlation functions at arbitrary
times [42]. This significantly reduces the computational overhead of
calculating the behaviour of the environment.

In a previous work [43] we showed how the displacement of a bosonic bath that
is linearly coupled to a spin system can be calculated from a simple
transformation of the system dynamics. We used this to gain insight into the
origin of some of the non-Markovian features in the dynamics. The approach
used there is limited to the calculation of first moments. In this paper we
reformulate the problem to show how we can calculate all higher order bath
correlation functions. In Sec. 2 we begin by extending the existing formalism
to show how any bath correlation function can be calculated in terms of equal
order system moments. An overview of the process tensor formulation [44, 42]
of open quantum systems is given in Sec. 3 where we demonstrate how it can be
used to efficiently calculate the system moments necessary for calculation of
higher order bath correlation functions. Finally, in Sec. 4 we use this
approach to study how the biased spin-boson model [45] can be tuned to
redistribute energy in its environment.

##  2 Bath Dynamics

Here we provide a method to calculate any correlation function of bath
operators in terms of those of system operators, for a system linearly coupled
to a Gaussian environment. We do this by expressing the relevant expectation
values of system operators as derivatives of a generating functional. By
adding source terms to the propagator that remain easy to integrate out for
the Gaussian bath, the derivatives required to calculate the relevant moments
can then be carried out explicitly.

The system-environment Hamiltonian we consider has the form:

$\displaystyle H$ $\displaystyle=H_{S}+H_{I}+H_{B}$
$\displaystyle=H_{S}+s\sum_{q}g_{q}(a_{q}+a_{q}^{\dagger})+\sum_{q}\omega_{q}a_{q}^{\dagger}a_{q},$
(1)

where $a_{q}^{(\dagger)}$ destroys (creates) an excitation in bosonic
environment mode $q$ with frequency $\omega_{q}$. Here, $s$ is a general
operator acting exclusively on the system Hilbert space and $g_{q}$ is the
coupling strength of mode $q$ to the system. Formally, the full dynamics, in
the interaction picture, can be written as [46]:

$\rho(t)=\overleftarrow{T}\exp\left[\int_{0}^{t}\mathcal{L}_{I}(t^{\prime})\,dt^{\prime}\right]\rho(0).$
(2)

The full density operator of the system and environment, is $\rho(t)$ and
$\mathcal{L}_{I}$ is the Liouvillian superoperator corresponding to
commutation with the interaction Hamiltonian and $\overleftarrow{T}$ indicates
that superoperators are to be time-ordered from right to left. Details of the
superoperator formalism are given in App. A.

To proceed we assume that the initial state is factorizable and can be written
as $\rho(0)=\rho_{S}(0)\otimes\rho_{B}(0)$. We can now trace out the bath and
arrive at an expression for the evolution of the system:

$\rho_{S}(t)=\left\langle\overleftarrow{T}\exp\left[\int_{0}^{t}\mathcal{L}_{I}(t^{\prime})\,dt^{\prime}\right]\right\rangle_{B}\rho_{S}(0),$
(3)

where $\left\langle\cdot\right\rangle_{B}$ denotes an expectation value with
respect to the initial state of the bath.

If we further assume that the initial state of the bath is Gaussian we can
apply the relation $\langle\text{e}^{X}\rangle=\text{e}^{\langle
X^{2}\rangle/2+\langle X\rangle}$ which is valid for expectations taken over
any $X$ that has a Gaussian distribution. This condition holds for
$\mathcal{L}_{I}$, which is linear in the $a_{q}^{(\dagger)}$, when the
expectation is taken with respect to $\rho_{B}$. We define the Hamiltonians
such that $\langle X\rangle=0$ and so the environment part of Eq. (3) can be
simplified to:

$\left\langle\overleftarrow{T}\exp\left[\int_{0}^{t}\mathcal{L}_{I}(t^{\prime})\,dt^{\prime}\right]\right\rangle_{B}=\\\
\exp\left[\frac{1}{2}\int_{0}^{t}\int_{0}^{t}\left\langle\overleftarrow{T}\mathcal{L}_{I}(t^{\prime})\mathcal{L}_{I}(t^{\prime\prime})\right\rangle_{B}\,dt^{\prime\prime}\,dt^{\prime}\right].$
(4)

The right-hand side of the above equation is exactly the Feynman-Vernon
influence functional [47] in superoperator form. Given that this relation
holds for any $\mathcal{L}_{I}(t)$ that is linear in $a_{q}$ and
$a_{q}^{\dagger}$, we can take
$\mathcal{L}_{I}(t)\rightarrow\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t)$
with the addition of source terms:

$\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t)=\mathcal{L}_{I}(t)+\sum_{q}\sum_{\alpha=L,R}\Big{[}\Lambda_{q}^{\alpha}(t)\mathcal{A}_{q}^{\alpha}(t)+\Lambda^{\alpha*}_{q}(t)\mathcal{A}_{q}^{\alpha\dagger}(t)\Big{]}+\phi(\Lambda,\Lambda^{*},t),$
(5)

where $\Lambda_{q}^{\alpha(*)}(t)$ are scalar valued functions associated with
superoperators $\mathcal{A}_{q}^{\alpha(\dagger)}(t^{\prime})$. These
superoperators can be expressed in terms of the action of $a^{(\dagger)}$ on
$\rho$:

$\displaystyle\mathcal{A}_{q}^{L(\dagger)}\rho=a_{q}^{(\dagger)}\rho,$
$\displaystyle\mathcal{A}_{q}^{R(\dagger)}\rho=\rho a_{q}^{(\dagger)}.$ (6)

We have included a scalar valued function $\phi$ in the definition of
$\overline{\mathcal{L}_{I}}$; the exact form of $\phi$ is dependent on the
order of the ladder operators in the quantity of interest. Using Eq. (5) we
can define a generating functional

$G(\Lambda,\Lambda^{*})=\left\langle\overleftarrow{T}\exp\left[\int_{0}^{t}\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t^{\prime})\,dt^{\prime}\right]\right\rangle,$
(7)

where the angular brackets now denote the expectation be taken with respect to
the total initial density matrix $\rho$. The expectation of any function of
ladder operators $f$ can be expressed as functional derivatives of Eq. (7):

$\operatorname{Tr}[f(\\{a_{q}(t_{m}),a_{q}^{\dagger}(t_{n})\\})\rho]=f\left(\left\\{\frac{\delta}{\delta\Lambda_{q}^{L}(t_{m})},\frac{\delta}{\delta\Lambda_{q}^{L*}(t_{n})}\right\\}\right)G(\Lambda,\Lambda^{*})\Bigg{\rvert}_{\Lambda,\Lambda^{*}=0}.$
(8)

In Eq. (8) we have implicitly performed the derivative before taking either of
the traces over the system or bath. These operations commute so we can either:
trace over both system and environment and then take a numerical derivative
(this is the counting field approach [33, 48, 35]); or take an analytic
derivative after tracing out the bath to get an expression purely in terms of
system correlation functions. We use the latter approach as it allows us to
calculate multiple different bath observables from a single influence
functional.

As an example, we can use the technique introduced here to derive the
following expression for the occupation of a specific mode (see App. B for
details):

$n_{q}(t)=n_{q}(0)+|g_{q}|^{2}\int_{0}^{t}\int_{0}^{t}\operatorname{Tr}[s(t^{\prime})s(t^{\prime\prime})\rho]\times\\\
\times\bigg{\\{}\cos[\omega_{q}(t^{\prime}-t^{\prime\prime})]-i\sin[\omega_{q}(t^{\prime}-t^{\prime\prime})]\coth(\frac{\omega_{q}}{2T})\bigg{\\}}\,dt^{\prime\prime}\,dt^{\prime},$
(9)

where $n_{q}(t)=\langle a_{q}^{\dagger}(t)a_{q}(t)\rangle$ and $T$ is the
temperature of the bath. We also include details of how this quantity can be
calculated for a simple model including just two modes in the environment in
App. C.

To calculate a single correlation function
$\operatorname{Tr}[s_{1}(t^{\prime})s_{2}(t^{\prime\prime})\rho]$ in the
integrand of Eq. (9) would typically require evolving $\rho$ to time
$t^{\prime\prime}$, acting with $s_{2}$, evolving to $t^{\prime}$ and finally
taking the expectation of $s_{1}$. To calculate the time evolution for each
required $t^{\prime}$ and $t^{\prime\prime}$ needed to evaluate the integral
in Eq. (9) is a computationally expensive procedure. To reduce this problem we
use the process tensor approach as this requires a single (albeit individually
more expensive) calculation from which any system observable can be determined
with little further cost. Below, we introduce the idea of a process tensor
conceptually and give an overview of how it can be calculated efficiently by
contracting a tensor network.

##  3 Process Tensors

The process tensor formalism is a general operational approach to non-
Markovian open quantum systems [44]. It supposes a finite set of control
operations act on the system at a sequence of time points, such as the
preparation of a particular state or a projective measurement for example. The
process tensor has two indices for each time at which a control operation can
be applied. It therefore allows the computation of any multi-time correlation
function of system operators by inserting control operations at the
corresponding times.

For Gaussian environments it has been shown that the process tensor can be
efficiently calculated by contracting a tensor network [41, 49]. This is done
by an alternative contraction order of the TEMPO network originally derived in
Ref. [40] as a representation of the quasi-adiabatic path integral (QUAPI)
method [50]. The process tensor fully captures the non-Markovian nature of the
interaction of the system with its environment and is therefore challenging to
compute. Once obtained, however, it can be used repeatedly to compute any
multi-time correlation at very little cost. With this, we have a tool at hand
to evaluate the occupation, Eq. (9), with moderate computational effort. In
the following, we briefly outline the form of the TEMPO tensor network and
then explain how it can be adapted to compute the process tensor. We refer the
reader to [40, 51] for more details on the derivation of the network and how
it can be efficiently contracted.

The TEMPO network can be understood in terms of superoperators beginning from
a discretisation of Eq. (3) with Eq. (4) substituted. We perform a symmetrised
Trotter splitting between the system and interaction parts:

$\rho_{S}(t_{N})=\overleftarrow{T}\mathcal{U}\prod_{k=0}^{N-1}\mathcal{U}^{2}\prod_{k=0}^{N-1}\prod_{k^{\prime}=0}^{k}\mathcal{I}_{k-k^{\prime}}(t_{k},t_{k^{\prime}})\mathcal{U}\rho_{S}(t_{0}),$
(10)

where we are now using the Schrödinger picture and have discretised the
evolution into $N$ timesteps of length $\Delta t$ and $t_{k}=k\Delta t$. Here
we assume $H_{S}$ is time independent such that
$\mathcal{U}=\text{e}^{\mathcal{L}_{S}\Delta t/2}$. For each timestep $t_{k}$
the terms
$\prod_{k^{\prime}=0}^{k}\mathcal{I}_{k-k^{\prime}}(t_{k},t_{k^{\prime}})$
account for the influence of the past system states on $t_{k}$ via the
environment. We can represent propagation forward 4 timesteps with the
following network:

(11)

Here the sum caps represent summing over the index of the tensor associated
with the capped leg i.e.

$\includegraphics[width=173.44534pt,clip={true},valign={c}]{tikzdiag3}.$ (12)

In practice the initial state is treated as a one site matrix product state
(MPS) and each timestep it is contracted with a matrix product operator (MPO)
corresponding to a row in the above network that grows the MPS by one site. In
other literature this MPS is referred to as the augmented density tensor
(ADT). The state is then read out at each timestep by summing over all but the
rightmost legs of the ADT. The process tensor can be identified as the above
network with the $\mathcal{U}$ propagators between each timestep disconnected
and the legs left open. Formally this can be represented as Eq. (10) without
performing the time-ordering. This network can then be contracted from left to
right and the process tensor is stored as an MPO. In practice this looks like:

(13)

Alternatively, one can view the system propagators $\mathcal{U}$ as control
operations and exclude them from the process tensor. Doing so can be utilised
to efficiently explore the dynamics for different system Hamiltonians [42]. In
this work, however, we include the system propagators in the process tensor
where arbitrary operations can be inserted between the open legs. In
particular, we can retrieve the system evolution by inserting identity
operations.

Direct contraction of Eq. (13), column-by-column, would result in an object
that grows exponentially with each column contracted. Many of the degrees of
freedom of this object, however, contribute negligibly and can be discarded
without significantly sacrificing accuracy [52, 53]. These irrelevant degrees
of freedom can be systematically identified by performing a singular value
decomposition (SVD) on each tensor in the chain and dispensing any components
with corresponding singular values below a threshold value. These SVD sweeps
are carried out once in each direction after contracting a column. For all of
the results in this paper we discarded components with singular values of
relative magnitude less than $10^{-8}$. The other main approximation we make
is a finite memory cutoff where the influence of the bath is only stored for
$K$ timesteps into the past [50, 54]. The value of $K$ depends on $\Delta t$:
for convergence we require that $K$ be increased until $K\Delta t$ exceeds the
correlation time of the bath. However $\Delta t$ must in turn be reduced to
minimise the error from the Trotter splitting. For all the results presented
here $\Omega\Delta t=0.05$ and $K=50$ were sufficient for convergence. These
parameter values lead to more intensive calculations than would typically
required for this form of coupling; this is because the calculation of bath
dynamics is highly sensitive to small errors in system expectation values.

The final process tensor is of the form shown in Eq. (13) which for $N$
timesteps is $N-1$ sites long; the bond dimension is related to the number of
relevant bath degrees of freedom [55]. From this a correlation such as
$\langle\mathcal{B}(t_{3})\mathcal{A}(t_{1})\rangle$, where $\mathcal{A}$ is
an arbitrary system superoperator, can be calculated by contracting the open
legs of the process tensor with $\mathcal{A}$ and $\mathcal{B}$ inserted at
the corresponding timesteps as follows:

(14)

For a product Hilbert space spanned by $\\{\ket{e_{i}}\otimes\ket{e_{j}}\\}$
the trace cap is defined as $\sum_{i}\ket{e_{i}}\otimes\ket{e_{i}}$, the
vectorised Hilbert space identity matrix and a null eigenvector of any
Liouvillian.

##  4 Biased Spin Boson Model

In this section we study a simple example system of a driven spin-1/2 coupled
to a continuum of bosonic modes, the biased spin-boson model [56, 45]. The
Hamiltonian of the system is:

$H=\epsilon s_{z}+\Omega
s_{x}+s_{z}\sum_{q}g_{q}(a_{q}+a_{q}^{\dagger})+\sum_{q}\omega_{q}a_{q}^{\dagger}a_{q},$
(15)

where we have introduced the spin-1/2 operators
$s_{z}=(\outerproduct{1}{1}-\outerproduct{0}{0})/2$ and
$s_{x}=(\outerproduct{1}{0}+\outerproduct{0}{1})/2$ which act on the states
$\\{\ket{0},\ket{1}\\}$. The transition between these states is driven
classically with strength $\Omega$ and bias $\epsilon$. The system is in turn
coupled with strength $g_{q}$ to a bath of bosonic modes of energy
$\omega_{q}$. The bosonic degrees of freedom are described by ladder operators
$a_{q}$ and $a_{q}^{\dagger}$ which satisfy the canonical commutation
relations: $[a_{q},a_{q^{\prime}}^{\dagger}]=\delta_{qq^{\prime}}$ and
$[a_{q},a_{q^{\prime}}]=[a_{q}^{\dagger},a_{q^{\prime}}^{\dagger}]=0$.

This kind of model can underpin a wide range of physical systems, for example
biological or molecular systems undergoing energy transport and interacting
with vibrational modes [57, 5], energy transfer in solid state systems [58],
superconducting qubits in microwave resonators [59], or quantum dots
interacting with a micromechanical resonator [60].

The environment in this model is characterised by its spectral density which
captures both the coupling to each mode $g_{q}$ and the density of states of
the environment. It is defined as:

$J(\omega)=\sum_{q}|g_{q}|^{2}\delta(\omega-\omega_{q}).$ (16)

This can be related to the auto-correlation function of the bath coupling
operator $B=\sum_{q}g_{q}(a_{q}+a_{q}^{\dagger})$,

$\displaystyle C(t)=\left\langle
B(t)B(0)\right\rangle=\int_{0}^{\infty}C(\omega,t)\,d\omega,$ (17)

where we have introduced

$C(\omega,t)=J(\omega)\bigg{[}\cos(\omega
t)\coth\left(\frac{\omega}{2T}\right)-i\sin(\omega t)\bigg{]}$ (18)

and $T$ is the temperature of the bath.

The spectral density is typically given a functional form informed by
experimental measurements. Here we use an Ohmic spectral density with an
exponential cutoff [45]:

$J(\omega)=2\alpha\omega\text{e}^{-\omega/\omega_{c}}$ (19)

where $\omega_{c}$ is the cut-off frequency and $\alpha$ a dimensionless
coupling constant.

###  4.1 Steady State Properties

Here we consider the change in energy of each bath mode, defined as:

$\Delta Q_{q}(t)=\omega_{q}\left[\langle
a_{q}^{\dagger}(t)a_{q}(t)\rangle-\langle
a_{q}^{\dagger}(0)a_{q}(0)\rangle\right].$ (20)

However, this is ill-defined in the continuum limit. We cannot refer to the
energy change of a specific mode of the environment as the coupling to any
single mode is infinitesimal. Instead we consider the change over a narrow
band of modes:

$\Delta
Q(\omega,t)=i\int_{0}^{t}\int_{0}^{t}\operatorname{Tr}[s(t^{\prime})s(t^{\prime\prime})\rho]\int_{\omega-\delta/2}^{\omega+\delta/2}C_{t}(\omega^{\prime},t^{\prime}-t^{\prime\prime})\,d\omega^{\prime}\,dt^{\prime\prime}\,dt^{\prime}.$
(21)

Here the frequency interval is of width $\delta$ and centered on $\omega$; the
notation $C_{t}$ signifies a derivative with respect to $t$.

For the results in this section we choose a relatively weak system-bath
coupling. This is not due to any inherent limitations of the methods used but
rather because the behaviour in this limit is more interesting. At stronger
coupling contributions from the re-organisation energy of the environment
dominate and the behavior reduces to that of the exactly solvable independent
boson model. This effect has been seen in previous calculations involving the
total heat of the bath [33].

Figure 1:  Frequency resolved steady-state variation of heat in bath as a
function of bias $\epsilon$ (main panel) and the total heat exchanged with the
bath compared with Eq. (30) (upper panel). Other parameters are $T=\Omega$,
$\alpha=0.05$, $\omega_{c}=10\Omega$ and $\delta=0.1\Omega$. Figure 2:  Steady
state change in bath heat as a function of temperature for varying initial
state where (a) $p=1.0$, (b) $p=0.9$ and (c) $p=0.8$. Here we fixed
$\epsilon=1.5\Omega$ and all other parameters are unchanged from Figure 1.

In the main panel of Fig. 1 we show the period-averaged steady state value for
$\Delta Q$ as a function of the bias $\epsilon$ with the two-level system
initialised in the low energy state, $\ket{0}$. The averaging is performed
over a single oscillation of each mode such that we calculate,

$\Delta Q(\omega)=\frac{\omega}{2\pi}\int_{t_{SS}}^{t_{SS}+2\pi/\omega}\Delta
Q(\omega,t^{\prime})\,dt^{\prime},$ (22)

where $t_{SS}$ is chosen to be long enough that the system density operator is
stationary $\dot{\rho}_{S}(t\geq t_{SS})\simeq 0$. We do this because even
after the system degrees of freedom reach a stationary state there are still
oscillations in the mode occupations [43]. This quantity is then able to probe
the way in which energy is removed from or added to the environment to
populate the spin-1/2 in its steady state configuration. We see that for small
bias the environment gains energy from the system while at large bias the
opposite is true.

To confirm the accuracy of this approach we also checked that the total energy
was conserved. This required calculation of the total change in $\langle
H_{S}\rangle$, $\langle H_{I}\rangle$ and $\langle H_{B}\rangle$ which we will
denote $\Delta Q_{S}$, $\Delta Q_{I}$ and $\Delta Q_{B}$ respectively. The
change in system energy, $\Delta Q_{S}$ is trivial to calculate from the
dynamics. The total energy exchanged with the bath, $\Delta Q_{B}$, is given
by:

$\Delta
Q_{B}=i\int_{0}^{t_{SS}}\int_{0}^{t_{SS}}\operatorname{Tr}[s(t^{\prime})s(t^{\prime\prime})\rho]\times\\\
\times C_{t}(t^{\prime}-t^{\prime\prime})\,dt^{\prime\prime}\,dt^{\prime}.$
(23)

The bath will generally have a finite memory time $t_{K}$ such that the auto-
correlation function satisfies $C(t>t_{K})=0\implies C_{t}(t>t_{K})=0$.
Therefore if the total energy exchanged is the relevant quantity then only
correlation functions that span this memory time are required.

A similar expression can be derived for $\Delta Q_{I}$ again using the
generating functional formalism. This is given by:

$\Delta
Q_{I}=2\imaginary\left[\int_{0}^{t_{SS}}\operatorname{Tr}[s(t)s(t^{\prime})\rho]\times\right.\\\
\left.\times C(t-t^{\prime})\,dt^{\prime}\vphantom{\int_{0}^{t_{SS}}}\right].$
(24)

All of the bath observables calculated here can be derived from the same set
of correlation functions. The only variation is in the form of the kernel used
for the integral transformation. Generally the dynamics of an $n^{\text{th}}$
order bath correlation function up to time $t$ depends on the set of
$n^{\text{th}}$ order system correlation functions up to that point.

In the upper panel of Fig. 1 we compare the numerical calculation of the heat
gained by the bath, $\Delta Q_{B}$ with the remaining contribution $-(\Delta
Q_{I}+\Delta Q_{S})$. We see that, up to small numerical inaccuracies, these
quantities are effectively equal as would be expected from energy
conservation.

In Ref. [33] it was shown that for weak coupling $\Delta Q_{B}$ can be well
approximated by

$\Delta Q_{B}\approx E_{r}-\Delta Q_{S}^{G},$ (25)

where $E_{r}$ is the reorganisation energy of the bath given by

$E_{r}=\int_{0}^{\infty}d\omega^{\prime}\frac{J(\omega^{\prime})}{\omega^{\prime}}=\alpha\omega_{c},$
(26)

and $\Delta Q_{S}^{G}$ is approximate change in system energy defined as

$\Delta Q_{S}^{G}=\operatorname{Tr}\\{H_{S}[\rho_{SS}(T)-\rho_{S}(0)]\\},$
(27)

where $\rho_{SS}(T)$ is the reduced Gibbs state of the system Hamiltonian,

$\rho_{S}(t_{SS})\approx\rho_{SS}(T)=\frac{\text{e}^{-H_{S}/T}}{\mathrm{Tr}[\text{e}^{-H_{S}/T}]}.$
(28)

Since the coupling to the environment is weak we expect there to be no
significant correlations between the system and environment such that the
system simply reaches a state in thermal equilibrium with respect to $H_{S}$
at the temperature of the environment.

From this we can then find an analytic expression for $\Delta Q_{S}^{G}$:

$\Delta
Q_{S}^{G}=\frac{1}{2}\left(\epsilon-\tilde{\Omega}\tanh\left(\frac{\tilde{\Omega}}{2T}\right)\right),$
(29)

where $\tilde{\Omega}=\sqrt{\Omega^{2}+\epsilon^{2}}$ is the generalised Rabi
frequency and we have used $\rho_{S}(0)=\ket{0}\bra{0}$. This in turn allows
us to derive an (approximate) analytic expression for $\Delta Q_{B}$ by
substituting (29) and (26) into (25) giving

$\Delta
Q_{B}\approx\frac{1}{2}\left(\alpha\omega_{c}-\epsilon+\tilde{\Omega}\tanh\left(\frac{\tilde{\Omega}}{2T}\right)\right).$
(30)

In the upper panel of Fig. 1 we compare the analytic approximation given by
Eq. (30) with the numerical results and find good agreement as would be
expected at weak system environment coupling.

From the form of Eq (29) we can identify a crossover temperature $T^{*}$,
above which the system energy increases and vice-versa. This is given by the
solution of

$\operatorname{Tr}_{S}\\{H_{S}[\rho_{SS}(T^{*})-\rho_{S}(0)]\\}=0.$ (31)

For initial states that are incoherent mixtures of the two basis states
$\rho_{S}(0)=p\ket{0}\bra{0}+(1-p)\ket{1}\bra{1}$ we can solve this to find

$T^{*}=\frac{\tilde{\Omega}}{2\tanh^{-1}\left[\epsilon(2p-1)/\tilde{\Omega}\right]}.$
(32)

In Fig. 2 we show the steady state heat change in the environment $\Delta
Q(\omega)$ as a function of temperature for three values of $p$. From this we
can see that the analytic expression is a good predictor of whether the system
has absorbed or emitted more energy into the bath at a specific frequency.

Figure 3:  Movement of heat around the bath over time using the driving
routine described in Eq. (34) with $\Omega t_{1}=100$, $\epsilon_{1}=2\Omega$,
$\Omega t_{2}=200$, $\epsilon_{2}=3\Omega$ and $\delta=0.01\Omega$. The other
parameters are the same as Fig. 1. The upper panel shows the system energy.

###  4.2 Dynamics

We now move on to show how the methods described above can also be used to
access the complete dynamics of the bath. To do this we implement a time-
dependent driving sequence with the goal of demonstrating the absorption of
energy from the bath over a range of modes followed by emission into the bath
over a different range of modes. This shows how, by controlling the dynamics
of the system, we can control the state of the environment. We initialize the
system in the $\ket{0}$ state at $t=0$ and allow the bias to vary over time
with the system Hamiltonian

$H_{S}(t)=\epsilon(t)s_{z}+\Omega s_{x}.$ (33)

We use a linear ramp for the bias and hence set

$\epsilon(t)=\begin{cases}\epsilon_{1}&0\leq t<t_{1}\\\
\epsilon_{1}+(\epsilon_{2}-\epsilon_{1})\frac{t-t_{1}}{t_{2}-t_{1}}&t_{1}\leq
t\leq t_{2}\\\ \epsilon_{2}&t>t_{2}\end{cases}.$ (34)

We choose $t_{1}$ to be large enough to allow the system to first reach
thermal equilibrium with the environment and $\epsilon_{1}$ is chosen such
that this equilibrium state is higher in energy than the initial state and
thus absorbs energy from the bath. The bias is then increased to push the
system into a higher temperature state such that energy is emitted into the
bath to re-equilibrate. As can be seen in the upper panel of Fig. 3 during the
linear ramp the energy of the system decreases, emitting energy into the bath.
This can be seen in the heat distribution of the bath as a band of modes which
increase in energy. The modes targeted correspond to frequencies close to the
(time-dependent) Rabi frequency
$\tilde{\Omega}(t)=\sqrt{\epsilon(t)^{2}+\Omega^{2}}$ which changes over the
course of the linear ramp. This shows that by using a time varying drive,
along with the intuition we have built up about the behaviour of the
environment, we are able to control the (frequency dependent) absorption and
emission of energy into and out of the bath.

##  5 Conclusions & Outlook

In this paper we have introduced a method for calculating any dynamical
observable of a Gaussian bath purely in terms of correlation functions of a
linearly coupled system. For systems coupled to infinite baths we demonstrated
that the process tensor formalism provides a framework to calculate, cheaply,
the large number of system correlation functions needed. In general it is
numerically costly to calculate a process tensor. However, once it exists it
can then be used to completely characterise any time dependent observable or
any multi-point correlation function of the system, bath, or a combination of
both.

Specifically, we were able to calculate how energy is moved between different
environment modes and a spin-1/2 particle in the biased spin boson model. We
showed that, for weak system environment coupling, the spin reaches thermal
equilibrium with the bath removing or adding energy to modes which correspond
to the system frequency. We also saw how the techniques described here can
allow us to track the flow of energy into and away from the system when it is
subject to a time dependent driving.

The framework we presented is not specific to any technique for calculating
the required correlation functions and opens up many avenues to analyze a
broad range of open quantum systems. In particular, when the environment
induces complex non-Markovian dynamics in the system this method has the
potential to reveal the origin of complicated dynamical effects, and to design
new systems that could control and harness these effects.

## Acknowledgements

DG and GF acknowledge studentship funding from EPSRC (EP/L015110/1). AS
acknowledges support from the Australian Research Council Centre of Excellence
for Engineered Quantum Systems (EQUS, CE170100009). We acknowledge support
from EPSRC (EP/T014032/1). This work used EPCC’s Cirrus HPC Service
(https://www.epcc.ed.ac.uk/cirrus).

## References

  * [1] H. P. Breuer and F. Petruccione  “The Theory of Open Quantum Systems”  Oxford University Press, 2002
  * [2] A. W. Chin et al.  “The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment–protein complexes”  In _Nat. Phy._ 9.2 Springer ScienceBusiness Media LLC, 2013, pp. 113–118  DOI: 10.1038/nphys2515
  * [3] M. Thorwart et al.  “Enhanced quantum entanglement in the non-Markovian dynamics of biomolecular excitons”  In _Chem. Phys. Lett._ 478.4-6 Elsevier BV, 2009, pp. 234–237  DOI: 10.1016/j.cplett.2009.07.053
  * [4] M. Mohseni, A. Shabani, S. Lloyd and H. Rabitz  “Energy-scales convergence for optimal and robust quantum transport in photosynthetic complexes”  In _J. Chem. Phys._ 140.3 AIP Publishing, 2014, pp. 035102  DOI: 10.1063/1.4856795
  * [5] Marco Rey, Alex W. Chin, Susana F. Huelga and Martin B. Plenio  “Exploiting Structured Environments for Efficient Energy Transfer: The Phonon Antenna Mechanism”  In _J. Phys. Chem. Lett._ 4.6 American Chemical Society (ACS), 2013, pp. 903–907  DOI: 10.1021/jz400058a
  * [6] Christine Maier et al.  “Environment-Assisted Quantum Transport in a 10-qubit Network”  In _Phys. Rev. Lett._ 122.5 American Physical Society (APS), 2019  URL: https://doi.org/10.1103%2Fphysrevlett.122.050501
  * [7] Dvira Segal and Bijay Kumar Agarwalla  “Vibrational Heat Transport in Molecular Junctions”  In _Ann. Rev. Phys. Chem._ 67.1 Annual Reviews, 2016, pp. 185–209  DOI: 10.1146/annurev-physchem-040215-112103
  * [8] S. Gröblacher et al.  “Observation of non-Markovian micromechanical Brownian motion”  In _Nat. Commun._ 6.1 Springer ScienceBusiness Media LLC, 2015  URL: https://doi.org/10.1038%2Fncomms8606
  * [9] Anton Potočnik et al.  “Studying light-harvesting models with superconducting circuits”  In _Nat. Commun._ 9.1 Springer Nature, 2018, pp. 904  URL: https://doi.org/10.1038%2Fs41467-018-03312-x
  * [10] Heng-Na Xiong, Yi Li, Yixiao Huang and Zichun Le  “Photonic Crystal with Infinite Cavity-array Structure”  In _Quant. Inf. & Comp._ 18, 2018, pp. 267–284  DOI: 10.26421/QIC18.3-4-6
  * [11] Inés Vega and Daniel Alonso  “Dynamics of non-Markovian open quantum systems”  In _Rev. Mod. Phys._ 89.1 American Physical Society (APS), 2017, pp. 015001  URL: https://doi.org/10.1103/revmodphys.89.015001
  * [12] H. P. Breuer, E. M. Laine, J. Piilo and B. Vacchini  “Colloquium: Non-Markovian dynamics in open quantum systems”  In _Rev. Mod. Phys._ 88 American Physical Society, 2016, pp. 021002  DOI: 10.1103/RevModPhys.88.021002
  * [13] Mark T Mitchison  “Quantum thermal absorption machines: refrigerators, engines and clocks”  In _Contemp. Phys._ 60.2 Taylor & Francis, 2019, pp. 164–187  DOI: 10.1080/00107514.2019.1631555
  * [14] “Thermodynamics in the Quantum Regime”  Springer, Cham, 2018  DOI: 10.1007/978-3-319-99046-0
  * [15] M. Brenes et al.  “Tensor-Network Method to Simulate Strongly Interacting Quantum Thermal Machines”  In _Phys. Rev. X_ 10.3 APS, 2020, pp. 031040  DOI: 10.1103/PhysRevX.10.031040
  * [16] Philipp Strasberg, Gernot Schaller, Neill Lambert and Tobias Brandes  “Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping”  In _New J. Phys._ 18.7 IOP Publishing, 2016, pp. 073007  DOI: 10.1088/1367-2630/18/7/073007
  * [17] Gerald D Mahan  “Many-particle physics”  Springer Science & Business Media, 2013
  * [18] Robert Silbey and Robert A. Harris  “Variational calculation of the dynamics of a two level system interacting with a bath”  In _J. Chem. Phys._ 80.6 AIP Publishing, 1984, pp. 2615–2617  DOI: 10.1063/1.447055
  * [19] R. A. Harris and R. Silbey  “Variational calculation of the tunneling system interacting with a heat bath. II. Dynamics of an asymmetric tunneling system”  In _J. Chem. Phys._ 83.3 AIP Publishing, 1985, pp. 1069–1074  DOI: 10.1063/1.449469
  * [20] Jozef T Devreese and Alexandre S Alexandrov  “Fröhlich polaron and bipolaron: recent developments”  In _Rep. Prog. Phys._ 72.6 IOP Publishing, 2009, pp. 066501  DOI: 10.1088/0034-4885/72/6/066501
  * [21] Alexey Kavokin and Guillaume Malpuech  “Cavity polaritons”  Elsevier, 2003
  * [22] Peter Kirton and Jonathan Keeling  “Superradiant and lasing states in driven-dissipative Dicke models”  In _New J. Phys._ 20.1 IOP Publishing, 2018, pp. 015009  DOI: 10.1088/1367-2630/aaa11d
  * [23] Alicia J Kollár et al.  “Supermode-density-wave-polariton condensation with a Bose–Einstein condensate in a multimode cavity”  In _Nat. commun._ 8.1 Nature Publishing Group, 2017, pp. 1–10  DOI: 10.1038/ncomms14386
  * [24] Martin Gärttner et al.  “Measuring out-of-time-order correlations and multiple quantum spectra in a trapped-ion quantum magnet”  In _Nat. Phy._ 13.8 Nature Publishing Group, 2017, pp. 781–786  DOI: 10.1038/nphys4119
  * [25] Jun Li et al.  “Measuring out-of-time-order correlators on a nuclear magnetic resonance quantum simulator”  In _Phys. Rev. X_ 7.3 APS, 2017, pp. 031011  DOI: 10.1103/PhysRevX.7.031011
  * [26] Mohamad Niknam, Lea F. Santos and David G. Cory  “Sensitivity of quantum information to environment perturbations measured with a nonlocal out-of-time-order correlation function”  In _Phys. Rev. Res._ 2 American Physical Society, 2020, pp. 013200  DOI: 10.1103/PhysRevResearch.2.013200
  * [27] Jake Iles-Smith, Neill Lambert and Ahsan Nazir  “Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems”  In _Phys. Rev. A_ 90.3 American Physical Society (APS), 2014, pp. 032114  URL: https://doi.org/10.1103/physreva.90.032114
  * [28] Jake Iles-Smith, Arend G. Dijkstra, Neill Lambert and Ahsan Nazir  “Energy transfer in structured and unstructured environments: Master equations beyond the Born-Markov approximations”  In _J. Chem. Phys._ 144.4 AIP Publishing, 2016, pp. 044110  DOI: 10.1063/1.4940218
  * [29] Neill Lambert, Shahnawaz Ahmed, Mauro Cirio and Franco Nori  “Modelling the ultra-strongly coupled spin-boson model with unphysical modes”  In _Nat. Commun._ 10 Nature Publishing Group, 2019, pp. 3721  DOI: 10.1038/s41467-019-11656-1
  * [30] Lili Zhu, Hao Liu, Weiwei Xie and Qiang Shi  “Explicit system-bath correlation calculated using the hierarchical equations of motion method”  In _J. Chem. Phys._ 137.19 American Institute of Physics, 2012, pp. 194106  DOI: 10.1063/1.4766358
  * [31] Linze Song and Qiang Shi  “Hierarchical equations of motion method applied to nonequilibrium heat transport in model molecular junctions: Transient heat current and high-order moments of the current operator”  In _Phys. Rev. B_ 95.6 APS, 2017, pp. 064308  DOI: 10.1103/PhysRevB.95.064308
  * [32] C Schinabeck, R Härtle and M Thoss  “Hierarchical quantum master equation approach to electronic-vibrational coupling in nonequilibrium transport through nanosystems: Reservoir formulation and application to vibrational instabilities”  In _Phys. Rev. B_ 97.23 APS, 2018, pp. 235429  DOI: 10.1103/PhysRevB.97.235429
  * [33] Maria Popovic et al.  “Non-equilibrium quantum thermodynamics with time-evolving matrix product operators”, 2021  arXiv:2008.06491 [quant-ph]
  * [34] Massimiliano Esposito, Upendra Harbola and Shaul Mukamel  “Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems”  In _Rev. Mod. Phys._ 81.4 APS, 2009, pp. 1665  DOI: 10.1103/RevModPhys.81.1665
  * [35] Michael Kilgour, Bijay Kumar Agarwalla and Dvira Segal  “Path-integral methodology and simulations of quantum thermal transport: Full counting statistics approach”  In _J. Chem. Phys._ 150.8 AIP Publishing LLC, 2019, pp. 084111  DOI: 10.1063/1.5084949
  * [36] Javier Prior, Alex W. Chin, Susana F. Huelga and Martin B. Plenio  “Efficient Simulation of Strong System-Environment Interactions”  In _Phys. Rev. Lett._ 105 American Physical Society, 2010, pp. 050404  DOI: 10.1103/PhysRevLett.105.050404
  * [37] Dario Tamascelli  “Excitation dynamics in chain-mapped environments”  In _Entropy_ 22.11 Multidisciplinary Digital Publishing Institute, 2020, pp. 1320  DOI: 10.3390/e22111320
  * [38] Florian AYN Schröder and Alex W Chin  “Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution”  In _Phys. Rev. B_ 93.7 APS, 2016, pp. 075105  DOI: 10.1103/PhysRevB.93.075105
  * [39] C Gonzalez-Ballestero, Florian AYN Schröder and Alex W Chin  “Uncovering nonperturbative dynamics of the biased sub-Ohmic spin-boson model with variational matrix product states”  In _Phys. Rev. B_ 96.11 APS, 2017, pp. 115427  DOI: 10.1103/PhysRevB.96.115427
  * [40] A. Strathearn et al.  “Efficient non-Markovian quantum dynamics using time-evolving matrix product operators”  In _Nat. Commun._ 9.1 Springer Nature, 2018, pp. 3322  URL: https://doi.org/10.1038%2Fs41467-018-05617-3
  * [41] Mathias R. Jørgensen and Felix A. Pollock  “Exploiting the Causal Tensor Network Structure of Quantum Processes to Efficiently Simulate Non-Markovian Path Integrals”  In _Phys. Rev. Lett._ 123 American Physical Society, 2019, pp. 240602  DOI: 10.1103/PhysRevLett.123.240602
  * [42] Gerald E. Fux et al.  “Efficient Exploration of Hamiltonian Parameter Space for Optimal Control of Non-Markovian Open Quantum Systems”  In _Phys. Rev. Lett._ 126 American Physical Society, 2021, pp. 200401  DOI: 10.1103/PhysRevLett.126.200401
  * [43] Dominic Gribben et al.  “Exact quantum dynamics in structured environments”  In _Phys. Rev. Res._ 2.1 APS, 2020, pp. 013265  DOI: 10.1103/PhysRevResearch.2.013265
  * [44] Felix A. Pollock et al.  “Non-Markovian quantum processes: Complete framework and efficient characterization”  In _Phys. Rev. A_ 97 American Physical Society, 2018, pp. 012127  DOI: 10.1103/PhysRevA.97.012127
  * [45] A. J. Leggett et al.  “Dynamics of the dissipative two-state system”  In _Rev. Mod. Phys._ 59 American Physical Society, 1987, pp. 1–85  DOI: 10.1103/RevModPhys.59.1
  * [46] Shaul Mukamel  “Principles of nonlinear optical spectroscopy”  Oxford University Press, New York, 1995
  * [47] R.P Feynman and F.L Vernon  “The theory of a general quantum system interacting with a linear dissipative system”  In _Ann. Phys._ 24, 1963, pp. 118–173  DOI: https://doi.org/10.1016/0003-4916(63)90068-X
  * [48] Mihail Silaev, Tero T. Heikkilä and Pauli Virtanen  “Lindblad-equation approach for the full counting statistics of work and heat in driven quantum systems”  In _Phys. Rev. E_ 90 American Physical Society, 2014, pp. 022103  DOI: 10.1103/PhysRevE.90.022103
  * [49] The TEMPO collaboration  “TimeEvolvingMPO: A Python 3 package to efficiently compute non-Markovian open quantum systems.”  GitHub, 2020  URL: {https://github.com/tempoCollaboration/TimeEvolvingMPO}
  * [50] Nancy Makri and Dmitrii E. Makarov  “Tensor propagator for iterative quantum time evolution of reduced density matrices. I. Theory”  In _J. Chem. Phys._ 102.11 AIP Publishing, 1995, pp. 4600–4610  DOI: 10.1063/1.469508
  * [51] Aidan Strathearn  “Modelling Non-Markovian Quantum Systems Using Tensor Networks”  Springer Nature, 2020
  * [52] Román Orús  “A practical introduction to tensor networks: Matrix product states and projected entangled pair states”  In _Ann. Phys._ 349 Elsevier BV, 2014, pp. 117–158  DOI: 10.1016/j.aop.2014.06.013
  * [53] Ulrich Schollwöck  “The density-matrix renormalization group in the age of matrix product states”  In _Ann. Phys._ 326.1 Elsevier BV, 2011, pp. 96–192  DOI: 10.1016/j.aop.2010.09.012
  * [54] Nancy Makri and Dmitrii E. Makarov  “Tensor propagator for iterative quantum time evolution of reduced density matrices. II. Numerical methodology”  In _J. Chem. Phys._ 102.11 AIP Publishing, 1995, pp. 4611–4618  DOI: 10.1063/1.469509
  * [55] M. Cygorek et al.  “Numerically-exact simulations of arbitrary open quantum systems using automated compression of environments”, 2021  arXiv:2101.01653 [quant-ph]
  * [56] A. J. Leggett  “Quantum tunneling in the presence of an arbitrary linear dissipation mechanism”  In _Phys. Rev. B_ 30.3 American Physical Society (APS), 1984, pp. 1208–1218  DOI: 10.1103/physrevb.30.1208
  * [57] S. Rackovsky and R. Silbey  “Electronic energy transfer in impure solids: I. Two molecules embedded in a lattice”  In _Mol. Phys._ 25.1 Taylor & Francis, 1973, pp. 61–72  DOI: 10.1080/00268977300100081
  * [58] Emil Rozbicki and Paweł Machnikowski  “Quantum Kinetic Theory of Phonon-Assisted Excitation Transfer in Quantum Dot Molecules”  In _Phys. Rev. Lett._ 100 American Physical Society, 2008, pp. 027401  DOI: 10.1103/PhysRevLett.100.027401
  * [59] Huaixiu Zheng and Harold U. Baranger  “Persistent Quantum Beats and Long-Distance Entanglement from Waveguide-Mediated Interactions”  In _Phys. Rev. Lett._ 110 American Physical Society, 2013, pp. 113601  DOI: 10.1103/PhysRevLett.110.113601
  * [60] I. Yeo et al.  “Strain-mediated coupling in a quantum dot–mechanical oscillator hybrid system”  In _Nat. Nanotech._ 9 Nature Publishing Group SN -, 2013, pp. 106 EP–  URL: https://doi.org/10.1038/nnano.2013.274

##  Appendix A  Superoperator formalism

To calculate the evolution of the total system-bath density matrix we start
with the von Neumann equation (in the interaction picture) and use the
superoperator formalism:

$\frac{d}{dt}\rho(t)=-i[H_{I}(t),\rho(t)]\rightarrow\frac{d}{dt}\ket{\rho(t)}\rangle=\mathcal{L}_{I}(t)\ket{\rho(t)}\rangle$
(35)

where $\mathcal{L}_{I}(t)$ is the superoperator associated with taking the
commutator with $-iH_{I}(t)$ i.e.
$\mathcal{L}_{I}(t)\ket{\rho}\rangle=-i\ket{[H_{I}(t),\rho]}\rangle$. We can
then formally integrate this and arrive at the expression:

$\ket{\rho(t)}\rangle=\overleftarrow{T}\text{e}^{\int_{0}^{t}dt^{\prime}\mathcal{L}_{I}(t^{\prime})}\ket{\rho(0)}\rangle.$
(36)

To find the dynamics of reduced system we trace over the bath to leave

$\ket{\rho_{S}(t)}\rangle=\left\langle\overleftarrow{T}\text{e}^{\int_{0}^{t}dt^{\prime}\mathcal{L}_{I}(t^{\prime})}\right\rangle_{B}\ket{\rho_{S}(0)}\rangle,$
(37)

where we have made the assumption of a separable initial state such that
$\ket{\rho(0)}\rangle=\ket{\rho_{S}(0)}\rangle\otimes\ket{\rho_{B}(0)}\rangle$
and $\langle\cdot\rangle_{B}$ represents an expectation taken with respect to
$\ket{\rho_{B}(0)}\rangle$.

Throughout this paper we represent operators in Hilbert space as lowercase
letters in normal italics (e.g. $s$ and $x$) and superoperators in Liouville
space superoperators as uppercase calligraphic letters (e.g. $\mathcal{L}$).
If a superoperator corresponds to simple extension of an operator then it will
be labelled as such, for example:

$[a,\rho]=a\rho-\rho a\rightarrow\mathcal{A}^{-}\ket{\rho}\rangle.$ (38)

In a similar vein we can define the left-acting superoperator
$a\rho\rightarrow\mathcal{A}^{L}\ket{\rho}\rangle$, the right-acting
superoperator $\rho a\rightarrow\mathcal{A}^{R}\ket{\rho}\rangle$ and the
anti-commutator superoperator
$\\{a,\rho\\}\rightarrow\mathcal{A}^{+}\ket{\rho}\rangle$. The only exception
is that when referring to superoperators which correspond to Hamiltonians in
Hilbert space we will use the typical $\mathcal{L}$ symbol. Although the time-
ordering is carried out on superoperators we reserve the calligraphic
$\mathcal{T}$ for the trace operation.

In Hilbert space the trace of an operator returns a scalar. In Liouville space
the density matrix is vectorised and thus to return a scalar the trace must
correspond to an inner product. To see how we define this we give a brief
overview of the formalism of superoperators.

Consider a Hilbert space spanned by the vectors $\\{\ket{e_{i}}\\}$ a density
matrix can be written as a vector in the corresponding Liouville space by
using the mapping

$\rho=\sum_{i,j}\rho_{ij}\outerproduct{e_{i}}{e_{j}}\rightarrow\ket{\rho}\rangle=\sum_{i,j}\rho_{ij}\ket{e_{i}}\otimes\ket{e_{j}}.$
(39)

An equivalent mapping can be carried out on any operator in Hilbert space. We
can then define the trace as a dual vector in the Liouville space which acts
as follows

$\operatorname{Tr}[\rho]=\sum_{k}\sum_{i,j}\rho_{ij}\braket{e_{k}\outerproduct{e_{i}}{e_{j}}e_{k}}\rightarrow\langle\braket{\mathcal{T}}{\rho}\rangle=\sum_{k}\sum_{i,j}\rho_{ij}\braket{e_{k}}{e_{i}}\otimes\braket{e_{k}}{e_{j}}.$
(40)

From this we can see that the correct dual vector is

$\langle\bra{\mathcal{T}}\equiv\sum_{k}\bra{e_{k}}\otimes\bra{e_{k}}.$ (41)

Now consider a bipartite system $\mathcal{H}_{A}\otimes\mathcal{H}_{B}$
spanned by $\\{\ket{f_{i}}\otimes\ket{h_{k}}\\}$. The corresponding Liouville
space is in turn spanned by
$\\{\ket{f_{i}}\otimes\ket{f_{j}}\otimes\ket{h_{k}}\otimes\ket{h_{l}}\\}$. The
partial trace over $\mathcal{H}_{B}$ is then

$\mathrm{Tr}_{B}[\cdot]\rightarrow\langle\bra{\mathcal{T}_{B}}\equiv\sum_{i,j,k}\outerproduct{f_{i}}{f_{i}}\otimes\outerproduct{f_{j}}{f_{j}}\otimes\bra{h_{k}}\otimes\bra{h_{k}}.$
(42)

We now define the partial trace over a series of superoperators as:

$\langle\bra{\mathcal{T}_{B}}\mathcal{A}(t_{A})\mathcal{B}(t_{B})\mathcal{C}(t_{C})\dots\ket{\rho_{B}}\rangle\equiv\braket{\mathcal{A}(t_{A})\mathcal{B}(t_{B})\mathcal{C}(t_{C})\dots}_{B}.$
(43)

We will reserve the use of angular brackets to refer to expectations of
superoperators. If there is no subscript after the angular brackets then the
trace is over the entire space.

If we consider the bath to initially be in a Gaussian state, e.g. thermal
equilibrium, then the trace over the bath can be simplified using Wick’s
theorem to give:

$\left\langle\text{e}^{\int_{0}^{t}dt^{\prime}\mathcal{L}_{I}(t^{\prime})}\right\rangle_{B}=\exp\left(\frac{1}{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}\left\langle\overleftarrow{T}\mathcal{L}_{I}(t^{\prime})\mathcal{L}_{I}(t^{\prime\prime})\right\rangle_{B}\right).$
(44)

The key result of this paper relies on an equivalent relation holding when
source terms linear in the ladder operators are added to the exponent.

##  Appendix B Number operator derivation

In this section we derive Eq. (9) of the main text. For simplicity we consider
a bath that consists of just a single mode such that $H_{I}=gs(a^{\dagger}+a)$
and $H_{B}=\omega a^{\dagger}a$, where $s$ is an arbitrary system operator.
The generalisation to a multimode bath is straightforward. We wish to
calculate $n(t)=\operatorname{Tr}[a^{\dagger}(t)a(t)\rho]$ and write this
quantity in terms of expectation values of superoperators as follows,

$n(t)=\left\langle\overleftarrow{T}\mathcal{A}^{L\dagger}(t)\mathcal{A}^{L}(t)\text{e}^{\int_{0}^{t}\mathcal{L}_{I}(t^{\prime})dt^{\prime}}\right\rangle.$
(45)

Now we can write

$\displaystyle\mathcal{A}^{L\dagger}(t)\mathcal{A}^{L}(t)$
$\displaystyle=\left.\frac{d^{2}}{d\Lambda
d\Lambda^{*}}\text{e}^{\Lambda^{*}\mathcal{A}^{L\dagger}(t)}\text{e}^{\Lambda\mathcal{A}^{L}(t)}\right|_{\Lambda,\Lambda^{*}=0}$
$\displaystyle=\left.\frac{d^{2}}{d\Lambda
d\Lambda^{*}}\text{e}^{\Lambda^{*}\mathcal{A}^{L\dagger}(t)+\Lambda\mathcal{A}^{L}(t)-\frac{|\Lambda|^{2}}{2}}\right|_{\Lambda,\Lambda^{*}=0},$
(46)

where we have combined the exponentials in the second line using the Baker-
Campbell-Hausdorff (BCH) formula and
$[\mathcal{A}^{L\dagger}(t)\mathcal{A}^{L}(t)]=[a^{\dagger}(t),a(t)]=-1$. Care
must be taken when explicitly calculating the derivatives above since the
operators involved do not commute. Inserting Eq. (B) into Eq. (45) gives

$n(t)=\left.\frac{d^{2}}{d\Lambda
d\Lambda^{*}}\left\langle\overleftarrow{T}\text{e}^{\Lambda\mathcal{A}^{L}(t)+\Lambda^{*}\mathcal{A}^{L\dagger}(t)-\frac{|\Lambda|^{2}}{2}}\text{e}^{\int_{0}^{t}\mathcal{L}_{I}(t^{\prime})dt^{\prime}}\right\rangle\right|_{\Lambda,\Lambda^{*}=0}.$
(47)

We can combine the exponentials in Eq. (47) without using the BCH formula by
virtue of the identity

$\overleftarrow{T}\left[\mathcal{B}(t),\int\mathcal{C}(t^{\prime})dt^{\prime}\right]=0,$
(48)

which holds for arbitrary $\mathcal{B}$, $\mathcal{C}$, and integral limits.
This is because the portion of the integral that contains $\mathcal{C}(t)$,
which in general does not commute with $\mathcal{B}(t)$ under time ordering,
is infinitesimal, i.e. has measure zero. Thus, Eq. (47) becomes

$\displaystyle n(t)=\left.\frac{d^{2}}{d\Lambda
d\Lambda^{*}}\left\langle\overleftarrow{T}\text{e}^{\int_{0}^{t}\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t^{\prime})dt^{\prime}}\right\rangle\right|_{\Lambda,\Lambda^{*}=0}=\left.\frac{d^{2}}{d\Lambda
d\Lambda^{*}}G(\Lambda,\Lambda^{*})\right|_{\Lambda,\Lambda^{*}=0}.$ (49)

with

$\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t^{\prime})=\mathcal{L}_{I}(t^{\prime})+\delta(t-t^{\prime})\left(\Lambda\mathcal{A}^{L}(t^{\prime})+\Lambda^{*}\mathcal{A}^{L\dagger}(t^{\prime})-\frac{|\Lambda|^{2}}{2}\right)$
(50)

To evaluate $G(\Lambda,\Lambda^{*})$ we now split the trace over the full
system in Eq. (49) into a partial trace over the bath, followed by a partial
trace over the system

$G(\Lambda,\Lambda^{*})=\left\langle\overleftarrow{T}\left\langle\overleftarrow{T}\text{e}^{\int_{0}^{t}\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t^{\prime})dt^{\prime}}\right\rangle_{B}\right\rangle_{S},$
(51)

where we have used idempotency of time ordering,
$\overleftarrow{T}=\overleftarrow{T}\overleftarrow{T}$. Since the operational
part of $\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t^{\prime})$ is linear
in bath operators we can use standard Gaussian integral results to evaluate
the trace over the bath analytically. This gives

$\left\langle\overleftarrow{T}\text{e}^{\int_{0}^{t}\overline{\mathcal{L}_{I}}(\Lambda,\Lambda^{*},t^{\prime})dt^{\prime}}\right\rangle_{B}\\\
=\exp\left(\Phi(t)+\Lambda^{*}\int_{0}^{t}dt^{\prime}\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{L}_{I}(t^{\prime})\right\rangle_{B}+\Lambda\int_{0}^{t}dt^{\prime}\left\langle\mathcal{A}^{L}(t)\mathcal{L}_{I}(t^{\prime})\right\rangle_{B}+|\Lambda|^{2}n(0)\right),$
(52)

where $\Phi(t)$ is the usual Feynman-Vernon influence phase (i.e. the exponent
in Eq. (44)) and $n(0)=\operatorname{Tr}[a^{\dagger}(0)a(0)\rho]$. Since this
is the only part of $G(\Lambda,\Lambda^{*})$ that carries dependence on
$\Lambda$ and $\Lambda^{*}$ we can now take the derivatives and set
$\Lambda=\Lambda^{*}=0$ to obtain

$n(t)\\\
=n(0)+\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}\left\langle\overleftarrow{T}\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{L}_{I}(t^{\prime})\right\rangle_{B}\left\langle\mathcal{A}^{L}(t)\mathcal{L}_{I}(t^{\prime\prime})\right\rangle_{B}\text{e}^{\Phi(t)}\right\rangle_{S}.$
(53)

To evaluate Eq. (53) we now use

$\mathcal{L}_{I}(t^{\prime})=-ig\mathcal{S}^{L}(t^{\prime})\left(\mathcal{A}^{L\dagger}(t^{\prime})+\mathcal{A}^{L}(t^{\prime})\right)+ig\mathcal{S}^{R}(t^{\prime})\left(\mathcal{A}^{R\dagger}(t^{\prime})+\mathcal{A}^{R}(t^{\prime})\right),$
(54)

which results in

$\displaystyle
n(t)=n(0)+g^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}$
$\displaystyle\Big{[}\left\langle\overleftarrow{T}\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{R}(t^{\prime\prime})\right\rangle_{H}\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{A}^{L}(t^{\prime})\right\rangle_{B}\left\langle\mathcal{A}^{L}(t)\mathcal{A}^{R\dagger}(t^{\prime\prime})\right\rangle_{B}$
$\displaystyle-\left\langle\overleftarrow{T}\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{L}(t^{\prime\prime})\right\rangle_{H}\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{A}^{L}(t^{\prime})\right\rangle_{B}\left\langle\mathcal{A}^{L}(t)\mathcal{A}^{L\dagger}(t^{\prime\prime})\right\rangle_{B}$
$\displaystyle+\left\langle\overleftarrow{T}\mathcal{S}^{R}(t^{\prime})\mathcal{S}^{L}(t^{\prime\prime})\right\rangle_{H}\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{A}^{R}(t^{\prime})\right\rangle_{B}\left\langle\mathcal{A}^{L}(t)\mathcal{A}^{L\dagger}(t^{\prime\prime})\right\rangle_{B}$
$\displaystyle-\left\langle\overleftarrow{T}\mathcal{S}^{R}(t^{\prime})\mathcal{S}^{R}(t^{\prime\prime})\right\rangle_{H}\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{A}^{R}(t^{\prime})\right\rangle_{B}\left\langle\mathcal{A}^{L}(t)\mathcal{A}^{R\dagger}(t^{\prime\prime})\right\rangle_{B}\Big{]}.$
(55)

where we have used the fact that

$\left\langle\mathcal{A}^{\alpha\dagger}(t)\mathcal{A}^{\alpha^{\prime}\dagger}(t^{\prime})\right\rangle_{B}=\left\langle\mathcal{A}^{\alpha}(t)\mathcal{A}^{\alpha^{\prime}}(t^{\prime})\right\rangle_{B}=0,$
(56)

for any $\alpha$ and $\alpha^{\prime}$. We have also defined

$\left\langle\overleftarrow{T}\mathcal{S}^{\alpha}(t^{\prime})\mathcal{S}^{\alpha^{\prime}}(t^{\prime\prime})\right\rangle_{H}=\left\langle\overleftarrow{T}\mathcal{S}^{\alpha}(t^{\prime})\mathcal{S}^{\alpha^{\prime}}(t^{\prime\prime})\text{e}^{\Phi(t)}\right\rangle_{S}=\left\langle\overleftarrow{T}\mathcal{S}^{\alpha}(t^{\prime})\mathcal{S}^{\alpha^{\prime}}(t^{\prime\prime})\text{e}^{\int_{0}^{t}\mathcal{L}_{I}(t^{\prime})dt^{\prime}}\right\rangle$
(57)

to indicate expectation values taken with respect to the full system density
matrix under full system evolution. To be clear, the subscript $H$ means
superoperators within the brackets are in the full system Heisenberg picture
and the trace is taken with respect to the full system density matrix, the
subscript $S$ means superoperators are in the interaction picture and the
trace is taken with respect to the reduced system density matrix, and the
brackets without subscript means superoperators are in the interaction picture
and the trace is taken with respect to the full system density matrix. For
instance, if $t^{\prime}>t^{\prime\prime}$ then

$\left\langle\overleftarrow{T}\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{L}(t^{\prime\prime})\right\rangle_{H}=\left\langle\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{L}(t^{\prime\prime})\right\rangle_{H}=\mathrm{Tr}[s(t^{\prime})s(t^{\prime\prime})\rho].$
(58)

Now we evaluate the bath expectation values in Eq. (B). For example, the in
the first line of Eq. (B) these evaluate to

$\left\langle\mathcal{A}^{L\dagger}(t)\mathcal{A}^{L}(t^{\prime})\right\rangle_{B}\left\langle\mathcal{A}^{L}(t)\mathcal{A}^{R\dagger}(t^{\prime\prime})\right\rangle_{B}=n(0)^{2}\text{e}^{-i\omega(t^{\prime}-t^{\prime\prime})},$
(59)

where we have used cyclicity of the trace,
$\left\langle\mathcal{A}^{L}(t)\mathcal{A}^{R\dagger}(t^{\prime\prime})\right\rangle_{B}=\left\langle\mathcal{A}^{L\dagger}(t^{\prime\prime})\mathcal{A}^{L}(t)\right\rangle_{B}$,
and made the interaction picture time evolution of bath superoperators
explicit, $\mathcal{A}^{\alpha}(t)=\mathcal{A}^{\alpha}\text{e}^{-i\omega t}$.
The overall result is

$\displaystyle
n(t)=n(0)+g^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}$
$\displaystyle\text{e}^{-i\omega(t^{\prime}-t^{\prime\prime})}\Big{[}\left\langle\overleftarrow{T}\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{R}(t^{\prime\prime})\right\rangle_{H}n(0)^{2}$
$\displaystyle-\left\langle\overleftarrow{T}\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{L}(t^{\prime\prime})\right\rangle_{H}n(0)(n(0)+1)$
$\displaystyle+\left\langle\overleftarrow{T}\mathcal{S}^{R}(t^{\prime})\mathcal{S}^{L}(t^{\prime\prime})\right\rangle_{H}(n(0)+1)^{2}$
$\displaystyle-\left\langle\overleftarrow{T}\mathcal{S}^{R}(t^{\prime})\mathcal{S}^{R}(t^{\prime\prime})\right\rangle_{H}n(0)(n(0)+1)\Big{]}.$
(60)

To help us enforce time-ordering upon system superoperators we split the
rectangular integral domain into two triangular domains,

$\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}=\int_{0}^{t}dt^{\prime}\int_{0}^{t^{\prime}}dt^{\prime\prime}+\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}dt^{\prime},$
(61)

where the first domain has $t^{\prime}>t^{\prime\prime}$ and the second had
$t^{\prime\prime}>t^{\prime}$. With this, and again using cyclicity of the
trace, e.g.
$\left\langle\mathcal{S}^{L}(t^{\prime})\mathcal{S}^{R}(t^{\prime\prime})\right\rangle=\left\langle\mathcal{S}^{L}(t^{\prime\prime})\mathcal{S}^{L}(t^{\prime})\right\rangle$,
we find that terms propotional to $n(0)^{2}$ cancel out and we are left with

$\displaystyle n(t)$
$\displaystyle=n(0)+g^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}\operatorname{Tr}[s(t^{\prime})s(t^{\prime\prime})\rho]\\{[n(0)+1]\text{e}^{-i\omega(t^{\prime}-t^{\prime\prime})}-n(0)\text{e}^{i\omega(t^{\prime}-t^{\prime\prime})}\\}$
(62)
$\displaystyle=n(0)+g^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}\operatorname{Tr}[s(t^{\prime})s(t^{\prime\prime})\rho](\cos(\omega(t^{\prime}-t^{\prime\prime}))-i\sin(\omega(t^{\prime}-t^{\prime\prime}))\coth(\omega/2T)),$
(63)

as used in the main text.

##  Appendix C Toy model

Figure 4:  Comparison of prediction by Eq. (62) (dots) with the exact solution
(lines) calculated by propagating the full system and environment representing
the bosons with a truncated basis of 4 levels each. The other parameters are
$\epsilon=0.1\Omega$, $g_{0}=0.1\Omega$, $g_{1}=0.2\Omega$,
$\omega_{0}=0.9\Omega$, $\omega_{1}=1.1\Omega$ and $T=0.1\Omega$.

In this section we benchmark our technique by calculating $\Delta Q_{q}$ for a
two-level system coupled to two bosonic modes. With so few degrees of freedom
we can simulate the dynamics of the modes explicitly and compare with the
prediction given by Eq. (9). The model we use is the two-mode Rabi model, the
Hamiltonian of which is given by

$\displaystyle H=\epsilon s_{z}+\Omega s_{x}$
$\displaystyle+s_{z}g_{0}(a_{0}+a_{0}^{\dagger})+\omega_{0}a_{0}^{\dagger}a_{0}$
$\displaystyle+s_{z}g_{1}(a_{1}+a_{1}^{\dagger})+\omega_{1}a_{1}^{\dagger}a_{1},$
(64)

and describes a two-level system driven with strength $\Omega$ and bias
$\epsilon$. The environment consists of bosons of frequency $\omega_{q}$ where
$a_{q}^{\dagger}$ creates an excitation in the mode indexed by
$q\in\\{0,1\\}$. The system is coupled to each of these modes with strength
$g_{q}$.

In Fig. 4 we plot the dynamics of each mode, calculated first by numerically
integrating the Schrödinger equation for the full system-environment dynamics,
and then from Eq. (62). We initialised the two-level system in the spin-up
state and each bosonic mode in a thermal state at temperature $T=0.1\Omega$.
It is clear that our new method correctly predicts the environment dynamics
for both modes.