††thanks<EMAIL_ADDRESS>
# Quantum phases of the biased two-chain-coupled Bose-Hubbard Ladder

Jingtao Fan  State Key Laboratory of Quantum Optics and Quantum Optics
Devices, Institute of Laser Spectroscopy, Shanxi University, Taiyuan 030006,
China  Collaborative Innovation Center of Extreme Optics, Shanxi University,
Taiyuan 030006, China  Xiaofan Zhou  State Key Laboratory of Quantum Optics
and Quantum Optics Devices, Institute of Laser Spectroscopy, Shanxi
University, Taiyuan 030006, China  Collaborative Innovation Center of Extreme
Optics, Shanxi University, Taiyuan 030006, China  Suotang Jia  State Key
Laboratory of Quantum Optics and Quantum Optics Devices, Institute of Laser
Spectroscopy, Shanxi University, Taiyuan 030006, China  Collaborative
Innovation Center of Extreme Optics, Shanxi University, Taiyuan 030006, China

###### Abstract

We investigate the quantum phases of bosons in a two-chain-coupled ladder.
This bosonic ladder is generally in a biased configuration, meaning that the
two chains of the ladder can have dramatically different on-site interactions
and potential energies. Adopting the numerical density-matrix renormalization-
group method, we analyze the phase transitions in various parameter spaces. We
find signatures of both insulating-to-superfluid and superfluid-to-insulating
quantum phase transitions as the interchain tunnelling is increased.
Interestingly, tunning the interaction to some intermediate values, the system
can exhibit a reentrant quantum phase transition between insulating and
superfluid phases. We show that for infinite interaction bias, the model is
amenable to some analytical treatments, whose prediction about the phase
boundary is in great agreement with the numerical results. We finally clarify
some critical parameters which separate the system into regimes with distinct
phase behaviours, and briefly compare typical properties of the biased and
unbiased bosonic ladder systems. Our work enriches the Bose-Hubbard physics.

###### pacs:

42.50.Pq

##  I Introduction

Strongly correlated bosons, especially those moving in the periodic
potentials, have always been the research interest for both experimentalists
and theorists, as they are related to a variety of quantum phenomenon Boson1 ;
Boson2 . The simplest model describing such systems is the Bose-Hubbard (BH)
model, which incorporates the contributions from the kinetic energy of
individual atoms and the repulsive interactions between them BHM1 ; BHM2 ;
BHM3 ; BHM4 ; BHM5 ; BHM6 . Although originally developed in the context of
4He liquid BHM1 , it has been demonstrated that the BH model can be feasibly
implemented with ultracold atoms trapped in optical lattices Book1 ; OL1 ; OL2
; OL3 . Utilizing the unprecedented degree of controllability of the laser
fields, all the characteristic parameters of the BH model can be tuned in the
optical lattice with high precision OL2 ; OL3 . Relying on this, the quantum
phase transition from a superfluid (SF) to a Mott insulator (MI), which is the
most important prediction of the BH model, has been experimentally realized in
one OL4 , two OL5  and three dimensions OL6 . Since then, lots of related
studies have been performed on the extensions of the BH model, by considering,
for example, diverse forms of interactions ExtendBH1 ; ExtendBH2 ; ExtendBH3
and gauge fields ExtendBH4 . These extensions stimulate the development of new
directions bridging condensed matter physics, stochastic physics and quantum
optics.

In this context, the bosons confined to low-dimensional lattices merit special
attention, since the correlations built up in these systems are considerably
enhanced by the interactions between atoms OneD . Among various low-
dimensional lattice models, the two-chain-coupled BH ladder is of particular
importance BHLtheory1 ; BHLtheory2 , since it serves as an intermediate
geometry intervening the one and two lattice systems Ladder . This provides
beneficial insights about the characteristics of the SF-to-MI transition in
going from one to two dimensions. As a matter of fact, the BH ladder has been
experimentally simulated in different artificial systems BHLexperiment1 ;
BHLexperiment2 ; BHLexperiment3 , stimulating immense interests of research
towards various aspects of this model, such as the chiral currents Current1 ;
Current2 ; Current3 ; Current4 ; Current5 ; Current6 , the quantum magnetism
Magnetism1 ; Magnetism2 ; Magnetism3 ; Magnetism4 ; Magnetism5 , and the
topological states Topo1 ; Topo2 ; Topo3 ; Topo4 ; Topo5 ; Topo6 . The BH
ladders considered in these studies, however, are mostly limited to the
symmetric case where both the on-site interactions and potential energies are
identical for the two chains. Notice that general ladder systems should also
involve the configurations where the two chains have distinctly different
system parameters. Actually, letting the two chains of the ladder asymmetric,
with respect to either the interactions or potential energies, may impose
major impacts on various properties of the system Biased1 ; Biased2 ; Biased3
; Biased4 ; Biased5 . More importantly, this “biased” ladder structure also
underlies the physics behind a large class of systems, such as the dressed
dipolar molecules Dipolar1 ; Dipolar2  or Rydberg gases Rydberg1 ; Rydberg2 ;
Rydberg3  in optical lattices and the low-dimensional magnetic materials under
external fields Material1 ; Material2 .

In this work, we investigate the ground-state properties of a biased bosonic
ladder at half filling, using state-of-the-art density-matrix renormalization-
group (DMRG) numerical methods dmrg1 ; dmrg2 . By saying “biased”, we mean
that the two chains constituting the whole ladder can have dramatically
different on-site interactions and potential energies. We first provide an
analysis of the quantum phases in the limit of infinite interaction bias,
where the on-site interaction is infinite for one chain and finite for the
other. It is found that, as the interchain tunnelling is increased, either the
MI-to-SF or the SF-to-MI quantum phase transition can occur depending on the
value of interactions. More interestingly, tuning the finite interaction to
some intermediate values, the system may even exhibit a reentrant quantum
phase transition between MI and SF. By mapping the finite interaction into an
effective canonical Kerr nonlinear form, we analytically derive the phase
boundary between MI and SF, which agrees well with the numerical results. With
the knowledge of the system under infinite interaction bias, we then discuss
the more general parameter regime where the interactions of both chains of the
ladder are finite. We map out the ground-state phase diagrams in various
parameter spaces, and characterize several critical parameters which separate
the system into regimes with distinct phase behaviours. Finally, we briefly
compare the typical properties of the biased and unbiased bosonic ladder
systems.

##  II Model and method

Figure 1: (a) Schematic picture of the ladder system. The two chains
constituting the ladder are designated as spin-up and spin-down, respectively.
The interspin and intraspin hopping amplitudes are denoted, respectively, as
$t$ and $h$. (b) Possible implementation of the bosonic ladder in optical
superlattices. The optical double well is generally tilted by an energy
difference $\Delta$. The boson tunnelling rates along different directions
simulate the hopping rates $t$ and $h$.

As illustrated in Fig. 1(a), the system in consideration is a bosonic ladder
with two coupled chains, which we denote as spin-up and spin-down,
respectively. The interspin tunnelling is allowed along the rung. We assume
the bosonic ladder is typically biased, i.e., atoms with different spins
experience different potential energies and local interactions. Such a
scenario can be effectively engineered in spin-dependent optical lattices SOL1
or optical superlattices Current3 ; DOL1 , where the spin index distinguishing
different chains can be represented by either the hyperfine sublevels or
optical wells, according to different experiment implementations [see Fig.
1(b) for illustration]. The Hamiltonian describing this system reads

$\displaystyle\hat{H}\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!-t\sum_{\left\langle
i,j\right\rangle,\sigma}\hat{b}_{i,\sigma}^{{\dagger}}\hat{b}_{j,\sigma}-h\sum_{j}(\hat{b}_{j,\uparrow}^{{\dagger}}\hat{b}_{j,\downarrow}+\text{H.c.})$
(1)
$\displaystyle\\!\\!+\Delta\sum_{j}(\hat{n}_{j,\uparrow}-\hat{n}_{j,\downarrow})+\sum_{j,\sigma}\frac{U_{\sigma}}{2}\hat{n}_{j,\sigma}\left(\hat{n}_{j,\sigma}-1\right)$

where the field operator $\hat{b}_{j,\sigma}$
($\hat{b}_{j,\sigma}^{{\dagger}}$) annihilates (creates) a bosonic atom with
spin $\sigma$ (=$\uparrow,\downarrow$) at the lattice site $j$. While atoms
with the same spin can hop between adjacent sites $\left\langle
i,j\right\rangle$ with the intraspin hopping rate $t$, an interspin field
along the rung of the ladder couples atoms with different spins at rate $h$.
The energy bias $\Delta$, which we assume to be positive in this work, tends
to polarize the atoms along the rung of the ladder and $U_{\sigma}$ denotes
the on-site repulsive interaction of atoms with spin $\sigma$
($=\uparrow,\downarrow$). In this work, we focus on the commensurate ladder
with total atomic density to be $\rho=N/2L=1/2$. Here,
$N=N_{\uparrow}+N_{\downarrow}$ is the total number of bosons on the two-chain
ladder, each of which has length $L$. This amounts to setting the total system
size to be $2L$. In the following discussion, we set the energy scale by
taking $t=1$, and we also take $\Delta=10$ unless otherwise specified.

The Hamiltonian (1) can be viewed as a natural extension of the single-
component BH model incorporating the spin degree of freedom, which is
controlled by both the transverse and longitude magnetic fields. Without the
interspin tunnelling $h$, the ladder decouples and reduces to two independent
BH chains. A finite energy bias $\Delta$ ($>0$) then polarizes the bosons to
the spin-down chain with commensurability of one boson per site, leaving the
spin-up chain empty. In this case, the physics is entirely governed by the
one-dimensional BH model, which has been extensively explored BHM1 ; BHM2 ;
BHM3 ; BHM4 ; BHM5 ; BHM6 . With a non-zero interspin tunnelling $h$, however,
the two BH chains are coupled together, meaning that the characteristic
parameters of each individual chain may impact the ground-state properties of
the composite system in a collective manner. This becomes especially
interesting if the on-site interactions of the two chains are tuned quite
different. Without loss of generality, let us assume
$U_{\uparrow}>U_{\downarrow}$ and take a preanalysis on the behavior of the BH
chain with spin-down. In this case, the on-site interaction may effectively
enhanced through high-order tunnelling process triggered by $h$ Rydberg4 ,
favoring the formation of MI, whereas the particle filling factor may deviate
from unity at the same time, which in turn promotes the SF character. The
seemingly opposite tendency of the ground-state property makes the roles of
the parameters in the biased bosonic ladder less intuitive and to be
quantitatively clarified.

The MI and SF phases can be directly identified by calculating the charge gap
$\delta_{L}$, defined as the difference between the energies needed to add and
remove one particle from the system, i.e.,

$\delta_{L}=\mu_{L}^{+}-\mu_{L}^{-},$ (2) $\mu_{L}^{+}=E_{L}(N+1)-E_{L}(N),$
(3) $\mu_{L}^{-}=E_{L}(N)-E_{L}(N-1),$ (4)

where $E_{L}(N)$ is the the ground-state energy for $L$ sites and $N$
particles, and the chemical potentials $\mu_{L}^{+}$ and $\mu_{L}^{-}$
characterize the energy cost to add and remove one particle, respectively. The
insulating phase is signaled by the opening up of $\delta_{L}$ in the
thermodynamical limit $N,L\rightarrow\infty$ with fixed density $\rho$,
consisting with the zero compressibility $\kappa$
($=\partial\rho/\partial\mu$) of an insulator CF1 . In the SF phase, however,
the charge gap $\delta_{L}$ closes and the system becomes compressible in the
thermodynamical limit. Since $\delta_{L}$ keeps finite for any finite systems,
in order to pinpoint the MI-to-SF transition in the parameter space, we should
extrapolate to the $L\rightarrow\infty$ limit by utilizing the standard
finite-size scaling. Here, we perform state-of-the-art DMRG calculations to
compute the many-body ground state of the system, with which various physical
observable can be obtained. In our numerical simulations, we set the cutoff of
the single-site atom number as $n_{\text{cutoff}}=6$. We set lattice size up
to $L=40$, for which we retain 600 truncated states per DMRG block and perform
20 sweeps with a maximum truncation error of $\sim 10^{-9}$.

