# Two-dimensional lattice with an imaginary magnetic field

Tomoki Ozawa  Advanced Institute for Materials Research (WPI-AIMR), Tohoku
University, Sendai 980-8577, Japan  Tomoya Hayata  Departments of Physics,
Keio University, 4-1-1 Hiyoshi, Kanagawa 223-8521, Japan  RIKEN iTHEMS, RIKEN,
Wako 351-0198, Japan

###### Abstract

We introduce a two-dimensional non-Hermitian lattice model with an imaginary
magnetic field and elucidate various unique features which are absent in
Hermitian lattice models with real magnetic fields. To describe the imaginary
magnetic field, we consider both the Landau gauge and the symmetric gauge,
which are related by a generalized gauge transformation, changing not only the
phase but also the amplitude of the wave function. We discuss the complex
energy spectrum and the non-Hermitian Aharonov-Bohm effect as examples of
properties which are due to the imaginary magnetic field independent of the
generalized gauge transformation. We show that the energy spectrum does not
converge as the lattice size is made larger, which comes from the intrinsic
nonperiodicity of the model. However, we have found that the energy spectrum
does converge if one fixes the length of one side and makes the other side
longer; this asymptotic behavior can be understood in the framework of the
non-Bloch band theory. We also find an analog of the Aharonov-Bohm effect; the
net change of the norm of the wave function upon adiabatically forming a
closed path is determined by the imaginary magnetic flux enclosed by the path,
which provides an experimentally observable feature of the imaginary magnetic
field.

††preprint: RIKEN-iTHEMS-Report-23

##  I Introduction

Physics of a charged particle in an external magnetic field has been of
fundamental importance in condensed matter physics. In two dimensions, a
charged particle in a magnetic field forms the equally spaced energy spectrum
called the Landau level, which is directly responsible for phenomena such as
the Landau diamagnetism, de Haas-van Alphen effect [1], and the integer and
fractional quantum Hall effects [2, 3]. Moreover, a charged particle on a two-
dimensional lattice under a magnetic field is described by the Harper-
Hofstadter model [4, 5], which is the paradigmatic model of the Chern
insulator [6].

Recently, there has been an increasing interest in non-Hermitian physics [7,
8]. Also in non-Hermitian quantum mechanics, the effect of vector potentials
has played significant roles. For example, the seminal Hatano-Nelson model is
the one-dimensional lattice model under an imaginary vector potential [9, 10],
and has been of fundamental importance showing the non-Hermitian skin effect
and nontrivial point gap topology [11, 12, 13, 14, 15]. The non-Hermitian skin
effect has been experimentally realized recently in a variety of systems
including mechanical metamaterials [16, 17], electrical circuits [18, 19],
photonics [20, 21, 22], and ultracold atomic gases [23]. The imaginary vector
potential has also been crucial in understanding the Landau-Zener transition
[24, 25].

With the recent experimental development of non-Hermitian quantum mechanics,
one can now realize a variety of non-Hermitian models under control, and there
is an increasing interest in experimentally realizing two or higher
dimensional non-Hermitian models [26, 27, 28, 29]. Despite this rapid progress
in non-Hermitian physics and the important role of magnetic fields played in
condensed matter physics, there has been little study on properties of the
imaginary magnetic fields in two dimensions, namely, properties of non-
Hermitian lattice systems where magnetic fields are imaginary, analogous to
the imaginary vector potential in the Hatano-Nelson model.

In this paper, we explore basic properties of two-dimensional lattices with an
imaginary magnetic field. We first elucidate the meaning of gauge invariance
in non-Hermitian settings, to distinguish between properties which are gauge
independent or gauge invariant, i.e., intrinsically due to the imaginary
magnetic field, and which are gauge dependent, i.e., dependent on specific
realizations and setups. We find that certain spectral properties are gauge
invariant, and discuss that the asymptotic energy spectrum as one makes the
length of the system larger, fixing the other length, can be nicely understood
within the framework of the non-Bloch band theory. Especially, even though the
system is non-Hermitian, the asymptotic spectrum under open boundary
conditions can be related to the spectrum under the periodic boundary
condition satisfying certain conditions. We also find an analog of the
Aharonov-Bohm effect for imaginary magnetic fields. Upon making a wave packet
move to form a closed trajectory in real space, the overall change of the norm
of the wave function is related to the imaginary magnetic flux enclosed by the
trajectory. Our paper lays a foundation to understand gauge invariant
properties in the setup of imaginary magnetic fields, generalizing the concept
of magnetic fields to two-dimensional non-Hermitian settings.

The structure of the paper is the following. In Sec. II, we introduce the
model we study, which is the two-dimensional non-Hermitian lattice model with
an imaginary magnetic field. In Sec. III, we discuss generalized gauge
transformations which are relevant to discussing properties intrisically due
to the imaginary magnetic field. In Sec. IV, we discuss various properties of
the complex energy spectrum of the model. We first discuss open boundary
conditions in both directions and then the effect of taking one direction to
be periodic, which we call the cylindrical configuration. In Sec. V, we
discuss the asymptotic limit of the energy spectrum when we fix the length of
one side and make the length of the other side infinitely long. We analyze the
asymptotic spectrum using the non-Bloch band theory. In Sec. VI, we discuss
the analog of the Aharonov-Bohm effect under the imaginary magnetic field. We
finally give the conclusion and future prospects in Sec. VII.

##  II Model

We consider a two-dimensional square lattice with an imaginary magnetic field.
We label the lattice sites by coordinates $(x,y)$, where $x$ and $y$ are both
integers. We let $\psi_{x,y}$ denote the amplitude of the wave function at
site $(x,y)$. The Schrödinger equation governing the dynamics of the system is

$\displaystyle i\frac{d\psi_{x,y}}{dt}=$ $\displaystyle
J\left(e^{i\theta_{X}(x-1,y)}\psi_{x-1,y}+e^{-i\theta_{X}(x,y)}\psi_{x+1,y}\right.$
$\displaystyle\left.+e^{i\theta_{Y}(x,y-1)}\psi_{x,y-1}+e^{-i\theta_{Y}(x,y)}\psi_{x,y+1}\right),$
(1)

where $t$, and $J$ are the time, and hopping parameter, respectively. In this
paper, we consider two gauge choices: the Landau gauge and the symmetric
gauge. The Landau gauge is defined by $(\theta_{X},\theta_{Y})=(0,Bx)$, and
the symmetric gauge is defined by $(\theta_{X},\theta_{Y})=(-By/2,Bx/2)$. When
$B$ is real, these gauges correspond to the ordinary Landau and symmetric
gauges with a real magnetic field. In this paper, however, we take $B$ to
represent a purely imaginary magnetic field, $B=i\mathcal{B}$, with
$\mathcal{B}$ being a real number. We note that, when $B$ is imaginary, the
factors $e^{i\theta_{X}(x,y)}$ and $e^{i\theta_{Y}(x,y)}$ can have a modulus
different from one, implying non-Hermiticity.

##  III Gauge transformation

Landau and symmetric gauges are equivalent for a real magnetic field because
they are related by a gauge transformation. We first review the gauge
transformation in Hermitian setups, and then extend the concept to non-
Hermitian settings. The gauge transformation is to consider a state which is
related to the original state by a position-dependent phase factor
$\psi^{\prime}_{x,y}=e^{i\chi(x,y)}\psi_{x,y}$, where $\chi(x,y)$ is a real
function. This transformation amounts to applying a local unitary
transformation to wave function. If we take the Hamiltonian $H$ in the
Hermitian Landau gauge, by choosing the gauge transformation
$e^{i\chi(x,y)}=e^{iBxy/2}$, the Hamiltonian $H^{\prime}$ transformed under
the gauge transformation is in the Hermitian symmetric gauge. Since the gauge
transformation is a unitary transformation, the energy spectrum is invariant
under the gauge transformation.

Now we extend the concept of gauge transformations to non-Hermitian setups of
imaginary magnetic fields. An important feature of imaginary magnetic fields
is that the Landau gauge and the symmetric gauge cannot be connected via an
ordinary gauge transformation $e^{i\chi(x,y)}$ determined by a real function
$\chi(x,y)$. Instead, it is more appropriate to consider a generalized gauge
transformation, $\psi^{\prime}_{x,y}=f(x,y)\psi_{x,y}$ with $f(x,y)$ being a
nonzero complex function, which does not just multiply a phase factor but also
allows scale change for the wavefunction. The Hamiltonian changes under this
generalized gauge transformation, not by a unitary transformation but by a
local (diagonal) similarity transformation. The Landau and symmetric gauges
are related via the generalized gauge transformation
$f(x,y)=e^{iBxy/2}=e^{-\mathcal{B}xy/2}$. Since this generalized gauge
transformation is a similarity transformation, the energy spectrum is
invariant. Furthermore, upon the generalized gauge transformation, the product
of hopping amplitudes as one goes around a plaquette of the square lattice
does not change, implying that the imaginary magnetic field is also invariant.

There are various realizations of non-Hermitian Hamiltonians, and what is
observable depends on the individual system that one works on. Upon studying
properties of imaginary magnetic fields, one should thus make a clear
distinction between what are universal properties of imaginary magnetic fields
and what are gauge- and system-specific features which depend on particular
realizations. We consider properties intrinsic to imaginary magnetic fields to
be those invariant under the generalized gauge transformation.

##  IV Energy spectrum

As is well-known in the study of non-Hermitian Hamiltonians, an energy
spectrum under periodic and open boundary conditions can take drastically
different values [9, 10, 30, 14]. We should therefore analyze the energy
spectrum together with the boundary conditions. We first note that, unlike the
case of real magnetic fields, lattice models with imaginary magnetic fields
cannot be made periodic in both $x$ and $y$ directions. For the Landau gauge,
we can make the lattice periodic in the $y$ direction, but not in the $x$
direction, and for the symmetric gauge we cannot make the Hamiltonian periodic
in either direction. In this paper, we call the Landau gauge with the periodic
boundary condition in the $y$ direction a cylindrical configuration. As we
see, the energy spectrum under open boundary conditions in both directions and
that under the cylindrical configuration are qualitatively different.

###  IV.1 Open boundary conditions

