# Dynamics simulation and numerical analysis of arbitrary time-dependent
$\mathcal{PT}$-symmetric system based on density operators

Xiaogang Li  State Key Laboratory of Low-Dimensional Quantum Physics and
Department of Physics, Tsinghua University, Beijing 100084, China  Chao Zheng
Department of Physics, College of Science, North China University of
Technology, Beijing 100144, China  Jiancun Gao  State Key Laboratory of Low-
Dimensional Quantum Physics and Department of Physics, Tsinghua University,
Beijing 100084, China  Frontier Science Center for Quantum Information,
Beijing 100084, China  Guilu Long<EMAIL_ADDRESS>State Key
Laboratory of Low-Dimensional Quantum Physics and Department of Physics,
Tsinghua University, Beijing 100084, China  Frontier Science Center for
Quantum Information, Beijing 100084, China  Beijing Academy of Quantum
Information Sciences, Beijing 100193, China  Beijing National Research Center
for Information Science and Technology and School of Information, Tsinghua
University, Beijing 100084, China

###### Abstract

$\mathcal{PT}$-symmetric system has attracted extensive attention in recent
years because of its unique properties and applications. How to simulate
$\mathcal{PT}$-symmetric system in traditional quantum mechanical system has
not only fundamental theoretical significance but also practical value. We
propose a dynamics simulation scheme of arbitrary time-dependent
$\mathcal{PT}$-symmetric system based on density operators, and the results
are compatible with previous methods based on pure-state vectors. Based on the
above, we are able to study the influence of quantum noises on the simulation
results with the technique of vectorization of density operators and
matrixization of superoperators (VDMS), and we show the depolarizing (Dep)
noise is the most fatal and should be avoided as much as possible. Meanwhile,
we also give a numerical analysis. We find that the problem of chronological
product usually has to be solved not only in the numerical calculation, but
also even in the experiment, because the dilated higher-dimensional
Hamiltonian is usually time-dependent. Through theoretical analysis and
numerical calculation, we find that on the premise of meeting the goal of
calculation accuracy and saving computing resources, the time step of
calculation and the cut-off term of Magnus series have to be carefully
balanced.

††preprint: APS/123-QED

##  I Introduction

That all physical observables, including Hamiltonians, must be Hermitian
operators has long been seen as one of the axioms in _conventional quantum
mechanics_ (CQM) [1], because Hermitian operators have real eigenspectrums as
we all know. However, Bender et al. found that some non-Hermitian
Hamiltonians, which are parity-time ($\mathcal{PT}$)-reversal symmetric, may
also have real eigenspectrums in 1998 [2], and then established the
_$\mathcal{PT}$ -symmetric quantum mechanics_ ($\mathcal{PT}$-QM) [3, 4, 5].
In last two decades, $\mathcal{PT}$-symmetry theory has developed rapidly [6,
7, 8, 9, 10, 11, 12, 13], and aroused wide attention [14, 15]. The phenomenon
related to $\mathcal{PT}$-symmetry exists widely, not only in classical
systems, such as optical systems[16, 17], microcavities [18] and circuits
[19], but also in quantum systems, such as strongly correlated many-body
systems [20], quantum critical spin chains [21], and ultracold atoms [22]. In
addition, $\mathcal{PT}$-symmetric system also shows practical value in
quantum sensors [23, 24, 25], which can use the sensitivity of
$\mathcal{PT}$-symmetric system near the exceptional points (EPs) to amplify
small signals. In particular, it is worth noting that, with the increasing
interest in $\mathcal{PT}$-symmetric systems, some new phenomena have emerged
[26, 27], which are impressive because they seem to conflict with theory of
conventional quantum mechanics or theory of relativity, such as the
instantaneous quantum brachistochrone problem [28, 29, 30, 31, 15, 32, 33, 34,
35], the discrimination of nonorthogonal quantum states [36, 37], and the
violation of no-signaling principle [38, 39, 40, 41, 42, 34, 43]. However,
some anomalies actually come from not knowing how to simulate
$\mathcal{PT}$-symmetric system in conventional quantum system, for instance,
if the discarded probability during the simulation of $\mathcal{PT}$-symmetric
system is considered, the no-signaling principle will still hold [44].
Therefore, finding a way to simulate $\mathcal{PT}$-symmetric system in
conventional quantum system has not only practical value, but also theoretical
significance.

At present, there are at least three technical routes to simulate
$\mathcal{PT}$-symmetric system, and all of them can be realized in
experiments [33, 25, 45]. The method of linear combination of unitaries (LCU)
[46], can be used to simulate various non-Hermitian $\mathcal{PT}$-symmetric
systems in discrete time, whether they are in unbroken or broken phase [46,
33, 47, 48], and even anti-$\mathcal{PT}$-symmetric systems [49]. The method
of weak measurement[50, 25], can also be used to simulate various time-
independent unbroken and broken $\mathcal{PT}$-symmetric systems in restricted
continuous time under the condition of weak interaction by a weak measurement.
The methods based on embedding [44, 51, 45], can be used to simulate the
dynamics of unbroken time-independent $\mathcal{PT}$-symmetric system (in
pure-states case [44] or mixed-states case [51]) in unrestricted continuous
time with only one qubit as an auxiliary system.

What deserves special attention is that in 2019, Wu et al. proposed a general
simulation scheme of dynamics of time-dependent (TD) arbitrary
$\mathcal{PT}$-symmetric system based on pure-state vectors through the
dilation method, and realized it with a single nitrogen-vacancy center in
diamond [45]. Specifically, their method is performed by dilating a general TD
$\mathcal{PT}$–symmetric Hamiltonian into a higher dimensional TD Hermitian
one with the help of an auxiliary qubit system, and evolving the state in the
dilated Hermitian system for a period of time, then performing a fixed
projection measurement on the auxiliary system, after that, the remained main
system is equivalent to going through the evolution process govern by the
$\mathcal{PT}$-symmetric Hamiltonian. However, the dilated Hamiltonian is
usually time-dependent, which means the system is actually an open quantum
system. From one point of view, as we all know, in an open quantum system, the
influence of quantum noises in the environment is usually inevitable so that
they have to be considered [52]. The time evolution of an open quantum system
interacting with memoryless environment can be described by the Lindblad
master equation [52, 53], which is usually based on the density operators
(matrices) rather than the state vectors. From another point of view, a pure
state will evolve to a mixed state under the quantum noises, while the theory
based on pure-state vectors can not conveniently deal with the question
related to mixed states, then the tool of pure-state vectors may be failed, in
this situation, the density operators will also be a better tool.

In this paper, we first generalize the outstanding work of Wu et al. based on
dilation method from the pure-state vectors case to the mixed-state density
operators case with the tool of density operators [45], and provide more
mathematical and physical completeness. It is worth emphasizing that it is not
a trivial process for generalizing the quantum state from pure-state vectors
to mixed-state density operators in the simulation of $\mathcal{PT}$-symmetric
system, and the difficulty mainly comes from the flexibility of characterizing
quantum states in $\mathcal{PT}$-symmetric system and the uncertainty of
mapping them to the high-dimensional quantum states in conventional quantum
system [51, 54]. Only based on this groundwork can we deal with problems in
open quantum systems, and we deal with these problems using the vectorization
of density operators and matrixization of superoperators (VDMS) technique.
Then we study the influence of quantum noises to the dynamics of time-
dependent arbitrary $\mathcal{PT}$-symmetric system, meanwhile, we also give a
numerical analysis. Through theoretical analysis and numerical calculation, we
find that on the premise of meeting the goal of calculation accuracy and
saving computing resources, the time step of calculation and the cut-off term
of Magnus series have to be carefully balanced. In the numerical analysis, we
also find that the time step $h$ of each numerical calculation step shall be
limited to its corresponding critical time $T_{c}$ of the convergence of
Magnus series, especially when the high-order terms of Magnus series are
considered. This phenomenon occurs because the dilated higher-dimensional
Hamiltonian is usually time-dependent, then the problem of chronological
product usually has to be dealt with, and the Magnus series may have to be
calculated [55, 56, 57], which may diverge when $t\rightarrow T_{c}$ so that
the error may be amplified after the Magnus series is truncated to a high-
order term in calculation. Meanwhile, the implemented duration of experimental
running is actually bounded by the critical time $T_{l}$ of the legitimacy of
dilation method. This phenomenon occurs because when $t\rightarrow T_{l}$, the
energy may diverge. In fact, the problem of chronological product may have to
be solved not only in the numerical calculation, but also even in the
experiment, because the dilated $\hat{H}_{AS}(t)$ has to be parameterized in
advance by numerically calculating the chronological product caused by
$H_{S}(t)$ needed to be dilated. In addition, when considering the influence
of quantum noises, we find that the depolarizing (Dep) noise (channel) is the
most fatal to the simulation of $\mathcal{PT}$-symmetric system among three
kinds of quantum noises we considered and should be avoided as much as
possible. It is worth noting that when the system considered is time-
independent and $\mathcal{PT}$-symmetry unbroken, the results of dynamics
simulation in this work are consistent with our previous results in Ref.[51],
and when the state considered is the pure state, the results of this work are
consistent with the theoretical results given in Ref.[45]. In summary, this
work provides a general theoretical framework based on density operators to
analytically and numerically analyze the dynamics of time-dependent arbitrary
$\mathcal{PT}$-symmetric system and the influence of quantum noises.

The rest of this paper is organized as follows. In Sec.II, we give some
necessary basic theories of $\mathcal{PT}$-symmetric system. In Sec.III, we
give a universal Hermitian dilation method of non-Hermitian Hamiltonians with
density operators, and based on that, we give a universal simulation scheme of
the dynamics of arbitrary TD $\mathcal{PT}$-symmetric system. To be able to
solve problems in open quantum system, we vectorize density operators and
matrixize the Liouvillian superoperators in Sec.IV. In addition, we specially
discuss the numerical calculation methods of time-dependent linear matrix
differential equations involved in this paper in Sec.V. In Sec.VI, we give an
example of two-dimensional $\mathcal{PT}$-symmetric system, and numerically
analyze its dynamics, meanwhile, we also consider the influence of three kinds
of quantum noises. In Sec.VII, we give conclusions and discussions. Further
more, we also make the Appendix A to show the details of the dilated
Hamiltonians, and the Appendix B to introduce the problem of chronological
product at the end of this paper.

##  II Theoretical preparations

Given a $n$-dimensional non-Hermitian $\mathcal{PT}$-symmetric Hamiltonian
$\mathcal{H}$, the parity operator $\mathcal{P}$ and time-reversal operator
$\mathcal{T}$, where $\mathcal{T}$ is an anti-linear operator, and $H$, $P$,
$T$ denote their matrix representation, respectively. They have the following
properties:

$\displaystyle P^{2}=I,T\overline{T}=$ $\displaystyle I,PT=T\overline{P},$
$\displaystyle PT\overline{H}=$ $\displaystyle HPT,$ (1)

where $\overline{[\cdot]}$ denotes complex conjugate of $[\cdot]$, and it
occurs because $\mathcal{T}$ is an anti-linear operator. If $H$ is similar to
a real diagonal matrix, $H$ will be $\mathcal{PT}$-symmetric unbroken,
otherwise, $H$ is called $PT$-symmetry broken if and only if it satisfies
either of these two conditions [7, 6, 58, 51]: (1) it cannot be diagonalized,
(2) it has complex eigenvalues that appears in complex conjugate pairs.

For a time-independent $\mathcal{PT}$-symmetry $H$, there is a time-
independent operator $\eta$ that satisfies:

$\eta H=H^{\dagger}\eta,$ (2)

where $\eta$ is called the metric operator of $H$, and it is a reversible
operator (usually Hermitian), when $H$ is $\mathcal{PT}$-symmetry unbroken, it
can be completely positive. The metric operator is usually not unique, for
instance, if $\eta$ is a metric operator of $H$, so is $r\eta$
($r\in\mathbb{R}$). The above Eq.(2) is also referred to as the pseudo-
Hermiticity relation, and $H$ is also referred to as pseudo-Hermitian
(Hamiltonian) [6, 59]. It is worth mentioning that in 2002, Mostafazadeh
pointed out that all the $\mathcal{PT}$-symmetric non-Hermitian Hamiltonians
belong to the class of pseudo-Hermitian Hamiltonians [6], recently, in 2020,
this conclusion was more strictly proved and strengthened by Ruili Zhang et
al., i.e., $\mathcal{PT}$-symmetry entails pseudo-Hermiticity regardless of
diagonalizability [13]. The theory of pseudo-Hermiticity provides more
convenience in dealing with questions related to mixed states, so we use the
tools usually used in the framework of pseudo-Hermiticity hereafter, such as
biorthonormal eigenbasis, metric operators , etc.

The elements of $\mathcal{PT}$-QM in the unbroken phase of $H$ can be
represented by biorthogonal basis of $H$:
{$|\chi_{n},a\rangle,|\phi_{n},a\rangle$}, and they have the properties as
follows [6]:

$\displaystyle\langle\chi_{m},a|\phi_{n},b\rangle=$
$\displaystyle\delta_{mn}\delta_{ab}$ (3a) $\displaystyle
H|\phi_{n},a\rangle=E_{n}|\phi_{n},a\rangle,\quad$ $\displaystyle
H^{\dagger}|\chi_{n},a\rangle=E_{n}|\chi_{n},a\rangle$ (3b)
$\displaystyle\sum\limits_{n}\sum\limits_{a=1}^{d_{n}}|\chi_{n},a\rangle\langle\phi_{n},a|=$
$\displaystyle\sum\limits_{n}\sum\limits_{a=1}^{d_{n}}|\phi_{n},a\rangle\langle\chi_{n},a|=I,$
(3c) $\displaystyle|\chi_{n},a\rangle=$ $\displaystyle\eta|\phi_{n},a\rangle,$
(3d) $\displaystyle\eta=\sum_{n}\sum_{a=1}^{d_{n}}$
$\displaystyle|\chi_{n},a\rangle\langle\chi_{n},a|,$ (3e)
$\displaystyle\eta^{-1}=\sum_{n}\sum_{a=1}^{d_{n}}$
$\displaystyle|\phi_{n},a\rangle\langle\phi_{n},a|,$ (3f)

where $d_{n}$ is the degree of degeneracy of the eigenvalue $E_{n}$ (in the
$\mathcal{PT}$-unbroken case, $E_{n}$ is real), and $a$ and $b$ are degeneracy
labels, and $|\phi_{k}\rangle$s ($|\chi_{k}\rangle$s) are usually not
orthogonal to each other. The possible real coefficients before $\eta$ have
been absorbed into the biorthogonal basis. It is worth noting that in the
extreme case, when $H$ becomes Hermitian, then the biorthogonal basis will
become orthogonal basis because of $H=H^{\dagger}$ and then
$\\{|\phi_{k}\rangle\\}=\\{|\chi_{k}\rangle\\}$. For convenience, we set
$d_{n}=1$ hereafter. If we recorded that
$\Phi=[|\phi_{1}\rangle,...|\phi_{i}\rangle,...|\phi_{n}\rangle]$,
$\Xi=[|\chi_{1}\rangle,...|\chi_{i}\rangle,...|\chi_{n}\rangle]$,
$E=$diag$(E_{1},...E_{i},...E_{n})$, then according to Eqs.(3) we will get
[44, 51]:

$\begin{split}&\Phi^{-1}H\Phi=E,\quad\Xi^{-1}H^{\dagger}\Xi=E\\\
&\eta=\Xi\Xi^{\dagger},\quad\eta^{-1}=\Phi\Phi^{\dagger}\\\
&\Xi=\eta\Phi,\quad\Xi^{\dagger}\Phi=I_{n}.\end{split}$ (4)

Through the positive Hermitian metric operator $\eta$, the representations of
a quantum observable $\mathcal{O}$ under the framework of CQM and the
framework of $\mathcal{PT}$-QM can be connected by a similarity
transformation, i.e., the Dyson map [60, 61]:

$\displaystyle O_{c}=\eta^{\frac{1}{2}}\cdot
O_{\mathcal{PT}}\cdot\eta^{-\frac{1}{2}},$ (5)

where $O_{c}$ is the observable in $\mathcal{PT}$-QM framework,
$O_{\mathcal{PT}}$ is the corresponding observable in CQM. Similarly, there is
a relation between the quantum state $\rho_{c}$ in CQM and the state
$\rho_{\mathcal{PT}}$ in $\mathcal{PT}$-QM:

$\rho_{c}=\sum_{mn}{\rho_{c}}_{mn}|m\rangle\langle
n|\Leftrightarrow\rho_{\mathcal{PT}}=\sum_{mn}{\rho_{c}}_{mn}|\phi_{m}\rangle\langle\chi_{n}|,$
(6)

where $\\{|n\rangle\\}$ is a mutual orthogonal basis in CQM;
${\rho_{c}}_{mn}$s are the matrix elements; $\rho_{c}$ and
$\rho_{\mathcal{PT}}$ are connected by the Dyson map given above in Eq.(5).
Therefore, there is a fixed relation between the quantum state (unnormalized)
$\rho_{S}$ in the framework of CQM and the quantum state $\rho_{\mathcal{PT}}$
in the framework of $\mathcal{PT}$-QM (for simplicity, we only use this
conclusion without more derivation, more details can be found in our previous
work in Ref.[51]):

$\displaystyle\rho_{S}$ $\displaystyle=\rho_{\mathcal{PT}}\cdot\eta^{-1}$
$\displaystyle=\sum_{mn}{\rho_{c}}_{mn}|\phi_{m}\rangle\langle\phi_{n}|,$ (7)

where $\rho_{S}$ can be normalized by $\rm{Tr}(\eta\rho_{S})=1$. Clarifying
the relation between the density operators in $\mathcal{PT}$-QM and the
density operators in CQM is an important step to extend the dynamics
simulation scheme of $\mathcal{PT}$-symmetric system from the pure-state
vectors case to the mixed-state density operators case [51]. From the above,
we know that this process is not trivial.

Now we introduce the concept of $\eta$-inner product [6, 51]:

$\displaystyle(|\psi_{1}\rangle,|\psi_{2}\rangle)_{\eta}\equiv\langle\psi_{1}\mid\psi_{2}\rangle_{\eta}:=\left\langle\psi_{1}|\eta|\psi_{2}\right\rangle$
$\displaystyle\forall\left|\psi_{1}\right\rangle,\left|\psi_{2}\right\rangle\in
L(\mathcal{H}),$ (8)

where $L(\mathcal{H})$ denotes Hilbert space, and $\eta$ is a reversible
Hermitian metric operator of this $\eta$-inner space, especially in the
unbroken phase of $\mathcal{PT}$-QM, it can be a positive operator [6].

Next we discuss the situation that non-Hermitian $\mathcal{PT}$-symmetric
system $H(t)$ is time-dependent. Considering two evolving states
$|\psi_{1}(t)\rangle$ and $|\psi_{2}(t)\rangle$, we assume that they satisfy
the Schrödinger-like equation:

$\displaystyle\frac{\mathrm{d}|\psi(t)\rangle}{\mathrm{d}t}=-iH(t)|\psi(t)\rangle,$
(9)

where we have set $\hbar=1$ here and after. According to the probability
conservation in the inner product space defined like in Eq.(II), we can obtain
that

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle\psi_{1}(t)|\psi_{2}(t)\rangle_{\eta(t)}$
$\displaystyle\equiv$
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\langle\psi_{1}(t)|\eta(t)|\psi_{2}(t)\rangle$
$\displaystyle=$
$\displaystyle\langle\psi_{1}(t)|-i\eta(t)H(t)+iH^{\dagger}(t)\eta(t)+\eta^{\prime}(t)|\psi_{2}(t)\rangle$
$\displaystyle=$ $\displaystyle 0,$ (10)

where we have recorded the differential operator
$\frac{\mathrm{d}}{\mathrm{d}t}$ as the symbol ”′”. Then we can get the
result:

$\displaystyle\eta^{\prime}(t)=i[\eta(t)H(t)-H^{\dagger}(t)\eta(t)].$ (11)

The above Eq.(11) is referred to as the time-dependent (TD) pseudo-Hermiticity
relation, $\eta(t)$ is the time-dependent (TD) metric operator in its
corresponding inner product space, i.e., $\eta(t)$-inner product space, which
leads to the probability conservation [60, 61]. Nothing that $\eta(t)$-inner
product space is time-dependent.

The solution of Eq.(11) can be obtained as:

$\displaystyle\eta(t)=\mathbb{T}e^{-i\int_{0}^{t}{H}^{\dagger}(\tau)d\tau}\eta(0)\mathbb{\overline{T}}e^{i\int_{0}^{t}H(\tau)d\tau},$
(12)

where $\mathbb{T}$ is time-ordering operator and $\mathbb{\overline{T}}$ is
the anti-time-ordering operator, moreover, $\eta(0)$ can be arbitrary
Hermitian operator. If we take $\eta(0)>0$, then there must exist a period of
time $T_{t}$ make $\eta(t)>0$ during $t\in[0,T_{t})$. It is worth noting that
when $H$ is time-independent, $\eta(t)$ may be time-independent so that
$\eta^{\prime}=0$ (for instance, when $H$ is $\mathcal{PT}$-symmetric
unbroken, $\eta(0)$ can be taken as the metric operator like the one in
Eq.(2)), and then the TD pseudo-Hermiticity relation given above in Eq.(11)
will be reduced to the pseudo-Hermiticity relation given in Eq.(2).

In addition, according to Eq.(II), if we set
$|\psi(t)\rangle_{1}=|\psi(t)\rangle_{2}=|\psi(t)\rangle$, and
$\rho_{S}(t)=|\psi(t)\rangle\langle\psi(t)|$, then we know:

$\displaystyle\frac{\mathrm{d}\mathrm{Tr}[\rho_{S}(t)\eta(t)]}{\mathrm{d}t}\equiv\frac{\mathrm{d}\mathrm{Tr}[\rho_{\mathcal{PT}}(t)]}{\mathrm{d}t}\equiv
0.$ (13)

Here $\rho_{\mathcal{PT}}(t)\equiv\rho_{S}(t)\eta(t)$ can be seen as a quantum
state in TD $\mathcal{PT}$-QM, and can be normalized by
$\mathrm{Tr}[\rho_{\mathcal{PT}}(t)]=\mathrm{Tr}[\rho_{S}(t)\eta(t)]\equiv 1$.
The form of $\rho_{\mathcal{PT}}(t)$ can be easily generalized, and can be
mapped to a quantum state $\rho_{c}(t)$ in (TD) CQM through a time-dependent
(TD) similarity transformation, i.e., the TD Dyson map similar to Eq.(5) [61]:

$\displaystyle O_{c}(t)=\eta^{\frac{1}{2}}(t)\cdot
O_{\mathcal{PT}}(t)\cdot\eta^{-\frac{1}{2}}(t),$ (14)

and $\rho_{c}(t)$ is similar to Eq.(6):

$\displaystyle\rho_{c}(t)=$
$\displaystyle\sum_{mn}{\rho_{c}}_{mn}(t)|m(t)\rangle\langle
n(t)|\Leftrightarrow$ $\displaystyle\rho_{\mathcal{PT}}(t)=$
$\displaystyle\sum_{mn}{\rho_{c}}_{mn}(t)|\phi_{m}(t)\rangle\langle\chi_{n}(t)|,$
$\displaystyle=$
$\displaystyle\sum_{mn}{\rho_{c}}_{mn}(t)|\phi_{m}(t)\rangle\langle\phi_{n}(t)|\cdot\eta(t),$
(15)

where $\rho_{c}(t)$ is a quantum state in CQM with a TD mutual orthogonal
basis $\\{|m(t)\rangle\\}$ in CQM, and $\\{|\phi(t)\rangle,|\chi(t)\rangle\\}$
is a TD biorthogonal basis in $\mathcal{PT}$-QM, and
$|\chi(t)\rangle=\eta(t)|\phi(t)\rangle$. Therefore, if we take $\eta(0)>0$
given in Eq.(11), similar to Eq.(6), we can obtain that

$\displaystyle\rho_{S}(t)$
$\displaystyle=\rho_{\mathcal{PT}}(t)\cdot\eta^{-1}(t)$
$\displaystyle=\sum_{mn}{\rho_{c}}_{mn}(t)|\phi_{m}(t)\rangle\langle\phi_{n}(t)|,$
(16)

which means on the promise that $\eta(t)>0$, we can always find an
unnormalized quantum state $\rho_{S}(t)$ in TD CQM related to a quantum state
$\rho_{\mathcal{PT}}(t)$ in TD $\mathcal{PT}$-QM (refer to Appendix B in
Ref.[51] for the proof that $\rho_{S}$ is actually an unnormalized state in
(TD) CQM). According to the relation between the unnormalized state of
$\rho_{S}(t)$ in TD CQM and quantum state $\rho_{\mathcal{PT}}(t)$ in TD
$\mathcal{PT}$-QM given above in Eq.(II), we know that once $\eta(t)$ is known
(given, or calculated), for the purpose of simulation, $\rho_{S}$ can be used
to represent $\rho_{\mathcal{PT}}$. It is worth noting that the two are
actually different in the physical sense, because they are not satisfied with
similarity transformation. Fortunately, for a simulation task, it is not
necessary to pursue the absolute equivalence of the two physical meanings, but
only to ensure that their form is appropriate and can be realized physically
[51, 31, 62]. The Eq.(II) is actually the prerequisite for the implementation
of the dilation method based on density operators we will discuss next.

##  III Universal Hermitian dilation method of non-Hermitian Hamiltonians
based on density operators

One of the methods to simulate the dynamics of $\mathcal{PT}$-symmetric system
is to find a dilated higher-dimensional Hermitian system (marked by ”$AS$”,
where ”$A$”, ”$S$”, ”$AS$” represents the auxiliary system, the main system
used to generate the dynamics of the non-Hermitian system, and the composite
system, respectively), which obey the von Neumann equation (it reduces to the
Schrödinger equation when restricted to pure-state vectors), to simulate the
dynamics of non-Hermitian system, which obey the von Neumann-like equation
(Schrödinger-like equation in pure-state vectors case). We assume the
evolution equation (the von Neumann-like equation) of unnormalized state
$\rho_{S}$ that has been mentioned in Eq.(II) is (hereafter, we set $\hbar=1$)
[54, 62, 27, 63]:

$\displaystyle\frac{\mathrm{d}{\rho_{S}(t)}}{\mathrm{d}t}$
$\displaystyle=-i[H_{S}(t),\rho_{S}(t)]_{{\dagger}}$
$\displaystyle\equiv-i[H_{S}(t)\rho_{S}(t)-\rho_{S}(t){H_{S}}^{\dagger}(t)],$
(17)

where $H_{S}(t)$ is non-Hermitian Hamiltonian in the system $S$ and can be
$\mathcal{PT}$-symmetric Hamiltonian, for generality, we assume it is time-
dependent, so the time-independent Hamiltonian can be regarded as a special
case of it. It is worth noting that because $H(t)$ is non-Hermitian, the
Eq.(III) cannot usually be realized directly in physics. $\rho_{S}(t)$ is the
unnormalized state in the system $S$ at time $t$s, its normalized form is

$\displaystyle\rho_{S_{\mathrm{N}}}(t)=\frac{\rho_{S}(t)}{\mathrm{Tr}[\rho_{S}(t)]},$
(18)

where $\rho_{S_{\mathrm{N}}}(t)$ is a legal quantum state (i.e., a positive-
semidefinite operator with unit trace) actually measured in the experiment,
and $P_{0}(t)\equiv\mathrm{Tr}[\rho_{S}(t)]$ represents the corresponding
measured probability of the state $\rho_{S_{\mathrm{N}}}(t)$ at the time
moment $t$. Because both $P_{0}(t)$ and $\rho_{S_{\mathrm{N}}}(t)$ can be
measured in an experiment, $\rho_{S}(t)$ can be used to represent
$\rho_{S_{\mathrm{N}}}(t)$, the probability $P_{0}(t)$ is actually absorbed in
the unnormalized density matrix $\rho_{S}(t)$ (see the Ref.[51] and its
Appendix B for details), and this practice meet the common usage[29, 30, 31,
45, 51]. Without misunderstanding, we assume that they are equivalent and are
strictly distinguished only when calculating probability. It should be noted
that since $H_{S}$ is non-Hermitian, $\mathrm{Tr}[\rho_{S}(t)]$ is usually not
constant, and may be greater than one so that loses its physical meaning.
Therefore, we set the constraint condition $\mathrm{Tr}[\rho_{S}(t)]\leqslant
1$, which can be satisfied by choosing the appropriate unnormalized initial
state $\rho_{S}(0)$ and duration time $T$ of the evolution.

However, because $H_{S}$ is non-Hermitian, the state evolution satisfying von
Neumann-like equation like Eq.(III) can not be realized directly in CQM.
Fortunately, we can map it into the following von Neumann equation. The von
Neumann equation related to the above von Neumann-like equation in the
Eq.(III) is:

$\displaystyle\frac{\mathrm{d}\rho_{AS}(t)}{\mathrm{d}t}=-i[\hat{H}_{AS}(t),\rho_{AS}(t)],$
(19)

where $\hat{H}_{AS}(t)$ is a Hermitian Hamiltonian in the system $AS$, so the
Eq.(19) can be realized directly in physics compared with Eq.(III), and
$\rho_{AS}(t)$ is the density operator in this system.

We assume that

$\displaystyle\rho_{AS}(t)$ $\displaystyle=\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho_{S}(t)\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\rho_{S}(t)&\rho_{S}(t)\cdot\xi^{\dagger}(t)\\\
\xi(t)\rho_{S}(t)&\xi(t)\cdot\rho_{S}(t)\cdot\xi^{\dagger}(t)\end{bmatrix},$
(20)

where $\begin{bmatrix}I_{S}\\\ \xi(t)\end{bmatrix}=|0\rangle_{A}\otimes
I_{S}+|1\rangle_{A}\otimes\xi(t)$, and we can see that the state
$\rho_{AS}(t)$ is usually an entanglement state.

Meanwhile, we assume

$\hat{H}_{AS}(t)=\begin{bmatrix}H_{1}(t)&H_{2}(t)\\\
H_{2}^{\dagger}(t)&H_{4}(t)\end{bmatrix},$ (21)

where $H_{1}(t),H_{2}(t),H_{4}(t)$ are all operators, and it is obvious that
$H_{1}^{\dagger}(t)=H_{1}(t),H_{4}^{\dagger}(t)=H_{4}(t)$, i.e., they are
Hermitian, because $\hat{H}_{AS}(t)$ is Hermitian.

Then according to the probability conservation principle, we know that

