# Artificial magnetic field for synthetic quantum matter without dynamical
modulation

Tomoki Ozawa  Advanced Institute for Materials Research, Tohoku University,
Sendai 980-8577, Japan  RIKEN Interdisciplinary Theoretical and Mathematical
Sciences Program, Wako 351-0198, Japan

###### Abstract

We propose an all-static method to realize an artificial magnetic field for
charge neutral particles without introducing any time modulation. Our proposal
consists of one-dimensional tubes subject to harmonic trapping potentials with
shifted centers. We show that this setup realizes an artificial magnetic field
in a hybrid real-momentum space. We discuss how characteristic features of
particles in a magnetic field, such as chiral edge states and the quantized
Hall response, can be observed in this setup. We find that the mean-field
ground state of bosons in this setup in the presence of long-range
interactions in physical real space can have quantized vortices in the hybrid
real-momentum space; such a state with vortices exhibits a supersolid
structure in the physical real space. Our method can be applied to a variety
of synthetic quantum matter, including ultracold atomic gases, coupled
photonic cavities, coupled waveguides, and exciton-polariton lattices.

Introduction Quantum simulation has been a central theme in the research of
synthetic quantum matter based on atomic, molecular, and optical systems, such
as ultracold atomic gases, trapped ions, polaritons, and photonics Bloch:2012
; Blatt:2012 ; Guzik:2012 ; Carusotto:2013 ; Georgescu:2014 ; Bloch:2018 ;
Boulier:2020 . In particular, simulation of topological phenomena using
synthetic quantum matter has attracted a considerable attention in the recent
years Hasan:2010 ; Qi:2011 ; Chiu:2018 ; Armitage:2018 ; Cooper:2019 ;
Ozawa:2019 . A prototypical example of topological phenomena is the quantum
Hall effect, which was first found in a two-dimensional electron gas under a
magnetic field Klitzing:1980 ; Tsui:1982 ; Yoshioka:Book . To simulate the
quantum Hall effect and other topological phenomena, one often needs to
simulate physics of charged particles in a magnetic field, which is not
straightforward because atoms in ultracold gases and photons in cavities and
waveguides are charge-neutral.

Various methods have been employed to simulate effects of a magnetic field,
namely to realize an artificial magnetic field Aidelsburger:2018 . In
ultracold atomic gases, the earliest example is to rotate the system,
exploiting the similarity between the Coriolis force and the Lorentz force
Madison:2000 ; Cooper:2008 ; Fetter:2009 . Other examples include Floquet
engineering GoldmanDalibard:2014 , light-induced gauge potential Goldman:2014
, and synthetic dimensions Ozawa:2019SD . Most of the existing proposals
involve adding fast time modulation and/or exploiting different energy and
time scales present in the system to explore adiabatic physics in a low energy
subspace. These methods often suffer from instability and heating as a result
of time modulation and spontaneous emissions. A notable exception is the
method using spinor condensates in a quadrupole magnetic field where spatially
dependent spin texture results in an artificial magnetic field Ho:1996 ;
Pietila:2009 ; Choi:2013 ; Ray:2014 . The situation is similar in photonics
and optics. For light close to optical frequencies, magneto-optical effects
are generally weak and thus, similar to ultracold atomics gases, one needs to
realize an artificial magnetic field to explore the physics of quantum Hall
effects Rechtsman:2013 ; Yuan:2018 ; Lustig:2019 . Static realizations of
topological models in photons are in the class of quantum spin-Hall insulators
where two degrees of freedom feel opposite magnetic fields and thus
topological protection relies upon decoupling of different degrees of freedom
Hafezi:2013 ; Schine:2016 .

Figure 1: Illustration of the scheme. (a) One-dimensional tubes elongated in
the $x$ direction are aligned along the $y$ direction. Harmonic trapping
potentials, depicted by parabolas, are applied to the tubes. The trapping
centers are on the dashed line, which is tilted with respect to the $x$ axis
by an angle $\theta$. (b) Setup viewed from top. Tubes are located at
$y=ny_{0}$, where $n$ is an integer. The trapping centers are on the dashed
line $y=x\tan\theta$. (c) In the $x$-$p_{y}$ plane, an artificial magnetic
field of strength $B=-\tan\theta$ is realized. We have the geometry of a
cylinder, where the $p_{y}$ direction is periodic with periodicity
$2\pi/y_{0}$.

In this paper, we propose a completely static way to realize an artificial
magnetic field regardless of the spin of the particles, which can be applied
to a variety of synthetic quantum matter. Our proposal is based on a well-
known fact that the Hamiltonian of a charged quantum particle in two
dimensions under a uniform magnetic field is equivalent to a set of harmonic
oscillators with shifted centers LL:QM . This fact implies, conversely, that a
set of harmonic oscillators with shifted centers can be viewed as a charged
particle in a uniform magnetic field. We pursue this analogy further and
consider a set of one-dimensional systems elongated along the $x$ direction
and align these one-dimensional tubes in the $y$-direction and place them
under harmonic trapping potentials with shifted centers (Fig. 1). We find that
the Hamiltonian of a charged particle in a magnetic field is realized in the
$x$-$p_{y}$ plane where $p_{y}$ is the momentum along the $y$ direction. Even
though the proposal looks deceptively simple, as we discuss in detail, we can
still observe most of the phenomena characteristic of charged particles in a
magnetic field, such as the chiral edge states and quantized Hall response. We
note that there is no dynamical component in the proposal, and, in particular,
our proposal does not break the time-reversal symmetry in the physical $x$-$y$
plane. Our proposal is thus expected to be more stable and free from heating
compared with existing proposals which rely on dynamical modulations.

We also explore the effect of interparticle interactions in the artificial
magnetic field. Weakly interacting bosons under a magnetic field form vortices
Madison:2000 ; Cooper:2008 ; Fetter:2009 ; Spielman:2009 ; Lin:2009 . However,
a contact interaction, common in ultracold atomic gases and photons with
optical nonlinearity, translates into extremely long-range interaction in the
$p_{y}$ direction, and we find that vortices do not form under ordinary
contact interactions. On the other hand, by including nearest-neighbor
interaction in $y$ direction we find that vortices can form in the $x$-$p_{y}$
plane. Such a vortex structure corresponds to supersolid structure in the
physical $x$-$y$ plane, where the phase coherence of a condensate is
maintained and the discrete translational symmetry is broken. Finally, we
discuss possible experimental platforms where the proposal can be implemented.

##  I Setup

We consider a set of uncoupled one-dimensional tubes elongated in the $x$
direction and aligned in the $y$ direction, where each tube is subject to a
harmonic trapping potential. The center of the trapping potential shifts
between tubes by a constant amount $x_{0}$. We assume that tubes are separated
by a constant amount $y_{0}$ in the $y$ direction. See Fig. 1 for an
illustration. We use the second-quantized formalism, where
$\hat{\psi}^{\dagger}_{n}(x)$ and $\hat{\psi}_{n}(x)$ are the creation and
annihilation operators, respectively, of a particle at coordinate
$(x,ny_{0})$, where $n$ indicates that the particle is in the $n$th tube. The
single-particle Hamiltonian is

$\displaystyle\hat{H}=$ $\displaystyle\sum_{n}\int
dx\left[\frac{1}{2m}\left(\partial_{x}\hat{\psi}_{n}^{\dagger}(x)\right)\left(\partial_{x}\hat{\psi}_{n}(x)\right)\right.$
$\displaystyle+\left.\frac{1}{2}m\omega^{2}(x-nx_{0})^{2}\hat{\psi}_{n}^{\dagger}(x)\hat{\psi}_{n}(x)\right],$
(1)

Note that there is no motion (coupling) between different tubes in the $y$
direction. By performing Fourier transformation in the $y$ direction,

$\displaystyle\hat{\psi}_{n}(x)=\sqrt{\frac{y_{0}}{2\pi}}\int_{-\pi/y_{0}}^{\pi/y_{0}}dp_{y}e^{iny_{0}p_{y}}\hat{\psi}(x,p_{y}),$
(2)

the Hamiltonian can be rewritten as

$\displaystyle\hat{H}=$ $\displaystyle\int
dx\int_{-\pi/y_{0}}^{\pi/y_{0}}dp_{y}\left[\frac{1}{2m}\left(\partial_{x}\hat{\psi}^{\dagger}(x,p_{y})\right)\left(\partial_{x}\hat{\psi}(x,p_{y})\right)\right.$
$\displaystyle+\left.\frac{1}{2}\frac{m\omega^{2}}{\tan^{2}\theta}\left\\{(i\partial_{p_{y}}+x\tan\theta)\hat{\psi}^{\dagger}(x,p_{y})\right\\}\right.$
$\displaystyle\hskip
56.9055pt\left.\left\\{(-i\partial_{p_{y}}+x\tan\theta)\hat{\psi}(x,p_{y})\right\\}\right].$
(3)

In the first-quantized form, this Hamiltonian is

$\displaystyle
H=-\frac{1}{2m}\partial_{x}^{2}+\frac{m\omega^{2}}{2\tan^{2}\theta}\left(-i\partial_{p_{y}}+x\tan\theta\right)^{2}$
(4)

This is nothing but the Hamiltonian of a charged particle in the $x$-$p_{y}$
plane under a magnetic vector potential $\mathbf{A}(x,p_{y})=(0,-x\tan\theta)$
in the Landau gauge, corresponding to the artificial magnetic field
$B=-\tan\theta$ penetrating the plane. (We have set the charge $e$, speed of
light $c$, and Plank constant $\hbar$ to be unity.) We have chosen to use
$-i\partial_{p_{y}}=-y$ as the “momentum” operator along the $p_{y}$
direction, regarding $p_{y}$ as the synthetic “position.” The mass is
anisotropic; the mass in the $x$ direction is the original mass $m_{x}=m$,
whereas the mass in the $p_{y}$ direction is
$m_{p_{y}}=\tan^{2}\theta/m\omega^{2}$. This anisotropy can be absorbed by
rescaling units of length in two directions and does not pose any problem in
the rest of the discussion. Another noticeable feature here is that the
$x$-$p_{y}$ plane is a cylinder, open along the $x$ direction and periodic
along the $p_{y}$ direction.

Let us comment on the time-reversal symmetry in our system. Our original
Hamiltonian (1) does not break the physical time-reversal symmetry. For a
single component (spinless) system, the time-reversal operation flips the sign
of the momentum $\mathbf{p}\to-\mathbf{p}$. Since we want to look at the
physics related to the artificial magnetic field in the $x$-$p_{y}$ plane, the
relevant effective time-reversal symmetry is now to flip their reciprocal
variables $\mathcal{T}:(p_{x},-y)\to(-p_{x},y)$; the effective time-reversal
symmetry is the simultaneous operation of flipping of $p_{x}$ and the parity
in $y$. Our Hamiltonian does break the parity in the $y$ direction due to the
shifted harmonic trapping potentials, which enables us to realize the
artificial magnetic field in the $x$-$p_{y}$ plane without breaking the
physical time-reversal symmetry.

Our method is related to the idea of the topological charge pumping, where a
momentum is replaced by an external parameter; by sweeping the parameter one
obtains topological properties of the original Hamiltonian Thouless:1983 ;
Nakajima:2016 ; Lohse:2016 . Our method, however, is crucially different from
charge pumping in that we replace a momentum $p_{y}$ by a physical real
coordinate $-y$, and in this way we succeed in recovering the full dynamics in
all dimensions, whereas no dynamics occurs in charge pumping in the direction
replaced by an external parameter.

Figure 2: Energy level and the Hall response. (a) The energy level of each
tube at position $n$, where length of each tube is $20a_{\mathrm{osc}}$ and we
employ an open boundary condition. We choose $\tan\theta=1$, and the energy is
in units of the cyclotron frequency $\omega$. Solid lines are to guide the
eye. (b) The center-of-mass motion of wave packets when a synthetic electric
field $E_{x}=0.1\omega/a_{\mathrm{osc}}$ is applied. Three solid lines show
the time evolution of the center of mass in the $p_{y}$ direction for
$\tan\theta=2/3$ (top, green), $\tan\theta=1$ (middle, blue), and
$\tan\theta=1/2$ (bottom, red). The dashed lines going through the solid lines
are theoretical predictions from Eq. (5) with $\mathcal{C}=1$. The units of
time are $1/\omega$, and the units of $p_{y}$ are $1/y_{0}$.

##  II Chiral edge states

We now analyze how characteristic features of charged particles in a magnetic
field can be observed in our setup in the $x$-$p_{y}$ plane. We first point
out that the cyclotron frequency of the artificial magnetic field is $\omega$,
as expected. This can be derived from the definition of the cyclotron
frequency $B/\tilde{m}$; as the mass $\tilde{m}$, we need to use a geometric
mean of anisotropic effective mass in two directions
$\tilde{m}=\sqrt{m_{x}m_{p_{y}}}$. Using $B=-\tan\theta$, we then obtain
$\omega=|B|/\tilde{m}$. This observation implies the well-known fact that the
energy level of harmonic oscillators in the $x$-$y$ plane, which is
$\omega(l+1/2)$, where $l=0,1,2,\cdots$ is a nonnegative integer, is nothing
but the Landau-level energy spectrum. High degeneracy of Landau levels come
from the many uncoupled tubes present in the system, which all have the same
energy spectrum.

In the presence of edges, charged particles in two dimensions under a magnetic
field should have chiral edge states. In Fig. 2(a), we plot the energy level
of the system in the presence of sharp edges in the $x$ direction as a
function of $n$, which is the reciprocal variable of $p_{y}$. We consider 21
tubes and the length of each tube in the $x$ direction is
$20a_{\mathrm{osc}}$, where $a_{\mathrm{osc}}\equiv 1/\sqrt{m\omega}$ is the
oscillator length. We can see that at large values of $|n|$, the dispersion
goes upward; these states are localized at the edges and they propagate in one
direction. The existence of chiral edge states can be probed from the dynamics
of wave packets localized at the edges as we show below. Since there is no
coupling between different tubes, the population of particles in each tube
does not change; in real space, the propagation of chiral edge states in the
$p_{y}$ direction appears as an appropriate evolution of the relative phases
in each tube, which can be probed, for example, in ultracold atomic gases
through time-of-flight imaging.

