Quantum simulation of parity-time symmetry breaking with a superconducting quantum processor

The observation of genuine quantum effects in systems governed by non-Hermitian Hamiltonians has been an outstanding challenge in the field. Here we simulate the evolution under such Hamiltonians in the quantum regime on a superconducting quantum processor by using a dilation procedure involving an ancillary qubit. We observe the parity-time ($\mathcal{PT}$)-symmetry breaking phase transition at the exceptional points, obtain the critical exponent, and show that this transition is associated with a loss of state distinguishability. In a two-qubit setting, we show that the entanglement can be modified by local operations.


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 shown [1] that non-Hermitian Hamiltonians endowed with parity and time (PT ) symmetry possess real positive eigenvalues and eigenvectors with positive norm. Experimental platforms where non-Hermitian Hamiltonians can be implemented comprise optical waveguides [2][3][4], polarized photons [5], nuclear spins [6,7], superconducting circuits [8,9], mechanical oscillators [10], nitrogen-vacancy centers in diamond [11,12], fiber-optics networks [13], and ultracold Fermi gases [14]. Open systems are a natural candidate for realizing these Hamiltonians, since non-Hermitian terms appear naturally as a consequence of energy being injected or lost [15].
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 losses [5], 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 noise [16]. Thus, in order to make progress one would need genuine realizations of non-Hermitian dynamics in the quantum regime [4,9,12], which maintain and allow the measurement of delicate quantum effects.
Here we show that the non-Hermitian dynamics can be simulated digitally [17] in a superconducting quantum processor by extending the Hilbert space with the use of * shruti.dogra@aalto.fi † Also at: Terra Quantum AG, St. Gallerstrasse 16A, 9400 Rorschach, Switzerland. ‡ sorin.paraoanu@aalto.fi an ancilla qubit and under the action of appropriatelydefined gates. To achieve this, we combine two techniques: a dilation method which 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 gates [18,19]. This combination enables us to observe and fully characterize the broken PT -symmetry transition and to settle decisively the relationship between non-Hermitian quantum mechanics and no-go theorems on state distinguishability and monotony of entanglement [20][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 IBM [24]. Although still imperfect, these devices have already enabled important results such as quantum error correction [25], fault-tolerant gates [26], proofs of violation of Mermin [27] and Leggett-Garg [28] inequalities, demonstrations of non-local parity measurements [29,30], simulations of paradigmatic models in open quantum systems [31], the creation of highly entangled graph states [32], the determination of molecular ground-state energies [33], the implementation of quantum witnesses [34], and quantumenhanced solutions to large systems of linear equations [35]. 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 Hamiltonian [36] with natural units = 1 where r is a real parameter and σ x and σ z are the Pauli matrices. The eigenvalues are ± √ 1 − r 2 , and the condition for non-Hermiticity is simply r = 0. The parity operator is P = σ x and the time-inversion operator is the complex conjugation operator T = . This Hamiltonian has an exceptional point at |r| = 1, where the two eigenvectors coalesce and the eigenvalues become paral-

Evolution
H a,q Initialization Readout lel. For |r| < 1 the eigenvalues are real, corresponding to distinct eigenvectors, and the Hamiltonian satisfies PTsymmetry [PT , H q ] = 0 (see Supplementary Note I); for |r| > 1, the eigenvalues become imaginary and the PTsymmetry 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 PT -symmetry [37]. 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 H a,q (t) acting on the total qubit-ancilla Hilbert space. The dynamics under H a,q (t) is determined by the Schrödinger equation whose solution is given by where |ψ(t) q is the solution of i d dt |ψ(t) q = H q |ψ(t) q .
In order to obtain a solution of the form Eq. (3), the ancilla and the qubit are initialized in the state |Ψ(0) a,q = (|0 a +η 0 |1 a )⊗|ψ(0) q , as shown in Fig. 1(a), by using a rotation R y (θ) on the ancilla qubit, where θ = 2 tan −1 η 0 and |ψ(0) q = |0 q . Thus, the dynamics of the system qubit in the subspace with the ancilla in state |0 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 a,q (t) = T exp[−i t 0 H a,q (τ )dτ ] is obtained numerically as described below.
Parity-time symmetry breaking in a single-qubit We implement the unitary operator U a,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. 1(b) the expected theoretical results for the ground-state population obtained under H q , i.e. generated by the non-unitary evolution operator exp(−iH q t) (see Methods). The breaking of PT -symmetry, as one crosses the exceptional point r = 1, is clearly visible. Indeed, for r < 1 the eigenvalues of H q are real and one observes Rabi oscillations. When |r| exceeds 1, the eigenvalues of H q become imaginary, the PT -symmetry is broken, and what one observes is the decay of the population. Fig. 1(c-e) present the results from the experimental realization of H q 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/ √ 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 V.
The details of the implementation are shown in Fig. 1(a). We start with the qubit and the ancilla both initialized in the state |0 , 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 U a,q (t) at any arbitrary time t is found by a numerical decomposition into single and twoqubit gates [18]. This decomposition U num (t) = U n . . . U 1 matches with the desired unitary operator U a,q (t) with a fidelity F U > 0.99, where the fidelity is defined as F U = 1 − ||U num (t) − U a,q (t)||/||U a,q (t)|| (also see Supplementary Note II). The quantum circuit which implements the decomposition of U num (t) (see inset of Fig. 1(a)) comprises a sequence of single qubit rotations U j q(a) , each of them having up to three degrees of freedom, and three two-qubit controlled-NOT gates [18,19]. The width of this circuit is 2 and the depth is 8. Specifically, the single-qubit gates are general rotations given by , where α, β, γ are the angles of rotations and the operators R y , R z correspond to the rotations generated by the Pauli operators σ y and σ z respectively. The operator U j q(a) (α, β, γ) 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. 1 After the U num (t) implementation, the post-selected subspace of our interest corresponds to the ancilla in state |0 a . At the end of the algorithm we measure the probabilities P kl of the qubit-ancilla state in the computational basis {|kl a,q ≡ |k a |l q } with k, l ∈ {0, 1}. Finally, the ground state population in the desired post-selected subspace of the system-qubit is given by, p 0 (t) = P 00 /(P 00 + P 01 ), which can be obtained directly from the experiments. These are shown with red dots in Fig. 1(c-e) and follow very closely the results for the population in the |0 q state of the qubit under the non-Hermitian Hamiltonian Eq. (1). The results demonstrate a high-fidelity simulation of the 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 orthogonal [20,22,23]. To observe this unusual feature of non-Hermitian dynamics we use the quantum circuit in Fig. 1(a) to evolve the system qubit, initialized respectively in the orthogonal states |0 q and |1 q . At various different instances of time, the state of the system qubit in the postselected sub-space with ancilla in state |0 a is obtained and the trace distance between the respective states is determined, where ρ diff (t) = ρ 1q (t) − ρ 2q (t) and ρ iq (t) = |ψ i (t) q ψ i (t)| q . For the given pair of initial states, the expected pattern for the variation of D with r and t is shown in Fig. 2(a). The characteristic recurrence time T R in the PT -symmetric phase and the decay time τ D in the broken-symmetry phase are plotted in Fig. 2(b) and compared to their analytical expressions (T R = π/ √ 1 − r 2 and τ D = 1/2 √ r 2 − 1 from Supplementary Equations (38) and (39)). Note that these times reflect the delicate balance between gain and loss, which is encoded in the structure of the Hamiltonian (see Supplementary Note IV). In Fig. 2(c)-(e) we show the experimentallyobtained 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 state [22] (see Supplementary Fig. 8). 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 D ∼ t −δ [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 ln D versus ln t starts to be approximately linear, with slope δ = 1.93 ± 0.08, see inset of Fig. 2(d). In this region, we obtain by fitting the experimental data δ = 1.75 ± 0.15 (dashed red line in the inset), a reasonably close value.
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 H q , for different values of r. Such scenarios have been studied theoretically [21,22], and it was shown that entanglement restoration and information recovery can happen in the 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 postselected subspace due to mere existence of a component of the total wavefunction outside this subspace [4,38,39].
To study this phenomenon, we consider a system consisting of two qubits q and q initialized in a maximally entangled Bell state, |Φ + q,q = (|00 q,q + |11 q,q )/ √ 2. One system qubit (say q) undergoes a non-Hermitian evolution by H q,q = H q ⊗ I q with the help of an ancillary qubit a, such that the total Hamiltonian including the dilation is H a,q,q = H a,q ⊗ I q leading to a unitary evolution, U a,q,q = U a,q ⊗ I q . Finally, the three-partite state of the system is measured and postselected subject to the state of the ancilla being |0 a .
The experimental implementation on the IBM quantum processor is carried out using three qubits, as shown in Fig. 3(a). 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 q , followed by σ z -measurements of all three qubits as shown in Fig. 3(a) -see Supplementary Note V for further details. At the end of each of these tomography measurements, the populations are obtained as p kl = P 0kl / 1 m,n=0 P 0mn , where k, l ∈ {0, 1}, from which the desired density operator of the system qubits ρ (0) q,q in the post-selected subspace is obtained. To study the entanglement dynamics we use the concurrence [40,41] as a measure of entanglement, given by C The change of concurrence with time is then observed for different values of r, as shown in Fig. 3 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. 3(b), while for r > 1 the Hamiltonian H q governing local evolution ceases to obey the PT symmetry, and the entanglement gradually decays with time, as shown in Fig. 3(d). For r = 1 we find the same theoretical asymptotic scaling as for distinguishability C (0) q,q ∼ t −δ , 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. 3(f) we present the corresponding theoretical curves for the time-variation of concurrence for a wider range of r parameters.  τ 123 Demonstration of an apparent violation of entanglement monotonicity. (a) Quantum circuit implemented on the superconducting quantum processor, where the qubit a serves as ancilla and q and q form a qubit-qubit bipartite system. In (b)-(d) 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 behaviour (continuous black lines). The inset in (c) presents the variation of the logarithm of concurrence versus 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 (b)-(d). 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 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 q under the action of the PT -symmetric non-Hermitian Hamiltonian. Specifically, at t = 0 we prepare the state cos(ϑ)|Φ − q,q − i sin(ϑ)|Ψ + q,q , where |Ψ ± q,q = (|01 q,q ± |10 q,q )/ √ 2 and |Φ ± q,q = (|00 q,q ± |11 q,q )/ √ 2 are the standard maximally entangled Bell states. The angle ϑ defines the concurrence | cos(2ϑ)| of this state. For the experimentshown in Fig. 3 (e) -we took ϑ = 59.185 0 , 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 |0 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. 3(e). In the first case, the entanglement increases up to approximately 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 tri-partite system (initialized in the state |0 ⊗ |Φ + ), where, as we have seen in Fig. 3(a), one of the system qubits q along with the ancillary qubit 'a' undergo the unitary evolution U a,q . The relevant correlations for the ensuing analysis are plotted in Fig.3(g) for r = 0.6. We define the single-qubit reduced states by tracing out the other qubits ρ i = Tr j,h [ρ i,j,h ], while the two-qubit reduced density operators are obtained by a single partial trace op-eration ρ i,j = Tr h [ρ i,j,h ], where the three qubits are labeled by i, j, h ∈ {a, q, q }. The concurrence associated with the state ρ i,j is denoted by 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 q are always in the permutation symmetric subspace of the two-qubit Hilbert space as one of the qubits evolves under H q (see Supplementary Note III). 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 ρ q remain maximally mixed all through the evolution, with a constant value of linear entropy s q = 1 − Tr(ρ 2 q ) = 0.5. Next, we observe that the concurrences C a,q and C a,q between the ancilla and the respective system qubits, i.e. a and q (or a and q ) always remain zero. This shows that the dynamics under U a,q,q 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 twoqubit reduced state of q, q and the ancilla. To quantify this tripartite correlation we use the three-tangle for pure states [41] τ a,q,q = C 2 a:q,q − C 2 a:q − C 2 a:q , or 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 C q:a ≡ C q,a and C q:q ≡ C q,q are the concurrences of the two-party reduced states ρ a,q and ρ q,q . The first term on the right hand side of Eq. (5) is the square of concurrence between the bi-partitions ρ q : ρ a,q , where one partition consists of the qubit q while the other partition is formed by the ancilla a and the system qubit q . For a pure three-qubit state ρ a,q,q , the quantity C q:a,q is effectively related to the mixedness of its bipartitions. More specifically, the square of concurrence between the partitions ρ q and ρ a,q is twice the linear entropy of the reduced density operator of either partition, given by 2(1 − Trρ 2 q ) or 2(1 − Trρ 2 a,q ). We now know from the simulated dynamics that the linear entropy s q(q ) = 1 − Trρ 2 q(q ) = 0.5, therefore at all times C 2 q:a,q = 1. Further, as shown in Fig. 3(g), there is no bipartite entanglement between the ancilla and the respective system qubits q (or q'), which implies that C q,a(a,q ) = 0. From Eq. (5), we obtain, Thus, the three-tangle among system qubits and ancilla and the concurrence between the system qubits are complementary to each other (see Supplementary Note III). By permuting the partitions in Eq. (5), it is easy to obtain τ a,q,q = 2s a , where s a = 1 − Tr(ρ 2 a ). The unitary U a,q,q , which induces a local non-Hermitian drive of qubit q in the post-selected subspace of the ancilla, is in fact a non-local operation on the system qubit q and the ancilla a. Under U a,q,q , 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 rdependent characteristic times, when postselected 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 Figure 6.

CONCLUSION
We have realized a quantum simulation of a singlequbit under a non-Hermitian Hamiltonian, observing the 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 η 2 + I from the dilation approach.
The simulation of phenomena governed by PTsymmetry at the single-quantum level open up several novel perspectives. Our scheme provides a systematic way of studying more complex non-Hermitian manyqubit 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 openconnected to a source of energy providing gain and reservoir for dumping this energy -they naturally will lead to new insights into quantum thermodynamics.

Simulating the non-Hermitian Hamiltonian in the dilated space
The single-qubit evolution under a general timedependent non-Hermitian Hamiltonian H q (t) is obtained in a certain subspace of an ancilla-qubit system undergoing a unitary evolution generated by H a,q (t). The Hamiltonian H a,q (t) in a four-dimensional Hilbert space can be obtained from H q by Naimark dilation [12]. Using this method, we can write the Hamiltonian H a,q (t) in the form: with where η(t) and M(t) are Hermitian operators; T and T are time-ordering and anti-time-ordering operators respectively, and I is the 2 × 2 identity operator. The Hamiltonian H a,q (t) can be obtained as follows. First, the equations Eq. (10,11) for η(t) and M (t) reflect the invariance of the norm of |Ψ a,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(d/dt)|Ψ a,q = H a,q (t)|Ψ a,q , and i(d/dt)|ψ(t) q = H q (t)|ψ(t) q together with Eq. (7) and Eq. (3) produce a linear system of equations in the unknown operators Λ(t) and Γ(t), 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,9). Note also that for Hermitian Hamiltonians H q the second term in Eq. (7), which is the qubit-ancilla interaction, becomes zero.

Initial conditions
To obtain an explicit form of H a,q (t), one should choose the operator M(t) at time t = 0 such that M(t)−I is positive for all t in the desired time interval. As a preliminary choice, we can take where m 0 > 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 labelled as µ 1 (t) and µ 2 (t), from where we numerically obtain µ min (t) = min{µ 1 (t), µ 2 (t)}. Interestingly, for r = 0, H q is Hermitian and with eigenvalues m 0 , which is the maximum value that µ min can assume. Therefore, for any arbitrary r and t, m 0 /µ min ≥ 1. Thus, at t = 0, M(t) is chosen to be, where f > 1, which ensures that M(t)−I remains positive for all t. From Eq. 10 we have The dynamics of the total ancilla-qubit system under H a,q (t) is obtained from the Schrödinger equation whose solution is given by where |ψ(t) q is the solution of i d dt |ψ(t) q = H q |ψ(t) q , and |ψ(t) q = η(t)|ψ(t) q . At t = 0 the state of the total system is which is a separable state of the ancilla |ψ(0) a and the system qubit |ψ(0) q . For preparing the initial state Eq. (19) the ancilla is taken in one of the eigenvectors of σ z , say |0 a . This is then subjected to a rotation by an angle θ around the y axis, R y (θ) = exp(−iθσ y /2), where, θ = 2 tan −1 η 0 . This leads to which is the initial state as defined by the protocol. On the other hand, the qubit q may be initialized in any arbitrary state |ψ(0) 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) |ψ(0) q = |0 and (ii) |ψ(0) q = |1 . 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 |ψ(0) q → |Φ + = 1 √ 2 (|00 + |11 ). For instance, in case i) we choose m 0 = 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, i.e. r = 1, in the same time range t ∈ [0, 8], we get η 0 = 16.1112 and θ = 3.0176 radians. Further, for r = 1.3, µ 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 product [23,36,42]. Here we make an explicit connection with this approach, showing that the Hermitian operator M (t) = η(t) 2 + 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 |Ψ(t) a,q , This norm has to be conserved during the time evolution. By taking the time-derivative of Eq. (21) and using Eq.
(2) we get This is the defining relation for the metric [23]. 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 |ψ(t) q in the Hilbert space we define the covector as q ψ(t)|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 H q studied in this work, the metric M (t) can be obtained analytically by employing the properties of 2 × 2 matrices (see also Supplementary Equation (30)). In the PTsymmetric phase we obtain an exact formula for the metric, 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) = M 0 . The result above Eq. (23) can be also obtained from the generic formula for the metric, as per [23], for parameters A = 0, B = −r/(1 − r 2 ), C = 1/(1 − r 2 ), 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 H a,q (t) in Eq. (7), where ψ(0) is the initial state at t = 0. This may also be written as 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 U a,q (t). This is done by observing the respective H a,q (t) evolutions of the complete set of two-qubit basis states. To find U a,q (t) we solve the Schrödinger equation numerically for different initial states, i d dt |ψ kl (t) a,q = H a,q (t)|ψ kl (t) a,q |ψ kl (0) = |kl a,q k, l = 0, 1.
where |00 a,q , |01 a,q , |10 a,q , and |11 a,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 H a,q (t) for a given time. Then, after solving this equation for |ψ 00 (t) , |ψ 01 (t) , |ψ 10 (t) , and |ψ 11 (t) we obtain the closed form of the unitary operator at an arbitrary time t, given by U a,q (t) = |ψ 00 (t) a,q 00| + |ψ 01 (t) a,q 01| +|ψ 10 (t) a,q 10| + |ψ 11 (t) a,q 11|. (26) For various different values of time, U a,q (t) is obtained, which is a general unitary operator in the four dimensional Hilbert space. Each U a,q (t) at a given time is then decomposed numerically in the form of single-qubit rotations and two-qubit controlled-NOT gates, as shown in Fig.(1e). This quantum circuit decomposition gives rise to U num (t), whose operation is very close to the theoretical U a,q (t). To characterize this, we calculate the error function err U (t) = ||U a,q (t) − U num (t)|| 2 /||U a,q (t)|| 2 , with the 2-norm defined by ||A|| 2 = √ λ max , where λ max is the largest eigenvalue of the matrix A * A. Here U num (t) is an unitary operator generated by the circuit in the inset of Fig. 1(a), where the parameters α, β, γ of U j q(a) (α, β, γ) are chosen to minimize the expression ||U a,q (t) − U num (t)|| 2 . Typically, we find err U (t) = ||U a,q (t) − U num (t)|| 2 /||U a,q (t)|| 2 to be of the order of 10 −4 , which demonstrates the high accuracy of our U a,q implementation. The accuracy with which our gate decomposition and the U a,q (t) operator match with each other is presented by an example data set in the Supplementary Table 1. U a,q (t) for arbitrary values of (r, t) and the corresponding U num (t) can be obtained from a GitHub code repository [43].

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 2 2n − 1 measurements corresponding to combinations of σ x , σ y , σ z and the identity matrix 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 |0 . 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 T 1 = I ⊗ I, T 2 = H ⊗ I, T 3 = R x (π/2) ⊗ I, T 4 = I ⊗ H, T 5 = I ⊗ R x (π/2), T 6 = (I ⊗ H)CNOT, T 7 = (I ⊗ R x (π/2))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 p i,j,k = |c i,j,k | 2 . Finally, the corresponding populations of the two-qubit reduced density operator in the postselected subspace with ancilla in state |0 a are given by Next, these populations are corrected for measurement errors, and the post-selected two-qubit density operators obtained further undergo convex optimization [44,45] (see Supplementary Note V for further details).