We first consider the open boundary conditions. As we have seen, the energy
spectrum under the Landau and symmetric gauges are the same because they are
related by the generalized gauge transformation. We also note that the energy
spectrum is invariant upon the change of the origin of the coordinate: the
spectrum is invariant under changing $x$ to $x+x_{0}$ and $y$ to $y+y_{0}$ in
the hopping factors $\theta_{X}$ and $\theta_{Y}$. This invariance can be
shown, for example, for the symmetric gauge, by noting that the shift $x\to
x+x_{0}$ can be realized by $f(x,y)=e^{-\mathcal{B}x_{0}y/2}$ and the shift
$y\to y+y_{0}$ by $f(x,y)=e^{-\mathcal{B}xy_{0}/2}$. Since we do not want the
imaginary magnetic fields and their properties to depend on the origin of the
coordinates, these transformation properties are desirable.

Figure 1: Energy spectrum for different system sizes and magnetic fields under
open boundary conditions for (a)-(d) $B=0.001i$ and (e)-(h) $B=0.01i$. The
system sizes are (a), (e) $N_{x}=N_{y}=20$; (b), (f) $N_{x}=N_{y}=40$; (c),
(g) $N_{x}=N_{y}=60$; and (d), (h) $N_{x}=N_{y}=80$, respectively. All axes
are in units of $J$.

In Fig. 1, we plot the energy spectrum in the complex plane for a lattice of
size $N_{x}\times N_{y}$ with $N_{x}=N_{y}=20,40,60$, and 80 for the imaginary
magnetic fields of $B=0.001i$, and $0.01i$. We observe that the the spread of
the energy spectrum along the real axis is from $-4J$ to $4J$, whereas the
spread of the energy along imaginary axis varies with the strength of the
imaginary magnetic fields as well as the size of the system. This spread of
the energy spectrum along real and imaginary axes can be explained assuming
$B=\mathcal{B}i$ is small, as we shall show.

To understand the spectral range of $H(B)$, we expand the Hamiltonian $H(B)$
in the Taylor series around $B=0$ and keep up to the first order in $B$
assuming $B$ is small. We write this expansion as $H(B)=H(0)+B\delta H$. When
$B$ is purely imaginary, the first-order term $B\delta H=i\mathcal{B}\delta H$
is an antisymmetric matrix. Since the zeroth order term $H(0)$ is a symmetric
matrix, $H(B)=H(0)+i\mathcal{B}\delta H$ is a decomposition of $H(B)$ into its
symmetric and antisymmetric components. Since $H(0)$ is a symmetric matrix,
its eigenvalues are all purely real, whereas the eigenvalues of
$i\mathcal{B}\delta H$, which is an anti-symmetric matrix, are all purely
imaginary. In the Appendix, we show that, for any diagonalizable real matrix
$M$, by writing $M=S+A$, where $S$ and $A$ are symmetric and antisymmetric
real matrices, the real part of the eigenvalues of $M$ are contained within
the minimum and maximum values of the eigenvalues of $S$, and the imaginary
part of the eigenvalues of $M$ are contained within the minimum and maximum
imaginary eigenvalues of $A$. The spectral range of
$H(B)=H(0)+i\mathcal{B}\delta H$ is thus contained, in the real-axis
direction, by the minimum and maximum eigenvalues of $H(0)$ and is contained
in the imaginary-axis direction by the minimum and maximum imaginary
eigenvalues of $i\mathcal{B}\delta H$. Now, $H(0)$ is the Hamiltonian of a
two-dimensional square lattice without any magnetic field. The eigenvalues in
the thermodynamic limit are given by the simple sinusoidal band structure
$2J\left(\cos(k_{x})+\cos(k_{y})\right)$, where $(k_{x},k_{y})$ is the two-
dimensional quasimomentum. Therefore, their maximum and minimum are given by
$\pm 4J$. The antisymmetric part, $i\mathcal{B}\delta H$, in the Landau gauge,
is the Hamiltonian of decoupled $N_{x}$ one-dimensional chains with
nonreciprocal hoppings; for a chain at position $x$, the hopping along the
chain is $\pm\mathcal{B}xJ$, where $x$ runs from some initial value $x_{0}$ to
$x_{0}+N_{x}-1$. (Note that, in $i\mathcal{B}\delta H$, there is no hopping
between chains.) For a chain at position $x$, its eigenvalues, as a function
of the quasimomentum $k_{y}$ along $y$ direction, are
$2iJ\mathcal{B}x\cos(k_{y})$. The range of eigenvalues of $i\mathcal{B}\delta
H$ thus depends on how $x_{0}$ is chosen. To obtain the tightest bound, we can
take $x_{0}$ to be $-(N_{x}-1)/2$, so $x$ takes values from $-(N_{x}-1)/2$ to
$(N_{x}-1)/2$. The minimum and maximum imaginary eigenvalues of
$i\mathcal{B}\delta H$ are then $\pm J(N_{x}-1)\mathcal{B}i$. The spread of
the energy spectrum of $H(B)$ obtained numerically in Fig. 1 along the
imaginary direction agrees with the bound $\pm J(N_{x}-1)\mathcal{B}i$; we
note that the actual spectral spread along the imaginary direction is roughly
half compared to the bound $\pm J(N_{x}-1)\mathcal{B}i$. The spectrum of
$i\mathcal{B}\delta H$ provides only the upper and lower bounds, and does not
necessarily give limits that can be saturated.

We note that, while the spread along the real axis obtained from the
eigenvalues of $H(0)$ does not depend on the system size, the spread along the
imaginary axis obtained from the eigenvalues of $i\mathcal{B}\delta H$ depends
linearly on $N_{x}=N_{y}=N$. Numerically calculated eigenvalues in Fig. 1 also
shows that, by fixing the magnetic field and increasing the system size, the
energy spectrum increases its spreading along the imaginary direction and does
not converge. This nonconvergence of the energy spectrum as $N\to\infty$ is in
stark contrast to ordinary two-dimensional Hermitian lattices in which
increasing the system size makes the energy spectrum converge to a continuous
band structure, except for possible edge-localized states. The origin of the
non-convergence of our energy spectrum is because, even though the imaginary
magnetic field is fixed and constant over the entire lattice, the hopping
strength such as $e^{-\mathcal{B}x}$ keeps increasing in the $x$ direction.

The spectral bounds we have obtained assume that the imaginary magnetic field
$B=i\mathcal{B}$ is small. As the magnetic field becomes larger, the spread of
the energy spectrum of $H(B)$ can go beyond the bounds obtained by assuming
small $B$. Even in such a case, the spectral bounds can be obtained by
decomposing $H(B)$ into the symmetric and asymmetric parts. This predicts that
the bounds on both the real and imaginary parts grow exponentially as the
imaginary magnetic field or the system size increases. However, the numerical
diagonalization of $H(B)$ becomes more and more difficult as $B$ becomes
larger because of the existence of matrix elements which significantly differ
in their magnitudes. In this paper, we focus our analysis on the case of small
$B$, where we can make definite statements, both numerically and analytically,
about the spectral structure and bounds.

Although the energy spectrum does not converge when both sides of the lattice
are taken to increase, we will see later that the energy spectrum does
converge when we fix the length of one side and take the other direction to
increase. We discuss the asymptotic spectrum when one side is made longer in
Sec. V.

###  IV.2 Cylindrical configuration

We next consider the case of the periodic boundary condition along the $y$
direction under the Landau gauge, namely, the cylindrical configuration. The
energy spectrum can be obtained by performing the Fourier transformation in
the $y$ direction and diagonalizing the Hamiltonian for each momentum
separately. Writing $\psi_{x,y}=\psi_{x}e^{ik}$, the equation to solve is

$\displaystyle
E\psi_{x}=J\left\\{\psi_{x-1}+\psi_{x+1}+\left(e^{-x\mathcal{B}-ik}+e^{x\mathcal{B}+ik}\right)\psi_{x}\right\\}.$
(2)

We note that this is an analog of the Harper equation for the imaginary
magnetic field [4]. In Fig. 2, we plot the energy spectrum in the cylindrical
configuration with two different values of the imaginary magnetic field and
two different ways to choose the origin of $x=0$. We find an unexpected
feature that the energy spectrum depends on the origin of the coordinates.
Under open boundary conditions, we saw that shifting of $x\to x+x_{0}$ is
achieved by the generalized gauge transformation of
$f(x,y)=e^{-\mathcal{B}x_{0}y}$. However, this gauge transformation is not
periodic in the $y$ direction and thus is not compatible with the cylindrical
configuration. The energy spectrum therefore depends on how the origin of $x$
is taken. This is a unique feature of the imaginary magnetic field; the energy
spectrum of the real magnetic field is independent of how the origin of $x$ is
taken, even when the cylindrical configuration is taken along the $y$
direction under the Landau gauge.

Figure 2: Energy spectrum under the periodic boundary condition along the $y$
direction (cylindrical configuration) with $N_{x}=40$ for (a), (b) $B=0.001i$
and (c), (d) $B=0.01i$. The $x$ coordinates are labeled as $x=0,1,2,\cdots,40$
for (a), (c) and $x=-20,-19,\cdots,19$ for (b), (d). All the axes are in units
of $J$.

##  V Asymptotic spectrum

Even though the spectrum does not converge keeping $N_{x}=N_{y}$, we find that
the spectrum does converge as one fixes the size of one side and makes the
other side become longer. In this section, we discuss the asymptotic behavior
of the energy spectrum as one side is made longer. In Fig. 3, we plot the
energy spectrum under open boundary conditions when $B=0.01i$, fixing
$N_{x}=40$ and choosing $N_{y}=50$, 100, 150. One sees that the overall shape
tends to stabilize as $N_{y}$ becomes large. We can understand this asymptotic
energy spectrum in the limit of large $N_{y}$ by means of the non-Bloch band
theory [31, 32]. The non-Bloch band theory is a formalism to obtain the
continuous energy spectrum of non-Hermitian systems under open boundary
conditions. To understand the asymptotic behavior of fixing $N_{x}$ and making
$N_{y}\to\infty$, we now regard the index $x$ to be an internal index of a
one-dimensional system elongated along the $y$ direction.

Figure 3: (a)-(c) Asymptotic energy spectrum fixing $N_{x}=40$ under open
boundary conditions for $B=0.01i$. (a) $N_{y}=50$. (b) $N_{y}=100$. (c)
$N_{y}=150$. (d) Asymptotic energy spectrum predicted from the non-Bloch band
theory. All axes are in units of $J$.

###  V.1 Non-Bloch band theory

