Quantum simulation of the pseudo-Hermitian

,


I. INTRODUCTION
Standard quantum mechanics is formulated in terms of Hermitian Hamiltonians, which guarantee real values for energies as well as unitary time evolution.However, it is possible to extend the theory to non-Hermitian Hamiltonians in order to model open systems that significantly interact with their environment [1].The study of open, effectively non-Hermitian quantum systems already began in the early days of quantum theory with works such as Gamow's use of a complex energy for an α-decaying radioactive nucleus in 1928 [2] and Majorana's theory of α-particle absorption by a nucleus [3].Subsequent early applications of non-Hermitian quantum theory included a derivation of the nuclear dispersion formula [4] and a unified theory of nuclear reactions [5,6].
Non-Hermitian features are also a well-established part of fields such as superconductivity [7,8] and X-ray absorption spectroscopy [9], where an imaginary term iΓ is often added to energies to model the phenomenon of lifetime broadening; i.e. the broadening of measured excitation peaks caused by the finite lifetime of the excited state.In the context of superconductivity, Γ is referred to as the Dynes parameter.
The above historical examples mainly feature non-Hermitian elements (such as the addition of imaginary terms to energies) as ad-hoc modifications of otherwise Hermitian standard quantum mechanics.However, in 1998 Bender and Boettcher showed [10] that non-Hermitian Hamiltonians (NHHs) satisfying parityinversion and time-reversal (PT ) symmetry can exhibit real eigenvalues (i.e.energies) despite their non-Hermiticity, which led to an increased interest in the study of explicitly non-Hermitian Hamiltonians.In 2002, Mostafazadeh [11] showed that PT symmetry being associated with real eigenvalues is a consequence of PT -symmetric Hamiltonians being a subset of the larger class of pseudo-Hermitian Hamiltonians, and that the pseudo-Hermiticity of a Hamiltonian is a necessary but not sufficient condition for it having real eigenvalues.
Interest into PT -symmetric NHHs in particular has been driven by the fact that they can be experimen-tally realized as effective Hamiltonians in systems that involve a balanced gain-loss interaction with the environment [1].These experimental realizations have been especially prominent in the field of optics [12,13], where the Maxwell equations governing the dynamics of an electromagnetic wave in a waveguide under certain conditions can be written in a way that is identical with the Schrödinger equation describing a PT -symmetric system, thus providing an analog to PT -symmetric quantum mechanics [14][15][16][17][18]. PT -symmetric optics have been used to demonstrate exciting and exotic phenomena such as unidirectional invisibility [18][19][20][21] and single-mode lasing [17,[22][23][24][25][26].Alongside optics, PT -symmetric behaviour has also been realized on a large variety of other experimental platforms, such as the spins of ultracold 6 Li atoms [27], multilevel superconducting transmons [28], coupled LRC circuits [29], electromechanical resonators [30], trapped ions [31,32], 13 C-labeled chloroform dissolved in acetone-d6 [33], optically-induced atomic lattices [34], and even macroscopic coupled pendula [35].
In addition to direct realizations, NHHs can also be studied by digital quantum simulation [36,37] by employing the Naimark dilation method [38].This has been successfuly applied to platforms such as nitrogenvacancy centres in diamond [39][40][41] and superconducting qubits [42].However, all of these examples utilize a time-independent non-Hermitian Hamiltonian.
In this work we perform a quantum simulation of a time-dependent NHH.We extend the use of the Naimark dilation method to a time-dependent non-PT -symmetric NHH in order to simulate the pseudo-Hermitian Landau-Zener-Stückelberg-Majorana (LZSM) model [43] on a superconducting quantum processor.The LZSM model finds use in the description of transitions in a large variety of two-level systems, such as the valence/conduction bands in graphene [44], molecular states in ultracold Cs 2 [45], electron transfer between two As or P donors in a silicon nanowire transistor [46], and even the in-plane and out-of-plane resonator modes of a classical nanomechanical resonator [47].The pseudo-Hermitian extension of the LZSM model can be applied to e.g. the description of two electromagnetic modes travelling to opposite directions in a waveguide [48], and boosting a weak signal with a stronger one in the sum-frequency generation process [49].
This article is organized as follows: The next section (Section II) presents a theoretical overview of the pseudo-Hermitian LZSM model.Section III then outlines the methods employed to simulate a non-Hermitian Hamiltonian on a quantum processor, and Section IV presents the results of the application of this methodology to the simulation of the LZSM model described in Section II.Section V ends the main article with a brief conclusion.After the main article there are two Appendices: Appendix A presents the mathematical definitions of PT symmetry and pseudo-Hermiticity, and Appendix B includes calculations showing the application of these definitions to the pseudo-Hermitian LZSM model.

II. THE PSEUDO-HERMITIAN LANDAU-ZENER-ST ÜCKELBERG-MAJORANA MODEL
The pseudo-Hermitian Landau-Zener-Stückelberg-Majorana (LZSM) model [43] features a Hamiltonian of the form where σ i are the Pauli matrices and ε is varied linearly in time at a rate of v: All parameters in the Hamiltonian are assumed to be real, with v being positive.The parameter k controls the degree of non-Hermiticity; k = 1 reduces to the standard, Hermitian LZSM model [50].For k ̸ = 1 the Hamiltonian is not Hermitian or PT -symmetric, but is pseudo-Hermitian.Mathematical definitions of PT symmetry and pseudo-Hermiticity and details on how they apply to this Hamiltonian are included in Appendices A and B.
At Ω 0 = 0, the eigenstates of the system are the diabatic states [50] with the energy levels Adding a nonzero Ω 0 gives rise to the adiabatic energy states [50] with the corresponding eigenvalues where is the difference between the energy levels.Note that the form given in Eq. ( 6) is not normalized.Fig. 1 presents the diabatic and adiabatic energy levels as a function of time.For t < 0, |1⟩ is the ground state and |0⟩ the excited state of the system, but at t = 0 these levels cross, and for t > 0, |1⟩ is the excited state.At large |t| the diabatic and adiabatic levels approach each other, but at small |t| the adiabatic energy levels diverge from the diabatic ones and exhibit an avoided crossing near t = 0 [50].
The question asked and answered by the LZSM model is the following: When time is scanned from t = −∞ to t = ∞ across the avoided crossing region, what is the probability P that these following, equivalent statements happen?[50] • A system that starts in the diabatic state |0⟩ or |1⟩ at t = −∞ will remain in the same state at t = ∞.
• A system that starts in the adiabatic state |E + ⟩ or |E − ⟩ at t = −∞ will have transitioned to the other adiabatic state at t = ∞.
The answer is given by the LZSM formula [43,50] As we see, increasing the value of the parameter Ω 0 decreases P (assuming k > 0 for simplicity).Ω 0 determines the smallest distance ∆E = √ k|Ω 0 | (found at t = 0) between the adiabatic energy levels, with a large Ω 0 meaning the levels are far apart and transitions are unlikely.For the diabatic levels, as the off-diagonal element of Ĥ, Ω 0 represents the coupling between the states |0⟩ and |1⟩, and a large Ω 0 makes transitions between them more likely.
Counteracting the influence of Ω 0 is the rate v of the change in ε(t).If ε changes slowly (hv/Ω 2 0 ≪ 1), we have P ≈ 0. The system evolves adiabatically; i.e. it mostly follows the adiabatic state |E + ⟩ or |E − ⟩.With a fast rate of change (hv/Ω 2 0 ≫ 1) we have P ≈ 1, and the system undergoes diabatic evolution where it likely transitions from |E + ⟩ to |E − ⟩ or vice versa by following the diabatic state |0⟩ or |1⟩ [50].
In the Hermitian case, the rates P 0→1 and P 1→0 of transitions between diabatic states are simply equal to 1 − P , but in the pseudo-Hermitian case probability is not conserved and these instead take the forms [43] The time-evolution of a system under the pseudo-Hermitian LZSM Hamiltonian presented in Eq. ( 1) has been analytically solved [43], giving transition probabilities P i→f (t) between states not just in the scenario where time evolves from t 0 = −∞ to t = ∞, but also for any finite length of time evolution.The dynamical invariants of the system are given by the formulas which replace the conservation of total probability that is present in time-evolution under a Hermitian Hamiltonian [43,51].
For the k = −1 case, the transition probabilities and the conservation law in Eq. ( 12) have alse been presented by Malla et al. [52], who investigated transition probabilities, invariants, and integrability in the context of anti-Hermitian LZSM systems.

A. Time evolution operators
The time evolution operator Û (t, t 0 ) evolves an initial state |ψ(t 0 )⟩ to a final state |ψ(t)⟩: Û (t, t 0 ) is determined by the Hamiltonian of the system as a time-ordered exponential where T is the time-ordering operator which rearranges a product between multiple operators (there can be more than two) in a decreasing order of time arguments when read from left to right.
If the Hamiltonian Ĥ(t) is Hermitian, Û (t, t 0 ) is a unitary operator ( Û † Û = 1), which means it preserves the inner product between states: A particular consequence of this is the conservation of probability: if the state |ψ⟩ has been normalized so that at some time t = t 0 , then Eq. ( 18) will also hold for all other values of t as well.
B. Simulating time evolution

Hermitian Hamiltonians
A Hermitian Hamiltonian can be simulated on a quantum processor by computing (either numerically or analytically) the time evolution operator Û (t i , t 0 ) from Eq. ( 15) for a range of time values t i .Each of these operators can then be implemented as a quantum circuit by decomposing it into a product of operators corresponding to quantum gates.In practice we have used the QuantumCircuit.unitarymethod from the Qiskit library [53] to perform this decomposition.
An example of a two-qubit quantum circuit is presented in Fig. 2. The circuit consists of U(θ, ϕ, λ) quantum gates with controlled-NOT (CNOT) gates between them.A U(θ, ϕ, λ) gate performs a three-dimensional one-qubit rotation by the Euler angles θ, ϕ, λ, and has a matrix representation of When a quantum circuit that implements a time evolution operator is executed, its effect on the initial state |ψ(t 0 )⟩ of the qubits is equivalent to an operation by Û (t i , t 0 ), so that the state of the qubits at the end of  29)).This particular circuit was used for the final data point (at t = 20 h/v) for the k = 0.5 case in the results presented below in Fig. 3.The circuit features two qubits (q0 as the system qubit and q1 as the ancilla) as well as a classical register (labeled "c") where the state of the system is measured to at the end of the circuit.The purple boxes are U(θ, ϕ, λ) quantum gates that represent three-dimensional single-qubit rotations by the three given Euler angles.The blue symbols are controlled-NOT gates, with the smaller circle denoting the control qubit and the larger one the target qubit.
measurement each qubit is collapsed into one of the basis states |0⟩ and |1⟩.Using a large number of shots (i.e.individual executions of a quantum circuit, each ending in a measurement) per circuit thus gives a good estimate of the probabilities of different measured states (which consist of |0⟩ and |1⟩ in a 1-qubit circuit, |00⟩, |01⟩, |10⟩, and |11⟩ in a 2-qubit one etc.).When this process is repeated for all values t i in the selected range, data series describing the time evolution of the system in the form of probability functions P state (t) are obtained.