$\displaystyle\frac{\mathrm{d}{\mathrm{Tr}[\rho_{AS}(t)]}}{\mathrm{d}t}$
$\displaystyle=\frac{\mathrm{d}\mathrm{Tr}[\rho_{S}(t)+\xi(t)\rho_{S}(t)\xi^{\dagger}(t)]}{\mathrm{d}t}$
$\displaystyle=\frac{\mathrm{d}\mathrm{Tr}[(\xi^{\dagger}(t)\xi(t)+I_{S})\rho_{S}(t)]}{\mathrm{d}t}$
$\displaystyle\equiv\frac{\mathrm{d}\mathrm{Tr}[M(t)\rho_{S}(t)]}{\mathrm{d}t}$
$\displaystyle=\mathrm{Tr}[M^{\prime}(t)\rho_{S}(t)+M(t)\rho_{S}^{\prime}(t)]$
$\displaystyle=\mathrm{Tr}\\{M^{\prime}(t)\rho_{S}(t)-iM(t)[H_{S}(t)\rho_{S}(t)-\rho_{S}(t){H_{S}}^{\dagger}(t)]\\}$
$\displaystyle=\mathrm{Tr}\\{[M^{\prime}(t)-i(M(t)H_{S}(t)-{H_{S}}^{\dagger}(t)M(t))]\rho_{S}(t)\\}$
$\displaystyle\equiv 0,$ (22)

where we have made

$M(t)=\xi^{\dagger}(t)\xi(t)+I_{S}.$ (23)

It is obvious that $M(t)$ is Hermitian and $M(t)>1$, so M(t) is reversible.
For the convenience, we set

$\displaystyle\mathrm{Tr}[\rho_{AS}(t)]=\mathrm{Tr}[M(t)\rho_{S}(t)]\equiv 1.$
(24)

Then according to Eq.(III) we get the result

$M^{\prime}(t)=-i[{H_{S}}^{\dagger}(t)M(t)-M(t)H_{S}(t)].$ (25)

It is worth noting that the above Eq.(25) has the same form with Eq.(11), so
this relation can also be a TD pseudo-Hermiticity relation, which can replace
the (time-independent) pseudo-Hermiticity relation like
${H_{S}}^{\dagger}M=MH_{S}$ [28]. In general, $M^{\prime}(t)\neq 0$, so
${H_{S}}^{\dagger}(t)M(t)\neq M(t)H_{S}(t)$, which means $M(t)$ is not the
metric operator of $H_{S}(t)$, but only the TD metric operator of inner
product space (i.e., the $M(t)$-inner product space), which leads to
probability conservation. We find that $M(t)$ is actually the metric operator
of TD pseudo-Hermitian Hamiltonian $h_{S}(t)$ that will be introduced next.

Then multiplying left and right sides of the above Eq.(25) by $M^{-1}(t)$, we
can get

$iM^{-1}(t)M^{\prime}(t)M^{-1}(t)=M^{-1}(t){H_{S}}^{\dagger}(t)-H_{S}(t)M^{-1}(t).$
(26)

If we introduce an operator:

$K(t)=H_{S}(t)M^{-1}(t)+\frac{i}{2}M^{-1}(t)M^{\prime}(t)M^{-1}(t),$ (27)

then according to the above Eq.(26), we will find $K^{\dagger}(t)=K(t)$, which
means $K(t)$ is Hermitian. The introduction of the operator $K(t)$ will be
very beneficial to our following work. We then introduce a new quantity
$h_{S}(t)$:

$\displaystyle h_{S}(t)$ $\displaystyle\equiv K(t)M(t)$
$\displaystyle=H_{S}(t)+\frac{i}{2}M^{-1}(t)M^{\prime}(t)$
$\displaystyle=\frac{1}{2}H_{S}(t)+\frac{1}{2}M^{-1}(t)H_{S}^{\dagger}(t)M(t),$
(28)

where the Eq.(25) is used in the derivation. Obviously, $h_{S}(t)$ is non-
Hermitian. At the same time, it is easy to verify that

$\displaystyle M(t)h_{S}(t)=h_{S}^{\dagger}(t)M(t),$ (29)

which is exactly pseudo-Hermitian relation mentioned in Eq.(2), and that means
$h_{S}(t)$ is actually a pseudo-Hermitian quantity, and $M(t)$ is actually its
metric operator. As we have known in Eq.(23), $M(t)>I_{S}$, so according to
above Eq.(29), there is a similarity transformation make

$\displaystyle h_{\mathrm{phys}}(t)$ $\displaystyle\equiv
M^{\frac{1}{2}}(t)\cdot h_{S}(t)\cdot M^{-\frac{1}{2}}(t)$
$\displaystyle=M^{-\frac{1}{2}}(t)\cdot M(t)h_{S}(t)\cdot M^{-\frac{1}{2}}(t)$
$\displaystyle=h_{\mathrm{phys}}^{\dagger}(t)$
$\displaystyle=\frac{1}{2}[M^{\frac{1}{2}}(t)H_{S}(t)M^{-\frac{1}{2}}(t)\\!+\\!M^{-\frac{1}{2}}(t)H_{S}^{\dagger}(t)M^{\frac{1}{2}}(t)],$
(30)

which means $h_{\mathrm{phys}}(t)$ is a Hermitian operator with real
eigenspectrum and can be a legal observable in physics, so $h_{S}(t)$ is also
a legal observable with real eigenspectrum [64, 14, 15, 60], and can be a
physically permissible Hamiltonian. We reveal the physical meaning of $M(t)$,
and find the physical observables $h_{S}(t)$ and $h_{\mathrm{phys}}(t)$
related to it, so provide more physical completeness [45].

By solving the Eq.(25), we obtain that

$\displaystyle
M(t)=\mathbb{T}e^{-i\int_{0}^{t}{H_{S}}^{\dagger}(\tau)d\tau}M(0)\mathbb{\overline{T}}e^{i\int_{0}^{t}H_{S}(\tau)d\tau},$
(31)

where $M(0)$ can be any Hermitian operator satisfying the condition
$M(0)>I_{S}$. The above Eq.(31) is also obtained in Ref.[45]. It is worth
noting that if we take $M(0)=\eta(0)>1$ given in Eq.(12), then $M(t)$ will be
$\eta(t)$. This means that $M(t)$ is actually an TD metric operator like
$\eta(t)$ given in Eq.(12). In addition, when $H_{S}$ is given, some
eigenvalues of $M(t)$ may be decreases with time in some cases, such as in the
$\mathcal{PT}$-symmetry broken phase of $H_{S}$, so we can defined a critical
time $T_{l}$ of the legitimacy that make $M(t)$ legal, i.e., $M(t)>I_{S}$ when
$t\in[0,T_{l})$, while when $t=T_{l}$, at least one of eigenvalues of
$M(T_{l})$ become one. It should be noted that $T_{l}$ may be infinite in some
cases, such as in the $\mathcal{PT}$-symmetry unbroken phase of $H_{S}$. In
addition, in the cases that $T_{l}\neq\infty$, $T_{l}$ will depend on the
initial setting of $M(0)$, in general, $T_{l}$ will increase with the increase
of eigenvalues of $M(0)$, a big enough $T_{l}$ can always be obtained by
scaling $M(0)$ to meet the requirement of experiment or numerical calculation
tasks.

It is worth noting that when the system is time-independent and
$\mathcal{PT}$-symmetry unbroken, there must be a metric operator $\eta$ that
satisfies $\eta>1$, then $M(0)$ can be chosen as $M(0)=\eta$, after that

$\displaystyle M(t)$
$\displaystyle=e^{-i\int_{0}^{t}H^{\dagger}\mathrm{d}\tau}\eta
e^{i\int_{0}^{t}H\mathrm{d}\tau}$ $\displaystyle=\eta\cdot
e^{-i\int_{0}^{t}H\mathrm{d}\tau}\cdot e^{i\int_{0}^{t}H\mathrm{d}\tau}$
$\displaystyle=\eta,$ (32)

which means $M(t)$ will be the metric operator and is time-independent.

By the above Eq.(III), we have established a map between $\rho_{S}(t)$ and
$\rho_{AS}(t)$, now we try to establish the map between $H_{S}(t)$ and
$\hat{H}_{AS}(t)$, specifically, we need to find the solutions of
$H_{1}(t),H_{2}(t),H_{4}(t)$ in Eq.(21).

Substituting the Eqs.(III) and (III) into the Eq.(19), we obtain that

$\displaystyle\frac{\mathrm{d}{\rho_{AS}(t)}}{\mathrm{d}t}$
$\displaystyle=\begin{bmatrix}0\\\
\xi^{\prime}(t)\end{bmatrix}\cdot\rho_{S}(t)\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}+\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho^{\prime}_{0}\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}+\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho_{S}\cdot\begin{bmatrix}0&\xi^{\prime{\dagger}}(t)\end{bmatrix}$
$\displaystyle=-i\left\\{\begin{bmatrix}H_{S}(t)\\\
\xi(t)H_{S}(t)+i\xi^{\prime}(t)\end{bmatrix}\cdot\rho_{S}(t)\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}-\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho_{S}(t)\cdot\begin{bmatrix}{H_{S}}^{\dagger}(t)&{H_{S}}^{\dagger}(t)\xi^{\dagger}(t)-i\xi^{\prime{\dagger}}(t)\end{bmatrix}\right\\}$
$\displaystyle=-i[\hat{H}_{S}(t)\rho_{AS}(t)-\rho_{AS}(t)\hat{H}_{S}(t)]$
$\displaystyle=-i\left\\{\begin{bmatrix}H_{1}(t)&H_{2}(t)\\\
H_{2}^{\dagger}(t)&H_{4}(t)\end{bmatrix}\cdot\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho_{S}\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}-\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho_{S}\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}\cdot\begin{bmatrix}H_{1}(t)&H_{2}(t)\\\
H_{2}^{\dagger}(t)&H_{4}(t)\end{bmatrix}\right\\}$
$\displaystyle=-i\left\\{\begin{bmatrix}H_{1}(t)+H_{2}(t)\xi(t)\\\
H_{2}^{\dagger}(t)+H_{4}(t)\xi(t)\end{bmatrix}\cdot\rho_{S}\cdot\begin{bmatrix}I_{S}&\xi^{\dagger}(t)\end{bmatrix}-\begin{bmatrix}I_{S}\\\
\xi(t)\end{bmatrix}\cdot\rho_{S}(t)\cdot\begin{bmatrix}H_{1}(t)+\xi^{\dagger}(t)H_{2}^{\dagger}(t)&H_{2}(t)+\xi^{\dagger}(t)H_{4}\end{bmatrix}\right\\}.$
(33)

According to the above Eq.(III), we can get the following relation:

$\displaystyle H_{1}(t)+H_{2}(t)\xi(t)$ $\displaystyle=H_{S}(t)$ (34a)
$\displaystyle H_{2}(t)+H_{4}(t)\xi(t)$
$\displaystyle=i\xi^{\prime}(t)+\xi(t)H_{S}(t).$ (34b)

Observing the above equations, we know that the solutions of
$H_{1}(t),H_{2}(t),H_{4}(t)$ are not unique, and $H_{1}$ can be chosen as the
unique variable. Therefore, we can add a gauge in order to obtain unique
solutions. However, the value of $H_{1}$ is artificially assigned in Ref.[45],
so some unique properties of the resulting dilated $\hat{H}_{AS}$ may be
masked, we will see that in our Sec.VI, and especially in Fig.1 we will see
the eigenspectrum of $\hat{H}_{AS}$ actually has symmetric property, and
$\hat{H}_{AS}$ obtained in Ref.[45] is actually the result of the application
of a symmetric gauge. We then provide more mathematical completeness.

By observing the Eq.(21), we know that the space of $\hat{H}_{AS}(t)$ is
$2n$-dimensional, however, the space of $\rho_{S}$ defined in the Eq.(III) is
just $n$-dimensional, consequently, similar to the Eq.(III), we can assume
that

$\displaystyle\rho_{AS}^{\bot}(t)$
$\displaystyle=\begin{bmatrix}-\xi^{\dagger}(t)\\\
I_{S}\end{bmatrix}\cdot\rho_{S}(t)\cdot\begin{bmatrix}-\xi(t)&I_{S}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\xi^{\dagger}(t)\cdot\rho_{S}(t)\cdot\xi(t)&-\xi^{\dagger}(t)\cdot\rho_{S}(t)\\\
-\rho_{S}(t)\xi(t)&\rho_{S}(t)\end{bmatrix}.$ (35)

It is easy to check that $\mathrm{Tr}[\rho_{AS}(t)\rho^{\bot}(t)]\equiv 0$,
which means that the space of $\rho$ and $\rho^{\bot}$ are mutually
orthogonal, then we know that $\rho$ and $\rho^{\bot}$ are located in two
different orthogonal subspaces in the space where $\hat{H}_{AS}(t)$ is
located. Therefore, we can adopt the following symmetric gauge (there are also
some other valid gauges[59, 50, 61]):

$\frac{\mathrm{d}{\rho_{AS}^{\bot}(t)}}{\mathrm{d}t}=-i[\hat{H}_{AS}(t),\rho_{AS}^{\bot}(t)].$
(36)

With this gauge, according to the probability conservation principle again,
similar to Eq.(III), we can also derive another relation between $M(t)$ and
$\xi(t)$ as follows:

$M(t)=\xi(t)\xi^{\dagger}(t)+I_{S},$ (37)

where $M(t)$ has been given in Eq.(31). Comparing the Eq.(37) with the
Eq.(23), we get

$\xi(t)\xi^{\dagger}(t)=\xi^{\dagger}(t)\xi(t),$ (38)

which means $\xi(t)$ will be a normal operator under the symmetric gauge in
Eq.(36), therefore, for convenience, we take $\xi(t)$ Hermitian, and then

$\xi(t)=[M(t)-I_{S}]^{\frac{1}{2}}.$ (39)

At the same time, by adopting the similar method with Eq.(III), we can obtain
the following relation:

$\displaystyle-H_{1}(t)\xi^{\dagger}(t)+H_{2}(t)$
$\displaystyle=-i\xi^{\dagger}{{}^{\prime}}(t)-\xi^{\dagger}(t)H_{S}(t)$ (40a)
$\displaystyle-H_{2}^{\dagger}(t)\xi^{\dagger}(t)+H_{4}(t)$
$\displaystyle=H_{S}(t).$ (40b)

Then multiplying Eq.(34a) right by $\xi(t)$ and substituting it to Eq.(40a),
we can obtain

$\displaystyle H_{2}(t)$
$\displaystyle=[-i\xi^{\prime}(t)+H_{S}(t)\xi(t)-\xi(t)H_{S}(t)]M^{-1}(t)$
$\displaystyle=K(t)\xi(t)-\xi(t)K(t)-\frac{i}{2}[\xi^{\prime}(t)M^{-1}(t)+$
$\displaystyle\quad M^{-1}(t)\xi^{\prime}(t)].$ (41)

It is obvious that $H_{2}(t)$ is anti-Hermitian. The details of the derivation
is given in Appendix A.

Then in a similar way like above, multiplying Eq.(34b) right by $\xi(t)$ and
substituting it into Eq.(40b), we can obtain

$\displaystyle H_{4}(t)=$
$\displaystyle[i\xi^{\prime}(t)\xi(t)+\xi(t)H_{S}(t)\xi(t)+H_{S}(t)]M^{-1}(t)$
$\displaystyle=$ $\displaystyle
K(t)\\!+\\!\xi(t)K(t)\xi(t)\\!+\\!\frac{i}{2}[\xi^{\prime}(t)\xi(t)M^{-1}(t)\\!-$
$\displaystyle M^{-1}(t)\xi(t)\xi^{\prime}(t)],$ (42)

