IBM Q Experience as a versatile experimental testbed for simulating open quantum systems

The advent of noisy intermediate-scale quantum (NISQ) technology is changing rapidly the landscape and modality of research in quantum physics. NISQ devices, such as the IBM Q Experience, have very recently proven their capability as experimental platforms accessible to everyone around the globe. Until now, IBM Q Experience processors have mostly been used for quantum computation and simulation of closed systems. Here, we show that these devices are also able to implement a great variety of paradigmatic open quantum systems models, hence providing a robust and flexible testbed for open quantum systems theory. During the last decade an increasing number of experiments have successfully tackled the task of simulating open quantum systems in different platforms, from linear optics to trapped ions, from nuclear magnetic resonance (NMR) to cavity quantum electrodynamics. Generally, each individual experiment demonstrates a specific open quantum system model, or at most a specific class. Our main result is to prove the great versatility of the IBM Q Experience processors. Indeed, we experimentally implement one and two-qubit open quantum systems, both unital and non-unital dynamics, Markovian and non-Markovian evolutions. Moreover, we realise proof-of-principle reservoir engineering for entangled state generation, demonstrate collisional models, and verify revivals of quantum channel capacity and extractable work, caused by memory effects. All these results are obtained using IBM Q Experience processors publicly available and remotely accessible online.


INTRODUCTION
The theory of open quantum systems studies the dynamics of quantum systems interacting with their surroundings. [1][2][3] In its most general formulation, it allows us to describe the out-of-equilibrium properties of quantum systems, it provides a theoretical framework to assess the quantum measurement problem, and it gives us the tools to investigate, understand and counter the deleterious effects of noise on quantum technologies. For these reasons its range of applicability is extremely wide, from solid state physics to quantum field theory, from quantum chemistry and biology to quantum thermodynamics, from foundations of quantum theory to quantum technologies.
Generally, the dynamics of open quantum systems are described in terms of a master equation, i.e., the equation of motion for the reduced density operator describing the quantum state of the system. Master equations are either phenomenologically postulated or derived microscopically from a Hamiltonian model of quantum system plus environment. Contrarily to the case of closed quantum systems, where the equation of motion describing the state dynamics is the Schrödinger equation, the general form of the master equation for an open quantum system is not known. Only under certain assumptions, known as the Born-Markov approximation, one can derive a general equation in the so called Lindblad form, able to describe the physical evolution of quantum states. [1][2][3][4][5] When these assumptions are not satisfied, e.g., for strong system-environment interaction and/ or long-living environmental correlations, we enter the intricate (and somewhat fuzzy) reign of non-Markovian dynamics. This consideration already illustrates how, despite its indubitable foundational nature, open quantum systems theory is still far from being completely understood and, in fact, it is peppered with unanswered questions of deep nature.
The increasing ability to coherently control an ever increasing number of individual quantum systems, together with the discovery of quantum coherence in complex biological systems, 6 has brought to light several scenarios in which the Markovian assumption fails. This has in turn given rise to a proliferation of results on the characterisation of memory effects and non-Markovian dynamics. 7-10 Interestingly, the cross-fertilisation of ideas from quantum information theory and open quantum systems has led to a new understanding of memory effects in terms of information backflow, [11][12][13][14][15][16] namely a partial return of quantum information previously lost from the open system due to the interaction with the environment.
Experimental results on quantum reservoir engineering, including the possibility to design desired forms of non-Markovian dynamics, [17][18][19][20][21][22][23][24][25] naturally lead to the question of whether or not memory effects are useful for quantum technologies, in the sense of constituting a resource for certain tasks. 15,[26][27][28][29][30][31] This question has not yet been satisfactorily answered. Even more remarkably, a complete understanding of non-Markovianity is still missing, as clearly illustrated in the insightful review of ref. 10 Given the considerations above, it is not surprising that in recent years a number of experiments have been proposed and realised to verify paradigmatic open quantum system models and test some of their predictions. Examples are numerous: the milestone experiment on the decoherence of a Schrödinger cat state with trapped ions is one of the first examples of engineered environment, 32 followed by the open system quantum simulator of ref. 33 The latter is also the first experimental realisation of an idea that shifted our perspective about environmental noise. Following a proposal by Vertsrate et al., 34,35 experimentalists proved that, by engineering certain types of Markovian master equation, one may actually create entangled states as stationary asymptotic states of the dynamics, therefore turning dissipation and decoherence from enemies to allies of quantum technologies. 33,36,37 In optical platforms, simulators of Markovian open quantum systems have been used to prove the existence of interesting phenomena, such as sudden death of entanglement 38,39 and sudden transition from quantum to classical decoherence. 40,41 In the same platform, experiments have shown how to engineer collisional models, 19 wherein the microscopic interaction between system and environment is obtained through a sequence of collisions between the open quantum system and one or more ancillae, the latter collectively describing the environment. 42 More recently, experiments in linear optics, [17][18][19]21,23,24 cavity QED, 43 NMR 22 and trapped ions 25 successfully demonstrated the engineering of non-Markovian open quantum systems and monitored the Markovian to non-Markovian crossover. Also complex quantum networks have been proposed as new systems for reservoir engineering of arbitrary spectral densities, 44 and a bosonic implementation with optical frequency combs has been presented. 45 Most of the experiments until now realised for simulating open quantum systems rely on the idea of analogue quantum simulator, that is a quantum system whose dynamics resemble those of another quantum system that we wish to study and understand. In contrast, a digital quantum simulator is a gate-based quantum computer which can be used to simulate any physical system, if suitably programmed. 46,47 Theoretical and experimental research on open systems digital quantum simulators is only now starting to flourish. [48][49][50][51][52][53] While, in principle, general algorithms for digital simulation of open quantum systems have been theoretically investigated, 48,50,53 their experimental implementation poses several challenges, since the physical quantum gates depend on the experimental platform and the circuit decomposition needs to be optimised in view of gates and measurement errors as well as qubit connectivity. Therefore, the existence of general algorithms for implementing theoretically universal open quantum system simulators does not guarantee the practical implementability in a realistic experiment on a given platform.
In this paper, we demonstrate that a careful circuit decomposition allows us to experimentally implement a vast number of fundamental open quantum systems models for one and two qubits. Not only are we able to generate different classes of open quantum dynamics, namely, unital (e.g., pure dephasing dynamical maps), non-unital (e.g., amplitude damping dynamical maps), phase covariant, and non-phase covariant (e.g., Pauli dynamical maps), but also we can explore the Markovian to non-Markovian crossover, including the recently discovered examples of essential 16 and eternal 54,55 non-Markovianity. We implement a recently proposed non-Markovianity witness 56 and we use our simulator to prove the non-monotonic behaviour of quantum channel capacity 15 and extractable work, 31 with implications to quantum communication and quantum thermodynamics.
Overall, our results clearly prove that even small quantum processors are versatile and robust testbeds for verifying a number of theoretical open quantum systems results and predictions, therefore paving the way to both new discoveries and a deeper understanding of one of the most fascinating and fundamental fields of quantum physics.