##  III Results

In the following, we systematically study the ground-state properties of the
biased bosonic ladder. Before providing the results of general parameters, we
first put the interaction bias to infinity, i.e., we consider the bosonic
ladder consisting of one chain with infinite on-site interaction and the other
with finite on-site interaction. The physics in this limit serves as a
beneficial starting point to understand the essential mechanism behind which
various phase behaviours occur.

###  III.1 Infinite interaction bias:
$U_{\uparrow}\\!-\\!U_{\downarrow}=\infty$

Without loss of generality, we fix the on-site interaction of the spin-up BH
chain to be infinity, $U_{\uparrow}\rightarrow\infty$, and that of the spin-
down chain to be finite. This amounts to imposing a hard-core constraint on
each site of the spin-up chain, where only one boson is allowed to occupy.

Before showing the full phase diagram, we can gain some useful insights into
the system by inspecting certain limit parameters. As mentioned in Sec. II,
the simplest limit is the zero interspin tunnelling $h=0$, under which the
ground state is fully described by the one-dimensional BH model with unit
filling. It is well known that the BH model with unit filling shows a SF-to-MI
transition at $U_{\text{c}}/t\sim 3.3$ BHM4 ; BHM5 ; BHM6 . The physics
becomes richer when $h$ is turned on. In this case, if we further set the on-
site interaction $U_{\downarrow}$ to be zero, the system closely resembles
that of two-level atoms inside cavity arrays for which the Jaynes-Cummings-
Hubbard (JCH) model works JCH1 ; JCH2 ; JCH3 ; JCH4 . This can be seen clearly
if we map the field operators of the hardcore bosons to those of quantum spins
by $\hat{b}_{j,\uparrow}\longrightarrow\sigma_{-}$ and
$\hat{b}_{j,\uparrow}^{{\dagger}}\longrightarrow\sigma_{+}$, where
$\sigma_{-}$ and $\sigma_{+}$ are spin-1/2 lowering and raising operators,
respectively. Therefore, the Hamiltonian describing the tunnelling process
between the two BH chains becomes Jaynes-Cummings type, in which the spin-down
bosons act as interaction-free photons bridging adjacent lattice sites. It
follows directly from the JCH physics that, by increasing $h$, the spin-down
bosons behave more localized and consequently undergoes a phase transition
from SF to MI at some critical tunnelling strength $h_{\text{c}}$ JCH1 ; JCH2
.

Anther interesting limit which is opposite to the JCH regime is
$U_{\downarrow}\rightarrow\infty$, implying a hard-core constraint on both
chains. It has been analytically demonstrated that, for the symmetric hard-
core BH ladder with $\Delta=0$, the critical tunneling strength $h_{\text{c}}$
decreases down to zero BHL1 . As a finite energy bias $\Delta$ usually tends
to increase the band gap, we do not expect the critical tunneling strength
$h_{\text{c}}$ to be shifted away from that of the symmetric case (i.e.,
$\Delta=0$).

Figure 2: The charge gap $\delta_{L}$, extroplated to $L\rightarrow\infty$, as
a function of $h$ with (a) $U_{\downarrow}=0.0$, (b) $U_{\downarrow}=2.0$, (c)
$U_{\downarrow}=3.6$ and (d) $U_{\downarrow}=20.0$. The other parameters are
$\Delta=10$ and $U_{\uparrow}=\infty$.

We now start to show the numerical results of general parameters. In order to
clarify the effect of the interspin tunnelling on the SF-to-MI transition, we
vary $h$ from zero to some large value and calculate the corresponding charge
gap $\delta_{L}$, which ought to be extrapolated to the $L\rightarrow\infty$
limit. From the knowledge of the one-dimensional BH model, the system stays at
the SF phase (MI phase) for $U_{\downarrow}\lesssim 3.3$
($U_{\downarrow}\gtrsim 3.3$) and $h=0$, and potential phase transitions can
take place when increasing $h$. Figure 2(a) plots the charge gap $\delta_{L}$
as a function of $h$ for $U_{\downarrow}=0$. It can be seen that, starting
from zero, the charge gap gradually opens up when $h$ exceeds the critical
value $h_{\text{c}}\approx 12.5$, evidencing a SF-to-MI transition, as can be
inferred from the JCH physics JCH2 . The critical tunnelling strength
$h_{\text{c}}$ decreases as we increase the on-site interaction
$U_{\downarrow}$, as shown in Fig. 2(b). In Fig. 2(d), we exemplify another
limit where the on-site interaction is considerably strong by setting
$U_{\downarrow}=20$. It can be found that, in this case, the gap $\delta_{L}$
keeps open irrespective of the value of $h$, meaning the system remains a MI.
Something interesting happens when the on-site interaction $U_{\downarrow}$ is
tuned slightly larger than $U_{\text{c}}$. As illustrated in Fig. 2(c), we
plot $\delta_{L}$ versus $h$ for $U_{\downarrow}=3.6$. With the increase of
$h$, the gap first closes at $h\approx 2.0$ and then reopens at $h\approx
5.1$, indicating that the phase transition appears twice. That is, the system
starts from the MI, and subsequently traverses the SF phase, and ending up in
the MI eventually. Figures 3(a) and 3(b) show finite-size scaling of the DMRG
data of the charge gap, by linear and quadratic fittings, for two
representative points located in the SF and MI phases, respectively. This
reentrant MI phase transition induced by the interspin tunnelling strength
$h$, does not exist in the symmetric bosonic ladder with
$U_{\uparrow}=U_{\downarrow}$ and $\Delta=0$, and is thus exclusive for the
biased ladder here.

Figure 3: The finite-size scaling of chemical potentials $\mu_{L}^{+}$ and
$\mu_{L}^{-}$ for (a) $h=3.0$ and (b) $h=8.0$. The other parameters are
$\Delta=10$, $U_{\downarrow}=3.6$ and $U_{\uparrow}=\infty$.

With the understanding above, we map out the phase diagram in the
$U_{\downarrow}-h$ plane in Fig. 4. The phase boundary has been extrapolated
to the $L\rightarrow\infty$ limit by the finite-size scaling. It is to be seen
clearly that, while increasing the on-site interaction $U_{\downarrow}$ always
drives the system to the MI phase, the role of the interspin tunnelling $h$
can be somehow opposite, i.e., it can trigger both the MI and SF phases,
depending on the value of $U_{\downarrow}$. Notice that the MI in the BH model
is essentially stabilized by the direct interaction between bosons, whereas
bosons with different spins are dressed together here forming composite
polaritons. We therefore expect that some effective interaction between
polaritons, which plays the key role in inducing different behaviours of the
phase transition, may emerge.

Figure 4: The phase diagram in the $U_{\downarrow}-h$ plane for $\Delta=10$
and $U_{\uparrow}=\infty$. The red solid line with square symbol (black solid
line) denotes The MI-SF phase boundary obtained by the DMRG calculation [Eq.
(7)]. The values of $h_{\text{m}}$ for each $U_{\downarrow}$ is also
pinpointed in the phase diagram by the blue solid line with circle symbol. For
comparation, location of the minimum of $U_{\text{eff}}$, determined by Eq.
(8), is plotted by the green dashed line. The shaded area, bounded by critical
interactions $U_{\text{c}}\approx 3.3$ and $U_{c1}\approx 3.9$, characterizes
the parameter region where the MI-to-SF-to-MI transition can occur.

To see this clearly, we map the local contribution of the Hamiltonian (1) into
an effective kerr nonlinearity by a simple energy mismatch argument. As
detailed in Appendix A, the local energy of Hamiltonian (1) consists of two
plaritonic modes whose eigenenergies are

$\displaystyle\omega_{n}^{\pm}$ $\displaystyle=$
$\displaystyle\frac{U_{\downarrow}}{2}n(n-2)+\frac{1}{2}(\Delta+U_{\downarrow})$
(5)
$\displaystyle\pm\frac{1}{2}\sqrt{(nU_{\downarrow}-U_{\downarrow}-\Delta)^{2}+4nh^{2}}$

where $n$ is the excitation number. Since we are only interested in the low-
energy physics, the focus in the following will be on the lower branch
$\omega_{n}^{-}$. We define the effective Hubbard interaction $U_{\text{eff}}$
as the energy cost incurred by forming a two-particle plaritonic excitation
(with energy $\omega_{2}^{-}$) from two single-particle plaritonic excitation
(with energy 2$\omega_{1}^{-}$) in neighboring lattice sites JCH3 ; JCH4 ,
i.e.,

$\displaystyle U_{\text{eff}}$ $\displaystyle=$
$\displaystyle\omega_{2}^{-}-2\omega_{1}^{-}$ (6) $\displaystyle=$
$\displaystyle\frac{1}{2}(U_{\downarrow}-\Delta)+\sqrt{\Delta^{2}+4h^{2}}$
$\displaystyle-\frac{\sqrt{(U_{\downarrow}-\Delta)^{2}+8h^{2}}}{2}$

With this understanding, we can obtain an analytical expression of the phase
boundary between the MI and SF phases by equating the effective interaction
$U_{\text{eff}}$ with the critical interaction strength of the BH model with
unit filling, namely

$U_{\text{eff}}=U_{\text{c}}\approx 3.3$ (7)

As shown in Fig. 4, the curve defined by Eq. (7), agrees well with the
numerical results obtained by the DMRG calculation. It should be emphasized
that, in deriving Eq. (6), we have implicitly assumed that the ground-state
property of the whole lattice system is mainly governed by its low-energy
local physics. This requires that (i) the energy scale owned by each local
lattice sites is considerably larger than the kinetic energy of bosons, namely
at least $\Delta\gg 1$ or $h\gg 1$, and (ii) the density fluctuations are weak
enough so that only the lowest-lying excitations of individual lattice sites
need to be taken into consideration. This guarantees the effectiveness of Eq.
(7) in predicting the MI-to-SF phase boundary, since the density fluctuations
are extremely suppressed in the MI. Eq. (6) provides further guidance to the
driving force inducing different phase transitions. An interesting finding is
that $U_{\text{eff}}$ exhibits nonmonotonic behaviour as $h$ increases from
zero. As illustrated in Fig. 5(a), with the increase of $h$, the effective
interaction $U_{\text{eff}}$ decreases first to a minimum and then increases
monotonically [see blue line], which explains the MI-to-SF-to-MI transition
found in Fig. 4. The location of the minimum of $U_{\text{eff}}$ can be easily
deduced by requiring $\partial U_{\text{eff}}/\partial h=0$, yielding a
trivial solution $h=0$ and a nontrivial solution,

$h=\frac{\sqrt{\Delta^{2}-(\Delta-U_{\downarrow})^{2}}}{2}.$ (8)

It is straightforward to show that Eq. (8), within its range of values,
minimizes $U_{\text{eff}}$. The curve obtained from Eq. (8) is depicted in
Fig. 4. Notice that, whereas $U_{\text{eff}}$ is minimized by $h=0$ when
$U_{\downarrow}=0$, consistent with the JCH physics JCH1 ; JCH2 , a nonzero
$U_{\downarrow}$ shifts the location of the interaction minimum (i.e., $h=0$)
to some finite value. Within this picture, an upper bound of $U_{\downarrow}$,
beyond which no SF phase would exist, can be obtained. This is immediately
achieved by substituting Eq. (8) into Eq. (7), which is then solved by

$U_{\downarrow}=U_{\text{c1}}\equiv\Delta+U_{\text{c}}-\sqrt{\Delta^{2}-U_{\text{c}}^{2}}.$
(9)