To apply the non-Bloch band theory, we perform the Fourier transformation
along the $y$ direction, as done in Eq. (2) above. In the non-Bloch band
theory, we replace $e^{ik}$ by a general complex number $\beta$, and solve the
above eigenvalue equation for a given value of $E$. Writing the above
eigenvalue equation as $E\vec{\psi}_{X}=H_{X}(\beta)\vec{\psi}_{X}$, where
$\vec{\psi}_{X}$ is a vector whose element is $\psi_{x}$, solutions to the
eigenvalue equation for a given value of $E$ are given by the solutions of
$\det[H_{X}(\beta)-E]=0$. This equation is an algebraic equation for $\beta$
with degree $2N_{x}$, and thus we generally have $2N_{x}$ solutions of
$\beta$. Writing $2N_{x}$ solutions of $\beta$ in the ascending order of their
magnitudes and labeling them as $\beta_{1}$, $\beta_{2}$, $\cdots$, the
eigenvalue $E$ belongs to a continuum of energy band if and only if
$|\beta_{N_{x}}|=|\beta_{N_{x}+1}|$ [31]. The corresponding values of
$\beta_{N_{x}}$ and $\beta_{N_{x}+1}$ form the generalized Brillouin zone in
the complex plane.

As we show below, the generalized Brillouin zone coincides with the ordinary
Brillouin, namely, $\beta=e^{ik}$ for real $k$, when the $x$ coordinate is
labeled so $x=0$ is in the center of the system. We want to look for the
generalized Brillouin zone and the continuum of energy bands by solving
$\det[H_{X}(\beta)-E]=0$. Remembering that the spectrum under open boundary
conditions is independent of the choice of the origin of $x=0$, we can choose
the coordinates most convenient for our purpose. It turns out to be
particularly useful to take a coordinate system where $x=0$ is in the center,
namely, if $N_{x}$ is an odd number, we write $N_{x}=2p+1$ with a positive
integer $p$ and take the $x$ coordinate to be
$x=-p,-p+1,\cdots,-1,0,1,\cdots,p-1,p$. If $N_{x}$ is an even number, we write
$N_{x}=2p$ with a positive integer $p$ and take
$x=-p+\frac{1}{2},\cdots,-\frac{3}{2},-\frac{1}{2},\frac{1}{2},\frac{3}{2},\cdots,p-\frac{1}{2}$.

Let us first consider the case where $N_{x}$ is an odd number: $N_{x}=2p+1$.
We take the $x$ coordinate label as above so $x=0$ is in the center. Then,

$\displaystyle
H_{X}(\beta)=\begin{pmatrix}e^{-p\mathcal{B}}\beta+e^{p\mathcal{B}}\frac{1}{\beta}&1\\\
1&e^{(-p+1)\mathcal{B}}\beta+e^{(p-1)\mathcal{B}}\frac{1}{\beta}&\ddots\\\
&\ddots&\ddots\\\ &&&e^{-\mathcal{B}}\beta+e^{\mathcal{B}\frac{1}{\beta}}&1\\\
&&&1&\beta+\frac{1}{\beta}&1\\\
&&&&1&e^{\mathcal{B}}\beta+e^{-\mathcal{B}}\frac{1}{\beta}&\ddots\\\
&&&&&\ddots&\ddots&1\\\
&&&&&&1&e^{p\mathcal{B}}\beta+e^{-p\mathcal{B}}\frac{1}{\beta}\end{pmatrix},$
(3)

in units of $J$. One thus sees that if $\beta$ is a solution of
$\det[H_{X}(\beta)-E]=0$, so is $1/\beta$, namely, $\det[H_{X}(1/\beta)-E]=0$.
The same holds for even $N_{x}$.

This implies that, upon writing $2N_{x}$ solutions of $\det[H_{X}(\beta)-E]=0$
in ascending order, we should have $\beta_{N_{x}+1}=1/\beta_{N_{x}}$. From
this condition, together with the condition that
$|\beta_{N_{x}}|=|\beta_{N_{x}+1}|$ for $E$ to be in the continuum of energy
bands [31], we arrive at the condition $|\beta_{N_{x}}|=|\beta_{N_{x}+1}|=1$.
This implies that $\beta$ belonging to the generalized Brillouin zone can be
written as $\beta=e^{ik}$ with a real number $k$, namely, the generalized
Brillouin zone coincides with the ordinary Brillouin zone with quasimomentum
$0\leq k<2\pi$.

Conversely, for a given value of $\beta=e^{ik}$ with a real $k$, all $N_{x}$
eigenvalues of $H_{X}(e^{ik})$ belong to the continuum of energy bands. This
is because any eigenvalue $E$ of $H_{X}(e^{ik})$ satisfies
$\det[H_{X}(e^{ik})-E]=0$, and when writing $2N_{x}$ solutions of $\beta$ in
an ascending order, $\beta=e^{\pm ik}$ must appear at $N_{x}$th and
$N_{x}+1$th positions. We have thus shown that the continuum of energy bands
obtained upon fixing $N_{x}$ and making $N_{y}$ large is nothing but the
energy spectrum in the cylindrical configuration where coordinates in the $x$
direction are taken so $x=0$ is placed at the center.

This implies that the solutions of Eq. (2) for real $k$, which are nothing but
the energy spectrum of the cylindrical configuration, are the asymptotic
spectrum when fixing $N_{x}$ and making $N_{y}$ large under open boundary
conditions. The fact that the generalized Brillouin zone coincides with the
ordinary Brillouin zone implies that there is no non-Hermitian skin effect.
This absence of the non-Hermitian skin effect is related to the
$\mathcal{PT}$-symmetry present in the system [33].

In Fig. 3(d), we show the continuum bands obtained from the energy spectrum of
a cylindrical configuration, taking $x=0$ to be at the center. We see that the
spectra in Figs. 3(a)-3(c) indeed approaches that of Fig. 3(d). With different
values of $\mathcal{B}$ and $N_{x}$, we find that there is a general structure
of a continuous spectrum along the real axis and several oval structures
spread along the imaginary direction, but the exact number of ovals and the
spread along the imaginary direction depend on specific values of the
parameters. We stress that the energy spectrum under open boundary conditions
does not depend on how the coordinates are chosen. Nevertheless, the
asymptotic spectrum coincides with the energy spectrum in the cylindrical
configuration where the coordinates are chosen in a symmetric manner.

##  VI Aharonov-Bohm effect for imaginary magnetic fields

We now discuss an effect analogous to the Aharonov-Bohm effect [34] for
imaginary magnetic fields. The non-Hermitian Aharonov-Bohm effect in a
parameter space due to the complex Berry phase has been experimentally
observed for synthetic mechanical metamaterials [35, 36], but, to our
knowledge, it has never been observed in real space. We first formulate the
non-Hermitian Aharonov-Bohm effect, which is the amplification and/or decay of
the magnitude of the wave function due to the imaginary magnetic flux enclosed
by a closed trajectory. Following the standard argument of the Aharonov-Bohm
effect in Hermitian systems [37], the wave function $\vec{\psi}$ of a wave
packet, which we assume to be well localized in real space, is transported
over a path $\mathcal{C}$ in the presence of a vector potential. Then, the
wave function $\vec{\psi}$ is related to the wave function $\vec{\psi_{0}}$
transported over the same path $\mathcal{C}$ but in the absence of a vector
potential by

$\displaystyle\vec{\psi}=\exp\left(i\int_{\mathbf{r}_{i},\mathcal{C}}^{\mathbf{r}_{f}}\mathbf{A}(\mathbf{r}^{\prime})\cdot
d\mathbf{r}^{\prime}\right)\vec{\psi_{0}},$ (4)

where the line integral is taken along the particular path $\mathcal{C}$,
connecting its initial position $\mathbf{r}_{i}$ and the final position
$\mathbf{r}_{f}$ of the wave packet, has been transported. This formula is
also valid in non-Hermitian systems, where $\mathbf{A}(\mathbf{r})$ becomes a
complex number. Taking a closed path $\mathcal{C}$, the wave function after
the transportation is related to the wavefunction in the absence of a magnetic
field by

$\displaystyle\vec{\psi}=\exp\left(i\oint_{\mathcal{C}}\mathbf{A}(\mathbf{r}^{\prime})\cdot
d\mathbf{r}^{\prime}\right)\vec{\psi_{0}}=e^{i\Phi}\vec{\psi_{0}},$ (5)

where $\Phi$ is the total flux through a surface that the path $\mathcal{C}$
encloses. In our case, since the flux is pure-imaginary, the prefactor
$e^{i\Phi}$ is a real number, giving rise to the amplifiaction and/or decay of
the norm of the wave function. Since the Hamiltonian we are considering is
Hermitian in the absence of the imaginary vector potential, the norm of the
wave function $\vec{\psi_{0}}$ remains one during the transportation.
Therefore, the effect of the Aharonov-Bohm factor $e^{i\Phi}$ can be obtained
directly by looking at the norm of the wave function in the final state.

To numerically demonstrate this non-Hermitian Aharonov-Bohm effect in real
space, we consider the setup where we start from a localized wave packet
around the center of the lattice, and then add external forces to make the
wave packet move. As the trajectory of the wave packet forms a closed path,
the change of the magnitude of the wave function is precisely related to the
imaginary magnetic flux enclosed by the path. We now numerically demonstrate
the effect.

Figure 4: Simulation of the Aharonov-Bohm effect under imaginary magnetic
fields for the lattice size of $50\times 50$ with $B=0.001i$ under open
boundary conditions. (a1) The counter-clockwise trajectory of the mean
position of the wave function, starting from the center of the lattice, which
we take to be $(x_{0},y_{0})$ = (0,0). (a2) $|\vec{\psi}|$ as a function of
time (in units of $1/J$) for the counter-clockwise trajectory. The horizontal
dotted line is the theoretical value of the Aharonov-Bohm factor
$e^{iBA_{\mathrm{Area}}}=e^{-\mathcal{B}A_{\mathrm{Area}}}\approx 0.965$. (a3)
$|\vec{\psi}|$ as a function of time for the clockwise trajectory. The
horizontal dotted line is at $e^{\mathcal{B}A_{\mathrm{Area}}}\approx 1.037$.
We also performed the same simulation with a different starting point located
at $(x_{0},y_{0})=(-5,-5)$, whose trajectory is plotted in (b1). The modulus
of the wave function $|\vec{\psi}|$ as a function of time for the counter-
clockwise trajectory and clockwise trajectory are plotted in (b2) and (b3),
respectively, for the wave packet starting from $(x_{0},y_{0})=(-5,-5)$. We
note that, although the $|\vec{\psi}|$ during the evolution is different for
(a2), (a3) and (b2), (b3), the final values, which are solely determined by
the enclosed imaginary magnetic flux, are common. In all plots, the solid
lines and curves are for the Landau gauge, whereas the dashed lines and curves
are for the symmetric gauge, respectively.