DATA AVAILABILITY
The data that support the findings of this study are available from authors upon reasonable request.

CODE AVAILABILITY
The codes used for the simulations can be found in the GitHub repository [43]. novation programme (grant agreement no. 862644, FET Open QUARTET), and from the Academy of Finland through project no. 328193 and through the "Finnish Center of Excellence in Quantum Technology QTF" project no. 312296. In addition, we would like to thank the Scientific Advisory Board for Defence (Finland) and Saab. We acknowledge the use of IBM Quantum services for this work. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team.

AUTHOR CONTRIBUTIONS
SD and GSP initiated the project and obtained the key analytical results. AM designed the quantum circuits and performed the simulations on the IBM Q machine with input from SD. All authors discussed the results. SD and GSP wrote the manuscript, with contributions also from AM.

CORRESPONDING AUTHORS
Correspondence to Shruti Dogra or Gheorghe Sorin Paraoanu.

Supplementary Material
where r is a real parameter and σ x and σ z are the Pauli matrices. Since this is essentially a spin-1/2 system, the parity operator is P = σ x and the time-inversion operator is the complex conjugation operator T = . PT -symmetry requires that H q and PT should share a common set of eigenvectors, which may be easily verified by operating PT on the eigenvectors of H q . We have three cases: Case I: 0 ≤ r < 1. In this situation √ 1 − r 2 is real and we obtain the eigenvalues of H q as ± √ 1 − r 2 , with the corresponding eigenvectors Acting with the PT operator we get Thus, |ψ ± is also an eigenvector of PT , which implies that, in the given range of the parameter r, the system Hamiltonian H q is PT -symmetric. In particular, at r = 0 the Hamiltonian reduces to σ x and the eigenvectors are real -the time-reversal in the PT operator does not play any role, and the only remaining action is that of σ x (identical to the Hamiltonian). Case II: r > 1. In this situation √ 1 − r 2 is purely imaginary; to put in evidence this we can write it as i √ r 2 − 1. The eigenvectors corresponding to the eigenvalues ±i √ r 2 − 1 are and acting with the time-reversal operator will result now in a change of sign of the square root term. Specifically, In this case the PT -symmetry is broken, and the action of the PT operator interchanges the eigenvalues. Case III: r = 1. In this case the two eigenvalues coalesce into a single one, zero, with corresponding eigenvector We get This is the exceptional point (EP).