Non-Hermitian Hamiltonians
For a non-Hermitian Hamiltonian, the simulation procedure described above cannot be used without alterations.The time evolution operators Û (t i , t 0 ) produced by a non-Hermitian Hamiltonian can be nonunitary, and thus cannot be decomposed into a product of quantum gates, which only produce unitary transformations.However, it is possible to simulate a Hermitian system which contains a subsystem that evolves in a non-Hermitian fashion.More specifically, we use the mathematical technique of Naimark dilation to transform Ĥq , a non-Hermitian Hamiltonian for one qubit, into Ĥaq , a Hermitian Hamiltonian for a qubit and an ancilla, so that the ancilla = 0 subspace of the system evolves according to Ĥq .

C. Naimark dilation
The Hermitian Hamiltonian Ĥaq , obtained from Ĥq as a result of Naimark dilation, can be expressed as [39,42] where 1 is a 2 × 2 identity matrix and σy is the Pauli y-matrix.The operators Λ(t) and Γ(t) are defined as and and Ĥq (τ )dτ .
In the last equation, T is the time-anti-ordering operator which arranges time arguments in an increasing instead of decreasing order.
The initial value M (t 0 ) = M0 must be set so that M (t) − 1 is positive for all t in the simulated interval.This is accomplished as follows: M0 is initially defined as where m 0 > 1 is an arbitrary constant.Then the eigenvalues of M ′ (t) are obtained for all t in the desired interval (in practice a discrete lattice of t-values is used).The smallest eigenvalue (minimized over both the time axis and over the two different eigenvalues per time point) is labeled µ min , and M ′ 0 is replaced by the final value where f > 1 is another arbitrary constant.Finally, in order to produce the desired time evolution for the qubit in the ancilla = 0 subspace, the qubitancilla system must be initialized in the correct initial state.The qubit can be initialized in an arbitrary state; this state serves as the initial state for the non-Hermitian time-evolution.However, the ancilla is initialized as |0⟩ and then rotated around the y-axis by an angle θ: The angle θ is given by with η(t 0 ) = η 0 × 1; i.e. η 0 is the scalar equivalent of the matrix η(t 0 ).The rotation Ry (θ) is combined with the time evolution operator Û (t i , t 0 ) before decomposition into a quantum circuit, so that the action of the quantum circuit upon the uninitialized state |00⟩ aq can be mathematically expressed as where Înit(|ψ 0 ⟩) is an operator that initializes the qubit to its given initial state |ψ 0 ⟩ q .The lower indices a and q indicate whether a state or an operator applies to the ancilla, the qubit, or both.The probabilities of the states |0⟩ q (P0→0, blue) and |1⟩ q (P0→1, orange) as functions of time (in units of h/v), with |0⟩ q as the initial state, Ω0 = √ hv, and various values of k.The invariant kP0→0 + P0→1 (green), which should be equal to k, is also included.Solid lines represent theoretical results, crosses are results obtained assuming a perfectly accurate quantum computer, and dots are results computed with a real quantum computer.

