# High-order series expansion of non-Hermitian quantum spin models

Lea Lenke  Lehrstuhl für Theoretische Physik I, Staudtstraße 7, Universität
Erlangen-Nürnberg, D-91058 Erlangen, Germany  Matthias Mühlhauser  Lehrstuhl
für Theoretische Physik I, Staudtstraße 7, Universität Erlangen-Nürnberg,
D-91058 Erlangen, Germany  Kai Phillip Schmidt  Lehrstuhl für Theoretische
Physik I, Staudtstraße 7, Universität Erlangen-Nürnberg, D-91058 Erlangen,
Germany

###### Abstract

We investigate the low-energy physics of non-Hermitian quantum spin models
with $PT$-symmetry. To this end we consider the one-dimensional Ising chain
and the two-dimensional toric code in a non-Hermitian staggered field. For
both systems dual descriptions in terms of non-Hermitian staggered Ising
interactions in a conventional transverse field exist. We perform high-order
series expansions about the high- and low-field limit for both systems to
determine the ground-state energy per site and the one-particle gap. The one-
dimensional non-Hermitian Ising chain is known to be exactly solvable. Its
ground-state phase diagram consists of second-order quantum phase transitions,
which can be characterized by logarithmic singularities of the second
derivative of the ground-state energy and, in the symmetry-broken phase, the
gap closing of the low-field gap. In contrast, the gap closing from the high-
field phase is not accessible perturbatively due to the complex energy and the
occurrence of exceptional lines in the high-field gap expression. For the two-
dimensional toric code in a non-Hermitian staggered field we study the quantum
robustness of the topologically ordered phase by the gap closing of the low-
field gap. We find that the well-known second-order quantum phase transition
of the toric code in a uniform field extends into a large portion of the non-
Hermitian parameter space. However, the series expansions become unreliable
for a dominant anti-Hermitian field. Interestingly, the analysis of the high-
field gap reveals the potential presence of an intermediate region.

##  I Introduction

Non-Hermitian Hamiltonians with $PT$-symmetry represent an interesting
extension of conventional quantum mechanics [1, 2], since such systems can
display spontaneous $PT$-symmetry breaking from a purely real to a complex
energy spectrum. Non-Hermitian operators typically arise when gain and loss
terms are appropriately adjusted in open physical systems [3]. As a
consequence, many experimental platforms are explored for non-Hermitian
quantum dynamics like optical waveguides, electronics, microwaves, acoustics
and single-spin systems [4]. On the theoretical side, many investigations have
focused either on single-particle or non-interacting non-Hermitian quantum
systems, e.g., fermionic or bosonic band insulators with topologically non-
trivial band structures [5] while non-Hermitian interacting quantum many-body
Hamiltonians are far less explored and understood.

One class of non-Hermitian interacting quantum many-body systems are solvable
one-dimensional quantum spin models [6, 7, 8, 9, 10, 11, 12] like the Ising
chain in a non-Hermitian staggered field [13], which can be mapped in many
cases to free fermions. However, higher-dimensional generalizations have not
been treated to the best of our knowledge. Another class of quantum many-body
systems are systems displaying intrinsic topological order which have a long-
range entangled ground state and exotic anyonic excitations with non-trivial
fractional statistics [14, 15]. One important question is the robustness of
intrinsic topological order in non-Hermitian quantum systems. Here an
attractive starting point are exactly solvable stabilizer codes like the toric
code [16], which has been extended to non-Hermitian operators with
$PT$-symmetry recently [17, 18, 19]. It has been found that the correspondence
between bulk quasi-particles and topologically protected degenerate ground
states breaks down. Further, a continuous quantum phase transition without gap
closing was explored that occurs in non-Hermitian topological orders [20].

Technically, one important tool to tackle quantum properties of interacting
quantum spin systems are high-order series expansions [21], which have been
used, among many other systems, to extract quantum critical properties in
transverse-field Ising models [22, 23, 24, 25, 26] as well as for topological
phase transitions in the toric code in the presence of external fields [27,
28, 29, 30]. However, these techniques have not been applied to non-Hermitian
quantum spin models which is exactly the main purpose of this article. We
extend high-order series expansions to investigate quantum critical properties
of non-Hermitian quantum spin models. In particular, we calculate high-order
series of energy gaps about different limits using the method of perturbative
continuous unitary transformations [31, 32, 25] or Takahashi’s perturbation
theory [33] and apply extrapolations to extract critical points and associated
critical exponents. In one dimension, we study the exactly solvable Ising
chain in a non-Hermitian staggered field [13] which allows a direct comparison
between the analytic solution and the findings from high-order series
expansions. In two dimensions, we replace the Ising interaction by the toric
code [16] in order to understand the breakdown of the topological phase by the
non-Hermitian staggered field.

The article is structured as follows. In Sec. II we introduce the one-
dimensional Ising chain and the two-dimensional toric code in a $PT$-symmetric
non-Hermitian staggered field and we describe duality transformations for both
systems to non-Hermitian staggered Ising interactions in a conventional
transverse field. Technical aspects of the high-order low- and high-field
series expansions and their extrapolations are given in Sec. III. All results
for the one-dimensional Ising chain in a non-Hermitian staggered field are
contained in Sec. IV. This includes the analytic solution as well as the
comparison to series expansion results. In Sec. V we discuss the breakdown of
the topological phase due to a non-Hermitian staggered Ising interaction. The
article is concluded in Sec. VI.

##  II Non-Hermitian models

We consider models with $N$ spin-$1/2$ particles on lattices that are coupled
by some interaction Hamiltonian $\mathcal{H}_{\textrm{I}}$ and subject to a
non-Hermitian staggered magnetic field. The prefactor of the magnetic field is
either given by $g\equiv\eta+\mathrm{i}\xi$ for spins on sublattice $\circ$ or
by $g^{*}\equiv\eta-\mathrm{i}\xi$ for spins on sublattice $\bullet$. The full
Hamiltonian reads

$\displaystyle\mathcal{H}=\mathcal{H}_{\textrm{I}}-g\sum_{j\in\circ}\sigma^{z}_{j}-g^{*}\sum_{j\in\bullet}\sigma^{z}_{j},$
(1)

where $\sigma^{\alpha}_{j}$ represents the Pauli matrix with flavor $\alpha$
on site $j$. The imaginary parts $\pm\xi$ of the magnetic field prefactors
correspond to the anti-Hermitian part of the Hamiltonian. The assignment of
spins to $\circ$ or $\bullet$ is done such that the full Hamiltonian is
$\mathcal{PT}$-symmetric.

For $\mathcal{H}_{\textrm{I}}=0$, the so-called high-field limit, the model
(1) is exactly solvable. The state with the lowest real eigenvalue $-|\eta|N$
is $\ket{\Uparrow}$. We will refer to this state as the ground state in the
high-field limit. If additionally $\eta=0$, the ground-state energy is
infinitely degenerate in the thermodynamic limit; we will not treat this case
with series expansions about the high-field limit. There are two types of
elementary spin-flip excitations, depending on whether the flipped spin is on
sublattice $\circ$ or $\bullet$. Their eigenvalues are given by
$-|\eta|(N-2)\pm 2\mathrm{i}\xi$ and are non-real for $\xi\neq 0$.

$g$
---
$g^{*}$
---
$J$
---
(a)
---
(b)
---
$A_{s}$
---
$B_{p}$
---
Figure 1: Ising chain (a) and toric code (b) in a staggered magnetic field.
The sites represent the spin-1/2 particles that are coupled by an Ising
interaction (a) or by the toric code operators (b). The colors of the sites
encode the imaginary part of the magnetic prefactor that is given by
$g=\eta+\mathrm{i}\xi$ for empty red sites, belonging to $\circ$, and
$g^{*}=\eta-\mathrm{i}\xi$ for filled green sites, belonging to $\bullet$. The
toric code star operators $A_{s}$ are depicted in yellow, the plaquette
operators $B_{p}$ in blue.

We study two interaction Hamiltonians $\mathcal{H}_{\textrm{I}}$: the one-
dimensional Ising chain and the 2D toric code [16] on a square lattice. Both
models are illustrated in Fig. 1. The full Hamiltonian of the Ising chain in a
staggered transverse magnetic field reads

$\displaystyle\mathcal{H}_{\textrm{TFIM}}=-J\sum_{j=1}^{N}\sigma^{x}_{j}\sigma^{x}_{j+1}-g\sum_{j\in\circ}\sigma^{z}_{j}-g^{*}\sum_{j\in\bullet}\sigma^{z}_{j},$
(2)

where we set $\sigma^{x}_{N+1}=\sigma^{x}_{1}$ and thus couple the chain
periodically. The Ising terms $\sigma^{x}_{j}\sigma^{x}_{j+1}$ have
eigenvalues $\pm 1$ depending on whether the respective spins are aligned or
anti-aligned. These eigenvalues are not conserved as they do not commute with
the Hamiltonian (2). The Hermitian version of this model is exactly solvable
[34, 35] using a Jordan-Wigner transformation [36] and a Bogoliubov
transformation in Fourier space. Li et al. [13] found a solution of the non-
Hermitian generalization that they obtain via a Jordan-Wigner transformation
and composite operators. In Sec. IV.1 we apply an alternative approach using a
generalized Bogoliubov transformation in momentum space to derive the exact
solution of the Ising chain in a non-Hermitian staggered magnetic field
confirming their findings.

The Hamiltonian of the toric code in a staggered parallel magnetic field reads

$\displaystyle\mathcal{H}_{\textrm{TCF}}=-J\sum_{s}A_{s}-J\sum_{p}B_{p}-g\sum_{j\in\circ}\sigma^{z}_{j}-g^{*}\sum_{j\in\bullet}\sigma^{z}_{j}.$
(3)

Every star operator $A_{s}\equiv\prod_{i\in s}\sigma^{x}_{i}$ commutes with
every plaquette operator $B_{p}\equiv\prod_{i\in p}\sigma^{z}_{i}$ because
they always share an even number of spins. Thus all of them commute with the
interaction Hamiltonian. The plaquette operators commute with the magnetic
field terms as well and therefore also with the full Hamiltonian (3). The star
and plaquette operators have eigenvalues $\pm 1$ whereof the eigenvalues of
the plaquette operators are conserved due to the afore mentioned commutation
relations. That is why the Hamiltonian decomposes into blocks for the
different combinations of plaquette eigenvalues.

The relevant low-energy physics takes place in the sector without plaquette
excitations, i.e., the block where all plaquette operators have eigenvalue
$\mathrm{sgn}(J)$. In this sector, the Hamiltonian simplifies to

$\displaystyle\mathcal{\widetilde{H}}_{\textrm{TCF}}=-|J|\frac{N}{2}-J\sum_{s}A_{s}-g\sum_{j\in\circ}\sigma^{z}_{j}-g^{*}\sum_{j\in\bullet}\sigma^{z}_{j}.$
(4)