As an initial state, we choose a normalized Gaussian wave packet
$\psi_{x,y}\propto e^{-\\{(x-x_{0})^{2}+(y-y_{0})^{2})\\}/(2\sigma^{2})}$
centered around the point $(x_{0},y_{0})$ with the spread $\sigma=5$. We apply
a force changing sinusoidally in time, created by a potential
$V_{x,y}=E_{x}\sin(2\pi t/T)x+E_{y}\sin(2\pi t/T)y$, so the wave packet makes
a rectangular trajectory either in the counter-clockwise or clockwise
direction. For the counter-clockwise trajectory, we apply

$\displaystyle(E_{x},E_{y})=\begin{cases}(1,0)&\text{for }0\leq t\leq T\\\
(0,1)&\text{for }T\leq t\leq 2T\\\ (-1,0)&\text{for }2T\leq t\leq 3T\\\
(0,-1)&\text{for }3T\leq t\leq 4T.\\\ \end{cases}$ (6)

For the clockwise trajectory, we make the wave packet move in the opposite
direction by applying

$\displaystyle(E_{x},E_{y})=\begin{cases}(0,1)&\text{for }0\leq t\leq T\\\
(1,0)&\text{for }T\leq t\leq 2T\\\ (0,-1)&\text{for }2T\leq t\leq 3T\\\
(-1,0)&\text{for }3T\leq t\leq 4T.\\\ \end{cases}$ (7)

We use $T=5/J$ in the numerical simulation.

In the following numerical simulation, we choose a lattice of size
$N_{x}=N_{y}=50$ with an imaginary magnetic field $B=0.001i$ under open
boundary conditions. We use coordinates so $(x,y)=(0,0)$ is located at the
center of the lattice. Starting from a wave packet centered around
$(x_{0},y_{0})=(0,0)$, and evolving in time until $t=4T$, the center of the
wave packet forms rectangles as plotted in Fig. 4(a1) for the counter-
clockwise trajectory. In Fig. 4(a2), we plot the modulus of the wave function,
$|\vec{\psi}|\equiv\sqrt{\sum_{x,y}|\psi_{x,y}|^{2}}$, as a function of time
for both the Landau and symmetric gauges for the counter-clockwise trajectory.
The same quantity for the clockwise trajectory is also plotted in Fig. 4(a3).
We see that during the time evolution, $|\vec{\psi}|$ is generally different
between the two gauges, but it is the same after the closed trajectory is
formed. The Aharonov-Bohm factor corresponding to this trajectory is
$e^{i\Phi}=e^{iBA_{\mathrm{Area}}}=e^{-\mathcal{B}A_{\mathrm{Area}}}\approx
0.965$ for the counter-clockwise trajectory, where the area enclosed by the
trajectory is $A_{\mathrm{Area}}\approx 36$ for both gauges. For the clockwise
trajectory,
$e^{i\Phi}=e^{-iBA_{\mathrm{Area}}}=e^{\mathcal{B}A_{\mathrm{Area}}}\approx
1.037$. Those values are also plotted in Figs. 4(a2) and 4(a3). Although, both
for counter-clockwise and clockwise trajectories, the time evolution of
$|\vec{\psi}|$ depends on the choice of gauge, the final values of Figs. 4(a2)
and 4(a3) agree well with the Aharonov-Bohm factor and are independent of the
gauge choice.

The detailed time-dependence of $|\vec{\psi}|$ can also be understood from the
Aharonov-Bohm factor during the evolution, which is
$\exp\left(i\int_{\mathbf{r}_{i},\mathcal{C}}^{\mathbf{r}_{f}}\mathbf{A}(\mathbf{r}^{\prime})\cdot
d\mathbf{r}^{\prime}\right)$. For example, during $0\leq t\leq T$,
$|\vec{\psi}|$ remains to be 1 both for the Landau and symmetric gauges in
Fig. 4(a2). This is because, both for the Landau gauge $\mathbf{A}=(0,Bx)$ and
the symmetric gauge $\mathbf{A}=(-By/2,Bx/2)$, the line integral along the
horizontal direction is zero if the wave packet starts from
$(x_{0},y_{0})=(0,0)$.

This observation also implies that the time evolution of $|\vec{\psi}|$
depends on the initial position of the wavepacket, whereas the final values of
$|\vec{\psi}|$ are independent of the initial position and the gauge one
chooses. This has been numerically demonstrated in Figs. 4(b1)-4(b3), where
the initial wavepacket is centered around $(x_{0},y_{0})=(-5,-5)$. Comparing
Fig. 4(a2) with Fig. 4(b2) for counter-clockwise trajectories, and Fig. 4(a3)
with Fig. 4(b3) for clockwise trajectories, we confirm that the final values
of $|\vec{\psi}|$ are $e^{\pm i\Phi}$, whereas the values during the evolution
depend on the initial position of the wave packet and the gauge choice.

##  VII Conclusion

We have studied spectral and geometrical properties of two-dimensional
lattices under a uniform imaginary magnetic field. Our results unveil features
of imaginary magnetic fields which are intrinsically different from real
magnetic fields, such as the impossibility to take periodic boundary
conditions in both directions and non-convergence of the energy spectrum in
the limit when both sides are taken large. On the other hand, there also are
similarities to the real magnetic field, such as description in terms of the
Harper equation and the analog of the Aharonov-Bohm effect. Although we
focused on the cases of purely imaginary magnetic fields, general results
presented in the paper, such as the non-Bloch band theory when increasing the
length of one direction and the non-Hermitian Aharonov-Bohm effect, should be
valid also for more general complex magnetic fields including both real and
imaginary components.

Experimental realization of the imaginary magnetic field we discussed in the
paper requires constructing an extended lattice with desired non-reciprocal
hopping amplitudes, which can be achieved in the platforms of mechanical
metamaterials [16, 17] and electrical circuits [18, 19, 35] in a
straightforward manner. The model, in the Landau gauge, can also be realized
by aligning a collection of one-dimensional Hatano-Nelson models and coupling
them; such a collection of Hatano-Nelson models can be realized through
Floquet methods, e.g., through the method described in Ref. [38], in various
photonic systems such as coupled waveguides or coupled resonators. In these
classical setups, change of the norm of the wavefunction due to the non-
Hermitian Aharonov-Bohm effect is observable through the amplification and/or
damping of corresponding signals. There is also a recent theoretical proposal
of an Aharonov-Bohm-type setting where the complex vector potential emerges
after post selection [39].

Our results provide a starting point toward the research field of non-
Hermitian magnetic fields. We have focused on physics on lattices;
understanding unique properties under an imaginary magnetic field in
continuous two-dimensional systems is also an open field of study. Moreover,
understanding properties under more general gauge fields such as complex
electromagnetic fields and non-Abelian gauge fields (e.g., spin-orbit
coupling) is left for future study.

###### Acknowledgements.

The authors thank Shuichi Murakami for helpful discussions on the non-Bloch
band theory. This work is supported by JSPS KAKENHI Grant No. JP20H01845,
Grant No. JP21H01007, Grant No. JP21H01084, JST PRESTO Grant No. JPMJPR2353,
and JST CREST Grant No.JPMJCR19T1.

*

##  Appendix A Spectral bounds of matrices

We show that the real part of the eigenvalues of a diagonalizable matrix $M$
is bounded by the eigenvalues of its symmetric component, and the imaginary
part of the eigenvalues is bounded by the eigenvalues of its anti-symmetric
components. For the discussion in the main text, we only need to consider the
case where $M$ is a real diagonalizable matrix, but here we provide a more
general property which holds even when $M$ is a complex diagonalizable matrix.

Assume that $M$ is a complex diagonalizable matrix. We then write

$\displaystyle M=\frac{M+M^{\dagger}}{2}+\frac{M-M^{\dagger}}{2},$ (8)

where $M^{\dagger}$ is the Hermitian conjugate of $M$. Defining the first term
as $S=\frac{M+M^{\dagger}}{2}$ and the second term as
$\frac{M-M^{\dagger}}{2}$, we can see that $S$ is a Hermitian matrix and $A$
is an anti-Hermitian matrix, namely, $S^{\dagger}=S$ and $A^{\dagger}=-A$. We
note that eigenvalues of a Hermitian matrix and an anti-Hermitian matrix are
real and pure imaginary, respectively.

For any normalized vector $|v\rangle$, we see the following property:

$\displaystyle\langle v|S|v\rangle$ $\displaystyle=\langle
v|\frac{M+M^{\dagger}}{2}|v\rangle=\frac{1}{2}\left(\langle
v|M|v\rangle+\langle v|M|v\rangle^{*}\right)$
$\displaystyle=\mathrm{Re}\left(\langle v|M|v\rangle\right).$ (9)

Since $S$ is a Hermitian matrix, the minimum value of the left hand side is
given when $|v\rangle$ is an eigenvector of the minimum eigenvalue
$\lambda_{\mathrm{min}}$ of $S$:

$\displaystyle\mathrm{min}\langle v|S|v\rangle=\lambda_{\mathrm{min}}.$ (10)

Similarly, the maximum value of the left-hand side is equal to the maximum
eigenvalue $\lambda_{\mathrm{max}}$ of $S$. Thus we obtain the following
inequality:

$\displaystyle\lambda_{\mathrm{min}}\leq\mathrm{Re}\left(\langle
v|M|v\rangle\right)\leq\lambda_{\mathrm{max}}.$ (11)

Now, taking $|v\rangle$ to be an eigenvector of $M$ with an eigenvalue $E$,
$M|v\rangle=E|v\rangle$, we finally obtain the desired inequality:

$\displaystyle\lambda_{\mathrm{min}}\leq\mathrm{Re}(E)\leq\lambda_{\mathrm{max}}.$
(12)

We note that the equality is not necessarily achieved because eigenvectors of
$M$ are not orthogonal and there can be ways to minimize
$\mathrm{Re}\left(\langle v|M|v\rangle\right)$ by choosing $|v\rangle$ to not
be an eigenvector of $M$.

We now prove a similar relation for the spread of the eigenvalues of $M$ along
the imaginary axis. The property of the imaginary part of the eigenvalues that
we want to show can be mapped to the problem of the real part by writing

$\displaystyle-iM=-iS-iA;$ (13)