It becomes clear that, the parameter region of $U_{\downarrow}$ within which
the MI-to-SF-to-MI transition can occur is $U_{\text{c}}<$
$U_{\downarrow}<U_{\text{c1}}$.

Figure 5: (a) The effective Hubbard interaction $U_{\text{eff}}$ as a function
of $h$ for different $U_{\downarrow}$ with $\Delta=10$ and
$U_{\uparrow}=\infty$. (b) The condensate fraction calculated for different
system sizes. Note that the result of $L\rightarrow\infty$ is obtained by
extrapolation using the finite-size scaling. $h_{\text{m}}$ specifies the
location of the maximum of the condensate fraction. The other parameters are
$\Delta=10$, $U_{\downarrow}=3.6$ and $U_{\uparrow}=\infty$.

An experimental measurable quantity that is able to mirror the effective
interaction is the condensate fraction (CF), defined as the number of bosons
in the condensate with respect to the total number of bosons CF1 ; CF2 . It
has been shown that the condensate fraction of Bose gases monotonically
decreases as the local interaction increases CF1 . For the bosonic ladder
considered here, the CF is defined as the largest eigenvalue of the matrix
$\left\langle\hat{b}_{i,\sigma}^{{\dagger}}\hat{b}_{j,\sigma^{\prime}}\right\rangle$
divided by the total number of bosons CF2 . Figure 5(b) shows the CF as a
function of $h$ for $U_{\downarrow}=3.6$ and different system sizes. It is
demonstrated that the CF increases first, reaching its maximum at $h\approx
3.4$, and then decreases. The location of the maximum of CF, designated as
$h_{\text{m}}$, depends sensitively on the value of $U_{\downarrow}$. As shown
in Fig. 4, we plot $h_{\text{m}}$ for varying $U_{\downarrow}$, which exhibits
the same behaviour as that obtained from Eq. (8). The agreement between
$h_{\text{m}}$ and Eq. (8) signals that the picture of the effective
interaction $U_{\text{eff}}$ works in a wide range of parameters, even inside
the SF phase where the density fluctuation is somehow enhanced.

###  III.2 Finite interaction bias: $U_{\uparrow}\\!-\\!U_{\downarrow}<\infty$

Figure 6: The phase diagram in the $U_{\uparrow}-h$ plane with $\Delta=10$ and
(a) $U_{\downarrow}=3.0$, (b) $U_{\downarrow}=3.6$, (c) $U_{\downarrow}=5.0$,
and (d) $U_{\downarrow}=15.0$.

Having understood the physics of bosonic ladder under the infinite interaction
bias, we are now in the stage to explore the more general parameter regime
where the interactions of both chain of the ladder are finite. Here we are
particularly interested in the influence of finite spin-up interaction on
various quantum phases. By calculating the charge gap $\delta_{L}$ with
extrapolation to the thermodynamic limit, we obtain the phase diagrams in the
$U_{\uparrow}-h$ plane in Figs. 6(a)-(d) with different $U_{\downarrow}$. As
shown in Fig. 6(a), in which the spin-down interaction is fixed as
$U_{\downarrow}=3.0\ $($<$ $U_{\text{c}}$), the SF region is confined by a
smooth phase boudary, which extends up to $U_{\uparrow}\rightarrow\infty$ and
$h\rightarrow\infty$. A phase transition from the SF to MI may occur when
increasing $U_{\uparrow}$ ($h$) for some fixed $h$ ($U_{\uparrow}$).
Increasing the spin-down interaction slightly larger than $U_{\text{c}}$, for
example $U_{\downarrow}=3.6$, the MI can emerge for small $h$, penetrating the
SF region, as illustrated in Fig. 6(b). Importantly, as $h$ approaches
infinity, the spin-up interaction delimiting different quantum phases
decreases and saturates to some critical value $U_{\text{c2}}$.

In fact, through an analysis of the plaritonic modes, the MI-to-SF phase
boundary in the $h\rightarrow\infty$ limit can be derived as
$U_{\uparrow}+U_{\downarrow}=4U_{\text{c}}\approx 13.2$ (see Appendix B for
details). Setting $U_{\uparrow}=U_{\downarrow}=U$, we immediately reproduce
the result of the symmetric case, i.e., $U=2U_{\text{c}}\approx 6.6$, obtained
by the bosonization method BHLtheory1 . Under this framework, the critical
interaction $U_{\text{c2}}$ is straightforwardly written as

$U_{\text{c2}}=4U_{\text{c}}-U_{\downarrow}.$ (10)

As marked in Fig. 6(b) by blue dashed line, the critical interaction
$U_{\text{c2}}$ defined above separates the phase diagram into two distinct
parameter regimes. For the $U_{\uparrow}>U_{\text{c2}}$ side, there exists the
interesting MI-to-SF-to-MI phase transition we explored in Subsection III.1,
whereas for $U_{\uparrow}<U_{\text{c2}}$, the MI-to-SF phase transition can
appear only once by monotonically varying $h$.

Adopting the description of the effective interaction introduced in Subsection
III.1, we anticipate that if $U_{\downarrow}$ is increased to be larger than
$U_{\text{c1}}$, a finite upper bound of $U_{\uparrow}$, beyond which the SF
phase disappears, can emerge. This is confirmed by the phase diagram in Fig.
6(c), where we take $U_{\downarrow}=5$ ($>U_{\text{c1}}$=3.9). As expected,
the SF phase is destroyed for $U_{\uparrow}\gtrsim 11$, in contrast to the
behaviour shown in Figs. 6(a) and 6(b). Increasing $U_{\downarrow}$ further
such that $U_{\downarrow}>4U_{\text{c}}$, the critical interaction
$U_{\text{c2}}$ touches zero, meaning that the SF disappears, at least when
$h$ is sufficiently large. The phase diagram with $U_{\downarrow}=15$
($>4U_{\text{c}}$) is plotted in Fig. 6(d), from which we find that the area
of SF completely vanishes.

Up to now, our focus are basically on the parameter regime where both the
interactions and potential energies of the two chains are asymmetric. The
individual effect incurred by one of the two asymmetric ingredients, i.e.,
either the interaction asymmetry or the potential energy asymmetry, has not
been elucidated. Here we complement this study by plotting two additional
phase diagrams, each of which has only one asymmetric ingredient. The phase
diagram in the $U_{\uparrow}-h$ plane with zero energy bias ($\Delta=0$) and
fixed spin-down interaction ($U_{\downarrow}=3.6$) is plotted in Fig. 7(a).
The phase diagram in this case shares the same structure with that in Fig.
6(a), albeit with shrunken SF area. It is also understood that no MI phase can
be found for sufficiently small $h$, contrasting the behaviour in Fig. 6(b),
since a zero $\Delta$ always closes the charge gap for lattices with non-
integer filling at $h=0$. As shown in Fig. 7(b), by requiring the interactions
of the two chains to be equal, saying $U=U_{\uparrow}=U_{\downarrow}$, we map
the phase diagram in the $U-h$ plane with finite energy bias $\Delta=10$. With
the increase of $h$, the critical interaction of the SF-to-MI transition
monotonically increases, asymptotically up to $U_{\text{c2}}\approx 6.6$,
showing a distinct behaviour compared to cases with asymmetric interactions.

Figure 7: (a) The phase diagram in the $U_{\uparrow}-h$ plane with
$U_{\downarrow}=3.6$ and $\Delta=0$. (b) The phase diagram in the $U-h$ plane
with $U=U_{\uparrow}=U_{\downarrow}$ and $\Delta=10$.

##  IV Discussion and conclusion

As mentioned in Sec. II, the considered model can be directly implemented with
ultracold atoms inside optical lattices under various experiment designs. For
example, the bosonic ladder can be prepared by growing optical superlattices,
which forms a double-well structure along one direction Current3 ; DOL1 . The
tunnelling strength $h$, on-site interactions $U_{\uparrow/\downarrow}$ and
energy bias $\Delta$ can be independently controlled by properly tuning the
geometry of the optical double well. Alternatively, one can employ the spin-
dependent optical lattices SOL1 , where atoms with different hyperfine states
experience different lattice potentials. In this scenario, $\Delta$ and $h$
are respectively controlled by the detuning and Rabi frequency of an
additional coupling laser, and the on-site interactions
$U_{\uparrow/\downarrow}$ can be tuned via Feshbach resonances or the lattice
depths experienced by atoms with different spins. In spite of the research
interests of the model in its own right, our results offer beneficial insights
into engineering effective interactions on demand by dressing different atomic
internal states Rydberg1 ; Rydberg4 . That said, our model constitutes only a
subset of the rich physics in the BH ladder with coupled chains, and many
interesting extensions are to be applied in the future. For example, with a
ladder structure, the hopping process of atoms may carry nontrivial peierls
phases, giving rise to synthetic gauge fields Gauge . These gauge fields may
not only affect the MI-to-SF transitions dramatically ExtendBH4  but also
induce various chiral currents Current1 ; Current2 ; Current3 ; Current4 ;
Current5 ; Current6 . Another direction is to fit the system into the grand-
canonical description by introducing a tunable chemical potential BHM1 ;
Biased1 . This may provide new perspectives on the magnetic or charge
correlations in Mott lobes with different filling factors.

In conclusion, we have theoretically studied the ground-state properties of
the BH ladder with half filling in a biased configuration by using state-of-
the-art DMRG numerical methods. It is found that the interchain tunnelling can
drive both the MI-to-SF and SF-to-MI quantum phase transitions, depending on
the value of interactions. A reentrant quantum phase transition between MI and
SF has also been predicted by setting the on-site interactions to some
intermediate values. Under appropriate conditions, the model is shown to be
amenable to some analytical treatment, whose predictions about the phase
boundary is in great agreement with the numerical results. Armed with these
knowledge, we have mapped out the full phase diagram and characterized some
critical parameters separating the system into regimes with distinct phase
behaviours.

## Acknowledgments

This work is supported by the National Key R&D Program of China under Grant
No. 2022YFA1404003, the National Natural Science Foundation of China (NSFC)
under Grant No. 12004230, 12174233 and 12034012, the Research Project
Supported by Shanxi Scholarship Council of China and Shanxi ’1331KSC’.

##  Appendix A Plaritonic modes of the local Hamiltonian

In this Appendix, we derive Eq. (5) in the main text. To that end, we
rearrange the Hamiltonian (1) as

$\hat{H}=\sum_{j}\hat{H}_{L}^{(j)}-t\sum_{\left\langle
i,j\right\rangle,\sigma}\hat{b}_{i,\sigma}^{{\dagger}}\hat{b}_{j,\sigma}$ (11)

where

$\displaystyle\hat{H}_{L}^{(j)}$ $\displaystyle=$
$\displaystyle-h(\hat{b}_{\uparrow}^{{\dagger}}\hat{b}_{\downarrow}+\text{H.c.})+\Delta(\hat{n}_{\uparrow}-\hat{n}_{\downarrow})$
(12)
$\displaystyle+\sum_{\sigma}\frac{U_{\sigma}}{2}\hat{n}_{\sigma}\left(\hat{n}_{\sigma}-1\right)$

describes the local physics at lattice site $j$. Note that we have omitted the
subscript $j$ in the right hand side of Eq. (12) for simplicity. The
Hamiltonian (12) can be spanned by the Fock basis
$\left|n_{\downarrow},n_{\uparrow}\right\rangle$, where $\hat{n}_{\sigma}$ is
the occupation number for bosons with spin $\sigma$($=\uparrow,\downarrow$).
In the $U_{\uparrow}\rightarrow\infty$ limit, a hardcore constraint on the
spin-up bosons can be imposed, meaning that we only need to retain the states
$\left|n,0\right\rangle$ and $\left|n-1,1\right\rangle$ with the total
occupation $n=n_{\downarrow}+n_{\uparrow}$. The plaritonic modes of
Hamiltonian (12) are therefore admixtures of $\left|n,0\right\rangle$ and
$\left|n-1,1\right\rangle$, and the eigenenergies are readily diagonalized as