For both systems (2) and (4) it suffices to study a restricted parameter
space. We may without loss of generality choose $\xi\geq 0$, because any sign
in front of $\xi$ can be absorbed into a reassignment of spins to $\circ$ and
$\bullet$. We may also choose $J,\eta\geq 0$ because any sign in front of $J$
or in front of both, $\eta$ and $\xi$, can be absorbed into a rotation in spin
space of some spins around the $z$-axis or all spins around the $x$-axis,
respectively.

Interestingly, both models are dual to the transverse field Ising model with
non-Hermitian staggered Ising interaction. The dual model of the Ising chain
(2) lives on a chain with effective sites centered on the links of the
original model. Similarly, for the toric code (4) in the subspace with all
plaquette operators having eigenvalues $+1$, the dual lattice is a square
lattice build by the centers of star operators. The effective sites are
labeled by $\mu$ and $\nu$. Because the interaction Hamiltonian contains in
both cases only operators with eigenvalues $\pm 1$, this mapping uses Pauli
matrices $\tau^{z}$ for pseudo-spins 1/2 that are defined to have the same
eigenvalues so that $\mathcal{H}_{\textrm{I}}$ is mapped to an effective
magnetic field in $z$-direction. Further, each local operator of the magnetic
field term changes the eigenvalues of two adjacent pseudo-spins on sites $\mu$
and $\nu$ which therefore gives an effective Ising interaction
$\tau^{x}_{\mu}\tau^{x}_{\nu}$. The resulting dual Hamiltonian reads then in
both cases

$\displaystyle\mathcal{H}=\bar{E}_{0}-J\sum_{\mu}\tau^{z}_{\mu}-g\sum_{\langle\mu,\nu\rangle\in\circ}\tau^{x}_{\mu}\tau^{x}_{\nu}-g^{*}\sum_{\langle\mu,\nu\rangle\in\bullet}\tau^{x}_{\mu}\tau^{x}_{\nu},$
(5)

where $\langle\mu,\nu\rangle$ represents bonds of the dual lattice. The
assignment of bonds to the former sublattices $\circ$ and $\bullet$ represent
the colors of the bonds. Note that the explicit value of the constant
$\bar{E}_{0}$ differs for the two cases, but plays no role physically.

##  III Series expansions

In order to explore the ground-state phase diagrams, we perform high-order
series expansions about the low-field limit, where $\eta=\xi=0$, and about the
high-field limit, where $J=0$. We will test the series expansion by comparing
to the exactly solvable Ising chain in a non-Hermitian staggered transverse
magnetic field given in Sec. IV.1. The goal of these methods is to transform
the Hamiltonians into effective quasi-particle (qp) conserving ones. Being
quasi-particle conserving, the effective Hamiltonians decompose into decoupled
blocks for different quasi-particle numbers, which enables us to obtain access
to the relevant low-energy physics, i.e., the ground-state energy via the zero
qp block and the one-particle gap via the one qp block.

###  III.1 Low-field expansions

In the following we study the low-field (lf) limit $g\ll J$ perturbatively. We
introduce hardcore-bosonic operators

$\displaystyle b^{\phantom{\dagger}}_{\mu}$
$\displaystyle=(\tau^{x}_{\mu}+\mathrm{i}\tau^{y}_{\mu})/2,$ $\displaystyle
b^{\dagger}_{\mu}$ $\displaystyle=(\tau^{x}_{\mu}-\mathrm{i}\tau^{y}_{\mu})/2$
(6)

to express Eq. (5) exactly as follows

$\displaystyle\mathcal{H}$
$\displaystyle=-JN+2J\sum_{\mu}b^{\dagger}_{\mu}b^{\phantom{\dagger}}_{\mu}-g\sum_{\langle\mu,\nu\rangle\in\circ}\left(b^{\dagger}_{\mu}+b^{\phantom{\dagger}}_{\mu}\right)\left(b^{\dagger}_{\nu}+b^{\phantom{\dagger}}_{\nu}\right)$
$\displaystyle\hphantom{{}=}-g^{*}\sum_{\langle\mu,\nu\rangle\in\bullet}\left(b^{\dagger}_{\mu}+b^{\phantom{\dagger}}_{\mu}\right)\left(b^{\dagger}_{\nu}+b^{\phantom{\dagger}}_{\nu}\right),$
(7)

where
$\mathcal{Q}\equiv\sum_{\mu}b^{\dagger}_{\mu}b^{\phantom{\dagger}}_{\mu}$ is a
counting operator reflecting the equidistant spectrum of the unperturbed
Hamiltonian $g=0$. The unperturbed ground state $\ket{0}_{\textrm{lf}}$ is
then defined by $\mathcal{Q}\ket{0}_{\textrm{lf}}=0$. Rescaling the
Hamiltonian by $2J$ yields

$\displaystyle\frac{\mathcal{H}}{2J}=-\frac{N}{2}+\mathcal{Q}-\frac{g}{2J}\sum_{n\in\\{-2,0,2\\}}T_{n}^{\circ}-\frac{g^{*}}{2J}\sum_{n\in\\{-2,0,2\\}}T_{n}^{\bullet},$
(8)

where $[\mathcal{Q},T_{n}^{\circ/\bullet}]=nT_{n}^{\circ/\bullet}$, i.e.
$T_{n}^{\circ/\bullet}$ creates $n$ elementary excitations. Every operator
used here is Hermitian, the anti-Hermitian parts arising only from the
imaginary part $\pm\xi$ of the magnetic field prefactors. This aspect is
essential for the straightforward generalization of the method of perturbative
continuous unitary transformations (pCUTs) [31, 32] to non-Hermitian
Hamiltonians (1).

For Hermitian Hamiltonians, the method of continuous unitary transformations
[37] works by introducing a unitary transformation $U(\ell)$ that depends on a
parameter $\ell$ flowing continuously from $0$ to $\infty$. The so called flow
Hamiltonian $\mathcal{H}(\ell)=U^{-1}(\ell)\mathcal{H}U(\ell)$ is also
$\ell$-dependent and becomes the desired effective Hamiltonian as $\ell$
approaches $\infty$. The unitary transformation is not explicitly needed. It
suffices to introduce its anti-Hermitian infinitesimal generator $\eta(\ell)$.
The effective Hamiltonian can be calculated using the flow equation

$\displaystyle\frac{\textrm{d}\mathcal{H}(\ell)}{\textrm{d}\ell}=[\eta(\ell),\mathcal{H}(\ell)]$
(9)

and taking the limit $\ell$ to $\infty$. If the Hermitian Hamiltonian consists
of a counting operator $\mathcal{Q}$ representing the unperturbed Hamiltonian
and a finite number of $T_{n}$ operators with $[\mathcal{Q},T_{n}]=nT_{n}$
connected to one or more perturbation parameters, the pCUT method is
applicable. For models with one perturbation parameter $\lambda$, the ansätze
for the infinitesimal generator and the flow Hamiltonian in the pCUT method
are

$\displaystyle\eta(\ell)$
$\displaystyle=\sum_{k=1}^{\infty}\lambda^{k}\sum_{(m_{1},\ldots,m_{k})}\mathrm{sgn}(m_{1}+\cdots+m_{k})F(\ell;(m_{1},\ldots,m_{k}))T_{m_{1}}\cdots
T_{m_{k}},$ (10) $\displaystyle\mathcal{H}(\ell)$
$\displaystyle=\mathcal{Q}+\sum_{k=1}^{\infty}\lambda^{k}\sum_{(m_{1},\ldots,m_{k})}F(\ell;(m_{1},\ldots,m_{k}))T_{m_{1}}\cdots
T_{m_{k}},$ (11)

where $F$ are real-valued coefficient functions. The generalization to
multiple perturbation parameters is straightforward. Inserting Eqs. (10) and
(11) into the flow equation (9), the functions $F$ can be calculated
recursively order by order. This can be done computer-aided up to high orders
resulting in an effective Hamiltonian $\mathcal{H}_{\textrm{eff}}$ for
$\ell\rightarrow\infty$ with rational coefficients which is blockdiagonal in
the qp number, i.e., $[\mathcal{H}_{\textrm{eff}},\mathcal{Q}]=0$ holds. One
advantage of this method is that it is not model specific, i.e., the solution
for the coefficient functions and therefore for the effective Hamiltonian does
only depend on the number and type of the $T_{n}$ operators. The specific
underlying model is only needed for explicitly calculating quantities like the
ground-state energy or the gap, which corresponds to normal-ordering the
effective pCUT Hamiltonian $\mathcal{H}_{\textrm{eff}}$.

If a Hamiltonian is not Hermitian but has a complex perturbation parameter
instead, the flow equation (9) still holds. Therefore the ansätze (10) and
(11) do not have to be altered and the calculated coefficient functions remain
unchanged. Note that in this case the infinitesimal generator $\eta(\ell)$ is
not anti-Hermitian and the transformation $U(\ell)$ not unitary but a
similarity transformation. Non-perturbative similarity transformations using
the quasi-particle generator were already applied successfully for non-
Hermitian Hamiltonians [38, 39] to describe magnetic excitations of the long-
range ordered Heisenberg model on the square lattice using the Dyson-Maleev
transformation. They are referred to as CST for continuous similarity
transformation. For our systems, the Hamiltonian can be written in both cases
as Eq. (8) with the two perturbation parameters $-g/(2J)$ and $-g^{*}/(2J)$.
In order to obtain the respective ground-state energy and one-particle gap, we
have calculated the zero and one qp block of
$\mathcal{H}_{\textrm{eff}}^{\textrm{lf}}$ for both models.

The unperturbed ground-state energy $-JN$ is non-degenerate in both cases,
thus the zero qp block is a number equal to the perturbed ground-state energy.
It is given by

$\displaystyle
E_{0}^{\textrm{lf}}=\tensor{*}[_{\textrm{lf}}]{\bra{0}}{}\mathcal{H}_{\textrm{eff}}^{\textrm{lf}}\ket{0}_{\textrm{lf}}.$
(12)

Its specific value has to be calculated for both models separately. The
results are given in the appendix A in Eq. (43) [Eq. (44)] up to 12 order for
the 1D [2D] case. For the 2D case we have applied a full graph decomposition.

The one qp block consists of the elements

$\displaystyle\bra{0}b^{\phantom{\dagger}}_{\mu}\mathcal{H}_{\textrm{eff}}^{\textrm{lf}}\,b^{\dagger}_{\nu}\ket{0}=E_{0}^{\textrm{lf}}\delta_{\mu\nu}+t_{\mu\nu},$
(13)

where $t_{\mu\nu}$ denotes the hopping amplitude from $\nu$ to $\mu$. After a
Fourier transform this block further decomposes into decoupled $2\times 2$
blocks for different momenta $k$. The two bands of the dispersion relation are
obtained by diagonalizing these blocks separately. The gap closes at $k=0$ for
both systems, so we only have to diagonalize the corresponding block and pick
the lower eigenvalue. The respective series are given in appendix A in Eq.
(45) [Eq. (46)] up to 10 order for the 1D [2D] case. The 2D results were again
obtained using a full graph decomposition.