the second term $-iA$ is Hermitian and the first term $-iS$ is anti-Hermitian.
This implies that the spread of the eigenvalues of $-iM$ along the real axis
is bounded by the real eigenvalues of $-iA$. Therefore, the spread of the
eigenvalues of $M$ along the imaginary axis is bounded by the eigenvalues of
$A$.

We have thus proved that, writing an arbitrary diagonalizable complex matrix
$M$ in terms of its Hermitian and anti-Hermitian parts by $M=S+A$, the
spectral range of $M$ along the real part is bounded by the eigenvalues of $S$
and the imaginary part is bounded by the eigenvalues of $A$. When $M$ is a
real matrix, the same result holds when writing $M=S+A$, with $S$ being a
symmetric matrix and $A$ being an anti-symmetric matrix.

## References

  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid State Physics_ (Saunders College Publishing, New York, 1976).
  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, New method for high-accuracy determination of the fine-structure constant based on quantized Hall resistance, Phys. Rev. Lett. 45, 494 (1980).
  * Tsui  _et al._ [1982] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Two-dimensional magnetotransport in the extreme quantum limit, Phys. Rev. Lett. 48, 1559 (1982).
  * Harper [1955] P. G. Harper, Single band motion of conduction electrons in a uniform magnetic field, Proceedings of the Physical Society. Section A 68, 874 (1955).
  * Hofstadter [1976] D. R. Hofstadter, Energy levels and wave functions of bloch electrons in rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized Hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).
  * Ashida  _et al._ [2020] Y. Ashida, Z. Gong, and M. Ueda, Non-Hermitian physics, Advances in Physics 69, 249 (2020).
  * Bergholtz  _et al._ [2021] E. J. Bergholtz, J. C. Budich, and F. K. Kunst, Exceptional topology of non-Hermitian systems, Rev. Mod. Phys. 93, 015005 (2021).
  * Hatano and Nelson [1996] N. Hatano and D. R. Nelson, Localization transitions in non-Hermitian quantum mechanics, Phys. Rev. Lett. 77, 570 (1996).
  * Hatano and Nelson [1997] N. Hatano and D. R. Nelson, Vortex pinning and non-Hermitian quantum mechanics, Phys. Rev. B 56, 8651 (1997).
  * Gong  _et al._ [2018] Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, and M. Ueda, Topological phases of non-Hermitian systems, Phys. Rev. X 8, 031079 (2018).
  * Kawabata  _et al._ [2019] K. Kawabata, K. Shiozaki, M. Ueda, and M. Sato, Symmetry and topology in non-Hermitian physics, Phys. Rev. X 9, 041015 (2019).
  * Borgnia  _et al._ [2020] D. S. Borgnia, A. J. Kruchkov, and R.-J. Slager, Non-Hermitian boundary modes and topology, Phys. Rev. Lett. 124, 056802 (2020).
  * Okuma  _et al._ [2020] N. Okuma, K. Kawabata, K. Shiozaki, and M. Sato, Topological origin of non-Hermitian skin effects, Phys. Rev. Lett. 124, 086801 (2020).
  * Lin  _et al._ [2023] R. Lin, T. Tai, L. Li, and C. H. Lee, Topological non-Hermitian skin effect, Frontiers of Physics 18, 53605 (2023).
  * Brandenbourger  _et al._ [2019] M. Brandenbourger, X. Locsin, E. Lerner, and C. Coulais, Non-reciprocal robotic metamaterials, Nature communications 10, 4608 (2019).
  * Ghatak  _et al._ [2020] A. Ghatak, M. Brandenbourger, J. Van Wezel, and C. Coulais, Observation of non-Hermitian topology and its bulk–edge correspondence in an active mechanical metamaterial, Proceedings of the National Academy of Sciences 117, 29561 (2020).
  * Helbig  _et al._ [2020] T. Helbig, T. Hofmann, S. Imhof, M. Abdelghany, T. Kiessling, L. Molenkamp, C. Lee, A. Szameit, M. Greiter, and R. Thomale, Generalized bulk–boundary correspondence in non-Hermitian topolectrical circuits, Nature Physics 16, 747 (2020).
  * Hofmann  _et al._ [2020] T. Hofmann, T. Helbig, F. Schindler, N. Salgo, M. Brzezińska, M. Greiter, T. Kiessling, D. Wolf, A. Vollhardt, A. Kabaši, C. H. Lee, A. Bilušić, R. Thomale, and T. Neupert, Reciprocal skin effect and its realization in a topolectrical circuit, Phys. Rev. Res. 2, 023265 (2020).
  * Weidemann  _et al._ [2020] S. Weidemann, M. Kremer, T. Helbig, T. Hofmann, A. Stegmaier, M. Greiter, R. Thomale, and A. Szameit, Topological funneling of light, Science 368, 311 (2020).
  * Xiao  _et al._ [2020] L. Xiao, T. Deng, K. Wang, G. Zhu, Z. Wang, W. Yi, and P. Xue, Non-Hermitian bulk–boundary correspondence in quantum dynamics, Nature Physics 16, 761 (2020).
  * Xiao  _et al._ [2021] L. Xiao, T. Deng, K. Wang, Z. Wang, W. Yi, and P. Xue, Observation of non-bloch parity-time symmetry and exceptional points, Phys. Rev. Lett. 126, 230402 (2021).
  * Liang  _et al._ [2022] Q. Liang, D. Xie, Z. Dong, H. Li, H. Li, B. Gadway, W. Yi, and B. Yan, Dynamic signatures of non-Hermitian skin effect and topology in ultracold atoms, Phys. Rev. Lett. 129, 070401 (2022).
  * Fukui and Kawakami [1998] T. Fukui and N. Kawakami, Breakdown of the Mott insulator: Exact solution of an asymmetric Hubbard model, Phys. Rev. B 58, 16051 (1998).
  * Oka and Aoki [2010] T. Oka and H. Aoki, Dielectric breakdown in a Mott insulator: Many-body Schwinger-Landau-Zener mechanism studied with a generalized Bethe ansatz, Phys. Rev. B 81, 033103 (2010).
  * Yamauchi  _et al._ [2020] L. Yamauchi, T. Hayata, M. Uwamichi, T. Ozawa, and K. Kawaguchi, Chirality-driven edge flow and non-Hermitian topology in active nematic cells, arXiv preprint arXiv:2008.10852 (2020).
  * Zou  _et al._ [2021] D. Zou, T. Chen, W. He, J. Bao, C. H. Lee, H. Sun, and X. Zhang, Observation of hybrid higher-order skin-topological effect in non-Hermitian topolectrical circuits, Nature Communications 12, 7201 (2021).
  * Palacios  _et al._ [2021] L. S. Palacios, S. Tchoumakov, M. Guix, I. Pagonabarraga, S. Sánchez, and A. G. Grushin, Guided accumulation of active particles by topological design of a second-order skin effect, Nature Communications 12, 4691 (2021).
  * Shang  _et al._ [2022] C. Shang, S. Liu, R. Shao, P. Han, X. Zang, X. Zhang, K. N. Salama, W. Gao, C. H. Lee, R. Thomale, _et al._ , Experimental identification of the second-order non-Hermitian skin effect with physics-graph-informed machine learning, Advanced Science 9, 2202922 (2022).
  * Yao and Wang [2018] S. Yao and Z. Wang, Edge states and topological invariants of non-Hermitian systems, Phys. Rev. Lett. 121, 086803 (2018).
  * Yokomizo and Murakami [2019] K. Yokomizo and S. Murakami, Non-Bloch band theory of non-Hermitian systems, Phys. Rev. Lett. 123, 066404 (2019).
  * Yokomizo and Murakami [2023] K. Yokomizo and S. Murakami, Non-Bloch bands in two-dimensional non-Hermitian systems, Phys. Rev. B 107, 195112 (2023).
  * Yi and Yang [2020] Y. Yi and Z. Yang, Non-Hermitian skin modes induced by on-site dissipations and chiral tunneling effect, Phys. Rev. Lett. 125, 186802 (2020).
  * Aharonov and Bohm [1959] Y. Aharonov and D. Bohm, Significance of electromagnetic potentials in the quantum theory, Phys. Rev. 115, 485 (1959).
  * Anandwade  _et al._ [2023] R. Anandwade, Y. Singhal, S. N. M. Paladugu, E. Martello, M. Castle, S. Agrawal, E. Carlson, C. Battle-McDonald, T. Ozawa, H. M. Price, and B. Gadway, Synthetic mechanical lattices with synthetic interactions, Phys. Rev. A 108, 012221 (2023).
  * Singhal  _et al._ [2023] Y. Singhal, E. Martello, S. Agrawal, T. Ozawa, H. Price, and B. Gadway, Measuring the adiabatic non-Hermitian berry phase in feedback-coupled oscillators, Phys. Rev. Res. 5, L032026 (2023).
  * Griffiths and Schroeter [2018] D. J. Griffiths and D. F. Schroeter,  _Introduction to Quantum Mechanics_ , 3rd ed. (Cambridge University Press, Cambridge, 2018).
  * Faugno and Ozawa [2022] W. N. Faugno and T. Ozawa, Interaction-induced non-Hermitian topological phases from a dynamical gauge field, Phys. Rev. Lett. 129, 180401 (2022).
  * Paiva  _et al._ [2023] I. L. Paiva, Y. Aharonov, J. Tollaksen, and M. Waegell, Aharonov–Bohm effect with an effective complex-valued vector potential, New Journal of Physics 25, 053017 (2023).


#  Constraining red supergiant mass-loss prescriptions through supernova radio
properties

Takashi J. Moriya1,2
1National Astronomical Observatory of Japan, National Institutes of Natural
Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan
2School of Physics and Astronomy, Faculty of Science, Monash University,
Clayton, Victoria 3800, Australia E-mail<EMAIL_ADDRESS>(TJM)

(Accepted 2021 February 13. Received 2021 February 13; in original form 2021
January 13)

###### Abstract

Supernova properties in radio strongly depend on their circumstellar
environment and they are an important probe to investigate the mass loss of
supernova progenitors. Recently, core-collapse supernova observations in radio
have been assembled and the rise time and peak luminosity distribution of
core-collapse supernovae at 8.4 GHz has been estimated. In this paper, we
constrain the mass-loss prescriptions for red supergiants by using the rise
time and peak luminosity distribution of Type II supernovae in radio. We take
the de Jager and van Loon mass-loss rates for red supergiants, calculate the
rise time and peak luminosity distribution based on them, and compare the
results with the observed distribution. We found that the de Jager mass-loss
rate explains the widely spread radio rise time and peak luminosity
distribution of Type II supernovae well, while the van Loon mass-loss rate
predicts a relatively narrow range for the rise time and peak luminosity. We
conclude that the mass-loss prescriptions of red supergiants should have
strong dependence on the luminosity as in the de Jager mass-loss rate to
reproduce the widely spread distribution of the rise time and peak luminosity
in radio observed in Type II supernovae.

