11institutetext: Vanderbilt University, Nashville TN 37235, USA
11email<EMAIL_ADDRESS>
# Cross-Modality Domain Adaptation for Vestibular Schwannoma and Cochlea
Segmentation

Han Liu  Yubo Fan  Can Cui  Dingjie Su  Andrew McNeil  Benoit M. Dawant

###### Abstract

Automatic methods to segment the vestibular schwannoma (VS) tumors and the
cochlea from magnetic resonance imaging (MRI) are critical to VS treatment
planning. Although supervised methods have achieved satisfactory performance
in VS segmentation, they require full annotations by experts, which is
laborious and time-consuming. In this work, we aim to tackle the VS and
cochlea segmentation problem in an unsupervised domain adaptation setting. Our
proposed method leverages both the image-level domain alignment to minimize
the domain divergence and semi-supervised training to further boost the
performance. Furthermore, we propose to fuse the labels predicted from
multiple models via noisy label correction. Our results on the challenge
validation leaderboard showed that our unsupervised method has achieved
promising VS and cochlea segmentation performance with mean dice score of
0.8261 ± 0.0416; The mean dice value for the tumor is 0.8302 ± 0.0772. This is
comparable to the weakly-supervised based method.

###### Keywords:

Vestibular schwannoma Cochlea Unsupervised domain adaptation

##  1 Introduction

Vestibular schwannoma (VS) is a benign tumor that arises from the Schwann
cells of the vestibular nerve, which connects the brain and the inner ear. To
facilitate the follow-up and treatment planning of VS, automatic methods to
segment the VS tumors and the cochlea from magnetic resonance imaging (MRI)
are proposed [1]. While the most commonly used modality for VS segmentation is
contrast-enhanced T1 (ceT1), high-resolution T2 (hrT2) imaging has been
demonstrated to be a possible alternative with less risk and lower cost.[2].

Supervised segmentation methods have shown the good performance in VS
segmentation [3], but they require to fully annotate image data which may not
be an option in practice. Weakly-supervised methods require less annotation
efforts, such as scribbles and bounding boxes, and sometimes they even achieve
comparable performance to the supervised ones [4]. In this work, we aim at
segmenting the VS tumor and the cochlea in hrT2 without any hrT2 annotations
during training. Specifically, we are provided with a dataset consisting of
ceT1 images and hrT2 images, but only the ceT1 images have the segmentation
labels. We consider the problem as an unsupervised domain adaptation (UDA)
problem. There are mainly two types of methods to tackle the UDA problem,
domain alignment and techniques based on semi-supervised learning (SSL).
Domain alignment focuses on reducing the distribution discrepancy by
optimizing some divergence metric [5, 6] or via adversarial learning [7, 8].
On the other hand, self-training [9], mean teacher [10], and other SSL-basd
techniques also offer competitive performance. In this work, we focus on
exploring methods that combines image-level domain alignment and SSL for UDA.

##  2 Methods

###  2.1 Problem formulation

For an unsupervised domain adaptation problem, we have access to a source
domain $D^{S}=\\{(x^{s}_{i},y^{s}_{i})|i=1,2,\cdots,n_{s}\\}$, and a target
domain $D^{T}=\\{x^{t}_{j}|j=1,2,\cdots,n_{t}\\}$, where $Y^{S}$ and $Y^{T}$
share the same $K$ classes. In our case, source and target domains correspond
to ceT1 and hrT2 respectively and $K=3$ representing background, VS and
cochlea. We aim to train a segmentation network $F_{t}$ that learns the
knowledge from the source domain and is capable to achieve robust and accurate
segmentation performance on the target domain, without accessing the target
domain labels $Y^{T}$.

###  2.2 Image-level Domain Alignment

Image-level domain alignment is a simple but effective method to tackle UDA
problem by reducing the distribution mismatch at the image-level, i.e., pseudo
image synthesis. Here, we propose to train the segmentation model $F_{t}$ with
the pseudo target domain images $\tilde{X}^{T}$, which are generated by
unpaired image-to-image translation. We explored both end-to-end training and
two-stage training. For end-to-end training, we rely on the Contrastive
Unpaired Translation (CUT) [11] as backbone for image synthesis and add an
extra segmentation module $F_{t}$ on top of the synthesized images. This
method will be referred as CutSeg. We select CUT for unpaired image-to-image
translation because it can be trained faster and is less memory-intensive,
allowing more flexibility when adding the 3D CNN-based segmentation module.
During training, we first train the CUT model alone till it achieves
reasonable synthesis performance. Then we train the CutSeg end-to-end with the
CUT subnetwork initialized with the pre-trained weights and the segmentation
module trained from scratch. For two-stage training, we used the CycleGAN [12]
to generate pseudo hrT2 images $\tilde{X}^{T}$. To improve the data diversity,
we trained both 2D and 3D CycleGANs and collected pseudo images from different
epochs. Lastly, we trained a segmentation module $F_{t}$ using
$\tilde{X}^{T}$.

###  2.3 Semi-supervised Training