###  III.2 High-field expansions

In the following we study the high-field (hf) limit $J\ll g$. The unperturbed
case $J=0$ corresponds to $\mathcal{H}_{\textrm{I}}=0$. This limit was already
brought up in Sec. II. Similar to before, we introduce hardcore-bosonic
operators

$\displaystyle b^{\phantom{\dagger}}_{j}$
$\displaystyle=(\sigma^{x}_{j}+\mathrm{i}\sigma^{y}_{j})/2$ $\displaystyle
b^{\dagger}_{j}$ $\displaystyle=(\sigma^{x}_{j}-\mathrm{i}\sigma^{y}_{j})/2$
(14)

and express the full Hamiltonians (2) and (4) as

$\displaystyle\begin{split}\mathcal{H}_{\textrm{TFIM}}&=-\eta
N+2g\sum_{j\in\circ}b^{\dagger}_{j}b^{\phantom{\dagger}}_{j}+2g^{*}\sum_{j\in\bullet}b^{\dagger}_{j}b^{\phantom{\dagger}}_{j}\\\
&\hphantom{{}=}-J\sum_{\langle
i,j\rangle}\left(b^{\dagger}_{i}+b^{\phantom{\dagger}}_{i}\right)\left(b^{\dagger}_{j}+b^{\phantom{\dagger}}_{j}\right),\end{split}$
(15) $\displaystyle\begin{split}\mathcal{\widetilde{H}}_{\textrm{TCF}}&=-\eta
N+2g\sum_{j\in\circ}b^{\dagger}_{j}b^{\phantom{\dagger}}_{j}+2g^{*}\sum_{j\in\bullet}b^{\dagger}_{j}b^{\phantom{\dagger}}_{j}\\\
&\hphantom{{}=}-J\frac{N}{2}-J\sum_{s}\prod_{j\in
s}\left(b^{\dagger}_{j}+b^{\phantom{\dagger}}_{j}\right).\end{split}$ (16)

When $\eta\neq 0$, the unperturbed ground state is unique in both cases,
defined by
$\sum_{j}b^{\dagger}_{j}b^{\phantom{\dagger}}_{j}\ket{0}_{\textrm{hf}}=0$ and
corresponds to the state $\ket{\Uparrow}$ introduced in Sec. II. It is not
possible to rescale these Hamiltonians such that the unperturbed parts are
expressed as counting operators. This implies that we cannot use the same
method as in the low-field limit. Instead, we use Takahashi’s perturbation
theory [33] that is based on [40].

This method is suited for Hamiltonians of the form $\mathcal{H}=H_{0}+\lambda
V$ with perturbation parameter $\lambda$. Assume that the unperturbed part has
an eigenvalue $E$ whose perturbative corrections we are interested in. Let
$P_{0}$ be a projector onto the space spanned by the unperturbed eigenvectors.
We can then construct a projector onto the space spanned by the perturbed
eigenvectors corresponding to the perturbed eigenvalue by

$\displaystyle
P=P_{0}-\sum_{k=1}^{\infty}\lambda^{k}\sum_{k_{1}+\cdots+k_{k+1}=k,k_{i}\geq
0}S^{k_{1}}V\cdots VS^{k_{k+1}}$ (17)

with $S=(1-P_{0})/(E-H_{0})$ being a projector onto the complement of the
space spanned by the unperturbed eigenvectors and $S^{0}=-P_{0}$. The
transformation that generates the effective Hamiltonian is obtained by
expanding $\Gamma=PP_{0}(P_{0}PP_{0})^{-1/2}$ in powers of $\lambda$.

The ground-state energy and gap for both models are obtained analogously to
Sec. III.1. Because of $\eta\neq 0$ the unperturbed ground-state energy is
non-degenerate and the perturbed ground-state energy can be calculated by

$\displaystyle
E_{0}^{\textrm{hf}}=\tensor{[}_{\textrm{hf}}]{\bra{0}}{}\mathcal{H}_{\textrm{eff}}^{\textrm{hf}}\ket{0}_{\textrm{hf}}.$
(18)

The series are given in appendix B in Eq. (47) [Eq. (48)] up to 10 [6] order
for the 1D [2D] case.

In the 1D case the one qp block consists of the elements

$\displaystyle\tensor{[}_{\textrm{hf}}]{\bra{0}}{}b^{\phantom{\dagger}}_{i}\mathcal{H}_{\textrm{eff}}^{\textrm{hf}}\,b^{\dagger}_{j}\ket{0}_{\textrm{hf}}=E_{0}^{\textrm{hf}}\delta_{ij}+t_{ij},$
(19)

where $t_{ij}$ denotes the hopping amplitude from $j$ to $i$. For the 2D case
the relevant excitation is composed of four spin flips on a star as all
excitations with less spin flips violate the constraint that all eigenvalues
of plaquette operators are $+1$. The respective qp block consists of the
elements

$\displaystyle\tensor{[}_{\textrm{hf}}]{\bra{0}}{}\left(\prod_{i\in
s}b^{\phantom{\dagger}}_{i}\right)\mathcal{H}_{\textrm{eff}}^{\textrm{hf}}\,\left(\prod_{j\in
s^{\prime}}b^{\dagger}_{j}\right)\ket{0}_{\textrm{hf}}=E_{0}^{\textrm{hf}}\delta_{ss^{\prime}}+t_{ss^{\prime}},$
(20)

where $t_{ss^{\prime}}$ denotes the hopping amplitude of the composite entity
from star $s^{\prime}$ to $s$. After a Fourier transform this block further
decomposes into decoupled $2\times 2$ blocks for different momenta $k$. The
two bands of the dispersion relation are obtained by diagonalizing these
blocks separately.

The energy of the elementary excitations with $k=0$ in the 1D case is given in
the appendix B in Eq. (49) up to 10 order. The gap for the 2D case is given up
to 4 order in Eq. (50). For the special case $\xi=0$ the gap is known up to
even higher orders [28] and the series up to order 10 is presented in Eq. (51)
for completeness. In the parameter region defined by $\eta=\xi$ we have
further calculated the gap explicitly up to 8 order in Eq. (52).

###  III.3 DLog Padé approximants

In the vicinity of a second-order quantum critical point $x_{\text{c}}$, the
gap typically vanishes as $\Delta\propto|x_{\text{c}}-x|^{z\nu}$ with the
dynamical critical exponent $z$ and the correlation length critical exponent
$\nu$. Similarly, the second derivative of the ground-state energy diverges
according to a power law
$\text{d}^{2}e_{0}/\text{d}x^{2}\propto|x_{\text{c}}-x|^{-\alpha}$ for
$\alpha\neq 0$ with the specific heat critical exponent $\alpha$.

Whenever we have a function $F$ with such a power law behavior, we can use the
logarithmic derivative in order to extract the value of the critical exponent
$\vartheta$ (see [41] for a general introduction of extrapolation techniques
like Padé and DLog Padé approximation). Assume that close to the critical
point $x_{\text{c}}$

$\displaystyle F(x)\approx|x_{\text{c}}-x|^{-\vartheta}A(x)$ (21)

with $A$ being analytic at $x_{\text{c}}$. This enables us to approximate
$A(x)\approx A(x_{\text{c}})\left(1+\mathcal{O}(x_{\text{c}}-x)\right)$ and to
calculate the logarithmic derivative

$\displaystyle
D(x)=\frac{\text{d}}{\text{d}x}\ln(F(x))\approx\frac{\vartheta}{x_{\text{c}}-x}\left(1+\mathcal{O}(x_{\text{c}}-x)\right)$
(22)

close to the critical point. The critical exponent is then given by the
residue of this derivative at the critical point.

If the function $F$ is not known but only its Taylor expansion about $x=0$ up
to some order $k$, we use DLog Padé approximants to estimate the critical
point and corresponding exponent. The DLog Padé approximant of $F$ of order
$(L,M)$ with $L+M=k-1$ is defined by

$\displaystyle P_{\text{DLog}}[L,M]_{F}(x)=P[L,M]_{D}(x),$ (23)

where

$\displaystyle
P[L,M]_{D}(x)=\frac{P_{L}(x)}{Q_{M}(x)}=\frac{p_{0}+p_{1}x+\cdots+p_{L}x^{L}}{1+q_{1}x+\cdots+q_{M}x^{M}}$
(24)

is the Padé approximant of order $(L,M)$ of the logarithmic derivative. The
coefficients are fixed by imposing that the Taylor series of the Padé
approximant coincides with the one of the logarithmic derivative up to order
$k-1$. A DLog Padé approximant indicates the critical point to be the root of
its denominator $Q_{M}(x)$ and the critical exponent to be its residue at the
critical point.

It is possible to obtain an approximation of the function $F$ by

$\displaystyle F(x)\approx
F(0)\exp\left(\int_{0}^{x}\mathop{}\\!\mathrm{d}y\,P_{\text{DLog}}[L,M]_{F}(y)\right).$
(25)

If the critical exponent $\vartheta$ is known we can perform a biased DLog
Padé approximation, i.e., we add the next highest order (the next even order
in case the series contains only even orders) to the Taylor series of $F$ and
take the respective prefactor $a$ such that the critical exponent is the
desired one. This prefactor can differ depending on the order $(L,M)$ of the
approximant we are interested in.

For our analysis, we calculate all possible DLog Padé approximants with
$L,M\geq 2$ and group them into families defined by constant $L-M$.
Approximants with removable root-pole pairs have to be excluded as they carry
the same information as their lower-ordered relatives and thus convey a false
sense of convergence. If the denominator of an approximant has additional
roots close to the expected one, the approximant has to be discarded as
defective as well, because the additional root may deform it. Here, we also
have to take complex roots into account and specify that additional roots are
considered to be close to the expected root $x_{\text{c}}$ if their distance
in the complex plane is smaller than $x_{\text{c}}/2$ . Finally, we might only
include families with at least two remaining members in our analysis,
depending on the total number of available approximants.

The critical exponents of the Hermitian limits of (5) are either known exactly
for the 1D TFIM or estimated with high precision for the 2D TFIM, which then
naturally translates to the dual toric code in a field. For the 1D chain, the
critical exponents are $\nu=z=1$ and the specific heat critical exponent is
$\alpha=0$, i.e., the second derivative of the ground-state energy diverges
logarithmically and not according to a power law. The best estimates for the
exponents of the 2D model are obtained via conformal bootstrap or Monte Carlo
simulations; the results given are $z=1$ [42, 43], $\nu=0.629971(4)$ [43], and
$\alpha=0.110087(12)$ [43].

##  IV Ising chain results

###  IV.1 Exact solution

The Ising chain in a Hermitian uniform transverse magnetic field is exactly
solvable [35, 34] by using a Jordan-Wigner transformation [36] and a
Bogoliubov transformation in Fourier space. The non-Hermitian generalization
(2) of this model is solvable by using the same steps, adapted to a larger
unit cell.