$\displaystyle\omega_{n}^{\pm}$ $\displaystyle=$
$\displaystyle\frac{U_{\downarrow}}{2}n(n-2)+\frac{1}{2}(\Delta+U_{\downarrow})$
(13)
$\displaystyle\pm\frac{1}{2}\sqrt{(nU_{\downarrow}-U_{\downarrow}-\Delta)^{2}+4nh^{2}}$

##  Appendix B Effective low-energy descritption in the $h\rightarrow\infty$
limit

Here we provide an effective low-energy decryption of the model in the
$h\rightarrow\infty$ limit. We first introduce two branches of quasi-modes
$\hat{b}_{j,+}$ and $\hat{b}_{j,-}$ defined as

$\hat{b}_{j,+}=\frac{1}{\sqrt{2}}(\hat{b}_{j,\uparrow}+\hat{b}_{j,\downarrow})$
(14)
$\hat{b}_{j,-}=\frac{1}{\sqrt{2}}(\hat{b}_{j,\uparrow}-\hat{b}_{j,\downarrow})$
(15)

Under the transformations of Eqs. (14) and (15), the Hamiltonian (1) is
rewritten as

$\hat{H}=\sum_{j}\hat{H}_{L}^{(j)}-t\sum_{\left\langle
i,j\right\rangle}\left(\hat{b}_{i,+}^{{\dagger}}\hat{b}_{j,+}+\hat{b}_{i,-}^{{\dagger}}\hat{b}_{j,-}\right)$
(16)

where the local Hamiltonian reads

$\displaystyle\hat{H}_{L}^{(j)}$ $\displaystyle=$
$\displaystyle(\frac{U_{\downarrow}}{8}+\frac{U_{\uparrow}}{8})[(\hat{n}_{j,+}+\hat{n}_{j,-})^{2}+(\hat{b}_{j,+}^{{\dagger}}\hat{b}_{j,-}+\hat{b}_{j,-}^{{\dagger}}\hat{b}_{j,+})^{2}-2(\hat{n}_{j,+}+\hat{n}_{j,-})]$
(17)
$\displaystyle+(\frac{U_{\downarrow}}{4}-\frac{U_{\uparrow}}{4})[(\hat{b}_{j,+}^{{\dagger}}\hat{b}_{j,-}+\hat{b}_{j,-}^{{\dagger}}\hat{b}_{j,+})(\hat{n}_{j,+}+\hat{n}_{j,-}-1)]$
$\displaystyle+\frac{\Delta}{2}(\hat{n}_{j,+}+\hat{n}_{j,-}-\hat{b}_{j,+}^{{\dagger}}\hat{b}_{j,-}-\hat{b}_{j,-}^{{\dagger}}\hat{b}_{j,+})+h(\hat{n}_{j,+}-\hat{n}_{j,-})$

It follows from Eq. (17) that the low-energy physics is dominated by bosons on
the “$-$” plaritonic branch in the $h\rightarrow\infty$ limit. We thus
anticipate an effective low-energy theory which is purely described by field
operators of the “$-$” plaritonic branch. The simplest way to achieve this is
to average the Hamiltonian (16) with respect to the vacuum state of the the
“$+$” plaritonic branch, yielding

$\displaystyle\hat{H}_{\text{eff}}$ $\displaystyle=$
$\displaystyle\sum_{j}\left[(\frac{\Delta}{2}-h)\hat{n}_{j,-}+\frac{\tilde{U}}{2}\hat{n}_{j,-}(\hat{n}_{j,-}-1)\right]$
(18) $\displaystyle-t\sum_{\left\langle
i,j\right\rangle}\hat{b}_{i,-}^{{\dagger}}\hat{b}_{j,-}$

where $\tilde{U}=(U_{\downarrow}+U_{\uparrow})/4$. Notice that the effective
description in the Hamiltonian (18) becomes accurate when $h$ approaches
infinity. More importantly, the Hamiltonian (18) is written in the same form
of the one-dimensional BH model with effective on-site interaction
$\tilde{U}$. It follows that the physics of our ladder system in this limit
can be effectively described by the one-dimensional BH model with simple
substitution of system parameters. Given this, the SF-to-MI phase boundary is
readily obtained as $\tilde{U}=U_{\text{c}}\approx 3.3$.

## References

  * (1) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).
  * (2) C. Orzel, A. K. Tuchman, M. L. Fenselau, M. Yasuda, and M. A. Kasevich, Squeezed States in a Bose-Einstein Condensate, Science 291, 2386 (2001).
  * (3) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Boson localization and the superfluid-insulator transition, Phys. Rev. B 40, 546 (1989).
  * (4) W. Krauth, M. Caffarel, and J.-P. Bouchaud, Gutzwiller wave function for a model of strongly interacting bosons, Phys. Rev. B 45, 3137 (1992).
  * (5) A. P. Kampf and G. T. Zimanyi, Superconductor-insulator phase transition in the boson Hubbard model, Phys. Rev. B 47, 279 (1993).
  * (6) J. K. Freericks and H. Monien, Strong-coupling expansions for the pure and disordered Bose-Hubbard model, Phys. Rev. B 53, 2691 (1996).
  * (7) V. A. Kashurnikov and B. V. Svistunov, Exact diagonalization plus renormalization-group theory: Accurate method for a one-dimensional superfluid-insulator-transition study, Phys. Rev. B 53, 11776 (1996).
  * (8) T. D. Kühner and H. Monien, Phases of the one-dimensional Bose-Hubbard model, Phys. Rev. B 58, R14741 (1998).
  * (9) M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices (Oxford University Press, United Kingdom, 2012).
  * (10) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Cold Bosonic Atoms in Optical Lattices, Phys. Rev. Lett. 81, 3108 (1998).
  * (11) E. Demler and F. Zhou, Spinor Bosonic Atoms in Optical Lattices: Symmetry Breaking and Fractionalization, Phys. Rev. Lett. 88, 163001 (2002).
  * (12) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (13) T. Stöferle, H. Moritz, C. Schori, M. Köhl, and T. Esslinger, Transition from a Strongly Interacting 1D Superfluid to a Mott Insulator, Phys. Rev. Lett. 92, 130403 (2004).
  * (14) I. B. Spielman, W. D. Phillips, and J. V. Porto, Mott-Insulator Transition in a Two-Dimensional Atomic Bose Gas, Phys. Rev. Lett. 98, 080404 (2007).
  * (15) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature 415, 39 (2002).
  * (16) A. Kuklov, N. Prokof’ev, and B. Svistunov, Commensurate Two-Component Bosons in an Optical Lattice:Ground State Phase Diagram, Phys. Rev. Lett. 92, 050402 (2004).
  * (17) V. W. Scarola and S. Das Sarma, Quantum Phases of the Extended Bose-Hubbard Hamiltonian: Possibility of a Supersolid State of Cold Atoms in Optical Lattices, Phys. Rev. Lett. 95, 033003 (2005).
  * (18) O. Mansikkamäki, S. Laine, and M. Silveri, Phases of the disordered Bose-Hubbard model with attractive interactions, Phys. Rev. B 103, L220202 (2021).
  * (19) R. Sachdeva, M. Singh, and T. Busch, Extended Bose-Hubbard model for two-leg ladder systems in artificial magnetic fields, Phys. Rev. A 95, 063601 (2017).
  * (20) M. A. Cazalilla, R. Citro, T. Giamarchi, E. Orignac, and M. Rigol, One dimensional bosons: From condensed matter systems to ultracold gases, Rev. Mod. Phys. 83, 1405 (2011).
  * (21) P. Donohue and T. Giamarchi, Mott-superfluid transition in bosonic ladders, Phys. Rev. B 63, 180508 (R) (2001).
  * (22) M. S. Luthra, T. Mishra, R. V. Pai, and B. P. Das, Phase diagram of a bosonic ladder with two coupled chains, Phys. Rev. B 78, 165104 (2008).
  * (23) E. Dagotto and T. M. Rice, Surprises on the Way from One- to Two-Dimensional Quantum Magnets: The Ladder Materials, Science 271, 618 (1996).
  * (24) M. Atala, M. Aidelsburger, Mi. Lohse, J. T. Barreiro, B. Paredes and I. Bloch, Observation of chiral currents with ultracold atoms in bosonic ladders, Nature Physics 10, 588 (2014).
  * (25) S. Greschner, M. Piraud, F. Heidrich-Meisner, I. P. McCulloch, U. Schollwőck, and T. Vekua, Spontaneous Increase of Magnetic Flux and Chiral-Current Reversal in Bosonic Ladders: Swimming against the Tide, Phys. Rev. Lett. 115, 190402 (2015).
  * (26) Q. Zhu, Z.-H. Sun, M. Gong, F. Chen, Y.-R. Zhang, Y. Wu, Y. Ye, C. Zha, S. Li, S. Guo, et al., Observation of Thermalization and Information Scrambling in a Superconducting Quantum Processor, Phys. Rev. Lett. 128, 160502 (2022).
  * (27) E. Orignac and T. Giamarchi, Meissner effect in a bosonic ladder, Phys. Rev. B 64, 144515 (2001).
  * (28) A. Petrescu and K. L. Hur, Bosonic Mott Insulator with Meissner Currents, Phys. Rev. Lett. 111, 150601 (2013).
  * (29) M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Observation of chiral currents with ultracold atoms in bosonic ladders, Nat. Phys. 10, 588 (2014).
  * (30) Y. Jian, X. Qiao, J.-C. Liang, Z.-F. Yu, A.-X. Zhang, and J.-K. Xue, Stability and superfluidity of the Bose-Einstein condensate in a two-leg ladder with magnetic field, Phys. Rev. A 104, 024212 (2021).
  * (31) C.-M. Halati, and T. Giamarchi, Bose-Hubbard triangular ladder in an artificial gauge field, Phys. Rev. Research 5, 013126 (2023).
  * (32) X. Zhou , S. Jia, and X.-W. Luo, Spin-tensor Meissner currents of ultracold bosonic gases in an optical lattice, Phys. Rev. A 108, 013304 (2023).
  * (33) S. Daul and D. J. Scalapino, Frustrated Hubbard ladders and superconductivity in $\kappa$-BEDT-TTF organic compounds, Phys. Rev. B 62, 8658 (2000).
  * (34) P.-B. He, Q. Sun, P. Li, S.-Q. Shen, and W. M. Liu, Magnetic quantum phase transition of cold atoms in an optical lattice, Phys. Rev. A 76, 043618 (2007).
  * (35) K. Totsuka, P. Lecheminant, and S. Capponi, Semiclassical approach to competing orders in a two-leg spin ladder with ring exchange, Phys. Rev. B 86, 014435 (2012).
  * (36) J. Fan, X. Zhou, W. Zheng, W. Yi, G. Chen, and S. Jia, Magnetic order in a Fermi gas induced by cavity-field fluctuations, Phys. Rev. A 98, 043613 (2018).
  * (37) X. Zhou, X.-W. Luo, G. Chen, S. Jia, and C. Zhang, Quantum spiral spin-tensor magnetism, Phys. Rev. B 101, 140412(R) (2020).
  * (38) M. Mancini, G. Pagano, G. Cappellini, L. Livi, M. Rider, J. Catani, C. Sias, P. Zoller, M. Inguscio, M. Dalmonte, L. Fallani, Observation of chiral edge states with neutral fermions in synthetic Hall ribbons, Science 349, 1510 (2015).
  * (39) B. K. Stuhl, H.-I. Lu, L. M. Aycock, D. Genkina, I. B. Spielman, Visualizing edge states with an atomic Bose gas in the quantum Hall regime, Science 349, 1514 (2015).
  * (40) S.-L. Zhang and Q. Zhou, Two-leg Su-Schrieffer-Heeger chain with glide reflection symmetry, Phys. Rev. A 95, 061601(R) (2017).
  * (41) K. Monkman and J. Sirker, Operational entanglement of symmetry-protected topological edge states, Phys. Rev. Research 2, 043191 (2020).
  * (42) P. Sompet, S. Hirthe, D. Bourgund, T. Chalopin, J. Bibo, J. Koepsell, P. Bojović, R. Verresen, F. Pollmann, G. Salomon, C. Gross, T. A. Hilker and I. Bloch, Realizing the symmetry-protected Haldane phase in Fermi–Hubbard ladders, Nature 606, 484 (2022).
  * (43) Y. Kuno and Y. Hatsugai, $Z_{2}\times Z_{2}$ symmetry and $Z_{4}$ Berry phase of bosonic ladders, Phys. Rev. A 108, 013319 (2023).
  * (44) I. Danshita, J. E. Williams, C. A. R. Sá de Melo, and C. W. Clark, Quantum phases of bosons in double-well optical lattices, Phys. Rev. A 76, 043606 (2007).
  * (45) I. Danshita, C. A. R. Sá de Melo, and C. W. Clark, Reentrant quantum phase transition in double-well optical lattices, Phys. Rev. A 77, 063609 (2008).
  * (46) H. Deng, H. Dai, J. Huang, X. Qin, J. Xu, H. Zhong, C. He, and C. Lee, Cluster Gutzwiller study of the Bose-Hubbard ladder: Ground-state phase diagram and many-body Landau-Zener dynamics, Phys. Rev. A 92, 023618 (2015).
  * (47) X. Qiao, X.-B. Zhang, Y. Jian, A.-X. Zhang, Z.-F. Yu, and J.-K. Xue, Quantum phases of interacting bosons on biased two-leg ladders with magnetic flux, Phys. Rev. A 104, 053323 (2021).
  * (48) A. Padhan, R. Parida, S. Lahiri, M. K. Giri and T. Mishra, Quantum phases of constrained bosons on a two-leg Bose-Hubbard ladder, Phys. Rev. A 108, 013316 (2023).
  * (49) G. Pupillo, A. Griessner, A. Micheli, M. Ortner, D.-W. Wang, and P. Zoller, Cold Atoms and Molecules in Self-Assembled Dipolar Lattices, Phys. Rev. Lett. 100, 050402 (2008).
  * (50) A. V. Gorshkov, S. R. Manmana, G. Chen, J. Ye, E. Demler, M. D. Lukin, and A. M. Rey, Tunable Superfluidity and Quantum Magnetism with Ultracold Polar Molecules, Phys. Rev. Lett. 107, 115301 (2011).
  * (51) L. I. R. Gil, R. Mukherjee, E. M. Bridge, M. P. A. Jones, and T. Pohl, Spin Squeezing in a Rydberg Lattice Clock, Phys. Rev. Lett. 112, 103601 (2014).
  * (52) Z. Lan, J. Minář, E. Levi, W. Li, and I. Lesanovsky, Emergent Devil’s Staircase without Particle-Hole Symmetry in Rydberg Quantum Gases with Competing Attractive and Repulsive Interactions, Phys. Rev. Lett. 115, 203001 (2015).
  * (53) Y. He, Z. Bai, Y. Jiao, J. Zhao, and W. Li, Superradiance-induced multistability in one-dimensional driven Rydberg lattice gases, Phys. Rev. A 106, 063319 (2022).
  * (54) M. Azuma, Z. Hiroi, M. Takano, K. Ishida, and Y. Kitaoka, Observation of a Spin Gap in SrCu2O3 Comprising Spin-$\frac{1}{2}$ Quasi-1D Two-Leg Ladders, Phys. Rev. Lett. 73, 3463 (1994).
  * (55) T. Hong, K. P. Schmidt, K. Coester, F. F. Awwadi, M. M. Turnbull, Y. Qiu, J. A. Rodriguez-Rivera, M. Zhu, X. Ke, C. P. Aoyama, et al., Magnetic ordering induced by interladder coupling in the spin-$\frac{1}{2}$ Heisenberg two-leg ladder antiferromagnet C9H18N2CuBr4, Phys. Rev. B 89, 174432 (2014).
  * (56) S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * (57) U. Schollwök, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * (58) G. Jotzu, M. Messer, F. Görg, D. Greif, R. Desbuquois, and T. Esslinger, Creating State-Dependent Lattices for Ultracold Fermions by Magnetic Gradient Modulation, Phys. Rev. Lett. 115, 073002 (2015).
  * (59) J. Sebby-Strabley, M. Anderlini, P. S. Jessen, and J. V. Porto, Lattice of double wells for manipulating pairs of cold atoms, Phys. Rev. A 73, 033605 (2006).
  * (60) N. Henkel, R. Nath, and T. Pohl, Three-Dimensional Roton Excitations and Supersolid Formation in Rydberg-Excited Bose-Einstein Condensates, Phys. Rev. Lett. 104, 195302 (2010).
  * (61) S. Ramanan, Tapan Mishra, Meetu Sethi Luthra, Ramesh V. Pai, and B. P. Das, Signatures of the superfluid–to–Mott-insulator transition in cold bosonic atoms in a one-dimensional optical lattice, Phys. Rev. A 79, 013625 (2009).
  * (62) A. D. Greentree, C. Tahan, J. H. Cole, and L. C. L. Hollenberg, Quantum phase transition of light, Nat. Phys. 2, 856 (2006).
  * (63) J. Koch and K. L. Hur, Superfluid–Mott-insulator transition of light in the Jaynes-Cummings lattice, Phys. Rev. A 80, 023811 (2009).
  * (64) C.-W. Wu, M. Gao, Z.-J. Deng, H.-Y. Dai, P.-X. Chen, and C.-Z. Li, Quantum phase transition of light in a one-dimensional photon-hopping-controllable resonator array, Phys. Rev. A 84, 043827 (2011).
  * (65) T. Grujic, S. R. Clark, D. Jaksch, and D. G. Angelakis, Non-equilibrium many-body effects in driven nonlinear resonator arrays, New. J. Phys. 14, 103025 (2012).
  * (66) F. Crépin,N. Laflorencie, G. Roux, and P. Simon, Phase diagram of hard-core bosons on clean and disordered two-leg ladders: Mott insulator-Luttinger liquid-Bose glass, Phys. Rev. B 84, 054517 (2011).
  * (67) A. J. Leggett, Bose-Einstein condensation in the alkali gases: Some fundamental concepts, Rev. Mod. Phys. 73, 307 (2001).
  * (68) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Artificial gauge potentials for neutral atoms, Rev. Mod. Phys. 83, 1523 (2011).