We now numerically demonstrate how the chiral edge mode in the $x$-$p_{y}$
plane can be observed through wave packet dynamics. We choose the strength of
the artificial magnetic field to be $B=-\tan\theta=-1$ for the numerical
simulation. As an initial state, we consider three wave packets, located at
the left edge, center, and the right edge of the system. To consider physics
relevant to a single band or band gap, we project the wave packet onto the
lowest level of each tube. Experimentally, in ultracold gases, such an initial
state can be prepared by creating a wave packet out of single-particle ground
states of each tube by starting from Bose-Einstein condensates in each tube.
In photonics, where light can be injected from the outside, we can insert
light with frequency corresponding to the lowest Landau level to prepare a
similar initial condition. The wave function of the initial state is plotted
in Fig. 3(a) and Fig. 3(b), in both the $x$-$y$ and $x$-$p_{y}$ planes. We
then let the wave packet evolve freely. The wave functions in the $x$-$p_{y}$
plane after time $t=1/\omega$, $t=2/\omega$, and $t=3/\omega$ are plotted in
Figs. 3(c)-3(e). We observe that the wave packet at the left edge goes upward,
whereas the wave packet at the right edge goes downward, showing that the edge
state goes clockwise in a chiral manner. This motion of the wave packet at the
edge is consistent with the artificial magnetic field in the $x$-$p_{y}$
plane; the strength of the magnetic field is $B=-\tan\theta<0$, which implies
that the cyclotron motion inside the bulk of the system is counterclockwise,
and the edge state goes clockwise. We note that the absolute value of the wave
function plotted in the $x$-$y$ plane looks identical during time evolution.
This is because our initial state is projected onto the lowest Landau level,
and thus the wave function in each tube is an eigenstate of the Hamiltonian
for each tube, and thus the only time evolution appears in the phase of the
wave function. The evolution of the phase in the $x$-$y$ plane appears as
chiral edge motion in the $x$-$p_{y}$ plane.

Figure 3: Dynamics of a wave packet and chiral edge modes. The units of length
in the $x$ direction are $a_{\mathrm{osc}}$, and the units of $p_{y}$ are
$1/y_{0}$. (a) The absolute value of the initial wave function is plotted in
the $x$-$y$ plane, where each $n$ represents a tube located at the position
$y=ny_{0}$. (b) The initial state in the $x$-$p_{y}$ plane. The lower panel
shows the wave function in the $x$-$p_{y}$ plane after time evolution for
duration (c) $t=1/\omega$, (d) $t=2\omega$, and (e) $t=3/\omega$. Chiral edge
motion is clearly visible.

We should also comment on the direction of the edge mode with respect to the
energy spectrum, see Fig. 2(a). One sees that, in the energy spectrum, the
energy goes up at large values of $n$. These states at large values of $n>0$
are localized at the right edge. Since the energy goes up, one might want to
conclude that the group velocity is positive at the right edge. However, as
one sees from Fig. 3, this is not the case; the group velocity is negative at
the right edge. This apparent discrepancy is because one needs to consider
$-n$ as the conjugate variable of $p_{y}$. Namely, the horizontal axis of Fig.
2(a) should be flipped when one wants to find the correct group velocity. By
flipping the horizontal axis of Fig. 2(a), one obtains the correct group
velocity, which is negative at the right edge.

##  III Quantized Hall response

Since Landau levels are topological, the Hall response should be quantized,
which can be a signature of nontrivial topology obtained from the bulk
property of the system. The topological Chern number $\mathcal{C}$ of Landau
levels with $B<0$ is one, $\mathcal{C}=1$; we now investigate if this unit
Chern number can be observed in our setup.

One way to see the topological nature of the system is through Laughlin charge
pumping Laughlin:1981 ; on a two-dimensional cylinder with a perpendicular
magnetic field, when a magnetic flux quantum is inserted through the cylinder,
the state on the cylinder moves along the cylinder by one unit. As we
naturally have a cylindrical structure in the $x$-$p_{y}$ plane, it is natural
to ask how Laughlin charge pumping appears in our setup. As we see below,
quantized pumping in our system occurs naturally due to construction. A
magnetic flux can be inserted through the cylinder in the $x$-$p_{y}$ plane by
shifting the center of the harmonic trapping potential in the $x$-$y$ plane.
If we shift the central position of the harmonic trapping potentials by
$\delta x$, the harmonic trapping potential acts as
$\frac{1}{2}m\omega^{2}(x-nx_{0}-\delta x)^{2}$ for a given tube labeled by an
integer $n$. In the $x$-$p_{y}$ plane, this shift corresponds to the shift of
the magnetic vector potential to $\mathbf{A}(x,p_{y})=(0,-x\tan\theta+\delta
x\tan\theta)$. The additional contribution of $\delta x\tan\theta$ is a
constant term in the magnetic vector potential, and thus if we integrate it
over the $p_{y}$ direction, we see that the cylinder in the $x$-$p_{y}$ plane
is penetrated by a magnetic flux $2\pi\delta x\tan\theta/y_{0}=2\pi\delta
x/x_{0}$. Therefore, inserting a magnetic flux quantum of $2\pi$ corresponds
to taking $\delta x=x_{0}$. This shift of $\delta x=x_{0}$ corresponds to no
change in the overall Hamiltonian except that a tube previously labeled by an
integer $n$ now acts as a tube with a label $n+1$. Thus, as in the Laughlin
pumping, the overall Hamiltonian comes back to the original one after the
shift of $\delta x=x_{0}$. In the $x$-$y$ plane, the shift of $\delta x=x_{0}$
corresponds to moving the entire system by an amount $\delta x=x_{0}$ to the
right, thus any state should also move horizontally by $x_{0}$. This implies
that any state on the cylinder of the $x$-$p_{y}$ plane should also move along
the cylinder by $x_{0}$, which is the manifestation of the Laughlin charge
pumping in our setup. Another way to see the topological nature of the system
is through dynamical response of the system. In ultracold gases, quantized
Hall response can be observed through the center-of-mass response upon
application of a force (synthetic electric field)
$\mathbf{E}=(E_{x},E_{p_{y}})$ Dauphin:2013 ; Aidelsburger:2015 ; Price:2016 .
When a synthetic electric field $E_{i}$ is applied along the $i$ direction,
the center-of-mass velocity is Price:2016
$v^{\mathrm{C.M.}}_{i}=-\epsilon_{ij}\frac{2\pi}{A_{\mathrm{BZ}}}E_{j}\mathcal{C}$,
where $A_{\mathrm{BZ}}$ is the area of the Brillouin zone and
$\epsilon_{xp_{y}}=-\epsilon_{p_{y}x}=1$. For our setup of a uniform magnetic
field in the $x$-$p_{y}$ plane, there is no translational symmetry in the $x$
direction, and therefore we need to consider a magnetic unit cell and the
magnetic Brillouin zone. A magnetic unit cell is an arbitrary rectangle in the
$x$-$p_{y}$ plane which encloses one magnetic flux quantum. Since the strength
of the artificial magnetic field is $|B|=\tan\theta$, the area of a magnetic
unit cell is $2\pi/\tan\theta$. Then the area of the magnetic Brillouin zone
is $(2\pi)^{2}$ divided by the area of a magnetic unit cell and thus
$A_{\mathrm{BZ}}=2\pi\tan\theta$. Therefore, the center-of-mass velocity
should follow

$\displaystyle
v^{\mathrm{C.M.}}_{i}=-\epsilon_{ij}E_{j}\mathcal{C}/\tan\theta=\epsilon_{ij}E_{j}\mathcal{C}/B.$
(5)

We should then be able to determine the Chern number of the Landau levels in
the $x$-$p_{y}$ plane by monitoring the center-of-mass motion. As we see
below, the quantized response upon adding $E_{p_{y}}$ is nothing but the
Laughlin charge pumping we saw above and turns out to be trivially satisfied,
whereas the response upon $E_{x}$ can provide observable signatures of the
quantized Hall response.

We first consider the Hall response upon adding $E_{p_{y}}$, which can be
applied by making the artificial magnetic vector potential time dependent
$\mathbf{A}=(0,-x\tan\theta-E_{p_{y}}t)$. Such a time dependence can again be
introduced by changing the position of the minima of the harmonic trapping
potentials. Changing the location of the trap minima from $nx_{0}$ to
$nx_{0}+\delta x$, in a $y$-independent manner, the magnetic vector potential
in the $x$-$p_{y}$ plane becomes $\mathbf{A}=(0,-x\tan\theta+\delta
x\tan\theta)$ as before. This implies that the desired magnetic vector
potential can be created by moving $\delta x$ as $\delta
x=-E_{p_{y}}t/\tan\theta$. Now, the quantized Hall response through the
center-of-mass velocity has a clear physical meaning. When the motion of the
trapping potential is slow enough, the center of mass should follow the motion
of the trap minima $\delta x(t)$. Therefore, the center-of-mass velocity in
the $x$ direction is $v^{\mathrm{C.M.}}_{x}=\partial_{t}\delta
x(t)=-E_{p_{y}}/\tan\theta$; from Eq. (5), this implies the correct relation
$\mathcal{C}=1$.

Adding $E_{x}$, the center-of-mass velocity follows
$v^{\mathrm{C.M.}}_{p_{y}}=E_{x}\mathcal{C}/\tan\theta$, which can provide a
more nontrivial observable signature of topology. A synthetic electric field
$E_{x}$ can be applied by adding a linear potential gradient $-E_{x}x$ to the
Hamiltonian. We numerically simulate the time evolution of a wave packet
located at the center, projected onto the lowest Landau level, under the field
$E_{x}$. In Fig. 2(b), we plot the center-of-mass position of the wave packet
in the $p_{y}$ direction under time evolution for different values of
$B=-\tan\theta$. We also plot the theoretical prediction from Eq. (5) assuming
$\mathcal{C}=1$, which agrees well with the numerical simulation of $\langle
p_{y}\rangle$. The numerical simulation also exhibits oscillation around the
theoretical line with oscillation period $2\pi/\omega$, which is a detailed
structure of the cyclotron motion. Thus, the quantized Hall response can be
observed from the center-of-mass motion in the $p_{y}$ direction under the
application of $E_{x}$, and the Chern number can be estimated from its slope.

##  IV Interaction and vortices

We now turn to interacting many-body cases. An important feature is that a
typical short-range interaction in the $x$-$y$ plane translates into long-
ranged and unconventional interaction in the $p_{y}$ direction. In particular,
the formation of vortices, which is a characteristic feature of weakly
interacting bosons in a magnetic field Madison:2000 ; Cooper:2008 ;
Fetter:2009 ; Spielman:2009 ; Lin:2009 , cannot be observed in the $x$-$p_{y}$
plane only with a contact interaction in the $x$-$y$ plane. Nonetheless, we
find that, by adopting a long-range interaction in the $y$ direction, vortices
can still form in our setup as we discuss below. We now focus on bosons with
the mean-field interaction and ask if vortices can form in $x$-$p_{y}$ plane.
We use description in terms of the Gross-Pitaevskii equation with the
condensate wave function $\Psi_{n}(x,t)$:

$\displaystyle i\partial_{t}\Psi_{n}(x,t)=H_{\mathrm{yhop}}+H_{\mathrm{int}}+$
$\displaystyle\left\\{-\frac{\partial_{x}^{2}}{2m}+\frac{1}{2}m\omega^{2}(x-nx_{0})^{2}+\frac{1}{2}m\tilde{\omega}^{2}x^{2}\right\\}\Psi_{n}(x,t),$
(6) $\displaystyle H_{\mathrm{yhop}}=J_{y}(\Psi_{n-1}(x,t)+\Psi_{n+1}(x,t)),$
(7) $\displaystyle H_{\mathrm{int}}=g_{s}|\Psi_{n}(x,t)|^{2}\Psi_{n}(x,t)$
$\displaystyle\hskip
28.45274pt+g_{1}(|\Psi_{n-1}(x,t)|^{2}+|\Psi_{n+1}(x,t)|^{2})\Psi_{n}(x,t).$
(8)

We additionally included an overall weak harmonic potential with the
oscillator frequency $\tilde{\omega}$ which confines the particles. We take
$\tilde{\omega}$ to be much smaller than $\omega$ so that its only effect is
to prevent particles from flying away. We also added a small intertube hopping
term, $H_{\mathrm{yhop}}$, which allows an exchange of particles between
different tubes. This hopping term allows the system to establish an overall
phase coherence. In the $x$-$p_{y}$ plane, such an intertube hopping term
enters as a weak periodic potential of the form $2J_{y}\cos(y_{0}p_{y})$. In
the numerical calculations below, we take $\tilde{a}_{\mathrm{osc}}\equiv
1/\sqrt{m\tilde{\omega}}=10a_{\mathrm{osc}}$ and $J_{y}=-0.01\omega$. The
interaction term $H_{\mathrm{int}}$ contains the contact interaction with
strength $g_{s}$ and also a nearest-neighbor (long-range) interaction between
tubes with strength $g_{1}$.

Although the interaction $H_{\mathrm{int}}$ looks innocuous in the $x$-$y$
plane, it translates into rather unusual infinite and long-ranged interactions
in the $x$-$p_{y}$ plane. Writing the Fourier transform of $\Psi_{n}(x)$ along
the $y$ direction as $\Psi(x,p_{y})$, and considering the corresponding
second-quantized creation and annihilation operators
$\hat{\Psi}^{\dagger}(x,p_{y})$ and $\hat{\Psi}(x,p_{y})$, the interaction
term written in $x$-$p_{y}$ operators, up to an overall normalization factor,
reads