B. Bloch sphere representation of PT -symmetry breaking
A simple geometric representation of PT -symmetry breaking can be obtained by using the Bloch sphere, see Supplementary Fig. 4. The eigenvalues derived above are represented by red and blue dashed lines, for r ∈ [0, 2]. At r = 0 the system Hamiltonian H q becomes σ x and its eigenvectors are orthogonal, lying along the ±x-axes on the Bloch sphere (shown with blue and red colored dots). As r increases from 0 to 1, the eigenvectors move closer to each other in the xy-plane and they are no longer orthogonal. Eventually these dashed lines join together at the exceptional point r = 1, corresponding to the coalescence of the eigenvectors of H s . This merging of the eigenvectors is shown in Supplementary Fig. 4 as the green colored dot along the y-axis, with coordinates (0, 1, 0). Going further to r > 1, the pair of eigenvectors separates again, moving in the yz-plane of the Bloch sphere. For very large values of r, the diagonal terms of H s begin to dominate and the corresponding eigenvectors approach the eigenvectors of σ z . Thus, the eigenvalues approach the ±z-axes asymptotically, which is shown by the black dots at (0, 0, ±1).

C. Single qubit evolution under Hq
In Fig. 1 b) we have presented the evolution of the ground-state population under the action of for various values of r. This can be obtained numerically as where we note that due to the non-Hermiticity the preservation of the norm is not guaranteed. This evolution can also be solved analytically. Indeed, we have On the single-qubit basis the action is If the initial state is |0 , this leads to the following expressions for the populations of the {|0 , |1 } basis, p 0 (t) = |α 0 (t)| 2 |α 0 (t)| 2 + |β 0 (t)| 2 and p 1 (t) = |β 0 (t)| 2 |α 0 (t)| 2 + |β 0 (t)| 2 .  I. Ground state population and errors at r = 0.6. We present the population of the ground state p0(t) calculated from the evolution under Hq and the corresponding numerical values p num 0 (t) obtained by using the quantum gate decomposition of Ua,q(t). We also show the error function errU which characterizes the mismatch between the two evolution operators.
For various different values of time, U a,q (t) is obtained, which is a general unitary operator in a four-dimensional Hilbert space. Each U a,q (t) at a given time is then decomposed numerically in the form of single-qubit rotations and two-qubit controlled-NOT gates as shown in Fig.(1e) of the maintext. This quantum circuit decomposition gives rise to U num (t), whose operation is very close to the theoretical U a,q (t). The accuracy with which our gate decomposition and the U a,q (t) operator match with each other is presented here by an example data set. Supplementary Table I shows a comparison between the theoretically expected values obtained from evolution with H q and the numerically obtained values obtained by simulating the results using the quantum circuit decomposition of U a,q (t). Here p 0 (t) and p num 0 (t) are the corresponding ground-state populations of the qubit at various different values of time t. To further characterize this, we calculate the error function err U (t) = ||U a,q (t) − U num (t)|| 2 /||U a,q (t)|| 2 , with the 2-norm defined by ||A|| 2 = √ λ max , where λ max is the largest eigenvalue of the matrix A * A. Here U num (t) is the unitary operator generated by the circuit in the inset of Fig.1(a) of the main text, where parameters α, β, γ of U j q(a) (α, β, γ) are chosen to minimize the expression ||U a,q (t) − U num (t)|| 2 . It is clearly apparent from the table that the theoretically expected and numerically obtained values of the ground state populations at various different times match exactly up to at least three decimal places, which shows the high accuracy level of our U a,q implementation.

