Introduction

The Hermiticity of physical observables is a fundamental tenant of standard quantum physics, guaranteeing real eigenspectra and leading to the generation of unitary dynamics in closed quantum systems. However, this is needlessly restrictive: it has been shown1 that non-Hermitian Hamiltonians endowed with parity and time (\({\mathcal{PT}}\)) symmetry possess real positive eigenvalues and eigenvectors with positive norm. Experimental platforms where non-Hermitian Hamiltonians can be implemented to comprise optical waveguides2,3,4, polarized photons5, nuclear spins6,7, superconducting circuits8,9, mechanical oscillators10, nitrogen-vacancy centers in diamond11,12, fiber-optic networks13, and ultracold Fermi gases14. Open systems are a natural candidate for realizing these Hamiltonians, since non-Hermitian terms appear naturally as a consequence of energy being injected or lost15.

However, a major drawback of open-system approaches is the need to control precisely the gain and the dissipation: these experiments require complicated setups with gain and alternating losses5, and yet only wave-like effects can be observed. Moreover, employing gain-loss systems for the study of quantum properties such as entropy, entanglement, and correlations is fundamentally impossible because gain inevitably adds noise16. Thus, in order to make progress one would need genuine realizations of non-Hermitian dynamics in the quantum regime4,9,12, which maintain and allow the measurement of delicate quantum effects.

Here we show that the non-Hermitian dynamics can be simulated digitally17 in a superconducting quantum processor by extending the Hilbert space with the use of an ancilla qubit and under the action of appropriately defined gates. To achieve this, we combine two techniques: a dilation method that is universal (applicable to any Hamiltonian)12 and an optimal method for generating any two-qubit gate with combinations of single-qubit gates and at most three CNOT (controlled-NOT) gates18,19. This combination enables us to observe and fully characterize the broken \({\mathcal{PT}}\)-symmetry transition and to settle decisively the relationship between non-Hermitian quantum mechanics and no-go theorems on state distinguishability and monotony of entanglement20,21,22,23. We achieve this by making use of the emergent technology of superconducting processors, on which significant technical progress has been shown in recent times by IBM24. Although still imperfect, these devices have already enabled important results such as quantum error correction25, fault-tolerant gates26, proofs of violation of Mermin27 and Leggett–Garg28 inequalities, demonstrations of non-local parity measurements29,30, simulations of paradigmatic models in open quantum systems31, the creation of highly entangled graph states32, the determination of molecular ground-state energies33, the implementation of quantum witnesses34, and quantum-enhanced solutions to large systems of linear equations35. The use of a superconducting quantum processor offers the possibility of extracting all relevant quantum correlations and of designing and programming efficiently the required gates, adapted to the particular topology of the machine.

We consider the generic system qubit non-Hermitian Hamiltonian36 with natural units  = 1

$${H}_{{\rm{q}}}={\sigma }_{x}+ir{\sigma }_{z},$$
(1)

where r is a real parameter and σx and σz are the Pauli matrices. The eigenvalues are \(\pm\! \sqrt{1-{r}^{2}}\), and the condition for non-Hermiticity is simply r ≠ 0. The parity operator is \({\mathcal{P}}={\sigma }_{x}\) and the time-inversion operator is the complex conjugation operator \({\mathcal{T}}=\star\). This Hamiltonian has an exceptional point at r = 1, where the two eigenvectors coalesce and the eigenvalues become parallel. For r < 1 the eigenvalues are real, corresponding to distinct eigenvectors, and the Hamiltonian satisfies \({\mathcal{PT}}\)-symmetry \([{\mathcal{PT}},{H}_{{\rm{q}}}]=0\) (see Supplementary Note 1); for r > 1, the eigenvalues become imaginary and the \({\mathcal{PT}}\) symmetry is broken. The Hamiltonian Eq. (1) can be understood as the standard non-Hermitian form providing equal coupling (off-diagonal terms) between the basis states as well as equal gain and loss via the complex diagonal terms required for \({\mathcal{PT}}\) symmetry37.

We realize the Hamiltonian Eq. (1) in a dilated space with the help of an ancilla, observing single-qubit dynamics under PT-symmetric Hamiltonians and witnessing the coalesce of eigenvectors at the exceptional points. Further, by allowing different quantum states to evolve under the same set of operations generated by non-Hermitian generators, we show that the trace distance between arbitrary states is modified—a task that is forbidden in Hermitian quantum mechanics. We extract the critical exponent of the transition, obtaining a value in agreement with theoretical predictions. We also observe an apparent violation of entanglement monotonicity in a two-qubit system, where one of the qubits is driven by a non-Hermitian Hamiltonian. Finally, we conclude by providing the complete dynamics of correlations developed between system qubits and the ancilla.

Results

To realize the Hamiltonian Eq. (1) we use a Naimark dilation procedure employing an additional ancilla qubit and a Hermitian operator \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) acting on the total qubit–ancilla Hilbert space. The dynamics under \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) is determined by the Schrödinger equation

$$i\frac{\mathrm{d}}{{\mathrm{d}}t}{\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}={{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t){\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}},$$
(2)

whose solution is given by

$${\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}={\left|0\right\rangle }_{{\rm{a}}}{\left|\psi (t)\right\rangle }_{{\rm{q}}}+{\left|1\right\rangle }_{{\rm{a}}}{\left|\tilde{\psi }(t)\right\rangle }_{{\rm{q}}},$$
(3)

where \({\left|\psi (t)\right\rangle }_{{\rm{q}}}\) is the solution of \(i\frac{\mathrm{d}}{{\mathrm{d}}t}{\left|\psi (t)\right\rangle }_{{\rm{q}}}={H}_{{\rm{q}}}{\left|\psi (t)\right\rangle }_{{\rm{q}}}\). Here \({\left|\tilde{\psi }(t)\right\rangle }_{{\rm{q}}}=\eta (t){\left|\psi (t)\right\rangle }_{{\rm{q}}}\), where η(t) is a positive linear operator given by \(\eta (t)={[(1+{\eta }_{0}^{2})\exp (-i{H}_{{\rm{q}}}^{\dagger }t)\exp (i{H}_{{\rm{q}}}t)-{\Bbb{I}}]}^{1/2}\)12, with \(\eta (0)={\eta }_{0}{\Bbb{I}}\) at the initial time t = 0 (see “Methods”).