$\displaystyle\sum_{p_{1}+p_{2}=p_{3}+p_{4}}$
$\displaystyle\left\\{\frac{g_{s}}{2}+g_{1}\cos\left(\frac{y_{0}(p_{1}-p_{2})}{2}\right)\cos\left(\frac{y_{0}(p_{3}-p_{4})}{2}\right)\right\\}$
$\displaystyle\int dx\,$
$\displaystyle\hat{\Psi}^{\dagger}(x,p_{1})\hat{\Psi}^{\dagger}(x,p_{2})\hat{\Psi}(x,p_{3})\hat{\Psi}(x,p_{4}).$
(9)

This is a position conserving interaction considering $p_{y}$ as a synthetic
position coordinate. The contact interaction $g_{s}$ is a uniform interaction
in the $x$-$p_{y}$ plane as long as the position is conserved, and thus
infinite-ranged. On the other hand, the intertube interaction $g_{1}$ has a
factor of $\cos(y_{0}\frac{p_{1}-p_{2}}{2})\cos(y_{0}\frac{p_{3}-p_{4}}{2})$;
the interaction is strongest when $p_{1}\approx p_{2}$ and $p_{3}\approx
p_{4}$, which corresponds to the situation where the interacting particles are
close to each other. Therefore, the intertube interaction introduces a
shorter-ranged interaction along the $p_{y}$ direction compared with the
contact interaction. We note, however, that, considering $p_{y}$ as a
synthetic position, this intertube interaction is still an unconventional form
which cannot be derived from a two-body potential depending only on the
relative positions between two interacting particles.

Figure 4: The interacting ground states. The ground-state wave function
$|\Psi_{n}(x)|$ when $g_{s}=a_{\mathrm{osc}}\omega$ and $g_{1}=0$, plotted in
(a) the $x$-$y$ plane, and (b) the corresponding wave function in the
$x$-$p_{y}$ plane, $|\Psi(x,p_{y})|$. The ground-state wave function when
$g_{s}=a_{\mathrm{osc}}\omega$ and $g_{1}=2a_{\mathrm{osc}}\omega$, plotted in
(c) the $x$-$y$ plane, and (d) the corresponding wave function in the
$x$-$p_{y}$ plane, where dips of the wave function can be observed. (e) The
integrated density for each tube, $\int dx|\Psi_{n}(x)|^{2}$, plotted as a
function of $n$; a density modulation is visible. (f) The phase profile of the
wave function $\Psi(x,p_{y})$. The phase singularities (vortices) at $p_{y}>0$
are marked by blue circles. The units of length in the $x$ direction are
$a_{\mathrm{osc}}$, and the units of $p_{y}$ are $1/y_{0}$.

We first consider the ground state when only the contact interaction $g_{s}$
is present. In Figs. 4(a) and 4(b), we plot the ground state, obtained by
imaginary-time propagation of the Gross-Pitaevskii equation (6), when
$g_{s}=a_{\mathrm{osc}}\omega$ and $g_{1}=0$, assuming the normalization
$\sum_{n}\int dx|\Psi_{n}(x,t)|^{2}=1$. Figure 4(a) shows $|\Psi_{n}(x)|$ in
the $x$-$y$ plane, whereas Fig. 4(b) shows $|\Psi(x,p_{y})|$. We see that the
particles are clumped at the center of the $x$-$p_{y}$ plane, and no vortex-
like structure is seen. This is because, in order to obtain vortices in the
$x$-$p_{y}$ plane, we need a short-range repulsive interaction in $x$-$p_{y}$,
but the interaction due to $g_{s}$ is infinite-ranged in the $p_{y}$
direction.

To obtain vortices in the $x$-$p_{y}$ plane, we need to have a long-range
interaction along the $y$ direction, which introduces a shorter-range
interaction in the $x$-$p_{y}$ plane. Technology to realize long-range
interactions in synthetic quantum matter has been rapidly developing;
representative setups include atoms with dipole-dipole interactions
Lahaye:2009 , systems involving Rydberg states Saffman:2010 ; Gorshkov:2011 ;
Sevincli:2011 , and atoms in a cavity with light-mediated interactions
Ritsch:2013 . Now we include a long-range interaction in the simplest form
through the nearest-neighbor interaction $g_{1}$, which can be realized in
ultracold atomic gases using dipole-dipole interactions with proper
orientation of the dipole field Lahaye:2009 . In Figs. 4(c) and 4(d), we plot
the ground state when $g_{s}=a_{\mathrm{osc}}\omega$ and
$g_{1}=2a_{\mathrm{osc}}\omega$. We see that $|\Psi(x,p_{y})|$ has six dips.
By looking at its phase profile in Fig. 4(f), we confirm that these dips are
the vortices formed in the $x$-$p_{y}$ plane. In the $x$-$y$ plane, the state
has density peaks at every other tube, as seen in the integrated density for
each tube plotted in Fig. 4(e). Vortices on a cylinder are known to have
dynamical properties different from those on a plane Ho:2015 ; Guenther:2017 ;
our proposal provides a unique setup where such phenomena can be studied.
Since the Hamiltonian has a discrete translational symmetry along the angle
$\theta$ from the $x$ axis, the state breaks this translational symmetry, thus
forming a supersolid state in the $x$-$y$ plane. Whether this supersolid
structure persists in the absence of the confining potential $\tilde{\omega}$
in a thermodynamic limit is a question to be addressed in a future work.

We note a crucial difference between previous studies relating long-range
interactions, such as a dipole-dipole interaction, and the formation of
vortices in the presence of an artificial magnetic field Lahaye:2009 ;
Zhang:2016 . In these earlier works, the dipole-dipole interaction, an
artificial magnetic field, and vortices all appear in the same (physical) two-
dimensional plane. Our result, on the other hand, considers a long-range
(intertube) interaction in the $x$-$y$ plane, but the vortices appear in the
reciprocal $x$-$p_{y}$ plane where the artificial magnetic field is realized.

Below, we give a more detailed analysis of the vortices. We first show the
wave function profile when a hypothetical contact interaction in the
$x$-$p_{y}$ plane is present. We see that the vortices corresponds to a
supersolid structure in the $x$-$y$ plane, confirming that the vortices we
observed above are indeed created by the same mechanism as the ordinary
vortices in the presence of a magnetic field. We then discuss how vortices
form as one increases the nearest-neighbor interaction $g_{1}$.

###  IV.1 Vortices under a contact interaction in the $x$-$p_{y}$ plane

We first analyze the bosonic ground state if the interaction is contact in the
$x$-$p_{y}$ plane. Note that this is the situation one often encounters in the
$x$-$y$ plane in ultracold atomic gases in the presence of an artificial
magnetic field Madison:2000 ; Cooper:2008 ; Fetter:2009 , and one expects the
formation of vortices. We include the contact interaction in the $x$-$p_{y}$
plane by adding a term $\int dx\int dp_{y}\tilde{g}_{s}|\Psi(x,p_{y})|^{4}$ in
the energy functional, where $\Psi(x,p_{y})$ is the normalized wave function
in the $x$-$p_{y}$ plane. We choose
$\tilde{g}_{s}=10a_{\mathrm{osc}}y_{0}\omega$, $g_{s}=g_{1}=0$ and we take all
the other parameters to be the same as before; in particular, we take
$\tan\theta=1$. In Fig. 5, we plot the ground state obtained by imaginary time
evolution of the Gross-Pitaevskii equation. The absolute values of the wave
function in the $x$-$y$ as well as the $x$-$p_{y}$ planes are plotted in Figs.
5(a) and 5(b), respectively. We can see two density peaks in the $x$-$y$
plane, similar to the structure seen in Fig. 4 when the nearest-neighbor
interaction is included. The total density of each tube is plotted in Fig.
5(c), where the vertical axis denotes the integral of $|\Psi_{n}(x)|^{2}$ over
$x$. The corresponding wave function in the $x$-$p_{y}$ plane shows the vortex
structure, as confirmed by the phase profile in Fig. 5(d). These vortices,
formed by an artificial magnetic field and a contact interaction in the
$x$-$p_{y}$ plane, are the ordinary vortices which have been observed in
ultracold gases in the presence of an artificial magnetic field in real
($x$-$y$) space. Comparing these results with the wave function found in Fig.
4, we conclude that the vortex structure found in Fig. 4 has the same origin
as the conventional vortices in the presence of an artificial magnetic field.

Figure 5: The mean-field ground state when we consider a hypothetical contact
interaction in the $x$-$p_{y}$ plane. The units of length in the $x$ direction
are $a_{\mathrm{osc}}$, and the units of $p_{y}$ are $1/y_{0}$. We take the
strength of the contact interaction in the $x$-$p_{y}$ plane to be
$\tilde{g}_{s}=10a_{\mathrm{osc}}y_{0}\omega$, and no other interparticle
interaction exists. Other parameters are $\tan\theta=1$,
$\tilde{\omega}=\omega/100$, and $J_{y}=-0.01\omega$. (a) The absolute value
of the wave function in the $x$-$y$ plane, where two density peaks are seen.
(b) The wave function in the $x$-$p_{y}$ plane, where two dips of the wave
function are seen. (c) The integrated wave function of each tube, $\int
dx|\Psi_{n}(x)|^{2}$, from which one can see the density peaks more clearly.
(d) The phase profile of the wave function in the $x$-$p_{y}$ plane; the two
blue circles indicate the position of the phase singularities, namely the
vortices.

###  IV.2 Vortex formation

Now we come back to the interaction we originally considered: the contact
interaction $g_{s}$ and the nearest-neighbor interaction $g_{1}$. We analyze
how vortices form by varying the value of $g_{1}$ and looking at the wave
function profile. In Fig. 6, we plot, for three different values of $g_{1}$,
the wave function in the $x$-$y$ plane and the $x$-$p_{y}$ plane, as well as
the phase profile of the wave function in the $x$-$p_{y}$ plane and the
integrated density for each tube. We use $\tan\theta=1$,
$g_{s}=a_{\mathrm{osc}}\omega$, $\tilde{\omega}=\omega/100$, and
$J_{y}=-0.01\omega$. We can observe that vortices form as one increases
$g_{1}$. We see that, when $g_{1}=0.8a_{\mathrm{osc}}\omega$, there is no
phase singularity in the $x$-$p_{y}$ plane, and the integrated density for
each tube has a peak at $n=0$ and monotonically decreases as one goes away
from $n=0$. The situation changes for $g_{1}=1.0a_{\mathrm{osc}}\omega$.
Although it is still difficult to see vortex formation in $|\Psi(x,p_{y})|$,
by looking at the phase profile of $\Psi(x,p_{y})$, we clearly see six phase
singularities, indicating the formation of vortices. Correspondingly, density
modulation starts to be seen in the integrated density for each tube. When
$g_{1}=1.2a_{\mathrm{osc}}\omega$, the vortices become easily visible in
$|\Psi(x,p_{y})|$. The phase singularity and the density modulation also
become more pronounced. The vortices are fragile against increase of the
intertube hopping amplitude $J_{y}$. If we use $J_{y}=-0.02\omega$, the
density dip in the $x$-$y$ plane, and hence the vortex structure in the
$x$-$p_{y}$ plane, similar to the one found for the middle column of Fig. 6,
will not appear until we increase $g_{1}$ until around $g_{1}\approx
1.3a_{\mathrm{osc}}\omega$. Thus, a stronger interaction is necessary to form
vortices. This tendency can be understood in the following way: Intertube
hopping $J_{y}$ introduces a sinusoidal potential energy in the $x$-$p_{y}$
plane along the $p_{y}$ direction. To create vortices in the $x$-$p_{y}$
plane, the interaction effect should be strong enough to overcome the
sinusoidal potential created by $J_{y}$. We note that the vortices can in
principle form also when there is no intertube hopping, $J_{y}=0$. In such a
case, however, there is no reason for the entire system to develop a phase
coherence, and thus the resulting structure in the $x$-$y$ plane is not a
supersolid anymore; it would rather be a phase incoherent density wave.

Figure 6: Formation of vortices in the $x$-$p_{y}$ plane. The three columns
represent the ground states for $g_{1}=0.8a_{\mathrm{osc}}\omega$,
$g_{1}=1.0a_{\mathrm{osc}}\omega$, and $g_{1}=1.2a_{\mathrm{osc}}\omega$,
respectively, from left to right. In each column, the top row shows
$|\Psi_{n}(x)|$ of the ground state, the second row shows $|\Psi(x,p_{y})|$,
the third row shows the phase profile of $\Psi(x,p_{y})$, and the bottom row
shows $\int dx|\Psi_{n}(x)|^{2}$ as a function of $n$. In the phase profile
panels, we marked phase singularities at $p_{y}>0$ by blue circles. (There is
also the same number of phase singularities at $p_{y}<0$, which are not
encircled.) The units of length in the $x$ direction are $a_{\mathrm{osc}}$,
and the units of $p_{y}$ are $1/y_{0}$.

##  V Experimental platforms

Our proposal is relevant to various platforms of synthetic quantum matter. In
ultracold atomic gases, decoupled arrays of one-dimensional tubes can be
realized with a two-dimensional optical lattice Bloch:NP . Our proposal can
then be realized by superposing another harmonic trapping potential with an
appropriate angle. The quantized Hall response and many-body physics are
directly relevant in ultracold atomic gases Dauphin:2013 ; Aidelsburger:2015 ;
Price:2016 . Various photonic platforms are also relevant to our proposal. For
example, exciton-polariton microcavity can be made into an array of one-
dimensional wires Boulier:2020 ; Amo:rv1 , and an additional harmonic
potential can be applied by varying the width of each wire. Since it is
extremely challenging to break the physical time-reversal symmetry in exciton-
polaritons to open a topological gap Klembt:2018 , our proposal provides a
valuable alternative to explore Chern insulator physics in exciton-polariton
microcavities. Photonic cavity arrays and propagating waveguide geometry are
also suited to realize our proposal, where harmonic potentials can be applied
by varying the cavity or waveguide size along the $x$ direction Hafezi:2013 ;
Rechtsman:2013 ; for these systems it is more natural to consider a lattice
structure in the $x$ direction, which results in a coupled wire setup in the
$x$-$p_{y}$ plane Ozawa:2017 . Edge-related physics are more easily accessible
in such photon-related platforms Hafezi:2013 ; Rechtsman:2013 ; Milicevic:2017
.