Using $\sigma^{\pm}_{j}=(\sigma^{x}_{j}\pm\sigma^{y}_{j})/2$, we define the
fermionic operators of the Jordan-Wigner transformation by

$\displaystyle c^{\phantom{\dagger}}_{j}$
$\displaystyle:=\left(\prod_{l=1}^{j-1}\sigma^{z}_{l}\right)\sigma^{+}_{j},$
$\displaystyle c^{\dagger}_{j}$
$\displaystyle:=\left(\prod_{l=1}^{j-1}\sigma^{z}_{l}\right)\sigma^{-}_{j}$
(26)

for $j\in\\{1,\ldots,N\\}$. Inverted, the Pauli operators are given by

$\displaystyle\sigma^{z}_{j}$
$\displaystyle=1-2\,c^{\dagger}_{j}c^{\phantom{\dagger}}_{j},$ (27)
$\displaystyle\sigma^{x}_{j}$
$\displaystyle=\prod_{l=1}^{j-1}\left(1-2\,c^{\dagger}_{l}c^{\phantom{\dagger}}_{l}\right)\left(c^{\dagger}_{j}+c^{\phantom{\dagger}}_{j}\right)$
(28)

for $j\in\\{1,\ldots,N\\}$. The Hamiltonian expressed in terms of the
fermionic operators (26) reads

$\displaystyle\begin{split}\mathcal{H}_{\textrm{TFIM}}&=-J\sum_{j=1}^{N}\left(c^{\dagger}_{j}c^{\phantom{\dagger}}_{j+1}+c^{\dagger}_{j}c^{\dagger}_{j+1}+\textrm{h.c.}\right)\\\
&\hphantom{{}=}-N\eta+2\sum_{j=1}^{N}g^{\phantom{*}}_{j}c^{\dagger}_{j}c^{\phantom{\dagger}}_{j},\end{split}$
(29)

where we have neglected a boundary term that vanishes in the thermodynamic
limit [34]. This term stems from transforming $\sigma^{x}_{N}\sigma^{x}_{1}$
and it arises from the fact that we impose $c_{N+1}=c_{1}$ which might differ
by a sign from inserting $j=N+1$ into Eq. (26).

We take the Fourier transform

$\displaystyle c^{\phantom{\dagger}}_{k,\circ}$
$\displaystyle=\sqrt{\frac{2}{N}}\sum_{j\in\circ}e^{-\mathrm{i}k(j-1)/2}c^{\phantom{\dagger}}_{j}$
(30) $\displaystyle c^{\phantom{\dagger}}_{k,\bullet}$
$\displaystyle=\sqrt{\frac{2}{N}}\sum_{j\in\bullet}e^{-\mathrm{i}kj/2}c^{\phantom{\dagger}}_{j},$
(31)

where $\circ$ contains all empty spins marked in red and $\bullet$ all filled
spins marked in green in Fig. 1(a). The Hamiltonian decomposes into a sum
$\mathcal{H}_{\textrm{TFIM}}=\sum_{k}H_{k}$ with

$\displaystyle H_{k}$
$\displaystyle=-J\left(c^{\dagger}_{k,\circ}c^{\phantom{\dagger}}_{k,\bullet}+e^{\mathrm{i}k}c^{\dagger}_{k,\bullet}c^{\phantom{\dagger}}_{k,\circ}+c^{\dagger}_{k,\circ}c^{\dagger}_{-k,\bullet}+e^{\mathrm{i}k}c^{\dagger}_{k,\bullet}c^{\dagger}_{-k,\circ}+\textrm{h.c.}\right)-2\eta+2(\eta+\mathrm{i}\xi)c^{\dagger}_{k,\circ}c^{\phantom{\dagger}}_{k,\circ}+2(\eta-\mathrm{i}\xi)c^{\dagger}_{k,\bullet}c^{\phantom{\dagger}}_{k,\bullet}.$
(32)

The transformed Hamiltonian is now block-diagonal with blocks
$h_{k}=H_{k}+H_{-k}$. These blocks are individually diagonalizable using a
generalized Bogoliubov transformation. We define vectors of fermionic
operators

$\displaystyle\mathbf{c}_{k}=\begin{pmatrix}c^{\phantom{\dagger}}_{k,\circ}&c^{\phantom{\dagger}}_{k,\bullet}&c^{\phantom{\dagger}}_{-k,\circ}&c^{\phantom{\dagger}}_{-k,\bullet}\end{pmatrix}$
(33)

to write the blocks as

$\displaystyle
h_{k}=\frac{1}{2}\begin{pmatrix}\mathbf{c}^{\dagger}_{k}&\mathbf{c}^{\phantom{\dagger}}_{k}\end{pmatrix}\begin{pmatrix}-\mathbf{A}_{k}&\mathbf{B}_{k}\\\
\mathbf{B}^{\dagger}_{k}&\mathbf{A}^{\textrm{T}}_{k}\end{pmatrix}\begin{pmatrix}\mathbf{c}^{\phantom{\dagger}}_{k}\\\
\mathbf{c}^{\dagger}_{k}\end{pmatrix},$ (34)

where the matrices $\mathbf{A}_{k}$ and $\mathbf{B}_{k}$ are the components of
the coefficient matrix. $\mathbf{A}_{k}$ is given by

$\displaystyle\begin{pmatrix}-2(\eta+\mathrm{i}\xi)&J(e^{\mathrm{i}k}+1)&0&0\\\
J(e^{-\mathrm{i}k}+1)&-2(\eta+\mathrm{i}\xi)&0&0\\\
0&0&2(\eta+\mathrm{i}\xi)&J(e^{-\mathrm{i}k}+1)\\\
0&0&J(e^{\mathrm{i}k}+1)&2(\eta-\mathrm{i}\xi)\end{pmatrix}$ (35)

and $\mathbf{B}_{k}$ by

$\displaystyle\begin{pmatrix}0&0&0&-J(e^{\mathrm{i}k}-1)\\\
0&0&J(e^{-\mathrm{i}k}-1)&0\\\ 0&-J(e^{-\mathrm{i}k}-1)&0\\\
J(e^{\mathrm{i}k}-1)&0&0&0\end{pmatrix}.$ (36)

For a Hermitian problem, the matrices $\mathbf{A}_{k}$ would be Hermitian and
the matrices $\mathbf{B}_{k}$ anti-symmetric. This would make the coefficient
matrix Hermitian and therefore unitary diagonalizable. Our matrix is not
Hermitian but normal and therefore still unitary diagonalizable. The
corresponding transformation maps the fermionic operators $c$ onto other
fermionic operators $\gamma$. The so diagonalized Hamiltonian takes the form

$\displaystyle\begin{split}\mathcal{H}_{\textrm{TFIM}}&=-\sum_{k}(\omega^{+}_{k}+\omega^{-}_{k})\\\
&+2\sum_{k}\left(\omega^{+}_{k}\gamma^{\dagger}_{k,+}\gamma^{\phantom{\dagger}}_{k,+}+\omega^{-}_{k}\gamma^{\dagger}_{k,-}\gamma^{\phantom{\dagger}}_{k,-}\right)\end{split}$
(37)

with the two bands

$\displaystyle\omega^{\pm}_{k}$
$\displaystyle=\sqrt{J^{2}+\eta^{2}-\xi^{2}\pm\sqrt{2J^{2}(\eta^{2}-\xi^{2})-4\eta^{2}\xi^{2}+2J^{2}(\eta^{2}+\xi^{2})\cos(k)}}$
(38) $\displaystyle=J\sqrt{1+r^{2}\cos(2\varphi)\pm
r\sqrt{2\cos(2\varphi)-r^{2}\sin^{2}(2\varphi)+2\cos(k)}},$ (39)

where Eq. (39) is obtained by inserting $\eta=Jr\cos(\varphi)$ and
$\xi=Jr\sin(\varphi)$ into Eq. (38). These bands can also be obtained through
linear combination of results by Li et al. [13].

###  IV.2 Phase diagram

$\displaystyle 0$$\displaystyle\pi/4$$\displaystyle\pi/2$$\displaystyle
3\pi/4$$\displaystyle\pi$$\displaystyle 5\pi/4$$\displaystyle
3\pi/2$$\displaystyle
7\pi/4$$\displaystyle\xi=0$$\displaystyle\eta=\xi$$\displaystyle\eta=0$$\displaystyle\xi=2J$$\displaystyle
r=0$$\displaystyle r=1$$\displaystyle r=\infty$

Figure 2: Phase diagram of the non-Hermitian Ising chain in a staggered
transverse magnetic field (2) as a function of $\arctan(r)$ and $\varphi$ in
polar coordinates. The ground-state energy diverges logarithmically at the
boundaries indicated by thick black lines. The yellow radial lines correspond
to paths with $\xi=0$, $\eta=\xi$, and $\eta=0$ displayed in Fig. 3 while the
blue line refers to the path with $\xi=2J$ displayed in Fig. 4.

The analytic solution to our model enables us to investigate the phase
transitions, as well as the points where energies become complex. There are
energies that are real in the whole parameter space and energies that are real
only in certain regions. The only region with real spectrum is the Hermitian
axis given by $\xi=0$.

The ground-state energy per spin is given by

$\displaystyle
e_{0}=-\frac{1}{4\pi}\int\limits_{0}^{2\pi}(\omega^{+}_{k}+\omega^{-}_{k})\mathop{}\\!\mathrm{d}k$
(40)

in the thermodynamic limit. This energy is real, because the imaginary parts
of the bands $\omega^{+}_{k}$ and $\omega^{-}_{k}$ cancel each other. Li et
al. [13] argue that the Laplacian of the ground-state energy diverges
logarithmically at $r=1$ and at $\eta=0$ for $\xi\geq J$. This indicates
second-order phase transitions at these boundaries, as depicted in Fig. 2.

$\displaystyle 0$$\displaystyle 2J$$\displaystyle\xi=0$$\displaystyle
0$$\displaystyle 2J$$\displaystyle\eta=\xi$$\displaystyle r=0$$\displaystyle
r=1$$\displaystyle r=\infty$$\displaystyle 0$$\displaystyle
2J$$\displaystyle\eta=0$$\displaystyle 0$$\displaystyle
2\eta$$\displaystyle-2\xi$$\displaystyle
2\xi$$\displaystyle-2\xi$$\displaystyle 0$$\displaystyle 2\xi$

Figure 3: Energies above the ground-state energy of elementary excitations
with momentum $k=0$ of the non-Hermitian Ising chain in a staggered transverse
magnetic field as a function of $\arctan(r)$ between $0$ and $\pi/2$. The
three plots correspond to the three yellow radial paths in Fig. 2 with $\xi=0$
(upper panel), $\eta=\xi$ (middle panel), and $\eta=0$ (lower panel). Real
parts are depicted using solid lines, imaginary parts using dotted lines. The
gap $\Delta=2\omega^{-}_{0}$ is depicted in black, the other elementary
excitation with the same momentum $2\omega^{+}_{0}$ in gray. The gap closes at
$r=1$, indicating a second-order phase transition. The energies are real if
and only if $\xi\leq J$.