A. Time evolution of probabilities
The dilation procedure described above was applied to simulating the non-Hermitian LZSM model, with the aim of replicating theoretically predicted results [43].Fig. 1 presents the results of the simulation of a qubit that starts in the state |0⟩ q at t 0 = −20 h/v and evolves until t = 20 h/v under the Hamiltonian from Eq. ( 1).The blue line and markers indicate the probability P 0→0 (t) of the qubit remaining in the initial state |0⟩ q at a given time value, while the orange line and markers indicate the probability P 0→1 (t) of the qubit having made the transition to |1⟩ q .
The solid lines represent a numerically computed theoretical prediction for the time evolution of the system, while the crosses represent the quantum circuits sent to be executed by the quantum processor.Any difference between the solid lines and the crosses is an error caused by numerical inaccuracies in the Naimark dilation process and/or the decomposition of a time evolution operator into a quantum circuit.The dots represent the results returned from the quantum processor; any difference between the crosses and the dots is an error introduced by the quantum processor, and is not present if a simulator emulating a perfectly accurate quantum processor is used in the place of a real quantum processor.
The simulation was performed on the ibm lagos processor, accessed through the IBM Quantum Platform cloud service.Qubit q 0 of ibm lagos was used as the simulated non-Hermitian qubit, while qubit q 1 served as the ancilla.Each data point was calculated with 10000 shots to ensure that the statistical uncertainties in the results are negligible.The qubit starts in the state |0⟩ q , while the initial state of the ancilla (represented as a0) is given by Eq. ( 27).
Running a quantum circuit on the processor returned values indicating how many times the system was found in each of the basis states |00⟩ aq , |01⟩ aq , |10⟩ aq , and |11⟩ aq of the qubit-ancilla system.This data corresponded to the time evolution of the Hermitian system with the Hamiltonian Ĥaq , and these basis states were equivalent to tensor products of the individual basis states of the qubit and the ancilla: Data corresponding to the time evolution of the non-Hermitian system with the Hamiltonian Ĥq was acquired through postselecting for the ancilla = 0 subspace.The results for the ancilla = 0 states |00⟩ aq and |01⟩ aq were used as those of the single-qubit subspace states |0⟩ q and |1⟩ q respectively, while the results for the ancilla = 1 states |10⟩ aq and |11⟩ aq were ignored.Results from before the postselection process are presented in Fig. 4.
Because the time evolution generated by Ĥq could be nonunitary, the probabilities P 0→0 (t) = | ⟨0|ψ(t)⟩ q | 2 and P 0→1 (t) = | ⟨1|ψ(t)⟩ q | 2 could be determined from the measured populations of the states |0⟩ q and |1⟩ q only up to a constant normalization factor N k , which could be different for each value of k.One possible way of determining N k was to note that in the theoretical model, P 0→0 (t 0 ) + P 0→1 (t 0 ) = 1 regardless of k, and set 1/N k equal to the sum of the measured populations of the states |0⟩ q and |1⟩ q at t 0 .However, with this method any errors in the measured populations at t 0 would affect the measured probabilities P 0→0 (t) and P 0→1 (t) at all values of t.This was particularly harmful in the k = −1 case, where the relative errors for the low-t data points were large due to the very small size of the ancilla = 0 subspace.To solve this problem, optimal values for N k were determined through performing a least-squares fit The probabilities of the states |0⟩ q (P1→0, blue) and |1⟩ q (P1→1, orange) as functions of time (in units of h/v), with |1⟩ q as the initial state, Ω0 = √ hv, and various values of k.The invariant kP1→0 + P1→1 (green), which should be equal to 1, is also included.
of the measured results (dots in the figures) to the postdecomposition predicted results (crosses in the figures).