We note that our proposal depends crucially on the fact that properties in
momentum space ($p_{y}$ direction) are observable in this synthetic quantum
matter. In ultracold gases, the momentum-space wave function is accessible
through time-of-flight imaging and, in photon-related systems, far-field
imaging provides information on momentum space.

##  VI Conclusion

We have shown that a set of tubes with harmonic trapping potentials whose
centers are shifted is equivalent to charged particles in a magnetic field in
the $x$-$p_{y}$ plane. The underlying idea of our proposal is to replace a
momentum in one direction by position to realize a simpler Hamiltonian, which
simulates the original Hamiltonian in a hybrid real-momentum space. This
method is quite versatile; for example, with a similar mechanism, a collection
of one-dimensional Aubry-André model Aubry:1980  (the Harper model) with
shifted modulation phases simulates a two-dimensional lattice with a magnetic
field, the Harper-Hofstadter model Harper:1955 ; Hofstadter:1976 , in real-
momentum space, without breaking the physical time-reversal symmetry. Our
proposal provides a simple and powerful method applicable to realize various
Hamiltonians, opening a possibility for quantum simulation in synthetic
quantum matter.

###### Acknowledgements.

I am grateful to Bryce Gadway, the discussion with whom led to the initial
idea of this paper. I also thank Iacopo Carusotto, Jean Dalibard, Nathan
Goldman, Maciej Lewenstein, and Hannah M. Price for useful discussions. I am
supported by JSPS KAKENHI Grant Number JP20H01845, JST PRESTO Grant Number
JPMJPR19L2, JST CREST Grant Number JPMJCR19T1, RIKEN Incentive Research
Project, and the Interdisciplinary Theoretical and Mathematical Sciences
Program (iTHEMS) at RIKEN.

## References

  * (1) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (2) I. Bloch, J. Dalibard, S. Nascimbène , Quantum simulations with ultracold quantum gases, Nature Physics 8, 267-276 (2012).
  * (3) R. Blatt, C. F. Roos, Quantum simulations with trapped ions, Nature Physics 8, 277-284 (2012).
  * (4) A. Aspuru-Guzik, P. Walther, Photonic quantum simulators, Nature Physics 8, 285-291 (2012).
  * (5) I. Carusotto and C. Ciuti, Quantum fluids of light, Rev. Mod. Phys. 85, 299 (2013).
  * (6) I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).
  * (7) T. Boulier, M. J. Jacquet, A. Maître, G. Lerario, F. Claude, S. Pigeon, Q. Glorieux, A. Bramati, E. Giacobino, A. Amo, J. Bloch, Microcavity Polaritons for Quantum simulation, arXiv:2005.12569 (2020).
  * (8) M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * (9) X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (10) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * (11) N.P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * (12) N. R. Cooper, J. Dalibard, and I. B. Spielman, Topological bands for ultracold atoms, Rev. Mod. Phys. 91, 015005 (2019).
  * (13) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 015006 (2019).
  * (14) 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).
  * (15) D. C. Tsui, H. L. Stormer, and A. C. Gossard, Two-Dimensional Magnetotransport in the Extreme Quantum Limit, Phys. Rev. Lett. 48, 1559 (1982).
  * (16) D. Yoshioka, The Quantum Hall Effect (Springer, Berlin, 2002).
  * (17) M. Aidelsburger, S. Nascimbene, and N. Goldman, Artificial gauge fields in materials and engineered systems, Comptes Rendus Physique 19, 394 (2018).
  * (18) K. W. Madison, F. Chevy, W. Wohlleben, and J. Dalibard, Vortex Formation in a Stirred Bose-Einstein Condensate, Phys. Rev. Lett. 84, 806 (2000).
  * (19) N. R. Cooper, Rapidly rotating atomic gases, Advances in Physics 57, 539 (2008).
  * (20) A. L. Fetter, Rotating trapped Bose-Einstein condensates, Rev. Mod. Phys. 81, 647 (2009).
  * (21) N. Goldman and J. Dalibard, Periodically Driven Quantum Systems: Effective Hamiltonians and Engineered Gauge Fields, Phys. Rev. X 4, 031027 (2014).
  * (22) N. Goldman, G. Juzeliūnas, P. Öhberg, and I. B. Spielman, Light-induced gauge fields for ultracold atoms, Rep. Prog. Phys. 77, 126401 (2014).
  * (23) T. Ozawa and H. M. Price, Topological quantum matter in synthetic dimensions, Nature Reviews Physics 1, 349-357 (2019).
  * (24) T.-L. Ho and V. B. Shenoy, Local Spin-Gauge Symmetry of the Bose-Einstein Condensates in Atomic Gases, Phys. Rev. Lett. 77, 2595 (1996).
  * (25) V. Pietilä and Mikko Möttönen, Creation of Dirac Monopoles in Spinor Bose-Einstein Condensates, Phys. Rev. Lett. 103, 030401 (2009).
  * (26) J.-y. Choi, S. Kang, S. W. Seo, W. J. Kwon, and Y.-i. Shin, Observation of a Geometric Hall Effect in a Spinor Bose-Einstein Condensate with a Skyrmion Spin Texture, Phys. Rev. Lett. 111, 245301 (2013).
  * (27) M. W. Ray, E. Ruokokoski, S. Kandel, M. Möttönen, and D. S. Hall, Observation of Dirac monopoles in a synthetic magnetic field, Nature 505, 657 (2014).
  * (28) M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Photonic Floquet topological insulators, Nature 496, 196 (2013).
  * (29) L. Yuan, Q. Lin, M. Xiao, and S. Fan, Synthetic dimension in photonics, Optica 5, 1396-1405 (2018).
  * (30) E. Lustig, S. Weimann, Y. Plotnik, Y. Lumer, M. A. Bandres, A. Szameit, and M. Segev, Photonic topological insulator in synthetic dimensions, Nature 567 356-360(2019).
  * (31) M. Hafezi, S. Mittal, J. Fan, A. Migdall, and J. M. Taylor, Imaging topological edge states in silicon photonics, Nature Photonics 7, 1001 (2013).
  * (32) N. Schine, A. Ryou, A. Gromov, A. Sommer, and J. Simon, Synthetic Landau levels for photons, Nature 534, 671-675 (2016).
  * (33) L. D. Landau and E. M. Lifshitz, Quantum Mechanics, Third edition, (Butterworth-Heinemann, Oxford, 1977).
  * (34) I. Spielman, Raman processes and effective gauge potentials, Phys. Rev. A 79, 063613 (2009).
  * (35) Y. J. Lin, R. L. Compton, K. Jiménez-García, J. V. Porto, and I. B. Spielman, Synthetic magnetic fields for ultracold neutral atoms, Nature 462, 628 (2009).
  * (36) D. J. Thouless, Quantization of particle transport, Phys. Rev. B 27, 6083 (1983).
  * (37) S. Nakajima, T. Tomita, S. Taie, T. Ichinose, H. Ozawa, L. Wang, M. Troyer, Y. Takahashi, Topological Thouless pumping of ultracold fermions, Nature Physics 12, 296-300 (2016).
  * (38) M. Lohse, C. Schweizer, O. Zilberberg, M. Aidelsburger, and I. Bloch, A Thouless quantum pump with ultracold bosonic atoms in an optical superlattice, Nature Physics 12, 350-354 (2016).
  * (39) R. Laughlin, Quantized Hall conductivity in two dimensions, Phys. Rev. B 23, 5632 (1981).
  * (40) A. Dauphin and N. Goldman, Extracting the Chern Number from the Dynamics of a Fermi Gas: Implementing a Quantum Hall Bar for Cold Atoms, Phys. Rev. Lett. 111, 135302 (2013).
  * (41) M. Aidelsburger, M. Lohse, C. Schweizer, M. Atala, J. T. Barreiro, S. Nascimbène, N. R. Cooper, I. Bloch, and N. Goldman, Measuring the Chern number of Hofstadter bands with ultracold bosonic atoms, Nature Physics 11, 162-166 (2015).
  * (42) H. M. Price, O. Zilberberg, T. Ozawa, I. Carusotto, and N. Goldman, Measurement of Chern numbers through center-of-mass responses, Phys. Rev. B 93, 245113 (2016).
  * (43) T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T Pfau, The physics of dipolar bosonic quantum gases, Rep. Prog. Phys. 72, 126401 (2009).
  * (44) M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with Rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * (45) A. V. Gorshkov, J. Otterbach, M. Fleischhauer, T. Pohl, and M. D. Lukin, Photon-Photon Interactions via Rydberg Blockade, Phys. Rev. Lett. 107, 133602 (2011).
  * (46) S. Sevinçli, N. Henkel, C. Ates, and T. Pohl, Nonlocal Nonlinear Optics in Cold Rydberg Gases, Phys. Rev. Lett. 107, 153001 (2011).
  * (47) H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Cold atoms in cavity-generated dynamical optical potentials, Rev. Mod. Phys. 85, 553 (2013).
  * (48) T.-L. Ho and B. Huang, Spinor Condensates on a Cylindrical Surface in Synthetic Gauge Fields, Phys. Rev. Lett. 115, 155304 (2015).
  * (49) N.-E. Guenther, P. Massignan, and A. L. Fetter, Quantized superfluid vortex dynamics on cylindrical surfaces and planar annuli, Phys. Rev. A 96, 063608 (2017).
  * (50) X.-F. Zhang, L. Wen, C.-Q. Dai, R.-F. Dong, H.-F. Jiang, H. Chang, and S.-G. Zhang, Exotic vortex lattices in a rotating binary dipolar Bose-Einstein condensate, Sci. Rep. 6, 19380 (2016).
  * (51) I. Bloch, Ultracold quantum gases in optical lattices, Nat. Phys. 1, 23 (2005).
  * (52) A. Amo and J. Bloch, Exciton-polaritons in lattices: A non-linear photonic simulator, Comp. Rend. Phys. 17, 805 (2016).
  * (53) S. Klembt, T. H. Harder, O. A. Egorov, K. Winkler, R. Ge, M. A. Bandres, M. Emmerling, L. Worschech, T. C. H. Liew, M. Segev, C. Schneider, and S. Höfling, Exciton-polariton topological insulator, Nature 562, 552 (2018).
  * (54) T. Ozawa and I. Carusotto, Synthetic Dimensions with Magnetic Fields and Local Interactions in Photonic Lattices, Phys. Rev. Lett. 118, 013601 (2017).
  * (55) M. Milićević, T. Ozawa, G. Montambaux, I. Carusotto, E. Galopin, A. Lemaître, L. Le Gratiet, I. Sagnes, J. Bloch, and A. Amo, Orbital edge states in a photonic honeycomb lattice, Phys. Rev. Lett. 118, 107403 (2017).
  * (56) S Aubry, G André, Analyticity breaking and Anderson localization in incommensurate lattices, Ann. Israel Phys. Soc 3, 133 (1980).
  * (57) P. G. Harper, Single Band Motion of Conduction Electrons in a Uniform Magnetic Field, Proc. Phys. Soc. Lond. A 68, 874 (1955).
  * (58) D. R. Hofstadter, Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields, Phys. Rev. B 14, 2239 (1976).


# A regular interior solution of Einstein field equations.

Gabino Estevez-Delgado
Facultad de Químico Farmacobiología de la UMSNH
Tzintzuntzan No. 173, Col. Matamoros,
Morelia Michoacán, C.P. 58240, México<EMAIL_ADDRESS>
Joaquin Estevez-Delgado
Facultad de Ciencias Físico Matemáticas de la UMSNH
Edificio B, Ciudad Universitaria,
Morelia Michoacán, C.P. 58040, México<EMAIL_ADDRESS>
Modesto Pineda Duran
Instituto Tecnológico Superior de Tacámbaro,
Av. Tecnológico No 201 Zona el Gigante, C. P. 61650,
Tacambaro Michoacán, México<EMAIL_ADDRESS>
and
Arthur Cleary-Balderas
Facultad de Ingeniería Eléctrica de la UMSNH
Edificio $\Omega$, Ciudad Universitaria, C.P. 58030,
Morelia Michoacán, México<EMAIL_ADDRESS>
###### Abstract

Starting from the solution of the Einstein field equations in a static and
spherically symmetric spacetime which contains an isotropic fluid, we
construct a model to represent the interior of compact objects with
compactness rate $u=\frac{GM}{c^{2}R}<0.23577$. The solution is obtained by
imposing the isotropy condition for the radial and tangential pressures, this
generates an ordinary differential equation of second order for the temporal
$g_{tt}$ and radial $g_{rr}$ metric potentials, which can be solved for a
specific function of $g_{tt}$. The graphic analysis of the solution shows that
it is physically acceptable, that is to say, the density, pressure and speed
of sound are positive, regular and monotonically decreasing functions, also,
the solution is stable due to meeting the criteria of the adiabatic index.
When taking the data of mass $M=1.44^{+0.15}_{-0.14}M_{\odot}$ and radius
$R=13.02^{+1.24}_{-1.06}km$ which corresponds to the estimations of the star
PSR J0030+045 we obtain values of central density $\rho_{c}=7.5125\times
10^{17}kg/m^{3}$ for the maximum compactness $u=0.19628$ and of
$\rho_{c}=2.8411\times 10^{17}kg/m^{3}$ for the minimum compactness
$u=0.13460$, which are consistent with those expected for this type of stars.