$\displaystyle\eta=-\infty$$\displaystyle\eta=0$$\displaystyle\eta=\infty$$\displaystyle
0$$\displaystyle 2\xi$$\displaystyle\xi=2J$

Figure 4: Energies above the ground-state energy of elementary excitations
with momentum $k=0$ of the non-Hermitian Ising chain in a staggered transverse
magnetic field as a function of $\arctan(\eta)$ between $-\pi/2$ and $\pi/2$.
The plot corresponds to the blue horizontal line indicated by $\xi=2J$ in Fig.
2. Real parts are depicted using solid lines, imaginary parts using dotted
lines. The gap $\Delta=2\omega^{-}_{0}$ is depicted in black, the other
elementary excitation with the same momentum $2\omega^{+}_{0}$ in gray. The
energies are complex conjugated. Their real parts vanish at $\eta=0$,
indicating a second-order phase transition.

The gap $\Delta=2\omega^{-}_{0}$ is one of the energies that are real in
certain extended regions of the parameter space, namely the region defined by
$\xi\leq J$. Outside of this region it is complex conjugated to
$2\omega^{+}_{0}$, the other elementary excitation with the same momentum
$k=0$. More specifically, it is given by

$\displaystyle\omega^{\pm}_{0}=\begin{cases}\left|\eta\pm\sqrt{J^{2}-\xi^{2}}\right|,&\xi\leq
J\\\ \eta\pm\mathrm{i}\sqrt{\xi^{2}-J^{2}},&\xi\geq J.\end{cases}$ (41)

The closing behavior of this gap indicates the same phase transitions as
indicated by the ground-state energy. At $r=1$ we have $\Delta=0$ as in the
Hermitian limit. This is depicted in Fig. 3 for different $\varphi$. At
$\eta=0$ we have

$\displaystyle\Delta\big{|}_{\eta=0}=2\sqrt{J^{2}-\xi^{2}},$ (42)

which has a vanishing real part for $\xi>J$, but a non-vanishing imaginary
part. This is depicted in Fig. 4. Since complex eigenvalues of $PT$-symmetric
Hamiltonians come in conjugate pairs, these imaginary parts cancel in the
ground-state energy and the vanishing real part plays the part of a vanishing
gap.

As a consequence, the gap vanishes with an exponent $z\nu=1/2$ for $\eta=0$.
In all other situations with $\eta\neq 0$ the gap closes linearly with
$z\nu=1$ including the well known Hermitian case where the gap is exactly
given by the first-order expression $|J-\eta|$.

Now we have that the ground-state energy is real in the whole parameter space
and the gap in the region defined by $|\xi|\leq J$. The question arises
whether we can find a region with a purely real spectrum therein. This
question is answered quickly by observing that e.g. $\omega^{-}_{\pi}$ is only
real on the Hermitian axis defined by $\xi=0$. Therefore we can conclude that
there exist states with broken $PT$-symmetry as soon as we turn on the non-
Hermitian perturbation.

###  IV.3 Gauging high-order series expansions

We can now use the analytic solution to gauge our series expansion results. We
derived the gap up to 10 order about the low-field limit (45) and about the
high-field limit (49). This order is similar to the maximal order for the 2D
case discussed below where the exact solution is not available. Note that in
principle we can reach much higher orders for the non-Hermitian Ising chain by
simply determining the Taylor series of the analytic expression.

$\displaystyle 0$$\displaystyle 2J$$\displaystyle\xi=0$(0, 1)$\displaystyle
0$$\displaystyle 2J$$\displaystyle\eta=\xi$$\displaystyle 0$$\displaystyle
1/10$$\displaystyle 1/6$$\displaystyle
1/\textrm{order}$$\displaystyle{1.000}$$\displaystyle{1.005}$$\displaystyle
r_{\textrm{c}}$(4, 4)$\displaystyle r=0$$\displaystyle r=1$$\displaystyle
r=\infty$$\displaystyle 0$$\displaystyle
2J$$\displaystyle\eta=0$$\displaystyle 0$$\displaystyle 1/10$$\displaystyle
1/6$$\displaystyle
1/\textrm{order}$$\displaystyle{1.0}$$\displaystyle{1.2}$$\displaystyle
r_{\textrm{c}}$(1, 2)

Figure 5: Real part of the gap for the non-Hermitian Ising chain in a
staggered transverse magnetic field as a function of $\arctan(r)$ between $0$
and $\pi/2$. The three plots correspond to the three yellow radial paths with
$\xi=0$ (upper panel), $\eta=\xi$ (middle panel), and $\eta=0$ (lower panel)
in Fig. 2. The analytic solution is depicted in black; the series expansions
(except for the 0 order) up to 10 order in red. Selected DLog Padé
approximants (25) are plotted in thick lines that are colored according to the
colors in Fig. 6. The roots of the series expansion (up to orders 4, 6, 8, 10
from right to left) are plotted as a function of $1/\textrm{order}$ in inset
plots in red. The black markers indicate the critical point $r_{\text{c}}=1$
in both cases.

Except for $\eta=0$ and $\xi=0$, the vanishing of the gap is only accessible
from the low-field limit due to the presence of an exceptional point at
$r>r_{\text{c}}$ as can be seen exemplary in the middle panel of FIG. 3. This
is different for $\eta=0$, where the two points coincide. The analytic low-
field gap, as well as its series expansion order by order, is plotted in Fig.
5 for exemplary regions of the parameter space. The Hermitian case $\xi=0$ is
exact already in first order, all other cases are only approximated by our
expansions. The critical point is predicted more accurately the closer we are
to the Hermitian limit. Using the bare series in tenth order, it is indicated
at $r_{\text{c}}=1.0004$ for $\eta=\xi$ and at $r_{\text{c}}=1.08$ for
$\eta=0$.

$\displaystyle{0.980}$$\displaystyle{0.985}$$\displaystyle{0.990}$$\displaystyle{0.995}$$\displaystyle{1.000}$$\displaystyle
r_{\textrm{c}}$$\displaystyle{0.995}$$\displaystyle{1.000}$$\displaystyle
r_{\textrm{c}}$01-12-23-345-5$\displaystyle{4}$$\displaystyle{5}$$\displaystyle{6}$$\displaystyle{7}$$\displaystyle{8}$$\displaystyle{9}$order$\displaystyle{0.88}$$\displaystyle{0.90}$$\displaystyle{0.92}$$\displaystyle{0.94}$$\displaystyle{0.96}$$\displaystyle{0.98}$$\displaystyle{1.00}$$\displaystyle
z\nu$$\displaystyle{0.96}$$\displaystyle{0.98}$$\displaystyle{1.00}$$\displaystyle
z\nu$

Figure 6: Critical point $r_{\text{c}}$ (upper panel) and critical exponent
$z\nu$ (lower panel) of the non-Hermitian Ising chain in a staggered
transverse magnetic field as indicated by DLog Padé approximants. The
parameter region is defined by $\eta=\xi$. Every family has a color assigned
and is labeled by their difference between numerator and denominator degree.
The averaged values and the respective sample standard deviations are
displayed in the inset plots. The gray lines are obtained by including single
family members, the black ones by considering them defective. In the latter
case, there is only one approximant left in 8 order, which is indicated by a
circle. In 9 order, the critical point [critical exponent] is found at
$r_{\text{c}}=0.9993(8)$ [$z\nu=0.991(9)$] when including single family
members and at $r_{\text{c}}=0.99967(19)$ [$z\nu=0.9955(23)$] when considering
them defective. The dashed black lines indicate the analytical findings
$r_{\text{c}}=1$ and $z\nu=1$.

The values are expected to become more accurate upon using DLog Padé
approximants. With Eq. (25), selected approximants are depicted in Fig. 5 as
well. Remarkably, not only the limiting case $\xi=0$ is described exactly by
these approximants, but also the purely imaginary case $\eta=0$ where the gap
closes as a square root. This is similar to the extrapolation of the mean-
field series in the limit of large spatial dimensions for the transverse-field
Ising model on the hypercubic lattice [25]. In fact, one can construct
$\Delta=2\sqrt{J^{2}-\xi^{2}}-2\eta$ exactly by adding the exact
extrapolations on the real and imaginary axis. The same holds for the
imaginary part of the gap about the high-field limit.

For any given order, there are multiple different DLog Padé approximants to
consider. For $\eta=\xi$, their respective critical points are plotted in the
upper panel of Fig. 6. Per order the predicted values are obtained by
averaging over all results by non-defective approximants. If single family
members are considered defective, the predicted value for the critical point
is given by $r_{\text{c}}=0.99967(19)$, if they are included it is given by
$r_{\text{c}}=0.9993(8)$. Both values are accurate to the same order of
magnitude as the value indicated by the bare 10 order series expansion
$r_{\text{c}}=1.0004$. DLog Padé approximants are especially well suited to
obtain critical exponents. For $\eta=\xi$, the critical exponent $z\nu$ is
plotted in the lower panel of Fig. 6. If single family members are considered
defective, the predicted value is given by $z\nu=0.9955(23)$, if they are
included it is given by $z\nu=0.991(9)$. The accuracy is lower than the
accuracy of the predicted critical point, which is expected. Nevertheless, it
is still convincing.

##  V Toric code results

$\displaystyle 0$$\displaystyle J$$\displaystyle 2J$$\displaystyle\xi=0$(5,
4)(4, 4)(3, 4)$\displaystyle 0$$\displaystyle J$$\displaystyle
2J$$\displaystyle\eta=\xi$$\displaystyle 0$$\displaystyle 1/8$$\displaystyle
1/4$$\displaystyle{0.45}$$\displaystyle{0.50}$$\displaystyle
r_{\textrm{c}}$$\displaystyle 1/\textrm{order}$$\displaystyle
r=0$$\displaystyle r=1$$\displaystyle r=\infty$$\displaystyle 0$$\displaystyle
J$$\displaystyle 2J$$\displaystyle\eta=0$$\displaystyle 0$$\displaystyle
4\eta$$\displaystyle 8\eta$$\displaystyle 0$$\displaystyle
4\eta$$\displaystyle 8\eta$

Figure 7: Gap of the non-Hermitian toric code in a staggered parallel magnetic
field as a function of $\arctan(r)$ between $0$ and $\pi/2$. The three plots
correspond to the three different parameter regions with $\xi=0$ (upper
panel), $\eta=\xi$ (middle panel), and $\eta=0$ (lower panel). The series
expansions (except for the 0 order) up to 10 order are depicted in red around
the low-field limit and in blue around the high-field limit, with larger
orders having darker colors. The high-field expansion for $\eta=\xi$ is only
available up to 8 order. Selected DLog Padé approximants (25) are plotted
using thick lines with colors matching the respective families. The roots of
the high-field series expansion for $\eta=\xi$ (up to orders 4, 6, 8 from
right to left) are plotted as a function of $1/\textrm{order}$ in the inset
plot in blue, those of the respective low-field expansion (up to orders 4, 5,
6, 7, 8, 9, 10 from right to left) in red.