###### keywords:

supernovae: general – radio continuum: transients – circumstellar matter –
stars: massive – stars: mass-loss – supergiants

††pubyear: 2021††pagerange:  Constraining red supergiant mass-loss
prescriptions through supernova radio properties – Constraining red supergiant
mass-loss prescriptions through supernova radio properties

##  1 Introduction

Red supergiants (RSGs) are progenitors of Type II supernovae (SNe, Smartt 2015
for a review) that are the most commonly observed core-collapse SNe (Shivvers
et al., 2017). While the upper mass limit of RSGs is found to be around 25
$\mathrm{M}_{\odot}$, the maximum mass of the RSG SN progenitors estimated
based on the progenitors observationally identified so far is around
$17~{}\mathrm{M}_{\odot}$ (e.g., Kochanek et al., 2008; Smartt et al., 2009).
The possible reasons for the discrepancy in the maximum masses of RSGs and RSG
SN progenitors have been actively discussed (e.g., Horiuchi et al., 2011;
Walmswell & Eldridge, 2012; Adams et al., 2017; Sukhbold & Adams, 2020),
although the exact maximum mass of the RSG SN progenitors could be higher than
17 $\mathrm{M}_{\odot}$ and the maximum mass difference may not be
statistically significant yet (Davies & Beasor, 2020). In particular, the
uncertainties in RSG mass loss have been suggested to be a possible reason for
the discrepancy (e.g., Yoon & Cantiello, 2010; Georgy, 2012). Several
empirical mass-loss prescriptions have been proposed for RSGs (e.g., Mauron &
Josselin 2011 for a summary) and they have a strong influence on the predicted
properties of SN progenitors (e.g., Meynet et al., 2015).

Matter lost from SN progenitors exists around the progenitors and forms
circumstellar matter (CSM). Thus, SN explosions occur within the CSM that
imprints the mass-loss history of the SN progenitors. Especially, radio
emission from SNe, which is dominated by synchrotron emission from the forward
shock, is strongly affected by the CSM properties (e.g., Chevalier, 1982b,
1998) and it has long been used to probe the CSM properties of the SN
progenitors (e.g., Weiler et al., 2002). In other words, radio observations
from SNe II can be used to constrain the mass-loss properties of RSGs
(Chevalier et al., 2006). In addition to radio observations, early SN
observations in optical and ultraviolet are starting to unveil the existence
of the dense CSM at the immediate vicinity of RSG SN progenitors (e.g., Yaron
et al., 2017; Förster et al., 2018). The dense CSM could be linked to unknown
RSG mass-loss mechanism connected to core collapse (e.g., Quataert & Shiode,
2012; Fuller, 2017). In this paper, however, we focus on the distant CSM
formed through the canonical RSG mass loss probed by the late SN radio
observations.

Although SN radio light curves (LCs) are known to have some modulations (e.g.,
Weiler et al., 1992; Soderberg et al., 2006), they are generally characterized
by a single-peak LC. The time of the radio LC peak corresponds to the time
when the synchrotron self-absorption or the free-free absorption optical depth
becomes around unity (Chevalier, 1998). Therefore, the rise time
($t_{\mathrm{rise}}$) and peak luminosity ($L_{\mathrm{peak}}$) of the SN
radio LCs depend strongly on the CSM density. Bietenholz et al. (2020)
recently compiled a large sample of core-collapse SN observations at
$2-10~{}\mathrm{GHz}$ and estimated the rise time and peak luminosity
distributions of SNe at around 8.4 GHz. They found that the rise time and peak
luminosity distribution of SNe II (excluding Type IIn SNe) is characterized by
lognormal distributions. Excluding SN 1987A, which is the explosion of a blue
supergiant, they found $\log(t_{\mathrm{rise}}/\mathrm{days})=1.7\pm 1.0$ and
$\log(L_{\mathrm{peak}}/\mathrm{erg~{}s^{-1}~{}Hz^{-1}})=25.4\pm 1.2$ for SNe
II.

In this paper, we aim to constrain the RSG mass-loss prescriptions by using
the radio rise time and peak luminosity relation of SNe II estimated by
Bietenholz et al. (2020). The rest of the paper is organized as follows. We
first formulate the radio emission from SNe in Section 2. We then indroduce
the RSG mass-loss prescriptions we investigate in this paper in Section 3. We
compare the rise time and peak luminosity distributions expected from these
RSG mass-loss prescriptions with that constrained from the observations and
discuss implications in Section 4. We summarize this paper in Section 5.

Figure 1:  Radio luminosity evolution from SNe in the
$\dot{M}=10^{-6}~{}\mathrm{M}_{\odot}~{}\mathrm{yr^{-1}}$ CSM with
$\varepsilon_{B}=0.01$ and SSA. The radio LCs with $M_{\mathrm{ej}}=5,10,15$,
and 20 $\mathrm{M}_{\odot}$ and $E_{\mathrm{ej}}=0.1,0.5,1,1.5,$ and 2 B are
presented. The red line presents the locations of the LC peaks from all the
combinations of $M_{\mathrm{ej}}=5-20~{}\mathrm{M}_{\odot}$ and
$E_{\mathrm{ej}}=0.1-2~{}\mathrm{B}$.

##  2 Radio luminosity evolution

We first formulate the SN radio luminosity evolution. We assume that
synchrotron emission from electrons accelerated at the forward shock by the
amplified magnetic field is the dominant radio luminosity source in SNe. The
synchrotron emission from the forward shock can be formulated as

$\nu L_{\nu}\sim\pi
r_{\mathrm{sh}}^{2}v_{\mathrm{sh}}n_{\mathrm{rel}}\left(\frac{\gamma_{\nu}}{\gamma_{\mathrm{min}}}\right)^{1-p}\gamma_{\nu}m_{e}c^{2}\left[1+\frac{t_{\mathrm{sync},\nu}}{t}\right]^{-1},$
(1)

where $r_{\mathrm{sh}}$ is the forward shock radius, $v_{\mathrm{sh}}$ is the
forward shock velocity, $n_{\mathrm{rel}}$ is the number density of
accelerated relativistic electrons, $\gamma_{\nu}=(2\pi m_{e}c\nu/eB)^{0.5}$
is the Lorentz factor of the relativistic electrons with the characteristic
frequency of $\nu$, $\gamma_{\mathrm{min}}$ is the minimum Lorentz factor of
the accelerated electrons, $m_{e}$ is the electron mass,
$t_{\mathrm{sync,\nu}}=6\pi m_{e}c/\sigma_{T}\gamma_{\nu}B^{2}$ is the
synchrotron cooling timescale, $e$ is the electron charge, $B$ is the magnetic
field strength, and $\sigma_{T}$ is the Thomson scattering cross section. We
assume $\gamma_{\mathrm{min}}\sim 1$ and $n_{\mathrm{rel}}\propto\gamma^{-p}$
with $p=3$ (e.g., Chevalier & Fransson, 2006; Maeda, 2013). Because we use the
radio LC properties at 8.4 GHz compiled by Bietenholz et al. (2020), all the
radio luminosity evolution we calculate in this paper is at 8.4 GHz.

Major uncertain parameters in estimating the radio luminosity are the fraction
of the shock kinetic energy used to the electron acceleration
($\varepsilon_{e}$) and the fraction converted to magnetic energy
($\varepsilon_{B}$). We adopt $\varepsilon_{e}=0.1$ and
$\varepsilon_{B}=0.01-0.001$ in this work following the constraints on these
parameters based on multi-frequency SN observations (e.g., Björnsson &
Fransson, 2004; Maeda, 2012; Kamble et al., 2016).

The synchrotron emission (Eq. 1) is affected by absorption processes. We first
assume that the synchrotron self-absorption (SSA) is the dominant absorption
process (Chevalier, 1998). The SSA optical depth is
$\tau_{\mathrm{SSA}}=(\nu/\nu_{\mathrm{SSA}})^{-(p+4)/2}$, where
$\nu_{\mathrm{SSA}}\simeq 3\times
10^{5}(r_{\mathrm{sh}}\epsilon_{e}/\epsilon_{B})^{2/7}B^{9/7}$ Hz in the cgs
unit for $p=3$. The free-free absorption can also dominate the absorption
process especially when the CSM density is high (e.g., Chevalier, 1982b, 1998;
Chevalier et al., 2006; Weiler et al., 2002; Bietenholz et al., 2020). We
discuss the effect of the free-free absorption at the end of Section 4.

The hydrodynamic properties of the forward shock ($r_{\mathrm{sh}}$ and
$v_{\mathrm{sh}}$) are estimated based on the formulation in Chevalier
(1982a). The SN ejecta are assumed to have the two-component power-law density
structure ($\propto r^{-n}$ outside and $\propto r^{-\delta}$ inside) with
$n=12$ and $\delta=1$ (Matzner & McKee, 1999). The SN explosion energy is
varied from 0.1 B to 2 B, where $1~{}\mathrm{B}\equiv 10^{51}~{}\mathrm{erg}$,
and the SN ejecta mass is set from 5 $\mathrm{M}_{\odot}$ to 20
$\mathrm{M}_{\odot}$ (e.g., Pejcha & Prieto, 2015).

Our interest in this paper is in the canonical RSG mass loss before the onset
of the possible mass-loss enhancement shortly before RSG SN explosions. Thus,
we assume that the CSM structure is formed through steady mass loss and has
the following density structure

$\rho_{\mathrm{CSM}}(r)=\frac{\dot{M}}{4\pi v_{\mathrm{CSM}}}r^{-2},$ (2)

where $\dot{M}$ is the mass-loss rate and $v_{\mathrm{CSM}}$ is the wind
velocity. We set $v_{\mathrm{CSM}}=20~{}\mathrm{km~{}s^{-1}}$, which is
typically observed in RSGs (Goldman et al., 2017).

