Keywords: physique quantique, systèmes quantiques ouverts, équation maîtresse de Lindblad, réseaux de tenseurs, circuits quantiques bruités, apprentissage automatique, états quantiques neuronaux, mécanisme de Kibble-Zurek, dynamique quantum physics, open quantum systems, Lindblad master equation, tensor networks, quantum circuits, machine learning, neural quantum states, Kibble-Zurek mechanism

Cette thèse est consacrée à l'étude de la dynamique hors-équilibre dans les systèmes quantiques fermés et ouverts. La question de l'obtention d'équations dynamiques numériquement efficaces pour simuler de tels systèmes, qui souffrent d'une malédiction de dimensionnalité lorsqu'ils sont de grande taille est au centre de la thèse, ainsi que les divers nouveaux résultats physiques qui peuvent être obtenus grâce à ces équations efficaces. La question de l'obtention de représentations efficaces des états quantiques est également centrale, à travers la construction d'Ansatz.

En particulier, nous étudions théoriquement la propagation de l'intrication dans les systèmes photoniques dissipatifs à forte interaction, motivés par la réalisation expérimentale récente de telles plateformes. Nos résultats montrent qu'en dépit du caractère dissipatif du système étudié, la propagation de l'intrication est ballistique avec une vitesse en accord avec l'image de quasiparticules impliquées dans la dynamique, des doublons et des holons photoniques. Notre analyse révèle que la dissipation est fortement asymétrique dans les deux cas étudiés de l'injection d'un doublon ou de son extraction, avec une influence plus forte dans le cas de l'extraction.

Ensuite, nous présentons la méthode du sous-espace dynamique, une nouvelle méthode numérique pour la simulation de l'évolution temporelle de systèmes quantiques ouverts avec une entropie modérée. En représentant la matrice densité par un sous-espace constitué de M états purs, la complexité de l'intégration de l'équation maîtresse de Lindblad est considérablement réduite, et l'algorithme obtenu est entièrement déterministe. Nous appliquons ensuite cette méthode pour modéliser des algorithmes quantiques bruités, en simulant un algorithme quantique bruité de transformation de Fourier qui est au centre de l'algorithme de Shor pour la factorisation de nombre premiers. Nous trouvons une loi d'échelle pour l'erreur en fonction du taux de dissipation, et le comportement de l'infidélité par rapport au circuit idéal est modélisé en fonction des états initiaux qui sont fournis en entrée à l'algorithme.

Enfin, l'application des méthodes d'apprentissage automatique à la mécanique quantique est présentée, à travers la classe d'Ansatz des états quantiques neuronaux (NQS). La construction générale des NQS est introduite, ainsi que sa relation avec d'autres Ansatz tels que les réseaux tensoriels. Ensuite, le problème de la simulation de la dynamique des systèmes fermés avec NQS est étudié, ainsi que les problèmes courants liés à la géométrie de l'espace des paramètres qui peuvent survenir. Un schéma alternatif est introduit et appliqué à des dynamiques de trempes dépendantes du temps dans des systèmes critiques, où le mécanisme de Kibble-Zurek peut être étudié. Enfin, différentes constructions d'Ansatz pour les systèmes ouverts sont présentées, notamment la combinaison de la méthode du sous-espace dynamique avec un Ansatz neuronal.

its mathematical structure, as has been the case with tensor network methods [START_REF] Orús | A practical introduction to tensor networks: matrix product states and projected entangled pair states[END_REF]. The development of such methods has led to the entanglement area law [START_REF] Eisert | Colloquium: area laws for the entanglement entropy[END_REF] for local unidimensional gapped systems, which has since then proved invaluable. Secondly, they enable researchers to study new problems, as they are more efficient to treat different regimes that were previously unattainable. A significant part of the thesis is dedicated to new physical results that have been obtained thanks to these effective descriptions. Such models for quantum mechanical systems have for now mostly been applied to equilibrium and ground state properties. Density functional theory [START_REF] Jones | Density functional theory: its origins, rise to prominence, and future[END_REF], a theoretically straightforward effective model useful for complex molecules, is perhaps the best example of the success of effective models for equilibrium properties, as it is nowadays ubiquitously used in research institutions and the private sector. Dynamical properties, as we will see in more detail in Chapter 1, are harder to pinpoint and involve a more complex mathematical structure than ground state properties in the general case. In particular, open quantum systems, a class of systems that describe quantum systems coupled to an external environment, hence being nonequilibrium almost by definition, are particularly elusive in that very few general results are known about them. Recently, general properties related to dissipative phase transitions [START_REF] Minganti | Spectral theory of Liouvillians for dissipative phase transitions[END_REF] have been put forward, but the exposition of dynamical properties for open systems remain rare for now.

In Chapter 1, we present the theoretical framework surrounding the study of nonequilibrium quantum systems. We start by presenting the problems that arise for closed quantum systems, where nonequilibrium dynamics involve superpositions of eigenstates of a many-body Hamiltonian, thus resulting in a highly complex evolution. We present results related to correlation propagation in such systems [START_REF] Cheneau | Experimental tests of Lieb-Robinson bounds[END_REF], as well as the Kibble-Zurek mechanism, that has been experimentally verified in the last decade [START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF][START_REF] Keesling | Quantum Kibble-Zurek mechanism and critical dynamics on a programmable rydberg simulator[END_REF]. We also present various numerical methods, from mean-field methods to tensor networks, that have yielded many important results in previous years. Then, we present the framework of open quantum systems. When a quantum system is coupled to an external environment, one can generally not compute the dynamics of the full composite system. Therefore, an efficient way to treat this problem is to treat the environment as a bath that is weakly coupled to the system, that is, remains weakly entangled with the system. By also supposing that the timescale related to the dynamics of the bath is much faster than that of the system, one can treat the system with a Lindblad master equation approach [START_REF] Breuer | The theory of open quantum systems[END_REF]. This equation closely resembles the Schrödinger equation, and has additional terms accounting for the coupling of the system with the environment. We present various numerical methods that were developed for such systems, as well as superconducting circuits [START_REF] Bishop | Circuit quantum electrodynamics[END_REF], being a prominent platform for the study of such systems. In Chapter 2, entanglement and its evolution in a strongly-correlated open quantum systems is studied. As mentioned, the study of correlation propagation as for now mostly been concentrated on closed systems [START_REF] Cheneau | Experimental tests of Lieb-Robinson bounds[END_REF][START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF], and the impact of the interplay of strong interactions and dissipation on entanglement spreading has remained relatively untouched [START_REF] Alba | Spreading of correlations in Markovian open quantum systems[END_REF]. We investigate the nonequilibrium dynamics in a photonic Mott insulator induced by adding or removing a particle from the center of the chain. As such, we observe a ballistic entanglement propagation, with a speed that is predicted by an approximate solution for the closed system in the strongly interacting regime. We also observe a drastic change in the behaviour of the system as a function of whether we start by removing or adding a particle, related to the form of the states involved in the dynamics. These results pave the way to more general results about entanglement dynamics in open systems.

In Chapter 3, the dynamical corner-space method is introduced. This method is based on the idea that there exists a class of systems whose von Neumann entropy is low enough that the density matrix describing the state of the open quantum system may be described as a collection of a limited number of pure states, thus rendering the integration of the Lindblad master equation more efficient. This idea is at the core of another method developed a few years ago, the corner-space renormalization method. The application of this method to the quantum Fourier transform will be presented. This is a key sub-routine in Shor's algorithm, and when subject to dissipation and/or decoherence, undergoes open system dynamics where a sequence of quantum gates are applied, entangling the qubits that form the register of the algorithm. Thanks to this method, new scaling behaviour for this algorithm as a function of dissipation and decoherence rates for this was possible to unveil. The impact of the initial state to the computation is also investigated, and a highly biased performance is found as a function of the loss channels. Chapter 4 is devoted to using machine learning techniques to study nonequilibrium dynamics. Having shown immense success to represent high-dimensional data, such as images, text and audio, the latest advances in machine learning are currently strongly impacting the quantum physics community, where such advances have allowed to make progress on a large number of problems [START_REF] Carleo | Machine learning and the physical sciences[END_REF][START_REF] Dawid | Modern applications of machine learning in quantum sciences[END_REF], from ground state search, to quantum tomography or the simulation of quantum circuits. Such methods will be presented in detail, and we will then focus on their use for the simulation of nonequilibrium dynamics. This endeavor has been plagued by several problems, to which we propose a solution by reformulating the problem at hand onto an optimization problem to be solved at each time step of the dynamics. We will also see how machine learning methods can be used to study open quantum systems dynamics, and that representing a density matrix with neural networks presents several problems. As such, we envision an extension of the dynamical corner-space method to the use of neural network states, which we describe towards the end of the chapter.

Nonequilibrium quantum systems I Introduction

As mentioned in the introduction, it quickly became clear that the Schrödinger equation could not be solved for more than a few particles in the presence of interactions [START_REF] Dirac | Quantum mechanics of many-electron systems[END_REF]. This means that one cannot know the exact dynamics of most physical systems. However, this does not mean that many features of these dynamics cannot be accessed. Luckily for us, one can, through various approximation methods and educated guesses on the form of the wavefunction, referred to as ansätze in this manuscript, make great progress towards this goal. In this section, we present the two main problems that are studied throughout the thesis. The first problem is understanding the nonequilibrium dynamics of ideal quantum systems, as described by the Schrödinger equation. Great progress has been made concerning quantum mechanics at equilibrium, however nonequilibrium properties are only starting to be properly investigated thanks to the rapid improvement of numerical methods and experimental platforms. Indeed, there is a much greater scarcity in the tools that can be used to study such systems. Many long standing problems have not been solved, such as the eigenstate thermalization hypothesis [START_REF] Srednicki | Chaos and quantum thermalization[END_REF], and the validity of the Kibble-Zurek mechanism [START_REF] Kibble | Topology of cosmic domains and strings[END_REF][START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF], to name a few. In ideal quantum systems, nonequilibrium dynamics occur when the initial state of a dynamical process is not an eigenstate of the applied Hamiltonian [START_REF] Mitra | Quantum quench dynamics[END_REF]. This can be done by locally adding a particle to the system, or to apply a quench, that can be either time-independent or time-dependent. In the first part of this chapter, the physical and numerical ramifications of this problem will be presented, as well as the Kibble-Zurek mechanism and the experimental works that have been devoted to studying these issues in recent years. The main numerical methods used to study these issues such as mean-field and tensor networks will be presented in section II. [START_REF] Donatella | Dynamics with autoregressive neural quantum states: application to critical quench dynamics[END_REF]. The second problem is understanding the dynamical properties of open quantum systems. Such systems are inherently out of equilibrium, as they can be viewed as ideal quantum systems coupled to an external environment. As explained in the introduction, studying open quantum systems is crucial as it enables one to understand quantum noise and its effect in quantum computational processes [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF] as well as thermodynamical tasks [START_REF] Carrega | Dissipative dynamics of an open quantum battery[END_REF]. In this context, dynamical properties can be studied by preparing the system in a certain state, and evolving the open quantum system in time. Most questions around this problem aim to understand how results for ideal quantum systems hold when an external environ-ment is present. For example, the spreading of correlations in open quantum systems is an active field of research [START_REF] Alba | Spreading of correlations in Markovian open quantum systems[END_REF], that will be discussed in more detail in Chapter 2. On the quantum information side, understanding how noise processes affect quantum algorithms is of crucial importance if we wish to use them one day [START_REF] Preskill | Quantum Computing in the NISQ era and beyond[END_REF]. The theoretical framework for treating open quantum systems will be presented, as well as various experimental platforms which can be described by this framework. In this context, systems that interact with a physical environment are considered, nontrivial dynamical phenomena will arise, as dissipation will act as a new player.

The many-body problem

Let us now be more specific about the many-body problem, in both closed and open quantum systems. According to the axioms of quantum mechanics, the state of an ideal quantum system is fully specified by its wavefunction |ψ⟩, that belongs to the Hilbert space associated to the system H. |ψ⟩ can be expanded in a given basis of this Hilbert space. A natural choice for this basis is a composite basis of local quantum numbers s i (for spin systems, these local quantum numbers are the local spins), since most quantum systems we wish to describe are composite: they consist of L interacting degrees of freedom. In this part, and in most of the manuscript, we suppose that the equations to be solved to describe the time evolution of a physical system cannot be solved. Therefore, the only way to go forward is to perform a numerical integration of this equation, which involves storing in memory the state of the system at each step, and incrementally calculating the evolution of the system. The basis expansion of |ψ⟩ reads: with c s 1 ,s 2 ,...,s L the amplitudes (complex numbers) corresponding to |s 1 , s 2 , . . . , s L ⟩, the basis elements of the Hilbert space (configurations) of the system. The Born rule tells us that by taking the modulus square of such amplitudes, we obtain the probabilities for finding the quantum state in such configuration after performing a measurement, i.e. p(s 1 , s 2 , . . . , s L ) = |c s 1 ,s 2 ,...,s L | 2 . A combinatorial explosion can be seen directly here, as there are l L amplitudes one must store (either on paper or on a computer) if one wishes to have a complete description of the system, with l the dimension of the local Hilbert space (which corresponds to the number of values the local quantum numbers s i can take, assuming they each can take the same). More concretely, to describe a system of just 50 two-level systems, 2 50 ∼ 10 15 complex numbers have to be described, or stored in memory. This is impossible even with modern day computers.

Regarding an open quantum system, the problem is similar but worse. As will be explained in this chapter, the state of an open quantum system can be fully described by a density matrix ρ, this time belonging to the bounded operator space B(H) associated to H. ρ can also be expanded over projectors, of which there will be l L : ρ =

II Dynamics of closed systems II.1 Sudden quenches

As mentioned, the problem of quench dynamics involves solving the time-dependent Schrödinger equation, starting from a certain state |ψ 0 ⟩. Throughout the manuscript, we fix ℏ = 1. The Schrödinger equation reads:

d |ψ(t)⟩ dt = -i Ĥ |ψ(t)⟩ . (1.3)
If |ψ(t = 0)⟩ = |ψ 0 ⟩ is an eigenstate of Ĥ, then the evolution as trivial: the phase of the state oscillates with frequency ω 0 corresponding to the energy of |ψ 0 ⟩ and the solution is:

|ψ(t)⟩ = exp(-iω 0 t) |ψ 0 ⟩ . (1.4)
However, if the initial state is not an eigenstate of Ĥ, the time evolution is nontrivial, and involves oscillations given by the frequencies of all the eigenstates of Ĥ. To see this, one can expand |ψ(t)⟩ as a superposition over eigenstates of Ĥ, denoted |ϵ i ⟩:

|ψ(t)⟩ = i c i (t) |ϵ i ⟩ .
where the coefficients c i (t) depend on time. Hence the right hand-side of Eq. (1.3) reads:

-i Ĥ |ψ(t)⟩ = -i i c i (t)ϵ i |ϵ i ⟩ = i ċi (t) |ϵ i ⟩ .
It is clear that if one wishes to integrate this equation and find the coefficients c i (t), one must know what the eigenvalues ϵ i of the Hamiltonian are. Computationally speaking, this requires to diagonalize Ĥ and obtain its full spectrum, in contrast to many problems in quantum physics and chemistry that require the ground state of the system and the first excited states only. Thus in general, if one wants to obtain a high-precision integration of Eq. (1.3), one requires to be able to describe all the excited states of the system, which can be of varying complexity. For instance, for local gapped (i.e., outside of the critical region) one-dimensional systems, the complexity of the ground state is significantly weaker than that of excited states. This fact will be explained and described in more detail when discussing tensor network methods, whose limits reside in such considerations. Note that a quench protocol instantly drives the state away from its previous equilibrium, thus cannot be treated perturbatively as is routinely done in the context of linear response theory, in which important results such as the fluctuation-dissipation have been found [START_REF] Kubo | The fluctuation-dissipation theorem[END_REF]. Since this protocol involves strong interactions, few general results have been obtained for it.

Correlation propagation

A set of results that has been obtained for sudden quench protocols in isolated quantum systems concerns correlation propagation. In a seminal work by Lieb and Robinson [START_REF] Lieb | The finite group velocity of quantum spin systems[END_REF], it was shown that in generic quantum systems with finite-range interactions and a finite local Hilbert state there exists a speed limit on information propagation. This is known as the Lieb-Robinson bound, which reads 1 :

||[ Â(t), B]|| ≤ c exp{-a[d(X, Y ) -v|t|]} (1.5)
with  and B two observables with finite supports X and Y (for lattice systems, the sites associated to the Hilbert spaces on which  and B act), and a, c, v positive constants.

As can be seen from the expression inside the exponential, v corresponds to the velocity of correlations2 : the commutator between two operators with distinct support will decay exponentially the further one goes out of the light cone d(X, Y ) -v|t|. From there an analogy with the light cone from special relativity was drawn, and such a limit was experimentally verified for various physical systems [START_REF] Cheneau | Experimental tests of Lieb-Robinson bounds[END_REF][START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF]. An intuitive explanation of this phenomenon was first provided by Calabrese and Cardy [START_REF] Calabrese | Time dependence of correlation functions following a quantum quench[END_REF]. When suddenly quenching some Hamiltonian parameters, excitations, also known as quasiparticles, will be created, since the equiblibrium quantum state now consists of a superposition of excitations in the new eigenbasis of the quenched Hamiltonian, as explained in section II.1. Only quasiparticles in close regions will be initially entangled. From there, quasiparticles spread baillistically in different directions with some group velocity, hence correlations between distant sites will appear as initially entangled quasiparticles reach them, as shown in Fig. 1.1. In this picture the Lieb-Robinson velocity is simply the maximal group velocity of quasiparticles. Such considerations will be explored further in Chapter 2, where the entanglement dynamics of dissipative quantum systems will be investigated, as well as the validity of the quasiparticle picture to explain such dynamics when coupling to an external environment is taken into account.

II.2 Time-dependent quenches and the Kibble-Zurek mechanism

Another equally interesting protocol consists in performing time-dependent quenches on critical systems that exhibit continuous phase transitions. In this context, a quantum state |ψ 0 ⟩ is prepared as the ground state of some time-dependent Hamiltonian Ĥ[λ(t)] at time t = 0, which depends on some control parameter λ(t). Depending on how the control parameter is tuned, various scenarios can take place. If the control parameter is tuned abruptly, we recover the case of sudden quenches. On the other hand, if the control parameter is tuned more slowly than the inverse of the gap of the system, the system remains in its ground state, which is known as adiabatic tuning. For physical systems that do not exhibit phase transitions, the system follows an adiabatic evolution and remains in the ground state of the Hamiltonian Ĥ[λ(t)] at all times (but may acquire a phase). Figure 1.1: Sketch of the quasiparticle picture, in this case for bosons trapped in an optical lattice, described by the Bose-Hubbard model [START_REF] Cheneau | Experimental tests of Lieb-Robinson bounds[END_REF]. By quenching the lattice depth (that is related to the interaction strength of such models), quasiparticle excitations are created and propagate in the lattice, which leads to correlation propagation that obeys the Lieb-Robinson bound. From [START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF]. Figure 1.2: Scheme of the Kibble-Zurek mechanism according to the adiabatic-impulse approximation. At times ± t, the relaxation time of the system t R crosses the time scale of the quench, which defines the limit between the adiabatic region, where the system essentially stays in the ground state of the system at that time, and the impulse regime, where nonequilibrium dynamics take place.

As an example, let us again consider the transverse-field Ising model, which will serve in the rest of the manuscript. The Hamiltonian is given by:

ĤTFI = -J ⟨m,n⟩ σz m σz n -g m σx m , (1.6)
with a nearest-neighbor coupling strength J and a magnetic field strength g. It exhibits a second-order quantum phase transition that separates a ferromagnetic (for J > 0) or antiferromagnetic (J < 0) phase for g < J from a paramagnetic phase for g > J. It faithfully describes many physical systems constituted of coupled two-level systems immersed in a magnetic field, such as neutral atom systems or magnetic nanomaterials. In one-dimension, the model is exactly solvable by using a Jordan-Wigner transformation. However, for higher dimensions, much less is known about the physics of the system. In 2D, the critical point is situated around g/J ≈ 3.044 as given by quantum Monte-Carlo simulations [START_REF] Blöte | Cluster Monte-Carlo simulation of the transverse ising model[END_REF], which are approximate methods. The Kibble-Zurek mechanism [START_REF] Kibble | Topology of cosmic domains and strings[END_REF][START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF] predicts the formation of topological defects quenched systems undergoing a second-order phase transition as the parameters are linearly swept across a critical point. This universal behavior stems from the fact that there exists a time, denoted t, at which the characteristic evolution time of the system (the relaxation time t R , related to the inverse of the energy gap) becomes larger than the characteristic quench time τ q , as schematically shown in Fig. 1.2. Before t, the dynamics remains quasiadiabatic, as the Hamiltonian parameters are tuned slowly with respect to the time scale of the system. After this time, in contrast, genuinely non-equilibrium dynamics takes place as critical slowing down sets in and the system dynamics becomes slower than the Hamiltonian parameters sweep. Ignoring specific details of this dynamics, one can derive a scaling law for the density of created defects3 at the end of the linear quench [START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF]. This scaling law is derived quite simply. Consider that for second-order phase transitions, the correlation length ξ and relaxation time τ scale as:

ξ ∼ g -ν , τ ∼ g -zν (1.7)
where we have taken the magnetic field g to be an external control parameter. The freeze-out time t is a solution of the equation

τ ( t) = g( t) -zν (1.8)
since we suppose these timescales cross at time t. With g(t) = g 0 (1 -t/τ q ), corresponding to a linear quench:

t = (1/τ q ) -zν/(zν+1) (1.9) ξ = (1/τ q ) -ν/(zν+1) (1.10)
with ξ the correlation length of the domains induced by the nonequilibrium dynamics.

The final defect density is inversely related to the this average correlation length:

⟨n d (τ q )⟩ ∼ ξ-d , (1.11)
and the final scaling law therefore reads

⟨n d (τ q )⟩ ∼ τ -dν/(zν+1) q .
(1.12)

The scaling law for the density of defects has a physical explanation in terms of the change in the different timescales during the quench. This is shown in Fig. 1.2. Far away from the critical point, the gap of the system is large, therefore the characteristic timescale of the system (related to the inverse of the energy gap) remains small; almost any parameter tuning will be adiabatic (corresponding to the far left and right of the figure).

As one approaches the critical point, precisely at the freeze-out timet, t R = ϵ(t): the characteristic time scales of the quench and of the system cross. From there, the evolution stops being adiabatic, entering the impulse regime. Correlations can therefore develop within domains of a limited size, related to the freeze-out time, which is time the system stays in this regime. Then, at + t, the evolution resumes into the adiabatic regime, and domains remain frozen.

II.3 Numerical methods

As we have seen in the previous sections, numerical simulations can be an important tool to study the dynamics of nonequilibrium quantum systems, both in sudden and time-dependent quench protocols. In this section we will present a few of the important methods used to study such protocols.

Mean-field approaches

A class of methods that enables one to treat large quantum systems in arbitrary dimensions is mean-field methods, also known as Gutzwiller ansätze [START_REF] Gutzwiller | Effect of correlation on the ferromagnetism of transition metals[END_REF]. These methods have previously mostly been used for electronic systems [START_REF] Hubbard | Electron correlations in narrow energy bands[END_REF][START_REF] Kanamori | Electron correlation and ferromagnetism of transition metals[END_REF][START_REF] Borghi | Self-consistent Gutzwiller study of bcc-Fe: interplay of ferromagnetic order and kinetic energy[END_REF] and variants of it have been used for the study of bosonic gases [START_REF] Jaksch | Cold bosonic atoms in optical lattices[END_REF]. The key idea is to write the wavefunction as a product of local wavefunctions:

|ψ⟩ = i |ϕ i ⟩ . (1.13)
As such, this enables one to separate the Schrödinger equation into L local equations.

One can also view this as a drastic reduction of the number of free parameters required to describe the many-body wavefunction. Indeed, the single-site wavefunctions are described by l parameters, with l the local Hilbert space dimension. They are subject to the orthogonality condition ⟨ϕ i |ϕ j ⟩ = δ ij , with δ ij denoting the Kronecker delta. In the case of spin systems, this results in |ϕ i ⟩ having only two coefficients corresponding to the probability amplitudes of the spin being up or down. This results in lL complex parameters for an L-body system of local Hilbert space dimension l, which is a drastic improvement upon the l L coefficients needed for the exact wavefunction. The coefficients of the many-body wavefunction are given by:

c s 1 ,s 2 ,...,s L = Π s L s i =s 1 a s i (1.14)
with a s i being vectors of dimension l. However, this approach will only work for weakly correlated quantum systems, since non-separable states cannot be described by such approaches.

An extension of this idea exists as cluster mean-field approaches [START_REF] Mcintosh | Multisite mean-field theory for cold bosonic atoms in optical lattices[END_REF]. The key idea is to consider more than one site per factorized wave function. The state is written as:

|ψ⟩ = Nc i |ϕ i ⟩ .
(1.15)

with N c the number of clusters. The size of a cluster is given by L/N c . This ansatz is able to capture states with an entanglement between sites belonging to the same cluster, thus having strictly more capacity to represent entanglement than the ansatz in Eq. (1.13).

Tensor networks and matrix product states

Beyond mean-field ansätze, powerful methods that have proved very effective in the study of one-dimensional systems are based on tensor networks [START_REF] Orús | A practical introduction to tensor networks: matrix product states and projected entangled pair states[END_REF]. The key idea behind such methods is to approximate the coefficients of the many-body wave function defined in Eq. (1.1) as a product of tensors:

c s 1 ,s 2 ,...,s L = i 1 ,i 2 ,...,i L A i 1 ,i 2 s 1 A i 2 ,i 3 s 2 . . . A i L ,i 1 s L . (1.16)
This is known as the matrix product state (MPS) ansatz. The tensors A i,j σ i are of rank 3, with one dimension corresponding to the physical degrees of freedom σ i , and two dimensions given by M 2 coefficients that can introduce correlations between sites. The number of parameters contained in an MPS ansatz is therefore lM 2 , which is no longer exponential in the system size. To see how the MPS ansatz is constructed, let us define the Schmidt decomposition of a quantum state |ψ⟩ belonging to a composite Hilbert space

H = H A ⊗ H B with dim(H A ) = m, dim(H B ) = n, m ≤ n.
|ψ⟩ = m i=1 Λ i |u i ⟩ ⊗ |v i ⟩ (1.17)
with Λ i the Schmidt coefficients, and m the Schmidt rank of the state. Note that for a separable state we have m = 1. Defining the reduced density matrix ρA = tr B |ψ⟩⟨ψ|, it can be shown that the Schmidt coefficients are related to the entanglement entropy

S A = -tr{ρ A log ρA } = - m i=1 |Λ i | 2 log |Λ i | 2 . (1.18)
For a one-dimensional system of L sites with a local Hilbert space dimension l, the coefficients c s 1 ,s 2 ,...,s L can be expressed as a function of the Schmidt coefficients, by successive Schmidt decompositions:

c s 1 ,s 2 ,...,s L = l α 1 =1 l 2 α 2 =1 • • • min(l n ,l L-n ) αn=1 • • • l α L =1 = Γ [1]s 1 1,α 1 Λ [1] α 1 Γ [2]s 2 α 1 ,α 2 Λ [2] α 2 • • • Λ [L-1] α L-1 Γ [L]s L α L ,1 . (1.19)
The Γ [i] matrices correspond to a basis transformation between the full basis and the Schmidt basis, and the Λ [i] vectors correspond to the Schmidt coefficients after tracing out the (L -i)th site. This encoding seems particularly inefficient, as the Γ [L/2] contains by itself l L elements. However, the beauty and central point of MPS ansätze is that the entanglement entropy of local, gapped many-body systems is limited. This is known as the entanglement area law [START_REF] Eisert | Icolloquium/i: area laws for the entanglement entropy[END_REF]. Formally, for such systems,

S A(L) ∼ ∂L ∼ 1D constant (1.20)
with A(L) a subsystem of size L, and ∂L the size of the boundary region between L and the rest of the system. This means that the number m of Schmidt coefficients Λ i is limited, and they in fact decay exponentially [START_REF] Orús | A practical introduction to tensor networks: matrix product states and projected entangled pair states[END_REF]. This enables one to safely truncate the size of the matrices in the previous equation to a certain bond dimension D. This leads to the construction presented in Eq.(1.16), with a number of parameters scaling as lLD 2 . These methods are powerful, as they enable to describe states with correlations beyond that of the mean-field ansatz while keeping the number of parameters polynomial in the system size. However, for more complex geometries, and in particular in 2D, constructing such an ansatz presents various problems, and in addition, even if the entanglement of a 2D system follows an area law, the area ∂L in 2D is not constant anymore, which means the bond dimension D will scale less favorably than in 1D.

II.4 Experimental systems

Sudden quenches

One of the main reasons to study the nonequilibrium dynamics of quantum systems is the recent progress of experimental platforms that enable researchers to perform controlled experiments, thus being able to directly test predictions. These platforms are known as quantum simulators [START_REF] Georgescu | Quantum simulation[END_REF], and consist in implementing a Hamiltonian with a physical analogue to that Hamiltonian. These allow one to tailor a specific Hamiltonian and the degree of control is such that quenches can be performed, particles can be added, and many different operations can be performed to study effects presented in the previous section. This is reminiscent of Feynman's vision of simulating quantum mechanics by using an analogous quantum mechanical system [START_REF] Feynman | Simulating physics with computers[END_REF].