Though image-level domain alignment can minimize the domain divergence, the
unlabeled target domain images $X^{T}$ are not directly involved in training
the segmentation model $F_{t}$. To overcome this limitation, we propose to
adapt a semi-supervised learning method named Mean Teacher (MT) [13] to make
better use of $X^{T}$. Specifically, a student model along with a teacher
model with the same network architecture are created and both models are
initialized with the best model weights obtained from Section 2.2. In our
semi-supervised setting, the labeled images are the pseudo hrT2 images while
the unlabeled images are the real hrT2 images. During training, the labeled
pseudo images are fed to student model and the segmentation loss $L_{seg}$ is
computed in a supervised manner. For unlabeled images, we first augment the
same image twice with different intensity transformation parameters. The
augmented images are then fed to the student model and the teacher model
separately and a consistency loss $L_{con}$ is computed. Both $L_{seg}$ and
$L_{con}$ are used to update the weights of the student model and the weights
of teacher model are updated as an exponential moving average of the student
weights. As suggested in [13], the teacher prediction is more likely to be
correct at the end of the training and thus the teacher model is taken as our
final $F_{t}$.

###  2.4 Noisy Label Correction as Label Fusion

In this challenge, we have obtained three models that were trained from
different strategies and each model alone has achieved satisfactory result on
the validation leaderboard. The first model is obtained by two-stage training
using the pseudo images from 2D CycleGAN, followed by semi-supervised training
by MT. The second model is fine-tuned based on the first model using the
pseudo images from 3D CycleGAN. The third model is a CutSeg model. Training
details can be found in Section 3.1.

Figure 1: The schematic diagram of our proposed method

Empirically, ensembles tend to yield better predictive performance when there
is a significant diversity among the models. Here, we propose to fuse the
labels from different models by treating the label fusion task as a noisy
label correction problem. We adapted a confident learning method called
cleanlab [14] which provides exact noise estimation and label error finding.
Note that we use cleanlab to directly fuse labels at the inference phase
rather than update the pseudo labels iteratively during training.
Specifically, we first obtain the softmax outputs of two models and convert
one output to an one-hot encoded label mask. The one-hot encoded mask was
considered as ’noisy label’ and was corrected by the softmax outputs from the
other model. Once the labels from the first two models were fused, the fused
labels are treated as noisy labels and fused again with the softmax outputs
from the remaining model. The labels fused from three models are used as our
final predictions.

##  3 Experiments and Results

###  3.1 Data and Implementation.

The dataset was released by the MICCAI challenge CrossMoDA 2021. All images
were obtained on a 32-channel Siemens Avanto 1.5T scanner using a Siemens
single-channel head coil [15]. ceT1 imaging was performed with an MPRAGE
sequence with in-plane resolution of 0.4×0.4 mm, in-plane matrix of 512×512,
and slice thickness of 1.0 to 1.5 mm. For hrT2 images, imaging was performed
with a 3D CISS or FIESTA sequence in-plane resolution of 0.5×0.5 mm, in-plane
matrix of 384×384 or 448×448, and slice thickness of 1.0 to 1.5 mm. The VS and
cochleas were manually segmented in consensus by the treating neurosurgeon and
physicist using both the ceT1 and hrT2 images. We randomly split the images
into 185 and 25 for training and validation respectively. Since the field of
views (FoV) of the source and target domain images vary significantly, we crop
each image into a cubic box, or ROI, using single-atlas registration [16]. The
ROI on the atlas image is manually cropped around the right side of the brain.
To obtain the ROI on the left side, we flip the volume left-to-right before
performing registration.

For preprocessing, in two-stage training, we resample the images to the most
common spacing in the target domain, i.e., (0.46875, 0.468975, 1.5) and
normalize the intensity to [0, 1]. In end-to-end training, we first train a
CutSeg models for 82 epochs with an auxiliary consistency loss, which is a
Mean Absolute Error (MAE) loss between the segmentation result of the real
hrT2 image and the prediction from the two-stage training model. Then the
CutSeg model was fine-tuned on the hrT2 images with in-plane resolution higher
than 0.5 mm. During inference, the fine-tuned CutSeg model was used to make
inference for the testing images with in-plane resolution above 0.5 mm. For
the segmentation module, we adapted the model architecture from [4] and used
dice + cross-entropy loss for training. For post-processing, we first reduce
the false positive VS prediction by removing the isolated components whose
center is greater than the adjacent cochlea center by 15 voxels along z-axis.
Then we take the largest connected components for both VS and cochlea within
each ROI.

For training, we used Adam optimizer with weight decay $10^{-4}$ and batch
size $1$. The learning rates were initialized to $5\times 10^{-4}$, $5\times
10^{-5}$ and $2\times 10^{-4}$ for two-stage training, MT and CutSeg
respectively. The hyperparameters were determined by grid-search within the
range of $10^{-2}$ to $10^{-6}$. The best hyperparameters were selected based
on the segmentation performance on our own validation set. The CNNs were
implemented in PyTorch [17] and MONAI on a Ubuntu desktop with an NVIDIA RTX
2080 Ti GPU. For quantitative evaluation, we measured the Dice score and
average symmetric surface distance (ASSD) between segmentation results and the
ground truth.

###  3.2 Experimental Results

Figure 2: Quantitative results. A to C and D to F show the best and worst VS
segmentation results. The image ID and the corresponding dice score are
displayed.

The following table shows the evaluation metrics of our proposed method on the
validation leaderboard. Due to the page length limit, the ablation studies and
other explored methods are not included in this paper.