In order to obtain a solution of the form Eq. (3), the ancilla and the qubit are initialized in the state \({\left|{{\Psi }}(0)\right\rangle }_{{\rm{a}},{\rm{q}}}=({\left|0\right\rangle }_{{\rm{a}}}+{\eta }_{0}{\left|1\right\rangle }_{{\rm{a}}})\otimes {\left|\psi (0)\right\rangle }_{{\rm{q}}}\), as shown in Fig. 1a, by using a rotation Ry(θ) on the ancilla qubit, where \(\theta =2\,{\tan }^{-1}\,{\eta }_{0}\) and \({\left|\psi (0)\right\rangle }_{{\rm{q}}}={\left|0\right\rangle }_{{\rm{q}}}\). Thus, the dynamics of the system qubit in the subspace with the ancilla in state \({\left|0\right\rangle }_{{\rm{a}}}\) satisfies the desired evolution by the non-Hermitian Hamiltonian Eq. (1). For an arbitrary time t and for any given r the corresponding unitary operator \({U}_{{\rm{a}},{\rm{q}}}(t)={\rm{T}}\exp [-i\mathop{\int}\nolimits_{0}^{t}{{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(\tau )\mathrm{d}\tau ]\) is obtained numerically as described below.

Fig. 1: Demonstration of \({\mathcal{PT}}\)-symmetry breaking for a single-qubit.
figure 1

a Quantum circuit implementing the dilated unitary operator Ua,q in the four-dimensional Hilbert space of system and ancilla, where Ry(θ) is a single-qubit rotation about the y-axis by an angle θ, and the subscript q(a) refers to the system qubit(ancilla). Details of the evolution under \({\mathcal{H}}_{{\rm{a}},{\rm{q}}}\) are given in the inset. b Dynamics of the population p0 of state \(\left|0\right\rangle\) under the Hamiltonian Hq for a range of parameters −1 ≤ r ≤ 1.5. As the eigenvectors of Hq coalesce at r = 1, the oscillations disappear, and beyond this point, \({\mathcal{PT}}\)-symmetry breaking is observed. The plots in ce present the results of the experiments (red dots) run on the IBM quantum experience (with 8192 repetitions) for the parameter r taking values r = 0.6, 1.0, and 1.3 respectively. The data matches well with theory (continuous black lines) corresponding to the three red lines from panel (b).

Parity–time symmetry breaking in a single qubit

We implement the unitary operator Ua,q(t) using the q[0] and q[1] qubits of a five-qubit IBM quantum processor for different values of the Hamiltonian parameter r. We start by presenting in Fig. 1b the expected theoretical results for the ground-state population obtained under Hq, that is, generated by the non-unitary evolution operator \(\exp (-i{H}_{{\rm{q}}}t)\) (see “Methods”). The breaking of \({\mathcal{PT}}\) symmetry, as one crosses the exceptional point r = 1, is clearly visible. Indeed, for r < 1 the eigenvalues of Hq are real and one observes Rabi oscillations. When r exceeds 1, the eigenvalues of Hq become imaginary, the \({\mathcal{PT}}\) symmetry is broken, and what one observes is the decay of the population. Figure 1c–e presents the results from the experimental realization of Hq on IBM quantum experience for three different values of r. We note that the agreement with the theoretical values is excellent. Each experiment is repeated 8192 times. Thus, the statistical errors here are of the order of \(1/\sqrt{8192}=0.01\), which lead to the error bars too small to be shown distinctly in the experimental plots. In addition, in various experiments presented here, we track the systematic errors in terms of measurement corrections and incorporate these corrections in respective experimental datasets as described in “Methods”, with further details given in the Supplementary Note 5.

The details of the implementation are shown in Fig. 1a. We start with the qubit and the ancilla both initialized in the state \(\left|0\right\rangle\), after which the ancilla alone is subjected to a rotation along the y-axis by an angle θ, which initializes the ancilla subspace θ12. The explicit form of the operator Ua,q(t) at any arbitrary time t is found by a numerical decomposition into single and two-qubit gates18. This decomposition Unum(t) = UnU1 matches with the desired unitary operator Ua,q(t) with a fidelity FU > 0.99, where the fidelity is defined as FU = 1 − Unum(t) − Ua,q(t)/Ua,q(t) (also see Supplementary Note 2). The quantum circuit that implements the decomposition of Unum(t) (see inset of Fig. 1a) comprises a sequence of single-qubit rotations \({{\mathcal{U}}}_{{\rm{q}}({\rm{a}})}^{j}\), each of them having up to three degrees of freedom, and three two-qubit CNOT gates18,19. The width of this circuit is 2 and the depth is 8. Specifically, the single-qubit gates are general rotations given by \({{\mathcal{U}}}_{{\rm{q}}({\rm{a}})}^{j}(\alpha ,\beta ,\gamma )={R}_{z}{(\alpha )}_{{\rm{q}}({\rm{a}})}^{j}{R}_{y}{(\beta )}_{{\rm{q}}({\rm{a}})}^{j}{R}_{z}{(\gamma )}_{{\rm{q}}({\rm{a}})}^{j}\), where α, β, and γ are the angles of rotations and the operators Ry, Rz correspond to the rotations generated by the Pauli operators σy and σz, respectively. The operator \({{\mathcal{U}}}_{{\rm{q}}({\rm{a}})}^{j}(\alpha ,\beta ,\gamma )\) has a direct correspondence with the single-qubit operator U3, as defined by IBM. For instance, given r = 0.6 and t = 0.5 (see Fig. 1a), we have the following set of operations: \({{\mathcal{U}}}_{{\rm{a}}}^{1}(2.83,0.55,3.72),\ {{\mathcal{U}}}_{{\rm{q}}}^{1}(0.51,-2.98,\,1.63)\), \({{\mathcal{U}}}_{{\rm{a}}}^{2}(-1.75,\,-3.34,\,-4.60),\ {{\mathcal{U}}}_{{\rm{q}}}^{2}(0.00,\,0.00,\,4.02)\), \({{\mathcal{U}}}_{{\rm{a}}}^{3}(4.81,3.08,-1.02),\ {{\mathcal{U}}}_{{\rm{q}}}^{3}(0.01,0.29,0.04)\), and \({{\mathcal{U}}}_{{\rm{a}}}^{4}(0.00,-5.19,0.50),\ {{\mathcal{U}}}_{{\rm{q}}}^{4}(0.46,-1.51,0.37)\). After the Unum(t) implementation, the post-selected subspace of our interest corresponds to the ancilla in state \({\left|0\right\rangle }_{a}\). At the end of the algorithm, we measure the probabilities Pkl of the qubit–ancilla state in the computational basis \(\{{\left|kl\right\rangle }_{{\rm{a}},{\rm{q}}}\equiv {\left|k\right\rangle }_{{\rm{a}}}{\left|l\right\rangle }_{{\rm{q}}}\}\) with k, l {0,1}. Finally, the ground-state population in the desired post-selected subspace of the system qubit is given by, p0(t) = P00/(P00 + P01), which can be obtained directly from the experiments. These are shown with red dots in Fig. 1c–e and follow very closely the results for the population in the \({\left|0\right\rangle }_{{\rm{q}}}\) state of the qubit under the non-Hermitian Hamiltonian Eq. (1). The results demonstrate a high-fidelity simulation of the \({\mathcal{PT}}\)-symmetry breaking in a single qubit.

Quantum state distinguishability

Next, we demonstrate an unexpected consequence of non-Hermitian dynamics concerning state distinguishability. Designing a general protocol to distinguish two (or more) arbitrary quantum states is a challenge in standard Hermitian quantum mechanics. On the other hand, the evolution of an arbitrary pair of states under a non-Hermitian operator can alter the distance between them, and may even make the arbitrary pair of quantum states orthogonal20,22,23. To observe this unusual feature of non-Hermitian dynamics, we use the quantum circuit in Fig. 1a to evolve the system qubit, initialized, respectively, in the orthogonal states \({\left|0\right\rangle }_{{\rm{q}}}\) and \({\left|1\right\rangle }_{{\rm{q}}}\). At various different instances of time, the state of the system qubit in the post-selected subspace with ancilla in state \({\left|0\right\rangle }_{{\rm{a}}}\) is obtained and the trace distance

$${\mathcal{D}}({\rho }_{1{\rm{q}}}({{t}}),{\rho }_{2{\rm{q}}}({{t}}))=\frac{1}{2}{\rm{tr}}\sqrt{{\rho }_{{\rm{diff}}}{({{t}})}^{\dagger }\ {\rho }_{{\rm{diff}}}({{t}})},$$
(4)

between the respective states is determined, where ρdiff(t) = ρ1q(t) − ρ2q(t) and \({\rho }_{i{\rm{q}}}(t)={\left|{\psi }_{i}(t)\right\rangle }_{{\rm{q}}}{\left\langle {\psi }_{i}(t)\right|}_{{\rm{q}}}\). For the given pair of initial states, the expected pattern for the variation of \({\mathcal{D}}\) with r and t is shown in Fig. 2a. The characteristic recurrence time TR in the \({\mathcal{PT}}\)-symmetric phase and the decay time \({\tau}_{\rm{D}}\) in the broken-symmetry phase are plotted in Fig. 2b and compared to their analytical expressions (\({T}_{\mathrm{R}}=\pi /\sqrt{1-{r}^{2}}\) and \(\tau_{\rm{D}} =1/2\sqrt{{r}^{2}-1}\) from Supplementary Eqs. (S11) and (S12)). Note that these times reflect the delicate balance between gain and loss, which is encoded in the structure of the Hamiltonian (see Supplementary Note 4). In Fig. 2c–e we show the experimentally obtained variation of the trace distance for three different values of r. An oscillating pattern in the trace distance is obtained for r < 1, which is a signature of information exchange between the system and the environment, while for r ≥ 1 we measure a decay pattern, which corresponds to loss of information to the environment. Interestingly, the oscillations in distinguishability correspond to oscillations in entanglement of qubit–ancilla state22 (see Supplementary Fig. 2). For r = 1 (exceptional point) these timescales diverge, and one cannot define anymore a characteristic time of the system. Instead, in close analogy with phase transitions in many-body systems, the distinguishability follows asymptotically a power-law \({\mathcal{D}} \sim {t}^{-\delta }\)22, where the critical exponent δ = 2 corresponds to two coalescing eigenstates. We have first checked numerically that for t 1 the distinguishability indeed displays this power-law behavior, with the critical exponent very close to 2. We can verify this scaling also experimentally, with the caveat that for t 3 the distingusihability becomes already smaller than the precision that we can reach on the IBM machine. Still, we can identify an interval t [1, 3] where the theoretical plot \(\mathrm{ln}\,{\mathcal{D}}\) versus \(\mathrm{ln}\,t\) starts to be approximately linear, with slope δ = 1.93 ± 0.08, see inset of Fig. 2d. In this region, we obtain by fitting the experimental data δ = 1.75 ± 0.15 (dashed red line in the inset), a reasonably close value.

Fig. 2: Quantum state distinguishability in a system driven by a non-Hermitian Hamiltonian.
figure 2

a Variation of trace distance Eq. (4) for single-qubit initial states \({\rho }_{1{\rm{q}}}(t=0)=\left|0\right\rangle \left\langle 0\right|\) and \({\rho }_{2{\rm{q}}}(t=0)=\left|1\right\rangle \left\langle 1\right|\), as they evolve under Hq. b Experimental (red dots) and theoretical recurrence (continuous black line) and decay times (continuous blue line) characterizing retrieval and loss of information in the two phases. The plots in ce present the theoretical (continuous black line) and experimental (red dots) curves of \({\mathcal{D}}({\rho }_{1{\rm{q}}}({\rm{t}}),{\rho }_{2{\rm{q}}}({\rm{t}}))\) for r = 0.6, 1.0, and 1.3, respectively. The inset in d shows \(\mathrm{ln}\,{\mathcal{D}}\) versus \(\mathrm{ln}\,t\) for t [1, 3], where the dotted red line is a linear interpolation of the data. Each experimental point is obtained by averaging over 8192 samples, such that the errors are below the size of the marker.

Evaluating the distinguishability requires a complete characterization of the single-qubit density matrices ρ1q(t) and ρ2q(t), which is done by a set of three operations that independently fetch the elements of the density matrix. Each of these experiments is repeated 8192 times, such that even after evaluating Eq. (4), the statistical error in the measure of distinguishability remains small.

Bipartite systems under non-Hermitian evolution

Next, we observe the dynamics of entanglement in a bipartite system when one of the parties undergo a local operation generated by Hq, for different values of r. Such scenarios have been studied theoretically21,22, and it was shown that entanglement restoration and information recovery can happen in the \({\mathcal{PT}}\)-symmetric phase. This breaks entanglement monotonicity, allowing the creation of entanglement by a local operation. This unusual effect is due to the modified evolution in the post-selected subspace due to mere existence of a component of the total wavefunction outside this subspace4,38,39.

To study this phenomenon, we consider a system consisting of two qubits q and \({\rm{q}}^{\prime}\) initialized in a maximally entangled Bell state, \({\left|{{{\Phi }}}^{+}\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }=({\left|00\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }+{\left|11\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} })/\sqrt{2}\). One system qubit (say q) undergoes a non-Hermitian evolution by \({{\Bbb{H}}}_{{\rm{q}},{\rm{q}}^{\prime} }={H}_{{\rm{q}}}\otimes {{\Bbb{I}}}_{{\rm{q}}^{\prime} }\) with the help of an ancillary qubit a, such that the total Hamiltonian including the dilation is \({{\Bbb{H}}}_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }={{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}\otimes {{\Bbb{I}}}_{{\rm{q}}^{\prime} }\) leading to a unitary evolution, \({{\Bbb{U}}}_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }={U}_{{\rm{a}},{\rm{q}}}\otimes {{\Bbb{I}}}_{{\rm{q}}^{\prime} }\). Finally, the three-partite state of the system is measured and post-selected subject to the state of the ancilla being \({\left|0\right\rangle }_{{\rm{a}}}\).