##  1 Introduction

Describing the interior of the stars and determining their average composition
requires many different complementary approaches as are: chemistry,
thermodynamics, nuclear physics, particle physics and gravitational physics.
And in the case that there are no instabilities generated, once all the
nuclear fuel has been used, the star shrinks and, depending of the mass and
the stability present, it may form a white dwarf, a neutron star or a quark
star [1, 2]. Of course the description of each one of these stages and their
stability is a far more delicate matter which involves a detailed analysis
that considers the type of predominant particles in the interior of the star,
whether they are electrons, neutrons, or quarks, even when the stars are
hybrids. According to the focus of this job, the matter in general is supposed
to be described in a satisfactory manner by a perfect fluid and it will not be
necessary to give a specific shape of a state equation. We know that,
depending on the value of the mass and radius of a star, it may be a white
dwarf, a neutron star or a quark star, this will also determine the orders of
magnitude of the density, for example densities in the order of
$10^{18}kg/m^{2}$ are typical for neutron stars. As such in this situation,
given the high density, it results adequate to describe the interior of the
stars by means of Einstein’s general relativity theory. The interior solutions
have been approached for over a century, the first of these was constructed
for a static and spherically symmetric spacetime with matter from a perfect
fluid and incompressible density, known as the interior Schwarzschild
solution. Although to start with, the density being constant is an unrealistic
requirement, its consequences revealed some differences with the treatment of
stellar models in the context of Newton’s theory of gravitation. One of these
is the compactness relation $u=GM/c^{2}R<4/9$, where $M$ is the mass and $R$
is the radius, this indicates that it is not possible to have stars with
arbitrary mass and radius. Afterwards it was shown that this relation is not
exclusive of this idealized model and that it is present for stars with a
monotonically decreasing density function for which the exterior geometry is
given by the Schwarzschild solution [3, 4]. And although more that 130
interior solutions with perfect fluid have been published, only a few met the
characteristics that makes them physically acceptable. From an analysis done
in 1998 on a total of 127 reported solutions, only 16 of these had their
density and pressure functions be positive, regular and monotonically
decreasing functions and also had a speed of sound that didn’t violate the
causality condition. And only 9 out of these 16 had a speed of sound which
decreases monotonically with the radius [5, 6, 7, 8, 9, 10, 11, 12, 13, 14,
15, 16]. Although this last requirement is debatable since, for example, for
the realistic MIT Bag state equation $P(\rho)=\frac{1}{3}(c^{2}\rho-4B_{g})$
associated to quark stars, the speed of sound is $v_{s}^{2}=\frac{c^{2}}{3}$,
which is not a monotonically decreasing function. The construction of stellar
solutions with perfect fluid is an active field although the difficulty in
obtaining physically acceptable solutions limits a great number of these. A
point which has been explored in relation to these is employing isotropic
coordinates [17, 18, 19, 20, 21] which has favored the integration of the
equations system and its application in the description of stars like HerX1,
4U1538-52, LMC X - 4, SAX J1808.4-3658 [22]. A particular class of solutions
constructed in Schwarzschild’s coordinates assumes a metric potential
$g_{tt}\\!=\\!-(1+ar^{2})^{n}$ to this group belong the Tolman IV [6] and
Durgapal [15] solutions, extensions for other values of $n$ a positive integer
or negative fractional value [23, 24, 25] have been done, showing that for
$n\geq 4$ the solutions that are generated are physically acceptable. Other
recent works have addressed the possibility of generating exact solutions with
metric potential $g_{tt}\\!=\\!-\frac{1+ar^{2}}{1+b^{2}}$ showing that this
functional form is adequate in obtaining physically acceptable solutions and
it’s consistent with the stars Her X-1 [26, 27], PSRJ0348+0432 [28], PSR
B0943+10 [29], PSR J0737 -3039A [30] and PSRJ1614 2230 [31]. Motivated by
these few last investigation works, in this report we present a new solution
to Einstein’s equations with perfect fluid with a metric function $g_{tt}$ and
its application to the star PSR J0030+045. The structure of this article will
be as follows: in the section 2 we present Einstein field equations for a
static and spherically symmetrical spacetime with a perfect fluid and assuming
the form of the metric function
$g_{tt}\\!=\\!-{{A^{2}\left(1+a{r}^{2}\right)^{2}}/\\!\left[{{1\\!+\\!(\frac{3}{\sqrt{2}}-1)a{r}^{2}}}\right]}$
we obtain the solution from the isotropy equation. In section 3 we determine
the hydrostatic functions and impose the coupling conditions between the
interior and the exterior solutions to determine the integration constants. In
section 4 graphic analysis of the solution is done, taking the observational
values of the mass and radius of the star PSRJ0030+045 we determine the
physical values of the pressure, density and speed of sound, starting from the
graphic analysis and from the data, we show that the solution is physically
acceptable. The conclusions and discussion of future works are presented in
the section 5.

##  2 The field equations and the solution

The interior geometry of a static and spherically symmetric spacetime can be
described through a line element [32, 33] :

$\displaystyle
ds^{2}=-y(r)^{2}dt^{2}+\frac{dr^{2}}{B(r)}+r^{2}(d\theta^{2}+\sin^{2}{\theta}d\phi^{2}),$
(1)

where $y,B$ are functions of the radial coordinate $r\leq R$. Einstein
equations $G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}=kT_{\mu\nu}$, (with
$R_{\mu\nu}$, $R$ and $g_{\mu\nu}$ the components of the Ricci tensor, the
Ricci scalar and the metric tensor respectively), have as source the
distribution of matter from a perfect fluid described by the energy-momentum
tensor:

$T_{\mu\nu}=c^{2}\rho u_{\mu}u_{\nu}+{\it P}(u_{\mu}u_{\nu}+g_{\mu\nu}),$ (2)

with $u^{\mu}$ the four velocity of the fluid, $\rho$ the energy density and
$P$ the pressure. The non zero components of the Einstein equations are:

$\displaystyle kc^{2}\rho$ $\displaystyle=$
$\displaystyle-\frac{B^{\prime}}{r}+\frac{1-B}{r^{2}},$ (3) $\displaystyle
k{\it P}$ $\displaystyle=$
$\displaystyle\frac{2By^{\prime}}{ry}-\frac{1-B}{{r}^{2}},$ (4) $\displaystyle
k{\it P}$ $\displaystyle=$
$\displaystyle\frac{(ry^{\prime\prime}+y^{\prime})B}{ry}+\frac{(ry^{\prime}+y)B^{\prime}}{2ry},\qquad$
(5)

with the derivative in relation to the radial coordinate $r$ denoted by ′.
Meanwhile the equation of conservation for the energy-momentum tensor
$\nabla_{\mu}T^{\mu}\,_{\nu}=0$ implies the Tolman-Oppenheimer-Volkov (TOV)
equation [32, 33]:

${\it P}^{\prime}=-{\frac{\left({\it P}+c^{2}\rho\right)y^{\prime}}{y}}.$ (6)

Although this last one is not an independent equation, since it can be
obtained from the system of equations (3) - (5). Being this the set of
equations for which we will obtain the solution starting from a function
$y(r)$.

###  2.1 The solution

For the integration of the system we propose a metric function
$g_{tt}=-y(r)^{2}$ with the form of $y(r)$ similar, but slightly different, to
one employed previously with which it was possible to integrate the system in
an adequate manner and which resulted useful to describe compact objects with
a compactness rate $u=GM/c^{2}R\leq 0.2660858316$ [28] . Specifically we have
that

$y\left(r\right)={\frac{A\left(1+a{r}^{2}\right)}{\sqrt{1+\left(\frac{3}{\sqrt{2}}-1\right)a{r}^{2}}}},$
(7)

where $A$ and $a$ are constants. One useful relation for the integration of
the system is obtained by replacing $y(r)$ in the difference of the equations
(4) and (5), leading to:

$B^{\prime}-{\frac{\sqrt{8}[3+\sqrt{2}+4(\sqrt{8}-1)a{r}^{2}+7{a}^{2}{r}^{4}]B}{(1+\sqrt{2}a{r}^{2})\left(2+3\sqrt{2}+7a{r}^{2}\right)r}}+{\frac{(1+a{r}^{2})(2+3\sqrt{2}+7a{r}^{2})\sqrt{2}}{(1+\sqrt{2}a{r}^{2})(3+\sqrt{2}+7a{r}^{2})r}}=0,$

after the integration we obtain:

$\displaystyle B\left(r\right)$ $\displaystyle=$ $\displaystyle
1+{\frac{\left(31-22\sqrt{2}\right)\left(3\,\sqrt{2}+2+7\,a{r}^{2}\right)^{3}a{r}^{2}}{343\,\left(1+a\sqrt{2}{r}^{2}\right)^{3}}\left[C+\ln\left({\frac{\sqrt{2}+3+7\,a{r}^{2}}{3\,\sqrt{2}+2+7\,a{r}^{2}}}\right)\right]}\qquad$
(8) $\displaystyle+$
$\displaystyle{\frac{\left(22-17\sqrt{2}\right)\left(732\sqrt{2}+1832+7\left(555\sqrt{2}+251\right)a{r}^{2}+4606{a}^{2}{r}^{4}\right)a{r}^{2}}{2303\left(\sqrt{2}+2a{r}^{2}\right)^{3}}},$

$C$ is the constant of integration.

##  3 Hydrostatic functions and physical conditions

Once we know the metric functions we will proceed to determine the hydrostatic
functions. Replacing the functions $y(r)$ and $B$ given by the equations (7)
and (8) in the equations (3) and (4) we determine the density and the
pressure:

$\displaystyle k{c}^{2}\rho\left(r\right)$ $\displaystyle=$
$\displaystyle{\frac{3\left(2\,\sqrt{2}+6+\left(-4+15\,\sqrt{2}\right)a{r}^{2}+14\,{a}^{2}{r}^{4}\right)\left(1-B\left(r\right)\right)}{\left(3\,\sqrt{2}+2+7\,a{r}^{2}\right)\left(\sqrt{2}+2\,a{r}^{2}\right){r}^{2}}}$
(9)
$\displaystyle-{\frac{14\left(12\,\sqrt{2}-13+7\,a{r}^{2}\right){a}^{2}{r}^{2}}{\left(3\,\sqrt{2}+2+7\,a{r}^{2}\right)\left(\sqrt{2}+2\,a{r}^{2}\right)\left(\sqrt{2}+3+7\,a{r}^{2}\right)}},$
$\displaystyle kP(r)$ $\displaystyle=$
$\displaystyle{\frac{2\left(6\,\sqrt{2}-3+7\,a{r}^{2}\right)a}{\left(1+a{r}^{2}\right)\left(3\,\sqrt{2}+2+7\,a{r}^{2}\right)}}$
(10)
$\displaystyle-{\frac{\left(3\,\sqrt{2}+2+3\,\left(5\,\sqrt{2}+1\right)a{r}^{2}+21\,{a}^{2}{r}^{4}\right)\left(1-B\left(r\right)\right)}{\left(1+a{r}^{2}\right)\left(3\,\sqrt{2}+2+7\,a{r}^{2}\right){r}^{2}}}.$

In these equations the expression $(1-B)/r^{2}$ appears, however, it is
regular when $r=0$ as it can be seen from the equation (8). Another important
relation to determine if the solution is physically acceptable is the speed of
sound, since it is required that the speed of sound in the model does not
violate the causality condition. In this case by means of the chain rule we
obtain the speed of sound:

$\frac{v^{2}(r)}{c^{2}}=\frac{1}{c^{2}}\frac{\partial
P(\rho)}{\partial\rho}={\frac{\left(S_{{1}}\left(r\right)B\left(r\right)+\left(3\,\sqrt{2}+2+7\,a{r}^{2}\right)^{2}\left(1+a{r}^{2}\right)^{2}\right)S_{{2}}\left(r\right)}{\left(S_{{3}}\left(r\right)\left(\sqrt{2}+3+7\,a{r}^{2}\right)^{2}B\left(r\right)+S_{{4}}\left(r\right)\right)\left(1+a{r}^{2}\right)^{2}}},$

where

$\displaystyle S_{{4}}(r)$ $\displaystyle=$
$\displaystyle(1+a{r}^{2})\left[15\sqrt{2}+150+7(34+9\sqrt{2})a{r}^{2}+98{a}^{2}{r}^{4}\right]\left[3\sqrt{2}+2+7a{r}^{2}\right]^{2}\\!,$
$\displaystyle S_{{1}}(r)$ $\displaystyle=$
$\displaystyle\left(\sqrt{2}-4\right)\left(\sqrt{2}+3+7\,a{r}^{2}\right)\left(2+\sqrt{2}+\left(5\,\sqrt{2}-1\right)a{r}^{2}+3\,{a}^{2}{r}^{4}\right),$
$\displaystyle S_{{3}}(r)$ $\displaystyle=$ $\displaystyle
3\,\left(\sqrt{2}-4\right)\left(10\,\sqrt{2}+30+3\,\left(-4+15\,\sqrt{2}\right)a{r}^{2}+14\,{a}^{2}{r}^{4}\right),$
$\displaystyle S_{{2}}(r)$ $\displaystyle=$
$\displaystyle\left(6\,\sqrt{2}-3+7\,a{r}^{2}\right)\left(\sqrt{2}+3+7\,a{r}^{2}\right)\left(\sqrt{2}+2\,a{r}^{2}\right).$

###  3.1 Criteria for physical acceptability

Obtaining a solution to Einstein’s equations is not a guarantee that said
solution is physically acceptable, there are many solutions that are not
physically acceptable [5] due to the fact that they do not comply with certain
properties. In the following we will mention the requirements that must be
satisfied, some of these will be applied directly and others will be shown in
a graphical manner in the following section [28].

