#  Entropy and entanglement in a bipartite
quasi-Hermitian system and its Hermitian counterparts

Abed Alsalam Abu Moise Graham Cox Marco<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Mathematics and Statistics
Memorial University of Newfoundland
St. John’s, NL A1C 5S7
Canada

###### Abstract

We consider a quantum oscillator coupled to a bath of $N$ other oscillators.
The total system evolves with a quasi-Hermitian Hamiltonian. Associated to it
is a family of Hermitian systems, parameterized by a unitary map $W$. Our main
goal is to find the influence of $W$ on the entropy and the entanglement in
the Hermitian systems. We calculate explicitly the reduced density matrix of
the single oscillator for all Hermitian systems and show that, regardless of
$W$, their von Neumann entropy oscillates with a common period which is twice
that of the non-Hermitian system. We show that generically, the oscillator and
the bath are entangled for almost all times. While the amount of entanglement
depends on the choice of $W$, it is independent of $W$ when averaged over a
period. These results describe some universality in the physical properties of
all Hermitian systems associated to a given non-Hermitian one.

##  1 Motivation & outline of main results

In recent years there has been much interest in extensions of quantum
mechanics that allow for non-Hermitian Hamiltonians. In $PT$-symmetric quantum
theory, for instance, non-Hermitian Hamiltonians with purely real spectrum are
often encountered. In many cases those Hamiltonians are quasi-Hermitian [1, 2,
3, 6, 20, 21, 22, 23, 26, 29]. For these one can produce associated Hermitian
Hamiltonians which can then be studied using ‘standard’ methods of quantum
theory.

The assignment of a Hermitian to a non-Hermitian system is called the Dyson
map [12, 8, 9]. The Dyson map is determined up to a unitary. To a given non-
Hermitian Hamiltonian $H$ are assigned Hermitian Hamiltonians $h$ which are
all unitarily equivalent. In this sense, physical properties of the original
$H$ are uniquely encoded by any associated $h$. However, if the system has a
subsystem structure, then the unitaries generally reshuffle local degrees of
freedom and the choice of the unitary plays a crucial physical role. This is
in particular the case for open systems, where system and bath degrees of
freedom are reshuffled. One should then ask what relevant information about
the original non-Hermitian system is encoded in ‘the’ associated Hermitian
one. In fact, if one allows the Dyson map to be time-dependent, as is often
the approach taken in the literature [11, 12, 14], then the above mentioned
non-uniqueness is compounded by additional freedom. It turns out that given
any quasi-Hermitian system $H$, one can chose a suitable time-dependent Dyson
map so that the associated Hermitian system has the trivial Hamiltonian $h=0$.
We explain this in Appendix A.3.

The ambiguity in the associated Hermitian systems is the subject of the
present paper. In the recent literature on $PT$-symmetric quantum theory it is
suggested that instead of studying the entropy of an open quasi-Hermitian
systems, one may study the entropy of the corresponding associated Hermitian
system [13, 17, 18, 5]. Concretely, in [13], the authors examine the entropy
of an oscillator coupled to an ‘environment’ of $N$ other oscillators via a
$PT$-symmetric Hamiltonian. As a proxy for the von Neumann entropy of the
single oscillator in the non-Hermitian system, they study the corresponding
quantity for one particular choice of an associated Hermitian system,
constructed by a time-dependent Dyson map. The authors find a different
qualitative behaviour of the entropy of the Hermitian system according to the
$PT$-symmetry phases of the non-Hermitian system. This interesting result is
found for one particular choice of the associated Hermitian system – the
choice may be regarded as a good one as it discerns the different symmetry
phases by the qualitative form of the dynamics of the entropy. For a different
choice however, the entropy dynamics would look entirely different (and in
particular, it would be time-independent upon chosing $h=0$, as explained
above).

In the current work we investigate all associated systems in the regime of
unbroken $PT$-symmetry – in which the non-Hermitian Hamiltonian is actually
quasi-Hermitian. We vary over all time-independent Dyson maps (excluding time
dependent ones for the reason mentioned above).

A Hamiltonian $H$ is said to be _quasi-Hermitian_ if

$H^{\dagger}=\eta H\eta^{-1}$

for some bounded, positive operator $\eta>0$, called a _metric operator_. One
obtains a Hermitian Hamiltonian by

$h=SHS^{-1},$

where $S$ is any operator with

$S^{\dagger}S=\eta.$

The assignment $H\mapsto h$ is called the Dyson map [8, 9], even though to a
single $H$ one can associate different $h$ for the following two reasons.

  * –

There are different metric operators $\eta$ for the same quasi-Hermitian $H$.

  * –

Once a metric $\eta$ is fixed, the general solution of the equation
$S^{\dagger}S=\eta$ is $S=W\sqrt{\eta}$, where $\sqrt{\eta}$ is the unique
positive operator that squares to $\eta$ and $W$ is an arbitrary unitary $W$.

The metric $\eta$ is uniquely determined if one fixes a sufficiently rich
(irreducible) set of quasi-Hermitian operators (including $H$) to be
observables; see [29] and Section 2. Once this choice is made, all associated
Hermitian systems are parameterized by the unitary $W$. Different choices of
$W$ give different Hermitian $h$, which are unitarily equivalent. However,
generally $W$ reshuffles the degrees of freedom and modifies local physical
properties, such as entanglement between subsystems.

In this work, we consider a concrete quasi-Hermitian system with Hamiltonian
$H$, describing the interaction of a single oscillator with a bath of $N$
other oscillators. This is an open quantum system of the ‘system–bath’ (${\rm
S}{\rm B}$) type, described by a bipartite Hilbert space ${\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}$. We take the metric $\eta$ to be general but
fixed, and then analyze all the associated Hermitian systems resulting from
varying over all choices of $W$. The dynamics of the non-Hermitian system and
associated Hermitian systems is given by the time-dependent wave functions

$|\psi(t)\rangle=e^{-itH}|\psi(0)\rangle\quad\mbox{and}\quad|\phi(t)\rangle=S|\psi(t)\rangle=W\sqrt{\eta}|\psi(t)\rangle,$

respectively. The vector $|\psi(t)\rangle$ is normalized with respect to the
metric induced by the inner product $\langle\cdot|\eta|\cdot\rangle$, while
$|\phi(t)\rangle$ is normalized relative to the ‘original’ inner product
$\langle\cdot|\cdot\rangle$. The reduced density matrices for the non-
Hermitian and the Hermitian systems are obtained by taking the partial trace
over the bath degrees of freedom (see Section 2.2),

$\bar{\rho}_{H}(t)={\rm tr}_{\rm
B}\big{(}|\psi(t)\rangle\langle\psi(t)|\eta\big{)}\quad\mbox{and}\quad\bar{\rho}_{h_{W}}(t)={\rm
tr}_{\rm B}\big{(}|\phi(t)\rangle\langle\phi(t)|\big{)}.$

Our main findings are summed up as follows:

  * 1.

Explicit evolution. We obtain explicit formulas for the states
$|\psi(t)\rangle$ and $|\phi(t)\rangle$ as well as the reduced states
$\bar{\rho}_{H}(t)$ and $\bar{\rho}_{h_{W}}(t)$. See Sections 3.2, 3.3.

  * 2.

Metric $\eta$. The operator $\bar{\rho}_{H}(t)$ is a density matrix exactly
when $\eta$ is of product form $\Lambda_{\rm S}\otimes\Lambda_{\rm B}$ –
otherwise $\bar{\rho}_{H}(t)$ has complex eigenvalues. See Section 3.1. We
thus take $\eta$ of product form in the further analysis.

  * 3.

Subsystem entropy. The reduced states $\bar{\rho}_{H}(t)$ and
$\bar{\rho}_{h_{W}}(t)$ are periodic222The whole ${\rm S}{\rm B}$ complex
consists of $N+1$ oscillators, so the energy spectrum of all the Hamiltonians
involved consists of discrete eigenvalues only, without continuous spectrum.
This explains the periodicity. in time, both having the same period regardless
of the choice of $W$. The von Neumann entropy $S=-{\rm tr}(\rho\ln\rho)$ of
$\bar{\rho}_{H}(t)$ and $\bar{\rho}_{h_{W}}(t)$ is periodic in time as well,
but for generic initial conditions333Only for specially tuned initial states
and $W$’s is the situation non-generic, see Section 3.3. and generic $W$, the
period of the entropy of the Hermitian system is double that of the non-
Hermitian system. See section 4.

  * 4.

${\rm S}{\rm B}$ entanglement. The non-Hermitian and the Hermitian ${\rm
S}{\rm B}$ states $|\psi(t)\rangle$, $|\phi(t)\rangle$ are entangled for all
times except at periodically reoccurring single instants.444For a class of
exceptional initial conditions the states are entangled for all times, see
Section 5. Given any entangled state $|\psi\rangle$, one can find $W$ such
that the associated $|\phi\rangle$ is disentangled, and for any disentangled
$|\psi\rangle$ there are $W$ such that $|\phi\rangle$ is entangled. However,
in an averaged sense, the choice of $W$ does not influence the entanglement at
all. Namely, the concurrence of the time-averaged density matrix,
$\langle\rho\rangle=\frac{1}{T}\int_{0}^{T}|\phi(t)\rangle\langle\phi(t)|\,dt$,
where $T$ is the period of $|\phi(t)\rangle\langle\phi(t)|$, is independent of
$W$. Its value is determined entirely by the initial condition and the choice
of the metric. We identify the initial states for which $\langle\rho\rangle$
is separable and for which it is maximally entangled. See Section 5.

##  2 Quasi-Hermitian systems

Let $\mathcal{H}$ be a finite-dimensional Hilbert space with inner product
$\langle\cdot|\cdot\rangle$. An operator $\eta$ is said to be positive,
denoted as $\eta>0$, if $\langle\psi|\eta\psi\rangle>0$ for all nonzero
$|\psi\rangle\in{\mathcal{H}}$. This is equivalent with saying that
$\eta^{\dagger}=\eta$ and all eigenvalues of $\eta$ are strictly positive.
Here $A^{\dagger}$ is the adjoint of the operator $A$, defined by
$\langle\psi|A\phi\rangle=\langle A^{\dagger}\psi|\phi\rangle$ for all
$|\phi\rangle,|\psi\rangle\in{\mathcal{H}}$. An operator $H$ on
${\mathcal{H}}$ is called ($\eta$-)quasi-Hermitian if there exists a positive
operator $\eta>0$ such that

$H^{\dagger}=\eta H\eta^{-1}.$ (1)

Quasi-hermiticity is a special case of pseudo-hermiticity, where (1) holds
with an invertible (but not necessarily positive) Hermitian operator $\eta$.
Pseudo- and quasi-hermitian Hamiltonians arise in $PT$-symmetric quantum
theory, see for instance [31] and references therein.

###  2.1 Hermitian counterparts

Let $H$ be a non-Hermitian operator on ${\mathcal{H}}$, a candidate for the
Hamiltonian of a physical system. In order to obtain a Hermitian quantum
theory, one could either:

  * (1)

Modify the inner product of ${\mathcal{H}}$ to
$\langle\cdot|\eta\,\cdot\rangle$ for some $\eta>0$ (called a metric
operator), such that $H$ becomes Hermitian in the Hilbert space
${\mathcal{H}}_{\eta}$ with this new inner product; or

  * (2)

Take a similarity transformation (invertible map) $S$ such that the
transformed $h=SHS^{-1}$ is Hermitian in the original Hilbert space
${\mathcal{H}}$.

If $H$ is quasi-Hermitian, then both options (1) and (2) are possible, but
neither the metric nor the similarity transform in options (1) and (2) are
unique. To explore this non-uniqueness, we first notice that any quasi-
Hermitian $H$ is diagonalizable [7, 10, 21]. More precisely,

$H=\sum_{n=1}^{N}E_{n}|\psi_{n}\rangle\langle\phi_{n}|,$ (2)