Table 1: Quantitative results on validation leaderboard | Dice | ASSD
---|---|---
VS | 0.8302 ± 0.0772 | 0.5686 ± 0.2675
Cochlea | 0.8220 ± 0.0310 | 0.1829 ± 0.0476

##  4 Conclusion

In this work, we exploited the image-level domain alignment and semi-
supervised training to tackle the unsupervised domain adaptation segmentation
problem. According to the validation leaderboard, our unsupervised method has
achieved a segmentation performance that is comparable to the performance of
the weakly-supervised method, demonstrating its effectiveness.

## References

  * [1] Vokurka, Elizabeth A., et al. ”Using Bayesian tissue classification to improve the accuracy of vestibular schwannoma volume and growth measurement.” American journal of neuroradiology 23.3 (2002): 459-467.
  * [2] Coelho, Daniel H., et al. ”MRI surveillance of vestibular schwannomas without contrast enhancement: clinical and economic evaluation.” The Laryngoscope 128.1 (2018): 202-209.
  * [3] Wang, Guotai, et al. ”Automatic segmentation of vestibular schwannoma from T2-weighted MRI by deep spatial attention with hardness-weighted loss.” International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Cham, 2019.
  * [4] Dorent, Reuben, et al. ”Scribble-based Domain Adaptation via Co-segmentation.” International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Cham, 2020.
  * [5] Lee, Chen-Yu, et al. ”Sliced wasserstein discrepancy for unsupervised domain adaptation.” Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2019.
  * [6] Long, Mingsheng, et al. ”Learning transferable features with deep adaptation networks.” International conference on machine learning. PMLR, 2015.
  * [7] Hoffman, Judy, et al. ”Cycada: Cycle-consistent adversarial domain adaptation.” International conference on machine learning. PMLR, 2018.
  * [8] Ganin, Yaroslav, et al. ”Domain-adversarial training of neural networks.” The journal of machine learning research 17.1 (2016): 2096-2030.
  * [9] Zou, Yang, et al. ”Unsupervised domain adaptation for semantic segmentation via class-balanced self-training.” Proceedings of the European conference on computer vision (ECCV). 2018.
  * [10] Perone, Christian S., et al. ”Unsupervised domain adaptation for medical imaging segmentation with self-ensembling.” NeuroImage 194 (2019): 1-11.
  * [11] Park, Taesung, et al. ”Contrastive learning for unpaired image-to-image translation.” European Conference on Computer Vision. Springer, Cham, 2020.
  * [12] Zhu, Jun-Yan, et al. ”Unpaired image-to-image translation using cycle-consistent adversarial networks.” Proceedings of the IEEE international conference on computer vision. 2017.
  * [13] Tarvainen, Antti, and Harri Valpola. ”Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results.” arXiv preprint arXiv:1703.01780 (2017).
  * [14] Northcutt, Curtis, Lu Jiang, and Isaac Chuang. ”Confident learning: Estimating uncertainty in dataset labels.” Journal of Artificial Intelligence Research 70 (2021): 1373-1411.
  * [15] Shapey, J., Kujawa, A., Dorent, R., Wang, G., Dimitriadis, A., Grishchuk, D., Pad-dick, I., Kitchen, N., Bradford, R., Saeed, S.R., Bisdas, S., Ourselin, S., Vercauteren,T.: Segmentation of vestibular schwannoma from mri — an open annotated datasetand baseline algorithm. Scientific Data (2021), in press. Preprint available athttps://doi.org/10.1101/2021.08.04.21261588medRXiv:10.1101/2021.08.04.21261588
  * [16] Avants, Brian B., et al. ”A reproducible evaluation of ANTs similarity metric performance in brain image registration.” Neuroimage 54.3 (2011): 2033-2044.
  * [17] Paszke, Adam, et al. ”Pytorch: An imperative style, high-performance deep learning library.” Advances in neural information processing systems 32 (2019): 8026-8037.


# Finite-temperature vibronic spectra from the split-operator thermofield
coherence dynamics

Zhan Tong Zhang  Jiří Vaníček<EMAIL_ADDRESS>Laboratory of Theoretical
Physical Chemistry, Institut des Sciences et Ingénierie Chimiques, Ecole
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

###### Abstract

The inclusion of temperature effects is important to properly simulate and
interpret experimentally observed vibrationally resolved electronic spectra.
We present a numerically exact approach for evaluating these spectra at finite
temperature using the thermofield coherence dynamics. In this method, which
avoids implementing an algorithm for solving the von Neumann equation for the
coherence, the thermal vibrational ensemble is first mapped to a pure-state
wavepacket in an augmented space, and this wavepacket is then propagated by
solving the standard, zero-temperature Schrödinger equation with the split-
operator Fourier method. We show that the finite-temperature spectra obtained
with the thermofield coherence dynamics in a Morse potential agree exactly
with those computed by Boltzmann-averaging the spectra of individual
vibrational levels. Because the split-operator thermofield dynamics on a full
tensor-product grid is restricted to low-dimensional systems, we briefly
discuss how the accessible dimensionality can be increased by various
techniques developed for the zero-temperature split-operator Fourier method.

##  I Introduction