RESULTS
We consider an open quantum system represented by a density operator ρ S . The dynamics of the open system are described by a family Φ t of completely positive and trace preserving (CPTP) maps, known as the dynamical map: ρ S ðtÞ ¼ Φ t ρ S ð0Þ, with ρ S ð0Þ the initial state. The equation of motion for the state of the system is the master equation and, if the dynamical map is invertible, can be written in a time-local form _ ρ S ðtÞ ¼ L t ρ S ðtÞ; (1) where L t is the time-dependent generator of the dynamics with T the chronological ordering operator, and Φ 0 ¼ I. Under rather general conditions, 1 the generator can be written in the form In the equation above, the first term on the r.h.s. describes the unitary dynamics, with H S the system Hamiltonian, and the second term, the dissipative dynamics induced by the interaction with the environment, with γ k ðtÞ and V k the decay rates and jump operators, respectively. If the decay rates are positive and constant, i.e., γ k ðtÞ γ ! 0, the dynamical map is a semigroup, and we refer to the dynamics as semigroup Markovian. Extending this definition, it is nowadays common to say that the dynamics are Markovian whenever all the decay rates γ k ðtÞ are positive at all times. In this case, the dynamical map satisfies the property of CP-divisibility, namely Φ t ¼ Φ t;s Φ s , with Φ t;s a two-parameter family of CPTP maps. Non-Markovian dynamics occur, instead, whenever at least one of the decay rates becomes negative for a certain interval of time. In this case, the intermediate map Φ t;s is not CP anymore and the dynamics is non-CP-divisible. In the following subsections, we present experiments run on the IBM Q Experience processors simulating different types of open quantum systems dynamics, both Markovian and non-Markovian. We begin by considering an example of Markovian semigroup master equation and dynamical evolution.
Markovian reservoir engineering For decades, noise induced by the environment has been considered the archetype enemy of quantum technologies. This is because very often the interaction between a quantum system and its surroundings leads to the fast disappearance of quantum properties, notoriously coherences and entanglement, playing a key role in providing quantum advantage. This point of view drastically changed as soon as physicists demonstrated that appropriate manipulation of an artificial environment (quantum reservoir engineering) would allow one to steer the open system towards, e.g., a maximally entangled state, 33,36 hence turning upside down the perspective of the environment as an enemy.
Following the lines of ref., 33 in this subsection we experimentally simulate a semigroup Markovian master equation for a twoqubit open system having as asymptotic stationary state the Bell state ψ À j i ¼ 1 ffiffi 2 p ð 01 j i À 10 j iÞ, where we indicate with 0 j i and 1 j i the computational basis of each qubit and we use the notation 01 j i ¼ 0 j i 1 1 j i 2 . This allows us to prepare a maximally entangled state as a result of the dissipative open system dynamics.
Each of the four Bell states is uniquely determined as an eigenstate with eigenvalues ± 1 with respect to σ The dissipative dynamics that pumps two qubits from an arbitrary initial state into the Bell state ψ À j i is realised by the composition of two channels that pump from the x . By changing the parameter 0 p 1 we simulate different types of open quantum system dynamics. For p ( 1, the repeated application of, e.g., Φ zz generates a master equation of Lindblad form with jump operator operator for any initial state. In Fig. 1, we show the action of the dissipative pumping maps Φ xx , Φ zz and their composition Φ xx Φ zz as a function of p, for a maximally mixed initial state. We compare the theoretical prediction (dashed lines) with the experimental data. Our results show a very good agreement between theory and experiment for the implementation of both the two families of maps and their composition. In the latter case, the results are more sensitive to errors because of the larger depth of the circuit implementing the composition of maps. Details on the circuit implementations and on their optimisation with respect to the specific qubit connectivity are given in Methods.
Collisional model and essential non-Markovianity Our second example of an open quantum system simulator deals with a class of models known as collisional models. These describe the interaction between a quantum system and its environment in terms of consecutive collisions between the system, in our case a qubit, and a sequence of environmental qubits (ancillae) in a given state. The system qubit and the nth environmental qubit interact pairwise during a time period τ. One assumes, as usual, a factorised initial state of system and environmental qubits. After the system has interacted with a given ancilla, one can trace out the ancillary degrees of freedom as it no longer affects the system's dynamics.
For a sufficiently large number of collisions, when the ancillae are in a thermal state, one can prove that the equation of motion describing the system's dynamics is of Gorini-Kossakowski-Sudarshan-Lindblad form. 4 Interestingly, contrarily to the microscopic derivation of the Markovian master equation, this approach does not rely on Born-Markov approximation, but it automatically leads to a CPTP dynamical semigroup. Also non-Markovian master equations can be introduced in the collisional model picture, for example by allowing for the ancillae to interact with each other in a well defined manner. [57][58][59] We implement experimentally an exemplary model of dynamics displaying the property of essential non-Markovianity, 16 following the lines of ref. 60 As research in the field of non-Markovian open quantum systems progressed, more refined definitions of memory effects started to appear in the literature. In ref. 16 a hierarchical classification of non-Markovianity was introduced, generalising the notion of CP-divisibility. In particular, the dynamical map is said to be P-divisible if the two-parameter family Φ t;s is positive. This is, of course, a weaker condition than CP-divisibility, since there exist maps that are P-divisible but not CP-divisible, while all CP-divisible maps are also P-divisible. We say that the dynamics is essentially non-Markovian if it is non-P-divisible, while it is weakly non-Markovian if it is non-CP-divisible but P-divisible.
We consider n ancillae initially prepared in the classically correlated state ρ corr ¼ 1 2 ð 0 j i n 0 h j n þ 1 j i n 1 h j n Þ. The collision between the system qubit and the kth environmental ancilla is described by the unitary operator U k ¼ e igτσz 0 j i k 0 h jþ e Àigτσz 1 j i k 1 h j, where g is the coupling strength. The dynamical map after n ¼ t=τ collisions is given by We compare this dynamics to the case in which the environmental ancillae are prepared in the uncorrelated state þ j i n . In this case, the dynamics is given by In the weak-coupling case (when gτ < π=4) the map in Eq. (6) gives rise to Markovian dynamics. In contrast, the map in Eq. (5) alternates intervals, with periodicity π=2, in which it is P-divisible (0 < ngτ < π=4) with intervals in which it is non-P-divisible (π=4 < ngτ < π=2), i.e., it is essentially non-Markovian. In Fig. 2, we compare, for both maps, the exact dynamics of the qubit coherence with the experimental results, for the initial state The experimental data are in good agreement with the theoretical predictions, and the oscillatory behaviour of the coherences in the case of correlated ancillae is a signature of the essential non-Markovianity of the dynamics. We note that the (much smaller) oscillations in the separable case are due to the repeated, imperfect application of the CNOT gates that Fig. 1 Simulation of the Markovian reservoir engineering. Every plot compares the measured overlap between the state of the system and the four Bell states (dots) with the corresponding analytical prediction for the ZZ pump (a), the XX pump (b) and their consecutive application (c). The pumps are applied to an initial maximally mixed state. The first two pumps clearly reduce the populations of the þ1 eigenspace of the corresponding stabiliser operator, and their composition results on a probability for the ψ À j i state around 0.94. Throughout the paper, the error bars indicating statistical errors are not shown, as they would be indistinguishable from the markers given the large number of experimental runs (see Methods).
are required to swap the ancillae that have interacted with the qubit with fresh ones: there is a small probability that the swap fails and the system interacts again with the same ancilla. This can introduce memory effects, as shown, e.g., in refs. 58,59 Markovian and non-Markovian dissipative dynamics The dynamical maps of Eqs. (5) and (6) are purely dephasing, since they affect only the coherences of the qubit. In this section we simulate an exactly solvable dissipative open quantum system known as the Jaynes-Cummings or generalised amplitude damping model. The microscopic derivation of the master equation can be found in textbooks (see, e.g., ref. 1 ). The master equation is given by where σ ± are the raising and lowering operators and the timedependent decay rate γðtÞ has the following analytical form with The equation above is obtained assuming that the environment is a bosonic zero-temperature reservoir with Lorentzian spectral density with ω 0 the qubit frequency, γ 0 an overall coupling strength, and λ the half height width of the Lorentzian profile. The timedependent coefficient c 1 ðtÞ defined in Eq. (9) crucially depends on the ratio R ¼ γ 0 =λ between the coupling strength and the width of the spectrum. Also, this coefficient fully determines the qubit dynamics as one can see from the following expression From Eq. (9) a straightforward calculation shows that the timedependent decay rate γðtÞ, defined in Eq. (8), takes negative values for certain time intervals whenever 2R ! 1. In this case the dynamical map is not CP-divisible, and therefore non-Markovian. In Fig. 3a, b we show, as example of Markovian and non-Markovian dynamics, the evolution of the excited state population of an initially excited state for two values of the parameter R, corresponding to Markovian (R ¼ 0:2) and non-Markovian (R ¼ 100) dynamics, respectively. The figures show the monotonic decay of the excited state population in the former case, while in the latter case the population oscillates since the qubit exchanges information and energy with the central mode of the Lorentzian peak, resonant with the qubit's Bohr frequency. As done for the other open quantum systems simulators, we compare the experimental data with theoretical predictions, finding a good agreement.
We also implement an experimentally friendly witness of non-Markovianity which was recently introduced in ref. 56 and stems from the spectral properties of the dynamical map. The witness is based on the behaviour of an initial maximally entangled state of qubit and auxiliary ancilla, and requires only the measurement of the expectation values of local observables such as σ i σ i , (i ¼ x; y; z). Non-monotonic behaviour of the witness as a function of time signals non-Markovianity. Other witnesses of non-Markovianity (such as the BLP measure introduced in ref. 11 ) could be implemented, and may be more efficient at detecting memory effects, but they generally require a larger number of measurements or optimisation over initial states. Fig. 2 Simulation of a collisional model. The red triangles show the collisions with the ancillae prepared in a classically correlated state (red triangles), while the blue dots show the case of a separable state. In both cases g ¼ 1, τ ¼ π=6 (weak-coupling regime). The red and blue dashed curves show, respectively, the theoretical predictions of Eqs. (5) and (6) for time t ¼ ngτ. Up to 7 collisions were simulated: the depth of the circuits grows with the number of collisions, causing increasing decoherence, but the oscillations due to the essential non-Markovianity are clearly visible in the case of correlated ancillae. In Fig. 3, we plot the dynamics of the entanglement witness for the amplitude damping channel. In Fig 3b, the witness clearly shows oscillatory behaviour, and therefore properly signals the presence of memory effects. The circuits implementing both the amplitude damping model and the non-Markovianity witness are presented and discussed in Methods.

Depolarising and Pauli channels
The depolarising channel is one of the most common models of qubit decoherence due to its nice symmetry properties. We can describe it by stating that, with probability 1 À p the qubit remains intact, while it becomes maximally mixed with probability p. This mixing can be modelled through a bit flip error, described by the action of σ x , a phase flip error, described by the action of σ z , or both, described by the action of σ y . The dynamical map of a Markovian open quantum system subjected to depolarising noise can be written as where i ¼ x; y; z and pðtÞ ¼ 1 À e Àγt , with γ the Markovian decay rate. In Fig. 4, we plot the qubit density matrix elements for various values of p, comparing the experimental data with the theoretical prediction. For exemplary purposes, we choose a specific initial state possessing non-zero coherences, but we have verified, by repeating the experiment for different initial states, that the agreement observed between experiment and theory is independent from the initial state chosen. As for the other simulated models, we postpone the discussion on the circuit implementation, the readout, and the error mitigation strategy to the Methods section.
Let us now introduce the most general single-qubit open quantum system model, namely the time-dependent Pauli channel. The master equation in this case takes the form We note that for γ i ðtÞ ¼ γ, we recover the Markovian depolarising channel. Generally, the dynamics described by the master equation above is not phase-covariant, 61 except for the case in which γ x ðtÞ ¼ γ y ðtÞ. Moreover, since the decay rates may take negative values, conditions for complete positivity must be imposed, and they are given in terms of a set of inequalities involving all the three decay rates, as one can see, e.g., from ref. 16 In the Discussion section we present the simulation of a specific form of time-dependent Pauli channel proposed in ref. 54 and used as an example of eternal non-Markovianity, i.e., an open quantum system dynamics for which the dynamical map is non-CP-divisible for all times t. More precisely, we use this experimental simulation to demonstrate a phenomenon predicted in ref., 31 namely the presence of oscillations in the extractable work. This shows an application of open quantum system simulation on the IBM Q Experience processors to fields other than quantum information theory, specifically quantum thermodynamics for the example here considered.

DISCUSSION
The results presented in the previous section highlight how fewqubit NISQ devices publicly available on the cloud already provide sufficient robustness and reliability to implement experimentally a number of open quantum systems dynamics, both Markovian and non-Markovian, both unital (dephasing collisional model) and non-unital (amplitude damping, depolarising and Pauli channels).
We have simulated all the paradigmatic open quantum systems models typically used to demonstrate physical phenomena induced by the presence of the environment, its consequences for quantum technological applications, but also possible benefits in the spirit of reservoir engineering. As an example of a fundamental study on open quantum systems, we have measured the dynamics of a non-Markovianity witness and showed that it correctly signals the presence of memory effects for time-local amplitude damping dynamics.
In this section, we build on these results to substantiate our claim of the IBM Q Experience being a versatile testbed for the experimental verification of physical effects due to the open character of the dynamics. Specifically, we focus on two applications: non-Markovian quantum channel capacity, of potential use in quantum communication, and extractable work dynamics, relevant in quantum thermodynamics.
Let us first consider a typical setting in quantum information processing and communication: Alice and Bob are at the opposite ends of a quantum channel, the former sending information (classical or quantum) and the latter receiving it. The maximum amount of information that can be reliably transmitted along a noisy quantum channel is known as the channel capacity. Specifically, we consider the quantum capacity Q defined as the limit to the rate at which quantum information can be reliably sent down a quantum channel.
In the following, we focus on the time-local amplitude damping model introduced in the previous section. In this case, the quantum channel capacity can be calculated exactly and takes the form 15 QðΦ t Þ ¼ max for jc 1 ðtÞj 2 > 1=2 and QðΦ t Þ ¼ 0 otherwise. In the equation above, H 2 is the binary Shannon entropy and c 1 ðtÞ is given by Eq. (9). One may be led into believing that the quantum channel capacity decreases, and in some case disappears, for increasing lengths of the noisy channel, due to the cumulative destructive effect of decoherence. This, however, only holds for Markovian noise, as theoretically demonstrated in ref. 15 Indeed, due to memory effects, the quantum channel capacity may take again non-zero values after disappearing for a certain finite length of the channel. In Fig. 5a, we experimentally demonstrate this effect, comparing the measured behaviour of QðΦ t Þ with the theoretical prediction for an amplitude damping model with R ¼ 100, R ¼ 200, and R ¼ 400.
We now consider an example of interest in the field of quantum thermodynamics, specifically related to the presence of memory effects in the open system dynamics. In ref., 31 the connection between information backflow and thermodynamic quantities was established by means of a non-Markovianity indicator based on the quantum mutual information. More specifically, it was proven that there exists a link between memory effects, as To understand this connection, let us recall that in order to link non-Markovianity with the evolution of thermodynamical quantities one needs to describe not only the system S, whose information content we are interested in, but also an observer or memory M. Interestingly, quantum mechanical correlations between system and observer may lead to the exciting possibility of extracting work while erasing information on the system. 62 Following ref., 31 we consider the case of a qubit S subjected to a Pauli dynamical map, see master Eq. (13), and an isolated qubit acting as memory, M. If the system and memory are prepared in a maximally entangled state, work can be extracted during the information erasure by using the initial entanglement, as predicted in ref. 62 In this framework, the extractable work takes the form 31 with n the number of qubits in S, HðΦ t ρ S Þ the von Neumann entropy of the evolved state of the system, IðΦ t ρ S : ρ M Þ the quantum mutual information quantifying the amount of total correlations between system S and memory M, k the Boltzmann constant and T the temperature of the reservoir used for the erasure. For the open system here considered, revivals of extractable work are due to the interplay between memory effects, witnessed by the quantum mutual information, and entropy dynamics, as dictated by Eq. (15).
We now further specify our analysis to the eternal non-Markovianity model, for which γ x ðtÞ ¼ γ y ðtÞ ¼ λ=2 and γ z ðtÞ ¼ Àω tanhðωtÞ=2. We compare the dynamics with the case in which γ x ðtÞ ¼ γ y ðtÞ ¼ λ=2 and γ z ðtÞ ¼ ω tanðωtÞ=2. The first dynamics is eternally non-Markovian, since γ z ðtÞ < 0 at all times, while the second one is non-CP-divisible but not eternally non-Markovian, since γ z ðtÞ takes both positive and negative values. This comparison allows us to distinguish between quantum and classical memory effects. 31 In Fig. 5b, we show the experimental results and the theoretical prediction for the dynamics of the extractable work for the two cases here considered. Despite the presence of experimental imperfections, clear oscillations of the extractable work are visible in one case, but completely absent in the other one, therefore showing for the first time not only the impact of memory effects on a thermodynamic quantity such as W ex , but also the subtle origin of these oscillations.
To conclude, it is worth mentioning that experiments on open quantum systems performed in dedicated experimental laboratories in different platforms such as, e.g., trapped ions, linear optical systems, superconducting qubits, NMR, NV centres in diamonds and cavity QED, generally reach higher precision and fidelity than those reported in our paper. However, this often occurs at the expense of increased specificity in terms of the models that can be simulated. Moreover, these experiments are accessible only to the researchers working in the given laboratory.
Open source small quantum processors already are, however, a reality and they have opened the door of experimental physics to quantum researchers who do not have either specific experimental training or sufficient economic resources to set up a laboratory. We therefore envisage a rapid change in the way in which research is conducted in a field which has been until now dominantly theoretical, namely the study of open quantum system dynamics. Perhaps, in a not too far future, this will blur the line separating theoretical and experimental research, and an increasing number of experiments will be remotely programmable by using a simple interface and language, in the true spirit of a quantum simulator.

METHODS
The results presented in this paper have been obtained on the IBM Q Experience processors freely available online: two 5-qubit machines, ibmqx2 and ibmq_vigo, and a 14-qubit machine, ibmq_16_melbourne. The circuits have been written using IBM's Qiskit, 63 a Python-based programming language for the IBM Q Experience. Each circuit has been run for 8192 shots (the maximum allowed by the devices) to gather statistics from the measurements. When needed, one-and two-qubit full state tomography has been obtained by using the tools provided in Qiskit, performing a maximum-likelihood reconstruction of the density matrix. 64 In the following, we give detailed explanations of the various circuits implemented to produce the results presented in the paper. We first make a few general considerations on the sources of error and on how to devise circuits with high fidelity on the IBM Q Experience devices.
Given the large number of experimental runs allowed by the devices, the errors due to statistical fluctuations, of order Oð1= ffiffiffi ffi N p Þ % 0:011, are too small for errorbars to be distinguishable from markers in Figs 1-5, and thus they are not shown. The discrepancy between the experimental points and the theoretical prediction, on the other hand, comes from systematic errors induced by various factors in the hardware implementation. First, the qubits undergo relaxation and decoherence due to external noise. This error becomes more relevant when increasing the depth of the circuit. Second, the gates have errors due to cross-talk and unwanted interactions between qubits when addressing them with pulses. In particular, the error rate of CNOT gates is about 10 times larger than single-qubit unitaries. Finally, there is a readout error affecting the quantum measurement, although this can be mitigated to some extent, as we discuss later on. The error rates and noise parameters on IBM Q devices are characterised on a daily basis using randomised benchmarking techniques. Since there are various, interdependent sources of noise, modelling and predicting the deviations in the experimental data is a non-trivial matter, beyond the scope of this paper. In this work, all quantum channels have been  14), where jc 1 ðtÞj 2 has been approximated by the ratio between the population of the excited state at time t and at t ¼ 0 in order to account for the deviations in the initial state preparation. b Rescaled extractable work W ex ðtÞ=kTln2 as a function of time for the eternally non-Markovian Pauli channel, with γ x ðtÞ ¼ γ y ðtÞ ¼ λ=2 and γ z ðtÞ ¼ Àω tanh ωt=2, with λ ¼ 1 and ω ¼ 1 2 (green circles) and for the non-CP-divisible Pauli channel γ x ðtÞ ¼ γ y ðtÞ ¼ λ=2 and γ z ðtÞ ¼ ω tan ωt=2, with λ ¼ 0:1 and ω ¼ 2. The extractable work has been evaluated using Eq. (15) from a full two-qubit tomography of system S and memory ancilla M.
decomposed into circuits bearing in mind some general guidelines targeted at minimising the aforementioned inaccuracies.
Since the IBM Q Experience devices are universal quantum computers, they enable the implementation of any unitary transformation of their constituent qubits; once a quantum circuit is provided for its execution, it is compiled into an equivalent circuit involving only the machine's basis gates, that is, those realisable experimentally. In the case of the IBM Q Experience devices, these include any single-qubit rotation, whereas the only multi-qubit gates are CNOTs. However, if the circuit requires a multiqubit gate among qubits that are not physically connected, the corresponding gate will be replaced with a longer circuit in which the states are swapped to neighbouring qubits in the compiled circuit. Since every swap gate includes a minimum of three CNOT gates, and these introduce considerable noise to the execution, it is crucial to assign the relevant qubits involved in the simulation (e.g., system, environment and ancillae) to the machine's qubits so that the number of CNOT gates between disconnected qubits is minimised. Furthermore, the devices are calibrated daily and the errors of the basis gates are reported. This information can also be taken into account in the qubit assignment, as using the CNOT gates with smaller errors is preferable.
In addition to the gate errors, the qubits' readout errors characterising the discrepancies between the qubit state and the measurement outcome probabilities are also provided. In the IBM Q Experience devices, there are considerable differences in the readout errors of the different qubits, so this information should also be taken into account in the qubit-assignment process: if possible, it is preferable to assign the system (and any auxiliary ancillae whose measurement is required) to low readout-error qubits, while qubits with large readout errors can still be used to simulate the environment or other ancillae that need not be measured. In any case, Qiskit provides post-processing error-mitigation tools that generally improve the experimental results under the package ignis. To do so, we first prepare all possible basis states 0 Á Á Á 0 j i ; 0 Á Á Á 1 j i ; ; 1 Á Á Á 1 j iand measure their corresponding outcome probabilities (notice that it is only necessary to include the qubits whose measurement outcomes need to be corrected). Once these are known, they can be used to correct any other experimental result by finding, via likelihood maximisation, the experimental outcome that is most congruent with the observed measurement deviations. All the data shown in the paper have been mitigated as described above, with the only exception of the collisional model with the ancillae prepared in the separable state; the reason why we have not mitigated those results is that, due to the qubit swapping, the measurements are performed on different physical qubits.