The experimental implementation on the IBM quantum processor is carried out using three qubits, as shown in Fig. 3a. As before, we average over 8192 realizations. We perform the complete quantum state tomography of the two-qubit system in the post-selected subspace at various different values of time t. This is done using a set of seven experiments on the system qubits q and \({\rm{q}}^{\prime}\), followed by σz measurements of all three qubits as shown in Fig. 3a—see Supplementary Note 5 for further details. At the end of each of these tomography measurements, the populations are obtained as \({p}_{kl}={P}_{0kl}/\mathop{\sum }\nolimits_{m,n = 0}^{1}{P}_{0mn}\), where k, l {0, 1}, from which the desired density operator of the system qubits \({\rho }_{{\rm{q}},{\rm{q}}^{\prime} }^{(0)}\) in the post-selected subspace is obtained. To study the entanglement dynamics, we use the concurrence40,41 as a measure of entanglement, given by \({{\mathcal{C}}}_{{\rm{q}},{\rm{q}}^{\prime} }^{(0)}=\max \{0,\sqrt{{\lambda }_{1}}-\sqrt{{\lambda }_{2}}-\sqrt{{\lambda }_{3}}-\sqrt{{\lambda }_{4}}\}\), where λi’s are the eigenvalues of the operator \({\rho }_{{\rm{q}},{\rm{q}}^{\prime} }^{(0)}({\sigma }_{y}\otimes {\sigma }_{y}){({\rho }_{{\rm{q}},{\rm{q}}^{\prime} }^{(0)})}^{\star }({\sigma }_{y}\otimes {\sigma }_{y})\) written in decreasing order.