B. Invariants
As can be seen in Fig. 3, in the non-Hermitian k ̸ = 1 cases the total probability P 0→0 + P 0→1 changes as a function of time instead of being restricted to 1.However, Eq. ( 12) holds for all k and t (up to the limits of the accuracy of the quantum computer), and replaces the total probability as the invariant of the non-Hermitian system as predicted [43,51].In order to show that the invariant from Eq. ( 13) holds as well, Fig. 3 was recreated with |1⟩ q as the initial state; the results are presented in Fig. 5.
Fig. 6 presents the dependence of the total probability P 0→0 +P 0→1 (left subplot) and the left side kP 0→0 +P 0→1 of the invariant relation given in Eq. ( 12) (right subplot) on the values of k and time.The three-dimensional surface represents numerically computed theoretical predictions, while the dots and crosses are based on the same data as the dots and crosses in Fig. 3.
C. Dependence on Ω0 Fig. 7 presents the dependence of the transition probabilities P 0→0 and P 0→1 on the parameter Ω 0 .The k = −1 case has been omitted, because for that case the theoretical results feature probabilities as high as 10 17 , which were too large for the numerical algorithms used in the simulation.The upper row features probabilities calculated with t 0 = −20 h/v and t = 20 h/v, so each data point there corresponds to the final data point of a plot like those presented in Fig. 3, but with varying values of Ω 0 .These data points have been measured and normalized in the manner explained in Section IV A. The . The probabilities of the states |0⟩ q (P0→0, blue) and |1⟩ q (P0→1, orange) as a function of Ω0 (in units of √ hv), with |0⟩ q as the initial state, t = 20 h/v, and various values of t0 and k.As Ω0 represents the coupling between the diabatic states, increasing its value increases the probability P0→1 of a transition.lower row presents similar results but with t 0 = 0 h/v and t = 20 h/v instead.

