# High-resolution [OI] line spectral mapping of TW Hya consistent with X-ray
driven photoevaporation.

Christian Rab University Observatory, Faculty of Physics, Ludwig-Maximilians-
Universität München, Scheinerstr. 1, 81679 Munich, Germany  Max-Planck-
Institut für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching,
Germany  Michael Weber University Observatory, Faculty of Physics, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 Munich, Germany
Exzellenzcluster ‘Origins’, Boltzmannstr. 2, 85748 Garching, Germany  Giovanni
Picogna University Observatory, Faculty of Physics, Ludwig-Maximilians-
Universität München, Scheinerstr. 1, 81679 Munich, Germany  Barbara Ercolano
University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität
München, Scheinerstr. 1, 81679 Munich, Germany  Exzellenzcluster ‘Origins’,
Boltzmannstr. 2, 85748 Garching, Germany  James Owen Astrophysics Group,
Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7
2AZ, UK

(Received August 11, 2023; Revised August 28, 2023; Accepted August 31, 2023)

###### Abstract

Theoretical models indicate that photoevaporative and magnetothermal winds
play a crucial role in the evolution and dispersal of protoplanetary disks and
affect the formation of planetary systems. However, it is still unclear what
wind-driving mechanism is dominant or if both are at work, perhaps at
different stages of disk evolution. Recent spatially resolved observations by
Fang et al. (2023) of the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$spectral
line, a common disk wind tracer, in TW Hya revealed that about 80% of the
emission is confined to the inner few au of the disk. In this work, we show
that state-of-the-art X-ray driven photoevaporation models can reproduce the
compact emission and the line profile of the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line. Furthermore, we show that
the models also simultaneously reproduce the observed line luminosities and
detailed spectral profiles of both the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and the
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$lines. While MHD wind models can
also reproduce the compact radial emission of the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line they fail to match the
observed spectral profile of the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line and underestimate the
luminosity of the $\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$line by a factor
of three. We conclude that, while we cannot exclude the presence of an MHD
wind component, the bulk of the wind structure of TW Hya is predominantly
shaped by a photoevaporative flow.

††journal: ApJL††software: This research made use of Astropy, a community-
developed core Python package for Astronomy (Astropy Collaboration et al.,
2013, 2018). This research made use of Photutils, an Astropy package for
detection and photometry of astronomical sources (Bradley et al., 2022).
matplotlib version 3.7.0 (Hunter, 2007); scipy (Virtanen et al., 2020); numpy
(Harris et al., 2020).

##  1 Introduction

The evolution and final dispersal of protoplanetary disks is thought to
strongly affect the formation and evolution of planetary systems. Disk winds
are considered to be significant contributors to the evolutionary processes
occurring within protoplanetary disks (Lesur et al., 2022; Pascucci et al.,
2022). Thermal winds can be launched through photoevaporation (PE) from the
central star (e.g. Gorti & Hollenbach, 2009; Nakatani et al., 2018; Picogna et
al., 2021; Ercolano et al., 2021) and are efficient at removing material at
rates comparable to the observed accretion rates of T-Tauri stars (e.g.
Ercolano & Pascucci, 2017). Thermal winds do not remove angular momentum from
the disk, and, when combined with viscous accretion models, they are
successful in reproducing the observed two-timescale behavior, evidenced by
the evolution of disk colors (e.g. Koepferl et al., 2013; Ercolano et al.,
2015) and several observational correlations such as the observed accretion
rates and the mass of the central star (Ercolano et al., 2014) or inner disk
life times (Picogna et al., 2021).

The inclusion of non-ideal magneto-hydrodynamical effects in weakly ionized
protoplanetary disk has shown that magnetorotational instability (MRI) (Balbus
& Hawley, 1991), hypothesized to drive viscosity in disks, is suppressed in
most regions of the disk (see e.g. Lesur et al., 2022, for a recent review),
except the very inner regions where thermionic emission from dust dominates
(Desch & Turner, 2015; Jankovic et al., 2021). Vigorous, magnetically
supported, disk winds (from now on MHD winds), are a solid prediction of most
recent simulations (e.g. Gressel et al., 2015; Bai et al., 2016; Wang et al.,
2019; Lesur, 2021; Gressel et al., 2020) and they replace MRI in most disk
regions by removing angular momentum from the disk, allowing material to flow
inward.

Which type of wind might dominate at different times and different locations
in a disk is an important question, which directly affects planet formation
models. The current picture emerging from the careful analysis of
spectroscopic diagnostics is that both types of winds operate in disks, with
MHD winds stronger in young objects and thermal winds dominating the final
evolution and eventual dispersal of disks (Ercolano & Pascucci, 2017; Weber et
al., 2020).

Currently used spectroscopic wind diagnostics, particularly the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$collisionally excited spectral
line, have complex line profiles (e.g. Simon et al., 2016; Fang et al., 2018;
Banzatti et al., 2019; Gangi et al., 2020), often preventing important wind
parameters, like the wind launching radius, to be directly determined (see
discussion in Weber et al. 2020). Rab et al. (2022) find that a combination of
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and molecular hydrogen
observations are consistent with thermal winds driven by X-ray
photoevaporation, but alternative models cannot be ruled out.

In order to break the degeneracies hidden in the interpretation of non-
spatially resolved line profiles, high-resolution spectral mapping of wind
diagnostics represent an attractive option. This has recently been done for
the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line from TW Hya by Fang et
al. (2023), using the multi-unit spectroscopic explorer (MUSE) at the Very
Large Telescope, who showed that about 80% of the [OI] emission is confined to
within 1 au radially from the star. In this paper we show that state-of-the-
art thermal wind models driven by X-ray photoevaporation (e.g. Picogna et al.,
2019; Ercolano et al., 2021; Picogna et al., 2021) are consistent with the
observations recently published by Fang et al. (2023). In Section 2 we briefly
describe the used photoevaporative disk wind models and our approach to
produce synthetic observables. In Section 3 we show our results in particular
the comparison to the observational data. We discuss our results in context to
previous works and MHD disk wind models and present our conclusions in Section
4.

##  2 Methods

In this section, we describe the physical models used in this work and how we
produce synthetic observables from those models that can be directly compared
to observational data.

###  2.1 Photoevaporative disks wind models

To model a photoevaporative disk wind we follow the approach by Picogna et al.
(2019, 2021); Ercolano et al. (2021)111X-ray PE models and data from Picogna
et al. (2021) are available here https://cutt.ly/lElY9JI, to which we refer
for details. This model uses a modified version of the PLUTO code (Mignone et
al., 2007; Picogna et al., 2019) to perform radiative-hydrodynamic simulations
of a disk irradiated by a central star. The temperatures in the wind and the
wind-launching regions, i.e. the upper layers of the disk, where the column
number density towards the central star is in the range between $5\times
10^{20}$ and $2.5\times 10^{22}$ cm-2, are determined by parameterizations
that are derived from detailed radiative transfer calculations with the gas
photo-ionization code MOCASSIN (Ercolano et al., 2003, 2005, 2008). For a
given column number density the respective parameterization yields the gas
temperature dependent on the ionization parameter $\xi=\frac{L_{X}}{nr^{2}}$,
where $L_{\mathrm{X}}$ is the X-ray luminosity of the star, $n$ the volume
number density and $r$ the spherical radius. In this work, we use a stellar
mass $M_{*}=0.7M_{\odot}$ and the parametrizations derived from the spectrum
labeled as Spec29 in Ercolano et al. (2021) with $L_{X}=2\times 10^{30}$ erg
s-1, which is appropriate considering observational constraints on TW Hya
(Robrade & Schmitt, 2006; Ercolano et al., 2017). The computational grid was
centered on the star. Spherical polar coordinates were adopted with 512
logarithmic spaced cells in the radial direction from $0.33$ au to $600$ au,
and 512 uniform spaced cells in the polar one from $0.005$ to $\pi/2$. Outflow
boundaries were adopted in the radial directions, while special reflective
boundaries were used to treat the regions close to the polar axis and the disk
mid-plane. A periodic boundary was assumed in the azimuthal direction. The
influence of the inner boundary was tested by decreasing the inner radial
boundary to $0.1$ au, while keeping the same radial resolution outside $0.33$
au.

###  2.2 Disk model without a wind