where the $E_{n}\in{\mathbb{R}}$ are the eigenvalues and the
$\\{|\psi_{n}\rangle,|\phi_{n}\rangle\\}_{n=1}^{N}$ form a complete bi-
orthonormal family, meaning that $\langle\psi_{k}|\phi_{l}\rangle=\delta_{kl}$
and $\sum|\psi_{n}\rangle\langle\phi_{n}|={\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$. Let us consider the
case where all eigenvalues $E_{n}$ are distinct for simplicity (a discussion
including degenerate eigenvalues can be done similarly, but this is not our
focus here). Then the decomposition (2) is unique, it is the spectral
representation of the operator $H$, and the
$P_{n}\equiv|\psi_{n}\rangle\langle\phi_{n}|$ are the uniquely defined
(generally not orthogonal) spectral projections. The vectors
$|\psi_{n}\rangle$ and $|\phi_{n}\rangle$, however, are determined only up to
a joint scaling $|\phi_{n}\rangle\mapsto z_{n}|\phi_{n}\rangle$ and
$|\psi_{n}\rangle\mapsto\frac{1}{\overline{z_{n}}}|\psi_{n}\rangle$, with
$0\neq z_{n}\in\mathbb{C}$ arbitrary.

$\bullet$ First let us explore the option (1). A metric $\eta$ is called a
metric for $H$ if $H$ is $\eta$-quasi-Hermitian. Let $A$ be a linear operator
on ${\mathcal{H}}$, with adjoint $A^{\dagger}$ as defined above. If $A$ is
viewed as an operator on ${\mathcal{H}}_{\eta}$, then $\langle\phi|\eta
A\psi\rangle=\langle\eta^{-1}A^{\dagger}\eta\phi|\eta\psi\rangle$, so the
adjoint of $A$ in ${\mathcal{H}}_{\eta}$ is
$A^{\ddagger}=\eta^{-1}A^{\dagger}\eta$. It follows that a given $\eta>0$ is a
metric for $H$ if and only if $H^{\ddagger}=H$, that is, if and only if $H$ is
Hermitian acting on ${\mathcal{H}}_{\eta}$. It is well known (see e.g. [26,
10]) that

$\eta\mbox{\ \, is a metric for
$H$}\quad\Longleftrightarrow\quad\eta=\sum_{n=1}^{N}x_{n}|\phi_{n}\rangle\langle\phi_{n}|\mbox{\quad
for some $x_{1},\ldots,x_{N}>0$},$ (3)

where the $|\phi_{n}\rangle$ are the vectors appearing in (3). The multitude
of metrics obtained by varying the $x_{j}$ in (2) naturally appears due to the
fact that $|\phi_{n}\rangle$ is only determined up to an arbitrary nonzero
scaling factor $z_{n}$ (as explained after (2)), which results in the scaling
$x_{n}\mapsto x_{n}|z_{n}|^{2}$. Given this non-uniqueness of the metric,
which one should be chosen to define the physical Hilbert space
${\mathcal{H}}_{\eta}$?

One answer is that the metric is fixed provided that instead of just $H$, one
chooses a whole irreducible family of operators to be Hermitian observables.
Namely, it is shown in [29] (see also [25] for the two-dimensional case) that
if there is a family of operators $\\{A_{i}\\}_{i}$ on $\mathcal{H}$, and
positive operators $\eta$, $\eta^{\prime}$ such that $A_{i}^{\dagger}=\eta
A_{i}\eta^{-1}$ and $A_{i}^{\dagger}=\eta^{\prime}A_{i}(\eta^{\prime})^{-1}$
for all $i$, then

$\mbox{$\eta^{\prime}$ is a scalar multiple of
$\eta$}\quad\Longleftrightarrow\quad\mbox{$\\{A_{i}\\}_{i}$ is an irreducible
family of operators on $\mathcal{H}$.}$

This means that for an irreducible family of quasi-Hermitian operators, there
is exactly one metric (up to a scalar multiple) that makes those operators
Hermitian. The chosen family555Examples of irreducible families are the Pauli
matrices for a spin, with the Euclidean inner product on ${\mathbb{C}}^{2}$,
or the position $\hat{x}$ and momentum $\hat{p}=-i\hbar\nabla_{x}$ for a
quantum particle (rather, the bounded Weyl operators generated by them) with
the inner product
$\langle\psi|\phi\rangle=\int_{{\mathbb{R}}^{3}}\bar{\psi}(x)\phi(x)\,d^{3}x$.
can then be viewed as the physical observables of the theory and the space of
pure states is $\mathcal{H}$ with inner product
$\langle\cdot|\cdot\rangle_{\eta}$.

On the other hand, if interested only in the single observable $H$ (the
Hamiltonian), one should keep the $x_{n}$ in (3) general.

$\bullet$ Next, let us investigate option (2) for $H$ of the form (2). Let
$\eta$ be a metric for $H$, so it is of the form (3). We find all invertible
$S$ such that the transformed $h\equiv SHS^{-1}$ is Hermitian,

$h=SHS^{-1}=\big{(}SHS^{-1}\big{)}^{\dagger}=h^{\dagger}.$ (4)

One readily sees that (4) is equivalent to $TH=HT$, where
$T=\eta^{-1}S^{\dagger}S$. That $T$ commutes with $H$, as in (2), is
equivalent to $T$ being diagonal in the same bi-orthonormal system as $H$,
that is, $T=\sum_{n=1}^{N}t_{n}|\psi_{n}\rangle\langle\phi_{n}|$ for some
$t_{n}\in\mathbb{C}$. Now

$S^{\dagger}S=\eta
T=\Big{(}\sum_{n=1}^{N}x_{n}|\phi_{n}\rangle\langle\phi_{n}|\Big{)}\Big{(}\sum_{k=1}^{N}t_{k}|\psi_{k}\rangle\langle\phi_{k}|\Big{)}=\sum_{n=1}^{N}x_{n}t_{n}|\phi_{n}\rangle\langle\phi_{n}|,$
(5)

and as $S^{\dagger}S>0$ and $x_{n}>0$, we have $t_{n}>0$ as well. It follows
from (3) and (5) that $\eta T$ is also a metric for $H$. In fact, (3) and (5)
show that given a fixed metric $\eta$ for $H$ and varying $\eta T$ over all
operators $T>0$ that commute with $H$, we obtain all of the metrics for $H$.
We conclude that given $\eta$, the $S$ we are looking for are the solutions of
$S^{\dagger}S=\eta T$, where $T>0$ is an operator that commutes with $H$
(equivalently, is diagonal in the same bi-orthonormal system as $H$). The
general solution is

$S=W\sqrt{\eta T},$ (6)

where $W$ is any unitary and where for a positive operator $A$, $\sqrt{A}$ is
the unique positive operator whose square equals $A$.

Once $W$ and $T$ are chosen, the associated Hermitian $h$ in (4) becomes

$h_{W,T}=W\sqrt{\eta T}\,H\frac{1}{\sqrt{\eta T}}W^{\dagger}.$ (7)

We stress with this notation that $h$ depends on the choice of $W$ and $T$.
The $h$ obtained from two different choices of unitaries, say $V$ and $W$, are
unitarily equivalent, with $h_{V,T}=Uh_{W,T}U^{\dagger}$ and $U=VW^{\dagger}$.
In this sense, the choice of $W$ is globally immaterial. However, if the
Hilbert space has a local structure, say is of bipartite nature
${\mathcal{H}}={\mathcal{H}}_{\rm S}\otimes{\mathcal{H}}_{\rm B}$, then the
global unitary $U$ may well change the local properties of the two local
subsystems, in which case the choice of $W$ will play a physically relevant
role. We also point out that the spectrum of $h_{W,T}$ does not depend on
either $W$ or $T$ (or $\eta$, for that matter).

In this work we take the following approach: We start with a given quasi-
Hermitian Hamiltonian $H$ and an arbitrary metric $\eta$ for $H$ and we view
${\mathcal{H}}_{\eta}$ as the physical Hilbert space. We analyze the class of
all associated Hermitian systems $h_{W,T}$, where $W$ and $T$ vary over all
unitaries and all positive operators commuting with $H$, respectively. As
explained above, considering all metrics $\eta$ is the same as considering all
metrics $\eta T$, so varying over $T$ is redundant if $\eta$ is kept
arbitrary. We may then set $T={\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$ and only consider

$S=W\sqrt{\eta},\qquad h_{W}=W\sqrt{\eta}\,H\frac{1}{\sqrt{\eta}}W^{\dagger}$
(8)

for all $W$ and $\eta$.

###  2.2 States, reduced states, von Neumann entropy

Consider now a fixed metric $\eta$, so that the physical Hilbert space is
${\mathcal{H}}_{\eta}$ and $H$ is Hermitian on ${\mathcal{H}}_{\eta}$,
$H^{\ddagger}=H$. Then $e^{-itH}$ is the unitary Schrödinger dynamics on
${\mathcal{H}}_{\eta}$. The average of an observable $A$ on
${\mathcal{H}}_{\eta}$ in the state $\psi\in{\mathcal{H}}_{\eta}$ is given by

$\langle\psi|A\psi\rangle_{\eta}=\langle\psi|\eta A\psi\rangle={\rm
tr}\big{(}|\psi\rangle\langle\psi|\eta A\big{)}={\rm
tr}\big{(}\widetilde{\rho}A\big{)},$ (9)

where

$\widetilde{\rho}=|\psi\rangle\langle\psi|\eta$ (10)

is a density matrix on ${\mathcal{H}}_{\eta}$ (a positive, trace-one
operator). This $\widetilde{\rho}$ is called the ‘generalized density matrix’
in [28]. It is important to point out that the trace in (9) is a purely
algebraic quantity: it is the sum of the eigenvalues of the operator, and
therefore does not depend on the choice of metric.

To arrive at a Hermitian Hamiltonian, it is necessary to make a choice for the
unitary $W$ in (8). The associated Hermitian Hamiltonian $h_{W}$ is then given
by (8). Let

$|\psi(t)\rangle=e^{-itH}|\psi(0)\rangle,\quad|\phi(t)\rangle=e^{-ith_{W}}|\phi(0)\rangle$
(11)

be the evolution of the initial states $|\psi(0)\rangle$, $|\phi(0)\rangle$
with respect to $H$ and $h_{W}$, respectively. The states are related by

$|\phi(t)\rangle=S|\psi(t)\rangle,\quad S=W\sqrt{\eta},$ (12)

and the density matrices associated to these vector states for the non-
Hermitian (see (10)) and the Hermitian systems are

$\rho_{H}(t)=|\psi(t)\rangle\langle\psi(t)|\eta\quad\mbox{and}\quad\rho_{h_{W}}(t)=|\phi(t)\rangle\langle\phi(t)|,$
(13)

respectively. (We adopt the notation $\rho_{h_{W}}$ and $\rho_{H}$ for the
density matrices on the Hermitian and non-Hermitian sides of the problem from
[13].) It is clear from (12) that

$\rho_{{h_{W}}}(t)=S|\psi(t)\rangle\langle\psi(t)|S^{\dagger}=S|\psi(t)\rangle\langle\psi(t)|(S^{\dagger}S)S^{-1}=S\rho_{H}(t)S^{-1}.$
(14)

It follows that $\rho_{h_{W}}(t)$ and $\rho_{H}(t)$ have the same eigenvalues,
and hence the same von Neumann entropy,
${\mathcal{E}}(\rho_{h_{W}}(t))={\mathcal{E}}(\rho_{H}(t))$, where

${\mathcal{E}}(\rho)=-{\rm
tr}\big{(}\rho\ln\rho\big{)}=-\sum_{i}\lambda_{i}\ln\lambda_{i}$ (15)

and $\\{\lambda_{i}\\}$ are the eigenvalues of $\rho$.

Consider now a bipartite system with ${\mathcal{H}}={\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}$ (‘system’ and ‘bath’). We consider the reduced
states (denoted by an overbar) defined by tracing out the degrees of freedom
of the subsystem ${\mathcal{H}}_{\rm B}$,

$\bar{\rho}_{H}(t)={\rm tr}_{{\mathcal{H}}_{\rm
B}}\big{(}\rho_{H}(t)\big{)},\qquad\bar{\rho}_{h_{W}}(t)={\rm
tr}_{{\mathcal{H}}_{\rm B}}\big{(}\rho_{h_{W}}(t)\big{)}.$ (16)

In some recent works [13, 17, 18, 5], the dynamics of a bipartite system
generated by a non-Hermitian Hamiltonian $H$ is studied, with particular focus
on the von Neumann entropy of the reduced density matrix $\bar{\rho}_{H}(t)$.
The strategy proposed in those works is to examine the entropy of
$\bar{\rho}_{h_{W}}(t)$ as a proxy for that of $\bar{\rho}_{H}(t)$. In this
respect, however, one should observe the following facts:

  * 1.

The operator $\bar{\rho}_{H}(t)$ always satisfies ${\rm
tr}_{{\mathcal{H}}_{\rm S}}(\bar{\rho}_{H}(t))=1$, but for some choices of
$\eta$ the eigenvalues of $\bar{\rho}_{H}(t)$ can be complex, in which case it
is not a valid density matrix.

  * 2.

Even if the metric $\eta$ is chosen such that $\bar{\rho}_{H}(t)$ is a density
matrix, for generic choices of $W$ the von Neumann entropies
${\mathcal{E}}(\bar{\rho}_{H}(t))$ and ${\mathcal{E}}(\bar{\rho}_{h_{W}}(t))$
are not the same. The latter in fact depends on the choice of $W$.

To understand the normalization of the trace mentioned in fact 1. above, we
observe (using ${\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{\rm S}$ as the system observable) that

${\rm tr}_{{\mathcal{H}}_{\rm S}}(\bar{\rho}_{H}(t))={\rm
tr}_{{\mathcal{H}}_{\rm S}}(\bar{\rho}_{H}(t){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\rm S})={\rm tr}_{{\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}}\big{(}\rho_{H}(t)({\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\rm S}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{\rm B})\big{)}={\rm
tr}_{{\mathcal{H}}_{\rm S}\otimes{\mathcal{H}}_{\rm B}}(\rho_{H}(t))=1.$

If $S=S_{\rm S}\otimes S_{\rm B}$, then $\bar{\rho}_{h}=S_{\rm
S}\bar{\rho}_{H}S_{\rm S}^{-1}$ and so the spectra and thus the von Neumann
entropies of $\bar{\rho}_{h}$ and $\bar{\rho}_{H}$ coincide. However, if $S$
is entangling (not of product form $S_{\rm S}\otimes S_{\rm B}$), then the
eigenvalues of the two reduced density matrices are not the same in general,
and neither are their entropies.

These difficulties are resolved in the next section, where we study the
concrete model used in [13]. In particular, we determine for which choices of
$\eta$ the reduced operator $\bar{\rho}_{H}(t)$ is indeed a density matrix,
and then we find the von Neumann entropy of $\bar{\rho}_{h_{W}}(t)$ for all
possible choices of the unitary $W$.

##  3 Model

An oscillator with creation and annihilation operators $a^{\dagger}$, $a$ is
coupled to a ‘bath’ of $N$ independent oscillators with creation and
annihilation operators $q_{i}^{\dagger}$, $q_{i}$, $i=1,\ldots,N$. The total
Hilbert space of the $N+1$ oscillators is

$\mathcal{H}=\mathcal{H}_{\rm S}\otimes\mathcal{H}_{\rm B},$ (17)

where $\mathcal{H}_{\rm S}$ is the space of a single oscillator and
$\mathcal{H}_{\rm B}$ is that of the other $N$. As in the previous section, we
denote the inner product by $\langle\cdot|\cdot\rangle$ and let $\dagger$
denote the adjoint in this inner product. The commutation relations are
$[a,a^{\dagger}]=1=[q_{i},q^{\dagger}_{i}]$, and all operators belonging to
different oscillators commute. This open quantum system, non-Hermitian model
was used in [13].

###  3.1 The quasi-Hermitian system

The coupled total system–bath Hamiltonian is

$H=\nu N_{\rm tot}+(g+\kappa)\sqrt{N}\,a^{\dagger}Q+(g-\kappa)\sqrt{N}\
aQ^{\dagger},$ (18)

where $\nu>0$ and $g,\kappa\in\mathbb{R}$ are parameters and

$N_{\rm tot}=a^{\dagger}a+\sum_{n=1}^{N}q_{n}^{\dagger}q_{n},\qquad
Q=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}q_{n}.$ (19)

Due to the different prefactors of $\kappa$ in the interaction term of (18),
$H$ is $\dagger$-Hermitian if and only if $\kappa=0$.

The ‘uncoupled’ ($g=\kappa=0$) Hamiltonian is simply $\nu N_{\rm tot}$, a
multiple of the total number operator $N_{\rm tot}$. As $H$ commutes with
$N_{\rm tot}$, each eigenspace of $N_{\rm tot}$, with a fixed number of
excitations (in the system plus the bath) is left invariant. Denote by
$|0_{\rm S}0_{\rm B}\rangle$ the ‘vacuum’ zero excitation state, where all
oscillators are in the ground state. The single excitation space is defined as

${\mathcal{E}}_{1}={\rm span}\big{\\{}|1_{\rm S}0_{\rm B}\rangle,|0_{\rm
S}1_{1}\rangle,|0_{\rm S}1_{2}\rangle,\ldots,|0_{\rm S}1_{N}\rangle\big{\\}},$
(20)

where $|1_{\rm S}0_{\rm B}\rangle=a^{\dagger}|0_{\rm S}0_{\rm B}\rangle$ and
$|0_{\rm S}1_{i}\rangle=q^{\dagger}_{i}|0_{\rm S}0_{\rm B}\rangle$ for
$i=1,\ldots,N$. When $H$ is applied to a vector in ${\mathcal{E}}_{1}$ the
result is again a vector in ${\mathcal{E}}_{1}$. Moreover, due to the
collective, symmetric nature of the system-bath interaction in (18), $H$
leaves the even smaller space

${\mathcal{H}}_{1}={\rm span}\big{\\{}|e_{\rm S}\rangle,\ |e_{\rm
B}\rangle\big{\\}}$ (21)

invariant, where

$|e_{\rm S}\rangle=|1_{\rm S}0_{\rm B}\rangle,\qquad|e_{\rm
B}\rangle=\frac{1}{\sqrt{N}}\sum_{n=1}^{N}|0_{\rm S}1_{n}\rangle.$ (22)

Those two vectors describe states in which a single excitation is either in
${\rm S}$ (the state $|e_{\rm S}\rangle$) or in ${\rm B}$, collectively spread
over the $N$ bath oscillators (the state $|e_{\rm B}\rangle$). Therefore, we
may view $H$ as an operator on ${\mathcal{H}}_{1}$. When we do this we denote
it by $H_{1}$, which has the form

$H_{1}=\nu{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}+(g-\kappa)\sqrt{N}\ |e_{\rm
B}\rangle\langle e_{\rm S}|+(g+\kappa)\sqrt{N}\ |e_{\rm S}\rangle\langle
e_{\rm B}|.$ (23)

The eigenvalues of $H_{1}$ are

$\omega_{\pm}=\nu\pm\omega,\qquad\omega=\sqrt{N}\sqrt{g^{2}-\kappa^{2}},$ (24)

which are real for $\kappa^{2}\leq g^{2}$ and (purely imaginary) complex
conjugates for $\kappa^{2}>g^{2}$. See [13] for a discussion of the $PT$
symmetry of $H$. The operator $H_{1}$ is diagonalizable except at the
transition points defined by $\kappa^{2}=g^{2}\neq 0$, where $H_{1}$ reduces
to a Jordan block. Note that increasing the number $N$ of oscillators in the
bath simply amounts to speeding up the dynamics (the frequency $\omega$) by a
factor $\sqrt{N}$.

We consider the ‘$PT$-symmetry unbroken regime’ $\kappa^{2}<g^{2}$, so that
$\omega_{\pm}\in{\mathbb{R}}$. For definiteness we take $g>0$ (the case $g<0$
can be dealt with in the same fashion), so

$0\leq|\kappa|<g,$ (25)

which is equivalent to $g+\kappa>0$ and $g-\kappa>0$. Then we have $\omega>0$
and

$a_{1}=\sqrt{g+\kappa}>0,\quad a_{2}=\sqrt{g-\kappa}>0,$ (26)

where the equalities in (26) define the quantities $a_{1}$, $a_{2}$. The two
linearly independent (not normalized) eigenvectors of $H_{1}$ and its adjoint
$H_{1}^{\dagger}$ are

$|v_{\pm}\rangle\propto a_{1}|e_{\rm S}\rangle\pm a_{2}|e_{\rm
B}\rangle\quad\mbox{and}\quad|v^{*}_{\pm}\rangle\propto a_{2}|e_{\rm
S}\rangle\pm a_{1}|e_{\rm B}\rangle,$

respectively. They satisfy $H_{1}|v_{\pm}\rangle=\omega_{\pm}|v_{\pm}\rangle$
and $H^{\dagger}_{1}|v^{*}_{\pm}\rangle=\omega_{\pm}|v^{*}_{\pm}\rangle$. Note
that $|v_{\pm}^{*}\rangle$ denote the eigenvectors of $H^{\dagger}$, not to be
confused with the complex conjugates of the eigenvectors $|v_{\pm}\rangle$ of
$H$. We normalize the vectors as

$|v_{\pm}\rangle=\tfrac{1}{\sqrt{2}}\Big{(}\sqrt{\tfrac{a_{1}}{a_{2}}}|e_{\rm
S}\rangle\pm\sqrt{\tfrac{a_{2}}{a_{1}}}|e_{\rm
B}\rangle\Big{)}\quad\mbox{and}\quad|v^{*}_{\pm}\rangle=\tfrac{1}{\sqrt{2}}\Big{(}\sqrt{\tfrac{a_{2}}{a_{1}}}|e_{\rm
S}\rangle\pm\sqrt{\tfrac{a_{1}}{a_{2}}}|e_{\rm B}\rangle\Big{)}.$ (27)

Then $\\{|v_{\pm}\rangle,|v^{*}_{\pm}\rangle\\}$ is a bi-orthonormal basis,
satisfying $\langle v^{*}_{\pm}|v_{\mp}\rangle=0$ and $\langle
v^{*}_{\pm}|v_{\pm}\rangle=1$, and the operator $H_{1}$ can be written as

$H_{1}=\omega_{+}|v_{+}\rangle\langle
v^{*}_{+}|+\omega_{-}|v_{-}\rangle\langle v^{*}_{-}|.$ (28)

Using this, one easily finds

$\displaystyle e^{-itH_{1}}$ $\displaystyle=$ $\displaystyle
e^{-it\omega_{+}}|v_{+}\rangle\langle
v^{*}_{+}|+e^{-it\omega_{-}}|v_{-}\rangle\langle v^{*}_{-}|$ (29)
$\displaystyle=$ $\displaystyle e^{-it\nu}\cos(\omega t){\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}-ie^{-it\nu}\sin(\omega t)\Big{(}\frac{a_{1}}{a_{2}}|e_{\rm
S}\rangle\langle e_{\rm B}|+\frac{a_{2}}{a_{1}}|e_{\rm B}\rangle\langle e_{\rm
S}|\Big{)}.$

We consider initial states which are vectors in ${\mathcal{H}}_{1}$, as
defined in (21), so the dynamics generated by $H$ is entirely given by the
operator $H_{1}$ from (23). We still consider the regime (25), so that the
spectrum of $H_{1}$ consists of two distinct real eigenvalues. Comparing (2),
(3) and (28), we see that $H_{1}$ is quasi-Hermitian and the set of all
associated metrics is

${\mathcal{M}}_{+}=\big{\\{}\eta=x_{1}|v_{+}^{*}\rangle\langle
v_{+}^{*}|+x_{2}|v_{-}^{*}\rangle\langle v_{-}^{*}|\ :\
x_{1},x_{2}>0\big{\\}}.$ (30)

Written as a matrix in the basis $\\{|e_{\rm S}\rangle,|e_{\rm B}\rangle\\}$,
we obtain from (27)

$\eta=\frac{1}{2}\begin{pmatrix}(x_{1}+x_{2})\,a_{2}/a_{1}&x_{1}-x_{2}\\\
x_{1}-x_{2}&(x_{1}+x_{2})\,a_{1}/a_{2}\end{pmatrix}.$ (31)

This is diagonal exactly when $x_{1}=x_{2}$. As we will see in Section A.1,
this is equivalent to $\eta$ being the restriction to ${\mathcal{H}}_{1}$ of a
product metric $\Lambda_{\rm S}\otimes\Lambda_{\rm B}$ on ${\mathcal{H}}$. And
as we discuss below after (39), this is also equivalent to the reduced system
state $\bar{\rho}_{H}(t)$ in (38) being a positive operator.

###  3.2 Reduced non-Hermitian system dynamics

Fix an $\eta\in\mathcal{M}_{+}$ and take an initial state of the form

$|\psi(0)\rangle=A|e_{\rm S}\rangle+B|e_{\rm B}\rangle$ (32)

for some $A,B\in\mathbb{C}$ normalized to have $\|\psi(0)\|^{2}_{\eta}=1$,
that is,

$1=\Big{(}\frac{x_{1}+x_{2}}{2}\Big{)}\Big{(}\frac{a_{2}}{a_{1}}|A|^{2}+\frac{a_{1}}{a_{2}}|B|^{2}\Big{)}+(x_{1}-x_{2}){\rm
Re}(AB^{*}).$ (33)

The dynamics is given by

$|\psi(t)\rangle=e^{-itH}|\psi(0)\rangle=e^{-it\nu}A(t)|e_{\rm
S}\rangle+e^{-it\nu}B(t)|e_{\rm B}\rangle,$ (34)

where

$\displaystyle A(t)$ $\displaystyle=A\cos(\omega
t)-iB\frac{a_{1}}{a_{2}}\sin(\omega t),$ (35) $\displaystyle B(t)$
$\displaystyle=B\cos(\omega t)-iA\frac{a_{2}}{a_{1}}\sin(\omega t).$

The normalization

$\|\psi(t)\|^{2}_{\eta}=x_{1}\big{|}\langle
v_{+}^{*}|\psi(t)\rangle\big{|}^{2}+x_{2}\big{|}\langle
v_{-}^{*}|\psi(t)\rangle\big{|}^{2}=1$ (36)

holds for all $t$, as $e^{-itH}$ acts unitarily on ${\mathcal{H}}_{1}$
equipped with the inner product $\langle\cdot|\cdot\rangle_{\eta}$. The
relation (36) is the same as (33) with $A$ and $B$ replaced by $A(t)$ and
$B(t)$.

We now introduce the reduction of the system to the single oscillator
($a^{\dagger}$, $a$). The average of a system observable $O_{\rm S}$
(observable of the single oscillator) in the state $|\psi(t)\rangle$ given by
(34) evolves according to

$\langle\psi(t)|\eta O_{\rm S}|\psi(t)\rangle={\rm tr}_{\rm
S}\big{(}\bar{\rho}_{H}(t)O_{\rm S}\big{)},$ (37)

where the reduced system state is

$\bar{\rho}_{H}(t)={\rm tr}_{\rm B}\,\rho_{H}(t)={\rm tr}_{\rm
B}\big{(}|\psi(t)\rangle\langle\psi(t)|\eta\big{)}.$ (38)

For the partial trace we have the identities ${\rm tr}_{\rm B}|e_{\rm
S}\rangle\langle e_{\rm S}|=|1_{\rm S}\rangle\langle 1_{\rm S}|$, ${\rm
tr}_{\rm B}|e_{\rm B}\rangle\langle e_{\rm B}|=|0_{\rm S}\rangle\langle 0_{\rm
S}|$ and ${\rm tr}_{\rm B}|e_{\rm S}\rangle\langle e_{\rm B}|=0={\rm tr}_{\rm
B}|e_{\rm B}\rangle\langle e_{\rm S}|$. Using (34) and $\eta$ of the form
(30), we obtain after a calculation

$\displaystyle\bar{\rho}_{H}(t)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{x_{1}+x_{2}}{2}\frac{a_{1}}{a_{2}}|B(t)|^{2}+\frac{x_{1}-x_{2}}{2}A(t)B(t)^{*}\Big{)}\,|0_{\rm
S}\rangle\langle 0_{\rm S}|$ (39)
$\displaystyle+\Big{(}\frac{x_{1}+x_{2}}{2}\frac{a_{2}}{a_{1}}|A(t)|^{2}+\frac{x_{1}-x_{2}}{2}A(t)^{*}B(t)\Big{)}\,|1_{\rm
S}\rangle\langle 1_{\rm S}|.$

This matrix is diagonal in the basis $\\{|0_{\rm S}\rangle,|1_{\rm
S}\rangle\\}$ and the two diagonal entries are its eigenvalues. One checks
directly that ${\rm tr}_{\rm S}(\bar{\rho}_{H}(t))=1$ (the sum of the diagonal
elements equals $\|\psi(t)\|_{\eta}^{2}=1$). However, the eigenvalues of
$\bar{\rho}_{H}(t)$ are complex, in general, unless the metric is chosen to
satisfy $x_{1}=x_{2}$. Indeed, the imaginary part of the first eigenvalue is
$\frac{x_{1}-x_{2}}{2}{\rm Im}(A(t)B(t)^{*})$. If $A,B$, the coefficients in
the initial state (32), are real then this quantity becomes666Use the
equations (35) and write $B$ as a function of $A$ according to the
normalization condition (33).
$-(x_{1}-x_{2})(\frac{1}{x_{1}+x_{2}}-\frac{a_{2}}{a_{1}}A^{2}-\frac{x_{1}-x_{2}}{x_{1}+x_{2}}AB)\cos(\omega
t)\sin(\omega t)$. Unless $x_{1}=x_{2}$ or the initial condition satisfies
$(x_{1}+x_{2})\frac{a_{2}}{a_{1}}A^{2}+(x_{1}-x_{2})AB=1$, the eigenvalues of
$\bar{\rho}_{H}(t)$ will not be real except at the discrete set of times $t$
when $\sin(\omega t)\cos(\omega t)=0$.

We require $\bar{\rho}_{H}(t)$ to be a density matrix (and in particular to
have non-negative eigenvalues) for all times. To do so with a metric that does
not depend on the initial conditions we therefore must choose $x_{1}=x_{2}$.
We thus take

$x_{1}=x_{2}=x>0$

for the remainder of the paper. In the basis $\\{|e_{\rm S}\rangle,|e_{\rm
B}\rangle\\}$ the metric $\eta$ is diagonal,

$\eta=x\begin{pmatrix}a_{2}/a_{1}&0\\\ 0&a_{1}/a_{2}\end{pmatrix},$ (40)

see (31). As explained after (31) above, this is equivalent to $\eta$ being of
product form. With this choice, $\bar{\rho}_{H}(t)$ given by (39) is
${\dagger}$-Hermitian. According to (39) and (35) we have

$\bar{\rho}_{H}(t)=p(t)\,|0_{\rm S}\rangle\langle 0_{\rm S}|+(1-p(t))\,|1_{\rm
S}\rangle\langle 1_{\rm S}|,$ (41)

where

$\displaystyle p(t)$ $\displaystyle=$ $\displaystyle
x\frac{a_{1}}{a_{2}}|B(t)|^{2}$ (42) $\displaystyle=$ $\displaystyle
x\Big{(}\frac{a_{1}}{a_{2}}|B|^{2}\cos^{2}(\omega
t)+\frac{a_{2}}{a_{1}}|A|^{2}\sin^{2}(\omega t)-2\sin(\omega t)\cos(\omega
t){\rm Im}(A^{*}B)\Big{)}$ $\displaystyle=$
$\displaystyle\frac{1}{2}+\Big{(}\frac{1}{2}-x\frac{a_{2}}{a_{1}}|A|^{2}\Big{)}\cos(2\omega
t)-x\sin(2\omega t){\rm Im}(A^{*}B).$

In the last step, we used the normalization condition (33), resulting in
$\frac{a_{1}}{a_{2}}|B|^{2}=\frac{1}{x}-\frac{a_{2}}{a_{1}}|A|^{2}$, and the
trigonometric identities $\sin(\omega t)\cos(\omega t)=\frac{1}{2}\sin(2\omega
t)$, $\cos^{2}(\omega t)=\frac{1}{2}(1+\cos(2\omega t))$ and $\sin^{2}(\omega
t)=\frac{1}{2}(1-\cos(2\omega t))$. In view of (42) it is natural to introduce
the parameter

$\alpha\equiv x\frac{a_{2}}{a_{1}}|A|^{2}\in[0,1].$ (43)

Equations (41) and (42) show the following.

Properties of $p(t)$:

  * 1.

$p(t)$ and $\bar{\rho}_{H}(t)$ depend on time unless $\alpha=\tfrac{1}{2}$ and
$A^{*}B\in\mathbb{R}$, in which case $p(t)=\tfrac{1}{2}$ and
$\bar{\rho}_{H}(t)=\tfrac{1}{2}{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$.

  * 2.

Otherwise $p(t)$ and $\bar{\rho}_{H}(t)$ are periodic in time, with period
$\pi/\omega$, and the mean value of $p(t)$ is

$p_{0}=\frac{\omega}{\pi}\int_{0}^{\pi/\omega}p(t)dt=\frac{1}{2}.$ (44)

###  3.3 Reduced Hermitian system dynamics

Next, we turn our attention to the density matrix of the Hermitian system,
which according to (14) is

$\rho_{h_{W}}(t)=S\rho_{H}(t)S^{-1}=W\sqrt{\eta}\rho_{H}(t)\tfrac{1}{\sqrt{\eta}}W^{\dagger}=W\sqrt{\eta}\,|\psi(t)\rangle\langle\psi(t)|\sqrt{\eta}\,W^{\dagger}.$
(45)

We keep $W$ in the notation $h_{W}$ to highlight that the choice of $h$
depends on $W$, see (8). Again choosing a metric $\eta$ of the form (30) with
$x_{1}=x_{2}=x>0$, we use (34) to obtain

$\sqrt{\eta}|\psi(t)\rangle=e^{-it\nu}\gamma(t)|e_{\rm
S}\rangle+e^{-it\nu}\delta(t)|e_{\rm B}\rangle,$ (46)

where

$\gamma(t)=\sqrt{x\frac{a_{2}}{a_{1}}}A(t),\qquad\delta(t)=\sqrt{x\frac{a_{1}}{a_{2}}}B(t).$
(47)

We then obtain

$\displaystyle\sqrt{\eta}\,|\psi(t)\rangle\langle\psi(t)|\sqrt{\eta}=\begin{pmatrix}|\gamma(t)|^{2}&\gamma(t)\delta(t)^{*}\\\
\gamma(t)^{*}\delta(t)&|\delta(t)|^{2}\end{pmatrix},$ (48)

written in matrix form in the ordered basis $\\{|e_{\rm S}\rangle,|e_{\rm
B}\rangle\\}$ of ${\mathcal{H}}_{1}$. Next, we take a general (time-
independent) unitary on ${\mathcal{H}}_{1}$, expressed in the same basis as

$W=\begin{pmatrix}a&b\\\ c&d\end{pmatrix},\qquad
ac^{*}+bd^{*}=0,\quad|a|^{2}+|b|^{2}=1=|c|^{2}+|d|^{2}.$ (49)

Using (45), (48), (49) and writing momentarily $\delta,\gamma$ for
$\delta(t),\gamma(t)$, we get

$\rho_{h_{W}}(t)=\begin{pmatrix}|a\gamma+b\delta|^{2}\quad&ac^{*}|\gamma|^{2}+bc^{*}\gamma^{*}\delta+ad^{*}\gamma\delta^{*}+bd^{*}|\delta|^{2}\\\
a^{*}c|\gamma|^{2}+b^{*}c\gamma\delta^{*}+a^{*}d\gamma^{*}\delta+b^{*}d|\delta|^{2}&|c\gamma+d\delta|^{2}\end{pmatrix}.$
(50)

We observe that $\rho_{h_{W}}(t)$ is periodic in time, with period
$\pi/\omega$. This follows from (47) together with (35), since
$\cos^{2}(\omega t)$, $\sin^{2}(\omega t)$ and $\sin(\omega t)\cos(\omega t)$
all have period $\pi/\omega$.

Next we calculate the reduced density matrix $\bar{\rho}_{h_{W}}(t)$ of ${\rm
S}$ by taking the partial trace of $\rho_{h_{W}}(t)$ over ${\rm B}$,

$\bar{\rho}_{h_{W}}(t)=q(t)\,|0_{\rm S}\rangle\langle 0_{\rm
S}|+\big{(}1-q(t)\big{)}\,|1_{\rm S}\rangle\langle 1_{\rm S}|,$ (51)

where $q(t)=\big{|}c\gamma(t)+d\delta(t)\big{|}^{2}$ and $\gamma(t),\delta(t)$
are given in (47). We get the expression

$q(t)=x\Big{|}c\sqrt{\frac{a_{2}}{a_{1}}}A(t)+d\sqrt{\frac{a_{1}}{a_{2}}}B(t)\Big{|}^{2},$
(52)

where $A(t)$ and $B(t)$ are given in (35). Expanding the square and using
$|d|^{2}=1-|c|^{2}$ as well as the normalization (33) (which is valid for
$A,B$ replaced by $A(t),B(t)$ at any time), we obtain

$\displaystyle q(t)$ $\displaystyle=$
$\displaystyle|c|^{2}+(1-2|c|^{2})x\frac{a_{1}}{a_{2}}|B(t)|^{2}+2x{\rm
Re}\big{(}cd^{*}A(t)B(t)^{*}\big{)}$ (53) $\displaystyle=$
$\displaystyle|c|^{2}+(1-2|c|^{2})p(t)+2x{\rm
Re}\big{(}cd^{*}A(t)B(t)^{*}\big{)},$

where $p(t)$ is the population of $\bar{\rho}_{H}(t)$ evaluated above in (42).
Expanding the real part term in (53) using (35), we arrive at

$\displaystyle q(t)$ $\displaystyle=$ $\displaystyle|c|^{2}+(1-2|c|^{2})p(t)$
(54) $\displaystyle-\sin(2\omega
t)\big{(}1-2x\frac{a_{2}}{a_{1}}|A|^{2}\big{)}{\rm Im}(cd^{*})-2x\cos(2\omega
t){\rm Im}(AB^{*}){\rm Im}(cd^{*})$ $\displaystyle+2x{\rm Re}(AB^{*}){\rm
Re}(cd^{*}).$

We take into account (42) to rewrite

$\displaystyle q(t)$ $\displaystyle=$ $\displaystyle\tfrac{1}{2}+2x{\rm
Re}(AB^{*}){\rm Re}(cd^{*})$ (55) $\displaystyle+\cos(2\omega
t)\Big{[}(\tfrac{1}{2}-\alpha)(1-2|c|^{2})-2x{\rm Im}(AB^{*}){\rm
Im}(cd^{*})\Big{]}$ $\displaystyle-\sin(2\omega t)\Big{[}(1-2\alpha){\rm
Im}(cd^{*})-x(1-2|c|^{2}){\rm Im}(AB^{*})\Big{]},$

where we recall that $\alpha$ is given in (43). Using the explicit form (55)
of $q(t)$, we obtain the following information.

Properties of $q(t)$:

  * 1.

$q(t)$ and $\bar{\rho}_{h_{W}}(t)$ depend on time unless both factors of the
cosine and sine terms in (55) vanish. Thus $q(t)$ is time-independent and
hence equal to

$q_{0}=\frac{\omega}{\pi}\int_{0}^{\pi/\omega}q(t)dt=\frac{1}{2}+2x{\rm
Re}(AB^{*}){\rm Re}(cd^{*})$ (56)

if and only if:

    * $\bullet$

$\alpha=1/2$, $|c|=1/\sqrt{2}$ and ${\rm Im}(AB^{*}){\rm Im}(cd^{*})=0$, or

    * $\bullet$

$\alpha=1/2$, $|c|\neq 1/\sqrt{2}$ and ${\rm Im}(AB^{*})=0$, or

    * $\bullet$

$\alpha\neq 1/2$, $|c|=1/\sqrt{2}$ and ${\rm Im}(cd^{*})=0$, or

    * $\bullet$

$\alpha\neq 1/2$, $|c|\neq 1/\sqrt{2}$ and ${\rm
Im}(A^{*}B)=\frac{1}{2x}|1-2\alpha|$ and ${\rm
Im}(cd^{*})=\frac{1}{2}(1-2|c|^{2}){\rm sgn}(1-2\alpha)$, where ${\rm
sgn}(x)=|x|/x$.

  * 2.

Otherwise $q(t)$ and $\bar{\rho}_{h_{W}}(t)$ are periodic in time, with period
$\pi/\omega$, and the mean value of $q(t)$ is $q_{0}$ from (56).

Generic initial states and unitaries.  As $x>0$, the average $q_{0}$ equals
$\frac{1}{2}$ exactly when ${\rm Re}(AB^{*}){\rm Re}(cd^{*})=0$. This is a
condition on the initial state (via $A,B$) and the unitary $W$ (via $c,d$). We
call the initial state and the unitary generic, respectively, when

$\mbox{$AB^{*}\not\in\mathbb{R}$ \quad and\quad$cd^{*}\not\in\mathbb{R}$}.$
(57)

In other words, for generic initial states and unitaries, the average $q_{0}$
of the population of $\bar{\rho}_{h_{W}}(t)$ differs from the average
$p_{0}=\frac{1}{2}$ of the population of $\bar{\rho}_{H}(t)$. As we show in
the next section, this deviation from the value $\frac{1}{2}$ causes the
entropy of $\bar{\rho}_{h_{W}}(t)$ to oscillate with exactly half the
frequency of the entropy of $\bar{\rho}_{H}(t)$.

##  4 Entropy

Recall that the states $\bar{\rho}_{H}(t)$ and $\bar{\rho}_{h_{W}}(t)$ are
given in (41) and (51), with associated populations $p(t)$, $q(t)$ evaluated
in (42) and (55). Their von Neumann entropy is given by

$\displaystyle\mathcal{E}(\bar{\rho}_{H}(t))$ $\displaystyle=$
$\displaystyle-p(t)\ln p(t)-(1-p(t))\ln(1-p(t)),$
$\displaystyle\mathcal{E}(\bar{\rho}_{h_{W}}(t))$ $\displaystyle=$
$\displaystyle-q(t)\ln q(t)-(1-q(t))\ln(1-q(t)).$ (58)

We show in Appendix A.2 that

$\mathcal{E}(\bar{\rho}_{H}(t))=\mathcal{E}(\bar{\rho}_{h_{W}}(t))\mbox{\ for
all $t\geq 0$ and initial conditions ($A,B$)\,}\quad\Longleftrightarrow\quad
cd=0.$

If ${\rm Re}(cd^{*})\neq 0$, then according to (44) and (56) the averages
$p_{0}$ and $q_{0}$ around which the populations $p(t)$ and $q(t)$ oscillate
are different for all generic initial conditions, i.e., all coefficients $A,B$
satisfying ${\rm Re}(AB^{*})\neq 0$. This translates into a modification of
the period of the entropy of $\bar{\rho}_{h_{W}}(t)$ as a function of time $t$
relative to that of $\bar{\rho}_{H}(t)$, as we explain now.

###  4.1 Period doubling of von Neumann entropy in the non-Hermitian versus
the Hermitian system

Consider the function

$\mathcal{E}(Q)=-Q\ln Q-(1-Q)\ln(1-Q),\quad Q\in[0,1].$

Suppose now that $Q=Q(t)$ depends periodically on time and has average
$Q_{0}$,

$Q(t)=Q_{0}+\Delta(t)\in[0,1],$ (59)

with $\Delta(t)$ having period $\pi/\omega$ and zero average. This setup
incorporates both cases $p(t)$ and $q(t)$ in one. As Figure 1 illustrates, if
$Q_{0}=\frac{1}{2}$, which is the value where $\mathcal{E}(Q)$ takes its
maximum, then as $Q(t)$ moves over one period, the entropy $\mathcal{E}(Q(t))$
moves over two periods.

Figure 1: Parameters $Q_{0}=0.5$, $\Delta\equiv\max_{t}\Delta(t)=0.2$.
According to (59), $Q(t)$ starts at $Q_{m}=0.3$ (at a time we take to be
$t=0$) and moves to $Q_{M}=0.7$ at time $\omega t=\pi/2$, and then back to
$Q_{m}$ at time $\omega t=\pi$ (left panel), so the value of the entropy
$\mathcal{E}(Q(t))$ evolves through two periods (right panel). In each period,
the entropy has two local minima (counting minima at the endpoints of the
considered intervals once).

On the other hand, if $Q_{0}\neq\frac{1}{2}$, then the period of the entropy
$\mathcal{E}(Q(t))$ is not doubled relative to that of $Q(t)$, as Figures 2
and 3 show.

Figure 2: Parameters $Q_{0}=0.6$, $\Delta=\max_{t}\Delta(t)=0.2$. In this case
$Q(t)$ starts at $Q_{m}=0.4$ when $t=0$ (upon a possible shift of the time
axis) and moves to $Q_{M}=0.8$ at time $\omega t=\pi/2$ and back to $Q_{m}$ at
time $\omega t=\pi$ (left panel). The value of the entropy $\mathcal{E}(Q(t))$
evolves through one single period (right panel). In each period, the entropy
has two local minima (counting minima at the endpoints once).

Figure 3: Parameters $Q_{0}=0.8$, $\Delta=\max_{t}\Delta(t)=0.2$. In this case
$Q(t)$ starts at $Q_{m}=0.6$ when $t=0$ (after possibly shift the time axis)
and moves to $Q_{M}=1.0$ at time $\omega t=\pi/2$ and back to $Q_{m}$ at time
$\omega t=\pi$ (left panel). The value of the entropy $\mathcal{E}(Q(t))$
evolves through one single period (right panel). Since $0.5$ is not in the
interval $(Q_{m},Q_{M})$, the graph of the entropy has only one local minimum
in each period, instead of two when the interval contains the value $0.5$ for
$Q$.

We draw the following conclusions:

  * $\bullet$

The period of the von Neumann entropy of the non-Hermitian system
$\mathcal{E}(\bar{\rho}_{H}(t))$ is $\frac{1}{2}\pi/\omega$, regardless of the
initial condition (except for the stationary state).

  * $\bullet$

Regardless of the metric (parameter $x$), the period of the von Neumann
entropy of the Hermitian system $\mathcal{E}(\bar{\rho}_{h_{W}}(t))$ is:

    * $\circ$

$\pi/\omega$, provided ${\rm Re}(AB^{*}){\rm Re}(cd^{*})\neq 0$ (generic
case),

    * $\circ$

$\frac{1}{2}\pi/\omega$, provided ${\rm Re}(AB^{*}){\rm Re}(cd^{*})=0$
(special case).

This means that for generic initial conditions (meaning ${\rm Re}(AB^{*})\neq
0$) and generic choices of the unitary $W$ (meaning ${\rm Re}(cd^{*})\neq 0$),
the period of the von Neumann entropy of the Hermitian system is double that
of the non-Hermitian system. That is, the entropy of the non-Hermitian system
oscillates faster. This is so even though the populations in both cases have
the same frequency $\pi/\omega$. The change of the period is due to the shift
of the average in the population induced by $W$, as given in (56).

###  4.2 Numerical illustration of the period doubling

We plot the populations and entropies for parameters in the regime

$\displaystyle x=x_{1}=x_{2}$ $\displaystyle>$ $\displaystyle
0\qquad\mbox{(metric $\eta$, \eqref{metricmatrix})}$ (60) $\displaystyle c,d$
$\displaystyle\geq$ $\displaystyle 0\qquad\mbox{(unitary $W$, \eqref{26})}$
(61) $\displaystyle A,B$ $\displaystyle\geq$ $\displaystyle
0\qquad\mbox{(initial state $|\psi(0)\rangle$, \eqref{ic})}$ (62)

According to (61) and the unitarity of $W$, we have $d=\sqrt{1-c^{2}}$.
Moreover, $x^{2}A^{2}B^{2}=\alpha(1-\alpha)$, where
$\alpha=x\frac{a_{2}}{a_{1}}A^{2}\in[0,1]$. The population $q(t)$ of the
Hermitian system reduced density matrix, given in (55), then becomes

$q(t)=q_{0}+\Delta\cos(2\omega t)$ (63)

with

$\displaystyle q_{0}$ $\displaystyle=$
$\displaystyle\tfrac{1}{2}+2\sqrt{c^{2}(1-c^{2})}\sqrt{\alpha(1-\alpha)},$
(64) $\displaystyle\Delta$ $\displaystyle=$
$\displaystyle\tfrac{1}{2}(1-2c^{2})(1-2\alpha).$ (65)

Here, $\alpha=x\frac{a_{2}}{a_{1}}A^{2}\in[0,1]$ and $c\in[0,1]$ can be chosen
freely. The population $p(t)$ of the non-Hermitian system, in (42), is simply
the expression (63) with $c=0$.

Note that the change $\alpha\mapsto 1-\alpha$ leaves $q_{0}$ invariant and
flips the sign of $\Delta$. It then suffices to plot graphs for
$\alpha\in[0,1/2]$. For $\alpha=1/2$ we get $\Delta=0$, which gives a
stationary state (for all $c$). The same invariance of $q_{0}$ and sign flip
of $\Delta$ is induced by $c^{2}\mapsto 1-c^{2}$.

In Figure 4, we compare the von Neumann entropies of the two density matrices
$\bar{\rho}_{H}(t)$ and $\bar{\rho}_{h_{W}}(t)$, directly seeing the doubling
of the period.

Figure 4: Comparing the entropies of $\bar{\rho}_{h_{W}}(t)$ and of
$\bar{\rho}_{H}(t)$ for the values $c=0.5$ and $\alpha=0,0.15,0.3,0.45$. The
doubling of the period for $\alpha\neq 0$ is manifest. The oscillations
decrease as $\alpha$ approaches $0.5$, which gives the stationary state.

##  5 Entanglement of system and bath oscillators

The total Hilbert space ${\mathcal{H}}={\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}$ in (17) is bipartite, one part being the
singled-out oscillator (system), the other being the remaining $N$ oscillators
(bath). We say that a nonzero vector $|\psi\rangle\in{\mathcal{H}}$ is of
product form, or disentangled, if $|\psi\rangle=|\psi_{\rm
S}\rangle\otimes|\psi_{\rm B}\rangle$ for some $|\psi_{\rm
S}\rangle\in{\mathcal{H}}_{\rm S}$ and some $|\psi_{\rm
B}\rangle\in{\mathcal{H}}_{\rm B}$. We call a nonzero
$|\psi\rangle\in{\mathcal{H}}$ entangled if it is not of product form. The
notion of being entangled or not does not depend on the metric determining the
inner product of ${\mathcal{H}}$. Nevertheless, the physical interpretation of
entanglement in terms of independence of the subsystems ${\rm S}$ and ${\rm
B}$ does depend on the metric. The physical manifestation of disentangled
states is the independence of the two subsystems ${\rm S}$ and ${\rm B}$.
Namely, if the inner product of ${\mathcal{H}}$ is given by a metric of the
form $\eta=\Lambda_{\rm S}\otimes\Lambda_{\rm B}$ (a particular example being
${\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}_{\rm S}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}_{\rm B}$), then
measurement outcomes of observables on either of the subsystems are
independent random variables. This follows because expectation values of
observables $O_{\rm S}\otimes O_{\rm B}$ in a state $|\psi_{\rm
S}\rangle\otimes|\psi_{\rm B}\rangle$ split into products,

$\langle\psi_{\rm S}\otimes\psi_{\rm B}|\eta(O_{\rm S}\otimes O_{\rm
B})\psi_{\rm S}\otimes\psi_{\rm B}\rangle=\langle\psi_{\rm S}|\Lambda_{\rm
S}O_{\rm S}\psi_{\rm S}\rangle\langle\psi_{\rm B}|\Lambda_{\rm B}O_{\rm
B}\psi_{\rm B}\rangle.$

However, those random variables become dependent (correlated) if $\eta$ is not
of product form, because then their average will not split into a product of a
system term times a bath term.

In the model defined in Section 3, the metrics $\eta$ we consider are
restrictions to the subspace ${\mathcal{H}}_{1}$ of product metrics
$\Lambda_{\rm S}\otimes\Lambda_{\rm B}$ of ${\mathcal{H}}$ (c.f. (31) and
Section A.1). The physical meaning of ${\rm S}{\rm B}$ entanglement in terms
of subsystem independence does therefore not depend on the choice of $\eta$
within this class. In other words, measurement outcomes of system and bath
observables in $|\psi\rangle\in{\mathcal{H}}_{1}$ are independent or
dependent, according to whether $|\psi\rangle\in{\mathcal{H}}_{1}$ is
disentangled or not, regardless of the choice of $\eta$. It is then sensible
to investigate the SB entanglement in pure states belonging to the subspace
${\mathcal{H}}_{1}$ for all $\eta$.

Any vector $|\psi\rangle\in{\mathcal{H}}_{1}$ is of the form

$|\psi\rangle=A|e_{\rm S}\rangle+B|e_{\rm B}\rangle,\qquad A,B\in\mathbb{C}.$
(66)

In accordance with (22), we may write $|e_{\rm
S}\rangle=|10\rangle=|1\rangle\otimes|0\rangle\in{\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}$ and $|e_{\rm
B}\rangle=|01\rangle=|0\rangle\otimes|1\rangle\in{\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}$. Explicitly,
$|0\rangle,|1\rangle\in{\mathcal{H}}_{\rm S}$ are the ground state and first
excited state of the system oscillator, and
$|0\rangle,|1\rangle\in{\mathcal{H}}_{\rm B}$ are the ground state of all the
$N$ bath oscillators and the distributed excitation state
$\frac{1}{\sqrt{N}}\sum_{n=1}^{N}|1_{n}\rangle$, respectively; see (22).

We now show that

$|\psi\rangle\mbox{\ of the form \eqref{2-73} is
disentangled}\quad\Longleftrightarrow\quad AB=0.$ (67)

To see that the implication $\Rightarrow$ in (67) holds, let $\rho_{\rm
S}\equiv{\rm tr}_{\rm B}|\psi\rangle\langle\psi|$ (partial trace over ${\rm
B}$). On the one hand, (66) gives $\rho_{\rm S}=|A|^{2}|1\rangle\langle
1|+|B|^{2}|0\rangle\langle 0|$. On the other hand, if $|\psi\rangle$ is
disentangled, then $\rho_{\rm S}$ must have rank one, since $\rho_{\rm S}={\rm
tr}_{\rm B}(|\psi_{\rm S}\rangle\langle\psi_{\rm S}|\otimes|\psi_{\rm
B}\rangle\langle\psi_{\rm B}|)=|\psi_{\rm S}\rangle\langle\psi_{\rm S}|\
\|\psi_{\rm B}\|^{2}$. This forces either $A=0$ or $B=0$. Conversely, to see
the implication $\Leftarrow$ in (67), we note that if either of $A$ or $B$
vanish, then $|\psi\rangle$ is is proportional to $|e_{\rm B}\rangle$ or
$|e_{\rm S}\rangle$, so $|\psi\rangle$ is disentangled.

Entanglement in the non-Hermitian system.  An initial state
$|\psi(0)\rangle=A|e_{\rm S}\rangle+B|e_{\rm B}\rangle$ evolves into
$|\psi(t)\rangle=e^{-itH}|\psi(0)\rangle=A(t)|e_{\rm S}\rangle+B(t)|e_{\rm
B}\rangle$, where the time dependent coefficients are given in (35). According
to (67), $|\psi(t)\rangle$ is disentangled exactly if $A(t)B(t)=0$. Let us
first analyze the condition $A(t)=0$. This equality is equivalent to the two
equations

$\displaystyle({\rm Re}A)\cos(\omega t)+({\rm
Im}B)\frac{a_{1}}{a_{2}}\sin(\omega t)$ $\displaystyle=0,$ (68)
$\displaystyle({\rm Im}A)\cos(\omega t)-({\rm
Re}B)\frac{a_{1}}{a_{2}}\sin(\omega t)$ $\displaystyle=0.$

The condition $B(t)=0$ is the same as (68) but with $A\leftrightarrow B$
swapped and $a_{1}\leftrightarrow a_{2}$ swapped.

Suppose $|\psi(0)\rangle$ is disentangled, so $AB=0$. Then exactly one of $A$
or $B$ vanish and the equations (68) are satisfied for $\omega
t\in\pi\mathbb{Z}$ (if $A=0$) or for $\omega t=\frac{\pi}{2}(2\mathbb{Z}+1)$
(if $B=0$). We conclude that $|\psi(t)\rangle$ is entangled except
periodically at discrete moments in time where it is disentangled.

On the other hand, if $|\psi(0)\rangle$ is entangled, then both $A$ and $B$ do
not vanish. If $A$ and $B$ are both real or both purely imaginary, then (68)
is not satisfied for any $t$. For all other $A$ and $B$ (68) is satisfied for
discrete, periodically repeating values of $t$.

We conclude:

  * (a)

If the initial state $|\psi(0)\rangle$ is entangled and both $A,B$ are either
purely real or purely imaginary, then $|\psi(t)\rangle$ is entangled for all
times $t$.

  * (b)

With the exception of case (a) and regardless of the entanglement in the
initial state $|\psi(0)\rangle$, the state $|\psi(t)\rangle$ is entangled
except at periodically repeating instants.

Entanglement in the Hermitian systems.  The Hermitian system pure state vector
is given by (see also (12))

$|\phi(t)\rangle=S|\psi(t)\rangle=W\sqrt{\eta}|\psi(t)\rangle=\widetilde{A}(t)|e_{\rm
S}\rangle+\widetilde{B}(t)|e_{\rm B}\rangle,$

a normalized vector in ${\mathcal{H}}$ (with the original inner product),
where

$\begin{pmatrix}\widetilde{A}(t)\\\
\widetilde{B}(t)\end{pmatrix}=T\begin{pmatrix}A(t)\\\ B(t)\end{pmatrix},\qquad
T=\sqrt{x\frac{a_{2}}{a_{1}}}\begin{pmatrix}a\ &ba_{1}/a_{2}\\\
c&da_{1}/a_{2}\end{pmatrix},\qquad\det T=x\det W\neq 0$ (69)

satisfy $|\widetilde{A}(t)|^{2}+|\widetilde{B}(t)|^{2}=1$; c.f. (47) and (49).
It follows from (67) that $|\phi(t)\rangle$ is entangled if and only if
$\widetilde{A}(t)\widetilde{B}(t)=0$. An analysis of the latter equality along
the lines of that carried out after (68) shows that $|\phi(t)\rangle$ is
entangled except at isolated, periodically reoccurring instants in time, just
like the state of the non-Hermitian system.

Effect of choice of $W$ on entanglement.  Given a state $|\psi\rangle=A|e_{\rm
S}\rangle+B|e_{\rm B}\rangle$ of the non-Hermitian system, the associated
Hermitian system state is $|\phi\rangle=W\sqrt{\eta}|\psi\rangle$. For the
choice $W={\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$ we have

$|\phi\rangle=\sqrt{\eta}|\psi\rangle=\sqrt{xa_{2}/a_{1}}A|e_{\rm
S}\rangle+\sqrt{xa_{1}/a_{2}}B|e_{\rm B}\rangle.$ (70)

Hence for $W={\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}$, $|\phi\rangle$ is entangled if and only
if $|\psi\rangle$ is entangled (recall (67)). The metric $\eta$ does not alter
the property of being entangled. Choosing a different $W$ to build
$|\phi\rangle$ from $|\psi\rangle$, however, changes this. It is not hard to
see that the unitaries $W$ that map every product state (that is $|e_{\rm
S}\rangle$ and $|e_{\rm B}\rangle$) into another product state are exactly the
diagonal and the off-diagonal $W$. Furthermore, given an entangled
$\sqrt{\eta}|\psi\rangle$ as in (70), one can always find unitaries $W$ such
that $W\sqrt{\eta}|\psi\rangle$ is not entangled. Those $W$ are precisely the
ones with $|a|=|d|=\sqrt{a_{1}/(xa_{2})}|B|$ and
$|b|=|c|=\sqrt{xa_{2}/a_{1}}|A|$.777This follows from the characterization
(67) of product states and the normalization
$\|\sqrt{\eta}|\psi\rangle\|_{\mathcal{H}}=1$.

We now examine the effect of $W$ on the time-averaged density matrix

$\langle\rho\rangle=\frac{\omega}{\pi}\int_{0}^{\pi/\omega}|\phi(t)\rangle\langle\phi(t)|\,dt$
(71)

where we integrate $|\phi(t)\rangle\langle\phi(t)|=\rho_{h_{W}}(t)$ over one
period, see (50). A direct calculation yields

$\langle\rho\rangle=\begin{pmatrix}q_{0}&z\\\ z^{*}&1-q_{0}\end{pmatrix},\quad
z=(bc^{*}+ad^{*})x{\rm Re}(AB^{*}),$ (72)

with $q_{0}=\tfrac{1}{2}+2x{\rm Re}(AB^{*}){\rm Re}(cd^{*})$, cf. (56). The
density matrix (72) is written in the basis $\\{|e_{\rm
S}\rangle\equiv|10\rangle,|e_{\rm B}\rangle\equiv|01\rangle\\}$ of
${\mathcal{H}}_{1}$, using the same notation as after (66). We view
${\mathcal{H}}_{1}$ as a subspace of the four dimensional space of two qubits,
spanned by the vectors $\\{|00\rangle,|01\rangle,|10\rangle,|11\rangle\\}$. In
this basis, the density matrix (72) takes the form

$\langle\rho\rangle=\begin{pmatrix}0&0&0&0\\\ 0&1-q_{0}&z^{*}&0\\\
0&z&q_{0}&0\\\ 0&0&0&0\end{pmatrix}.$ (73)

We calculate the concurrence [30, 19] of $\langle\rho\rangle$ to be888In the
present case, the square of the concurrence is the difference between the two
non-zero eigenvalues of the squared matrix $\langle\rho\rangle^{2}$.

$\mathcal{C}\big{(}\langle\rho\rangle\big{)}=2x|{\rm Re}(AB^{*})|.$ (74)

The concurrence of any two qubit density matrix is bounded below by $0$
(separable state) and above by $1$ (maximally entangled state). Both
inequalities in

$2x|{\rm Re}(AB^{*})|\leq 2x|A||B|\leq
x\Big{(}\frac{a_{2}}{a_{1}}|A|^{2}+\frac{a_{1}}{a_{2}}|B|^{2}\Big{)}=1$

(the last equality is the normalization (33), with $x=x_{1}=x_{2}>0$) are
saturated exactly if $AB^{*}\in\mathbb{R}$ and $|A|=\frac{a_{1}}{a_{2}}|B|$.

We conclude that the concurrence of $\langle\rho\rangle$ is the same for all
choices $W$, so it only depends on the initial state. The state
$\langle\rho\rangle$ is separable if and only if ${\rm Re}(AB^{*})=0$, and is
maximally entangled if and only if ${\rm Im}(AB^{*})=0$ and
$|A|=\frac{a_{1}}{a_{2}}|B|$.

##  6 Conclusion

The Dyson map assigns to a given quasi-Hermitian quantum system an associated
Hermitian system in a non-unique way. We quantify the non-uniqueness by means
of a metric operator $\eta$ and a unitary map $W$. The physical properties of
the Hermitian systems depend on the choice of $W$, and it is not obvious how
to capture the dynamics of the original quasi-Hermitian system in its
Hermitian counterparts – unless there happens to be some universality
throughout the Hermitian family. We describe an aspect of universality for a
quasi-Hermitian open system consisting of a single oscillator coupled to a
bath of $N$ oscillators. We show that there is a unique metric operator for
which the reduced state of the system (single oscillator) is a well-defined
density matrix. Using this metric, we construct all Hermitian systems obtained
from the quasi-Hermitian one by varying $W$. We find that the entropy of the
single oscillator in the Hermitian system evolves periodically in time with
exactly double the period of the corresponding entropy of the quasi-Hermitian
system, independently of $W$. We further show that averaged over one time
period, the entanglement between the oscillator and the bath is independent of
the choice of $W$.

##  Appendix A Appendix

###  A.1 Diagonal form of $\eta$ is equivalent to product form

We have seen above in (40) that $\eta$ must be diagonal in order for the
populations of $\bar{\rho}_{H}(t)$ to be non-negative and that conversely if
$\eta$ is diagonal, then the populations of $\bar{\rho}_{H}(t)$ are positive.
As it turns out, $\eta$ being diagonal is also equivalent to $\eta$ being of
product form. More precisely, the following two statements are equivalent:

  * (1)

$\eta$ is diagonal in the basis $|e_{\rm S}\rangle,|e_{\rm B}\rangle$ of
${\mathcal{H}}_{1}$.

  * (2)

There are metrics $\Lambda_{\rm S}$ and $\Lambda_{B}$ on ${\mathcal{H}}_{\rm
S}$ and ${\mathcal{H}}_{\rm B}$, respectively, such that $\Lambda_{\rm
S}\otimes\Lambda_{\rm B}$ leaves ${\mathcal{H}}_{1}$ invariant and
$\eta=\Lambda_{\rm S}\otimes\Lambda_{\rm
B}\\!\upharpoonright_{{\mathcal{H}}_{1}}$ is the restriction of this product
to ${\mathcal{H}}_{1}$.

Given $\eta$, the $\Lambda_{\rm S}$ and $\Lambda_{\rm B}$ are not unique.

Proof of (1) $\Leftrightarrow$ (2).

Consider a bipartite Hilbert space ${\mathcal{H}}={\mathcal{H}}_{\rm
S}\otimes{\mathcal{H}}_{\rm B}$ with an orthonormal basis
$|v_{ij}\rangle=|e_{i}\rangle\otimes|f_{j}\rangle$, the $|e_{i}\rangle$ and
$|f_{j}\rangle$ being orthonormal bases of ${\mathcal{H}}_{\rm S}$ and
${\mathcal{H}}_{\rm B}$, respectively. Let ${\mathcal{H}}_{1}$ be the two-
dimensional subspace ${\mathcal{H}}_{1}={\rm
span}\\{|v_{11}\rangle,|v_{22}\rangle\\}$. In this setup, $|v_{11}\rangle$ is
identified with $|e_{\rm S}\rangle$ and $|v_{22}\rangle$ with $|e_{\rm
B}\rangle$. Let $\eta>0$ be a strictly positive operator on
${\mathcal{H}}_{1}$.

We first show (1) $\Rightarrow$ (2). Assume that $\eta$ is diagonal, that is,
$\eta=a|v_{11}\rangle\langle v_{11}|+b|v_{22}\rangle\langle v_{22}|$, where
$a,b>0$. Set

$\displaystyle\Lambda_{\rm S}$ $\displaystyle=$
$\displaystyle\alpha_{1}|e_{1}\rangle\langle
e_{1}|+\alpha_{2}|e_{2}\rangle\langle e_{2}|+\Lambda_{\rm S}^{\perp}$
$\displaystyle\Lambda_{\rm B}$ $\displaystyle=$
$\displaystyle\beta_{1}|f_{1}\rangle\langle
f_{1}|+\beta_{2}|f_{2}\rangle\langle f_{2}|+\Lambda_{\rm B}^{\perp},$

where $\alpha_{1},\alpha_{2},\beta_{1},\beta_{2}>0$ satisfy
$\alpha_{1}\beta_{1}=a$, $\alpha_{2}\beta_{2}=b$ and $\Lambda_{\rm
S}^{\perp}$, $\Lambda_{\rm B}^{\perp}$ are arbitrary positive operators on the
orthogonal complements of ${\rm span}\\{|e_{1}\rangle,|e_{2}\rangle\\}$ and
${\rm span}\\{|f_{1}\rangle,|f_{2}\rangle\\}$ in ${\mathcal{H}}_{\rm S}$ and
${\mathcal{H}}_{\rm B}$, respectively. Then $\Lambda_{\rm
S}\otimes\Lambda_{\rm B}$ is a metric on ${\mathcal{H}}$ which leaves
${\mathcal{H}}_{1}$ invariant and satisfies $(\Lambda_{\rm
S}\otimes\Lambda_{\rm B})|v_{jj}\rangle=\eta|v_{jj}\rangle$ for $j=1,2$.

Next we prove that (2) $\Rightarrow$ (1). The orthogonal projection onto
${\mathcal{H}}_{1}$ is given by

$\pi=|v_{11}\rangle\langle v_{11}|+|v_{22}\rangle\langle v_{22}|=p_{11}\otimes
q_{11}+p_{22}\otimes q_{22},$

where $p_{ij}=|e_{i}\rangle\langle e_{j}|$ and $q_{ij}=|f_{i}\rangle\langle
f_{j}|$. Since $\Lambda_{\rm S}\otimes\Lambda_{\rm B}$ leaves
${\mathcal{H}}_{1}$ invariant, we have $(\Lambda_{\rm S}\otimes\Lambda_{\rm
B})\pi=\pi(\Lambda_{\rm S}\otimes\Lambda_{\rm B})\pi$. Now

$(\Lambda_{\rm S}\otimes\Lambda_{\rm B})\pi=\Lambda_{\rm
S}p_{11}\otimes\Lambda_{\rm B}q_{11}+\Lambda_{\rm S}p_{22}\otimes\Lambda_{\rm
B}q_{22}$ (75)

and, with $[\Lambda_{\rm S}]_{ij}=\langle e_{i},\Lambda_{\rm S}e_{j}\rangle$
and similarly for $\Lambda_{\rm B}$,

$\displaystyle\pi(\Lambda_{\rm S}\otimes\Lambda_{\rm B})\pi$ $\displaystyle=$
$\displaystyle[\Lambda_{\rm S}]_{11}[\Lambda_{\rm B}]_{11}\ p_{11}\otimes
q_{11}+[\Lambda_{\rm S}]_{12}[\Lambda_{\rm B}]_{12}\ p_{12}\otimes q_{12}$
(76) $\displaystyle=$ $\displaystyle[\Lambda_{\rm S}]_{21}[\Lambda_{\rm
B}]_{21}\ p_{21}\otimes q_{21}+[\Lambda_{\rm S}]_{22}[\Lambda_{\rm B}]_{22}\
p_{22}\otimes q_{22}.$

Taking the partial trace over ${\rm B}$ in (75) and (76) and equating the two
results gives

$[\Lambda_{\rm B}]_{11}\Lambda_{\rm S}p_{11}+[\Lambda_{\rm
B}]_{22}\Lambda_{\rm S}p_{22}=[\Lambda_{\rm S}]_{11}[\Lambda_{\rm B}]_{11}\
p_{11}+[\Lambda_{\rm S}]_{22}[\Lambda_{\rm B}]_{22}\ p_{22}.$ (77)

Since $[\Lambda_{\rm B}]_{11},[\Lambda_{\rm B}]_{22}>0$ we get from (77) that
$\Lambda_{\rm S}p_{jj}=[\Lambda_{\rm S}]_{jj}p_{jj}$ for $j=1,2$, so
$\Lambda_{\rm S}|e_{j}\rangle=[\Lambda_{\rm S}]_{jj}|e_{j}\rangle$. Hence the
restriction of $\Lambda_{\rm S}$ to ${\rm
span}\\{|e_{1}\rangle,|e_{2}\rangle\\}$ is diagonal, $\Lambda_{\rm
S}=[\Lambda_{1}]_{11}p_{11}+[\Lambda_{\rm S}]_{22}p_{22}+\Lambda^{\perp}_{\rm
S}$, where $\Lambda_{\rm S}^{\perp}$ is the block acting on the orthogonal
complement of that span. By taking the partial trace over ${\rm S}$ in (75)
and (76) and proceeding analogously, we see that $\Lambda_{\rm
B}=[\Lambda_{\rm B}]_{11}q_{11}+[\Lambda_{\rm B}]_{22}q_{22}+\Lambda_{\rm
B}^{\perp}$. It follows that $(\Lambda_{\rm S}\otimes\Lambda_{\rm
B})|v_{jj}\rangle=[\Lambda_{\rm S}]_{jj}[\Lambda_{\rm
B}]_{jj}\,|v_{jj}\rangle$ for $j=1,2$, so

$\eta=[\Lambda_{\rm S}]_{11}[\Lambda_{\rm B}]_{11}\,|v_{11}\rangle\langle
v_{11}|+[\Lambda_{\rm S}]_{22}[\Lambda_{\rm B}]_{22}\,|v_{22}\rangle\langle
v_{22}|$

is diagonal.

###  A.2 Conditions for $\bar{\rho}_{H}(t)=\bar{\rho}_{h_{W}}(t)$ and
$\mathcal{E}(\bar{\rho}_{H}(t))=\mathcal{E}(\bar{\rho}_{h_{W}}(t))$

In this section we assume the metric $\eta$ is of the form (30) with
$x_{1}=x_{2}=x>0$. Recall the formulas (41) and (51) for the quasi-Hermitian
and Hermitian density matrices.

First we ask when the two reduced density matrices coincide. We show that the
following statements are equivalent:

  * 1.

$\bar{\rho}_{H}(t)=\bar{\rho}_{h_{W}}(t)$ for all $t$ in an open interval
$I\subset{\mathbb{R}}$ and all $A,B\in\mathbb{C}$;

  * 2.

$\bar{\rho}_{H}(t)=\bar{\rho}_{h_{W}}(t)$ for all $t\in{\mathbb{R}}$ and all
$A,B\in\mathbb{C}$;

  * 3.

There are two real phases $\Phi_{1}$, $\Phi_{2}$, such that

$W=\begin{pmatrix}e^{i\Phi_{1}}&0\\\ 0&e^{i\Phi_{1}}\end{pmatrix}.$

1\. $\Rightarrow$ 3\. : Assume that $\bar{\rho}_{H}(t)=\bar{\rho}_{h_{W}}(t)$
for $t\in I$. Then $p(t)=q(t)$ for $t\in I$, where these quantities are given
in (42) and (55), respectively. Their equality is equivalent with

$\xi_{1}+\xi_{2}\cos(2\omega t)+\xi_{3}\sin(2\omega t)=0\mbox{\quad for all
$t\in I$},$ (78)

with $\xi_{1}=2x{\rm Re}(AB^{*})\,{\rm Re}(cd^{*})$,
$\xi_{2}=-|c|^{2}(1-2\alpha)-2x{\rm Im}(AB^{*}){\rm Im}(cd^{*})$ and
$\xi_{3}=-(1-2\alpha){\rm Im}(cd^{*})+2x(1-|c|^{2}){\rm Im}(AB^{*})$. As the
constant function and the sine and cosine are three independent functions, and
since (78) holds for all $t$ in an interval, we conclude that
$\xi_{1}=\xi_{2}=\xi_{3}=0$. Since $\xi_{1}=0$ for all $A,B$ and $x>0$, we
have ${\rm Re}(cd^{*})=0$, that is, ${\rm Im}(cd^{*})=cd^{*}$. The
coefficients $A,B$ and $\alpha$ are related by (33) and (43), resulting in
$A=\sqrt{\alpha}\sqrt{\frac{a_{1}}{xa_{2}}}e^{if_{1}}$ and
$B=\sqrt{1-\alpha}\sqrt{\frac{a_{2}}{xa_{1}}}e^{if_{2}}$, where
$f_{1},f_{2}\in\mathbb{R}$ are phases. This gives ${\rm
Im}(AB^{*})=\sqrt{\alpha(1-\alpha)}\frac{1}{x}\,{\rm Im}\,e^{i(f_{1}-f_{2})}$.
Then $\xi_{2}=0$ for all $A,B$ implies that

$|c|^{2}(1-2\alpha)=-2\sqrt{\alpha(1-\alpha)}cd^{*}\,{\rm
Im}\,e^{i(f_{1}-f_{2})}\quad\mbox{for all $f_{2},f_{2}\in\mathbb{R}$,
$\alpha\in[0,1]$.}$

This forces $|c|^{2}(1-2\alpha)=0=\alpha(1-\alpha)cd^{*}$ for all
$\alpha\in[0,1]$. Hence $c=0$. Then due to (49), $|d|=1$ and $bd^{*}=0$, so
$b=0$, and statement 3. holds.

3\. $\Rightarrow$ 2\. : Suppose $c=0$. Then $|d|=1$ and from (52) we have
$q(t)=x\frac{a_{1}}{a_{2}}|B(t)|^{2}$, which equals the population of $|0_{\rm
S}\rangle$ in $\bar{\rho}_{H}(t)$, see (39). Therefore 2. holds.

2\. $\Rightarrow$ 1\. : Obvious.

This completes the proof of the equivalence of the three statements 1.-3.

Next we ask when the entropies of the two density matrices coincide. We show
that the following statements 4.-6. are equivalent:

  * 4.

$\mathcal{E}(\bar{\rho}_{H}(t))=\mathcal{E}(\bar{\rho}_{h_{W}}(t))$ for all
$t$ in an open interval $I\subset{\mathbb{R}}$ and all $A,B\in\mathbb{C}$;

  * 5.

$\mathcal{E}(\bar{\rho}_{H}(t))=\mathcal{E}(\bar{\rho}_{h_{W}}(t))$ for all
$t\in{\mathbb{R}}$ and all $A,B\in\mathbb{C}$;

  * 6.

There are two real phases $\Phi_{1}$, $\Phi_{2}$ such that $W$ is of either of
the two forms

$W=\begin{pmatrix}e^{i\Phi_{1}}&0\\\ 0&e^{i\Phi_{1}}\end{pmatrix}\qquad\text{
or }\qquad W=\begin{pmatrix}0&e^{i\Phi_{1}}\\\ e^{i\Phi_{2}}&0\end{pmatrix}.$

4\. $\Rightarrow$ 6\. : Start by looking at the function
$\mathcal{E}(q)=-q\ln(q)-(1-q)\ln(1-q)$, for $q\in[0,1]$. It is clear from the
graph of $\mathcal{E}(q)$ (see the left panel of Figure 1) that
$\mathcal{E}(q)=\mathcal{E}(q^{\prime})$ exactly if either $q=q^{\prime}$ or
$q=1-q^{\prime}$. Consequently, if
${\mathcal{E}}(\bar{\rho}_{h_{W}}(t))={\mathcal{E}}(\bar{\rho}_{H}(t))$ for
all $t\in I$, then for each $t\in I$ individually, we have either $p(t)=q(t)$
or $p(t)=1-q(t)$. We now show that the same alternative must happen for all
$t\in I$.

Suppose first that $p(t_{0})\neq 1-q(t_{0})$ for some $t_{0}\in I$. Then by
the continuity of $p(t)$ and $q(t)$, we have $p(t)\neq 1-q(t)$ for all $t$ in
an open interval $I_{0}\subset I$ around $t_{0}$, so we must have $p(t)=q(t)$
for $t\in I_{0}$. But this means that
$\bar{\rho}_{H}(t)=\bar{\rho}_{h_{W}}(t)$ for all $t\in I_{0}$. Hence, as 2.
and 3. are equivalent, $W$ is of the diagonal form as given in point 3. above.
Similarly, if $p(t_{0})\neq q(t_{0})$ for some $t_{0}\in I$, we obtain
$p(t)=1-q(t)$ an an interval around $t_{0}$. Proceeding as in the proof of the
implication 1. $\Rightarrow$ 3\. above, this implies that $a=d=0$, so $W$ is
of the off-diagonal form given in statement 6. above.

6\. $\Rightarrow$ 5\. : If $W$ is of the diagonal form, then we already showed
that $p(t)=q(t)$ when we proved 3. $\Rightarrow$ 2\. In the same way, if $W$
is off-diagonal, then one sees that $p(t)=1-q(t)$. In either case,
$\mathcal{E}(\bar{\rho}_{H}(t))=\mathcal{E}(\bar{\rho}_{h_{W}}(t))$.

5\. $\Rightarrow$ 4\. : Obvious.

###  A.3 About the Dyson map

The idea of mapping a non-Hermitian Hamiltonian to a Hermitian one was
originally presented by Dyson in the context of the theory of magnetization
[8, 9]. Let ${\mathcal{H}}$ be a Hilbert space with inner product
$\langle\cdot|\cdot\rangle$ and let $H$ be an operator on ${\mathcal{H}}$ that
is not necessarily Hermitian with respect to $\langle\cdot|\cdot\rangle$.
Denote by $|\psi(t)\rangle=e^{-itH}|\psi(0)\rangle$ the solution of the
evolution equation

$i\partial_{t}|\psi(t)\rangle=H|\psi(t)\rangle.$ (79)

Next, let $S(t)$ be a differentiable family of operators on ${\mathcal{H}}$
such that $S(t)$ is invertible for each $t$, and set

$|\varphi(t)\rangle=S(t)|\psi(t)\rangle.$ (80)

Passing from $|\psi(t)\rangle$ to $|\varphi(t)\rangle$ represents a (possibly
time-dependent) change of variables. The evolution equation for
$|\varphi(t)\rangle$ is

$i\partial_{t}|\varphi(t)\rangle=h(t)|\varphi(t)\rangle,$ (81)

with

$h(t)=S(t)HS(t)^{-1}+i\dot{S}(t)S(t)^{-1},$ (82)

the dot being the time derivative. Conversely, if $|\varphi(t)\rangle$ solves
(81) then $|\psi(t)\rangle$ solves (79). Equation (82) is called the time-
dependent Dyson equation [13]. By means of $S(t)$, one may thus equivalently
solve (79) or (81). If $H$ is not Hermitian, one can look for $S(t)$ such that
the resulting $h(t)$ is Hermitian, hence trading a non-Hermitian problem with
constant Hamiltonian $H$ for a Hermitian problem with time-dependent
Hamiltonian $h(t)$. One readily sees that

$h(t)^{\dagger}=h(t)\quad\Longleftrightarrow\quad
i\partial_{t}\big{(}S^{\dagger}(t)S(t)\big{)}=H^{\dagger}\big{(}S^{\dagger}(t)S(t)\big{)}-\big{(}S^{\dagger}(t)S(t)\big{)}H.$
(83)

The operator $\eta(t)=S(t)^{\dagger}S(t)$ is automatically positive, so
$\eta(t)$ is a family of metrics. The equation for $\eta(t)$, according to
(83), is

$i\partial_{t}\eta(t)=H^{\dagger}\eta(t)-\eta(t)H.$ (84)

This is called the quasi-Hermiticity relation in [11]; note that it simplifies
to (1) if $\eta$ does not depend on time. It is clear that (84) has a unique
solution for any initial condition $\eta(0)$, namely

$\eta(t)=e^{-itH^{\dagger}}\eta(0)e^{itH},$ (85)

and that $\eta(t)$ is positive for all times if and only if it is positive at
some $t_{0}$.

A strategy to study the dynamics generated by a non-Hermitian $H$ is to find a
transformation $S(t)$ such that $h(t)$, as given by (82), is Hermitian, and
then analyze the dynamics of this Hermitian system using usual quantum
theoretical methods. Finding $S(t)$ for a specific Hamiltonian $H$ is not
easy, however. It often involves making a judicious ansatz containing
parameters that must solve rather complicated differential equations, which
are obtained by imposing the self-adjointness of $h(t)$. This can be done
explicitly for some models [4, 11, 13, 24, 27, 15, 16].

Given $H$, we seek _all_ possible $S(t)$, and the resulting Hermitian
Hamiltonians $h(t)$, with the sole requirement that
$\eta(t)=S(t)^{\dagger}S(t)$ is positive and satisfies the quasi-Hermiticity
relation (84). The solution $\eta(t)$ is uniquely determined by the initial
condition $\eta(0)$, which we may choose to be any positive, invertible
operator. The most general form of $S(t)$ is thus

$S(t)=W(t)\sqrt{\eta(0)}e^{itH},$ (86)

where $W(t)$ is any unitary family and $\sqrt{\eta(0)}$ denotes the unique
positive operator squaring to $\eta(0)$. The $h(t)$ associated to (86) by (82)
is

$h(t)=i\dot{W}(t)W(t)^{\dagger}.$ (87)

Note that we are entirely free to choose $W(t)$. For instance, given an
arbitrary $A=A^{\dagger}$, the choice $W(t)=e^{-itA}$ yields $h(t)=A$. This
means any time-independent Hermitian $h$ can be obtained from a suitable
choice of $S(t)$. A particularly simple choice is $S(t)=e^{itH}$, which
results from undoing the dynamics $e^{-itH}$ (going backwards in time) and has
$h=0$.

More generally, suppose $A(t)$ is a continuous family of operators, and let
$W(t)$ solve the differential equation

$i\dot{W}(t)=A(t)W(t).$ (88)

It is easily shown that if $A(t)=A(t)^{\dagger}$ for all $t$ and the initial
condition $W(0)$ is unitary, then the solution $W(t)$ is unitary for all $t$.
Choosing this $W(t)$, we find from (87) the Hermitian Hamiltonian
$h(t)=i\dot{W}(t)W(t)^{\dagger}=i\dot{W}(t)W(t)^{-1}=A(t)$. This means any
time-_dependent_ Hermitian $h(t)$ can also be obtained from a suitable choice
of $S(t)$.

Acknowledgements.  The authors are grateful to Andreas Fring for graciously
explaining the ideas and results of [13]. All authors are supported by
Discovery Grants from the Natural Sciences and Engineering Research Council of
Canada (NSERC).

## References

  * [1] C.M. Bender, S. Boettcher: Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry Phys. Rev. Lett. 80, 5243–5246 (1998)
  * [2] C.M. Bender, S. Boettcher, S. Meisinger: Real $PT$-Symmetric Quantum Mechanics, J. of Mat. Phys. 40, 2201-2229 (1999)
  * [3] C. M. Bender, D. C. Brody, H. F. Jones: Must a Hamiltonian be Hermitian?, Am. J. Phys. 71, 1095 (2003)
  * [4] O. A. Castro-Alvaredo, A. Fring: A spin chain model with non-Hermitian interaction: the Ising quantum spin chain in an imaginary field, J. Phys. A: Math. Theor. 42 (2009)
  * [5] J. Cen, A. Saxena: Anti-$PT$-symmetric qubit: Decoherence and entanglement entropy, Phys. Rev. A 105, 022404 (2022)
  * [6] J. Dieudonneé: Quasi-Hermitian operators Proceedings Of The International Symposium on Linear Spaces, 115–123 (1961)
  * [7] M. P. Drazin, E. V. Haynsworth: Criteria for the reality of matrix eigenvalues, Mathematische Zeitschrift 78, 449-452 (1962)
  * [8] F. J. Dyson: General Theory of Spin-Wave Interactions, Phys. Rev. 102(5), 1217-2030 (1956)
  * [9] F. J. Dyson: Thermodynamic Behavior of an Ideal Ferromagnet, Phys. Rev. 102(5), 1230-1244 (1956)
  * [10] J. Feinberg, M. Znojil: Which metrics are consistent with a given pseudo-Hermitian matrix? J. Math. Phys. 63, 013505 (2022)
  * [11] A. Fring, T. Frith: Exact analytical solutions for time-dependent Hermitian Hamiltonian systems from static unobservable non-Hermitian Hamiltonians, Phys. Rev. A 95, 010102(R) (2017)
  * [12] A. Fring, T. Frith: Mending the broken $PT$-regime via an explicit time-dependent Dyson map, Phys. Lett. A 381 2318-2323 (2017)
  * [13] A. Fring, T. Frith: Eternal life of entropy in non-Hermitian quantum systems, Phys. Rev. A 100, 010102(R) (2019)
  * [14] A. Fring, A., M.H.Y. Moussa: Unitary quantum evolution for time-dependent quasi-Hermitian systems with nonobservable Hamiltonians, Phys. Rev. A 93(4), 042114 (2016)
  * [15] A. Fring, R. Tenney: Spectrally equivalent time-dependent double wells and unstable anharmonic oscillators, Phys. Lett. A 384, 126530 (2020)
  * [16] A. Fring, R. Tenney: Exactly solvable time-dependent non-Hermitian quantum systems from point transformations, Phys. Lett. A 410, 127548 (2021)
  * [17] T. Frith: Time-dependence in non-Hermitian quantum systems , Ph.D. thesis, arXiv quant-ph, arXiv:2002.01977
  * [18] T. Frith: Exotic entanglement for non-Hermitian Jaynes–Cummings Hamiltonians, J. Phys. A: Math. Theor. 53 485303 (14pp) (2020)
  * [19] R. Horodecki, P. Horodecki, M. Horodecki, K. Horodecki: Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009)
  * [20] R. Kretschmer, L. Szymanowski: Quasi-Hermiticity in infinite-dimensional Hilbert spaces, Physics Letters A 325, 112 (2004)
  * [21] A. Mostafazadeh: Pseudo-Hermiticity versus $PT$-symmetry I. The necessary condition for the reality of the spectrum of a non-Hermitian Hamiltonian, J. Math. Phys. 43, 205-214 (2002)
  * [22] A. Mostafazadeh: Pseudo-Hermiticity versus $PT$-symmetry II. A complete characterization of non-Hermitian Hamiltonians with a real spectrum, J. Math. Phys. 43, 2814-2816 (2002)
  * [23] A. Mostafazadeh: Pseudo-Hermiticity versus $PT$-symmetry III. Equivalence of pseudo-Hermiticity and the presence of antilinear symmetries, J. Math. Phys. 43, 3944-3951 (2002)
  * [24] A. Mostafazadeh: $PT$-symmetric cubic anharmonic oscilator as a physical model, J. Phys. A: Math. Gen. 38, 6557-6570 (2005)
  * [25] A. Mostafazadeh, Seher Özçelik: Explicit Realization of Pseudo-Hermitian and Quasi-Hermitian Quantum Mechanics for Two-Level Systems, Turk. J. Phys. 30, 437-443 (2006)
  * [26] A. Mostafazadeh: Pseudo-Hermitian representation of quantum mechanics, Int. J. Geom. Methods in Mod. Phys. 7(7), 1191-1306 (2010)
  * [27] D.P. Musumbu, H.B. Geyer, W.D. Heiss: Choice of a metric for the non-Hermitian oscillator, J. Phys. A: Math. Gen. 40 F75-80 (2007)
  * [28] G. Scolarici, L, Solombrino: Time evolution of non-Hermitian quantum systems and generalized master equations, Czechoslovak J. of Phys. 56, No. 9 (2006)
  * [29] F.G. Scholtz, H.B. Geyer, F.J.W. Hahne: Quasi-Hermitian Operators in Quantum Mechanics and the Variational Principle, Ann. Phys. 213, 74-101 (1992)
  * [30] W. K. Wooters: Entanglement of Formation of an Arbitrary State of Two Qubits, Phys. Rev. Lett. 80, 2245 (1998)
  * [31] R. Zhang, H. Qin, J. Xiao: PT-symmetry entails pseudo-Hermiticity regardless of diagonalizability, J. Math. Phys. 61, 012101 (2020)


11institutetext: Computer Science and Engineering Department, Sichuan
University Jinjiang College, 620860 Pengshan, China
greenhat1016@gmail.com,mr.l172586418<EMAIL_ADDRESS>
# Hybrid Parallel Bidirectional Sieve based on SMP Cluster

Gang Liao  Lian Luo  Lei Liu

###### Abstract

In this article, hybrid parallel bidirectional sieve is implemented by SMP
Cluster, the individual computational units joined together by the
communication network, are usually shared-memory systems with one or more
multicore processor. To high-efficiency optimization, we propose average
divide data into nodes, generating double-ended queues (deque) for sieve that
are able to exploit dual-cores simultaneously start sifting out primes from
the head and tail.And each node create a FIFO queue as dynamic data buffer to
ache temporary data from another nodes send to. The approach obtains huge
speedup and efficiency on SMP Cluster. ybrid parallel, HPC, SMP Cluster, sieve

###### Keywords:

h

##  1 Introduction

Research into questions involving primes continues today, partly driven by the
importance of primes in modern cryptography. As our computational power
increases, researcher often pays more attention to Data analysis, Climate
modeling, Protein folding, Drug discovery etc. We can also exploit multicores
to efficiency solve some problem in the field of number theory.

M.Aigner and G.M.Ziegler [1] presented six quite different proofs of the
infinitude of primes. Mills[2] has shown that there is a constant $\Theta$
such that the function $f(n)=[{\Theta^{3}}^{n}]$ generates only primes. The
sieve of Eratosthenes-Legendre [3] [4] is an ancient algorithm for finding all
prime numbers up to any given limit. In number theory, tests distinguishing
between primes and composite integers will be crucial. The most basic
primality test is trial division, which tells us that integer $n$ is prime if
and only if it is not divisible by any prime not exceeding $\sqrt{n}$.

The computational complexity of algorithms for determining whether an integer
is prime is measured in terms of the number of binary digits in the integer.
The algorithm using trial divisions to determine whether an integer $n$ is
prime is exponential in terms of the number of binary digits of $n$, or in
terms of $\log_{2}n$ ,because $\sqrt{n}={2}^{log_{2}{n/2}}$.

As n gets large, an algorithm with exponential complexity quickly becomes
impractical. Leonard Adleman, Carl Pomerance, and Robert Rumely [5] [6]
developed an algorithm that can prove an integer is prime using $(\log
n)^{clogloglogn}$ nit operations, where c is a constant. In 2002, M. Agrawal,
N. Kayal, and N. Saxena [7], announced that they had found an algorithm PRIMES
is in P that can produce a certificate of primality for an integer n using
$O((logn)^{12})$ bit operations.

Karl Friedrich Gauss conjectured that $\pi(x)$ increases at the same rate as
the functions $\frac{x}{logx}$ and $Li(x)=\int_{2}^{x}\frac{dt}{logt}$. And
the Prime Number Theorem that the ratio of $\pi(x)$ to $\frac{x}{logx}$
approaches 1 as $x$ grows without bound. One way [11] to evaluate $\pi(x)$
only $O({x}^{\frac{3}{5}+e})$ bit operations without finding all the primes
less than $x$ is to use a counting argument based on the sieve of
Eratosthenes.

In this paper, Hybrid parallel bidirectional sieve based on SMP Cluster is
proposed to improve efficient and speedup. The result is proved to be
effective by MPI and OpenMP [8] [9] [10]. With Hybrid parallel, it has far-
reaching significance in cryptography.

##  2 Communication and Optimization

ILP and TLP provide parallelism at a very low level, they are typically
controlled by the processor and the operating system, and isn’t directly
controlled by the programmer. Parallel hardware is often classified using
Flynn’s taxonomy, which distinguished between the number of instruction
streams and the number of data streams a system can handle. A von Neumann
system is classified as SISD. Vector processors and graphics processing units
(GPU) are often classified as SIMD. MIMD execute multiple independent
instruction streams, each of which can have its own data stream. Shared-memory
or distributed-memory is typically MIMD. And most of the lager MIMD systems
are hybrid systems (Fig.1) in which a number of relatively small share-memory
are connected by an interconnection network. In such systems, the individual
shared-memory systems are sometimes called nodes.

Figure 1: SMP Cluster Architecture.

###  2.1 Interconnection networks

Currently the two most widely used interconnects on shared-memory systems are
buses and crossbars [15]. The key characteristic of a bus is that the
communication wires are shared by the devices that are connected to it. Buses
have the virtue of low cost and flexibility. Crossbars (Fig.2) allow
simultaneous communication among different devices, so they are much faster
than uses. But the cost of the switches and links is relatively high.

Distributed-memory interconnects are often divided into two groups: direct
interconnects and indirect interconnects. One measure of ”number of
simultaneous communications” or ”connectivity” is bisection width. To
understand this measure, imagine that the parallel system is divided into two
halves, and each half contains half of the processors or nodes. An alternate
way of computing the bisection width is to remove the minimum number of links
needed to split the set of nodes into two equal halves.

Figure 2: Shared-memory system simultaneous memory access

The hypercube (Fig.3) is a highly connected direct interconnect that has been
used in actual system. A hypercube of dimension d has $p={2}^{d}$ nodes, and a
switch in a d-dimensional hypercube is directly connected to a processor and d
switches. The bisection width of a hypercube is $\frac{p}{2}$.The switches
support $1+d=1+\log_{2}p$ wires. The hypercube is more powerful and expensive
to construct.

Figure 3: (a) two-dimensional hypercube (b) three-dimensional hypercube (c)
four-dimensional hypercube

The crossbar and the omega network are relatively simple examples of indirect
networks. The omega network (Fig.4) is less expensive than crossbar. The omega
network uses $\frac{1}{2}plog_{2}(p)$ of the 2 x 2 crossbar switches, so it
uses a total of ${2}plog_{2}(p)$ switches, while the crossbar users $p^{2}$.

###  2.2 Hybrid Parallelism

We define the speedup of a parallel program to be
$S=\frac{T_{serial}}{T_{parallel}}$ . Then linear speedup has $S=Pcores$, this
value, $\frac{S}{P}$, is sometimes called the efficiency of the parallel
program as follows:

$E=\frac{S}{P}=\frac{\frac{T_{serial}}{T_{parallel}}}{P}$ (1)

Back in the 1960s, Gene Amdahl [13] that’s become as Amdahl’s Law:

$S_{overall}=\frac{1}{(1-f)+\frac{f}{s}}$ (2)

It means that unless virtually all of a serial program is parallelized, the
possible speedup is going to be very limited-regardless of the number of cores
available. A more mathematical version of this statement is known as
Gustafson’s Law [14].

Unfortunately, there are several mismatch problem between the (hybrid)
programming schemes and the hybrid hardware architecture. Often, one can see
in publications, that applications may or may not benefit from hybrid
programming depending on some application parameters, e.g., in [16][17][18]
[19].

Polf Rabenseifner analyses strategies to overcome typical drawbacks of this
easily usable programming scheme on systems with weaker inter-connects [20].
Best performance can be achieved with overlapping communication and
computation, but this scheme is lacking in ease of use. Often, hybrid MPI $+$
OpenMP programming denotes a programming style with OpenMP shared memory
parallelization inside the MPI processes (i.e., each MPI process itself has
several OpenMP threads) and communication with MPI between the MPI processes,
but only outside of parallel regions.

This hybrid programming scheme will be named materonly in the following
classification, which is based on the question, when and by which thread(s)
the messages are sent between the MPI processes:

  * .

Pure MPI

  * .

Hybrid MPI $+$ OpenMP

  * .

Overlapping communication and computation

  * .

Pure OpenMP

Overlapping of communication and computation is a chance for an optimal usage
of the application itself, in the OpenMP parallelization and in the load
balancing. It requires a coarse-grained and thread-rank-based OpenMP
parallelization, the separation of halo-based computation from the computation
that can be overlapped with communication, and the threads with different
tasks must be load balanced. Advantages of the overlapping scheme are:

  * .

the problem that one CPU may not achieve the inter-node bandwidth is no longer
relevant as long as there is enough computational work that can be overlapped
with the communication

  * .

the saturation problem is solved as long as not more CPUs communicate in
parallel than necessary to achieve the inter-node bandwidth

  * .

the sleeping threads problem is solved as long as all computation and
communication is load balanced among the threads.

Figure 4:  (a) Crossbar (b) omega network

##  3 Bidirectional Sieve Model

Foster’s methodology [12] provides an outline of steps include

  * .

Partitioning.

  * .

Communication.

  * .

Agglomeration or aggregation

  * .

Mapping for parallel programming

###  3.1 Algorithm Design

The sieve of Eratosthenes does so by iteratively marking as composite the
multiples of each prime, starting with the multiples of 2 [4]. We can exploit
and improve the sieve of Eratosthenes based on SMP Cluster (Fig. 5). Assume
that there are some disorder integers which the scale of $n$, and when each
node sieve the integers in the block that the scale of $k$, it could achieve
high-efficiency optimization. We conjectured that the SMP Cluster requires at
least N nodes.The formula as follows:

$N=\frac{n}{k}+(n\bmod k)\And{1}$ (3)

And each node generate one deque and do with dual-cores. One core is located
in the head of the deque. On the contrary, the other one is located in the
tail of the deque. It’s easy to deduction the formula about the amount of
cores($C_{cores}$) and deques($D_{deques}$):

$C_{cores}=D_{deques}={2}N$ (4)

There is another point that’s worth considering. In most cases, the scale of
node $N$ is not exactly equal $k$. We can deal with the state as follows
Alg.1:

Algorithm 1 the scale of node $N^{th}$

0: $K$ denote that the currency scale of node $N^{th}$

0: $k$ denote that the general scale of node

if ${0}\leq K\leq\frac{k}{2}$ then

Node N assign single core to right or left sieve

else

Node N assign dual-cores to simultaneous bidirectional sieve

end if

Figure 5: Construct Bidirectional Sieve

And its flow diagram is shown in Fig.6.

Figure 6: High-level flow diagram of hybrid parallel bidirectional Sieve

###  3.2 Primality Testing : Non-deterministic

Primality testing of a number is perhaps the most common problem concerning
number theory.The problem of detecting whether a given number is a prime
number has been studied extensively but nonetheless,it turns out that all the
deterministic algorithms for this problem are too slow to be used in real life
situations and the better ones amongst them are tedious to code.But,there are
some probabilistic methods which are very fast and very easy to
code.Moreover,the probability of getting a wrong result with these algorithms
is so slow that it can be neglected in normal situations.

All the algorithms which we are going to discuss will require you to
efficiently compute $(a^{b})\bmod c$ (where a,b,c are non-negative integers).
A straightforward algorithm to do the task can be to iteratively multiply the
result with $a$ and take the remainder with $c$ at each step,this algorithm
takes $O(b)$ time and is not very useful in practice. We can do it
$O(\log{b})$ by using what is called as exponentiation by squaring as follows:

$f(n)=$ $\begin{cases}(a^{2})^{\frac{b}{2}},&\mbox{if }b\mbox{ is even and b
$>$ 0}\\\ a(a^{2})^{\frac{b-1}{2}},&\mbox{if }b\mbox{ is odd}\\\ $1$,&\mbox{if
}b\mbox{ $=0$}\end{cases}\\\ $

Algorithm 2 modulo(a,b,c) : Exponentiating by squaring to $(a^{b})\bmod c$

0: $x=1,y=a$

0: $(a^{b})\bmod c$

while $b>0$ do

if $b\And 1$ then

$x=(x*y)\bmod c$

end if

$y=(y*y)\bmod c$

$b>>=1$

end while

return $x\bmod c$

Pierre de Fermat first stated the Fermat’s Little Theorem in a letter dated
October 18, 1640, to his friend and confidant Fr$\acute{e}$nicle de Bessy as
the following [7]:

$a^{p}=a\pmod{p}$ (5)

or alternatively:

$a^{p-1}=1\pmod{p}$ (6)

According to Fermat’s Little Theorem[7], if $p$ is a prime number and a is
positive integer less than $p$ ($a<p$),and then calculate $a^{p-1}\bmod p$. If
the result is not 1, then by Fermat’s Little Theorem p cannot be prime.The
more iterations we do, the higher is the probability that our result is
correct.

Algorithm 3 Fermat(p,iterations) : Fermat′s primality test

if $p=1$ then

return $false$

end if

for $i:=1$ to $iterations$ do

$a=rand()\bmod(p-1)+1$

if modulo(a,p-1,p)!=1 (Alg.2) then

return $false$

end if

end for

return $true$

Though Fermat is highly accurate in practice there are certain composite
numbers $p$ known as Carmichael numbers for which all values of $a<p$ for
which $gcd(a,p)=1$,$(a^{p-1})\bmod p=1$.And in that case,the Fermat’s test
will return wrong result with very high probability.Out of the Carmichael
numbers less than ${10}^{16}$,about $95\%$ of them are divisible by primes
$<1000$.However,there are other improved primality tests which don’t have this
flaw as Fermat’s(e.g.Rabin-Miller test[21][22],Solovay-Strassen test [23]).

##  4 Performance Analysis

Figure 7: statistics and analysis hybrid parallel bidirectional sieve with
general method

Different programming schemes on clusters of SMPs show different performance
benefits or penalties in this paper. Fig.7 summarizes the result of hybrid
parallel bidirectional sieve .It’s obvious that nodes communication would
waste most of time when data scale is tiny.Even its slower than general
method.However, if there are hyper-data scale,hybrid parallel show huge
efficiency and optimization.Indeed,sometimes the waste of communication could
be neglected.In that case,multicores parallelism is effective approach to
solve some problem in number theory.

To achieve an optimal usage of the hardware,one can also try to use the idling
CPU’s for other applications,especially low-priority single-threaded or multi-
threaded non-MPI application if the parallel high-priority hybrid application
does not use the total memory of the SMP nodes.

##  5 Conclusion

In this study we haven shown that hybrid parallel on SMP cluster is an
applicable method to implement bidirectional sieve . The analysis demonstrated
that even hybrid parallel bidirectional sieve is efficiency and optimization
solution.

As our computational power increases,Most HPC system are clusters of shared
memory nodes.Parallel programming must combine the distributed memory
parallelization on the node inter-connect with shared memory parallelization
inside of each node.And Each parallel programming schema on hybrid
architecture has one or more significant drawbacks(e.g. sleeping-thread and
saturation problem). However,Hybrid parallel also has far-reaching
significance in many fields(e.g.Cryptography,Data analysis, Climate modeling,
Protein folding, Drug discovery).

We believe that hybrid parallel bidirectional sieve can be properly modeled
using techniques form number theory and this article is just an early trial of
using hybrid parallelism to improve speedup and efficiency.

## References

  * [1] M. Aigner , G. M. Ziegler: Proofs from THE BOOK,3rd ed.,Springer-Verlag,Berlin,(2003)
  * [2] W.H. Mills: A prime-representing function, Bulletin of the American Mathematical Society, Volume 53,604,(1947)
  * [3] H.Halberstam , H.-E.Richert: Sieve Methods, Academic Press, London,(1974)
  * [4] Horsley, Rev. Samuel, F. R. S.: The Sieve of Eratosthenes. Being an Account of His Method of Finding All the Prime Numbers, Philosophical Transactions (1683 C1775), Vol. 62., pp. 327 C347,(1772)
  * [5] L.M. Adleman, C. Pomerance, R.S. Rumly, On distinguishing prime numbers from composite numbers, Annals of Mathematics, Volume 117 (1983).
  * [6] R. Rumely,:Recent advances in primality testing, Notices of American Mathmatical Society, Volume 30 , 475-477,(1983)
  * [7] M.A Agrawal, N. Kayal, N. Saxena : PRIMES is in P, Department of Computer Science & Engineering, Indian Institute of Technology, Kanpur, India,(2002)
  * [8] R. Chandra, et al.: Parallel Programming in OpenMP, Morgan Kaufmann, San Francisco,(2001)
  * [9] P. Pacheco: Parallel Programming with MPI, Morgan Kaufmann, San Francisco,(1997)
  * [10] M.Quinn: Parallel Programming in C with MPI and OpenMP, McGraw-Hill Higher Education, Boston,(2004).
  * [11] J.C. Lagarias , A.M. Odlyzko: New algorithm for computing PI(x), Bell Laboratories Technical Memorandum TM-82-11218-57.
  * [12] I. Foster: Designing and Building Parallel Programs, Addison-Wesley, Reading, MA, 1995. Also available from http://www.mcs.anl.gov/ itf/dbpp/ (accessed 21.09.10)
  * [13] G.M. Amdahl: Validity of the single processor approach to achieving large scale computing capabilities, in: Proceedings of the American Federation of Information Processing Societies Conference, vol. 30, issue 2, Atlantic City, NJ, 1967, pp. 483-485
  * [14] J.L. Gustafson: Reevaluating Amdahl’s law, Commun. ACM 31 (5) (1988) 532-533
  * [15] Peter S. Pacheco: An introduction to PARALLEL PROGRAMMING, Elsevier (Singapore) Pte Ltd,(2011)
  * [16] Georg Hager, Frank Deserno, Gerhaed Wellein: Pseudo-Vectorization and RISC Optimization Techniques for the Hitachi SR8000 Architecture, in High Performance Computing in Science and Engineering in Munich ’02, Springer-Verlag Berlin Heidelberg, (2003)
  * [17] D. S. Henty: Performance of hybrid message-passingand shared-memory parallelism for discrete element modeling, in Proc. Supercomputing’00, Dallas, TX, (2000).
  * [18] Richard D. Loft, Stephen J. Thomas, John M. Dennis: Terascale spectral element dynamical core for atmospheric general circulation models, in proceedings, SC 2001, NOW. 2001, Nov. 2001, Denver, USA. www.sc2001.org/papers/pap.pap189.pdf
  * [19] Gerhard Wellein, GeorgHager, Achim Basermann, Holger Fehske: Fast sparse matrix-vector multiplication for TeraFlop/s computers, in proceedings of VECPAR’2002, 5th Int’l Conference on High Performance Computing and Computational Science, Porto, Portugal, June 26-28, 2002, part I, pp 57-70. http://vecpar.fe.up.pt/
  * [20] Rolf Rabenseifner: Hybrid Parallel Programming: Performance Problems and Chances, in proceeding of the 45th CUG Conference 2003, Columbus, Ohio, USA, May 12-16,2003, www.cug.org
  * [21] Miller, Gary L.:Riemann’s Hypothesis and Tests for Primality, Journal of Computer and System Sciences 13 (3): 300 C317, doi:10.1145/800116.803773,(1976)
  * [22] Rabin, Michael O.:Probabilistic algorithm for testing primality, Journal of Number Theory 12 (1): 128 C138, doi:10.1016/0022-314X(80)90084-0,(1980)
  * [23] Solovay, Robert M.Strassen, Volker.:A fast Monte-Carlo test for primality”. SIAM Journal on Computing 6 (1): 84 C85. doi:10.1137/0206006,(1977)