V. CONCLUSION
We have realized a quantum simulation of a single qubit under a time-dependent pseudo-Hermitian LZSM Hamiltonian, and used the resulting data to validate previous theoretical predictions.We have observed transition rates and the nonconservation of probability, the dependence of the results on time and the parameters k and Ω 0 , and the presence of the dynamical invariants of the system.The quantum simulation employs the Naimark dilation technique; by extending the Hilbert space with the use of an ancilla qubit and postselecting on the ancilla state, we were able to implement specific non-Hermitian dynamics.
This work demonstrates a useful application of current quantum computing technology, and shows that simulation of time-dependent non-Hermitian Hamiltonians is possible on a superconducting quantum processor.Future research will be focused on extending the methodology presented in this article to systems of more than one qubit, which will yield insights into multi-qubit phenomena such as quantum entanglement in a non-Hermitian context.The pseudo-Hermitian LZSM Hamiltonian from Eq. ( 1) is not PT -symmetric, as it does not match the generic forms given in Eqs.(33) and (34).An explicit calculation of the commutator [ Ĥ, P T ] yields [ Ĥ, P T ] = Ĥ P T − P T Ĥ

ACKNOWLEDGMENTS
where * is an operator that performs complex conjugation on everything to its right side.The result further confirms the lack of PT symmetry of the Hamiltonian.In the special case of k = 1, t = 0, the Hamiltonian is reduced to Ω 0 σx /2, which is PT -symmetric.

