<figure><img src="image_1.png" /></figure>

HAL
open science

# Beyond Markovian Dissipation at the Nanoscale :
towards finding quantum design rules for bio-organic
nanodevices

Thibaut Lacroix

► To cite this version:

Thibaut Lacroix. Beyond Markovian Dissipation at the Nanoscale: towards finding quantum design rules for bio-organic nanodevices. Quantum Physics [quant-ph]. Sorbonne Université; University of Saint Andrews, 2023. English. NNT: 2023SORUS106. tel-04135149

HAL Id: tel-04135149

https://theses.hal.science/tel-04135149

Submitted on 20 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Thèse de Doctorat

en vue de l'obtention du grade de **Docteur** délivré par Sorbonne Université
in partial fulfilment for the degree of **Doctor of Philosophy (PhD)** at the
University of St Andrews

Discipline : **Physique**

École Doctorale n°564 Physique en Île de France

Quantum Materials Centre for Doctoral Training

Réalisée à l'**Institut des NanoSciences de Paris & School of Physics and Astronomy**
sous la direction de Alex CHIN et Brendon LOVETT

Présentée et soutenue publiquement le 31 mars 2023 par

**Thibaut LACROIX**

**Beyond Markovian Dissipation at the Nanoscale**

*Towards finding quantum design rules for bio-organic nanodevices*

**Au delà de la dissipation markovienne à l'échelle nanométrique**

*Vers la découverte de règles quantiques pour la conception de nano-dispositifs bio-organiques*

Devant le jury composé de

Pr Fabienne

MICHELINI

Professeur

Aix-Marseille Université

Rapporteuse

Dr Stephen

CLARK

Associate Professor

University of Bristol

Rapporteur

Pr Natalia

KOROLKOVA

Professor

University of St Andrews

Examinatrice

Dr Riccardo

SPEZIA

Directeur de recherche

LCT | CNRS

Président du jury

Dr Alex

CHIN

Chargé de recherche

INSP | CNRS

Co-directeur de thèse

Pr Brendon

LOVETT

Professor

University of St Andrews

Co-directeur de thèse
<figure><img src="image_2.png" /></figure>
# Table of Contents

<table><tr><td>Acknowledgements</td><td>ix</td></tr><tr><td>Publications</td><td>xi</td></tr><tr><td>Abstract</td><td>xii</td></tr><tr><td>1 Introduction</td><td>1</td></tr><tr><td>1.1 Opening quantum physics</td><td>1</td></tr><tr><td>1.2 Inspirations from biology</td><td>3</td></tr><tr><td>1.3 Outline</td><td>6</td></tr><tr><td>BACKGROUND</td><td>9</td></tr><tr><td>2 Theory of Open Quantum Systems</td><td>11</td></tr><tr><td>2.1 Isolated, closed and open systems</td><td>11</td></tr><tr><td>2.2 Environment, bath and reservoir</td><td>13</td></tr><tr><td>2.3 Time evolution</td><td>13</td></tr><tr><td>2.4 Harmonic Bath, correlation function and spectral density</td><td>15</td></tr><tr><td>2.5 (Non-)Markovianity</td><td>19</td></tr><tr><td>METHODS</td><td>23</td></tr><tr><td>3 Tensor Networks</td><td>25</td></tr><tr><td>3.1 Quantum states, tensors and diagrammatic representation</td><td>25</td></tr><tr><td>3.2 Matrix Product States</td><td>26</td></tr><tr><td>3.3 Singular Value Decomposition</td><td>28</td></tr><tr><td>Restricted rank approximation of quantum states</td><td>30</td></tr><tr><td>Gauge freedom and canonical forms</td><td>30</td></tr><tr><td>3.4 Matrix Product Operators</td><td>33</td></tr><tr><td>Constructing a MPO for a Hamiltonian</td><td>34</td></tr><tr><td>3.5 Time-evolution methods</td><td>36</td></tr><tr><td>Time-Dependent Variational Principle</td><td>38</td></tr><tr><td>4 Chain Mapping of Bosonic Environments</td><td>43</td></tr><tr><td>4.1 Principle of the chain mapping</td><td>44</td></tr><tr><td>4.2 Extension to finite temperature</td><td>47</td></tr><tr><td>4.3 Spatially extended systems</td><td>50</td></tr><tr><td>Zero temperature</td><td>51</td></tr><tr><td>Finite temperature</td><td>53</td></tr><tr><td>4.A Finite-temperature extended bath and thermofields</td><td>56</td></tr><tr><td>Main results of thermofields</td><td>56</td></tr><tr><td>Alternative derivation from the bath correlation function</td><td>57</td></tr><tr><td>Relation with the extended bath</td><td>58</td></tr></table>
<table><tr><td>5</td><td>Influence of the ordering of the chains modes</td><td>59</td></tr><tr><td>5.1</td><td>How to manage multiple environments?</td><td>59</td></tr><tr><td>5.2</td><td>Two environments: Mode arrangements</td><td>60</td></tr><tr><td></td><td>Left-Right arrangement</td><td>62</td></tr><tr><td></td><td>Successive arrangement</td><td>62</td></tr><tr><td></td><td>Interleaved arrangement</td><td>63</td></tr><tr><td></td><td>Comparison</td><td>63</td></tr><tr><td>5.3</td><td>Three environments: Connectivity</td><td>68</td></tr><tr><td>5.4</td><td>Discussion</td><td>70</td></tr></table>

APPLICATIONS 73

<table><tr><td>6</td><td>Unveiling non-Markovian Signalling</td><td>75</td></tr><tr><td>6.1</td><td>Model</td><td>76</td></tr><tr><td></td><td>Hamiltonian MPO formulation</td><td>78</td></tr><tr><td></td><td>Long- and short-distance limits</td><td>79</td></tr><tr><td></td><td>Mapping to SBM for N = 2</td><td>80</td></tr><tr><td>6.2</td><td>Non-Markovian recurrences and information backflow</td><td>81</td></tr><tr><td></td><td>Zero temperature</td><td>81</td></tr><tr><td></td><td>Finite temperature</td><td>86</td></tr><tr><td>6.3</td><td>Discussion</td><td>88</td></tr></table>

<table><tr><td>7</td><td>Dissipation in multi-component Nanodevices</td><td>91</td></tr><tr><td>7.1</td><td>Quantum switch</td><td>94</td></tr><tr><td></td><td>Displaced bath in the chain representation</td><td>96</td></tr><tr><td>7.2</td><td>Reorganisation dynamics</td><td>97</td></tr><tr><td>7.3</td><td>Remote transient activation</td><td>100</td></tr><tr><td></td><td>Sign of the interaction Hamiltonian</td><td>101</td></tr><tr><td>7.4</td><td>Remote permanent activation</td><td>102</td></tr><tr><td>7.5</td><td>Energy landscape texture</td><td>104</td></tr><tr><td>7.6</td><td>Analogies with allosteric regulation</td><td>105</td></tr><tr><td>7.7</td><td>Discussion</td><td>107</td></tr><tr><td>7.A</td><td>Definition of the reorganisation energy with a toy example</td><td>108</td></tr></table>

8 Conclusions and Prospects 109

Bibliography 113
**Candidate's declaration**

I, Thibaut François Marie Lacroix, do hereby certify that this thesis, submitted for the degree of PhD, which is approximately 29,529 words in length, has been written by me, and that it is the record of work carried out by me, or principally by myself in collaboration with others as acknowledged, and that it has not been submitted in any previous application for any degree. I confirm that any appendices included in my thesis contain only material permitted by the 'Assessment of Postgraduate Research Students' policy.

I was admitted as a research student at the University of St Andrews in September 2019.

I received funding from an organisation or institution and have acknowledged the funder(s) in the full text of my thesis.

Date Signature of candidate

03/02/2023

<figure><img src="image_3.png" /></figure>

**Supervisor's declaration**

I hereby certify that the candidate has fulfilled the conditions of the Resolution and Regulations appropriate for the degree of PhD in the University of St Andrews and that the candidate is qualified to submit this thesis in application for that degree. I confirm that any appendices included in the thesis contain only material permitted by the 'Assessment of Postgraduate Research Students' policy.

Date Signature of supervisor

9/2/23

<figure><img src="image_4.png" /></figure>

**Permission for publication**

In submitting this thesis to the University of St Andrews we understand that we are giving permission for it to be made available for use in accordance with the regulations of the University Library for the time being in force, subject to any copyright vested in the work not being affected thereby. We also understand, unless exempt by an award of an embargo as requested below, that the title and the abstract will be published, and that a copy of the work may be made and supplied to any bona fide library or research worker, that this thesis will be electronically accessible for personal or research use and that the library has the right to migrate this thesis into new electronic forms as required to ensure continued access to the thesis.

I, Thibaut François Marie Lacroix, confirm that my thesis does not contain any third-party material that requires copyright clearance.

The following is an agreed request by candidate and supervisor regarding the publication of this thesis:
**Electronic copy**

No embargo on electronic copy.

Date Signature of candidate

03/02/2023

<figure><img src="image_5.png" /></figure>

Date Signature of supervisor

<figure><img src="image_6.png" /></figure>

9/2/23

**Underpinning Research Data or Digital Outputs**

**Candidate's declaration**

I, Thibaut Francois Marie Lacroix, hereby certify that no requirements to deposit original research data or digital outputs apply to this thesis and that, where appropriate, secondary data used have been referenced in the full text of my thesis.

Date Signature of candidate

03/02/2023

<figure><img src="image_7.png" /></figure>
Mon dessein n'est pas
d'amuser un public hautain,
mais d'attirer les suffrages de
l'amitié seule.
J'aurais voulu t'offrir une
œuvre plus digne de toi, plus
digne de ta belle âme,
étrangère à tout ce qui se passe
ici-bas, et remplie d'une poésie
vivante et pure, de pensées
simples et grandes!
Mais je n'ai à présenter à ton
accueil bienveillant et amical
que ces chants bizarres, moitié
plaisants et moitié tristes,
vulgaires et fantasques, fruits
tardifs de mes loisirs, de mes
légères inspirations, de mes
nuits sans sommeil, de mes
années jeunes et déjà flétries;
froides observations de mon
esprit et tristes remarques de
mon cœur.

Alexandre POUCHKINE,
*Eugène Onéguine*
<figure><img src="image_8.png" /></figure>
# Acknowledgements

Les remerciements occupent une place à part dans une thèse. Ceux-ci regroupent pêle-mêle parents, amis et collègues; toutes personnes ayant accompagné l'aspirant docteur lors de ses années de thèse. Alors que le reste du manuscrit relate un travail scientifique dans un style académique, les remerciements peuvent quant à eux aborder des moments de la vie du doctorant contemporain à sa thèse mais distincts. C'est l'unique section d'un ouvrage consacré à un travail de recherche qui laisse entrevoir l'épaisseur de l'existence de son auteur. Cette incongruité est peut-être l'unique trace visible de l'extérieur du fait qu'une thèse est un acte de compagnonnage lors duquel un apprenti réalise ses premières œuvres guidé par des pairs, des maîtres comme l'on disait autrefois.

This is why I want to first thank my supervisors Brendon W. Lovett and Alex W. Chin. You have both been wonderful supervisors on the scientific as well as on the personal level, and you are an inspiration for my scientific career. You guided me when I was taking my first steps into the real scientific world, helped me understand the science and regain vision when I was lost in my work, you also trusted me and gave me the freedom to explore my own research path, and you let me take the time necessary to take care of myself in troubled times. For all of that I am grateful.

Aline, je t'ai soutenue lorsque tu en avais besoin, tu m'as soutenu lorsque j'en avais besoin; je retiens les rires et les larmes, Paris, Saint Cloud, Saint Pétersbourg et aussi la Normandie, et je ne regrette rien.

À mes amis qui m'ont vus partir et revenir.

Les MécaQ avec qui j'ai grandi : Amaury, Félix, Clara, Fern, Enzo, Victor et Pierre, je suis finalement le seul à faire de la mécanique quantique ! Les agrégés : Tommy, je ne suis pas sûr d'être un futur collègue, Sylvain mon binôme et Christelle aka Kiristelle Boit Ce Lait, merci d'avoir toujours été là, quand ça allait bien et quand ça n'allait pas. Raphaelle, tant de choses depuis Lyon, Agoranov, le club de débat et jusqu'à aujourd'hui. Ta force de caractère et ton ambition (dans le bon sens du terme) m'inspirent. Laetitia, ma binôme ! Nous ne serons jamais d'accord sur rien à part l'humour mais nous serons toujours amis. Élise, qui sait naviguer dans un monde absurde comme personne, je retiens que lorsqu'on souhaite savoir quelque chose il suffit de poser la question. Louis, un autre "Lyonnais" irremplaçable, sémillant et humble, nous nous engageons ensemble sur les chemins accidentés de la recherche. Sophie, avec qui j'ai passé tant de temps à m'interroger sur les voies à emprunter, tu es la seule personne que je connaisse qui a plus de diplômes que moi. Marine qui aura toujours une place à part dans ma vie, même si nous nous voyons moins souvent je sais que tu seras toujours là.

Lorine, j'ai de la chance de t'avoir rencontrée et j'apprécie tes photos de lamas et de cactus.

Bien évidemment les SuperNanos de l'INSP ! Compagnons de thèse, de sorties et de mots fléchés. Thomas qui connaît tout le monde, Mathis toujours de bonne humeur et légèrement moqueur, Sarah et son énergie inépuisable, Marine retrouvée par hasard après l'agrég, Angeline toujours prête à en découdre, Thibault lui aussi toujours en chemise mais pas les mêmes, Benoît toujours combatif et là pour les autres, Lucille qui m'aide à préparer mes blagues, Brieuc nouveau compagnon de bureau, chimiste perdu chez les physiciens et théoricien perdu chez les expérimentateurs, Camille la plus ancienne d'entre nous toujours partante pour une bière bien méritée, Violette, partie quand je suis arrivé mais qu'on a néanmoins souvent revu, Davide who taught me a lot about coffee and cooking, Angela, I like you even though I don't understand how you can enjoy brutalist architecture; I hope this thesis won't look (too) pretentious!, and last but not least, Angus with whom I shared an office and many interesting discussions about science
and everything else. Je remercie également tous les autres doctorants, postdocs et permanents de l'INSP –
avec une mention spéciale pour mes voisins mitoyens Thierry et Laurent – qui ont fait partie de ma vie
pendant ces années, et rompu le calme du couloir 22-32.

To my colleagues and friends in St Andrews with whom I didn't spend enough time! Gerald, I had so many interesting discussions about science, politics and life with you! It's a chance that you are so opinionated. Dom Gribben and Dom Rouse, you guided me during my first year. Kristin, you are a role model for PhD candidates. Piper, Zia and Roos, it's a shame that we mostly saw one another through Teams, and I really enjoyed the few times we met in person. Chris Hooley, always cheerful and so knowledgeable about everything. I'm still impressed by your on-the-spot recitations of French poetry! Jonathan Keeling, so scientifically precise, always constructive and kind.

I would like to thank all the other PhD candidates, postdocs and permanent researchers in the School of Physics and Astronomy for the good times spent together during tea breaks, theory cakes, theory drinks and many other occasions. I also express my gratitude to the administrative staff who always did a wonderful job and made everything so easy.

I am grateful to my technical partner at dstl, Nigel Colley who followed my work thoroughly from the start.

À ma famille qui, comme de nombreuses familles de doctorants, ne comprend pas vraiment ce à quoi je passe mes journées. Papa, Maman, Martin, Anne et Guillaume, vous m'avez toujours supportés.

Three cheers for les escrimeurs de la rue d'Ulm and the fencers of St Andrews! Hip Hip Hoorroway! Hip Hip Hoorroway! Hip Hip Hoorroway!

Finally, I thank His Majesty's Government for its generosity.

Funding

This work was supported by the Defence Science and Technology Laboratory and Direction Générale de l'Armement through the Anglo-French Ph.D. scheme, and by the University of St Andrews (School of Physics and Astronomy).
Publications

Thibaut, tell me, why do French
scientists' names always start
with 'La': Laplace, Lagrange,
Laguerre, Lavoisier...

Discussion at a coffee break of
Polaron Day that made me
realise that I am, by accident,
part of a prestigious lineage.

Part of the work done during my PhD was published or submitted as articles in peer-reviewed journals

[1] Thibaut Lacroix, Angus Dunnett, Dominic Gribben, Brendon W. Lovett, and Alex Chin. *Unveiling non-Markovian spacetime signaling in open quantum systems with long-range tensor network dynamics*, Phys. Rev. A 104, 052204 – Published 4 November 2021.

[2] Thibaut Lacroix, Brendon W. Lovett, Alex W. Chin, *From Non-Markovian Dissipation to Spatiotemporal Control of Quantum Nanodevices*, arXiv:2205.11247, Submitted to Quantum.

[3] Thibaut Lacroix, *Un modèle simple de décohérence*, Under review at Bulletin de l'Union des Physiciens.

And other parts lead to an article currently in preparation

Thibaut Lacroix, Brendon W. Lovett, Alex W. Chin, *Connectivity matters: Impact of bath modes ordering and geometry in Open Quantum System simulation with Matrix Product State*, In Preparation.
Abstract

A better understanding of dissipation is crucial for understanding real-world quantum systems. Indeed, all quantum systems experience interactions with an (often) uncontrollable outside environment that can lead to a decay of excited state populations and a loss of quantum coherences. The study of dissipation is timely as the development of next-generation nanoscale quantum technologies is on its way, and the existence of non-trivial quantum effects in biological systems is being seriously investigated. However, descriptions of dissipation in quantum systems are reduced (most of the time) to time-local approaches and (everywhere) to space-local independent environments. These simplifying assumptions do render analytic and numerical calculations possible, yet they get rid of a breadth of physical processes that can alter radically the quantum systems’ dynamics. In this thesis, building on a numerically exact tensor networks method, we developed a technique able to handle spatio-temporal correlations between a quantum system and bosonic (i.e. vibrational, electromagnetic, magnons, etc.) environments. With this method we studied the signalling process - a form of information backflow - in quantum systems, and uncovered how it can induce non-trivial dynamics, and be leveraged to populate otherwise inaccessible excited states. We also evidenced the ability of ‘non-local’ environment reorganisation, induced by system-environment interactions, to radically change the nature of the thermodynamically favoured system ground state. The new phenomenology of physical processes, resulting from considering quantum systems interacting with a common environment, has important consequences for the design of nanodevices as it gives access to new control, sensing and cross-talk mechanisms. In another vein, these results might also give us a new framework to study and interpret (quantum?) effects in the biological realm.

Résumé

Une meilleure compréhension de la dissipation est cruciale pour décrire les systèmes quantiques de façon
réaliste. En effet, tous les systèmes quantiques interagissent avec un environnement extérieur (souvent)
incontrôlable qui peut conduire à une décroissance des populations d'états excités et à une perte de
cohérence quantique. L'étude de la dissipation est opportune car le développement des technologies
quantiques de nouvelle génération à l'échelle nanométrique est en cours, et l'existence d'effets quantiques
non triviaux dans les systèmes biologiques est sérieusement étudiée. Cependant, les descriptions de la
dissipation dans ces systèmes sont réduites (souvent) à des approches locales en temps et (partout) à des
environnements locaux et indépendants. Ces hypothèses simplificatrices rendent les calculs analytiques
et numériques possibles, mais elles font abstraction d'un grand nombre de processus physiques qui
peuvent modifier radicalement la dynamique des systèmes quantiques. Dans cette thèse, en se basant
sur une méthode de réseaux de tenseurs numériquement exacte, nous avons développé une technique
capable de traiter les corrélations spatio-temporelles entre un système quantique et des environnements
bosoniques (c.-à-d. vibrationnel, électromagnétique, de magnons, etc.). Avec cette méthode, nous avons
étudié le processus de signalisation - une forme de retour d'information (information backflow) - dans un
système quantique, et nous avons découvert comment il peut induire une dynamique non triviale, et être
utilisé pour peupler des états excités autrement inaccessibles. Nous avons également mis en évidence la
capacité de la réorganisation "non locale" de l'environnement, induite par les interactions avec le système,
à changer radicalement la nature de l'état fondamental thermodynamiquement favorisé du système. La
nouvelle phénoménologie de processus physiques, résultant de la description de systèmes quantiques
interagissant avec un environnement commun, a des conséquences importantes pour la conception de
nano-dispositifs, car elle donne accès à de nouveaux mécanismes de contrôle, de détection et de diaphonie.
Dans une autre direction, ces résultats pourraient également nous donner un nouveau cadre pour étudier
et interpréter des effets (quantiques ?) dans le domaine biologique.
# Introduction

> It hurt a little, scattering their
> hard-won revelation to the wind
> for anyone to use — perhaps even
> to beat them to some far greater
> prize — but they'd relied on the
> generosity of their predecessors
> from the moment they'd arrived
> on Nazdeek, and the sheer scale of
> the overall problem made it
> utterly perverse to cling selfishly
> to their own small triumph.

Greg EGAN,
Riding the Crocodile

## 1.1 Opening quantum physics

Since its discovery slightly more than a hundred years ago, quantum theory radically changed our understanding of physics at small length-scales. Nevertheless its *actual* technological outputs – even though they had a huge impact on modern society – mostly remain limited to cases that can be argued to be within the reach of semi-classical descriptions: lasers, semiconductors, Nuclear Magnetic Resonance, etc. Applications that necessitate a fundamentally quantum description (such as quantum communication, computation or sensing) exist only at the stage of proof of concepts on noisy implementations on lab benches. This state of affairs is the direct consequence of the fragile nature of quantum coherence which disappears when a quantum system is interacting with a macroscopic uncontrolled environment. This loss of coherence resulting from interactions with the outside world is the main reason why, for instance, quantum computers are so difficult to manufacture and physical qubits need to be kept near absolute zero temperature. Indeed, insulating a quantum system from its environment enables to treat it as a closed system, as described by textbook quantum physics where the system is associated with a state in a complex Hilbert space, and this state evolves in an unitary fashion through the Schrödinger equation.

However, as we are used to in classical physics, every system does

1.1 Opening quantum physics . 1

1.2 Inspirations from biology . . 3

1.3 Outline . . . . . . . . . . . . . . 6
interact with the rest of the world and in doing so dissipates energy. Whether this dissipation has to be taken into account or not is (simply) a question of time-and-length scales. These dissipative processes can 'live' at the 'margin' of the theory – like friction in mechanics – or be at its very core – like baths in thermodynamics. It is even argued that the open systems view should be *fundamentally* considered on an equal footing with the closed systems view [4]. In any cases, considering the interaction of a system of interest with its environment is *central* for the goal of understanding real life systems, even when they are as simple as the wheels of a car, or complex like a four-stroke engine.

Given that we live in a time – the advent of the so called second quantum revolution – when the wish to design (nano)devices that make use of quantum theory most fundamental features is vocal; a comprehensive description of dissipative (i.e. open) quantum system is needed. The foundations of the theory of open quantum systems have been laid down in the 60’s during the development of quantum optics [5, 6] but remained limited to the – perfectly justified in this context – simpler cases where interactions are weak and the environment ‘memoryless’. At the same time, the theoretical development of quantum master equations and their formalisation took place between the 60’s and the 80’s. The pivotal results of this period were the famous GKSL and Redfield master equations [7, 8] that allow to describe the dissipative dynamics of a quantum system without describing the microscopic degrees of freedom of the environment, nor *a fortiori* their dynamics. All these approaches rely on the so called Markovian approximation – the memoryless environment – that will be discussed in the next chapter. However, models originating from condensed matter require to go beyond weak coupling and the Markovian approximation [9]. The first investigations in this direction have been performed during the 80’s and the 90’s, especially by studying quantum Brownian motion [10] and the Spin Boson Model [11]. Concurrently, a different formalism relying on path integrals was developed by Vernon and Feynman to generally describe a quantum system interacting linearly with an environment [12]. This influence functional formalism was extremely insightful and is still nowadays the source of numerous perturbative and numerical methods (especially in quantum chemistry). In spite of major advances in the field, the applications of the theory of open quantum systems beyond the Markovian case remained limited because the computational tools needed to accurately describe the dynamics of non-analytical models were not available until the present day.

During the last decade, a broad range of numerically exact methods has been developed to tackle this challenging problem. Among these methods, a subset relies on *tensor networks*; a tool initially applied in the context of quantum many-body physics to describe the quantum state of multipartite systems. At the heart of tensor network methods is the idea – that will be developed further in a dedicated chapter – that physically relevant states inhabit a ‘small’ corner of the total Hilbert space, and that it is possible to restrict the dynamics to this corner to render complex computations tractable. Furthermore, tensor

A numerically exact method is a method where convergence to the actual dynamics is ensured by the tuning of some simulation parameters. Said differently, it is a method for which the numerical error can be estimated exactly in terms of the convergence parameters.
networks give us important insights about the entanglement geometry of quantum states, and are sometimes boldly considered to be a new language for (finite-dimensional) quantum physics [13, 14].