One of the conditions that must be met, is the regularity of the geometry when
approaching the center. Which can be expressed in an algebraic manner, through
the Kretschmann scalar, given its extension, it is enough with showing that
the metric coefficients around $r=0$ are of the form $\alpha+\beta
r^{2}+O(r^{4})$. The expansion of $B(r)$ and $y(r)$ in the proximity of $r=0$
gives us:

$y\left(r\right)=A\left((1-{\frac{\left(3\,\sqrt{2}-6\right)a}{4}}{r}^{2}-{\frac{\left(30\,\sqrt{2}-41\right){a}^{2}}{16}}{r}^{4}+O\left({r}^{6}\right))\right),$
$B\left(r\right)=1+{\frac{2\,a\left(\left(17\,\sqrt{2}-26\right)\left(C+\ln\left({\frac{\sqrt{2}+3}{3\,\sqrt{2}+2}}\right)\right)+41\,\sqrt{2}-80\right){r}^{2}}{49}}+O\left({r}^{4}\right),$

besides the regularity, the geometry must be absent of any event horizon, this
property is easier to demonstrate through a graphic analysis and it will be
analysed in the following section.
The density and pressure must be finite, positive and monotonically decreasing
as functions of the radial coordinate. That is to say, for $r\in(0,R)$,
$\rho^{\prime}<0$ and $P^{\prime}<0$ (condition that will be analysed
graphically) and in the center they must have their maximum value. which
implies the following set of inequalities:

$k{c}^{2}\rho(0)={\frac{6\,a\left(\left(132\,\sqrt{2}+193\right)\left(2\,C-\ln\left(2\right)\right)+3006\,\sqrt{2}+4286\right)}{6713\,\sqrt{2}+9506}}>0,$
(11)
$kP(0)=1/49\,a\left(\left(17\,\sqrt{2}-26\right)\left(2\,C-\ln\left(2\right)\right)-65\,\sqrt{2}+134\right)>0,$
(12)

$\rho^{\prime\prime}(0)=-{\frac{5{a}^{2}\left[6\left(195-103\sqrt{2}\right)\left(2C-\ln
2\right)+3819+6257\sqrt{2}\right]}{49\,\left(\sqrt{2}+3\right)^{3}}}<0,$ (13)
$\\!\\!P^{\prime\prime}(0)=-{\frac{3{a}^{2}\left[2\left(113-72\,\sqrt{2}\right)\left(2\,C-\ln
2\right)+2759-1248\,\sqrt{2}\right]}{49\,\left(\sqrt{2}+3\right)^{2}}}<0.$
(14)

In addition to these inequalities the solution satisfies $\rho^{\prime}(0)=$
and $P^{\prime}(0)=0$, that together with the inequalities (11)-(14) implies
that $r=0$ is a maximum for the functions $\rho$ and $P$. The causality
condition in the center of the star requires that it satisfies

$\\!\\!0\leq\frac{v(0)^{2}}{c^{2}}={\frac{3[24\,\sqrt{2}+4\,C-2\,\ln\,2+\,55\,]}{5[96\sqrt{2}+12C-6\ln
2+121]}}\leq 1.$ (15)

Combining the previous equations we can determine inequalities for the
constants $(C,a)$, in particular forming
$k\left(\rho\left(0\right){c}^{2}+3\,\Pr\left(0\right)\right)=9\,\left(2-\sqrt{2}\right)a>0$,
from where we obtain that $a>0$.
The constants $C$ and $W$ which appear in the metric functions are determined
by imposing that the interior and exterior geometry on the surface of the star
$r=R$ are joined in a continual manner and that the pressure is zero on the
surface. The exterior geometry is described by the exterior Schwarzschild
solution:

$\displaystyle ds^{2}\\!$ $\displaystyle\\!=\\!$
$\displaystyle-\left(1-\frac{2GM}{c^{2}r}\right)\,dt^{2}+\left(1-\frac{2GM}{c^{2}r}\right)^{-1}dr^{2}+\,r^{2}(d\theta^{2}+\sin^{2}{\theta}\,d\phi^{2}),\quad
r\geq R,$

where $M$ represents the total mass inside the fluid sphere. When we impose
$P(R)=0$, from the equation (10) we obtain $C$:

$\displaystyle C$ $\displaystyle=$
$\displaystyle-\ln\left({\frac{\sqrt{2}+3+7\,w}{3\,\sqrt{2}+2+7\,w}}\right)+{\frac{\left(25\,\sqrt{2}+47\right)W_{1}}{W_{2}\left(3\,\sqrt{2}+2+7\,w\right)^{2}}},$
(16)

where $w=aR^{2}$,
$W_{2}=822\,\sqrt{2}+548+822\,\left(5\,\sqrt{2}+1\right)w+5754\,{w}^{2}$ and

$W_{1}=3286\sqrt{2}+4048+(8149\sqrt{2}+18623)w+7(2757\sqrt{2}+1075){w}^{2}+13426{w}^{3}.$

Meanwhile from the continuity of the component $g_{tt}$ in $r=R$ it results:

${A}^{2}={\frac{\left(3\,\sqrt{2}-2\right)\left(3\,\sqrt{2}+2+7\,a{R}^{2}\right)^{2}}{14\left[3\,\sqrt{2}+2+3\,\left(5\,\sqrt{2}+1\right)a{R}^{2}+21\,{a}^{2}{R}^{4}\right]\left(1+a{R}^{2}\right)}}.$
(17)

The continuity of $g_{rr}$ in $r=R$ determines the value of the compactness as
function of $w$:

$u(w)=\frac{GM}{c^{2}R}=\frac{1}{2}(1-B(R))={\frac{\left(6\,\sqrt{2}-3+7\,w\right)w}{3\,\sqrt{2}+2+3\,\left(5\,\sqrt{2}+1\right)w+21\,{w}^{2}}}.$
(18)

The rest of the conditions require of a graphical analysis and these
correspond to the Energy conditions:

\- The Strong Energy Condition: $c^{2}\rho+3P\geq 0$, $c^{2}\rho+P\geq 0$ or

\- The Dominant Energy Condition: $\rho\geq 0$ and $c^{2}\rho\geq|P|$
And the Stability condition, a configuration of static and spherically
symmetrical matter is stable if it satisfies the relativistic condition for
the adiabatic index:

$\Gamma=\frac{P+c^{2}\rho}{c^{2}P}\frac{dP}{d\rho}>\frac{4}{3}\qquad\forall\;r\in[0,R]$

##  4 Graphic representation of the solution

From the graphic analysis of the functions of density, pressure, speed of
sound and adiabatic index we obtain that the function which restricts the
values of the parameter $w\leq w_{0}=0.90378$ is that of the adiabatic index,
specifically for values of $w>w_{0}$ the adiabatic index $\gamma(0)<4/3$,
which implies that the solution will be unstable. This maximum value $w_{0}$
through the equation (18) allows us to obtain the maximum permissible
compactness value for the compactness $u\leq u_{0}=0.23577$. Although the
solution is physically acceptable for the compactness values $u\leq u_{0}$ in
the graphic analysis we will focus on the particular case of the star PSR
J0030+0451 with estimates of mass $M=1.44^{+0.15}_{-0.14}M_{\odot}$ and radius
$R=13.02^{+1.24}_{-1.06}km$, obtained through the study of the X-ray emission
by means of the NICER (Neutron star Interior Composition Explorer) telescope
from the international space station [34]. The graphic representation will be
done in terms of the dimensionless variable $x=r/R$ and the dimensionless
functions associated to the physical quantities of density $kc^{2}R^{2}\rho$,
pressure $kR^{2}P$ and speed of sound $v^{2}/c^{2}$. The values of compactness
that were chosen for the graphic analysis are $u_{max}=0.19628$, $u=0.18086$,
$u=0.16545$, $u=0.15003$ and $u_{min}=0.13460$, where $u_{max}$ it’s
associated with the maximum mass $M=1.59M_{\odot}$ and the minimum radius
$R=11.96km$ $u=u_{max}=$, meanwhile $u_{min}$ is obtained by taking the
minimum mass $M=1.3M_{\odot}$ and the maximum radius $R=14.26km$.

Figure 1: Graphic representation of the density and the pressure for the
different values of compactness from the star PSR J0030+0451.

In the figure 1 we show the behaviour of the density and pressure for
different values of compactness which is obtained from the estimates in base
to observations. The graphics show that the density and pressure are positive
and monotonically decreasing functions, their values diminish as the
compactness value decreases, appearing in a more noticeable manner the
difference between the values of the density or the pressure in the center of
the star, we also observe that the pressure is zero on the surface. From the
figure 1 we observe that the Strong Energy Condition is satisfied, since both
the density and the pressure are positive. Also we have that for a specific
value of $u$ the value of the density is much greater than that of the
pressure ($c^{2}\rho>P$), which implies that the Dominant Energy Condition is
also satisfied. From the figure 2, graphic on the right, we observe that the
causality condition is met, since $0.2c^{2}<v^{2}<0.34c^{2}$ and that the
speed of sound is lower for lower compactness values, with maximum values on
the surface. The stability of the solution is guaranteed by the adiabatic
index, the left graph in the figure 2, with $\gamma$ being a monotonically
increasing function, the lowest value of the adiabatic index occurs in the
center of the star for the maximum compactness $u_{max}$ as such the set of
compactness values that is being analysed satisfies $\gamma>1.6843>4/3$.

Figure 2: Graphic representation of the speed of sound and the adiabatic index
for the different compactness values of the star PSR J0030+0451.

The absence of an event horizon and the continuity of the geometry on the
surface of the star is shown in the figure 3. In addition,in the graph on the
right side of the figure 3, we graph the forces present (the gravitational
force $F_{g}$ and the hydrostatic force $F_{h}$), identified by means of the
Tolman-Oppenheimer-Volkoff (TOV) equation

$-{\frac{\left({\it P_{r}}+c^{2}\rho\right)y^{\prime}}{y}}-{\it
P_{r}}^{\prime}=0,\qquad\Rightarrow\qquad F_{g}(r)=-{\frac{\left({\it
P_{r}}+c^{2}\rho\right)y^{\prime}}{y}},\quad F_{h}(r)=-{\it
P_{r}}^{\prime}.\qquad$

In the figure 3 we can observe the attractive effect of the gravitational
force $F_{g}$ countered by the hydrostatic repulsive force.

Figure 3: In the graph of the left side we present the behaviour of the metric
coefficients from the interior and exterior geometry, meanwhile in the graph
of the right side it´s shown the behaviour of the forces in the interior of
the star.

##  5 Discussion and conclusions

The graphic analysis has allowed us to show that the solution presented
satisfies every requirement which makes it physically acceptable and although
the graphic analysis was done considering estimated values of mass and radius
for the star PSR J0030+045, a similar behaviour is present for other values of
compactness, as long as $u\leq 0.23577$. To confirm that the behaviour of the
solution is not only graphically compatible but also that the orders of
magnitude which are obtained from the model are compatible with those
expected, on the tables 1 and 2 we report the values of density, pressure,
speed of sound and adiabatic index in the interior for the case of maximum 1
and minimum compactness 2.

Table 1: Interior behavior of the physical values of the density, pressure, speed of the sound and adiabatic index for the PSR J0030+0451, with $R=11.96km$ and $M=1.59M_{\odot}$, $u_{max}=0.19628$. $r(km)$ | $\\!\\!\rho\left(10^{17}\frac{kg}{m^{3}}\right)\\!\\!\\!$ | $\\!P(10^{33}Pa)$ | $v^{2}(c^{2})$ | $\gamma$
---|---|---|---|---
0. | 7.5125 | 9.4371 | 0.20654 | 1.6843
1.1960 | 7.4155 | 9.2563 | 0.20795 | 1.7052
2.3920 | 7.1398 | 8.7360 | 0.21214 | 1.7702
3.5880 | 6.7211 | 7.9258 | 0.21905 | 1.8887
4.7840 | 6.2131 | 6.9019 | 0.22857 | 2.0775
5.9800 | 5.6621 | 5.7429 | 0.24055 | 2.3719
7.1760 | 5.1146 | 4.5238 | 0.25474 | 2.8428
8.3720 | 4.5970 | 3.3040 | 0.27094 | 3.6581
9.5680 | 4.1283 | 2.1264 | 0.28887 | 5.3282
10.764 | 3.7155 | 1.0202 | 0.30824 | 10.395
11.960 | 3.3547 | 0 | 0.32880 | $\infty$
Table 2: Interior behavior of the physical values of the density, pressure, speed of sound and adiabatic index for the PSR J0030+0451 star, with $R=14.26km$ and $M=1.3M_{\odot}$, $u_{min}=0.13460$. $r(km)$ | $\\!\\!\rho\left(10^{17}\frac{kg}{m^{3}}\right)\\!\\!\\!$ | $\\!P(10^{33}Pa)$ | $v^{2}(c^{2})$ | $\gamma$
---|---|---|---|---
0 | 2.8411 | 2.1586 | 0.20605 | 2.6431
1.4260 | 2.8246 | 2.1275 | 0.20669 | 2.6723
2.8520 | 2.7760 | 2.0362 | 0.20873 | 2.7656
4.2780 | 2.6978 | 1.8894 | 0.21203 | 2.9321
5.7040 | 2.5965 | 1.6944 | 0.21667 | 3.2001
7.1300 | 2.4774 | 1.4588 | 0.22249 | 3.6176
8.5560 | 2.3462 | 1.1923 | 0.22956 | 4.2876
9.9820 | 2.2093 | 0.90520 | 0.23777 | 5.4501
11.408 | 2.0715 | 0.60586 | 0.24708 | 7.8340
12.834 | 1.9374 | 0.30200 | 0.25734 | 15.081
14.260 | 1.8055 | 0 | 0.26859 | $\infty$