Fig. 3: Demonstration of an apparent violation of entanglement monotonicity.
figure 3

a Quantum circuit implemented on the superconducting quantum processor, where the qubit a serves as ancilla and q and \({\rm{q}}^{\prime}\) form a qubit–qubit bipartite system. In bd we present the results from the experiments, where the variation of concurrence with time is shown for different values of r, with each experiment being repeated 8192 times. Experimental data (red dots) very closely follow the theoretically expected behavior (continuous black lines). The inset in c presents the variation of the logarithm of concurrence versus \(\mathrm{ln}\,t\) in the interval t [1, 3], where the dashed red line is a linear fit to the experimental red circles. e We prepare the qubits q and q′ in a state with concurrence 0.475 and we evolve the system under the local non-Hermitian evolution with up to t = 1.75, for different values of r. Solid blue, red, and black curves present the theoretically expected results, while blue circle, red square, and black diamond markers correspond to the respective experimental measurements. In f we show the theoretical plots for the variation of concurrence with time, with one of the qubits undergoing a local non-Hermitian evolution, for various parameters r [0, 1.5]; the horizontal green lines mark the values of r corresponding to the experimental results in panels (bd). The correlations among the system qubits and ancilla, simulated based on the quantum circuit in a are presented in g for r = 0.6, where the concurrence between the system qubits and between the system qubit and ancilla are shown with dotted and continuous red lines, respectively, the dot-dashed and continuous black lines present the linear entropy of the system qubit and of the ancilla, and the continuous blue line corresponds to the three-tangle entanglement.

The change of concurrence with time is then observed for different values of r, as shown in Fig. 3b–d. For r = 0 we have checked that the dynamics is unitary and there is no variation in the entanglement values. In this case, the standard result that the entanglement is not changed by local operations is confirmed. However, for 0 < r < 1, the concurrence is found to be oscillating, which is clearly seen in Fig. 3b, while for r  > 1 the Hamiltonian Hq governing local evolution ceases to obey the symmetry, and the entanglement gradually decays with time, as shown in Fig. 3d. For r = 1 we find the same theoretical asymptotic scaling as for distinguishability \({{\mathcal{C}}}_{{\rm{q}},{\rm{q}}^{\prime} }^{(0)} \sim {t}^{-\delta }\), where δ = 2. To compare with the experiment, again we restrict the time to t [1,  3], and find δ = −1.71 ± 0.01 from the theoretical curve and δ = −1.93 ± 0.27 from the measured data (see inset in Fig. 3c). In Fig. 3f we present the corresponding theoretical curves for the time variation of concurrence for a wider range of r parameters.

The obtained variation in concurrence under a local operation contradicts at first sight the well-known property of monotonicity of entanglement. To make this effect even more striking, we have performed another experiment where we observe the increase in entanglement between the qubits q and \({\rm{q}}^{\prime}\) under the action of the \({\mathcal{PT}}\)-symmetric non-Hermitian Hamiltonian. Specifically, at t = 0 we prepare the state \(\cos (\vartheta ){\left|{{{\Phi }}}^{-}\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }-i\sin (\vartheta ){\left|{{{\Psi }}}^{+}\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }\), where \({\left|{{{\Psi }}}^{\pm }\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }=({\left|01\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }\pm {\left|10\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} })/\sqrt{2}\) and \({\left|{{{\Phi }}}^{\pm }\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }=({\left|00\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} }\pm {\left|11\right\rangle }_{{\rm{q}},{\rm{q}}^{\prime} })/\sqrt{2}\) are the standard maximally entangled Bell states. The angle ϑ defines the concurrence \(| \cos (2\vartheta )|\) of this state. For the experiment—shown in Fig. 3e—we took ϑ = 59.185°, yielding a concurrence of 0.475 at t = 0. The preparation of this state is done by single- and two-qubit gates acting on the two qubits; the ancilla is not involved and remains separate in the state \({\left|0\right\rangle }_{{\rm{a}}}\). Next, we simulate the action of the non-Hermitian Hamiltonian for 0 ≥ t < 2 and r = 0.3, r = 1, r = 1.3, see Fig. 3e. In the first case, the entanglement increases up to ~0.8 (and would continue to oscillate at longer times), while for r = 1, r = 1.3 it decreases monotonously.

Entanglement correlations between system and ancilla

The simulation of non-Hermiticity by the dilation method allows us to get an in-depth understanding of this phenomenon. Let us look at the complete picture in the eight-dimensional Hilbert space of this tripartite system (initialized in the state \(\left|0\right\rangle \otimes \left|{{{\Phi }}}^{+}\right\rangle\)), where, as we have seen in Fig. 3a, one of the system qubits q along with the ancillary qubit “a” undergo the unitary evolution Ua,q. The relevant correlations for the ensuing analysis are plotted in Fig. 3g for r = 0.6. We define the single-qubit reduced states by tracing out the other qubits \({\rho }_{i}={{\rm{Tr}}}_{j,h}[{\rho }_{i,j,h}]\), while the two-qubit reduced density operators are obtained by a single partial trace operation \({\rho }_{i,j}={{\rm{Tr}}}_{h}[{\rho }_{i,j,h}]\), where the three qubits are labeled by \(i,j,h\in \{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} \}\). The concurrence associated with the state ρi,j is denoted by \({{\mathcal{C}}}_{i,j}\) and it is calculated using the formula for mixed two-qubit states mentioned earlier. It is interesting to note that q and \({\rm{q}}^{\prime}\) are always in the permutation symmetric subspace of the two-qubit Hilbert space as one of the qubits evolves under Hq (see Supplementary Note 3). Therefore, it is enough to observe any one of the system qubits. Analyzing first the single-qubit states, we find that the single-qubit reduced density operators ρq and \({\rho }_{{\rm{q}}^{\prime} }\) remain maximally mixed all through the evolution, with a constant value of linear entropy \({s}_{{\rm{q}}}=1-{\rm{Tr}}({\rho }_{q}^{2})=0.5\).