In contrast to the non-Hermitian Ising chain, there is no analytic solution
for the toric code in a non-Hermitian staggered field whose low-energy physics
corresponds to the dual non-Hermitian transverse-field Ising model on the
square lattice with non-Hermitian Ising interactions. However, the Hermitian
case is well studied by different techniques including high-order series
expansions [44, 27, 45, 29] (for the dual transverse-field Ising model on the
square lattice see [22, 23, 46, 43]). One finds a second-order phase
transition in the 3D* [47] universality class with $r_{\text{c}}=0.32847(4)$
[46] and $z\nu=0.629971(4)$ [43]. Here we are mainly interested in how the
non-Hermitian staggered field influence this quantum critical breakdown of the
topological phase. To this end we use the series expansion of the gap. The
expansion about the low-field limit is given up to 10 order in Eq. (46). The
expansion about the high-field limit (50) is given up to 10 order in Eq. (51)
for $\xi=0$ and up to 8 order in Eq. (52) for $\eta=\xi$. We note that the
series in the Hermitian limits agree with the ones known from the literature
[22, 23]. All of these series expansions are plotted in Fig. 7 order by order
together with some selected DLog Padé approximants about both limits.

The Hermitian case $\xi=0$ is well studied by series expansions. We
nevertheless summarize the main findings with the current perturbative order
which allows a good comparison to the non-Hermitian case. DLog Padé
approximation of the low-field gap results in 8 order in
$r_{\text{c}}=0.328905(19)$ if single family members are included [if they are
considered defective one finds $r_{\text{c}}=0.3288978(10)$]. Here, we use the
8 order results as there is only one non-defective approximant in 9 order. The
accuracy is expected to even increase for higher orders, as the families
converge towards the literature value $r_{\text{c}}=0.32847(4)$ [46].
Extracting the corresponding critical exponent of the gap closing yields
$z\nu=0.6471(4)$ [$z\nu=0.64698(18)$] when including [excluding] single family
members. Turning to the gap in the high-field limit, one observes that the gap
series is alternating and that it consists only of even powers in $r^{-1}$,
which is equivalent to a series in $r^{-2}$. As a consequence, only three of
the DLog Padé approximants are not defective (only two in ninth order). Taking
the average and sample standard deviation yields a critical point of
$r_{\text{c}}=0.33394(18)$ and a critical exponent of $z\nu=0.6354(23)$ in 9
order. The result for the critical point fits the literature value, even
though it is less accurate than the low-field results (the better agreement
for the critical exponent is taken as a coincidence). In conclusion, the high-
order series expansion does therefore well capture the second-order quantum
phase transition between the low-field topological and the high-field
polarized phase in a quantitative fashion. This is in particular true for the
low-field expansion while the accuracy of the high-field expansion is less
good.

$\displaystyle{0.400}$$\displaystyle{0.405}$$\displaystyle{0.410}$$\displaystyle{0.415}$$\displaystyle{0.420}$$\displaystyle{0.425}$$\displaystyle
r_{\textrm{c}}$$\displaystyle{7}$$\displaystyle{8}$$\displaystyle{9}$$\displaystyle{0.405}$$\displaystyle{0.410}$$\displaystyle
r_{\textrm{c}}$01-12-2$\displaystyle{4}$$\displaystyle{5}$$\displaystyle{6}$$\displaystyle{7}$$\displaystyle{8}$$\displaystyle{9}$order$\displaystyle{0.55}$$\displaystyle{0.60}$$\displaystyle{0.65}$$\displaystyle{0.70}$$\displaystyle{0.75}$$\displaystyle{0.80}$$\displaystyle
z\nu$$\displaystyle{7}$$\displaystyle{8}$$\displaystyle{9}$$\displaystyle{0.60}$$\displaystyle{0.65}$$\displaystyle
z\nu$3-345-5

Figure 8: Critical point $r_{\text{c}}$ (upper panel) and critical exponent $z\nu$ (lower panel) of the non-Hermitian toric code in a non-Hermitian staggered parallel magnetic field for $\eta=\xi$ obtained by DLog Padé approximants of the low-field gap as a function of the order. Every family of approximants has a distinct color and is labeled by their difference between numerator and denominator degree. The averaged values and the respective sample standard deviation are displayed in the inset plots. The gray line is obtained by including single family members, the black one by considering them defective. In the latter case, there is only one approximants left in 7 order, which is indicated by a circle. In 9 order, the critical point [exponent] is found at $r_{\text{c}}=0.405(3)$ [$z\nu=0.64(4)$] when including single family members and $r_{\text{c}}=0.40559(5)$ [$z\nu=0.6466(6)$] when considering them defective. The dashed black lines indicate the literature value $z\nu=0.629971(4)$ [43] for the exponent of the Hermitian model. DLog Padé | Critical point | Critical exponent
---|---|---
$(3,2)$ | $r_{\text{c}}=0.38$ | $z\nu=5.2$
$(3,4)$ | $r_{\text{c}}=0.51$ | $z\nu=1.1$
$(5,2)$ | $r_{\text{c}}=0.58$ | $z\nu=0.4$
TAB. 1: Critical point and exponent of the non-Hermitian toric code in a staggered parallel magnetic field for $\eta=\xi$ as indicated by DLog Padé approximants about the high-field limit. Biased DLog Padé | Critical point | Prefactor
---|---|---
$(3,6)$ | $r_{\text{c}}=0.5511$ | $a=0.011453466$
$(5,4)$ | $r_{\text{c}}=0.5460$ | $a=0.00536516$
$(7,2)$ | $r_{\text{c}}=0.5469$ | $a=0.00609564$
TAB. 2: Critical point and exponent of the non-Hermitian toric code in a
staggered parallel magnetic field for $\eta=\xi$ as indicated by biased DLog
Padé approximants about the high-field limit. The bias is taken such that the
critical exponent is $z\nu=0.629971$. The critical point is found at
$r_{\text{c}}=0.5480(23)$.

Next we turn on the non-Hermitian staggered field and focus on $\eta=\xi$.
Estimates of the critical point $r_{\text{c}}$ and critical exponent $z\nu$
obtained by DLog Padé approximation of the low-field gap are shown in Fig. 8.
One finds $r_{\text{c}}=0.405(3)$ with the highest order if single family
members are included (at $r_{\text{c}}=0.40559(5)$ if they are considered
defective). The critical value is therefore shifted to larger values of $r$ in
comparison to the Hermitian case. The series expansion about the high-field
limit consists again only of even powers in $r^{-1}$. Except for the one
corresponding to $r^{-2}$, all coefficients of this series are negative. The
bare order 8 series indicates a critical point at $r_{\text{c}}\approx 0.51$.
As for the Hermitian case, there are only three non-defective DLog Padé
approximants; their respectively indicated critical points and exponents are
listed in Tab. 1. Taking average and sample standard deviation yields a
critical point of $r_{\text{c}}=0.55(4)$ using the highest order. This result
suggests that there could be at least two critical points in contrast to the
Hermitian case. However, this result is to be handled with care, because the
extrapolation of the high-field gap series does not work as well as the one
for the low-field gap. This is also apparent because the critical exponents
differ largely. We can, however, gain further insights by performing DLog Padé
approximation biased with the 3D Ising critical exponent $z\nu=0.629971(4)$,
i.e., assuming no change in the universality class for finite non-Hermitian
staggered fields. Although this must not be the correct exponent, the quality
of the extrapolation is convincing. Because the series consists only of even
powers in $r^{-1}$, the next higher order we add has to be even as well. More
specifically, we add a term of the form $ar^{-10}$. The respective prefactors
$a$ as well as corresponding critical points are stated in Tab. 2. The listed
averaged critical point is given by $r_{\text{c}}=0.5480(23)$. The critical
points of the biased approximants match far better than those of the unbiased
approximants, which supports our assumed critical exponent and the potential
presence of an intermediate regime between the low-field and the high-field
phase. If we instead perform a DLog Padé approximation biased with the
critical point of the low-field expansion $r_{\text{c}}=0.405$, i.e., assuming
no intermediate phase, we obtain unphysical critical exponents. This further
supports the potential presence of such an intermediate phase.

DLog Padé | Critical point | Critical exponent
---|---|---
$(3,2)$ | — | —
$(3,4)$ | $r_{\text{c}}=1.17$ | $z\nu=1.14$
$(5,2)$ | — | —
$(3,6)$ | $r_{\text{c}}=0.89$ | $z\nu=0.42$
$(5,4)$ | $r_{\text{c}}=0.72$ | $z\nu=0.12$
$(7,2)$ | — | —
TAB. 3: Critical point and exponent of the non-Hermitian toric code in a
staggered parallel magnetic field for $\eta=0$ as indicated by DLog Padé
approximants about the low-field limit. The approximants that do not close are
indicated by dashes.

For the purely anti-Hermitian staggered field ($\eta=0$), the series expansion
of the low-field gap consists only of even powers in $r$ and is alternating.
As a consequence, the extrapolation does not work in a convincing manner. Most
of the DLog Padé approximants are defective and half of those that are not
defective do not close. The associated critical points and exponents are
listed in Tab. 3. If one only considers the approximants that close, the
critical point is indicated at $r_{\text{c}}=0.80(9)$ using the highest order.
If one considers all approximants, it is not certain that the gap closes at
all. Furthermore, no large-field expansion is available for this case.

$\displaystyle\varphi=0$$\displaystyle\varphi=\pi/4$$\displaystyle\varphi=\pi/2$$\displaystyle
r=0$$\displaystyle r=1$$\displaystyle r=\infty$$\displaystyle
0$$\displaystyle\pi/4$$\displaystyle\pi/2$$\displaystyle\varphi$$\displaystyle{0.0}$$\displaystyle{0.5}$$\displaystyle{1.0}$$\displaystyle
z\nu$

Figure 9: Phase diagram of the non-Hermitian toric code in a staggered
parallel magnetic field (3) obtained from DLog Padé approximants of the low-
field gap. The radial plot displays the indicated critical point as a function
of $\arctan(r)$ and $\varphi$ in polar coordinates. The inset plot displays
the corresponding critical exponent as a function of $\varphi$. The points and
errorbars are obtained by averaging over results by selected approximants and
taking the corresponding sample standard deviation. For every point, selected
are either all non-defective 9 order approximants or, if there is only one in
9 order, all non-defective 8 order approximants. The critical point shifts to
higher values of $r$ if the anti-Hermitian part increases. For large anti-
Hermitian parts it is not certain that the gap closes at all. The displayed
crosses are the values from Tab. 3. The critical exponents are close to the
3D* Ising critical exponent $z\nu=0.629971(4)$ [43] that is indicated by the
dashed line.