Reservoir engineering
In ref., 33 the authors provide the circuits for the implementation of the Bell-state pumping. However, these are composed of gates that are natural to the trapped-ions platform used in that work, so their direct implementation on the IBM Q Experience devices would result in far too long circuits. Therefore, we propose a different set of circuits that follow the same basic working principles, but have been designed specifically keeping in mind the characteristics of the IBM Q Experience platform.
The pumping circuits proposed in ref. 33 are composed of four parts. First, the relevant information regarding the state of the system (that is, whether the system is in the þ1 or the À1 eigenspaces of the stabiliser operators) is mapped into an ancilla. Second, the state of the system is modified depending on the state of the ancilla. Third, the mapping circuit is reversed. At this stage, the system has been pumped, but if the ancilla is to be used again for a new pumping cycle, it needs to be reset, which is the fourth step. We follow these same lines, designing circuits that perform these same steps while minimising the number of gates involved. Before we explain the resulting circuits, let us mention that, since the IBM Q Experience devices are not equipped with the reset operation, we must use a different ancilla for every pump.
The way we map the eigenspace information into an ancilla is by first applying a CNOT gate between the system qubits. Suppose that qubits s 1 and s 2 are initially in some Bell state, for instance, ϕ ± j i ¼ ð 00 j i± 11 j iÞ= ffiffi ffi 2 p . A CNOT gate controlled by s 1 transforms the state into ± j i 0 j i. Instead, ψ ± j i would be transformed into ± j i 1 j i. Hence, we see that the information regarding the σ  Fig. 6a. To map the eigenspace information into the environment ancilla a ZZ , we apply a CNOT controlled by the relevant qubit, s 2 . After these two gates (and considering that the initial state of the ancilla is 1 j i), a ZZ will be in state 1 j i if the initial state of the system is ϕ ± j i and 0 j i if it is ψ ± j i. Therefore, the conditional rotation gate only acts in the former case, while it does not modify the state in the latter. The angle of the controlled rotation, in turn, controls the efficiency of the pump p via the relation θ ¼ 2 arcsin ffiffi ffi p p .
Finally, the last two CNOT gates simply revert the mapping part of the circuit. The working principle of the σ ð1Þ x σ ð2Þ x pump (Fig. 6b) is essentially the same. However, we need to add an extra Hadamard gate to transform the state of s 1 before mapping the information to the ancilla a XX . As for the composite pump, we can simply concatenate the two circuits. Notice that in the direct concatenation there would be two consecutive CNOTs between the system qubits, which can be removed.
Regarding the measurement process, the IBM Q Experience platform only enables measurement in the computational basis. Hence, to assess the probabilities for each of the Bell states, we need to change basis by applying again a combination of CNOT and Hadamard to the system qubits, so ϕ þ j i is mapped into 00 j i, etc. Again, notice that this would result in repeated consecutive gates, the effect of which amounts to identity, so they can be removed from the circuits. Finally, in Fig. 1, we show the results starting from the maximally mixed state. To do so, we simulated the circuits preparing the system in four different initial conditions, 00 j i, 01 j i, 10 j i and 11 j i.

