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 \({\mathcal{O}}(1/\sqrt{N})\approx 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 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 multi-qubit 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 \(\left|0\cdots 0\right\rangle ,\left|0\cdots 1\right\rangle ,\ldots ,\left|1\cdots 1\right\rangle\) and 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, \(\left|{\phi }^{\pm }\right\rangle =(\left|00\right\rangle \pm \left|11\right\rangle )/\sqrt{2}\). A CNOT gate controlled by \({s}_{1}\) transforms the state into \(\left|\pm \right\rangle \left|0\right\rangle\). Instead, \(\left|{\psi }^{\pm }\right\rangle\) would be transformed into \(\left|\pm \right\rangle \left|1\right\rangle\). Hence, we see that the information regarding the \({\sigma }_{x}^{(1)}\otimes {\sigma }_{x}^{(2)}\) eigenspace (namely, the sign) is contained in the state of \({s}_{1}\) after the transformation, whereas the one corresponding to the \({\sigma }_{z}^{(1)}\otimes {\sigma }_{z}^{(2)}\) is in qubit \({s}_{2}\). Now, let us consider the circuit implementing the \({\sigma }_{z}^{(1)}\otimes {\sigma }_{z}^{(2)}\) pump in Fig. 6a. To map the eigenspace information into the environment ancilla \({a}_{{\rm{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 \(\left|1\right\rangle\)), \({a}_{{\rm{ZZ}}}\) will be in state \(\left|1\right\rangle\) if the initial state of the system is \(\left|{\phi }^{\pm }\right\rangle\) and \(\left|0\right\rangle\) if it is \(\left|{\psi }^{\pm }\right\rangle\). 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 \(\theta =2\arcsin \ \sqrt{p}\). Finally, the last two CNOT gates simply revert the mapping part of the circuit. The working principle of the \({\sigma }_{x}^{(1)}\otimes {\sigma }_{x}^{(2)}\) 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}_{{\rm{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.

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}_{{\rm{XX}}}\) and \({a}_{{\rm{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. Full size image

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 \(\left|{\phi }^{+}\right\rangle\) is mapped into \(\left|00\right\rangle\), 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, \(\left|00\right\rangle\), \(\left|01\right\rangle\), \(\left|10\right\rangle\) and \(\left|11\right\rangle\).

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 \({\sigma }_{x}\) basis. The collision unitary \(U={e}^{ig\tau {\sigma }_{z}}\otimes \left|0\right\rangle \left\langle 0\right|+{e}^{-ig\tau {\sigma }_{z}}\otimes \left|1\right\rangle \left\langle 1\right|\) can be implemented by means of a rotation \({R}_{z}(g\tau )\) around the \(Z\)-axis between two CNOT gates.

Fig. 7: Circuits implementing the collisional model. The top circuit shows the case of the collision with three ancillae in the separable \(\left|+++\right\rangle\) state, 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 \((\left|000\right\rangle +\left|111\right\rangle )/\sqrt{2}\), 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. Full size image

When the ancillae are prepared in the state \({\rho }_{{\rm{corr}}}=({\left|0\right\rangle }^{\otimes n}{\left\langle 0\right|}^{\otimes n}+{\left|1\right\rangle }^{\otimes n}{\left\langle 1\right|}^{\otimes 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 \(\left|{\psi }_{{\rm{GHZ}}}\right\rangle =(\left|000\right\rangle +\left|111\right\rangle )/\sqrt{2}\): by tracing out the third one we are left with the two-qubit state \({\rho }_{E}^{(2)}=(\left|00\right\rangle \left\langle 00\right|+\left|11\right\rangle \left\langle 11\right|)/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).

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 \({\left|\psi \right\rangle }_{s}=\alpha {\left|0\right\rangle }_{s}+\beta {\left|1\right\rangle }_{s}\), and setting the state of the environment to vacuum \({\left|0\right\rangle }_{e}\), the two gates between the system and environment qubits transform the joint state into \(\alpha {\left|0\right\rangle }_{s}{\left|0\right\rangle }_{e}+\beta \cos \theta {\left|1\right\rangle }_{s}{\left|0\right\rangle }_{e}+\beta \sin \theta {\left|0\right\rangle }_{s}{\left|1\right\rangle }_{e}\). Therefore, identifying the states \({\left|0\right\rangle }_{s}\) and \({\left|1\right\rangle }_{s}\) with the ground and excited states respectively, and by choosing \(\theta =\arccos {c}_{1}(t)\), the reduced state of the system becomes Eq. (11).

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 \(\theta =\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 \({\sigma }_{x}\otimes {\sigma }_{x}\). To measure \({\sigma }_{y}\otimes {\sigma }_{y}\), the combination \({S}^{\dagger }H\) can be used instead, whereas no gates are required for \({\sigma }_{z}\otimes {\sigma }_{z}\). Full size image

The non-Markovianity witness for a channel \({\Phi }_{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}_{\Phi }=\langle {\phi }^{+}| {\mathbb{I}}\otimes {\Phi }_{t}[| {\phi }^{+}\rangle \langle {\phi }^{+}| ]| {\phi }^{+}\rangle\). Hence, implementing this quantity requires preparing the state \(\left|{\phi }^{+}\right\rangle\) 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 \(\left|{\phi }^{+}\right\rangle\). 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 \(\left|{\phi }^{+}\right\rangle \left\langle {\phi }^{+}\right|=({\mathbb{I}}\otimes {\mathbb{I}}+{\sigma }_{x}\otimes {\sigma }_{x}-{\sigma }_{y}\otimes {\sigma }_{y}+{\sigma }_{z}\otimes {\sigma }_{z})/4\) and measure the corresponding local observables. Figure 8 shows the circuit corresponding to the measurement of \({\sigma }_{x}\otimes {\sigma }_{x}\).

Depolarising channel

The depolarising channel defined in Eq. (12) can be implemented, for any value of \(p\equiv p(t)\in [0,1]\), with the circuit shown in Fig. 9. Three ancillary qubits are prepared in a state \(\left|{\psi }_{\theta }\right\rangle =\cos \theta /2\left|0\right\rangle +\sin \theta /2\left|1\right\rangle\), 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}\theta /2\).

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 \(\left|0\right\rangle\) and rotated around the \(y\) direction by an angle \(\theta (p)=\frac{1}{2}\arccos (1-2p)\). They then act as controls for controlled-\(X\), -\(Y\) and -\(Z\) gates, respectively. Full size image

The rotation angle \(\theta\) 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 \(\theta\) to \(p\) is thus

$${\sin }^{2}\frac{\theta }{2}{\cos }^{4}\frac{\theta }{2}+{\sin }^{4}\frac{\theta }{2}{\cos }^{2}\frac{\theta }{2}=\frac{p}{4},$$ (16)

with solution \(\theta (p)=\frac{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

$${\mathcal{E}}(\rho )=\sum _{i=0}^{3}{p}_{i}{\sigma }_{i}\rho {\sigma }_{i},$$ (17)

with \(0\le {p}_{i}\le 1\) and \({\sum }_{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 \(\left|\psi \right\rangle\) of the ancillae needed for the Pauli channel can be implemented by the circuit in Fig. 10, parametrised by the three angles \({\theta }_{1},{\theta }_{2},{\theta }_{3}\)

$$\begin{array}{l}\left|\psi \right\rangle =({c}_{1}{c}_{2}{c}_{3}+{s}_{1}{s}_{2}{s}_{3})\left|00\right\rangle +({c}_{1}{c}_{2}{s}_{3}-{s}_{1}{s}_{2}{c}_{3})\left|01\right\rangle\\ \qquad \quad +\,({c}_{1}{s}_{2}{c}_{3}-{s}_{1}{c}_{2}{s}_{3})\left|10\right\rangle + ({s}_{1}{c}_{2}{c}_{3}+{c}_{1}{s}_{2}{s}_{3})\left|11\right\rangle \end{array}$$ (18)

where \({c}_{i}\equiv \cos {\theta }_{i}\) and \({s}_{i}\equiv \sin {\theta }_{i}\).

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 \(\left|0\right\rangle\) and rotated around \(y\) with the angles \({\theta }_{1},{\theta }_{2}\) and \({\theta }_{3}\) found by solving the system in Eq. (19), and they act as controls for a controlled \(X\) and a controlled \(Y\). Full size image

The angles \({\theta }_{i}\) can be found by solving the following system of equations:

$$\left\{\begin{array}{ll}{p}_{0}=| \langle 00| \psi \rangle {| }^{2}={({c}_{1}{c}_{2}{c}_{3}+{s}_{1}{s}_{2}{s}_{3})}^{2}&\\ {p}_{1}=| \langle 01| \psi \rangle {| }^{2}={({c}_{1}{c}_{2}{s}_{3}-{s}_{1}{s}_{2}{c}_{3})}^{2}&\\ {p}_{2}=| \langle 10| \psi \rangle {| }^{2}={({c}_{1}{s}_{2}{c}_{3}-{s}_{1}{c}_{2}{s}_{3})}^{2}&\\ {p}_{3}=| \langle 11| \psi \rangle {| }^{2}={({s}_{1}{c}_{2}{c}_{3}+{c}_{1}{s}_{2}{s}_{3})}^{2}&\end{array}\right..$$ (19)

The above system of equations allows for multiple analytical solutions whose expressions are too cumbersome to be reported here. The choice of the solution to use in each case depends on a number of factors, such as the gate fidelity for the specific values of the parameters.

Notice that the circuit for the Pauli channel can be used also for the depolarising channel, instead of the three-qubit circuit of Fig. 9, but it proves to be less accurate on the IBM Q Experience devices, presumably because of the required entanglement between the two qubits.