### 1.2 Inspirations from biology

Renewed interest about whether quantum physics is necessary to understand some biological processes has grown in the previous decade. In a trivial sense, all of biology relies on quantum physics at the biochemical level as it is the theory governing chemical bonds, and thus the properties of molecules and their possible chemical reactions. Therefore, when we are asking the question of the role of quantum effects in biological processes, we are speaking of effects beyond these trivial ones. In the last decades new results have pointed to quantum mechanics being decisive to explain key features of some biological functions [15, 16]. Photosynthesis is the canonical example of such systems. Photosynthesis is the light-harvesting process that transforms energy contained in photons into chemical energy that can be used by plants, algae and cyanobacteria to grow and perform vital functions. This process relies on a large number of so called pigment-protein complexes forming light-harvesting complexes that collect electromagnetic energy from sunlight and funnel it to a reaction centre. Schematically, the biological energy conversion consists in a photon being absorbed by the pigment complexes, this absorption creates an exciton, i.e. a bounded pair of electron and hole, that is transported thanks to dissipative processes to the reaction centre where charge separation occurs. The extracted electron will then be used for the splitting of water molecules that produces hydrogen. The obtained hydrogen is then used to form the NADPH and ATP molecules that chemically store the harvested energy, and will fuel the reactions leading to the synthesis of sugar and starch. The quantum efficiency of the photon/exciton conversion is estimated around 95% and the yield of charge separation event per formed exciton is close to 100% [15]. The main pigments involved in light harvesting are (bacterio-)chlorophylls and carotenoids, and the protein scaffolds around them tune their excited state energies and the inter-pigment interactions. Figure 1.1 shows a light-harvesting complex and a reaction centre found in Thermochromatium tepidum bacteria.

The green sulfur bacteria's photosystem has been studied extensively by biologists and spectroscopists. It possesses an additional type of pigment-protein complex in between the antennas and the reaction centre – an excitonic wire – called the Fenna-Mathews-Olson complex (FMO), depicted in Fig. 1.2. The FMO is made of seven bacteriochlorophylls pigments held together non-covalently by a protein scaffold. The FMO is often considered as the prototypical example of biological systems evidencing the existence of non-trivial quantum effects in biology. Experimental results of beatings in a cross-peak of 2D electronic spectra of the FMO at low and physiological temperature [17, 18] –

Light-harvesting complexes are often referred to as 'antennas', as they absorb electromagnetic radiation.

The other product of the water splitting reaction is the dioxygen released by photosynthetic life forms.
**Figure 1.1:** Structure of photosynthetic LH1-RC super-complex in Thermochromatium tepidum. The reaction centre is made of the protein depicted in green and orange in the foreground and the light-harvesting complex 1 (LH1) is in the background. The diameter of LH1 is of the order of 10 nm and the height of the super-complex is around 20 nm.

<figure><img src="image_9.png" /></figure>

which have been associated with the existence of quantum coherences – have opened a fascinating and animated debate on the interpretation of these signals. Similar results have been reported on spectra of light-harvesting complexes in marine algae [19].

These results open the door to the exciting question of the possibility of maintaining the existence of quantum features in a complex and highly noisy system at room temperature. Previous works have studied *Environment-Assisted Quantum Transport* to describe energy transport in light-harvesting complexes, where interactions with Markovian baths (i.e. memoryless) can enhance excitonic current, enable specific energy transitions, or make the system insensitive to structural disorder [20–25]. However, the highly structured protein environment of the photosynthetic pigments that constitute the light-harvesting complexes and reaction centres requires to go beyond the Markovian treatment [21, 23].

This possibility resonates strongly with the aim of the current quantum 'revolution' and the current engineering capabilities of quantum nanodevices. Therefore, investigating the dynamics of open quantum systems in structured environments at finite temperature is interesting and useful for at least two reasons. Firstly, it might help us gain a new understanding of possible bio-physical mechanisms, and secondly – even if it turns out that biological systems do not rely on quantum mechanical effects – adding new finite-temperature quantum mechan-
<figure><img src="image_10.png" /><figcaption>**Figure 1.2:** Fenna-Matthews-Olson complex (FMO) from Chlorobaculum Tepidum bacteria. Seven bacteriochlorophylls pigments are held inside a protein backbone.</figcaption></figure>

ical processes to our toolbox is valuable. Hence, pushing further the investigation of out-of-equilibrium open quantum dynamics, with state-of-the-art methods, is necessary for both the development of quantum technologies and the study of fundamental biological processes. Moreover, these two endeavours are not disjoint as the results (positive or negative) obtained in quantum biology might be useful to propose novel organic quantum nanodevices. One of the key technical issues to the description of these complex multi-components quantum systems is the inclusion of environmental spatial correlations, as stated in the monograph of reference in the field of quantum biology [15]

*It should be noted however, that at present the inclusion of spatial quantum correlations between environments is challenging and restricted to simple cases. This may be a drawback for certain applications in which spatial correlations play a significant role*

This problem will be addressed directly in the work presented in this thesis.

More generally, the richness of nanoscale biological processes that are taking place inside living organisms, and are relying on the concerted actions of a few fermions or excitons – over very diverse timescales

In addition, these nanoscale biological processes are cyclic processes, and can thus be compared to heat engines. This comparison opens the questions of their study in terms of (quantum) thermodynamics [26, 27].
and lengthscales – to perform energy harvesting, energy transduction, charge separation, catalysis, etc., make them a wonderful playground to explore many-body dissipative quantum physics.

## 1.3 Outline

This thesis is divided into three parts. The first part, **Background**, is made of chapter 2 and is concerned with the theory of open quantum systems, how they are related to the usual closed systems described by the Schrödinger equation and what are the consequences of the different approximations one might perform to study such systems. Notably, the concepts of Markovian and non-Markovian environments are introduced. Some material of this chapter lead to a pedagogical article currently under review [3].

The second part, **Methods**, composed of chapters 3, 4 and 5, deals with the methods I have used and developed to study the dynamics of non-Markovian quantum systems. The fundamental concepts of tensor networks and why they are a tool of choice to describe one-dimensional many-body quantum states and their dynamics are introduced in chapter 3. In chapter 4 we present the environment chain-mapping technique at the core of the versatile TEDOPA method that enables a description of zero and finite temperature open quantum systems states as simple tensor networks called Matrix Product States. We then derive an original extension of this technique to spatially extended systems [1], that enable the resolution of the problem quoted above. Finally, chapter 5 compares different ways of representing the state of an open quantum system as a tensor network and shows that the ‘intuitive’ one usually used to perform simulation is not the most efficient one. Furthermore, I show that simpler representations seem to be the most efficient ones and relate this to the connectivity of the network.

In the last part, **Applications**, I introduce applications of the framework and tools developed in the previous chapters [1, 2]. Chapter 6 – motivated by questions arising from biological systems – explores the concept of environmentally mediated spatiotemporal signalling in open systems, and chapter 7 applies it to a generic model of dissipation in multi-components quantum system that can describe a wide diversity of (bio)chemical and man-made systems. I also discuss in this chapter the biological problem of allostery and how the new framework introduced in this work could be useful to improve its understanding.

The final chapter, chapter 8, aims at summarising the main results of this thesis and discusses the next steps of this work (some of them already under investigation).
In this thesis, alongside the scientific and research materials, are small
biographical notes on the scientists explicitly cited in the main text.
I decided to add them as a reminder that science is done by actual
human beings and that our work is built on top of theirs – *Nanos
gigantum umeris insidentes*. These biographical sidenotes become rarer
as one advances in the thesis because most recent results are not named
(yet?) upon their discoverer.
<figure><img src="image_11.png" /></figure>
BACKGROUND
<figure><img src="image_12.png" /></figure>
> If we wish to remain human, then there is only one way, the way into the open society. We must go into the unknown, the uncertain and the insecure [...]

Karl POPPER, The Open Society and its Enemies

## 2.1 Isolated, closed and open systems

Physicists draw imaginary boundaries between the object they are interested in and the rest of the world in order to focus only on the features that are characteristic of the object. This distinction between the object of interest – called the physical system or *system* – and the rest of the world – called the environment of the system or *environment* – leads to a classification of physical systems depending on the nature of their interactions with their environment. This classification is relevant at every scale regardless of the nature of the system itself: it can be a collection of point-like particles as in mechanics, a macroscopic system as in thermodynamics or quantum fields as in condensed matter. An illustration in the context of quantum systems is shown in Fig. 2.1.

**Definition 2.1.1** A physical system is *isolated* if it can't exchange energy or matter with the outside world.

The prototypical example of an isolated system is the Thermos bottle and the trivial one is the entire universe.

**Definition 2.1.2** A physical system is *closed* if it can only exchange energy with the outside world.

An everyday example of a closed system would be a hermetically sealed room where no matter could get out but heat could flow through the windows, walls, etc.

**Definition 2.1.3** A physical system is *open* if it can exchange energy and matter with the outside world.

2.1 Isolated, closed and open systems 11

2.2 Environment, bath and reservoir 13

2.3 Time evolution 13

2.4 Harmonic Bath, correlation function and spectral density 15

2.5 (Non-)Markovianity 19

<figure><img src="image_13.png" /><figcaption>**Figure 2.1:** Schematic representation of the separation of a physical object *S* + *E* into a *System S* whose free dynamics is induced by the Hamiltonian $\hat{H}_S$ and an *Environment E* with a Hamiltonian $\hat{H}_E$. The black double arrows represent the possible interactions between these two parts through the Hamiltonian $\hat{H}_{int}$.</figcaption></figure>
The reader, as well as every other living creature, is an open system who needs to exchange energy (e.g heat) and matter (e.g food) with the outside world.

Quantum theory, as developed in the first half of the twentieth cen-
tury, is formulated for a closed system which can be described by a
wave-vector (if it is a pure state) or a density matrix (if it is a statistical
mixture). However, realistic quantum objects are never isolated sys-
tems but are interacting with 'external' (e.g. electromagnetic, electronic,
vibrational...) degrees of freedom which constitute their environments.
These interactions have consequences for the dynamics of the system
as they often lead to the decay of excited states populations toward the
ground state of the system (i.e. dissipation) and in the loss of quantum
coherences in the eigenbasis (i.e. decoherence). The study of the result-
ing dynamics is the subject of the theory of open quantum systems
(TOQS).

Understanding the processes leading to dissipation and decoherence
is crucial both from a theoretical and an experimental point of view.
It allows us to understand the observed dynamics of quantum sys-
tems (physical, chemical or biological ones) and to invent strategies
to suppress or harness dissipation and decoherence. This knowledge
becomes of paramount importance as we are entering a so called sec-
ond quantum revolution [16, 28] in which quantum superposition and
entanglement are being used as resources to reach a desired function
(e.g. communication, computation, energy transport...). As the scale
of these devices reduces for integration purposes, they become more
exposed to the very dissipative phenomena described by TOQS. En-
vironments are often complex and composed of a large number of
degrees of freedom. Hence, the dynamics of the system is not an easy
thing to predict.

Another motivation for the study of OQS comes, surprisingly, from biology. Spectroscopic experiments performed on photosynthetic complexes show the presence of quantum coherences during energy transport in the form of beatings [17, 19]. A lot of debate emerged following these measurements on the nature of the observed signal [29], the functional role of quantum coherences [25] and the (absence of) selection pressure associated with the existence of coherence [30]. These results and debates around photosynthetic complexes rejuvenated interest in other biological processes where quantum physics might play a non-trivial role [20, 31]. Among the other processes where quantum effects are *hypothesised* we can cite avian navigation (which might be the most promising one), olfaction, the physiological mechanism behind anaesthetics, or proton tunnelling in enzymatic reaction and DNA mutations [20, 31–33]. All these biological processes have in common the strong interplay between a quantum system (excitons, spins...) with a few degrees of freedom and a very structured environment at room temperature (solvent and protein complexes). This is a striking difference from the man-made quantum devices that need to be shielded from their environment.
## 2.2 Environment, bath and reservoir

The environment of the system of interest could be anything: vibrational modes of molecules, a gas of conduction electrons in a metal, electronic spins, etc. Nevertheless, big picture characterizations that don't take into account the actual nature of the degrees of freedom can be made. Two types of environments are particularly relevant in the study of open quantum systems: reservoir and heat bath (commonly abridged bath). These definitions are in accordance with the use of the same terms in classical thermodynamics. When the environment has a larger number of degrees of freedom than the reduced system, it is often possible to describe it with infinitely many degrees of freedom. Sometimes the environment really has an infinite number of degrees of freedom. For example, this is the case of the electromagnetic field.

**Definition 2.2.1** A reservoir is an environment with an infinite number of degrees of freedom.

The continuous nature of the reservoir's spectrum implies that the Poincaré recurrence time of the system's state is in principle infinite; meaning that the dynamics of a system coupled to a reservoir is irreversible. If the reservoir has a well-defined temperature, it is called a bath.

**Definition 2.2.2** A heat bath (or bath) is a reservoir in a thermal equilibrium state.

<figure><img src="image_14.png" /><figcaption>**Figure 2.2:** John von Neumann (December 28, 1903 – February 8, 1957) was a Hungarian-American mathematician, physicist, computer scientist and economist.</figcaption></figure>

More generally, a bath is a reservoir with a well defined and constant conjugate thermodynamic variable: temperature for a heat bath (or thermostat), pressure for a pressurestat, chemical potential for a matter reservoir...

## 2.3 Time evolution

The quantum theory of closed systems describes the evolution between two times $t > t_0$ of a quantum system associated with a density operator $\hat{\rho}(t)$ via the unitary time-evolution operator $\hat{U}(t, t_0)$ defined with the Hamiltonian operator $\hat{H}^1$ [34]

$$ \hat{\rho}(t) = \hat{U}(t, t_0) \hat{\rho}(t_0) \hat{U}^{\dagger}(t, t_0). \qquad (2.1) $$

The equation of motion (EOM) of the density matrix $\hat{\rho}$ is given by the von Neumann equation

$$ \frac{d}{dt} \hat{\rho}(t) = -\frac{i}{\hbar} [\hat{H}, \hat{\rho}(t)] \stackrel{\text{def.}}{=} \mathcal{L} \hat{\rho}(t). \qquad (2.2) $$

In Eq. (2.2) the Liouvillian super-operator $\mathcal{L}$ has been defined. A super-operator $\mathcal{S}$ is an application that acts on the space of operators to the space of operators. With the Liouvillian, we can formally solve the equation of motion

$$ \hat{\rho}(t) = \overleftarrow{T} \exp \left( \int_{t_0}^{t} dt' \mathcal{L}(t') \right) \hat{\rho}(t_0) \qquad (2.3) $$

¹ In the general case

$$ \hat{U}(t, t_0) = \overleftarrow{T} \exp \left( -\frac{i}{\hbar} \int_{t_0}^{t} dt' \hat{H}(t') \right) $$

where $\overleftarrow{T}$ is time-ordering operators chronologically from right to left, and $\hat{H}(t)$ depends parametrically on time $t$.

All the operators introduced so far are in the so called Schrödinger picture.
<figure><img src="image_15.png" /><figcaption>**Figure 2.3:** Joseph Liouville (March 24, 1809 – September 8, 1882) was a French mathematician who worked on number theory, complex analysis, differential geometry and topology, but also mathematical physics and astronomy.</figcaption></figure>

where, here, $\hat{T}$ is time-ordering super-operators chronologically from right to left. The goal of the TOQS is to be able to describe the time-evolution of the density matrix of the system of interest given its interactions with its environment [9, 35, 36].

The main approach is to consider the system and the environment as a larger isolated system, where the evolution of the joint density matrix of the {System + Environment} is described by the von Neumann equation (2.2), and then to get rid of the degrees of freedom of the environment. Several approximations are then possible to obtain a tractable formulation of the evolution of the system.

The total Hilbert space $\mathcal{H}$ of the isolated {System + Environment} is identified as a tensor product $\mathcal{H} = \mathcal{H}_S \otimes \mathcal{H}_E$ where $\mathcal{H}_S$ is associated with the system, and $\mathcal{H}_E$ associated with the environment.

We can thus decompose the Hamiltonian into three different terms:

* the system Hamiltonian $\hat{H}_S$ which contains only terms acting on $\mathcal{H}_S$ – it describes the free evolution of the reduced system.

* the environment Hamiltonian $\hat{H}_E$ which contains only terms acting on $\mathcal{H}_E$ – it describes the free evolution of the environment.

* the interaction Hamiltonian $\hat{H}_{\text{int}}$ which acts on the full Hilbert space $\mathcal{H}$.

The density matrix of the system, often referred to as the reduced density matrix, can be defined by taking the partial trace over the environment degrees of freedom of the full density matrix

$$ \hat{\rho}_S(t) = \text{tr}_E[\hat{\rho}(t)]. \qquad (2.4) $$

Expectation values of observables acting on the system's Hilbert space are defined as follows

$$ \langle \hat{A} \rangle = \text{tr}_S [\hat{A} \hat{\rho}_S(t)] = \text{tr} [\hat{A} \otimes \hat{1}_E \hat{\rho}(t)], \qquad (2.5) $$

where $\hat{1}_E$ is the identity on $\mathcal{H}_E$.

Therefore taking the trace of Eq. (2.2) over the environment formally gives the equation of motion of the system density matrix

$$ \frac{d}{dt} \hat{\rho}_S(t) = -\frac{i}{\hbar} \text{tr}_E [\hat{H}, \hat{\rho}(t)]. \qquad (2.6) $$

The equation of motion Eq. (2.6) is called the master equation.

**Definition 2.3.1** A master equation (ME) is an equation of evolution of the density matrix $\hat{\rho}_S(t)$ of the reduced system S. It has the form

$$ \frac{d}{dt} \hat{\rho}_S(t) = L \hat{\rho}_S \qquad (2.7) $$

where L is a linear operation on $\hat{\rho}_S$ that can be non-local in time.
The formal solution of Eq. (2.7) gives the so called dynamical map $V(t, t_0)$ giving the evolution of $\hat{\rho}_S$ between $t_0$ and $t > t_0$

$$ \hat{\rho}_S(t) = V(t, t_0) \hat{\rho}_S(t_0). \qquad (2.8) $$

We can define as follow the so called interaction picture in which the interaction Hamiltonian $\hat{H}_{\text{int}}$ gives the evolution of the density matrix and the free Hamiltonian $\hat{H}_0 = \hat{H}_S + \hat{H}_E$ defines a unitary transformation for the other operators

$$ \hat{\rho}^I(t) = \hat{U}_I(t, t_0) \hat{\rho}(t_0) \hat{U}_I^\dagger(t, t_0) \qquad (2.9) $$

$$ \text{with } \hat{U}_I(t, t_0) = \overleftarrow{T} \exp \left( -\frac{i}{\hbar} \int_{t_0}^{t} dt' \hat{H}_{\text{int}}(t') \right), \qquad (2.10) $$

$$ \hat{O}^I(t) = \hat{U}_0^\dagger(t, t_0) \hat{O}(t) \hat{U}_0(t, t_0) \qquad (2.11) $$

$$ \text{with } \hat{U}_0(t, t_0) = \overleftarrow{T} \exp \left( -\frac{i}{\hbar} \int_{t_0}^{t} dt' \hat{H}_0(t') \right). \qquad (2.12) $$

The main point of this new picture is that the evolution of the density matrix only depends on the interaction Hamiltonian. In the interaction picture the von Neumann equation and its formal solution become

$$ \frac{d}{dt} \hat{\rho}^I(t) = -\frac{i}{\hbar} [\hat{H}_{\text{int}}^I(t), \hat{\rho}^I(t)] \stackrel{\text{def.}}{=} \mathcal{L}^I(t) \hat{\rho}^I(t), \qquad (2.13) $$

$$ \hat{\rho}^I(t) = \overleftarrow{T} \exp \left( \int_{t_0}^{t} dt' \mathcal{L}^I(t') \right) \hat{\rho}(t_0). \qquad (2.14) $$

From there, the theory of open quantum systems could be (crudely) summarised as a variety of cases and the corresponding approximations that can be made to write down an explicit form of the dynamical map $V(t, t_0)$.

We will now focus on a specific type of bath that will be of particular interest in the rest of this manuscript: the harmonic bath.

For the rest of the manuscript, the I superscript on the interaction picture density matrix will be dropped.

## 2.4 Harmonic Bath, correlation function and spectral density

Equation (2.14) gives the formal solution of the von Neumann equation for the total density matrix $\hat{\rho}(t)$. Tracing over the environment, we obtain a formal solution for the reduced density matrix $\hat{\rho}_S(t)$ which can then be tractable with further simplification coming from the specific physical problem under consideration. If there are no initial correlations between the system of interest and its environment, the initial state is a product states $\hat{\rho}(t_0) = \hat{\rho}_S(t_0) \otimes \hat{\rho}_E(t_0)$. In that case, the trace in Eq. (2.14) becomes an expectation value on the environment

$$ \hat{\rho}_S(t) = \overleftarrow{T}_S \left\langle \exp \left( \int_{t_0}^{t} dt' \mathcal{L}^I(t') \right) \right\rangle_E \hat{\rho}_S(t_0). \qquad (2.15) $$

Which also implies that we can separate the time-ordering $\overleftarrow{T} = \overleftarrow{T}_S \overleftarrow{T}_E$.
where the super-operator valued expectation is defined as

$$ \langle S \rangle_E = \text{tr}_E \left[ \overleftarrow{T}_E S \hat{\rho}_E(t_0) \right] . \qquad (2.16) $$

In Eq. (2.15), we can identify $V(t, t_0) = \overleftarrow{T}_S \langle \exp \left( \int_{t_0}^t dt' \mathcal{L}^I(t') \right) \rangle_E$ as being a time-propagator between times $t_0$ and $t$ for the system density matrix. We can pause for a few lines to illustrate an important property of dynamical maps. Contrary to isolated systems' evolution operator, open systems' time-propagator are not in general divisible

$$ \text{For } t_0 < u < t, V(t, t_0) \neq V(t, u) \cdot V(u, t_0) \text{ in general.} \qquad (2.17) $$

Indeed, as can be noticed in Eq. (2.15), in general

$$ \langle \exp \left( \int_{t_0}^{t} dt' \mathcal{L}^{I}(t') \right) \rangle_E \neq \langle \exp \left( \int_{u}^{t} dt' \mathcal{L}^{I}(t') \right) \rangle_E \langle \exp \left( \int_{t_0}^{u} dt' \mathcal{L}^{I}(t') \right) \rangle_E, \qquad (2.18) $$

because the expectation value of a product is different from the product of expectation values.

We are now going to specify an explicit Liouvillian $\mathcal{L}^I(t)$ by choosing a specific Hamiltonian that describes the environment as a bosonic bath with a linear coupling to the system

$$ \hat{H}_B = \sum_i \hbar \omega_i \hat{a}_i^\dagger \hat{a}_i, \qquad (2.19) $$

$$ \hat{H}_{\text{int}} = \hat{O} \sum_{i} (g_{i} \hat{a}_{i} + g_{i}^{*} \hat{a}_{i}^{\dagger}) \qquad (2.20) $$

where $\hat{O}$ is an operator acting on the reduced system Hilbert space ($\hat{O} = \hat{O} \otimes \hat{1}_B$), $\omega_i$ is the angular frequency associated with the $i^{th}$ bath mode and $g_i$ the coupling strength between the system and the $i^{th}$ bath mode. Furthermore, if we consider that the bath initial state is a Gibbs state$^2$, we can use the following relation [37]

$$ \langle \exp(\mathcal{X}) \rangle_B = \exp \left( \langle \mathcal{X} \rangle_B + \frac{1}{2} \langle \mathcal{X}^2 \rangle_B \right), \qquad (2.21) $$

where, here, $\mathcal{X}$ is a super-operator.

Hence, choosing $\langle \mathcal{X} \rangle_B = 0$, we can cast Eq. (2.15) in the following form$^3$

## System Density Matrix with Time Propagator (Bosonic Bath)

$$ \hat{\rho}_S(t) = \overleftarrow{T}_S \exp \left( \int_{t_0}^{t} dt' dt'' \langle \mathcal{L}^I(t') \mathcal{L}^I(t'') \rangle_B \right) \hat{\rho}_S(t_0) \qquad (2.22) $$

