# Multi-level Encoder-Decoder Architectures for Image Restoration

Indra Deep Mastan and Shanmuganathan Raman
Indian Institute of Technology Gandhinagar
Gandhinagar, Gujarat, India
{indra.mastan<EMAIL_ADDRESS>
###### Abstract

Many real-world solutions for image restoration are learning-free and based on
handcrafted image priors such as self-similarity. Recently, deep-learning
methods that use training data, have achieved state-of-the-art results in
various image restoration tasks (e.g., super-resolution and inpainting).
Ulyanov et al. bridge the gap between these two families of methods in [29].
They have shown that learning-free methods perform close to the state-of-the-
art learning-based methods ($\approx$ 1 PSNR). Their approach benefits from
the encoder-decoder network ($ed$).

In this paper, we propose a framework based on the multi-level extensions of
the encoder-decoder network ($med$) to investigate interesting aspects of the
relationship between image restoration and network construction independent of
learning. Our framework allows various network structures by modifying the
following network components: skip links, cascading of the network input into
intermediate layers, a composition of the encoder-decoder subnetworks, and
network depth. These handcrafted network structures illustrate how the
construction of untrained networks influence the following image restoration
tasks: denoising, super-resolution, and inpainting. We also demonstrate image
reconstruction using flash and no-flash image pairs. We provide performance
comparisons with the state-of-the-art methods for all the restoration tasks
above. The additional results are available on the project page.

##  1 Introduction

Image restoration is an ill-posed problem which aims to recover an image given
its corrupted observation (e.g., denoising [39, 5, 31], super-resolution [14,
23, 2], and inpainting [36, 34, 33]). Corruption may occur due to noise,
camera shake, and due to the fact that the picture was taken in rain or
underwater [16]. Image restoration methods could be mainly classified into two
types - _traditional_ methods and _deep-learning_ (DL) methods. Traditional
methods include spatial filtering methods (e.g., bilateral filters [28], non-
local means [4]), wavelet transform based methods [6], and dictionary learning
and sparse coding [17, 37]. DL methods generally include a neural network to
learn image prior from the training samples (learning-based111The learning
refers to training the network on the collection of images and learning-free
refers to the methods which do not use training data.) for restoration, where
the training samples contain paired examples of corrupted and high-quality
images.

Traditional methods are generally faster and comparatively less cumbersome to
implement, e.g., filtering approaches [8]. Whereas DL methods could be tricky
to implement. For example, methods based on adversarial loss require training
of two separate networks, namely a generator and a discriminator [14].
Moreover, DL methods output photo-realistic images with finer details of
features due to the image prior being captured by feature learning on a
collection of images [14, 38, 3].

Representation learning from images gives insight into the image statistics
captured by the network. The main idea is to perform various image restoration
tasks to learn a better image prior [12]. However, it is focused on the
learning-based setting [1]. There are fewer studies that directly investigate
the image prior captured by the neural network without using any training
datasets. Ulyanov et al. first conducted the studies to achieve image
restoration without using a training sample (learning-free) [29]. This paper
focuses on the research thread mentioned above. Our work combine the ideas of
traditional methods and the DL approaches similar to [10, 25, 15, 35, 29].

Our ablation study shows how the structure of the untrained network influences
the quality of image restoration achieved by them. For example, inpainting of
a large missing region is qualitatively better-achieved using an encoder-
decoder network _without_ skip connections, whereas super-resolution is
better-achieved _with_ skip links (Fig. 11 and Fig. 12).

We have performed extensive experiments on various handcrafted network
architectures obtained by modifying the network components. We focus on the
following network components: depth of the network, skip connections,
cascading of the network input into intermediate layers (cascade), and
composition of the encoder-decoder subnetworks (composition). We show how each
of the above network components affects image restoration. For example, we
show how the performance of denoising gets affected when we increase the depth
of the network (Fig. 5).

We have formulated a framework called multi-level encoder decoder ($med$) that
models various handcrafted network architectures. An instance from our
framework $med$ is a composition of three encoder-decoder networks (Fig. 2).
The multi-level extension of encoder-decoder is motivated to exploit the re-
occurrence of an image patch at different resolutions. We show our analysis
using six different network instances of $med$ (Table 1). These handcrafted
network architectures help us develop insight into how the network
construction influences image restoration (Fig. 5, Fig. 8, Fig. 9, Fig. 11,
and Fig. 12). The key idea is to iteratively minimize the loss between the
network output and the corrupted image to implicitly capture the image prior
in the network.

There is an inherent contrast in our objectives. On the one hand, we aim to
experiment with various high capacity networks to show the relation between
image restoration and network construction. The higher depth allows more
network components and various network structures for the analysis of the
image prior. On the other hand, the high capacity network should not
negatively influence the quality of image restoration. This is due to the fact
that the higher depth network suffers from the _vanishing gradients_ problem
[18, 26]. One option is to use skip links to propagate the gradients and feed
the image features from the intermediate layers to the last layers of the
network [18]. Our main contributions are summarized as follows.

  * •

To the best of our knowledge, this is the first study of a multi-level
encoder-decoder framework ($med$) designed to illustrate the relationship
between image restoration and network construction, independent of training
data and using DL. The $med$ framework allows analysis of the deep prior by
using four networks components (depth, skip connections, composition, cascade)
whereas DIP [29] includes the investigation based on the two network
components (depth and skip connections). The $med$ framework provides a more
rigorous evaluation of the usefulness of skip connections compared to [29].

  * •

We also perform various image restoration tasks to show the quality of the
image prior captured by the multi-level network architectures. We have
achieved results comparable to the state-of-the-art methods for denoising,
super-resolution, and inpainting with $x\%$ pixels drop despite experimenting
with various high-capacity networks. We also observe a better flash no-flash
based image construction when compared to [29].

##  2 Related work

Image restoration aims to recover a good quality image from a corrupted
observation. It is a useful preprocessing step for other problems, e.g.,
classification [30]. Mao et al. have shown image restoration using an $ed$
network with symmetric skip links between the layers of encoder and decoder
[18]. There are various proposals for the loss functions for the image
restoration tasks, e.g., adversarial loss [14], perceptual loss [14], or
contextual loss [20, 19]. In addition, Chang et al. have proposed a single
generic network for various image restoration tasks [22]. However, the
drawback to this line of work is that the restoration output could be biased
toward a particular training dataset.

Ulyanov et al. showed that a randomly-initialized $ed$ network works as a
hand-crafted prior for restoring images without training data [29]. Motivated
by their approach, our learning-free framework only uses the handcrafted
structure of the network for image restoration. However, unlike [29], we
explore how the network components directly influence various image
restoration tasks.

##  3 Multi-level Encoder-Decoder Framework

In this section, we explain the multi-level encoder-decoder framework ($med$)
and its major components. We shall also discuss an example construction of a
multi-level encoder-decoder network and then provide a classification of the
networks useful for our experiments.

The $med$ is one of the general class of networks, where each network is a
composition of encoder-decoder blocks as subnetworks. We address $med$ as a
network $\mathcal{F}$ for devising a simpler explanation. The $med$ network
$\mathcal{F}$ is a composition of two subnetworks, namely a generator $G$ and
an enhancer $E$. The image restoration network $\mathcal{F}$ is defined in Eq.
1.

$\mathcal{F}=E\circ G$ (1)

Here, the generator and the enhancer are either an encoder-decoder network
($ed$) or a composition of $ed$ networks. The encoders determine the abstract
representation of the image features, which are used by the decoder for the
reconstruction of the image. The composition of networks allows multiple sub-
networks to learn image features from the down-sampled versions of the
corrupted image. This would enforce the output of the generator to be
consistent across the multiple scales of the target image222Target image
refers to the high-quality image whose corrupted observation $\hat{I}$ is
given for restoration. to improvise the quality of the image restoration.

The multi-level encoder-decoder framework is motivated to model various
network architectures by modifying the network components described in
Subsection 3.1. For example, let’s suppose the generator is a depth-$k$ $ed$
network. There are five network configurations obtained by modifying the skip
connections, namely, Intra-skip, Inter-skip encoder-encoder, Inter-skip
decoder-encoder, No-skip, and Full-skip connections333In the supplementary
material we have provided the details of different types of skip connections..
There are two network configurations based on the cascading of the network
input, i.e., network with cascade or network without cascade. There could be
$(k-1)$ different generator-enhancer compositions for a depth-$k$ generator
network. We do not consider depth-$k$ enhancer to reduce the model capacity.
Finally, given a depth-$k$ $ed$ network as the generator, the $med$ framework
will allow $1\times 5\times 2\times(k-1)=10(k-1)$ different network
structures. On the other hand, [29] will allow only two different network
configurations (network with skip connections and without skip connections).
Therefore, the generalization $med$ provides various networks to analyze the
effects of network components on the quality of the image restoration.
Technically, the $med$ is a general framework to explore the nature of the
mapping between the network parameter space and the natural image space.

###  3.1 Network Components

We focus on the following components to show how the network structure affects
the image restoration output. (a) skip connections, (b) depth of the network,
(c) cascading of the network input into the intermediate layers, and (d)
composition of two $ed$ networks. We describe each of these components as
follows.

(a) Skip connections. The skip link between the layers $L_{i}$ and $L_{j}$,
where $i$ and $j$ are the indices of the network layers with $i<j$, is made by
concatenating the output of the layer $L_{j-1}$ with the output of the layer
$L_{i}$ and then feeding into the layer $L_{j}$. We have provided the detailed
classification of the skip connections in supplementary material. In Fig. 1(a)
and Fig. 1(b), we have pictorially shown useful skip link configurations for
the paper.

(b) Depth of the network. It is measured by the number of layers present in
the network. Higher depth networks capture finer feature details. However, a
very high depth could negatively influence the performance (Fig. 4). There are
two ways to increase network depth. First, by introducing a new layer into the
encoder-decoder ($ed$) network. Second, by performing a composition of the two
$ed$ networks.

(c) Cascading of network input (cascade). It is a procedure to successively
down-sample the network input and then feed it into the intermediate layers of
the network. Formally, to provide the network input at the intermediate layer
$L$, we resize the network input and then concatenate it with the layer $L-1$.
Next, we feed the resulting tensor into the layer $L$. Cascading of network
inputs was also utilized by Chen et al. [7]. We use it to provide the image
features into the enhancer network (Fig. 1(c)).

(d) Composition of $ed$ networks (composition). The composition of two
encoder-decoder networks is achieved by feeding the output of the first $ed$
network into the second $ed$ network. The composition of two $ed$ networks
increases the network depth and the number of skip connections. The main
objective of performing the network composition is to learn image features
from the downsampled versions of the corrupted image.

(a) Intra-skip.

(b) Full-skip.

(c) Cascading of network input.

Figure 1: Network components. Layers of the encoder are in red and layers of
the decoder are in blue. (a) Intra-skip: the skip connections within EDS
network. (b) Full-skip: both the Intra-skip connections and Inter-skip
connections are present. (c) Cascading of the network input.

###  3.2 Multi-level Encoder-Decoder Network

Here, we give an example construction of $med$ network $\mathcal{F}$. It is a
three-level $ed$ network where the generator is the first $ed$ and the
enhancer is a composition of the other two $ed$ (Fig. 2).

$\mathcal{F}=E^{2}\circ E^{1}\circ G$ (2)

In Eq. 2, the subnetwork $G$ is the generator and subnetwork $E^{2}\circ
E^{1}$ is the enhancer $E$.

The networks $G$, $E^{1}$, and $E^{2}$ are defined as follows.
$G:\mathbb{R}^{m\times n\times c}\rightarrow\mathbb{R}^{m\times n\times c}$,
$E^{1}:\mathbb{R}^{\frac{m}{2}\times\frac{n}{2}\times
c}\rightarrow\mathbb{R}^{\frac{m}{2}\times\frac{n}{2}\times c}$, and
$E^{2}:\mathbb{R}^{\frac{m}{4}\times\frac{n}{4}\times
c}\rightarrow\mathbb{R}^{\frac{m}{4}\times\frac{n}{4}\times c}$.

Here, $c$ is the number of channels (c is 3 for RGB images). The generator $G$
operates at $2\times$ the resolution of $E^{1}$ and $4\times$ the resolution
of $E^{2}$. A resize operator $R$ is used to down-sample the output of $G$ to
feed into $E^{1}$ and down-sample the output of $E^{1}$ to feed into $E^{2}$.
We have abstracted out $R$ in Eq. 2 for devising a simpler explanation. As
described earlier, the enhancer $E=E^{2}\circ E^{1}$ is mainly used to
improvise the output of the generator $G$ by making it consistent across
different resolutions of the target images.

Figure 2: Multi-level encoder-decoder network architecture. An example
construction of a three-level $med$ network. The generator $G$ is an $ed$
network and enhancer $E=E^{1}\circ E^{2}$ is the composition of two $ed$
networks. There are skip connections within each $ed$ subnetwork. The layers
are shown using colors as follows:  Convolutional layer with stride =1,
Convolutional layer with stride=2,  Batch Normalization, and  Upsampling. The
subnetwork $G$ is a depth-5 $ed$ network, $E^{1}$ is a depth-4 $ed$ network,
and $E^{2}$ is a depth-3 $ed$ network.

###  3.3 Network Classification

We have provided an example construction of a multi-level encoder-decoder
network in Fig. 2. Similarly, there are various other network architectures we
can get by modifying the network components. We give a classification of $med$
networks useful for _our methods_ to analyze these network architectures. The
$med$ network is classified based on skip links and cascading of the network
input, as shown in Table 1. The network MED has no skip connections and MEDS
has Intra-skip connections (the character “S” in MEDS denotes the presence of
skip connections). The network MEDSF has Full-skip connections. Similarly,
MEDC has cascading of network input without skip connections (the character
“C” in MEDC denotes the cascading of network input). MEDSFC has cascading of
network input with Full-skip connections. We will use the networks given in
Table 1 for our experiments. For example, to see the effects of the decreasing
skip links, one could perform image restoration with MEDSF, MEDS, and MED
networks.

|  No skip |  Intra-skip |  Full-Skip
---|---|---|---

Cascade |  MEDC |  MEDSC |  MEDSFC
---|---|---|---
No Cascade |  MED |  MEDS |  MEDSF

Table 1: Classification of $med$ networks. The classification is based on the
following network components: skip-links and cascading of network input at
intermediate layers. The graphical representations of the above network
components are shown in Fig. 1.

##  4 Applications

In this section, we show the performance on the following image restoration
tasks: super-resolution, denoising, inpainting, and flash no-flash. We provide
the technical details of the experiments in the supplementary material.

The aim of image restoration is to reconstruct the image features given a
corrupted image $\hat{I}$. The image $\hat{I}$ is computed by adding noise or
blur or downsampling the target image $I$. Ulyanov et al. formulated the image
restoration problem to the setting of DL based learning-free framework [29].
The image restoration framework is as follows.

$\displaystyle\theta^{*}=$
$\displaystyle{\underset{\theta}{\operatorname{arg\,\text{min}}}}\;\mathcal{L}\big{(}\mathcal{F}_{\theta}(\hat{z}),\hat{I}\big{)};$
(3)

Here, $\mathcal{L}$ is the loss function and $\mathcal{F}$ is a network with
parameters denoted by $\theta$ and the network input $z$ is prepared from the
corrupted image $\hat{I}$. The loss function in the Eq. 3 is a general
definition. We now discuss how to perform various image restoration tasks.

Denoising. Denoising aims to reduce noise and recover the clean image where
the learning process is assisted only by the corrupted image. Consider a noisy
image $\hat{I}$. Let $d_{1}=\mathcal{D}(\frac{1}{2},\hat{I})$ and
$d_{2}=\mathcal{D}(\frac{1}{4},\hat{I})$ be the down-sampled versions of the
image $\hat{I}$. Our approaches are based on the following property of a
natural image: patch recurrence within and across multiple scales. Using this
property, one could say that the down-sampled corrupted image contains some of
the image features. To make the best use of the property above, our multi-
scale loss $\mathcal{L}\big{(}\mathcal{F}_{\theta}(z),\hat{I}\big{)}$ (Eq. 3)
for denoising is defined in Eq. 4.

$\begin{split}\theta^{*}=&\;{\underset{\theta}{\operatorname{arg\,\text{min}}}}\;\lambda_{1}\|G_{\theta}(z)-\hat{I}\|\\\
&+\lambda_{2}\|E^{1}_{\theta}(z)-d_{1}\|+\lambda_{3}\|E^{2}_{\theta}(z)-d_{2}\|\end{split}$
(4)

Here, $\mathcal{F}=E^{2}\circ E^{1}\circ G$ (Eq. 2). Loss function in Eq. 4
enforces the output of the generator to be consistent across the multiple
resolutions of the target image. Stated differently, the network performs
image restoration at multiple resolutions. Intuitively, achieving restoration
at multiple scales is more challenging than at a single scale. Therefore, we
expect that solving a harder problem could help in learning a better image
prior [12]. The image prior is implicitly captured by the network which is
required to restore the image features [29].

Denoising using our MEDSF is shown in Fig. 4. Our MEDSF achieves SSIM=0.72
whereas the baseline DIP [29] outputs a SSIM of 0.71 for a noise strength of
$\sigma=100$. The PSNR values for our MEDSF is 20.95 and DIP outputs a PSNR of
21.36. In Fig. 6, we can observe that a higher PSNR value do not imply higher
perceptual quality. We emphasize that the learning-free methods are sensitive
to hyper-parameters444The learning-free methods are sensitive to hyper-
parameters shown in Fig. 4 of the supplementary material and DIP [29]..
Therefore, the performance of DIP and our MEDSF could probably be further
maximized by changing the hyper-parameters.

In Fig. 5, we can observe the effects of network depth on denoising. The
network initially learns the global features from the corrupted image by
minimizing the loss function defined in Eq. 4. Later, the network starts
learning fine feature details which includes noise. Therefore, due to over
learning, it produces noisy spots similar to the ones contained in the
corrupted image. For example, MEDSF intermediate output at around 1000
iterations is the desired noise free image because it achieves the maximum
PSNR.

(a) Original,
image

(b) Noisy,
image

(c) DIP,
(0.479, 18.65)

(d) Ours,
(0.496, 18.39)

Figure 4: Denoising. A comparison between DIP [29] and our MEDSF for denoising
with noise strength of $\sigma=100$ using the performance metric (SSIM, PSNR).

Figure 5: Network depth effects on denoising. EDS5 is a depth-5 $ed$ network
with skip connections (similarly for EDS6 and EDS7). The highest-depth network
MEDSF converges faster. EDS5 network (lower depth) achieves the highest PSNR
value but converges the slowest. This shows that a higher model capacity does
not necessarily lead to improved performance.

Super-resolution. Given a low-resolution (LR) image
$\hat{I}\in\mathbb{R}^{m\times n\times 3}$, and a scaling factor $t$, super-
resolution aims to enhance the image quality and generate a high-resolution
(HR) image $I^{H}\in\mathbb{R}^{mt\times nt\times 3}$. We feed network input
$z$ into $med$ network $\mathcal{F}=E^{2}\circ E^{1}\circ G$ and solve the
following minimization problem given in Eq. 5.

$\begin{split}\theta^{*}=\;&{\underset{\theta}{\operatorname{arg\,\text{min}}}}\;\lambda_{1}\|G_{\theta}(z)-u_{0}\|&\\\
&+\lambda_{2}\|E^{1}_{\theta}(z)-u_{1}\|+\lambda_{3}\|E^{2}_{\theta}(z)-\hat{I}\|&\end{split}$
(5)

Here, $u_{0}=\mathcal{U}(\hat{I},4)$ and $u_{1}=\mathcal{U}(\hat{I},2)$ are
the up-sampled versions of the corrupted LR image $\hat{I}$. Eq. 5 determines
the network parameter $\theta^{*}$ which minimizes the loss
$\mathcal{L}\big{(}\mathcal{F}_{\theta}(z),\hat{I}\big{)}$.

Super-resolution achieved by Ulyanov et al. in Deep Image Prior (DIP) is the
state-of-the-art in DL-based learning-free methods to the best of our
knowledge [29]. DIP does not use training samples to learn the image prior in
contrast to the learning-based methods which benefit from the training data
and adversarial loss or perceptual loss [23, 14]. Thus, it lacks local level
features in the output image. However, it is shown to output better images
than various learning-free methods such as bicubic upsampling [29].

We achieved an average SSIM of 0.80, whereas DIP [29] achieved an average SSIM
of 0.81 for $4\times$super-resolution. We obtained 24.48 as the average PSNR.
Whereas DIP achieved an average PSNR of 25.14555RGB images in Set14 dataset
had three channels and our $med$ network also outputs RGB images having three
channels. However, super-resolution output of DIP [29] have images with four
channels (including the alpha channel). Therefore, to get a _fair comparison_
, we reproduced the DIP output before drawing the comparison.. The perceptual
quality of the generated images by the proposed approach is observed to be
comparable to that of DIP (Fig. 6).

(a) High resolution image.

(b) Low resolution image.

(c) DIP [29], (0.88, 28.2).

(d) MEDSF, (0.88, 25.43).

Figure 6: $4\times$ Image super-resolution. A qualitative comparision using
performance metric (SSIM, PSNR). We can observe that a higher PSNR value does
not imply a higher perceptual quality.

Image inpainting. It involves computing missing pixel values in the corrupted
image $\hat{I}$ using the corresponding binary mask $m\in\\{0,1\\}^{k\times
l}$. Inpainting has various applications such as removing undesirable objects
and text in an image, restoring damaged paintings, and computing missing
pixels lost during transmission.

Suppose $\mathcal{I}$ is the target image and the corrupted image $\hat{I}$ is
obtained using the mask $m$ as follows $\hat{I}=\mathcal{I}\odot m$, where
$\odot$ is the Hadamard product. Let $d_{1}=\mathcal{D}(\frac{1}{2},\hat{I})$
and $d_{2}=\mathcal{D}(\frac{1}{4},\hat{I})$ be the down-sampled versions of
the corrupted image $\hat{I}$, and $m_{1}=\mathcal{D}(\frac{1}{2},m)$ and
$m_{2}=\mathcal{D}(\frac{1}{4},m)$ be the down-sampled versions of the mask
$m$. We solve the following minimization problem given in Eq. 6.

$\displaystyle\theta^{*}=$
$\displaystyle{\underset{\theta}{\operatorname{arg\,\text{min}}}}\;\lambda_{1}\|(G_{\theta}(z)-\hat{I})\odot
m\|$ (6) $\displaystyle+\lambda_{2}\|(E^{1}_{\theta}(z)-d_{1})\odot
m_{2}\|+\lambda_{3}\|(E^{2}_{\theta}(z)-d_{2})\odot m_{3}\|$

We show the following three inpainting tasks. (1) _restoring missing pixels_
lost by masking the target image with a randomly generated binary mask (Fig.
7), (2) _region-inpainting_ which includes painting a large region (Fig. 8 and
Fig. 11), and (3) _removing text_ superimposed on an image (Fig. 9).

Inpainting requires understanding the global context and the local structure
of the target image [34]. We believe that region-inpainting is the most
challenging task because the information from the nearby pixels might not
always be sufficient to complete the scene.

We obtained 24.62 as the average PSNR and 0.86 as the average SSIM for
inpainting with $90\%$ missing pixels. Whereas DIP [29] achieved an average
PSNR of 25.05 and an average SSIM of 0.86. The perceptual quality of the
generated images by the proposed approach is observed to be comparable to that
of the other methods (Fig. 7).

(a) Original image

(b) Corrupted image

(c) DIP [29], (0.85, 25.48)

(d) Ours, (0.85, 26.18)

Figure 7: Inpainting. A comparision for restoration of $90\%$ missing pixels
using performance metric (SSIM, PSNR).

(a) Original
image

(b) MEDSFC,
24.87

(c) MEDSF,
23.20

(d) MEDSC,
24.28

(e) MEDS,
24.04

(f) MEDC,
26.80

(g) MED,
25.83

Figure 8: Cascading of network input. Effects of cascading of the network
input in the intermediate layers of the network on removing an object from an
image given in (a). The vase present in (a) is removed using a white mask, and
then inpainting is performed. Networks in (b) and (c) have the same set of
skip links. Similarly, (d) and (e) have the same collection of skip links, and
(f) and (g) do not have skip links. Cascading of network input is performed in
(b), (d) and (f). (b) and (c) shows that cascading of network input into the
intermediate layers of the network improves the performance. Similarly, we can
observe that the cascading of network input performed better for other
networks: (d) and (e), and (f) and (g).

(a) Masked
image

(b) MEDSF*,
30.92

(c) MEDSF,
30.40

(d) MEDS*,
31.05

(e) MEDS,
30.35

(f) MED*,
30.18

(g) MED,
29.35

Figure 9: Composition of networks. Effects of the composition of the $ed$
networks. MEDSF∗, MEDS∗, and MED∗ are two levels $ed$ networks. MED, MEDS, and
MEDSF are three level $ed$ networks (enhancer is a composition of two $ed$
networks). (b) and (c) shows that the two-level full-skip network performed
better than three levels of the full-skip network. Similarly, we can observe
that the two-level $med$ network performed better for other networks: (d) and
(e), and (f) and (g).

(a) Flash image

(b) No-flash

(c) DIP, 17.03

(d) Ours, 18.54

Figure 10: Flash-no flash reconstruction. (a) Flash image. (b) No flash image.
(c) DIP. (d) Ours MEDS.

Flash No-flash. Given a pair of flash and no-flash images, the objective is to
get a single high-quality image which incorporates details of the scene from
the flash image and ambient illumination from the no-flash image [21, 9]. The
combined image helps to achieve denoising, white balancing, red-eye correction
[21], foreground extraction [27], and saliency detection [11].

Consider a pair $(I^{F},I^{NF})$, where $I^{F}$ is a flash image and $I^{NF}$
is a no-flash image. The network input $z$ is prepared by concatenating
$I^{F}$ and $I^{NF}$. Let $f_{1}=\mathcal{D}(\frac{1}{2},I^{NF})$ and
$f_{2}=\mathcal{D}(\frac{1}{4},I^{NF})$ be the down-sampled versions of
$I^{NF}$. We solve the optimization problem given in Eq. 7.

$\displaystyle\theta^{*}=$
$\displaystyle{\underset{\theta}{\operatorname{arg\,\text{min}}}}\;\lambda_{1}\Big{(}\|G_{\theta}(z)-I^{NF}\|+\|E^{1}_{\theta}(z)-f_{1}\|$
(7)
$\displaystyle+\|E^{2}_{\theta}(z)-f_{2}\|\Big{)}+\lambda_{2}\|G_{\theta}(z)-I^{F}\|$

Here, $\lambda_{1}$ and $\lambda_{2}$ are the coefficients to control the
image features from $I^{NF}$ and $I^{F}$. The flash no-flash output is shown
in Fig. 10. It is worth noting that our implementation of flash no-flash is
more flexible in providing features from both flash and no-flash images using
coefficients $\lambda_{1}$ and $\lambda_{2}$, unlike [29] (Fig. 12 of the
supplementary material).

##  5 Network Structures Effects on Restoration

(a) Masked image.

(b) MEDSF.

(c) MED.

Figure 11: Skip connections (I). The network with skip links (MEDSF) does not
perform well for region inpainting compared to the network without skip
connections.

(a) LR image.

(b) MED, 0.63.

(c) MEDSF, 0.70.

Figure 12: Skip connections (II). Skip links (MEDSF) improves 4$\times$ super-
resolution as shown by SSIM.

Here, we discuss the various aspects of the relation between the network
construction and image restoration using the $med$ framework. Our choice of
the multi-level architecture (a high capacity network) is motivated to
illustrate the behavior of various network components (Sec. 3). We emphasize
that the image restoration quality from untrained networks is sensitive to
hyper-parameters search [29]. We now discuss the results of the ablation
studies that we have conducted.

Effects of Skip links. Skip connections have shown adverse effects on
inpainting, see Fig. 8, Fig. 9, and Fig. 11 (the number of skip connections in
the above figures decreases from left to right). Our interpretation of the
adverse effects is as follows. The layers of encoder have under-developed
regions and their pixel values are close to that of the mask (either zero or
one). The skip connections pass such intermediate representation to the
decoder, which leads to reconstruction bias. Therefore, output images have
pixel values that are close to the mask.

Effects of Depth.  In Fig. 5, we observe a higher the depth network converges
faster because it has a large number of parameters. However, a lower depth
network EDS5 could achieve better restoration than the higher depth network
MEDSF. There could be two major factors for the above result. First, higher
depth network suffer from the vanishing gradient problem which negatively
influences the performance [18, 26]. Second, the increase in the number skip
connections due to higher depth, influence the performance positively [18]. We
believe that the decrease in the PSNR value indicates that the negative
influence of the network depth could have more impact compared to the
performance enhancement we get from skip connections.

Effects of Cascading of network input (cascade). The cascade and the skip
connection looks similar because they both provide image features to the
intermediate layers of the network. However, they provide a different type of
image features. Cascade provides image features from the corrupted image.
Whereas, skip connections pass the image features from the intermediate layers
of the network. Object removal (inpainting) is better achieved using cascade
(Fig. 8). Whereas providing image features using skip connections have shown
adverse effects for inpainting (Fig. 11). This could be because of the image
features captured at the intermediate layers of the network are less
interpretable than the features present in the corrupted image.

|  DIP [29] |  $med$ (Ours)
---|---|---

Depth |  ✔ |  ✔
---|---|---
Skip-links |  ✔ |  ✔
Composition of $ed$ |  ✖ |  ✔
Cascading of input |  ✖ |  ✔

Table 2: Network components to investigate the influence of the network
structures for image restoration tasks.

|  PSNR |  SSIM
---|---|---

|  DIP [29] |  Ours |  DIP [29] |  Ours
---|---|---|---|---
Denoising |  21.36 |  20.95 |  0.71 |  0.72
Inpainting |  25.05 |  24.62 |  0.86 |  0.86
SISR |  25.14 |  24.48 |  0.81 |  0.80

Table 3: A quantitative comparison for denoising, inpainting, and single image
super-resolution (SISR) using average PSNR and SSIM. We provide the visual
comparison of generated images in the supplementary material. The perceptual
quality of the generated images is comparable to DIP [29] despite the $med$
network has a higher capacity to accommodate various network components (Table
2).

Effects of Composition of $ed$ networks. The two-level $med$ network performed
better than a three level $med$ network for text-removal from an image (Fig.
9). However, the performance difference is not very significant (less than one
PSNR). A network composition increases the network depth and the number of
skip connections. Therefore, a three-level $med$ could have more influence on
restoration from skip connections compared to a two-level $med$ network.
Similarly, a three-level $med$ could also increase the effects of vanishing
gradients due to the higher depth. The composition of the networks shows the
combined effects of increasing depth and skip connections.

##  6 Conclusion

We have shown interesting aspects of the relationship between image
restoration and network construction. Our methods are unsupervised and they
only use the corrupted image for restoration instead of using any training
data. Therefore, we believe that it does not produce a biased output unlike
learning-based methods, e.g., model collapse [24]. We feel it is a challenging
experimental setup compared to supervised learning setup because the network
is not learning image features by the pairs of low and high-quality images.
The challenge is the limited contextual understanding due to the lack of
feature learning from the training data.

Our $med$ framework is a generalization of DIP [29]. This generalization is
novel because it incorporates various network components and an enhancer
network. The $med$ framework is more expressive in terms of casting different
network structures to perform the ablation studies for various aspects of the
network (Table 2). We also discuss image restoration task specific network
structures that perform comparably to the state-of-the-art methods (Table 3).

The major components of the restoration framework are the network and the loss
function (Eq. 3). We have shown analysis using various network structures and
MSE loss666In Fig. 13 of the supplementary material, we show that MSE
performed better than contextual loss [20].. The study of MSE loss is useful
as it is used in other image restoration methods. For example, MSE with
adversarial loss in [13, 25] and MSE with contextual loss in [19].

We observed that some network components do not enhance the restoration
quality. For example, a network with skip links does not perform well for
inpainting. Therefore, the experiments on a network with skip connections for
inpainting will not be efficient. Wang et al. have used skip connections for
video inpainting [32]. However, their approach is in the supervised learning
setup, unlike our unsupervised setup. We believe that there are similarities
in both setups. For example, if a network component is negatively influencing
the image prior learning from the corrupted image (unsupervised setup), then
it should also negatively influence the learning from the multiple images of
training data (supervised setup). We propose as future work to study our
restoration framework in the supervised learning setup.

## References

  * [1] D. Bau, B. Zhou, A. Khosla, A. Oliva, and A. Torralba.  Network dissection: Quantifying interpretability of deep visual representations.  CVPR, 2017.
  * [2] A. Bhowmik, S. Shit, and C. S. Seelamantula.  Training-free, single-image super-resolution using a dynamic convolutional network.  IEEE Signal Processing Letters, 25(1):85–89, 2018.
  * [3] S. A. Bigdeli and M. Zwicker.  Image restoration using autoencoding priors.  arXiv preprint arXiv:1703.09964, 2017.
  * [4] A. Buades, B. Coll, and J.-M. Morel.  A non-local algorithm for image denoising.  In Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, volume 2, pages 60–65. IEEE, 2005.
  * [5] H. C. Burger, C. J. Schuler, and S. Harmeling.  Image denoising: Can plain neural networks compete with bm3d?  In Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on, pages 2392–2399. IEEE, 2012.
  * [6] S. G. Chang, B. Yu, and M. Vetterli.  Adaptive wavelet thresholding for image denoising and compression.  IEEE transactions on image processing, 9(9):1532–1546, 2000.
  * [7] Q. Chen and V. Koltun.  Photographic image synthesis with cascaded refinement networks.  In The IEEE International Conference on Computer Vision (ICCV), volume 1, 2017.
  * [8] C. E. Duchon.  Lanczos filtering in one and two dimensions.  Journal of applied meteorology, 18(8):1016–1022, 1979.
  * [9] E. Eisemann and F. Durand.  Flash photography enhancement via intrinsic relighting.  In ACM transactions on graphics (TOG), volume 23, pages 673–678. ACM, 2004.
  * [10] Y. Gandelsman, A. Shocher, and M. Irani.  Double-dip: Unsupervised image decomposition via coupled deep-image-priors.  arXiv preprint arXiv:1812.00467, 2018.
  * [11] S. He and R. W. Lau.  Saliency detection with flash and no-flash image pairs.  In European Conference on Computer Vision, pages 110–124. Springer, 2014.
  * [12] D. Kim, D. Cho, D. Yoo, and I. S. Kweon.  Learning image representations by completing damaged jigsaw puzzles.  WACV, 2018.
  * [13] C. Ledig, L. Theis, F. Huszár, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang, et al.  Photo-realistic single image super-resolution using a generative adversarial network.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 4681–4690, 2017.
  * [14] C. Ledig, L. Theis, F. Huszar, J. Caballero, A. Cunningham, A. Acosta, A. Aitken, A. Tejani, J. Totz, Z. Wang, and W. Shi.  Photo-realistic single image super-resolution using a generative adversarial network.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), July 2017.
  * [15] S. Lefkimmiatis.  Non-local color image denoising with convolutional neural networks.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), July 2017.
  * [16] Z. Li, Z. Murez, D. Kriegman, R. Ramamoorthi, and M. Chandraker.  Learning to see through turbulent water.  In 2018 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 512–520. IEEE, 2018.
  * [17] J. Mairal, F. Bach, J. Ponce, and G. Sapiro.  Online learning for matrix factorization and sparse coding.  Journal of Machine Learning Research, 11(Jan):19–60, 2010.
  * [18] X. Mao, C. Shen, and Y.-B. Yang.  Image restoration using very deep convolutional encoder-decoder networks with symmetric skip connections.  In Advances in neural information processing systems, pages 2802–2810, 2016.
  * [19] R. Mechrez, I. Talmi, F. Shama, and L. Zelnik-Manor.  Learning to maintain natural image statistics.  arXiv preprint arXiv:1803.04626, 2018.
  * [20] R. Mechrez, I. Talmi, and L. Zelnik-Manor.  The contextual loss for image transformation with non-aligned data.  European Conference on Computer Vision (ECCV), 2018.
  * [21] G. Petschnigg, R. Szeliski, M. Agrawala, M. Cohen, H. Hoppe, and K. Toyama.  Digital photography with flash and no-flash image pairs.  In ACM transactions on graphics (TOG), volume 23, pages 664–672. ACM, 2004.
  * [22] J. H. Rick Chang, C.-L. Li, B. Poczos, B. V. K. Vijaya Kumar, and A. C. Sankaranarayanan.  One network to solve them all – solving linear inverse problems using deep projection models.  In The IEEE International Conference on Computer Vision (ICCV), Oct 2017.
  * [23] M. S. M. Sajjadi, B. Scholkopf, and M. Hirsch.  Enhancenet: Single image super-resolution through automated texture synthesis.  In The IEEE International Conference on Computer Vision (ICCV), Oct 2017.
  * [24] T. Salimans, I. Goodfellow, W. Zaremba, V. Cheung, A. Radford, and X. Chen.  Improved techniques for training gans.  In Advances in neural information processing systems, pages 2234–2242, 2016.
  * [25] A. Shocher, S. Bagon, P. Isola, and M. Irani.  Internal distribution matching for natural image retargeting.  arXiv preprint arXiv:1812.00231, 2018.
  * [26] R. K. Srivastava, K. Greff, and J. Schmidhuber.  Training very deep networks.  In Advances in neural information processing systems, pages 2377–2385, 2015.
  * [27] J. Sun, S. B. Kang, Z.-B. Xu, X. Tang, and H.-Y. Shum.  Flash cut: Foreground extraction with flash and no-flash image pairs.  In Computer Vision and Pattern Recognition, 2007. CVPR’07. IEEE Conference on, pages 1–8. IEEE, 2007.
  * [28] C. Tomasi and R. Manduchi.  Bilateral filtering for gray and color images.  In Computer Vision, 1998. Sixth International Conference on, pages 839–846. IEEE, 1998.
  * [29] D. Ulyanov, A. Vedaldi, and V. Lempitsky.  Deep image prior.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), June 2018.
  * [30] R. G. Vidal, S. Banerjee, K. Grm, V. Struc, and W. J. Scheirer.  Ug^ 2: A video benchmark for assessing the impact of image restoration and enhancement on automatic visual recognition.  In 2018 IEEE Winter Conference on Applications of Computer Vision (WACV), pages 1597–1606. IEEE, 2018.
  * [31] P. Vincent, H. Larochelle, Y. Bengio, and P.-A. Manzagol.  Extracting and composing robust features with denoising autoencoders.  In Proceedings of the 25th international conference on Machine learning, pages 1096–1103. ACM, 2008.
  * [32] C. Wang, H. Huang, X. Han, and J. Wang.  Video inpainting by jointly learning temporal structure and spatial details.  arXiv preprint arXiv:1806.08482, 2018.
  * [33] J. Xie, L. Xu, and E. Chen.  Image denoising and inpainting with deep neural networks.  In Advances in neural information processing systems, pages 341–349, 2012.
  * [34] C. Yang, X. Lu, Z. Lin, E. Shechtman, O. Wang, and H. Li.  High-resolution image inpainting using multi-scale neural patch synthesis.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), volume 1, page 3, 2017.
  * [35] D. Yang and J. Sun.  Bm3d-net: A convolutional neural network for transform-domain collaborative filtering.  IEEE Signal Processing Letters, 25(1):55–59, 2018.
  * [36] R. A. Yeh, C. Chen, T. Y. Lim, A. G. Schwing, M. Hasegawa-Johnson, and M. N. Do.  Semantic image inpainting with deep generative models.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5485–5493, 2017.
  * [37] R. Zeyde, M. Elad, and M. Protter.  On single image scale-up using sparse-representations.  In International conference on curves and surfaces, pages 711–730. Springer, 2010.
  * [38] K. Zhang, W. Zuo, S. Gu, and L. Zhang.  Learning deep cnn denoiser prior for image restoration.  In The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), July 2017.
  * [39] F. Zhu, G. Chen, and P.-A. Heng.  From noise modeling to blind image denoising.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 420–429, 2016.