Combining the assumptions presented so far, we can estimate the radio
luminosity evolution for a given $\dot{M}$. Fig. 1 shows the radio luminosity
evolution for $\dot{M}=10^{-6}~{}\mathrm{M}_{\odot}~{}\mathrm{yr^{-1}}$ with
some combinations of $M_{\mathrm{ej}}$ and $E_{\mathrm{ej}}$ adopting only
SSA. For a given $\dot{M}$, the rise time and peak luminosity make a single
line when we plot the radio LCs with all the combinations of
$M_{\mathrm{ej}}=5-20~{}\mathrm{M}_{\odot}$ and
$E_{\mathrm{ej}}=0.1-2~{}\mathrm{B}$. An example is shown in Fig. 1 as the red
line in the case of $\dot{M}=10^{-6}~{}\mathrm{M}_{\odot}~{}\mathrm{yr^{-1}}$.

Figure 2:  Empirical mass-loss rates of RSGs as a function of their bolometric
luminosity. The effective temperature of 3500 K is assumed. We adopt the two
representative mass-loss prescriptions of van Loon et al. (2005) and de Jager
et al. (1988).

Figure 3:  Rise time and peak luminosity of SNe II expected from the de Jager
mass-loss rate (top) and the van Loon mass-loss rate (bottom).
$\varepsilon_{B}=0.01$ and SSA are adopted. The colored lines show those for
the RSG SN progenitor luminosity range ($\log L/\mathrm{L}_{\odot}=4.4-5.2$).
For a given luminosity, each line shows the rise time and peak luminosity
range for $M_{\mathrm{ej}}=5-20~{}\mathrm{M}_{\odot}$ and
$E_{\mathrm{ej}}=0.1-2~{}\mathrm{B}$ (cf. Fig. 1). The gray lines show RSGs
with $\log L/\mathrm{L}_{\odot}=5.3,5.4,$ and 5.5 which may not explode as SNe
II. The red squares show the rise time and peak luminosity of SNe II obtained
in Bietenholz et al. (2020). The blue ellipse shows the region where the 68%
of SNe II are estimated to exist based on all the observational information
available including upper limits (Bietenholz et al., 2020).

Figure 4:  Same as Fig. 3, but for $\varepsilon_{B}=0.001$.

##  3 RSG mass-loss prescriptions

Many empirical prescriptions for the RSG mass-loss rate have been proposed
(e.g., Mauron & Josselin, 2011, for a summary). They can be broadly classified
to two types; those with relatively high mass-loss rates as represented by van
Loon et al. (2005) and those with relatively low mass-loss rates as
represented by de Jager et al. (1988). We show these RSG mass-loss rates in
Fig. 2. The empirical mass-loss rates depend on the bolometric luminosity
($L$) and effective temperature of RSGs. Because we investigate the RSG mass
loss through SNe, the ranges of the luminosity and effective temperature we
can constrain here are those of RSG SN progenitors. The effective temperature
of the RSG SN progenitors is around 3500 K (Smartt, 2015) and we assume this
temperature for all the RSG SN progenitors. The lowest luminosity of the RSG
SN progenitors is observationally constrained to be
$\log(L/\mathrm{L}_{\odot})\simeq 4.4$ (Davies & Beasor, 2020), which is the
lowest luminosity we adopt in this paper. The maximum luminosity of RSGs is
$\log(L/\mathrm{L}_{\odot})\simeq 5.5$ (Davies et al., 2018; Davies & Beasor,
2020), while the estimated maximum luminosity of the RSG SN progenitors based
on the current RSG SN progenitor identification is
$\log(L/\mathrm{L}_{\odot})\simeq 5.2$ (e.g., Smartt et al., 2009). We adopt
the RSG luminosity of up to $\log(L/\mathrm{L}_{\odot})\simeq 5.5$.

Figure 5:  Rise times and peak luminosities estimated by changing the mass-
loss rates by a factor of 10 and 0.1. We exclude RSGs with $\log
L/\mathrm{L}_{\odot}=5.3,5.4,$ and 5.5 in this figure. The top panel is for
the de Jager mass-loss rate and the bottom panel is for the van Loon mass-loss
rate. $\varepsilon_{B}=0.001$ and SSA are assumed.

Figure 6:  Same as Fig. 5 but the free-free absorption is taken into account.

##  4 Constraining RSG mass-loss prescriptions with radio SN observations

In this section, we estimate the rise time and peak luminosity of SNe II in
radio expected from the RSG mass-loss prescriptions introduced in the previous
section and compare the results to the observed distribution.

The top panel of Fig. 3 shows the comparison using the de Jager RSG mass-loss
rate. Each line shows the rise time and peak luminosity range for a given RSG
luminosity with $M_{\mathrm{ej}}=5-20~{}\mathrm{M}_{\odot}$ and
$E_{\mathrm{ej}}=0.1-2~{}\mathrm{B}$ (cf. Fig. 1). The lines for the less
luminous RSGs, which have smaller mass-loss rates, locate more left in the
rise time and peak luminosity plane. The three gray lines in Fig. 3 are for
$\log(L/\mathrm{L}_{\odot})=5.3$, 5.4, and 5.5. As presented before, it is not
clear if RSGs with these luminosities explode as SNe. The red squares show the
rise times and peak luminosities of SNe II observed so far (Bietenholz et al.,
2020). The blue ellipse shows the region where 68% of SNe II are estimated to
locate in the rise time - peak luminosity plane based on all the observational
information including upper limits (Bietenholz et al., 2020). Most of the
observational points (red squares) locate at the regions that are expected
from the de Jager rate, but the theoretical prediction does not reach the low
peak luminosity regions included in the blue ellipse.

The bottom panel in Fig. 3 shows the same comparison as the top panel but for
the van Loon RSG mass-loss rate. The van Loon mass-loss rate is less dependent
on the luminosity than the de Jager mass-loss rate (Fig. 2). Thus, the
expected area to cover in the rise time and peak luminosity plane for the van
Loon mass-loss rate is smaller than that covered by the de Jagear mass-loss
rate. We can also find that the expected rise time and peak luminosity from
the van Loon mass-loss rate only cover the largest rise time and peak
luminosity ranges found from the observations. Overall, the van Loon mass-loss
rate does not explain the widely spread rise time and peak luminosity relation
in SNe II, and it does not explain the low peak luminosity that dominates SNe
II.

In both de Jager and van Loon mass-loss rates, we found that the theoretical
predictions do not reach the low peak luminosity expected from the
observations shown as the blue ellipse in Fig. 3. This could be partly because
of our assumption of $\varepsilon_{B}=0.01$. As discussed before,
$\varepsilon_{B}$ is uncertain and a smaller $\varepsilon_{B}$ can make the
expected peak luminosity lower. Fig. 4 shows the same estimates as in Fig. 3
but for $\varepsilon_{B}=0.001$. We can find that the rise time and peak
luminosity expected from the de Jager rate starts to cover some low peak
luminosity region, while the van Loon mass-loss rate still does not explain
the low peak luminosities.

In addition to the uncertainties in physics, we need to take the diversity in
the RSG mass-loss rates into account. Although the mass-loss prescriptions
provide one mass-loss rate for one luminosity, dispersion within a factor of
around 10 in the RSG mass-loss rate is known to exist (e.g., Mauron &
Josselin, 2011; Goldman et al., 2017). Fig. 5 shows the expected rise time and
peak luminosity regions with the mass-loss rate within a factor of 10 from the
de Jager and van Loon prescriptions in the case of $\varepsilon_{B}=0.001$.
The de Jager mass-loss rate reproduces the diverse rise time and peak
luminosity estimated from the observations, while the van Loon mass-loss rate
predict larger rise time and higher luminosity than those observed in SNe II.
The region covered in the rise time and peak luminosity plane is not as widely
spread as that covered by the de Jager mass-loss rate. As discussed before,
this is because the van Loon mass-loss rate does not strongly depend on
luminosity and the CSM density around RSG SN progenitors is not expected to
have much diversity.

We have not taken the effect of the free-free absorption into account so far.
The optical depth to the free-free absorption at 8.4 GHz becomes unity at
(Chevalier et al., 2006)

$t_{\mathrm{ff}}\simeq
6\left(\frac{\dot{M}/10^{-6}~{}\mathrm{M}_{\odot}~{}\mathrm{yr^{-1}}}{v_{\mathrm{CSM}}/10~{}\mathrm{km~{}s^{-1}}}\right)^{2/3}\left(\frac{T_{\mathrm{CSM}}}{10^{5}~{}\mathrm{K}}\right)^{-1/2}\left(\frac{v_{\mathrm{sh}}}{10^{4}~{}\mathrm{km~{}s^{-1}}}\right)^{-1}~{}\mathrm{days},$
(3)

where $T_{\mathrm{CSM}}$ is the CSM temperature. $T_{\mathrm{CSM}}$ is
uncertain, but we adopt $T_{\mathrm{CSM}}\simeq 10^{5}~{}\mathrm{K}$ in this
work (Lundqvist & Fransson, 1988; Chevalier et al., 2006). If the time when
$\tau_{\mathrm{SSA}}$ becomes unity is smaller than $t_{\mathrm{ff}}$, the
free-free absorption determines the radio LC peak. Fig. 6 shows the rise time
- peak luminosity relations where the free-free absorption is taken into
account. We can find that the models with large $\dot{M}$ are strongly
affected by free-free absorption. They show a break in the rise time - peak
luminosity relation and the region with the long rise time and low peak
luminosity (the right bottom part of the rise time - peak luminosity plane)
that is not covered by the SSA models is covered by taking the free-free
absorption into account. Our conclusion that the de Jager rate explains the
widely spread rise time and peak luminosity relation better, however, remains
unchanged.

##  5 Summary

We have compared the rise time and peak luminosity distribution of SNe II in
radio with those expected from the RSG mass-loss prescriptions of de Jager et
al. (1988) and van Loon et al. (2005). We found that the de Jager mass-loss
rate explains the rise time and peak luminosity of SNe II in radio well
considering the dispersion in the RSG mass-loss rates (Fig. 6). The van Loon
mass-loss rate generally predicts longer rise times and higher peak
luminosities than those found in SNe II. This is because of the relatively
high mass-loss rates that the van Loon prescription provides. The uncertainty
in physics and the dispersion in the RSG mass-loss rates are not likely to
fully resolve this discrepancy.

Another superiority of the de Jager mass-loss rate in explaining the radio
properties of SNe II is that it can explain the wide range of the rise times
and peak luminosities found in SNe II. The de Jager mass-loss rate can
reproduce the wide range because of its steep dependence on luminosity that
allows it to produce a variety of mass-loss rates with the RSG SN progenitor
luminosity range (Fig. 2). The de Jager mass-loss rate changes by a factor of
around 100 from the lowest luminosity RSG SN progenitors to the highest, while
the van Loon mass-loss rate changes only by a factor of around 10. The RSG
mass-loss prescriptions having as strong luminosity dependence as the de Jager
prescription is preferred to explain the diversity in SN II radio properties.