Next, we observe that the concurrences \({{\mathcal{C}}}_{{\rm{a}},{\rm{q}}}\) and \({{\mathcal{C}}}_{{\rm{a}},{\rm{q}}^{\prime} }\) between the ancilla and the respective system qubits, that is, a and q (or a and \({\rm{q}}^{\prime}\)) always remain zero. This shows that the dynamics under \({{\mathbb{U}}}_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }\) does not develop bipartite correlations between the respective system qubits and the ancilla qubit. Therefore, the creation of a tripartite correlation between the system and the ancilla can happen only through entangling correlations between the two-qubit reduced state of \({\rm{q}},{\rm{q}}^{\prime}\) and the ancilla. To quantify this tripartite correlation, we use the three-tangle for pure states41

$${\tau }_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }= \, {{\mathcal{C}}}_{{\rm{a}}:{\rm{q}},{\rm{q}}^{\prime} }^{2}-{{\mathcal{C}}}_{{\rm{a}}:{\rm{q}}}^{2}-{{\mathcal{C}}}_{{\rm{a}}:{\rm{q}}^{\prime} }^{2},\, {\mathrm{or}}\\ {\tau }_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }= \, {{\mathcal{C}}}_{{\rm{q}}:{\rm{a}},{\rm{q}}^{\prime} }^{2}-{{\mathcal{C}}}_{{\rm{q}}:{\rm{a}}}^{2}-{{\mathcal{C}}}_{{\rm{q}}:{\rm{q}}^{\prime} }^{2},$$
(5)

where in the last equation we used the invariance of the tangle under permutations. As shown by Eq. (5), the maximum value of the three tangle is obtained in the absence of concurrence between the individual components. Here \({{\mathcal{C}}}_{{\rm{q}}:{\rm{a}}}\equiv {{\mathcal{C}}}_{{\rm{q}},{\rm{a}}}\) and \({{\mathcal{C}}}_{{\rm{q}}:{\rm{q}}^{\prime} }\equiv {{\mathcal{C}}}_{{\rm{q}},{\rm{q}}^{\prime} }\) are the concurrences of the two-party reduced states ρa,q and \({\rho }_{{\rm{q}},{\rm{q}}^{\prime} }\). The first term on the right-hand side of Eq. (5) is the square of concurrence between the bipartitions \({\rho }_{{\rm{q}}}:{\rho }_{{\rm{a}},{\rm{q}}^{\prime} }\), where one partition consists of the qubit q, while the other partition is formed by the ancilla a and the system qubit \({\rm{q}}^{\prime}\). For a pure three-qubit state \({\rho }_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }\), the quantity \({{\mathcal{C}}}_{{\rm{q}}:{\rm{a}},{\rm{q}}^{\prime} }\) is effectively related to the mixedness of its bipartitions. More specifically, the square of concurrence between the partitions ρq and \({\rho }_{{\rm{a}},{\rm{q}}^{\prime} }\) is twice the linear entropy of the reduced density operator of either partition, given by \(2(1-{\rm{Tr}}{\rho }_{q}^{2})\) or \(2(1-{\rm{Tr}}{\rho }_{{\rm{a}},{\rm{q}}^{\prime} }^{2})\). We now know from the simulated dynamics that the linear entropy \({s}_{{\rm{q}}({\rm{q}}^{\prime} )}=1-{\rm{Tr}}{\rho }_{{\rm{q}}({\rm{q}}^{\prime} )}^{2}=0.5\); therefore, at all times \({{\mathcal{C}}}_{{\rm{q}}:{\rm{a}},{\rm{q}}^{\prime} }^{2}=1\). Further, as shown in Fig. 3g, there is no bipartite entanglement between the ancilla and the respective system qubits q (or q′), which implies that \({{\mathcal{C}}}_{{\rm{q}},{\rm{a}}({\rm{a}},{\rm{q}}^{\prime} )}=0\). From Eq. (5), we obtain,

$${\tau }_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }=1-{{\mathcal{C}}}_{{\rm{q}},{\rm{q}}^{\prime} }^{2}.$$
(6)

Thus, the three tangle among system qubits and ancilla and the concurrence between the system qubits are complementary to each other (see Supplementary Note 3). By permuting the partitions in Eq. (5), it is easy to obtain \({\tau }_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }=2{s}_{{\rm{a}}}\), where \({s}_{{\rm{a}}}=1-{\rm{Tr}}({\rho }_{{\rm{a}}}^{2})\). The unitary \({{\Bbb{U}}}_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }\), which induces a local non-Hermitian drive of qubit q in the post-selected subspace of the ancilla, is in fact a nonlocal operation on the system qubit q and the ancilla a. Under \({{\Bbb{U}}}_{{\rm{a}},{\rm{q}},{\rm{q}}^{\prime} }\), as the ancilla entangles and dis-entangles with the joint state of the system qubits, we see the resulting oscillations of various correlations in time. These oscillating correlations with r-dependent characteristic times, when post-selected in the ancilla subspace, produce an apparent violation of entanglement monotonicity. While we observe experimentally the variation in entanglement under local operations in only one post-selected subspace, other subspaces of the same system also witnesses similar patterns for the variation of entanglement under local operations as shown in Supplementary Fig. 3.

Conclusion

We have realized a quantum simulation of a single-qubit under a non-Hermitian Hamiltonian, observing the \({\mathcal{PT}}\)-symmetry breaking as the exceptional point is crossed and the associated change in distinguishability. The use of a quantum processor for the simulation has the advantage that more complex scenarios can be studied, such as a bipartite system with one of the qubits driven by a non-Hermitian Hamiltonian. In this case, we observe the violation of the entanglement monotonicity no-go result from standard quantum mechanics. We also note that, while our method relies on dilation by the use of an ancilla, another approach to non-Hermitian evolution exists, where the dimension of the Hilbert space remains the same but the standard inner product is modified (see “Methods”). These two methods can be put in an exact correspondence—the metric used in the latter approach can be identified as the operator \({\eta }^{2}+{\Bbb{I}}\) from the dilation approach.