A_{1}}{\delta{\tilde{\Upsilon}}_{ab}}=-\frac{\hbar^{2}}{8}\sum_{hk}\tilde{\Upsilon}_{ch}\tilde{\Upsilon}_{ck}\frac{\partial{\tilde{\Upsilon}}^{-1}_{hk}}{\partial\tilde{\Upsilon}_{ab}}+2\hbar^{2}\sum_{cde}\tilde{C}_{cd}\frac{\partial{\tilde{\Upsilon}}^{-1}_{de}}{\partial\tilde{\Upsilon}_{ab}}\tilde{C}_{ec}$
(199)

In the same way $A_{2}$:

$\frac{\delta
A_{2}}{\delta\Upsilon_{ab}}=\frac{1}{2}\sum_{cd}\frac{\partial\Upsilon^{-1}_{cd}}{\partial\Upsilon_{ab}}\left\langle\frac{\partial^{2}{V^{(\text{tot})}}}{\partial
R_{c}\partial R_{d}}\right\rangle_{\rho}$ (200) $\frac{\delta
A_{2}}{\delta\mathcal{R}_{a}}=-\left\langle f_{a}^{(tot)}\right\rangle_{\rho}$
(201) $\frac{\delta A_{2}}{\delta Q_{a}}=\frac{\delta A_{2}}{\delta C_{ab}}=0$
(202)