Vibrationally resolved electronic spectroscopy has proven to be a powerful
tool for probing molecular systems and understanding photochemical processes.
These experiments are typically conducted at ambient or higher temperatures,
which are relevant for understanding biochemical and atmospheric processes as
well as developing new materials.[1, 2, 3, 4] Other than the appearance of hot
bands and the thermal Doppler broadening effects, non-Condon effects can also
noticeably alter spectral features at finite temperatures.[5, 6, 7, 8] These
thermal effects may become more pronounced in two-dimensional spectroscopy
experiments.[9] Including temperature effects in theoretical simulations is
thus essential to accurately interpret and predict the results of vibronic
spectroscopy.[10, 11, 3, 12, 13, 14, 15, 9, 16, 4]

Given the large number of vibrational levels that may be occupied, time-
independent sum-over-state expressions become impractical at higher
temperatures for more than a few vibrational degrees of freedom. Time-
dependent approaches offer a promising alternative, in which the spectrum is
evaluated as the Fourier transform of an appropriate autocorrelation function
without the need to select relevant transitions and compute their Franck-
Condon factors.[17, 11] However, the most widely used quantum mechanical[18,
19, 20, 21, 22] and semiclassical[23, 17, 24, 25, 26, 27, 28] methods were
developed only for pure states and cannot directly handle a thermal ensemble
of states.

Standard approaches to include temperature effects either directly evolve the
thermal density matrix or employ statistical averaging over pure states. An
alternative way to represent quantum dynamics at finite temperature, called
thermofield dynamics,[29, 30] maps the thermal ensemble to a pure state in an
augmented space, which can be propagated by conventional wavepacket methods
for solving the time-dependent Schrödinger equation. Thermofield dynamics has
been applied to problems in electronic structure[31] and molecular quantum
dynamics,[32, 33, 34, 35, 36, 37, 38] including the simulation of vibronic
spectra.[39, 40, 12, 9, 41, 42, 43]

In vibronic spectroscopy applications, the evolution of the thermofield
wavepacket replaces the dynamics of the coherence between two electronic
states. Due to the doubling of dimensions, the thermofield coherence dynamics
obviously benefits from a combination with semiclassical or other approximate
propagation methods that scale favorably with the number of degrees of
freedom. Combining the thermofield dynamics with the extended thawed Gaussian
approximation[44, 45, 46] has enabled the simultaneous description of finite-
temperature, non-Condon,[40] and even anharmonicity effects.[12, 9, 41] Yet,
one would sometimes like to avoid dynamical approximations and obtain the
exact quantum solution.

For small systems with a known potential energy surface, the dynamical Fourier
method may be used to solve the time-dependent Schrödinger equation on a grid
with the split-operator algorithm.[18, 47, 48, 22, 26] Although limited to low
dimensions, the split-operator Fourier method offers a numerically exact
solution, which can validate semiclassical and other more approximate
approaches.

Here, we combine the thermofield coherence dynamics with the split-operator
Fourier method in order to simulate vibronic spectra at finite temperature.
This approach is validated on a one-dimensional anharmonic Morse model by
comparing the spectra computed with the split-operator thermofield coherence
dynamics to the results of the Boltzmann averaging of spectral contributions
from different initial vibrational levels. We also examine the thermofield
wavepacket’s evolution in the augmented configuration space consisting of both
the physical and fictitious degrees of freedom.

##  II Theory

Let us consider a molecular system with $D$ nuclear degrees of freedom and two
uncoupled adiabatic electronic states: the ground state $|g\rangle$ and the
excited state $|e\rangle$. The dimensionless spectrum for vibronic transitions
from $|g\rangle$ to $|e\rangle$ is obtained as the Fourier transform

$\sigma(\omega)=\operatorname{Re}\int_{0}^{\infty}C(t)e^{i\omega t}\,dt,$ (1)

of the dipole-dipole correlation function

$C(t)=\operatorname{Tr}(\hat{\mu}^{\dagger}e^{-i\hat{H}_{e}t/\hbar}\hat{\mu}\hat{\rho}e^{i\hat{H}_{g}t/\hbar}),$
(2)

where $\hat{\mu}$ represents the transition dipole moment operator,
$\hat{H}_{g}$ and $\hat{H}_{e}$ are the vibrational Hamiltonians associated
with the ground and excited electronic states, respectively, and $\hat{\rho}$
is the vibrational density operator of the ground electronic state.[17, 22]

At finite temperature $T$, the density of occupied vibrational levels in the
initial population follows the Boltzmann distribution

$\hat{\rho}=e^{-\beta\hat{H}_{g}}/\operatorname{Tr}(e^{-\beta\hat{H}_{g}}),$
(3)

with the inverse temperature $\beta=1/(k_{B}T)$. We can then rewrite $C(t)$ as

$C(t)=\sum_{n}p_{n}e^{i\omega_{n,g}t}\langle
n,g|\hat{\mu}^{\dagger}e^{-i\hat{H}_{e}t/\hbar}\hat{\mu}|n,g\rangle,$ (4)

where $p_{n}=e^{-\beta\hbar\omega_{n,g}}/\sum_{k}e^{-\beta\hbar\omega_{k,g}}$
is the Boltzmann probability, $|n,g\rangle$ denotes the $n$-th vibrational
eigenfunction of $\hat{H}_{g}$, and $\omega_{n,g}$ is the corresponding
vibrational frequency.[11] The thermal spectrum can be obtained from

$C(t)=\sum_{n}p_{n}e^{i\omega_{n,g}t}\langle\varphi_{n}(0)|\varphi_{n}(t)\rangle$
(5)