The simulation of phenomena governed by \({\mathcal{PT}}\)symmetry at the single-quantum level open up several novel perspectives. Our scheme provides a systematic way of studying more complex non-Hermitian many-qubit systems. It is important to realize that for a system of N qubits the overhead in the width of the circuit is just one ancilla qubit. For example, it would be straightforward to generalize to the study of entanglement that we have performed to one non-Hermitian qubit and N − 1 Hermitian ones, in which case the depth of the circuit remains equal to 8. Furthermore, because we have access to the quantum regime, our scheme enables the study of quantum fluctuations. Since these systems are open—connected to a source of energy providing gain and reservoir for dumping this energy—they naturally will lead to new insights into quantum thermodynamics.

Methods

Simulating the non-Hermitian Hamiltonian in the dilated space

The single-qubit evolution under a general time-dependent non-Hermitian Hamiltonian Hq(t) is obtained in a certain subspace of an ancilla–qubit system undergoing a unitary evolution generated by \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\). The Hamiltonian \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) in a four-dimensional Hilbert space can be obtained from Hq by Naimark dilation12. Using this method, we can write the Hamiltonian \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) in the form:

$${{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)={\Bbb{I}}\otimes {{\Lambda }}(t)+{\sigma }_{y}\otimes {{\Gamma }}(t),$$
(7)

with

$${{\Lambda }}(t)=\left[{H}_{{\rm{q}}}(t)+i\frac{{{\mathrm{d}}}\eta (t)}{{\mathrm{d}}t}\eta (t)+\eta (t){H}_{{\rm{q}}}(t)\eta (t)\right]{{{M}}}^{-1}(t),$$
(8)
$${{\Gamma }}(t)=i\left[{H}_{{\rm{q}}}(t)\eta (t)-\eta (t){H}_{{\rm{q}}}(t)-i\frac{{\mathrm{d}}\eta (t)}{{\mathrm{d}}t}\right]{{{M}}}^{-1}(t),$$
(9)
$$\eta (t)={({{M}}(t)-{\Bbb{I}})}^{\frac{1}{2}},\quad {\rm{and}}$$
(10)
$${{M}}(t)={{T}}\exp \left[-i\int_{0}^{t}{\mathrm{d}}\tau {H}_{\mathrm{q}}^{\dagger }(\tau )\right]{{M}}(0)\widetilde{{{T}}}\exp \left[i\int_{0}^{t}{\mathrm{d}}\tau {H}_{\mathrm{q}}(\tau )\right],$$
(11)

where η(t) and M(t) are Hermitian operators; T and \(\widetilde{{{T}}}\) are time-ordering and anti-time-ordering operators, respectively, and \({\Bbb{I}}\) is the 2 × 2 identity operator. The Hamiltonian \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) can be obtained as follows. First, the Eqs. (10) and (11)) for η(t) and M(t) reflect the invariance of the norm of \({\left|{{\Psi }}\right\rangle }_{{\rm{a}},{\rm{q}}}(t)\) in the form Eq. (3) on the dilated space under evolution (see also the discussion about metric below). Then, the Schrödinger equations \(i({\mathrm{d}}/{\mathrm{d}}t){\left|{{\Psi }}\right\rangle }_{{\rm{a}},{\rm{q}}}={{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t){\left|{{\Psi }}\right\rangle }_{{\rm{a}},{\rm{q}}}\), and \(i({\mathrm{d}}/{\mathrm{d}}t){\left|\psi (t)\right\rangle }_{{\rm{q}}}={H}_{{\rm{q}}}(t){\left|\psi (t)\right\rangle }_{{\rm{q}}}\) together with Eqs. (7) and (3) produce a linear system of equations in the unknown operators Λ(t) and Γ(t),

$${{\Lambda }}(t)-i{{\Gamma }}(t)\eta (t)={H}_{{\rm{q}}}(t),$$
(12)
$${{\Lambda }}(t)\eta (t)+i{{\Gamma }}(t)=i\frac{{\mathrm{d}}\eta (t)}{{\mathrm{d}}t}+\eta (t){H}_{{\rm{q}}}(t).$$
(13)

By multiplying the second equation to the right with η(t) and with −η−1(t) and adding the results to the first equation, we obtain immediately the solution as given by Eqs. (8) and (9). Note also that for Hermitian Hamiltonians Hq the second term in Eq. (7), which is the qubit–ancilla interaction, becomes zero.

Initial conditions

To obtain an explicit form of \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\), one should choose the operator M(t) at time t = 0 such that \({{M}}(t)-{\Bbb{I}}\) is positive for all t in the desired time interval. As a preliminary choice, we can take

$${{M}}(t=0)={{{M}}}_{0}={m}_{0}\times {\Bbb{I}},$$
(14)

where m0 > 1, may be chosen arbitrarily. Further, we obtain the eigenvalues of M(t) in the desired time interval. Fixing the value of r, at any arbitrary time t, the eigenvalues of M(t) are labeled as μ1(t) and μ2(t), from where we numerically obtain \({\mu }_{\min }(t)=\min \{{\mu }_{1}(t),{\mu }_{2}(t)\}\). Interestingly, for r = 0, Hq is Hermitian and

$${{M}}(t)={m}_{0}\times {\Bbb{I}}\quad \quad \forall \ t,$$

with eigenvalues m0, which is the maximum value that μmin can assume. Therefore, for any arbitrary r and t, \({m}_{0}/{\mu }_{\min }\ge 1\). Thus, at t = 0, M(t) is chosen to be,

$${{M}}(t=0)={{{M}}}_{0}=\frac{{m}_{0}}{{\mu }_{\min }}f\times {\Bbb{I}},$$
(15)

where f > 1, which ensures that \({{M}}(t)-{\Bbb{I}}\) remains positive for all t. From Eq. (10) we have

$$\eta (t=0)={\left(\frac{{m}_{0}}{{\mu }_{\min }}f-1\right)}^{\frac{1}{2}}\times {\Bbb{I}}={\eta }_{0}\times {\Bbb{I}}.$$
(16)

The dynamics of the total ancilla–qubit system under \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) is obtained from the Schrödinger equation

$$i\frac{\mathrm{d}}{{\mathrm{d}}t}{\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}={{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t){\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}},$$
(17)

whose solution is given by

$${\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}={\left|0\right\rangle }_{{\rm{a}}}{\left|\psi (t)\right\rangle }_{{\rm{q}}}+{\left|1\right\rangle }_{{\rm{a}}}{\left|\tilde{\psi }(t)\right\rangle }_{{\rm{q}}},$$
(18)

where \({\left|\psi (t)\right\rangle }_{{\rm{q}}}\) is the solution of \(i\frac{\mathrm{d}}{{\mathrm{d}}t}{\left|\psi (t)\right\rangle }_{{\rm{q}}}={H}_{{\rm{q}}}{\left|\psi (t)\right\rangle }_{{\rm{q}}}\), and \({\left|\tilde{\psi }(t)\right\rangle }_{{\rm{q}}}=\eta (t){\left|\psi (t)\right\rangle }_{{\rm{q}}}\). At t = 0 the state of the total system is