# Efficiently Enumerating Scaled Copies of Point Set Patterns

Aya Bernstine111 School of Computer Science and Engineering , The Hebrew
University, Jerusalem, Israel<EMAIL_ADDRESS>and Yehonatan
Mizrahi222 School of Computer Science and Engineering , The Hebrew University,
Jerusalem, Israel<EMAIL_ADDRESS>
###### Abstract

Problems on repeated geometric patterns in finite point sets in Euclidean
space are extensively studied in the literature of combinatorial and
computational geometry. Such problems trace their inspiration to Erdős’
original work on that topic. In this paper, we investigate the particular case
of finding scaled copies of any pattern within a set of $n$ points, that is,
the algorithmic task of efficiently enumerating all such copies. We initially
focus on one particularly simple pattern of axis-parallel squares, and present
an algorithm with an $O(n\sqrt{n})$ running time and $O(n)$ space for this
task, involving various bucket-based and sweep-line techniques. Our algorithm
is worst-case optimal, as it matches the known lower bound of
$\Omega(n\sqrt{n})$ on the maximum number of axis-parallel squares determined
by $n$ points in the plane, thereby solving an open question for more than
three decades of realizing that bound for this pattern. We extend our result
to an algorithm that enumerates all copies, up to scaling, of any full-
dimensional fixed set of points in $d$-dimensional Euclidean space, that works
in time $O(n^{1+1/d})$ and space $O(n)$, also matching the corresponding lower
bound due to Elekes and Erdős.

##  1 Introduction

The problems of geometric point pattern matching and the identification of
repeated geometric patterns are fundamental computational problems with a
myriad of applications, ranging from computer vision [13, 11], image and video
compression [2], model-based object recognition [16], structural biology [12]
and even computational chemistry [9]. Such problems were originally motivated
by questions from extremal combinatorics regarding the maximal number of
occurrences of a given pattern determined by a set of points, a field
historically inspired by Erdős’ well-known Unit Distance Problem (1946)
regarding the maximal number of unit distance pairs induced by such sets [8].
Our paper approaches the computational problems of identifying patterns using
tools and techniques encountered in the framework of computational geometry,
ensuring exact, provably correct and efficient solutions. For a detailed
survey that discusses the applicability of these methods and their relevance
to such problems, see [1].

In this paper, we analyze a particular case of interest – the problem of
identifying shifted and scaled copies of any point set pattern in Euclidean
space, termed homothetic copies of this set. For brevity, we refer to those as
copies of the pattern. Note that each copy is certainly determined by the
images of any two points in the set. We begin with focusing on one of the
simplest forms of these problems, due to its historical significance, where
the goal is to identify repeated patterns of axis-parallel squares in the
plane. Simply put, we strive for an efficient algorithm that lists all subsets
of four points that constitute the vertices of a square with axis-parallel
edges, an algorithm that could be easily generalized to treat the
identification of any pattern in any dimension. Note that coping with the
challenge of enumerating all copies of a given pattern trivially encompasses a
solution for the task of finding one single copy of that pattern in a
background set of points, namely the task of pattern matching. This manifests
the importance of dealing with the more general of those two challenges,
thereby solving both at one fell swoop.

This problem was extensively studied in the literature of combinatorial and
computational geometry. As articulated in 1990 by van Kreveld and de Berg
[15], the maximum possible number of axis-parallel squares determined by $n$
points in the plane is $\Theta(n\sqrt{n})$, and given any set of $n$ points,
all such squares can be enumerated in time $O(n\sqrt{n}\log n)$333The analysis
given throughout this paper of time and space complexities is based on the
conventional and very standard word-RAM model of computation [10]. and space
$O(n)$. They consolidate their result by showing an extension of their
algorithm to one that enumerates all full-dimensional axis-parallel hypercubes
in $d$-dimensional Euclidean space in time $O(n^{1+1/d}\log n)$, with a
logarithmic-factor gap separating this computational result from the
corresponding combinatorial geometric bound of a maximum of
$\Theta(n^{1+1/d})$ possible hypercubes, raising the challenge of overcoming
that gap as an open question. The latter combinatorial result was further
extended by Elekes and Erdős [7]. They established a bound of
$\Theta(n^{1+1/d})$ on the maximum number of copies of any pattern in
$\mathbb{R}^{d}$, assuming the pattern is full-dimensional and has rational
coordinates. The computational aspect of this generalization occurs in [3],
providing an algorithm that works in time $O(n^{1+1/d}\log n)$ for the task of
enumerating all copies of any fixed full-dimensional set consisting of a
constant number of points in $d$-dimensional Euclidean space, exhibiting the
same logarithmic-factor gap between the two results. We remark that the
complexity of finding copies in these algorithms becomes smaller with
increasing dimension, a rare occurrence.

###  1.1 Our Results

Our main result of this paper is an efficient deterministic algorithm that
enumerates all scaled copies of any given pattern consisting of a constant
number of points in $d$-dimensional Euclidean space. As earlier stressed, we
first handle the case of listing all axis-parallel squares in the plane.
Although being a special case that follows immediately from our main theorem,
there are several motivations for its individual treatment. First, it arguably
forms the simplest symmetric full-dimensional pattern in the plane.
Consequently, it earned its right of having the first attempt, as we are aware
of, among all problems of listing scaled copies of some pattern, of admitting
an improved performance compared to the trivial algorithm, in terms of running
time analysis. Adaptations to other patterns appeared, e.g., in [3], but [15]
were the first to raise the question of whether it is computationally feasible
to realize the combinatorial bound of $\Theta(n\sqrt{n})$ possible axis-
parallel squares, thereby improving their result having the running time of
$O(n\sqrt{n}\log n)$. Our algorithm which is, in some sense, an amplified
version of their original algorithm, is therefore a worst-case optimal
solution for this problem, and fully addresses this question which was open
for more than three decades.