by propagating the initial vibrational wavepackets
$|\varphi_{n}(0)\rangle=\hat{\mu}|n,g\rangle$ individually as
$|\varphi_{n}(t)\rangle=e^{-i\hat{H}_{e}t/\hbar}|\varphi_{n}(0)\rangle$. In
high-dimensional systems, propagating wavepackets associated with all occupied
vibrational levels becomes impractical, and a statistical sampling approach is
needed.[49, 50, 51, 15, 16]

Using the thermofield coherence dynamics, the thermal ensemble is instead
mapped to a pure state in a doubled configuration space. The thermal dipole-
dipole correlation function (2) can be rewritten exactly[12] as a wavepacket
autocorrelation function

$C(t)=\langle\bar{\psi}|e^{-i\hat{\bar{H}}t/\hbar}|\bar{\psi}\rangle$ (6)

of an initial thermofield wavepacket, represented in the augmented thermofield
space (denoted by a bar) as

$\bar{\psi}(\bar{q})=\langle q|\hat{\mu}\hat{\rho}^{1/2}|q^{\prime}\rangle,$
(7)

where $\bar{q}=(q,q^{\prime})^{T}$ is a $2D$-dimensional coordinate vector in
the augmented configuration space and $|q^{\prime}\rangle$ represents the
position states in the “fictitious” Hilbert space (denoted by a prime). The
wavepacket is propagated according to the time-dependent Schrödinger equation
with the “augmented” Hamiltonian $\hat{\bar{H}}$ expressed in the position
representation as

$\bar{H}(\bar{q})=H_{e}(q)-H_{g}(q^{\prime}),$ (8)

with a $2D\times 2D$ mass matrix

$\bar{m}=\begin{pmatrix}m&0\\\ 0&-m\end{pmatrix}$ (9)

replacing the $D\times D$ mass matrix $m$. In contrast to standard thermofield
dynamics, where the two Hamiltonians in the right-hand side of Eq. (8) are
identical, in vibronic spectroscopy the thermofield Hamiltonian (8) depends on
the two Hamiltonians associated with the two electronic states involved in the
transition. In essence, the thermofield wavepacket encodes the coherence
between the dynamics on the two potential energy surfaces.

If the ground surface

$V_{g}(q^{\prime})=v_{0}+(q^{\prime}-q_{\text{ref}}^{\prime})^{T}\cdot\kappa\cdot(q^{\prime}-q_{\text{ref}}^{\prime})/2$
(10)

is harmonic, the initial thermofield wavepacket[12]

$\psi(\bar{q})=e^{(i/\hbar)[(\bar{q}-\bar{q}_{0})^{T}\cdot\bar{A}_{0}\cdot(\bar{q}-\bar{q}_{0})/2+\bar{p}_{0}^{T}\cdot(\bar{q}-\bar{q}_{0})+\bar{\gamma}_{0}]}$
(11)

is a Gaussian centered at

$\bar{q}_{0}=\begin{pmatrix}q_{\text{ref}}^{\prime}\\\
q_{\text{ref}}^{\prime}\end{pmatrix},\qquad\bar{p}_{0}=\begin{pmatrix}0\\\
0\end{pmatrix}$ (12)

in the phase space. This Gaussian has a temperature-dependent width matrix[12]

$\bar{A}_{0}=i\begin{pmatrix}\mathrm{A}&\mathrm{B}\\\
\mathrm{B}&\mathrm{A}\end{pmatrix},$ (13)

where $\mathrm{A}$ and $\mathrm{B}$ are $D\times D$ block matrices

$\displaystyle\mathrm{A}=m^{1/2}\cdot\Omega\cdot\operatorname{coth}(\beta\hbar\Omega/2)\cdot
m^{1/2},$ (14)
$\displaystyle\mathrm{B}=-m^{1/2}\cdot\Omega\cdot\sinh(\beta\hbar\Omega/2)^{-1}\cdot
m^{1/2},$ (15)

with $\Omega=\sqrt{m^{-1/2}\cdot\kappa\cdot m^{1/2}}$. The phase factor,

$\bar{\gamma}_{0}=(-i\hbar/2)\ln[\operatorname{det}(m\cdot\Omega/\pi\hbar)],$
(16)

serves to ensure normalization.

The spectrum at finite temperature can now be obtained by applying standard
techniques for solving the time-dependent Schrödinger equation to the
thermofield wavepacket with a doubled number of coordinates. Although the
increase in dimensionality implies that an exact thermofield solution does not
save computational time compared to the direct propagation of the density
matrix, the wavepacket approach permits the use of zero-temperature numerical
methods and avoids implementing methods for solving the von Neumann equation.

Among various numerical methods for propagating wavepackets, we choose to
combine the thermofield dynamics with the second-order split-operator
algorithm, a popular choice for separable Hamiltonians, because it is
explicit, numerically stable, easy to implement and because it preserves
several geometric properties of the exact evolution operator, such as the
unitarity, symplecticity, and time reversibility.

##  III Numerical examples

In the following, we use natural units with $m=\hbar=k_{B}=\hat{\mu}=1$. As in
Ref. 12, the excited-state surface is a one-dimensional Morse potential

$V_{e}(q)=V_{e}(q_{\text{ref}})+\frac{\omega_{e}}{4\chi}\left[1-e^{-\sqrt{2m\omega_{e}\chi}(q-q_{\text{ref}})}\right]^{2},$
(17)