Additionally to the PE disk wind models, we use an existing radiation thermo-
chemical disk model for TW Hya from the DIANA (DIsc
ANAlysis222https://diana.iwf.oeaw.ac.at) project presented in Woitke et al.
(2019). This model does not include a wind component but was made to reproduce
existing (mostly spatially unresolved) observational data including the
spectral energy distribution and about 50 spectral lines (i.e. line fluxes are
matched within a factor of two to three). With this model we show how a pure
disk model compares with the spatially resolved
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$observations. As at the time of
the publication of this model no spatially resolved observables for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$were produced, we rerun the model
with a more recent version of the radiation thermo-chemical code PRODIMO
(PROtoplanetary DIsk MOdel333https://prodimo.iwf.oeaw.ac.at revision: 66efbd75
2023/06/27, Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Woitke et
al. 2016) to produce line cubes and images that can be compared to the
spatially resolved data.

###  2.3 Synthetic observables

To produce synthetic observables we use two different approaches to post-
process the PE disk wind model. Similar to Weber et al. (2020) we use the
MOCASSIN Monte Carlo radiative transfer code that allows to model spectral
line emission for atomic species in the optical and infrared (Ercolano et al.,
2003, 2005, 2008). Furthermore we use the radiation thermo-chemical code
PRODIMO that was recently applied in Rab et al. (2022) to produce synthetic
observables for atomic and molecular species that are supposed to trace disk
winds. We use both approaches to show that our results, in particular the
spatial extent of the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$, are robust
and do not strongly depend on the details of the post-processing method (e.g.
chemistry, heating-cooling, line excitation).

For both approaches, we use the same physical structure (density and velocity
field), the same X-ray/stellar spectrum, and the same dust properties. To
account for the accretion luminosity we add to the X-ray spectrum a black body
spectrum with $T_{\mathrm{eff}}=12000\,\mathrm{K}$ normalized to
$L_{\mathrm{acc}}=2.95\times 10^{-2}~{}L_{\odot}$ (Fang et al., 2018). For the
dust we assume a gas-to-dust ratio of 100 and an ISM size distribution (see
Weber et al. 2020 for details). Quantities such as temperatures, line
populations, chemical abundances and the synthetic observables are self-
consistently calculated within each post-processing framework. For a more
detailed discussion on these different post-processing approaches and a
comparison see Rab et al. (2022).

At first we use the line radiative transfer modules of MOCASSIN and PRODIMO to
produce synthetic model images for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$emission assuming a distance of
$60\,\mathrm{pc}$ (Gaia Collaboration et al., 2021) and a disk inclination of
$7^{\circ}$ (Qi et al., 2004), consistent with Fang et al. (2023). For the
produced model images we use an oversampling factor of seven
($0.003615\arcsec$) compared to the pixel scale of the observations. Following
Fang et al. (2023) we downsample the model images to the pixel scale of the
MUSE data before convolving it with the PSF. For these model images and for
the observed image we produce azimuthally averaged radial profiles using the
corresponding pixel scale as the width for the radial bins. Additionally we
also produce radial profiles for the unconvolved model images to better show
the real extend of the emission. We also produce synthetic images in the same
way for the toy model (power-law model) presented in Fang et al. (2023) and
the thermo-chemical model without a wind. We note that the used oversampling
factor is not enough to fully resolve the toy model of Fang et al. (2023) but
we found that it is sufficient to reproduce the observation after convolution
of the toy model image (mainly because a downsampling to the pixel size of the
observations is anyway required). We also note that a higher resolution
spatial grid is used for the radiative transfer step, but as we are not aiming
for a detailed fitting of the observational data we use the oversampling
factor of all model images for consistency (e.g. between the two different
post-processing methods) and efficiency.

We also compare our model with the observations of the
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$line presented in Pascucci et al.
(2011). For this line we simply produce spectral line profiles again using
both post-processing approaches.

##  3 Results

Here, we compare our modeling results to the observational data. As we do not
present newly developed PE models we focus only on the observables. In Sect.
3.1 we present normalized azimuthally averaged radial profiles for
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$in a similar fashion as Fang et
al. (2023) and in Sect. 3.2 we compare our PE wind models to the observed
spectral profiles for the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$spectral lines.

###  3.1 Radial profiles

Figure 1: The toy (power-law) model (green) from Fang et al. (2023) compared
to results from a radiation thermo-chemical disk model (purple) without a wind
that fits a large set of observational data for TW Hya (Woitke et al., 2019).
The dashed and solid colored lines show the radial profiles for the
unconvolved and convolved images, respectively. The gray solid line shows the
radial profile of the PSF and the black solid line for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$MUSE observations (Fang et al.,
2023). The inset shows a zoom-in to the region marked by the brown box
($r=1.3-6.4\,\mathrm{au}$) in the main panel.

####  3.1.1 Disk only model

In Fig. 1, we show the toy model from Fang et al. (2023) in comparison to the
disk-only model from Woitke et al. (2019). The thermo-chemical disk model
shows a similar radial profile for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$emission, in particular the steep
slope, but is slightly more compact and hence does not match the data (i.e. it
remains unresolved). We note that the model of Woitke et al. (2019) uses a
parameterized disk structure, which is quite different to the disk wind models
used in this work or in Fang et al. (2023), in particular, it has a dust and
gas depleted (optically thin but not empty) inner hole extending up to almost
$3\,\mathrm{au}$. Although it is possible to adapt this model to achieve a
better match to the spatially resolved data, such a disk only model, by
construction, cannot match the blue-shifts seen in the observed line profiles
of the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$lines. Nevertheless, it is
interesting to see that such a model is almost in agreement with the spatially
resolved data, although such constraints were not included in the modeling.
Furthermore, this model indicates that a disk only solution seems to produce
an even more compact emission region for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$compared to the wind models
presented here or in Fang et al. (2023) and also that the contribution from a
disk itself can be significant in the case of TW Hya.

Figure 2: Comparison of the photoevaporative disk wind model with the observed
normalized radial intensity profiles of the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line. The black solid line shows
the observation, the gray solid line the PSF (both from Fang et al. 2023. The
colored solid and dashed lines show the model results, convolved with the PSF
and at the resolution of the model image, respectively. In red, we show the
results from the MOCASSIN and in blue from the PRODIMO post-processing. The
inset shows a zoom-in to the region marked by the brown box
($r=1.3-6.4\,\mathrm{au}$) in the main panel. Figure 3:
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$(top) and
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$(bottom) line flux of our
photoevaporation disk wind model obtained by multiplying the emissivity with
$2\pi R$ in order to illustrate the origin of the emission. White solid
contour lines indicate the regions where the top 80 % of emission originates.
White arrows represent the velocity field in the PE model.

####  3.1.2 Photoevaporative disk wind models

In Fig. 2 we compare the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$radial
intensity profiles of our PE wind model at an inclination of $7^{\circ}$ to
the observation by Fang et al. (2023). As can be seen in the unconvolved
radial profiles, both, the MOCASSIN and PRODIMO models yield very similar
results, with the PRODIMO model having slightly more emission at very low
($\lesssim 0.4$ au) and at extended radii, and MOCASSIN showing enhanced
emission at intermediate radii between $\approx 0.5$ and 3 au. In both models,
the emission peaks well inside of 1 au with a steep decrease in intensity at
larger radii. Computing the cumulative integral, we find that 80% of the
emission originates inside 2 au of the central star, compared to 1 au for the
toy model that Fang et al. (2023) derived as a fit to the data. This is also
visible in Fig. 3, where we show 2D emission maps overlain by contours showing
the 80% regions. It is worth pointing out that although the 80% regions extend
to $\approx 2$ au and $\approx 14$ au for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$lines, respectively, the emission
inside this region is not uniform but has a strong gradient with the peak
close to the star. Comparing the profiles after convolution with the PSF, it
can be seen that the radial profiles of the PE wind model are in very good
agreement with the observations. This shows that current state-of-the-art
photoevaporative disk models produce compact emission consistent with the
spatially resolved observations of TW Hya.

Figure 3 shows significant emission originating from the inner 1 au, which is
inside the gravitational radius of the X-ray photoevaporation models shown
here ($\sim 3.5$ au). The emitting material close to the star is thus bound to
the inner disk and not affected by photoevaporation. This implies that our
main result would not qualitatively change in the case of a ”gapped” inner
disk as suggested by the observations of a dark annulus in mm-wave dust
emission at $\approx$ 1 au (Andrews et al., 2016). Indeed, the presence of gas
close to the central star of TW Hya, which shows clear accretion signatures,
justifies the employment of primordial disk models in this case. While we do
not expect significant changes to the wind structure due to the gap and hence
on the main picture of having compact emission, a detailed hydrodynamical
model would be required to assess potential differences in the detailed
spectroscopic line profiles emitted from this region (see also Sect. 4). Such
a model is, however, outside the scope of this study and would also require
higher spatial resolution in the observations of the inner disc, which is not
possible with current instrumentation.

Figure 4: Comparison of observed
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$(Simon et al., 2016; Fang et al.,
2018, 2023) and $\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$(Pascucci et al.,
2011) spectral line profiles to synthetic line profiles of photoevaporative
disk wind models. The [OI] line profiles from the models and that observed by
Fang et al. (2023) have been degraded to a resolving power $R=40000$,
comparable to the resolution in the observation by Simon et al. (2016). For
the Fang et al. (2023) data the error bars correspond to $3\sigma$
uncertainties (for visibility). The Gaussian distribution shown in green
represents the derived mean centroid shift of $-0.8\,\mathrm{kms^{-1}}$ with a
standard deviation of $0.4\,\mathrm{kms^{-1}}$ as derived by Fang et al.
(2023) from various observations of
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$. The [NeII] model profile has
been degraded to $R=30000$ and the error bars in the observed profile indicate
$1\sigma$ uncertainties.

###  3.2 Spectral line profiles

In Fig. 4 we compare the modelled spectral line profiles to the observations
of the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$spectral lines. For
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$we show three different observed
spectra representing the scatter observed in the measured centroid velocities
$v_{\mathrm{c}}$ with a mean of $v_{\mathrm{c}}=-0.8\,\mathrm{kms^{-1}}$ and a
standard deviation of $0.4\,\mathrm{kms^{-1}}$ (Fang et al., 2023). Both
models are in good agreement with the observations. In particular the MOCASSIN
($v_{c}\approx-0.8\,\mathrm{kms^{-1}}$) model matches the profile of (Simon et
al., 2016), that shows $v_{c}\approx-0.8\,\mathrm{kms^{-1}}$, exceptionally
well. The PRODIMO profile is slightly broader in the blue part of the spectrum
and therefore appears more blue-shifted
($v_{c}\approx-1.5\,\mathrm{kms^{-1}}$). Nevertheless this is still consistent
with observations as for example the profile of Fang et al. (2023) shows a
similar behavior with $v_{c}\approx-1.5\,\mathrm{kms^{-1}}$, for the spectrum
downgraded to R=40000. The $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line
fluxes from the models are $1.2\times 10^{-5}\,L_{\sun}$ (MOCASSIN) and
$5.5\times 10^{-6}\,L_{\sun}$ (PRODIMO ), in good agreement with the observed
values of $1.0-1.5\times 10^{-5}\,L_{\sun}$ (Simon et al., 2016; Fang et al.,
2018, 2023).

The difference in the shape of the two model spectra can be explained by the
slightly more extended $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$emission of
the PRODIMO model (see Fig.2), which traces slightly faster velocities of the
PE wind compared to the MOCASSIN model. We tested this by simply removing all
emission for $r>4\,\mathrm{au}$ in the synthetic observables of the PRODIMO
model and find that for this case the MOCASSIN and PRODIMO line profiles
become almost identical. As the density structure and velocity fields in both
models are the same; the differences arise from different radial temperature
gradients and differences in the line excitation calculations. However, our
results indicate that detailed models for TW Hya, fully considering the
stellar properties and possibly also details of the disk structure (i.e. a gap
in the disk structure; see Owen 2011) are required for a comprehensive
interpretation of the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line
profile.

For $\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$the modeled profiles are very
similar and match very well the observed spectral profile. As
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$traces regions further out and
higher up (up to r $\approx 10\,\mathrm{au})$ in the disk/wind with respect to
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$, it traces higher velocities of
the photoevaporative flow (see Fig. 3), consistent with the observed
$v_{\mathrm{c}}\approx-5\,\mathrm{kms^{-1}}$. The line luminosity of the
models is $3.4\times 10^{-6}\,L_{\sun}$ (MOCASSIN) and $5.1\times
10^{-6}\,L_{\sun}$ (PRODIMO ) which are in excellent agreement with the
observed range of luminosities of $\approx\\!3.5-6.2\times 10^{-6}\,L_{\sun}$
(Pascucci et al., 2011; Pascucci & Sterzik, 2009; Najita et al., 2010).

##  4 Discussion and Conclusions

Our results show that current state-of-the-art PE disk wind models are
consistent with the observational data presented in Fang et al. (2023). We
want to emphasize that for this work, we did not use any new developments or
adaptations to the PE wind modeling approach already presented and used in
Picogna et al. (2019); Weber et al. (2020); Picogna et al. (2021); Ercolano et
al. (2021); Rab et al. (2022); Weber et al. (2022). The only difference
between this work and previously published work is the choice of an
appropriate X-ray spectrum for TW Hya.

Fang et al. (2023) compared the spatially resolved
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$emission to very low resolution
early photoevaporation models (Owen et al., 2010; Ercolano & Owen, 2010) which
present a more extended emission profile for the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$. Based on this comparison, they
concluded that a magneto-thermal driven wind is necessary to explain the
spatially resolved $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line data. We
show here, however, that, compared to the fiducial magneto-thermal wind model
presented in Fang et al. (2023), modern PE models produce only very slightly
more extend $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$emission but are still
fully consistent with the data. Apart from the higher inner grid resolution,
an important difference between the ”old” PE models (Owen et al., 2010;
Ercolano & Owen, 2010) and the newer PE models (Picogna et al., 2019; Weber et
al., 2020; Picogna et al., 2021; Ercolano et al., 2021) is the temperature
parameterization. The new models take into account the detailed column density
distribution (attenuation) in the simulation regions, yielding a more accurate
temperature (and thus density) profile. More specifically, the density in the
inner disk is higher for the new models, resulting in a more compact emission
region.

The new PE models can match the observed line luminosity and the shape of the
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$profile very well, in particular
the small observed blue-shift of $v_{c}\approx-0.8\,\mathrm{km/s}$. As noted
by Fang et al. (2023) the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$spectral
profile of their MHD model is too blue-shifted compared to observations, which
is likely caused by the higher wind velocities in the inner regions compared
to our PE model. Fang et al. (2023) argue that the discrepancy could be caused
by the lack of an inner hole in their primordial disk model and that the
presence of such a hole would allow for red-shifted emission from the back
side of the disk to contribute to the line profile, reducing its blue-shift.
However, in order to zero the velocity center by this effect, the line would
then result broadened by the blue-shifted value (Ercolano & Owen, 2010, 2016),
which would then again be in tension with the observations. Additionally, as
TW Hya is seen almost face-on, the contribution from the back-side of the disk
will be limited unless an unsuitable large inner hole (several au) is assumed.
More quantitatively speaking, our model has an inner radius of
$r_{\mathrm{in}}=0.33\,\mathrm{au}$ (i.e. larger than the MHD model), but
still the back side of the disk contributes to less than 2% to the total flux
and does not have any significant impact on the centroid velocity. We also
tested a PE wind model with $r_{\mathrm{in}}=0.1\,\mathrm{au}$ and found no
significant differences. In any case, for a more thorough interpretation of
the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$spectral profile of TW Hya, a
more detailed disk structure model for the inner few au is required. In
particular constraints from ALMA observations (gap at
$\approx\,1\,\mathrm{au}$ and unresolved emission from $r<0.5\,\mathrm{au}$;
Andrews et al. 2016$)$ and VLT/SPHERE (marginal detection of emission from the
inner few au; van Boekel et al. 2017) indicate that there is likely still some
dust in the inner $1-2\,\mathrm{au}$ (see also Ercolano et al., 2017) that
would reduce $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$emission from the
backside of the disk.

As already noted in Pascucci et al. (2011) and also discussed in Fang et al.
(2023) the $\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$line observation points
towards a thermally driven wind. This is supported by the models presented in
this work, as the agreement of the PE wind model with the spectral line
profile and the observed line fluxes (which is under-predicted by a factor of
three by the MHD model of Fang et al. 2023) is excellent, indicating that at
least for the radii larger than a few au, traced by the
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$line, the disk wind structure of TW
Hya is predominantly shaped by a PE flow.

We conclude that, while the currently available spatially resolved data does
not allow to clearly distinguish a pure thermally driven wind from a magneto-
thermal wind, the simultaneous agreement of the PE wind models with the
spatially resolved data and the spectral profiles of
$\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$and
$\mathrm{[NeII]}\,12.8\,\mathrm{\micron}\,$strongly indicates that at least
large parts of the disk wind seen in TW Hya are driven by X-ray
photoevaporation.

Acknowledgements: We thank the anonymous referee for a quick and constructive
report. We thank Min Fang for providing their reduced observational data for
the $\mathrm{[OI]}\,6300\,\mathrm{\mbox{\AA}}\,$line including the PSF. We
acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) Research Unit “Transition discs” - 325594231. This
research was supported by the Excellence Cluster ORIGINS which is funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy - EXC-2094 - 390783311. CHR is grateful for
support from the Max Planck Society. This research has made use of NASA’s
Astrophysics Data System.

## References

  * Andrews et al. (2016) Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJ, 820, L40, doi: 10.3847/2041-8205/820/2/L40
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152, doi: 10.3847/0004-637X/818/2/152
  * Balbus & Hawley (1991) Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214, doi: 10.1086/170270
  * Banzatti et al. (2019) Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76, doi: 10.3847/1538-4357/aaf1aa
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.5.0, 1.5.0, Zenodo, doi: 10.5281/zenodo.6825092
  * Desch & Turner (2015) Desch, S. J., & Turner, N. J. 2015, ApJ, 811, 156, doi: 10.1088/0004-637X/811/2/156
  * Ercolano et al. (2005) Ercolano, B., Barlow, M. J., & Storey, P. J. 2005, MNRAS, 362, 1038, doi: 10.1111/j.1365-2966.2005.09381.x
  * Ercolano et al. (2003) Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X.-W. 2003, MNRAS, 340, 1136, doi: 10.1046/j.1365-8711.2003.06371.x
  * Ercolano et al. (2015) Ercolano, B., Koepferl, C., Owen, J., & Robitaille, T. 2015, MNRAS, 452, 3689, doi: 10.1093/mnras/stv1528
  * Ercolano et al. (2014) Ercolano, B., Mayr, D., Owen, J. E., Rosotti, G., & Manara, C. F. 2014, MNRAS, 439, 256, doi: 10.1093/mnras/stt2405
  * Ercolano & Owen (2010) Ercolano, B., & Owen, J. E. 2010, MNRAS, 406, 1553, doi: 10.1111/j.1365-2966.2010.16798.x
  * Ercolano & Owen (2016) —. 2016, MNRAS, 460, 3472, doi: 10.1093/mnras/stw1179
  * Ercolano & Pascucci (2017) Ercolano, B., & Pascucci, I. 2017, Royal Society Open Science, 4, 170114, doi: 10.1098/rsos.170114
  * Ercolano et al. (2021) Ercolano, B., Picogna, G., Monsch, K., Drake, J. J., & Preibisch, T. 2021, MNRAS, 508, 1675, doi: 10.1093/mnras/stab2590
  * Ercolano et al. (2017) Ercolano, B., Rosotti, G. P., Picogna, G., & Testi, L. 2017, MNRAS, 464, L95, doi: 10.1093/mnrasl/slw188
  * Ercolano et al. (2008) Ercolano, B., Young, P. R., Drake, J. J., & Raymond, J. C. 2008, ApJS, 175, 534, doi: 10.1086/524378
  * Fang et al. (2018) Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28, doi: 10.3847/1538-4357/aae780
  * Fang et al. (2023) Fang, M., Wang, L., Herczeg, G. J., et al. 2023, Nature Astronomy, doi: 10.1038/s41550-023-02004-x
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Smart, R. L., Sarro, L. M., et al. 2021, A&A, 649, A6, doi: 10.1051/0004-6361/202039498
  * Gangi et al. (2020) Gangi, M., Nisini, B., Antoniucci, S., et al. 2020, A&A, 643, A32, doi: 10.1051/0004-6361/202038534
  * Gorti & Hollenbach (2009) Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539, doi: 10.1088/0004-637X/690/2/1539
  * Gressel et al. (2020) Gressel, O., Ramsey, J. P., Brinch, C., et al. 2020, ApJ, 896, 126, doi: 10.3847/1538-4357/ab91b7
  * Gressel et al. (2015) Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84, doi: 10.1088/0004-637X/801/2/84
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hunter (2007) Hunter, J. D. 2007, Computing In Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jankovic et al. (2021) Jankovic, M. R., Owen, J. E., Mohanty, S., & Tan, J. C. 2021, MNRAS, 504, 280, doi: 10.1093/mnras/stab920
  * Kamp et al. (2010) Kamp, I., Tilling, I., Woitke, P., Thi, W.-F., & Hogerheijde, M. 2010, A&A, 510, A18, doi: 10.1051/0004-6361/200913076
  * Koepferl et al. (2013) Koepferl, C. M., Ercolano, B., Dale, J., et al. 2013, MNRAS, 428, 3327, doi: 10.1093/mnras/sts276
  * Lesur et al. (2022) Lesur, G., Ercolano, B., Flock, M., et al. 2022, arXiv e-prints, arXiv:2203.09821, doi: 10.48550/arXiv.2203.09821
  * Lesur (2021) Lesur, G. R. J. 2021, A&A, 650, A35, doi: 10.1051/0004-6361/202040109
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, The Astrophysical Journal Supplement Series, 170, 228, doi: 10.1086/513316
  * Najita et al. (2010) Najita, J. R., Carr, J. S., Strom, S. E., et al. 2010, ApJ, 712, 274, doi: 10.1088/0004-637X/712/1/274
  * Nakatani et al. (2018) Nakatani, R., Hosokawa, T., Yoshida, N., Nomura, H., & Kuiper, R. 2018, ApJ, 857, 57, doi: 10.3847/1538-4357/aab70b
  * Owen (2011) Owen, J. 2011, PhD thesis, University of Cambridge, UK
  * Owen et al. (2010) Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415, doi: 10.1111/j.1365-2966.2009.15771.x
  * Pascucci et al. (2022) Pascucci, I., Cabrit, S., Edwards, S., et al. 2022, arXiv e-prints, arXiv:2203.10068, doi: 10.48550/arXiv.2203.10068
  * Pascucci & Sterzik (2009) Pascucci, I., & Sterzik, M. 2009, ApJ, 702, 724, doi: 10.1088/0004-637X/702/1/724
  * Pascucci et al. (2011) Pascucci, I., Sterzik, M., Alexander, R. D., et al. 2011, ApJ, 736, 13, doi: 10.1088/0004-637X/736/1/13
  * Picogna et al. (2021) Picogna, G., Ercolano, B., & Espaillat, C. C. 2021, MNRAS, 508, 3611, doi: 10.1093/mnras/stab2883
  * Picogna et al. (2019) Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691, doi: 10.1093/mnras/stz1166
  * Qi et al. (2004) Qi, C., Ho, P. T. P., Wilner, D. J., et al. 2004, ApJ, 616, L11, doi: 10.1086/421063
  * Rab et al. (2022) Rab, C., Weber, M., Grassi, T., et al. 2022, A&A, 668, A154, doi: 10.1051/0004-6361/202244362
  * Robrade & Schmitt (2006) Robrade, J., & Schmitt, J. H. M. M. 2006, A&A, 449, 737, doi: 10.1051/0004-6361:20054247
  * Simon et al. (2016) Simon, M. N., Pascucci, I., Edwards, S., et al. 2016, ApJ, 831, 169, doi: 10.3847/0004-637X/831/2/169
  * Thi et al. (2011) Thi, W.-F., Woitke, P., & Kamp, I. 2011, MNRAS, 412, 711, doi: 10.1111/j.1365-2966.2010.17741.x
  * van Boekel et al. (2017) van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132, doi: 10.3847/1538-4357/aa5d68
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wang et al. (2019) Wang, L., Bai, X.-N., & Goodman, J. 2019, ApJ, 874, 90, doi: 10.3847/1538-4357/ab06fd
  * Weber et al. (2020) Weber, M. L., Ercolano, B., Picogna, G., Hartmann, L., & Rodenkirch, P. J. 2020, MNRAS, 496, 223, doi: 10.1093/mnras/staa1549
  * Weber et al. (2022) Weber, M. L., Ercolano, B., Picogna, G., & Rab, C. 2022, MNRAS, 517, 3598, doi: 10.1093/mnras/stac2954
  * Woitke et al. (2009) Woitke, P., Kamp, I., & Thi, W.-F. 2009, A&A, 501, 383, doi: 10.1051/0004-6361/200911821
  * Woitke et al. (2016) Woitke, P., Min, M., Pinte, C., et al. 2016, A&A, 586, A103, doi: 10.1051/0004-6361/201526538
  * Woitke et al. (2019) Woitke, P., Kamp, I., Antonellini, S., et al. 2019, PASP, 131, 064301, doi: 10.1088/1538-3873/aaf4e5


# Thermalization of open quantum systems using multiple-Davydov $\text{D}_{2}$
variational approach

Mantas Jakučionis, Darius Abramavičius Institute of Chemical Physics, Vilnius
University, Sauletekio Ave. 9-III, LT-10222 Vilnius, Lithuania

###### Abstract

Numerical implementations of explicit phonon bath require a large number of
oscillator modes in order to maintain oscillators at the initial temperature
when modelling energy relaxation processes. An additional thermalization
algorithm may be useful in controlling the local temperature. In this paper we
extend our previously proposed thermalization algorithm [1] to be used with
the numerically exact multiple-Davydov $\text{D}_{2}$ trial wavefunction for
simulation of relaxation dynamics and spectroscopic signals of open quantum
systems using the time-dependent Dirac-Frenkel variational principle. By
applying it to the molecular aggregate model, we demonstrate how the
thermalization approach significantly reduces the numerical cost of
simulations by decreasing the number of oscillators needed to explicitly
simulate the aggregate’s environment fluctuations, while maintaining
correspondence to the exact population relaxation dynamics. Additionally, we
show how the thermalization can be used to find the equilibrium state of the
excited molecular aggregate, which is necessary for simulation of the
fluorescence and other spectroscopic signals. The presented thermalization
algorithm open possibilities to investigate larger system-bath models than was
previously possible using the multiple-Davydov $\text{D}_{2}$ trial
wavefunction and local heating effects in molecular complexes.

##  I Introduction

Open quantum system models are widely used to describe properties of molecular
aggregates [2, 3]. The system usually consists of molecular electronic states.
Intramolecular vibrational degrees of freedom (DOFs), which play the major
role in relaxation process of the systems of interest can also be included
into the quantum system model. The rest of the DOFs are treated as an
environment of a constant temperature – the bath. The bath is modelled as a
collection of quantum harmonic oscillators (QHO) and is characterized by a
continuous fluctuation spectral density function [3, 4, 5, 6]. Separation into
the system and the bath parts is mostly formal as the system-bath coupling has
to be included to account for molecular environment-induced decoherence and
temperature effects, hence the quantum dynamics penetrates into the bath, and
bath also changes its state.

When using wavefunction-based simulation approaches, it can be challenging to
maintain a precise representation of the bath as a constant temperature
thermostat, because energy exchange between the system and the bath can alter
thermal properties of the bath. Generally, a large number of explicitly
modelled QHO modes have to be included to minimize the negative effects of
thermal energy accumulation in the bath, but this is numerically expensive.
Therefore, one always has to balance between the size of the model, accuracy
of the chosen numerical method and its numerical cost. Alternatively, one
could numerically change the wavefunction variables during its time evolution
in a way as to prevent accumulation of the thermal energy in the bath and to
maintain it at a desired temperature, i.e., perform thermalization.

It is challenging to accurately simulate dynamics of quantum systems that
exchange energy and (quasi-) particles with their surroundings, i.e., of the
open quantum systems [7, 8], because a numerical cost needed to propagate the
corresponding dynamical equations in time increases exponentially with the
number of DOFs. The wavefunction approach based on the multiple-Davydov
$\text{D}_{2}$ trial wavefunction ($\text{mD}_{2}$ ansatz) [9, 10, 11, 12,
13], along with the time-dependent variational principle, has been shown to be
an excellent tool for accurately simulating the dynamics of system-bath models
[9, 14, 15, 16, 17, 18, 19, 20] and spectroscopic signals [21, 22, 19, 20,
23]. Despite relying on an adaptive, time-dependent state basis set, the
problem of rapidly growing numerical costs remain.

In a previous study, we proposed the thermalization algorithm [1] to be used
with the Davydov $\text{D}_{2}$ ansatz [24, 25, 14, 26, 27, 28, 29], which
restricts QHOs to their lowest uncertainty states – coherent states [30, 31]
with Gaussian wavepackets in their coordinate-momentum phase space. We
demonstrated how the thermalization algorithm regulates the temperature of the
environment and enables the electronically excited molecular system to relax
into its equilibrium state at a given temperature [32, 33, 34] even when using
a reduced number of bath oscillators, which greatly reduces numerical costs.
The characteristics of the resulting equilibrium state are essential for
modeling fluorescence, excited state emission, excited state absorption and
other spectroscopic signals [2]. However, the $\text{D}_{2}$ ansatz is a crude
approximation of the actual system-bath model eigenstates and thus is unable
to completely capture electronic population relaxation dynamics [22].

Meanwhile, the system-bath dynamics obtained using the multiple-Davydov
ansätze are consistent with the results from other state-of-the-art methods,
such as hierarchical equations of motion [22, 10, 12], quasiadiabatic path
integral [15], and multi-configuration time-dependent Hartree [16, 35], even
when the number of bath oscillators is large. Due to the more complicated
wavefunction structure of the $\text{mD}_{2}$ ansatz, straightforward
application of the $\text{D}_{2}$ ansatz thermalization algorithm is not
possible. In this work, we extend the thermalization algorithm for the
$\text{mD}_{2}$ ansatz by introducing an additional state projection algorithm
and adopting the coarse-grained scattering approximation.

In Section (II) we describe the thermalization algorithm for the
$\text{mD}_{2}$ ansatz, and in Section (III) we provide a theoretical
description of its application to simulating the fluorescence spectra. Then,
in Section (IV) we demonstrate its capabilities by simulating excitation
relaxation dynamics of an H-type molecular aggregate and its fluorescence
spectrum. Lastly, in Section (V) we discuss changes made to adapt the
thermalization algorithm of the $\text{D}_{2}$ ansatz for the $\text{mD}_{2}$
ansatz.

##  II Thermalization of the $\text{mD}_{2}$ ansatz

We consider a molecular aggregate model, where each molecule $n=1,2,\ldots,N$
couples to its own local reservoir $k=1,2,\ldots,N$, each of which consists of
$q=1,2,\ldots,Q$ QHO modes. The model is given by the Hamiltonian
$\hat{H}=\hat{H}_{\text{S}}+\hat{H}_{\text{B}}+\hat{H}_{\text{SB}}$ with the
system, the bath and the system-bath coupling terms defined as:

$\displaystyle\hat{H}_{\text{S}}=$
$\displaystyle\sum_{n}^{N}\varepsilon_{n}\hat{a}_{n}^{\dagger}\hat{a}_{n}+\sum_{n,m}^{n\neq
m}J_{nm}\hat{a}_{n}^{\dagger}\hat{a}_{m},$ (1)
$\displaystyle\hat{H}_{\text{B}}=$
$\displaystyle\sum_{k,q}^{N,Q}\omega_{kq}\hat{b}_{kq}^{\dagger}\hat{b}_{kq},$
(2) $\displaystyle\hat{H}_{\text{SB}}=$
$\displaystyle-\sum_{n}^{N}\hat{a}_{n}^{\dagger}\hat{a}_{n}\sum_{q}^{Q}\omega_{nq}g_{nq}\left(\hat{b}_{nq}^{\dagger}+\hat{b}_{nq}\right),$
(3)

with the reduced Planck’s constant set to $\hbar=1$. Here $\varepsilon_{n}$ is
the $n$-th molecule electronic excitation energy, $J_{nm}$ denotes the
resonant coupling between the $n$-th and $m$-th molecules, $\omega_{nq}$
denotes the frequency of the $q$-th QHO in the $k$-th local reservoir, while
$g_{nq}$ is the coupling strength between the $q$-th oscillator in the $n$-th
local reservoir to the $n$-th molecule. The operators $\hat{a}_{n}^{\dagger}$
and $\hat{a}_{n}$ represent the creation and annihilation operators for
electronic excitations, respectively, while $\hat{b}_{nq}^{\dagger}$ and
$\hat{b}_{nq}$ represent the creation and annihilation bosonic operators for
QHOs.

In addition, we implicitly couple the system-bath model to the secondary bath
characterized by a fixed temperature, $T.$ The coupling between the secondary
and primary baths occurs via the scattering events that allow the system-bath
model to exchange energy with the secondary bath and thermalize local
reservoirs as is described below.

A state of the system-bath model is given by the $\text{mD}_{2}$ wavefunction

$|\Psi\left(t\right)\rangle=\sum_{i,n}^{M,N}\alpha_{i,n}\left(t\right)|n\rangle\otimes|\boldsymbol{\lambda}_{i}\left(t\right)\rangle,$
(4)

where $\alpha_{i,n}\left(t\right)$ is the $i$-th multiple complex amplitude
associated with a singly excited electronic state $|n\rangle$ localized on the
$n$-th molecule, $|n\rangle=\hat{a}_{n}^{\dagger}|0\rangle_{el}$.
$|0\rangle_{el}$ is the electronic ground state. The complexity and accuracy
of the $\text{mD}_{2}$ ansatz can be adjusted by varying the multiplicity
number, $M$. The states of QHO modes are represented by multi-dimensional
coherent states

$\displaystyle|\boldsymbol{\lambda}_{i}\left(t\right)\rangle$
$\displaystyle=\exp\sum_{k,q}^{N,Q}\left(\lambda_{i,kq}\left(t\right)\hat{b}_{kq}^{\dagger}-\lambda_{i,kq}^{\star}\left(t\right)\hat{b}_{kq}\right)|0\rangle_{\text{vib}},$
(5)

with $\lambda_{i,kq}\left(t\right)$ being the $i$-th multiple complex
displacement parameter and
$|0\rangle_{\text{vib}}=\otimes_{k,q}|0\rangle_{k,q}$ is the global
vibrational ground state of all QHOs.

The $\text{mD}_{2}$ wavefunction describes a state of the system-bath model as
a superposition of $M$ multi-dimensional coherent state terms, which allows it
to represent a wide range of system-bath model states beyond the Born-
Oppenheimer and Gaussian approximations.

The thermalization algorithm for the $\text{mD}_{2}$ ansatz is realized by
stochastic scattering events [36, 37] during time evolution of the
wavefunction. These events change momenta, $p_{kq}$, of all $q$-th QHO modes
of the $k$-th local reservoir at once. We assume that the scattering
probability, $P_{k}\left(\theta,\tau_{\text{sc}}\right)$, of $\theta$
scattering events occurring per time interval, $\tau_{\text{sc}}$, with a
scattering rate, $\nu_{k}$, is given by a Poisson distribution

$P_{k}\left(\theta,\tau_{\text{sc}}\right)=\frac{1}{\theta!}\left(\tau_{\text{sc}}\nu_{k}\right)^{\theta}\text{e}^{-\tau_{\text{sc}}\nu_{k}}.$
(6)

Numerically, Poisson statistics are realised by simulating Bernoulli processes
[38, 39] in the limit of $\tau_{\text{sc}}\rightarrow 0$, while maintaining
the condition that $\nu_{k}\tau_{\text{sc}}\ll 1$. To simulate the scattering
events we divide wavefunction propagation into equal time length,
$\tau_{\text{sc}}$, intervals

$t_{i}=\left(i\tau_{\text{sc}},\left(i+1\right)\tau_{\text{sc}}\right],\
i=0,1,\ldots.$ (7)

At the end of each interval, for each local reservoir, we flip a biased coin
with the probability $\nu_{k}\tau_{\text{sc}}$ of landing heads for all local
reservours. If a $k$-th coin lands heads, we change momenta of _all oscillator
modes of that $k$-th reservoir,_ otherwise, no changes are made. A list of
scattering moments at which the numerical simulation is paused to perform the
scatterings can be precomputed prior to starting the simulation by drawing
probabilities for all local reservoir and all time intervals $t_{i}$ from Eq.
(6).

We assume that during the scattering event the local bath, which experiences
the scattering, acquires thermal-equilibrium kinetic energy. Such state is
given by a single coherent state for one specific QHO. In order to set the new
momenta values of the scattered reservoir oscillator modes, we first project
the $\text{mD}_{2}$ wavefunction of Eq. (4) to its single-multiple Davydov
$\text{D}_{2}$ form

$|\psi\left(t\right)\rangle=\sum_{n}^{N}\beta_{n}\left(t\right)|n\rangle\otimes|\tilde{\boldsymbol{\lambda}}\left(t\right)\rangle,$
(8)

where $\beta_{n}$ are the projected complex electronic amplitudes and
$|\tilde{\boldsymbol{\lambda}}\left(t\right)\rangle$ is the projected multi-
dimensional coherent state, which is defined later. This follows the
decoherence idea [40], where the macroscopic environment perform a collapse of
the wavefunction into a set of preferred states, in our case, the electronic-
vibrational states,
$|n\rangle\otimes|\tilde{\boldsymbol{\lambda}}\left(t\right)\rangle$. The
projected complex electronic amplitudes are equal to

$\beta_{n}\left(t\right)=\sum_{i}^{M}\alpha_{i,n}\left(t\right)\langle\tilde{\boldsymbol{\lambda}}\left(t\right)|\boldsymbol{\lambda}_{i}\left(t\right)\rangle,$
(9)

while the projected multi-dimensional coherent state

$\displaystyle|\tilde{\boldsymbol{\lambda}}\left(t\right)\rangle$
$\displaystyle=\exp\sum_{k,q}^{N,Q}\left(\tilde{\lambda}_{kq}\left(t\right)\hat{b}_{kq}^{\dagger}-\tilde{\lambda}_{kq}^{\star}\left(t\right)\hat{b}_{kq}\right)|0\rangle_{\text{vib}}$
(10)

is defined in terms of the complex displacements

$\displaystyle\tilde{\lambda}_{kq}\left(t\right)$
$\displaystyle=\frac{1}{\sqrt{2}}\left(x_{kq}\left(t\right)+\text{i}p_{kq}\left(t\right)\right),$
(11)

where $x_{kq}\left(t\right)$ and $p_{kq}\left(t\right)$ are QHO coordinate and
momentum expectation values

$\displaystyle x_{kq}=$
$\displaystyle\frac{1}{\sqrt{2}}\sum_{i,j,n}^{M,M,N}\alpha_{i,n}^{\star}\alpha_{j,n}\langle\boldsymbol{\lambda}_{i}|\boldsymbol{\lambda}_{j}\rangle\sum_{k,q}^{N,Q}\left(\lambda_{i,kq}^{\star}+\lambda_{j,kq}\right),$
(12) $\displaystyle p_{kq}=$
$\displaystyle\frac{\text{i}}{\sqrt{2}}\sum_{i,j,n}^{M,M,N}\alpha_{i,n}^{\star}\alpha_{j,n}\langle\boldsymbol{\lambda}_{i}|\boldsymbol{\lambda}_{j}\rangle\sum_{k,q}^{N,Q}\left(\lambda_{i,kq}^{\star}-\lambda_{j,kq}\right),$
(13)

calculated from the $\text{mD}_{2}$ ansatz, where
$\langle\boldsymbol{\lambda}_{i}|\boldsymbol{\lambda}_{j}\rangle$ is the
overlap of two coherent states

$\langle\boldsymbol{\lambda}_{i}|\boldsymbol{\lambda}_{j}\rangle=\exp\sum_{k,q}^{N,Q}\left(\lambda_{i,kq}^{\star}\lambda_{j,kq}-\frac{1}{2}\left(\left|\lambda_{i,kq}\right|^{2}+\left|\lambda_{j,kq}\right|^{2}\right)\right).$
(14)

This completes the projection operation of the $\text{mD}_{2}$ state, given by
Eq. (4), into its simplified $\text{D}_{2}$ form in Eq. (8).

Once the projected wavefunction is deduced, we modify momenta of the scattered
oscillators by sampling the QHO diagonal density operator distribution in the
coherent state representation at temperature $T$, known as the Glauber-
Sudarshan distribution [41, 42, 43, 44]

$\mathcal{P}\left(\tilde{\lambda}_{kq}\right)=\mathcal{Z}_{kq}^{-1}\exp\left(-\left|\tilde{\lambda}_{kq}\right|^{2}\left(\text{e}^{\frac{\omega_{kq}}{k_{\text{B}}T}}-1\right)\right).$
(15)

For scattered modes, we set momenta values in Eq. (11) to

$p_{kq}\left(t\right)=\sqrt{2}\text{Im}\left(\tilde{\lambda}_{kq}^{\mathcal{P}}\right),$
(16)

where $\tilde{\lambda}_{kq}^{\mathcal{P}}$ is a sample drawn from the Glauber-
Sudarshan distribution. $\mathcal{Z}_{kq}^{-1}$ and $\omega_{kq}$ are
partition functions and frequencies of QHO, respectively, while $k_{\text{B}}$
is the Boltzmann constant. During the scattering events, coordinates,
$x_{kq}$, of both the scattered and non-scattered modes remain unchanged.
Notice, that the local baths, which do not experience scattering, remain
unaffected by the scattering of other modes.

Now that the wavefunction of the system-bath model after scattering is known
(given by Eq. (8)) we rewrite it in the $\text{mD}_{2}$ wavefunction form of
Eq. (4) by populating amplitudes and displacements of the first multiple,
$i=1$, as

$\displaystyle\alpha_{1,n}\left(t\right)$
$\displaystyle=\beta_{n}\left(t\right),$ (17)
$\displaystyle\lambda_{1,kq}\left(t\right)$
$\displaystyle=\tilde{\lambda}_{kq}\left(t\right).$ (18)

Amplitudes of the unpopulated multiples are set to
$\alpha_{j=2,\ldots,M,n}\left(t\right)=0$, while the unpopulated displacements
are positioned in a layered hexagonal pattern around the populated coherent
state [19]

$\displaystyle\lambda_{j=2,\ldots,M,kq}\left(t\right)$
$\displaystyle=\lambda_{1,kq}\left(\tau\right)$
$\displaystyle+\frac{1}{4}\left(1+\left\lfloor\beta\left(j\right)\right\rfloor\right)e^{\text{i}2\pi\left(\beta\left(j\right)+\frac{1}{12}\left\lfloor\beta\left(j\right)\right\rfloor\right)},$
(19)

where $\beta\left(j\right)=\left(j-2\right)/6$ is the coordination function
and $\left\lfloor x\right\rfloor$ is the floor function. The exact arrangement
of displacements of the unpopulated multiples is not critical, as long as the
distance in the phase space to the populated multiple coherent state is not
too large, otherwise, the initially unpopulated multiples will not contribute
to further dynamics [17, 45].

Once the scattered $\text{mD}_{2}$ wavefunction is determined and the
scattering event is finalised, further simulation of $\text{mD}_{2}$ dynamics
according to equations of motion proceeds. This procedure generates a
stochastic wavefunction trajectory, where the system-bath model at each time
moment is described by a pure state, which is a single member of a statistical
ensemble [42, 43]. The thermalized model dynamics are obtained by averaging
observables over an ensemble of wavefunction trajectories, $\gamma$, which
differ by their: initial amplitudes, $\alpha_{i,n}\left(0\right)$, initial
coherent state displacements, $\lambda_{i,kq}\left(0\right)$, and a sequence
of scattering events. Ensemble averaging is performed in a parallelized Monte
Carlo scheme.

##  III Thermalized fluorescence spectra

Wavefunction trajectories allows calculation of an arbitrary observable.
Calculation of equilibrium fluorescence spectrum requires to know thermally
equilibrated state of the excited model. The presented thermalization
procedure allows to obtain such state and calculate fluorescence spectrum.

In general, the frequency-domain spectrum of a quantum system can be written
as a Fourier transform

$F\left(\omega\right)=\text{Re}\int_{0}^{\infty}\text{d}t\text{e}^{i\omega
t}S\left(t\right),$ (20)

of the corresponding time-domain response function, $S\left(t\right)$.

The fluorescence (FL) response function, $S_{\text{fl}}\left(t\right)$, is a
specific case of the more general time-resolved fluorescence (TRF) response
function [2, 46]

$\displaystyle S_{\text{trf}}\left(\tau_{\text{eq}},t\right)$
$\displaystyle=\frac{1}{\Gamma}\sum_{\gamma=1}^{\Gamma}\langle\Psi_{\text{G}}\left(0\right)|_{\gamma}\hat{\mu}_{-}\hat{\mathcal{V}}_{\text{E}}^{\dagger}\left(\tau_{\text{eq}}+t\right)\hat{\mu}_{+}$
$\displaystyle\times\hat{\mathcal{V}}_{\text{G}}\left(t\right)\hat{\mu}_{-}\hat{\mathcal{V}}_{\text{E}}\left(\tau_{\text{eq}}\right)\hat{\mu}_{+}|\Psi_{\text{G}}\left(0\right)\rangle_{\gamma},$
(21)

where $\hat{\mathcal{V}}_{\text{E}}$ and $\hat{\mathcal{V}}_{\text{G}}$ are
the excited and ground state system-bath propagators

$\displaystyle\hat{\mathcal{V}}_{A}\left(t_{1}\right)|\Psi_{A}\left(t_{2}\right)\rangle$
$\displaystyle=|\Psi_{A}\left(t_{1}+t_{2}\right)\rangle,$ (22)

while
$\hat{\mu}_{+}=\sum_{n}^{N}\left(\boldsymbol{e}\cdot\boldsymbol{\mu}_{n}\right)\hat{a}_{n}^{\dagger}$
and
$\hat{\mu}_{-}=\sum_{n}^{N}\left(\boldsymbol{e}\cdot\boldsymbol{\mu}_{n}\right)\hat{a}_{n}$
are the excitation creation and annihilation operators of the system [19],
$\boldsymbol{\mu}_{n}$ is the electronic transition dipole moment vector,
$\boldsymbol{e}$ is the external field polarization vector.
$|\Psi_{\text{G}}\left(0\right)\rangle_{\gamma}$ is a model ground state with
an initial condition of the $\gamma$-th trajectory. The EOMs for propagating
the $\text{mD}_{2}$ wavefunction, as well as the approach to solving them, are
described in detail in Refs. [45, 19].

$S_{\text{trf}}\left(\tau_{\text{eq}},t\right)$ is a function of two times:
the equilibration time, $\tau_{\text{eq}}$, and the coherence time, $t$.
During the equilibration time, the system evolves in its excited state and,
due to the system-bath interaction, relaxes to an equilibrium state. After
this, during the coherence time, spontaneous emission occurs.

We will apply thermalization during the equilibration time to facilitate the
relaxation of the system-bath model into the lowest energy equilibrium state
by removing excess thermal energy from local reservoirs. We denote
$\hat{\mathcal{G}}_{\text{E},\gamma}$ as the excited state propagator
$\hat{\mathcal{V}}_{\text{E}}$, but with the thermalization. Then the
thermalized TRF (tTRF) response function can be written as

$\displaystyle\tilde{S}_{\text{trf}}\left(\tau_{\text{eq}},t\right)$
$\displaystyle=\frac{1}{\Gamma}\sum_{\gamma=1}^{\Gamma}\langle\Psi_{\text{G}}\left(0\right)|_{\gamma}\hat{\mu}_{-}\hat{\mathcal{G}}_{\text{E},\gamma}^{\dagger}\left(\tau_{\text{eq}}\right)\hat{\mathcal{V}}_{\text{G}}^{\dagger}\left(t\right)\hat{\mu}_{+}$
$\displaystyle\times\hat{\mathcal{V}}_{\text{G}}\left(t\right)\hat{\mu}_{-}\hat{\mathcal{G}}_{\text{E},\gamma}\left(\tau_{\text{eq}}\right)\hat{\mu}_{+}|\Psi_{\text{G}}\left(0\right)\rangle_{\gamma}.$
(23)

By considering the equilibration time to be long enough to reach the
equilibrium state of the system-bath model, we define the FL response function
to be

$S_{\text{fl}}\left(t\right)=\lim_{\tau_{\text{eq}}\rightarrow\infty}S_{\text{trf}}\left(\tau_{\text{eq}},t\right),$
(24)

and the thermalized fluorescence (tFL) response function as

$\tilde{S}_{\text{fl}}\left(t\right)=\lim_{\tau_{\text{eq}}\rightarrow\infty}\tilde{S}_{\text{trf}}\left(\tau_{\text{eq}},t\right).$
(25)

The spectra obtained using the fluorescence response function without and with
thermalization will be compared in the next section. For the numerical
simulation, the required equilibration time interval has to be deduced by
increasing $\tau_{\text{eq}}$ until the resulting fluorescence spectra
converges.

##  IV Results

To investigate the thermalization algorithm for the $\text{mD}_{2}$ ansatz, we
will analyse the linear trimer model, which we previously used to study
thermalization of the $\text{D}_{2}$ ansatz [1]. The model consists of $N=3$
coupled molecules, with excited state energies $\varepsilon_{n}$ being equal
to $0,\ 250,\ 500\ \text{cm}^{-1}$, forming an energy funnel. The nearest
neighbour couplings are set to $J_{1,2}=J_{2,3}=100\ \text{cm}^{-1}$,
$J_{3,1}=0$. The electronic dipole moment vectors of molecules are
$\boldsymbol{\mu}_{n}=(1,0,0)$ in Cartesian coordinate system. This classifies
the trimer as the H-type molecular aggregate [47].

QHOs of local molecular reservoirs are characterized by the super-Ohmic [48]
spectral density function
$C^{\prime\prime}\left(\omega\right)=\omega\left(\omega/\omega_{c}\right)^{s-1}\exp\left(-\omega/\omega_{c}\right)$
with an order parameter $s=2$ and a cut-off frequency $\omega_{\text{c}}=100\
\text{cm}^{-1}$. The QHO frequencies are
$\omega_{kq}=\omega_{0}+\left(q-1\right)\Delta\omega$, where the frequency
off-set is $\omega_{0}=0.01\ \text{cm}^{-1}$. The reorganization energy of
each local reservoir is $\Lambda_{k}=\sum_{q}\omega_{kq}g_{kq}^{2}=100\
\text{cm}^{-1}$. The scattering time step-size is set to
$\tau_{\text{sc}}=0.01\ \text{ps}$. Finally, the ensemble consists of $900$
wavefunctions trajectories, which we found to be sufficient to obtain the
converged model dynamics. The $\text{mD}_{2}$ ansatz multiplicity $M=5$ is
used as the results with a higher multiplicity quantitatively match the $M=5$
case.

We will be considering three bath models: the dense bath model, where the
spectral density function $C^{\prime\prime}\left(\omega\right)$ is discretized
into $Q=75$ oscillators per local reservoir with a step-size of
$\Delta\omega=10\ \text{cm}^{-1}$; the sparse bath model, where the number of
modes is reduced by a factor of 5 to just $Q=15$ oscillators per local
reservoir with $\Delta\omega=50\ \text{cm}^{-1}$; and the sparse bath with
thermalization model, where $C^{\prime\prime}\left(\omega\right)$ is
discretized according to the sparse bath model and thermalization is used.

In the absence of the bath, the system has three single-excitation stationary
exciton states with energies: $E_{1}^{\text{exc}}\approx-37.23\
\text{cm}^{-1}$, $E_{2}^{\text{exc}}=250\ \text{cm}^{-1}$,
$E_{3}^{\text{exc}}\approx 537.23\ \text{cm}^{-1}$, satisfying the time-
independent Schrödinger equation

$\hat{H}_{\text{S}}\Phi_{n}^{\text{exc}}=E_{n}^{\text{exc}}\Phi_{n}^{\text{exc}},$
(26)

with the system Hamiltonian given by Eq. (1). The exciton eigenstates [3, 5],
$\Phi_{n}^{\text{exc}}$, have their excitations delocalized over multiple
molecules [42]. Therefore, it is convenient to analyse molecular aggregate
excitation relaxation dynamics in terms of excitons as they are natural quasi-
particles of the aggregate. We denote the probability of finding the aggregate
in its $n$-th excitonic state as the population, given by

$\rho_{n}^{\text{exc}}\left(t\right)=\text{{\hbox{\sum_{k,l,i,j}}}}\left(\Phi_{k}^{\text{exc}}\right)_{n}^{\star}\left\langle\alpha_{i,k}^{\star}\left(t\right)\alpha_{j,l}\left(t\right)S_{i,j}\left(t\right)\right\rangle_{\text{th}}\left(\Phi_{l}^{\text{exc}}\right)_{n},$
(27)

where $\left\langle\ldots\right\rangle_{\text{th}}$ is the averaging over an
ensemble of wavefunction trajectories.

Using $\text{mD}_{2}$ ansatz we have

Figure 1: The exciton state populations,
$\rho_{n}^{\text{exc}}\left(t\right)$, and the average temperatures,
$T_{k}\left(t\right)$, of local reservoirs of the trimer with (a, b) the dense
bath model, (c, d) the sparse bath model and (e, f) the sparse bath model with
thermalization.

First, we study the electronic excitation dynamics. The initial excitonic
state populations correspond to the optically excited highest energy state:
$\rho_{3}^{\text{exc}}=1$, $\rho_{1,2}^{\text{exc}}=0$, while the initial QHO
displacements, $\lambda_{i,kq}\left(0\right)$, are sampled from the Glauber-
Sudarshan distribution in Eq. (15) to account for the initial temperatures of
$T_{k}=77\ \text{K}$.

In Fig. (1) we display the trimer model exciton state populations
$\rho_{n}^{\text{exc}}\left(t\right)$ and average temperatures [49]
$T_{k}\left(t\right)$ of local reservoirs for all three bath models. The
aggregate environment causes dephasing between excitonic states and induces
irreversible population relaxation [2, 3]. The population dynamics of the
dense bath models exhibits a sequential relaxation from the initially
populated highest energy excitonic state to the lowest energy state via the
intermediate state. Eventually, the population distribution reaches the
equilibrium state. The majority of the excitation energy is transferred to
oscillators of local reservoirs. We observe an increase of temperatures [50,
49, 51] due to the finite number of oscillators in local reservoirs. An
infinite number of oscillators would have to be included to maintain
temperature constant at the initial value. The initial rapid rise in
temperature is due to oscillator reorganization in the aggregate’s electronic
excited state manifold, while the following slow rise is due to energy
transfer from the system to local reservoirs.

In the sparse bath model, we observe that if the number of vibrational modes
is reduced, the population dynamics become skewed due to insufficiently dense
representation of the spectral density function. Furthermore, the temperature
increase is higher than compared to the dense bath model, which further
changes the characteristics of the resulting equilibrium state.

When the thermalization algorithm is applied to the sparse bath model with a
scattering rate $\nu_{k}=1.25\ \text{ps}^{-1}$, the population dynamics are
restored and qualitatively match those of the dense model. Although the
initial temperatures of the local reservoirs exceed those of the dense bath
model, they gradually decrease due to thermalization, and this rate can be
adjusted by changing the scattering rate.

Figure 2: (a) the TRF and (b) tTRF spectra of the trimer with the dense bath
model, simulated with an increasing equilibration time $\tau_{\text{eq}}$. The
absorption spectrum is also shown. Vertical dashed lines show energies
$E^{\text{exc}}$ of the excitonic states.

Next, we turn our attention to simulating the FL spectrum of the linear trimer
with the dense bath model with scattering rate $\nu_{k}=1\ \text{ps}^{-1}$.
The initial excitonic state population distribution is now calculated in terms
of the system-field interaction, as described in Ref. [19].

In Fig. (2) we compare the TRF and tTRF spectra with increasing equilibration
times $\tau_{\text{eq}}$. When $\tau_{\text{eq}}=0$, both the TRF and tTRF
spectra are equivalent and exactly match the absorption spectrum, which
consists of three peaks due to transition involving the combined excitonic-
vibrational (vibronic) states and can not be regarded as purely excitonic. For
reference, vertical dotted lines indicate energies $E^{\text{exc}}$ of
excitonic states. These do not match the three peak energies exactly due to
the system being coupled to the environment.

By allowing equilibration to occur, $\tau_{\text{eq}}>0$, both the TRF and
tTRF spectra show peak intensity shift towards lower energies as excitation
relaxes towards the equilibrated state during the equilibration time. After
equilibrating for $\tau_{\text{eq}}=2\ \text{ps}$, we find that both spectra
have converged and do not change with longer $\tau_{\text{eq}}$. Therefore,
the TRF and tTRF spectra at $\tau_{\text{eq}}=2\ \text{ps}$ can be considered
as the FL and tFL spectra of the trimer model as defined in Eqs. (24), (25).

Both spectra exhibit their highest intensities at the energies of the lowest
vibronic states. However, the FL spectrum also has considerable intensities at
energies of the intermediate and highest vibronic states. Surprisingly, the
higher energy peak is more intense than the intermediate peak. The tFL
spectrum intensities at these energies are negligible, which indicate that the
thermalization allows the trimer model to reach a lower energy equilibrium
state, which is no longer hindered by the excess thermal energy accumulation
in QHOs of local reservoirs.

In Fig. (3), we also compare the obtained FL and tFL spectra with the FL
spectrum simulated using a previously proposed excited state numerical
optimization approach [20, 52, 53]. It relies on finding the model’s lowest
energy excitonic state in terms of the $\text{mD}_{2}$ ansatz parameters and
then applying thermal fluctuations to effectively generate the model in a
lowest energy equilibrium state at the temperature of $77\ \text{K}$. We see
that all three methods produce a similar lowest vibronic peak, but the tFL
spectrum has a higher intensity tail towards the low energy side and almost no
intensities at energies of the intermediate and the highest vibronic states,
while the FL spectrum simulated using the optimization approach has a small
intensity at the energy of the intermediate vibronic states. The optimization
approach spectrum more closely resembles that of the thermalized model than
the non-thermalized spectrum.

Figure 3: The fluorescence spectra comparison of the trimer with the dense
bath model obtained without thermalization, with thermalization and using the
optimization approach. The equilibration time is $\tau_{\text{eq}}=2\
\text{ps}$. The absorption spectrum is also shown. Vertical dashed lines show
energies $E^{\text{exc}}$ of the excitonic states.

##  V Discussion

Starting from an arbitrary non-equilibrium initial condition a closed quantum
system will not equilibrate due to energy conservation. Thermalization
procedure is necessary to guarantee proper thermal equilibrium in the long run
for _all_ bath oscillators. This requires introducing the concept of primary
and secondary baths. In our model the primary bath is a part of explicit
quantum DOFs, while the secondary bath is a thermal reservour with infinite
thermal capacity, i. e. it keeps constant temperature in any energy exchange
process. In this case, the secondary bath cannot be described by mechanical
equations – only statistical or thermodynamical conepts apply. Our statistical
algorithm performs energy exchange between the primary and secondary baths
using the statistical scattering idea: the primary bath state is being reset
to the thermally equilibrated state, thus giving up excess energy to or
drawing additional energy from the secondary bath. This is a major extension
of the explicit quantum TDVP theory – the extended model covers a broader
range of phenomena: local heating and cooling, as well as bath oscillator
dynamic localization, which are not available in the standard TDVP theory.

In order to adapt the $\text{D}_{2}$ ansatz thermalization algorithm for the
$\text{mD}_{2}$ ansatz, several extensions were made. During the time
evolution of the system-bath model, the $\text{mD}_{2}$ ansatz multiples
become correlated, leading to a non-Gaussian bath wavefunction. It becomes
impossible to represent a new Gaussian wavefunction of scattered QHO modes,
sampled from Eq. (15), without changing the wavefunction of all the rest non-
scattered oscillators at the same time. Therefore, we chose to project the
$\text{mD}_{2}$ into the $\text{D}_{2}$ ansatz whenever scatterings occurs,
allowing to correctly represent the newly sampled Gaussian wavefunction of
scattered oscillators. This idea requires the consideration of a few aspects.

The projected $\text{D}_{2}$ wavefunction accurately maintains average
coordinates and momenta of the $\text{mD}_{2}$ ansatz QHO states, while
variances and higher-order moments become affected. This causes variation of
excitation relaxation dynamics compared to the standard $\text{mD}_{2}$
ansatz. However, system-bath models mostly rely only on the linear coupling
between the system and average coordinates of QHO modes, therefore, as seen in
Fig. (1), the discrepancy is minimal. The higher-order couplings become
necessary when anharmonic vibrational modes or changes to their frequencies
upon excitation are considered [54, 20].

To maintain the close correspondence to the standard $\text{mD}_{2}$ ansatz
the projection should not occur too often. This is because it takes time for
the wavefunction after scattering to again become correlated between its many
multiples, i.e., to take advantage of the unpopulated $\text{mD}_{2}$ ansatz
multiples after projection. If the repopulation time is shorter than the time
between projection operations, the model population dynamics become similar to
those of the $\text{D}_{2}$ ansatz, even though the $\text{mD}_{2}$ ansatz is
being used. The average time interval between projection operations is
determined by the scattering rate $\nu_{k}$, a property of the physical
system, while the scattering time $\tau_{\text{sc}}$ is a parameter of the
model and must be as small as necessary to ensure the Bernoulli-to-Poisson
statistics transition condition, $\nu_{k}\tau_{\text{sc}}\ll 1$.

To increase the average time between projection operations, we adopt a coarser
scattering approach for the $\text{mD}_{2}$ ansatz as compared to the
$\text{D}_{2}$ ansatz. Instead of considering scattering events of individual
oscillators, we consider events, where all oscillators of certain local
reservoirs are scattered at once, requiring only a single projection operation
to scatter many oscillators at once. This approach allows the $\text{mD}_{2}$
ansatz to continue utilizing all its multiples for the improved accuracy over
the $\text{D}_{2}$ ansatz, while reducing the number of explicitly modelled
oscillators needed to maintain local reservoirs’ temperatures close to initial
values, thereby reducing the numerical cost.

Using the $\text{mD}_{2}$ ansatz to simulate population dynamics of the trimer
model, it took an average of 166 minutes per trajectory using the dense bath
model, but only 1.3 minutes using the sparse bath model and 2 minutes using
the sparse bath model with thermalization. The computational overhead of
thermalization is small compared to the overall time savings when switching
from using the dense bath to the sparse bath. The numerical cost reduction is
also greater for the $\text{mD}_{2}$ ansatz than for the $\text{D}_{2}$ ansatz
in Ref. [1], because the $\text{mD}_{2}$ ansatz EOMs constitute an implicit
system of differential equations, which require a more involved, two-step
numerical approach to find a solution [19, 45]. By considering fewer
oscillators in each local reservoir, simulations of the dynamics and
spectroscopic signals of aggregates made up of more molecules becomes
possible.

Computing a single trajectory of the tTRF response function in Eq. (23) with
an equilibration time of $\tau_{\text{eq}}=2\ \text{ps}$ took an average of 79
minutes. The previously proposed optimization method [20] for simulating FL
spectra does not require propagation during the equilibration time interval of
the TRF response function and has to be computed only once, but it took 193
minutes. In general, we find that the computation of tTRF is more reliable and
numerically stable. The optimization approach struggles to consistently find
the lowest energy excitonic state of the model due to its heuristic nature,
requiring many attempts to find the solution and eventually having to choose
the lowest energy one. This is particularly apparent when a wide range of
oscillator frequencies are included.

For elementary system-bath models without Hamiltonian parameter disorder, the
optimization approach can be a good starting point for FL spectra simulation.
However, a more accurate spectra most likely will be obtained using the tTRF
approach. For models with Hamiltonian disorder, e.g., static molecule
excitation energy disorder [55, 56, 57], the optimization approach would
require finding model’s lowest energy excitonic state for each realization of
the Hamiltonian, negating its advantage of having to perform optimization
procedure only once.

In conclusion, the presented thermalization algorithm for the numerically
exact $\text{mD}_{2}$ ansatz allows to reduce numerical cost of system-bath
model simulations by having to explicitly include fewer bath oscillators,
while maintaining correspondence with the exact relaxation dynamics. The
thermalization algorithm efficiently controls molecular heating effects due to
the reduced number of oscillators. Furthermore, the application of
thermalization to simulation of fluorescence spectra demonstrates lower
computation time, greater numerically stability and higher accuracy compared
to the numerical optimization approach.

## Conflicts of interest

There are no conflicts of interest to declare.

###### Acknowledgements.

We thank the Research Council of Lithuania for financial support (grant No:
S-MIP-23-48). Computations were performed on resources at the High Performance
Computing Center, “HPC Sauletekis” in Vilnius University Faculty of Physics.

## References

  * Jakučionis and Abramavičius [2021] M. Jakučionis and D. Abramavičius, Temperature-controlled open-quantum-system dynamics using a time-dependent variational method, Physical Review A 103, 032202 (2021).
  * Mukamel [1995] S. Mukamel,  _Principles of nonlinear optical spectroscopy_ (Oxford University Press, 1995).
  * Valkunas  _et al._ [2013] L. Valkunas, D. Abramavicius, and T. Mančal,  _Molecular Excitation Dynamics and Relaxation_ (Wiley-VCH Verlag GmbH, 2013).
  * Bardeen [2014] C. J. Bardeen, The Structure and Dynamics of Molecular Excitons, Annual Review of Physical Chemistry 65, 127 (2014).
  * van Amerongen  _et al._ [2000] H. van Amerongen, R. van Grondelle, and L. Valkunas,  _Photosynthetic Excitons_ (World Scientific, 2000).
  * Schröter  _et al._ [2015] M. Schröter, S. Ivanov, J. Schulze, S. Polyutov, Y. Yan, T. Pullerits, and O. Kühn, Exciton-vibrational coupling in the dynamics and spectroscopy of Frenkel excitons in molecular aggregates, Physics Reports 567, 1 (2015).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Weiss [2012] U. Weiss,  _Quantum Dissipative Systems_ (WORLD SCIENTIFIC, 2012).
  * Zhou  _et al._ [2015] N. Zhou, Z. Huang, J. Zhu, V. Chernyak, and Y. Zhao, Polaron dynamics with a multitude of Davydov D 2 trial states, The Journal of Chemical Physics 143, 014113 (2015).
  * Wang  _et al._ [2016] L. Wang, L. Chen, N. Zhou, and Y. Zhao, Variational dynamics of the sub-Ohmic spin-boson model on the basis of multiple Davydov D 1 states, The Journal of Chemical Physics 144, 024101 (2016).
  * Zhou  _et al._ [2016a] N. Zhou, L. Chen, Z. Huang, K. Sun, Y. Tanimura, and Y. Zhao, Fast, Accurate Simulation of Polaron Dynamics and Multidimensional Spectroscopy by Multiple Davydov Trial States, Journal of Physical Chemistry A 120, 1562 (2016a).
  * Chen  _et al._ [2018] L. Chen, M. Gelin, and Y. Zhao, Dynamics of the spin-boson model: A comparison of the multiple Davydov D1,D1.5,D2 Ansätze, Chemical Physics 515, 108 (2018).
  * Zhao  _et al._ [2022] Y. Zhao, K. Sun, L. Chen, and M. Gelin, The hierarchy of Davydov’s Ansätze and its applications, WIREs Computational Molecular Science 12, e1589 (2022).
  * Huang  _et al._ [2017] Z. Huang, L. Chen, N. Zhou, and Y. Zhao, Transient dynamics of a one-dimensional Holstein polaron under the influence of an external electric field, Annalen der Physik 529, 1600367 (2017).
  * Wang  _et al._ [2017a] L. Wang, Y. Fujihashi, L. Chen, and Y. Zhao, Finite-temperature time-dependent variation with multiple Davydov states, J. Chem. Phys. 146, 124127 (2017a).
  * Chen  _et al._ [2019] L. Chen, M. F. Gelin, and W. Domcke, Multimode quantum dynamics with multiple Davydov D2 trial states: Application to a 24-dimensional conical intersection model, Journal of Chemical Physics 150, 24101 (2019).
  * Jakučionis  _et al._ [2020] M. Jakučionis, T. Mancal, and D. Abramavičius, Modeling irreversible molecular internal conversion using the time-dependent variational approach with sD2 ansatz, Physical Chemistry Chemical Physics 22, 8952 (2020).
  * Wang  _et al._ [2021] L. Wang, F. Zheng, J. Wang, F. Großmann, and Y. Zhao, Schrödinger-Cat States in Landau-Zener-Stückelberg-Majorana Interferometry: A Multiple Davydov Ansatz Approach, Journal of Physical Chemistry B 125, 3184 (2021).
  * Jakučionis  _et al._ [2022] M. Jakučionis, A. Žukas, and D. Abramavičius, Modeling molecular J and H aggregates using multiple-Davydov D2 ansatz, Physical Chemistry Chemical Physics 24, 17665 (2022).
  * Jakučionis  _et al._ [2023] M. Jakučionis, A. Žukas, and D. Abramavičius, Inspecting molecular aggregate quadratic vibronic coupling effects using squeezed coherent states, Physical Chemistry Chemical Physics 25, 1705 (2023).
  * Sun  _et al._ [2015] K. W. Sun, M. F. Gelin, V. Y. Chernyak, and Y. Zhao, Davydov Ansatz as an efficient tool for the simulation of nonlinear optical response of molecular aggregates, Journal of Chemical Physics 142, 212448 (2015).
  * Zhou  _et al._ [2016b] N. Zhou, L. Chen, Z. Huang, K. Sun, Y. Tanimura, and Y. Zhao, Fast, Accurate Simulation of Polaron Dynamics and Multidimensional Spectroscopy by Multiple Davydov Trial States, Journal of Physical Chemistry A 120, 1562 (2016b).
  * Zhao [2023] Y. Zhao, The hierarchy of Davydov s Ansätze: From guesswork to numerically "exact" many-body wave functions, Journal of Chemical Physics 158, 80901 (2023).
  * Sun  _et al._ [2010] J. Sun, B. Luo, and Y. Zhao, Dynamics of a one-dimensional Holstein polaron with the Davydov ansätze, Physical Review B - Condensed Matter and Materials Physics 82, 014305 (2010).
  * Luo  _et al._ [2010] B. Luo, J. Ye, C. Guan, and Y. Zhao, Validity of time-dependent trial states for the Holstein polaron, Physical Chemistry Chemical Physics 12, 15073 (2010).
  * Jakučionis  _et al._ [2018] M. Jakučionis, V. Chorošajev, and D. Abramavičius, Vibrational damping effects on electronic energy relaxation in molecular aggregates, Chemical Physics 515, 193 (2018).
  * Jakucionis  _et al._ [2022] M. Jakucionis, I. Gaiziunas, J. Sulskus, and D. Abramavicius, Simulation of Ab Initio Optical Absorption Spectrum of $\beta$-Carotene with Fully Resolved S0 and S2 Vibrational Normal Modes, The Journal of Physical Chemistry A 126, 180 (2022).
  * Davydov [1979] A. S. Davydov, Solitons in molecular systems, Physica Scripta 20, 387 (1979).
  * Scott [1991] A. C. Scott, Davydov’s soliton revisited, Physica D: Nonlinear Phenomena 51, 333 (1991).
  * Zhang  _et al._ [1990] W.-M. Zhang, D. H. Feng, and R. Gilmore, Coherent states: Theory and some applications, Rev. Mod. Phys 62, 867 (1990).
  * Kais and Levine [1990] S. Kais and R. D. Levine, Coherent states for the Morse oscillator, Phys. Rev. A 41, 2301 (1990).
  * Moix  _et al._ [2012] J. M. Moix, Y. Zhao, and J. Cao, Equilibrium-reduced density matrix formulation: Influence of noise, disorder, and temperature on localization in excitonic systems, Phys. Rev. B 85, 115412 (2012).
  * Subasi  _et al._ [2012] Y. Subasi, C. H. Fleming, J. M. Taylor, and B. L. Hu, Equilibrium states of open quantum systems in the strong coupling regime, Phys. Rev. E 86, 061132 (2012).
  * Gelzinis and Valkunas [2020] A. Gelzinis and L. Valkunas, Analytical derivation of equilibrium state for open quantum system, The Journal of Chemical Physics 152, 051103 (2020).
  * Zeng and Yao [2022] J. Zeng and Y. Yao, Variational Squeezed Davydov Ansatz for Realistic Chemical Systems with Nonlinear Vibronic Coupling, Journal of Chemical Theory and Computation 18, 1255 (2022).
  * Plenio and Knight [1998] M. B. Plenio and P. L. Knight, The quantum-jump approach to dissipative dynamics in quantum optics, Reviews of Modern Physics 70, 101 (1998).
  * Luoma  _et al._ [2020] K. Luoma, W. T. Strunz, and J. Piilo, Diffusive Limit of Non-Markovian Quantum Jumps, Physical Review Letters 125, 150403 (2020).
  * Kampen [2007] V. N. Kampen,  _Stochastic Processes in Physics and Chemistry_ (Elsevier, 2007).
  * Bertsekas and Tsitsiklis [2008] D. Bertsekas and J. Tsitsiklis, Introduction to probability, in  _Introductory Medical Statistics_ (IOP Publishing Ltd, 2008).
  * Schlosshauer [2007] M. Schlosshauer,  _Decoherence and the Quantum-To-Classical Transition_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007).
  * Glauber [1963] R. J. Glauber, Coherent and incoherent states of the radiation field, Physical Review 131, 2766 (1963).
  * Chorošajev  _et al._ [2016] V. Chorošajev, O. Rancova, and D. Abramavicius, Polaronic effects at finite temperatures in the B850 ring of the LH2 complex, Physical Chemistry Chemical Physics 18, 7966 (2016).
  * Wang  _et al._ [2017b] L. Wang, Y. Fujihashi, L. Chen, and Y. Zhao, Finite-temperature time-dependent variation with multiple Davydov states, The Journal of Chemical Physics 146, 124127 (2017b).
  * Xie  _et al._ [2017] Q. Xie, H. Zhong, M. T. Batchelor, and al, Including temperature in a wavefunction description of the dynamics of the quantum Rabi model, Journal of Physics A: Mathematical and Theoretical 51, 014001 (2017).
  * Werther and Großmann [2020] M. Werther and F. Großmann, Apoptosis of moving nonorthogonal basis functions in many-particle quantum dynamics, Physical Review B 101, 174315 (2020).
  * Balevičius  _et al._ [2015] V. Balevičius, L. Valkunas, and D. Abramavicius, Modeling of ultrafast time-resolved fluorescence applied to a weakly coupled chromophore pair, The Journal of Chemical Physics 143, 074101 (2015).
  * Hestand and Spano [2018] N. J. Hestand and F. C. Spano, Expanded Theory of H- and J-Molecular Aggregates: The Effects of Vibronic Coupling and Intermolecular Charge Transfer, Chemical Reviews 118, 7069 (2018).
  * Kell  _et al._ [2013] A. Kell, X. Feng, M. Reppert, and R. Jankowiak, On the shape of the phonon spectral density in photosynthetic complexes, Journal of Physical Chemistry B 117, 7317 (2013).
  * Abramavicius  _et al._ [2018] D. Abramavicius, V. Chorošajev, and L. Valkunas, Tracing feed-back driven exciton dynamics in molecular aggregates, Physical Chemistry Chemical Physics 20, 21225 (2018).
  * Scarani  _et al._ [2002] V. Scarani, M. Ziman, P. Štelmachovič, N. Gisin, and V. Bužek, Thermalizing Quantum Machines: Dissipation and Entanglement, Physical Review Letters 88, 097905 (2002).
  * Chin  _et al._ [2013] A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio, The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes, Nature Physics 9, 113 (2013).
  * K Mogensen and N Riseth [2018] P. K Mogensen and A. N Riseth, Optim: A mathematical optimization package for Julia, Journal of Open Source Software 3, 615 (2018).
  * Zhan  _et al._ [2009] Z. H. Zhan, J. Zhang, Y. Li, and H. S. Chung, Adaptive particle swarm optimization, IEEE transactions on systems, man, and cybernetics. Part B, Cybernetics : a publication of the IEEE Systems, Man, and Cybernetics Society 39, 1362 (2009).
  * Chorošajev  _et al._ [2017] V. Chorošajev, T. Marčiulionis, and D. Abramavicius, Temporal dynamics of excitonic states with nonlinear electron-vibrational coupling, The Journal of Chemical Physics 147, 074114 (2017).
  * Abramavicius and Valkunas [2003] D. Abramavicius and L. Valkunas, Geminate pair recombination in molecular systems with correlated disorder, Physical Review B 68, 245203 (2003).
  * Eisfeld  _et al._ [2010] A. Eisfeld, S. M. Vlaming, V. A. Malyshev, and J. Knoester, Excitons in molecular aggregates with lévy-type disorder: Anomalous localization and exchange broadening of optical spectra, Physical Review Letters 105, 137402 (2010).
  * Rancova  _et al._ [2017] O. Rancova, M. Jakučionis, L. Valkunas, and D. Abramavicius, Origin of non-Gaussian site energy disorder in molecular aggregates, Chemical Physics Letters 674, 120 (2017).