A second motivation for separately handling this case, is that it results in a
relatively natural algorithm having less obfuscating technicalities as
compared to the general case of treating arbitrary patterns. In other words,
it crystallizes the additional techniques we incorporate within the original
solution of [15], yielding the required improvement. Briefly sketched, the
gist of these techniques includes the following ideas. First, we use a
reduction from arbitrary input points to points having ”compressed”
coordinates, that is, we relabel their values, allowing the use of linear non-
comparative sorting methods on the result. Second, we deploy a sweep-line sub-
procedure that marks points forming a square, instead of searching those in a
set, avoiding the corresponding logarithmic factor without involving
probabilistic tools as hashing. Third, in order to carry out the mentioned
sweep-line process correctly, we also relabel the sum and the difference of
the input coordinates, in addition to the relabeling of the coordinates
themselves. We show why the last step is crucial for the algorithm to succeed
in the full elaboration of our solution in Section 2.

Theorem: Given a planar set of points $P$ of size $n$, all axis-parallel
squares defined by points from $P$ can be enumerated in time $O(n\sqrt{n})$
and $O(n)$ space.

As previously stressed, this algorithm matches, in the worst-case, the known
bound of $\Theta(n\sqrt{n})$ on the maximum possible number of axis-parallel
squares. We further extend this result to a companion proposition, that
discusses the task of reporting all full-dimensional hypercubes in
$d$-dimensional Euclidean space. A solution for the latter task that works in
time $O(n^{1+1/d})$ and space $O(n)$, which also matches the corresponding
known lower bound and addresses a similar open question, is provided.

Establishing our main result, which is concerned with finding scaled copies of
any arbitrary and full-dimensional pattern in $d$-dimensions, relies on the
majority of the ideas behind the proof of the previous theorem. Assuming the
pattern consists of a constant number of points, we provide an algorithm with
a running time of $O(n^{1+1/d})$. Similarly to the proof idea of the previous
theorem mentioned earlier, to accomplish this task we relabel some particular
and carefully chosen affine transformations of the input coordinates along
with the coordinates themselves, a relabeling that creates a convenient
representation of the points, for the purpose of scanning them using a sweep-
line technique.

Theorem: Given a fixed set of points $Q$ of full dimension in the
$d$-dimensional Euclidean space, and a set of points $P$ of size $n$, all
scaled copies of $Q$ defined by points from $P$ can be enumerated in time
$O(n^{1+1/d})$ and $O(n)$ space.

The running time in this theorem matches the corresponding combinatorial bound
of the same magnitude, and improves the best known running time for this task
of $O(n^{1+1/d}\log n)$ [3]. We provide the proof of the last theorem in
Section 3.

Note that although the improvement suggested is by a logarithmic factor, the
upshot is a worst-case optimal algorithm in terms of running time analysis,
even for the most general case of arbitrary patterns. This can be compared
with a work by de Rezende and Lee [6], who studied the problem of enumerating
all rotated copies of a given pattern, improving the running time of the
trivial algorithm for this task by a logarithmic factor as well. Other than
that, aside from the worst-case optimality of our results, the techniques
deployed form a rather general scheme, and may therefore be potentially useful
to treat other variants of the problem studied.

###  1.2 More Related Work

An excellent survey that covers a variant of our problem, that deals with
rotated copies instead of scaled ones, along with combinations of rotations
and scaling, affine transformations of patterns, Ramsey-type graph-theoretic
related problems and more can be found in [4].

A result related to this variant, namely that of rotated copies, by Cho and
Mount [5] considers the problem of approximated pattern matching allowing only
rotations. They showed an efficient algorithm that admits a
$3.13$-approximation to the optimum. Namely, it returns a transformation whose
Hausdorff distance is at most a factor $3.13$ larger than the optimal
Hausdorff distance from the pattern.

In [15], the problem of counting the number of all axis-parallel rectangles
(boxes in $d$-dimensions) is also considered, providing an algorithm with a
running time of $O(n^{2-1/d})$. Note that as opposed to enumerating axis-
parallel squares, it is an increasing function of $d$.

###  1.3 Paper Organization

The rest of the paper is organized as follows. In Section 2 we provide an
algorithm that efficiently enumerates all axis-parallel squares. An
elaboration on the techniques used for the speeding-up previous results is
also given. Section 3 generalizes the result from Section 2 to full-
dimensional hypercubes, and ultimately to an efficient enumeration scheme of
all copies of any arbitrary point set pattern. Concluding remarks and
directions for further research are given in Section 4.

##  2 Axis-Parallel Squares

The purpose of the current section is to present an efficient algorithm that
reports all axis-parallel squares defined by a planar set of points. That is,
given $n$ points on the plane, the algorithm addresses the task of listing all
subsets of four points that form the vertices of an axis-parallel square.
Obviously, this task can be tackled by a trivial algorithm that sorts the
input points, say, according to their $x$ coordinate, and then scans every
pair of points residing in the same column and determines whether that pair
can be complemented to a square from the right, by searching in a set, for
instance. However, the performance of this algorithm is rather poor in terms
of complexity analysis, as it results in a running time of $O(n^{2}\log n)$,
e.g., when all points share the same $x$ coordinate.

A more efficient algorithm, devised by van Kreveld and de Berg [15], manages
to improve upon this by first treating only columns with at most $\sqrt{n}$
points, in an identical fashion to that described in the trivial algorithm,
then deleting those columns, and finally applying a similar procedure on the
remainder, only considering rows instead of columns. In this section, we refer
to columns with at most $\sqrt{n}$ points as ”short” columns. The full
description of their algorithm is as follows:

Squares-Listing$(p_{1},...,p_{n})$:

  1. 1.

Build a balanced search tree $T$ and an array $A$ on the input points, both
sorted by the $x$ coordinate.

  2. 2.

For every pair of points $p$ and $q$ in $A$ residing in a short column, search
in $T$ whether they can be complemented to a square from the right or from the
left. Report each square found, unless the other two vertices defining it are
on a short column to the left of $p$ and $q$, to avoid reporting the same
square more than once.

  3. 3.

Delete all short columns from $T$ and $A$, and convert each remaining point
$(x,y)$ to $(y,x)$.

  4. 4.

Apply step 2 on the remaining converted points.

The correctness of the algorithm is proven in their work, along with its
complexity analysis. They show that it operates in a running time of
$O(n\sqrt{n}\log n)$ and $O(n)$ space. As we integrate most of its ideas in
our amplified version, we present the analysis of the relevant steps in detail
in the proof of our theorem in the current section, later on.

We strive for an algorithm that does not include the above logarithmic factor,
i.e., one with a running time of $O(n\sqrt{n})$ and space $O(n)$. As shown in
[15], there are configurations of $n$ points that define $\Theta(n\sqrt{n})$
squares, or even more generally:

###### Theorem 2.1

(van Kreveld, de Berg) For a set $P$ of $n$ points in d-dimensional space, the
maximum possible number of $2^{d}$ points that form the vertices of an axis-
parallel hypercube is $\Theta(n^{1+1/d})$.

This theorem obviously induces a lower bound on the running time of the
optimal algorithm for the enumeration task. Thus, our result in this section
bridges the logarithmic-factor gap between this lower bound, and the
previously best known upper bound for this objective, attained by Squares-
Listing$(p_{1},...,p_{n})$.

###  2.1 Main Ideas Towards an Improvement

The logarithmic factor in the running time of Squares-
Listing$(p_{1},...,p_{n})$ is present due to the step in which we search
points in a set. There are several steps we deploy in order to overcome this.
Assume first, for simplicity, that all input points have coordinates in
$\\{1,...,n\\}$. Then, instead of searching in a set the two query pairs that
complement the pair $(x,y),(x,y+\delta)$ to a square, i.e., the pair
$(x+\delta,y),(x+\delta,y+\delta)$ and the pair
$(x-\delta,y),(x-\delta,y+\delta)$, we put all query points (points being part
of some query pair) along with the original input points in an array, apply
radix sort on that array, treating each point as a two-digit number, and mark
all positive queries. That is, we mark each query point adjacent to an
existing input point sharing the same coordinates, or to an already marked
identical query point. The resulting marked pairs define the appropriate
existing squares. Using this non-comparative sorting algorithm, along with its
linear scanning that follows, we obtain a process that does not involve the
logarithmic factor of marking a query point.

Generally, however, we cannot assume that all coordinates are taken from
$\\{1,...,n\\}$. We address this issue by ”shrinking” the coordinates of all
input points. In other words, we label each $x$ coordinate in the input with a
value in $\\{1,...,n\\}$, and likewise for the $y$ coordinate, so effectively
we reduce the problem to the aforementioned restricted version. The main
caveat of this idea, though, is that arithmetical considerations regarding
labels are invalid. For example, if we map the pair of points $(10,20)$ and
$(10,40)$ to the labeled points $(1,2)$ and $(1,4)$ respectively, the pre-
translated pair may be complemented to a square using the pair $(30,20)$ and
$(30,40)$, but there is no guarantee that $(30,20)$, if exists in the input,
was mapped to $(3,2)$.

To overcome the above issue, we avoid using arithmetic considerations when
defining a query pair of points $q_{1},q_{2}$ that complement the post-labeled
pair (i.e., the points after being renamed) $p_{1}=(x,y),p_{2}=(x,y+\delta)$
to a square (from the right, assuming that $\delta>0$, for example). That is,
instead of using the invalid label $x+\delta$ as a coordinate, we define the
pair $q_{1},q_{2}$ using identical labels as those of $p_{1},p_{2}$ in the
following manner: Note that $q_{2}$, aside from sharing the same $y$
coordinate as $p_{2}$ (assuming $q_{2}$ lies above $q_{1}$, for instance), it
also shares the same Manhattan distance from the origin as $p_{1}$444This is
true assuming all points have positive coordinates, but we can safely assume
this, as otherwise we can shift the input so it maintains that form, without
affecting our result. , namely they have the same sum of coordinates. The full
elaboration of this consideration appears in the course of the correctness
proof of our algorithm. Therefore, if we label also each sum of coordinates in
the input points, similarly to the labeling of the coordinates themselves, we
can search for that exact label and avoid the arithmetic considerations. In
short – the query point $q_{2}$ is defined having the same $y$ label as
$p_{2}$ and the same $x+y$ label as $p_{1}$. Similarly, we label also the
differences of the coordinates, so $q_{1}$ has the same $y$ label as $p_{1}$
and the same $y-x$ label as $p_{2}$, as we show in the proof of our following
theorem. Note that the label of $x+y$ need not be the sum of the labels of $x$
and $y$ in this approach. Similar considerations apply for the pair that
complements to a square from the left.

Another observation related to the idea described in the last paragraph, which
will be useful in further developments of the current section’s algorithm, is
that rotating any vector $(x,y)$ in the plane by $45^{o}$ and stretching it by
a factor of $\sqrt{2}$ yields the vector $(x+y,y-x)$, as illustrated in Figure
1. This technique can be shown to establish, for example, the famous geometric
equivalence between the $l_{1}$ and $l_{\infty}$ norms (see, e.g., [14]),
namely that the $l_{1}$ distance between any two pre-rotated points equals the
$l_{\infty}$ distance between the post-rotated points, after scaling. For our
purposes, we note that instead of viewing the process as labeling the $x+y$
and $y-x$ values of the coordinates, as suggested above, one can think of this
process as the labeling of the post-rotated points’ coordinates – in the case
of axis-parallel squares, a rotation by the aforementioned factor.