III. SUPPLEMENTARY NOTE 3: DYNAMICS OF ENTANGLEMENT
In this section we focus on the dynamics of entanglement in the two experimental setups (single qubit and ancilla and two qubits and ancilla).

A. One-qubit system
We observed non-Hermitian PT -symmetric (0 < r < 1) and non-PT -symmetric (r ≥ 1) dynamics of the system qubit as shown in the main text. Considering the system qubit and the ancilla undergoing the unitary dynamics generated by H a,q (t) in the dilated space, non-local correlations between the system and the ancilla are developed in time. By performing quantum tomography we obtain the concurrence between the system and the ancilla shown in Supplementary Fig. 5 with red dots. The theoretically expected dependence is shown with continuous black line. Interestingly, the concurrence oscillate in time for 0 < r < 1, which may be interpreted as the information retrieval of the system.

B. Two-qubit system
Permutation symmetry of the two-qubit state We show here that the permutation symmetry (interchange of qubits q ↔ q ) of the initial Bell state |Φ + = 1 √ 2 (|00 + |11 ) is preserved under time evolution. This is not a priori obvious, since H q acts locally only on qubit q. The analytical form of the two-qubit state at an arbitrary time t is given by where N (t) = 1 i=0 (|α i | 2 + |β i | 2 ) is the normalization constant and coefficients α i (t) and β i (t) are defined in Eqs. (31). Since α 1 (t) = β 0 (t), we have, It is now easy to see that interchanging the qubits does not alter |Φ(t) .