Collisional model
The circuit used to implement the collisional model described in the section "Collisional model and essential non-Markovianity" is depicted in Fig. 7. Initially, the system and ancillae are prepared in the desired initial state, then each collision is applied, and finally the qubit is measured in the σ x basis. The collision unitary U ¼ e igτσz 0 j i 0 h j þ e Àigτσz 1 j i 1 h j can be implemented by means of a rotation R z ðgτÞ around the Z-axis between two CNOT gates.
When the ancillae are prepared in the state ρ corr ¼ ð 0 j i n 0 h j n þ 1 j i n 1 h j n Þ=2, we can simplify the circuit by noticing that each ancilla is in the maximally mixed state, and the action of U does not affect its state. We thus only need three ancillae prepared in the GHZ state ψ GHZ j i¼ ð 000 j iþ 111 j iÞ= ffiffi ffi 2 p : by tracing out the third one we are left with the two-qubit state ρ ð2Þ E ¼ ð 00 j i 00 h j þ 11 j i 11 h jÞ=2, and then we can keep colliding with either of the two ancillae.
The circuit with the ancillae in the separable state was implemented and run on the device ibmq_16_melbourne (14 qubits), in order to have a larger number of ancillae for the collisions. The system qubit was chosen to have the highest connectivity (three qubits) and smallest readout error. The connectivity layout of the computer does not allow for direct collisions with more than three ancillae, and swaps between qubits are required, increasing the errors in the simulation and possibly introducing memory effects (as discussed in the Results section). Fig. 6 Circuits implementing the reservoir engineering protocol. a ZZ pump, b XX pump and c ZZ and XX pump. The circuits were run on ibmqx2. Qubits s 1 and s 2 (corresponding to q 2 and q 1 on the device) are the system qubits, while qubits a XX and a ZZ (q 4 and q 0 ) are the environment ancillae for the two maps. State preparation and measurement are not shown in the circuits above.