And for $A_{3}$:

$\frac{\delta A_{3}}{\delta Q_{a}}=\hbar\frac{d\mathcal{R}_{a}}{dt}$ (203)
$\frac{\delta A_{3}}{\delta\mathcal{R}_{a}}=-\hbar\frac{dQ_{a}}{dt}$ (204)
$\frac{\delta A_{3}}{\delta
C_{ab}}=-\hbar\sum_{cd}\frac{\partial\Upsilon^{-1}_{ab}}{\partial\Upsilon_{cd}}\frac{d\Upsilon_{cd}}{dt}$
(205) $\frac{\delta
A_{3}}{\delta\Upsilon_{ab}}=\hbar\sum_{cd}\frac{\partial\Upsilon^{-1}_{cd}}{\partial\Upsilon_{ab}}\frac{dC_{cd}}{dt}$
(206)

By imposing that:

$\delta A=\delta A_{1}+\delta A_{2}+\delta A_{3}=0$

we get the equation of motion.

$\frac{\delta A}{\delta Q}=0$ $\frac{d\mathcal{R}_{a}}{dt}=\frac{\hbar
Q_{a}}{m_{a}}$ (207) $\frac{\delta A}{\delta\mathcal{R}}=0$
$\frac{dQ_{a}}{dt}=\frac{1}{\hbar}\left\langle
f_{a}^{(tot)}\right\rangle_{\rho(t)}$ (208)