Rotated Hamiltonian
While the Hamiltonian from Eq. ( 1) is not PT -symmetric, it can be transformed into a PT -symmetric one by a π/2 rotation around the x-axis.This rotation transforms the Pauli matrices as σx → σx , σy → σz , and σz → −σ y , so that Eq. ( 1) is transformed into which matches the forms given in Eq. ( 33) and Eq. ( 34).The commutator is which explicitly shows the PT symmetry.

B Eigenvalues and eigenvectors 1 Original Hamiltonian
The following calculation shows that the expressions given in Eqs. ( 6) and ( 7) are indeed the eigenvectors and eigenvalues of the Hamiltonian Ĥ from Eq. (1): 2 Rotated Hamiltonian For the rotated Hamiltonian Ĥrot from Eq. ( 37), the eigenvalues are the same as given in Eq. ( 7), but the eigenvectors change to since The eigenvalues are real when kΩ 2 0 + ε 2 ≥ 0, and appear as imaginary complex conjugate pairs when kΩ 2 0 + ε 2 < 0. kΩ 2 0 + ε 2 = 0 is an exceptional point where the eigenvalues and eigenvectors of both Ĥ and Ĥrot become degenerate.Based on the realness of eigenvalues, kΩ 2 0 + ε 2 > 0 can be assumed to be the region of unbroken PT symmetry for Ĥrot , with kΩ 2 0 + ε 2 < 0 corresponding to broken PT symmetry.This can be confirmed by computing the effect of the PT operator on the eigenstates |E ± ⟩ rot of Ĥrot : Here we can transform the fraction in the first element of the vector into since In the region where ∆E ∈ R ⇔ kΩ 2 0 +ε 2 ≥ 0, the complex conjugation in the second element of the vector in Eq. ( 42) has no effect, and we have This shows that the eigenstates of Ĥrot are also eigenstates of P T , which is the condition for unbroken PT symmetry.

C Pseudo-Hermiticity
The Hamiltonian Ĥ from Eq. ( 1) can be easily seen to be η-pseudo-Hermitian with where a is a real constant, since η Ĥ η−1 = a k 0 0 1 Note that this choice of η is a particularly simple solution, but not a unique one.

Figure 1 .
Figure 1.A plot of the energies E (in units of Ω0) of the diabatic and adiabatic states |0⟩ (cyan), |1⟩ (yellow), |E+⟩ (red), and |E−⟩ (blue) as a function of time (in units of Ω0/v), with various values of the parameter k.A system starting in the state |E+⟩ = |0⟩ (or |E−⟩ = |1⟩) at t = −∞ and traversing the avoided crossing region towards t = ∞ will follow the diabatic state |0⟩ (|1⟩) and transition to |E−⟩ (|E+⟩) with a probability of P (black arrows).Solid lines indicate purely real energy values, while the dashed lines in the k = −1 case represent the imaginary parts of purely imaginary energy values.

Figure 3 .
Figure 3.The probabilities of the states |0⟩ q (P0→0, blue) and |1⟩ q (P0→1, orange) as functions of time (in units of h/v), with |0⟩ q as the initial state, Ω0 = √ hv, and various values of k.The invariant kP0→0 + P0→1 (green), which should be equal to k, is also included.Solid lines represent theoretical results, crosses are results obtained assuming a perfectly accurate quantum computer, and dots are results computed with a real quantum computer.

FK
acknowledges funding from the Magnus Ehrnrooth foundation and the Nokia Industrial Doctoral School in Quantum Technology.SD and GSP acknowledge support from the Academy of Finland (Research Council of Finland) through the Centre of Excellence program (project 352925).Appendix B MATHEMATICAL DETAILS ON THE PSEUDO-HERMITIAN LZSM HAMILTONIAN A PT symmetry 1 Original Hamiltonian Figure 2.An example of a quantum circuit obtained by decomposing a time evolution operator (including state initialization as described in Eq. (