Amplitude damping channel
The circuit in Fig. 8 implements the amplitude damping channel with the non-Markovianity witness. For an arbitrary pure state of the system ψ j i s ¼ α 0 j i s þ β 1 j i s , and setting the state of the environment to vacuum 0 j i e , the two gates between the system and environment qubits transform the joint state into α 0 j i s 0 j i e þ β cos θ 1 j i s 0 j i e þ β sin θ 0 j i s 1 j i e . Therefore, identifying the states 0 j i s and 1 j i s with the ground and excited states respectively, and by choosing θ ¼ arccos c 1 ðtÞ, the reduced state of the system becomes Eq. (11).
The non-Markovianity witness for a channel Φ t proposed in ref. 56 is based on the dynamical behaviour of the entanglement between the system and an ancilla initially prepared in a maximally entangled state. Namely, this quantity is defined as f Φ ¼ hϕ þ jI Φ t ½jϕ þ ihϕ þ jjϕ þ i. Hence, implementing this quantity requires preparing the state ϕ þ j i between the system and an ancilla, which can be achieved using a Hadamard and a CNOT gate, applying the dynamical map to the system and measuring the probability of finding the joint system and ancilla state in ϕ þ j i. As suggested in ref., 56 we need not use an extra CNOT gate to project on the Bell basis; instead, we can take advantage of the fact that ϕ þ j i ϕ þ h j ¼ ðI I þ σ x σ x À σ y σ y þ σ z σ z Þ=4 and measure the corresponding local observables. Figure 8 shows the circuit corresponding to the measurement of σ x σ x .