where $\mathcal{L}^I(t) \bullet = -\frac{i}{\hbar} [\hat{H}_{\text{int}}(t), \bullet]$. In general, the propagator is not time-local because of the time dependence of $\langle \mathcal{L}^I(t') \mathcal{L}^I(t'') \rangle_B$.

The next step is then to evaluate $\langle \mathcal{L}^I(t') \mathcal{L}^I(t'') \rangle_B$ in order to be able

When $\hat{H}_B$ describes vibrational modes and $\hat{H}_{\text{int}}$ their coupling to an electronic system, the model they define is called Linear Vibrronic Coupling.

$$ \hat{\rho}(t) = \frac{1}{Z} \exp \left( -\frac{\hat{H}}{k_B T} \right) $$

with $Z = \text{tr} [\exp(-\frac{\hat{H}}{k_B T})]$ the partition function.

<footnote>$^2$ A Gibbs state is a state of thermal equilibrium where</footnote>

<footnote>$^3$ Because $\hat{H}_{\text{int}}^I(t)$ is linear in bath operators.</footnote>
to have an explicit form for the propagator acting on $\hat{\rho}_S(t_0)$.

## Super-operator notation

Let us introduce a useful set of notation defining super-operators. Given an operator $\hat{O}$ we can define the four following super-operators

$$O^L \hat{\rho} \stackrel{\text{def.}}{=} \hat{O} \hat{\rho}, \qquad (2.23)$$

$$O^R \hat{\rho} \stackrel{\text{def.}}{=} \hat{\rho} \hat{O}, \qquad (2.24)$$

$$O^+ \hat{\rho} \stackrel{\text{def.}}{=} \{\hat{O}, \hat{\rho}\}, \qquad (2.25)$$

$$O^- \hat{\rho} \stackrel{\text{def.}}{=} [\hat{O}, \hat{\rho}]. \qquad (2.26)$$

Moreover, one can notice that $O^\pm = O^L \pm O^R$.

Therefore, we can rewrite the interaction Liouvillian as

$$\mathcal{L}^I(t) = -\frac{i}{\hbar} H_{\text{int}}^{I-}(t). \qquad (2.27)$$

The interaction Hamiltonian for a linear coupling with a bosonic bath (Eq. (2.20)) is a product of a system operator and a bath operator which can be written $\hat{H}_{\text{int}}^I(t) = \hat{O}(t) \otimes \hat{B}(t)$.

Thus, with our new super-operator notation

$$\mathcal{L}^I(t) = -\frac{i}{\hbar} (H_{\text{int}}^{I L}(t) - H_{\text{int}}^{I R}(t)) \qquad (2.28)$$

$$= -\frac{i}{\hbar} (O^L(t) B^L(t) - O^R(t) B^R(t)) \qquad (2.29)$$

$$= -\frac{i}{4\hbar} \left( (O^{-}(t) + O^{+}(t))(B^{-}(t) + B^{+}(t)) - (O^{-}(t) - O^{+}(t))(B^{-}(t) - B^{+}(t)) \right) \qquad (2.30)$$

$$\mathcal{L}^I(t) = -\frac{i}{2\hbar} (O^+(t)B^-(t) + O^-(t)B^+(t)). \qquad (2.31)$$

This leads us to

$$\langle \mathcal{L}^I(t') \mathcal{L}^I(t'') \rangle_B = -\frac{1}{4\hbar^2} \left\langle \begin{aligned} & (O^+(t') B^-(t') + O^-(t') B^+(t')) \\ & \times (O^+(t'') B^-(t'') + O^-(t'') B^+(t'')) \end{aligned} \right\rangle_B . \qquad (2.32)$$

One can show easily that for a Gibbs state and an interaction Hamiltonian given by Eq. (2.20), $\text{tr}[B^{-}B^{-}\hat{\rho}_{B}] = \text{tr}[B^{-}B^{+}\hat{\rho}_{B}] = 0$. Hence

Because $B^{\pm}\hat{\rho}_{B} \propto g_{i}\hat{a}_{i} + g_{i}^{*}\hat{a}_{i}^{\dagger}$.

$$\langle \mathcal{L}^I(t') \mathcal{L}^I(t'') \rangle_B = -\frac{1}{4\hbar^2} O^{-}(t') (O^{+}(t'') \langle B^{+}(t') B^{-}(t'') \rangle_B + O^{-}(t'') \langle B^{+}(t') B^{+}(t'') \rangle_B). \qquad (2.33)$$

Defining the two-time bath correlation function $C(t', t'')$

$$C(t', t'') = \text{tr}_B [\hat{B}(t') \hat{B}(t'') \hat{\rho}_B(t_0)], \qquad (2.34)$$
we have $\text{tr}_B [B^+(t')B^-(t'')\hat{\rho}_B(t_0)] = 4\text{Im}[C(t',t'')] \text{ and } \text{tr}_B [B^+(t')B^+(t'')\hat{\rho}_B(t_0)] = 4\text{Re}[C(t',t'')].$

Finally, we can express $\langle \mathcal{L}^I(t')\mathcal{L}^I(t'') \rangle_B$ in terms of super-operators acting only on the system and the bath two-times correlation function

### Super-Operator Valued Expectation of $\mathcal{L}^I(t')\mathcal{L}^I(t'')$

$$ \langle \mathcal{L}^I(t')\mathcal{L}^I(t'') \rangle_B = -\frac{1}{\hbar^2} O^-(t') (O^-(t'')\text{Re}[C(t',t'')] + iO^+(t'')\text{Im}[C(t',t'')]) \quad (2.35) $$

The two-time correlation function $C(t,t')$ can be made explicit using the expression of $\hat{H}_{\text{int}}^I(t)$ given in Eq. (2.20)

$$ \begin{aligned} C(t, t') &= \text{tr}_B [\hat{B}(t) \hat{B}(t') \hat{\rho}_B(t_0)] \\ &= \sum_{i,j} \delta_{i,j} (g_i^* g_j e^{i\omega_i(t-t')} n_\beta(\omega_i) + g_i g_j^* e^{-i\omega_i(t-t')} (n_\beta(\omega_i) + 1)) \\ &= \sum_i |g_i|^2 (e^{i\omega_i(t-t')} n_\beta(\omega_i) + e^{-i\omega_i(t-t')} (n_\beta(\omega_i) + 1)) \\ C(t-t') &= \sum_i |g_i|^2 \left[ \coth\left(\frac{\beta\hbar\omega_i}{2}\right) \cos(\omega_i(t-t')) - i \sin(\omega_i(t-t')) \right]. \end{aligned} \quad (2.36) $$

If we define a new quantity called the spectral density of the bath, Eq. (2.36) can be cast in an integral form.

**Definition 2.4.1** The bath spectral density contains all the information about the interaction between the system and the bath

$$ J(\omega) = \sum_i |g_i|^2 \delta(\omega - \omega_i) \text{ where } \omega \ge 0. \quad (2.37) $$

For a bosonic bath, the bath Power Spectrum $S(\omega)$ is proportional to the spectral density

$$ S(\omega) = 2\pi n_B(\omega) J(\omega) \text{ where } \omega \ge 0. \quad (2.38) $$

Finally, the general form of the correlation function of a bosonic bath is

### Bosonic Bath Auto-Correlation Function

$$ C(\tau) = \int_0^\infty d\omega J(\omega) \left[ \coth\left(\frac{\beta\hbar\omega}{2}\right) \cos(\omega\tau) - i \sin(\omega\tau) \right] \quad (2.39) $$

Hence, for a bosonic bath with a linear coupling, the time propagator of the reduced system density matrix is defined by the super-operator valued expectation $\langle \mathcal{L}^I(t')\mathcal{L}^I(t'') \rangle_B$, which is entirely defined by the

with $n_\beta(\omega_i) = \frac{1}{\exp(\beta\hbar\omega_i) - 1}$, where $\beta = (k_B T)^{-1}$, the Bose-Einstein distribution.

We can notice that, due to the stationarity of the initial bath state, the bath correlation function became a function $C(t-t')$ of a single time argument $\tau = t-t'$.

To approximate the limit case of a perfectly flat spectral density at $T = 0$ K, we look at the bath correlation function for a rectangular spectral density of height $J_0$ centered on $\omega_0$ with a width $\Delta\omega$

$$ C(\tau) = J_0 \Delta\omega \text{sinc}\left(\frac{\Delta\omega\tau}{2}\right) e^{-i\omega_0\tau}. $$

The characteristic correlation time of the bath correlation function (corresponding to the first zero of the sinc) is $\tau_B = \frac{2\pi}{\Delta\omega}$. Taking the infinite bandwidth limit we obtain an instantaneous correlation $\lim_{\Delta\omega \to \infty} C(\tau) \propto \delta(\tau)$.

If the spectral density is a Lorentzian of the same height, bandwidth and peaked at the same frequency, the correlation function becomes

$$ C(\tau) = \pi \frac{\Delta\omega}{2} J_0 e^{-\frac{\Delta\omega\tau}{2}} e^{-i\omega_0\tau} $$

which has the similar characteristic time $\tau_B$ and the same infinite bandwidth limit.
bath auto-correlation function $C(\tau)$ or equivalently by the bath spectral density $J(\omega)$. Here we recover the well-known result that for a Gaussian environment, the dynamics of the system does not depend on the microscopic details of the environment but only on the bath correlation function [9, 35].

We can already see that for different spectral densities $J(\omega)$ the bath correlation function is going to have different behaviours, and especially it will decay on different timescales $\tau_B$. We call this timescale $\tau_B$ the bath correlation time.

For example, if only one bath mode couples to the system the correlation function has a purely oscillatory behaviour and thus an infinite correlation time.

Generally, if the system couples to a range of frequencies $\Delta\omega$ the correlation time is going to be of the order $\tau_B \approx 2\pi/\Delta\omega$.

## 2.5 (Non-)Markovianity

If all the bath modes couple to the system the same way, i.e. $J(\omega) = J_0$, at zero temperature the correlation function is a Dirac delta function. In that case, there are only instantaneous correlations into the bath. Such a 'white noise' spectral density is a limit case of what is called a Markovian behaviour.

**Definition 2.5.1** An environment is called Markovian if its timescale of relaxation $\tau_B$ is much shorter than the timescales of the processes happening in the system $\tau_S$

$$ \tau_B \ll \tau_S. \qquad (2.40) $$

Being in a Markovian regime allows one to make important simplifications to the equations of motion as the dynamic of the system doesn't depend on past interactions with the bath. With the additional assumptions of separability of the system and the bath states $\hat{\rho}(t) \approx \hat{\rho}_S(t) \otimes \hat{\rho}_B(t)$ because of weak coupling (Born approximation) and the secular approximation, one can write down the famous Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) equation.

**Definition 2.5.2** Assuming the Born, Markov and secular approximations, the most general form of a ME is the GKSL master equation

$$ \frac{d}{dt} \hat{\rho}_S(t) = -\frac{i}{\hbar} [\hat{H}_{\text{eff}}, \hat{\rho}_S(t)] + \sum_{k=1}^{d^2-1} \gamma_k \left( \hat{A}_k \hat{\rho}_S(t) \hat{A}_k^\dagger - \frac{1}{2} \{ \hat{A}_k^\dagger \hat{A}_k, \hat{\rho}_S(t) \} \right), \qquad (2.41) $$

where $\hat{H}_{\text{eff}} = \hat{H} + \hat{H}_{\text{LS}}$ is an effective Hamiltonian with $\hat{H}_{\text{LS}}$ a Lamb shift term renormalising the system's energy levels due to the interaction with the bath, $\gamma_k > 0$ are relaxation rates, $\hat{A}_k \in \mathcal{B}(\mathcal{H}_S)$ are traceless so called jump operators and $d = \text{dim}(\mathcal{H}_S)$.

<figure><img src="image_16.png" /><figcaption>**Figure 2.4:** Andrei Markov (June 14, 1856 – July 20, 1922) was a Russian mathematician best known for his work on stochastic processes. A primary subject of his research later became known as Markov chains and Markov processes.</figcaption></figure>

The secular approximation considers that the infinitesimal timespan $dt$ is larger than the fastest transition in the system.

Vittorio Gorini (born 4 December 1940) is an Italian physicist and professor at the Università degli Studi dell'Insubria. His work include research in the structure of quantum dynamical semigroups and formulation of de Sitter relativity.
⁴ The adjective ‘non-Markovian’ can also be used in a general sense of ‘not relying on a Markovian approximation’ without implying that this is of particular importance for the dynamics of the system.

<figure><img src="image_17.png" /></figure>

**Figure 2.5:** Andrzej Kossakowski (20 February 1938 – 31 January 2021) was a Polish theoretical physicist and a professor at the Nicolaus Copernicus University. He was best known for his work on open quantum systems.

**Figure 2.6:** When the environment is Markovian, energy injected by the system is simply dissipated away, whereas in the non-Markovian case it can be injected back into the system at a latter time and a different place.

Equation (2.41) is a master equation describing the evolution of $\hat{\rho}_S$ with a unitary contribution given by the effective Hamiltonian and a dissipative part generated by the action of the jump operators on the instantaneous density matrix. Hence, by construction, Eq. (2.41) is a time-local equation. The GKSL equation of motion is often used in the context of quantum optics [6].

We call a bath *non-Markovian* if we cannot neglect its past interactions with a system to describe the dynamics of the system⁴. Whereas the Markovian case is well-studied, the non-Markovian one is still poorly understood as it is more complex. A consequence of non-Markovianity is that energy that has been dissipated into the environment can be re-injected into the system at later times (as illustrated in Fig. 2.6) and thus possibly altering drastically its dynamics. This corresponds to the so called *information backflow* from the environment to the system which is often responsible for recurrences in the system’s state [38].

<figure><img src="image_18.png" /></figure>

Note that the backflow can happen at a different position in space than the one where the energy was dissipated initially. This is of strong importance for multi-component systems like machines designed by humans or emerging from natural selection which often rely on different pieces being arranged in space to realise a desired function. For instance, an engine is composed of valves, pistons, etc. that act in an ordered manner, such as a four-stroke cycle, to perform energy conversion; and a light harvesting complex is made of pigments held together by a protein scaffold in such a way to efficiently transfer energy. Hence, to describe accurately such concerted actions in a extended device under dissipation, it is needed to go beyond the usual assumption of Markovian bath dynamics.

This opens the possibility that injecting an excitation in an extended system causes the local structure to relax to a new equilibrium position, but in doing so key system properties such as energy gaps or couplings to other systems can be strongly modified in the new conformation [39, 40]. It thus becomes possible for local reorganization dynamics to *propagate* and effect dramatic changes at distal locations at later times [41].
This phenomenology will be of central interest in this thesis as we will show in subsequent chapters examples of the richness of non-Markovian dynamics. However, the necessity to take into account the interaction history in order to accurately describe non-Markovian dynamics is the very source of the difficulty of this description. Major advances on this problem happened recently as modern computation techniques became available and as new numerical approaches were uncovered. For instances, the reaction-coordinate mapping incorporate the ‘non-Markovian part’ of the environment in the definition of the system and the remaining modes are described with a Markovian EOM [42]; at zero temperature, the pseudo-modes method is able to access the system dynamics in the ultra-strong coupling regime by replacing the environment by a few unphysical pseudo-modes [43]; and the reduced hierarchical equations of motion (HEOM) method maps non-Markovian effects into new terms in the density matrix EOM that, as suggested by the name, belongs to a hierarchy of equations [44]. However, these methods are either not generally applicable or numerically expensive. New techniques relying on so called tensor networks have emerged to represent either the system propagator [45, 46], the bath influence functional in the form of a process tensor [47] or the total wave function of the system and the environment together [48]. They allow a general, numerically exact and efficient description of the system and the bath dynamics. Chapters 3 and 4 treat in more detail the basics of tensor networks and a specific method, the Time-Dependent Density operator with Orthonormal Polynomials (TEDOPA), that will be used in the rest of the thesis.

<figure><img src="image_19.png" /><figcaption>**Figure 2.7:** Ennackal Chandy George Sudarshan (16 September 1931 – 13 May 2018) was an Indian theoretical physicist and a professor at the University of Texas. Sudarshan has been credited with, inter alia, Glauber–Sudarshan P representation, V-A theory, tachyons and quantum Zeno effect.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 2.8:** Göran Lindblad (9 July 1940 - 1 December 2022) was a Swedish theoretical physicist and a professor emeritus at the KTH Royal Institute of Technology, Stockholm.</figcaption></figure>
<figure><img src="image_21.png" /></figure>
METHODS
<figure><img src="image_22.png" /></figure>
The identity matrix $1$ is usually represented as a simple wire $\delta_{i,j} = i \text{---} j$.

<figure><img src="image_23.png" /><figcaption>**Figure 3.3:** Sir Roger Penrose (born 8 August 1931) is an English mathematician and mathematical physicist who was awarded the Nobel Prize in Physics "for the discovery that black hole formation is a robust prediction of the general theory of relativity".</figcaption></figure>

The action of a matrix on a vector is represented by joining the legs that correspond to the summed over index during contraction. For instance the equation $\vec{v} = A\vec{u}$, where $A$ is a $n \times m$-matrix and $\vec{u}$ is a $m$-dimensional vector and $\vec{v}$ is a $n$-dimensional vector, can be written

$$ \text{---} n \text{---} \mathbf{v} \text{---} = \text{---} n \text{---} \mathbf{A} \text{---} m \text{---} \mathbf{u} \text{---}, \quad (3.2) $$

where, to render the different dimensions of the matrix and vectors explicit, we have written them in extenso inside of the legs themselves. A matrix multiplication $C = AB$ can be written in a similar fashion

$$ \text{---} n \text{---} \mathbf{C} \text{---} p \text{---} = \text{---} n \text{---} \mathbf{A} \text{---} m \text{---} \mathbf{B} \text{---} p \text{---}. \quad (3.3) $$

More complex linear operations can also be represented, such as the trace of a $n \times n$-matrix $M$ which is depicted by contracting the two legs of the matrix together

$$ \text{tr}[M] = \text{---} \stackrel{\text{---} n}{\text{---}} \stackrel{\text{---} M}{\text{---}}. \quad (3.4) $$

The diagrammatic notation of tensors is a useful tool to simply represent complex operations involving several contractions at the same time. It gives a representation of equations and mathematical operations that is easy to read and intuitive to manipulate. Such notation was first introduced by Roger Penrose [49] and gained traction in the last twenty years when one of the leading methods to determine the ground state of 1d-system, the Density Matrix Renormalization Group (DMRG) [50, 51], was reformulated in terms of tensor networks [52], and when tensor-network-based methods applicable to other geometries [53–56] to find ground states or time-evolve quantum states were found [57–60]. Simultaneously, a reformulation of (finite dimensional) quantum theory in terms of abstract diagrams representing states, effects and processes was also developed stemming from category theory [14].

## 3.2 Matrix Product States

The tensor $c$ of a quantum state $|\psi\rangle$ corresponding to a one-dimensional system can be decomposed into a product of $N$ smaller rank-3 tensors $T_k$ (except for the first and last sites where the tensors will have a rank-2)

$$ c_{i_1 \dots i_N} = \sum_{\{\alpha\}} T_{i_1}^{\alpha_1} T_{i_2}^{\alpha_2} T_{i_3}^{\alpha_3} \dots T_{i_N}^{\alpha_{N-1}}. \quad (3.5) $$

In this form, the local tensor $T_k$ contains the information on the quantum state on site $k$ and its relation (especially the entanglement) with the neighbouring sites.
9 1 tebibit = 2⁴⁰ bits ≈ 1.0995 Tbit. In the Julia programming language or in Python, unsigned integers are coded on 8 bits, hence storing the state of a 40-site long spin chain requires around 8 Tbit.

The SVD can be interpreted geometrically. The unitary matrices $U$ and $V^\dagger$ correspond to a rotation in the complex plane and the singular matrix $S$ corresponds to a "squeezing" of the basis vectors.

There exists another version of the SVD called 'compact' where $S$ is a $p \times p$-matrix with $p = \min\{n, m\}$, and $U$ and $V^\dagger$ are respectively $n \times p$ and $p \times m$ isometries (whose definition will be given latter on). The compact version is often used in practice for numerical calculations.

reduces to a point, thus the entropy will be independent of the size of the system $S_{VN} \sim$ constant. The MPSs are states that satisfy this area law, as we will show in the next section. Hence, MPSs are the relevant states to describe 1d physical systems that are interacting through local interactions. Figure 3.4 shows a drawing of the hierarchy of manifolds inside a 1d many-body Hilbert space. Moreover, the main asset of this Ansatz – that follows from the reduction to a corner of the Hilbert space – is its practical numerical advantage over generic formulations of many-body wave functions. In full generality, for a tensor $c$ representing a $N$-component many-body state with local Hilbert spaces of dimension $d$, the number of elements that need to be stored is $d^N$. The number of coefficients thus grows exponentially with the size of the system, rapidly preventing any reasonable size calculations. Storing the state of one of the simplest many-body systems, a spin chain, of 40 sites already requires 8 tebibits.⁹ It is therefore unimaginable to practically perform any simulation this way even on such a simple system. Whereas the number of coefficients describing a MPS of the order of ≈ $NdD^2$, which grows linearly with the size of the system.

## 3.3 Singular Value Decomposition

The Singular Value Decomposition (SVD) of a tensor is a factorisation that generalises the eigen-decomposition of square matrices to rectangular matrices and higher rank tensors. As in an eigen-decomposition the goal is to represent a tensor in a diagonal form (along a given partition of its indices), denoted $S$, whose entries are called the Singular Values of the tensor.

**Definition 3.3.1** *The Singular Value Decomposition of a complex $n \times m$ rectangular matrix $T$ is a decomposition such that*

$$T = USV^\dagger \quad (3.8)$$

where $U$ and $V$ are respectively $n \times n$ and $m \times m$ square unitary matrices and $S$ is a $n \times m$ rectangular diagonal matrix whose entries are positive numbers ordered in decreasing order and called the singular values of the matrix.

This definition can be extended to higher-rank tensors by partitioning their indices into two groups forming two super-indices. The tensor can then be represented as a rectangular matrix on which the SVD can be applied. In that case the result of the decomposition depends on the partition of the indices.

Equation (3.9) shows a diagrammatic representation of this decomposition for a rank-3 tensor $T$, along a partition of its indices given by the red dotted line, into the product $USV^\dagger$ where $U^\dagger U = 1$ and $V^\dagger V = 1$ and $S$ is a diagonal matrix whose elements (the singular values of the tensor along the given partition) are positive and arrange in a
attached to the apex of the triangle

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node (X) at (0,0) {$X$}; \node at (-0.5,0.5) {$d_1$}; \node at (0.5,0.5) {$d_2$}; \node at (0,-0.5) {$d_1$}; \draw[thick] (X) -- (-0.5,0.5); \draw[thick] (X) -- (0.5,0.5); \draw[thick] (X) -- (0,-0.5); \end{tikzpicture} = \hat{P}, \qquad (3.13) $$

$$ \begin{tikzpicture}[baseline=(current bounding box.center)] \node (X) at (0,0) {$X$}; \node at (-0.5,0.5) {$d_2$}; \node at (0.5,0.5) {$d_1$}; \node at (0,-0.5) {$d_2$}; \draw[thick] (X) -- (-0.5,0.5); \draw[thick] (X) -- (0.5,0.5); \draw[thick] (X) -- (0,-0.5); \end{tikzpicture} = \hat{\mathcal{I}}. \qquad (3.14) $$

Hence, upon contraction with the corresponding tensors of the associated bra $\langle\psi|$ they will simply reduce to the identity

$$ \langle\psi|\psi\rangle = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (X) at (0,0) {$\cdots$}; \node at (-0.5,0.5) {$\cdots$}; \node at (0.5,0.5) {$\cdots$}; \node at (0,-0.5) {$\cdots$}; \draw[thick] (X) -- (-0.5,0.5); \draw[thick] (X) -- (0.5,0.5); \draw[thick] (X) -- (0,-0.5); \end{tikzpicture} = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (X) at (0,0) {$\cdots$}; \node at (-0.5,0.5) {$\cdots$}; \node at (0.5,0.5) {$\cdots$}; \node at (0,-0.5) {$\cdots$}; \draw[thick] (X) -- (-0.5,0.5); \draw[thick] (X) -- (0.5,0.5); \draw[thick] (X) -- (0,-0.5); \end{tikzpicture} = \dots = 1. \qquad (3.15) $$

By performing a SVD on the left-most tensor, and then contracting the singular matrix $S$ and the right unitary $V^\dagger$ with the tensor on the right, the left-most tensor of the MPS becomes a unitary itself. By repeating this process with the tensor on its right, and then with the next tensor, etc., we arrive at the so called *left canonical gauge* where all the tensors are unitaries

### SVD

$$ |\psi\rangle = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (T1) at (0,0) {$T_1$}; \node at (-0.5,0.5) {$i_1$}; \node at (0.5,0.5) {$i_2$}; \node at (0,-0.5) {$i_3$}; \node at (0,0) {$\cdots$}; \draw[thick] (T1) -- (-0.5,0.5); \draw[thick] (T1) -- (0.5,0.5); \draw[thick] (T1) -- (0,-0.5); \end{tikzpicture} \qquad (3.16) $$

### Contract

$$ = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (U1) at (0,0) {$U_1$}; \node at (-0.5,0.5) {$i_1$}; \node at (0.5,0.5) {$i_2$}; \node at (0,-0.5) {$i_3$}; \node at (0,0) {$\cdots$}; \draw[thick] (U1) -- (-0.5,0.5); \draw[thick] (U1) -- (0.5,0.5); \draw[thick] (U1) -- (0,-0.5); \end{tikzpicture} \qquad (3.17) $$

### SVD

$$ = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (U1) at (0,0) {$U_1$}; \node at (-0.5,0.5) {$i_1$}; \node at (0.5,0.5) {$i_2$}; \node at (0,-0.5) {$i_3$}; \node at (0,0) {$\cdots$}; \draw[thick] (U1) -- (-0.5,0.5); \draw[thick] (U1) -- (0.5,0.5); \draw[thick] (U1) -- (0,-0.5); \end{tikzpicture} \qquad (3.18) $$

$$ = \dots \qquad (3.19) $$

$$ |\psi\rangle = \begin{tikzpicture}[baseline=(current bounding box.center)] \node (U1) at (0,0) {$U_1$}; \node at (-0.5,0.5) {$i_1$}; \node at (0.5,0.5) {$i_2$}; \node at (0,-0.5) {$i_3$}; \node at (0,0) {$\cdots$}; \draw[thick] (U1) -- (-0.5,0.5); \draw[thick] (U1) -- (0.5,0.5); \draw[thick] (U1) -- (0,-0.5); \end{tikzpicture} \qquad (3.20) $$

Of course, we can apply the same procedure starting from the right-

Because of the normalisation of the wave function $\langle\psi|\psi\rangle = 1$, we can show that the last, right-most, tensor is also a unitary.
where *D* is the bond dimension of the MPS. Hence, because *D* is independent of the size *N* of the system, MPSs belong to the class of area law states.

The SVD is not the only way to decompose a tensor into a product of other tensors. For instance, if it were possible to write the considered tensor *T* as a square matrix, we could have simply used a diagonalization *T* = *UAU*†. For rectangular matrices, the QR decomposition can also be used where a *n* × *m*-matrix *M* (*n* > *m*) is decomposed into a *n* × *m* isometry *Q* and a *m* × *m* upper-triangular matrix *R*

$$ \text{--- } n \text{--- } \mathbf{M} \text{--- } m \text{---} = \text{--- } n \text{--- } \mathbf{Q} \text{--- } m \text{--- } \mathbf{R} \text{--- } m \text{---} \quad (3.27) $$

The QR decomposition is useful when one is mainly interested in the unitarity properties and not in the Schmidt coefficients, as this decomposition is computationally faster than the SVD.

## 3.4 Matrix Product Operators

In order to compute expectation values of observables or apply unitary transformations to a quantum state, we need a TN representation of operators. In the same fashion as a one-dimensional quantum state can be represented as a MPS, operators acting on those states can be represented as *Matrix Product Operators* (MPO). For an operator $\hat{O}$, its MPO can be defined as follows

**Definition 3.4.1** The decomposition of the tensor of the coefficients of an operator acting on a quantum state into a product of smaller rank tensors is called a *Matrix Product Operator* decomposition.

$$ \hat{O} = \sum_{\{i_k\} \{i'_k\} \{w\}} W_{1 w_0 w_1}^{i_1 i'_1} \cdots W_{N w_{N-1} w_N}^{i_N i'_N} |\phi_{i'_1} \cdots \phi_{i'_N}\rangle \langle \phi_{i_1} \cdots \phi_{i_N}| \quad (3.28) $$

The contracted indices between the tensors are called virtual indices. The free indices are called physical indices and correspond to the different input and output local quantum states. They can take *d* values (with *d* the dimension of the local Hilbert space).

A diagrammatic formulation of a MPO representation of an operator is

$$ \hat{O} = -1 - \begin{array}{c|c|c|c} \multicolumn{4}{c}{\begin{array}{c} i_1 \\ \hline i'_1 \end{array}} \\ \multicolumn{4}{c}{\begin{array}{c} i_2 \\ \hline i'_2 \end{array}} \\ \multicolumn{4}{c}{\begin{array}{c} i_3 \\ \hline i'_3 \end{array}} \\ \multicolumn{4}{c}{\begin{array}{c} \vdots \\ \hline \vdots \end{array}} \\ \multicolumn{4}{c}{\begin{array}{c} i_N \\ \hline i'_N \end{array}} \end{array} - 1 - \begin{array}{c} w_1 \\ \hline w_2 \end{array} - \begin{array}{c} w_2 \\ \hline w_3 \end{array} - \begin{array}{c} \vdots \\ \hline \vdots \end{array} - \begin{array}{c} w_{N-1} \\ \hline w_N \end{array} - \begin{array}{c} w_N \end{array} $$
where the convention to put free legs of dimension 1 to the first and last tensors have been used to have only rank-4 tensors. These extra free legs are equivalent to no legs.

The action of an operator on a quantum state can then be represented as the action of a MPO on a MPS and the resulting state is given by the contraction of the local tensors

<figure><img src="image_24.png" /></figure>

Note that during the contraction of the physical legs of the MPS and MPO, the dimensions of their respective virtual indices multiply. This is the reason why we need the restricted rank approximation introduced above. Without it, the bond dimension of a MPS will grow exponentially during the time evolution performed with a MPO. In Sec. 3.5 we will show how to write a MPO for the time evolution. To do so, we first need to be able to write the Hamiltonian as a MPO. Hence, we present a general method to construct the MPO representation of a Hamiltonian.

## Constructing a MPO for a Hamiltonian

To construct the MPO representation of a Hamiltonian $\hat{H}$ which is made of a sum of local terms, we use a method based on the recurrence relation presented in Ref. [60].

To define the $k^{th}$ tensor of the MPO, we have to decompose the Hamiltonian into

* a part $\hat{H}_{k-1}^L$ that describes what happens before the bond $k$

* a part $\hat{H}_{k+1}^R$ describing what happens after the bond $k$,

* a part $\sum_a \hat{h}_{k a}^L \otimes \hat{h}_{k a}^R$ which is an interaction Hamiltonian between the part of the system on the left of bond $k$ and the one on the right of bond $k$.

The Hamiltonian thus reads

$$ \hat{H} = \hat{H}_{k-1}^L \otimes \hat{1}_k^R + \hat{1}_k^L \otimes \hat{H}_{k+1}^R + \sum_a \hat{h}_{k a}^L \otimes \hat{h}_{k a}^R \quad (3.31) $$
where $\hat{\mathbf{1}}_k^R = \underbrace{\hat{\mathbf{1}} \otimes \dots \otimes \hat{\mathbf{1}}}_{N-k+1 \text{ times}}$ and $\hat{\mathbf{1}}_k^L = \underbrace{\hat{\mathbf{1}} \otimes \dots \otimes \hat{\mathbf{1}}}_{k \text{ times}}$. Hence $\hat{h}_k^L$ contains an

operator defined on the left of $k$ and $\hat{h}_k^R$ an operator defined on the right of $k$. For instance, if we consider a XYZ-Hamiltonian with nearest neighbours couplings, we could have $\hat{h}_k^L = J_a \hat{S}_k^a$ and $\hat{h}_k^R = \hat{S}_{k+1}^a$ with $a \in \{x, y, z\}$. A recurrence relation between the 'right parts' of the Hamiltonian at two consecutive sites can be defined

$$ \begin{pmatrix} \hat{H}_k^R \\ \hat{h}_k^R \\ \hat{\mathbf{1}}_k^R \end{pmatrix} = W_{k+1} \begin{pmatrix} \hat{H}_{k+1}^R \\ \hat{h}_{k+1}^R \\ \hat{\mathbf{1}}_{k+1}^R \end{pmatrix}, \qquad (3.32) $$

with the matrices $W_k$ defining the Hamiltonian MPO

$$ \hat{H} = \sum_{\sigma, \sigma', w} W_{1 w_1}^{\sigma_1 \sigma'_1} W_{2 w_2}^{\sigma_2 \sigma'_2} \dots W_{N w_{N-1}}^{\sigma_N \sigma'_N} |\sigma_1 \dots \sigma_N\rangle \langle \sigma'_1 \dots \sigma'_N| . \quad (3.33) $$

In Eq. (3.33) the $\sigma$ and $\sigma'$ indices refer to the local Hilbert spaces of the different parts of the system (e.g. spins in the XYZ-Hamiltonian) whereas the $w$ indices relate to virtual bonds between the different parts of the system.

The on-site tensors $W_k$ have a general structure

$$ W_k = \begin{pmatrix} \hat{\mathbf{1}} & C_k & D_k \\ 0 & \cdots & \cdots \\ \vdots & A_k & B_k \\ 0 & \cdots & \cdots \\ 0 & \cdots & \cdots \\ \vdots & \ddots & \ddots \\ 0 & \cdots & \cdots \\ \hat{\mathbf{1}} \end{pmatrix}, \qquad (3.34) $$

The first and last tensors will be respectively a row-vector and a column-vector, as can be seen in the example below.

where the blocks $A_k, \dots, D_k$ are operator-valued, $\hat{\mathbf{1}}$ is an identity operator on the local Hilbert space of site $k$ and $0$ is a notation for tensor elements equal to zero. The number of columns in $A_k$ and $C_k$ is equal to size of the set ran through by the index $a$ in Eq. (3.31), which also corresponds to the number of lines of $A_k$ and $B_k$. The different blocks can be interpreted as follow:

$D_k$ is the on-site energy.

$C_k$ corresponds to the local contribution of the coupling between site $k$ and sites to its right.

$B_k$ corresponds to the local contribution of the coupling between site $k$ and sites to its left.

$A_k$ corresponds to long-range coupling terms (as we will show in chapter 6).

For example, a XYZ-Hamiltonian with a nearest neighbours interaction on a 1d lattice with an external field $\vec{h} = (h_x, 0, h_z)$

$$ \hat{H}_{XYZ} = \sum_i (J_x \hat{S}_i^x \hat{S}_{i+1}^x + J_y \hat{S}_i^y \hat{S}_{i+1}^y + J_z \hat{S}_i^z \hat{S}_{i+1}^z + h_x \hat{S}_i^x + h_z \hat{S}_i^z), \quad (3.35) $$
will have on-site MPO tensors of the form

$$W_i = \begin{pmatrix} \hat{1}_{2 \times 2} & J_x \hat{S}_i^x & J_y \hat{S}_i^y & J_z \hat{S}_i^z & h_x \hat{S}_i^x + h_z \hat{S}_i^z \\ 0 & 0 & 0 & 0 & \hat{S}_i^x \\ 0 & 0 & 0 & 0 & \hat{S}_i^y \\ 0 & 0 & 0 & 0 & \hat{S}_i^z \\ 0 & 0 & 0 & 0 & \hat{1}_{2 \times 2} \end{pmatrix}, \quad (3.36)$$

that has a bond dimension $D = 5$ and a physical dimension (dimension of the local Hilbert space) $d = 2$.

We can see how the contraction of the {$W_i$}$_{i=1...N}$ gives back the Hamiltonian by looking at a three-spins XYZ-chain. Starting from the left, the on-site tensors contract as follow

$$W_1 W_2 = (\hat{1}_{2 \times 2} \quad J_x \hat{S}_1^x \quad J_y \hat{S}_1^y \quad J_z \hat{S}_1^z \quad h_x \hat{S}_1^x + h_z \hat{S}_1^z) \begin{pmatrix} \hat{1}_{2 \times 2} & J_x \hat{S}_2^x & J_y \hat{S}_2^y & J_z \hat{S}_2^z & h_x \hat{S}_2^x + h_z \hat{S}_2^z \\ 0 & 0 & 0 & 0 & \hat{S}_2^x \\ 0 & 0 & 0 & 0 & \hat{S}_2^y \\ 0 & 0 & 0 & 0 & \hat{S}_2^z \\ 0 & 0 & 0 & 0 & \hat{1}_{2 \times 2} \end{pmatrix} \quad (3.37)$$

$$= (\hat{1}_{2 \times 2} \quad J_x \hat{S}_2^x \quad J_y \hat{S}_2^y \quad J_z \hat{S}_2^z \quad h_x \hat{S}_2^x + h_z \hat{S}_2^z + J_x \hat{S}_1^x \hat{S}_2^x + J_y \hat{S}_1^y \hat{S}_2^y + J_z \hat{S}_1^z \hat{S}_2^z + h_x \hat{S}_1^x + h_z \hat{S}_1^z). \quad (3.38)$$

The resulting tensor has the same shape as $W_1$ with its C-block being $C_2$ and its D-block corresponding to both the on-site energies of the first and second spins and their interaction term. Contraction with the last tensor $W_3$ is simply a dot product of two operator-valued vectors

$$W_1 W_2 W_3 = (\hat{1}_{2 \times 2} \quad J_x \hat{S}_2^x \quad J_y \hat{S}_2^y \quad J_z \hat{S}_2^z \quad h_x \hat{S}_2^x + h_z \hat{S}_2^z + J_x \hat{S}_1^x \hat{S}_2^x + J_y \hat{S}_1^y \hat{S}_2^y + J_z \hat{S}_1^z \hat{S}_2^z + h_x \hat{S}_1^x + h_z \hat{S}_1^z) \begin{pmatrix} h_x \hat{S}_3^x + h_z \hat{S}_3^z \\ \hat{S}_3^x \\ \hat{S}_3^y \\ \hat{S}_3^z \\ \hat{1}_{2 \times 2} \end{pmatrix} \quad (3.39)$$

$$= h_x \hat{S}_3^x + h_z \hat{S}_3^z + J_x \hat{S}_2^x \hat{S}_3^x + J_y \hat{S}_2^y \hat{S}_3^y + J_z \hat{S}_2^z \hat{S}_3^z + h_x \hat{S}_2^x + h_z \hat{S}_2^z + J_x \hat{S}_1^x \hat{S}_2^x + J_y \hat{S}_1^y \hat{S}_2^y + J_z \hat{S}_1^z \hat{S}_2^z + h_x \hat{S}_1^x + h_z \hat{S}_1^z \quad (3.40)$$

$$W_1 W_2 W_3 = \hat{H}_{XYZ}. \quad (3.41)$$

The operator recovered in Eq. (3.40) is the Hamiltonian in Eq. (3.35) for three spins.

The interpretation of the blocks inside of the site tensors of the MPO in Eq. (3.34) allows them to be constructed 'by hand' for local interactions. We will make use of this construction in chapters 5, 6 and 7.

## 3.5 Time-evolution methods

To study the dynamics of a quantum state $|\psi(t)\rangle$ starting from a given initial condition $|\psi(0)\rangle$, one has to be able to perform time-evolution
on a MPS. The general solution of the Schrödinger equation for a time-independent Hamiltonian is ($\hbar = 1$)

$$|\psi(t)\rangle = e^{-i\hat{H}t} |\psi(0)\rangle. \qquad (3.42)$$

The straightforward solution would thus be to write down the evolution operator $\hat{U}(t) = \exp(-i\hat{H}t)$ as a MPO. However, this is unpractical for long time simulations as the size of the bond dimension would become too large. A solution is to split up $\hat{U}(t)$ into a product of N one-time-step evolution operators $\hat{U}(\delta t)$, with $\delta t = t/N$, and truncate the MPS bond dimension after each application

$$\hat{U}(t) = e^{-i\hat{H}N\delta t} = \left(e^{-i\hat{H}\delta t}\right)^N = \hat{U}(\delta t)^N. \qquad (3.43)$$

The *Time-Evolving Block Decimation* (TEBD) method is based on this approach [53, 57]. It is adapted for Hamiltonians with short-range interactions<fnref n="10" /> that can be split into internally commuting parts. For example, let us say that the Hamiltonian $\hat{H}$ can be decomposed into two terms $\hat{H} = \hat{H}_1 + \hat{H}_2$ where the terms forming each $\hat{H}_i$ commute with one another but $\hat{H}_1$ and $\hat{H}_2$ don't necessarily commute. Then, the single time step evolution operator can be split up using a Suzuki-Trotter decomposition [66]

$$\hat{U}(\delta t) = e^{-i\hat{H}_1\delta t} e^{-i\hat{H}_2\delta t} + \mathcal{O}(\delta t^2). \qquad (3.44)$$

Each exponential can then be written in the form of a MPO. The error induced by the Trotterization of the exponential can be decreased to $\mathcal{O}(\delta t^3)$ by symmetrising the decomposition. The impact of this time step error on the unitarity of the time evolution can be controlled by choosing a sufficiently small $\delta t$. The other source of error originates from the truncation of the bond dimension performed after each time step. This truncation affects directly the unitarity of the time-evolution but convergence can be easily controlled by varying the number of kept singular values.

Other methods share a similar philosophy: find an approximate formulation of the single-time-step evolution operator $\hat{U}(\delta t)$ and write it in the form of an MPO. For example, the $W^I$ method relies on a second order Taylor expansion of $\hat{U}(\delta t)$ where terms acting on overlapping supports are being discarded, i.e. products of the type $\hat{H}_i\hat{H}_j$, with $\hat{H} = \sum_i \hat{H}_i$, are set to zero if $\hat{H}_i$ and $\hat{H}_j$ act on a common part of the system [59]. This method is simple and has the advantage of having an implementation that is agnostic of the particular problem at hand; it only requires the knowledge of the MPO representation of $\hat{H}$. However, it suffers from strong limitations and is not able to accurately compute the evolution generated by single-site Hamiltonians [60]. An improvement named $W^{II}$ takes into account single-site operator overlapping but requires a more complex implementation [59].

Another set of time evolution methods takes a different perspective. Instead of finding an efficient MPO representation $\hat{U}(\delta t)$ they aim at approximating the *action* of the time evolution operator on a quantum state. This means that they offer a way to evolve $|\psi(t)\rangle$ into $|\psi(t + \delta t)\rangle$ without constructing an explicit representation of $\hat{U}(\delta t)$. The two main

<figure><img src="image_25.png" /><figcaption>Figure 3.5: Matsuo Suzuki (born on March 3, 1937) is an emeritus professor at the University of Tokyo. He has been working on critical phenomena and phase transitions in spin systems.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 3.6: Hale Freeman Trotter (30 May 1931 – 17 January 2022) was a Canadian-American mathematician known for the Lie–Trotter product formula.</figcaption></figure>

<footnote n="10">¹⁰ Long-range interaction can also be described with TEBD using a succession of swap gates to exchange local terms together in order to effectively write the long-range interaction as a local one [65]. Nonetheless, this increases the complexity of the method.</footnote>
If $N = \dim(\mathcal{H})$ the Krylov subspace encompasses the whole Hilbert space.

The term 'variational' in the name of the method comes from the fact that in practice one aims at minimising the right-hand side of Eq. (3.45) to find $|\varphi\rangle$.

<figure><img src="image_27.png" /><figcaption>**Figure 3.7:** Nikolay Mitrofanovich Krylov (29 November 1879 – 11 May 1955) was a Soviet mathematician known for his works in mathematical physics, numerical methods and mechanics. He was the PhD supervisor of Nikolay Bogoliubov.</figcaption></figure>

methods representative of this approach are the global Krylov method [60, 67] and the time-dependent variational principle (TDVP). The global Krylov method approximates the action of $\hat{U}(\delta t)$ by restricting oneself to the so called Krylov subspace $\mathcal{H}_N$ which is the space spanned by the set of vectors $\{|\psi(0)\rangle, \hat{H} |\psi(0)\rangle, \dots, \hat{H}^N |\psi(0)\rangle\}$. Defining the Krylov vectors $\{|k_1\rangle, \dots, |k_N\rangle\}$, which form an orthonormal basis of $\mathcal{H}_N$, the method consists in searching for the vector $|k_i\rangle$ that is the closest approximation to $\hat{U}(\delta t) |\psi(t)\rangle$.

We now focus more specifically on the TDVP method as it is the one we have been using in our work.

## Time-Dependent Variational Principle

The original idea behind TDVP goes back to Dirac [68] and Frenkel [69]. The main point, in the modern tensor networks formulation, is that instead of solving the Schrödinger equation and then truncating the MPS representation of the quantum state, one can solve the equations of motion projected into a space of restricted bond dimension [58, 70]. The general formulation of the Dirac-Frenkel Variational Principle [71] is that one looks for a solution $|\varphi\rangle \in \mathcal{M}$ of the Schrödinger equation where $\mathcal{M} \subset \mathcal{H}$ is a manifold of the total Hilbert space $\mathcal{H}$ in which we think that the relevant physical states 'live'. We define $T_{|\varphi\rangle} \mathcal{M}$ the tangent space of $\mathcal{M}$ around the state $|\varphi\rangle$. The criterion to find $|\varphi\rangle$ is that for every state $|\chi\rangle \in T_{|\varphi\rangle} \mathcal{M}$

$$ \langle \chi | \left( \frac{d}{dt} - \frac{1}{i\hbar} \hat{H} \right) | \varphi \rangle = 0, \qquad (3.45) $$

which can be interpreted as saying that the time evolution procedure should keep $|\varphi\rangle$ inside of the manifold $\mathcal{M}$. Introducing $\hat{P}_{T_{|\varphi\rangle} \mathcal{M}}$ the projector onto the tangent space $T_{|\varphi\rangle} \mathcal{M}$, we can write the state $|\chi\rangle = \hat{P}_{T_{|\varphi\rangle} \mathcal{M}} |\varphi\rangle$ with $|\varphi\rangle$ a state in $\mathcal{H}$. Leading to

$$ \forall |\phi\rangle \in \mathcal{H}, \quad \langle \phi | \hat{P}_{T_{|\phi\rangle} \mathcal{M}} \left( \frac{d}{dt} - \frac{1}{i\hbar} \hat{H} \right) |\phi\rangle = 0. \qquad (3.46) $$

Because the time derivation and the projector commute, we have

$$ \forall |\phi\rangle \in \mathcal{H}, \quad \langle \phi | \left( \frac{d}{dt} - \frac{1}{i\hbar} \hat{P}_{T_{|\phi\rangle} \mathcal{M}} \hat{H} \right) |\phi\rangle = 0. \qquad (3.47) $$

This equation must be true for any $|\phi\rangle \in \mathcal{H}$, Eq. (3.45) can thus be written

$$ \left( \frac{d}{dt} - \frac{1}{i\hbar} \hat{P}_{T_{|\phi\rangle} \mathcal{M}} \hat{H} \right) |\phi\rangle = 0. \qquad (3.48) $$

In the context of MPS, the manifold $\mathcal{M}$ will correspond to the space of full-ranked MPS of a given bond dimension $D$, and the tangent space will be the space spanned by variations of single MPS tensors. The major advantage of this method is that it naturally preserves the unitarity of the time evolution and conserves the energy.
At first glance, Eq. (3.48) is quite a complex equation as the projector itself depends on $|\varphi\rangle$. However it has a simple decomposition [70] that can be written easily in a tensor network form

$$ \hat{P}_{T_{|\varphi\rangle}} = \sum_{i=1}^{N} \hat{P}_{|\varphi\rangle,i-1}^{L} \otimes \hat{1}_{i} \otimes \hat{P}_{|\varphi\rangle,i+1}^{R} - \sum_{i=1}^{N-1} \hat{P}_{|\varphi\rangle,i}^{L} \otimes \hat{P}_{|\varphi\rangle,i+1}^{R}, \quad (3.49) $$

where $\hat{P}_{|\varphi\rangle,i}^{L/R}$ projects on all the sites left (right) of site i included. The first terms of Eq. (3.49) neutralises states that differ from $|\varphi\rangle$ by more than one tensor. The second term is there to ensure normalisation of the projected states by removing contributions in excess that correspond exactly to $|\varphi\rangle$. These left and right projectors can be written as TN in the left and right canonical gauges

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

The projection of the application of the Hamiltonian on $|\varphi\rangle$ results in $2N-1$ terms. They can be decomposed for each site i and bond connecting sites i and i + 1, hence giving a set of local equations for orthogonality centres $T_i^C$ at site i and the corresponding bond-centred matrix $S_i$

$$ \frac{d}{dt} T_i^C = -i \hat{H}_{\text{eff}}^i T_i^C \quad \text{for } i \in [1, N], \qquad (3.52) $$

$$ \frac{d}{dt} S_i = +i \hat{K}_{\text{eff}}^i S_i \quad \text{for } i \in [1, N-1], \qquad (3.53) $$

with $\hat{H}_{\text{eff}}^i$ and $\hat{K}_{\text{eff}}^i$ that can be obtained by contractions of the dashed boxes in Eqs. (3.54) & (3.55)

<figure><img src="image_30.png" /><figcaption>Figure 3.8: Paul Adrien Maurice Dirac (8 August 1902 – 20 October 1984) was an English theoretical physicist who made fundamental contributions to the early development of both quantum mechanics and quantum electrodynamics. He was awarded the Nobel prize in Physics with Schrödinger "for the discovery of new productive forms of atomic theory".</figcaption></figure>

In our implementation of TDVP, we actually perform the gauge transformations with a QR decomposition instead of a SVD, hence the unitaries of Eqs. (3.54) & (3.55) are replaced by isometries and the bond-centred matrices $S_i$ are not diagonal.
$$ \hat{P}_{|\varphi\rangle,i-1}^L \otimes \hat{1}_i \otimes \hat{P}_{|\varphi\rangle,i+1}^R \hat{H} |\varphi\rangle = \left. \begin{array}{c} \text{def.} \\ \equiv \hat{H}_{\text{eff}}^i \end{array} \right\} , \quad (3.54) $$

$$ \hat{P}_{|\varphi\rangle,i}^L \otimes \hat{P}_{|\varphi\rangle,i+1}^R \hat{H} |\varphi\rangle = \left. \begin{array}{c} \text{def.} \\ \equiv \hat{K}_{\text{eff}}^i \end{array} \right\} . \quad (3.55) $$

Formally, these $2N-1$ equations are coupled as $\hat{H}_{\text{eff}}$ and $\hat{K}_{\text{eff}}$ depend on $|\varphi\rangle$. However, we can still solve these equations sequentially and update $\hat{H}_{\text{eff}}$ and $\hat{K}_{\text{eff}}$ after each site evolution. This procedure is known as one-site TDVP (1TDVP) [60].

The error induced in this method has three contribution:

(1) the projection error resulting from restricting the dynamics to the manifold $\mathcal{M}$. However, this error is easily controllable by increasing the bond dimension of the MPS until convergence is reached.

(2) the finite time step error which is of order $\mathcal{O}(\delta t^3)$.

(3) the error coming from considering the $2N-1$ equations as independent.

Alternative formulations of TDVP that evolve several sites at once at each time step also exist. The manifold considered in those cases include MPS of larger bond dimension than the initial one. The most famous several site TDVP is the two-site implementation (2TDVP). This implementation of TDVP has the 'advantage' of allowing the bond dimension to evolve dynamically during the time evolution but its complexity scales poorly with the physical dimension of the MPS tensors.

Recently, a new implementation of TDVP combining the simplicity of 1TDVP with the flexibility of the bond dimension of 2TDVP has been developed under the name Adaptive 1TDVP (A1TDVP) [72]. This new method is able to increase the bond dimension of the MPS 'on the fly' while performing one-site evolution of the MPS. It relies on

<figure><img src="image_31.png" /><figcaption>**Figure 3.9:** Yakov Il'ich Frenkel (10 February 1894 – 23 January 1952) was a Soviet physicist who worked in condensed matter.</figcaption></figure>
performing sub-space expansion (i.e. adding basis vectors to the space
under consideration) in the one-site projectors $\hat{P}_{|\phi\rangle,i}^{L/R}$. These projectors
now project to a manifold of higher bond dimension.

In the next chapter, we will see how to naturally give a one-dimensional
geometry to a joint {System + Environment} state. This will enable the
construction of a MPS representation of the wave function of an open
quantum system.
<figure><img src="image_32.png" /></figure>
This form for $\hat{H}_{\text{int}}$ is the simplest form of linear coupling. More complex interactions can have the system's operator depend on the bath or the bath's operator depend on the system as we will see later.

**Figure 4.1:** Geometry of the system interacting with a bosonic environment. Black lines represent interactions. (a) For a continuum of normal modes depicted as a green shade. (b) For a finite set of modes, depicted as green discs. This structure where the system couples to all the modes and the modes only to the system is called the *star geometry*.

## 4.1 Principle of the chain mapping

Let us consider a bosonic bath at zero temperature coupled to a system via an operator $\hat{O}$. The bath free Hamiltonian and the interaction Hamiltonian are

$$ \hat{H}_B = \int_0^{+\infty} \hbar\omega \hat{a}_\omega^\dagger \hat{a}_\omega \mathrm{d}\omega, \quad (4.1) $$

$$ \hat{H}_{\text{int}} = \hat{O} \int_{0}^{+\infty} g_{\omega}(\hat{a}_{\omega} + \hat{a}_{\omega}^{\dagger}) \mathrm{d}\omega, \quad (4.2) $$

where $\hat{a}_\omega^\dagger$ creates an excitation in the bath mode of energy $\hbar\omega$ and $g_\omega$ is the coupling strength between the system and the bath mode of angular frequency $\omega$. In practice, there exists a maximal frequency of the bath $\omega_c$ that can be excited through its interaction with the system. These Hamiltonians define a geometry of the system and the bath interaction. The bath modes are independent (normal) modes and the system interacts with all of them. Figure 4.1 shows schematics of this interaction geometry in the cases of a continuum of modes and a finite number of modes. This type of geometry is called a *star geometry*.

<figure><img src="image_33.png" /></figure>

We can introduce the following unitary transformation of the bath creation operators

$$ \hat{a}_{\omega}^{\dagger} = \sum_{n=0}^{\infty} U_{n}(\omega) \hat{b}_{n}^{\dagger} \quad (4.3) $$

where $U_n(\omega)$ is defined with the real orthonormal polynomials $P_n$

$$ U_n(\omega) = g_\omega P_n(\omega). \quad (4.4) $$

Equation (4.3) expresses the decomposition of the continuum of independent bosonic modes labelled by $\omega \in \mathbb{R}^+$ onto a new infinite set of discrete modes labelled by $n \in \mathbb{N}$. This is a reflection of the separability of the bath Hilbert space [37]. Because the transformation from $\{\hat{a}_\omega\}_{\omega \in \mathbb{R}^+}$ to $\{\hat{b}_n\}_{n \in \mathbb{N}}$ is unitary, the bosonic canonical commutation relations are satisfied by $\hat{b}_n$ and $\hat{b}_n^\dagger$. The unitarity of the transformation imposes an orthogonality relation for the polynomials

$$ \int_{0}^{+\infty} d\omega U_{n}(\omega) U_{m}(\omega) = \int_{0}^{+\infty} d\omega g_{\omega}^{2} P_{n}(\omega) P_{m}(\omega) = \delta_{n,m}. \quad (4.5) $$

This orthogonality relation defines the family of polynomials used for the transformation. Thus, the chosen polynomials depend on the bath spectral density $J(\omega) = g_\omega^2$. Another useful property of these
polynomials is that they obey a recurrence relation

$$P_n(\omega) = (C_{n-1}\omega - A_{n-1})P_{n-1}(\omega) + B_{n-1}P_{n-2}(\omega), \quad (4.6)$$

where $A_n$ is related to the first moment of $P_n$, $B_n$ and $C_n$ to the norms of $P_n$ and $P_{n-1}$ [37, 73]. This recurrence relation can be used to construct

the polynomials with the conditions that $P_0(\omega) = \|p_0\|^{-1} = \left(\int_{\mathbb{R}^+} J(\omega)d\omega\right)^{-\frac{1}{2}}$ and $P_{-1}(\omega) = 0$, with $\|\cdot\|$ the norm of $\cdot$ with respect to the measure $J(\omega)$, and $P_n(\omega) = p_n(\omega)\|p_n\|^{-1}$.

If we apply this unitary transformation to the interaction Hamiltonian

The polynomials {$p_n$}$_{n \in \mathbb{N}}$ are the so called monic polynomials where the factor in front of $\omega^n$ is equal to 1.

$$\hat{H}_{\text{int}} = \hat{O} \int_{0}^{+\infty} g_{\omega} \sum_{n} U_{n}(\omega)(\hat{b}_{n} + \hat{b}_{n}^{\dagger})d\omega \quad (4.7)$$

$$= \sum_{n} \hat{O} \int_{0}^{+\infty} J(\omega) P_{n}(\omega)(\hat{b}_{n} + \hat{b}_{n}^{\dagger})d\omega \quad (4.8)$$

$$= \sum_{n} \hat{O} \underbrace{\left( \int_{0}^{+\infty} J(\omega) P_{n}(\omega) P_{0}(\omega) d\omega \right)}_{\delta_{n,0}} \|p_{0}\|(\hat{b}_{n} + \hat{b}_{n}^{\dagger}) \quad (4.9)$$

$$\hat{H}_{\text{int}} = \|p_0\| \hat{O}(\hat{b}_0 + \hat{b}_0^\dagger) \quad (4.10)$$

we obtain a new expression where the system couples only to the first mode with the coupling strength $\|p_0\| \stackrel{\text{def.}}{=} c_0$.

The same transformation applied to the bath Hamiltonian yields to the following nearest neighbours hopping Hamiltonian where $\omega_n = \frac{A_n}{C_n}$ is the energy of the chain mode $n$ and $t_n = C_n^{-1}$ is the coupling between mode $n$ and $n+1$ [73]

The chain coefficients $\omega_n$, $t_n$ and $c_0$ can sometimes by calculated analytically, for instance at $\beta = \infty$ for the family of Ohmic SD defined below, otherwise they are computed numerically using an implementation of the ORTHPOL package [75] in Julia [76].

$$\hat{H}_B = \sum_n \omega_n \hat{b}_n^\dagger \hat{b}_n + t_n (\hat{b}_n^\dagger \hat{b}_{n+1} + \hat{b}_{n+1}^\dagger \hat{b}_n). \quad (4.11)$$

From the new bath and interaction Hamiltonians of Eqs. (4.10) and (4.11) we can see that the unitary transformation $U_n(\omega)$ transforms the bosonic environment composed of a continuum of independent modes — the star environment — into a semi-infinite chain of interacting modes (see Fig. 4.2) [23].

<figure><img src="image_34.png" /></figure>

The new chain coefficients $c_0$, $\omega_n$ and $t_n$ are all defined by the bath spectral density. However, the asymptotic behaviours of the energy $\omega_n$ and the hopping term $t_n$ for $n \to \infty$ are universal and only depend on the support of the measure $J(\omega)$

<caption>**Figure 4.2:** Schematic drawing of the geometry of the system interacting with the new semi-infinite chain environment. Black lines represent interactions. This new geometry is one-dimensional.</caption>

$$\lim_{n \to \infty} \omega_n = \frac{\omega_c}{2} = \lim_{n \to \infty} 2t_n, \quad (4.12)$$

where $\omega_c$ is the maximal bath frequency coupled to the system. Figure 4.3 shows the distribution of the chain coefficients for an Ohmic SD.
For an Ohmic SD ($s = 1$) the coupling strength $\alpha$ is a dimensionless parameter.

**Definition 4.1.1** The family of Ohmic spectral densities is defined as follow

$$J(\omega) = 2\alpha \frac{\omega^s}{\omega_c^{s-1}} f_{\omega_c}(\omega), \qquad (4.13)$$

where $\alpha$, called the Kondo parameter, is the strength of the system-bath coupling, $\omega_c$ is the bath cut-off frequency, $f_{\omega_c}(\omega)$ is the a cut-off function, and $s$ is the Ohmicity. When $s = 1$, the SD is called Ohmic; for $s < 1$ it is sub-Ohmic; and for $s > 1$ it is super-Ohmic.

We can see in Fig. 4.3 that the chain coefficients quickly become uniform, defining the so called translationally invariant part of the chain. This also means that the bath-specific physics of the system-bath interaction happens in the first few modes of the chain.

<figure><img src="image_35.png" /><figcaption>**Figure 4.3:** Chain modes energy $\omega_n$ and tunneling $t_n$ for an Ohmic spectral density with a hard cut-off. These chain coefficients are independent of the Kondo parameter $\alpha$.</figcaption></figure>

The nature of the orthonormal polynomials, and hence of the chain coefficient, depends on the cut-off function $f_{\omega_c}(\omega)$. For a Heaviside step function (hard cut-off), the corresponding polynomials will be Jacobi polynomials. Whereas for an exponential cut-off the polynomials will be Laguerre polynomials.

One can show that for Ohmic spectral densities, the chain modes' on-site and hopping energies are independent of the coupling strength $\alpha$. Initially the system was coupled to all the normal bath modes, and after this transformation it is only coupled to the first chain mode $n = 0$. An excitation injected into this mode, i.e. the system dissipating energy into the environment, can then travel along the chain as a wavefront.

At zero temperature, this chain-mapped environment is well-suited for a representation of the joint {System + Bath} wave function as a MPS because the bath is now made of discrete modes, and all the couplings of the joint system are local. One just needs to truncate the chain keeping only the first $N_m$ modes. The only constraint is that this
number should be large enough so that the wavefront generated by the interaction with the system does not reach the end of the chain. Retaining only a finite number of modes in the chain representation also corresponds to considering a finite set of normal modes in the original bath. Truncating the chain can thus be seen as an optimal sampling procedure where the bath modes are not sampled uniformly. These sampled bath modes can be recovered by diagonalizing the tri-diagonal symmetric $N_m \times N_m$-matrix $h_B$

$$h_B = \begin{pmatrix}
\omega_0 & t_0 & 0 & 0 & \dots & 0 \\
t_0 & \omega_1 & t_1 & 0 & \dots & 0 \\
0 & t_1 & \omega_2 & t_2 & \dots & 0 \\
\vdots & & \ddots & & \ddots & \vdots \\
0 & 0 & 0 & 0 & t_{Nm-2} & \omega_{Nm-1}
\end{pmatrix}, \quad (4.14)$$

defined such that

$$\hat{H}_B^{\text{truncated}} = \hat{b}^{\dagger} h_B \hat{b} \quad \text{with} \quad \hat{b} = \begin{pmatrix} \hat{b}_0 \\ \hat{b}_1 \\ \vdots \\ \hat{b}_{N_m-1} \end{pmatrix}. \quad (4.15)$$

The combination of the chain mapping and a MPS representation gives the method known as Time Evolving Density Operator with Orthonormal Polynomial (TEDOPA) [48, 74]. This method is very efficient for simulating OQS non-perturbatively and in the non-Markovian regime as it allows the simulation of the evolution of the full wave-function of the system and its environment.

## 4.2 Extension to finite temperature

The chain mapping allows the quantum state to have a natural representation as a MPS at $\beta = \infty$ when the initial state of the environment is a vacuum state. However, when one wants to describe an environment at finite temperature ($\beta \neq \infty$), the joint state is not a wave function anymore and must be written as a density matrix. Especially, we are interested in initial states of the following form

$$\hat{\rho}(0) = |\psi_S(0)\rangle \langle \psi_S(0)| \otimes \frac{1}{Z} \exp(-\beta \hat{H}_B), \quad (4.16)$$

where $|\psi_S(0)\rangle$ is an initial pure state for the system, and $Z = \text{tr}[\exp(-\beta \hat{H}_B)]$ is the partition function of the bath. This initial state $\hat{\rho}(0)$ is a product state between the system and its environment where the latter is in a Gibbs state at the inverse temperature $\beta$. A natural approach there could be to write the density matrix as a MPO [77] and use methods originating from DMRG to time-evolve it but the scaling of these methods with the local Hilbert's space dimensions is fairly poor. Other methods based on an expression of the bath influence functional as
a tensor network are able to handle (reduced) density matrices and finite temperatures [46, 47].

Nevertheless, in some cases TEDOPA can be extended to finite temperature [78]. As we have seen in Sec. 2.4, for harmonic baths the dynamics of the reduced system only depends on the bath auto-correlation function $C(\tau)$. The exact nature of the environment is not essential to calculate the system's observables, thus the real microscopic bath can be replaced by a fictitious one that is more convenient. Hence, if we are able to find a zero temperature bath with the same correlation function as the finite temperature one, then it will be possible to use a MPS representation of the quantum state and use TEDOPA. The finite temperature bath correlation function is

$$C_{\beta}(\tau) = \int_{0}^{\omega_{c}} d\omega J(\omega) \left( n_{\beta}(\omega)e^{i\omega\tau} + (n_{\beta}(\omega) + 1)e^{-i\omega\tau} \right). \quad (4.17)$$

Whereas, for $\beta = \infty$, the correlation function reduces to

$$C_{\infty}(\tau) = \int_{0}^{\omega_{c}} d\omega J(\omega)e^{-i\omega\tau}. \quad (4.18)$$

We want to rewrite $C_{\beta}(\tau)$ in the same form as Eq. (4.18), i.e. the integral of a spectral density multiplied by a time-dependent phase factor. In other words, we want to find a bath at zero $T$ with a different spectral density but with the same system dynamics as the finite $T$ bath.

We recast the first term of $C_{\beta}(\tau)$ in Eq. (4.17) such that the argument of the exponential is the same as the second term by sending $\omega \to -\omega$, allowing for negative frequencies and using the identity

$$n_{\beta}(-|\omega|) = -(n_{\beta}(|\omega|) + 1). \quad (4.19)$$

The correlation function thus reads

$$C_{\beta}(\tau) = \int_{-\omega_c}^{0} d\omega \operatorname{sign}(\omega) J(|\omega|) (n_{\beta}(|\omega|) + 1) e^{-i\omega\tau} + \int_{0}^{\omega_c} d\omega J(\omega) (n_{\beta}(\omega) + 1) e^{-i\omega\tau}. \quad (4.20)$$

We define the extended SD $J_{\text{ext}}(\omega) = \operatorname{sign}(\omega)J(|\omega|)$ with a domain extended to negative frequencies. Finally, the bath correlation function can be written

$$C_{\beta}(\tau) = \int_{-\omega_c}^{+\omega_c} d\omega J_{\text{ext}}(\omega) (n_{\beta}(\omega) + 1) e^{-i\omega t} \quad (4.21)$$

where the new bath has an effective temperature-dependent SD

$$J_{\beta}(\omega) = J_{\text{ext}}(\omega)(n_{\beta}(\omega) + 1). \quad (4.22)$$

The behaviour of the thermal SD $J_{\beta}(\omega)$ when varying $\beta$ is shown in Fig. 4.4. We can see that it is always positive and converges to the zero temperature limit over a couple orders of magnitude.
<figure><img src="image_36.png" /><figcaption>**Figure 4.4:** Thermal spectral density $J_β$ in the Ohmic case with $α = 0.01$ and a hard cut-off for different values of $β$. The thermal SD is always positive and its contribution in the negative part of the spectrum decreases rapidly with the inverse temperature $β$.</figcaption></figure>

With this transformation we have, in a sense, double the number of bath modes. The new corresponding free bath Hamiltonian and interaction Hamiltonian are

$$ \hat{H}_B = \int_{-\omega_c}^{\omega_c} d\omega \omega \hat{a}_{\omega}^{\dagger} \hat{a}_{\omega}, \qquad (4.23) $$

$$ \hat{H}_{\text{int}} = \hat{O} \int_{-\omega_c}^{\omega_c} d\omega \sqrt{J_{\beta}(\omega)} (\hat{a}_{\omega} + \hat{a}_{\omega}^{\dagger}) . \qquad (4.24) $$

The new initial state of the bath now corresponds to the vacuum state [78, 79]

$$ \hat{\rho}(0) = (|\psi_S(0)\rangle \otimes |\{0\}\rangle)(\langle\psi_S(0)| \otimes \langle\{0\}|) . \qquad (4.25) $$

This new fictitious bath at zero temperature replaces emission of energy at frequency $ω$ from the bath into the system by an absorption in the bath mode of frequency $-ω$.

We can now define a set of orthonormal polynomials with the associated measure $J_β(ω)$, and follow again the procedure outlined in Sec. 4.1 to map the extended environment to a chain. The combination of the chain mapping technique, the description of a finite temperature environment as a zero temperature one and the description of the joint wave function as a MPS is called Thermalised-TEDOPA (TEDOPA) [78] and enables the efficient description of OQS at finite temperature.

The new temperature-dependent chain parameters have an asymptotic behaviour given by

$$ \lim_{n \to \infty} \omega_n^\beta = \frac{\omega_{\text{Max}} + \omega_{\text{min}}}{2} = 0, \qquad (4.26) $$

$$ \lim_{n \to \infty} t_n^\beta = \frac{\omega_{\text{Max}} - \omega_{\text{min}}}{4} = \frac{\omega_c}{2}, \qquad (4.27) $$
<figure><img src="image_37.png" /><figcaption>**Figure 4.5:** Chain mode energy $\omega_n^\beta$ and tunneling $t_n^\beta$ for an Ohmic spectral density at the inverse temperature $\beta = 5$. These chain coefficients are independent of the Kondo parameter $\alpha$.</figcaption></figure>

where $\omega_{\min/\max}$ is the lowest (largest) frequency that couples to the system. The asymptotic speed of excitations $t_n^\beta$ on the finite-temperature chain is two times larger than the speed for zero-temperature. As a consequence, to simulate the same evolution time, the temperature-dependent chain must be at least twice as long – which is consistent with the fact that the number of normal modes in the environment has been doubled. The distribution of the chain coefficients for an Ohmic SD at $\beta = 5$ can be seen in Fig. 4.5.

## 4.3 Spatially extended systems

The chain mapping procedures presented in the two previous sections were for baths 'indifferent' to the potential spatial structure of the system, i.e. different parts of the system being at different positions in space. This can be seen by the fact that the bath part of the interaction Hamiltonian $\hat{H}_{\text{int}}$ in Eq. (4.2) does not depend on a spatial coordinate. From the point of view of the environment, the system acts as a whole and there is no propagation of bath excitations in space. However, in actual multipartite quantum systems, the different parts are at a given distance from one another and information takes time to travel in space. Indeed, in many cases (for instance, concerted dynamics in biological systems [15]) it is important to be able to take into account the spatial quantum correlations inside of the environment. This means that we have to describe space-dependent interactions and thus adapt our chain mapping technique. Most of the original results to be presented here have been published in Ref. [1].

We consider a one-dimensional bosonic bath with modes characterised by the wave-vectors $k \in [-k_c, +k_c]$, where $k_c$ is the environment cut-off wave-vector. The environment linear dispersion relation is given by $\omega_k = |k|c$ with $c$ the speed of the phonons in the bath. The bath and
interaction Hamiltonians are

$$ \hat{H}_E + \hat{H}_{\text{int}} = \int_{-k_c}^{+k_c} \omega_k \hat{a}_k^\dagger \hat{a}_k dk + \sum_\alpha \hat{O}_\alpha \int_{-k_c}^{+k_c} (g_k^\alpha \hat{a}_k + \text{h.c.}) dk, \quad (4.28) $$

where $\hat{a}_k$ is the annihilation operator of a bath mode of wave-vector $k$, $g_k^\alpha = g_k e^{ikr_\alpha}$, with $g_k = g_{-k} \in \mathbb{R}$, is the coupling strengths between the system and the bath and $r_\alpha$ is the position of the system site $\alpha$,

$\hat{O}_\alpha$ is a system operator acting locally on the system at the position $r_\alpha$. The phase factor $e^{ikr_\alpha}$ in the interaction Hamiltonian takes into account the phase difference between a plane wave generated at a position $r_1$ and one generated at another position $r_2$.

## Zero temperature

We separate positive and negative wave-vector modes and apply to them two different chain mappings, and we introduce $\hat{b}_k \stackrel{\text{def.}}{=} \hat{a}_{-k}$. Two chain mappings are required here because of the bath dispersion relation $\omega_k = |k|c$ which prevents Eq. (4.6) to be successfully applied on the interval $[-k_c, +k_c]$. The bath and interaction Hamiltonians become

$$ \hat{H}_E + \hat{H}_{\text{int}} = \int_0^{+k_c} dk \omega_k (\hat{a}_k^\dagger \hat{a}_k + \hat{b}_k^\dagger \hat{b}_k) \\ + \sum_\alpha \hat{O}_\alpha \int_0^{+k_c} dk g_k (e^{ikr_\alpha} (\hat{a}_k + \hat{b}_k^\dagger) + \text{h.c.}). \quad (4.29) $$

We now introduce two unitary transformations

$$ \hat{a}_{k \ge 0} = \sum_n U_n(k) \hat{c}_n, \quad (4.30) $$

$$ \hat{b}_{k \ge 0} = \sum_m V_m(k) \hat{d}_m, \quad (4.31) $$

where the matrix elements are

$$ U_n(k) = V_n(k) = g_k P_n(k) \quad (4.32) $$

where $\{P_n\}_{n \in \mathbb{N}}$ are orthonormal polynomials with respect to the measure $\mu(k) = |g_k^\alpha|^2 = g_k^2 \stackrel{\text{def.}}{=} J(k)$ (which is the bath spectral density) such that

$$ \int_0^{+k_c} P_n(k) P_m(k) J(k) dk = \delta_{n,m}. \quad (4.33) $$

We can then map the bath Hamiltonian using the unitary transformations from Eqs. (4.30)-(4.31) to two independent tight-binding chains with the same on-site energies $\omega_n$ and hopping energies $t_n$:

$$ \hat{H}_E = \sum_n \omega_n (\hat{c}_n^\dagger \hat{c}_n + \hat{d}_n^\dagger \hat{d}_n) + t_n (\hat{c}_n^\dagger \hat{c}_{n+1} + \hat{c}_{n+1}^\dagger \hat{c}_n + \hat{d}_n^\dagger \hat{d}_{n+1} + \hat{d}_{n+1}^\dagger \hat{d}_n). \quad (4.34) $$
**Figure 4.6:** The unitary transformation $U_n(k)$ transforms a continuous environment of normal $k$-modes to semi-infinite discrete tight-binding chains with interacting $n$-modes. A sub-system at position $r$ is now coupled to the modes of the chains with the coupling strength $\gamma_n(r)$. This new chain-mapping preserves the 1d geometry of the system and the environment but introduces long-range couplings.

<figure><img src="image_38.png" /></figure>

For the interaction Hamiltonian, we apply the same procedure and make use of Eq. (4.6) and find that the chains couple to the system with coupling coefficients $\gamma_n(r_\alpha)$ and $\gamma_n(r_\alpha)^*$

$$ \hat{H}_{\text{int}} = \sum_{\alpha} \hat{O}_{\alpha} \sum_{n} \left( \gamma_{n}(r_{\alpha}) (\hat{c}_{n} + \hat{d}_{n}^{\dagger}) + \text{h.c.} \right) \quad (4.35) $$

where

$$ \gamma_n(r_\alpha) = \int_0^{+k_c} dk J(k) e^{ikr_\alpha} P_n(k) . \quad (4.36) $$

In preceding works, TEDOPA resulted in the system being connected only to the first site of the chain. By contrast, here the system is generally coupled to all the sites of the chain, as represented in Fig. 5.1. Nevertheless, the description of the Hamiltonian as a MPO is not compromised by this change of geometry, as we will show in chapter 6. Another difference between this chain mapping and the previously introduced ones, is that after the mapping the system is now formally coupled to two environments, as can clearly be seen in Fig. 5.1.

We can study the behaviour of these new long-range coupling coefficients for the widely used Ohmic spectral density with a hard cut-off at $k_c$. Figure 4.7 shows the absolute values of the system-chain coupling for zero temperature for an Ohmic spectral density. A sub-system at $k_c R = 0$ couples only to the first site of the chain. However sub-systems at other positions $R$ couple to a range of modes with a maximum strength for the mode $n \sim k_c R/2c$. We see there that a sub-system at a given position couples mostly to a specific portion of the chain and not to the whole chain. Hence, it is still possible to perform numerically exact simulations by truncating the chain. Because of this correspondence between a part of the chain and a region in space, we call this newly introduced chain-mapped environment a "correlated environment".
<figure><img src="image_39.png" /><figcaption>**Figure 4.7:** Absolute value of the system-chain coupling constants, for a bosonic bath with a hard cut-off Ohmic spectral density, as a function of the chain modes *n* and the sites separations *R*. Note that the main peak is centred around *R*/2*c*. Here *α* = 0.12, *c* = 1 and *k*<sub>c</sub> = 1.</figcaption></figure>

Finite temperature

To describe finite temperature systems in a statistical mixture as equiv-
alent zero temperature state vectors, we again apply the same strategy
as the one developed for T-TEDOPA by allowing the bath to have neg-
ative frequency modes to describe thermal fluctuations and using an
alternative bath spectral density that captures the temperature depen-
dence. To identify this new effective spectral density, we put the finite
temperature bath auto-correlation functions C<sub>β</sub>(r,t) for propagating
and counter-propagating modes in the form of a zero temperature
auto-correlation C<sub>∞</sub>(r,t).

The interaction Hamiltonian in interaction picture is

$$
\begin{align}
\hat{H}_{\text{int}}^{I} &= \sum_{\alpha} \hat{O}_{\alpha} \int_{0}^{+k_{c}} dk g_{k} \left( e^{i(kr_{\alpha} - \omega_{k}t)} \hat{a}_{k} + \text{h.c.} \right) \nonumber \\
&\quad + \sum_{\alpha} \hat{O}_{\alpha} \int_{0}^{+k_{c}} dk g_{k} \left( e^{-i(kr_{\alpha} + \omega_{k}t)} \hat{b}_{k} + \text{h.c.} \right) \tag{4.37}
\end{align}
$$

$$
= \sum_{\alpha} \hat{O}_{\alpha} (\hat{B}_{r_{\alpha}}^{1}(t) + \hat{B}_{r_{\alpha}}^{2}(t)) . \qquad (4.38)
$$

Hence the bath correlation function for the propagating modes is

$$
C_{\beta}(r-r',t) = \langle \hat{B}_{r}^{1}(t) \hat{B}_{r'}^{1}(0) \rangle_{B} \quad (4.39)
$$

$$
= \int_{0}^{+k_c} dk J(\omega_k) \left( n_{\beta}(\omega_k) e^{-i(k(r-r')-\omega_k t)} \right. \\
\left. \qquad + (n_{\beta}(\omega_k) + 1) e^{i(k(r-r')-\omega_k t)} \right) . \quad (4.40)
$$

We could write the corresponding correlation function for the counter-
propagating modes which would be the same except for the sign of
the wave-number *k*. These two correlation functions have, in addition
to the usual temperature-dependence *β* and time-dependence *t*, a
spatial-dependence *r* − *r*′ originating from the spatial-dependence of
the coupling coefficients {*g*<sub>*k*</sub><sup>*α*</sup>} between the system and the bosonic bath.
Hence these correlation functions contain more information (i.e. about
space and time) than the correlation functions usually encountered. For zero-temperature, the correlation function reduces to

$$C_{\infty}(r-r',t) = \int_{0}^{+k_c} dk J(\omega_k) e^{i(k(r-r')-\omega_k t)}. \quad (4.41)$$

We want to rewrite $C_{\beta}(r-r',t)$ in the same form as Eq. (4.41), i.e. the integral of a spectral density times a plane-wave phase factor. In other words, we want to find a bath at zero $T$ with a different spectral density but with the same system dynamics as the finite $T$ bath.

We recast the first term of $C_{\beta}(r-r',t)$ in Eq. (4.40) such that the argument of the exponential is the same as the second term by sending $k \to -k$, allowing for negative frequencies (hence, $\omega_{-k} = -\omega_k$) and using the identity in Eq. (4.19). With this transformation we have, in a sense, double the number of propagating modes. There are the propagating positive $k$ modes with positive energies and the propagating negative $k$ modes with negative energies (coming from the second term of the correlation function).

Finally, the bath correlation function for propagating modes can be written

$$C_{\beta}(r-r',t) = \int_{-k_c}^{+k_c} dk \, J_{\text{ext}}(\omega_k)(n_{\beta}(\omega_k) + 1)e^{i(k(r-r')-\omega_k t)} \quad (4.42)$$

with $J_{\text{ext}}$ is the spectral density with a domain extended to negative frequencies and antisymmetrized such that $J_{\text{ext}}(-|\omega_k|) = -J_{\text{ext}}(|\omega_k|)$. The same procedure can be applied to the counter-propagating modes. We can thus define orthonormal polynomials with the finite-temperature spectral density

$$J_{\beta}(k) = J_{\text{ext}}(\omega_k)(n_{\beta}(\omega_k) + 1), \quad (4.43)$$

which is always positive and continuously differentiable. We define the unitary transformation to chain modes

$$\hat{a}_k = \sum_n U_n^\beta(k) \hat{c}_n \quad \text{for } k \in [-k_c, +k_c], \quad (4.44)$$

$$\hat{b}_k = \sum_n U_n^\beta(k) \hat{d}_n \quad \text{for } k \in [-k_c, +k_c] \quad (4.45)$$

where $U_n^\beta(k) = \sqrt{J_\beta(k)} P_n^\beta(k)$ and $P_n^\beta(k)$ is a polynomial of order $n$ from a family of orthonormal polynomials with respect to the measure $d\mu(k) = J_\beta(k)dk$, i.e.

$$\int_{-k_c}^{+k_c} P_n^\beta(k) P_m^\beta(k) d\mu(k) = \delta_{n,m}. \quad (4.46)$$

With this set of orthogonal polynomials, we can map the environment to two tight binding chains and a coupling coefficient

$$\gamma_n(r) = \int_{-k_c}^{+k_c} dk \, J_{\beta}(\omega_k) e^{ikr} P_n^{\beta}(k) \quad (4.47)$$
<figure><img src="image_40.png" /><figcaption>**Figure 4.8:** Absolute value of the system-chain coupling constants at finite temperature, for a bosonic bath with a hard cut-off Ohmic spectral density, as a function of the chain modes *n* and the site separations *R*. The peaks are centred around *n* = *R*/*c*. Here *α* = 0.12, *β* = 0.5, *c* = 1 and *k*<sub>c</sub> = 1.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Figure 4.9:** Absolute value of the system-chain coupling constants at finite temperature, for a bosonic bath with a hard cut-off Ohmic spectral density, as a function of the chain mode number *n* for a fixed *R* = 5 and several temperatures (*α* = 0.12 and *k*<sub>c</sub> = 1).</figcaption></figure>

between the system and the $\hat{c}_n$ and $\hat{d}_n^\dagger$ operators.

The finite temperature coupling constants between the system and
the chain for a hard cut-off Ohmic SD keep broadly the same form
as the zero temperature ones. A coupling profile as a function of
system site separations is displayed in Fig. 4.8. The differences are
that the amplitudes increase with temperature, and the peak value is
no longer centred around the mode *n* = *R*/2*c* but rather *n* = *R*/*c*. For
*β* = 0.5 the amplitude of the coupling is doubled compared to the zero
temperature case. We also note that the tail before the peak presents
more oscillations than the zero-temperature one which is smoother.
The change in the coupling profile as a function of temperature is
shown in Fig. 4.9. For high and moderately high temperatures, the
couplings decrease in amplitude as *β* increases but are still centred
around *n* ≈ *R*/*c*. For high values of *β*, the amplitude stays constant but
the maximum swaps to *n* ≈ *R*/2*c* as we recover the zero temperature
value.
# Appendix

## 4.A Finite-temperature extended bath and thermofields

The thermofields method is a purification procedure that maps a system in interaction with a bath at temperature $T$ to a system in interaction with two baths at $T = 0$ in a vacuum state. It thus enables the use of a wave function representation at finite temperature, instead of a density matrix, at the cost of doubling the number of environmental modes. The standard derivation of thermofields relies on the bijection between a Hilbert space and the linear forms on this Hilbert space (i.e. the bra/ket correspondence) and Bogoliubov transformations. Clearly, there are parallels between this method and the construction of the extended environment presented above in the derivation of the finite-temperature chain mapping. Here, we give an alternative derivation of the thermofields solely using bath correlation functions and present its relationship with the extended bath method.

### Main results of thermofields

We outline the principle results of thermofields, i.e. the characteristics of the two new baths. Detailed derivations of the thermofields approach can be found in Refs. [80, 81]. We start from the following OQS Hamiltonian

$$ \hat{H} = \hat{H}_S + \underbrace{\int_{\mathbb{R}^+} d\omega \, \hbar\omega \hat{a}_\omega^\dagger \hat{a}_\omega}_{\hat{H}_B} + \sum_\alpha \underbrace{\hat{Q}_\alpha \int_{\mathbb{R}^+} d\omega \, (g_\omega \hat{a}_\omega + \text{h.c.})}_{\hat{H}_{\text{int}}} \quad (4.48) $$

where $\hat{Q}_\alpha$ is a system operator coupling to the bath and $\hat{a}_\omega^\dagger$ is the bath creation operator for the mode of energy $\hbar\omega$.

A new set of environmental modes $\{\tilde{a}_\omega\}_{\omega \in \mathbb{R}^+}$ independent of the original one is introduced, and a new bath term is added to the Hamiltonian

$$ \hat{H}_B = \int_{\mathbb{R}^+} d\omega \, (-\omega) \tilde{a}_\omega^\dagger \tilde{a}_\omega . \quad (4.49) $$

Starting from a bath with coupling coefficients $\{g_\omega\}_{\omega \in \mathbb{R}^+}$, the two new baths have coupling strengths

$$ g_\omega^1 = g_\omega \cosh(\theta_\omega), \quad (4.50) $$

$$ g_\omega^2 = g_\omega \sinh(\theta_\omega), \quad (4.51) $$

where $g_\omega^1$ couples to the initial bath modes $\hat{a}_\omega$ and $g_\omega^2$ to the additional bath modes $\tilde{a}_\omega$; with the mixing angle $\theta_\omega$ being defined by

$$ \tanh(\theta_\omega) = \exp\left(-\frac{\beta\omega}{2}\right). \quad (4.52) $$
Alternative derivation from the bath correlation function

From the Hamiltonian in Eq. (7.1) we define the bath correlation function

$$C_{\beta}(\tau) = \text{tr}_B [\hat{B}(\tau)\hat{B}(0)\hat{\rho}_B] = \int_{\mathbb{R}^+} |g_{\omega}|^2 (n_{\beta}(\omega)e^{i\omega\tau} + (n_{\beta}(\omega)+1)e^{-i\omega\tau}), \quad (4.53)$$

where $\hat{B} = \int_{\mathbb{R}^+} d\omega (g_{\omega}\hat{a}_{\omega} + \text{h.c.})$ and $\hat{\rho}_B$ is a Gibbs state at the inverse temperature $\beta$. At zero-temperature, the correlation function becomes

$$C_{\infty}(\tau) = \int_{\mathbb{R}^{+}} |g_{\omega}|^{2} e^{-i\omega\tau}. \qquad (4.54)$$

We can see the correlation function $C_{\beta}(\tau)$ as the sum of two zero T correlation functions with coupling constants $g'_{\omega} = g_{\omega}\sqrt{n_{\beta}(\omega)}$ and $g''_{\omega} = g_{\omega}\sqrt{n_{\beta}(\omega)+1}$. Furthermore, one can notice that

$$\left(\sqrt{n_{\beta}(\omega) + 1}\right)^{2} - \left(\sqrt{n_{\beta}(\omega)}\right)^{2} = 1. \qquad (4.55)$$

Thus, we can define

$$\cosh(\theta_{\omega}) \stackrel{\text{def.}}{=} \sqrt{n_{\beta}(\omega) + 1}, \qquad (4.56)$$

$$\sinh(\theta_{\omega}) \stackrel{\text{def.}}{=} \sqrt{n_{\beta}(\omega)}. \qquad (4.57)$$

Hence, we have

$$\tanh(\theta_{\omega}) = \sqrt{\frac{n_{\beta}}{n_{\beta} + 1}} = \exp\left(-\frac{\beta\omega}{2}\right). \qquad (4.58)$$

For these two terms to appear in the correlation function we could have $\hat{B} = \hat{B}' + \hat{B}''$ with $\langle \hat{B}'\hat{B}'' \rangle_B = \langle \hat{B}''\hat{B}' \rangle_B = 0$, which is possible at $T = 0$ if $\hat{B}'$ and $\hat{B}''$ act on different bath modes. Hence, we introduce

$$\hat{B}' \stackrel{\text{def.}}{=} \int_{\mathbb{R}^+} d\omega (g'_{\omega}\hat{a}'_{\omega} + \text{h.c.}), \qquad (4.59)$$

$$\hat{B}'' \stackrel{\text{def.}}{=} \int_{\mathbb{R}^+} d\omega (g''_{\omega}\hat{a}''_{\omega} + \text{h.c.}) \qquad (4.60)$$

with $[\hat{a}'_{\omega_1}, \hat{a}''_{\omega_2}] = [\hat{a}'_{\omega_1}^{\dagger}, \hat{a}''_{\omega_2}] = 0$ and $[\hat{a}'_{\omega_1}, \hat{a}'_{\omega_2}^{\dagger}] = \delta(\omega_1 - \omega_2)$ (idem for $\hat{a}'$).

The two terms $C_{\beta}(\tau)$ have opposite phases, thus we also require that in the interaction picture these new bath $\hat{B}'$ and $\hat{B}''$ operators acquire an opposite phase

$$U_I^{\dagger}(t) \hat{a}''_{\omega} U_I(t) = \hat{a}''_{\omega} e^{-i\omega t}, \qquad (4.61)$$

$$U_I^{\dagger}(t) \hat{a}'_{\omega} U_I(t) = \hat{a}'_{\omega} e^{i\omega t} \qquad (4.62)$$
where $U_I(t)$ is the interaction picture evolution operator, which implies that the bath Hamiltonian is

$$ \hat{H}_B = \hat{H}'_B + \hat{H}''_B \stackrel{\text{def.}}{=} \int_{\mathbb{R}^+} (-\omega) \hat{a}_\omega^{\dagger} \hat{a}_\omega' d\omega + \int_{\mathbb{R}^+} \omega \hat{a}''_\omega^{\dagger} \hat{a}''_\omega d\omega . \quad (4.63) $$

Because the dynamics of the system is entirely determined (given an initial condition) by the system Hamiltonian $\hat{H}_S$ and the bath correlation function $C_\beta(\tau)$, there is an equivalence between a bosonic bath at temperature $T$ (described by a density matrix) coupled to the system with the coupling strength $g_\omega$, and two baths at $T=0$ (described with a wave-function) coupled to the system with strengths $g'_\omega = g_\omega \sqrt{n_\beta(\omega)}$ and $g''_\omega = g_\omega \sqrt{n_\beta(\omega) + 1}$.

## Relation with the extended bath

The extended bath introduced in the T-TEDOPA [78] method merges the two types of effective zero $T$ modes into a single bath where negative energies are allowed by noticing that

$$ n_{\beta}(-|\omega|) = -(n_{\beta}(|\omega|) + 1). \quad (4.64) $$

This identity can be 'propagated' to the two sets of thermofields annihilation and creation operators by imposing

$$ \hat{a}'_{\omega} = \hat{a}''_{-\omega}. \quad (4.65) $$

This corresponds to merging together the two Hilbert spaces $\mathcal{H}'$ and $\mathcal{H}''$ of the 'primed' and 'seconded' modes respectively into a single Hilbert space $\mathcal{H}_B^{\text{ext}} = \mathcal{H}' \oplus \mathcal{H}''$. The bath Hamiltonian can now be written

$$ \hat{H}_B = \int_{\mathbb{R}^+} (-\omega) \hat{a}''^{\dagger}_{-\omega} \hat{a}''_{-\omega} d\omega + \int_{\mathbb{R}^+} \omega \hat{a}''^{\dagger}_{\omega} \hat{a}''_{\omega} d\omega \quad (4.66) $$

$$ = \int_{\mathbb{R}} \omega \hat{a}''^{\dagger}_{\omega} \hat{a}''_{\omega} d\omega . \quad (4.67) $$

It is then natural to interpret $\hat{a}''^{\dagger}_{-|\omega|}$ as creating a bath excitation of negative energy $-\hbar|\omega|$. Similarly, following the consequences of Eq. (4.64), we can define a single set of temperature-dependent coupling constant for the new bath modes and the system

$$ g''_{\omega} = g_{\omega} \sqrt{n_{\beta}(|\omega|) + 1}. \quad (4.68) $$

These new zero $T$ independent bath modes with positive and negative energies form the extended bath on which we apply the chain mapping.
# 5

## Influence of the ordering of the chains modes

> I think the problem is not to find
> the best or most efficient method
> to proceed to a discovery, but to
> find any method at all.

Richard FEYNMAN,
Nobel Lecture

### 5.1 How to manage multiple environments?

In the previous chapter, we introduced chain-mapping techniques of bosonic environments which enable us to represent the joint state of an OQS as MPS. The current chapter focuses on how the obtained MPS is going to be written in practice when several environments are interacting with the system. This is the case, for instance, with the 'correlated environment' we introduced previously. There, the environment is split in two baths corresponding to different directions of propagation for the environmental excitations. Another example, of high interest for quantum technologies, comes from quantum thermodynamics where one might be interested in studying the heat exchange of a quantum system interacting with a hot and a cold reservoirs. A last example comes from molecular systems where electronic excitations can effectively couple to several different vibrational environments [82, 83]. Having several environments interacting with the system leads to an increase of the required computational resources needed to perform an accurate simulation, because the amount of correlation between the system and the environmental degrees of freedom grows. A few methods exist to reduce the computational cost of such simulations, either based on entanglement renormalisation (i.e. changing the structure of the TNS Ansatz to lower locally the entanglement) [82] or on 'clustering' of environmental modes (i.e. changing the order of modes). The entanglement renormalisation approach relies on expressing the joint quantum state as a tree tensor network (TTN) state and then analysing the entanglement properties of system/environments and environments/environments partitions. From this analysis, a new layer of so called entanglement renormalising (ER) tensors is constructed and placed in between the system and the environments to lower the entanglement but keep correlations between environments. This approach necessitates a thorough analysis of the entanglement structure between

5.1 How to manage multiple environments? 59

5.2 Two environments: Mode arrangements 60
Left-Right arrangement 62
Successive arrangement 62
Interleaved arrangement 63
Comparison 63

5.3 Three environments: Connectivity 68

5.4 Discussion 70

The ER tensors don't have physical legs and are there to 'compress' correlations. They are to be related to the disentanglers that can be found in another TNS the Multiscale Entanglement Renormalisation Ansatz (MERA) [84].
the system and the environments and between the environments them-
selves before being able to write the TTN state.

On the other hand, the methods relying solely on the reordering of the
environmental modes preserve a simple TTN state or MPS structure.
How environmental modes are ordered in the MPS representation of
the quantum state is a priori arbitrary. One just has to make sure that
the operators in the MPOs applied to the state follow the same order-
ing (which might make their expression more complex). Nevertheless,
it does not mean that all orderings of the environmental modes are
equivalent. The computational cost associated with MPSs is related to
their bond dimension (as introduced in chapter 2). As explained be-
fore, these bond dimensions are related to the amount of correlations,
or the degree of entanglement, between bi-partitions of the quantum
state. It is known that reordering the modes in some specific cases of
single bosonic or fermionic environments leads to reducing the degree
of entanglement, and thus to an improved efficiency of the computa-
tions. In the case of bosonic modes, the On the Fly Swapping (OFS)
is a procedure applied during the sweeps of the (real or imaginary)
time evolution where neighbouring tensor sites are swapped if it re-
sults in decreasing a given cost function (entanglement entropy or the
truncation error) [85]. Even though this method is a practical way of
improving computational performance, it needs to be applied at each
time step of the state evolution. In the case of fermionic bath modes,
the reordering is particularly simple and consists of alternating bath
modes corresponding to filled orbitals with modes corresponding to
empty ones [86, 87]. Whereas the reordering procedure for bosonic de-
grees of freedom is a dynamic procedure, this fermionic reordering is
static and done once and for all when writing down the initial state.

A natural question then emerges: is there a similar efficiency-increasing
static reordering of bath modes in the bosonic case? We study here the
influence of different bosonic modes arrangements on the efficiency of
open quantum systems simulations. We initially conjecture that entan-
glement in MPS states can be lowered by diminishing the ‘correlation
length’ between two environmental excitations created at the same
time. Investigating the convergence of different arrangements of baths
modes we are able to give a negative answer to this hypothesis. This
exploration brings us to study another property of tensor network
states that we show to be related to entanglement: the connectivity of
the system with its environments. We bring evidence that adopting a
geometry where there is a single ‘interface’ between the system and the
several environments it connects reduces the amount of entanglement
in the state.

## 5.2 Two environments: Mode arrangements

In several cases the quantum system of interest is interacting with
two baths. A first instance comes from the study of thermodynamic
properties and dynamics of quantum systems. For example, when
<figure><img src="image_42.png" /></figure>

**Figure 5.1:** Schematic drawing of the unitary transformation $U$ – the chain mapping – changing two baths of independent bosonic modes into two chains with nearest neighbour interactions. The correlation length $\xi(t)$ measures the distance in between two dynamically correlated environmental excitations.

studying heat flows when a quantum system is between a cold bath and a hot bath or the system thermalization behaviour [72, 88]. This category of situations are important in the context of quantum thermal machines. Another example is the two chains of our correlated environment for spatially extended systems. More generally, the two bosonic baths could be of different nature, for example one representing a vibrational environment and the other an electromagnetic environment [89]. There, the two environments could play a different role, and their effects can simply add up, or new features that don't come from one environment alone could emerge (such behaviour is called non-additive). In any cases, the problem of a quantum system interacting with two bosonic environments is more than just a case study and has physical relevance.

When the two environments have similar characteristic timescales, energy dissipation from the system results in the generation of environmental excitations that are then propagating in their respective chains in a dynamically correlated way. It thus appears ‘natural’ to define a ‘correlation length’ $\xi$ between these excitations, namely the number of MPS’s sites in between them. The hypothesis that directly follows is that the degree of entanglement, as described by the bond dimension, grows with this correlation length. Though not formulated in this way, this hypothesis forms a main conclusion of the work previously performed on modes reordering in the context of fermionic baths [87].

Hence, in this section we introduce three different orderings of bath
**Figure 5.2:** Diagram of the different mode arrangements in the MPS for a system interacting with two bosonic baths. The dotted lines represent the long-range couplings between the system and chain modes. (a) 'Left-Right' arrangement of the chain-mapped environment. Chain modes are symmetrically placed on both sides of the system. The correlation length between the excitations of the two chains grows linearly in time. (b) 'Successive' arrangement of the chain-mapped environment. One chain connects directly to the system and the second one is appended to the first one. The correlation length between the excitations of the two chains is a constant of the order of the length of the chain. (c) 'Interleaved' arrangement of the chain-mapped environment. Chain modes are alternating and thus coupling to their next nearest neighbours. The correlation length is constant of the order of 1.

The Left-Right MPS is 'isomorphic' to a TTN state with two branches: one for each environment.

<figure><img src="image_43.png" /></figure>

modes (depicted in Fig. 5.2) that display three different behaviours for the correlation length $\xi$, and we compare the minimum bond dimension they require in order to reach convergence in 1TDVP simulations.

## Left-Right arrangement

Usually, the MPS and MPO representation of the joint system-environment state and Hamiltonian follows what we call a 'Left-Right' ordering where the system is placed in between the two chains. This arrangement is 'intuitive' because it reproduces the structure of the Hamiltonian and, thus, often used when describing a quantum system interacting with two independent baths [72, 86]. However, this ordering of the environmental modes might lead to a non-optimal maximal bond dimension. Indeed, the initial environmental excitations created by the environment propagate towards the end of their respective chains but are dynamically correlated. Hence, as time passes, these excitations that are highly correlated move further apart thus a priori requiring a higher bond dimension for the MPS as the correlation length grows linearly with time $\xi \propto t$. Figure 5.1 shows a representation of this correlation length.

## Successive arrangement

An alternative ordering of the chain modes corresponds to placing the system at one end of a chain and the other chain at the other end (i.e. concatenate the two chains), as shown in Fig. 5.2 (b). We coin this ordering of the bath modes 'Successive'. This might seem counterintuitive as the initial bath excitations thus induce immediately a correlation over a chain-long region. However, the correlation length stays approximately constant $\xi \approx N_m$ as the excitations propagate along the chain. Moreover, the arrangement of the chain modes enables
a small optimization of the size of the system's tensor in the MPO representation of the Hamiltonian.

## Interleaved arrangement

To diminish the bond dimension a solution would be to bring closer the modes that are highly correlated with one another. In the Successive arrangement case, compared to the Left-Right one, this is done by having a fixed distance of $\xi \approx N_m$ between the correlated excitations. An improvement could be achieved by interleaving the two chains together into a single one. This is done by alternating the mode of one chain with the corresponding mode of the other one, as shown in Fig. 5.2 (c). Thus, correlated excitations are now separated by a fixed distance of $\xi \approx 1$ thus diminishing the correlation length. This type of mode ordering has already been applied to fermionic problems, namely Anderson impurity, where one bath describes empty orbitals and the other one filled orbitals.

## Comparison

To compare these three different mode arrangements, we are first going to apply them to a well known and studied OQS model, the Independent Boson Model (IBM), and later on to the Spin Boson Model (SBM).

**Definition 5.2.1** The Spin Boson Model is a paradigmatic model of OQS where a single two-level system interacts linearly with a bosonic environment

$$ \hat{H} = \frac{\epsilon}{2} \hat{\sigma}_z + \frac{\Delta}{2} \hat{\sigma}_x + \int_0^{+\infty} \omega \hat{a}_\omega^\dagger \hat{a}_\omega d\omega + \frac{\hat{\sigma}_z}{2} \int_0^{+\infty} (g_\omega \hat{a}_\omega + \text{h.c.}) d\omega, \quad (5.1) $$

with $\epsilon$ the TLS energy gap and $\Delta$ its tunnelling term. When the tunnelling is null the SBM reduces to the so called **Independent Boson Model** where the system and interaction Hamiltonians commute. Contrary to the SBM, the IBM is analytically solvable.

In the IBM, because $[\hat{H}_S, \hat{H}_{\text{int}}] = 0$, the system's populations are invariant. The system only experiences pure dephasing of its initial coherences and is not subjected to decay.

## Independent Boson Model

We first show that all arrangements are able to accurately describe system's dynamics at zero and finite temperature by considering an analytically solvable model. The IBM describes a spin interacting with a bosonic environment where the system coupling operator commutes with the system's Hamiltonian $\hat{H}_S$. The Hamiltonian considered here is rotated by $\pi/2$ compared to Def. 5.2.1, i.e the coupling operator is
Because the coupling constants $g_k$ are the same for both baths, this model could also be written as a single bath IBM with a fourfold Kondo parameter $\alpha$.

$\hat{\sigma}_x$. The Hamiltonian we consider has two bosonic baths

$$ \hat{H} = \frac{\Delta}{2} \hat{\sigma}_x + \sum_{i=1}^{2} \int_{\mathbb{R}} \omega_k \hat{a}_k^{i\dagger} \hat{a}_k^i dk + \frac{\hat{\sigma}_x}{2} \sum_{i=1}^{2} \int_{\mathbb{R}} (g_k \hat{a}_k^i + \text{h.c.}) dk, \quad (5.2) $$

where $\hat{a}_k^{i\dagger}$ creates an excitation in bath $i$ with energy $\omega_k$, and $g_k$ is the coupling strength between bath mode $k$ and the system for both baths. The expectation value $\langle \sigma_z \rangle$ is given by [9, 35]

$$ \langle \sigma_z \rangle (t) = \cos(\Delta t) \text{Re} \left[ \rho_{\uparrow\downarrow}(0) \exp \left( - \int_0^\infty \frac{J(\omega)}{\omega^2} (1 - \cos(\omega t)) \coth \left( \frac{\beta\omega}{2} \right) \right) \right], \quad (5.3) $$

where $\rho_{\uparrow\downarrow}(0)$ is the initial coherence of the system, $J(\omega) = \sum_k |g_k|^2 \delta(\omega - \omega_k) = 2\alpha\omega H(\omega_c - \omega)$ is the bath spectral density and $\beta = (k_B T)^{-1}$ is the bath inverse temperature. Figures 5.3 (a) and 5.3 (b) show that the Left-Right, Successive and Interleaved arrangements all recover the analytical expression for the smallest non-trivial bond dimension $D=2$ for widely different values for $\Delta$ and $\alpha$. This small bond dimension is expected as the analytical solution of state of the system (for the initial state given below) is known to be an entangled state between the system and a displaced environment

$$ |\psi(t)\rangle = \frac{e^{-i\theta(t)}}{\sqrt{2}} \left( |\uparrow_x\rangle \otimes \hat{D}(\{\alpha_k(t)\}) |\{0\}_k\rangle - |\downarrow_x\rangle \otimes \hat{D}(\{-\alpha_k(t)\}) |\{0\}_k\rangle \right) \quad (5.4) $$

where $\theta(t)$ and $\alpha_k(t)$ depend solely on the SD $J(\omega)$ and the energies of bath modes $\omega_k$ [35] and $\hat{D}(\{\alpha_k\})$ is a multi-mode displacement operator. Given that the environment state is simply a displaced state, manifesting no entanglement of its own, the information about the joint state can be stored with a single bit of information, hence the value of the bond dimension $D=2$. Similarly, Fig. 5.3 (c) shows that the analytical behaviour is also recovered at finite temperature for the first non-trivial bond dimension. The initial joint state is a product state

$$ \left\{ \begin{array}{ll} |\downarrow_z\rangle \otimes |\{0\}_k\rangle & \text{if } \beta = \infty \\ |\downarrow_z\rangle \langle\downarrow_z| \otimes \exp(-\beta\hat{H}_B)/Z & \text{if } \beta \neq \infty \end{array} \right. \quad (5.5) $$

where $|\downarrow_z\rangle$ is the eigenstate of $\hat{\sigma}_z$ associated with the eigenvalue $-1$, and $Z$ is the bath partition function.

## Correlated Environment

Now that we have established that our three different arrangements are able to accurately describe OQS dynamics, we are going to study the convergence behaviour of these different orderings with a non-trivial model: a correlated environment model. We consider a 1D system composed of two sites labelled by $\gamma$ placed in space at the position $r_\gamma$. These two sites are interacting with a common bosonic bath of plane waves labelled by the wave-vector $k \in \mathbb{R}$. The corresponding
<figure><img src="image_44.png" /><figcaption>**Figure 5.3:** Dynamics of the $\langle \sigma_z \rangle$ for the Interleaved, Successive and Left-Right chain arrangements for several maximal bond dimensions $D$. (a) The model parameters are $\Delta = 0.8\omega_c$, $\beta = \infty$, $\alpha = 0.2$. (b) The model parameters are $\Delta = 0.2\omega_c$, $\beta = \infty$, $\alpha = 0.1$. (c) The model parameters are $\Delta = 0.4\omega_c$, $\omega_c\beta = 10$, $\alpha = 0.1$. The three different arrangements are able to recover the analytical results at zero and finite temperature for the first non trivial bond dimension $D = 2$.</figcaption></figure>

Hamiltonian is

$$ \hat{H} = \sum_{\gamma} E_{\gamma} \hat{f}_{\gamma}^{\dagger} \hat{f}_{\gamma} + \omega_0 (\hat{f}_{\gamma} \hat{f}_{\gamma+1}^{\dagger} + \text{h.c.}) + \int_{\mathbb{R}} \omega_k \hat{a}_k^{\dagger} \hat{a}_k dk \\ + \sum_{\gamma} \hat{f}_{\gamma}^{\dagger} \hat{f}_{\gamma} \int_{\mathbb{R}} (g_k e^{ikr_{\gamma}} \hat{a}_k + \text{h.c.}) dk, \qquad (5.6) $$

where $\hat{f}_{\gamma}^{\dagger}$ creates an excitation on site $\gamma$, $\hat{a}_k^{\dagger}$ creates a plane wave of energy $\omega_k = |k|c$ with $c$ the phonon speed (and $\hbar = 1$), and $g_k e^{ikr_{\gamma}}$ is a coupling strength between bath mode $k$ and the site $\gamma$. As introduced in chapter 4, the bosonic environment is mapped to two chains, one for the propagating bath modes and the other one for the counter-propagating modes, with long-ranged couplings with the system. Figure 5.1 shows a schematic diagram of the chain mapping procedure that transforms a continuous environment into two interacting chains. We comment on the fact that this is a particular instance of a system coupled to two bosonic baths (they happen to be described by the same parameters up to complex conjugation). We consider this specific type of environment because it is general and can be easily extended to the case of two bosonic baths with different parameters. Furthermore, in that case, the dynamics of the two environments are highly correlated. The results shown hereafter are for degenerate sites $E_{\gamma} = 0$ at positions $k_c r_1 = 0$ and $k_c r_2 = 5$ for different values of coupling strengths to the bath $\alpha$ (the SD is again taken to be Ohmic) and coherent coupling $\omega_0$. Initially the system state is localised on the first site

$$ \left\{ \begin{array}{ll} |r_1\rangle \otimes |\{0\}_k\rangle & \text{if } \beta = \infty \\ |r_1\rangle \langle r_1| \otimes \exp(-\beta \hat{H}_B)/Z & \text{if } \beta \neq \infty \end{array} \right. \qquad (5.7) $$

where $|r_1\rangle$ is the state corresponding to site 1 being occupied, i.e. $\langle r_1 | \hat{f}_1^\dagger \hat{f}_1 | r_1 \rangle = 1$. Figure 5.4 (a) shows the population of the second site at zero temperature for a strong system-bath coupling $\alpha = 0.1$ and a weak tunnelling energy $\omega_0 = 0.1\omega_c$. The Left-Right arrangement converges at $D = 5$ whereas successive and interleaved are already converged at $D = 4$. One can notice that for the first non-trivial bond dimension $D = 2$, the Interleaved ordering is slightly better than the Successive one. Increasing the system-bath coupling and the tunnelling
<figure><img src="image_45.png" /><figcaption>**Figure 5.4:** Population of the second site in the correlated environment model. Generally the Left-Right arrangement converges for higher bond dimensions than the Interleaved and Successive arrangements that have a similar convergence behaviour. The amount of correlations in the joint state is increasing throughout the three parameter regimes (a), (b) and (c). In the legends, Int., Succ. and L.-R. respectively stand for an Interleaved, a Successive, and a Left-Right chain arrangement. The model parameters are (a) $E_Y = 0$, $\alpha = 0.1$, $\omega_0 = 0.1\omega_c$, $\beta = \infty$, (b) $E_Y = 0$, $\alpha = 0.2$, $\omega_0 = 0.4\omega_c$, $\beta = \infty$, and (c) $E_Y = 0$, $\alpha = 0.1$, $\omega_0 = 0.1\omega_c$, $\beta = 1$.</figcaption></figure>

energy requires a larger bond dimension to describe accurately the dynamics of the system, as shown in Fig. 5.4 (b). The system-bath coupling has been doubled $\alpha = 0.2$ and the tunnelling energy multiplied by four $\omega_0 = 0.4\omega_c$ compared to the previous case. The different mode arrangements are all converged for $D = 10$ but only Interleaved one is represented for readability. We can see that Left-Right is generally performing worse than the Successive and Interleaved arrangements. For unconverged values of $D$, it exhibits dynamics that are further away from the converged one than the Successive and Interleaved arrangements. These two arrangements look similar but interleaved seems to catch the coherent oscillations better for $D = 6$. In order to require a higher bond dimension for the state, we now consider a finite temperature $\omega_c\beta = 1$ case. The couplings are the same as in Fig. 5.4 (a). We do not show the converged results in Fig. 5.4 (c) (which are for $D = 15$) as we are only interested in the differences between the different arrangements, and want to prevent the graph to be too cluttered. We can see in Fig. 5.4 (c) that the Left-Right ordering is again less accurate than the other two for smaller bond dimensions. It exhibits a dynamic similar to the other orderings only when reaching $D = 10$. The Interleaved arrangement performs better than the Successive one for $D = 5$, but both become similar for larger values.

From these examples, we can rule out the explanation in terms of correlation length because the Interleaved and Successive arrangements have very similar convergence properties despite having an order of magnitude difference in their correlation lengths. Nevertheless, these two orderings are still converging faster than the usual Left-Right one. This result might be surprising as these arrangements don't reproduce the structure of the underlying Hamiltonian. This leads us to consider another hypothesis related to the connectivity of the system-environment couplings. In the Left-Right case the system is connected twice to the environment, doubling – so to speak – its surface of exchange with it, whereas the Successive and Interleaved ones are coupled only once to a composite environment. A second piece of evidence that the correlation length is not a valid explanation for the efficiency of the computation is given by a fourth mode arrangement which has a constant correlation length and two points
<figure><img src="image_46.png" /></figure>

**Figure 5.5:** 'Reverse Left-Right' arrangement of the chain-mapped environment. Chains modes are symmetrically placed on both side of the system but on one chain they are in the reverse order. The correlation length between the excitations of the two chains remains constant and is the same as in the 'Successive' ordering.

**Figure 5.6:** Comparison of the Left-Right and Reverse Left-Right ordering of chains modes for the same parameters as in Fig. 5.4 (b) ($E_γ = 0, α = 0.2, ω_0 = 0.4ω_c, β = ∞$). Even though the correlation length $ξ$ of correlated environmental excitations are different by an order of magnitude, the convergence behaviour of the dynamics with respect to $D$ is the same.

of connection – 'interfaces' – between the system and the environment. Keeping a Left-Right arrangement of the mode but taking one of the chain in the reverse order (see Fig. 5.5) gives a situation where the correlation length is fixed (and similar to the Successive arrangement) as the excitations are now propagating in the same direction. Hence, this new Reverse Left-Right arrangement has the same topology as the Left-Right ordering but a different correlation length. Figure 5.6 shows the second site population for the same choice of parameters as in Fig. 5.4 (b). The two versions of the Left-Right arrangement have a similar convergence behaviour despite having radically different correlation lengths. Hence, we can rule out the hypothesis that the important metric for the MPS bond dimension is the correlation length between excitations in both environments.

Now that the hypothesis of the connection between the correlation length and the bond dimension has been refuted, it is important to notice that the intuitive relation between the dynamics of correlated excitations and the growth of bond dimension – often invoked to justify the feasibility of a given TNS representation – is hand-wavy. As convincing as it may sound initially, there is no firm ground on which this assertion is set. Indeed, one could argue that dynamically correlated excitations carry very low entanglement in the MPS as the only
information they share is their distance from the system which is updated by the action of the time evolution operator. Therefore, to better understand the absence of connection between the correlation length $\xi$ and the growth of the bond dimension, a rigorous mathematical analysis for simple local Hamiltonian should be performed.

## 5.3 Three environments: Connectivity

To test the hypothesis of the importance of connectivity on the growth of bond dimensions, we study a spin interacting with three identical bosonic baths. We name these baths ‘up’, ‘left’ and ‘right’. The Hamiltonian of this three baths Spin Boson Model (S3BM) is

$$ \hat{H} = \frac{\epsilon}{2}\hat{\sigma}_z + \frac{\Delta}{2}\hat{\sigma}_x + \sum_{i=1}^{3} \hat{\sigma}_x \int_{0}^{\infty} d\omega \sqrt{J(\omega)} (\hat{a}_{\omega}^{i} + \hat{a}_{\omega}^{\dagger i}), \quad (5.8) $$

where $\hat{\sigma}_j$ are Pauli matrices, $J(\omega)$ is the bath spectral density and $\hat{a}_{\omega}^{i}$ creates a bosonic excitation of energy $\omega$ in the $i^{th}$ bath. We keep an Ohmic spectral density with a hard cut-off $J(\omega) = 2\alpha\omega_c H(\omega_c - \omega)$.

Three different configurations are looked at to test the hypothesis of the role of connectivity in the growth of maximal bond dimension: (1) the tree configuration where the system is connected directly to three baths; (2) the left-right configuration where the system is connected to the left and (right + up) chains; and (3) the successive configuration where the system is connected to a single chain (left + right + up). In order to study in real time the evolution of the required bond dimension necessary for the convergence of the MPS state, we time-evolve the state using a 2-site variant of the TDVP method (2TDVP) and the adaptive variant of the 1-site TDVP (DTDVP) [72]. These two variants are able to update in real time the bond dimensions of the MPS state. Figure 5.7 (a) shows the evolution of the bond dimension of the MPS with a Left-(Successive Right+Up) arrangement during a time evolution performed with 2TDVP [60]. The chains all have $N_m = 30$ modes. We can see that the bond dimension grows quickly and reaches the cut-off bond dimension of the simulation ($D_{Max} = 50$) at $\omega_c t \approx 6$. The other noticeable element is that the third chain which was appended to the right chain was not updated by 2TDVP. This is consistent with the known fact that 2TDVP has difficulties taking into account long-ranged interaction. Hence, this case has also been studied with an Interleaved arrangement of the (right + up) chain to make it tractable with 2TDVP as shown in Fig. 5.7 (b). The main conclusion on the growth of the bond dimensions is unchanged. These heat maps directly show that the growth of the bond dimensions is localised around the system. The DTDVP algorithm which enables a one-site update of the MPS and dynamically evolving bond dimensions was also tried on this arrangement but always got stuck in the initial manifold.

However, the DTDVP method was used to study the one chain (Left + Right + Up) case shown in Fig. 5.7 (c). This chain was in a Successive

Again, this simple model can be mapped to a SBM with a single environment, and is used here as a first investigation of our hypothesis.
<figure><img src="image_47.png" /><figcaption>**Figure 5.7:** Evolution for a S3BM at $\beta = \infty$ of the bond dimensions of the MPS with: (a) a Left-S-(Right-Up) configuration obtained with 2TDVP. The cut-off value for the bond dimensions $D_{Max} = 50$ is quickly reached and 2TDVP is not able to handle the long range coupling induced by appending the right and up chains. (b) a Left-S-(Right interleaved with Up) configuration obtained with 2TDVP. The cut-off value for the bond dimensions $D_{Max} = 50$ is quickly reached. (c) a S-(Left+Right+Up) configuration obtained with DTDVP. The bond dimensions evolve more slowly.</figcaption></figure>

arrangement and the long range coupling were handled correctly by the method. Here, the bond dimensions grow at a slower pace and reach maximal allowed bond dimension $D_{Max} = 50$ at $\omega_c t \approx 25$, which corresponds to a fourfold improvement on the case Left-(Right+Up) geometry that has a connectivity of 2. The growth of the bond dimension with this arrangement is not localised around the system per se but at the beginning of each chain.

The 2TDVP and DTDVP algorithms, as currently implemented in our MPSDynamics package [76], are not yet straightforward to apply to TTN. In that case we went back to the 1TDVP method and chose the relative error of the expectation value $\langle \sigma_z \rangle$ with respect to the $D = 50$ results to be our metric for the growth of the bond dimension as we are not interested in the dynamics of this observable in itself. Figure 5.8 shows the relative error for $D \in \{10, 20, 30\}$ for a selection of parameters spanning a variety of dynamics. In most cases a residual error still persists in the $D = 30$ case whose value is on the order of a percent. Most importantly, these errors consistently manifest at early times $\omega_c t \lesssim 5$ which does not contradict our hypothesis and seems to point in the same direction as the results obtained for the Left-Right and Successive arrangements. However, in order to formulate a firm conclusion we would need to have access to the dynamics of the bond dimension for the TTN as we did in Fig. 5.7. To do so, we are currently developing an implementation of the DTDVP method to TTN.

With these examples we can see that connectivity seems to play an important role in limiting the growth of the bond dimensions needed to accurately describe the state in time.
<figure><img src="image_48.png" /><figcaption>**Figure 5.8:** Relative error between the dynamics of $\langle \sigma_z \rangle$ for different values of $D$ and $D = 50$ obtained with 1TDVP for different simulation parameters with a TTN state. In every case small errors of the order of a percent are always present for $D = 30$ at early time. The simulation parameters are (a) $\alpha = 0.1$, $\epsilon = 0.2\omega_c$ and $\Delta = 0.2\omega_c$, (b) $\alpha = 0.5$, $\epsilon = 0.2\omega_c$ and $\Delta = 0.2\omega_c$, and (c) $\alpha = 0.5$, $\epsilon = 0.2\omega_c$ and $\Delta = 0.5\omega_c$.</figcaption></figure>

## 5.4 Discussion

In this chapter, we have shown that, when a system is interacting with several environments, different arrangements of the environmental modes – and crucially, not only the ‘intuitive’ TTN structures – can lead to well-converged results. The hand-wavy argument often given against arrangements where independent environments are not coupled directly to the system, namely that the growth of bond dimension in the multi-environments MPS states is related to the correlation length $\xi(t)$ of dynamically correlated environmental excitations, has been disproved for bosonic environments. We came to this conclusion by

(1) categorising (and introducing new) environmental mode arrangements in the chain representation,

(2) showing that they are all able to accurately simulate the OQS dynamics,

(3) comparing their convergence behaviour with respect to the joint state bond dimension,

(4) therefore showing that their convergence behaviours are independent of the correlation length $\xi$.

This first result has a consequence of practical importance as it implies that joint system-environments state can always be written as MPS which are slightly easier to implement than TTN.

Showing that the usual Left-Right ordering of bath modes is less efficient than the less ‘natural’ Successive and Interleaved arrangements lead us to consider an alternative hypothesis to the correlation length one. For bosonic environments, the number of bonds connecting the system to the environments (‘interfaces’) is a quantity associated with the growth of the bond dimension. The lower the number of interfaces, the lower the required bond dimension. For simple models that can be mapped easily to single environment problems, our results agree with this hypothesis. Nevertheless, investigations must continue to test this new conjecture in more complex settings and – if it holds true – gain a deeper understanding of its origin.
These results have a compounding impact on the efficiency of TN
simulations as the number of elements in the system's tensor scales
proportionally to $D^N$ with $N$ the number of interfaces between the sys-
tem and its environments. For instance, in the S3BM the scaling is pro-
portional to $D^3$ for TTN states, $D^2$ for Left-Right MPS, and only $D$ for
Successive and Interleaved MPS. The results presented here thus im-
ply that simple TTN states can always be recast as Successive-ordered
MPS, thus a priori removing the need for entanglement renormalisa-
tion of environmental couplings. These results add to the diverse body
of knowledge around the importance of selecting the appropriate ge-
ometry when using TNS and highlight that the appropriate geometries
are not always the ones reproducing the structure of the interaction
Hamiltonian.

The following chapters will focus on applying the tools introduced
and developed so far to study non-Markovian effects in OQS.
<figure><img src="image_49.png" /></figure>
APPLICATIONS
<figure><img src="image_50.png" /></figure>
*Someone I say, will remember us
in another time*

Sappho, FRAGMENT 147

Building on the concepts and tools we introduced and developed in the previous parts, in this chapter, we formulate a fully quantum mechanical model that allows us to explore strong spatio-temporal correlations. The existence of such correlations necessitates a manifestly non-Markovian description of the dynamics. The majority of the original results we are presenting in this chapter have been published in Ref. [1]. This model opens a route to establishing the phenomenology of non-Markovian dissipation in the regime where all the relevant timescales – system dynamics, relaxation transitions and environmental signalling – are similar. Our objective is to identify and understand the underlying microscopic phenomena behind the non-Markovian phenomena emerging from our model, in order to create a conceptual toolkit that could be used to exploit these effects, including any explicitly non-classical effects, in artificial nanoscale devices. The motivation for this exploration arises from protein-based ‘nanomachines’ selected through natural selection to perform the key optoelectronic tasks of photosynthesis as presented in chapter 1. For example, the pigment-protein complexes (PPC) that perform the electron transfers at the core of photosynthesis are composed of photoactive pigments in interaction with a highly structured environment made of a protein scaffold that tunes the electronic and vibrational properties of the molecular network. Coordinating multiple charge dynamics in structures with poor dielectric screening and typical lateral sizes of only 5–6 nm requires exquisite spatio-temporal control of energy transfer and electron transport, including mechanisms of feedback to ensure the processes occur in the correct order without waste of excited state energies. While the role of the structured environments found PPCs has been widely discussed in terms of transport efficiency and the possible support of coherent electronic dynamics in light-harvesting [17, 19, 90, 91], the signalling and potential efficiency gains from spatio-temporal feedback and heralding feedforward processes in the environment has received rather scant attention. However, first principles methods based on crystal structures do show that the large secondary protein elements that span the electron transport chain in the Reaction Center (RC) could ‘communicate’ the initial and final sites of the electron trans-

<table><tr><td>6.1 Model</td><td>76</td></tr><tr><td>Hamiltonian MPO formula-<br>tion</td><td>78</td></tr><tr><td>Long- and short-distance lim-<br>its</td><td>79</td></tr><tr><td>Mapping to SBM for N = 2</td><td>80</td></tr></table>

<table><tr><td>6.2 Non-Markovian recurrences<br>and information backflow</td><td>81</td></tr><tr><td>Zero temperature</td><td>81</td></tr><tr><td>Finite temperature</td><td>86</td></tr></table>

<table><tr><td>6.3 Discussion</td><td>88</td></tr></table>
**Figure 6.1:** Biological inspiration for our correlated bath model. (Left) The protein structure of a nanoscale photosynthetic reaction centre. Photoactive pigments are held rigidly by non-covalent protein interactions that also tune their electronic overlaps, interactions and excited state energies. The coordination of multiple cofactors by extended structures, such as quasi-1d alpha helices, allows vibrational fluctuations to act on different cofactors in a spatio-temporally correlated manner. (Right) Structure of the cofactors active in charge separation through quantum electron transport (ET). The oxidation of water in photosynthesis requires four successful electron transports, and this multi-fermion process is regulated through feed-forward (FF) and feedback (FB) mechanisms induced by strong electron-hole interactions with the dissipative protein scaffold.

<figure><img src="image_51.png" /></figure>

port, and may act to prevent accumulation of further charges [92]. Elsewhere in biology, the idea of dynamical structural changes as a way to regulate processes is well established, especially in the field of allosteric regulation [93, 94]. Considered as an open quantum system problem, the existence of strong spatio-temporal correlations necessitates a manifestly non-Markovian description of the dynamics, as the key physics is encoded in the retarded 'action at a distance' that results from previous system-bath interactions, energy exchange, etc.

## 6.1 Model

We consider a one-dimensional chain of $N$ sites labelled by $\{\alpha\}_{\alpha=1...N}$ in a common one-dimensional bosonic bath with modes characterised by the wave-vectors $k \in [-k_c, +k_c]$, where $k_c$ is the environment cut-off wave-vector. The environment dispersion relation is chosen as the usual linear dispersion $\omega_k = |k|c$ with $c$ the speed of the phonons in the bath. We restrict ourselves to the single excitation subspace of the system described by a Hamiltonian $\hat{H}_S$ with nearest neighbour hopping. The total Hamiltonian has the form

$$ \hat{H} = \hat{H}_S + \hat{H}_E + \hat{H}_{\text{int}} \qquad (6.1) $$

$$ \begin{aligned} = & \sum_{\alpha=1}^{N} E_{\alpha} \hat{f}_{\alpha}^{\dagger} \hat{f}_{\alpha} + \sum_{\alpha=1}^{N-1} J (\hat{f}_{\alpha}^{\dagger} \hat{f}_{\alpha+1} + \text{h.c.}) \\ & + \int_{-k_c}^{+k_c} \omega_k \hat{a}_k^{\dagger} \hat{a}_k dk + \sum_{\alpha=1}^{N} \hat{f}_{\alpha}^{\dagger} \hat{f}_{\alpha} \int_{-k_c}^{+k_c} (g_k^{\alpha} \hat{a}_k + \text{h.c.}) dk \end{aligned} \qquad (6.2) $$

where $\hat{f}_{\alpha}$ annihilates a system excitation on site $\alpha$, $\hat{a}_k$ is the annihilation operator of a bath mode of wave-vector $k$, $g_k^{\alpha} = g_k e^{ikr_{\alpha}}$, with $g_k = g_{-k} \in \mathbb{R}$, is the coupling strengths between the system and the bath and $r_{\alpha}$ is the position of the site $\alpha$. A schematic of the model is presented in Fig. 6.2.

Because we restrict ourselves to the system's single excitation subspace, we don't have to deal with the fermionic nature of the system excitation.
<figure><img src="image_52.png" /><figcaption>**Figure 6.2:** Schematic diagram of the model under study. A system composed of interacting sites is embedded into a single bosonic environment. Each site couples differently to the environment.</figcaption></figure>

In the rest of the chapter, we consider a system made of two degenerate sites with an initial state where the system and its environment are decoupled and the bath is empty

$$|\psi(t=0)\rangle = |S(0)\rangle \bigotimes_{k \in [-k_c, k_c]} |0_k\rangle = |S(0)\rangle \bigotimes_{n \in \mathbb{Z}} |0_n\rangle, \quad (6.3)$$

where $|S(0)\rangle$ is the initial state of the system, $|0_k\rangle$ represents the vacuum state of the mode $k$ of the bath, and $|0_n\rangle$ the vacuum state of the mode $n$ of the chains (as presented in chapter 4). The restriction to two degenerate sites is made for simplicity and ease of results interpretation. As we have already shown in Sec. 4.3, the empty bath can also be used to effectively describe a Gibbs state at temperature $T$, and so we are also presenting results at non-zero temperature. Even though the method here presented works for any initial state of the system, the initial state of the system is chosen to be the highest energy eigenstate (i.e. upper eigenstate) of the system Hamiltonian $\hat{H}_S$.

The resulting {System + Chains} joint wave function is represented as a MPS in a successive ordering (as presented in chapter 5) where the first couple tensors correspond to the system sites and the following ones to the modes of the two chains

Here we introduce a convention where one chain is labelled with positive integers and the other one with negative integers.

$$|\psi(t)\rangle = \underset{\substack{\text{S1} \\ 2}}{\bigcirc} - \underset{\substack{\text{S2} \\ 2}}{\bigcirc} - \underset{\substack{\text{n=1} \\ d}}{\bigcirc} - \underset{\substack{\text{2} \\ d}}{\bigcirc} - \cdots - \underset{\substack{-1 \\ d}}{\bigcirc} - \underset{\substack{D \\ d}}{\bigcirc} - \underset{\substack{-2 \\ d}}{\bigcirc} - \cdots \quad (6.4)$$

where $d < \infty$ is the truncated Hilbert space dimension of the chain modes, i.e. the maximal number of excitations of a mode. The number of modes of the two semi-infinite chains are truncated at large enough values $N_m$, such that an excitation on the chain does not have the possibility to reach the end of the chain during the time evolution.
Hamiltonian MPO formulation

Following the structure of the MPS presented above, we want to write
the corresponding MPO representation of the Hamiltonian in Eq. (6.2)
in the chain picture using the method presented in chapter 3. The
on-site tensor has a bond dimension D = 2(α + 2) for the α-th site and a
physical dimension d_S = 2

$$
W_1 = (\hat{1} \quad J \hat{f}_1 \quad J \hat{f}_1^\dagger \quad \hat{f}_1^\dagger \hat{f}_1 \quad \hat{f}_1^\dagger \hat{f}_1 \quad E_1 \hat{f}_1^\dagger \hat{f}_1), \quad (6.5)
$$

and (in the general case of *N* sites for completeness)

$$
W_{1<\alpha\le N} =
\begin{pmatrix}
\hat{1} & J\hat{f}_\alpha & J\hat{f}_\alpha^\dagger & \underbrace{0 \ 0 \ \dots}_{2(\alpha-2)} & \hat{f}_\alpha^\dagger\hat{f}_\alpha & \hat{f}_\alpha^\dagger\hat{f}_\alpha & E_\alpha\hat{f}_\alpha^\dagger\hat{f}_\alpha \\
0 & 0 & \cdots & \cdots & \hat{f}_\alpha^\dagger & \hat{f}_\alpha & \hat{0} \\
0 & 0 & \cdots & \cdots & \hat{f}_\alpha & \hat{f}_\alpha & \hat{0} \\
\hat{1} & \hat{1} & \cdots & \cdots & \hat{0} & \hat{0} & \hat{1}
\end{pmatrix}
\tag{6.6}
$$

with $J \to 0$ for the last system tensor, and empty elements correspond-
ing to zeros. The first chain on-site tensor has a similar structure, but
with a constant bond dimension $D = 2(N+2)$ for each mode, and a
physical dimension $d < \infty$

$$
W_{0 \le n \le N_m - 1} =
\begin{pmatrix}
\hat{1} & t_n \hat{c}_n^\dagger & t_n \hat{c}_n & 0 & 0 & \dots & 0 & \omega_n \hat{c}_n^\dagger \hat{c}_n \\
& 0 & & & & & & \hat{c}_n \\
& 0 & & & & & & \hat{c}_n^\dagger \\
& \hat{1} & & & & & & \gamma_n^1 \hat{c}_n \\
& & \hat{1} & & & & & \gamma_n^{1*} \hat{c}_n^\dagger \\
& \dots & & & & & & \vdots \\
& & & \hat{1} & & & & \gamma_n^{N*} \hat{c}_n^\dagger \\
& & & & \hat{1} & & & \hat{1}
\end{pmatrix},
\quad (6.7)
$$

with $t_{N_m-1} = 0$. The second chain tensors are identical with $\hat{d}_n$ and $\gamma_n(r)\hat{d}_n^\dagger$ instead of $\hat{c}_n$ and $\gamma_n(r)\hat{c}_n$. The last tensor is

$$
W_{-(N_m-1)} = \begin{pmatrix}
\omega_{N_m-1} \hat{d}_{N_m-1}^{\dagger} \hat{d}_{N_m-1} \\
\hat{d}_{N_m-1} \\
\hat{d}_{N_m-1}^{\dagger} \\
\gamma_{N_m-1}^{*} \hat{d}_{N_m-1}^{\dagger} \\
\gamma_{N_m-1}^{1*} \hat{d}_{N_m-1} \\
\vdots \\
\gamma_{N_m-1}^{N*} \hat{d}_{N_m-1} \\
\hat{1}
\end{pmatrix} . \quad (6.8)
$$

One might notice that the chains sites tensors have a bond dimension
*D* that is fixed by the number of sites in the system *N*. This means that
having a large environment only increases the number of individual tensors one needs but not their size. We would like to stress that this result is central for the tractability of this approach. We have introduce identity operators on the diagonals of the $A_k$ blocks to carry out along the chain the long range coupling coefficients such that they are associated with the corresponding system site operator $\hat{f}_a^\dagger \hat{f}_a$. Hence, they allow a local representation of the Hamiltonian as a MPO even though the interactions are long range across the chains.

To illustrate how the Hamiltonian is recovered from these tensors, we perform the calculation in the case where there is only one site in the system and two modes on a unique chain. In that case there are only three tensors:

$$W_1 = (\hat{1} \quad \hat{f}_1^\dagger \hat{f}_1 \quad \hat{f}_1^\dagger \hat{f}_1 \quad E_1 \hat{f}_1^\dagger \hat{f}_1), \qquad (6.9)$$

$$W_2 = \begin{pmatrix} \hat{1} & t_1 \hat{c}_1^\dagger & t_1 \hat{c}_1 & 0 & 0 & \omega_1 \hat{c}_1^\dagger \hat{c}_1 \\ 0 & 0 & 0 & \hat{1} & 0 & \gamma_1^1 \hat{c}_1 \\ 0 & 0 & 0 & 0 & \hat{1} & \gamma_1^{1*} \hat{c}_1^\dagger \\ 0 & 0 & 0 & 0 & 0 & \hat{1} \end{pmatrix}, \qquad (6.10)$$

$$W_3 = \begin{pmatrix} \omega_2 \hat{c}_2^\dagger \hat{c}_2 \\ \hat{c}_2 \\ \hat{c}_2^\dagger \\ \gamma_2^1 \hat{c}_2 \\ \gamma_2^{1*} \hat{c}_2^\dagger \\ \hat{1} \end{pmatrix}. \qquad (6.11)$$

The contraction of $W_2$ and $W_3$ gives a $5 \times 1$ tensor – the same shape as $W_3$ and the transpose of the shape of $W_1$

$$W_2 \cdot W_3 = \begin{pmatrix} \omega_2 \hat{c}_2^\dagger \hat{c}_2 + t_1 (\hat{c}_1^\dagger \hat{c}_2 + \hat{c}_1 \hat{c}_2^\dagger) + \omega_1 \hat{c}_1^\dagger \hat{c}_1 \\ \gamma_2^1 \hat{c}_2 + \gamma_1^1 \hat{c}_1 \\ \gamma_2^{1*} \hat{c}_2^\dagger + \gamma_1^{1*} \hat{c}_1^\dagger \\ \hat{1} \end{pmatrix}. \qquad (6.12)$$

The $\gamma_2^1$ and $\gamma_2^{1*}$ couplings have been ‘captured’ by the identity operators and added to the $\gamma_1^1$ and $\gamma_1^{1*}$ terms. Further contraction with $W_1$ gives a ‘scalar’ corresponding to the example Hamiltonian

$$W_1 \cdot W_2 \cdot W_3 = \omega_2 \hat{c}_2^\dagger \hat{c}_2 + t_1 (\hat{c}_1^\dagger \hat{c}_2 + \hat{c}_1 \hat{c}_2^\dagger) + \omega_1 \hat{c}_1^\dagger \hat{c}_1 \\ + \gamma_2^1 \hat{f}_1^\dagger \hat{f}_1 \hat{c}_2 + \gamma_1^1 \hat{f}_1^\dagger \hat{f}_1 \hat{c}_1 + \gamma_2^1 \hat{f}_1^\dagger \hat{f}_1 \hat{c}_2 \\ + \gamma_1^{1*} \hat{f}_1^\dagger \hat{f}_1 \hat{c}_1^\dagger + \gamma_2^{1*} \hat{f}_1^\dagger \hat{f}_1 \hat{c}_2^\dagger \\ + E_1 \hat{f}_1^\dagger \hat{f}_1 \qquad (6.13)$$

$$W_1 \cdot W_2 \cdot W_3 = \hat{H}. \qquad (6.14)$$

## Long- and short-distance limits

Because of the dependence of the system-bath coupling strengths on the spatial configuration of the system, the system-chain couplings are long-ranged; thus the system can create excitations on different regions
of the chain. We consider an Ohmic SD as presented in chapter 4. As we have seen in Fig. 4.7, the amplitude of the coupling before the peak decreases with the position of the peak. Said differently, the larger distance between the two sites, the less the second site interacts with the beginning of the chain. Thus, we can expect that for infinite separation when $R \to \infty$ this system will solely interact with the first mode of each chain, and hence behave like a SBM. The Hamiltonian in Eq. (6.2) can be mapped onto the SBM Hamiltonian in Eq. (5.1) with $\epsilon = 0, \Delta = 2J$, and $\hat{f}_{\alpha}^{\dagger}\hat{f}_{\alpha} = (\hat{\sigma}_{z} + \hat{1})/2$, for large $R$. Figure 6.3 shows the comparison between the SBM and the infinite separation case. As expected, in this limit we recover the well-studied SBM. Looking at

<figure><img src="image_53.png" /><figcaption>**Figure 6.3:** Dynamics of the up-state $|\uparrow_z\rangle$ of a Spin Boson Model (SBM) compared with the dynamics of the upper eigenstate of the Correlated Environment model for corresponding parameters ($c = 1, J = 0.25\omega_c$ and $\alpha = 0.2$) with a large separation $k_c R = 200$ between the two sites of the system. The two dynamics are the same.</figcaption></figure>

the opposite limit, when the separation between the two system sites vanishes, Eq. (6.2) tells us that the system completely decouples from the environment because in the single excitation subspace $\sum_{\alpha} \hat{f}_{\alpha}^{\dagger}\hat{f}_{\alpha} = \hat{1}$.

### Mapping to SBM for $N = 2$

The $N = 2$ case is a specific case where, because of the symmetry around the midpoint between the two sites, the problem presented here can be written in the form of a SBM with an effective spectral density that depends explicitly on the sites separation. Consider the interaction Hamiltonian $\hat{H}_{\text{int}}$ in Eq. (6.2) in the case of a two-site system with intersite distance $R$, we have

$$ \hat{H}_{\text{int}} = \sum_{\alpha=1}^{2} \hat{f}_{\alpha}^{\dagger} \hat{f}_{\alpha} \int_{-k_c}^{+k_c} (g_k^{\alpha} \hat{a}_k + g_k^{\alpha*} \hat{a}_k^{\dagger}) dk \quad (6.15) $$

$$ = \sum_{\alpha=1}^{2} \hat{f}_{\alpha}^{\dagger} \hat{f}_{\alpha} \int_{0}^{+k_c} \left( g_k^{\alpha}(\hat{a}_k + \hat{a}_{-k}^{\dagger}) + g_k^{\alpha*}(\hat{a}_k^{\dagger} + \hat{a}_{-k}) \right) dk. \quad (6.16) $$
We can introduce a new set of vibrational modes, the symmetric mode $\hat{c}_k$ and the antisymmetric mode $\hat{d}_k$

$$ \hat{c}_k = \frac{\hat{a}_k + \hat{a}_{-k}}{\sqrt{2}}, \qquad (6.17) $$

$$ \hat{d}_k = \frac{\hat{a}_k - \hat{a}_{-k}}{\sqrt{2}}. \qquad (6.18) $$

Hence, the interaction Hamiltonian becomes

$$ \hat{H}_{\text{int}} = \sum_{\alpha=1}^{2} \hat{f}_{\alpha}^{\dagger} \hat{f}_{\alpha} \int_{0}^{+k_c} \left[ g_k^{\alpha} \sqrt{2}(\hat{c}_k + \hat{c}_k^{\dagger}) + g_k^{\alpha*} \sqrt{2}(\hat{c}_k + \hat{c}_k^{\dagger}) \right. \\ \left. + g_k^{\alpha} \sqrt{2}(\hat{d}_k - \hat{d}_k^{\dagger}) - g_k^{\alpha*} \sqrt{2}(\hat{d}_k - \hat{d}_k^{\dagger}) \right] dk \qquad (6.19) $$

We choose the origin of position at the midpoint between the two sites so that (in the single excitation subspace)

$$ \hat{H}_{\text{int}} = (\hat{f}_1^\dagger \hat{f}_1 + \hat{f}_2^\dagger \hat{f}_2) \int_0^{+k_c} 2\sqrt{2}g_k \cos\left(\frac{kR}{2}\right) (\hat{c}_k + \hat{c}_k^\dagger) dk \\ + (\hat{f}_1^\dagger \hat{f}_1 - \hat{f}_2^\dagger \hat{f}_2) \int_0^{+k_c} 2\sqrt{2}ig_k \sin\left(\frac{kR}{2}\right) (\hat{d}_k^\dagger - \hat{d}_k) dk \qquad (6.20) $$

$$ \hat{H}_{\text{int}} = \hat{1}_{S} \text{ const} + \hat{\sigma}_{z} \int_{0}^{+k_c} 2\sqrt{2}ig_k \sin\left(\frac{kR}{2}\right) (\hat{d}_k^{\dagger} - \hat{d}_k) dk \qquad (6.21) $$

Therefore the system only couples to the antisymmetric vibration modes and thus corresponds to a SBM with an effective spectral density $J_{\text{eff}}(k) = 8|g_k|^2 \sin^2(\frac{kR}{2})$. However for larger values of $N$ it is no longer possible to map the system to a SBM. This is similar to the spin-mapping presented in [46].

Now, we study the system and bath dynamics for intermediate values of the separation between the two sites. In this distance range we expect that the new phenomenology of non-Markovian effects associated with spatial correlations between the system and the environment will be manifest.

## 6.2 Non-Markovian recurrences and information backflow

### Zero temperature

At zero temperature, the dynamics of the two-level system in a bosonic environment is well known and described by the SBM [35]. In the system's eigen-basis, the population of the upper state (high energy state) should spontaneously decay to the lower state on a time-scale given by the intensity of the coupling between the system and the bath (as illustrated in Fig. 6.3). The right panel of Figure 6.4 shows the evolution of the eigen-populations with an initial state of the system being the upper eigenstate. We clearly see that the upper level population
<figure><img src="image_54.png" /><figcaption>**Figure 6.4:** System and bath dynamics. (Left) A heatmap of the chain occupation in time showing the propagation of bath excitations along the chains. (Right) Upper eigenstate population. An eigenstate revival and a site localisation are associated with a chain excitation reaching the beginning of the chain. The separation between the two sites is $k_c R = 40$, their coupling is $J = 0.25 \omega_c$, the speed of sound is $c = 1$, $\alpha = 0.12$ and $k_c = 1$.</figcaption></figure>

decays as expected until $t \approx R/c$ when a revival happens. This revival corresponds to an increased localisation of the excitation on the second site of the system after following an evolution in a spatial superposition. With the same conditions, the SBM exhibits the same dynamics except for the revival. The study of the bath in the chain representation allows us to have a spatial interpretation of the interaction between the system and its environment as the maximum coupling between a system's site and the chain is localised around $n = k_c R/2$. The left part of Fig. 6.4 shows a heatmap of the occupation of the modes of the chains as a function of time. The positive and negative chain modes each correspond to one of the two chains necessary to take into account propagating and counter-propagating $k$-modes. The corresponding initial system state is an excitation delocalised on the two sites with a separation $k_c R = 40$.

We can see that the chain modes around $n = \pm k_c R/2 = \pm 20$ get populated first and that the corresponding bath's excitations then propagate on the chains. At $\omega_c t \approx 20$ an excitation propagating from the mode $n = 0$, coupled mostly to the first site, and an excitation propagating from the mode $n = 20$ constructively interfere around $n = 10$. The former continues to propagate on the chain and traces a ray in the diagram. The latter reaches $n = 0$ at $\omega_c t \approx 40$ and is reflected. We can see from this diagram that revivals happen when the excitation emitted along the chain by one site reaches the part of the chain interacting with the other site. We thus have a feedback effect of the environment on the system. This transfer of bath's excitations from one site to the
<figure><img src="image_55.png" /><figcaption>**Figure 6.5:** (Left) A heatmap of the chain occupation in time showing the propagation of bath excitations along the chains. Arrows have been added to represent the trajectories of chains' excitations. (Right) System eigenstate population for an initial state in the upper eigenstate. The separation between the two sites is $k_c R = 20$, their coupling is $J = 0.25\omega_c$, the speed of sound is $c = 2$ and $\alpha = 0.12$. We can definitely see a revival of population at a time consistent with the amount of time needed for a bosonic excitation to travel into the bath from one system's site to the other.</figcaption></figure>

other is a direct manifestation of the information backflow we mentioned in chapter 2 as a source of non-Markovian effects.

The dynamics of the chain with negative modes is not the reflection of the dynamics of the chain with positive modes. Indeed the negative chain modes correspond to the propagating *k*-modes, hence the excitations created by the second site move away from the origin of the chain (which is coupled to the first site). On the contrary, bath's excitations created by the second site on the positive modes chain correspond to the counter-propagating *k*-modes and move toward the origin of the chain. On both chains the excitations created by the first system site propagate toward the end of the chain as they move away in real space from the first site. This explains the apparent 'asymmetry' between the two chains.

The dynamics of the system, all other parameters being the same, only depends on the ratio $R/c$. This is also true for the chain dynamics, for example the ($R = 40, c = 1$) and ($R = 20, c = 0.5$) cases have the same heatmaps. This was expected as the system's sites couple in both cases to the same parts of the chain and the bath's excitations travel on the chain at the same speed. Increasing the propagation speed of the bath excitations we can generate several revivals with something like an echo between the two sites, as shown in Fig. 6.5 where revivals with decreasing amplitudes can be observed with a periodicity of $R/c$. All the parameters are the same as in Fig. 6.4 except the speed of the bath's excitations that has been doubled. The left panel of Fig. 6.5