RSG mass-loss rates are an essential ingredient to understand the evolution of
massive stars towards SN explosions and SN properties. We have shown in this
paper that radio SN observations can provide a way to constrain the uncertain
RSG mass-loss prescriptions. The radio SN observations are still limited and
many more observations are required to constrain the rise time and peak
luminosity of SNe in radio. Further radio observations would lead to a better
understanding of mass loss in massive stars.

## Acknowledgements

I would like to thank the anonymous referee for the helpful comments that
improved this work. This work is supported by the Grants-in-Aid for Scientific
Research of the Japan Society for the Promotion of Science (JP18K13585,
JP20H00174).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Adams et al. (2017) Adams S. M., Kochanek C. S., Gerke J. R., Stanek K. Z., 2017, MNRAS, 469, 1445
  * Bietenholz et al. (2020) Bietenholz M. F., Bartel N., Argo M., Dua R., Ryder S., Soderberg A., 2020, arXiv e-prints, p. arXiv:2011.11737
  * Björnsson & Fransson (2004) Björnsson C.-I., Fransson C., 2004, ApJ, 605, 823
  * Chevalier (1982a) Chevalier R. A., 1982a, ApJ, 258, 790
  * Chevalier (1982b) Chevalier R. A., 1982b, ApJ, 259, 302
  * Chevalier (1998) Chevalier R. A., 1998, ApJ, 499, 810
  * Chevalier & Fransson (2006) Chevalier R. A., Fransson C., 2006, ApJ, 651, 381
  * Chevalier et al. (2006) Chevalier R. A., Fransson C., Nymark T. K., 2006, ApJ, 641, 1029
  * Davies & Beasor (2020) Davies B., Beasor E. R., 2020, MNRAS, 493, 468
  * Davies et al. (2018) Davies B., Crowther P. A., Beasor E. R., 2018, MNRAS, 478, 3138
  * Förster et al. (2018) Förster F., et al., 2018, Nature Astronomy, 2, 808
  * Fuller (2017) Fuller J., 2017, MNRAS, 470, 1642
  * Georgy (2012) Georgy C., 2012, A&A, 538, L8
  * Goldman et al. (2017) Goldman S. R., et al., 2017, MNRAS, 465, 403
  * Horiuchi et al. (2011) Horiuchi S., Beacom J. F., Kochanek C. S., Prieto J. L., Stanek K. Z., Thompson T. A., 2011, ApJ, 738, 154
  * Kamble et al. (2016) Kamble A., et al., 2016, ApJ, 818, 111
  * Kochanek et al. (2008) Kochanek C. S., Beacom J. F., Kistler M. D., Prieto J. L., Stanek K. Z., Thompson T. A., Yüksel H., 2008, ApJ, 684, 1336
  * Lundqvist & Fransson (1988) Lundqvist P., Fransson C., 1988, A&A, 192, 221
  * Maeda (2012) Maeda K., 2012, ApJ, 758, 81
  * Maeda (2013) Maeda K., 2013, ApJ, 762, 14
  * Matzner & McKee (1999) Matzner C. D., McKee C. F., 1999, ApJ, 510, 379
  * Mauron & Josselin (2011) Mauron N., Josselin E., 2011, A&A, 526, A156
  * Meynet et al. (2015) Meynet G., et al., 2015, A&A, 575, A60
  * Pejcha & Prieto (2015) Pejcha O., Prieto J. L., 2015, ApJ, 806, 225
  * Quataert & Shiode (2012) Quataert E., Shiode J., 2012, MNRAS, 423, L92
  * Shivvers et al. (2017) Shivvers I., et al., 2017, PASP, 129, 054201
  * Smartt (2015) Smartt S. J., 2015, Publ. Astron. Soc. Australia, 32, e016
  * Smartt et al. (2009) Smartt S. J., Eldridge J. J., Crockett R. M., Maund J. R., 2009, MNRAS, 395, 1409
  * Soderberg et al. (2006) Soderberg A. M., Chevalier R. A., Kulkarni S. R., Frail D. A., 2006, ApJ, 651, 1005
  * Sukhbold & Adams (2020) Sukhbold T., Adams S., 2020, MNRAS, 492, 2578
  * Walmswell & Eldridge (2012) Walmswell J. J., Eldridge J. J., 2012, MNRAS, 419, 2054
  * Weiler et al. (1992) Weiler K. W., van Dyk S. D., Pringle J. E., Panagia N., 1992, ApJ, 399, 672
  * Weiler et al. (2002) Weiler K. W., Panagia N., Montes M. J., Sramek R. A., 2002, ARA&A, 40, 387
  * Yaron et al. (2017) Yaron O., et al., 2017, Nature Physics, 13, 510
  * Yoon & Cantiello (2010) Yoon S.-C., Cantiello M., 2010, ApJ, 717, L62
  * de Jager et al. (1988) de Jager C., Nieuwenhuijzen H., van der Hucht K. A., 1988, A&AS, 72, 259
  * van Loon et al. (2005) van Loon J. T., Cioni M. R. L., Zijlstra A. A., Loup C., 2005, A&A, 438, 273


We study novel simulation-like preorders for quotienting nondeterministic Büchi automata.

We define fixed-word delayed simulation,
a new preorder coarser than delayed simulation.
We argue that fixed-word simulation is the coarsest forward simulation-like preorder which can be used for quotienting Büchi automata,
thus improving our understanding of the limits of quotienting.
Also, we show that computing fixed-word simulation is PSPACE-complete.

On the practical side,
we introduce proxy simulations,
which are novel polynomial-time computable preorders sound for quotienting.
In particular, delayed proxy simulation induce quotients that can be smaller by an arbitrarily large factor than direct backward simulation.
We derive proxy simulations as the product of a theory of refinement transformers:
A refinement transformer maps preorders nondecreasingly,
preserving certain properties.
We study under which general conditions refinement transformers are sound for quotienting.

We thank Richard Mayr and Patrick Totzke for helpful discussions,
and two anonymous reviewers for their valuable feedback.

[1]
Abdulla, P., Chen, Y.F., Holik, L., Vojnar, T.: Mediating for
Reduction. In: FSTTCS. pp. 1–12. Schloss
Dagstuhl–Leibniz-Zentrum fuer Informatik (2009)

[2]
Aziz, A., Singhal, V., Swamy, G.M., Brayton, R.K.: Minimizing Interacting
Finite State Machines. Tech. Rep. UCB/ERL M93/68, UoC, Berkeley (1993)

[3]
Clemente, L., Mayr, R.: Multipebble Simulations for Alternating Automata -
(Extended Abstract). In: CONCUR. LNCS, vol. 6269, pp. 297–312.
Springer-Verlag (2010), <http://dx.doi.org/10.1007/978-3-642-15375-4_21>

[4]
Dill, D.L., Hu, A.J., Wont-Toi, H.: Checking for Language Inclusion Using
Simulation Preorders. In: CAV. LNCS, vol. 575. Springer-Verlag (1991),

[5]
Etessami, K.: A Hierarchy of Polynomial-Time Computable Simulations for
Automata. In: CONCUR. LNCS, vol. 2421, pp. 131–144. Springer-Verlag (2002),

[6]
Etessami, K., Wilke, T., Schuller, R.A.: Fair Simulation Relations, Parity
Games, and State Space Reduction for Büchi Automata. SIAM J. Comput.
34(5),  1159–1175 (2005),

[7]
Fritz, C., Wilke, T.: Simulation Relations for Alternating Büchi
Automata. Theor. Comput. Sci.  338(1-3),  275–314 (2005),

[8]
Gramlich, G., Schnitger, G.: Minimizing NFA's and Regular Expressions.
Journal of Computer and System Sciences  73(6),  908–923 (2007),

[9]
Gurumurthy, S., Bloem, R., Somenzi, F.: Fair Simulation Minimization. In:
CAV. LNCS, vol. 2404, pp. 610–624. Springer-Verlag (2002),

[10]
Henzinger, T.A., Kupferman, O., Rajamani, S.K.: Fair Simulation. Information
and Computation  173,  64–81 (2002),

[11]
Henzinger, T.A., Rajamani, S.K.: Fair Bisimulation. In: TACAS. LNCS, vol.
1785, pp. 299–314. Springer-Verlag (2000),

[12]
Jiang, T., Ravikumar, B.: Minimal NFA Problems are Hard. In: Albert, J.,
Monien, B., Artalejo, M. (eds.) Automata, Languages and Programming, Lecture
Notes in Computer Science, vol. 510, pp. 629–640. Springer Berlin /
Heidelberg (1991), <http://dx.doi.org/10.1007/3-540-54233-7_169>

[13]
Juvekar, S., Piterman, N.: Minimizing Generalized Büchi Automata. In:
CAV. LNCS, vol. 4414, pp. 45–58. Springer-Verlag (2006),

[14]
Kupferman, O., Vardi, M.: Verification of Fair Transition Systems. In: CAV,
LNCS, vol. 1102, pp. 372–382. Springer-Verlag (1996),

[15]
Kupferman, O., Vardi, M.: Weak Alternating Automata Are Not That Weak. ACM
Trans. Comput. Logic  2,  408–429 (Jul 2001),

[16]
Lynch, N.A., Vaandrager, F.W.: Forward and Backward Simulations. Part I:
Untimed Systems. Information and Computation  121(2),  214–233 (1995),

[17]
Milner, R.: Communication and Concurrency. Prentice-Hall (1989)

[18]
Raimi, R.S.: Environment Modeling and Efficient State Reachability Checking.
Ph.D. thesis, The University of Texas at Austin (1999)

[19]
Schewe, S.: Beyond Hyper-Minimisation—Minimising DBAs and DPAs is
NP-Complete. In: Lodaya, K., Mahajan, M. (eds.) FSTTCS. LIPIcs, vol. 8, pp.
400–411. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl,
Germany (2010), <http://drops.dagstuhl.de/opus/volltexte/2010/2881>

[20]
Somenzi, F., Bloem, R.: Efficient Büchi Automata from LTL Formulae. In:
CAV, LNCS, vol. 1855, pp. 248–263. Springer-Verlag (2000),

[21]
Vardi, M.: Alternating Automata and Program Verification. In: Computer
Science Today, LNCS, vol. 1000, pp. 471–485. Springer-Verlag (1995),