$${\left|{{\Psi }}(0)\right\rangle }_{{\rm{a}},{\rm{q}}}= \;{\left|0\right\rangle }_{{\rm{a}}}{\left|\psi (0)\right\rangle }_{{\rm{q}}}+{\left|1\right\rangle }_{{\rm{a}}}{\left|\tilde{\psi }(0)\right\rangle }_{{\rm{q}}},\\ = \;{\left|0\right\rangle }_{{\rm{a}}}{\left|\psi (0)\right\rangle }_{{\rm{q}}}+{\left|1\right\rangle }_{{\rm{a}}}\eta (0){\left|\psi (0)\right\rangle }_{{\rm{q}}},\\ = \;\left({\left|0\right\rangle }_{{\rm{a}}}+{\eta }_{0}{\left|1\right\rangle }_{{\rm{a}}}\right)\otimes {\left|\psi (0)\right\rangle }_{{\rm{q}}}={\left|\psi (0)\right\rangle }_{{\rm{a}}}\otimes {\left|\psi (0)\right\rangle }_{{\rm{q}}},$$
(19)

which is a separable state of the ancilla \({\left|\psi (0)\right\rangle }_{{\rm{a}}}\) and the system qubit \({\left|\psi (0)\right\rangle }_{{\rm{q}}}\). For preparing the initial state Eq. (19) the ancilla is taken in one of the eigenvectors of σz, say \({\left|0\right\rangle }_{{\rm{a}}}\). This is then subjected to a rotation by an angle θ around the y-axis, \({R}_{y}(\theta )=\exp (-i\theta {\sigma }_{y}/2)\), where, \(\theta =2{\tan }^{-1}{\eta }_{0}\). This leads to

$${\left|\psi (0)\right\rangle }_{{\rm{a}}}= \;\cos \frac{\theta }{2}{\left|0\right\rangle }_{{\rm{a}}}+\sin \frac{\theta }{2}{\left|1\right\rangle }_{{\rm{a}}}=\cos \frac{\theta }{2}\left({\left|0\right\rangle }_{{\rm{a}}}+\tan \frac{\theta }{2}{\left|1\right\rangle }_{{\rm{a}}}\right),\\ = \;\frac{1}{\sqrt{{\eta }_{0}^{2}+1}}\left({\left|0\right\rangle }_{{\rm{a}}}+{\eta }_{0}{\left|1\right\rangle }_{{\rm{a}}}\right),$$
(20)

which is the initial state as defined by the protocol. On the other hand, the qubit q may be initialized in any arbitrary state \({\left|\psi (0)\right\rangle }_{{\rm{q}}}\). For the case of a single qubit, as discussed in the first part of the paper, we considered two different values of the state of the qubit: (i) \({\left|\psi (0)\right\rangle }_{{\rm{q}}}=\left|0\right\rangle\) and (ii) \({\left|\psi (0)\right\rangle }_{{\rm{q}}}=\left|1\right\rangle\). The same formalism applies also in the case of the two-qubit system discussed in the second part of the paper, in which case the qubit–qubit system is initialized in a maximally entangled Bell state \({\left|\psi (0)\right\rangle }_{{\rm{q}}}\to \left|{{{\Phi }}}^{+}\right\rangle =\frac{1}{\sqrt{2}}(\left|00\right\rangle +\left|11\right\rangle )\).

For instance, in case (i) we choose m0 = 2 and f = 1.01. At r = 0.6 in the time range t [0, 8], we obtain η0 = 1.7436 and θ = 2.1001 (radians). At the exceptional point, that is, r = 1, in the same time range t [0, 8], we get η0 = 16.1112 and θ = 3.0176 radians. Further, for r = 1.3, \({\mu }_{\min }\) is obtained separately for various time intervals to increase the probability of success. This then led to different values of η0 and hence θ for each time point.

The metric

Non-Hermitian quantum dynamics can be alternatively formulated by using Hilbert spaces with a modified bra vector, resulting in a redefinition of the inner product23,36,42. Here we make an explicit connection with this approach, showing that the Hermitian operator \(M(t)=\eta {(t)}^{2}+{\Bbb{I}}\) can be identified as the metric that plays a key role in this formalism.

Indeed, from Eq. (18) we can calculate the norm of the dilated vector \({\left|{{\Psi }}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}\),

$${\,}_{{\rm{a}},{\rm{q}}}{\langle {{\Psi }}(t)| {{\Psi }}(t)\rangle }_{{\rm{a}},{\rm{q}}}={\,}_{{\rm{q}}}{\langle \psi (t)| M(t)| \psi (t)\rangle }_{{\rm{q}}}$$
(21)

This norm has to be conserved during the time evolution. By taking the time derivative of Eq. (21) and using Eq. (2) we get

$$i\frac{\mathrm{d}}{\mathrm{d}t}M(t)={H}_{\mathrm{q}}^{\dagger }(t)M(t)-M(t){H}_{\mathrm{q}}(t).$$
(22)

This is the defining relation for the metric23. Note that this can also be obtained in a straigthforward way from Eq. (11). Thus, in this approach to non-Hermitian quantum mechanics for every vector \({\left|\psi (t)\right\rangle }_{{\rm{q}}}\) in the Hilbert space we define the covector as \({\,}_{{\rm{q}}}\left\langle \psi (t)\right|M(t)\), which ensures that the inner product qψ(t)M(t)ψ(t)〉q from Eq. (21) has the meaning of a conserved probability.

For the particular case of the Hamiltonian Hq studied in this work, the metric M(t) can be obtained analytically by employing the properties of 2 × 2 matrices (see also Supplementary Eq. (3)). In the \({\mathcal{PT}}\)-symmetric phase we obtain an exact formula for the metric,

$$M(t)/{M}_{0}= \,\frac{1}{1-{r}^{2}}\left[1-{r}^{2}\cos \left(2\sqrt{1-{r}^{2}}t\right)\right]{\mathbb{I}}\\ +\frac{r}{\sqrt{1-{r}^{2}}}\sin \left(2\sqrt{1-{r}^{2}}t\right){\sigma }_{z}\\ +\frac{r}{1-{r}^{2}}\left[1-\cos 2\sqrt{1-{r}^{2}}t\right]{\sigma }_{y}$$
(23)

One can check also that M(t) is positively defined for r < 1, while for r = 0 we recover the standard Hermitian quantum mechanics with M(t) = M0. The result above Eq. (23) can be also obtained from the generic formula for the metric, as per ref. 23, for parameters A = 0, B = − r/(1 − r2), C = 1/(1 − r2), and D = 0.

It is interesting to remark how the main problem of non-Hermitian quantum mechanics, that of nonconservation of probability, has been dealt with in completely different ways: either by the introduction of a metric and modifying the inner product, or, in the dilation method, by adding an ancilla that absorbs the excess population.

Evolution and measurement in the dilated space

Let us consider the evolution of an arbitrary state of a two-qubit system under the Hamiltonian \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) in Eq. (7),