Figure 1: Illustrating the rotation by $45^{o}$ and the stretch by a factor of
$\sqrt{2}$ applied on four points in the plane. Each point $(x,y)$ was
converted to the point $(x+y,y-x)$ as a result.

###  2.2 The Efficient Solution

These ideas from the previous subsection lead to our main theorem of this
section:

###### Theorem 2.2

Given a planar set of points $P$ of size $n$, all axis-parallel squares
defined by points from $P$ can be enumerated in time $O(n\sqrt{n})$ and $O(n)$
space.

Proof:  Assume w.l.o.g. that all points have positive coordinates, as
otherwise we can shift those without affecting the result. The theorem’s
statement is accomplished using the following algorithm:

Amplified-Squares-Listing$(p_{1},...,p_{n})$:

  1. 1.

Change the representation of each point $p=(x,y)$ to the representation
$(x,y,x+y,y-x)$.
Map each $x$ coordinate in the input to a value in $\\{1,...,n\\}$. Perform a
similar procedure for the $y$ coordinates, the $x+y$ coordinates and the $y-x$
coordinates.

  2. 2.

Build an array $A$ on the input points, sorted by the $x$ coordinate.

  3. 3.

For each pair of post-labeled points $p_{1}=(x,y_{1},x+y_{1},y_{1}-x)$ and
$p_{2}=(x,y_{2},x+y_{2},y_{2}-x)$ with $y_{2}>y_{1}$, out of the first $n$
pairs of points in $A$ that reside in a short column – define the query pair
$q_{1}=(*,y_{1},x+y_{2},*),q_{2}=(*,y_{2},*,y_{1}-x)$ that complements
$p_{1},p_{2}$ to a square from the right, and the query pair
$q_{1}^{\prime}=(*,y_{1},*,y_{2}-x),q_{2}^{\prime}=(*,y_{2},x+y_{1},*)$ that
complements to a square from the left. The wildcards replace the two unknown
coordinates, which are uniquely defined anyway, given the other two.

  4. 4.

Place each query point defined by its $y$ and $x+y$ coordinates, i.e., points
of the form of $q_{1}$ and $q_{2}^{\prime}$ from step 2, in an array $B_{1}$
along with all input points, and apply radix sort on $B_{1}$ based on those
two coordinates. Perform a similar procedure for points of the form of $q_{2}$
and $q_{1}^{\prime}$ in another array $B_{2}$.

  5. 5.

Scan $B_{1}$ and mark each query point adjacent to an input point sharing the
same coordinates, or to an already marked identical query point. Act similarly
on $B_{2}$. Report each square found during this scanning, unless the other
two vertices defining it are on a short column and complement to a square from
the left, to avoid reporting the same square more than once.

  6. 6.

Perform steps $3$-$5$ iteratively on each subsequent $n$ pairs of points in
$A$ that reside in a short column.

  7. 7.

Delete all points that are on short columns from $A$. Convert each remaining
point $(x,y)$ to $(y,x)$. Apply steps $1$-$6$ on the remaining converted
points.

Correctness: Most of the main ideas behind the algorithm’s correctness were
described earlier, in the previous subsection. We now elaborate on the rest of
the details. All squares having at least one edge on a short column are
reported in step 5 of the algorithm, before applying step 6. The rest have
both edges on long columns, so they are reported in step 7.

Given a pair $p_{1}=(x_{1},y_{1})$ and $p_{2}=(x_{1},y_{2})$ with
$y_{2}>y_{1}$, the pair $q_{1}=(x_{2},y_{1})$ and $q_{2}=(x_{2},y_{2})$ that
complement to a square from the right (i.e., $x_{2}>x_{1}$) maintains that
$x_{2}+y_{1}=x_{1}+y_{2}$ since they lie on a diagonal with the same Manhattan
distance from the origin. Putting if differently,

$x_{2}=x_{1}+(y_{2}-y_{1})\Longrightarrow
x_{2}+y_{1}=x_{1}+(y_{2}-y_{1})+y_{1}=x_{1}+y_{2}$

As for $q_{2}$, the latter equality also shows that $x_{2}-y_{2}=x_{1}-y_{1}$
by subtracting $y_{1}+y_{2}$ from both sides. These are exactly the query
points defined by the algorithm, up to the labeling that maintains those
properties. The analysis for the pair that complements to a square from the
left is symmetric.

Complexity: As for the running time, the first two steps of the algorithm cost
$O(n\log n)$ using some standard sorting algorithm, e.g., merging sort. Each
time step 3 is performed, at most $2n$ query points are defined in $O(n)$
time. Each time steps 4-5 are performed, two arrays, each of size at most $n$,
are sorted and then scanned in a linear time. Defining the obtained squares in
step 5, based on the marked queries, can also be carried out in $O(n)$ time.
The total number of query points defined after finishing step 6 is

$O\left(\sum_{i}s_{i}^{2}\right)\leq
O\left(\sum_{i}s_{i}\sqrt{n}\right)=O\left(\sqrt{n}\cdot\sum_{i}s_{i}\right)\leq
O\left(n\sqrt{n}\right)$

where $s_{i}$ denotes the length of the $i$’th short column, i.e., the number
of points in it. As mentioned, each batch of $O(n)$ queries is handled in
$O(n)$ time, so the total running time analysis for steps 1-6 of this
algorithm is $O(n\sqrt{n})$. The analysis for the converted points in step 7
is symmetric. It only remains to notice that the number of pairs, this time,
is

$O\left(\sum_{i}d_{i}^{2}\right)\leq
O\left(\sum_{i}d_{i}\sqrt{n}\right)=O\left(\sqrt{n}\cdot\sum_{i}d_{i}\right)\leq
O\left(n\sqrt{n}\right)$

where $d_{i}$ denotes the length of the $i$’th row out of the remaining rows,
after deleting the short columns. We used the fact that there are at most
$\frac{n}{\sqrt{n}}$ long columns, as otherwise there are more than $n$ input
points. Therefore, the length of each remaining row, after deletion, is at
most $\frac{n}{\sqrt{n}}=\sqrt{n}$, and all points are treated. We remark that
parts of this analysis occur in [15], but we provide those in our proof for
completeness.

As for the space complexity, note that each of the data structures defined in
the above algorithm is of size $O(n)$, and that each step involving those
structures does not cost more than $O(n)$ space.

##  3 The General Case

In this section, we treat the generalized version of the problem from the
previous section. We describe an algorithm that enumerates all scaled copies
of any arbitrary full-dimensional pattern in $d$-dimensions, assuming the
pattern is of a fixed size. Since such copies are determined by the images of
any two points in the pattern, a trivial listing algorithm defined similarly
to the trivial algorithm from the previous section can be easily obtained. We
can safely assume w.l.o.g. that there are two pattern points that share a
coordinate, as otherwise we can rotate the input. The trivial algorithm merely
chooses these two pattern points and calculates in $O(1)$ the transformations
that produce each of the other pattern points. Then it scans each pair of the
background points that share this coordinate, and verifies whether they could
be complemented to that pattern using the computed transformations. The
running time of this algorithm is therefore $O(n^{2}\log n)$.

In [15], the authors managed to generalize their axis-parallel squares listing
algorithm to an algorithm that enumerates all axis-parallel full-dimensional
hypercubes in $d$-dimensions, outputting the list of all hypercubes in time
$O(n^{1+1/d}\log n)$. Their algorithm relies on similar ideas that guided the
treatment of the planar case, namely the partitioning to short axis-parallel
lines (columns, in the planar case), handling those first. In [3], the latter
result was further generalized to any full-dimensional pattern in
$d$-dimensions, providing an algorithm that lists all copies of it, and that
admits the same running time complexity.

###  3.1 Warm Up - Axis-Parallel Hypercubes

We first treat the case of enumerating all full-dimensional axis-parallel
hypercubes. We provide an algorithm with a running time complexity of
$O(n^{1+1/d})$ and with a linear space complexity, addressing the open
question of matching the lower bound from Theorem 1. Our algorithm builds on
the techniques and the observations from the previous section, with some
additions and adjustments so it complies with the properties of the
$d$-dimensional space.

One observation that is true for the $d$-dimensional case, is that any two
points with all but one equal coordinate, determine $2^{d-1}$ full-dimensional
possible hypercubes. Denote such a pair of points by $t$ and $r$. Each of
these hypercubes is uniquely associated with a vector $e\in\\{-1,1\\}^{d-1}$,
in which the $j$’th coordinate determines the direction of progress from $t$
and $r$ along the $j$’th axis, where $j$ is any coordinate except the one in
which they differ. In a similar fashion to the planar case, we would like to
relabel the input coordinates, their sums and their differences, place them in
an array, radix sort it and mark the correct vertices that complement to a
hypercube in an efficient manner, while scanning this array. Moreover, we scan
only pairs of points lying on short axis-parallel lines, similarly to the
planar case, only that this time, by ”short” we mean having not more than
$n^{1/d}$ point on it.

Keeping these ideas in mind, we are ready to present the generalized version
of Theorem 2:

###### Proposition 3.1

Given a set of points $P$ of size $n$ in $d$-dimensional Euclidean space, all
axis-parallel full-dimensional hypercubes defined by points from $P$ can be
enumerated in time $O(n^{1+1/d})$ and $O(n)$ space.

Proof:  Assume w.l.o.g. that all points have positive coordinates, as
otherwise we can shift those without affecting the result. The following
algorithm establishes the proposition’s statement:

Amplified-Hypercubes-Listing$(p_{1},...,p_{n})$:

  1. 1.

For each input point $p=(x_{1},x_{2},...,x_{d})$, add the following additional
list of coordinates:

$\left((x_{i}-x_{j}),(x_{i}+x_{j})~{}~{}|~{}~{}\forall 1\leq i<j\leq d\right)$

Map each of those augmented coordinates, including the original ones, to a
label in $\\{1,...,n\\}$.

  2. 2.

Build an array $A$ on the input points, sorted by each of their coordinates
based on the coordinates’ order, except for the last original coordinate
(i.e., $x_{d}$).

  3. 3.

For each pair of points $t,r$ that lie in the same short axis-parallel line,
having the same coordinates except for the last, out of the first $n$ pairs
with this property, add $2^{d}-2$ query points which define together a
hypercube. Do this for all $2^{d-1}$ possible hypercubes in the following
manner. First, any axis-parallel hypercube having $t$ and $r$ as its vertices,
is defined using one additional vertex

$r^{\prime}=(r_{1}+e_{1}\cdot\delta,...,r_{d-1}+{e_{d-1}}\cdot\delta,t_{d})$

where $r_{i}$ is the $i$’th coordinate in $r$ (and similarly for $t$),
$\delta=r_{d}-t_{d}$ and $e\in\\{-1,1\\}^{d-1}$. The rest of the vertices in
each such hypercube are defined similarly, except for replacing all subsets of
the coordinates in the vector $e$ by zeros, and using $t$ instead of $r$.

Now, define the coordinates that are to be searched – not in the
aforementioned arithmetic manner, but using the labels from step 1 instead.
That is, translate $r_{i}+e_{i}\cdot\delta$ to the label of $r_{i}+r_{d}$ if
$e_{i}=1$, and to that of $r_{i}-r_{d}$ otherwise. Fill in the unknown
coordinates using wildcards, as those are uniquely defined anyway, given the
others.

  4. 4.

Place all query points defined by the same coordinates in an array along with
all input points. Apply radix sort on each of those arrays, according to the
known coordinates in it.

  5. 5.

Scan each array from step 4, and mark each query point adjacent to an input
point sharing the same coordinates, or to an already marked identical query
point. Report all hypercubes that were found (by checking that all vertices
are present for each hypercube), except for hypercubes that have two vertices
on a short axis-parallel line of the same type, only with a smaller index.

  6. 6.