Depolarising channel
The depolarising channel defined in Eq. (12) can be implemented, for any value of p pðtÞ 2 ½0; 1, with the circuit shown in Fig. 9. Three ancillary qubits are prepared in a state ψ θ j i ¼ cos θ=2 0 j i þ sin θ=2 1 j i, and are used as controls for, respectively, a controlled-X (CNOT), a controlled-Y and a controlled-Z rotation. This way, each gate will be applied with a probability sin 2 θ=2.
The rotation angle θ must be chosen so that each of the gates is applied with probability p. Notice that applying X and then Y, but not applying Z is equivalent (up to global phases) to just applying Z, and so on. The resulting equation that binds θ to p is thus with solution θðpÞ ¼ 1 2 arccosð1 À 2pÞ.

Pauli channel
At a specific time instant t, the Pauli channel described by the master Eq. (13) can be written as with 0 p i 1 and P i p i ¼ 1. The depolarising channel is a special case of the Pauli channel where p 1 ¼ p 2 ¼ p 3 ¼ p=4.
One could think to implement the channel by generalising the circuit of Fig. 9, specifically, by allowing for the three ancillae to be prepared in different states. One can see, however, that this cannot be done in the most general case. For example, the eternally non-Markovian channel presented in the Discussion section is not implementable in this way.
It is possible to implement the general Pauli channel with just two qubits, if we allow them to be in an entangled state. The first qubit acts as the control for a controlled-X (CNOT) gate, and the second one for a controlled-Y. Notice that, as remarked in the previous section, applying both X and Y is effectively equivalent to applying a Z gate.
The state ψ j i of the ancillae needed for the Pauli channel can be implemented by the circuit in Fig. 10, parametrised by the three angles θ 1 ; θ 2 ; θ 3 ψ j i ¼ ðc 1 c 2 c 3 þ s 1 s 2 s 3 Þ 00 j iþ ðc 1 c 2 s 3 À s 1 s 2 c 3 Þ 01 j i þ ðc 1 s 2 c 3 À s 1 c 2 s 3 Þ 10 j iþ ðs 1 c 2 c 3 þ c 1 s 2 s 3 Þ 11 j i where c i cos θ i and s i sin θ i . The angles θ i can be found by solving the following system of equations:  Fig. 7 Circuits implementing the collisional model. The top circuit shows the case of the collision with three ancillae in the separable þ þ þ j istate, prepared by means of Hadamard gates. Each collision consists of a rotation around Z between two CNOTs. A final measurement in the X basis is performed in order to measure the coherence of the system qubit. The bottom circuit shows how, in the case of ancillae prepared in the correlated state, we can effectively simulate the map with three ancillae prepared in the GHZ state ð 000 j iþ 111 j iÞ= ffiffi ffi 2 p , by colliding alternately with two of them. Notice that, in principle, we could have always a collision with the same ancilla. This, however, would cause the compiler to remove consecutive CNOTs and join the rotations into a single gate, making the circuit trivial. The top circuit was run on ibmq_16_melbourne, while the bottom circit was run on ibmqx2. Fig. 8 Circuit for the amplitude damping channel with non-Markovianity witness. The circuit was run on ibmq_vigo. Qubits q 1 , q 2 and q 3 were used as system, environment and witness ancilla, respectively. Setting θ ¼ arccos c 1 ðtÞ, the channel acting on the system qubit simulates the amplitude damping dynamics. The two Hadamard gates rotate the state of the qubits in order to measure the observable σ x σ x . To measure σ y σ y , the combination S y H can be used instead, whereas no gates are required for σ z σ z . Fig. 9 Circuit implementation of the depolarising channel. The circuit was run on ibmqx2. Qubit q 2 is used as the system qubit, as it is the only one that can be targeted with three CNOTs. The ancillae a 1 ; a 2 and a 3 are initially in the state 0 j i and rotated around the y direction by an angle θðpÞ ¼ 1 2 arccosð1 À 2pÞ. They then act as controls for controlled-X, -Y and -Z gates, respectively. Fig. 10 Circuit implementation of a general Pauli channel on ibmqx2. Qubit q 2 is used as the system qubit because it can be targeted with two CNOTs and has the smallest readout error. The ancillae are initially in the state 0 j i and rotated around y with the angles θ 1 ; θ 2 and θ 3 found by solving the system in Eq. (19), and they act as controls for a controlled X and a controlled Y.