We finally discuss the extension of the topological phase in the full
parameter space which is illustrated in Fig. 9 using DLog Padé approximation
of the low-field gap series. The latter works reliably except for the regime
of small $\eta$ as already discussed for the purely anti-Hermitian staggered
field. We observe that the critical point shifts to higher values of $r$ for
increasing $\xi$. The robustness of the topological phase therefore increases
continuously with the non-Hermiticity of the staggered field. The
corresponding averaged critical exponents are illustrated in Fig. 9 as well.
As long as the approximants are reliable, the deduced critical exponents are
almost independent of $\varphi$ and close to the 3D* Ising critical exponent
$z\nu=0.629971(4)$ [43]. We therefore conjecture that the critical breakdown
of the topological phase remains in the 3D* Ising universality class in this
parameter regime and can be described by the condensation of anyonic quasi-
particles.

##  VI Conclusions

In this work we have extended high-order linked-cluster expansions like
perturbative continuous unitary transformations and Takahashi perturbation
theory to non-Hermitian quantum spin models with $PT$-symmetry. In practice,
we have reached similar maximal perturbative orders for the ground-state
energy per site and the energy gap of elementary excitations as for the more
conventional Hermitian counterparts. Here we have considered the one-
dimensional Ising chain and the two-dimensional toric code in a non-Hermitian
staggered field. For both systems we have exploited dual descriptions in terms
of non-Hermitian staggered Ising interactions in an uniform transverse field.

For the one-dimensional non-Hermitian Ising chain we have demonstrated an
alternative analytic solution in full agreement with the one given in Ref. 13.
Its ground-state phase diagram can be characterized by logarithmic
singularities of the Laplacian of the ground-state energy and the gap closing
of the low-field gap in the symmetry-broken phase. In particular, the gap
closing of the low-field gap can be quantitatively described by extrapolating
the high-order gap series.

For the two-dimensional toric code in a non-Hermitian staggered field we study
the quantum robustness of the topologically ordered phase by the gap closing
of the low-field gap. We find that the well-known second-order quantum phase
transition of the toric code in an uniform field extends into a large portion
of the non-Hermitian parameter space. However, the series expansions become
unreliable for a dominant anti-Hermitian field. Interestingly, the analysis of
the high-field gap reveals the potential presence of an intermediate region,
which, however, deserves further studies in the future. It would be further
interesting to use high-order series expansions also for other perturbed non-
Hermitian toric codes [17, 18, 19] where the correspondence between bulk
quasi-particles and topologically protected degenerate ground states breaks
down. At this point we stress that we have concentrated here on applying high-
order series expansions to extract spectral properties like the ground-state
energy and elementary gaps in non-Hermitian quantum spin systems. Extensions
towards other physical quantities are certainly possible and interesting. In
particular, entanglement measures like the entanglement entropy would directly
probe topological quantum order.

## Acknowledgements

KPS acknowledges financial support by the German Science Foundation (DFG)
through the grant SCHM 2511/11-1. LL and KPS acknowledge support by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project-ID
429529648 – TRR 306 QuCoLiMa (”Quantum Cooperativity of Light and Matter”).

##  Appendix A Low-field expansions

The ground-state energy per spin of the Ising chain in a non-Hermitian
staggered transverse field (2) is given by

$\displaystyle e_{0}^{(12,\eta,\xi)}$
$\displaystyle=-J-\frac{1}{4}\frac{\eta^{2}-\xi^{2}}{J}-\frac{1}{64}\frac{\eta^{4}+\xi^{4}}{J^{3}}-\frac{5}{32}\frac{\eta^{2}\xi^{2}}{J^{3}}-\frac{1}{256}\frac{\eta^{6}-\xi^{6}}{J^{5}}+\frac{7}{256}\frac{\eta^{4}\xi^{2}-\eta^{2}\xi^{4}}{J^{5}}-\frac{25}{16384}\frac{\eta^{8}+\xi^{8}}{J^{7}}$
$\displaystyle\hphantom{{}=}-\frac{129}{4096}\frac{\eta^{6}\xi^{2}+\eta^{2}\xi^{6}}{J^{7}}-\frac{171}{8192}\frac{\eta^{4}\xi^{4}}{J^{7}}-\frac{49}{65536}\frac{\eta^{10}-\xi^{10}}{J^{9}}+\frac{781}{65536}\frac{\eta^{8}\xi^{2}-\eta^{2}\xi^{8}}{J^{9}}-\frac{33}{32768}\frac{\eta^{6}\xi^{4}-\eta^{4}\xi^{6}}{J^{9}}$
(43)
$\displaystyle\hphantom{{}=}-\frac{441}{1048576}\frac{\eta^{12}+\xi^{12}}{J^{11}}-\frac{7631}{524288}\frac{\eta^{10}\xi^{2}+\eta^{2}\xi^{10}}{J^{11}}-\frac{18551}{1048576}\frac{\eta^{8}\xi^{4}+\eta^{4}\xi^{8}}{J^{11}}-\frac{7241}{262144}\frac{\eta^{6}\xi^{6}}{J^{11}}.$

The ground-state energy per spin of the toric code in a non-Hermitian
staggered field (4) is given by

$\displaystyle e_{0}^{(12,\eta,\xi)}$
$\displaystyle=-J-\frac{1}{4}\frac{\eta^{2}-\xi^{2}}{J}-\frac{15}{64}\frac{\eta^{4}+\xi^{4}}{J^{3}}-\frac{11}{32}\frac{\eta^{2}\xi^{2}}{J^{3}}-\frac{147}{256}\frac{\eta^{6}-\xi^{6}}{J^{5}}-\frac{131}{256}\frac{\eta^{4}\xi^{2}-\eta^{2}\xi^{4}}{J^{5}}-\frac{18003}{8192}\frac{\eta^{8}+\xi^{8}}{J^{7}}$
(44)
$\displaystyle\hphantom{{}=}-\frac{1699}{512}\frac{\eta^{6}\xi^{2}+\eta^{2}\xi^{6}}{J^{7}}-\frac{9981}{4096}\frac{\eta^{4}\xi^{4}}{J^{7}}-\frac{5420775}{524288}\frac{\eta^{10}-\xi^{10}}{J^{9}}-\frac{29368399}{1572864}\frac{\eta^{8}\xi^{2}-\eta^{2}\xi^{8}}{J^{9}}-\frac{6703565}{786432}\frac{\eta^{6}\xi^{4}-\eta^{4}\xi^{6}}{J^{9}}$
$\displaystyle\hphantom{{}=}-\frac{33446240377}{603979776}\frac{\eta^{12}+\xi^{12}}{J^{11}}-\frac{108757489433}{905969664}\frac{\eta^{10}\xi^{2}+\eta^{2}\xi^{10}}{J^{11}}-\frac{6105915487}{67108864}\frac{\eta^{8}\xi^{4}+\eta^{4}\xi^{8}}{J^{11}}-\frac{23451561631}{452984832}\frac{\eta^{6}\xi^{6}}{J^{11}}.$

The gap of the Ising chain in a non-Hermitian staggered transverse field (2)
is given by

$\displaystyle\Delta^{(10,\eta,\xi)}$
$\displaystyle=2J-2\eta-\frac{\xi^{2}}{J}-\frac{1}{4}\frac{\xi^{4}}{J^{3}}-\frac{1}{8}\frac{\xi^{6}}{J^{5}}-\frac{5}{64}\frac{\xi^{8}}{J^{7}}-\frac{7}{128}\frac{\xi^{10}}{J^{9}}.$
(45)

The gap of the toric code in a non-Hermitian staggered field (4) is given by

$\displaystyle\Delta^{(10,\eta,\xi)}$
$\displaystyle=2J-4\eta-2\frac{\eta^{2}}{J}-2\frac{\xi^{2}}{J}-3\frac{\eta^{3}}{J^{2}}-3\frac{\eta\xi^{2}}{J^{2}}-\frac{9}{2}\frac{\eta^{4}}{J^{3}}-3\frac{\eta^{2}\xi^{2}}{J^{3}}+\frac{3}{2}\frac{\xi^{4}}{J^{3}}-11\frac{\eta^{5}}{J^{4}}-\frac{23}{2}\frac{\eta^{3}\xi^{2}}{J^{4}}-\frac{1}{2}\frac{\eta\xi^{4}}{J^{4}}-\frac{2625}{128}\frac{\eta^{6}}{J^{5}}$
$\displaystyle\hphantom{{}=}-\frac{3325}{128}\frac{\eta^{4}\xi^{2}}{J^{5}}-\frac{1463}{128}\frac{\eta^{2}\xi^{4}}{J^{5}}-\frac{763}{128}\frac{\xi^{6}}{J^{5}}-\frac{14771}{256}\frac{\eta^{7}}{J^{6}}-\frac{21595}{256}\frac{\eta^{5}\xi^{2}}{J^{6}}-\frac{8733}{256}\frac{\eta^{3}\xi^{4}}{J^{6}}-\frac{1909}{256}\frac{\eta\xi^{6}}{J^{6}}-\frac{940739}{8192}\frac{\eta^{8}}{J^{7}}$
$\displaystyle\hphantom{{}=}-\frac{180361}{1024}\frac{\eta^{6}\xi^{2}}{J^{7}}-\frac{227827}{4096}\frac{\eta^{4}\xi^{4}}{J^{7}}+\frac{557}{16}\frac{\eta^{2}\xi^{6}}{J^{7}}+\frac{238689}{8192}\frac{\xi^{8}}{J^{7}}-\frac{11472297}{32768}\frac{\eta^{9}}{J^{8}}-\frac{1274661}{2048}\frac{\eta^{7}\xi^{2}}{J^{8}}-\frac{4894637}{16384}\frac{\eta^{5}\xi^{4}}{J^{8}}$
$\displaystyle\hphantom{{}=}+\frac{55859}{4096}\frac{\eta^{3}\xi^{6}}{J^{8}}+\frac{1313867}{32768}\frac{\eta\xi^{8}}{J^{8}}-\frac{287258435}{393216}\frac{\eta^{10}}{J^{9}}-\frac{183060165}{131072}\frac{\eta^{8}\xi^{2}}{J^{9}}-\frac{179805973}{196608}\frac{\eta^{6}\xi^{4}}{J^{9}}-\frac{86807677}{196608}\frac{\eta^{4}\xi^{6}}{J^{9}}$
$\displaystyle\hphantom{{}=}-\frac{47455265}{131072}\frac{\eta^{2}\xi^{8}}{J^{9}}-\frac{66440327}{393216}\frac{\xi^{10}}{J^{9}}.$
(46)

##  Appendix B High-field expansions

The ground-state energy per spin of the Ising chain in a non-Hermitian
staggered transverse field (2) is given by