Perform steps 3-5 on each subsequent $n$ pairs of points in $A$ on a short
axis-parallel line of that type.

  7. 7.

Delete all points that are on a short axis-parallel line of a currently
analyzed type from $A$, and convert each remaining point
$(x_{1},x_{2},...,x_{d})$ to $(x_{d},x_{1},...,x_{d-1})$. Apply steps 1-6 on
the remaining converted points. This step is carried out $d-1$ times.

Correctness: Almost all details regarding the analysis of the correctness of
this algorithm already appeared in that of Amplified-Squares-
Listing$(p_{1},...,p_{n})$. As for the phase of searching by labels, note that
if $r_{j}^{\prime}$ is some unknown coordinate, for which we only have an
undesired arithmetic definition based on the coordinate $r_{j}$ and on
$\delta$, then

$r_{j}^{\prime}=r_{j}+\delta=r_{j}+(r_{d}-t_{d})\Longrightarrow
r_{j}^{\prime}+t_{d}=r_{j}+r_{d}$

which exactly corresponds to the labels that the algorithm searches (the
treatment of positive or negative values of $\delta$ is symmetric), and
similarly for subtraction.

Complexity: The running time analysis is similar to that of Amplified-Squares-
Listing$(p_{1},...,p_{n})$ with the following differences. The relabeling in
step 1 and the sorting of $A$ in step 2; the definitions of $2^{d-1}$
hypercubes, each consisting of additional $2^{d}-2$ points in step 3; applying
radix sort on $2^{O(d)}$ arrays in step 4; the linear scanning of those arrays
in step 5 and applying step 7 for $d-1$ times – all of those involve
multiplying the complexity of the planar case by at most a constant factor of
$2^{O(d)}$. The only major difference concerns the total number of query
points defined each time step 7 in the algorithm is invoked. In the treatment
of the first $d-1$ axis-parallel lines, only short lines are considered, so
the cost for each line is

$O\left(\sum_{i}s_{i}^{2}\right)\leq O\left(\sum_{i}s_{i}\cdot
n^{1/d}\right)=O\left(n^{1/d}\cdot\sum_{i}s_{i}\right)\leq
O\left(n^{1+1/d}\right)$

where $s_{i}$ denotes the length of the $i$’th short axis-parallel line of
that form. Regarding the treatment of the last axis-parallel line, we note
that all remaining such lines are short. Assume that there exists a long
remaining line. Then all points on it are on long axis-parallel lines, with
respect to some axis, which induces more points that are on long lines with
respect to another axis, yielding that there are more than
$\left(\left(n^{1/d}\right)\right)^{d}$ input points in total, which is
obviously a contradiction. Thus, all input points are treated with the
mentioned running time, so the total running time analysis is indeed
$O(n^{1+1/d})$. The analysis of the space complexity is similar, with a
constant multiplicative factor of $2^{O(d)}$ compared to the analysis of the
planar case.

###  3.2 General Scaled Copies of a Pattern

At this point, we have all the necessary building blocks to present our final
result of this paper, that deals with enumerating all full-dimensional copies
of a given fixed-size pattern. We show an algorithm for this task that works
in time $O(n^{1+1/d})$ with a linear space complexity. This answers the open
question of realizing the lower bound of Elekes and Erdős [7] regarding the
maximum possible number of occurrences of such patterns, raised in [3]. In the
following result, we rely on an observation from Section 2, by which we can
generate identical labels by means of rotating the input points, which is
equivalent to using the sum (or difference) of the coordinates in the simple
case of hypercubes.

The full statement of our theorem is as follows:

###### Theorem 3.2

Given a fixed set of points $Q$ of full dimension in the $d$-dimensional
Euclidean space, and a set of points $P$ of size $n$, all scaled copies of $Q$
defined by points from $P$ can be enumerated in time $O(n^{1+1/d})$ and $O(n)$
space.

Proof:  Assume w.l.o.g. that all points have positive coordinates, as
otherwise we can shift those without affecting the result. Assume also that
the pattern has more than two points, as the other case may be easily handled.
The following algorithm is used to prove the theorem:

Amplified-Scaled-Copies-Listing$(p_{1},...,p_{n})$:

  1. 1.

Rotate the input points and the pattern points simultaneously, such that two
of the pattern points, $p$ and $q$, lie in the same column afterwards (in
$d$-dimensions – in some axis-parallel line). Calculate the transformations
that produce the other pattern points out of $p$ and $q$.

  2. 2.

For a point $r\neq p,q$ in $Q$, calculate the rotation that places $r$ and $p$
on the same column. Thus, $p$ and $r$ share the same value of some linear
transformation of their coordinates. Label the original input points’
coordinates, along with that linear transformation value for those, using
labels in $\\{1,...,n\\}$ (in $d$-dimensions – perform this labeling procedure
in all dimensions, as the pattern is full-dimensional).

  3. 3.

Build an array $A$ on the input points, sorted by their $x$ coordinate (in
$d$-dimensions – a sorting based on all coordinates except for one).

  4. 4.

For each pair $r$ and $t$ residing in a short column (in $d$-dimensions –
short axis parallel line that corresponds to the sorting from step 3),
calculate all values of the relevant pattern points according to the
transformations from step 1. Express those using the corresponding labels
obtained from step 2, until defining $n$ such query set. If there are other
pairs among the pattern points determining a line parallel to that through $r$
and $t$, let $r$ and $t$ have also the role of that pair when defining the
other points that complement to the pattern, in a similar manner to the
treatment of the hypercubes case.

  5. 5.

Place all query points defined by the same (augmented) coordinates in an array
with all input points. Apply radix sort on each such array.

  6. 6.

Scan each array from step 5, and mark each query point adjacent to an input
point sharing the same coordinates, or to an already marked identical query
point. Report all copies that were found, except for copies that are duplicate
(may be generated by two pairs in the pattern having the same slope).

  7. 7.

Perform steps 4-6 on each subsequent $n$ pairs of points in $A$ on a short
column (in $d$-dimensions – axis-parallel line of that type).

  8. 8.

Delete all points that are on a short column (in $d$-dimensions – axis-
parallel line of a currently analyzed type) from $A$. Apply steps 1-6 on the
remaining points, choosing this time a different pair which is not parallel to
the previous pairs chosen. This step is carried out $d-1$ times.

Analysis: The correctness and the complexity analysis of this algorithm are
essentially the same as those of Theorem 2 and Proposition 3, with some minor
exceptions. The labeling of the points, and searching by those labels, are
both legitimate, as justified earlier. The first two steps cost $O(1)$, as the
pattern is of fixed size. The analysis from the proof of Proposition 3, by
which at each iteration (dimension), the number of treated points is
$O(n^{1+1/d})$ is similar here, as we did not use there the fact that the
lines are axis-parallel, but rather only the fact that the corresponding
vectors form an independent set. Thus, the argument by which there would be
more than $\left(n^{1/d}\right)^{d}$ points, if the last iteration includes a
long axis-parallel line, still holds.

##  4 Conclusion and Further Work

In this paper, we analyzed the problem of efficiently enumerating all axis-
parallel squares within a set of $n$ points in the plane. We provide an
algorithm with a running time of $O(n\sqrt{n})$ for this task, thereby solving
an open question for more than three decades [15] of matching the
corresponding lower bound on the running time of the optimal algorithm for
this task, as the maximum possible number of such squares is
$\Theta(n\sqrt{n})$. We consolidate our result by treating the $d$-dimensional
case of listing all full-dimensional hypercubes in $\mathbb{R}^{d}$ in time
$O(n^{1+1/d})$. Next, we show that our framework can even solve the general
case of reporting all scaled copies of a given pattern of fixed size in time
$O(n^{1+1/d})$, answering an open question from [3] by matching a similar
lower bound due to Elekes and Erdős [7]. All of our algorithms admit a linear
space complexity.

As for the techniques used, we relied on some existing algorithms for those
tasks, amplifying and speeding up their running time using bucket-based
methods (non-comparative sorting, labeling), sweep-line scanning and more. As
far as we are aware of, the combinations of these techniques in the way we
deployed was not noted in the literature so far for dealing with tasks having
a similar flavor.

One immediate open question is whether these can be applied or adjusted for
other pattern matching problems, such as copies obtained from rotating a
patters, combining rotations and scaling and more, solving some of the open
problems in this field. Another possibility is applying them even for
different computational geometric problems. Another questions is whether the
task of finding a single copy of a general pattern is an easier computational
task than enumerating all copies of it. For the case where the pattern is not
of a constant size, but rather input dependent, the best known algorithm [3]
includes a multiplicative factor of the magnitude of the pattern size, as
compared to the combinatorial bound on the maximum possible number of its
copies, that does not include that factor. Thus, also shaving this factor from
the computational results becomes another important open question for further
research.

#### Acknowledgements

We would to thank the students in the Algorithmic Problem Solving class
(2019/2020) at the Hebrew University of Jerusalem for bringing to our
attention that this was a known open problem in the field of computational
geometry.

## References

  * [1] Helmut Alt and Leonidas Guibas.  Discrete geometric shapes: Matching, interpolation, and approximation: A survey.  Handbook of Computational Geometry, 03 1997.
  * [2] Marc Alzina, Wojciech Szpankowski, and Ananth Grama.  2d-pattern matching image and video compression: Theory, algorithms, and experiments.  IEEE transactions on image processing : a publication of the IEEE Signal Processing Society, 11:318–31, 02 2002.
  * [3] P. Braß.  Combinatorial geometry problems in pattern recognition.  Discrete & Computational Geometry, 28:495–510, 2002.
  * [4] P. Braß and J. Pach.  Chapter 2 problems and results on geometric patterns.  2005\.
  * [5] Minkyoung Cho and David Mount.  Improved approximation bounds for planar point pattern matching.  Algorithmica, 50:175–207, 02 2008.
  * [6] Pedro De Rezende and D. Lee.  Point set pattern matching in d-dimensions.  Algorithmica, 13:387–404, 04 1995.
  * [7] György Elekes and Paul Erdos.  Similar configurations and pseudo grids.  Intuitive geometry, 63:85–104, 1994.
  * [8] P. Erdös.  On sets of distances of n points.  American Mathematical Monthly, 53:248–250, 1946.
  * [9] Paul Finn, Lydia Kavraki, Jean-Claude Latombe, Rajeev Motwani, Christian Shelton, Suresh Venkatasubramanian, and A. Yao.  Rapid: Randomized pharmacophore identification for drug design.  volume 10, pages 324–333, 01 1997.
  * [10] Michael L Fredman and Dan E Willard.  Blasting through the information theoretic barrier with fusion trees.  In Proceedings of the twenty-second annual ACM symposium on Theory of Computing, pages 1–7, 1990.
  * [11] David M Mount, Nathan S Netanyahu, and Jacqueline Le Moigne.  Efficient algorithms for robust feature matching.  Pattern Recognition, 32(1):17–38, 1999.
  * [12] Raquel Norel, Daniel Fischer, Haim Wolfson, and Ruth Nussinov.  Molecular surface-recognition by a computer vision-based technique.  Protein engineering, 7:39–46, 02 1994.
  * [13] Grant Schindler, Panchapagesan Krishnamurthy, Roberto Lublinerman, Yanxi Liu, and Frank Dellaert.  Detecting and matching repeated patterns for automatic geo-tagging in urban environments.  pages 1–7, 07 2008.
  * [14] Meenakshi Sharma.  Analysis of distance measures in content based image retrieval.  09 2018.
  * [15] Marc J. van Kreveld and Mark T. de Berg.  Finding squares and rectangles in sets of points.  In Manfred Nagl, editor, Graph-Theoretic Concepts in Computer Science, pages 341–355, Berlin, Heidelberg, 1990. Springer Berlin Heidelberg.
  * [16] Vasileios Zografos and Bernard Buxton.  Affine invariant, model-based object recognition using robust metrics and bayesian statistics.  volume 3656, 09 2005.