Entanglement restoration
It has been shown in Fig. (3) of the main text that the concurrence between the system qubits varies with time when one of the system qubits undergoes a local evolution generated by a non-Hermitian Hamiltonian. This situation is experimentally realized in a dilated space of two system qubits and one ancilla, such that the dynamics of our interest lies in the post-selected subspace corresponding to ancilla in state |0 . The results from the simulation in Supplementary Fig. 6(a,b) depict the oscillations of the concurrence between the system qubits for r < 1 and decay for r > 1 respectively in the desired post-selected subspace. A similar pattern of variation in concurrence is also observed in the post-selected subspace with ancilla in state |1 , as shown in parts (c,d) of Supplementary Fig. 6. However, in the complete eight dimensional space of system qubits and ancilla, these three-qubits undergo a unitary dynamics and entanglement does not vary under local operations.
While we observe the variation in entanglement under local operations in one of the post-selected subspaces, other subspace of the same system also witnesses similar patterns for the variation of entanglement under local operations. As expected from standard Hermitian quantum mechanics, in the combined picture of system and ancilla the apparent violation of entanglement monotonicity no longer exists. A complete picture of correlations in this three-qubit system is presented in Supplementary Fig. 7.

IV. SUPPLEMENTARY NOTE 4: INFORMATION RETRIEVAL AND LOSS
Recalling the single-qubit dynamics from Section I C, the time-evolved state of the system ρ(t) with the non-Hermitian Hamiltonian H q is given by [22] where ρ(0) is the initial state of the system. Using the spectral decomposition, ρ(0) = mn ρ mn |ψ m ψ n |; in the eigenbasis of the Hamiltonian (H q ) we obtain In the case of a single qubit system, m, n ∈ {+, −} such that E ± are the eigenvalues of H q with corresponding eigenvectors |ψ ± and the difference in the eigenvalues, E m − E n = 2 |1 − r 2 |. For r = 0, the eigenvectors |ψ ± are not mutually orthogonal, which results in a time-dependent denominator (from now on, denoted by N (t)) in Supplementary Eq. 37. Supplementary Fig. 8 presents the variation of the normalization 1/N (t) with time t for different values of r. If r ∈ (0, 1), both the numerator and denominator of Supplementary Eq. 37 oscillate in time at the same rate. Thus there is information exchange between the system and the environment. The system retrieves the information from the environment in a time T R , the so-called recurrence time. Therefore for any arbitrary t, N (t + T R ) = N (t), leading to which gives For r > 1, the Hamiltonian H q does not exhibit PT symmetry. In this regime, the quantum system looses its information exponentially to the environment and never regains it back. The rate of loss of information may be quantified in terms of the decay time (τ D ) of the system. This may be interpreted as the time in which the norm 1/N (t) of the state ρ(t) diminishes by e. Therefore we have These properties of recurrence and decay for r ∈ (0, 1) and r > 1 respectively may be observed in the measurement of various quantities, as shown in the main text.
The values of recurrence time (T R ) and decay time (τ D ) are obtained experimentally by fitting the data using Supplementary Eq. (32). In the case r < 1, √ 1 − r 2 in Eqs. (31) is substituted by π/T R , while for r > 1, √ 1 − r 2 is replaced by −i/2τ D . These characteristic times are the properties of the Hamiltonian H q as they emerge solely from the quantum state evolution. Interestingly, we encounter these characteristic time scales in the study of various quantities. For instance, the measure of distinguishability, which is based on the trace distance between two arbitrary pairs of single-qubit states, oscillates at T R for the case of r < 1 and decays at τ D for r > 1. This may be readily seen in the simulation in Fig. 2(a) of the main text, and also in the experimental curves for r = 0.6 and r = 1.3 in Fig.  2(c,e) of the main text. The same T R and τ D emerge in the study of entanglement dynamics between two-qubits when one of these qubits is evolved under the non-Hermitian Hamiltonian H q . Indeed, these characteristic times appear in the experimental data from the main text, and they are well supported by the theoretical results of Fig. 3(f) in the main text and of Supplementary Fig. 7. An interesting scenario is presented in Supplementary Fig. 6, where in the subspaces post-selected by ancilla in state |0 and |1 respectively, oscillations in the measure of concurrence for 0 < r < 1 and decay for r > 1 exhibit different patterns. The underlying cause is the evolution under different non-Hermitian Hamiltonians, effective in the respective post-selected subspaces. In general, these oscillating and decaying patterns signify the information retrieval or information loss at r-dependent characteristic rates. These may be further interpreted to study the memory effects as in a non-Markovian dynamics and information decay in a Markovian dynamics respectively.