Indeed, Eq. (208) is equal to Eq. (54a). We can proceed to get the equations
for $\bm{C}$ and $\bm{\Upsilon}$. To proceed, we need to remove the
$\partial\bm{\Upsilon}^{-1}/\partial\bm{\Upsilon}$ tensor:

$\frac{\partial\Upsilon^{-1}_{cd}}{\partial\Upsilon_{ab}}=-\frac{1}{2}\Upsilon_{ac}^{-1}\Upsilon_{bd}^{-1}-\frac{1}{2}\Upsilon_{cb}^{-1}\Upsilon_{ad}^{-1}$
(209) $\frac{\delta A}{\delta C_{ab}}=0$
$\frac{d\tilde{\bm{\Upsilon}}}{dt}=2\hbar\left({\tilde{\bm{\Upsilon}}}\tilde{\bm{C}}+\tilde{\bm{C}}{\tilde{\bm{\Upsilon}}}\right)$
(210)

Eq. (210) is equal to Eq. (54b) if we set $\bm{A}=0$.

Lastly,

$\frac{\delta A}{\delta\Upsilon_{ab}}=0$
$\frac{d\tilde{\bm{C}}}{dt}=\frac{1}{2\hbar}\left\langle\frac{\partial^{2}V}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}\right\rangle_{\rho}-\frac{\hbar}{8}\tilde{\bm{\Upsilon}}\tilde{\bm{\Upsilon}}+2\hbar\tilde{\bm{C}}\tilde{\bm{C}}$
(211)

Indeed, also in this case, Eq. (211) is equal to Eq. (54e) if we set
$\bm{A}=0$.

Therefore, we found that the dynamical equations for obtained with the least
action principle coincides with those obtained from the TD-SCHA.

##  Appendix D Equilibrium

Here we show that, among all possible stationary solutions, the one that
minimizes the functional:

$F=E-TS$

is the SCHA solution, i.e. the one where temperature is uniform on all the
degrees of freedom.

Since we are in a stationary solution, we can write everything in the basis of
$\bm{e_{\mu}}$ vectors that diagonalizes the $\bm{\Upsilon}$, $\bm{A}$ and
$\left\langle\partial^{2}V/\partial\tilde{\bm{R}}^{2}\right\rangle_{\rho}$.

The $F$ functional is:

$F[\hat{\rho}]=\left\langle\hat{K}+\hat{V}\right\rangle_{\hat{\rho}}-TS[\hat{\rho}]$
(212)

The average of the kinetic energy (Eq. 63) is:

$\left\langle\hat{K}\right\rangle_{\hat{\rho}}=\frac{\hbar^{2}}{2}\sum_{\mu}\left(\frac{\tilde{\Upsilon}_{\mu}}{4}+\tilde{A}_{\mu}\right)$
(213)

Substituting Eq. (77b), we get:

$\left\langle\hat{K}\right\rangle_{\hat{\rho}}=\sum_{\mu}\frac{\hbar\omega_{\mu}(2n_{\mu}+1)}{4}$
(214)

The entropy $S[\hat{\rho}]$ can be obtained as the sum of the entropies of the
Harmonic oscillators on the $\mu$ modes:

$S=\frac{k_{b}}{2}\sum_{\mu}\left[\frac{\beta_{\mu}\hbar\omega_{\mu}}{\tanh\frac{\beta_{\mu}\hbar\omega}{2}}-2\log\left(\sinh\frac{\beta_{\mu}\hbar\omega_{\mu}}{2}\right)\right]$
(215)

We must impose that

$\frac{\partial F}{\partial\beta_{\mu}}=0$ (216)
$\frac{\partial\braket{\hat{K}}}{\partial\beta_{\mu}}+\frac{\partial\braket{\hat{V}}}{\partial\beta_{\mu}}-T\frac{\partial
S}{\partial\beta_{\mu}}=0$ (217)
$\frac{\partial\braket{\hat{K}}}{\partial\beta_{\mu}}=-\frac{\hbar^{2}\omega_{\mu}^{2}}{8\sinh^{2}\left(\frac{\beta_{\mu}\hbar\omega_{\mu}}{2}\right)}$
(218)

The derivative of the average of the potential can be computed exploiting the
formalism introduced by Bianco et al.[33]. In particular, they showed that the
derivative of any observable with respect to a the quantum fluctuations is:

$\frac{\partial\left\langle
V(\bm{R})\right\rangle_{\rho}}{\partial\Upsilon_{\mu}}=\frac{1}{2}\sum_{ab}\frac{\partial{{\tilde{\Upsilon}}^{-1}}_{ab}}{\partial\Upsilon_{\mu}}\left\langle\frac{\partial^{2}V}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}}\right\rangle_{\rho}$
(219)

We can write it in the polarization basis (exploiting Eq. (77b):

$\frac{\partial\left\langle
V(\bm{R})\right\rangle_{\rho}}{\partial\Upsilon_{\mu}}=-\frac{1}{2}{\tilde{\Upsilon}}_{\mu}^{-2}\omega_{\mu}^{2}$
(220)

From which we get the average of the potential:

$\frac{\partial\left\langle
V(\bm{R})\right\rangle_{\rho}}{\partial\beta_{\mu}}=-\frac{\hbar^{2}\omega_{\mu}^{2}}{8\sinh^{2}\left(\frac{\beta_{\mu}\hbar\omega_{\mu}}{2}\right)}$
(221)

While the derivative of the entropy is:

$\frac{\partial
S}{\partial\beta_{\mu}}=-\frac{\beta_{\mu}\hbar^{2}\omega_{\mu}^{2}k_{b}}{4\sinh^{2}\left(\frac{\beta_{\mu}\hbar\omega_{\mu}}{2}\right)}$
(222)

Putting all together we get:

$\frac{2\hbar^{2}\omega_{\mu}^{2}-k_{b}T2\beta_{\mu}\hbar^{2}\omega_{\mu}^{2}}{8\sinh^{2}\left(\frac{\beta_{\mu}\hbar\omega_{\mu}}{2}\right)}=0$
(223)

And we get the condition that $\beta_{\mu}$ must satisfy to minimize the free
energy:

$\beta_{\mu}=\frac{1}{k_{b}T}$ (224)

##  Appendix E Derivation of the linear response system

Here, we apply perturbation theory on the TD-SCHA equations for small
perturbations. Since we are expanding around equilibrium solutions, we have:

$\bm{C}=0\qquad\bm{Q}=0\qquad\mathrm{Im}\bm{A}=0$ (225)

When $t\geq t_{0}$, we add an external time dependent perturbation to the non
interacting potential:

$V(\bm{R},t)=V^{(0)}(\bm{R})+V^{(1)}(\bm{R},t)$ (226)

And we want to study the dynamics of the system close to equilibrium. This
perturbation affects the dynamics (Eq. 54e) only in the two averages:

$\left\langle
f_{a}\right\rangle_{\rho}\qquad\left\langle\frac{d^{2}V}{dR_{a}dR_{b}}\right\rangle_{\rho}$
(227)

We indicate with (1) the perturbed quantities. If omitted, we refer to static
quantities.

Since the perturbation changes the parameters of the density matrix, they also
change the ionic probability distribution:

$\rho(\bm{R},t)=\rho(\bm{R})+\rho^{(1)}(\bm{R},t)$ (228)

This affect also the generic average over an observable:

$\left\langle O(\bm{R})\right\rangle_{\rho(t)}=\left\langle
O(\bm{R})\right\rangle_{{\rho}}+\left\langle
O^{(0)}(\bm{R})\right\rangle_{\rho^{(1)}}$ (229)

By expanding the probability distribution $\rho^{(1)}(\bm{R},t)$ as a function
of the perturbed centroid position $\bm{\mathcal{R}}^{(1)}$ and the perturbed
fluctuations $\bm{\Upsilon}^{(1)}$, we get:

$\displaystyle\left\langle
O\right\rangle_{\rho^{(1)}}=\frac{1}{2}\sum_{ab}\bigg{[}$
$\displaystyle\Upsilon_{ab}\left(\left\langle
Ou_{a}\right\rangle_{{\rho}}\mathcal{R}_{b}^{(1)}+\mathcal{R}_{a}^{(1)}\left\langle
Ou_{b}\right\rangle_{{\rho}}\right)+$
$\displaystyle-\Upsilon_{ab}^{(1)}\left\langle
u_{a}u_{b}O\right\rangle_{{\rho}}\bigg{]}$ (230)

where $u$ is the displacement with respect to the equilibrium centroid. In
particular, we are interested in the averages of forces and second derivative
of the BO potential:

$\displaystyle\left\langle\frac{dV}{dR_{a}}\right\rangle_{\rho^{(1)}}$
$\displaystyle=\frac{1}{2}\sum_{hk}\bigg{[}\Upsilon_{hk}\left\langle\frac{dV}{dR_{a}}u_{h}\right\rangle_{{\rho}}\mathcal{R}^{(1)}_{k}+$
$\displaystyle+\Upsilon_{hk}\mathcal{R}_{h}^{(1)}\left\langle\frac{dV}{dR_{a}}u_{k}\right\rangle_{{\rho}}-\Upsilon_{hk}^{(1)}\left\langle
u_{h}u_{k}\frac{dV}{dR_{a}}\right\rangle_{{\rho}}\bigg{]}$ (231)

Using the definition of the 3 and 4 phonon scattering tensor
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{\Phi}}$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{\Phi}}$ (ref.[33]), we have:

$\Phi_{ab}=\sum_{p}\Upsilon_{ap}\left\langle
u_{p}\frac{dV}{dR_{b}}\right\rangle_{{\rho}}$ (232)
$\overset{\scriptscriptstyle{\text{(3)}}}{\Phi}_{abc}=\sum_{pq}\Upsilon_{ap}\Upsilon_{bq}\left\langle
u_{p}u_{q}\frac{dV}{dR_{c}}\right\rangle_{{\rho}}$ (233)
$\overset{\scriptscriptstyle{\text{(4)}}}{\Phi}_{abcd}=\sum_{pqr}\Upsilon_{ap}\Upsilon_{bq}\Upsilon_{cr}\left\langle
u_{p}u_{q}u_{r}\frac{dV}{dR_{d}}\right\rangle_{{\rho}}$ (234)

Thus we have:

$\displaystyle\left\langle\frac{dV}{dR_{a}}\right\rangle_{\rho^{(1)}}=-\frac{1}{2}$
$\displaystyle\sum_{hkpqrs}\Upsilon^{(1)}_{hk}\Upsilon^{-1}_{kp}\Upsilon_{ps}\Upsilon^{-1}_{hq}\Upsilon_{qr}\left\langle
u_{r}u_{s}\frac{dV}{dR_{a}}\right\rangle_{{\rho}}+$
$\displaystyle+\sum_{h}\Phi_{ah}\mathcal{R}_{h}^{(1)}$ (235)

$\left\langle\frac{dV}{dR_{a}}\right\rangle_{\rho^{(1)}}=\sum_{h}\Phi_{ah}\mathcal{R}_{h}^{(1)}-\frac{1}{2}\sum_{hkpq}\Upsilon^{(1)}_{hk}\Upsilon^{-1}_{kp}\Upsilon^{-1}_{hq}\overset{\scriptscriptstyle{\text{(3)}}}{\Phi}_{qpa}$
(236)

To simplify notation, in the rest of this section we use the convention:

$\braket{\cdot}=\left\langle\cdot\right\rangle_{{\rho}}$

in an analogous way we can get the other term of the perturbation:

$\displaystyle\left\langle\frac{d^{2}V}{dR_{a}dR_{b}}\right\rangle_{\rho^{(1)}}=\frac{1}{2}\sum_{hk}$
$\displaystyle\bigg{[}\Upsilon_{hk}\braket{\frac{d^{2}V}{dR_{a}dR_{b}}u_{h}}\mathcal{R}^{(1)}_{k}+$
$\displaystyle+\Upsilon_{hk}\mathcal{R}_{h}^{(1)}\braket{\frac{d^{2}V}{dR_{a}dR_{b}}u_{k}}+$
$\displaystyle-\Upsilon_{hk}^{(1)}\braket{u_{h}u_{k}\frac{d^{2}V}{dR_{a}dR_{b}}}\bigg{]}$
(237)

Using the Bianco formalism introduced in Eq. (184), we get:

$\displaystyle\left\langle\frac{d^{2}V}{dR_{a}dR_{b}}\right\rangle_{\rho^{(1)}}=-\frac{1}{2}$
$\displaystyle\sum_{hkpq}\Upsilon^{-1}_{kp}\Upsilon^{-1}_{hq}\stackrel{{\scriptstyle(4)}}{{\Phi}}_{abpq}\Upsilon^{(1)}_{hk}+$
$\displaystyle+\sum_{h}\stackrel{{\scriptstyle(3)}}{{\Phi}}_{abh}\mathcal{R}_{h}^{(1)}$
(238)

We can derive now the explicit expression of all the perturbed equation of
motion:

$\displaystyle m_{a}\ddot{\mathcal{R}_{a}^{(1)}}$
$\displaystyle=\braket{f_{a}^{(1)}}-\sum_{h}\Phi_{ah}\mathcal{R}_{h}^{(1)}+$
$\displaystyle+\frac{1}{2}\sum_{hkpq}\Upsilon^{(1)}_{hk}\Upsilon^{-1}_{kp}\Upsilon^{-1}_{hq}\stackrel{{\scriptstyle(3)}}{{\Phi_{qpa}}}$
(239)

$\displaystyle\dot{\Upsilon^{(1)}}_{ab}$
$\displaystyle=\sum_{c}\frac{\hbar}{m_{c}}\bigg{[}\Upsilon_{cb}(2C_{ac}^{(1)}-\mathrm{Im}A_{ac}^{(1)})+\Upsilon_{ac}(2C_{cb}^{(1)}+$
$\displaystyle+\mathrm{Im}A_{cb}^{(1)})+\mathrm{Im}A_{cb}^{(1)}\mathrm{Re}A_{ac}-\mathrm{Im}A_{ac}^{(1)}\mathrm{Re}A_{cb}\bigg{]}$
(240)

$\displaystyle\mathrm{Im}\dot{A_{ab}^{(1)}}$
$\displaystyle=\sum_{c}\frac{\hbar}{2m_{c}}\Big{(}\Theta_{ac}^{(1)}\mathrm{Re}A_{cb}+$
$\displaystyle-\Theta_{bc}^{(1)}\mathrm{Re}A_{ac}+\Theta_{ac}\mathrm{Re}A_{cb}^{(1)}-\Theta_{bc}\mathrm{Re}A_{ac}^{(1)}\Big{)}$
(241)

$\displaystyle\mathrm{Re}\dot{A_{ab}^{(1)}}$
$\displaystyle=\sum_{c}\frac{\hbar}{2m_{c}}\Big{(}4C_{ac}^{(1)}\mathrm{Re}A_{cb}+4C_{bc}^{(1)}\mathrm{Re}A_{ca}+$
$\displaystyle-\Theta_{ac}\mathrm{Im}A_{cb}^{(1)}+\Theta_{bc}\mathrm{Im}A_{ac}^{(1)}\Big{)}$
(242)

$\displaystyle\dot{C_{ab}^{(1)}}$
$\displaystyle=\frac{1}{2\hbar}\braket{\frac{d^{2}V^{(1)}}{dR_{a}dR_{b}}}+\frac{1}{2\hbar}\sum_{h}\stackrel{{\scriptstyle(3)}}{{\Phi}}_{abh}\mathcal{R}_{h}^{(1)}+$
$\displaystyle-\frac{1}{4\hbar}\sum_{hkpq}\Upsilon^{-1}_{kp}\Upsilon^{-1}_{hq}\stackrel{{\scriptstyle(4)}}{{\Phi}}_{abpq}\Upsilon^{(1)}_{hk}+$
$\displaystyle+\sum_{c}\frac{\hbar}{2m_{c}}\bigg{(}-\frac{1}{4}\Theta_{ac}^{(1)}\Theta_{cb}-\frac{1}{4}\Theta_{ac}\Theta_{cb}^{(1)}+$
$\displaystyle+\mathrm{Re}A_{ac}\mathrm{Re}A_{bc}^{(1)}+\mathrm{Re}A_{ac}^{(1)}\mathrm{Re}A_{bc}\bigg{)}$
(243)

$\Theta^{(1)}_{ab}=\Upsilon_{ab}^{(1)}+2\mathrm{Re}A_{ab}^{(1)}$ (244)

We can remove the theta dependence on ${\mathrm{Im}\dot{\bm{A}}}$ and
$\dot{\bm{C}}$:

$\displaystyle\mathrm{Im}\dot{A_{ab}^{(1)}}$
$\displaystyle=\sum_{c}\frac{\hbar}{2m_{c}}\Big{(}\mathrm{Re}A_{cb}\Upsilon^{(1)}_{ac}-\mathrm{Re}A_{ac}\Upsilon^{(1)}_{bc}+$
$\displaystyle+\mathrm{Re}A_{cb}^{(1)}\Upsilon_{ac}-\mathrm{Re}A_{ac}^{(1)}\Upsilon_{bc}\Big{)}$
(245) $\displaystyle\dot{C}^{(1)}_{ab}$
$\displaystyle=\frac{1}{2\hbar}\braket{\frac{d^{2}V^{(1)}}{dR_{a}dR_{b}}}+\frac{1}{2\hbar}\sum_{h}\stackrel{{\scriptstyle(3)}}{{\Phi}}_{abh}\mathcal{R}_{h}^{(1)}+$
$\displaystyle-\frac{1}{4\hbar}\sum_{hkpq}\Upsilon^{-1}_{kp}\Upsilon^{-1}_{hq}\stackrel{{\scriptstyle(4)}}{{\Phi}}_{abpq}\Upsilon^{(1)}_{hk}+$
$\displaystyle+\sum_{c}\frac{\hbar}{8m_{c}}\Big{(}2\mathrm{Re}A_{ac}\Upsilon_{cb}^{(1)}+2\mathrm{Re}A_{cb}\Upsilon_{ac}^{(1)}+$
$\displaystyle+2\mathrm{Re}A_{ac}^{(1)}\Upsilon_{cb}+2\mathrm{Re}A_{cb}^{(1)}\Upsilon_{ac}-\Upsilon_{ac}\Upsilon^{(1)}_{cb}-\Upsilon_{cb}\Upsilon^{(1)}_{ac}\Big{)}$
(246)

By further deriving $\dot{\bm{\Upsilon}}^{(1)}$ and
$\mathrm{Re}\dot{\bm{A}^{(1)}}$ we delete two variables from the equations
($\bm{C}$ and $\mathrm{Im}\bm{A}$). And arrive to the final system of Eq.
(81c).

The $\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ terms are obtained from the
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{\Phi}}$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{\Phi}}$ when transforming in the
mass-rescaled variables.

$\overset{\scriptscriptstyle{\text{(3)}}}{D}_{abc}=\frac{\overset{\scriptscriptstyle{\text{(3)}}}{\Phi}_{abc}}{\sqrt{m_{a}m_{b}m_{c}}}$
(247)
$\overset{\scriptscriptstyle{\text{(4)}}}{D}_{abcd}=\frac{\overset{\scriptscriptstyle{\text{(4)}}}{\Phi}_{abcd}}{\sqrt{m_{a}m_{b}m_{c}m_{d}}}$
(248)

##  Appendix F Full expression of the Linear Perturbation system

Here we report the full expressions of the 2,3,4-rank tensors that define the
linear response system of the TD-SCHA. These are obtained writing the full
expression of the system as derived in appendix E. These are expressed in the
polarization basis, i.e. the basis that diagonalizes the SCHA dynamical matrix
at equilibrium:

$\displaystyle{\bar{X}}_{\mu\nu\eta\lambda}$
$\displaystyle=-\frac{\hbar(2n_{\eta}+1)(2n_{\lambda}+1)(2\omega_{\mu}n_{\nu}+2\omega_{\nu}n_{\mu}+\omega_{\mu}+\omega_{\nu})}{4(2n_{\mu}+1)(2n_{\nu}+1)\omega_{\eta}\omega_{\lambda}}{\overset{\scriptscriptstyle{\text{(4)}}}{D}}_{\mu\nu\eta\lambda}-\frac{\delta_{\mu\eta}\delta_{\nu\lambda}+\delta_{\mu\lambda}\delta_{\nu\eta}}{2}\left(\omega_{\mu}^{2}+\omega_{\nu}^{2}+\frac{2\omega_{\mu}\omega_{\nu}}{(2n_{\mu}+1)(2n_{\nu}+1)}\right)$
(249) $\displaystyle{\bar{X}^{\prime}}_{\mu\nu\eta\lambda}$
$\displaystyle=-\frac{\hbar\left[\omega_{\mu}n_{\mu}(n_{\mu}+1)(2n_{\nu}+1)+\omega_{\nu}n_{\nu}(n_{\nu}+1)(2n_{\mu}+1)\right](2n_{\eta}+1)(2n_{\lambda}+1)}{4(2n_{\mu}+1)(2n_{\nu}+1)\omega_{\eta}\omega_{\lambda}}\overset{\scriptscriptstyle{\text{(4)}}}{D}_{\mu\nu\eta\lambda}+$
$\displaystyle-\left(\frac{\delta_{\mu\eta}\delta_{\nu\lambda}+\delta_{\mu\lambda}\delta_{\nu\eta}}{2}\right)\frac{(2n_{\mu}n_{\nu}+n_{\mu}+n_{\nu})(2n_{\mu}n_{\nu}+n_{\mu}+n_{\nu}+1)2\omega_{\mu}\omega_{\nu}}{(2n_{\mu}+1)(2n_{\nu}+1)}$
(250)
${\bar{Y}}_{\mu\nu\eta\lambda}=-\frac{4(\delta_{\mu\eta}\delta_{\nu\lambda}+\delta_{\mu\lambda}\delta_{\nu\eta})\omega_{\mu}\omega_{\nu}}{(2n_{\mu}+1)(2n_{\nu}+1)}$
(251) $\displaystyle{\bar{Y}^{\prime}}_{\mu\nu\eta\lambda}$
$\displaystyle=\frac{\delta_{\mu\eta}\delta_{\nu\lambda}+\delta_{\mu\lambda}\delta_{\nu\eta}}{2}\cdot$
$\displaystyle\cdot\left(\frac{2\omega_{\mu}\omega_{\nu}}{(2n_{\mu}+1)(2n_{\nu}+1)}-\omega_{\mu}^{2}-\omega_{\nu}^{2}\right)$
(252)
${\bar{Z}}_{\mu\nu\eta}=\frac{2[(2n_{\mu}+1)\omega_{\nu}+(2n_{\nu}+1)\omega_{\mu}]}{\hbar(2n_{\mu}+1)(2n_{\nu}+1)}\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\mu\nu\eta}$
(253) $\displaystyle{\bar{Z}^{\prime}}_{\mu\nu\eta}$
$\displaystyle=\frac{2\omega_{\mu}n_{\mu}(n_{\mu}+1)(2n_{\nu}+1)}{\hbar(2n_{\mu}+1)(2n_{\nu}+1)}\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\mu\nu\eta}+$
$\displaystyle+\frac{2\omega_{\nu}n_{\nu}(n_{\nu}+1)(2n_{\mu}+1)}{\hbar(2n_{\mu}+1)(2n_{\nu}+1)}\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\mu\nu\eta}$
(254)
${\bar{X}^{\prime\prime}}_{\mu\nu\eta}=\frac{\hbar^{2}(2n_{\nu}+1)(2n_{\eta}+1)}{8\omega_{\nu}\omega_{\eta}}\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\mu\nu\eta}$
(255) ${\bar{Z}^{\prime\prime}}_{\mu\nu}=-\delta_{\mu\nu}\omega_{\mu}^{2}$
(256)

##  Appendix G Proof of the dynamical ansatz

In this appendix we proof that the self-energy calculated from the TD-SCHA in
Eq. (118) correspond to the dynamical _ansatz_ proposed by Bianco et al.[33],
reported in Eq. (119).

Our expression for the self-energy is:

$\displaystyle\bm{\Sigma}(\omega)=-\bar{\bm{X}^{\prime\prime}}$
$\displaystyle\left[(\bar{\bm{X}}+\omega^{2})-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{X}^{\prime}}\right]^{-1}\cdot$
$\displaystyle\cdot\left[\bar{\bm{Z}}-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{Z}^{\prime}}\right]$
(257)

To simplify the expression, we define the tensor $\bm{U}$, $\bm{P}$ and
$\bm{T}$, that in the polarization basis are:

$\bar{X^{\prime\prime}}_{\alpha\beta\gamma}=U_{\alpha\beta}\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\alpha\beta\gamma}$
(258)
$\bm{P}=(\bar{\bm{X}}+\omega^{2})-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{X}^{\prime}}$
(259)
$T_{\alpha\beta}=\frac{\bar{Z}_{\alpha\beta\gamma}}{\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\alpha\beta\gamma}}-{\bar{Y}}_{\alpha\beta}\left(\bar{Y^{\prime}}_{\alpha\beta}+\omega^{2}\right)^{-1}\frac{Z^{\prime}_{\alpha\beta\gamma}}{\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\alpha\beta\gamma}}$
(260)

Both $\bm{U}$ and $\bm{T}$ are 2-rank tensors, while $\bm{P}$ is in general a
4-rank tensor. However, if
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$ also $\bm{P}$ becomes a
2-rank tensor. For this reason it is easier to split $\bm{P}$ in two
contribution:

$\bm{P}={\bm{P}^{(0)}}+{\bm{P}^{(1)}}$ (261)

where ${\bm{P}^{(0)}}=\bm{P}$ if
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$.

With this expression, the self-energy becomes:

$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}{\bm{U}}{\bm{P}}^{-1}\bm{T}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
(262)

Here, the inversion of the 4-rank $\bm{P}$ tensor is equal to the inversion of
a big rank 2 tensor in which the first two and last two indices are grouped
together. Moreover, $\bm{T}$ has only two indices even if in the right-hand
expression of Eq. (260) three indices appears, as both
$\bar{Z}_{\alpha\beta\gamma}$ and $\bar{Z^{\prime}}_{\alpha\beta\gamma}$ loose
the dependence on the $\gamma$ index if divided by
$\overset{\scriptscriptstyle{\text{(3)}}}{D}_{\alpha\beta\gamma}$.

First of all, lets consider the simple case for which
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$ but
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\neq 0$. In this case the
Bianco self-energy reduces to the bubble diagram:

$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\bm{\Lambda}(\omega)\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
(263)

If we set $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$ in our Eq.
(262) we get:

$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}{\bm{U}}{\bm{P}^{(0)}}^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
(264)

The two equations are equal if:

$\bm{\Lambda}(\omega)={\bm{U}}{\bm{P}^{(0)}}^{-1}{\bm{T}}$ (265)

This equality can be proved by a simple algebric calculations. In fact,
$\bm{P}^{(0)}$ is a diagonal 4-rank tensor in the polarization basis, and can
be inverted by inverting its elements.

The explicit expression in the polarization basis of the $\bm{T}$,
$\bm{P}^{(0)}$ and $\bm{U}$ are:

$U_{\alpha\beta}=-\hbar^{2}\frac{(2n_{\alpha}+1)(2n_{\beta}+1)}{8\omega_{\alpha}\omega_{\beta}}$
(266)

$\displaystyle T_{\gamma\delta}$
$\displaystyle=\bigg{[}\frac{(4n_{\gamma}+2)\omega_{\delta}+(4n_{\delta}+2)\omega_{\gamma}}{\hbar(2n_{\gamma}+1)(2n_{\delta}+1)}+\frac{8\omega_{\gamma}\omega_{\delta}}{(2n_{\gamma}+1)(2n_{\delta}+1)}\frac{(2n_{\gamma}+1)(2n_{\delta}+1)}{(\omega^{2}-\omega_{\gamma}^{2}-\omega_{\delta}^{2})(2n_{\gamma}+1)(2n_{\delta}+1)+2\omega_{\gamma}\omega_{\delta}}\cdot$
$\displaystyle\cdot\left(\frac{2\omega_{\gamma}n_{\gamma}(n_{\gamma}+1)}{\hbar(2n_{\gamma}+1)}+\frac{2\omega_{\delta}n_{\delta}(n_{\delta}+1)}{\hbar(2n_{\delta}+1)}\right)\bigg{]}$
(267) $\displaystyle T_{\gamma\delta}$
$\displaystyle=\bigg{[}\frac{(4n_{\gamma}+2)\omega_{\delta}+(4n_{\delta}+2)\omega_{\gamma}}{\hbar(2n_{\gamma}+1)(2n_{\delta}+1)}+\frac{8\omega_{\gamma}\omega_{\delta}}{(\omega^{2}-\omega_{\gamma}^{2}-\omega_{\delta}^{2})(2n_{\gamma}+1)(2n_{\delta}+1)+2\omega_{\gamma}\omega_{\delta}}\left(\frac{2\omega_{\gamma}n_{\gamma}(n_{\gamma}+1)}{\hbar(2n_{\gamma}+1)}+\frac{2\omega_{\delta}n_{\delta}(n_{\delta}+1)}{\hbar(2n_{\delta}+1)}\right)\bigg{]}$
(268) $\displaystyle P_{\alpha\beta}^{(0)}$
$\displaystyle=\omega^{2}-\omega_{\alpha}^{2}-\omega_{\beta}^{2}-\frac{2\omega_{\alpha}\omega_{\beta}}{(2n_{\alpha}+1)(2n_{\beta}+1)}-\frac{8\omega_{\alpha}\omega_{\beta}}{(\omega^{2}-\omega_{\alpha}^{2}-\omega_{\beta}^{2})(2n_{\alpha}+1)(2n_{\beta}+1)+2\omega_{\alpha}\omega_{\beta}}\cdot$
$\displaystyle\cdot\frac{2\omega_{\alpha}\omega_{\beta}(2n_{\alpha}n_{\beta}+n_{\alpha}+n_{\beta})(2n_{\alpha}n_{\beta}+n_{\alpha}+n_{\beta}+1)}{(2n_{\alpha}+1)(2n_{\beta}+1)}$
(269) $\displaystyle P_{\alpha\beta}^{(0)}$
$\displaystyle=\omega^{2}-\omega_{\alpha}^{2}-\omega_{\beta}^{2}-\frac{2\omega_{\alpha}\omega_{\beta}}{(2n_{\alpha}+1)(2n_{\beta}+1)}\left[1+\frac{8\omega_{\alpha}\omega_{\beta}(2n_{\alpha}n_{\beta}+n_{\alpha}+n_{\beta})(2n_{\alpha}n_{\beta}+n_{\alpha}+n_{\beta}+1)}{(\omega^{2}-\omega_{\alpha}^{2}-\omega_{\beta}^{2})(2n_{\alpha}+1)(2n_{\beta}+1)+2\omega_{\alpha}\omega_{\beta}}\right]$
(270)

With some algebric manipulation, it is straightfoward to show that

$\displaystyle\frac{U_{\alpha\beta}T_{\alpha\beta}}{P^{(0)}_{\alpha\beta}}=-\frac{\hbar^{2}}{4\omega_{\alpha}\omega_{\beta}}\bigg{[}$
$\displaystyle\frac{(\omega_{\alpha}+\omega_{\beta})(n_{\alpha}+n_{\beta}+1)}{(\omega_{\alpha}+\omega_{\beta})^{2}-\omega^{2}}+$
$\displaystyle-\frac{(\omega_{\alpha}-\omega_{\beta})(n_{\alpha}-n_{\beta})}{(\omega_{\alpha}-\omega_{\beta})^{2}-\omega^{2}}\bigg{]}$
(271)

That is exactly the expression in the polarization basis of the
$\bm{\Lambda}(\omega)$ tensor. Therefore we proved Eq. (265), and the
dynamical ansatz of Bianco in the case
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$. To proceed with the case
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}\neq 0$ we must add the
$\bm{P}^{(1)}$.

$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}{\bm{U}}\left[{\bm{P}^{(0)}}+\bar{\bm{P}^{(1)}}\right]^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
(272)
$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}{\bm{U}}\left[\mathbbm{1}+{\bm{P}^{(0)}}^{-1}\bar{\bm{P}^{(1)}}\right]^{-1}{\bm{P}^{(0)}}^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
(273)

We already proved that

$\bm{\Lambda}(\omega)={\bm{U}}{\bm{P}^{(0)}}^{-1}{\bm{T}}$ (274)

Therefore, it is trivial to show that Eq. (273) is the Bianco self-energy (Eq.
119) if we prove that:

${\bm{P}^{(0)}}^{-1}\bar{\bm{P}^{(1)}}={\bm{P}^{(0)}}^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}{\bm{U}}.$
(275)

In fact, if we substitute Eq. (275) into Eq. (273), and we perform the Taylor
expansion for small $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$, we
get:

$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}{\bm{U}}{\bm{P}^{(0)}}^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}+\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}{\bm{U}}{\bm{P}^{(0)}}^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}{\bm{U}}{\bm{P}^{(0)}}^{-1}{\bm{T}}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}+\cdots$
(276)
$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\bm{\Lambda}(\omega)\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}+\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\bm{\Lambda}(\omega)\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}\bm{\Lambda}(\omega)\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}+\cdots$
(277)

That is the correct diagrammatic expansion of the self-energy, as defined in
ref.[33].

We give here the explicit expression of the $\bm{P}^{(1)}$ 4-rank tensor.