$H_{4}(t)$ is Hermitian, obviously. The details of the derivation is also
given in Appendix A.

Next, according to Eq.(III) and Eq.(34a), the result as follows will be
obtained:

$\displaystyle H_{1}(t)=$ $\displaystyle H_{S}(t)-H_{2}(t)\xi(t)$
$\displaystyle=$ $\displaystyle
H_{S}(t)-[-i\xi^{\prime}(t)+H_{S}(t)\xi(t)-\xi(t)H_{S}(t)]\cdot$
$\displaystyle M^{-1}(t)\cdot\xi(t)$ $\displaystyle=$
$\displaystyle[i\xi^{\prime}(t)\xi(t)+\xi(t)H_{S}(t)\xi(t)+H_{S}(t)]M^{-1}(t)$
$\displaystyle=$ $\displaystyle
K(t)+\xi(t)K(t)\xi(t)+\frac{i}{2}[\xi^{\prime}(t)\xi(t)M^{-1}-$ $\displaystyle
M^{-1}\xi(t)\xi^{\prime}(t)]$ $\displaystyle=$ $\displaystyle H_{4}(t).$ (43)

Consequently, according to Eq.(III) and Eq.(III) we can easily get

$\displaystyle H_{1}(t)+iH_{2}(t)$
$\displaystyle=[I_{S}-i\xi(t)]H_{S}(t)M^{-1}(t)[I_{S}+i\xi(t)]+\xi^{\prime}(t)M^{-1}(t)$
$\displaystyle=[I_{S}-i\xi(t)]K(t)[I_{S}+i\xi(t)]+\frac{1}{2}[\xi^{\prime}(t)M(t)^{-1}[I_{S}+i\xi(t)]+[I_{S}-i\xi(t)]M(t)^{-1}\xi^{\prime}(t)]$
$\displaystyle=[I_{S}-i\xi(t)]\left\\{K(t)+\frac{1}{2}[I_{S}-i\xi(t)]^{-1}[\xi^{\prime}(t)[I_{S}-i\xi(t)]^{-1}+[I_{S}+i\xi(t)]^{-1}\xi^{\prime}(t)][I_{S}+i\xi(t)]^{-1}\right\\}[I_{S}+i\xi(t)],$
(44a) $\displaystyle H_{1}(t)-iH_{2}(t)$
$\displaystyle=[I_{S}+i\xi(t)]H_{S}(t)M^{-1}(t)[I_{S}-i\xi(t)]-\xi^{\prime}(t)M^{-1}(t)$
$\displaystyle=[I_{S}+i\xi(t)]K(t)[I_{S}-i\xi(t)]-\frac{1}{2}[\xi^{\prime}(t)M(t)^{-1}[I_{S}-i\xi(t)]+[I_{S}+i\xi(t)]M(t)^{-1}\xi^{\prime}(t)]$
$\displaystyle=[I_{S}+i\xi(t)]\left\\{K(t)-\frac{1}{2}[I_{S}+i\xi(t)]^{-1}[\xi^{\prime}(t)[I_{S}+i\xi(t)]^{-1}+[I_{S}-i\xi(t)]^{-1}\xi^{\prime}(t)][I_{S}-i\xi(t)]^{-1}\right\\}[I_{S}-i\xi(t)],$
(44b)

where $M(t)=\xi^{2}(t)+I_{S}=[I_{S}\pm i\xi(t)]\cdot[I_{S}\mp i\xi(t)]$. Both
$H_{1}(t)+iH_{2}(t)$ and $H_{1}(t)-iH_{2}(t)$ are both Hermitian.

Finally, according to Eq.(21), we obtain the final form of $\hat{H}_{AS}(t)$:

$\displaystyle\hat{H}_{AS}(t)=$ $\displaystyle I_{S}\otimes
H_{1}(t)+i\sigma_{y}\otimes H_{2}(t)$ $\displaystyle=$
$\displaystyle|+_{y}\rangle\langle+_{y}|\otimes[H_{1}(t)+iH_{2}(t)]$
$\displaystyle+$
$\displaystyle|-_{y}\rangle\langle-_{y}|\otimes[H_{1}(t)-iH_{2}(t)],$ (45)

where
$|+_{y}\rangle=\frac{1}{\sqrt{2}}(|0\rangle_{A}+i|1\rangle_{A}),|-_{y}\rangle=\frac{1}{\sqrt{2}}(|0\rangle_{A}-i|1\rangle_{A})$
are the eigenstates of $\sigma_{y}$ corresponding to its eigenvalues of
$+1,-1$, respectively. The equation above has a similar form to the result of
pure-state case given in Ref.[45] (see Eqs.(19)-(21) in their Supplementary
Materials), while the differences are mainly caused by they actually use the
basis $\\{|-_{y}\rangle,-i|+_{y}\rangle\\}$, while we use the basis
$\\{|0\rangle,|1\rangle\\}$. According to Eq.(III), we can find $\hat{H}(t)$
is highly symmetric under the symmetric gauge of Eq.(36), and use only one
single qubit as its auxiliary system. We have to point out that, whether in
the experiments or numerical calculations, as long as we use a time-dependent
$\hat{H}_{AS}(t)$ given above in Eq.(III), it may be inevitable to solve the
problem of chronological product. The reason for the numerical calculation
situation is obvious, while in an experiment, $\hat{H}_{AS}(t)$ usually has to
be parameterized in advance by numerically calculating $M(t)$ given in
Eq.(31), which needs to numerically computing the chronological product caused
by $H_{S}(t)$, so it usually needs to deal with the problem of chronological
product unless any two moments of $\hat{H}_{AS}(t)$ are commute to each other
(such as the case of $H_{S}$ is time-independent and $\mathcal{PT}$-symmetry
unbroken, then $M(0)$ can be chosen as its metric operator $\eta$ just like
the case in Eq.(III) [45]. We will see the impact of the chronological product
on simulation accuracy more clearly in Fig.3 and Fig.4 in our example given in
Sec.VI.

For the convenience of expression, we define the ”$\circ$” operation as
$C\circ[\cdot]=C[\cdot]C^{\dagger}$. After that, looking back Eqs.(III) and
(19), their solution can be expressed as:

$\displaystyle\rho_{S}(t)$
$\displaystyle=\mathcal{U}_{S}\circ\rho_{S}(0)\equiv\mathbb{T}e^{-i\int^{t}_{0}H_{S}(\tau)\mathrm{d}\tau}\circ\rho_{S}(0),$
(46a) $\displaystyle\rho_{AS}(t)$
$\displaystyle=U_{AS}\circ\rho_{AS}(0)\equiv\mathbb{T}e^{-i\int^{t}_{0}\hat{H}_{AS}(\tau)\mathrm{d}\tau}\circ\rho_{AS}(0),$
(46b)

where $\hat{H}_{AS}(t)$ has been obtained in Eq.(III), and $\mathcal{U}_{S}$
is the non-unitary evolution operator related to the non-Hermitian Hamiltonian
$H_{S}$, while $U_{AS}$ is the unitary evolution operator related to the
Hermitian Hamiltonian $H_{AS}$. We will call the method of obtaining
$\rho_{AS}(t)$ through Eq.(46b) the dilation method, while call the method of
first obtaining $\rho_{S}(t)$ according to Eq.(46a) and then combining it into
$\rho_{AS}(t)$ by Eq.(III) the combination method. The difference between two
methods is that the delated Hamiltonian $\hat{H}_{AS}$ containing the term
$\xi^{\prime}(t)$ is avoid to be calculated for latter, which may increase the
error in numerical calculation, while only need to calculate $M(t)$ given in
Eq.(31) and $\xi(t)$ given in Eq.(39). The error of the numerical calculation
between the dilation method and the combination method can be defined as
follows:

$\displaystyle\Delta_{\rho}(t)=\|{\rho_{AS}}_{\mathrm{dilation}}(t)-{\rho_{AS}}_{\mathrm{combination}}(t)\|_{F},$
(47)

where ${\rho_{AS}}_{\mathrm{dilation}}(t)$ denotes the result calculated by
the dilation method, while ${\rho_{AS}}_{\mathrm{combination}}(t)$ denotes the
result calculated by the combination method, and the symbol ”$F$” denotes
Frobenius norm.

We define the measurement operator: $\Pi_{k}=|k\rangle_{A}\langle k|\otimes
I_{S},k\in\\{0,1\\}$, and the map $\mathcal{M}_{k}$:
$\mathcal{M}_{k}[\rho_{AS}]=\mathrm{Tr}_{A}[\Pi_{k}\circ\rho_{AS}]$.
Therefore, the Eq.(46b) can be mapped to the Eq.(46a) by the map
$\mathcal{M}_{0}$, experimentally, by performing a projection measurement
$\Pi_{0}\equiv|0\rangle_{A}\langle 0|$ on the auxiliary qubit:

$\displaystyle\mathcal{M}_{0}[\rho_{AS}(t)]=\rho_{S}(t)=\mathbb{T}e^{-i\int^{t}_{0}H_{S}(\tau)\mathrm{d}\tau}\circ\rho_{S}(0),$
(48)

which means we can simulate the non-unitary evolution (the dynamics) of non-
Hermitian system $H_{S}$ in the higher-dimensional system $\hat{H}_{AS}$. It
is worth noting that the map $\mathcal{M}_{0}$ is realized by a fixed
projection measurement (post-selection) $\Pi_{0}\equiv|0\rangle_{A}\langle
0|$, which is time-independent, so it will be easy to be realized in
experiment [45]. However, the success of this process $\mathcal{M}_{0}$ is
probabilistic, and the corresponding success probability $P_{0}(t)$ is

$\displaystyle P_{0}(t)=\mathrm{Tr}[\rho_{S}(t)].$ (49)

Obviously, in general, $P_{0}(t)<1$. After the measurement $\Pi_{0}$, the
state in main system $S$ will be ${\rho_{S}}_{\mathrm{N}}(t)$ in Eq.(18) with
success probability $P_{0}(t)$ given above. It is worth mentioning that the
success probability may be optimized by technical means, such as using the
local-operations-and-classical-communication (LOCC) protocol scheme proposed
in Ref.[51], which may be of significance for experiment.

In particular, in the situation that $H_{S}$ is $\mathcal{PT}$-unbroken and
time-independent, according to the result of Eq.(III), we can take the
operator $M(t)\equiv\eta$, where $\eta$ is a positive metric operator and
$\eta>1$, and then according to the Eqs.(27), (III), (III) and (III),
$K,H_{1},H_{2},H_{4}$ all will be time-independent, specifically, according to
Eqs.(4), they will become [44, 50, 51]:

$\displaystyle M(t)$ $\displaystyle\equiv\eta\Rightarrow\xi=(\eta-
I_{S})^{\frac{1}{2}},$ (50a) $\displaystyle K$ $\displaystyle\equiv
H_{S}\cdot\eta^{-1}=\Phi\cdot E_{S}\cdot\Phi^{\dagger},$ (50b) $\displaystyle
H_{1}$ $\displaystyle=H_{S}\eta^{-1}+\xi H_{S}\eta^{-1}\xi$
$\displaystyle=\Phi E_{S}\Phi^{\dagger}+\xi\cdot\Phi
E_{S}\Phi^{\dagger}\cdot\xi$ (50c) $\displaystyle H_{2}$
$\displaystyle=H_{S}\eta^{-1}\xi-\xi H_{S}\eta^{-1},$ $\displaystyle=\Phi
E_{S}\Phi^{\dagger}\cdot\xi-\xi\cdot\Phi E_{S}\Phi^{\dagger}$
$\displaystyle=-H_{2}^{\dagger},$ (50d) $\displaystyle H_{4}$
$\displaystyle=H_{S}\eta^{-1}+\xi H_{S}\eta^{-1}\xi$
$\displaystyle=\xi^{-1}(H_{S}\eta^{-1}+\eta
H_{S}-H_{S}-H_{S}^{\dagger})\xi^{-1}+H_{S}\eta^{-1}$ $\displaystyle=H_{1}.$
(50e)

They are the same as the results we obtained in Ref.[51]. Therefore, according
to Eq.(44), $H_{1}\pm iH_{2}$ will also become time-independent:

$\displaystyle H_{1}+iH_{2}$ $\displaystyle=(I_{S}-i\xi)\Phi\circ
E_{S}=\Phi_{+}\circ E_{S},$ (51a) $\displaystyle H_{1}-iH_{2}$
$\displaystyle=(I_{S}+i\xi)\Phi\circ E_{S}=\Phi_{-}\circ E_{S},$ (51b)

where $\Phi_{\pm}=(I_{S}\mp i\xi)$. At the same time, the dilated higher-
dimensional system $\hat{H}_{AS}$ will also become time-independent, and
according to the result of Eq.(III):

$\displaystyle\hat{H}_{AS}$ $\displaystyle=I_{A}\otimes
H_{1}+i\sigma_{y}\otimes H_{2}$
$\displaystyle=\\!I_{A}\\!\\!\otimes\\!(\\!H_{S}\eta^{-1}\\!\\!+\\!\\!\xi
H_{S}\eta^{-1}\xi)\\!+\\!i\sigma_{y}\\!\\!\otimes\\!(H_{S}\eta^{-1}\xi\\!-\\!\xi
H_{S}\eta^{-1}\\!)$ $\displaystyle=V_{AS}\circ[I_{A}\otimes E_{S}],$ (52)

where

$\begin{split}V_{AS}&=\frac{1}{\sqrt{2}}[(I_{A}+i\sigma_{x})\otimes
I_{S}-i(\sigma_{y}+\sigma_{z})\otimes\xi]\cdot I_{A}\otimes\Phi\\\
&=\frac{1}{\sqrt{2}}\begin{bmatrix}\Phi_{+}&i\Phi_{-}\\\
i\Phi_{+}&\Phi_{-}\end{bmatrix}\\\
&=\frac{1}{\sqrt{2}}\begin{bmatrix}(I_{S}-i\xi)\Phi&i(I_{S}+i\xi)\Phi\\\
i(I_{S}-i\xi)\Phi&(I_{S}+i\xi)\Phi\end{bmatrix}\end{split}$

is an unitary operator, while $V_{AS}$ is not unique. This result is also the
same as that in Ref.[51]. From the above Eq.(III), it can be found that the
degeneracy of higher-dimensional dilated system $\hat{H}_{AS}$ is twice that
of lower-dimensional $\mathcal{PT}$-symmetric system $H_{S}$ in this
situation.

##  IV Vectorization of density operators and matrixization of Liouvillian
superoperators in open quantum system

The evolution equation of an open quantum system with a Markovian
approximation (i.e., memoryless) can be expressed by Lindblad master equation
[53]:

$\frac{\mathrm{d}\rho_{AS}(t)}{\mathrm{d}t}=\mathcal{L}\rho_{AS}(t)=-i[\hat{H}_{AS}(t),\rho_{AS}(t)]+\sum_{\mu}\mathcal{D}[\Gamma_{\mu}]\rho_{AS}(t),$
(53)

where $\rho_{AS}(t)$ is the density operator of the system, $\Gamma_{\mu}$ is
the jump operator, $\mathcal{L}$ is the Liouvillian superoperator, and
$\mathcal{D}[\Gamma_{\mu}]$ is the dissipator related to the $\Gamma_{\mu}$,
which is used to describe the dissipation:

$\mathcal{D}\left[\Gamma_{\mu}\right]\rho_{AS}(t)=\Gamma_{\mu}\rho_{AS}(t)\Gamma_{\mu}^{\dagger}-\frac{\Gamma_{\mu}^{\dagger}\Gamma_{\mu}}{2}\rho_{AS}(t)-\rho_{AS}(t)\frac{\Gamma_{\mu}^{\dagger}\Gamma_{\mu}}{2}.$
(54)

The operator-sum representation is a convenient tool to describe the open
system, various models of decoherence and dissipation in open quantum systems
have been widely studied, such as amplitude damping (AD) channel model, phase
damping (PD) channel model, and depolarizing (Dep) channel. In some cases of
simple decoherence (such as the case of Hamiltonian $H=0$), these models have
a concise form in the framework of operator-sum representation, however, when
the Hamiltonian gets complicated just like the one in Eq.(III), the
application of the operator-sum representation will be indirect and
inconvenient. Therefore, the vectorization of density operators and
matrixization of superoperators (VDMS) technique may be a more convenient tool
(see the Appendix of Ref.[53] for details).

We adopt VDMS technique to carry out Kraus decomposition of density operators
in matrix basis $\\{|\alpha\rangle\langle\beta|,\alpha,\beta=1,2,\cdots,n\\}$
[65]. In this VDMS technique, a matrix $A$ can be mapped to a vector $\vec{A}$
by stacking all the rows of the matrix $A$ to a column in order:

$\displaystyle A=\begin{bmatrix}a_{1}&a_{2}\\\
a_{3}&a_{4}\end{bmatrix}\rightarrow\vec{A}=\begin{bmatrix}a_{1}\\\ a_{2}\\\
a_{3}\\\ a_{4}\end{bmatrix},$

In the similar way, the density operator $\rho$ will be mapped to the vector
$\vec{\rho}$, and the superoperators $\mathcal{A}[\cdot]=A[\cdot]B$ will be
mapped to a matrix, which can be recorded as
$\overline{\overline{\mathcal{A}}}=A\otimes B^{\mathrm{T}}$. so
$\overrightarrow{{\cal
A}[\cdot]}=\overrightarrow{A[\cdot]B}=\overline{\overline{\mathcal{A}}}\vec{[\cdot]}=A\otimes{B^{\rm{T}}}\vec{[\cdot]}$,
where the superscript ”$\mathrm{T}$” denotes the transpose operation..

Meanwhile, give a matrix

$\displaystyle B=\begin{bmatrix}b_{1}&b_{2}\\\ b_{3}&b_{4}\end{bmatrix},$

in the VDMS technique, the Hilbert-Schmidt inner product can be introduced as
follows:

$\displaystyle\langle
B|A\rangle=\begin{bmatrix}b_{1}^{*}&b_{2}^{*}&b_{3}^{*}&b_{4}^{*}\end{bmatrix}\begin{bmatrix}a_{1}\\\
a_{2}\\\ a_{3}\\\ a_{4}\end{bmatrix}=\mathrm{Tr}[B^{\dagger}A],$ (55)

where $|A\rangle,\langle B|$ are actually the $\vec{A},\vec{B}^{\dagger}$
mentioned above. From this point of view, it is easy to find that the vector
of $\rho_{AS}$ given in Eq.(III) and $\rho^{\bot}_{AS}$ given in Eq.(III) will
still be mutually orthogonal.

Therefore, the Lindbald superoperator $\mathcal{L}$ in Eq.(53) will be mapped
to

$\displaystyle\overline{\overline{\mathcal{L}}}(t)=$
$\displaystyle-i[\hat{H}_{AS}(t)\otimes
I_{AS}-I_{AS}\otimes\hat{H}^{\mathrm{T}}_{AS}(t)]$
$\displaystyle+\sum_{\mu}{[\Gamma_{\mu}\otimes({\Gamma^{\dagger}_{\mu}})^{\mathrm{T}}-\frac{\Gamma^{\dagger}_{\mu}\Gamma_{\mu}\otimes
I_{AS}}{2}-\frac{I_{AS}\otimes({\Gamma^{\dagger}_{\mu}\Gamma_{\mu}})^{\mathrm{T}}}{2}]},$
(56)

After that the Eq.(19) will be mapped to

$\frac{\mathrm{d}\vec{\rho}_{AS}(t)}{\mathrm{d}t}=\overline{\overline{\mathcal{L}}}(t)\vec{\rho}_{AS}(t),$
(57)

and then we can get the solution of the above equation:

$\vec{\rho}_{AS}(t)=\mathbb{T}e^{\int_{0}^{t}\overline{\overline{\mathcal{L}}}(\tau)\mathrm{d}\tau}\vec{\rho}_{AS}(0),$
(58)

where $\mathbb{T}$ is the time-ordering operator mentioned above in Eq.(12),
and $\vec{\rho}_{AS}(0)$ is the vector representation of the initial density
operator $\rho_{AS}(0)$. In general, the dilated $\hat{H}_{AS}(t)$ is time-
dependent, so $\overline{\overline{\mathcal{L}}}$ will also be time-dependent,
then the problem of chronological product may also have to be dealt with (see
more details in Appendix B and in our example given in Sec.VI).

##  V Numerical calculation of the relevant linear time-dependent matrix
differential equations

As we all know, when the system $H_{S}$ and $H_{AS}$ is time-independent, the
related calculations including Eq.(31), Eq.(III), Eqs.(46) and et.al. are
trivial. However, when they are time-dependent, the problem of chronological
product may have to be considered. According to Magnus’s theory [55, 57], the
solution of a linear matrix differential equation (include time-dependent
Schrödinger equation):

$\displaystyle Y^{\prime}(t)=A(t)Y(t),\quad Y\left(t_{0}\right)=Y_{0},$ (59)

can be expressed by the form of exponential matrix like:

$\displaystyle Y(t)=\exp\left(\Omega\left(t,t_{0}\right)\right)Y_{0},$ (60)

where $Y_{0}$ is the initial vector (state), and

$\displaystyle\Omega(t,t_{0})=\sum_{k=1}^{\infty}\Omega_{k}(t,t_{0})$ (61)

is the matrix series: so-called Magnus expansion or Magnus series, and
$\Omega_{n}(t,t_{0})$ is the $n$-th term (see more details in Appendix B).

However, for an arbitrary TD linear operator (matrix) $A(t)$, the above Magnus
series may diverge [57], and a sufficient but unnecessary condition for it to
converge in $t\in[t_{0},t_{f})$ is:

$\displaystyle\int_{t_{0}}^{t_{f}}\|A(s)\|_{2}\mathrm{d}s<\pi,$ (62)

where $\|A\|_{2}$ denotes 2-norm of $A$. The above integral converges
completely in the interval $t\in[t_{0},t_{0}+T_{c})$ (see more details in
Appendix B), where $T_{c}$ is defined as the critical time of convergence that
makes the above integral take $\pi$, and it only depends on the operator $A$
itself. This brings some restrictions to the time step $h$ in numerical
calculation, because when $h<T_{c}$, the results can be completely trusted,
while when $h>T_{c}$, the results becomes untrustworthy and the error may be
amplified after Magnus series is truncated to a high-order term $\Omega_{n}$.
The higher the order of items $\Omega_{n}$, the more obvious it is.

An alternative method to overcome the above finite convergence interval is to
divide the interval into $N$ segments so that the Magnus series in each
segment $[t_{k},t_{k+1}]$ meets the above convergence conditions Eq.(62),
i.e., $t_{k+1}-t_{k}<T_{c}$, then the Eq.(60) can be replaced by (see Eq.(240)
in Ref.[57]):

$\displaystyle
Y\left(t_{N}\right)=\prod_{k=0}^{N-1}\exp\left(\Omega\left(t_{k+1},t_{k}\right)\right)Y_{0},$
(63)

where $\Omega\left(t_{k+1},t_{k}\right)$ can be appropriately truncated to
$\Omega_{n}\left(t_{k+1},t_{k}\right)$ in practical use.

Suppose that a fixed step size $h$ is adopted in computation, according to the
analysis in Ref.[57]) (see around Eq.(66) and Eq.(242) for details), the error
of computations is $O(h^{3})$ when the computations up to first term of Magnus
series $\Omega_{1}$ are carried out, and $O(h^{5})$ when the computations up
to second term $\Omega_{2}$ are carried out. It is worth noting that, for the
purpose of numerical calculation, we usually have to consider the
computational complexity of the exponential matrix
$\exp\left(\Omega\left(t,t_{0}\right)\right)$, because the computation cost of
the exponential matrix is usually very expensive, especially when the matrix
becomes large, so in this case, we should try to avoid its frequent
computation. Therefore, on the premise of meeting the goal of calculation
accuracy and saving computing resources, the time step of calculation and the
cut-off term of Magnus series $\Omega_{n}$ have to be carefully balanced,
specifically, we can reduce the calculation times of the exponential matrix by
selecting a big enough $h$, and compensate for the loss of calculation
accuracy by calculating up to higher order terms of Magnus series in Eq.(63),
and we will see this compensation in the comparison between Fig.3 and Fig.4 in
next section.