with the equilibrium position $q_{\text{ref}}=1.5$, the minimum energy
$V_{e}(q_{\text{ref}})=10$, the fundamental frequency $\omega_{e}=0.9$ at
$q_{\text{ref}}$, and the anharmonicity measure $\chi=0.02$. The initial,
ground-state potential energy surface is harmonic and given by Eq. (10) with
$q_{\text{ref}}^{\prime}=0$ and $\kappa=1$.

Spectra are computed either by Boltzmann averaging [Eq. (4)] or with
thermofield dynamics [Eq. (6)] at scaled temperatures
$T_{\omega}=k_{B}T/(\hbar\omega_{g})=0$, $0.5$, $1$, and $2$, where
$\omega_{g}=\sqrt{\kappa/m}$ is the ground-state vibrational frequency. Here,
as a demonstration, we apply the standard second-order TVT algorithm on a grid
with 256 points for each degree of freedom. The wavepacket is propagated for a
total time of 1000 with a time step of 0.1. A Gaussian damping function with a
half-width at half-maximum of 15 is applied to the autocorrelation functions
before the spectra are evaluated with Eq. (1).

Figure 1: Comparison between spectra obtained via thermofield coherence
dynamics (solid black line) and Boltzmann averaging (red dashed line) at four
different scaled temperatures $T_{\omega}$. Figure 2: Contributions from
different vibrational levels, scaled by their Boltzmann probability of
occupation, to the Boltzmann-averaged spectrum at $T_{\omega}=1$.

As shown in Fig. 1, the thermofield split-operator Fourier method gives the
same results as the Boltzmann-weighted average. As the temperature increases,
the spectral range broadens and additional peaks (the “hot bands”) appear. In
Ref. 12, a comparison was made between the thermal spectra obtained using
thermofield thawed Gaussian approximation and “exact” spectra obtained by
computing individual Franck-Condon factors by numerical integration. With the
time-dependent split-operator approach, we avoid the need to preselect the
relevant transitions and can give numerically exact results for arbitrary
global potentials, even if their vibrational eigenfunctions do not have
analytical expressions. This approach thus has the potential to validate
semiclassical methods in a wider range of situations, including those with a
fitted global potential energy surface.

In Fig. 2, we show the decomposition of the thermal spectrum at $T_{\omega}=1$
into the contributions from different initial vibrational levels. Only spectra
from levels $n=0$, $1$, $2$, and $3$ are shown since the contributions from
$n>3$ are negligible here (the Boltzmann weight is 0.0002 for $n=4$). As
expected, the higher vibrational levels are the source of hot bands and expand
the spectral range. Whereas the Boltzmann averaging requires propagating each
initial vibrational level separately, a single thermofield wavepacket
trajectory takes into account all thermal contributions.

With a grid-based method, we can easily visualize the wavepacket and its
evolution. Figure 3 shows the initial thermofield wavepacket at different
temperatures in the augmented coordinate space. As the temperature increases,
the off-diagonal block $\mathrm{B}$ of the initial width matrix (13) increases
and the initial thermofield wavepacket rotates and stretches as a result.

Figure 3: Initial thermofield wavepackets at different scaled temperatures
$T_{\omega}$.

In Fig. 4, we show the evolution of the wavepacket at $T_{\omega}=1$. Whereas
the extent of the wavepacket over the fictitious degree of freedom
$q^{\prime}$ remains approximately constant, the initial excitation of the
coordinate $q$ on the anharmonic excited-state surface leads to significant
evolution of the wavepacket in the physical dimension.

Figure 4: Evolution of the thermofield wavepacket (for $T_{\omega}=1$) in the
augmented configuration space at different times $t$; $q$ is the physical
degree of freedom, whereas $q^{\prime}$ is fictitious.

##  IV Discussion and conclusion

We have described a very simple method for computing vibronic spectra at
finite temperature. This method applies the Fourier split-operator propagation
algorithm to the thermofield wavepacket representation of the electronic
coherence, avoiding both Boltzmann averaging and solving the von Neumann
equation. Instead, we solve the zero-temperature Schrödinger equation in an
augmented space. The method yields quantum mechanically exact spectra of
arbitrary global potentials, making it a valuable tool for validating
semiclassical and other approximate techniques for finite-temperature spectra
calculations. If highly accurate results are desired, the simple approach
described here can be made both more accurate and efficient by implementing
high-order geometric integrators obtained by symmetric composition of the
elementary second-order algorithm, in the same way as for the zero-temperature
split-operator calculations.[52, 53, 54, 55]

The method presented here is based on the full tensor-product grid and,
therefore, is limited to systems with very few degrees of freedom. To extend
its applicability to higher-dimensional systems, one can combine the split-
operator thermofield dynamics with adaptive moving grids,[56, 57, 58]
matching-pursuit coherent-state basis,[59] tensor-train representation,[60] or
any other technique designed for the zero-temperature case. In fact, tensor-
train representations have already been used extensively[32, 36] for the
thermofield dynamics, although not in conjunction with the split-operator
Fourier method. Here, we have intentionally restrained ourselves to the
simplest, full-grid method, because we did not want to obscure the main idea
by combining it with various strategies for increasing the accessible
dimensionality. We believe that even the full-grid version can be interesting
for researchers who would like to quickly check the thermal effects on their
zero-temperature spectra of low-dimensional systems without implementing new
codes.