As an example, let us consider the quantum simulation of a spin Hamiltonian, as presented in ref. [START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF]. An important platform that has been perfected in the last decades has been optical lattices, where atoms can be trapped by optical tweezers and one can use the strong nonlinearity of the energy spectrum of atoms to effectively realize two-level systems. The atoms are often put in highly excited states, such as Rydberg states [START_REF] Browaeys | Many-body physics with individually controlled rydberg atoms[END_REF] that correspond to a high principal quantum number n. This enables a larger degree of control and coupling to electromagnetic fields. Such a platform can be faithfully described by the following Hamiltonian, where the two eigenstates of the σz operator correspond to the ground |g⟩ and excited states |e⟩ of the Rydberg atoms:

Ĥ = - i Ωσ x i - i (I i -∆)σ z i - i̸ =j V ij 2 σz i σz j (1.21)
where V ij = C 6 /|r ir j | 6 is an atomic potential (C 6 is a constant given by experimental details), Ω is the Rabi frequency of an applied transverse-field that drives transitions between |e⟩ and |g⟩, and (I i -∆) determines the frequency of the longitudinal field. ∆ is the laser detuning, and I i = j,i̸ =j V ij /2 represents a site-independent detuning. Here V ij is taken to be positive, corresponding to attractive interactions, and in the absence of transverse and longitudinal fields the ground state of this Hamiltonian is an antiferromagnet 4 . Note that by taking a homogeneous spacing between atoms, one can take V ij to be constant, and simply consider V ij ≈ J. Therefore, one obtains a slightly modified antiferromagnetic Ising Hamiltonian, and quench dynamics can be induced by suddenly tuning the transverse and longitudinal field. In Fig. 1.4, results of the correlation functions C(i, r) = 4 ⟨σ z i σz i+r ⟩ -⟨σ z i ⟩⟨σ z i+r ⟩ are shown for a quench from a paramagnetic phase, with ∆ ≫ J ≫ Ω, to an antiferromagnetic phase with J ≫ ∆, Ω. In the figure ,C(i,i) is related to the local magnetization as

C(i, i) = 1 -4⟨σ z i ⟩ 2 .
One can clearly see an inversion of the sign of C(0, 1) and C(1, 0) as a function of the detuning, transitioning from ferromagnetic (∆ < 0) to antiferromagnetic (∆ > 0) correlations. Note that for the upper panels, corresponding to sudden quenches, the correlations will remain short-ranged, thus the results could be compared to data from exact diagionalization for a 4 × 4 lattice, which are a good fit to the experiment. Another One is mostly interested in the antiferromagnetic region (AFM), whose eigenstates can be probed by suddenly quenching the detuning from deep in the paramagnetic phase, with ∆ ≪ J. (b) Scheme of the experimental setup: a 2D array of atoms trapped in an optical lattice whose states |e⟩ (large red spheres) |g⟩ (small blue spheres) can be coupled with an infrared laser. The top part of the scheme represents the high-resolution objective that can resolve individual sites, enabling precise analyses of experimental data. (c): Typical protocol. The system is prepared in a trivial phase, with a small number of defects (left), and is then quenched into an excited configuration (center). Then, increasing the lattice depth causes the atoms in the excited state to be lost (right), hence only the atoms in the ground states remain to simplify detection. Reproduced from [START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF]. c) and (d) correspond to time ΩT = 0.5π, which is fast enough to be considered a sudden quench, while lower panels (e, (f), (g) and (h) correspond to ΩT = 2.97π, where long-range correlations develop. Reproduced from [START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF]. group has performed a similar study the same year [START_REF] Lienhard | Observing the space-and time-dependent growth of correlations in dynamically tuned synthetic ising models with antiferromagnetic interactions[END_REF], and more results on more complex models, for which the phase diagram is unknown, are expected in the coming years.

Time-dependent quenches

Time-dependent quenches can also be implemented in quantum simulators, and in recent years many works have investigated this, using Rydberg atoms [START_REF] Keesling | Quantum Kibble-Zurek mechanism and critical dynamics on a programmable rydberg simulator[END_REF], in a setup close to what was presented previously, trapped ions [START_REF] Cui | Experimentally testing quantum critical dynamics beyond the Kibble-Zurek mechanism[END_REF], Bose-Einstein condensates [START_REF] Damski | Dynamics of a quantum phase transition in a ferromagnetic Bose-Einstein condensate[END_REF], and even digital quantum simulators such as IBM's machines [START_REF] Higuera-Quintero | Experimental validation of the Kibble-Zurek Mechanism on a Digital Quantum Computer[END_REF]. To implement a time-dependent quench on Rydberg atom platforms, the transverse and longitudinal fields have to be tuned linearly in time.

Here we will focus on the implementation of a linear quench to examine the validity of the Kibble-Zurek mechanism on Rybderg simulator, close to what was presented previously. In Fig. 1.5, experimental results are shown for a linear quench performed on a system that undergoes an Ising-like quantum phase transition, breaking Z 2 symmetry. In the two-level basis spanned by the ground state |g⟩ and the excited state |e⟩, the Hamiltonian describing the system is given by:

Ĥ = Ω 2 i σx i -∆ i P |g⟩ i + V ij ⟨i,j⟩ P |g⟩ i P |g⟩ j (1.22)
where Ω is the Rabi coupling frequency between states |g⟩ and |e⟩, ∆ is the laser detuning, V ij is the interaction strength and P |g⟩ i = 1 2 ( 1i + σz i ) is the projector onto state |g⟩ on site i. compared with results from an MPS numerical integration. Reproduced from [START_REF] Keesling | Quantum Kibble-Zurek mechanism and critical dynamics on a programmable rydberg simulator[END_REF].

The system is initialized in a homogeneous state with all atoms in the |g⟩ state, and after the quench, some domains with atoms all in state |e⟩ appear, whose density is predicted by the Kibble-Zurek mechanism. As such, the critical detuning can be estimated, and the Kibble-Zurek mechanism can be verified. In the same work [START_REF] Keesling | Quantum Kibble-Zurek mechanism and critical dynamics on a programmable rydberg simulator[END_REF], authors investigate other phases, as their machine is programmable, and other critical phases of the system can be investigated. Once again, it is really the degree of control of such platforms that make them ideal to study quench dynamics.

III Open quantum systems III.1 General setting

Many other interesting cases of nonequilibrium quantum phenomena can be found within the framework of open quantum systems. To describe this class of systems, let us consider a quantum system S, described by a density matrix ρ. After performing a few experiments, one can quickly notice that S interacts with its environment E: for example, if S describes a photonic system, photons can leak out in the environment. A straightforward way to understand how the leakage of photons takes place would be to describe the composite system S + E, made of the system S and the environment E. The composite system is described by a Hamiltonian ĤS+E that describes the degrees of freedom of both the system and the environment. It can be written in a general form

ĤS+E = Ĥ { Ŝk } + ĤE { Γk } , (1.23)
with Ĥ the system Hamiltonian, ĤE the environment Hamiltonian, { Ŝk } the set of system operators and { Γk } the set of environment (or bath) operators. It is clear that the manybody problem appears once again, only worse, as realistic environments can be comprised of many bosonic and fermionic degrees of freedom that describe the electromagnetic field and the molecules composing the air and the experimental apparatus. Neglecting these last contributions and supposing that the system interacts with only a few bosonic modes is still not enough for us to compute the dynamics of most composite systems, its Hilbert space still being too large.

A way around this is to search for an effective description of ρS and obtain a time-evolution equation for it. Note that if the reduced system interacts with its environment, the wave function of the composite system will become entangled (non-separable):

|Ψ⟩ SE ̸ = |ψ⟩ S ⊗ |χ⟩ E . (1.24)
This means that if one wants to obtain a description of S, by performing a partial trace on |Ψ⟩ SE one will obtain a density matrix ρS , and cannot use the Schrödinger equation anymore. Hence one needs to change the formalism, and see how a time-evolution equation can be obtained for the density matrix when taking into account the interaction between system and environment. Let us first properly define the objects we will need to use to obtain an effective dynamical equation on the system.

III.2 Density matrices

A density matrix ρ is an operator in the bounded operator space B(H) associated to a Hilbert space H and should satisfy the following properties to describe a physical system [START_REF] Breuer | The theory of open quantum systems[END_REF]:

1. Trace one: Trρ = 1, 2. Positive-semidefinite: p i ≥ 0, with p i its eigenvalues, 3. Hermiticity: ρ † = ρ.

It also admits infinitely many decompositions of the form

ρ = i p i |ϕ i ⟩ ⟨ϕ i | , ( 1.25)
with p i its eigenvalues and |ϕ i ⟩ the corresponding eigenstates. The physical interpretation for such a decomposition is that the density matrix represents a collection of states that each have a probability p i . According to this view, the density matrix tells us which states will be measured with which frequency if we perform a large number of experiments. One can also stay agnostic about the physical interpretation of this picture and consider that a density matrix describes a reduced quantum state, with uncertainty coming from its entropy that stems from entanglement with the environment.

III.3 Quantum operations

Let us now see how the a density matrix can change, meaning how it evolves under quantum operations. A quantum operation is a completely-positive, trace-preserving (CPTP) map, that one applies to a density matrix. The Kraus theorem [START_REF]States, effects, and operations fundamental notions of quantum theory[END_REF] states that this map can be decomposed as the application of N 2 so-called Kraus operators:

χ = K[ρ] = i Ki ρ K † i (1.26)
with χ the transformed density matrix. A quantum operation corresponds physically to any transformation that has happened to the physical system, such as time evolution, a measurement or a certain experimental protocol. The Kraus operators must satisfy the condition i K † i Ki = 1 for Eq. (1.26) to be valid. K must also be linear

K[αρ 1 + β ρ2 ] = αK[ρ 1 ] + βK[ρ 2 ] (1.27)
and conserve the Hermiticity of ρ. All these conditions ensure that the transformed quantum state χ conserves the properties of ρ. K is a superoperator, meaning that it can be compactly represented by a N 2 × N 2 matrix acting on operators, much like operators act on vectors. One can already see the inherent complexity to numericaly study open quantum systems, as objects of size N 2 × N 2 must be manipulated.

As an example of a quantum operation, let us consider a single qubit, whose levels are written |↑⟩ and |↓⟩. We assume that the state |↑⟩ has been prepared, and encodes some kind of information. The simplest model to describe energy relaxation is through the amplitude damping channel N p , that is a quantum operation with Kraus operators given by:

K 0 = 1 0 0 √ 1 -p , K 1 = 0 √ p 0 0 . (1.28)
The strength of the channel is given by p. Applying this channel on the state |↑⟩, we have:

N p [|↑⟩⟨↑|] = p 0 0 1 -p = p |↓⟩⟨↓| + (1 -p) |↑⟩⟨↑| (1.29)
which is a statistical mixture of states |↓⟩ and |↑⟩ with probabilities p and 1 -p, respectively. Through this simple example, we also see how entropy has been created through the process, since we go from a pure state to a mixed state. We have kept things very simple here, with only a single qubit, therefore no coupling between qubits, and have supposed no Hamiltonian inducing some time-evolution. In the general case, these ingredients come into play and the dynamics is much more complex.

III.4 The Lindblad master equation

Having described what a general quantum operation looks like, recall that we are interested in obtaining an equation for the time evolution of the quantum system represented by ρ. Specifically, we are interested in finding a general form for the map M defined by

ρ(t + δt) = Mρ(t) (1.30)
for vanishing δt. By expanding up to second-order, one obtains:

ρ(t + δt) = i Mi (t)ρ(t) Mi (t) = ρ(t) + δt∂ t ρ(t) + O(δt 2 ) (1.31)
By assuming that the environment of the system is much larger than the system itself, one can consider time-independent Kraus operators Mi . Indeed, in this case, the environment itself acts as a bath and remains at equilibrium at all times, hence the action of the environment on the system is the same at all times. This is known as the Born approximation. Now, some care must be taken with respect to the time increment δt, by considering the physical timescales at hand. δt must be small with respect to the system timescale ∆t S (that can be taken to be the smallest characteristic evolution time) for the limit to make sense. However, δt must be much larger than ∆t E , the environment timescale, so that no coherent effects remain in the system-environment interaction. This happens when the spectrum of the environment is flat, and the environment operators are then delta-correlated:

⟨ Γk (t) Γk (t ′ )⟩ = δ(t -t ′ ). (1.32)
This is known as the Markov approximation. We can now take the appropriate limit δt → 0, and obtain a general form for the Kraus operators Mi . We take one of them to be of the order of unity, and others of the order of √ δt in order to obtain a second-order update on the density matrix:

   M0 = 1 -iδt K Mi>0 = √ δt Ĵi . (1.33)
One can split the operator K into a Hermitian part and an anti-Hermitian part, given respectively by Î = ( K + K † )/2 and Ĝ = i( K -K)/2, with K = Î -i Ĝ. One has

M † 0 ρ(t) M0 = ρ(t) -iδt[ Î, ρ(t)] -δt{ Ĝ, ρ(t)} + O(δt 2 ). (1.34)
Using the normalization condition on the Kraus operators, we have:

1 = i M † i Mi = M † 0 M0 + δt i>0 Ĵ † i Ĵi + O(δt 2 ) = 1 -2δt Ĝ + δt i>0 Ĵ † i Ĵi + O(δt 2 ). (1.35)
By rearranging terms, we obtain Ĝ = i>0 Ĵ † i Ĵi /2. The dynamics of the density matrix ρ(t) is therefore given by an equation of the form

ρ(t + δt) = ρ(t) -iδt[ Î, ρ(t)] + δt i>0 Ĵi ρ(t) Ĵ † i - 1 2 { Ĵ † i Ĵi , ρ(t)} + O(δt 2 ). (1.36)
Recall that the evolution of a closed quantum system is given by the von Neumann equation:

ρ(t + δt) = ρ(t) -iδt Ĥ, ρ + O(δt 2 ) (1.37)
we identify the Hamiltonian in Eq. (1.36), and we finally obtain the Lindblad master equation:

∂ t ρ = -i Ĥ, ρ + 1 2 i>0 2 Ĵi ρ Ĵ † i -Ĵ † i Ĵi , ρ (1.38)
where we have dropped the time dependence for readability. By comparing the von Neumann equation and the Lindblad master equation, one can see that the operators Ĵi must vanish when the evolution is that of a closed system: therefore we see that they must describe an effect of the environment. The Lindblad master equation can be derived in many other ways more formal ways [START_REF] Breuer | The theory of open quantum systems[END_REF][START_REF] Carmichael | An open systems approach to quantum optics[END_REF]. Let us see how the operators Ĵi can be physically interpreted as jump operators, that induce a collapse of the system wavefunction due to a continuous weak measurement process.

III.5 The stochastic Schrodinger equation and quantum trajectories

Another way to derive the Lindblad master equation is by considering a physical system under continuous measurement. By doing this, we obtain what is known as an unraveling of the master equation in terms of a stochastic process. The picture here is that of a physical system that stochastically undergoes jumps from one state to another, with a collapse of the wavefunction occuring at random times. For example, consider an atomic system at room temperature. At each instant, there is a nonzero probability that excitations are lost to the thermal environment. This can be seen as a stochastic process with two outcomes: either the excitation is lost or it is not. In the general case, this picture is well described by the generalized measurement picture: the quantum state after a generalized measurement is given by:

ρ(t) → ρ(t + δt) = nµ µ Mµ ρ(t) M † µ p µ (t) , ( 1.39)
with Mµ the measurement operators. That is, by performing a measurement of the system described by ρ(t), the state has collapsed to a certain state, but we are uncertain about which state it is, of which there are n µ possible outcomes. p µ (t) = ⟨ M † µ Mµ ⟩(t) is therefore the probability of measuring the µth outcome if a projective measurement is carried out. Let us consider the following choice of measurement operators

   M0 = 1 -(i Ĥ + Ĵ † Ĵ/2)δt, M1 = Ĵ√ δt, ( 1.40)
where we have supposed that there is only one jump operator Ĵ. This closely corresponds to what was given in Eq. (1.33): the density matrix gets updated in time as under the action of the Lindblad master equation. The dynamics of an open system can thus be interpreted as the result of a continuous weak measurement of the jump operator Ĵ by the environment. The time at which a jump happens is non-deterministic in this view. and its probability is given by p

1 = ⟨ M † 1 M1 ⟩(t) = ⟨ Ĵ † Ĵ⟩(t)δt.
We refer to a specific realisation of the stochastic process as a quantum trajectory. In this picture, the density matrix is interpreted statistically, and represents our ignorance of the particular realisation of the outcome record r(t). When M1 is measured (r(t) = 1), the state vector becomes .41) This corresponds to the collapse of the wave function upon the environment having successfully measured the outcome of L. If instead M0 is measured (r(t) = 0), corresponding to no jump having occured with associated probability p 0 = 1 -p 1 , the system's state vector becomes