##  VI An example: 2-dimensional $\mathcal{PT}$-symmetric system

In this section, we analyze an example: 2-dimensional $\mathcal{PT}$-symmetric
system:

$H_{S}=\begin{bmatrix}re^{i\theta}&s\\\
s&re^{-i\theta}\end{bmatrix},r,s\in\mathbb{R},\theta\in[-\pi/2,\pi/2],$ (64)

$\theta$ can be understood as the parameter representing the degree of non-
Hermiticity of the Hamiltonian $H_{S}$, and the degree of non-Hermiticity will
increase with $|\theta|$ (when $\theta=0$, $H_{S}$ will be Hermitian, when
$\theta=\pi/2$, $H_{S}$ will be anti-Hermitian, see Sec.V in Ref.[51] for
details). The eigenvalues of $H_{S}$ are
$E_{\pm}=r\cos{\theta}\pm\sqrt{s^{2}-r^{2}\sin^{2}{\theta}}$, and when
$s^{2}-r^{2}\sin^{2}{\theta}>0$, $H$ is $\mathcal{PT}$-symmetry unbroken,
otherwise, when $s^{2}-r^{2}\sin^{2}{\theta}<0$, $H$ is
$\mathcal{PT}$-symmetry broken, then the two eigenvalues become complex
conjugate.

The results of pure-state vectors case has been given in Ref.[45], in order to
facilitate comparison with the results in Ref.[45], we set $\theta=\pi/2,s=1$,
which makes $H_{S}$ have the same form as the one in Ref.[45], and under this
parameter configuration, the system is in the $\mathcal{PT}$-symmetry unbroken
phase, then the eigenvalues of $H_{S}$ will be $E_{\pm}=\pm\sqrt{1-r^{2}}$. At
the same time, we set $M(0)=5I_{S}$ given in Eq.(31). We set the initial
density operator of pure-state case as
$\rho_{S}(0)=\frac{1}{5}|0\rangle_{S}\langle 0|$, which correspondings to the
initial state $|0\rangle_{S}$ in Ref.[45], and the coefficient $1/5$ is
required by $\mathrm{Tr}[\rho_{AS}]\equiv 1$ according to Eq.(24). We set the
initial density operator of mixed-state case as
${\rho_{S}}_{\mathrm{mixed}}(0)=\frac{1}{30}\begin{pmatrix}4&1\\\
1&2\end{pmatrix}$, which can not be described by the pure-state vectors used
in Ref.[45].

###  VI.1 The effectiveness of the density operators tool and the
eigenspectrum of Hamiltonian before and after dilation

Figure 1: The eigenvalues of $\mathcal{PT}$-symmetric Hamiltonian $H_{S}$
(green curves with small circles or boxes, $E_{\pm}=\pm 0.8$) and the
eigenvalues of its dilated Hermitian Hamiltonian $\hat{H}_{AS}$ (curves with
other colors, $E_{1},E_{2},E_{3},E_{4}$). (a) $\mathcal{PT}$-symmetry unbroken
phase ($r=0.6$, $T_{l}=\infty$). (b) $\mathcal{PT}$-symmetry broken phase,
$E_{\pm}=\pm 0.98i$ ($r=1.4$, $T_{l}\approx 0.604$).

The eigenvalues of $\mathcal{PT}$-symmetric Hamiltonian $H_{S}$ ($E_{\pm}$)
and the eigenvalues of its dilated Hermitian Hamiltonian $\hat{H}_{AS}$
($E_{1},E_{2},E_{3},E_{4}$, arranged in descending order) are given in Fig.1.
Fig.1(a) is about the $\mathcal{PT}$-symmetric unbroken phase ($r=0.6$,
$E_{\pm}=\pm\sqrt{1-0.6^{2}}=\pm 0.8$), while Fig.1(b) is about the broken
phase ($r=1.4,E_{\pm}=\pm\sqrt{1-1.4^{2}}=\pm 0.98i$ will be conjugate complex
numbers, which is exactly what $\mathcal{PT}$-symmetry theory predicts [6]).
From Fig.1(a), we can see $E_{\pm}$ remain unchanged, while $E_{1}$ and
$E_{2}$ oscillate around $E_{+}$, $E_{3}$ and $E_{4}$ oscillate around $E_{-}$
and change periodically with time $t$ (in this case, the critical time $T_{l}$
of the legitimacy given below Eq.(62) under the setting $r=0.6$ is infinite,
i.e., $T_{l}=\infty$). Meanwhile, $E_{1}$ and $E_{4}$, $E_{2}$ and $E_{3}$ are
symmetric about $E=0$ (black dashed line). From Fig.1(b), we can see that
$E_{k}$s ($k=1,2,3,4$) are also symmetric about $E=0$ (black dashed line),
which is similar as the case of unbroken phase. In fact, the symmetry in the
case of unbroken and broken phase is caused by the symmetric gauge adopted in
Eq.(36). However, the periodicity oscillation of $E_{k}$s are destroyed, and
$E_{k}$s increase (decrease) monotonically with time $t$. Especially, when
$t\rightarrow T_{l}$ (the critical time $T_{l}$ of the legitimacy under the
setting $r=1.4$ is about 0.604), $E_{1}$ ($E_{4}$) will increases (decreases)
sharply to infinity, which is caused by one of eigenvalues of $M(t)$ given in
Eq.(31) tend to one. At this moment, the energy of system $\hat{H}_{AS}$ may
diverge, so can not be realized by an experiment. In a summary, the critical
time $T_{l}$ of the legitimacy actually bounds the duration of implemented
experimental running.

We also characterize the evolution using the renormalized population
${P_{\mathrm{N}}}_{0}(t)$ of state $\rho_{S}(0)$ in main system $S$ [45]. In
this situation, the renormalized population ${P_{\mathrm{N}}}_{0}(t)$ can be
obtained by

$\displaystyle{P_{\mathrm{N}}}_{0}(t)=\frac{\mathrm{Tr}[|0\rangle_{S}\langle
0|\cdot\rho_{S}(t)]}{\mathrm{Tr}[\rho_{S}(t)]}=\frac{\langle
0|\rho_{S}(t)|0\rangle_{S}}{P_{0}(t)},$ (65)

where $P_{0}(t)=\mathrm{Tr}[\rho_{S}(t)]$ has been given below Eq.(18). The
result based on the pure-state vectors has been given analytically in Ref.[45]
as follows:

$\displaystyle{P_{\mathrm{N}}}_{0}(t)=\begin{cases}\frac{\left|e^{t\sqrt{r^{2}-1}}\left(r+\sqrt{r^{2}-1}\right)-e^{-t\sqrt{r^{2}-1}}\left(r-\sqrt{r^{2}-1}\right)\right|^{2}}{\left|e^{t\sqrt{r^{2}-1}}\left(r+\sqrt{r^{2}-1}\right)-e^{-t\sqrt{r^{2}-1}}\left(r-\sqrt{r^{2}-1}\right)\right|^{2}+\left|ie^{-t\sqrt{r^{2}-1}}-ie^{t\sqrt{r^{2}-1}}\right|^{2}},&r\neq
1\\\ \frac{(t+1)^{2}}{(t+1)^{2}+t^{2}},&r=1.\end{cases}$ (66)

According to Eq.(46a) and above Eq.(65), we can calculate the renormalized
population ${P_{\mathrm{N}}}_{0}(t)$ under the initial density operator
$\rho_{S}(0)$. To verify the effectiveness of density operators method in
Sec.III, we take $r=0.6$ and draw the Fig.2. The green curve is drawn
according to the analytical results based on the pure-state vector method
given in the above Eq.(66), while the blue dotted curve is drawn according to
Eqs.(46a) and (65), and we can find their results are completely the same,
which means the density operator method are compatible with pure-state vector
method, just as we intuitively think. It should be pointed out that the
Hamiltonian $H_{S}$ in this example is time-independent, so the problem of
chronological product can be avoid in the calculation of $M(t)$ according to
Eq.(31).

Figure 2: Comparisons of the results for 2-dimensional
$\mathcal{PT}$-symmetric dynamics simulation (unbroken phase) between the
method based on pure-state vectors (green curve) and the method based on
density operator (blue dotted curve) ($\theta=\pi/2,s=1,r=0.6$). Both $X$ and
$Y$ axes are dimensionless. In this case, the system is in
$\mathcal{PT}$-symmetry unbroken phase.

###  VI.2 The influence of three kinds of quantum noises generated in main
system

In the experimental simulation of $\mathcal{PT}$-symmetric system dynamics,
the quantum state will inevitably be disturbed by quantum noises, especially
when the state is entangled just like the state in Eq.(III), the damage of
quantum noise to entanglement may be so fatal that has to be studied
carefully. In this section, we introduce three common quantum noise (channel)
models and using the VDMS technique mentioned in Sec.IV to characterize the
Lindblad master equation corresponding to them. We only consider the situation
where the noises act on the main system space (similarly, the situation where
the noises act on the auxiliary system can also be studied). Based on that, we
analysis their effects on the simulation of $\mathcal{PT}$-symmetric system
dynamics.

In AD channel and PD channel mentioned blow Eq.(54), there is only one jump
operator (for convenience, all parameters of decay rate have been set to
$\gamma$):
$\Gamma_{\mathrm{AD}}=\sqrt{\gamma}\sigma^{S}_{-}=\sqrt{\gamma}|0\rangle_{S}\langle
1|$, $\Gamma_{\mathrm{PD}}=\sqrt{\gamma}\sigma^{S}_{z}$; while in Dep channel,
there are three jump operators:
${\Gamma_{\mathrm{Dep}}}_{k}=\sqrt{\gamma}\sigma^{S}_{k}$, where $k=x,y,z$,
$\\{\sigma^{S}_{k}\\}$ are Pauli operators, and the subscript ”$S$” indicates
it acts on the main system, and the complete form of any operator $X$ with
different superscripts is $X^{S}=I_{A}\otimes X^{S},X^{A}=X^{A}\otimes I_{S}$.

Considering the system $\hat{H}_{AS}(t)$ used to simulate the
$\mathcal{PT}$-symmetric system dynamics in Eq.(III), and substituting jump
operators into Eq.(IV) we obtain that

$\displaystyle\overline{\overline{\mathcal{L}}}_{\mathrm{AD}}(t)$
$\displaystyle=-i[\hat{H}_{AS}(t)\otimes
I_{AS}-I_{AS}\otimes\hat{H}^{\mathrm{T}}_{AS}(t)]+\gamma[|0\rangle_{S}\langle
1|\otimes(|1\rangle_{S}\langle 0|)^{\mathrm{T}}-\frac{|1\rangle_{S}\langle
1|}{2}\otimes I_{S}-I_{S}\otimes\frac{(|1\rangle_{S}\langle
1|)^{\mathrm{T}}}{2}],$ (67a)
$\displaystyle\overline{\overline{\mathcal{L}}}_{\mathrm{PD}}(t)$
$\displaystyle=-i[\hat{H}_{AS}(t)\otimes
I_{AS}-I_{AS}\otimes\hat{H}^{\mathrm{T}}_{AS}(t)]+\gamma[\sigma^{S}_{z}\otimes({\sigma^{S}_{z}})^{\mathrm{T}}-I_{S}\otimes
I_{S}],$ (67b)
$\displaystyle\overline{\overline{\mathcal{L}}}_{\mathrm{Dep}}(t)$
$\displaystyle=-i[\hat{H}_{AS}(t)\otimes
I_{AS}-I_{AS}\otimes\hat{H}^{\mathrm{T}}_{AS}(t)]+\gamma\sum_{k\in\\{x,y,z\\}}[\sigma^{S}_{k}\otimes({\sigma}^{S}_{k})^{\mathrm{T}}-I_{S}\otimes
I_{S}],$ (67c)

where $\hat{H}_{AS}(t)$ is the dilated Hamiltonian in Eq.(III), $I_{AS}$ and
$I_{S}$ represent the identity operator in the composite system $AS$ and the
main system $S$ respectively, and they are the same when written in complete
form. We calculate Eqs.(67) and Eqs.(46) by Magnus series given in Eq.(60),
and only the first two terms of Magnus series, i.e., $\Omega_{1},\Omega_{2}$
are considered. Specifically, in the calculation, the operator $A(t)$ given in
Eq.(59) can be replaced by $-i\hat{H}_{AS}(t)$ obtained in Eq.(III), and
$\overline{\overline{\mathcal{L}}}_{\mathrm{AD}}(t),\overline{\overline{\mathcal{L}}}_{\mathrm{PD}}(t),\overline{\overline{\mathcal{L}}}_{\mathrm{Dep}}(t)$
given in Eqs.(67).

Based on the density operators, we are able to consider the effect of quantum
noise to the simulation of the dynamics of $\mathcal{PT}$-symmetric system
$H_{S}$. According to Eq.(48), Eqs.(67) and Eq.(65), we can calculate the
renormalized population ${P_{\mathrm{N}}}_{0}(t)$ under AD, PD, Dep channel
and no noise, respectively. Under the parameter settings of Fig.3, Fig.4 and
Fig.5, i.e., $\gamma=0.25,\theta=\pi/2,s=1,r=0.6$, the critical time of the
legitimacy mentioned below Eq.(31) $T_{l}=\infty$ because in the
$\mathcal{PT}$-symmetry unbroken phase of $H_{S}$, all eigenvalues of $H_{S}$
are real. In addition, the $\max_{t\in[0,8]}\|H_{AS}(t)\|_{2}\approx 1.08$,
$\max_{t\in[0,8]}\|\overline{\overline{\mathcal{L}}}_{AD}(t)\|_{2}\approx
2.19$,
$\max_{t\in[0,8]}\|\overline{\overline{\mathcal{L}}}_{PD}(t)\|_{2}\approx
2.34$,
$\max_{t\in[0,8]}\|\overline{\overline{\mathcal{L}}}_{Dep}(t)\|_{2}\approx
2.38$, so when we set time step $h=0.02$ (Fig.3(b)(d) and Fig.5) or $h=0.2$
(Fig.3(a)(c) and Fig.4), the convergence condition, i.e., Eq.(62), will always
be satisfied in every step because $h<T_{c}$.

The relations between the renormalized population ${P_{\mathrm{N}}}_{0}$ and
time $t$ in the case of quantum noise and no noise are given in Fig.3 and
Fig.4 respectively. At first, we focus on the cases of no noise, which are
related to the green curves (the combination method) and blue dotted curves
(the dilation method) in Fig.3. In Fig.3, the subfigures (a) and (b) are
related to the relation between ${P_{\mathrm{N}}}_{0}$ and $t$ with Magnus
series are calculated to the first term $\Omega_{1}$ in time steps $h=0.2$ and
$h=0.02$ respectively according to Eqs.(46), Eqs.(67) and Eq.(63), while in
Fig.4(a), Magnus series is calculated to the second term $\Omega_{2}$ with
time step $h=0.2$. In Fig.3(a-b) and Fig.4(c), the green lines represent the
results of no noise computed by the combination method (act as a theoretical
result), which only involves calculations involving $H_{S}(t)$, not
calculations involving $H_{AS}(t)$, while the blue doted lines represent the
results of no noise computed by the dilation method (act as a simulation
result), which involves calculations involving both $H_{S}(t)$ and
$H_{AS}(t)$. The Fig.3(c-d) and Fig.4(b) are errors between the green lines
and the blue lines. Compared Fig.3(c) with (d), we can find that when the time
step is increased from $h=0.2$ to $h=0.02$, the error is reduced by two orders
of magnitude, which is duce to the error is $O(h^{3})$ when the Magnus series
is computed to the first term $\Omega_{1}$ according to Eq.(63) [57].
Meanwhile, compared Fig.3(c) with Fig.4(b), by calculating Magnus series to
the second term $\Omega_{2}$, we can find the error in Fig.4(b) is also
reduced by two orders of magnitude although they are obtained with the same
time step $h=0.2$, which is due to the error is $O(h^{5})$ when the Magnus
series is computed to the second term $\Omega_{2}$ [57]. It is worth pointing
that although computing Magnus series to high-order terms
$\Omega_{n}(n\geqslant 2)$ may lead to higher accuracy, $\Omega_{n}$ are
usually difficult to be computed especially when $H_{AS}$ is high-dimensional
multivariable symbolic matrix (see Eq.(79) in Appendix B). On the contrary,
the advantage of only computing $\Omega_{1}$ is very convenient, however, the
improvement of accuracy has to be achieved by reducing the time step $h$,
which means that more exponential matrices have to be calculated, which is
usually computing-expensive especially when the matrix is high-dimensional.
Therefore, considering the computational cost of achieving a specific
accuracy, we usually have to make a compromise between the fewer computed
terms of Magnus series and the bigger step size.

Figure 3: The effects of quantum noises on the dynamics simulation of
$\mathcal{PT}$-symmetry unbroken system under the initial state of pure-state
density operator ($\gamma=0.25,\theta=\pi/2,s=1,r=0.6$). Both $X$ and $Y$ axes
are dimensionless. (a-b) Renormalized population ${P_{\mathrm{N}}}_{0}$ with
Magnus series are calculated to the first term $\Omega_{1}$ in time steps
$h=0.2$ and $h=0.02$ respectively. The linestyles and colors related to each
noise (including no noise) are represented differently. (c-d) The errors
between the green curves and the blue dotted curves corresponding to
subfigures (a-b) respectively.

Figure 4: The effects of quantum noises on the dynamics simulation of
$\mathcal{PT}$-symmetry unbroken system under the initial state of pure-state
density operator ($\gamma=0.25,\theta=\pi/2,s=1,r=0.6$). Both $X$ and $Y$ axes
are dimensionless. (a) Renormalized population ${P_{\mathrm{N}}}_{0}$ with
Magnus series are calculated to the second term $\Omega_{2}$ in time step
$h=0.2$. The linestyles and colors related to each noise (including no noise)
are represented differently. (b) The errors between the green curves and the
blue dotted curves corresponding to subfigure (a).