V. SUPPLEMENTARY NOTE 5: ADDITIONAL EXPERIMENTAL INFORMATION
In this work we used the ibmq ourense computer from IBM Q Experience. In experiments with 2 qubits we employed the zeroth and the first qubits. In experiments with 3 qubits we used the zeroth, the first and the second qubits. In Supplementary Table II  In this scheme, the first operator (T 1 ) provides the populations (4 − 1 = 3 variables). The operators T 2 and T 3 project the real and imaginary parts of the single-quantum coherences of the second qubit (4 variables) on the diagonal positions of the density operator, which is then read by σ z measurements. In a similar way, T 4 and T 5 provide the single-quantum coherences of the first qubit (4 variables). In order to obtain the terms corresponding to the zero quantum coherence (ρ exp 2,3 ) and double-quantum coherence (ρ exp 1,4 ) of the density operator, we use two-party correlations, which is a faster way to obtain these terms. Thus, the remaining four variables are obtained via operations T 6 and T 7 .
A. Experimentally obtained state of the system Multiple implementations (here 8192) of these seven tomography operations are performed on the system qubits (q and q'), while at the end of the sequence, all three qubits: ancilla (a) and system qubits are measured. These measurements actually provide the relative occupancy of the eight bases states of a three-qubit system at the end of each measurement, which we call here 'populations'. There is a possiblity of measurement errors associated with one or more of these populations, such as, a single-qubit state, which is in state |0 might be wrongly measured as |1 or vice versa. For the three-qubit case, a correction matrix [46,47] can be associated with each qubit, where f are the measurement fidelities of the bases states |0 and |1 respectively of the i th qubit (i ∈ {a, q, q }). Thus the correction matrix for a three-qubit system is given by, In order to compensate the effect of errors induced by the measurements, the corrected populations (p corr ) are obtained by, p corr = F −1 p m , where p m stand for the experimentally measured populations of the eight bases states. Measurement fidelities of each of these qubits are measured independent of each other by measuring their respective bases states (|0 and |1 ). We obtain the measurement fidelities: f = 0.98. However, we need only to examine the post-selected subspace of the total system with the ancilla in state |0 . 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. This is elaborated in the next sub-section. Using these intrinsic population values (p corr ), post-selected two-qubit density operators are obtained, which further undergo the convex optimization.

B. Postselection
Let us denote by ρ a,q,q the complete three-qubit density operator of the ancilla and the system qubits ρ a,q,q , ρ a,q,q = 1 i,j,k,l,m,n=0 c i,j,k c * l,m,n (|i a |j q |k q )( l| a m| q n| q ), where |i, j, k a,q,q (or |l, m, n a,q,q ) form the three-qubit computational basis, and c i,j,k (or c l,m,n ) are the corresponding coefficients of our three-qubit state, subject to the respective basis vectors. A two-qubit reduced density operator of (say) system qubits is obtained by the partial trace operation, which eliminates the information about the ancilla qubit. The two-qubit reduced density operator is thus given by ρ q,q = Tr a (ρ a,q,q ) = 1 j,k,m,n=0 c 0,j,k c * 0,m,n (|j q |k q )( m| q n| q ) + c 1,j,k c * 1,m,n (|j q |k q )( m| q n| q ) .
Most often, we post-select the subspace with ancilla qubit in state |0 a , and correspondingly, the two-qubit postselected reduced density operator is given by, c 0,j,k c * 0,m,n (|j q |k q )( m| q n| q ), where N (0) qq is responsible for the normalization of this post-selected subspace. To be more explicit, the two-qubit reduced subspace of our interest lies in a 4 × 4 block of the three-qubit state (ρ a,q,q ), which is shown colored in Supplementary Fig. 9. Interestingly, ρ q,q = ρ (0) q,q + ρ (1) q,q . Further, the single-qubit reduced density operators in respective subspaces is obtained by partial trace of one of the qubits. For instance, q,q ) and ρ q = Tr q (ρ q,q ).
As we described in the main text, a set of seven experiments are used to completely determine ρ (0) q,q . At the end of each of these operators (I ⊗ T ι , ι ∈ {1, 7}), elements at the diagonal positions are directly accessed via σ z measurements. Thus, in each of these experiments, we measure all three qubits, and obtain eight diagonal elements p i,j,k = |c i,j,k | 2 . Finally, the corresponding populations of the two-qubit reduced density operator in the postselected subspace, with ancilla in state |0 a is given by ρ a,q,q ′ = |0 a |j, k q,q ′ 0| a m, n| q,q ′ ρ (0) q,q ′ |0 a |j, k q,q ′ 1| a m, n| q,q ′ |1 a |j, k q,q ′ 0| a m, n| q,q ′ |1 a |j, k q,q ′ 1| a m, n| q,q ′ ρ (1) q,q ′ FIG. 9. Abstract respresentation of the three-qubit density operator, with the desired post-selected subspace, which is a 4 × 4 block highlighted with pink color.