$$i\frac{\mathrm{d}}{{\mathrm{d}}t}{\left|\psi (t)\right\rangle }_{{\rm{a}},{\rm{q}}}={{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t){\left|\psi (0)\right\rangle }_{{\rm{a}},{\rm{q}}},$$

where ψ(0) is the initial state at t = 0. This may also be written as

$${\left|\psi (t)\right\rangle }_{{\rm{a}},{\rm{q}}}= \;{{T}}\exp \left[i\int_{0}^{t}{{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(\tau ){\mathrm{d}}\tau \right]{\left|\psi (0)\right\rangle }_{{\rm{a}},{\rm{q}}}\\ = \;{U}_{{\rm{a}},{\rm{q}}}(t){\left|\psi (0)\right\rangle }_{{\rm{a}},{\rm{q}}},$$
(24)

where T is the time-ordering operator. For a given set of values of r and t, it is useful to obtain an explicit form of the unitary operator Ua,q(t). This is done by observing the respective \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) evolutions of the complete set of two-qubit basis states. To find Ua,q(t) we solve the Schrödinger equation numerically for different initial states,

$$\left\{\begin{array}{ll}i\frac{\mathrm{d}}{{\mathrm{d}}t}{\left|{\psi }_{kl}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}={{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t){\left|{\psi }_{kl}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}&\\ \left|{\psi }_{kl}(0)\right\rangle ={\left|kl\right\rangle }_{{\rm{a}},{\rm{q}}}\hfill&k,l=0,1.\end{array}\right.$$
(25)

where \({\left|00\right\rangle }_{{\rm{a}},{\rm{q}}}\), \({\left|01\right\rangle }_{{\rm{a}},{\rm{q}}}\), \({\left|10\right\rangle }_{{\rm{a}},{\rm{q}}}\), and \({\left|11\right\rangle }_{{\rm{a}},{\rm{q}}}\) correspond to the complete set of basis vectors in the four-dimensional Hilbert space. The system qubit and the ancilla are initialized in all four bases states, respectively, and then evolved numerically under \({{\mathcal{H}}}_{{\rm{a}},{\rm{q}}}(t)\) for a given time. Then, after solving this equation for \(\left|{\psi }_{00}(t)\right\rangle\), \(\left|{\psi }_{01}(t)\right\rangle\), \(\left|{\psi }_{10}(t)\right\rangle\), and \(\left|{\psi }_{11}(t)\right\rangle\), we obtain the closed form of the unitary operator at an arbitrary time t, given by

$${U}_{{\rm{a}},{\rm{q}}}(t)= \, {\left|{\psi }_{00}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}\left\langle 00\right|+{\left|{\psi }_{01}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}\left\langle 01\right|\\ \, +{\left|{\psi }_{10}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}\left\langle 10\right|+{\left|{\psi }_{11}(t)\right\rangle }_{{\rm{a}},{\rm{q}}}\left\langle 11\right|.$$
(26)

For different values of time, Ua,q(t) is obtained, which is a general unitary operator in the four-dimensional Hilbert space. Each Ua,q(t) at a given time is then decomposed numerically in the form of single-qubit rotations and two-qubit CNOT gates, as shown in Fig. 1e.

This quantum circuit decomposition gives rise to Unum(t), whose operation is very close to the theoretical Ua,q(t). To characterize this, we calculate the error function errU(t) = Ua,q(t) − Unum(t)2/Ua,q(t)2, with the 2-norm defined by \(| | A| {| }_{2}=\sqrt{{\lambda }_{\max }}\), where \({\lambda }_{\max }\) is the largest eigenvalue of the matrix A*A. Here Unum(t) is an unitary operator generated by the circuit in the inset of Fig. 1a, where the parameters α, β, γ of \({{\mathcal{U}}}_{{\rm{q}}({\rm{a}})}^{j}(\alpha ,\beta ,\gamma )\) are chosen to minimize the expression Ua,q(t) − Unum(t)2. Typically, we find errU(t) = Ua,q(t) − Unum(t)2/Ua,q(t)2 to be of the order of 10−4, which demonstrates the high accuracy of our Ua,q implementation. The accuracy with which our gate decomposition and the Ua,q(t) operator match with each other is presented by an example data set in the Supplementary Table 1. Ua,q(t) for arbitrary values of (r, t) and the corresponding Unum(t) can be obtained from a GitHub code repository43.

Quantum state reconstruction

For single-qubit tomography we take 4 − 1 = 3 measurements, corresponding to the set of Pauli operators σx, σy, σz. For higher-dimensional quantum systems of n qubits, we need 22n − 1 measurements corresponding to combinations of σx, σy, σz and the identity matrix \({\Bbb{I}}\) of the two qubits. Therefore, a complete quantum state tomography of a two-qubit system requires a set of (16 − 1) experiments, which correspond to determining the expectation values of all the two-qubit operators formed by products of Pauli operators and the identity. In the present work, we need only to examine the post-selected subspace of the total system with the ancilla in state \(\left|0\right\rangle\). Therefore, we circumvent the complexities of three-qubit tomography by restricting our measurement to a 4 × 4 block of the complete 8 × 8 three-qubit density operator.

We perform a complete quantum state tomography of the system qubits by applying the following seven operators, namely \({{\Bbb{T}}}_{1} = {\Bbb{I}}\otimes {\Bbb{I}}, {{\Bbb{T}}}_{2} = H\otimes {\Bbb{I}}, {{\Bbb{T}}}_{3} = {R}_{x}(\pi /2)\otimes {\Bbb{I}},{{\Bbb{T}}}_{4} = {\Bbb{I}}\otimes H,{{\Bbb{T}}}_{5}={\Bbb{I}}\otimes {R}_{x}(\pi /2),{{\Bbb{T}}}_{6}=({\Bbb{I}}\otimes H){\rm{CNOT}}, {{\Bbb{T}}}_{7}=({\Bbb{I}}\otimes {R}_{x}(\pi /2)){\rm{CNOT}}\). The application of each of these operators is followed by the measurement in the σz bases and post-selection of the desired subspace. Thus, in each of these experiments, we measure all three qubits, and obtain eight diagonal elements pi,j,k = ci,j,k2. Finally, the corresponding populations of the two-qubit reduced density operator in the post-selected subspace with ancilla in state \({\left|0\right\rangle }_{{\rm{a}}}\) are given by

$${p}_{{\rm{j}},{\rm{k}}}^{(0)}=\frac{{p}_{0,{\rm{j}},{\rm{k}}}}{\mathop{\sum }\nolimits_{{\rm{j}},{\rm{k}} = 0}^{1}{p}_{0,{\rm{j}},{\rm{k}}}}.$$
(27)

Next, these populations are corrected for measurement errors, and the post-selected two-qubit density operators obtained further undergo convex optimization44,45 (see Supplementary Note 5 for further details).