$\displaystyle P^{(1)}_{\alpha\beta\gamma\delta}$
$\displaystyle=-\frac{\hbar(2n_{\gamma}+1)(2n_{\delta}+1)(2\omega_{\alpha}n_{\beta}+2\omega_{\beta}n_{\alpha}+\omega_{\alpha}+\omega_{\beta})}{4(2n_{\alpha}+1)(2n_{\beta}+1)\omega_{\gamma}\omega_{\delta}}\overset{\scriptscriptstyle{\text{(4)}}}{D}_{\alpha\beta\gamma\delta}-\frac{8\omega_{\alpha}\omega_{\beta}}{(\omega^{2}-\omega_{\alpha}^{2}-\omega_{\beta}^{2})(2n_{\alpha}+1)(2n_{\beta}+1)+2\omega_{\alpha}\omega_{\beta}}\cdot$
$\displaystyle\cdot\frac{\hbar\left[\omega_{\alpha}n_{\alpha}(n_{\alpha}+1)(2n_{\beta}+1)+\omega_{\beta}n_{\beta}(n_{\beta}+1)(2n_{\alpha}+1)\right](2n_{\gamma}+1)(2n_{\delta}+1)}{4(2n_{\alpha}+1)(2n_{\beta}+1)\omega_{\gamma}\omega_{\delta}}\overset{\scriptscriptstyle{\text{(4)}}}{D}_{\alpha\beta\gamma\delta}$
(278)
$P^{(1)}_{\alpha\beta\gamma\delta}=-\frac{\hbar(2n_{\gamma}+1)(2n_{\delta}+1)\overset{\scriptscriptstyle{\text{(4)}}}{D}_{\alpha\beta\gamma\delta}}{4(2n_{\alpha}+1)(2n_{\beta}+1)\omega_{\gamma}\omega_{\delta}}\left[2\omega_{\alpha}n_{\beta}+2\omega_{\beta}n_{\alpha}+\omega_{\alpha}+\omega_{\beta}+\frac{8\omega_{\alpha}\omega_{\beta}[\omega_{\alpha}n_{\alpha}(n_{\alpha}+1)(2n_{\beta}+1)+\omega_{\beta}n_{\beta}(n_{\beta}+1)(2n_{\alpha}+1)]}{(\omega^{2}-\omega_{\alpha}^{2}-\omega_{\beta}^{2})(2n_{\alpha}+1)(2n_{\beta}+1)+2\omega_{\alpha}\omega_{\beta}}\right]$
(279)

Again, with straightfoward algebra we can prove that:

$\frac{P^{(1)}_{\alpha\beta\gamma\delta}}{P^{(0)}_{\alpha\beta}}=\frac{T_{\alpha\beta}\overset{\scriptscriptstyle{\text{(4)}}}{D}_{\alpha\beta\gamma\delta}U_{\gamma\delta}}{P^{(0)}_{\alpha\beta}}$
(280)

This concludes the proof that the one-phonon self-energy proposed as _ansatz_
by Bianco et al.[33] can be formally derived in a full dynamical contex within
the TD-SCHA theory.

##  Appendix H Additional details on the Lanczos algorithm

The Lanczos matrix $\mathcal{L}$ depends on the static SCHA self-consistent
frequencies $\omega_{\mu}$ that appears in the expressions of the Lanczos
matrix $\mathcal{L}$ (Eq. 137, 249,250,251,252,253,254,255,256) and on the
anharmonic scattering tensors
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ defined in Eq. (84) and
(85). The frequencies $\omega_{\mu}$ are obtained through the SCHA
minimization[20, 37], while $\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
and $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ are computed using only
first derivatives of the BO potential (i.e. _ab initio_ forces) with the
method described in ref.[33].

In particular they are evaluated with a Monte Carlo average on an ensemble of
$N_{\text{conf}}$ configurations (see Appendix F and ref.[33]).

For the simulation of hydrogen reported in Sec. VII.2, we employed an ensemble
of 40000 configurations.

The Lanczos algorithm works computing the tri-diagonal elements of the
$\mathcal{T}$ matrix ($\alpha_{i},\beta_{i},\gamma_{i}$). Each step we apply
the $\mathcal{L}$ matrix (Eq. 137) to a vector, and we obtain the next values
of $\alpha_{i}$,$\beta_{i}$ and $\gamma_{i}$. These values are used to obtain
the Green function with Eq. (152). As can be seen from the expression of the
Green function, each new iteration introduces a new pole in the Green
function. This is well represented in Figure 6 of the main text.

To see the signal coming from the poles of the spectral function, we need to
introduce an imaginary part on the Green function. This is done in two ways.
The standard way is introducing a finite smearing $\eta$ by changing
$\omega\rightarrow\omega+i\eta$. Also, we can employ the “terminator” of Eq.
(155). It introduces a branch cut in the Green function that smooth the
overall result non uniformly, allowing to converge with a much lower number of
iterations.

We present in Figure 7 the result of the Raman signal of high-pressure
hydrogen presented in Sec. VII.2, with and without the terminator
($N_{\text{steps}}=200$, $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$,
$\eta=$5\text{\,}{\mathrm{cm}}^{-1}$$).

As shown, the low region of the spectrum is completely unaffected by the
terminator. In fact, each new iterations adds a new pole on the vibron, whose
intensity is distributed in several poles. For this reason, converging a
complex non Lorentzian shape with the continued fraction without terminator
requires a lot of iterations. The terminator fixes this behaviour providing a
smooth lineshape for the whole vibron, without affecting the rest of the
spectrum. A similar effect can be obtained by increasing the smearing $\eta$
to a value larger than the distance between two poles. However, in this way,
$\eta$ will affect the full spectrum, even the low frequency region.

##  Appendix I Computational details for the hydrogen simulation

The simulations are performed on a 2x2x1 supercell to sample the Brillouin
zone for phonons (96 atoms). The _ab initio_ energy landscape is simulated
through DFT, with the BLYP[62] functional. For all the DFT calculations, we
employed the Quantum ESPRESSO[63, 64] package, with a plane wave basis set and
a norm conserving pseudo-potential from the Pseudo Dojo library[65]. The
energy cutoff for the basis of the wave-functions was 60 Ry (240 Ry for the
electronic density). The Brillouin zone for electrons is sampled on a 4x4x4
mesh in the phonon displaced supercell (2x2x1). Additional details can be
found in ref.[26], where we published the Raman and IR signal for several
pressures and focused the discussion on the high-pressure hydrogen physics.

##  Appendix J Harmonic two-phonon propagator

In this section we derive the Harmonic two-phonon propagator from the response
system of Eq. (81c).

The harmonic limit is obtained as
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}=0$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$. In this limit, the one-
phonon and two-phonon Green functions are decoupled. The two-phonon
propagation is described by the variables $\mathrm{Re}\bm{A}$ and
$\bm{\Upsilon}$, so we isolate these contributions in the green function (Eq.
(88)).

The two phonon green function is obtained inverting the matrix:

$G^{(2ph)}(\omega)=\begin{pmatrix}-\bar{\bm{X}}-\omega^{2}&-\bar{\bm{Y}}\\\
-{\bar{\bm{X}^{\prime}}}&-\bar{\bm{Y}^{\prime}}-\omega^{2}\end{pmatrix}^{-1}$
(281)

In the polarization basis, the tensor $\bar{\bm{X}}$, $\bar{\bm{Y}}$,
$\bar{\bm{X}^{\prime}}$ and $\bar{\bm{Y}^{\prime}}$ are diagonal. Therefore,
once we identify a couple of modes $\mu\nu$, we only need to invert a 2x2
matrix in Eq. (281).

Since we are interested in the fluctuations-fluctuations correlation function,
we only need the response in the $\bm{\Upsilon}$ block (the fluctuations are
the convariance matrix $\bm{\Upsilon}^{-1}$). This term is:

$G^{(\Upsilon)}_{\mu\nu}(\omega)=\frac{-\bar{Y^{\prime}}_{\mu\nu\mu\nu}-\omega^{2}}{\bar{X^{\prime}}_{\mu\nu\mu\nu}\bar{Y}_{\mu\nu\mu\nu}+(\bar{X}_{\mu\nu\mu\nu}+\omega^{2})(\bar{Y^{\prime}}_{\mu\nu\mu\nu}+\omega^{2})}$
(282)

Subsittuting the expression defined in F (taking care of the symmetry
exchanging $\mu\leftrightarrow\nu$) we get:

$\displaystyle G^{(\Upsilon)}_{ab}(\omega)$
$\displaystyle=\frac{2}{\hbar\left(4{n}_{a}{n}_{b}+2{n}_{a}+2{n}_{b}+1\right)}\cdot$
$\displaystyle\cdot\bigg{[}\frac{-2\omega^{2}{\omega}_{a}{n}_{b}-\omega^{2}{\omega}_{a}-2\omega^{2}{\omega}_{b}{n}_{a}-\omega^{2}{\omega}_{b}}{[(\omega_{a}-\omega_{b})^{2}-\omega^{2}][(\omega_{a}+\omega_{b})^{2}-\omega^{2}]}+$
$\displaystyle+\frac{2{\omega}_{a}^{3}{n}_{b}+{\omega}_{a}^{3}-2{\omega}_{a}^{2}{\omega}_{b}{n}_{a}-{\omega}_{a}^{2}{\omega}_{b}}{[(\omega_{a}-\omega_{b})^{2}-\omega^{2}][(\omega_{a}+\omega_{b})^{2}-\omega^{2}]}+$
$\displaystyle+\frac{-2{\omega}_{a}{\omega}_{b}^{2}{n}_{b}-{\omega}_{a}{\omega}_{b}^{2}+2{\omega}_{b}^{3}{n}_{a}+{\omega}_{b}^{3}}{[(\omega_{a}-\omega_{b})^{2}-\omega^{2}][(\omega_{a}+\omega_{b})^{2}-\omega^{2}]}\bigg{]}$
(283)

This green function has the poles in the correct position, when:

$\omega^{2}=(\omega_{a}\pm\omega_{b})^{2}$

Since we want the fluctuation-fluctuation correlation function, we need to
change variable for the Green function from $\bm{\Upsilon}$ to
$\bm{\Upsilon}^{-1}$. This is achieved knowing how the perturbation in
$(\bm{\Upsilon}^{-1})^{(1)}$ depend on those on $\bm{\Upsilon}^{(1)}$:

$(\bm{\Upsilon}^{-1})^{(1)}_{ab}=-\sum_{cd}\Upsilon^{-1}_{ac}\Upsilon^{-1}_{bd}\Upsilon^{(1)}_{cd}$
(284)

From this expression we can compute the fluctuation-fluctuation response
function:

$\chi_{abcd}(\omega)=G^{(\Upsilon^{-1})}_{abcd}(\omega)=-\sum_{ef}\Upsilon^{-1}_{ae}\Upsilon^{-1}_{bf}G^{(\Upsilon)}_{efcd}(\omega)$
(285)

If we substitute the expressions in the polarization basis, we obtain:

$\displaystyle\chi_{\mu\nu}(\omega)$
$\displaystyle=-\frac{\hbar\left({\omega}_{\mu}-{\omega}_{\nu}\right)\left({n}_{\mu}-{n}_{\nu}\right)}{2\left(\omega-{\omega}_{\mu}+{\omega}_{\nu}\right)\left(\omega+{\omega}_{\mu}-{\omega}_{\nu}\right){\omega}_{\mu}{\omega}_{\nu}}+$
$\displaystyle+\frac{\hbar\left({\omega}_{\mu}+{\omega}_{\nu}\right)\left({n}_{\mu}+{n}_{\nu}+1\right)}{2\left(\omega^{2}-\left({\omega}_{\mu}+{\omega}_{\nu}\right)^{2}\right){\omega}_{\mu}{\omega}_{\nu}}$
(286)

That is the standard two-phonon propagator.