From the tables 1 and 2 it can be noticed that the orders of magnitude of the
density and pressure are also those expected for the star PSR J0030+0451. With
which we can conclude that the model obtained is physically acceptable and
useful to represent stars with compactness $u\leq 0.23577$. Another relevance
of the solution constructed is that it can be useful as seed for obtaining new
physically acceptable solutions [35] in which we consider the contribution of
the presence of electric charge [36] or from an anisotropy factor [37], as
well as in the determination of new solutions in alternative gravitational
theories [38], investigations that could be developed in future works.

## Acknowledgments

We appreciate the facilities provided by the Universidad Michoacana de San
Nicolás de Hidalgo and the CIC -UMSNH during the realization of this
investigation as well as the CONACYT for the support given.

## References

  * [1] S. L. Shapiro and S. A. Teukolsk, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects, WILEY-VCH Verlag GmbH and Co. KGaA. (2004)
  * [2] R. Kippenhahn, A. Weigert and A. Weiss Stellar Structure and Evolution Secon Edition Springer-Verlag Berlin Heidelberg, (2012)
  * [3] R. M. Wald , General Relativity. University of Chicago Press, 1984.
  * [4] H.A. Buchdahl, Physical Review, 116, 1027 (1959).
  * [5] M. S. R. Delgaty, K. Lake, Comput. Phys. Commun. 115 (1998) 395.
  * [6] R. C. Tolman, Phys. Rev. 55, 364 (1939)
  * [7] G. K. Patwardhan and P. C. Vaidya, J. Univ. Bombay 12, 23 (1943).
  * [8] S. Nariai, Sci. Rep. Tohoku Univ., 34, 160 (1950)
  * [9] A. L. Mehra, J. Aust. Math. Soc. 6, 153 (1966)
  * [10] B. Kuchowicz, Acta Phys. Pol. 34, 131 (1968)
  * [11] H. Heintzmann, Z. Phys. 228, 489 (1969)
  * [12] R. J. Adler, J. Math. Phys. 15, 727 (1974)
  * [13] S. P. Goldman, Astrophys. J. 226, 1079 (1978)
  * [14] J. J. Matese and P. G. Whitman, Phys. Rev. D 22, 1270 (1980)
  * [15] M. C. Durgapal, J. Phys. A, Math. Gen. 15, 2637 (1982)
  * [16] N. R. Finch and J.E.F. Skea, Class and Quantum Gravity, 6, 467 (1989)
  * [17] N. Pant, R. N. Mehta and M. J. Pant, Astrophys Space Sci 330, 353 (2010).
  * [18] N. Pant, P. Fuloria and B. C. Tewari, Astrophys Space Sci 340, 407 (2012).
  * [19] N. Pant, p. Fuloria and N Pradhan, Int J Theor Phys 53, 993 (2014).
  * [20] S. Molla, M. Murshid and M. Kalam, Astrophys Space Sci 367, 4 (2022).
  * [21] M. H. Murad and N. Pant Astrophys Space Sci 350, 349 (2014).
  * [22] T. Gangopadhyay et al  M.N.R.A.S., 431, 3216 (2013),
  * [23] N. Pant, Astrophys Space Sci 331, 633 (2011).
  * [24] S.K. Maurya · Y.K. Gupta Astrophys Space Sci 334, 145 (2011)
  * [25] S. K. Maurya and Y. K. Gupta, Astrophys Space Sci 337, 151 (2012).
  * [26] G Estevez-Delgado, J Estevez-Delgado, N Montelongo Garcia and M. Pineda Duran, Mod. Phys. Lett. A 33, 1850237 (2018)
  * [27] G Estevez-Delgado, J Estevez-Delgado, J M Paulin-Fuentes, N Montelongo Garcia and M Pineda Duran, Mod. Phys. Lett. A 34, 1950115 (2019)
  * [28] G Estevez-Delgado, J Estevez-Delgado, N Montelongo Garcia, M Pineda Duran Can. J. Phys. 97, 988 (2019)
  * [29] G Estevez-Delgado, J Estevez-Delgado, J M Paulin-Fuentes, M. Pineda Duran, N. Montelongo García, Rev. Mex. Fís. 65, 392 (2019)
  * [30] J Estevez-Delgado, J V Cabrera, J. A. Rodriguez Ceballos, A Cleary- Balderas, M Paulin-Fuentes, Phys. Lett. A 35, 2050141 (2020)
  * [31] J Estevez-Delgado, N. Enrique Rodriguez Maya, J Martinez Peña and D. Rivera Rangel and Nancy Cambron Muñoz, Mod. Phys. Lett. A 36, 2150068 (2021)
  * [32] S. Weinberg, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity. (John Wiley and Sons 1972), pp 299-304
  * [33] Schutz B F. A First Course in General Relativity, 2nd ed. (Cambridge University Press 2009), pp 258-269.
  * [34] M. C. Miller et al ApJL 887 L24 (2019)
  * [35] P. Boonserm, M. Visser and S. Weinfurtner, Phys. Rev. D  71 (2005) 124037,
  * [36] B. V. Ivanov, Eur. Phys. J. C 81 (2021) 227,
  * [37] J Ovalle, R Casadio, R da Rocha, A Sotomayor, Eur. Phys. J. C 78 (2018) 122,
  * [38] A. M. Oliveira, H. E. S. Velten, J. C. Fabris, L. Casarini, Phys. Rev. D 92 (2015) 044020,


#  POLYIAMONDS ATTAINING EXTREMAL TOPOLOGICAL PROPERTIES, PART I

Greg Malen
Department of Mathematics  Duke University
e-mail address<EMAIL_ADDRESS>Érika Roldán
Department of Mathematics  The Ohio State University
e-mail address<EMAIL_ADDRESS>
###### Abstract

We consider two optimization questions with respect to polyiamonds. What is
the maximum number of holes that a polyiamond with $n$ tiles can enclose, and
what is the minimum number of tiles required to construct a polyiamond with
$h$ holes? These numbers will be given by the sequences $f_{\triangle}(n)$ and
$g_{\triangle}(h)$, respectively. In this paper, we construct a sequence of
polyiamonds with $h_{k}=(3/2)(k^{2}-k)$ holes and $n_{k}=(1/2)(9k^{2}+3k-4)$
tiles such that $g_{\triangle}(n_{k})=h_{k}$. Furthermore, these polyiamonds
all attain a specific set of efficient geometric and topological properties.

##  1 Introduction and main results

A polyiamond is a planar shape formed by gluing together finitely many
congruent equilateral triangles along their edges. If two triangles of a
polyiamond intersect, then their intersection is an entire edge, and the
gluing requirement implies that a polyiamond must have connected interior. We
only consider _free_ polyiamonds; that is two polyiamonds are equivalent if
they agree up to isometry. We refer to the triangles of a polyiamond as either
triangles or tiles.

One longstanding open problem in enumerative combinatorics is determining the
number of polyiamonds with exactly $n$ tiles. We denote the set of all
polyiamonds with $n$ tiles by $\mathcal{T}_{n}$, the cardinality of
$\mathcal{T}_{n}$ by $a(n)$, and the number of tiles in a polyiamond $A$ by
$|A|$. The values of the sequence $a(n)$ are known up to $n\leq 75$ [4]. It is
known that the sequence grows exponentialy fast and that its asymptotic growth
constant is equal to $\lim_{n\to\infty}a(n+1)/a(n)$ (see [1] and [2] for the
most recent improvements on the upper and lower bounds for this constant).

In the present paper, we study the extremal topological combinatorics problems
of maximizing the number of holes for a polyiamond with $n$ tiles, and
conversely, of minimizing the number of tiles required to create a polyiamond
with $h$ holes.

To be precise about the topology, we consider the tiles of a polyiamond to be
closed. Polyiamonds are finite unions of these closed tiles, so they are
compact. By Alexander duality [6], the number of holes of a polyiamond (the
rank of the first homology group) is the number of finite bounded, connected
components of the polyiamond’s complement in the plane. The area of a hole is
defined to be the number of tiles needed to fill it in.

Given a polyiamond $A$, we denote by $h(A)$ the number of holes in $A$. In
Figure 1, we show polyiamonds with 1, 2, and 3 holes, respectively. These
polyiamonds are optimal in the sense that it is impossible to create more
holes with the same number of tiles, and it is impossible to create these
numbers of holes with fewer tiles. These two topological and geometric
properties are precisely defined in the next two definitions.

######  Definition 1.1.

For $n\geq 1$, we define the sequence of the maximum number of holes that a
polyiamond with $n$ tiles can have as

$f_{\triangle}(n):=\max_{A\in\mathcal{T}_{n}}h(A).$ (1)

######  Definition 1.2.

For $h\geq 1$, we define the sequence of the minimum number of tiles needed
for constructing a polyiamond with $h$ holes as

$g_{\triangle}(h):=\min_{h(A)=h}|A|.$ (2)

As far as we know, this is the first time that the sequences
$f_{\triangle}(n)$ and $g_{\triangle}(h)$ are being defined and studied.

(a) $T_{1}$

(b) $T_{2}$

(c) $T_{3}$

Figure 1: Polyiamonds with the minimum number of tiles for up to three holes.

Here and throughout,

$h_{k}=3\binom{k}{2}=\frac{3}{2}k^{2}-\frac{3}{2}k,\qquad
n_{k}=\frac{9}{2}k^{2}+\frac{3}{2}k-2.$ (3)

Our main results are the following.

######  Theorem 1.1.

For $k>1$, $g_{\triangle}(h_{k})=n_{k}$

These values are realized by the sequence of spiral polyiamonds $Spir_{k}$,
with $n_{k}$ tiles and $h_{k}$ holes. The elements of this sequence also
exhibit efficient geometric and topological properties.

######  Theorem 1.2.

For all $k>1$ the dual graph of $Spir_{k}$ is a tree, the outer perimeter of
$Spir_{k}$ coincides with that of a shape of minimum perimeter for
$Spir_{k}+h(Spir_{k})$ many tiles, and every hole in $Spir_{k}$ has an area of
1.

Analogous sequences to $f_{\triangle}$ and $g_{\triangle}$ can be defined for
polyominoes instead of polyiamonds. A polyiomino is a planar shape formed by
gluing together a finitely many congruent squares along their edges. With
exhaustive computational work in [3] in 2015, Tomás Oliveira e Silva
enumerated all polyominoes with less than $28$ tiles and determined
$f_{\square}(n)$ for $n\leq 28$ and $g_{\square}(h)$ for $h\leq 8$.

In 2019, Kahle and Roldán were the first to explicitly study these sequences
for polyominoes in [7], establishing the first and second order asymptotics of
$f_{\square}(n)$ and obtaining the values of $g_{\square}(h)$ for a sequence
with exponential growth. In particular, they show that given
$C_{1}>\sqrt{5/2}$ and $C_{2}<\sqrt{3/2}$, there exists an
$n_{0}=n_{0}(C_{1},C_{2})$ such that

$\frac{1}{2}n-C_{1}\sqrt{n}\leq f(n)\leq\frac{1}{2}n-C_{2}\sqrt{n},$

for $n>n_{0}$.

They also prove, in an analogous statement as our result in Theorem 1.1, that
for

$s_{k}=\frac{1}{3}\left(2^{2k}-1\right),\qquad
l_{k}=\frac{1}{3}\left(2^{2k+1}+3\cdot 2^{k+1}+3\right),$

it is the case that $g_{\square}(s_{k})=l_{k}$ for all integers $k\geq 1$. In
proving this, they construct highly symmetric polyominoes with $s_{k}$ tiles
and $l_{k}$ holes (Figures 1 and 5 in [7]), which also satisfy the geometric
and topological structural requirements analogous to Theorem 1.2.

##  2 Preliminary Results

Define the perimeter $p(A)$ of a polyiamond to be the number of edges that are
on the topological boundary of $A$. For all integers $n\geq 1$ we denote by
$p_{\min}(n)$, the minimum perimeter possible for polyiamonds with $n$ tiles.
In 1976 [5], Harary and Harborth proved that

$p_{min}(n)=2\left\lceil\frac{1}{2}(n+\sqrt{6n})\right\rceil-n.$ (4)

Furthermore, they showed that the polyiamond in $\mathcal{T}_{n}$ in which
tiles are placed in a minimal hexagonal spiral attains $p_{min}(n)$ for every
$n$ (See Figure 3 in [5] and Figure 13 in [9]).

Although $p_{min}$ is not a monotone function, the increments are always $+1$
or $-1$. This fact appears as a comment along with Figure 13 in [9], and we
include it as part of the following lemma.

######  Lemma 2.1.

For any natural number $n$, $p_{min}(n+1)=p_{min}(n)\pm 1$. Let $k\geq 1$.
Then $p_{min}(n)-1\leq p_{min}(n+k)\leq p_{min}(n)+k$, and if
$p_{min}(n+1)=p_{min}(n)+1$, then $p_{min}(n)\leq p_{min}(n+k)$.

###### Proof.

As noted in [9], in the hexagonal spiral the $(n+1)$-th tile shares one edge
with the $n$-th tile, and if it is pointing outward from the center then it
shares a second edge with a tile that was placed in the previous layer of the
spiral. So if the $(n+1)$-th tile is pointing outward, it closes two edges of
the perimeter while adding one, and

$p_{min}(n+1)=p_{min}(n)-1.$

Conversely, if the $(n+1)$-th tile is pointing inward, it only shares the one
connecting edge with the $n$-th tile. So it closes one while adding two, and

$p_{min}(n+1)=p_{min}(n)+1.$

Thus $p_{min}(n+k)-p_{min}(n)\leq k$, since each additional tile can add at
most one edge to the perimeter.