Now we focus on the effects of quantum noises on the renormalized population
${P_{\mathrm{N}}}_{0}$. From Fig.3(a-b) and Fig.4, we can find that the curves
related to Dep channel drops rapidly, and faster than the cases of AD, PD
channels, which is caused by the tendency of Dep channel that transforms any
quantum state into the maximal mixed state, where ${P_{\mathrm{N}}}_{0}$ will
keep $1/2$. The similar phenomenon can also be seen in mixed-state density
operator case in Fig.5, where we have set the initial density operator of the
mixed state as
${\rho_{S}}_{\mathrm{mixed}}(0)=\frac{1}{30}\begin{pmatrix}4&1\\\
1&2\end{pmatrix}$. Meanwhile, in the initial short time, the red curves
related to AD channel almost coincides with the cyan curves related to PD
channel, while after a long time, they separate from each other. The similar
phenomenon can also be seen in Fig.5. However, this phenomenon is accidental
because the initial pure-state density operator
$\rho_{S}(0)=\frac{1}{5}|0\rangle_{S}\langle 0|$ is just located on the
eigenstate (steady state, or the fixed point of superoperator of the noise
channel) of the dissipation terms of AD channel given in Eq.(67a) and PD
channel given in Eq.(67b) [53, 65], so their roles can be ignored in a short
time, and only the left terms containing $\hat{H}_{AS}(t)$ play the roles;
while in the long run, the evolution state $\rho_{S}(t)$ is far away from the
initial pure state $|0\rangle_{S}$, so it can be affected by the dissipation
terms where their roles cannot be ignored. Compared with the case of mixed
state, we can understand this phenomenon more clearly. From Fig.5 we can see
that in a short time, all the curves including noise and no noise cases rise,
which are the results driven by Hamiltonian $\hat{H}_{AS}(t)$. However, the
case of AD channel rises faster than all other curves, because the AD channel
has a tendency to change all states to the state $|0\rangle_{S}$, which will
contribute the renormalized population ${P_{\mathrm{N}}}_{0}$ (more strictly,
$|0\rangle_{S}$ is the fixed point (steady state) of AD channel [53, 65].

Figure 5: The effects of quantum noises on the dynamics simulation of
$\mathcal{PT}$-symmetry unbroken system ($\gamma=0.25,\theta=\pi/2,s=1,r=0.6$,
and set time step =0.02) under the initial state of mixed-state density
operator. Both $X$ and $Y$ axes are dimensionless. (a) Renormalized population
${P_{\mathrm{N}}}_{0}$ under the first term of Magnus series $\Omega_{1}$
according to Eq.(63). The linestyles (dotted curves) and colors related to
each noise (including no noise) are represented differently. (b) The error
between the green curve and the blue dotted curve corresponding to (a).

The effects of quantum noises on the dynamics simulation of
$\mathcal{PT}$-symmetry broken system under the initial state of pure-state
density operator $\rho_{S}=\frac{1}{5}|0\rangle_{S}\langle 0|$ are given in
Fig.6, and the renormalized population ${P_{\mathrm{N}}}_{0}$ is calculated to
the first term of Magnus series $\Omega_{1}$. In this figure, we set the
parameter as $\gamma=0.25,\theta=\pi/2,s=1,r=1.4$, then the legitimacy will be
lost after the critical time $T_{l}\approx 0.604$. Under this parameter
configuration, $\max_{t\in[0,0.6]}\|H_{AS}(t)\|_{2}\approx 14.20$,
$\max_{t\in[0,0,6]}\|\overline{\overline{\mathcal{L}}}_{AD}(t)\|_{2}\approx
28.40$,
$\max_{t\in[0,0.6]}\|\overline{\overline{\mathcal{L}}}_{PD}(t)\|_{2}\approx
28.40$,
$\max_{t\in[0,0.6]}\|\overline{\overline{\mathcal{L}}}_{Dep}(t)\|_{2}\approx
28.41$, so when we set time step $h=0.02$, the convergence condition, i.e.,
Eq.(62), will always be satisfied in every step because $h<T_{c}$ (including
noises case and no-noise case). From Fig.6(a) we can see that all curves
decrease monotonically because the eigenvalues of $H_{S}$ will be complex
numbers, which can be seen from Fig.1, so the evolution operator
$\mathbb{T}e^{-i\int_{0}^{t}H_{S}(\tau)d\tau}$ may cause the decay of model
$|0\rangle_{S}$. In addition, the blue curve and the green curve almost
completely coincide in the whole interval $t\in[0,T_{l})$, which can be seen
more clearly in the error diagram given in Fig.6(b), which shows high
accuracy.

Figure 6: The effects of quantum noises on the dynamics simulation of
$\mathcal{PT}$-symmetry broken system ($\gamma=0.25,\theta=\pi/2,s=1,r=1.4$,
and set time step $h=0.02$) under the initial state of pure-state density
operator. Both $X$ and $Y$ axes are dimensionless. (a) Renormalized population
${P_{\mathrm{N}}}_{0}$ under the first term $\Omega_{1}$ of Magnus series
according to Eq.(63). (b) The error between the green curve and the blue
dotted curve corresponding to (a).

##  VII Conclusions and discussions

In this work, we generalized the results of Wu et al. in Ref.[45], which are
based on the dilation method, from the pure-state vectors case to the mixed-
state case with the help of density operators, and provided a general
theoretical framework based on density operators to analytically and
numerically analyze the dynamics of TD arbitrary $\mathcal{PT}$-symmetric
system and the influence of quantum noises. We make conclusions from the
perspective of analytical analysis and numerical analysis, respectively.

At first, from the perspective of analytical analysis, more physical
completeness was provided. In the process of derivation, we discussed the
physical meaning of $M(t)$ ignored in Ref.[45]. Specifically, we proved that
$M(t)$ is not the metric operator of $\mathcal{PT}$-symmetric system $H_{S}$,
but the TD metric operator of $M(t)$-inner product space, which satisfies
probability conservation. Meanwhile, we also gave a quantity $h_{S}(t)$
related to $M(t)$, and proved that it actually is a physical observable,
because it can be mapped to the Hermitian quantity $H_{\mathrm{phys}}$, which
has a real eigenspectrum, through a TD similarity transformation, i.e., the TD
Dyson map. In addition, more mathematical completeness was also provided by
us. Specifically, in the derivation, we obtained the dilated Hamiltonian
$H_{AS}(t)$ by attaching a symmetric gauge instead of artificially assigning a
quantity to the free variable as in Ref.[45]. As a result, the hidden symmetry
of the eigenspectrum of the dilated Hamiltonian $H_{AS}(t)$ was able to be
revealed. It is worth noting that when the system considered is time-
independent and $\mathcal{PT}$-symmetry unbroken, the results of dynamics
simulation in this work are consistent with our previous results in Ref.[51],
and when the state considered is pure state, the results of this work are
consistent with the results given in Ref.[45]. Because the dilated system
$\hat{H}_{AS}(t)$ is actually located in an open quantum system, the influence
of environment will be inevitable, we introduced the tool of VDMS to solve the
Lindblad master equation under three kinds of quantum noises, then the
influence of quantum noises can be studied in the dynamics simulation of
$\mathcal{PT}$-symmetric system.

Then, from the perspective of numerical analysis, we pointed out that on the
premise of meeting the goal of calculation accuracy and saving computing
resources, the time step $h$ of calculation and the cut-off term of Magnus
series have to be carefully balanced. In addition, the time step $h$ of the
numerical calculation should be restricted by the critical time $T_{c}$ of
convergence of Magnus series in every step, because the dilated higher-
dimensional Hamiltonian $\hat{H}_{AS}(t)$ is usually time-dependent, the
problem of chronological product has to be solved, and the Magnus series has
to be calculated, which may diverge when $t>T_{c}$ so that the error may be
amplified after Magnus series is truncated to high-order terms in calculation.
Meanwhile, the implemented duration of experimental running is actually
bounded by the critical time $T_{l}$ of the legitimacy of dilation method.
This phenomenon occurs because when $t\rightarrow T_{l}$, at least one of
eigenvalues of $M(t)$ given in Eq.(31) will be close to one, then the
corresponding eigenvalue of $\xi(t)$ given in Eq.(39) will be close to zero,
so the energy may diverge (see Fig.1(b)). In fact, the problem of
chronological product has to be solved not only in the numerical calculation,
but also even in the experiment, because $\hat{H}_{AS}(t)$ has to be
parameterized in advance by calculating it. When considering the influence of
quantum noises, according to the results of the numerical calculation in our
example, we know the depolarizing noise is the most fatal to the dynamics
simulation of $\mathcal{PT}$-symmetric system among three kinds of quantum
noises we considered and should be avoided as much as possible.

Finally, we have to mention that in Sec.VI, we actually give an example of
time-independent Hamiltonian rather that TD Hamiltonian, just like in
Ref.[45], for the purpose of facilitate display and comparison. However, the
method of analysis and calculation is universal, and what only needs
additional attention is to be careful about the convergence of chronological
product when calculating $M(t)$ given in Eq.(31), which has been avoid in our
example.

##  VIII Acknowledgements

We acknowledge the National Key Research and Development Program of China
(Grant No. 2017YFA0303700), Beijing Advanced Innovation Center for Future Chip
(ICFC), Tsinghua University Initiative Scientific Research Program, and the
National Natural Science Foundation of China (Grant No. 11974205). C.Z. thanks
the National Natural Science Foundation of China (Grants No. 12175002 and No.
11705004), the Beijing Natural Science Foundation (Grant No. 1222020), and
NCUT Talents Project and Special Fund.

##  Appendix A Derivations of $H_{2}(t)$ and $H_{4}(t)$

Here we drive $H_{2}(t)$, $H_{4}(t)$ expressed by the operator $K(t)$ in the
main text. First of all, we know

$M^{\prime}(t)=(\xi^{2}(t)+I)^{\prime}=\xi^{\prime}(t)\xi(t)+\xi(t)\xi^{\prime}(t),$
(68)

then according to Eq.(III),

$\displaystyle H_{2}(t)=$
$\displaystyle[-i\xi^{\prime}(t)+H_{S}(t)\xi(t)-\xi(t)H_{S}(t)]M^{-1}(t)\xi(t)$
$\displaystyle=$
$\displaystyle-i\xi^{\prime}(t)M^{-1}(t)+H_{S}(t)M^{-1}(t)\cdot\xi(t)-\xi(t)H_{S}(t)M^{-1}(t)$
$\displaystyle=$
$\displaystyle-i\xi^{\prime}(t)M^{-1}(t)+[K(t)-\frac{i}{2}M^{-1}M^{\prime}(t)M^{-1}(t)]\xi(t)-\xi(t)[K(t)-\frac{i}{2}M^{-1}M^{\prime}(t)M^{-1}(t)]$
$\displaystyle=$ $\displaystyle
K(t)\xi(t)-\xi(t)K(t)-\frac{i}{2}[2\xi^{\prime}(t)M^{-1}(t)+M^{-1}(t)M^{\prime}(t)\xi(t)M^{-1}(t)-M^{-1}\xi(t)M^{\prime}(t)M^{-1}(t)]$
$\displaystyle=$ $\displaystyle
K(t)\xi(t)\\!-\\!\xi(t)K(t)\\!-\\!\frac{i}{2}\\{2\xi^{\prime}(t)M^{-1}(t)\\!+\\!M^{-1}(t)[\xi^{\prime}(t)\xi(t)\\!+\\!\xi(t)\xi^{\prime}(t)]\xi(t)M^{-1}(t)-$
$\displaystyle
M^{-1}\xi(t)[\xi^{\prime}(t)\xi(t)\\!+\\!\xi(t)\xi^{\prime}(t)]M^{-1}(t)\\}$
$\displaystyle=$ $\displaystyle
K(t)\xi(t)-\xi(t)K(t)-\frac{i}{2}\\{2\xi^{\prime}(t)M^{-1}(t)+M^{-1}(t)[\xi^{\prime}(t)(M(t)-I)M^{-1}(t)-M^{-1}(t)(M(t)-I)\xi^{\prime}(t)M^{-1}\\}$
$\displaystyle=$ $\displaystyle
K(t)\xi(t)-\xi(t)K(t)-\frac{i}{2}[\xi^{\prime}(t)M^{-1}(t)-M^{-1}(t)\xi^{\prime}(t)].$
(69)

In the similar way above, according to Eq.(III),

$\displaystyle H_{4}(t)=$
$\displaystyle[i\xi^{\prime}(t)\xi(t)+\xi(t)H_{S}(t)\xi(t)+H_{S}(t)]M^{-1}(t)$
$\displaystyle=$ $\displaystyle
i\xi^{\prime}(t)\xi(t)M^{-1}(t)+H_{S}(t)M^{-1}(t)+\xi(t)H_{S}(t)M^{-1}(t)\xi(t)$
$\displaystyle=$ $\displaystyle
i\xi^{\prime}(t)\xi(t)M^{-1}(t)+K(t)-\frac{i}{2}M^{-1}M^{\prime}(t)M^{-1}(t)+\xi(t)[K(t)-\frac{i}{2}M^{-1}M^{\prime}(t)M^{-1}(t)]\xi(t)$
$\displaystyle=$ $\displaystyle
K(t)+\xi(t)K(t)\xi(t)+\frac{i}{2}M^{-1}(t)[2M(t)\xi^{\prime}(t)\xi(t)-M^{\prime}(t)-\xi(t)M^{\prime}(t)\xi(t)]M^{-1}(t)$
$\displaystyle=$ $\displaystyle
K(t)+\xi(t)K(t)\xi(t)+\frac{i}{2}M^{-1}(t)\\{2M(t)\xi^{\prime}(t)\xi(t)-[\xi^{\prime}(t)\xi(t)+\xi(t)\xi^{\prime}(t)]-\xi(t)[\xi^{\prime}(t)\xi(t)+\xi(t)\xi^{\prime}(t)]\xi(t)\\}M^{-1}(t)$
$\displaystyle=$ $\displaystyle
K(t)+\xi(t)K(t)\xi(t)\\!+\\!\frac{i}{2}M^{-1}(t)\\{2M(t)\xi^{\prime}(t)\xi(t)\\!-\\![\xi^{\prime}(t)\xi(t)+\xi(t)\xi^{\prime}(t)]\\!-\\!\xi(t)\xi^{\prime}(t)[M(t)\\!-\\!I]-$
$\displaystyle[M(t)\\!-\\!I]\xi^{\prime}(t)\xi(t)\\}M^{-1}(t)$
$\displaystyle=$ $\displaystyle
K(t)+\xi(t)K(t)\xi(t)+\frac{i}{2}M^{-1}(t)[M(t)\xi^{\prime}(t)\xi(t)-\xi(t)\xi^{\prime}(t)M(t)]M^{-1}(t)$
$\displaystyle=$ $\displaystyle
K(t)+\xi(t)K(t)\xi(t)+\frac{i}{2}[\xi^{\prime}(t)\xi(t)M^{-1}-M^{-1}\xi(t)\xi^{\prime}(t)].$
(70)

##  Appendix B Problem of chronological product

Considering a matrix differential equation [55]:

$\displaystyle Y^{\prime}(t)=A(t)Y(t),\quad Y\left(t_{0}\right)=Y_{0},$ (71)

where $A(t)$ is a known time-dependent matrix, $Y_{0}$ is the initial value of
$Y(t)$, and $Y(t)$ is the matrix to be solved. The formal solution of the
above equation is [56, 57]:

$\displaystyle Y(t)=\mathbb{T}\exp\left(\int_{t_{0}}^{t}A(s)ds\right)Y_{0},$
(72)

where $\mathbb{T}$ is the time-ordering operator. For arbitrary two time
$t_{1}$ and $t_{2}$ ($t_{1}\neq t_{2}$), in general, $[A(t_{1}),A(t_{2})]\neq
0$, then $e^{A(t_{1})+A(t_{2})}\neq e^{A(t_{1})}\cdot e^{A(t_{2})}$, the
symbol $\mathbb{T}$ can not be ignored. When $[A(t_{1}),A(t_{2})]=0$ for
arbitrary two time $t_{1}$ and $t_{2}$, especially when $A$ is time-
independent, the symbol $\mathbb{T}$ can be ignored.

The Eq.(72) can be expressed as [57]:

$\displaystyle Y(t)=\exp\left(\Omega\left(t,t_{0}\right)\right)Y_{0},$ (73)

where $\Omega(t)$ can be written as the sum of series:

$\displaystyle\Omega(t)=\sum_{k=1}^{\infty}\Omega_{k}(t),$ (74)

and $\Omega_{n}(t)$ is the $n$-th term of Magnus series. Magnus points out
that the differential of $\Omega$ with respect to $t$ can be written as:

$\displaystyle\Omega^{\prime}=\frac{\operatorname{ad}_{\Omega}}{\exp\left(\operatorname{ad}_{\Omega}\right)-1}A,$
(75)

so the solutions of the above equation constitute Magnus expansion, or Magnus
series.

The term $\Omega_{n}$ can be obtained by $S_{n}^{(j)}$, which can be obtained
by the following recursive formula:

$\displaystyle S_{n}^{(j)}$
$\displaystyle=\sum_{m=1}^{n-j}\left[\Omega_{m},S_{n-m}^{(j-1)}\right],\quad
2\leq j\leq n-1$ $\displaystyle S_{n}^{(1)}$
$\displaystyle=\left[\Omega_{n-1},A\right],\quad
S_{n}^{(n-1)}=\operatorname{ad}_{\Omega_{1}}^{n-1}(A),$ (76)

where $\operatorname{ad}_{\Omega}^{n}$ is a shorthand for an iterated
commutator, and

$\displaystyle\operatorname{ad}_{\Omega}^{0}A=A,\quad\operatorname{ad}_{\Omega}^{k+1}A=\left[\Omega,\operatorname{ad}_{\Omega}^{k}A\right].$
(77)

For convenience, we set $t_{0}=0$. Therefore, we can get

$\displaystyle\Omega_{1}$ $\displaystyle=\int_{0}^{t}A(\tau)d\tau$
$\displaystyle\Omega_{n}$
$\displaystyle=\sum_{j=1}^{n-1}\frac{B_{j}}{j!}\int_{0}^{t}S_{n}^{(j)}(\tau)d\tau,\quad
n\geq 2,$ (78)

where $B_{j}$ are the Bernoulli numbers, and $B_{1}=-1/2$. For the convenience
of using, we write out the first four terms of $\Omega_{n}$ as follows:

$\displaystyle\begin{aligned}
\Omega_{1}(t)=&\int_{0}^{t}A\left(t_{1}\right)dt_{1}\\\
\Omega_{2}(t)=&\frac{1}{2}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\left[A\left(t_{1}\right),A\left(t_{2}\right)\right]\\\
\Omega_{3}(t)=&\frac{1}{6}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}\cdot\\\
&\left(\left[A\left(t_{1}\right),\left[A\left(t_{2}\right),A\left(t_{3}\right)\right]\right]+\left[A\left(t_{3}\right),\left[A\left(t_{2}\right),A\left(t_{1}\right)\right]\right]\right),\\\
\Omega_{4}(t)=&\frac{1}{12}\int_{0}^{t}dt_{1}\int_{0}^{t_{1}}dt_{2}\int_{0}^{t_{2}}dt_{3}\int_{0}^{t_{3}}dt_{4}\left(\left[\left[\left[A_{1},A_{2}\right],A_{3}\right],A_{4}\right]\right.\\\
&+\left[A_{1},\left[\left[A_{2},A_{3}\right],A_{4}\right]\right]+\left[A_{1},\left[A_{2},\left[A_{3},A_{4}\right]\right]\right]+\\\
&\left[A_{2},\left[A_{3},\left[A_{4},A_{1}\right]\right]\right]).\end{aligned}$
(79)

It is worth noting that the Magnus series in Eq.(74) may diverge [57], and a
sufficient condition for it to converge in $t\in[0,T)$ is :

$\displaystyle\int_{0}^{T}\|A(s)\|_{2}\mathrm{d}s<\pi,$ (80)

where $\|A\|_{2}$ denotes 2-norm of $A$.

## References

  * Griffiths and Schroeter [2018] D. J. Griffiths and D. F. Schroeter,  _Introduction to quantum mechanics_ (Cambridge university press, 2018).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Real spectra in non-hermitian hamiltonians having $\mathcal{PT}$ symmetry, Phys. Rev. Lett. 80, 5243 (1998).
  * Bender  _et al._ [1999] C. M. Bender, S. Boettcher, and P. N. Meisinger, Pt-symmetric quantum mechanics, J. Math. Phys. 40, 2201 (1999).
  * Bender  _et al._ [2002] C. M. Bender, D. C. Brody, and H. F. Jones, Complex extension of quantum mechanics, Phys. Rev. Lett. 89, 270401 (2002).
  * Bender  _et al._ [2004] C. M. Bender, D. C. Brody, and H. F. Jones, Erratum: Complex extension of quantum mechanics [phys. rev. lett. 89, 270401 (2002)], Phys. Rev. Lett. 92, 119902 (2004).
  * Mostafazadeh [2002a] A. Mostafazadeh, Pseudo-hermiticity versus pt symmetry: The necessary condition for the reality of the spectrum of a non-hermitian hamiltonian, J. Math. Phys. 43, 205 (2002a).
  * Mostafazadeh [2002b] A. Mostafazadeh, Pseudo-hermiticity versus pt-symmetry. ii. a complete characterization of non-hermitian hamiltonians with a real spectrum, J. Math. Phys. 43, 2814 (2002b).
  * Mostafazadeh [2002c] A. Mostafazadeh, Pseudo-hermiticity versus pt-symmetry iii: Equivalence of pseudo-hermiticity and the presence of antilinear symmetries, J. Math. Phys. 43, 3944 (2002c).
  * Mostafazadeh [2003] A. Mostafazadeh, Pseudo-hermiticity and generalized pt- and cpt-symmetries, J. Math. Phys. 44, 974 (2003).
  * Mostafazadeh [2007a] A. Mostafazadeh, Time-dependent pseudo-hermitian hamiltonians defining a unitary quantum system and uniqueness of the metric operator, Phys. Lett. B 650, 208 (2007a).
  * Curtright and Mezincescu [2007] T. Curtright and L. Mezincescu, Biorthogonal quantum systems, J. Math. Phys. 48, 092106 (2007).
  * Brody [2013] D. C. Brody, Biorthogonal quantum mechanics, J. Phys. A: Math. Theor. 47, 035305 (2013).
  * Zhang  _et al._ [2020] R. Zhang, H. Qin, and J. Xiao, Pt-symmetry entails pseudo-hermiticity regardless of diagonalizability, J. Math. Phys. 61, 012101 (2020).
  * Bender [2007] C. M. Bender, Making sense of non-hermitian hamiltonians, Rep. Prog. Phys. 70, 947 (2007).
  * Mostafazadeh [2010] A. Mostafazadeh, Pseudo-hermitian representation of quantum mechanics, Int J Geom Methods Mod Phys 07, 1191 (2010).
  * Makris  _et al._ [2008] K. G. Makris, R. El-Ganainy, D. N. Christodoulides, and Z. H. Musslimani, Beam dynamics in $\mathcal{P}\mathcal{T}$ symmetric optical lattices, Phys. Rev. Lett. 100, 103904 (2008).
  * Rüter  _et al._ [2010] C. E. Rüter, K. G. Makris, R. El-Ganainy, D. N. Christodoulides, M. Segev, and D. Kip, Observation of parity–time symmetry in optics, Nature physics 6, 192 (2010).
  * Peng  _et al._ [2014] B. Peng, S. K. Ozdemir, F. Lei, F. Monifi, M. Gianfreda, G. L. Long, S. Fan, F. Nori, C. M. Bender, and L. Yang, Parity-time-symmetric whispering-gallery microcavities, Nature Physics 10, 394 (2014).
  * Assawaworrarit  _et al._ [2017] S. Assawaworrarit, X. Yu, and S. Fan, Robust wireless power transfer using a nonlinear parity-time-symmetric circuit, Nature 546, 387 (2017).
  * Ashida  _et al._ [2017] Y. Ashida, S. Furukawa, and M. Ueda, Parity-time-symmetric quantum critical phenomena, Nature Communications 8, 15791 (2017).
  * Couvreur  _et al._ [2017] R. Couvreur, J. L. Jacobsen, and H. Saleur, Entanglement in nonunitary quantum critical spin chains, Phys. Rev. Lett. 119, 040601 (2017).
  * Li  _et al._ [2019] J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar, and L. Luo, Observation of parity-time symmetry breaking transitions in a dissipative floquet system of ultracold atoms, Nature communications 10, 1 (2019).
  * Zhang  _et al._ [2019a] M. Zhang, W. Sweeney, C. W. Hsu, L. Yang, A. D. Stone, and L. Jiang, Quantum noise theory of exceptional point amplifying sensors, Phys. Rev. Lett. 123, 180501 (2019a).
  * Chu  _et al._ [2020] Y. Chu, Y. Liu, H. Liu, and J. Cai, Quantum sensing with a single-qubit pseudo-hermitian system, Phys. Rev. Lett. 124, 020501 (2020).
  * Yu  _et al._ [2020] S. Yu, Y. Meng, J.-S. Tang, X.-Y. Xu, Y.-T. Wang, P. Yin, Z.-J. Ke, W. Liu, Z.-P. Li, Y.-Z. Yang, G. Chen, Y.-J. Han, C.-F. Li, and G.-C. Guo, Experimental investigation of quantum $\mathcal{P}\mathcal{T}$-enhanced sensor, Phys. Rev. Lett. 125, 240506 (2020).
  * Croke [2015] S. Croke, $\mathcal{PT}$-symmetric hamiltonians and their application in quantum information, Phys. Rev. A 91, 052113 (2015).
  * Kawabata  _et al._ [2017] K. Kawabata, Y. Ashida, and M. Ueda, Information retrieval and criticality in parity-time-symmetric systems, Phys. Rev. Lett. 119, 190401 (2017).
  * Mostafazadeh [2007b] A. Mostafazadeh, Quantum brachistochrone problem and the geometry of the state space in pseudo-hermitian quantum mechanics, Phys. Rev. Lett. 99, 130502 (2007b).
  * Bender  _et al._ [2007] C. M. Bender, D. C. Brody, H. F. Jones, and B. K. Meister, Faster than hermitian quantum mechanics, Phys. Rev. Lett. 98, 040403 (2007).
  * Günther and Samsonov [2008a] U. Günther and B. F. Samsonov, $\mathcal{P}\mathcal{T}$-symmetric brachistochrone problem, lorentz boosts, and nonunitary operator equivalence classes, Phys. Rev. A 78, 042115 (2008a).
  * Günther and Samsonov [2008b] U. Günther and B. F. Samsonov, Naimark-dilated $\mathcal{P}\mathcal{T}$-symmetric brachistochrone, Phys. Rev. Lett. 101, 230404 (2008b).
  * Ramezani  _et al._ [2012] H. Ramezani, J. Schindler, F. M. Ellis, U. Günther, and T. Kottos, Bypassing the bandwidth theorem with $\mathcal{PT}$ symmetry, Phys. Rev. A 85, 062122 (2012).
  * Zheng  _et al._ [2013] C. Zheng, L. Hao, and G. L. Long, Observation of a fast evolution in a parity-time-symmetric system, Phil. Trans. R. Soc. A 371, 20120053 (2013).
  * Beygi and Klevansky [2018] A. Beygi and S. P. Klevansky, No-signaling principle and quantum brachistochrone problem in $\mathcal{PT}$-symmetric fermionic two- and four-dimensional models, Phys. Rev. A 98, 022105 (2018).
  * Brody [2021] D. C. Brody, Pt symmetry and the evolution speed in open quantum systems, in  _Journal of Physics: Conference Series_, Vol. 2038 (IOP Publishing, 2021) p. 012005.
  * Bender  _et al._ [2013] C. M. Bender, D. C. Brody, J. Caldeira, U. Günther, B. K. Meister, and B. F. Samsonov, Pt-symmetric quantum state discrimination, Phil. Trans. R. Soc. A 371, 20120160 (2013).
  * Wang  _et al._ [2020] Y.-T. Wang, Z.-P. Li, S. Yu, Z.-J. Ke, W. Liu, Y. Meng, Y.-Z. Yang, J.-S. Tang, C.-F. Li, and G.-C. Guo, Experimental investigation of state distinguishability in parity-time symmetric quantum dynamics, Phys. Rev. Lett. 124, 230402 (2020).
  * Barnett and Andersson [2002] S. M. Barnett and E. Andersson, Bound on measurement based on the no-signaling condition, Phys. Rev. A 65, 044307 (2002).
  * Lee  _et al._ [2014] Y.-C. Lee, M.-H. Hsieh, S. T. Flammia, and R.-K. Lee, Local $\mathcal{P}\mathcal{T}$ symmetry violates the no-signaling principle, Phys. Rev. Lett. 112, 130404 (2014).
  * Brody [2016] D. C. Brody, Consistency of PT-symmetric quantum mechanics, J. Phys. A: Math. Theor. 49, 10LT03 (2016).
  * Tang  _et al._ [2016] J.-S. Tang, Y.-T. Wang, S. Yu, D.-Y. He, J.-S. Xu, B.-H. Liu, G. Chen, Y.-N. Sun, K. Sun, Y.-J. Han, _et al._ , Experimental investigation of the no-signalling principle in parity–time symmetric theory using an open quantum system, Nat. Photonics 10, 642 (2016).
  * Feng  _et al._ [2017] L. Feng, R. El-Ganainy, and L. Ge, Non-hermitian photonics based on parity–time symmetry, Nat. Photonics 11, 752 (2017).
  * Bagchi and Barik [2020] B. Bagchi and S. Barik, Remarks on the preservation of no-signaling principle in parity-time-symmetric quantum mechanics, Mod. Phys. Lett. A 35, 2050090 (2020).
  * Huang  _et al._ [2018] M. Huang, A. Kumar, and J. Wu, Embedding, simulation and consistency of pt-symmetric quantum theory, Phys. Lett. A 382, 2578 (2018).
  * Wu  _et al._ [2019] Y. Wu, W. Liu, J. Geng, X. Song, X. Ye, C.-K. Duan, X. Rong, and J. Du, Observation of parity-time symmetry breaking in a single-spin system, Science 364, 878 (2019).
  * Gui-Lu [2006] L. Gui-Lu, General quantum interference principle and duality computer, Commun. Theor. Phys. 45, 825 (2006).
  * Zheng [2018] C. Zheng, Duality quantum simulation of a general parity-time-symmetric two-level system, EPL (Europhysics Letters) 123, 40002 (2018).
  * Gao  _et al._ [2021] W.-C. Gao, C. Zheng, L. Liu, T.-J. Wang, and C. Wang, Experimental simulation of the parity-time symmetric dynamics using photonic qubits, Opt. Express 29, 517 (2021).
  * Zheng [2019] C. Zheng, Duality quantum simulation of a generalized anti-PT-symmetric two-level system, Europhys. Lett. 126, 30005 (2019).
  * Huang  _et al._ [2019] M. Huang, R.-K. Lee, L. Zhang, S.-M. Fei, and J. Wu, Simulating broken $\mathcal{PT}$-symmetric hamiltonian systems by weak measurement, Phys. Rev. Lett. 123, 080404 (2019).
  * Li  _et al._ [2022] X. Li, C. Zheng, J. Gao, and G. Long, Efficient simulation of the dynamics of an $n$-dimensional $\mathcal{PT}$-symmetric system with a local-operations-and-classical-communication protocol based on an embedding scheme, Phys. Rev. A 105, 032405 (2022).
  * Nielsen and Chuang [2002] M. A. Nielsen and I. Chuang, Quantum computation and quantum information (2002).
  * Minganti  _et al._ [2019] F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Quantum exceptional points of non-hermitian hamiltonians and liouvillians: The effects of quantum jumps, Phys. Rev. A 100, 062131 (2019).
  * Ohlsson and Zhou [2021] T. Ohlsson and S. Zhou, Density-matrix formalism for $\mathcal{PT}$-symmetric non-hermitian hamiltonians with the lindblad equation, Phys. Rev. A 103, 022218 (2021).
  * Magnus [1954] W. Magnus, On the exponential solution of differential equations for a linear operator, Commun. Pure Appl. Math. 7, 649 (1954).
  * Blanes  _et al._ [1998] S. Blanes, F. Casas, J. Oteo, and J. Ros, Magnus and fer expansions for matrix differential equations: the convergence problem, J. Phys. A: Math. Gen. 31, 259 (1998).
  * Blanes  _et al._ [2009] S. Blanes, F. Casas, J.-A. Oteo, and J. Ros, The magnus expansion and some of its applications, Physics reports 470, 151 (2009).
  * Huang  _et al._ [2021] M. Huang, R.-K. Lee, G.-Q. Zhang, and J. Wu, A solvable dilation model of pt-symmetric systems, arXiv:2104.05039 [quant-ph] (2021).
  * Zhang  _et al._ [2019b] D.-J. Zhang, Q.-h. Wang, and J. Gong, Time-dependent $\mathcal{PT}$-symmetric quantum mechanics in generic non-hermitian systems, Phys. Rev. A 100, 062121 (2019b).
  * Fring and Moussa [2016] A. Fring and M. H. Y. Moussa, Unitary quantum evolution for time-dependent quasi-hermitian systems with nonobservable hamiltonians, Phys. Rev. A 93, 042114 (2016).
  * Luiz  _et al._ [2020] F. S. Luiz, M. A. de Ponte, and M. H. Y. Moussa, Unitarity of the time-evolution and observability of non-hermitian hamiltonians for time-dependent dyson maps, Physica Scripta 95, 065211 (2020).
  * Brody and Graefe [2012] D. C. Brody and E.-M. Graefe, Mixed-state evolution in the presence of gain and loss, Phys. Rev. Lett. 109, 230405 (2012).
  * Xiao  _et al._ [2019] L. Xiao, K. Wang, X. Zhan, Z. Bian, K. Kawabata, M. Ueda, W. Yi, and P. Xue, Observation of critical phenomena in parity-time-symmetric quantum dynamics, Phys. Rev. Lett. 123, 230401 (2019).
  * Bender  _et al._ [2003] C. M. Bender, D. C. Brody, and H. F. Jones, Must a hamiltonian be hermitian?, Am. J. Phys. 71, 1095 (2003).
  * Andersson  _et al._ [2007] E. Andersson, J. D. Cresser, and M. J. Hall, Finding the kraus decomposition from a master equation and vice versa, J. Mod. Opt. 54, 1695 (2007).