## References

  * Li  _et al._ [2014] Y. Li, J. Hao, H. Liu, Y. Li, and Y. Ma, The metallization and superconductivity of dense hydrogen sulfide, The Journal of Chemical Physics 140, 174712 (2014).
  * Peng  _et al._ [2017] F. Peng, Y. Sun, C. J. Pickard, R. J. Needs, Q. Wu, and Y. Ma, Hydrogen clathrate structures in rare earth hydrides at high pressures: Possible route to room-temperature superconductivity, Physical Review Letters 119, 10.1103/physrevlett.119.107001 (2017).
  * Drozdov  _et al._ [2015] A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system, Nature 525, 73 (2015).
  * Somayazulu  _et al._ [2019] M. Somayazulu, M. Ahart, A. K. Mishra, Z. M. Geballe, M. Baldini, Y. Meng, V. V. Struzhkin, and R. J. Hemley, Evidence for superconductivity above 260 k in lanthanum superhydride at megabar pressures, Physical Review Letters 122, 10.1103/physrevlett.122.027001 (2019).
  * Drozdov  _et al._ [2019] A. Drozdov, P. Kong, V. Minkov, S. Besedin, M. Kuzovnikov, S. Mozaffari, L. Balicas, F. Balakirev, D. Graf, V. Prakapenka, _et al._ , Superconductivity at 250 k in lanthanum hydride under high pressures, Nature 569, 528 (2019).
  * Wang  _et al._ [2010] Y. Wang, J. Lv, L. Zhu, and Y. Ma, Crystal structure prediction via particle-swarm optimization, Physical Review B 82, 10.1103/physrevb.82.094116 (2010).
  * Oganov and Glass [2006] A. R. Oganov and C. W. Glass, Crystal structure prediction usingab initioevolutionary techniques: Principles and applications, The Journal of Chemical Physics 124, 244704 (2006).
  * Pickard and Needs [2011] C. J. Pickard and R. J. Needs, Ab initiorandom structure searching, Journal of Physics: Condensed Matter 23, 053201 (2011).
  * Runge and Gross [1984] E. Runge and E. K. U. Gross, Density-functional theory for time-dependent systems, Physical Review Letters 52, 997 (1984).
  * Marques  _et al._ [2006] M. A. Marques, C. A. Ullrich, F. Nogueira, K. Burke, A. Rubio, and E. K. U. Gross,  _Time-Dependent Density Functional Theory_ (Springer Science & Business Media, 2006).
  * Car and Parrinello [1985] R. Car and M. Parrinello, Unified approach for molecular dynamics and density-functional theory, Physical review letters 55, 2471 (1985).
  * Cupo  _et al._ [2019] A. Cupo, D. Tristant, K. Rego, and V. Meunier, Theoretical analysis of spectral lineshapes from molecular dynamics, npj Computational Materials 5, 10.1038/s41524-019-0220-1 (2019).
  * Baym and Mermin [1961] G. Baym and N. D. Mermin, Determination of thermodynamic green's functions, Journal of Mathematical Physics 2, 232 (1961).
  * Krilov and Berne [1999] G. Krilov and B. J. Berne, Real time quantum correlation functions. II. maximum entropy numerical analytic continuation of path integral monte carlo and centroid molecular dynamics data, The Journal of Chemical Physics 111, 9147 (1999).
  * Poulsen  _et al._ [2003] J. A. Poulsen, G. Nyman, and P. J. Rossky, Practical evaluation of condensed phase quantum correlation functions: A feynman–kleinert variational linearized path integral method, The Journal of Chemical Physics 119, 12179 (2003).
  * Hernandez and Voth [1998] R. Hernandez and G. A. Voth, Quantum time correlation functions and classical coherence, Chemical Physics 233, 243 (1998).
  * Sun  _et al._ [1998] X. Sun, H. Wang, and W. H. Miller, Semiclassical theory of electronically nonadiabatic dynamics: Results of a linearized approximation to the initial value representation, The Journal of Chemical Physics 109, 7064 (1998).
  * Bonella  _et al._ [2010] S. Bonella, M. Monteferrante, C. Pierleoni, and G. Ciccotti, Path integral based calculations of symmetrized time correlation functions. II, The Journal of Chemical Physics 133, 164105 (2010).
  * Snider  _et al._ [2020] E. Snider, N. Dasenbrock-Gammon, R. McBride, M. Debessai, H. Vindana, K. Vencatasamy, K. V. Lawler, A. Salamat, and R. P. Dias, Room-temperature superconductivity in a carbonaceous sulfur hydride, Nature 586, 373 (2020).
  * Errea  _et al._ [2014] I. Errea, M. Calandra, and F. Mauri, Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides, Physical Review B 89, 10.1103/physrevb.89.064302 (2014).
  * Errea  _et al._ [2016] I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Quantum hydrogen-bond symmetrization in the superconducting hydrogen sulfide system, Nature 532, 81 (2016).
  * Bianco  _et al._ [2018] R. Bianco, I. Errea, M. Calandra, and F. Mauri, High-pressure phase diagram of hydrogen and deuterium sulfides from first principles: Structural and vibrational properties including quantum and anharmonic effects, Physical Review B 97, 10.1103/physrevb.97.214101 (2018).
  * Errea  _et al._ [2020] I. Errea, F. Belli, L. Monacelli, A. Sanna, T. Koretsune, T. Tadano, R. Bianco, M. Calandra, R. Arita, F. Mauri, and J. A. Flores-Livas, Quantum crystal structure in the 250-kelvin superconducting lanthanum hydride, Nature 578, 66 (2020).
  * Borinaga  _et al._ [2016a] M. Borinaga, P. Riego, A. Leonardo, M. Calandra, F. Mauri, A. Bergara, and I. Errea, Anharmonic enhancement of superconductivity in metallic molecularCmca - 4 hydrogen at high pressure: a first-principles study, Journal of Physics: Condensed Matter 28, 494001 (2016a).
  * Borinaga  _et al._ [2016b] M. Borinaga, I. Errea, M. Calandra, F. Mauri, and A. Bergara, Anharmonic effects in atomic hydrogen: Superconductivity and lattice dynamical stability, Phys. Rev. B 93, 174308 (2016b).
  * Monacelli  _et al._ [2020] L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Black metal hydrogen above 360 GPa driven by proton quantum fluctuations, Nature Physics 10.1038/s41567-020-1009-3 (2020).
  * Bianco  _et al._ [2019] R. Bianco, I. Errea, L. Monacelli, M. Calandra, and F. Mauri, Quantum enhancement of charge density wave in NbS2 in the two-dimensional limit, Nano Letters 19, 3098 (2019).
  * Bianco  _et al._ [2020] R. Bianco, L. Monacelli, M. Calandra, F. Mauri, and I. Errea, Weak dimensionality dependence and dominant role of ionic fluctuations in the charge-density-wave transition of NbSe2, Physical Review Letters 125, 10.1103/physrevlett.125.106101 (2020).
  * Zhou  _et al._ [2020] J. S. Zhou, L. Monacelli, R. Bianco, I. Errea, F. Mauri, and M. Calandra, Anharmonicity and doping melt the charge density wave in single-layer TiSe2, Nano Letters 20, 4809 (2020).
  * Aseginolaza  _et al._ [2019a] U. Aseginolaza, R. Bianco, L. Monacelli, L. Paulatto, M. Calandra, F. Mauri, A. Bergara, and I. Errea, Phonon collapse and second-order phase transition in thermoelectric SnSe, Physical Review Letters 122, 10.1103/physrevlett.122.075901 (2019a).
  * Aseginolaza  _et al._ [2019b] U. Aseginolaza, R. Bianco, L. Monacelli, L. Paulatto, M. Calandra, F. Mauri, A. Bergara, and I. Errea, Strong anharmonicity and high thermoelectric efficiency in high temperature sns from first-principles, arXiv preprint arXiv:1906.02047 (2019b).
  * Ribeiro  _et al._ [2018] G. A. S. Ribeiro, L. Paulatto, R. Bianco, I. Errea, F. Mauri, and M. Calandra, Strong anharmonicity in the phonon spectra of PbTe and SnTe from first principles, Physical Review B 97, 10.1103/physrevb.97.014306 (2018).
  * Bianco  _et al._ [2017] R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Mauri, Second-order structural phase transitions, free energy curvature, and temperature-dependent anharmonic phonons in the self-consistent harmonic approximation: Theory and stochastic implementation, Physical Review B 96, 10.1103/physrevb.96.014111 (2017).
  * Lazzeri and Mauri [2003] M. Lazzeri and F. Mauri, First-principles calculation of vibrational raman spectra in large systems: Signature of small rings in CrystallineSiO2, Physical Review Letters 90, 10.1103/physrevlett.90.036401 (2003).
  * Simoncelli  _et al._ [2019] M. Simoncelli, N. Marzari, and F. Mauri, Unified theory of thermal transport in crystals and glasses, Nature Physics 15, 809 (2019).
  * Monacelli  _et al._ [2017] L. Monacelli, G. Batignani, G. Fumero, C. Ferrante, S. Mukamel, and T. Scopigno, Manipulating impulsive stimulated raman spectroscopy with a chirped probe pulse, The Journal of Physical Chemistry Letters 8, 966 (2017).
  * Monacelli  _et al._ [2018] L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Pressure and stress tensor of complex anharmonic crystals within the stochastic self-consistent harmonic approximation, Phys. Rev. B 98, 024106 (2018).
  * cheng Li and qing Tong [1985] T. cheng Li and P. qing Tong, Hohenberg-kohn theorem for time-dependent ensembles, Physical Review A 31, 1950 (1985).
  * Huber and Heller [1987] D. Huber and E. J. Heller, Generalized gaussian wave packet dynamics, The Journal of Chemical Physics 87, 5302 (1987).
  * Huber  _et al._ [1988] D. Huber, E. J. Heller, and R. G. Littlejohn, Generalized gaussian wave packet dynamics, schrödinger equation, and stationary phase approximation, The Journal of Chemical Physics 89, 2003 (1988).
  * Pal  _et al._ [2016] H. Pal, M. Vyas, and S. Tomsovic, Generalized gaussian wave packet dynamics: Integrable and chaotic systems, Physical Review E 93, 10.1103/physreve.93.012213 (2016).
  * Haegeman  _et al._ [2011] J. Haegeman, J. I. Cirac, T. J. Osborne, I. Pižorn, H. Verschelde, and F. Verstraete, Time-dependent variational principle for quantum lattices, Physical Review Letters 107, 10.1103/physrevlett.107.070601 (2011).
  * Guaita  _et al._ [2019] T. Guaita, L. Hackl, T. Shi, C. Hubig, E. Demler, and J. I. Cirac, Gaussian time-dependent variational principle for the bose-hubbard model, Physical Review B 100, 10.1103/physrevb.100.094529 (2019).
  * Hackl  _et al._ [2020] L. Hackl, T. Guaita, T. Shi, J. Haegeman, E. Demler, and J. I. Cirac, Geometry of variational methods: dynamics of closed quantum systems, SciPost Phys. 9, 48 (2020).
  * Meyer and Worth [2003] H.-D. Meyer and G. A. Worth, Quantum molecular dynamics: propagating wavepackets and density operators using the multiconfiguration time-dependent hartree method, Theoretical Chemistry Accounts 109, 251 (2003).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2007).
  * Kohn and Luttinger [1957] W. Kohn and J. M. Luttinger, Quantum theory of electrical transport phenomena, Physical Review 108, 590 (1957).
  * Pribram-Jones  _et al._ [2016] A. Pribram-Jones, P. E. Grabowski, and K. Burke, Thermal density functional theory: Time-dependent linear response and approximate functionals from the fluctuation-dissipation theorem, Physical Review Letters 116, 10.1103/physrevlett.116.233001 (2016).
  * Errea  _et al._ [2012] I. Errea, B. Rousseau, and A. Bergara, Isotope effect in the superconducting high-pressure simple cubic phase of calcium from first principles, Journal of Applied Physics 111, 112604 (2012).
  * Hellman  _et al._ [2011] O. Hellman, I. A. Abrikosov, and S. I. Simak, Lattice dynamics of anharmonic solids from first principles, Phys. Rev. B 84, 180301 (2011).
  * Tadano  _et al._ [2014] T. Tadano, Y. Gohda, and S. Tsuneyuki, Anharmonic force constants extracted from first-principles molecular dynamics: applications to heat transfer simulations, Journal of Physics: Condensed Matter 26, 225402 (2014).
  * Lihm and Park [2020] J.-M. Lihm and C.-H. Park, Gaussian time-dependent variational principle for the finite-temperature anharmonic lattice dynamics, arXiv preprint arXiv:2010.15725 (2020).
  * Putrino and Parrinello [2002] A. Putrino and M. Parrinello, Anharmonic raman spectra in high-pressure ice fromAb InitioSimulations, Physical Review Letters 88, 10.1103/physrevlett.88.176401 (2002).
  * Silvestrelli  _et al._ [1997] P. L. Silvestrelli, M. Bernasconi, and M. Parrinello, Ab initio infrared spectrum of liquid water, Chemical Physics Letters 277, 478 (1997).
  * Windl  _et al._ [1993] W. Windl, P. Pavone, K. Karch, O. Schütt, D. Strauch, P. Giannozzi, and S. Baroni, Second-order raman spectra of diamond fromab initiophonon calculations, Physical Review B 48, 3164 (1993).
  * Press [2017] W. H. Press,  _Numerical Recipes: The Art of Scientific Computing, Third Edition_ , 3rd ed. (Cambridge University Press, 2017).
  * Rocca  _et al._ [2008] D. Rocca, R. Gebauer, Y. Saad, and S. Baroni, Turbo charging time-dependent density-functional theory with lanczos chains, The Journal of Chemical Physics 128, 154105 (2008).
  * Errea  _et al._ [2015] I. Errea, M. Calandra, C. J. Pickard, J. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, High-pressure hydrogen sulfide from first principles: A strongly anharmonic phonon-mediated superconductor, Physical Review Letters 114, 10.1103/physrevlett.114.157004 (2015).
  * Errea [2016] I. Errea, Approaching the strongly anharmonic limit with ab initio calculations of materials’ vibrational properties – a colloquium, The European Physical Journal B 89 (2016).
  * Magdău and Ackland [2013] I. B. Magdău and G. J. Ackland, Identification of high-pressure phases III and IV in hydrogen: Simulating raman spectra using molecular dynamics, Physical Review B 87, 10.1103/physrevb.87.174110 (2013).
  * Zhang  _et al._ [2018] C. Zhang, C. Zhang, M. Chen, W. Kang, Z. Gu, J. Zhao, C. Liu, C. Sun, and P. Zhang, Finite-temperature infrared and raman spectra of high-pressure hydrogen from first-principles molecular dynamics, Physical Review B 98, 10.1103/physrevb.98.144301 (2018).
  * Miehlich  _et al._ [1989] B. Miehlich, A. Savin, H. Stoll, and H. Preuss, Results obtained with the correlation energy density functionals of becke and lee, yang and parr, Chemical Physics Letters 157, 200 (1989).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Advanced capabilities for materials modelling with quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving vanderbilt pseudopotentials, Phys. Rev. B 88, 085117 (2013).