Along a given side of this hexagonal shape, the tiles alternate between
outward and inward, except at the corners where two inward facing tiles are
placed one after the other (see Figure 13 in [9]). Hence along a given side,
$p_{min}$ alternates between values $i-1$ and $i$, and along the next side it
alternates between $i$ and $i+1$. Therefore once the threshold $i$ is reached
by $p_{min}(n)$, it can never drop below $i-1$ as $n$ increases. Furthermore,
the outward facing tiles achieve the smaller of the two perimeter values along
a side, and $p_{min}(n+1)=p_{min}(n)+1$ if and only if the $n$-th tile is
outward facing, so for such an $n$ we have $p_{min}(n)\leq p_{min}(n+k)$ for
any $k\geq 1$. ∎

The number of edges which are on the boundary of two tiles of a polyiamond $A$
will be denoted by $b(A)$. These are the edges contained in the interior of
$A$. Observe that all the edges of the individual tiles of $A$ necessarily
either belong to the perimeter or to the interior of $A$. This means that

$3n=p(A)+2b(A).$ (5)

Let $b_{min}(n)$ be the minimum number of interior edges that a polyiamond
with $n$ tiles can have. We associate a dual graph to a polyiamond by
considering each polygon as a vertex and by connecting any two of these
vertices if they share an edge, so that $b(A)$ is precisely the number of
edges in this graph.

Since polyiamonds have connected interiors by definition, the associated dual
graph of a polyiamond is connected and has a spanning tree. Then the dual
graph of a polyiamond with $n$ tiles has $n$ vertices, and its spanning tree
has $n-1$ edges. Hence $b_{\min}(n)\geq n-1$. Observing that the polyiamond
$C$ in which $n$ tiles are placed in a row has exactly $b(C)=n-1$, we conclude
that $b_{min}(n)=n-1$.

To distinguish which edges on the perimeter bound holes in a polyiamond $A$,
we define the hole perimeter $p_{h}(A)$ to be the number of edges bounding
holes, and the remaining outer perimeter $p_{out}(A)$ to be the number of
edges which do not bound holes, so that

$p(A)=p_{out}(A)+p_{h}(A).$ (6)

If a polyiamond $A$ is simply connected, then $p(A)=p_{out}(A)$. Using these
parameters, we formulate a necessary condition for the existence of a
polyiamond with $n$ tiles and $h$ holes.

###  2.1 A Topological Isoperimetric Inequality

######  Lemma 2.2.

If there exists a polyiamond with $n$ tiles and $h$ holes, then

$h\leq\frac{3n-2b_{min}(n)-p_{min}(n+h)}{3}.$ (7)

###### Proof.

Let $A$ be a polyiamond with $n$ tiles and $h$ holes, and total hole area
$h+j$ for some $j\geq 0$. Then by (5) and (6),

$p_{h}(A)=3n-2b(A)-p_{out}(A).\\\ $ (8)

Since $b(A)\geq b_{min}(n)$, and by filling in each hole we observe that
$p_{out}(A)\geq p_{min}(n+h+j)$, we have that

$p_{h}(A)\leq 3n-2b_{min}(n)-p_{min}(n+h+j).\\\ $ (9)

Holes in polyiamonds have at least three edges, so let $p_{h}(A)=3h+i$ for
some $i=i(j)\geq 0$. Then

$\begin{split}&3h+j\leq 3n-2b_{min}(n)-p_{min}(n+h+j)\\\ \Longrightarrow\hskip
17.07182pt&3h\leq 3n-2b_{min}(n)-p_{min}(n+h+j)-i.\end{split}$

If $j=i=0$, then trivially $p_{min}(n+h)=p_{min}(n+h+j)-i$. Otherwise $j,i>0$,
and by Lemma 2.1,

$\begin{split}&p_{min}(n+h)-1\leq p_{min}(n+h+j)\\\ \Longrightarrow\hskip
17.07182pt&p_{min}(n+h)\leq p_{min}(n+h+j)+1\leq p_{min}(n+h+j)+i.\end{split}$

Plugging this into the right side of the preceding inequality and dividing by
three, we have the desired result. ∎

For convenience, we define notation for this upper bound function,

$\begin{split}M(n,h)&=\frac{3n-2b_{min}(n)-p_{min}(n+h)}{3}\\\\[8.0pt]
&=\frac{n+2-p_{min}(n+h)}{3}.\end{split}$ (10)

######  Corollary 2.2.1.

If a polyiamond exists with $n$ tiles and $h$ holes, then $M(n,h)\geq h$. And
conversely if $M(n,h)<h$, then no such polyiamond can exist. Thus if a
polyiamond exists with $n$ tiles and $h$ holes such that
$n=\min\\{n:M(n,h)=h\\}$, then $g_{\triangle}(h)=n$.

Next we examine $M(n,h)$ as a function of $n$, and consider the set of $n$ for
which $M(n,h)\geq h$.

######  Lemma 2.3.

The function $M(n,h)$ is non decreasing in $n$.

###### Proof.

Let $i\geq 1$ be an integer. Then applying Lemma 2.1 to $M(n+i,h)$,

$\displaystyle M(n+i,h)$ $\displaystyle=\frac{n+i+2+2-p_{min}(n+h+i)}{3}$
$\displaystyle\geq\frac{n+i+2-p_{min}(n+h)-i}{3}$
$\displaystyle\geq\frac{n+2-p_{min}(n+h)}{3}$ $\displaystyle=M(n,h).$

∎

In Section 3, we construct a sequence of polyominoes $Spir_{k}$ and prove
Theorem 1.1 by showing that

$|Spir_{k}|=\min\\{n:M(n,h_{k})=h_{k}\\}.$ (11)

###  2.2 Basic Shapes in the Triangular Lattice

In the rest of this section we introduce some terminology to describe
arrangements of polyiamonds and to assist in counting. We take the equilateral
triangular lattice in the plane, oriented so that each triangle has a
horizontal edge. We call a triangle an _up_ -triangle if the horizontal edge
is at the base, and a _down_ -triangle if the horizontal edge is on top. Let
$\text{Hex}_{k}$ be the regular hexagon in this lattice with side length $k$,
and refer to the concentric hexagonal layers as
$L_{k}=\text{Hex}_{k}-\text{Hex}_{k-1}$. $\text{Hex}_{k}$ and $L_{k}$ are
polyiamonds themselves, but will primarily be used as metric benchmarks in
this lattice to describe the positions of the tiles of a given polyiamond.
Since we only consider free polyiamonds, these are all assumed to be centered
at the same origin.

(a) $\text{Hex}_{3}$

(b) $\text{Hex}_{3}=\text{Hex}_{2}\cup L_{3}$

Figure 2: Central hexagons and layers, with $L_{3}$ divided into trapezoids.

For every $k\geq 1$, $L_{k}$ decomposes into 6 trapezoidal row polyiamonds
each containing $k$ up-triangles and $k-1$ down-triangles, or vice versa (see
Figure 2). Thus,

$\begin{split}|L_{k}|&=6(k+k-1)=12k-6,\\\\[5.0pt]
|\text{Hex}_{k}|&=\sum_{i=1}^{k}(12i-6)=12\binom{k+1}{2}-6k=6k^{2}.\end{split}$
(12)

The sequences $h_{k}$ and $n_{k}$ from (3) are such that
$n_{k}+h_{k}=6k^{2}-2$. These sequences will be, respectively, the maximum
number of holes that can be created within $\text{Hex}_{k}$ and the minimum
number of tiles necessary to do so. The difference terms in these sequences
are

$h_{k+1}-h_{k}=3k.$ (13)

And

$n_{k+1}-n_{k}=9k+6.$ (14)

That is, in the $(k+1)$-th hexagonal layer $L_{k+1}$, we will create $3k$
holes by adding $9k+6$ triangles.

##  3 Proof of Theorems 1.1 and 1.2

We begin by constructing a sequence of spiral polyiamonds
$\\{Spir_{k}\\}_{k\geq 2}$. It is constructed inductively from $Spir_{2}$ and
$Spir_{3}$ in Figure 3.

(a) $Spir_{2}$

(b) $Spir_{3}$

Figure 3: The first two elements of the sequence $Spir_{k}$.

The $k$-th element in the sequence is built from $Spir_{k-2}$ by first adding
the last two tiles in the bottow right corner of $L_{k-2}$. Then add all down-
triangles in $L_{k-1}$, and all of $L_{k}$ except the two tiles in its bottom
right corner. Next, to ensure the interior is connected, take the penultimate
up-triangle, from left to right, in the bottom row of $L_{k-2}$ and shift it
down and to the right, to the penultimate up-space in the bottom row of
$L_{k-1}$ (see Figure 4). This creates $Spir_{k}$, a spiral of holes with
central $Spir_{2}$ and $Spir_{3}$ configurations, depending on whether $k$ is
even or odd, respectively. We depict $Spir_{15}$ in Figure 5.

$t$

$t$

Figure 4: From left to right: Fill in bottom corner of $L_{2}$ in $Spir_{2}$;
add all down-triangles in $L_{3}$ and all but the bottom right two tiles in
$L_{4}$; finally move the tile $t$ from $L_{2}$ to $L_{3}$ to get $Spir_{4}$.

######  Lemma 3.1.

For all $k>1$, we have that $|Spir_{k}|=n_{k}$ and $h(Spir_{k})=h_{k}$.

###### Proof of Lemma 3.1.

All elements of $\\{Spir_{k}$}, have holes outside of the central
configuration are up-triangles. Triangles in a layer are evenly split between
up and down, so $L_{k-1}$ contains $6k-9$ up-triangles by equation (12). In
building $Spir_{k}$ from $Spir_{k-2}$, every up-triangle in $L_{k-1}$ becomes
a hole except for the one which is filled by the tile from $L_{k-2}$ to
connect the spiral, which of course creates one up-triangle hole in $L_{k-2}$
that was not there before. So the number of holes added in this process is
precisely the number of up-triangles in $L_{k-1}$. Thus
$h(Spir_{k})=h(Spir_{k-2})+6k-9$.

From Figure 3 it is clear that $h(Spir_{2})=3=h_{2}$ and
$h(Spir_{3})=9=h_{3}$. Fix $k\geq 4$ and assume that the formula holds for all
values less than $k$. Then $h(Spir_{k-2})=h_{k-2}$, and

$\displaystyle h(Spir_{k})$ $\displaystyle=h_{k-2}+6k-9$
$\displaystyle=\frac{3}{2}k^{2}-\frac{3}{2}k.$

By induction, $h(Spir_{k})=h_{k}$ for all integers $k\geq 2$. And these
polyiamonds fill out all of $\text{Hex}_{k}$ except for the bottom right two
triangles, so

$|Spir_{k}|=6k^{2}-2-h_{k}=n_{k}.$

∎

###### Proof of Theorem 1.1.

In Lemma 3.1, we proved that for all $k>1$, we have that $|Spir_{k}|=n_{k}$
and $h(Spir_{k})=h_{k}$. Thus, to finish the proof we need to show that
$g_{\triangle}\left(h\left(Spir_{k}\right)\right)=|Spir_{k}|.$ It is easy to
check by a straight forward computation that $M(n_{k},h_{k})=h_{k}$ and that
$M(n_{k}-1,h_{k})<h_{k}$. Thus, by Corollary 2.2.1 we conclude that

$g_{\triangle}\left(h\left(Spir_{k}\right)\right)=|Spir_{k}|.$

∎

###### Proof of Theorem 1.2.

Observe first that for each $k>1$, every hole of the polyiamond $Spir_{k}$ has
an area of one. Also, by definition of $M(n,h)$ we have that
$b(Spir_{k})=n_{k}-1$ which implies that the dual graph of $Spir_{k}$ is a
tree. Finally, by (4) and the fact that $Spir_{k}$ has only holes of area 1 we
conclude that $p_{out}(Spir_{k})=p_{min}(n_{k}+h_{k})$. ∎

##  4 Concluding Remarks

In Part II of this work we completely answer this problem, determining the
value of $g_{\triangle}(h)$ for every $h\geq 1$. We call a polyiamond $A$ with
$|A|=g_{\triangle}\left(h(A)\right)$ crystallized, and prove that the set of
structural conditions in Theorem 1.2 is in fact equivalent to crystallization
for polyiamonds. This is in contrast to the situation for polyominoes [8],
where these optimal structural conditions imply crystallization but are not
equivalent.

Figure 5: $Spir_{15}$

## Acknowledgments.

This work was supported by HFSP RGP0051/2017, NSF DMS 17-13012, and NSF
DMS-1352386.

## References

  * [1] Gill Barequet and Mira Shalah.  Improved bounds on the growth constant of polyiamonds.  In Proceedings of the 32nd European Workshop on Computational Geometry, Lugano, Switzerland, pages 67–70, 2016.
  * [2] Gill Barequet, Mira Shalah, and Yufei Zheng.  An improved lower bound on the growth constant of polyiamonds.  J. Comb. Optim., 37(2):424–438, 2019.
  * [3] Tomás Oliveira e Silva.  Animal enumerations on the $\\{4,4\\}$ Euclidian tiling.  http://sweet.ua.pt/tos/animals/a44.html, December 2015.
  * [4] Anthony J. Guttmann and Iwan Jensen.  Appendix: Series data and growth constant, amplitude and exponent estimates.  In Polygons, Polyominoes and Polycubes, pages 469–482. Springer, 2009.
  * [5] Frank Harary and Heiko Harborth.  Extremal animals.  J. Combinatorics Information Syst. Sci., 1(1):1–8, 1976.
  * [6] Allen Hatcher.  Algebraic topology.  2005\.
  * [7] Matthew Kahle and Érika Roldán.  Polyominoes with maximally many holes.  Geombinatorics, Vol XXIX(Issue I):5–20, 2019.
  * [8] Greg Malen and Érika Roldán.  Topology and geometry of crystallized polyominoes.  arXiv preprint arXiv:1910.10342, 2019.
  * [9] Robert R. Meyer and Winston C. Yang.  Maximal and minimal polyiamonds.  Technical report, University of Wisconsin-Madison, 05 2002.