$\displaystyle\begin{split}e_{0}^{(10,J)}&=-\eta-\frac{J^{2}}{4\eta}-\frac{J^{4}}{64\eta^{3}}\frac{\eta^{2}-3\xi^{2}}{\eta^{2}+\xi^{2}}-\frac{J^{6}}{256\eta^{5}}\frac{\eta^{2}+5\xi^{2}}{\eta^{2}+\xi^{2}}-\frac{J^{8}}{16384\eta^{7}}\frac{25\eta^{6}-269\eta^{4}\xi^{2}-405\eta^{2}\xi^{4}-175\xi^{6}}{(\eta^{2}+\xi^{2})^{3}}\\\
&\hphantom{{}=}-\frac{J^{10}}{65536\eta^{9}}\frac{49\eta^{6}+715\eta^{4}\xi^{2}+1043\eta^{2}\xi^{4}+441\xi^{6}}{(\eta^{2}+\xi^{2})^{3}}.\end{split}$
(47)

The ground-state energy per spin of the toric code in a non-Hermitian
staggered field (4) is given by

$\displaystyle e_{0}^{(6,J)}$
$\displaystyle=-\eta-\frac{J}{2}-\frac{J^{2}}{16\eta}-\frac{J^{4}}{1024\eta^{3}}\frac{3\eta^{2}-5\xi^{2}}{9\eta^{2}+\xi^{2}}-\frac{J^{6}}{32768\eta^{5}}\frac{54\eta^{6}+9\eta^{4}\xi^{2}+368\eta^{2}\xi^{4}+29\xi^{6}}{(9\eta^{2}+\xi^{2})^{2}(4\eta^{2}+\xi^{2})}.$
(48)

The energy of the elementary excitations with $k=0$ of the Ising chain in a
non-Hermitian staggered transverse field (2) is given by

$\displaystyle\Delta_{\pm}^{(10,J)}$
$\displaystyle=2(\eta\pm\textrm{i}\xi)\pm\frac{1}{\textrm{i}\xi}J^{2}\pm\frac{1}{4\textrm{i}\xi^{3}}J^{4}\pm\frac{1}{8\textrm{i}\xi^{5}}J^{6}\pm\frac{5}{64\textrm{i}\xi^{7}}J^{8}\pm\frac{7}{128\textrm{i}\xi^{9}}J^{10}$
(49)

for $\xi\neq 0$. The gap of the toric code in a non-Hermitian staggered field
(4) is given by

$\displaystyle\Delta^{(4,J)}=\begin{cases}\displaystyle
8\eta+\frac{J^{2}}{4\eta}\frac{-3\eta^{2}+5\xi^{2}}{\eta^{2}+\xi^{2}}+\frac{J^{4}}{768\eta^{3}}\frac{387\eta^{8}-5232\eta^{6}\xi^{2}-258\eta^{4}\xi^{4}-936\eta^{2}\xi^{6}-153\xi^{8}}{(\eta^{2}+\xi^{2})^{3}(9\eta^{2}+\xi^{2})},&\eta\geq\xi\\\
\displaystyle
8\eta+\frac{J^{2}}{4\eta}+\frac{J^{4}}{768\eta^{3}}\frac{-45\eta^{8}-96\eta^{6}\xi^{2}-114\eta^{4}\xi^{4}-120\eta^{2}\xi^{6}-57\xi^{8}}{(\eta^{2}+\xi^{2})^{3}(9\eta^{2}+\xi^{2})},&\eta<\xi.\end{cases}$
(50)

For $\xi=0$, this gap is given up to higher orders by

$\displaystyle\Delta^{(10,J)}\big{|}_{\xi=0}=8\eta-\frac{3}{4}\frac{J^{2}}{\eta}+\frac{43}{768}\frac{J^{4}}{\eta^{3}}-\frac{19993}{884736}\frac{J^{6}}{\eta^{5}}+\frac{82873487}{10192158720}\frac{J^{8}}{\eta^{7}}-\frac{1901437203257}{587068342272000}\frac{J^{10}}{\eta^{9}}$
(51)

and for $\eta=\xi$ by

$\displaystyle\Delta^{(8,J)}\big{|}_{\eta=\xi}=8\eta+\frac{1}{4}\frac{J^{2}}{\eta}-\frac{129}{1280}\frac{J^{4}}{\eta^{3}}-\frac{100511}{12288000}\frac{J^{6}}{\eta^{5}}-\frac{7866337483}{52140441600000}\frac{J^{8}}{\eta^{7}}.$
(52)

## References

  * Bender [2007] C. M. Bender, Rep. Prog. Phys. 70, 947 (2007).
  * Bender and Boettcher [1998] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998).
  * El-Ganainy  _et al._ [2018] R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nature Phys. 14, 11 (2018).
  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys 69, 249 (2020).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Rev. Mod. Phys. 93, 015005 (2021).
  * Wang  _et al._ [2020] C. Wang, M.-L. Yang, C.-X. Guo, X.-M. Zhao, and S.-P. Kou, EPL 128, 41001 (2020).
  * Yang  _et al._ [2020] M.-L. Yang, H. Wang, C.-X. Guo, X.-R. Wang, G. Sun, and S.-P. Kou, arXiv:2006.10278 [cond-mat.str-el] (2020).
  * Korff [2008] C. Korff, J. Phys. A: Math. Theor. 41, 295206 (2008).
  * Castro-Alvaredo and Fring [2009] O. A. Castro-Alvaredo and A. Fring, J. Phys. A: Math. Theor. 42, 465211 (2009).
  * Albertini  _et al._ [1996] G. Albertini, S. R. Dahmen, and B. Wehefritz, Journal of Physics A: Mathematical and General 29, L369 (1996).
  * Liu  _et al._ [2021] Y.-G. Liu, L. Xu, and Z. Li, J. Phys.: Condens. Matter 33, 295401 (2021).
  * Zhang and Song [2013] X. Z. Zhang and Z. Song, Phys. Rev. A 88, 042108 (2013).
  * Li  _et al._ [2014] C. Li, G. Zhang, X. Z. Zhang, and Z. Song, Phys. Rev. A 90, 012103 (2014).
  * Leinaas and Myrheim [1977] J. M. Leinaas and J. Myrheim, Nuovo Cim. B 37, 1 (1977).
  * Wilczek [1982] F. Wilczek, Phys. Rev. Lett. 48, 1144 (1982).
  * Kitaev [2003] A. Kitaev, Ann. Phys. 303, 2 (2003).
  * Guo  _et al._ [2020a] C.-X. Guo, X.-R. Wang, C. Wang, and S.-P. Kou, Phys. Rev. B 101, 144439 (2020a).
  * Guo  _et al._ [2020b] C.-X. Guo, X.-R. Wang, and S.-P. Kou, EPL 131, 27002 (2020b).
  * Shackleton and Scheurer [2020] H. Shackleton and M. S. Scheurer, Phys. Rev. Research 2, 033022 (2020).
  * Matsumoto  _et al._ [2020] N. Matsumoto, K. Kawabata, Y. Ashida, S. Furukawa, and M. Ueda, Phys. Rev. Lett. 125, 260601 (2020).
  * Oitmaa  _et al._ [2006] J. Oitmaa, C. Hamer, and W. Zheng,  _Series Expansion Methods for Strongly Interacting Lattice Models: Series expansions for lattice gauge models_ (Cambridge University Press, Cambridge, 2006).
  * He  _et al._ [1990] H. X. He, C. J. Hamer, and J. Oitmaa, J. Phys. A: Math. Gen. 23, 1775 (1990).
  * Weihong  _et al._ [1994] Z. Weihong, J. Oitmaa, and C. J. Hamer, J. Phys. A: Math. Gen. 27, 5425 (1994).
  * Powalski  _et al._ [2013] M. Powalski, K. Coester, R. Moessner, and K. P. Schmidt, Phys. Rev. B 87, 054404 (2013).
  * Coester  _et al._ [2016] K. Coester, D. G. Joshi, M. Vojta, and K. P. Schmidt, Phys. Rev. B 94, 125109 (2016).
  * Röchner  _et al._ [2016] J. Röchner, L. Balents, and K. P. Schmidt, Phys. Rev. B 94, 201111 (2016).
  * Vidal  _et al._ [2009] J. Vidal, S. Dusuel, and K. P. Schmidt, Phys. Rev. B 79, 033109 (2009).
  * Dusuel  _et al._ [2010] S. Dusuel, M. Kamfor, K. P. Schmidt, R. Thomale, and J. Vidal, Phys. Rev. B 81, 064412 (2010).
  * Dusuel  _et al._ [2011] S. Dusuel, M. Kamfor, R. Orús, K. P. Schmidt, and J. Vidal, Phys. Rev. Lett. 106, 107203 (2011).
  * Schmidt [2013] K. P. Schmidt, Phys. Rev. B 88, 035118 (2013).
  * Knetter and Uhrig [2000] C. Knetter and G. Uhrig, Eur. Phys. J. B 13, 209 (2000).
  * Knetter  _et al._ [2003] C. Knetter, K. P. Schmidt, and G. S. Uhrig, J. Phys. A: Math. Gen. 36, 7889 (2003).
  * Takahashi [1977] M. Takahashi, J. Phys. C: Solid State Phys. 10, 1289 (1977).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 16, 407 (1961).
  * Pfeuty [1970] P. Pfeuty, Ann. Phys. 57, 79 (1970).
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Z. Physik 47, 631 (1928).
  * Wegner [1994] F. Wegner, Ann. Phys. 506, 77 (1994).
  * Powalski  _et al._ [2015] M. Powalski, G. S. Uhrig, and K. P. Schmidt, Phys. Rev. Lett. 115, 207202 (2015).
  * Powalski  _et al._ [2018] M. Powalski, K. P. Schmidt, and G. S. Uhrig, SciPost Phys. 4, 001 (2018).
  * Kato [1949] T. Kato, Prog. Theor. Phys. 4, 514 (1949).
  * Guttmann [1989] A. C. Guttmann,  _Phase Transitions and Critical Phenomena_ , edited by C. Domb and J. Lebowitz, Vol. 13 (Academic Press, New York, 1989).
  * Pfeuty and Elliott [1971] P. Pfeuty and R. J. Elliott, Jour. Phys. C: Solid State Phys. 4, 2370 (1971).
  * Kos  _et al._ [2016] F. Kos, D. Poland, D. Simmons-Duffin, and A. Vichi, J. High Energ. Phys. 2016, 36.
  * Trebst  _et al._ [2007] S. Trebst, P. Werner, M. Troyer, K. Shtengel, and C. Nayak, Phys. Rev. Lett. 98, 070602 (2007).
  * Tupitsyn  _et al._ [2010] I. S. Tupitsyn, A. Kitaev, N. V. Prokof’ev, and P. C. E. Stamp, Phys. Rev. B 82, 085114 (2010).
  * Blöte and Deng [2002] H. W. J. Blöte and Y. Deng, Phys. Rev. E 66, 066110 (2002).
  * Schuler  _et al._ [2016] M. Schuler, S. Whitsitt, L.-P. Henry, S. Sachdev, and A. M. Läuchli, Phys. Rev. Lett. 117, 210401 (2016).