The method is not restricted to linear spectra. In the same way as the
Gaussian,[9] multiconfigurational Ehrenfest,[42] or optimized mean-
trajectory[43] semiclassical thermofield dynamics, the method can be applied
to more sophisticated experimental techniques, such as pump-probe and two-
dimensional spectroscopies. Beyond vibronic spectra, this approach may also
find applications in other finite-temperature dynamics problems such as the
calculation of the rates of internal conversion.[27, 37, 61, 62]

###### Acknowledgements.

The authors thank Tomislav Begušić for discussions and acknowledge the
financial support from the EPFL.

## Author declarations

### Conflict of interest

The authors have no conflicts to disclose.

## Data availability

The data that support the findings of this study are available within the
article.

## References

  * Peterman  _et al._ [1997] E. J. G. Peterman, T. o. Pullerits, R. Van Grondelle, and H. Van Amerongen, J. Phys. Chem. B 101, 4448 (1997).
  * Gaiduk  _et al._ [2010] A. Gaiduk, M. Yorulmaz, P. V. Ruijgrok, and M. Orrit, Science 330, 353 (2010).
  * Qian  _et al._ [2020] Y. Qian, X. Li, A. R. Harutyunyan, G. Chen, Y. Rao, and H. Chen, J. Phys. Chem. A 124, 9156 (2020).
  * Ndengué  _et al._ [2023] S. Ndengué, E. Quintas-Sánchez, R. Dawes, C. C. Blackstone, and D. L. Osborn, J. Phys. Chem. A 127, 6051 (2023).
  * Orlandi [1973] G. Orlandi, Chem. Phys. Lett. 19, 459 (1973).
  * Tanimura and Mukamel [1993] Y. Tanimura and S. Mukamel, J. Opt. Soc. Am. B 10, 2263 (1993).
  * Kundu  _et al._ [2022] S. Kundu, P. P. Roy, G. R. Fleming, and N. Makri, J. Phys. Chem. B 126, 2899 (2022).
  * Roy  _et al._ [2022] P. P. Roy, S. Kundu, N. Makri, and G. R. Fleming, J. Phys. Chem. Lett. 13, 7413 (2022).
  * Begušić and Vaníček [2021a] T. Begušić and J. Vaníček, J. Phys. Chem. Lett. 12, 2997 (2021a).
  * Zucchelli  _et al._ [2002] G. Zucchelli, R. C. Jennings, F. M. Garlaschi, G. Cinque, R. Bassi, and O. Cremonesi, Biophys. J. 82, 378 (2002).
  * Baiardi, Bloino, and Barone [2013] A. Baiardi, J. Bloino, and V. Barone, J. Chem. Theory Comput. 9, 4097 (2013).
  * Begušić and Vaníček [2020] T. Begušić and J. Vaníček, J. Chem. Phys. 153, 024105 (2020).
  * Xun, Deng, and He [2020] J. Xun, J. Deng, and R. He, Chem. Phys. Lett. 759, 138045 (2020).
  * Sršeň  _et al._ [2020] v. Sršeň, J. Sita, P. Slavíček, V. Ladányi, and D. Heger, J. Chem. Theory Comput. 16, 6428 (2020).
  * Wernbacher and González [2021] A. M. Wernbacher and L. González, Phys. Chem. Chem. Phys. 23, 17724 (2021).
  * Prlj  _et al._ [2022] A. Prlj, E. Marsili, L. Hutton, D. Hollas, D. Shchepanovska, D. R. Glowacki, P. Slavíček, and B. F. E. Curchod, ACS Earth Space Chem. 6, 207 (2022).
  * Heller [1981a] E. J. Heller, Acc. Chem. Res. 14, 368 (1981a).
  * Feit, Fleck, and Steiger [1982] M. D. Feit, J. A. Fleck, Jr., and A. Steiger, J. Comp. Phys. 47, 412 (1982).
  * Meyer, Manthe, and Cederbaum [1990] H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990).
  * Ben-Nun and Martínez [1999] M. Ben-Nun and T. J. Martínez, J. Phys. Chem. A 103, 10517 (1999).
  * Worth, Robb, and Burghardt [2004] G. A. Worth, M. A. Robb, and I. Burghardt, Faraday Discuss. 127, 307 (2004).
  * Tannor [2007] D. J. Tannor,  _Introduction to Quantum Mechanics: A Time-Dependent Perspective_ (University Science Books, Sausalito, 2007).
  * Heller [1981b] E. J. Heller, J. Chem. Phys. 75, 2923 (1981b).
  * Herman and Kluk [1984] M. F. Herman and E. Kluk, Chem. Phys. 91, 27 (1984).
  * Hagedorn [1998] G. A. Hagedorn, Ann. Phys. (NY) 269, 77 (1998).
  * Lubich [2008] C. Lubich,  _From Quantum to Classical Molecular Dynamics: Reduced Models and Numerical Analysis_ , 12th ed. (European Mathematical Society, Zürich, 2008).
  * Heller [2018] E. J. Heller,  _The semiclassical way to dynamics and spectroscopy_ (Princeton University Press, Princeton, NJ, 2018).
  * Bonfanti  _et al._ [2018] M. Bonfanti, J. Petersen, P. Eisenbrandt, I. Burghardt, and E. Pollak, J. Chem. Theory Comput. 14, 5310 (2018).
  * Suzuki [1985] M. Suzuki, J. Phys. Soc. Jap. 54, 4483 (1985).
  * Takahashi and Umezawa [1996] Y. Takahashi and H. Umezawa, Int. J. Mod. Phys. B 10, 1755 (1996).
  * Harsha, Henderson, and Scuseria [2019] G. Harsha, T. M. Henderson, and G. E. Scuseria, J. Chem. Theory Comput. 15, 6127 (2019).
  * Borrelli and Gelin [2016] R. Borrelli and M. F. Gelin, J. Chem. Phys. 145, 224101 (2016).
  * Borrelli and Gelin [2017] R. Borrelli and M. F. Gelin, Sci. Rep. 7, 9127 (2017).
  * Gelin and Borrelli [2017] M. F. Gelin and R. Borrelli, Annalen der Physik 529, 1700200 (2017).
  * Chen and Zhao [2017] L. Chen and Y. Zhao, J. Chem. Phys. 147, 214102 (2017).
  * Borrelli and Gelin [2021] R. Borrelli and M. F. Gelin, WIREs Comput. Mol. Sci. 11, e1539 (2021).
  * Fischer and Saalfrank [2021] E. W. Fischer and P. Saalfrank, J. Chem. Phys. 155, 134109 (2021).
  * Gelin and Borrelli [2023] M. F. Gelin and R. Borrelli, J. Chem. Theory Comput. 19, 6402 (2023).
  * Reddy and Prasad [2015] C. S. Reddy and M. D. Prasad, Mol. Phys. 113, 3023 (2015).
  * Reddy and Prasad [2016] C. S. Reddy and M. D. Prasad, J. Phys. Chem. A 120, 2583 (2016).
  * Begušić and Vaníček [2021b] T. Begušić and J. Vaníček, Chimia 75, 261 (2021b).
  * Chen  _et al._ [2021] L. Chen, R. Borrelli, D. V. Shalashilin, Y. Zhao, and M. F. Gelin, J. Chem. Theory Comput. 17, 4359 (2021).
  * Polley and Loring [2022] K. Polley and R. F. Loring, J. Chem. Phys. 156, 124108 (2022).
  * Lee and Heller [1982] S.-Y. Lee and E. J. Heller, J. Chem. Phys. 76, 3035 (1982).
  * Patoz, Begušić, and Vaníček [2018] A. Patoz, T. Begušić, and J. Vaníček, J. Phys. Chem. Lett. 9, 2367 (2018).
  * Prlj  _et al._ [2020] A. Prlj, T. Begušić, Z. T. Zhang, G. C. Fish, M. Wehrle, T. Zimmermann, S. Choi, J. Roulet, J.-E. Moser, and J. Vaníček, J. Chem. Theory Comput. 16, 2617 (2020).
  * Kosloff and Kosloff [1983a] D. Kosloff and R. Kosloff, J. Comp. Phys. 52, 35 (1983a).
  * Kosloff and Kosloff [1983b] R. Kosloff and D. Kosloff, J. Chem. Phys. 79, 1823 (1983b).
  * Manthe and Huarte-Larrañaga [2001] U. Manthe and F. Huarte-Larrañaga, Chem. Phys. Lett. 349, 321 (2001).
  * Meyer, Gatti, and Worth [2009] H.-D. Meyer, F. Gatti, and G. A. Worth, eds.,  _Multidimensional Quantum Dynamics: MCTDH Theory and Applications_ , 1st ed. (Wiley-VCH, Weinheim, 2009).
  * Barbatti and Sen [2016] M. Barbatti and K. Sen, Int. J. Quant. Chem. 116, 762 (2016).
  * Wehrle, Šulc, and Vaníček [2011] M. Wehrle, M. Šulc, and J. Vaníček, Chimia 65, 334 (2011).
  * Choi and Vaníček [2019] S. Choi and J. Vaníček, J. Chem. Phys. 150, 204112 (2019).
  * Roulet, Choi, and Vaníček [2019] J. Roulet, S. Choi, and J. Vaníček, J. Chem. Phys. 150, 204113 (2019).
  * Roulet and Vaníček [2021] J. Roulet and J. Vaníček, J. Chem. Phys. 155, 204109 (2021).
  * Thompson, Warsi, and Mastin [1997] J. F. Thompson, Z. U. Warsi, and C. W. Mastin,  _Numerical grid generation: foundations and applications_ , Vol. 45 (North-holland Amsterdam, 1997) Chap. 11.
  * Pettey and Wyatt [2006] L. R. Pettey and R. E. Wyatt, Chem. Phys. Lett. 424, 443 (2006).
  * Choi and Vaníček [2019] S. Choi and J. Vaníček, J. Chem. Phys. 151, 234102 (2019).
  * Chen and Batista [2006] X. Chen and V. S. Batista, J. Chem. Phys. 125, 124313 (2006).
  * Greene and Batista [2017] S. M. Greene and V. S. Batista, J. Chem. Theory Comput. 13, 4034 (2017).
  * Wang, Ren, and Shuai [2021] Y. Wang, J. Ren, and Z. Shuai, J. Chem. Phys. 154, 214109 (2021).
  * Wenzel and Mitric [2023] M. Wenzel and R. Mitric, J. Chem. Phys. 158, 034105 (2023).