|ψ 1 (t + δt)⟩ = M1 |ψ(t)⟩ p 1 (t) = Ĵ ⟨ Ĵ † Ĵ⟩(t) |ψ(t)⟩ . ( 1
|ψ 0 (t + δt)⟩ = M0 |ψ(t)⟩ √ p 0 = 1 -dt i Ĥ - 1 2 Ĵ † Ĵ + 1 2 ⟨ Ĵ † Ĵ⟩(t) |ψ(t)⟩ . (1.42)
Interestingly, in this picture, the coupling of the system to its environment, even in the event of no jump, influences the time-evolution of the system. By combining Eqs. (1.41) and 1.42, and considering infinitesimal time intervals dt and outcome records dr(t) 5 , the evolution of the wave function is given by the following stochastic differential equation [START_REF] Wiseman | Quantum Measurement and Control[END_REF]:

d |ψ(t)⟩ = (1 -dr(t)) |ψ 0 (t + dt)⟩ + dr(t) |ψ 1 (t + dt)⟩ -|ψ(t)⟩ = dt ⟨ Ĵ † Ĵ⟩(t) -Ĵ † Ĵ 2 -i Ĥ |ψ(t)⟩ + dr(t) Ĵ/ ⟨ Ĵ † Ĵ⟩(t) -1 |ψ(t)⟩ (1.43)
5 dr(t) is an infinitesimal stochastic variable, also known as a Wiener process.

which is known as the stochastic Schrödinger equation. This leads to the following picture: the wavefunction evolves in time according to a non-Hermitian Hamiltonian, with sudden quantum jumps occuring at random times. Particular realisations of this process correspond to conditional single trajectories {|ψ (i) ⟩} i . It can be shown that one recovers the density matrix ρ(t) at a given time by averaging over infinitely many trajectories ρ(t) = lim

N traj →+∞ 1 N traj N traj i=1 |ψ (i) (t)⟩⟨ψ (i) (t)| . (1.44)
thus the map given in Eq. (1.43) is equivalent to that of the Lindblad master equation, Eq. (1.38). It follows that any observable can be obtained analogously:

Tr[ρ(t) Ô] = lim N traj →+∞ 1 N traj N traj i=1 ⟨ψ (i) (t)| Ô|ψ (i) (t)⟩ . (1.45)
In the numerical methods section, we will see how this translates to a straightforward numerical method to reduce the complexity of integrating the Lindblad master equation in time.

III.6 Numerical methods

Now that the main theoretical framework to treat open quantum systems has been presented, let us examine in more detail some of the numerical methods that have been developed over the years to study such systems. There are two main objects of study with open quantum systems: (i) finding the steady-state and/or the Liouvillian gap, which is the distance between the two smallest eigenvalues of the Liouvillian, important in the study of dissipative phase transitions, and (ii): simulating the time-evolution to investigate transient regimes. This section and most results in this manuscript mostly focus on (ii).

Monte-Carlo wave function (MCWF)

A widely used numerical method that naturally appears from the stochastic Schrödinger equation is the Monte-Carlo wave function (MCWF) method [START_REF] Mølmer | Monte-Carlo wave-function method in quantum optics[END_REF]. This method consists in evolving a quantum state according to a non-Hermitian Hamiltonian, and applying quantum jumps by randomly choosing from the set of jump operators Ĵi describing the system-environment interaction, according to some probability. The algorithm to simulate the dynamics of an open quantum system can be decomposed in the following steps:

1. Draw a random number a uniformly between 0 and 1.

2. Time-evolve the wavefunction |ψ(t)⟩ according to the non-Hermitian Hamiltonian H until time t ′ , defined by ⟨ψ(t ′ )|ψ(t ′ )⟩ = a.

3. Compute the jump probabilities p i associated to each jump operator Ĵi , given by

p i = ⟨ψ(t ′ )| Ĵ † i Ĵi |ψ(t ′ )⟩ = ⟨ Ĵ † i Ĵi ⟩(t ′ ) (1.46)
4. Select the ith jump operator with probability p i and apply it to |ψ⟩ (t ′ ) along with renormalizing it as:

ψ(t ′ ) = Ĵi ⟨ Ĵ † i Ĵi ⟩(t ′ ) |ψ(t ′ )⟩ . (1.47)
5. Repeat until the desired time is reached by restarting step 1 with the state ψ(t ′ ) .

This algorithm has many variants, but this is computationally the most efficient [START_REF] Krämer | QuantumOptics.jl: a julia framework for simulating open quantum systems[END_REF] with respect to other implementations. By repeating this N traj times, one obtains dynamics that are closer to the exact solution. However, how many trajectories are needed in practice highly depends on the problem, and no general procedure exists to determine N traj [START_REF] Daley | Quantum trajectories and open many-body quantum systems[END_REF]. The complexity of this method is given by O(N traj l L ), an exponential advantage with respect to a full integration of the Lindblad master equation. This is particularly useful for highly dissipative systems and when the number of dissipative channels is high, as we will see in more detail in chapter 3.

Mean-field and cluster approaches

The ideas presented in section II.3 for mean-field ansätze can be naturally extended to open quantum systems. Rather than taking the wavefunction to be a product of local wavefunction, we take the density matrix to be a product of local density matrices:

ρ = Nc i=1 ρi , (1.48)
with N c (N c = L in the Gutzwiller mean-field approach) the number of clusters and ρi ∈ B(H i ). As for closed systems, this approach limits quantum correlations to be inside the cluster, meaning it will work best for weakly and/or locally interacting systems. This has the same value as for closed systems, reducing the number of free parameters required to describe the system density matrix ρ from l L to N c l Nc/N , with N c /N the size of the cluster. This leads to N c /N equations to solve to integrate the master equation in time, leading to a more tractable solution. This approach has been quite successful in capturing a certain number of features of dissipative phase transitions [START_REF] Jin | Cluster mean-field approach to the steady-state phase diagram of dissipative spin systems[END_REF] and nonequilibrium properties of strongly-correlated open quantum systems [START_REF] Boité | Steady-state phases and tunneling-induced instabilities in the driven dissipative Bose-Hubbard model[END_REF].

Matrix Product Operators

Tensor networks can also be applied to open quantum systems. The problem amounts to finding an expression for the many-body density matrix coefficients in terms of local tensors, much like Eq. (1.16). An additional issue when dealing with density matrices is that positivity must be enforced, without ever constructing the full density matrix 6 . A simple ansatz can be constructed by using the Choi isomorphism [START_REF] Choi | Completely positive linear maps on complex matrices[END_REF] that associates a density matrix to a pure state in an enlarged Hilbert space of dimension l 2L . The then vectorised density matrix can be described like an MPS in this Hilbert space, and the density matrix coefficients defined in Eq.(1.2) can be written as a tensor product:

c s ′ 1 ,s ′ 2 ,...,s ′ L s 1 ,s 2 ,...,s L = A i 1 ,i 2 s 1 ,s ′ 1 A i 2 ,i 3 s 2 ,s ′ 2 • • • A i N ,i 1 s L ,s ′ L (1.49)
which is equivalent to the MPS ansatz with an additional dimension for each tensor that corresponds to the set of spin indices {s i }. This is known as the matrix product operator (MPO) ansatz. While it has proved effective and yielded interesting results for 1D systems, one of the main issues of such a construction is that the positivity of the density matrix is not ensured when performing time evolution of the ansatz. Note that this problem has been proved to be computationally hard [START_REF] Kliesch | Matrix-product operators and states: nphardness and undecidability[END_REF], hence no approach can be a perfect solution 7 .

A partial solution to this issue is to built what is known as locally-purified density operators (LPDO) [START_REF] Werner | Positive tensor network approach for simulating open quantum many-body systems[END_REF]. This is done by constructing a locally enlarged Hilbert space, built as

H ⊗ A = H 1 ⊗ A 1 ⊗ H 2 ⊗ A 2 ⊗ . . . H L ⊗ A L
with A the ancilla space. A wavefunction is then written as an MPS in this basis, whose coefficients are given by:

⟨s 1 , a 1 s 2 , a 2 , . . . , s L , a L |ψ⟩ = c s 1 ,a 1 s 2 ,a 2 ,...,s L ,a L = A i 1 ,i 2 s 1 ,a 1 A i 2 ,i 3 s 2 ,a 2 • • • A i N ,i 1 s L ,a L . (1.50)
The density matrix coefficients are then obtained by tracing out the ancilla sites on |ψ⟩⟨ψ|:

c s ′ 1 ,s ′ 2 ,...,s ′ L s 1 ,s 2 ,...,s L = B i 1 ,i 2 s ′ 1 ,s ′ 1 B i 2 ,i 3 s ′ 2 ,s ′ 2 • • • B i N ,i 1 s ′ L ,s ′ L . ( 1.51)
Note that this expression is not equivalent to the MPO ansatz, as the B tensors are obtained by doing a partial trace on the wave function |ψ⟩, and the MPO simply corresponds to a vectorized density matrix. The LPDO ansatz ensures the positivity of the density matrix, as it is obtained by a partial trace on a quantum state. However, some works have shown that the required bond dimension to obtain a similar precision is higher with LPDOs than MPOs [START_REF] Cuevas | Purifications of multipartite states: limitations and constructive methods[END_REF][START_REF] Cuevas | Fundamental limitations in the purifications of tensor networks[END_REF]. This hints to a lower capacity of LPDO to represent highly entangled states. Finally, we note that most works employ the MPO ansatz without worrying about the positivity of the resulting density matrix while obtaining highly accurate results [START_REF] Biella | Photon transport in a dissipative chain of nonlinear cavities[END_REF][START_REF] Cui | Variational matrix product operators for the steady state of dissipative quantum systems[END_REF].

Corner-space renormalization

A method that is of particular interest for strongly-correlated systems is the corner-space renormalization method [START_REF] Finazzi | Corner-space renormalization method for driven-dissipative two-dimensional correlated systems[END_REF]. This method is rather different from mean-field or MPO methods, as it is constructed from the idea that there exists a class of physical systems, low-entropy quantum systems, for which the density matrix is accurately described by only a restricted set of M states. Mathematically, this is written as

ρ = N k p k |ϕ k ⟩⟨ϕ k | = M k p k |ϕ k ⟩⟨ϕ k | + O(ϵ M ) (1.52) with ϵ M = 1 -M k p k ≪ 1
. This means that one can truncate the eigenvalues to keep only M states, and the truncated density matrix will still faithfully describe the physical system at hand. The corner-space renormalization method aims to solve the steady-state ρ ∈ H of an open quantum system, described by Eq. (1.38). The basic idea behind this method is shown in Fig. 1.6. It can be decomposed in three steps, as pictured by the three light rectangles forming the figure :   1. Take two-subsystems A and B whose steady-state can be found exactly, and obtain the steady-state by brute-force integration. Diagonalize the two corresponding steady-states ρ(A) , ρ(B) to obtain the form given in Eq.(1.52). 3. Repeat until the size of the full system is reached.

Note that this method allows for the states in the solution to be entangled between subsystems A and B: product states are constructed as the basis for the solution of the steady-state, hence this steady-state can comprise any superpositions of these product states. In this sense, this method is more powerful than mean-field or tensor network methods. Its limitation mainly resides in the cut-off of the M states, which will be limiting for high-temperature states with a large entropy. Note that the complexity of the algorithm for this method resides in the diagonalization of the density matrix at the end of the first step. This complexity of this step is O(max[M 3 , l 3N A,B ]), since the first step is done with exact diagonalization of subsystems A, B with Hilbert space dimension l N A,B , and subsequent steps concern density matrices of size M × M .

A method to simulate the dynamics of open quantum systems, based on the same principle that only M states may suffice to describe low-entropy systems, will be introduced in chapter 3.

III.7 Experimental systems

As mentioned, any quantum system that couples to an external environment, provided the Born-Markov approximation is satisfied, can be depicted as an open quantum system.

Superconducting circuits

A class of experimental setups that can be described as open systems are superconducting circuits. These circuits have proved immensely useful for the development of quantum information in the past two decades [START_REF] Wendin | Quantum information processing with superconducting circuits: a review[END_REF], as well as for exploring many-body physics [START_REF] Blais | Quantum information processing and quantum optics with circuit quantum electrodynamics[END_REF]. One of the building blocks of such systems are electromagnetic LC resonators, composed of a capacitor C and an inductor L. One can describe a linear resonator with the Hamiltonian

Ĥ = 1 2C Q2 + 1 2L φ2 (1.53)
where Q is the quantized charge of the capacitor, and φ is the quantized magnetic flux flowing through the inductor. These operators are described by the same commutation rules as for position and momentum, i.e. [ Q, φ] = -i, hence the Hamiltonian can be diagonalized in exactly the same way as for a harmonic oscillator, by introducing the bosonic annihilation (and creation) operator, defined as:

â = 1 √ ω   φ √ 2L -i Q √ 2C   (1.54)
and the Hamiltonian is written as

Ĥ = ω â † â + 1 2 . (1.55)
Such a Hamiltonian has a linear spectrum of n levels corresponding to n bosonic excitations, and its energy levels are given by

E n = ω n + 1 2 , (1.56)
with ω the frequency of the resonator. What one really wants to investigate issues related to nonequilbrium dynamics is nonlinear terms in the Hamiltonian, so that the states involved in the dynamics are strongly correlated. The key ingredient to creating a nonlinear superconducting circuit is the Josephson junction [START_REF] Josephson | Possible new effects in superconductive tunnelling[END_REF], that adds nonlinearities to resonators. These are hybrid structures, composed of two superconductors separated by an insulating barrier, that behave as nonlinear inductors. The Hamiltonian term describing this element is given by ĤJ = E J cos φJ (1.57) with E J = I 0 ϕ 0 /2π (I 0 is the current passing through the junction, ϕ 0 is the flux quantum) and φJ is the flux flowing through the junction, that is again quantized. The cosine term is at the origin of the nonlinearity. By building a circuit made of a capacitor and a Josephson junction, one has the following Hamiltonian, that can be Taylor-expanded by supposing that the system is in a weakly nonlinear regime:

Ĥ = 1 2C Q2 + E J cos φJ (1.58) Ĥ = 1 2C Q2 + E J 1 - 1 2 φ2 J + 1 24 φ4 J + O φ6 J . (1.59)
The spectrum is now slightly modified, and the energy separation between levels is not homogeneous. This leads to the following Hamiltonian, in terms of creation and annihilation operators b, b † :

Ĥ = ω b † b + 1 2 -E J -E C b † + b 4 (1.60) with ω = √ 8E C E J , E C = e 2 /2C
, e being the charge of the electron [START_REF] Bishop | Circuit quantum electrodynamics[END_REF]. This leads to being able to realize two-level systems known as transmon qubits, crucial for quantum information provided the nonlinear terms are strong enough, such that a resonant drive at energy ω only couples the first two levels. Alternatively, one can study nonlinear lattice models with such platforms, such as the Bose-Hubbard model. This model is described by the following Hamiltonian:

Ĥ = L i=1 ω c b † i bi + U 2 b † i b † i bi bi -J L-1 i=1 b † i bi+1 + H.c. , (1.61)
for a 1D lattice with open boundary conditions, whose features will be studied in more detail in Chapter 2. ω c is the cavity frequency, J the nearest-neighbour coupling, and U the boson-boson interaction strength. The coupling J can be engineered by capacitive couplings, that are highly tunable elements [START_REF] Wendin | Quantum information processing with superconducting circuits: a review[END_REF]. Such systems exhibit natural loss channels, due to their coupling to the external environment. A simple model for such loss channels is that of local single-particle losses, that are described by jump operators

Ĵi = Γ 1 bi . (1.62)
Another ingredient that can be added is an external driving, that can be coherent (entering the Hamiltonian) or incoherent (entering the dissipative part of the Liouvillian). For a coherent drive in the rotating wave approximation, the following terms may be added to the Hamiltonian

ĤD = i F i ( b † i + bi ). (1.63)
This leads to richer physics, as the system consists of driven-dissipative coupled nonlinear elements. Many interesting phases can be engineered, and the physics of such systems has only begun to be explored. In Fig. 1.7, a driven-dissipative superconducting platform is depicted. The full system is about 10 mm 2 , and contains eight coupled transmon resonators, whose observables can be measured thanks to readout resonators. The chain of transmons is driven by a coherent driving on an external site, and another site is coupled to the system, acting as an energy dump. This enabled authors of ref. [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF] to stabilize strongly correlated phases, coined a Mott insulator of photons, that was never accessed before.

Entanglement dynamics in strongly-correlated open quantum systems I Introduction

After being a subject of early intense debate at the dawn of quantum mechanics [START_REF] Schrödinger | Die gegenwärtige situation in der quantenmechanik[END_REF][START_REF] Einstein | Can quantum-mechanical description of physical reality be considered complete?[END_REF], entanglement is now recognized as a key feature of quantum physics [START_REF] Horodecki | Quantum entanglement[END_REF]. The efforts towards building a complete mathematical description of this notion were instrumental in the development of quantum information. In this context, the core of the theory is centered around three main tasks: detecting [START_REF] Gühne | Entanglement detection[END_REF], quantifying [START_REF] Plenio | An introduction to entanglement measures[END_REF] and manipulating entanglement [START_REF] Bennett | Purification of noisy entanglement and faithful teleportation via noisy channels[END_REF]. The progress made on these three fronts would allow to outperform classical methods in the fields of metrology [START_REF] Giovannetti | Advances in quantum metrology[END_REF], cryptography [START_REF] Gisin | Quantum cryptography[END_REF] and computation [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF].

In addition to providing sound foundations to the field of quantum information, entanglement theory has also paved the way to new discoveries in other areas of physics. As anticipated at the beginning of the millenium [START_REF] Preskill | Quantum information and physics: some future directions[END_REF], quantities such as the entanglement entropy have proved to be very valuable tools for characterizing the ground-state wave function of many-body quantum systems [START_REF] Amico | Entanglement in many-body systems[END_REF][START_REF] Chiara | Genuine quantum correlations in quantum manybody systems: a review of recent progress[END_REF][START_REF] Islam | Measuring entanglement entropy in a quantum many-body system[END_REF]. The study of entanglement in manybody systems has not been restricted to their ground state properties: entanglement dynamics and its propagation in space in quantum systems has also been the subject of intense research activities for spin chains [START_REF] Amico | Dynamics of entanglement in one-dimensional spin systems[END_REF][START_REF] Jurcevic | Quasiparticle engineering and entanglement propagation in a quantum many-body system[END_REF][START_REF] Kastner | Entanglement-enhanced spreading of correlations[END_REF], fermionic [START_REF] Alba | Spreading of correlations in Markovian open quantum systems[END_REF] and bosonic systems [START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF][START_REF] Läuchli | Spreading of correlations and entanglement after a quench in the one-dimensional Bose-Hubbard model[END_REF][START_REF] Daley | Measuring entanglement growth in quench dynamics of bosons in an optical lattice[END_REF][START_REF] Bernier | Light-Cone and Diffusive Propagation of Correlations in a Many-Body Dissipative System[END_REF]. These works were important in inspecting the validity and limits of predictions about Lieb-Robinson bounds in lattice systems, as well as providing information about the properties and excitations of complex many-body systems.

Extending these investigations to open quantum systems represents a timely frontier of research [START_REF] Aolita | Open-system dynamics of entanglement:a key issues review[END_REF] that is of fundamental importance because much less is known with respect to the state-of-the-art in isolated quantum many-body systems at thermal equilibrium or exhibiting unitary Hamiltonian dynamics. Whereas in general experimentalists try to protect their system from interacting with its environment, other approaches based on the general concept of "reservoir engineering" try to exploit the openness of a system and take advantage of judiciously designed dissipation to reach non-trivial quantum states in the transient regime [START_REF] Leghtas | Confining the state of light to a quantum manifold by engineered two-photon loss[END_REF] or in the steady state [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF]. In recent years, experimental progress in tailoring effective photon-photon interactions in cavity and circuit quantum electrodynamics (QED) devices has lead to the emergence of controllable quantum optical many-body systems [START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Schmidt | Circuit qed lattices: towards quantum simulation with superconducting circuits[END_REF][START_REF] Roy | Colloquium: strongly interacting photons in one-dimensional continuum[END_REF]. Unlike most condensed matter setups where the system is close to thermal equilibrium, this new class of systems are open quantum platforms in which intrinsic losses, due to the photon finite lifetime, have to be compensated by an external coherent or incoherent driving.

In the past years, several works have been devoted to transport properties of stronglycorrelated photonic platforms [START_REF] Biella | Photon transport in a dissipative chain of nonlinear cavities[END_REF][START_REF] Lee | Few-photon transport in manybody photonic systems: a scattering approach[END_REF][START_REF] Mertz | Photonic currents in driven and dissipative resonator lattices[END_REF][START_REF] Debnath | Nonequilibrium photonic transport and phase transition in an array of optical cavities[END_REF]. However, only a few exist regarding entanglement and correlation propagation in driven-dissipative systems, which mainly focus on free fermion systems [START_REF] Alba | Spreading of correlations in Markovian open quantum systems[END_REF][START_REF] Maity | Growth of mutual information in a quenched one-dimensional open quantum many-body system[END_REF]. The recent experimental demonstration of dissipatively stabilized photonic Mott insulators [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF] in chains of superconducting microwave resonators paves the way to the exploration of such an exciting frontier.

In this chapter, after presenting the physics of strongly-interacting bosonic systems, we will explore entanglement propagation in photonic Mott insulators, showing genuine physical effects associated to the openness of such systems. In contrast to most works about correlation propagation in interacting bosonic systems, here we do not consider global quenches of the system that typically consist in abruptly changing the value of the interaction strength in all the lattice [START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF][START_REF] Bernier | Light-Cone and Diffusive Propagation of Correlations in a Many-Body Dissipative System[END_REF]. Instead, we consider two configurations where one photon is injected or removed from one cavity in the middle of a chain and investigate the propagation of entanglement that is produced between distant sites as a function of time and of their spatial separation. Such a study is achieved by monitoring the negativity of two-site reduced density matrices, that witnesses entanglement. We show a strinking different role of photon dissipation in the two configurations.

II The Bose-Hubbard model

Let us now present a few features of the paradigmatic Bose-Hubbard model. This model consists of L interacting bosonic modes, described by the following Hamiltonian:

Ĥ = L i=1 ω c b † i b i + U 2 b † i b † i b i b i -J L-1 i=1 (b † i b i+1 + h.c.), ( 2.1)
with ω c the cavity mode frequency, U the photon-photon (Kerr) on-site interaction, J the nearest-neighbor photon hopping coupling, and b i (b † i ) the annihilation (creation) photon operators for each site. It presents a Mott phase in the strongly interacting (U ≫ J) limit and a superfluid phase in the weakly interacting (U ≪ J) regime. In this section we restrict the discussion to one-dimensional systems for simplicity. The physical systems described by the Bose-Hubbard Hamiltonian include, but are not limited to, lattices of microwave resonators in circuit QED platforms as described in Chapter 1 [START_REF] Houck | On-chip quantum simulation with superconducting circuits[END_REF][START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Fitzpatrick | Observation of a dissipative phase transition in a one-dimensional circuit qed lattice[END_REF][START_REF] Hartmann | Quantum simulation with interacting photons[END_REF][START_REF] Noh | Quantum simulations and many-body physics with light[END_REF][START_REF] Carusotto | Photonic materials in circuit quantum electrodynamics[END_REF], semiconductor microcavities [START_REF] Carusotto | Quantum fluids of light[END_REF][START_REF] Goblot | Nonlinear polariton fluids in a flatband reveal discrete gap solitons[END_REF] and ultracold gases in optical lattices [START_REF] Greiner | Quantum phase transition from a superfluid to a mott insulator in a gas of ultracold atoms[END_REF][START_REF] Jaksch | Cold bosonic atoms in optical lattices[END_REF]. These platforms, as most realistic many-body systems, exhibit dissipation and dephasing due to the coupling to the environment. In cold atom systems dephasing is dominant [START_REF] Bernier | Light-Cone and Diffusive Propagation of Correlations in a Many-Body Dissipative System[END_REF] while for microwave photons in circuit QED platforms particle loss is typically the most important channel [START_REF] Pichler | Nonequilibrium dynamics of bosonic atoms in optical lattices: decoherence of many-body states due to spontaneous emission[END_REF][START_REF] Hartmann | Quantum simulation with interacting photons[END_REF].

II.1 Strongly correlated regime

One can start by considering the limit U/J → ∞. In this limit, one can exactly perform a fermionization procedure in 1D, which enables one to solve the model [START_REF] Girardeau | Relationship between systems of impenetrable bosons and fermions in one dimension[END_REF]. In such a limit, bosons become impenetrable, and the Hamiltonian reads:

H = ω c L i=1 b † i b i -J L-1 i=1 (b † i b i+1 + h.c.), (2.2)
which corresponds to a tight-biding model with a local energy ω c . By using a Jordan-Wigner transform, this model is exactly solvable [START_REF] Girardeau | Relationship between systems of impenetrable bosons and fermions in one dimension[END_REF][START_REF] Paredes | Tonks-girardeau gas of ultracold atoms in an optical lattice[END_REF].

In the following we will focus on a more interesting and practical case, where U ≫ J (but not so much larger). This does not completely wash out interactions, hence there will be some interplay between interactions and dissipative processes. In such a regime, in order to describe the physics of a photonic Mott insulator with one photon per site on average, we can safely truncate the local Hilbert space to a maximum of two photons per site by retaining only the |0⟩, |1⟩, and |2⟩ Fock number states. The validity of this assumption was carefully tested numerically by increasing the local Hilbert space cutoff and checking that the results were unchanged. A Mott insulator phase corresponding to one photon per site for U ≫ J is approximately described by the factorized state

|Ψ Mott ⟩ = |1⟩ 1 ⊗ |1⟩ 2 ⊗ ... ⊗ |1⟩ L = |11...1⟩ . ( 2.3)
In the regime of strong interactions, the Hamiltonian can be diagonalized by using generalized Jordan-Wigner and Bogoliubov transformations [START_REF] Cheneau | Light-cone-like spreading of correlations in a quantum many-body system[END_REF], via a mapping to a spin-1 model. By considering L coupled three level systems, one can rewrite the Hamiltonian as interacting doublon and holon excitations. These fermionic quasiparticles are described by local creation operators d † i and h † i for doublons and holons respectively, such that

d † i |1⟩ i = |2⟩ i , h † i |1⟩ i = |0⟩ i .
We therefore have:

Ĥ = L j=1 P -2J d † j dj+1 -J ĥ † j ĥj+1 - √ 2J d † j ĥ † j+1 -ĥj dj+1 + H.c. + U 2 (n d,j + nh,j ) P
(2.4) with P = Π j (1 -nd,j nh,j ). One sees that the complexity of the model resided in the projector P , which contains non-quadratic terms. However, one can reasonably approximate P ≈ 1, since the probability of having both a doublon and a holon on the same site |⟨n d,j nh,j ⟩| 2 should be 0. By doing this, one obtains a quadratic Hamiltonian that we can diagonalize using a Bogoliubov transform. The eigenmodes of the system are given by:

γ † d,k = u(k)d † k + v(k)h † -k (2.5) γ † h,-k = u(k)h † -k -v(k)h † k (2.6) with u(k) = cos(θ k /2), v(k) = i sin(θ k /2), (2.7
)

θ k = arctan √ 32J sin k U -6J cos k (2.8) dk = 2 L + 1 i sin(ki) di (2.9) ĥ-k = - 2 L + 1 i sin(ki) ĥi (2.10) k = πn L + 1 , n ∈ N for open boundary conditions. (2.11)
The Hamiltonian can now be written in a diagonal form:

Ĥ = k ϵ d (k)γ † d,k γd,k + ϵ d (-k)γ † h,-k γh,-k (2.12) with ϵ d (k) = -J cos k + 1 2 (U -6J cos k) 2 + 32J 2 sin 2 k (2.13) ϵ h (-k) = J cos k + 1 2 (U -6J cos k) 2 + 32J 2 sin 2 k. ( 2.14)
For a closed system, we can extract the velocity for the corresponding doublons and holons as the derivative of the eigenenergies with respect to k:

v d,h (k) = ∂ k ϵ d,h (k) = ±J sin k + 1 4 12JU sin k -8J 2 sin k cos k (U -6J cos k) 2 + 32J 2 sin 2 k (2.15)
with ± respectively corresponding to doublons and holons. This enables one to obtain bounds on the quasiparticle speed, which for doublons reads

v D max = 4J 1 -4J 2 U 2 +O J 3 U 4 and v H max = 2J 1 + 17J 2 2U 2 + O J 3 U 4
for holons. These results will serve as a basis to understand what happens in terms of entanglement propagation and the effect of dissipation and decoherence in the case of an open system.

III Entanglement propagation in dissipative photonic Mott insulators

We now consider the system described by the Hamiltonian of Eq. (B.1) coupled to an external environment. Within an open quantum systems approach, the time evolution of the system density matrix ρ can be described by the following Lindblad master equation:

dρ dt = -i[ Ĥ, ρ] + 1 2 L i=1 C 2J (C) i ρJ (C) † i -{J (C) † i J (C) i , ρ}, (2.16) J i c i 1 i c+1 i c-1 i L γ U J γ γ γ γ γ J γ γ γ γ U U J Time t = 0 0.0 0.1 0.2 0.3 0.4 0.5

III.1 Entanglement generation protocol

Since |Ψ Mott ⟩ is a factorized state, an interesting question is how to perturb such a photonic Mott insulator in order to create entanglement in a simple way and study its propagation in a direct fashion. In the following we will show that this is possible by injecting (or removing) one photon from an occupied site. As shown in the upper panels of Figs.

2.1 and 2.2, we will consider such manipulation on the central site of a linear chain of resonators. In the case of a photonic insulator with a large U , this can be achieved simply by applying a coherent π-pulse drive on the central site that induces a Rabi rotation from the |1⟩ to the |2⟩ (or to the |0⟩) Fock number state in the considered site. We have explicitly verified that such operation can be performed with fidelity close to 1 thanks to the strong anharmonicity produced by the large on-site interaction U . This way, it is possible to prepare the state

|Ψ D ⟩ = |1...2...1⟩ (|Ψ H ⟩ = |1...0...1⟩)
where D (H) stands for doublon (holon), corresponding to the injection of a single localized excitation on top of the quasiparticle vacuum. In the doublon case, the excess photon in the central site can hop to the right nearest-neighbor site or, with the same probability, to the left site. Due to the symmetry of the chain with respect to the central site and the lack of which-path information, such propagation creates an entangled state that can propagate along the chain. For circuit QED platforms, a Mott insulator can be prepared and maintained through an active stabilization process [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF]. In the following, we will consider the dissipative dynamics of the system after the creation of the localized doublon (holon) in the absence of stabilization.

III.2 Entanglement detection

In order to witness bipartite entanglement between two partitions A and B of the system, we have considered the negativity function [START_REF] Życzkowski | Volume of the set of separable states[END_REF], defined as N (ρ AB ) = λ<0 |λ|, where the sum is taken over the negative eigenvalues λ of ρ T A AB , which is the partial transpose with respect to subsystem A of the joint density matrix ρ AB . In the following, A and B will be two resonators at a distance r from the central cavity (i = i c ). The time-dependent negativity of the reduced density matrix for these two sites at the positions i c-r and i c+r will be denoted by N r (t). For systems with a relatively moderate Hilbert space dimension, we can compute the time evolution of the full density matrix of the system via an exact integration of the master equation. Once we get the full density matrix, we can trace out with respect to the degrees of freedom of all the sites except the two sites under study.

In the regime where we can consider only the |0⟩, |1⟩, and |2⟩ states as local basis for a given site, we can reconstruct a two-site reduced density matrix by exploiting the fact that any Hermitian operator can be decomposed over the generators of the group associated to its Hilbert space [START_REF] Thew | Qudit quantum-state tomography[END_REF]. As we truncated the local Hilbert space to states with up to 2 photons, the generators of the SU(3) ⊗ SU(3) group allows us to reconstruct the reduced density matrix as:

ρ (2) 3 = 1 9 8 i 1 ,i 2 =0 r i 1 i 2 Λ (i 1 ) ⊗ Λ (i 2 ) ,
(2.17) with Λ (i) the generators 1 of the SU(3) group and

J i c i 1 i c+1 i c-1 i L γ J γ γ γ γ γ J γ γ γ γ J Time t = 0 0.0 0.2 0.4 0.6 0.
r i 1 i 2 = 9⟨Λ (i 1 ) ⊗ Λ (i 2 ) ⟩ tr((Λ (i 1 ) ⊗ Λ (i 1 ) ) 2 )
.

Using this method we can reconstruct the two-site density matrix for all times and obtain the entanglement negativity. Since the Hermitian operators Λ (i) can also be expressed as a function of the bosonic creation and annihilation operators up to their third-order power 2 , this tomographic method can be used in experiments to measure the two-site entanglement negativity. Indeed, the measurement of expectation values of moments of the photon fields has become a rather standard procedure in circuit QED platforms (see, 1 The nonzero matrix elements of the Hermitian Λ (i) matrices, such that Λ

(i) r,s = (Λ (i) s,r ) ⋆ , are the following ones: Λ (0) 11 = Λ (0) 22 = Λ (0) 33 = 1, Λ (1) 12 = 1, Λ (2) 12 = -i, Λ (3) 11 = -Λ (3) 22 = 1, Λ (4) 13 = 1,Λ (5) 13 = -i, Λ (6) 23 = 1, Λ (7) 23 = -i, Λ (8) 22 = 1/ √ 3 and Λ (8) 33 = -2/ √ 3. 2
The eight generators Λ (i) can be expressed in terms of the bosonic annihilation and creation operators. Namely:

Λ (1) = 1 2 (bb †2 + b 2 b † ), Λ (2) = i 2 (b 2 b † -bb †2 ), Λ (3) = 3 4 b 2 b †2 -1 2 bb † , Λ (4) = 1 √ 2 (b †2 + b 2 ), Λ (5) = i √ 2 (b 2 -b †2 ), Λ (6) = 1 √ 2 (b †2 b + b † b 2 ), Λ (7) = i √ 2 (b † b 2 -b †2 b), and Λ (8) = 1 √ 3 (bb † -b † b).
e.g. [START_REF] Besse | Parity detection of propagating microwave fields[END_REF]). Note that we have conveniently developed and used this approach for numerical simulations based on the Matrix Product Operator (MPO) technique [START_REF] Biella | Photon transport in a dissipative chain of nonlinear cavities[END_REF][START_REF] Zwolak | Mixed-state dynamics in one-dimensional quantum lattice systems: a time-dependent superoperator renormalization algorithm[END_REF][START_REF] Verstraete | Matrix product density operators: simulation of finite-temperature and dissipative systems[END_REF]. Indeed, MPO simulations are effective to simulate longer chains of cavities but do not allow for a direct access to the full system density matrix, an issue that we bypassed with the procedure described above.

IV Results and discussion

In Fig. 2.1, we report results for the negativity N r (t) for different values of r and of the photon loss rate γ (here no pure dephasing is considered, Γ d = 0). The negativity shows a well resolved peak for most values of the spatial separation r and γ/J: increasing r delays the negativity peak, showing a clear entanglement propagation. A revival peak of entanglement is visible in the r = 1 curve at longer times. Previous studies of two-qubit systems with non-Markovian environments have revealed entanglement revival effects [START_REF] Bellomo | Non-markovian effects on the dynamics of entanglement[END_REF][START_REF] Ficek | Dark periods and revivals of entanglement in a two-qubit system[END_REF][START_REF] Galve | Entanglement dynamics of nonidentical oscillators under decohering environments[END_REF]. In our system, the two-site dynamics is non-unitary and non-Markovian even for γ = 0 since the other sites of the chain have been traced out for the calculation of the negativity. The value of the negativity peaks decreases with increasing dissipation γ. However, it is remarkable that the entanglement propagation speed is negligibly influenced by dissipation and remains essentially ballistic. In Fig. 2.2, we report the analogous dynamics of negativity for the other configuration where the holon state |Ψ H ⟩ is prepared. With the same parameters as in Fig. 2.1, in the holon case not only the propagation speed is slower, but the role of dissipation is more dramatic, as we do not see any peaks for the chosen values of γ/J as soon as r > 1.

IV.1 Propagation speed

In Fig. 2.3, we report the calculated entanglement propagation speeds versus J/U both for the case of photon injection (doublon excitation) and extraction (holon excitation). In the same plot, we have also reported the maximal propagation speed of doublons and that of holons for a closed system, namely

v D max = 4J 1 -4J 2 U 2 + O J 3 U 4 , and v H max = 2J 1 + 17J 2 2U 2 + O J 3 U 4 .
The fact that the predicted entanglement speed for a closed system remains close to what we observed in the dissipative case is striking. This means that even in the presence of losses, the properties of the quasiparticle picture remain valid.

IV.2 Particle-hole asymmetry

The effect of dissipation for the two considered configurations is presented in Fig. 2.4 in which the peak value of N r=1 is plotted. To compare the genuine effect of dissipation and dephasing, we considered a holon propagation in a chain with hopping coupling 2J and a doublon propagation in a chain with hopping J in order to have the same speed (for U/J ≫ 1 the speeds differ by a factor 2). From the peak value of the negativity (occurring at the same time), we see that pure dephasing acts on the two cases in the same identical way (dashed lines), with an exponential decay of the negativity. Indeed, pure dephasing conserves the total number of particles and does not break the particlehole symmetry. On the other hand, in the presence of photon losses, our investigation reveals a striking asymmetry between the doublon (thick solid line) and holon (thin solid line) cases. Indeed, the negativity vanishes much faster for the holon configuration even when the speed is the same. The asymmetry can be qualitatively understood if we consider the quantum jump picture [START_REF] Mølmer | Monte-Carlo wave-function method in quantum optics[END_REF] for two sites that are entangled. In the quantum jump picture, the density matrix is seen as the statistical mixture of time-dependent quantum trajectories. These quantum trajectories represent the time-evolution of the wavefunction of the system along with sudden changes (quantum jumps) in the state of the wavefunction due to physical processes such as particle losses, that are inherently stochastic. In the holon case, by diagonalizing the reduced density matrix associated to two entangled sites, we have found that the entanglement is mostly due to the Bell state |ψ

H,+ ⟩ = 1 √ 2 |0⟩ i c-r |1⟩ i c+r + |1⟩ i c-r |0⟩ i c+r .
In a single quantum trajectory picture, a single quantum jump due to photon decay in one of the two sites transforms such state into the factorized state |0⟩ i c-r ⊗ |0⟩ i c+r . By contrast, in the doublon configuration, the entanglement is mostly due to the state |ψ

D,+ ⟩ = 1 √ 2 |2⟩ i c-r |1⟩ i c+r + |1⟩ i c-r |2⟩ i c+r .
In this case, the quantum jump produced by a photon loss in site i c-r produces the (normalized) state 1 dissipation rate. The holon negativity becomes exactly zero for a finite value of the dissipation γ whereas the decay in the doublon case slows down for increasing values of γ. Indeed, the two-site reduced density matrix satisfies an effective master equation, which is in general non-markovian as obtained by tracing out the other degrees of freedom of the chain, whose dynamics is non-trivial.

√ 3 √ 2|1⟩ i c-r |1⟩ i c+r + |0⟩ i c-r |2⟩ i c+r ,

V Conclusion

We have theoretically explored the physics of entanglement propagation in photonic Mott insulators in the presence of photon particle losses and dephasing. We have investigated a scheme where the entanglement is generated by injecting (or extracting) a photon from a site of a photonic Mott insulator. We have introduced a relatively simple quantum state tomography protocol, valid in the limit of strong photon-photon interactions, to study the bipartite entanglement properties. Our scheme is particularly suited for circuit QED platforms exhibiting strong photon-photon interactions and allowing the measurement of quantum optical correlation functions between distant sites. In spite of the losses, the propagation of the negativity peak exhibits a speed, which is close to the doublon (holon) quasiparticle propagation speed respectively in the case of the injection (extraction) of a photon. Remarkably, the impact of particle losses is highly asymmetric for these two configurations, while pure dephasing does not break the doublon-holon symmetry. Our work paves the way to new investigations on the entanglement propagation in open quantum systems, where multipartite entanglement spreading could also be investigated [START_REF] Horodecki | Quantum entanglement[END_REF][START_REF] Giorgi | Genuine quantum and classical correlations in multipartite systems[END_REF].

A future interesting research direction is the characterization of the entanglement dynamics at long times (diffusive vs. ballistic) and the quest for universal features underlying the dynamics of correlations in this class of systems. Another challenging problem to be investigated in the future, given the recent experimental success in the dissipative stabilization of photonic Mott insulators [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF], is the search for protocols to stabilize entanglement propagation in open quantum systems.

3

The dynamical corner-space method

In this chapter, we will introduce a numerical method to simulate the time-evolution of an open quantum system. This method is based on the same idea as the corner-space renormalization method, presented in Chapter 1, section III.6, which is that there exists a class of systems for which M pure states suffices to faithfully represent the density matrix, with M ≪ N , N being the Hilbert space dimension. This new method, the dynamical corner-space method, efficiently computes the time evolution of open quantum systems that have a moderate entropy with controllable accuracy, depending on the effective rank M of the density matrix. It will here be applied to the continuous-time dynamics of noisy quantum circuits in the presence of dissipation and decoherence, beyond digital error models that are commonly used in the quantum information community [START_REF] Arute | Quantum supremacy using a programmable superconducting processor[END_REF].

I Introduction

As we have seen in the previous chapter, strongly-correlated open quantum systems have remained relatively unexplored due to the lack of existing numerical methods to tackle them. In particular, we are interested in methods that put no restriction on the amount of entanglement that quantum states exhibits, unlike tensor network and mean-field methods.

A standard method that does this job is the Monte-Carlo wave function method [START_REF] Mølmer | Monte-Carlo wave-function method in quantum optics[END_REF][START_REF] Daley | Measuring entanglement growth in quench dynamics of bosons in an optical lattice[END_REF].

As explained previously, this method reduces the problem to evolving N traj wavefunctions, and in the N traj → ∞ limit one recovers the exact evolution. However, the number of required trajectories N traj to obtain accurate results is not known in general [START_REF] Daley | Quantum trajectories and open many-body quantum systems[END_REF] and, in the case of weak dissipation, the method can quickly become equivalent to a full integration of the master equation as a greater amount of trajectories are needed to reach convergence. In recent years, there has been a growing interest in the idea that for a certain class of low-entropy systems, a limited number of states, belonging to a so-called "corner" subspace, can efficiently and faithfully represent the density matrix [START_REF] Finazzi | Corner-space renormalization method for driven-dissipative two-dimensional correlated systems[END_REF][START_REF] Bris | Low-rank numerical approximations for high-dimensional lindblad equations[END_REF][START_REF] Chen | Low-rank density-matrix evolution for noisy quantum circuits[END_REF][START_REF] Mccaul | Fast computation of dissipative quantum systems with ensemble rank truncation[END_REF]. As presented earlier, such a method, the corner-space renormalization method, exists and has proved to be extremely useful in the context of dissipative phase transitions. However, this method is restricted to finding steady-states. An alternative method is therefore needed to simulate the time-evolution of low-entropy strongly-correlated open quantum systems.

In the field of quantum information, the tremendous advances on the control of artificial quantum systems, such as superconducting Josephson qubits [START_REF] Blais | Quantum information processing and quantum optics with circuit quantum electrodynamics[END_REF] and trapped ions [START_REF] Bruzewicz | Trapped-ion quantum computing: Progress and challenges[END_REF], are allowing dramatic progress towards the realisation of devices for quantum computation [START_REF] Arute | Quantum supremacy using a programmable superconducting processor[END_REF][START_REF] Zhong | Quantum computational advantage using photons[END_REF]. We are currently in the noisy intermediate-scale quantum (NISQ) era [START_REF] Preskill | Quantum Computing in the NISQ era and beyond[END_REF], where error correction is not yet possible due to daunting overheads [START_REF] Devitt | Quantum error correction for beginners[END_REF], but where quantum advantage might be already exploited for applications in quantum chemistry [START_REF] Lanyon | Towards quantum chemistry on a quantum computer[END_REF], optimisation [START_REF] Biamonte | Quantum machine learning[END_REF] and even finance [START_REF] Orús | Quantum computing for finance: Overview and prospects[END_REF]. Scaling up quantum circuits and designing practically efficient quantum correction protocols make it of crucial importance to precisely understand the effects of both incoherent and coherent sources of noise on quantum algorithms [START_REF] Martinis | Qubit metrology for building a fault-tolerant quantum computer[END_REF][START_REF] Harper | Efficient learning of quantum noise[END_REF][START_REF] Deutsch | Harnessing the Power of the Second Quantum Revolution[END_REF][START_REF] Cross | Validating quantum computers using randomized model circuits[END_REF].

To meet these challenges, there is a strong need for accurate numerical simulations of quantum hardware on classical computers [START_REF] Aharonov | Polynomial simulations of decohered quantum computers[END_REF][START_REF] Aaronson | Boson sampling with lost photons[END_REF][START_REF] Harrow | Robustness of quantum gates in the presence of noise[END_REF][START_REF] Koch | Demonstrating NISQ era challenges in algorithm design on IBM's 20 qubit quantum computer[END_REF][START_REF] Larose | Overview and Comparison of Gate Level Quantum Software Platforms[END_REF]. In this perspective, the application of tensor-network methods to quantum-circuit simulation has been shown to be effective to model circuits, although with restrictions on the amount of entanglement building up through the circuit [START_REF] Vidal | Efficient Classical Simulation of Slightly Entangled Quantum Computations[END_REF][START_REF] Zhou | What Limits the Simulation of Quantum Computers?[END_REF][START_REF] Noh | Efficient classical simulation of noisy random quantum circuits in one dimension[END_REF][START_REF] Napp | Efficient classical simulation of random shallow 2D quantum circuits[END_REF][START_REF] Oh | Classical simulation of lossy boson sampling using matrix product operators[END_REF]. In addition, most existing simulators of quantum hardware consider local and digital error models [START_REF] Arute | Quantum supremacy using a programmable superconducting processor[END_REF][START_REF] Devitt | Quantum error correction for beginners[END_REF][START_REF] Jones | QuEST and high performance simulation of quantum computers[END_REF][START_REF] Gottesman | Fault-tolerant quantum computation with constant overhead[END_REF], that consist in extending the quantum circuit model to noisy algorithms by including noise gates applied after each unitary gate. While in close proximity with classical error models, these two approximations do not necessarily hold, especially for highly-entangling circuits [START_REF] Weinstein | Quantum process tomography of the quantum Fourier transform[END_REF], and remain a challenge in quantum error correction [START_REF] Klesse | Quantum error correction in spatially correlated quantum noise[END_REF]. To take into account realistic sources of noise for highly-entangling circuits one should revert to a continuous-time description, where noise is taken into account continuously during the dynamics associated to the quantum gates. If one neglects non-Markovian effects, this can be done in the framework of the Lindblad master equation [START_REF] Breuer | The theory of open quantum systems[END_REF]. However, such a description is numerically expensive; for a chain of L qubits with Hilbert space dimension N = 2 L , a density matrix of size N 2 must be evolved. Since quantum processors are conceived to be weakly dissipative and with low entropy, they belong to the class of systems that have a low-entropy. Stabilized arrays [START_REF] Ma | A dissipatively stabilized Mott insulator of photons[END_REF][START_REF] Lebreuilly | Stabilizing strongly correlated photon fluids with non-Markovian reservoirs[END_REF], cat qubit systems [START_REF] Guillaud | Repetition cat qubits for fault-tolerant quantum computation[END_REF] and quantum hardware with state-of-the-art dissipation rates [START_REF] Arute | Quantum supremacy using a programmable superconducting processor[END_REF][START_REF] Preskill | Quantum computing in the NISQ era and beyond[END_REF] also belong to this category. In this chapter, we present the dynamical-corner space method that is particularly suited to simulate quantum processors in the NISQ era. We start by presenting the method and the steps that compose the corresponding algorithm in section II. In the next section, various computational details, such as the complexity of the algorithm, and the question of stiffness for time-integration are discussed in section III. Then, in section IV, we present an application of the method to the simulation of a noisy quantum Fourier transform processor. We discuss the applicability of our method to this quantum algorithm, and find scaling laws for the error as a function of the number of qubits, as well as the impact of dissipative channels on the performance of the algorithm as a function of the considered initial state. Finally, digital error models are discussed in more detail in section V, in which key differences between our approach to simulate quantum circuits and such models are discussed.

II The dynamical corner-space method

Let us consider an open quantum system whose dynamics is governed by the Lindblad master equation [START_REF] Breuer | The theory of open quantum systems[END_REF]:

∂ t ρ = Lρ = -i[ Ĥ, ρ] + D i=1 Ĵi ρ Ĵ † i - 1 2 { Ĵ † i Ĵi , ρ} , (3.1)
where Ĥ is the system Hamiltonian acting on a Hilbert space H of dimension N , and Ĵi is the ith jump operator. At any time t, the solution ρ may be approximated by

ρ(t) ≃ M (t) k=1 p k (t) |ϕ k (t)⟩⟨ϕ k (t)| , p k (t) ≥ p k+1 (t), ∀k, (3.2)
where p k (t) are the M (t) largest eigenvalues of ρ at the time t and |ϕ k (t)⟩ their associated eigenvectors. By construction, the controlled truncation error introduced by such an approximation is monotonically decreasing with M and quantified by ϵ M = 1 -M k=1 p k , so that the decomposition becomes exact for M (t) = r(t), with r(t) denoting the rank of ρ(t), equivalent to the α = 0 Rényi entropy [START_REF] Beck | Thermodynamics of Chaotic Systems: An Introduction[END_REF]. Therefore, in a wide class of lowentropy systems including most platforms relevant for quantum computing, ρ is very well approximated by M ≪ N basis vectors, and even by M ≳ 1 for close to pure states.

M will be referred to as the corner dimension. The accuracy of the calculations will be controlled by a fixed maximum error ϵ with ϵ M ≤ ϵ enforced at any time.

It follows from the above that all the information of the density matrix is carried by a set of weighted corner base vectors of the form √ p k |ϕ k ⟩. In some arbitrary computational basis {|n⟩} N n=1 , these can be represented by a N × M matrix with elements C nk (t) = p k (t) ⟨n|ϕ k (t)⟩. We have:

ρ(t) = M (t) k=1 p k (t) |ϕ k (t)⟩⟨ϕ k (t)| = Ĉ(t) Ĉ † (t).
(3.

3)

The idea behind the dynamical corner-space method is then to evolve Ĉ in time, without ever reconstructing ρ. The evolution Ĉ(t) → Ĉ(t + δt) over a small time step δt, schematically illustrated in Fig. 3.1, involves two computational operations: (i) the construction of the transition basis and (ii) the projection onto new principal components. Below we detail these two steps, respectively represented as the left and right parts of the figure .  Step (i): The weighted corner basis Ĉ(t) evolves into the transition basis T (t + δt) as

ρ(t + δt) = D i=0 Mi ρ(t) M † i = M (D+1) m=1 |ψ m (t + δt)⟩⟨ψ m (t + δt)| = T (t + δt) T † (t + δt). (3.4)
Here, the following Kraus map [START_REF] Wiseman | Quantum Measurement and Control[END_REF], equivalent to Eq. (3.1), is used to time-evolve the density matrix:

e δtL ρ = D i=0 Mi ρ M † i , M0 = exp(-iδt H), Mi≥1 = √ δt Li , ( 3.5
)

with H = Ĥ -i 2 D i=1
L † i Li a non-Hermitian operator depending on the Hamiltonian and the quantum jump operators, which is reminiscent of the stochastic Schrödinger equation presented in section III.5. By construction, the transition basis T (t+δt) is a N ×[M (t)(D+ 1)] rectangular matrix, where D is the number of dissipation channels. Its mth column is given by |ψ m (t + δt)⟩ = √ p µ Mν |ϕ µ (t)⟩, with ν = (m -1) ÷ M (t) and µ = (m -1) mod M (t) + 1. In the quantum trajectories picture, a Kraus operator is drawn at random at every time step according to the probability for an external observer to measure its outcome. This generated single trajectories conditioned by a specific record of the history of the outcomes, which by repeating the procedure could be averaged to recover a correct density matrix. We here employ a completely deterministic approach: at every time step, all the (D + 1) possible trajectories branching from each of the corner's M base state vectors are generated, naturally weighted by their likelihood. Although exact, no further processing would result in an exponential growth of the corner dimension with time, growing as D n δt , with n δt the number of timesteps. The second step of the algorithm solves this problem.

Step (ii): The transition matrix is now projected to a new weighted corner basis Ĉ(t + δt) of lower dimension M (t + δt) ≤ M (t)(D + 1) via a new truncated eigendecomposition P of the form of Eq. (3.3). Crucially, this is possible without ever reconstructing the full density matrix. Indeed, the N × N matrix ρ(t + δt) = T (t + δt) T † (t + δt) and the much smaller [M (t)(D +1)]×[M (t)(D +1)] matrix σ(t+δt) = T † (t+δt) T (t+δt) share the same non-zero eigenvalues p k . The associated eigenvectors, |ϕ k,ρ (t + δt)⟩ and |ϕ k,σ (t + δt)⟩ are related by the identity [START_REF] Chen | Low-rank density-matrix evolution for noisy quantum circuits[END_REF][START_REF] Gentle | Matrix Algebra: Theory, Computations, and Applications in Statistics[END_REF]:

√ p k |ϕ k,ρ (t + δt)⟩ = T (t + δt) |ϕ k,σ (t + δt)⟩ . (3.6)
The components of the new decomposition can then be judiciously truncated to only retain the leading M (t + δt) eigenvalues p k , yielding an updated weighted corner basis Ĉ(t + δt), with the same structure as the initial one Ĉ(t). This procedure can be iterated for an arbitrary time with the possibility of tuning the length of the time step δt and the tolerance parameter ϵ to control the accuracy at the desired level. The time evolution of the corner basis involves no explicit reconstruction of the full density matrix. Indeed, the largest representation of the state of the system involved in the process, the transition basis T , is comparable in size to that of M (D + 1) closed systems. As we will see, this method is comparable to MCWF in terms of complexity, where the the entropy is the essential quantity to check for: if the required corner dimension M is too large, the method fails terribly. For systems with moderate entropy instead, it can be very efficient.

III Computational details III.1 Complexity

The method presented above amounts to evolving M (t) closed systems, with M (t) the corner dimension. For a quantum computation, the initial state is pure, M (t = 0) = 1. As shown in Fig. 3.1(c), the dimension M (t) grows moderately in time. At every time step, the most demanding operation is related to the construction of the matrix σ = T † T . This involves a number of operations of order O(M 2 [D + 1] 2 2 L ). Indeed, the transition matrix T has a size M (D + 1) × N , where N = 2 L denotes the dimension of the Hilbert space and L the size of the system under consideration. Note that the number of jump operators D typically scales as L for most relevant quantum computing platforms consisting of L coupled units. The complexity of the method is thus of order O(M 2 L 2 2 L ). This represents an exponential reduction of the complexity with respect to that of a brute-force master equation integration, which is of order O(4 L ).

III.2 Memory use

The Kraus operators Ki are N ×N matrices, the same size as the density matrix. However, their memory requirement is much smaller, thanks to their extreme sparsity for most. Indeed, jump operators corresponding to dissipative processes are typically single-body operators and thus as memory-consuming as state vectors1 . The largest Kraus operator One can see that the gain in memory requirement from our method scales exponentially. Note tha the Hamiltonian, being very sparse, is not the limitation. The largest stored object is the corner Ĉ, that is composed of M states. Here, as an example, we considered M = 100. (b) Infidelity between the dissipative output density matrix ρ(c) γ (calculated via the corner-space method) and the exact integration density matrix ργ as a function of the maximum corner dimension M max , for an L-qubit circuit and γT QFT = 1.42 × 10 -2 . The initial state is

|ψ 0 ⟩ = QFT -1 |GHZ⟩.
is K0 , which roughly corresponds to Ĥ. The RAM needed to store it is shown in Fig. 3.2 together with that associated to other numerically relevant objects as a function of the number of qubits. Another feature of our implementation is that the dimension of the corner basis M (t) is dynamically adapted to match the maximum allowed error that we impose, thus optimizing the computing and memory resources. This is done by means of memory-contiguous dynamically resizable arrays. This saves considerable computation time, as M (t) grows in time when starting from a pure state (M (0) = 1).

III.3 Efficient evaluation of relevant observables

The evaluation of relevant metrics in quantum information is challenging with quantumtrajectory approaches, such as the Monte Carlo wave function method [START_REF] Carmichael | An open systems approach to quantum optics[END_REF][START_REF] Mølmer | Monte-Carlo wave-function method in quantum optics[END_REF][START_REF] Plenio | The quantum-jump approach to dissipative dynamics in quantum optics[END_REF] (MCWF). Such an approach relies on the evolution of n stochastic trajectories {|ψ i (t)⟩} n i=1 . The density matrix can then be reconstructed as ρ

(t) = (1/n) n i=1 |ψ i (t)⟩⟨ψ i (t)|.
To simulate a quantum computation, this has two downsides: firstly, for weakly dissipative systems (γT ≪ 1), most trajectories will experience no quantum jump and thus be identical. This means that most of the computing resources are wasted in performing a redundant task. Secondly, most relevant metrics involved in quantum-information processing, such as fidelity and entanglement measures, namely concurrence, negativity, or entanglement entropy [START_REF] Horodecki | Quantum entanglement[END_REF], require constructing explicitly the (dense) density matrix of the system and diagonalizing it. In practice, the latter, of complexity O(N 3 ), is not feasible for systems larger than ∼ 15 qubits. In contrast, the presented method yields explicitly both the eigenvalues {p k (t)} M k=1 and the eigenvectors {|ϕ k (t)⟩} M k=1 at every time step, with no need for additional calculations. To give a concrete example, let us consider the evaluation of the fidelity between two arbitrary mixed states ρ and ρ′ with rank M and M ′ respectively, as given by

F(ρ, ρ′ ) = Tr[ ρρ ′ ρ] = M m=1 ϕ m M k,m=1 |ϕ k ⟩ A km ⟨ϕ m | 1/2 ϕ m , with A km = M ′ k ′ =1 p ′ k ′ √ p k p m ⟨ϕ k |ϕ ′ k ′ ⟩ ⟨ϕ ′ k ′ |ϕ m ⟩ , ( 3.7)
where p

(′)
k and ϕ

(′)
k correspond to the kth eigenvalue and eigenvector of ρ(′) . It appears from the above that the total complexity of this evaluation is of order

O(2N 3 ) + O(2M ′ N ) + O(M 3 ) + O(M N ) , ( 3.8)
where the first term accounts for the eigendecomposition of the two density matrices, the second for the construction of the matrix A km , the third for its diagonalization and the last for the trace. An additional subleading complexity of order O(n×N ) is to be added if the density matrices are obtained via Monte Carlo wave function calculations in order to account for their construction. In contrast, the order O(2N 3 ) is to be discarded when using the dynamical corner-space method, as the eigendecompositions are known explicitly.

Then, for each method, one finally has to leading order and for M (′) ≪ N the following scaling figures of merit: In practice, the inconvenient scaling of the complexity for the MCWF Dynamical corner-space

F(ρ, ρ′ ) O(2N 3 ) O(max(M, 2M ′ )N ) F(ρ, |ϕ⟩) O(nN ) O(M N ) S(ρ) O(N 3 ) O(M )
Monte Carlo wave function approach, stemming from the two density-matrix diagonalizations, combined with the necessity of storing dense matrices well beyond the realistically available RAM makes it impossible to compute the fidelity between two mixed states from trajectories for systems larger than ∼ 13 sites on a desktop computer. A similar discussion can be made for the entropy. We here see the superiority of the dynamical corner-space method, in particular in the context of quantum information where such metrics are often calculated.

III.4 Scaling of the corner dimension

In our method, the tolerance on the precision of the density matrix is tunable by design and the convergence versus M is ensured for moderately dissipative systems with low entropy. This is contrast to the case of quantum trajectories for which the number of needed trajectories n for a given problem is currently unknown a priori [START_REF] Daley | Quantum trajectories and open many-body quantum systems[END_REF]. In Fig. 3.2(b), one can see that for γT QFT = 1.42 × 10 -2 , an infidelity below 10 -3 can be reached by choosing M max ∼ L ln L. Hence, the corner dimension grows only polynomially with the size L, in particular sub-quadratically. Note that this value of γT QFT is slightly higher than the state-of-the-art, hence the method is well suited to simulate experimental platforms in the coming years.

III.5 Integration method and stiffness

An issue with the map Ĉ(t) → Ĉ(t + δt) that updates the corner basis at each time step is that it involves terms proportional to √ δt in the Kraus operators. This seems to restrict the method to a first order explicit integration scheme, which could result in a poor stability of the method when dealing with stiff dynamics ensued from nonlinear Hamiltonians. However, we circumvent this limitation by separating the pseudo-unitary evolution generated by K0 from that of the rest of the Kraus operators. Instead of computing T im = k K 0,ik C km = k (δ ik -iδt Hik )C km , we perform an exact numerical integration of the differential equation ∂ t Ĉ = -i H Ĉ over the time interval [t, t + δt] via an ordinary differential equation (ODE) solver well-adapted to the level of stiffness induced by the Hamiltonian. This allows us to treat stiff problems, via adapted implicit methods, and to use adaptive time-stepping.

To illustrate this point, let us consider a Kerr-cat qubit [START_REF] Grimm | Stabilization and operation of a kerr-cat qubit[END_REF] as described by the following Hamiltonian:

Ĥ = Kâ † â † ââ + ω c â † â + G(â 2 + â †2 ) , ( 3.9)
where â (â † ) is a cavity mode annihilation (creation) operator, ω c is the frequency of the cavity, K its Kerr nonlinearity, and G the two-photon driving frequency. If one considers jump operators Ĵ1 = √ γâ, Ĵ2 = √ κâ 2 that respectively describe single-and two-photon losses, a logical qubit can be conceived by considering logical states

|+⟩ = |C + α ⟩ = (|α⟩ + i |-α⟩)/ √ 2 and |-⟩ = |C - α ⟩ = (|α⟩ -i |-α⟩)/ √
2 as these are steady states of the system in the γ → 0 limit. Numerically, the simulation of such systems cannot be efficiently treated with explicit first-order methods as the differential equation corresponding to their time evolution is stiff because of the Kerr nonlinearity. Thanks to the numerical integrator used for the coherent part of the evolution of the corner Ĉ, our method is capable of describing such systems, whose entropy is limited when κ ≫ γ. In Fig. 3.3, the Wigner functions of each of the 4 most populated states of the corner are shown with their associated probabilities p i , after having evolved the system for a time t = 10/γ via the corner-space method, setting a photon cutoff n ph = 20. One sees that the low-dimensional basis found by the corner indeed closely matches that of a qubit with Schrödinger-cat logical states while keeping track of the effects of the single-body loss on the lowly probable states. By tuning the tolerance parameter of the method ϵ, such dissipative effects can be captured to any desired order. Our method therefore enables one to investigate the dynamics of such systems in an efficient way, and could be used to understand how single-photon dissipation impacts quantum operations in multi-Kerr-catqubit systems, among other applications involving bosonic qubits.

|x 1 H R 2 R L-1 R L |y 2 |x 2 H R L-2 R L-1 |y 2 |x L-1 H R 2 |y L-1

IV Application to the noisy QFT IV.1 The quantum Fourier transform

We will now apply the method to the numerical simulation of noisy quantum circuits.

A quantum circuit consists in a set of gates G, that correspond to unitary operators Ûg , with g ∈ G. These operators are applied to a register of L qubits that have as input state |x⟩, that are equivalent to the random-access memory of a classical computer. The output of a quantum circuit is given by state |y⟩ = f (|x⟩), where f denotes any recursive function [START_REF] Deutsch | Quantum theory, the Church-Turing principle and the universal quantum computer[END_REF], and |y⟩ corresponds to the answer we want from the algorithm. We consider the noisy quantum Fourier transform (QFT) circuit, that is the cornerstone of the Shor algorithm and many algorithms related to the hidden subgroup problem [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF].

In doing so, we will consider the architecture depicted in Fig. 3.4 (a). The input state to be Fourier transformed through the circuit is encoded in the state of L qubits. Given an orthonormal computing basis {|n⟩} N n=0 spanning the Hilbert space of the system, of dimension N = 2 L , this circuit linearly transforms any input of the form

|x⟩ = N -1 j=0 x j |j⟩ , (3.10)
where x is a vector of complex amplitudes, into

|y⟩ = QFT(|x⟩) = N -1 k=0 y k |k⟩ , with y k = 1 √ N j e 2iπjk/N x j . ( 3.11)
This therefore corresponds to the discrete Fourier transform of the N wave-function amplitudes. This operation can be reduced to subsequent applications of only two different types of quantum gates [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF], the Hadamard and the controlled-phase gates. In the local basis:

ÛH | i = 1 √ 2    1 1 1 -1    , ÛR | jk =            1 1 1 e 2πi/2 k            . ( 3
.12)

⟨0 i | ⟨1 i | |0 i ⟩ |1 i ⟩ ⟨0 j 0 k | ⟨0 j 1 k | ⟨1 j 0 k | ⟨1 j 1 k | |0 j 0 k ⟩ |0 j 1 k ⟩ |1 j 0 k ⟩ |1 j 1 k ⟩
The number of involved such elementary gates, L(L + 1)/2, is only polynomial in L = log 2 (N ); the complexity of the QFT algorithm is hence of order O(L 2 ). In contrast, that of the usual fast Fourier transform (FFT) is exponential: O(N log N ) ∼ O(L2 L ). This quantum advantage makes the QFT a central building block of many other algorithms.

In the following, the quantum gates introduced above will be executed via a continuoustime evolution defined by an appropriate master equation taking the form of Eq. (3.1). The effect of two types of noise will be addressed. We will first consider dissipative processes induced by a weak coupling to a zero-temperature environment, as described by jump operators of the form:

Ĵi = √ γ σ- i , ( 3.13)
ubiquitous in superconducting platforms. This corresponds to local decay processes from the excited qubit state |↑⟩ j ≡ |0⟩ to the lower energy qubit state |↓⟩ j ≡ |1⟩. We will treat the effect of decoherence as well. This will be described by local pure-dephasing jump operators:

Ĵi = √ γ σz i , ( 3.14)
whose microscopic origin can be related to different sources such as disorder and phase fluctuations in the drive [START_REF] Marquardt | Introduction to dissipation and decoherence in quantum systems[END_REF]. We consider a drive rotating at the qubit Rabi frequency. The Hadamard gate acting on the ith qubit can be realized via the subsequent application of the two following Hamiltonians:

Ĥ1 i = δ 2 σy i , Ĥ2 i = δ 2 σz i , (3.15)
for a time ∆t H 1 = π/2δ and ∆t H 2 = π/δ, respectively. The controlled-phase gates with control qubit j and target qubit k can be performed through the Hamiltonian:

Ĥj,k = δ 2 σz j + δ 2 σz k - δ 2 σz j σz k + 1 , ( 3.16)
applied for a time ∆t R,k = π/δ2 k . For simplicity, we will assume sudden switching between gate Hamiltonians and will not include coherent errors (related to some switching time,or to some unwanted detuning between the drive and the qubit frequency), although both effects could be accurately described by the dynamical corner-space method.

A first example of the dynamics of such a circuit is presented in Fig. 3.4 (b) and Fig. 3.4 (c). There, the initial state is chosen as

|ψ 0 ⟩ = QFT -1 (|GHZ⟩) = 1 √ 2N N -1 n=0
(1 + e 2iπn/N ) |n⟩ , (3.17) so as to be the inverse of the Greenberger-Horne-Zeilinger (GHZ) state [START_REF] Greenberger | Going Beyond Bell's Theorem[END_REF][START_REF] Bouwmeester | Observation of Three-Photon Greenberger-Horne-Zeilinger Entanglement[END_REF][START_REF] Pan | Experimental test of quantum nonlocality in three-photon Greenberger-Horne-Zeilinger entanglement[END_REF]:

|GHZ⟩ = 1 √ 2 (|00...0⟩ + |11...1⟩) . (3.18)
The latter is known to be a maximally entangled state. Through this choice, the output state is ensured to be highly entangled, therefore demonstrating that entanglement is not a limiting factor for the corner-space method, as the intermediate states build up long-range entanglement. In panel (b) of Fig. that evaluates entanglement between bipartitions of the form {{1, . . . , n}, {n+1, . . . , L}}. While this is a rigorous measure of entanglement only for pure states [START_REF] Horodecki | Quantum entanglement[END_REF], it still gives a valid qualitative description of the entanglement temporal build-up for states that are γ calculated with the corner-space method and the output ργ obtained via an exact integration as a function of the maximum corner dimension M max for different values of γT QFT and a fixed number of qubits L = 10. (b) Fidelity as a function of M max and L for γT QFT = 2.5 × 10 -2 . The parameter N denotes the dimension of the total Hilbert space. Values corresponding to M max = L are highlighted by hollow markers (we find a fidelity F ≳ 0.997 for M max ∼ L ln L). The initial state is |ψ 0 ⟩ as in Fig. 3.1 for all three panels. (c) Simulation time of the QFT algorithm in the presence of dissipation versus the number L of qubits for the exact solution of the master equation (squared markers) and for the time-dependent corner-space method (circles) for two different values of the control parameter ϵ. The dissipation rate is set by γT QFT = 2.5 × 10 -2 . close to pure. Note that S ent (n) > S is a sufficient condition 2 for showing quantum entanglement between the two subsystems 1 ⊕ . . . ⊕ n and n + 1 ⊕ . . . ⊕ L. One sees that the initial localized entanglement spreads through the system in close relation to the architecture of the circuit. As clearly follows from the shading in panel (b), this takes place when the (entangling) controlled-phase gates are applied.

IV.2 Benchmarking the method

Let us benchmark the accuracy of our calculations against the results of an exact integration of the master equation for small values of L, the numbers of qubits. In what follows, ργ and ρ(c) γ denote the output density matrices of the noisy QFT obtained via the exact integration and via the corner method, respectively. Instead, ρ0 denotes the ideal outcome of the noiseless QFT, which is a pure state. The results of this benchmarking process are presented in Fig. 3.5 for fixed values of γT QFT , where T QFT denotes the physical duration of the QFT operation. This ensures that the output infidelity with respect to ρ0 remains constant as the circuit size is increased. In particular, Fig. 3.5 (a) shows the infidelity of the method 1 -F(ρ (c) γ , ργ ) as a function of the rescaled maximum corner dimension M max /L for L = 10. One sees that for M max ∼ 10L the exact results are excellently approximated by the time-dependent corner-space method for all the considered values of γT QFT . The method still performs reasonably well for noise rates as high as γT QFT = 1.5 × 10 -1 , where the fidelity to the output of the noiseless circuit is as low as F(ρ 0 , ργ ) = 0.758. Fig. 3.5 (b) shows the fidelity of the method F(ρ (c) γ , ργ ) for different numbers of qubits L as a function of the fraction M max /N between the dimension of the basis spanning the corner and that of the basis of the entire Hilbert space. These results show that the advantage of the dynamical corner-space method over exact integration of the master equation increases with L. In particular, to match a desired fidelity, the required corner dimension M is found to grow as L ln L with the system size. For L = 12, an excellent agreement of the corner-space method with the exact integration is already obtained for M max /N = 10 -2 . Finally, in Fig. 3.5 (c), we compare the computation time of the corner-space method to the exact integration, for two different values of ϵ3 . The corner-space method achieves an exponential speed-up with respect to the master equation integration. This leads to simulations faster by more than three orders of magnitude for L ∼ 15. Moreover, tuning the tolerance ϵ from 1 × 10 -4 down to 5 × 10 -5 preserves the scaling of the simulation time with L. Simulations of up to 21 qubits are presented. This represents a Hilbert-space dimension (N = 2 L ) above 2 million states, corresponding to (dense) density matrices weighting 64 TiB that could never have been handled with a brute-force integration of the master equation (3.1). The method presented above is thus capable of efficiently evolving such high-dimensional objects by means of parsimonious corner representations of the density matrix. The calculations were performed on a single six-core Intel Xeon processor at a clock speed of 1.9GHz.

IV.3 Scaling laws

One can now evaluate the impact of incoherent processes on intermediate-scale devices via a continuous-time description and determine the scaling of errors. In Fig. 3.6, the fidelity F(ρ (c) γ , ρ0 ) is shown for up to L = 21 qubits, for different values of γ/δ. Here, we consider dissipation channels described by the jump operators Ĵi = √ γ σi . Remarkably, the infidelity scales only quadratically as a function of the number of qubits L. This scaling dependence allows one to precisely estimate the impact of γ = 1/T 1 on the QFT algorithm, T 1 being the energy relaxation time of the considered system.

IV.4 Impact of initial states

Another key property is the dependence of the fidelity on the initial state, crucial to redesign algorithms that rely preferentially on a certain class of states. In Fig. 3.7, we address this question for the QFT by sampling initial states. Either energy relaxation produced by the jump operators Ĵi = √ γ σi or pure dephasing described by Ĵi = √ γ σz i are considered. Only two simple parameters that characterise the initial state are found to be crucial for the considered architecture: the total number n S of spins up and the spin-up "barycenter"

B(ρ) = 1 n S ℓ ℓ × Tr |↑⟩⟨↑| ℓ ρ . (3.20)
Our findings show that, in presence of energy relaxation, the fidelity of the noisy QFT decreases linearly with the number of spins up in the initial state. This is in stark contrast to the case of pure dephasing, which shows no significant dependence on n S . The fidelity also exhibits a strong dependence on the spin barycenter. Indeed, energy relaxation only affects excited states and the circuit's Hadamard gates are applied one qubit at a time starting from the beginning of the chain. As a result, excited qubit states (spin up) close to the end of the chain are rotated down to the Bloch-sphere equator by the Hadamard gates later than those on the opposite end. Thus, they are globally more affected by dissipation.

V Quantum errors

A fundamental question that quickly appeared in the debates on quantum computing in the early 90s is the question of how quantum errors occur. In classical computing, this was a problem for a long time and it took some time for classical error correction to be developed, leading to fault-tolerant classical computers. Quantum error correction is trickier, as for instance quantum information cannot simply be copied because of the no-cloning theorem, so all repetition codes developed for classical computing cannot be reused. Quantum error correction codes are based on digital error models, that are a simplification of how quantum errors microscopically occur. Here we will examine the differences between such a model and a full master equation approach.

As presented in section IV, a convenient and widely used model for quantum computation is the gate-based model [START_REF] Nielsen | Quantum Computation and Quantum Information: 10th Anniversary Edition[END_REF]. For closed systems, this model is strictly equivalent to the successive application of unitary time-evolution operators of the form ÛG = e -i ĤG τ , with ĤG the gate Hamiltonian and τ the gate time. In most current classical simulations of noisy quantum processors, errors are accounted for by extending the gate-based model to what has been coined digital error models [START_REF] Arute | Quantum supremacy using a programmable superconducting processor[END_REF]. These consist in applying noise gates after each unitary gate, expressed as Kraus operators, in analogy to error models for classical processors. However, in general, a quantum system is subject to a completely-positive, trace-preserving (CPTP) map acting on the system density matrix ρ in continuous time.

The generator L of such a map can always be expressed in the Lindblad form as a Liouvillian whose action on the density matrix takes the form

L[ρ] = U[ρ] + D[ρ],
with U and D the unitary and dissipative contributions to the time evolution of ρ. Explicitly, for a quantum gate G, one has:

U[ρ] = -i[ ĤG , ρ], (3.21)
D[ρ] = D i=1 Ĵi ρ Ĵ † i - 1 2 { Ĵ † i Ĵi , ρ} , (3.22)
with Ĵi the jump operators that describe the dissipative channels, which take a simple form when the environment can be treated within the Born-Markov approximation. Given an initial density matrix ρ, after a time interval τ , the density matrix at the output of the gate ρ(τ ) is given by ρ(τ ) = e τ L ρ. One can always separate the unitary and dissipative parts of L to obtain

ρ(τ ) = E G ÛG (τ )ρ Û † G (τ ) ≡ E G (ρ ÛG ), (3.24)
where we have defined the density matrix after the ideal unitary process ρ ÛG = ÛG (τ )ρ Û † G (τ ), with ÛG (τ ) the time evolution operator that corresponds to the application of Hamiltonian ĤG for time τ , and E G a non-unitary map. In cases where D and U commute, one explicitly has ÛG (τ ) = e -i ĤG τ and E G = e τ D , for the unitary and the error processes, respectively. However, in general, obtaining E G is non-trivial, and experimentally one would need to perform a tomography on each gate to determine its exact error process: this is known as quantum process tomography [START_REF] Mohseni | Quantum-process tomography: resource analysis of different strategies[END_REF].

Here we numerically simulate the quantum process tomography of a controlled phase gate of duration τ = π/2δ, δ being the Rabi frequency of the qubits. To do so, we decompose the error process E in the Pauli basis as:

ρ(τ ) = E(ρ CR π/2 ) = m,n χ err mn Pm ρCR π/2 P † n , ( 3.25)
with ρCR π/2 the output density matrix corresponding to an ideal controlled phase gate, Pn ∈ {I, X, Y, Z} ⊗2 the generators of the Pauli group on 2 qubits and χ err the error matrix that completely characterizes the error process and that we aim at numerically determining.

In Fig. 3.8, the error matrix corresponding to a noisy CR π/2 gate is shown for two different cases. Panel (b) corresponds to a Lindblad evolution with jump operators Ĵi = √ γ σi , and panel (a) corresponds to a Lindblad evolution with jump operators Ĵi = √ γ σz i . By using a Choi decomposition [START_REF] O'brien | Quantum process tomography of a controlled-not gate[END_REF], one can numerically estimate the error matrix that recovers the density matrix ρ(τ ) of the realistic Lindblad evolution. In the first case, Ĵi = √ γ σi was considered, which does not commute with the Hamiltonian. This leads to a complex error process accounting for the spatial propagation of local errors upon (non-local) Hamiltonian time evolution. In this case, this manifests through the presence of single-and two-qubit error events with different magnitudes in the components of the error matrix. Similar error matrices have also been found experimentally [START_REF] Weinstein | Quantum process tomography of the quantum Fourier transform[END_REF]. Each non-negligible element corresponds to a Kraus operator to be applied after the ideal gate in the digital-error approach. While two local jump operators suffice to describe the qubit-environment interaction, up to 16 Kraus operators could be necessary in such an approach. Note that the noise process acting on two qubits at the same time generally has a smaller magnitude, which can justify that one can neglect them for low-depth circuits (the quantum Fourier transform has L(L + 1)/2 gates, meaning these errors can accumulate in the long run). By contrast, in panel (b) pure-dephasing jump operators were considered, Ĵi = √ γ σz i , which commute with the controlled phase Hamiltonian. Hence, the error process only involves Z errors in this case, although not strictly local (notice the ZZ component).

As appears from the above discussion, modeling errors in continuous time is of fundamental and applied importance for the following reasons:

(1) The errors induced by the presence of local dissipative events, as captured by local jump operators, cannot be accounted for via local noise gates, and are affected by the applied Hamiltonian, an effect that is generally not described by digital error models. For d-qubit gates, the number of noise gates to be applied goes up to 4 d in a digital-error approach.

(2) For circuits such as the quantum Fourier transform, there are L controlled rotation gates that are applied for different times, hence one would need to perform tomography on each of them to recover their error matrices (as well as for the Hadamard gate). For multi-qubit gates such as the Toffoli gate, tomography becomes even more expensive.

(3) Characterizing the jump operators of an experimental platform is much easier than performing tomography to obtain the error process for each gate. For example, in the case of superconducting circuits, measuring the T 1 and T 2 relaxation times is enough.

(4) Our method is also able to treat collective dissipative processes described by jump operators such as Ĵz = (1/ √ L) i σz i [START_REF] Nissen | Collective suppression of linewidths in circuit QED[END_REF][START_REF] Shammah | Open quantum systems with local and collective incoherent processes: efficient numerical simulations using permutational invariance[END_REF]. This corresponds to a single Kraus operator in our method, and is hence inexpensive. With a digital error model, describing collective effects would require to perform tomography on a L-qubit system, which is intractable since it requires O(4 L ) measurements.

VI Conclusion

We investigated the role of dissipation and decoherence in noisy intermediate-scale quantum circuits. Focusing on the key algorithm of the QFT, we revealed the scaling behavior for the fidelity and explored its dependence on the initial state. To achieve this goal, we have introduced and demonstrated a numerical time-dependent corner-space method that performs a judicious compression of the Hilbert space to faithfully represent the system density matrix. The method is not limited by entanglement and is suitable for systems with moderate entropy. Furthermore, our approach could be combined with efficient representations of the corner-space wavefunctions, such as neural-network ansätze [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF][START_REF] Carleo | Constructing exact representations of quantum many-body systems with deep neural networks[END_REF][START_REF] Schmitt | Quantum Many-Body Dynamics in Two Dimensions with Artificial Neural Networks[END_REF][START_REF] Sharir | Deep Autoregressive Models for the Efficient Variational Simulation of Many-Body Quantum Systems[END_REF]. These qualities make our approach ideally tailored for the NISQ era, where the noise processes will become weaker as progress is made. This provides a new tool to improve our understanding of quantum hardware. The presented method can indeed be applied in many contexts related to quantum information: algorithm design for quantum feedback [START_REF] Fösel | Reinforcement learning with neural networks for quantum feedback[END_REF], machine learning for quantum control [START_REF] Zeng | Quantum control based on machine learning in an open quantum system[END_REF] and quantum error mitigation [START_REF] Sun | Mitigating Realistic Noise in Practical Noisy Intermediate-Scale Quantum Devices[END_REF].

Neural-network quantum states

In section 1, we have seen that there exists a plethora of methods to simulate quantum systems, both in closed and open settings. Each methods has its strengths and weaknesses; for instance mean-field methods are computationally extremely efficient while they poorly represent entanglement. Tensor networks slightly improve on this, being more resourceintensive but accurately representing area-law entangled states. Finally, corner methods have less limitations on entanglement, but deal poorly with high-entropy quantum systems. In this chapter, we present neural quantum states (NQS), a new class of methods that enable to treat systems in the strongly correlated regime, with a higher representative power than tensor networks. These methods are based on machine learning algorithms, more precisely on generative models.

I Introduction

Since LeCun's revolutionary paper on image classification in the 1990s [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], the broad field of machine learning has exploded. It can perhaps be considered as the latest technological revolution, nowadays permeating every field of science and having immense commercial impact. Its impact on physics has already been considerable [START_REF] Carleo | Machine learning and the physical sciences[END_REF]. In the braodest sense, machine learning consists in building algorithms that are not symbolically written: for instance, in the AlphaGo program [START_REF] Silver | Mastering the game of go with deep neural networks and tree search[END_REF], the instructions of how to play Go are not explicitly programmed, only a general purpose architecture is built-in, and it is only through some training procedure that the performance of the algorithm can become satisfying.

Generative models have recently been at the forefront of machine learning research. These models enable to reconsitute objects of a certain class from having observed a large number of instances of objects of this class. For instance, after observing a large number of human faces, new realistic faces can be constructed [START_REF] Brock | Large Scale GAN Training for High Fidelity Natural Image Synthesis[END_REF]. As such, one can claim the function corresponding to generating a realistic face will have been learned. The same idea has been applied to language processing, where models are now able to complete sentences in a high-fidelity and realistic fashion, thus being able to write, e.g., newspaper columns [START_REF] Brown | Language Models are Few-Shot Learners[END_REF].

It is therefore a somewhat natural question to ask if such powerful models could help for our understanding of complex many-body systems. Wavefunctions are indeed highdimensional objects that must carry some structure, much like images and language. In this context, neural quantum states (NQS), a variant of generative models, that consist in finding faithful low-dimensional representations of many-body wavefunctions by a neural network were introduced [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF]. Many important applications of NQS have been explored, including finding the ground state [START_REF] Sharir | Deep autoregressive models for the efficient variational simulation of many-body quantum systems[END_REF][START_REF] Choo | Two-dimensional frustrated J 1 -J 2 model studied with neural network quantum states[END_REF], investigating the dynamics of strongly correlated systems [START_REF] Schmitt | Quantum many-body dynamics in two dimensions with artificial neural networks[END_REF][START_REF] Gutiérrez | Real time evolution with neural-network quantum states[END_REF], quantum tomography [START_REF] Torlai | Neuralnetwork quantum state tomography[END_REF], open quantum systems [START_REF] Vicentini | Variational neural-network ansatz for steady states in open quantum systems[END_REF][START_REF] Hartmann | Neural-network approach to dissipative quantum many-body dynamics[END_REF][START_REF] Yoshioka | Constructing neural stationary states for open quantum many-body systems[END_REF][START_REF] Nagy | Variational quantum Monte Carlo method with a neuralnetwork ansatz for open quantum systems[END_REF] and classical simulation of quantum circuits [START_REF] Jónsson | Neural-network states for the classical simulation of quantum computing[END_REF][START_REF] Medvidović | Classical variational simulation of the quantum approximate optimization algorithm[END_REF]. In fact, NQS have proved to be the most accurate variational method in approximating the J 1 -J 2 model's ground state in the frustrated regime [START_REF] Choo | Two-dimensional frustrated J 1 -J 2 model studied with neural network quantum states[END_REF]. In addition, several works have recently demonstrated the superior capacity of some neural-network architectures over tensor-network states in representing volume-law entangled states [START_REF] Levine | Quantum entanglement in deep learning architectures[END_REF][START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF] or area-law entangled states in 2D [START_REF] Wu | From Tensor Network Quantum States to Tensorial Recurrent Neural Networks[END_REF]. All these results indicate that NQS will become a standard tool to numerically study many-body systems in the coming years. In this chapter, we introduce neural quantum states as well as applications to simulating dynamics of closed and open quantum systems. In section II, we introduce the task of generative modeling, and how its building blocks can be applied to learning quantum states. In section III, we present present Monte-Carlo Markov chains and their application to NQS. Then, we show how many-body wavefunctions are represented in practice in section IV, and present a few of the main architectures for NQS that were introduced in recent years. The question of their expressive power, or capacity, with respect to other approximation methods will be examined in section V. We then turn to results obtained for the dynamics of many-body systems using NQS in section VI. We present some of the issues that can occur when trying to propagate an NQS in time, and provide an alternative solution for these problems. As such, we are able to obtain results for the Kibble-Zurek mechanism in large two-dimensional lattices, and can compare these results to state-of-the-art numerical methods. Finally, in section VII, we present the application of NQS to open quantum systems. We explain why it is general trickier to construct an NQS for a density matrix, and present some of the solutions to these issues that have been developed recently. One of these solutions, that will work well for low-entropy systems, is to use a low-rank representation of the density matrix and to represent each state of its corner basis by an NQS.

II Generative models

Let us start by introducing the basic ingredients of generative modelling. A generative model is a statistical model f θ , that depends on a set of variational parameters θ, of a function f that operates on elements x of a configuration space X . For simplicity, we will consider a discrete configuration space in what follows, with N the total number of configurations. As such, the goal of generative modelling is to find the statistical model f θ that approximates the function f as best as possible. This means that at the end of the execution of some training algorithm

D (f θ , f ) = ϵ, (4.1)
with ϵ → 0 if the number of training steps goes to infinity with a carefully selected distance measure D, provided the generative model f θ has enough representative power1 . For instance, a commonly used distance measure in machine learning is the mean square error [START_REF] Goodfellow | Deep Learning[END_REF]:

D (f θ , f ) = x N x i (f θ (x i ) -f (x i )) 2 (4.2)
which we see approaches 0 as f θ → f . To minimize this, one can use many different algorithms, that are all variants of gradient descent. Gradient descent consists in calculating the gradient of the distance, and updating the variational parameters so as to follow the steepest descent of the distance in the variational space. The update rule for the parameters θ is given by:

θ n+1 = θ n -η∇ θ D(f, f θ ), (4.3)
where η, known as the learning rate, is a hyperparameter that is adjusted for better convergence, and n denotes the nth iteration of the training procedure.

In practice, the number N of possible configurations may be exponentially large, hence the sum in Eq. ( 4.2) can become intractable. One must therefore approximate the distance D as well as its gradient. We can therefore identify two major requirements to perform generative modelling:

1. The evaluation of f θ (x) must be tractable.

2. The sums appearing in D(f, f θ ) and ∇ θ D(f, f θ ) must be tractable or made tractable by approximation.

Before presenting how these requirements are fulfilled in practice, let us specialize the discussion to quantum states. In the general case, one wishes to approximate a certain quantum state |ψ⟩. The generative model, also known as neural quantum state (NQS) will therefore be denoted by |ψ θ ⟩ in the rest of the chapter. Throughout the chapter we will restrict the discussion to spin systems, where spin configurations s = (σ 1 , σ 2 , . . . , σ L ) are associated to states |s⟩ that belong to a Hilbert space H. |ψ θ ⟩ and |ψ⟩ can be expanded in the basis spanned by these states as:

|ψ⟩ = s ψ(s) |s⟩ (4.4) |ψ θ ⟩ = s ψ θ (s) |s⟩ . (4.5)
The functions ψ and ψ θ operate on vectors s and the resulting amplitudes ψ(s) and ψ θ (s) are complex numbers. Thus the generative model is fully specified by its amplitudes ψ θ (s), that form a 2 L vector of complex numbers. Note that the requirement that ψ θ (s) be tractable means that the number of variational parameters θ must be polynomial in the system size, or else computing amplitudes will include a super-polynomial number of operations.

Originally, these ideas were applied to finding the ground state of a many-body Hamiltonian Ĥ. In this case, the cost function to minimize is the energy, given by

E θ = ⟨ψ θ | Ĥ|ψ θ ⟩ ⟨ψ θ |ψ θ ⟩ . (4.6)
As the energy converges to its minimum, |ψ θ ⟩ → |ψ 0 ⟩, with the ground state defined by Ĥ |ψ 0 ⟩ = E 0 |ψ 0 ⟩. By doing this one effectively learns the ground state. The expression for the expectation value of any operator  can be written in the following form, for a given NQS |ψ θ ⟩:

⟨ Â⟩ = ⟨ψ θ | Â|ψ θ ⟩ ⟨ψ θ |ψ θ ⟩ (4.7) = s,s ′ ⟨ψ θ |s⟩⟨s| Â|s ′ ⟩⟨s ′ |ψ θ ⟩ s |⟨ψ θ |s⟩| 2 (4.8) = s ⟨ψ θ |s⟩ ⟨s|ψ θ ⟩ ⟨s|ψ θ ⟩ s ′ ⟨s| Â|s ′ ⟩⟨s ′ |ψ θ ⟩ s |⟨ψ θ |s⟩| 2 (4.9) = s |⟨ψ θ |s⟩| 2 s ′ ⟨s| Â|s ′ ⟩ ⟨s ′ |ψ θ ⟩ ⟨s|ψ θ ⟩ s |⟨ψ θ |s⟩| 2 . ( 4.10)
Two main terms can be identified:

P (s) = |⟨ψ θ |s⟩| 2 s |⟨ψ θ |s⟩| 2 (4.11) A loc (s) = s ′ ⟨s| Â|s ′ ⟩ ⟨s ′ |ψ θ ⟩ ⟨s|ψ θ ⟩ . (4.12)
where A loc (s) is the so-called local estimator of Â. Therefore, the expectation value of an observable  can be written as the expectation value of its local estimator A loc over the probability distribution P (s):

⟨ Â⟩ = s P (s)A loc (s) = ⟨A loc ⟩ P . (4.13)
Let us stress that these calculations are feasible for operators with the property that the number of states s ′ such that the elements |⟨s| Â|s ′ ⟩| ̸ = 0, for arbitrary s is at most polynomial in the number of spins. This is known as k-locality, with k the maximum number of sites a single term in the Hamiltonian acts on (for the transverse-field Ising model, k = 2). Conversely, evaluating A loc (s) would not be tractable, given that the sum over s ′ in 4.12 would be over an exponential number of elements. By defining the logarithmic derivatives, O k (s) = ∂ θ k ψ θ (s)/ψ θ (s), the components of the gradient can be calculated straightforwardly:

∂ θ k ⟨ Â⟩ = ∂ θ k s,s ′ ψ * θ (s) ⟨s| Â|s ′ ⟩ ψ θ (s ′ ) s |ψ θ (s)| 2 = s,s ′ ψ * θ (s) ⟨s| Â|s ′ ⟩ O k (s ′ )ψ θ (s ′ ) s |ψ θ (s)| 2 + s,s ′ ψ * θ (s)O * k (s) ⟨s| Â|s ′ ⟩ ψ θ (s ′ ) s |ψ θ (s)| 2 -⟨ Â⟩ s |ψ θ (s)| 2 (O k (s) + O * k (s)) s |ψ θ (s)| 2 , which implies ∂ θ k ⟨ Â⟩ = ⟨A loc O * k ⟩ -⟨A loc ⟩⟨O k ⟩ + c.c. (4.14)

III Monte-Carlo Markov Chains

To calculate these quantities in order to perform gradient descent, one needs to truncate the sum over s in Eq. ( 4.13) for it to become tractable. A way to do this is to obtain samples according to P , so as to have a correct estimation of the averages in the expression of the gradient. This will lead to averages given by

⟨ Â⟩ ≈ 1 M M i=1 A loc (s (i) ), (4.15)
where s (i) are samples drawn according to P . The statistical error associated with such an estimate scales as ε ∝ σ2 /M , and it is bounded as long as the variance σ 2 of A loc is finite 2 . For instance, when  is a k-local spin operator with bounded coefficients, its variance is strictly finite. 3 Therefore, the error in the estimate of expectation values decreases as ε ∼ 1/ √ M , which allows us to reach arbitrary accuracy in the estimation by increasing the number of samples M , given that lim M →∞ ε = 0. However, generating a set of samples according to the Born distribution, {s (i) } ∼ P (s), is in general a non-trivial computational task in the case where the variational ansatz, Ψ θ (s), is parameterized by an efficiently computable, yet arbitrary function. One of the most commonly adopted strategies to sample from P (s) is through Markov chain Monte Carlo [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] methods, including the Metropolis-Hastings method [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF], which generate a sequence of correctly distributed samples s (i) .

Metropolis-Hastings methods construct a markovian stochastic process which satisfies the detailed balance condition for the target probability distribution

P (s)T (s → s ′ ) = P (s ′ )T (s ′ → s) , (4.16)
where T (s (i) → s (i+1) ) is the probability that the state s (i) at step i transitions to the state s (i+1) at the following step. As the process is Markovian, the transition probability at every step depends exclusively on the current configuration. Assuming that the process is ergodic, the detailed balance condition ensures that, regardless of the initial configuration s (0) , the sequence eventually converges to the correct distribution P (s) in the long time limit.

One possible choice of the transition probability T is given by the Metropolis-Hastings algorithm [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. The main idea is to express T in terms of a local transition kernel T and an acceptance probability A such that

T (s → s ′ ) = T (s → s ′ )A(s → s ′ ) . (4.17)
This way, we split the global stochastic process into the product of two local subprocesses that we can compute efficiently. For instance, it is very easy to find a normalized local transition kernel that allows us to modify only a few degrees of freedom, like flipping a single spin in a given configuration. Conversely, it is hard to find a normalized global kernel that would act on all spins. The acceptance probability to go from a configuration s to s ′ through a local transition is defined as

A(s → s ′ ) = min 1, P (s ′ )T (s ′ → s) P (s)T (s → s ′ ) . (4.18)
Notice that the normalization of the Born probabilities cancel out, giving the tractable expression

P (s ′ ) P (s) = ⟨s ′ |ψ θ ⟩ ⟨s|ψ θ ⟩ 2 , ( 4.19)
which allows us to consider unnormalized variational ansätze. Therefore if the variational state is computationally tractable, the transition probability is also tractable. Choosing a valid transition rule T (s → s ′ ) is important, particularly in the case of systems whose states obey certain symmetries. For example, if the total magnetization along the direction of the computational basis is known, it may be fixed during the sampling procedure, and we can use a transition rule that does not project the Markov chain outside of a certain region in Hilbert space. For ground state search, a somewhat expensive yet effective choice for the transition kernel is to use the Hamiltonian itself:

T (s → s ′ ) = ⟨s| Ĥ|s ′ ⟩ (1 -δ s,s ′ ) s ′ ̸ =s ⟨s| Ĥ|s ′ ⟩ , (4.20)
which is known as the Hamiltonian transition rule [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF].

Finally, note that obtaining samples through Markov chains yields issues, as the samples are generally correlated. In addition, since the proposed moves are generally local, this can fail terribly for states that have multimodal distributions in Hilbert space. We will come back to these issues while discussing direct sampling in the next section.

IV Representing the wavefunction

Now that we have seen how to compute the quantities of interest using parametrized quantum states, let us dive into how to devise expressive variational states in practice. The main idea is that we need to represent high-dimensional functions with a parametrization that is flexible and general enough to describe physical systems, while involving only a polynomial amount of parameters.

Traditionally, researchers have relied on physically-inspired variational ansätze. An extension of the mean-field ansätze presented in the introduction is the Jastrow wave function [START_REF] Jastrow | Many-body problem with strong forces[END_REF][START_REF] Manousakis | The spin-½ Heisenberg antiferromagnet on a square lattice and its application to the cuprous oxides[END_REF], which stands out as one of the most successful and widely used ones for quantum chemistry. It is based on the assumption that two-body interactions are the most physically relevant, and it assigns a trainable potential to every interacting pair. Formally, the Jastrow wave function is given by:

ψ J θ (s) = e -1 2 i̸ =j θ ij σ i σ j , (4.21)
where the sum runs over all possible spin pairs, and θ ij are the parameters encoding pairwise spin correlations. Therefore, for a system of L spins, the resulting wave function has L 2 parameters. Moreover, in translationally invariant systems, the parameters θ ij can be made depend exclusively on the distance between i and j, resulting in a reduced number of parameters L.

To design more general variational ansätze, one can use neural networks. In this case, ψ θ (s) corresponds to the output of a neural network, that takes the configuration s as input in the form of a vector. In a feed-forward neural network of depth D, every layer l consists of a nonlinear activation function g (l) that acts, component-wise, on a vector resulting from applying the weight matrix W (l) to the output of the previous layer. This way, it is possible to write the variational state as the composition of operations g (l) ⊙W (l) , where ⊙ indicates point-wise operation, such that ψ θ (s) = g (D) ⊙ W (D) . . . g (2) ⊙ W (2) g (1) ⊙ W (1) s .

(4.22)

Hence, the output is a complex scalar representing the probability amplitude of configuration s.

From a mathematical perspective, these ansätze are of great interest given that neural networks are subject to universal representation theorems [START_REF] Kolmogorov | On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition[END_REF]. Consider a neural network without nonlinear activation functions. The function realized by such a network is a simple affine map that consists of multiplying the input by a weight matrix and adding to it an additional bias vector. This is evidently not enough to represent nonlinear functions. Thus, the addition of nonlinear activation functions is crucial for neural networks to be able to represent a larger class of functions. Kolmogorov and Arnold [START_REF] Kolmogorov | On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition[END_REF] have shown that any arbitrary continuous high-dimensional function can be expressed as a linear combination of the composition of a set of nonlinear functions

f (x) = 2m i=0 ζ i   m j=1 ς i,j (x j )   , ( 4.23)
where ζ i , ς i,j are nonlinear functions that act on the individual components of the input x ∈ R m . This means that we could represent any function f (x) with a polynomial number O(m 2 ) of one-dimensional nonlinear functions. This strongly resembles the structure of a neural network with two hidden layers. Note, however, that the nonlinear functions must be carefully chosen depending on the target function. In most machine learning architectures, the nonlinearities are typically fixed, i.e. ζ i = ς i,j ∀i, j. It turns out that fully-connected neural network composed of a single hidden layer and nonlinear activation functions are also universal function approximators. That is, given that the target function is reasonably well-behaved it can be approximated to any desired accuracy given that its hidden layer contains enough nodes [START_REF] Goodfellow | Deep Learning[END_REF][START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF][START_REF] Hornik | Approximation capabilities of multilayer feedforward networks[END_REF]. Note that this may still require a hidden layer that is exponentially large in the number of nodes. This raises the question what one can achieve with neural networks that have multiple hidden layers. The universal approximation theorem guarantees that there exists a neural network, i.e., choice of the architecture, which approximates the given target function arbitrarily well. However, it does not guarantee that we are able to find this choice. It turns out that, in practice, deep neural networks are capable of solving many problems with much less nodes, i.e., trainable parameters, compared to shallow neural networks. However, these results hold for arbitrary nonlinear functions, ζ i , ς i,j in 4.23, that must be appropriately found in order to represent the target function. In practice, neural networks use a fixed nonlinear activation, and we can only adjust the number of operations. In these cases, the number of neurons does not have a strict polynomial scaling and it can be, in the worst case, exponential in L [START_REF] Cybenko | Approximation by superpositions of a sigmoidal function[END_REF]. Nevertheless, the state-of-the-art results in computer vision and natural language processing [START_REF] Brown | Language Models are Few-Shot Learners[END_REF][START_REF] Caron | Emerging Properties in Self-Supervised Vision Transformers[END_REF][START_REF] Nichol | GLIDE: Towards Photorealistic Image Generation and Editing with Text-Guided Diffusion Models[END_REF] should be sufficient motivation to employ similar techniques to represent quantum states. In recent years, a large number of variational ansätze based on neural networks have also been developed for fermionic [START_REF] Ruggeri | Nonlinear network description for manybody quantum systems in continuous space[END_REF][START_REF] Hermann | Deep-neural-network solution of the electronic schrödinger equation[END_REF] and bosonic systems [START_REF] Saito | Solving the Bose-Hubbard Model with Machine Learning[END_REF][START_REF] Pescia | Neural-network quantum states for periodic systems in continuous space[END_REF].

IV.1 Restricted Boltzmann machines

NQS were first introduced using restricted Boltzmann machines (RBMs) [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF]. RBMs are shallow models featuring two fully-connected layers: a visible layer, consisting of L units, and a hidden layer, consisting of M units. A scheme of an RBM architecture is presented in 4.1. The wave function amplitudes of an RBM ansatz are given by:

Ψ θ (s) = h e b † v s+b † h h+h † W s . ( 4.24)
where s, h represent the visible and hidden units, respectively, and the parameters θ = {b v , b h , W } represent the visible and hidden biases and the weight matrix, respectively. In the NN picture, the RBM is a single-layer nonlinear feed-forward NN, with the visible units serving as inputs and the exponential serving as the activation function. By construction, RBMs are designed in such a way that computing the summation over hidden units, as in Eq. (4.24), can be done analytically. To see this, we can rewrite Eq. (4.24) in a tractable form considering binary hidden units h i ∈ {-1, 1}, leading to

ψ θ (s) = e b † v s M i=1 2 cosh (b h,i + W i• s) , ( 4.25)
where b h,i and W i• denote the i-th hidden bias and weight matrix row, respectively. To treat spin systems, the visible units will represent the N physical spins, thus the input of Pictorial representation of a RBM that represents the wave function of an Lspin system, with s = (σ 1 , σ 2 , . . . , σ L ) and h = (h 1 , h 2 , . . . , h d ) the hidden units. Adapted from [START_REF] Dawid | Modern applications of machine learning in quantum sciences[END_REF].

the RBM is simply the spin configuration s. In this way we obtain an analytical expression to evaluate the amplitude for a given spin configuration, and thus represent the full wave function with this ansatz. One can also interpret the hidden units as M hidden spins, and in this picture the RBM can be thought of as an interacting spin model with interaction strengths W ij . Moreover, we can treat an RBM as a model with an associated energy depending on its parameters, input, and hidden spin values. This is known as an energybased model, and explains why input biases are present in Eq. (4.24). In fact, the RBM is equivalent to a Hopfield network, a type of spin glass [START_REF] Barra | On the equivalence of Hopfield networks and Boltzmann machines[END_REF].

Being the first to be introduced in this context, most of the early works about NQS employ RBMs, but other architectures have been explored in more recent years, that we present in the next section. The capacity of RBMs and its relationship to quantum entanglement has been examined in various works [START_REF] Deng | Quantum entanglement in neural network states[END_REF][START_REF] Chen | Equivalence of restricted Boltzmann machines and tensor network states[END_REF]. An extension of this architecture, the deep RBM, has also been introduced to solve more complex problems [START_REF] Gao | Efficient representation of quantum many-body states with deep neural networks[END_REF], which consists in stacking more than two fully connected layers.

IV.2 Autoregressive neural networks

Autoregressive neural networks can also be used for constructing NQS, as introduced in [START_REF] Sharir | Deep autoregressive models for the efficient variational simulation of many-body quantum systems[END_REF] and applied to both quantum [START_REF] Luo | Gauge equivariant neural networks for quantum lattice gauge theories[END_REF] and classical problems [START_REF] Wu | Solving statistical mechanics using variational autoregressive networks[END_REF]. Their main advantage is that the wave function outputs of such models are normalized, allowing for direct sampling, which is easier to parallelize than Markov chain Monte Carlo and leads to independent samples. A pictorial representation of both the network and the sampling algorithm is presented in 4.2. Here, the many-body wave function is expressed in terms One samples consecutive spins using direct sampling on the conditional probabilities at each step. From [START_REF] Dawid | Modern applications of machine learning in quantum sciences[END_REF]. of a product of conditional complex amplitudes:

ψ θ (s) = L i=1 ϕ i (σ i |σ i-1 , . . . , σ 1 ) , (4.26)
which is subject to the normalization condition σ ϕ i (σ|σ i-1 , . . . , σ 1 ) 2 = 1. With this kind of architecture, we can compute expectation values by directly sampling state configurations instead of building a Markov Chain through the Metropolis-Hastings algorithm. This is done by iteratively sampling one spin after the other: we start sampling the first spin σ 1 from the reduced probability distribution |ϕ(σ 1 )| 2 . Then, we sample the second one σ 2 according to the conditional probability distribution |ϕ 2 (σ 2 |σ 1 )| 2 , then the next one |ϕ(σ 3 |σ 2 , σ 1 )| 2 , and so on until σ N . This sampling procedure is embarrassingly parallel. 4This sampling procedure yields independent identically distributed samples. Conversely, Markov chain Monte Carlo methods may suffer from highly correlated consecutive samples, which is problematic for complex probability distributions, e.g. that are far from Gaussian. Markov chain Monte Carlo methods such as the Metropolis-Hastings algorithm, generally rely on performing modifications to the spin configurations to sample subsequent states. Therefore, this process could yield highly correlated consecutive samples that may have a negative impact on the results. In order to compute expectation values, we need to estimate the autocorrelation time to draw uncorrelated samples from the resulting chain. Moreover, when approaching a phase transition, such methods suffer from critical slowing down, making the sampling of uncorrelated configurations unfeasible in many situations.

Consider a quantum state that spans several separated regions in the Hilbert space, where the probability is concentrated. In this case, Markov chains generally remain stuck in one of the regions, given that it must take several penalizing steps to travel from one to another, resulting into a highly inaccurate sampling. In contrast, the direct sampling procedure can seamlessly draw spin configurations belonging to all the regions according to the probability distribution, yielding much better samples.

Let us finish by mentioning that a plethora of different architectures have been implemented as NQS in recent years, such as convolutional neural network wave functions [START_REF] Schmitt | Quantum many-body dynamics in two dimensions with artificial neural networks[END_REF], and, more recently, group convolutional neural networks [START_REF] Roth | Group Convolutional Neural Networks Improve Quantum State Accuracy[END_REF], which are especially convenient to implement certain symmetries directly built-in to the wave function.

IV.3 The gated recurrent unit ansatz

An ansatz that we will employ in section VI is the gated recurrent unit (GRU) ansatz, a type of recurrent neural network (RNN). This ansatz is autoregressive, and its conditional amplitudes are obtained through the following transformation, which depends on s i , the local spin, and h i-1 , the hidden unit coming out of the previous recurrent cell:

ϕ i (s i , h i-1 ) := GRU(s i , h i-1 ) = A(s i ) exp{iφ(s i )}, (4.27)
A(s i ) = ς(U A h i + b A ), (4.28)
φ(s i ) = ςς(U φ h i + b φ ), (4.29)
where ς and ςς denote a softmax and softsign transformation respectively, and where we have defined the conditional modulus A(s i ) and phase φ(s i ) that depend on matrices U A and U φ respectively, which each contain variational parameters. The hidden vector h i is given by

h i = (1 -z j ) ⊙ h i-1 + z j ⊙ hj , ( 4.30)
with the latent hidden vector hj given by

hj = tanh W[r j ⊙ h j-1 , s j ] + b (4.31)
and the update gate z j and reset gate r j :

z j = sig (W z [h j-1 ; s j-1 ] -b z ) , ( 4.32
)

r j = sig (W r [h j-1 ; s j-1 ] -b r ) , ( 4.33)
where sig denotes the sigmoid funtion and where we have defined the matrices W , W r , W z and bias vectors b, b r , b z which are all variational parameters. [h j-1 ; s j-1 ] denotes a concatenation of vectors h j-1 and s j-1 , the latter corresponding to the one-hot encoding of the local spin configuration s j-1 . The important point here is that the total number of variational parameters scales quadratically with the dimension of the hidden vector d h , since the W matrices each contain d h × (d h + 2) variational parameters for a one-hot encoding of the local spins s i . One can see from this transformation that the GRU ansatz naturally respects the autoregressive property. Panel (a) is for real-valued wave functions, which can be relevant for a certain class of problems, and panel (b) is for complex-valued wave functions. In both schemes, a local spin configuration s i and a hidden vector h i are fed into an RNN cell, which performs a nonlinear transformation. Then an activation function (ς, for softmax and/or ςς, for softsign) is applied to obtain the final probability and/or phase corresponding to the configuration. At the end, the probabilities (and phases) are combined to obtain the final wave function amplitudes ψ(s).

s 2 s 1 s 0 s 2 s 1 s 0 RNN P 1 h 1 RNN P 2 h 2 h 0 h 1 RNN P 3 h 3 ψ ( s ) h 2 (a) RNN P 1 ς h 1 RNN h 0 h 1 ψ ( s ) (b) ςς Φ 1 P 2 h 2 Φ 2 RNN P 3 h 3 Φ 3 h 2 ς ςς ς ςς ς ς ς

V Capacity and entanglement

A natural question arises regarding the expressive capacity of NQS and how different model architectures compare to each other in this respect. In particular, understanding how they compare to tensor network methods is crucial, as it would be natural that there is some limit, perhaps stemming from entanglement structure as for tensor networks. However, as is common in machine learning, this question is tricky to answer. A point we do not cover here but that is also important is that trainability is also a common issue for machine learning methods: an ansatz may be highly expressive but not so trainable, hence in practice it will not be useful.

Tensor networks have been a recurrent tool to study many-body systems. They are now well established and characterized, and their limits are well understood, hence they constitute a conceptual framework to study the complexity of quantum states. For this reason, there has been a significant community effort to study the relationship between tensor network states and NQS [START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF][START_REF] Chen | Equivalence of restricted Boltzmann machines and tensor network states[END_REF][START_REF] Glasser | Neural-network quantum states, string-bond states, and chiral topological states[END_REF] , which provides insight about the expressive capacity of NQS [START_REF] Levine | Quantum entanglement in deep learning architectures[END_REF]. Following the first introduction of NQS implementing RBM [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF], early works focused on finding direct relationships between various kinds of RBM-based states and tensor network states [START_REF] Chen | Equivalence of restricted Boltzmann machines and tensor network states[END_REF][START_REF] Glasser | Neural-network quantum states, string-bond states, and chiral topological states[END_REF]. Recently, it has been proven that NQS can efficiently approximate, in logarithmic space, all efficiently contractible tensor network states with arbitrary precision [START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF]. Therefore, for every tensor network state there exists an equivalent NQS of polynomial size. Conversely, there are quantum states that can be efficiently described by NQS, whose representation in terms of tensor networks requires an exponential amount of parameters. Hence, tensor network states are a subset [START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF] by mapping tensor networks to NQS. PEPS* refers to a sub-class of projected entangled pair states, a generalization of MPS for two-dimensional systems.

Adapted from [START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF].

of NQS [START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF], as depicted in 4.4.

More precisely, one study the entanglement scaling captured by the different ansätze. In a generic quantum many-body system with density matrix ρ, the entanglement entropy is defined as

S(ρ) = -Tr [ρ log 2 ρ] , (4.34)
which is zero for any pure state. Let us consider a partition of the system in two subsets: I and its complementary O, as well as the reduced density matrix ρ I = Tr O [ρ]. In general, ρ I represents a mixed state, which can have nonzero von Neumann entanglement entropy.

For a generic quantum state, the entanglement entropy of ρI grows with the volume of the cut. Thus, it corresponds to a volume-law scaling. NQS can efficiently capture such scaling with architectures ranging from very basic shallow ones, such as RBM [START_REF] Deng | Quantum entanglement in neural network states[END_REF], to more modern and deeper approaches, such as convolutional neural networks or RNNs [START_REF] Levine | Quantum entanglement in deep learning architectures[END_REF].

We can understand most differences between the ansätze at an intuitive level by, simply, looking at how they are built. In 4.5, we provide a pictorial representation of the different connections that some ansätze can draw in a bi-dimensional system. Clearly, the MPS ansatz, depicted in 4.5(b), is the most locally restricted one, as it can only account for nearest neighbour connections in a snake-like pattern. This effectively limits [START_REF] Hibat-Allah | Variational neural annealing[END_REF].

(b) (a) (d) (c)
the entanglement that MPS can capture. The RNN ansatz, illustrated in 4.5(d), while it is limited to parse the state in the same pattern as the MPS, it has the freedom to account for additional information, allowing it to capture richer correlations.

In contrast, other ansätze such as the Jastrow or RBM wave functions, respectively illustrated in 4.5(a) and (c), can draw connections between arbitrary sites. The Jastrow ansatz can account for all possible pairs in the system, regardless of the distance. Then, the RBM anstaz is a generalization of the Jastrow by means of an auxiliary hidden layer of variable size. Through the hidden neurons, the ansatz is no longer limited to pairs, and it can actually consider up to all-to-all connections. This non-local character allows them to capture volume-law entanglement.

VI Dynamics of closed systems with NQS

We now turn to the application of NQS to dynamics of nonequilibrium quantum systems. Impressive results for the transverse-field Ising model on large lattices have been obtained using convolutional neural networks [START_REF] Schmitt | Quantum many-body dynamics in two dimensions with artificial neural networks[END_REF]. It is therefore of particular interest to continue investigating such problems that could yield important insights into nonequilibrium phenomena. In particular, studying the quantum Kibble-Zurek mechanism in higher dimensions has remained relatively untouched [START_REF] Schmitt | Quantum phase transition dynamics in the two-dimensional transverse-field Ising model[END_REF], and NQS is deemed to be employed for the dynamics of time-dependent finite-size systems. While most previous works focus on the restricted Boltzmann machine (RBM) ansatz, more recent works on ground state search have employed networks that are closer to the state-of-the-art in machine learning such as autoregressive convolutional models [START_REF] Sharir | Deep autoregressive models for the efficient variational simulation of many-body quantum systems[END_REF] or recurrent neural networks [START_REF] Hibat-Allah | Recurrent neural network wave functions[END_REF], for which the accuracy of the variational ground state energy was significantly improved. These networks have a so-called autoregressive structure, which means that one can perform direct sampling of uncorrelated configurations for arbitrarily large system sizes, thereby reducing the number of required samples. It is therefore of crucial importance to investigate the use of similar networks for quantum dynamics.

However, the widespread application of more complex and autoregressive neural networks for quantum dynamics has been held back by issues arising with time-dependent variational Monte-Carlo (t-VMC) [START_REF] Becca | Quantum Monte Carlo approaches for correlated systems[END_REF]. This procedure, used to propagate an NQS in time, involves the inversion of a stochastically constructed singular matrix, which makes it particularly prone to noise. While some regularization techniques have helped improving the accuracy of the method [START_REF] Schmitt | Quantum many-body dynamics in two dimensions with artificial neural networks[END_REF][START_REF] Hofmann | Role of stochastic noise and generalization error in the time propagation of neural-network quantum states[END_REF], accessing all regimes at long times via t-VMC remains a challenge [START_REF] Dawid | Modern applications of machine learning in quantum sciences[END_REF][START_REF] Hofmann | Role of stochastic noise and generalization error in the time propagation of neural-network quantum states[END_REF][START_REF] Park | Geometry of learning neural quantum states[END_REF][START_REF] Vicentini | NetKet 3: Machine Learning Toolbox for Many-Body Quantum Systems[END_REF].

In this section, we show that the stability of t-VMC strongly depends on the chosen ansatz, and that, in particular, it fails when applied to those based upon recurrent neural networks (RNN). To circumvent this issue, we propose an alternative scheme to numerically solve the dynamics of quantum systems. The scheme consists in casting an arbitrary Runge-Kutta integration scheme of any order into minimizing a variational distance at each time step, while only involving a polynomial overhead. This enables one to employ recurrent neural networks for quantum dynamics, which leads to a drastic reduction the the number of required samples. Our scheme can be implemented for any order of a chosen integration method, for a polynomial overhead in memory. We then apply our scheme to both time-dependent and sudden quenches, enabling one to recover Kibble-Zurek scaling laws for large system sizes and high precision on the dynamics during the full quench.

VI.1 Issues with t-VMC

In order to propagate a NQS in time according to some Hamiltonian Ĥ, one should minimize the following variational distance:

D( θ) = dist |ψ θ⟩ , exp -i Ĥτ |ψ θ ⟩ (4.35)
over the variational parameters θ after each time step τ . By expanding |ψ θ⟩ as a Taylor expansion in the variational parameters to obtain the following equation to be solved at with open boundary conditions. The system is subject to the TFI Hamiltonian, and configurations are fed into GRU in an autoregressive fashion, each cell being fed with a hidden vector h i,j and a spin σ i . After the GRU cell, a softmax (ς) and a softsign (ςς) transformation are applied to obtain the amplitudes and the phases of the conditional amplitudes ϕ i (σ i |σ i-1 , .., σ 1 ). (b) Spectra of the quantum geometric tensors S for different ansätze optimized to the ground state of the TFI model, for various values of g/J. The left panel is for a complex RBM ansatz with α = 2, the central panel for a two three layer convolutional autoregressive network with complex parameters (ARNN), and the right panel for a GRU ansatz with real parameters with d h = 10, for which is plotted the spectrum of the imaginary part of the S matrix, as this is the quantity to invert for ansätze with real parameters. each time step, so as to obtain θk := ∂ t θ k [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF]:

S k,k ′ θk ′ = -iF k , ( 4.36)
with

F k = ⟨O * k E loc ⟩ -⟨O * k ⟩⟨E loc ⟩, and S k,k ′ = ⟨O * k O k ′ ⟩ -⟨O * k ⟩⟨O k ′ ⟩
when keeping only second-order terms. The O k terms are the log-derivatives of the NQS ansatz, as given by

O k (σ) = ∂ θ k ln ψ θ (σ). (4.37)
With most machine learning algorithms, these are obtained through automatic differentiation. The local energy E loc is local estimator which evaluates as

E loc (σ) = σ ′ ψ θ (σ ′ ) ψ θ (σ) ⟨σ| Ĥ|σ ′ ⟩ . (4.38)
The parameters θ are then updated at each time step according to their derivatives θk using a numerical solver. This procedure is known as time-dependent variational Monte Carlo (t-VMC), and is related to natural gradient methods in machine learning [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF].

Several issues can arise: (i) The S matrix to be inverted is in general singular, which makes the process extremely sensitive to stochastic variations coming from sampling. Hence, regularization must be used to obtain a non-diverging parameters derivative θk [START_REF] Hofmann | Role of stochastic noise and generalization error in the time propagation of neural-network quantum states[END_REF]. This increases the stability of the method, although impacting its accuracy, making it challenging to obtain accurate long-time dynamics. (ii) State-of-the-art regularization techniques involve the inversion of the S matrix via its singular-value decomposition, which makes the complexity of the method of order O(N 3 par ), with N par the number of parameters for the NQS. 5 This is one of the reasons why natural gradient methods is rarely used for other machine-learning applications, which can involve models with up to billions of parameters. (iii) For NQS with real parameters, one may split the real and imaginary parts of Eq. (4.36) and solve either of the two resulting equations [START_REF] Yuan | Theory of variational quantum simulation[END_REF]:

Re(S kk ′ ) θk ′ = 0 Im(S kk ′ ) θk ′ = -F k .
Solving the first forces one to impose a condition on θk to obtain a nontrivial solution ( θ ̸ = 0), and solving the second one is difficult in general since the diagonal elements of Im(S) are zero (S is Hermitian) and remaining off-diagonal elements are close to 0 for various ansätze with real parameters we have tested. When the minimum and maximum eigenvalues of a matrix are highly separated, λ min /λ max ∼ ϵ, ϵ → 0, this is known as poor conditioning [START_REF] Goodfellow | Deep Learning[END_REF]. It makes regularization schemes even harder to implement and yields large parameter time derivatives, thereby requiring impractically small time steps for t-VMC to work.

To further investigate such issues, let us consider the transverse-field Ising (TFI) Hamiltonian:

ĤTFI = -J ⟨m,n⟩ σz m σz n + g m σx m , (4.39)
with J the nearest-neighbor coupling strength, and g the transverse field strength. As mentioned in section II.1, at zero temperature this model exhibits a quantum phase transition for g c = J in one dimension and for g c ≈ 3.044J in two dimensions [START_REF] Blöte | Cluster Monte-Carlo simulation of the transverse ising model[END_REF]. For J > 0, the transition separates a ferromagnetic6 phase from a paramagnetic phase, where the spins tend to align with the transverse-field. In the former, when g ≪ J, the ground state is degenerate and spanned by |↑, ↑, . . . , ↑⟩ and |↓, ↓, . . . ↓⟩, while in the latter, when g ≫ J the ground state is |→, →, . . . , →⟩, with |→⟩ the eigenstate of σx . This model serves as a convenient benchmark for the dynamics of NQS [START_REF] Carleo | Solving the quantum many-body problem with artificial neural networks[END_REF][START_REF] Schmitt | Quantum many-body dynamics in two dimensions with artificial neural networks[END_REF][START_REF] Gutiérrez | Real time evolution with neural-network quantum states[END_REF]. In Fig. 4.6(b), the spectra of S, computed at the ground state of the TFI Hamiltonian is shown for three different ansätze: a restricted Boltzmann machine (RBM), an autoregressive convolutional network with complex parameters (ARNN), and the GRU ansatz presented in the previous section, for different values of g/g c (we have fixed J = 1) for a two-dimensional array of N = 16 spins. Two key observations can be made from this figure: first, all spectra become more singular as g/J decreases, something that makes t-VMC less applicable to ferromagnetic-like states, whether it is for ground state search or for dynamics, which involves the full spectrum of the Hamiltonian. One can also see that the RBM spectrum is the less singular, and presents a shell-like structure, making it simpler to regularize. The ARNN ansatz is similar to this, and presents many eigenvalues of order O(1). In stark contrast, the GRU spectrum indicates that the variational space that is extremely flat for all values of g/J, as very few eigenvalues are of a high enough value, and the magnitude of the eigenvalues decrease extremely rapidly. This inherent poor conditioning implies that if such states are the initial states of a dynamical evolution, the inversion of S will yield a parameter update θ whose norm is large with respect to ∥θ∥, forcing one to dramatically reduce the time step which makes calculations impractical. This last spectrum is tricky to regularize, as there is no clear separation of the eigenvalues.

VI.2 Variational Runge-Kutta algorithms

As identified in the previous section, t-VMC suffers from a number of limitations that make it unsuitable for use with certain ansätze, in particular for the GRU ansatz. We propose to circumvent these issues by resorting back to the original Dirac-Frenkel variational principle (Eq. (4.35)) and solving an optimization problem at each time step. The general form of this problem is the minimization of the distance

D( θ) = dist |ψ θ⟩ , T |ψ θ ⟩ , ( 4.40)
where θ denotes the set of variational parameters to optimize, θ those at the previous time step t, and T a propagator evolving the state of the system between times t and t + τ under the action of the Hamiltonian of interest. While expanding the propagator to first order in τ yields a valid first-order update of the variational state, we instead propose to build a variational principle upon a s-order Runge-Kutta approximant, as generated by a propagator Ts such that:

|ψ(t + τ )⟩ = Ts (t) |ψ(t)⟩ + O(τ s+1 ). (4.41)
General explicit expressions for this operator are derived in Appendix A for arbitrary diagonally implicit Runge-Kutta schemes. In particular, for the standard (second-order) Heun method used throughout the paper, one has:

THeun (t) = 1 -iτ Ĥ(t) - τ 2 2 Ĥ(t + τ ) Ĥ(t). (4.42)
The variational method induced by this scheme involves an integration error per time step of third order in τ . In practice, rather than the Fubini-Study distance, we use the following numerically well-behaved metric based upon the quantum fidelity: At times ± t, the relaxation time of the system t R crosses the time scale of the quench, which defines the limit between the adiabatic region, where the system essentially stays in the ground state of the system at that time, and the impulse regime, where nonequilibrium dynamics take place. (b) Scaling of the injected energy Ĥ(τ q ) -E 0 as a function of the quench time τ q for different system sizes N = L × L. One can see that for fast quenches, the NQS dynamics predict very well the injected energy. The data for system sizes with more than L = 5 cannot be simulated via exact diagionalization, hence they are compared to the τ 0.72 q scaling obtained in ref. [START_REF] Schmitt | Quantum phase transition dynamics in the two-dimensional transverse-field Ising model[END_REF]. (c) Average energy for a quench of total time Jτ q = 0.31, for a 6 × 6 system and a GRU ansatz with where Tloc is an operator acting as a local estimator and whose non-zero entries are given by

dist |ψ⟩ , |ϕ⟩ = 1 - |⟨ψ|ϕ⟩| 2 ⟨ψ|ψ⟩ ⟨ϕ|ϕ⟩ , ( 4
T loc (σ) := ⟨σ| Tloc |σ⟩ = σ ′ ψ θ (σ ′ ) ψ θ(σ) ⟨σ| Ts |σ ′ ⟩ . (4.47)
It thus appears from the above that any s-order Runge-Kutta update can be cast into a variational problem consisting in maximizing the squared expectation value of an observable. Furthermore, this expectation value can be efficiently sampled at every step of the optimization process as the average

⟨ψ θ| Tloc |ψ θ⟩ = E σ∼|ψ θ | 2 [T loc (σ)]. (4.48)
The error of an integration using the Heun propagator is of the order of τ 3 , which will be the lower bound for the error at each time step for the variational method we present. Note that this involves summing over the connected elements of T as can be seen in Eq. (4.47). This will in general become more expensive as the order s increases, as powers of the Hamiltonian will be contained in T . The number of connected elements will depend Numerical results obtained for a linear quench for various values of the total quench time τ q . On the upper panels, the average correlation function is shown as a function of time for both the exact (dashed line) and the GRU simulations (solid lines). The error as a function of time ϵ(t) is shown in the lower panels. For panels (b) and (c), the time step was chosen to be linear and non-constant, as the characteristic time scale changes as a function of time. Parameters for each quench are are given in the appendix. on the sparsity of the Hamiltonian; for the transverse-field Ising model for instance, the number of connected elements in Ĥs scales as N s (which means N 2 for a second-order explicit integration scheme).

VI.3 Application: critical quench dynamics The Kibble-Zurek mechanism

The Kibble-Zurek mechanism [START_REF] Kibble | Topology of cosmic domains and strings[END_REF][START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF] predicts the formation of topological defects in quenched systems undergoing a second-order phase transition as the parameters are linearly swept across a critical point. This universal behavior stems from the fact that there exists a time, denoted t, at which the characteristic evolution time of the system (the relaxation time t R , related to the inverse of the energy gap) becomes larger than the characteristic quench time τ q , as schematically shown in Fig. 4.7(a). Before t, the dynamics remains quasi-adiabatic, as the Hamiltonian parameters are tuned slowly with respect to the time scale of the system. After this time, non-equilibrium dynamics takes place as critical slowing down sets in and the system dynamics becomes slower than the Hamiltonian parameters sweep. Ignoring specific details of this dynamics, one can derive a scaling law for the density of created defects at the end of the linear quench [START_REF] Zurek | Cosmological experiments in superfluid helium?[END_REF], which reads ⟨n d (τ q )⟩ ∼ τ -dν/(zν+1) q (4.49) with d the dimensionality of the system and z, ν universal critical exponents. Note that nd will depend on the geometry of the system, and is not always easily accessible experimentally. For this reason some works consider other observables, such as the number of domain walls [START_REF] Puebla | Quantum kibble-zurek physics in longrange transverse-field ising models[END_REF], or the injected energy density [START_REF] Schmitt | Quantum phase transition dynamics in the two-dimensional transverse-field Ising model[END_REF], defined as

E = 1 L 2 ⟨ Ĥ(t)⟩ -E 0 (t) . (4.50)
This quantity is a witness of the injected defects, which are excitations of the system with respect to its ground state at time t, hence E must scale as the density of defects. This quantity is convenient as it does not depend on the geometry of the physical system under consideration. In Ref. [START_REF] Schmitt | Quantum phase transition dynamics in the two-dimensional transverse-field Ising model[END_REF], authors extract a scaling law in τ 0.72 q for the 2D transversefield Ising model based upon various numerical methods that they used to estimate the gap of the system at different times to extract a correlation length. We consider a similar protocol in this paper, in which both J and g are varied in time. The Hamiltonian reads

ĤTFI = -J(t) ⟨m,n⟩ σz m σz n + g(t) m σx m , ( 4.51)
with J(t) = J(1 -t/τ q ), g(t) = g 0 (1 + t/τ q ), with g 0 = g c and J = 1, which will be used as a reference unit in the rest of the paper. Note that this type of quench involves all regimes (0 ≤ g/J < +∞), hence stochastic reconfiguration does not enable one to access the full dynamics [START_REF] Schmitt | Quantum phase transition dynamics in the two-dimensional transverse-field Ising model[END_REF]. Our approach does not suffer from such issues, which is one of its major strengths. In Fig. 4.7(b), we show the injected energy (not rescaled for clarity) ⟨ Ĥ(t)⟩ -E 0 (t) for various system sizes, for different values of the quench time. One can Curves collapse, which is what is predicted by a direct sampling scheme. This enables us to estimate the required number of samples to obtain a given error. Note that as the dynamics progress, the variance of the residual error rises, indicating a rising difficulty to find the best minima in time.

see that the exact results are closely matched for fast quenches, and precision drops for slower quenches (Jτ q > 1). This is related to the fact that dynamics become almost fully adiabatic, hence the NQS dynamics is not precise enough to capture the small amount of injected energy. The presented method most accurately reproduces the dynamics of the system in the non-adiabatic region, where the dynamics presents universal features; this is the regime of relevance when studying the Kibble-Zurek mechanism. For lattices above 5 × 5, one cannot compare our results with exact results, but the scaling predicted in [START_REF] Schmitt | Quantum phase transition dynamics in the two-dimensional transverse-field Ising model[END_REF] is recovered, indicating that the results are reliable. We also show correlation functions for various quench times, as well as the residual error corresponding to the full dynamics, as given by

ϵ(t) = dist |ψ θ⟩ , T |ψ θ ⟩ (4.52) = 1 -|⟨ψ θ| T |ψ θ ⟩ | 2 . (4.53)

VI.4 Sudden quenches

One can also investigate the nonequilibrium dynamics of many-body systems without relying on regularization hyperparameters and instabilities stemming from stochastic reconfiguration. As a benchmark of the presented method, we consider the two-dimensional time-independent transverse-field Ising model, defined in Eq. (4.39). We prepare the system in the ground state of the TFI Hamiltonian for g ≫ J, which is |ψ 0 ⟩ = |→, →, . . . , →⟩, and quench the magnetic field to the values g/g c = 2, 1, and 1/10; results are displayed in panels (a), (b), and (c) of Fig. 4.9 respectively. By doing this, one spontaneously creates excitations of all the eigenstates of the Hamiltonian corresponding to the final value of the parameters. This is of high interest experimentally, as it can be used to probe the properties of non-integrable systems close to the critical point. It is expected that a quench to the critical point is the most difficult to simulate, since it involves states that are correlated at all scales in the thermodynamical limit. In contrast to other approaches, the use of direct sampling here ensures that the dynamics will be accurate for a lower number of samples. Here, we have imposed the reflection and Z 2 symmetries for panels (a) and (b), and no symmetry for panel (c), where, surprisingly, imposing symmetries was found to degrade the accuracy.

VI.5 Error analysis

The residual error ϵ depends on the chosen time step δt. In Fig. 4.10, this is shown as a function of the chosen time step δt for different integration schemes, namely Euler, (implicit) midpoint, and Heun methods. The propagators T for the Euler and the implicit midpoint method are respectively given by TEul = 1 -iτ Ĥ (4.54)

Tmid = 1 -iτ Ĥ + τ 2 4 Ĥ2 , ( 4.55)
and the Heun propagator was given by Eq. (4.42). Euler is a first-order method, while the midpoint and Heun methods are both second-order methods. The midpoint method conserves a symplectic symmetry in the exact case, which for the Schrödinger equation corresponds to the energy. Note, however, that propagating an NQS in time, this is no longer guaranteed, as the time evolution is approximated stochastically. Using Tmid as a propagator is equivalent to minimizing the distance proposed in Ref. [START_REF] Gutiérrez | Real time evolution with neural-network quantum states[END_REF]. Interestingly though, when using the midpoint method, the norm of the denominator of Eq. (4.44) contains terms of order 2, whereas they vanish for the Heun method. Indeed, the implicit midpoint method yields an update of the form

|ψ(t + δt)⟩ = Tmid (t) |ψ(t)⟩ + O(δt 3 ), (4.56)
and as such may in principle be regarded as being of second order, provided one is able to approximate |ψ(t + δt)⟩ with an NQS up to an error of order O(δt 3 ). However, this condition cannot be satisfied with a normalized ansatz such as those used in this work. Indeed, the norm of the ideal updated state to be matched is given by:

⟨ψ(t)| T † mid Tmid |ψ(t)⟩ = ⟨ψ(t)| 1 + 3δt 2 2 Ĥ2 |ψ(t)⟩ + O(δt 3 ) = 1 + O(δt 2 ) . ( 4

.57)

This must depart from 1 to second order in δt, implying that the implicit midpoint method effectively yields a first-order update when using ansätze normalized by design.

In contrast to the midpoint method, the Heun method does not suffer from this issue, indeed:

T † Heun THeun = 1 + iδt Ĥ - δt 2 2 Ĥ2 1 -iδt Ĥ - δt 2 2 Ĥ2 = 1 + O(δt 4 ) . ( 4

.58)

This property allowed us to simplify the expression of the distance in Eq. (4.44), considerably reducing the complexity of evaluating the gradients of the loss function.

These considerations generalize to the time-dependent scenario, where for the Heun method we have

T † Heun (t) THeun (t) = 1 + iδt Ĥ(t) - δt 2 2 Ĥ(t) Ĥ(t + δt) 1 -iδt Ĥ(t) - δt 2 2 Ĥ(t + δt) Ĥ(t) = 1 + δt 2 Ĥ(t) 2 - δt 2 2 { Ĥ(t), Ĥ(t + δt)} + O(δt 3 ). ( 4.59)
Provided one can expand Ĥ(t + δt) = Ĥ(t) + δt∂ t Ĥ(t) + O(δt 2 ), as is the case for any analytic quench, the second-order terms in Eq. (4.59) cancel out exactly, implying that the Heun method remains of second order when using normalized ansätze. Therefore, in this context, this method will effectively be of first-order. One can clearly observe in Fig. 4.10(a) a δt 4 scaling for the Heun method as expected from the form of the distance we minimize, and a δt 2 error for the two other methods. Hence the midpoint method carries the disadvantage of having an error comparable to an Euler integration while having the computational complexity of the Heun method. One can also see that for the Heun method there is a cutoff around gδt = 0.01, for which the error cannot be reduced by decreasing the time step. This is due to the fact that the method displays a systematic error coming from the optimization for a given set of hyperparameters (number of samples, learning rate, number of training step). Nevertheless, the error can always be decreased by increasing the number of samples. Note that in the constant error regime, reducing the time step will be penalizing as the total number of time steps at fixed total simulation time. Hence there is an optimal time step to be chosen, which we find to be gδt ≈ 0.01 for various hyperparameter choices. In Fig. 4.10(b), the standard deviation of the error is shown for each time step of a sudden quench from a paramagnetic phase to g/g c = 2. One can distinguish two regimes, separated by gt ≈ 1, where before this time the standard deviation of the residual error grows, and then saturates after this time. This can be attributed to a harder optimization error after a given time, where the minimum of the distance is found but with a greater variance. These results are reminiscent of those presented in a recent work where authors demonstrate that quantum states become harder to fit with an NQS at lated times after a sudden quench, which is not always related to the growth of entanglement [START_REF] Lin | Scaling of Neural-Network Quantum States for Time Evolution[END_REF]. Standard deviations are shown for various values of N s , the number of samples considered for each optimisation. One can see that σ ϵ (t) decreases as √ N s , as predicted by a direct sampling scheme.

VII Dynamics of open quantum systems with NQS VII.1 Representing density matrices

Extensions of neural network techniques to the time evolution of open quantum systems have been proposed. A first approach relies on purification [START_REF] Torlai | Neuralnetwork quantum state tomography[END_REF], which obliges one to use Restricted Boltzmann Machine (RBM) as the ansatz [START_REF] Vicentini | Variational neural-network ansatz for steady states in open quantum systems[END_REF][START_REF] Hartmann | Neural-network approach to dissipative quantum many-body dynamics[END_REF][START_REF] Yoshioka | Constructing neural stationary states for open quantum many-body systems[END_REF][START_REF] Nagy | Variational quantum Monte Carlo method with a neuralnetwork ansatz for open quantum systems[END_REF] for the density matrix, a requirement that bounds the representative power of such methods. More recently, new approaches based on autoregressive networks were proposed for the time evolution of open quantum systems. These are based on the Positive-Operator-Value-Measurement (POVM) formalism, in which the density matrix is represented as an average over outcomes of measurements [START_REF] Luo | Autoregressive neural network for simulating open quantum systems via a probabilistic formulation[END_REF][START_REF] Reh | Time-dependent variational principle for open quantum systems with artificial neural networks[END_REF]. By recasting the Lindblad master equation into a differential equation on a probability distribution, a time-dependent variational principle can be obtained. An issue with this approach is that the space of density matrices (positive semi-definite) is more restricted than that of arbitrary probability distributions, therefore the time evolving map may project distributions out of the space of density matrices, for example by violating positivity. Finally, a more recent approach consists in building an autoregressive and positive density matrix that relies on a construction based on Gram matrices [START_REF] Vicentini | Positive-definite parametrization of mixed quantum states with deep neural networks[END_REF]. We also present the extension of the dynamical-corner space method, which consists in representing each state of the corner basis by an NQS. This restricts the method to lowentropy quantum systems, therefore will not be practical to study, for instance, steadystates of critical systems (where the entropy can be high), but ensures positivity and avoids enforcing a certain structure for the ansatz. In principle, one could even combine different neural network architectures with such an approach, but this is generally computationally impractical.

The RBM approach

One wants to obtain a generative model for the coefficients of the density matrix ρ θ (s, s ′ ), defined as:

ρθ = s,s ′ ρ θ (s, s ′ ) |s⟩⟨s ′ | . ( 4.60)
The first approach consists in using an RBM to represent a pure state |ψ θ ⟩ that belongs to an enlarged Hilbert space H = H s ⊗ H a , with H s describing the physical spins and H a describing the ancilla spins. Provided that there are at least as many ancilla spins as physical spins, one can always write:

ρ θ (s, s ′ ) = a ψ θ (s, a)ψ * θ (s ′ , a). (4.61)
However, in the general case, the sum over a is intractable, as there are 2 L elements to sum over. The key point is that since the RBM ansatz has an analytical form, this sum can be performed analytically and a closed expression for ρ θ (s, s ′ ) can be obtained [START_REF] Vicentini | Variational neural-network ansatz for steady states in open quantum systems[END_REF]. This is a calculation much like that which led to Eq. (4.25), where a trace over hidden spins was performed analytically. This leads to a positive and Hermitian ansatz for the density matrix, with derivatives of the density matrix coefficients also having an analytical form. However, as we have seen, RBMs are restricted in terms of capacity, hence other approaches were quickly proposed to obtain more freedom on the choice of the neural network.

The Gram-Hadamard matrix approach

The Gram-Hadamard approach proposes to re-express the coefficients of a density matrix by using a Gram matrix form and useful properties of the Hadamard product [START_REF] Vicentini | Positive-definite parametrization of mixed quantum states with deep neural networks[END_REF]. The where the element-wise product is taken in the same basis for each Gram matrix. The elements of the density matrix are therefore given by ρθ (s, s

′ ) = K h=1 R a=1 ψ (a,h) θ (s)ψ * (a,h) θ (s ′ ). (4.65)
where ψ (a,h) θ (s) are complex amplitudes. Note that this ansatz closely resembles the purification of Eq. 4.61. A diagram of the ansatz is shown in Fig. 4.11. The power of such a construction is that by restricting the rank R of the Gram matrices to a small number, the rank of the final density matrix is obtained as R K . This form is reminiscent of the locallypurified density operator presented in section III.6 of Chapter 1, where a partial trace is taken on finite-rank states, therefore is tractable. However, there is no local purification here, as ψ (h) are global states. A natural question that arises is whether this construction could also be used with tensor networks architectures as an alternative to LPDOs, hence creating an alternative positive tensor network ansatz. Note that Gram-Hadamard construction may constrain the entanglement structure of ρθ , and further studies are needed to verify this, in which case the full power of machine learning methods could be deployed to study open quantum systems. Finally, in the same work [START_REF] Vicentini | Positive-definite parametrization of mixed quantum states with deep neural networks[END_REF], a similar construction has been proposed that explictily enforces the autoregressive property, to enable for direct sampling and normalized wavefunctions. However, again this construction may constrain the entanglement structure of the ansatz.

The POVM approach

Another approach that does not enforce positivity but allows absolute freedom on the choice of the ansatz will be referred to as the POVM approach. Here, we decompose the density matrix as

ρθ = b Π b θ M b , ( 4.66)
where { M b } is a set of tomographically complete measurement operators, with 4 L elements. A simple choice for this is the Pauli set, composed of strings of Pauli operators, i.e.

{ M b } = {1, σx , σy , σz } ⊗L . Π b θ = Tr{ρ θ M b }
is the probability for the operator string to be measured. As such, Π θ is a parametrized probability distribution over 4 L outcomes, and can be parametrized by a neural network, with a freedom on the choice of architecture. In [START_REF] Luo | Gauge equivariant neural networks for quantum lattice gauge theories[END_REF] and [START_REF] Reh | Time-dependent variational principle for open quantum systems with artificial neural networks[END_REF], the authors use a Transformer architecture and a recurrent neural network, respectively, to represent Π θ . Violating the positivity of the ansatz, although unphysical, has not proved to be such an issue in practical applications, much like for MPOs. However, we expect such approaches to work best for systems with strong dissipation. Consider a singe qubit system, whose state is described by a density matrix ρ which lives inside a Bloch sphere. If dissipation is strong, the state will tend to occupy central regions of the Bloch sphere, thus remaining in the space of valid quantum states. For states close to pure, thus living on the edge of the Bloch sphere, if positivity is not enforced the reconstructed state from Π θ may reside outside the Bloch sphere, i.e. outside the space of physical quantum states.

VII.2 Combining the dynamical corner-space method with NQS

As a perspective to be further investigated in the future, an alternative to the POVM and Gram-Hadamard matrix approaches is to employ the low-rank representation of the density matrix detailed in chapter 3. By parametrizing each state in the corner-space by a NQS and propagating the corner in time, that is now parametrized by a set of NQS, the time-evolution of the open quantum system is recovered. This approach is particularly suited to weakly dissipative systems, thus being complementary to POVM approaches. Its advantage with respect to Gram-Hadamard approaches is that it bears no constraint on the form of the ansatz, hence the entanglement structure of the final density matrix is the same as that of the states composing it. In this representation, the corner C(t) at each time t fully specifies the density matrix, and is defined by its elements as C sk = √ p k ψ θ k (s), with θ k the weights vector of the network that parametrises the k-th corner state. In the NQS language, the algorithm consists in four important steps: with updated weights η l (that form the transition basis).

(ii) Probabilities of each new state are obtained by using probability conservation; the probability π l of each state ϕ η l of the initial corner basis has to lead to D + 1 states whose probabilities sum to p l .

(iii) Now that the parameters η l represents the M (D + 1) states of the transition basis and π l their respective probabilities, one can construct T † T through a sampling procedure, diagonalize, and truncate the eigenvalues below ϵ to obtain the new probabilities pm that will correspond to M new states.

(iv) Obtain the new corner basis C from T . In the neural network picture, obtain the parameters θm that correctly represent the states of C from the eigenvalues and eigenvectors of T † T .

Below each step is detailed. Note that a quantum trajectory approach is also possible, where a single trajectory is propagated and subject to random jumps, with the same optimization procedure as the one described below.

Evaluation of corner states generated by jump operators

This parts amounts to minimizing, at each time step t, for every jump operator Ĵi , the distance

C η l := dist ϕ η l , Ĵi ψ k θ , ( 4.67)
with η l the weights of the ansätz to optimize (once optimized, these will be the parameters of the lth state of the transition basis), and θ k the weights of the states of the kth state of the corner basis at time step t. For ease of notation, we drop the superscripts l and k in the calculations that follow. For normalized ansätze, the cost function is simply the fidelity. We denote a configuration. For local jump operators Ĵi = 11 ⊗ 12 ⊗ . . . ⊗ ȷi ⊗ . . . ⊗ 1N (with ȷi a local operator acting on a single spin) we have:

C (η) = |⟨ϕ η | Ĵi |ψ θ ⟩| 2 ⟨ϕ η | Ĵ † i Ĵi |ϕ η ⟩ , ⟨ϕ η | Ĵi |ψ θ ⟩ = s,s ′ ϕ * η (s)ψ θ k (s ′ )⟨s| Ĵi |s ′ ⟩ = s |ϕ η (s)| 2 s ′ ψ θ k (s ′ ) ϕ η (s) ⟨s| Ĵi |s ′ ⟩ = s |ϕ η (s)| 2 s ′ i ψ θ (s ′ i |s\s i ) ϕ η (s) ⟨s i |ȷ i |s ′ i ⟩, ⟨ϕ η | Ĵi |ψ θ ⟩ = s |ϕ η (s)| 2 Ĵψ ϕ (s) = ⟨ Ĵψ ϕ ⟩ |ϕη| 2
where we have defined

Ĵψ ϕ (s) = σ ′ i ψ θ (σ ′ i |s\σ i ) ϕ η (s) ⟨σ i |ȷ i |σ ′ i ⟩,
and ⟨...⟩ |ϕ| 2 denotes an average over the probability distribution given by |ϕ(s)| 2 . The point here is that Ĵψ ϕ (s) can be evaluated exactly (in fact, for a k-local operator, it will be a sum of 2 k terms) and its average can therefore be approximated by sampling the probability distribution over which it is averaged. Therefore, one has:

C (η) = |⟨ϕ η | Ĵi |ψ θ ⟩| 2 ⟨ϕ η | Ĵ † i Ĵi |ϕ η ⟩ = ⟨ Ĵψ ϕ ⟩ |ϕη| 2 ⟨ Ĵϕ ψ ⟩ |ψ θ | 2 ⟨ Ĵ † i Ĵi ⟩ |ϕη| 2 . ( 4
∂ k C (η) = O k -⟨O k ⟩ |ψ θ | 2 Ĵϕ ψ |ψ θ | 2 ⟨ Ĵψ ϕ ⟩ |ϕη| 2 /⟨ Ĵ † i Ĵi ⟩ |ϕη| 2 (4.69)
which has the same computational complexity as that of minimizing the energy as presented at the beginning of the chapter.

Evaluate the probabilities

After step (i), one obtains the correct parameters that represent the states ϕ η l of the transition basis. We can almost build T † T to diagonalize it, but we still need the associated probabilities π l . In a regular corner approach, these would naturally arise from step (i). This is not the case here, so we have to obtain them in another way from the corner basis C. Using the facts that each state in Ĉ generates D + 1 states and that probability is conserved, we have:

D+1 i=1 π (k) i = p k ,
and we know that the states generated by the jumps (i ̸ = 1) will have a probability π

(k) i = δt⟨ψ θ k | Ĵ † i Ĵi |ψ θ k ⟩. It naturally follows that π (k) 1 = p k -D+1 i=2 π (k)
i , and we can therefore obtain the complete set of probabilities {π l } M (D+1) l=1 . This step is computationally for free since the averages ⟨ Ĵ † k Ĵk ⟩ were computed at the previous step. At this step, one can discard all the probabilities that are below some truncation error ϵ, and renormalize the probability distribution given by π l .

Evaluation of T † T

From the previous steps one obtains probabilities π l and the parameters η l corresponding to each state of the transition basis, ϕ η l . Since we cannot explicitly reconstruct T † T , we can sample it. Explicitly, the elements of T † T are

[T † T ] ll ′ = s ϕ * η l (s)ϕ η l ′ (s) √ π l π l ′ (4.70) = s |ϕ η l (s)| 2 ϕ η l ′ (s) ϕ η l (s) √ π l π l ′ (4.71) = ϕ η l ′ ϕ η l π l π ′ l |ϕ η l | 2 . ( 4.72)
We once again find the form of an expectation value. A difference here is that these expectation values are complex, hence may be more difficult to estimate due to a sign problem that comes from adding phases together. However there is no optimization here so it is possible to use a very large number of samples to reconstruct T † T with a low computational cost. After reconstructing T † T , one can diagonalize it and obtain M , the updated corner dimension, as well as pm , the corresponding probabilities. However, we do not obtain the set of parameters θm corresponding to these states, thus another optimization is needed. This step implies the diagonalization of an M (D + 1) square matrix, T † T , with no dependence on N .

Obtention of the updated corner basis C

In the original dynamical corner-space algorithm, one easily obtains C from T by multiplying the eigenstates of T by T itself. By denoting |t m ⟩ the m-th eigenvector of T † T (which was explicitly diagonalized in step 3) and t m q its q-th component, the elements of C are obtained as: Each state of the updated corner is therefore a particular superposition of the states of the transition basis for which we know the corresponding amplitudes √ πq √ pk t m q . By considering normalized ansätze, we are able to evaluate the correct amplitudes Ψ tar (s) and evaluate this distance and its gradient.

Cm = T |t m ⟩ , Csm = q T sq t m q , |ψ θ m ⟩ = q √ π q √ pm t m q |ϕ η q ⟩ .

Computational complexity

The computational complexity of such a method bears no dependence on N , the Hilbert space dimension. Combining the four steps of the algorithm, the complexity of a single time step is dominated by the first and third steps. The first step consists in M D optimizations that have a number of optimization steps n steps , where a gradient of size N par (the number of variaitonal parameters) is calculated by estimating averages over N samples samples. The third step consists in diagonalizing a M (D + 1) × M (D + 1) matrix. The complexity of a single time step is therefore given by O max [M (D + 1)] 3 , M Dn steps N par N samples .

(4.74)

VIII Conclusion

We have seen how neural-network quantum states presents itself as a powerful numerical toolbox to study the dynamics of complex quantum systems.

We have presented an alternative scheme for the real time evolution of quantum manybody systems with an NQS ansatz. This scheme does not rely on hyperparameter tuning for regularization and is found to be stable, and is not limited by the regime of the physical system. We therefore expect it to be useful when stochastic reconfiguration fails, when one uses autoregressive models or for when the network contains a large number of variational parameters. This scheme could be used to further investigate the limits of the Kibble-Zurek mechanism in nontrivial quantum systems, such as frustrated systems in triangular lattices or in systems with multiple critical points. As other applications, quantum control and pulse optimization problems could also be considered, where accurate methods to simulate time-dependent dynamics are required, as well as the quench dynamics of bosonic systems [START_REF] Carleo | Unitary dynamics of strongly interacting bose gases with the time-dependent variational monte carlo method in continuous space[END_REF]. As an outlook, the scheme could further be improved by considering different cost functions that lead to weaker optimization errors and variance. We also expect the scheme to be helpful for the dynamics of open quantum systems [START_REF] Hartmann | Neural-network approach to dissipative quantum many-body dynamics[END_REF][START_REF] Reh | Time-dependent variational principle for open quantum systems with artificial neural networks[END_REF][START_REF] Vicentini | Positive-definite parametrization of mixed quantum states with deep neural networks[END_REF], where t-VMC can present similar difficulties to those exposed. Finally, we expect the various approaches for the simulation of open quantum systems to be applied to noisy quantum circuits. Indeed, ideal quantum circuits have already been shown to be simulatable by NQS more efficiently than with tensor network methods [START_REF] Medvidović | Classical variational simulation of the quantum approximate optimization algorithm[END_REF]. In the regime of low noise, therefore weak entropy, a hybrid dynamical corner-NQS approach would be particularly suited. We also expect the Gram-Hadamard approach to be effective for the study of dissipative phase transitions, as well as the study of an extended form of the Kibble-Zurek mechanism to dissipative critical systems.

General conclusion and perspectives

In this manuscript, we explored the nonequilibrium dynamics of closed and open quantum systems and presented some novel methods to deal with the difficult problem of simulating such dynamics, along with new physical results for various interesting problems entering this class of systems.

After explaining the key questions and the main numerical methods surrounding the topic, we investigated the dynamics of entanglement in a strongly-interacting dissipative quantum system was presented. This chapter was motivated by recent experimental progress towards large chains of dissipatively stabilized Mott insulator of photons, which open up a new playground for many-body physics, where solid-state phases can be achieved in photonic systems. We introduced a new simple protocol to monitor entanglement dynamics, where bipartite entanglement can be investigated by looking at the negativity of reduced density matrices describing pairs of distant sites. We observe a ballistic propagation of entanglement, in agreement with the quasiparticle picture provided by a closed system analysis. However, a strikingly asymmetric behaviour is observed as a function of the chosen protocol, which depends on whether we initially inject or extract a particle.

Then, we presented the dynamical corner-space method. This new method to simulate the dynamics of open quantum systems fills up a gap in the literature that existed for low-entropy, highly-entangled states. If the rank of the density matrix describing the system of interest is small with respect to the Hilbert space dimension, the method can outperform other methods that do not restrict entanglement, such as Monte-Carlo wave function methods. We applied this new method to the simulation of a noisy quantum Fourier transform processor by considering a piecewise continuous-time implementation of the algorithm, so that microscopic dissipative channels were taken to account as opposed to digital error models. By doing so, we were able to derive scaling laws for the error as a function of the dissipation rate. In addition, the impact of the initial state on the performance of the algorithm was unveiled, revealing a highly biased algorithm. We believe this method will enable the study of dissipatively stabilized quantum systems such as cat qubits, and error mitigation protocols that are strongly needed in the quantum information community until the era of fault-tolerant quantum computation is reached.

Finally, the general class of ansätze of neural-network quantum states (NQS) was presented. These were introduced in the context of the general task of generative modeling in standard machine learning. After presenting Monte-Carlo sampling algorithms, we introduced the various architectures used to represent quantum states. In particular, autoregressive models, which enable for normalized quantum states and efficient sampling algorithms were exposed. We believe such architectures, that are becoming standard in state-of-the-art machine learning routines, will become a key element of the numerical physics toolbox to study complex quantum systems. We presented various issues that arise when trying to simulate quantum dynamics with NQS, as well as an alternative scheme to perform this task that we developed. This scheme allows one to study timedependent quenches, and the transition from adiabatic to genuine nonequilibrium dynamics that occur during such quenches, as explained by the Kibble-Zurek mechanism. We also presented the issue of representing the state of an open quantum system with NQS, whose difficulty resides in preserving the positivity of the density matrix. We displayed key ingredients of a new hybrid method, which uses ideas from the dynamical corner-space method combined with NQS. We believe such a method will be particularly useful for the simulation of noisy quantum circuits, to reach sizes well beyond what was achieved using the original dynamical corner-space method and examine, for example, the entanglement dynamics and spreading that occur in noisy quantum circuits.

Very recently, a new efficient algorithm for matrix multiplications was discovered by a machine learning model, coined AlphaTensor [START_REF] Fawzi | Discovering faster matrix multiplication algorithms with reinforcement learning[END_REF]. Prior to this discovery, the last algorithm for such a task was about fifty years old. We believe such phenomena will only amplify with the increase of computing power along with the improvement of algorithms. The work presented in this thesis is part of the efforts to design and improve effective models, to ultimately achieve an era of computationally-assisted scientific discovery [START_REF] Krenn | On scientific understanding with artificial intelligence[END_REF], where the task of the scientist becomes that of formulating a problem as precisely as possible, and algorithms provide solutions to these problems, as complex as they may be.

A High-order integration methods for (A.4)

NQS I Generalization to higher-order methods

Note that the entries are nonzero only on the lowest triangular matrix for explicit methods.

We will here consider methods whose Butcher coefficients satisfy either a i,j>i = 0 (explicit) or a i,j≥i = 0 (implicit).

For the considered tableaus, the system of Eqs. (A.2) and (A.3) can be solved explicitly by using the following recurrence identity:

|κ (n) i ⟩ = Π(n) i |ψ (n) ⟩ + h n i-1 j=1 Π(n) i a ij |κ (n) j ⟩, (A.5) with Π(n) i = 1 -a ii h n Φt n,i -1 Φt n,i . (A.6)
Indeed, one has:

|κ (n) i ⟩ = Π(n) i + h n Π(n) i a ij Π(n) j + h 2 n Π(n) i a ij Π(n) j a jk Π(n) k + . . . |ψ (n) ⟩, (A.7)
and thus |ψ (n+1) ⟩ = T (n) s |ψ (n) ⟩, (A.8)

with

T (n) s = 1 + s i=1 b i h n Π(n) i + h 2 n Π(n) i a ij Π(n) j + h 3 n Π(n) i a ij Π(n) j a jk Π(n) k + . . . . (A.9)
Representing the wavefunctions above with a variational ansatz, namely σ ψ η n = ψ η n (σ), the generic update can finally be recast into the following optimization process:

|ψ η n+1 ⟩ = arg min |ψη⟩ dist |ψ η ⟩, T (n)
s |ψ η n ⟩ . (A.10)

Time-independent explicit case

The equations above considerably simplify when considering time-independent Hamiltonians and an explicit integration method. Indeed, to any order s, we have For these, we have:

T (n) s = s-1 m=0 λ m h m n Π(n)m =
λ 0 λ 1 λ 2 λ 3 1 1 1/2 1/6
This corresponds to the usual factor λ m = 1/m! of the truncated Taylor expansion of the propagator. Note, however, that this is generally no longer the case in a time-dependent scenario, as will appear below.

B Résumé substantiel

Cette thèse est consacrée à l'étude de la dynamique hors équilibre dans les systèmes quantiques idéaux (fermés, isolés de leur environnement direct) et dissipatifs (ouverts, couplés à leur environnement direct). Après une présentation générale de questions liées à la dynamique de systèmes physique ainsi que de leur Histoire, une introduction détaillée du sujet ainsi qu'une exposition des diverses difficultés apparaissant lors de son étude sont faites dans le chapitre 1. Ensuite, un premier objet d'étude, la dynamique de l'intrication dans des systèmes photoniques dissipatifs fortement corrélés, est exposé dans le chapitre 2 ainsi que les divers nouveaux résultats physiques obtenus par une approche de réseaux de tenseurs. Puis, la méthode du sous-espace dynamique, méthode générale pour la dynamique de systèmes ouverts à entropie faible est présentée dans le chapitre 3. Son application à l'algorithme de la transformée de Fourier quantique est présentée, encore une fois accompagné de nouveaux résultats obtenus grâce à cette nouvelle méthode. Enfin, la simulation de la dynamique de systèmes isolés grâce à une représentation en états quantiques neuronaux est détaillée, ainsi qu'une nouvelle méthode développée pour obtenir des résultats plus fiables pour les systèmes avec un Hamiltonien dépendant du temps sont présentés dans le chapitre 4. Divers résultats concernant le mécanisme de Kibble-Zurek dans des systèmes d'Ising avec champ transverse dans la région critique ont ainsi pu être obtenus.

Un système photonique fortement corrélé composé de L résonateurs électromagnétiques couplés peut être décrit par l'Hamiltonien suivant l'opérateur de saut pour le site i et le canal de dissipation C. Nous considérons deux canaux : le canal de dissipation d'énergie, et le canal de déphasage. En utilisant une approche par réseaux de tenseurs (MPO), nous avons obtenus de nombreux résultats sur la dynamique de l'intrication. Nous considérons deux cas : dans le premier, nous ajoutons une particule au site central, au dessus d'un état type Mott fortement corrélé. Dans le deuxième cas, nous considérons le retrait d'une particule (donc l'ajout d'un trou) sur le site central. Nous pouvons estimer l'intrication grâce à la négativité, qui est une mesure de l'intrication bipartite. Nous pouvons donc quantifier l'intrication entre deux sites i et j. En faisant une trace partielle sur l'ensemble du système privé des sites i, j, nous obtenons une matrice densité réduite, pour laquelle nous calculons la négativité. Ainsi, nous obtenons que la propagation de l'intrication est ballistique aux temps courts, et est correctement prédite par les calculs pour le système fermé. Cela indique que l'image de quasiparticules reste pertinente même lorsque la dissipation est forte. Remarquablement, la dissipation affecte beaucoup plus la propagation de l'intrication dans le cas du retrait d'une particule au temps initial. Cela s'explique par une image de trajectoires quantiques, qui rend compte de la fragilité de la propagation du trou.

Ĥ = L i=1 ω c b † i b i + U 2 b † i b † i b i b i -J L-
Le chapitre 3 est dédié à la simulation numérique de systèmes quantiques ouverts en général, et en particulier des systèmes ayant une entropie limitée. Ce type de systèmes est celui généralement construit pour les tâches de traitement de l'information quantique, ou les systèmes sont isolés le plus possible de leur environnement direct. Cette technique est basée sur la décomposition de la solution de l'équation maîtresse ρ(t) sous la forme, pour tout instant t :

ρ(t) ≃ M (t

|ψ⟩ = s 1

 1 ,s 2 ,...,s L c s 1 ,s 2 ,...,s L |s 1 , s 2 , . . . , s L ⟩ (1.1)

Figure 1 . 3 :

 13 Figure 1.3: (a): Phase diagram of the model described by the Hamiltonian in Eq.(1.21).One is mostly interested in the antiferromagnetic region (AFM), whose eigenstates can be probed by suddenly quenching the detuning from deep in the paramagnetic phase, with ∆ ≪ J. (b) Scheme of the experimental setup: a 2D array of atoms trapped in an optical lattice whose states |e⟩ (large red spheres) |g⟩ (small blue spheres) can be coupled with an infrared laser. The top part of the scheme represents the high-resolution objective that can resolve individual sites, enabling precise analyses of experimental data. (c): Typical protocol. The system is prepared in a trivial phase, with a small number of defects (left), and is then quenched into an excited configuration (center). Then, increasing the lattice depth causes the atoms in the excited state to be lost (right), hence only the atoms in the ground states remain to simplify detection. Reproduced from[START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF].

Figure 1 . 4 :

 14 Figure 1.4: Correlation functions after a sudden quench in the two-dimensional antiferromagnetic Ising model as a function of the detuning ∆, for different values of the quench time T . Dashed lines correspond to linked-cluster calculations, while solid lines correspond to exact diagonalization data. Upper panels (a), (b), (c) and (d) correspond to time ΩT = 0.5π, which is fast enough to be considered a sudden quench, while lower panels (e, (f), (g) and (h) correspond to ΩT = 2.97π, where long-range correlations develop. Reproduced from[START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF].

Figure 1 . 5 :

 15 Figure 1.5: Experimental results for a time-dependent linear quench on a Rydberg quantum simulator. (a):Single-shot images of the system before the quench (top row), after a fast quench (medium row) and after a slow quench (bottom row). The correlated domains are larger in the bottom row, as more time is spent in the impulse region hence domains have more time to develop. A higher density of defects is observed for fast quenches. (b) Correlation length as a function of the final detuning. The grey dashed line indicates the critical detuning: if the quench is stopped before this, the correlation length remains small, as the dynamics remain quasi-adiabatic. (c) Scaling of the correlation length as a function of the sweep rate, which is the inverse of the quench time τ q , compared with results from an MPS numerical integration. Reproduced from[START_REF] Keesling | Quantum Kibble-Zurek mechanism and critical dynamics on a programmable rydberg simulator[END_REF].

Figure 1 . 6 :

 16 Figure 1.6: Sketch of the corner-space renormalization method. Taken from [65].

2 .

 2 Merge the two sub-systems into a larger system A ∪ B. Construct a basis C(M ), spanned by product states of the form ϕ (A) r ϕ (B) r ′ , and retain only the M most probable ones by keeping only the largest joint probabilities p (A) r p (B) r ′ . Solve the Lindblad master equation in this subspace of dimension M .

Figure 1 . 7 :

 17 Figure 1.7: (a): Optical image of a driven-dissipative superconducting circuit where a chain of transmons, labeled from 1 to 8, are coupled to a reservoir and to a stabilizing drive. Observables can be measured by the readout resonators. The inset to the figure is a scan-electron microsopy image of a single transmon qubit. (b) Energy scales in such a system, where the interaction dominates all other energy scales. (c) Sketch of the effective model obtained, where a strongly-interacting Bose-Hubbard model is realized, coupled to two sites that act as a stabilizer for Mott insulator states. Reproduced from [69].

33 Figure 2 . 1 :

 3321 Figure2.1: Upper panel: sketch of the considered system, a chain of coupled photonic resonators with on-site photon-photon interaction U and nearest-neighbor hopping coupling J. The top chain depicts the initial time configuration with a Mott insulator of photons (one photon per cavity) where a double occupation (doublon) has been created in the central site i = i c . The bottom chain depicts the configuration at a later observation time t, with entanglement existing between distant sites. The photonic modes are subject to losses and dephasing. Lower panel: entanglement negativity N r (t) between sites i c-r and i c+r as a function of time t for different values of the spatial separation r = {1, 2, 3, 4} from the central site i c . The shade of the lines gradually decreases going from r = 1 to r = 4. Calculations were performed with an MPO ansatz with bond link dimension χ = 200 on a chain of L = 20 cavity sites. For each value of r, results for different values of the photon loss rate γ in units of the hopping J are shown. The initial state at t = 0 is |Ψ D ⟩ (see the text) corresponding to a doublon excitation localized in the central site. In this figure, the pure dephasing rate Γ d is 0. The on-site interaction for all the cavities is U/J = 33.3.

33 Figure 2 . 2 :

 3322 Figure 2.2: Upper panel: sketch like in Fig. 2.1, but where the initial state has an empty central site (holon). Lower panel: temporal dynamics of the negativity with same parameters as in Fig. 2.1, but with an initial state |Ψ H ⟩ corresponding to a holon localized in the central site.

Figure 2 . 3 :

 23 Figure2.3: Entanglement propagation speed (units of J) versus the normalized hopping J/U . Filled (blue) circles: entanglement speed for the configuration corresponding to the injection of one additional photon (doublon injection). Empty (red) circles: entanglement speed for the holon injection. Solid (blue) thick line: maximal speed for doublon quasiparticles in a closed Hamiltonian system. Thin (red) solid line: maximal holon speed. Dashed lines are polynomial fits with linear and quadratic terms in J/U . Error bars were estimated taking into account uncertainty due to time discretization and to the finite bond link dimension in the MPO calculations. Parameters: γ/J = 0.1 and for the doublon (holon) configuration U/γ = 100 (1000).

Figure 2 . 4 :

 24 Figure 2.4: Peak value of N r=1 (log scale) as a function of the photon dissipation rate γ or dephasing Γ d . Thick (thin) solid line: curve for the doublon (holon) configuration versus γ when only photon particle losses are present (Γ d = 0). The holon injection case has been calculated for a chain with a hopping coupling 2J to have the same maximal speed of a doublon in a chain with hopping coupling J (see text). Dashed (doublon) and dotted (holon) lines: negativity peak versus Γ d when there is only the pure dephasing channel (γ = 0). Parameters: L = 5, U/J = 100.

Figure 3 . 1 :

 31 Figure 3.1: (a) Sketch depicting one iteration of the time-dependent corner-space method.

Figure 3 . 2 :

 32 Figure 3.2: (a) Memory use of different quantities.One can see that the gain in memory requirement from our method scales exponentially. Note tha the Hamiltonian, being very sparse, is not the limitation. The largest stored object is the corner Ĉ, that is composed of M states. Here, as an example, we considered M = 100. (b) Infidelity between the dissipative output density matrix ρ(c) γ (calculated via the corner-space method) and the exact integration density matrix ργ as a function of the maximum corner dimension M max , for an L-qubit circuit and γT QFT = 1.42 × 10 -2 . The initial state is |ψ 0 ⟩ = QFT -1 |GHZ⟩.

Figure 3 . 3 :

 33 Figure 3.3: Wigner function of the 4 leading components of the corner-space describing a cavity subject to a Lindblad evolution as described in section III.5, after time γt = 10. The first two correspond to states |C ± α ⟩ with opposite parities. Parameters here are κ/γ = 2, ω c /γ = 1, K/γ = 10, G/γ = 50.

Figure 3 . 4 :

 34 Figure 3.4: (a) Quantum circuit representing the quantum Fourier transform (QFT) in the presence of dissipation. (b) Continuous-time evolution of the exponential of the von Neumann entropy S (dash-dotted) with the input state |ψ 0 ⟩ defined in Eq. (3.17), for L = 8 qubits and γ/δ = 1 × 10 -3 . The corner-space dimension M (t) is also plotted (plain line). Temporal intervals corresponding to Hadamard and controlled-phase gates are indicated by lighter (red) and darker (blue) background colors, respectively. (c) Temporal build-up of the entanglement entropy S ent (n), as defined in Eq. (3.19), as the different gates of the QFT are performed.

3 . 4 ,

 34 the time evolution of the corner-space dimension M as well as the exponential exp(S) of the von Neumann entropy S(ρ) = Tr[ρ ln ρ] are shown, exhibiting similar trends. Panel (c) depicts the spatial entanglement propagation as the circuit's gates are progressively applied from the first to the last qubit. This is quantified by the entanglement entropy S ent (n) ≡ S(Tr 1,...,n [ρ]), (3.19)

γTQFT = 1 . 5 •Figure 3 . 5 :

 1535 Figure 3.5: (a) Infidelity between the output density matrix ρ(c)γ calculated with the corner-space method and the output ργ obtained via an exact integration as a function of the maximum corner dimension M max for different values of γT QFT and a fixed number of qubits L = 10. (b) Fidelity as a function of M max and L for γT QFT = 2.5 × 10 -2 . The parameter N denotes the dimension of the total Hilbert space. Values corresponding to M max = L are highlighted by hollow markers (we find a fidelity F ≳ 0.997 for M max ∼ L ln L). The initial state is |ψ 0 ⟩ as in Fig.3.1 for all three panels. (c) Simulation time of the QFT algorithm in the presence of dissipation versus the number L of qubits for the exact solution of the master equation (squared markers) and for the time-dependent corner-space method (circles) for two different values of the control parameter ϵ. The dissipation rate is set by γT QFT = 2.5 × 10 -2 .

1 Figure 3 . 6 :Figure 3 . 7 :

 13637 Figure 3.6: (a) Infidelity between the noisy output density matrix ρ (c) γ and the noiseless output ρ0 = |GHZ⟩⟨GHZ| obtained from the initial state |ψ 0 ⟩ for different values of γ/δ. The dashed line is a guide to the eye showing a growth ∝ L 2 . (b) Infidelity for three different initial states, described in the legend, for the same dissipation rate γ/δ = 2.5 × 10 -4 .

(3. 23 )Figure 3 . 8 :

 2338 Figure 3.8: Absolute value of the elements of the error matrices χ err (normalized by the characteristic error magnitude γτ ) after the continuous-time evolution described in section V. In panel (a), the system is subject to dissipation ( Ĵi = √ γσ - i ) and in panel (b) to decoherence channels ( Ĵi = √ γσ z i ). Here γ/δ = 10 -3 . For each error matrix the element χ 11 has been set to zero for readability (for such weak values of γ, χ 11 ∼ 1).

Figure 4 . 1 :

 41 Figure 4.1: Pictorial representation of a RBM that represents the wave function of an Lspin system, with s = (σ 1 , σ 2 , . . . , σ L ) and h = (h 1 , h 2 , . . . , h d ) the hidden units. Adapted from[START_REF] Dawid | Modern applications of machine learning in quantum sciences[END_REF].

Figure 4 . 2 :

 42 Figure 4.2: Example of an ARNN quantum state for four spins. (a) Pictorial representation of the network. The arrows representing the weights of the model are skewed in order not to break the conditional structure of the output probability distribution. These layers are "masked", due to some connections being deleted. (b) Sampling algorithm.One samples consecutive spins using direct sampling on the conditional probabilities at each step. From[START_REF] Dawid | Modern applications of machine learning in quantum sciences[END_REF].

Figure 4 . 3 :

 43 Figure 4.3: Pictorial representation of a recurrent neural network architecture for NQS.Panel (a) is for real-valued wave functions, which can be relevant for a certain class of problems, and panel (b) is for complex-valued wave functions. In both schemes, a local spin configuration s i and a hidden vector h i are fed into an RNN cell, which performs a nonlinear transformation. Then an activation function (ς, for softmax and/or ςς, for softsign) is applied to obtain the final probability and/or phase corresponding to the configuration. At the end, the probabilities (and phases) are combined to obtain the final wave function amplitudes ψ(s).

Figure 4 . 4 :

 44 Figure 4.4: Expressive capacity of different classes of variational states, as explicitly proven in Ref.[START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF] by mapping tensor networks to NQS. PEPS* refers to a sub-class of projected entangled pair states, a generalization of MPS for two-dimensional systems. Adapted from[START_REF] Sharir | Neural tensor contractions and the expressive power of deep neural quantum states[END_REF].

Figure 4 . 5 :

 45 Figure 4.5: Schematic representation of various ansätze inspired by [217]. (a) The Jastrow ansatz draws connections between all possible pairs of sites. (b) The MPS ansatz draws connections between nearest-neighbour sites along a line. (c) The RBM ansatz connects all the sites to every hidden neuron, illustrated in different colours. (d) The RNN ansatz processes the state sequentially, following the green arrows. The dark blue arrows indicate the flow of information within the model. Arrows without a starting site correspond to free parameters. Sketch inspired by [218].

Figure 4 . 6 :

 46 Figure 4.6: (a): Scheme of the GRU ansatz for a two-dimensional system of L = 4 spinswith open boundary conditions. The system is subject to the TFI Hamiltonian, and configurations are fed into GRU in an autoregressive fashion, each cell being fed with a hidden vector h i,j and a spin σ i . After the GRU cell, a softmax (ς) and a softsign (ςς) transformation are applied to obtain the amplitudes and the phases of the conditional amplitudes ϕ i (σ i |σ i-1 , .., σ 1 ). (b) Spectra of the quantum geometric tensors S for different ansätze optimized to the ground state of the TFI model, for various values of g/J. The left panel is for a complex RBM ansatz with α = 2, the central panel for a two three layer convolutional autoregressive network with complex parameters (ARNN), and the right panel for a GRU ansatz with real parameters with d h = 10, for which is plotted the spectrum of the imaginary part of the S matrix, as this is the quantity to invert for ansätze with real parameters.

. 43 )Figure 4 . 7 :

 4347 Figure 4.7: (a) Scheme of the Kibble-Zurek mechanism according to the adiabatic-impulse approximation. At times ± t, the relaxation time of the system t R crosses the time scale of the quench, which defines the limit between the adiabatic region, where the system essentially stays in the ground state of the system at that time, and the impulse regime, where nonequilibrium dynamics take place. (b) Scaling of the injected energy Ĥ(τ q ) -E 0 as a function of the quench time τ q for different system sizes N = L × L. One can see that for fast quenches, the NQS dynamics predict very well the injected energy. The data for system sizes with more than L = 5 cannot be simulated via exact diagionalization, hence they are compared to the τ 0.72

  d h = 8, with N s = 200 samples. (d) Average correlation function along the z-axis for the same quench as in (c). with ⟨ψ θ| Ts |ψ θ ⟩ ≡ ⟨ψ θ| Tloc |ψ θ⟩ , (4.46)

Figure 4

 4 Figure 4.8: Numerical results obtained for a linear quench for various values of the total quench time τ q . On the upper panels, the average correlation function is shown as a function of time for both the exact (dashed line) and the GRU simulations (solid lines). The error as a function of time ϵ(t) is shown in the lower panels. For panels (b) and (c), the time step was chosen to be linear and non-constant, as the characteristic time scale changes as a function of time. Parameters for each quench are are given in the appendix.

Figure 4 . 9 :

 49 Figure 4.9: Sudden quench dynamics induced by quenching the Hamiltonian parameters to g/g c = 2, g/g c = 1 and g/g c = 1/10 for various hidden vector dimensions d h , starting from an initial paramagnetic state. Upper panels show the average magnetization along the x axis ⟨ Ŝx ⟩ = i σx i /L 2 and lower panels show the cumulative error as a function of time.

Figure 4 .

 4 Figure 4.10: (a): residual error ϵ as a function of the time step δt for the Euler method, the midpoint method, and the Heun integration methods averaged over the 10 first time steps of quench dynamics of the TFI model to g/g c = 2. (b) Standard deviation of the error as a function of time rescaled by √ N s .Curves collapse, which is what is predicted by a direct sampling scheme. This enables us to estimate the required number of samples to obtain a given error. Note that as the dynamics progress, the variance of the residual error rises, indicating a rising difficulty to find the best minima in time.

Figure 4 . 11 :

 411 Figure 4.11: Tensor network diagram of the Gram-Hadamard ansatz for the density matrix.The plaquette structure of the ansatz is shown: with K plaquettes of rank R, the obtained density matrix has rank R K , which enables for states that have a high entropy to be efficiently represented. Adapted from[START_REF] Vicentini | Positive-definite parametrization of mixed quantum states with deep neural networks[END_REF].

  (i) Evaluate the M (D + 1) states generated by the Kraus operators Ki acting on the M states of Ĉ. The state generated by Kraus operator K0 can be obtained as for closed systems, i.e. using t-VMC with a non-Hermitian Hamiltonian or an optimisation scheme as presented in the previous section. This step consists of O(M D) optimizations. This leads to M (D + 1) states ϕ η l M (D+1) l=1

Ψ tar ψ θ m |ψ θ m | 2 ,tar ψ θ m |ψ θ m | 2 ψ θ m Ψ tar |Ψtar| 2 .

 222 By noting |Ψ tar ⟩ the target state, one can therefore minimizedist (|ψ θ m ⟩ , |Ψ tar ⟩) = |⟨ψ θ m |Ψ tar ⟩| 2 .In analogy with the previous steps, one has:⟨ψ θ m |Ψ tar ⟩ = s |ψ θ m (s)| 2 Ψ tar (s) ψ θ m (s) = hence dist (|ψ θ m ⟩ , |Ψ tar ⟩) = Ψ

  b T a m-2 c, else. (A.12) Two common Butcher tableaus corresponding to the fourth-order Runge-Kutta method are

1 i=1

 1 (b † i b i+1 + h.c.), (B.1)avec ω c la fréquence du mode de la cavité, U la force d'interaction photon-photon (Kerr), J le couplage entre sites voisins, et b i (b † i ) l'opérateur d'annihiliation (création) sur chaque site. Un tel système présente une transition de Mott entre un système fortement corrélé, lorsque U ≫ J et une phase superfluide dans le régime à interactions faibles U ≪ J). Lorsque les résonateurs sont dissipatifs, le système peut être décrit par l'équation de Lindblad suivante, sous les approximations de Born et Markov :

  It appears from this expansion that the complete description consists in l 2L complex numbers to be stored. Computationally, this means that a naive approach would enable to describe open quantum systems of half the size of the largest closed systems we can simulate. Conveniently, many approaches exist that are tailored for open quantum systems to describe the properties of open quantum systems without storing the full density matrix.

	s ′ 1 ,s ′ 2 ,...,s ′ L s 1 ,s 2 ,...,s L	c s ′ 1 ,s ′ 2 ,...,s ′ L s 1 ,s 2 ,...,s L |s 1 , s 2 , . . . , s L ⟩ ⟨s ′ 1 , s ′ 2 , . . . , s L | .	(1.2)

  There exists vectors {|u 1 ⟩ , . . . |u m ⟩} ∈ H A and {|v 1 ⟩ , . . . |v n ⟩} ∈ H B which are pairwise orthogonal ⟨u i |u j ⟩ = δ ij and ⟨v i |v j ⟩ = δ ij and a real numbers Λ 1 ≥ Λ 2 ≥ . . . Λ m ≥ 0 such that |ψ⟩ can be written as a Schmidt decomposition:

  .[START_REF] Josephson | Possible new effects in superconductive tunnelling[END_REF]) C(θ) can be minimized via a stochastic gradient descent, exactly like the cost functions considered in the chapter. The components of the gradient ∇C(θ) are again given by the log-derivative trick, with O k = ∂ θ k ln ψ

θ (s):

  Let us consider some linear ordinary differential equation of the form∂ t |ψ(t)⟩ = Φt |ψ(t)⟩ , (A.1)whose solution |ψ⟩ is discretized over a set of well-defined times {t n } N n=0 such that |ψ (n) ⟩ ≡ |ψ(t n )⟩. Linear Runge-Kutta (RK) methods provide us with the following prescription for approximately updating |ψ (n) ⟩: |ψ (n+1) ⟩ = |ψ (n) ⟩ + h n denotes the number of stages of the method, h n = t n+1 -t n , t n,i = t n -c i h n and the coefficients {a ij } i,j , {b i } i and {c i } i are completely determined by the Butcher tableau of the considered method:c 1 a 11 a 12 • • • a 1s c 2 a 21 a 22 • • • a 2s s a s1 a s2 • • • a ss b 1 b 2 • • • b s .

			s
				b i |κ (n) i ⟩,	(A.2)
			i=1
	|κ (n)
	. . .	. . .	. . .	. . .	. . .

i ⟩ = Φt n,i |ψ (n) ⟩ + h n s j=1 a ij |κ (n) j ⟩ , (A.3)

where s c

)

  k=1 p k (t) |ϕ k (t)⟩⟨ϕ k (t)| , p k (t) ≥ p k+1 (t), ∀k, (B.3) où les probabilités {p k (t)} k sont les M (t) valeurs propres les plus élevées à l'instant t et {|ϕ(t)⟩ k } k sont leurs vecteurs propres associés. Lorsque l'entropie du système reste faible au cours de son évolution temporelle, cette décomposition peut être tronquée à un nombre de composantes M bien plus faible que la dimension de l'espace de Hilbert associé au système. La méthode du sous-espace dynamique permet ainsi de limiter l'évolution temporelle à ce nombre très limité de composantes sans jamais reconstruire la matrice densité complète. S'ensuit un énorme avantage sur le plan de la complexité numérique de la simulation. À titre d'exemple, pour un ensemble de L spins 1/2, la complexité est réduite de O(2 L × 2 L ) à O(M × 2 L ). Cette méthode est appliquée à un algorithme réalisant la transformée de Fourier quantique en présence de dissipation ou de déphasage. Cet algorithme est au centre de nombreux algorithmes quantiques, comme l'algorithme de Shor, et en comprendre les mécanismes d'erreur avec une modélisation réaliste est ainsi de grande importance. Nos simulations montrent que l'erreur introduite par le couplage du circuit à un environnement extérieur croît quadratiquement avec la taille du système L. Nous montrons aussi que l'erreur dépend sensiblement du type d'état à partir duquel l'opération est menée. Nous envisageons aussi d'étendre ce type de méthode afin de pouvoir utiliser des états variationnels comme états composant le sous-espace, ce qui en augmentera considérablement sa compléxité numérique. Enfin, nous présentons le cadre théorique de la représentation d'états quantiques par des réseaux de neurones. Pour une fonction d'onde ψ dans un espace de Hilbert discret, l'amplitude de chacune des configurations s peut être représenté par un réseau de neurone, auquel l'on donne en entrée la configuration, et en sortie nous obtenons ψ(s). L'application canonique d'une telle méthode est trouver l'état fondamental d'un système décrit par un Hamiltonien Ĥ. En minimisant la quantité ⟨ψ θ | Ĥ|ψ θ ⟩ ⟨ψ θ |ψ θ ⟩ (B.4) sur l'ensemble des valeurs possible de θ, nous pouvons obtenir une estimation de l'état fondamental et de son énergie, en vertu du principe variationnel. En utilisant des méthodes standard d'apprentissage profond, cette quantité peut être minimisée et efficacement, et ainsi ces méthodes deviennent très utiles pour les problèmes fortement corrélés, ou les autres méthodes numériques à notre disposition sont très limitées. En particulier, nous étudions le mécanisme de Kibble-Zurek pour des systèmes d'Ising avec champ transverse. Ce système présente une transition de phase entre une phase paramagnétique et ferromagnétique (ou antiferromagnétique). Proche de la transition de phase, il est notoirement difficile d'étudier ce système, notamment lorsque l'Hamiltonien dépend du temps comme lors de dynamique type Kibble-Zurek. En considérant une intégration numérique alternative avec des états quantiques neuronaux, nous sommes capable d'obtenir de résultats pour ce mécanisme cohérent avec d'autres méthodes numériques, pour ainsi corroborer de tels résultats. Nous montrons aussi les limites d'autres méthodes d'intégration numériques pour les états neuronaux, comme les méthodes de Monte Carlo variationnelles en temps (t-VMC). Ces méthodes pourront donc être utilisées à l'avenir pour investiguer ce mécanisme dans des systèmes plus complexes, voire dissipatifs.

Note that there exists a whole class of Lieb-Robinson bounds, this is just the most generic.

Correlation functions appear in the expansion of the norm of the commutator [ Â(t), B].

with respect to the ground state at the end of the quench. In a 1D classical spin system, this will be proportional to the number of domains with spins flipped up if all spins point down in the ground state.

We have added a minus sign to the Hamiltonian in ref.[START_REF] Guardado-Sanchez | Probing the quench dynamics of antiferromagnetic correlations in a 2D quantum Ising spin system[END_REF] for the sake of simplicity.

Naively, enforcing positivity can be done by diagonalizing the full matrix, but the method then stops to be computationally efficient.

unless P=NP.

Chapter 2. Entanglement dynamics in strongly-correlated open quantum systems

For fully-connected operators of the form K = â1 ⊗ b2 ⊗ • • • ⊗ ẑL , this would be an issue, but such error models are quite unrealistic.

This corresponds to the violation of the classical identity S A ≤ S A⊕B .

In evaluating these execution times, both the exact integration and corner-space calculations were carried out on a single six-core Intel Xeon E5-2609 v3 processor at 1.9 GHz.

We will come back to question in section V.

If A loc is Gaussian, then the error is equal to σ 2 /M .

When |ψ θ ⟩ approaches an eigenstate of Â, it can be shown that the variance vanishes. Consequently, considering  = Ĥ, the statistical error vanishes as we approach the ground (or any excited) state. This is known as the zero-variance property, and is convenient for ground state search as it means the variance of the error will vanish as the solution is approached, thus ensuring a well-behaved training procedure.

We can use the intermediate conditional probabilities to draw samples for a low computational cost, e.g., use the probabilities for N -1 spins and sample from the last one, to obtain new samples; with Markov chain Monte Carlo we cannot do this.

Alternatively, one can also perform the inversion with iterative methods such as conjugate gradients, thereby reducing the complexity. This, however, does not enable all forms of regularization.

or anti-ferromagnetic if J < 0. For simplicity we will consider J > 0 throughout the paper.

Acknowledgements

Love and work are the cornerstones of our humanness.

-Sigmund Freud, apparently misattributed, ca. 1910

Time-dependent implicit midpoint method

The implicit midpoint method is characterized by the tableau

and yields

Time-dependent Heun method

The Butcher tableau of this explicit second-order Runge-Kutta method reads

and yields the following propagator:

(A.17)

Time-dependent Ralston method

The Butcher tableau of this explicit second-order Runge-Kutta method reads 0 2/3 2/3 1/4 3/4 (A.18)

and yields the following propagator: .19)