11institutetext: Department of Physics, National University of Singapore

# Photometric Redshift Analysis using Supervised Learning Algorithms and Deep
Learning

Chong De Wei, Kenny 11<EMAIL_ADDRESS>Abel Yang 11<EMAIL_ADDRESS>
###### Abstract

We present a catalogue of galaxy photometric redshifts for the Sloan Digital
Sky Survey (SDSS) Data Release 12. We use various supervised learning
algorithms to calculate redshifts using photometric attributes on a
spectroscopic training set. Two training sets are analysed in this paper. The
first training set consists of 995,498 galaxies with redshifts up to $z\approx
0.8$. On the first training set, we achieve a cost function of 0.00501 and a
root mean squared error value of 0.0707 using the XGBoost algorithm. We
achieved an outlier rate of 2.1% and 86.81%, 95.83%, 97.90% of our data points
lie within one, two, and three standard deviation of the mean respectively.
The second training set consists of 163,140 galaxies with redshifts up to
$z\approx 0.2$ and is merged with the Galaxy Zoo 2 full catalog. We also
experimented on convolutional neural networks to predict five morphological
features (Smooth, Features/Disk, Star, Edge-on, Spiral). We achieve a root
mean squared error of 0.117 when validated against an unseen dataset with over
200 epochs. Morphological features from the Galaxy Zoo, trained with
photometric features are found to consistently improve the accuracy of
photometric redshifts.

##  1 Overview and Methodology

Redshifts of celestial objects have been a vital component in the field of
astronomy and we use them to measure various attributes such as the rotation
of the galaxy and the distance from us. Traditionally, they have been measured
by spectroscopy. While spectroscopy is effective in determining redshifts of
galaxies, it is time consuming and expensive and therefore not scalable to map
a spectroscopic redshift for every celestial object. Spectroscopy data is
limited in the SDSS database with only 0.36% mapped 2011AJ….142…72E . There
are two methods in the literature currently deployed to determine photometric
redshifts: template-fitting and the statistical/machine learning method. We
will focus on the machine learning method and also use deep learning methods
to increase predictive power.

The SDSS is a major multi-filter imaging and spectroscopic redshift survey
using a 2.5m wide angle optical telescope. We queried out a data size of
995,498 galaxies with redshifts up to $z\approx 0.8$ for our training and
validation set. We first split the data points and photometric features using
the random subspace method. For each subset, we train them on an individual
tree and use standard deviation reduction to determine the quality of every
split. The process is repeated when building the tree down to the leaf node.
We use this method to find the more important features.

At the same time, understanding the physics behind the various photometric
features allow us to feature engineer attributes that perform stronger than
the native features. In another paper investigating photometric redshifts
2008ApJ…674..768O , concentration indices are found to break the degeneracy
between colour band and photometric redshifts. Indeed, concentration indices
contain some information on morphological classification and concentration
indices separate galaxies into early and late types 2003AJ….125.1682N , with
the early type galaxies having lower values of the concentration indices. This
motivates the need to include the concentration indices into the neural
network.

Since redshifts are determined by looking at transition lines and these lines
form very sharp peaks at certain wavelengths, we are able to extract some
photometric information on these peaks by taking the difference between the
colour bands. The colour bands and colour indexes interact with each other and
we use the multilayer perceptron to encourage these interactions.

For each neural network, we perform a 5-fold validation and the same neural
network is trained 5 times. The average results of the cost function are then
recorded.

##  2 Results

Model Features | Cost Function
---|---
modelmag color band | 0.0013276
(u, g, r, i, z) |
dered color band | 0.0013093
(u, g, r, i, z) |
dered color indexes | 0.0014040
(u-g, g-r, r-i, i-z) |
dered color band + dered color indexes | 0.0012729
dered color band + dered color indexes + concentration indices | 0.0012559
Table 1: Cost Function of Neural Networks with various features

As shown in Table 1, the best model consists of 14 features with a cost
function of 0.0012559 with the second training set. The main dataset is then
trained with the best model features across 4 algorithms (Linear Regression,
Random Forests, XGBoost and Neural Networks [MLP]) and we achieve a cost
function of 0.00501 using the XGBoost model Chen:2016:XST:2939672.2939785 .

In Figure 1, we plot photometric determined redshifts versus spectroscopic
redshifts in a gaussian density plot and the average redshift bias is
determined to be $2.59\times 10^{-3}$ with a standard deviation of 0.0559. As
seen in Figure 1, photometric redshift precision is found to decrease as the
wavelength increases. The cost function is found to reduce to 0.00183 from
0.00501 for $r<20$ spectroscopic objects whereas the cost function increase to
0.0106 for $r>20$ spectroscopic objects.

Figure 1: XGBoost model of first training set

Outliers are defined to be larger than three standard deviations of the mean.
We achieved an outlier rate of 2.10% and calculated 86.81%, 95.83%, 97.90% of
our data points lie within one, two, and three standard deviation of the mean
respectively.

##  3 Creating Scalable Features beyond SDSS photometric data

###  3.1 Performance of Galaxy Zoo 2 Features

One major limitation with current photometric features is the degeneracy
between spectroscopic determined redshifts and 5 colour SDSS photometric data.
There is a considerable number of celestial objects with the same colour band
range yet they give rise to different redshift values. No matter what
supervised learning algorithm you use, it won’t be able to differentiate these
galaxies apart.

We attempt to use morphological features to break the degeneracy between
redshift and colour band. The Galaxy Zoo 2013MNRAS.435.2835W  is a crowd
sourcing project where online citizens can contribute by classifying galaxies
through a decision tree process into various categories. With a sufficient
sample size, the end result takes on a value that is analogous to probability,
that is if 60% of online citizens voted a spiral pattern, the probability of a
galaxy being a spiral is 0.6. Together with features from the SDSS, we merge
the two databases together and deploy the same methodology to determine the
performance of Galaxy Zoo features.

Model | Cost Function
---|---
Models without Galaxy Zoo 2 Features | 0.001256
Models with Galaxy Zoo 2 Features | 0.001198
Table 2: Performance of Galaxy Zoo 2 Features

The selected Galaxy Zoo 2 features are the probability of galaxy being smooth,
the probability of galaxy having features/disk, probability of galaxy,
probability of galaxy being an edge-on galaxy and probability of galaxy having
a spiral. Similarly, we train neural networks with the same subset of features
in Table 1 together with Galaxy Zoo 2 features and determine their performance
using the cost function. Neural networks with Galaxy Zoo 2 features are found
to consistently outperform neural networks without Galaxy Zoo 2 features. This
demonstrates these five Galaxy Zoo 2 features provided morphological
information to break some degeneracy between redshift and information from the
SDSS.

###  3.2 Deep Learning

The objective of this component is to build a model that predicts
morphological shapes labelled by online citizens. We use a convolutional
neural network architecture to train a multi-label image classifier. We down
sample the input into 3 normalized channels (RGB) with $150\times 150$ inputs.
After several iterations, the final model consists of a total of 10 layers: 8
convolution layers, 4 max pooling layers, 2 full connected (with dropout)
layers. In all convolution layers, we apply the ‘same’ padding to retain the
output shape of the matrix. A ‘ReLu’ activation is applied, and the number of
kernels used is 32, 32, 64, 128 in that order. Every two convolution layers is
followed by a max pooling layer of size ($2\times 2$). The output is then
flattened into a vector and then fed into the fully connected layer with 1024
nodes, followed by an activation function. 50% of the features are then
dropped out. Similar to the MLP, the cost function used is the mean squared
error and the ‘adam’ optimizer 2014arXiv1412.6980K  is used with learning rate
= 0.001, beta1 = 0.9, beta2 = 0.999, fuzz factor = $10^{-8}$ and learning rate
decay = 0.

We achieve a root mean squared error of 0.117 when validated against an unseen
dataset with over 200 epochs. In addition, we can also map other morphological
classifications into the model as well, and deploy transfer learning or use a
more extensive architecture to predict the morphological features of any
arbitrary galaxy to a higher accuracy. Deep learning algorithms allow us to
scale morphological features without the need of manual labelling and we can
use the morphological features from the galaxy zoo to determine photometric
redshifts to a higher accuracy.

## References

  * (1) D.J. Eisenstein, D.H. Weinberg, E. Agol, H. Aihara, C. Allende Prieto, S.F. Anderson, J.A. Arns, É. Aubourg, S. Bailey, E. Balbinot et al., Astron. J. 142, 72 (2011), 1101.1529
  * (2) H. Oyaizu, M. Lima, C.E. Cunha, H. Lin, J. Frieman, E.S. Sheldon, Astrophys. J. 674, 768 (2008), 0708.0030
  * (3) O. Nakamura, M. Fukugita, N. Yasuda, J. Loveday, J. Brinkmann, D.P. Schneider, K. Shimasaku, M. SubbaRao, Astron. J. 125, 1682 (2003), astro-ph/0212405
  * (4) T. Chen, C. Guestrin, _XGBoost: A Scalable Tree Boosting System_ , in _Proceedings of the 22Nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ (ACM, New York, NY, USA, 2016), KDD ’16, pp. 785–794, ISBN 978-1-4503-4232-2, 1603.02754
  * (5) K.W. Willett, C.J. Lintott, S.P. Bamford, K.L. Masters, B.D. Simmons, K.R.V. Casteels, E.M. Edmondson, L.F. Fortson, S. Kaviraj, W.C. Keel et al., MNRAS 435, 2835 (2013), 1308.3496
  * (6) D.P. Kingma, J. Ba, ArXiv e-prints (2014), 1412.6980


# Impact of dark states on the stationary properties of quantum particles with
off-centered interactions in one dimension

G. Bougas  Center for Optical Quantum Technologies, Department of Physics,
University of Hamburg, Luruper Chaussee 149, 22761 Hamburg Germany  N. L.
Harshman  Physics Department, American University, Washington, DC 20016, USA
P. Schmelcher  Center for Optical Quantum Technologies, Department of Physics,
University of Hamburg, Luruper Chaussee 149, 22761 Hamburg Germany  The
Hamburg Centre for Ultrafast Imaging, University of Hamburg, Luruper Chaussee
149, 22761 Hamburg, Germany

###### Abstract

We present a generalization of the two-body contact interaction for non-
relativistic particles trapped in one dimension. The particles interact only
when they are a distance $c$ apart. The competition of the interaction length
scale with the oscillator length leads to three regimes identified from the
energy spectra. When $c$ is less than the oscillator length, particles avoid
each other, whereas in the opposite case bunching occurs. In the intermediate
region where the oscillator length is comparable to $c$, both exclusion and
bunching are manifested. All of these regions are separated by dark states,
i.e. bosonic or fermionic states which are not affected by the interactions.

##  I Introduction

A paradigmatic model for understanding interacting quantum systems is provided
by the two-body contact potential. This zero-range interaction occurs
naturally in effective field theories for bosons and multi-component fermions,
and it has proven remarkably useful for describing the physics of trapped
ultracold atomic gases, where the typical length scales are all much larger
than the effective range of the potential Pethick and Smith (2008); Bloch  _et
al._ (2008). The experimental control possible for ultracold atoms in
effectively one dimensional (1D) traps has been remarkably productive for
studying the dynamics of single and multi-species quantum gases Mistakidis
_et al._ (2023); Wenz  _et al._ (2013); Serwane  _et al._ (2011); Sowiński and
Ángel García-March (2019), in part because the theoretical description of 1D
contact interactions does not require regularization or renormalization to
achieve physically meaningful results Busch  _et al._ (1998); Farrell and van
Zyl (2009).

For single-component bosons, the $N$-body Hamiltonian with two-boy contact
potential

$H=\sum_{j=1}^{N}\left(-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{\partial
x_{j}^{2}}+V(x_{j})\right)+g\sum_{j<k}\delta(x_{j}-x_{k})$ (1)

has been exhaustively studied for various trapping potentials $V(x)$. In the
case of an infinite square well or no trapping potential, the model (1) is the
Lieb-Liniger model Lieb and Liniger (1963) and solvable in the bosonic sector
by the Bethe-ansatz for any interaction strength $g$ Gaudin (2014). Similarly,
for $N=2$ and a harmonic trapping potential the model is also solvable for any
$g$ Busch  _et al._ (1998). The model (1) can also be generalized to
multicomponent boson and fermion models Sowiński and Ángel García-March
(2019); García-March  _et al._ (2014); Harshman (2016a, b), for which the
previous special cases are also solvable Sutherland (2004).

One of the most remarkable results (valid for any trapping potential) is the
Bose-Fermi mapping due to Girardeau Girardeau (1960): in the hard core limit
$g\to\infty$ and for any trapping potential, bosonic solutions of (1) lie in
one-to-one correspondence with non-interacting fermionic solutions, sharing
the same energy and spatial probability density. As the interaction strength
is tuned from $g=0$ to $g\to\infty$, each bosonic energy level shifts to the
corresponding fermionic level. The contact interaction at $g\to\infty$
excludes the bosons from two-body coincidences, and the bosons are said to be
‘fermionized’.

A key feature of the contact interaction is that single-component fermions are
not at all impacted by the interaction, i.e. they are ‘dark’ to the contact
potential. Single-component fermionic wave functions are totally antisymmetric
under exchange, and this antisymmetry forces nodes in the wave function
precisely at the two-body coincidences where the contact interaction has
support. Note that in higher dimensions, two single-component bosons are also
dark to contact interactions when they have non-zero relative angular
momentum, which also forces a node at the two-body coincidences. In other
words, the contact interaction creates scattering only in the $s$-wave
channel. However in 1D, relative angular momentum is the same as relative
parity for two particles, and there are only two two-body channels: even or
odd relative parity. Odd relative parity coincides with fermionic antisymmetry
in 1D and even relative parity with bosonic symmetry (see discussion of
symmetry below), and so the contact interaction is felt by states with even
relative parity and can mimic statistical exclusion.

Several proposals to modify the contact interaction already exist. The most
general form for a point defect in one dimension gives a three parameter set
of contact interactions Albeverio  _et al._ (2012). This set includes so-
called $p$-wave interactions that single-component fermions feel but which are
dark to single-component bosons Cheon and Shigehara (1998, 1999).
Alternatively, if the two-body interaction range is taken to be a Gaussian or
some other regularized function that limits to a Dirac delta, then there is a
finite range to the interaction, and fermions are again no longer dark to the
interaction.

The goal of this article is to consider a simple generalization of the contact
interaction in 1D. The interaction we propose introduces a range to the two-
body interaction in a different way: two particles interact only when they are
exactly a distance $c$ apart. The two-body interaction describing such a
system is given by

$U(g,c)=g\delta(x_{1}-x_{2}+c)+g\delta(x_{1}-x_{2}-c).$ (2)

The two delta functions guarantee that particle exchange symmetry
$x_{1}\leftrightarrow x_{2}$ holds, and therefore states can be classified as
bosonic or fermionic. In the limit that $c\to 0$, $U(g,c)$ becomes the
standard contact interaction again, but for finite $c$ it introduces a new
length scale into the physical system that can compete with the other length
scales.

We consider the specific case of two particles in a harmonic trap
$H(g,c)=H_{0}+U(g,c)$, where $U(g,c)$ is defined above and

$H_{0}=-\frac{\hbar^{2}}{2m}\left(\partial_{x_{1}}^{2}+\partial_{x_{2}}^{2}\right)+\frac{1}{2}m\omega^{2}(x_{1}^{2}+x_{2}^{2}).$
(3)

The trap introduces another length scale to the problem, the harmonic
oscillator length $a_{\rm{ho}}=\sqrt{\hbar/(m\omega)}$. We show below how the
competition between the length scale $c$ and $a_{\rm{ho}}$ separates the
energy spectrum into three regimes: (1) the _exclusion regime_ for
$c<a_{\rm{ho}}$ and low energies, where the model behaves similar to the
contact interaction model, including two-body exclusion of particles for large
$g$; (2) the _truncation regime_ for $c>a_{\rm{ho}}$ and low energy, where the
interaction length scale $c$ is so large that it suppresses the tail of the
wave functions and creates a bunching effect; and (3) the _crossover regime_
for $c\sim a_{\rm{ho}}$ where exclusion and bunching compete, and small
variations of $c$ can lead to dramatic changes in the wave function variance.
The interfaces between these regions are defined by the appearance of dark
states to the interaction $U(g,c)$. Unlike the contact interaction, these dark
states are both fermionic and bosonic and are sporadically distributed
throughout the spectrum for specific values of $c$ corresponding to zeros of
Hermite polynomials. At these points triple degeneracy occurs at infinite
interactions, where two bosonic (fermionic) states cluster with another
fermionic (bosonic) one. Moreover, dark states determine the competition
between bunching and exclusion in the crossover regime.

The outline is as follows: in Sect. II, we analyze the symmetries of the model
and identify the bosonic and fermionic sectors. In Sect. III, we provide
analytic solutions for arbitrary $g$ and $c$ and analyze how the spectrum
varies with the parameters. In Sect. IV, we construct the solutions at
infinite interactions, analyze the corresponding energy spectrum, and
subsequently discuss the ensuing symmetries. Finally, in Sect. V we summarize
our results and examine possible realizations of the presented model.

##  II Structure and Symmetries of the Model

As a first step towards solving the model $H(g,c)$, we make a transformation
to center-of-mass and relative coordinates,

$\displaystyle X$ $\displaystyle=$ $\displaystyle\frac{1}{2}(x_{1}+x_{2})$
$\displaystyle x$ $\displaystyle=$ $\displaystyle x_{1}-x_{2},$ (4)

so that the Hamiltonian separates as
$H(g,c)=H_{\mathrm{com}}+H_{\mathrm{rel}}(g,c)$:

$\displaystyle H_{\mathrm{com}}$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{4m}\partial_{X}^{2}+m\omega^{2}X^{2}$
$\displaystyle H_{\mathrm{rel}}(g,c)$ $\displaystyle=$
$\displaystyle-\frac{\hbar^{2}}{m}\partial_{x}^{2}+\frac{m\omega^{2}}{4}x^{2}+g\delta(x+c)+g\delta(x-c).$

The model (3) therefore is equivalent to two one-dimensional quantum systems.
It is trivially integrable, with the role of the integral of motion for each
degree of freedom played by each sub-Hamiltonian energy. Another consequence
is that, like all 1D quantum systems with non-singular potentials, the
spectrum of each sub-Hamiltonian should be singly-degenerate except in the
case $g\to\infty$ Loudon (1959); Harshman (2017a). In what follows, we employ
harmonic oscillator units such that $\hbar$, $m$, and $\omega$ all equal $1$.

The center-of-mass sub-Hamiltonian $H_{\mathrm{com}}$ is the familiar 1D
harmonic oscillator, independent of interaction strength $g$ and interaction
displacement $c$. The energy eigenstates have wave functions

$\Phi_{N}(X)=\left(\frac{2}{\pi}\right)^{\frac{1}{4}}\frac{1}{\sqrt{2^{N}N!}}H_{N}(\sqrt{2}X)e^{-X^{2}},$
(6)

where $H_{N}(\sqrt{2}X)$ is the $N$-order Hermite polynomial. The center-of-
mass separates out and therefore we ignore it for the rest of the calculation.

For the relative sub-Hamiltonian $H_{\mathrm{rel}}(g,c)$, we denote the
eigenstates as $\phi_{n}(x)$ with eigenenergies $\epsilon_{n}$. The quantum
number $n\in\\{0,1,2,\ldots\\}$ counts the number of nodes (for finite $g$).
In the special case $g=0$, the solutions are the non-interacting relative wave
functions

$\phi^{0}_{n}(x)=\left(\frac{1}{2\pi}\right)^{\frac{1}{4}}\frac{1}{\sqrt{2^{n}n!}}H_{n}(x/\sqrt{2})e^{-x^{2}/4},$
(7)

with energy $(n+1/2)$. As we show in the next section, solutions for arbitrary
values of $g$ and $c$ can be found analytically by solving the transcendental
equation given by matching boundary conditions for parabolic cylindrical
functions at $x=\pm c$.

Figure 1: (a) Energy levels $\epsilon_{n}$ corresponding to even
($n=0,2,\ldots$) and odd ($n=1,3,\ldots$) parity eigenstates for $c=0.75$. The
panels on the right depict the wavefunctions of the ground bosonic (b) and
fermionic state (c) at various interaction strengths, according to the legend.

Before constructing these solutions, we first analyze the kinematic symmetries
of $H(g,c)$, by which we mean the group of symmetry transformations that
commute with the Hamiltonian $H(g,c)$ for any value of $g$ and $c$. Two
symmetries always hold: spatial inversion $\Pi$ and particle permutation
$\Sigma$. These transform the particle coordinates $(x_{1},x_{2})$ and the
COM-relative coordinates $(X,x)$ in the following manner:

$\displaystyle\Pi\cdot(x_{1},x_{2})\to(-x_{1},-x_{2}),$
$\displaystyle\Pi\cdot(X,x)\to(-X,-x)\ $ (8) and
$\displaystyle\Sigma\cdot(x_{1},x_{2})\to(x_{2},x_{1}),$
$\displaystyle\Sigma\cdot(X,x)\to(X,-x).$ (9)

The combined symmetry group $G=\\{e,\Pi,\Sigma,\Pi\Sigma\\}$ generated by
$\Pi$ and $\Sigma$ is isomorphic to the Klein four-group $G\sim V_{4}$. This
group is Abelian and has four one-dimensional unitary irreducible
representations corresponding to the possible signs $\pm 1$ representing $\Pi$
and $\Sigma$.

Defining $\hat{\Pi}$ and $\hat{\Sigma}$ as operators that represent the
generators $\Pi$ and $\Sigma$ acting on the Hilbert space of wave functions,
the states $\Phi_{N}(X)$ and $\phi_{n}(x)$ transform as:

$\displaystyle\hat{\Pi}\Phi_{N}(X)=(-1)^{N}\Phi_{N}(X);$
$\displaystyle\hat{\Sigma}\Phi_{N}(X)=\Phi_{N}(X)$
$\displaystyle\hat{\Pi}\phi_{n}(x)=(-1)^{n}\phi_{n}(x);$
$\displaystyle\hat{\Sigma}\phi_{n}(x)=(-1)^{n}\phi_{n}(x).$

The center-of-mass wave functions are invariant under particle exchange, but
the relative wave functions can be split into two sectors: bosonic wave
functions that are even under exchange and fermionic wave functions that are
odd under exchange. For the case of two particles, we also see that particle
exchange eigenstates coincide with relative parity eigenstates: bosonic wave
functions are necessarily even under (relative) parity and fermionic are odd
under relative parity. These sectors decouple and can be considered
separately.

For the sake of completeness, we note that the Hamiltonian $H_{0}+U(g,c)$ has
additional symmetries in the limits $g=0$ and $g\to\infty$. When $g=0$, the
Hamiltonian $H_{0}$ is the familiar two-dimensional isotropic harmonic
oscillator. Its _kinematic symmetry_ (i.e, the group of all symmetry
transformations that commute with the Hamiltonian) is given by
$\mathrm{U}(2)$. The relevance of this group can be understood in two ways: 1)
As the set of all unitary transformations of the ladder operators $a_{1}$ and
$a_{2}$ into $a_{1}^{\prime}=u_{11}a_{1}+u_{12}a_{2}$ and
$a_{2}^{\prime}=u_{21}a_{1}+u_{22}a_{2}$, or 2) as the set of all orthogonal
and symplectic transformations of four-dimensional classical phase space,
where the intersection
$\mathrm{O}(4)\cap\mathrm{Sp}(4,\mathbb{R})\sim\mathrm{U}(2)$. Either of these
sets of transformations are equivalent and preserve $H_{0}$. The dimensions of
the irreducible representations of $\mathrm{U}(2)$ lie in one-to-one
correspondence with the degeneracies of energy levels of $H_{0}$ Louck (1965).
Further, the uniform spacing between adjacent energy levels of $H_{0}$ can be
understood by considering the _dynamical symmetry_ group, i.e., the group of
all space-time transformations that preserve the action of the 2D harmonic
oscillator. This group is called the harmonic oscillator group
$\mathrm{HO}(2)$ and includes the kinematic symmetry group $\mathrm{U}(2)$ as
a subgroup Niederer (1973).

The symmetries in the case $g\to\infty$ will be discussed below in Sect. IV.

##  III Generic Interaction Strength

In this section, we provide expressions for the solutions $\phi_{n}(x)$ of the
relative sub-Hamiltonian $H_{\mathrm{rel}}(g,c)$ when the interaction strength
$g$ is finite. We will also identify two critical features of the model: $(i)$
the role played by dark states, i.e., states of the non-interacting
Hamiltonian whose wave function nodes coincide with the interaction
displacement $c$ and thus remain unperturbed; and $(ii)$ the role played by
the competition between the harmonic oscillator length scale and the
interaction offset scale $c$.

###  III.1 Solutions of the Relative Sub-Hamiltonian

The first step towards finding the eigenspectrum of the relative sub-
Hamiltonian Eq. (II) for generic $g$ is to divide the relative $x$ coordinate
into three regions $I=(-\infty,-c)$, $I\\!I=(-c,c)$ and
$I\\!I\\!I=(c,\infty)$. The solutions in these intervals are known: they are
linear combinations of parabolic cylinder functions $D_{Q}(x)$ Abramowitz and
Stegun (1948); Aouadi  _et al._ (2016); Avakian  _et al._ (1987),

$\displaystyle\phi_{n}(x)=\begin{cases}\phi^{(I)}_{n}(x),&x\in I\\\
\phi^{(I\\!I)}_{n}(x),&x\in I\\!I\\\ \phi^{(I\\!I\\!I)}_{n}(x),&x\in
I\\!I\\!I\end{cases},$
$\displaystyle\quad\phi^{(i)}_{n}(x)=\alpha_{i}D_{Q_{n}}(x)+\beta_{i}D_{Q_{n}}(-x),\quad
i=I,I\\!I,I\\!I\\!I,$ (11)

where $Q_{n}\equiv\epsilon_{n}-1/2$ and the coefficients are yet to be
determined. In the case where $Q_{n}=n$, an integer, these functions reduce to
the $n$-order Hermite polynomials Abramowitz and Stegun (1948),

$D_{n}(x)=2^{-n/2}H_{n}(x/\sqrt{2})e^{-1/4x^{2}}=(2\pi)^{1/4}\sqrt{n!}\,\phi^{0}_{n}(x).$
(12)

Note that since we are interested in bound state solutions due to the harmonic
trap, we demand that the eigenstates vanish at $x\to\pm\infty$. From the
asymptotic expansions Abramowitz and Stegun (1948); Aouadi  _et al._ (2016)
$D_{Q_{n}}(x\to\infty)\sim e^{-x^{2}/4}x^{Q_{n}}$ and
$D_{Q_{n}}(x\to-\infty)\sim\sqrt{2\pi}/\Gamma(-Q_{n})e^{x^{2}/4}x^{-1-Q_{n}}$
we infer that $\alpha_{I}=\beta_{I\\!I\\!I}=0$.

Figure 2: (a) Energy levels $\epsilon_{n}$ corresponding to even
($n=0,2,\ldots$) and odd ($n=1,3,\ldots$) parity eigenstates for $c=1.5$. The
panels on the right depict the wavefunctions of the ground bosonic (b) and
fermionic eigenstate (c) at various interaction strengths, according to the
legend.

To determine the remaining coefficients and the energy levels, boundary
conditions are imposed at the interval intersections, namely at $x=\pm c$.
These consist of continuity conditions for the wave functions and
discontinuity conditions for the first derivatives due to the delta
interaction potentials Belloni and Robinett (2014),

$\displaystyle\phi^{(I\\!I)}_{n}(-c)=\phi^{(I)}_{n}(-c),$ (13a)
$\displaystyle\phi^{(I\\!I)}_{n}(c)=\phi^{(I\\!I\\!I)}_{n}(c),$ (13b)
$\displaystyle\frac{d\phi^{(I\\!I)}_{n}(x)}{dx}\Big{|}_{x=-c}-\frac{d\phi^{(I)}_{n}(x)}{dx}\Big{|}_{x=-c}=g\phi^{(I)}_{n}(-c),$
(13c)
$\displaystyle\frac{d\phi^{(I\\!I\\!I)}_{n}(x)}{dx}\Big{|}_{x=c}-\frac{d\phi^{(I\\!I)}_{n}(x)}{dx}\Big{|}_{x=c}=g\phi^{(I\\!I\\!I)}_{n}(c).$
(13d)

Apart from the above boundary conditions, one needs to take into account the
symmetry of the relative sub-Hamiltonian under particle exchange, i.e. the
$\Sigma$ operation. The wave functions can be either even or odd upon particle
exchange, $\hat{\Sigma}\phi_{n}(x)=(-1)^{n}\phi_{n}(x)$, describing identical
bosons or fermions respectively. Therefore, an additional condition can be
imposed for the coefficients, namely $\alpha_{I\\!I\\!I}=(-1)^{n}\beta_{I}$.

Having at hand the above boundary conditions and constraints one can determine
the energy levels and coefficients up to a normalization constant. The
discontinuity conditions (13c), (13d) can be simplified by employing the
recurrence relations that the parabolic cylinder functions enjoy Abramowitz
and Stegun (1948),

$\frac{dD_{Q_{n}}(\pm x)}{dx}=\frac{1}{2}xD_{Q_{n}}(\pm x)\mp D_{Q_{n}+1}(\pm
x).$ (14)

Applying these conditions and properties, the following transcendental
equation for the energy levels $\epsilon_{n}$ is established,

$\displaystyle\frac{D_{Q_{n}}(c)D_{Q_{n}+1}(-c)}{D_{Q_{n}}(-c)+(-1)^{n}D_{Q_{n}}(c)}+\frac{D_{Q_{n}}(-c)D_{Q_{n}+1}(c)}{D_{Q_{n}}(-c)+(-1)^{n}D_{Q_{n}}(c)}$
$\displaystyle=-gD_{Q_{n}}(c).$ (15)

Note that the above relation is general and the energy levels
$\epsilon_{n}=Q_{n}+1/2$ are determined for arbitrary finite interaction
strength $g$. The equation holds also for arbitrary displacement $c$ except
from the special cases where $D_{Q_{n}}(\pm c)=0$, which are treated
separately. These occur whenever $Q_{n}$ is an integer, and the zeros
correspond to the roots of Hermite polynomials [Eq. (12)], and hence to non-
interacting eigenstates $\phi^{0}_{n}(x)$. These states do not ‘feel’ the
interaction, i.e. they are dark states. This behavior is similar to that of
fermionic states, which do not ‘feel’ a contact interaction; see below. We
will see evidence of these dark states in the spectrum below.

Apart from the energy levels, the coefficients are also determined from the
boundary conditions. An analytical expression is thus established for the
eigenstates,

$\displaystyle\phi_{n}(x)=\begin{cases}\beta_{I}D_{Q_{n}}(-x),&x\in I\\\ \\\
\beta_{I}\frac{D_{Q_{n}}(c)\left[D_{Q_{n}}(x)+(-1)^{n}D_{Q_{n}}(-x)\right]}{D_{Q_{n}}(-c)+(-1)^{n}D_{Q_{n}}(c)},&x\in
I\\!I\\\ \\\ (-1)^{n}\beta_{I}D_{Q_{n}}(x),&x\in I\\!I\\!I\end{cases},$ (16)

where $\beta_{I}$ is a normalization constant. As predicted from symmetry
analysis above, the description of bosons (fermions) requires even (odd)
quantum number $n$. Note that the eigenstates related to a particular relative
parity form an orthonormal basis.

###  III.2 Competition between scales

Now we analyze how the competition between the oscillator length scale ($1$ in
our units) and the interaction displacement $c$ is imprinted on the energy
levels and eigenstates. The dark states described above also leave evidence in
the spectrum.

First, consider the displacement $c=0.75<1$, smaller than the trap length
scale. The corresponding energy level structure is presented in Fig. 1. As
expected due to the non-zero displacement, fermionic energy levels are also
affected by the interactions [dashed lines in Fig. 1(a)], contrary to the case
of zero-range interactions Busch  _et al._ (1998); Budewig  _et al._ (2019).
There is however one energy level ($n=4$) which barely depends on the
interaction strength, i.e., a dark state. This occurs because the first root
of the corresponding Hermite polynomial $H_{4}(x/\sqrt{2})$ lies very close to
the displacement $c$ Salzer  _et al._ (1952). As $g\to\infty$, the energy
levels of the adjacent bosonic ($n=2$) and fermionic ($n=3$) state tend
asymptotically to the energy of the dark state, signalling a triple
degeneracy. The latter is present at attractive interactions as well, where
the $n=4$ dark state now clusters with the $n=5,n=6$ levels.

Except from the triple degeneracy induced by the dark state, all other low
energy levels are doubly degenerate as $g\to\infty$. Adjacent bosonic and
fermionic energy levels cluster together in a way reminiscent of the Tonks-
Girardeau regime, manifested in systems with strongly repulsive short-range
interactions Tonks (1936); Girardeau (1960); Kehrberger  _et al._ (2018). In
that regime, bosons turn into hardcore particles avoiding each other and their
energy levels become degenerate with those of non-interacting fermions. By
analogy, we call this displacement parameter range ($c<1$) the _exclusion
regime_.

The exclusion is demonstrated in the transition of the profiles of the first
bosonic and fermionic eigenstates ($n=0,1$) from attractive to strongly
repulsive interactions [Fig. 1(b), (c)]. The two particles tend to avoid the
narrow region $I\\!I$ and delocalize away from their interaction centers in
regions $I$ and $I\\!I\\!I$. For strong repulsion ($g=20$), the probability to
find the two particles inside the interval $(-c,c)$ is very small for both the
bosonic and the fermionic state. Note that as the interactions become even
more repulsive ($g=100$), the corresponding probability tends to zero.

At the attractive side ($g=-1$), the two particles tend to bunch, both in the
bosonic and the fermionic state. In fact, when $g\to-\infty$ the energies of
these two levels tend to $-\infty$, corresponding to two deeply bound two-body
bound states Frost (1954, 1956); Scott  _et al._ (1993). This feature that the
lowest two states when $g\to-\infty$ are the symmetric and antisymmetric
combinations of identical bound states is independent of $c$.

Considering now the energy levels for $c=1.5>1$ [Fig. 2(a)], their asymptotic
behavior at $\left|g\right|\gg 1$ hints at different degeneracies compared to
the case where $c<1$. The absence of any triple degeneracy in the displayed
energy range is attributed to the fact that no roots of low-lying Hermite
polynomials exist close to $c=1.5$ and therefore dark states are not present
for the low energy levels. High-level states with zeros near $c=1.5$ certainly
exist and would be dark (or nearly dark) to this interaction, but we do not
depict them.

Aside from the two-level clustering exhibited by a few eigenstates (e.g.
$n=1,2$), some low-energy levels are singly-degenerate (e.g. $n=0,3$ for
$g>0$). When a singly-degenerate eigenstate is occupied, the two particles
exhibit a bunching behavior inside the interval $(-c,c)$, becoming even more
pronounced at strong interaction. A characteristic example is the ground
bosonic state [Fig. 2(b)]. For very strong interaction ($g=100$), the two
particles are almost entirely localized within the interval $(-c,c)$. However,
when the particles reside in bosonic or fermionic states becoming doubly
degenerate at large $g$, they display the opposite behavior. Namely, they are
expelled further away from region $I\\!I$ as $g$ further increases [first
fermionic state in Fig. 2(c)]. This pattern is consistent with the one already
encountered in the doubly-degenerate eigenstates at $c<1$ [Fig. 1(b), (c)].

To understand the bunching manifested in singly-degenerate eigenstates, we
first focus on the transition between single and double degeneracy. In that
regard, the energy spectrum is investigated for variable displacement $c$,
while keeping the interaction strength fixed at a large value, $g=10$ [Fig.
3]. When $c\ll 1$, bosonic and fermionic states with adjacent quantum numbers
are doubly degenerate (_exclusion regime_). Their energies tend to the levels
of non-interacting fermions at $c=0$ [blue dash-dotted lines on the left of
Fig. 3]. For slightly larger $c$, the doubly degenerate energy levels increase
since the particles are pushed away from the interaction centers, to the edges
of the harmonic trap where the potential energy is higher [Fig. 1 (b), (c)].

As $c$ further increases, there are particular points where the energy levels
of bosonic states with quantum numbers $n$ and $n+2$ approach the energy of a
fermionic eigenstate with number $n+1$ [e.g. dashed circle in Fig. 3]. These
points mark a triple degeneracy, which will be clearly manifested as
$g\to\infty$ [see Fig. 1 (a)]. The corresponding $c$ values lie very close to
the first root of even Hermite polynomials, thus further establishing the link
between triple degeneracy and dark states. When the displacement is tuned a
bit further away from these points, the energy levels of the adjacent bosonic
states anticross, and one of them becomes singly-degenerate [e.g. $n=0$ at
$c\gtrsim 1$ or $n=2$ at $c\lesssim 1$]. This implies that dark states mark
the onset of singly-degenerate eigenstates. Moreover, they set the boundaries
between single and double degeneracy.

Apart from bosonic eigenstates, the energy levels of fermionic states can also
be singly-degenerate. The transition from double to single degeneracy takes
place close to the roots of odd Hermite polynomials [e.g. dash-dotted circle
in Fig. 3]. At these values, two adjacent fermionic levels with quantum
numbers $n$ and $n+2$ become degenerate with a bosonic state corresponding to
$n+1$.

Close to the triple degeneracy points, doubly- and singly-degenerate energy
levels coexist, in what we call the _crossover regime_. As the displacement
$c$ however is further tuned to higher values, singly-degenerate eigenstates
outnumber any other kind at low energies. These levels tend towards the non-
interacting bosonic and fermionic values as $c\gg 1$ [blue dash-dotted lines
on the right of Fig. 3]. In this region, the interaction centers are located
at the edges of the harmonic trap, and the two particles barely feel any
interaction. The oscillator length is the only relevant length scale,
dictating the exponential decay of the eigenstates at large separations. Since
the energy spectrum resembles that of an harmonic oscillator, truncated at the
edges due to the interaction centers, the regime of large interaction
displacement is called the _truncation region_.

Figure 3: Energy spectrum for varying displacement $c$ at $g=10$. The energy
levels correspond to even ($n=0,2,\ldots$) and odd ($n=1,3,\ldots$) parity
solutions. The dash and dash-dotted circles mark the onset of triple
degeneracy and the existence of bosonic and fermionic dark states
respectively. The horizontal blue dash-dotted lines on the left correspond to
the non-interacting fermionic energy levels. The respective lines on the
right-hand side correspond to the non-interacting energy spectrum.

Despite the smooth single-to-double degeneracy transition of the energy levels
[Fig. 3], the two particles display a completely different structure when
residing in a singly- or doubly-degenerate eigenstate [Figs. 1(b), 2(b)]. What
is the connection between these two different patterns? To answer that
question we probe the widths of the eigenstates, $\braket{x_{n}^{2}}=\int
dx\,x^{2}\left|\phi_{n}(x)\right|^{2}$, with respect to the displacement $c$.
Focusing on the first bosonic and fermionic eigenstates [$n=0,1$ in Fig. 4],
we observe that their extent is almost identical for $c<1$. This is a
manifestation of the double degeneracy, where the exclusion of the particles
inside the interval $(-c,c)$ results in identical density patterns for
adjacent bosonic and fermionic states [Fig. 1(b), (c)].

Close to the triple degeneracy mark however [dashed circle in Fig. 4], small
variations in the displacement $c$ result in a substantial drop of the spatial
extent of the ground bosonic state, almost an order of magnitude. The
particles are now found at very small separations, a manifestation of the
bunching effect for singly-degenerate eigenstates [see also Fig. 2 (b)]. As
$c$ is tuned to further larger values, the extent of the ground state
asymptotes to the non-interacting width, $\braket{x_{0}^{2}}=1$ Sakurai (1967)
[horizontal blue dash-dotted line in Fig. 4]. Note that $\braket{x_{0}^{2}}$
approaches unity from below, indicating that the bunching is strong, and the
particles localize within smaller distances than the oscillator length. The
same abrupt drop occurs in the spatial extent of the ground fermionic state as
well, in the vicinity of the first root of the $n=3$ Hermite polynomial [dash-
dotted circle in Fig. 4]. From this point on, this fermionic state becomes
singly-degenerate, and subsequently $\braket{x_{1}^{2}}$ saturates to the non-
interacting value, 3 [black dash-dotted horizontal line in Fig. 4].

The spatial extent of higher excited states features a very interesting
pattern. For the first excited bosonic state for example, $\braket{x_{2}^{2}}$
goes through a series of sudden drops and increases [$n=2$ in Fig. 4], as $c$
tends to larger values. These abrupt changes occur nearby triple degeneracy
points [see also Fig. 3]. The two drops [$c\simeq 0.75,2.3$] are attributed to
the single degeneracy and hence to the bunching effect. The one increase in
between [$c\simeq 1$] takes place because from this point on the $n=2$
eigenstate becomes doubly degenerate with the $n=1$ fermionic state. Due to
these series of abrupt transitions, the spatial extent $\braket{x_{2}^{2}}$
assumes its asymptotic value 5 [yellow dash-dotted horizontal line] at a
larger displacement in comparison to the other two eigenstates. Let us also
note that for even stronger repulsions, the changes in $\braket{x_{n}^{2}}$
near the triple degeneracy points grow increasingly abrupt.

##  IV Infinite interactions

Despite the fact that the onset of single, double and triple degeneracy can be
clearly observed already from large interactions [Figs. 1, 2 and 3], the
degeneracies occur only at infinite interaction strengths. Even though this is
a special regime requiring separate treatment, further insight can be gained
regarding the structure of eigenstates with different degree of degeneracy.

###  IV.1 Energy levels

Figure 4: Spatial extents $\langle x_{n}^{2}\rangle$ of the first two bosonic
eigenstates ($n=0,2$) and the fermionic ground state ($n=1$) with respect to
the displacement $c$ at $g=+10$. The dash and dash-dotted circles mark the
onset of triple degeneracy points [Fig. 3]. The horizontal dash-dotted lines
correspond to the widths of the non-interacting eigenstates $\langle
x^{2}_{n}\rangle=(2n+1)$.

To find the energy level structure at $g=+\infty$, we proceed as in Sec. III.
Namely, the relative $x$ coordinate is separated into three regions, $I$,
$I\\!I$ and $I\\!I\\!I$. The boundary conditions that the wave functions and
their derivatives satisfy at the interval intersections are the same as for
finite $g$ [Eqs. (13a)-(13d)]. The most important information that we get from
Eqs. (13c), (13d) is that the wavefunctions have to vanish at the interaction
centers, $x=\pm c$. This guarantees that the first derivatives remain finite.

This leads us to update the setup into three disjoint regions, given that
there is effectively a hard wall at the interaction centers. The relative
Hamiltonian in regions $I$ and $I\\!I\\!I$ corresponds to the Hamiltonian of a
single particle confined in a harmonic oscillator, truncated by a hard-wall
boundary at $x=-c$ ($I$) or $x=c$ ($I\\!I\\!I$). In region $I\\!I$,
$H_{\rm{rel}}$ corresponds to that of a single particle confined in a box
potential superimposed with a harmonic oscillator.

We first focus on the latter region. The relative wave functions can be
written as a superposition of parabolic cylinder functions [see also Eq.
(11)],

$\phi_{n}^{(in)}(x)=\alpha^{(in)}D_{Q^{(in)}_{n}}(x)+\beta^{(in)}D_{Q^{(in)}_{n}}(-x),$
(17)

where $Q^{(in)}_{n}\equiv\epsilon^{(in)}_{n}-1/2$. The superscript denotes
that the two particles are located inside the interval determined by the
interaction centers, $(-c,c)$. The hard-wall boundary conditions imposed by
the box potential imply that $\phi^{(in)}_{n}(-c)=\phi^{(in)}_{n}(c)=0$. These
two equations can be cast in the form of a matrix eigenvalue problem with zero
eigenvalue Consortini and Frieden (1976),

$\begin{pmatrix}D_{Q^{(in)}_{n}}(c)&D_{Q^{(in)}_{n}}(-c)\\\
D_{Q^{(in)}_{n}}(-c)&D_{Q^{(in)}_{n}}(c)\end{pmatrix}\begin{pmatrix}\alpha^{(in)}\\\
\beta^{(in)}\end{pmatrix}=0.$ (18)

Demanding that the determinant of the matrix be zero, we end up with the
following relation for the energy levels $\epsilon^{(in)}_{n}$ Grosche (1993);
Jafarov  _et al._ (2020); Aguilera-Navarro  _et al._ (1980); Aquino and Cruz
(2017),

$D^{2}_{Q^{(in)}_{n}}(c)-D^{2}_{Q^{(in)}_{n}}(-c)=0.$ (19)

For the other two symmetric regions it suffices to solve the relative
Hamiltonian in only one of them. This is due to the fact that $H_{\rm{rel}}$
is invariant under spatial inversion ($\Pi$), even for infinite interactions
[see also Sec. II]. As a result, the energy levels in regions $I$ and
$I\\!I\\!I$ coincide, denoted as $\epsilon^{(out)}_{\nu}$. Focusing on region
$I\\!I\\!I$ for example, the wave functions read
$\chi^{(out,I\\!I\\!I)}_{\nu}(x)=\alpha^{(out,I\\!I\\!I)}D_{Q^{(out)}_{\nu}}(x)$,
where $Q^{(out)}_{\nu}\equiv\epsilon^{(out)}_{\nu}-1/2$. The superscript
denotes that both particles are located outside of the interaction centers
interval $(-c,c)$, and in particular in region $I\\!I\\!I$. Imposing the hard-
wall boundary condition at $x=c$, the energy levels $\epsilon^{(out)}_{\nu}$
are determined Grosche (1993),

$D_{Q^{(out)}_{\nu}}(c)=0.$ (20)

The associated eigenstates $\chi^{(out,I\\!I\\!I)}_{\nu}(x)$ however do not
transform properly under the $\Pi$ and $\Sigma$ operations, since they are
mapped to interval $I$, e.g.
$\hat{\Sigma}\chi^{(out,I\\!I\\!I)}_{\nu}(x)=\chi^{(out,I)}_{\nu}(x)$. To
resolve this issue adequate superpositions of the eigenstates in the two
regions are constructed,

$\displaystyle\phi^{(out)}_{n}(x)=\frac{1}{\sqrt{2}}\left[\chi^{(out,I)}_{\nu}(x)+(-1)^{n}\chi^{(out,I\\!I\\!I)}_{\nu}(x)\right],$
(21) $\displaystyle
n=\begin{cases}2\nu+\theta\left[(n\mod\nu)-1/2\right],&\nu>0\\\
\theta[n-1/2],&\nu=0\end{cases},$ (22)

where $\theta(\cdot)$ is the Heaviside step function. The
$\phi^{(out)}_{n}(x)$ wave functions now span the entire region outside of the
interval $(-c,c)$. Moreover, according to Eq. (22), every quantum number $\nu$
corresponds to a pair of adjacent even and odd quantum numbers $n=2\nu$ and
$n=2\nu+1$. In this way, the $\phi^{(out)}_{n}(x)$ eigenstates can be
classified as even or odd parity under the action of the spatial inversion and
particle permutation operations, similar to finite $g$. The energies of even
and odd parity states with adjacent quantum numbers $2\nu$ and $2\nu+1$ are
degenerate and both equal to $\epsilon^{(out)}_{\nu}$.

Figure 5: (a) Energy levels at $g=+\infty$ versus the displacement $c$. The
positions of the exact crossings between the doubly-degenerate
$\epsilon^{(out)}_{n}$ and singly-degenerate $\epsilon^{(in)}_{n}$ levels
designate three regions, the exclusion (E), crossover (C), and truncation (T)
region.

From the above analysis it becomes clear that the bosonic and fermionic levels
$\epsilon^{(out)}_{n}$ are doubly-degenerate, and that particles occupying
these eigenstates are localized outside of the interval $(-c,c)$. On the other
hand, the two particles are strictly found within $(-c,c)$ when the
$\phi^{(in)}_{n}$ eigenstates are populated. In this case the bosonic and
fermionic levels are distinct [Eq. (19)], they are singly-degenerate. Varying
the displacement $c$ is equivalent to moving the hard walls, and thus shifting
the energy levels $\epsilon^{(in)}_{n}$ and $\epsilon^{(out)}_{n}$ [Fig. 5].
What distinguishes the spectra from those with large finite $g$ [Fig. 3] is
that the hard walls at $x=\pm c$ prevent any transition between singly- and
doubly-degenerate eigenstates. This is captured in Fig. 5, where there are
only exact crossings between $\epsilon^{(in)}_{n}$ and $\epsilon^{(out)}_{n}$
for variable displacement $c$. The positions of these crossings correspond to
roots of Hermite polynomials, marking the presence of dark states, and hence
triple degeneracy.

At large finite interactions, the exact crossings become avoided [Fig. 3]. In
this case, the double-to-single degeneracy transition manifests, identified in
Figs. 3 and 4. In this regime, there are no hard walls at $x=\pm c$ [or
additional symmetry, see the follow-up Section] preventing the particles from
changing their spatial configuration. When $c$ lies close to the roots of
Hermite polynomials, a transition occurs between a spatial structure
resembling the one for $\phi^{(in)}_{n}$ (single-degeneracy) and the one
similar to $\phi^{(out)}_{n}$ (double-degeneracy). From this vantage point,
the bunching occurs because the density profiles of singly-degenerate
eigenstates at $g\gg 1$ resemble the structure of $\phi^{(in)}_{n}$ at
$g=+\infty$.

Apart from understanding the bunching effect, the energy levels at $g=+\infty$
[Fig. 5] help us also to neatly classify the density profiles of the two-
particle system, depending on the value of $c$. In particular,
$\epsilon^{(in)}_{0}$ and $\epsilon^{(out)}_{0,1}$ delineate the three regimes
that were identified before [see also Sec. III]. On the left of
$\epsilon^{(in)}_{0}$, there is the _exclusion regime_ (E). All eigenstates
are doubly-degenerate, and particles are excluded from the interval $(-c,c)$.
On the right of $\epsilon^{(out)}_{0,1}$, all eigenstates are singly-
degenerate, marking the _truncation regime_ (T). Particles are localized
within the interaction center interval $(-c,c)$, and the energy levels
saturate to their non-interacting values at large $c$. In-between these two
regions, doubly- and singly-degenerate eigenstates coexist, denoting the
_crossover regime_ (C). In this regime dark states occur at the positions of
the exact crossings, a sign of triple degeneracy. Of course, such a
classification scheme serves only as a guide when tackling large finite
interactions. The boundaries between the regions slightly change due to the
occurrence of avoided crossings [Fig. 3].

###  IV.2 Symmetry analysis

For finite $g$ and $c$, all energy levels of the relative sub-Hamiltonian are
singly-degenerate, as is true for any solutions to the 1D Schrödinger equation
defined on a path-connected interval. More generally, any system whose maximal
kinematic symmetry is Abelian should only be singly-degenerate. How then can
we understand the additional double and triple degeneracies that occur in the
$g\to\infty$ limit?

These can be understood by recognizing that the domain of the relative
coordinate $x$ becomes effectively disconnected for infinite $g$. Each of the
three regions $I$, $I\\!I$ and $I\\!I\\!I$ acts like an independent quantum
system, each experiencing independent time evolution. Energy eigenstates are
localized to each of these three independent regions in the $g\to\infty$
limit. To see this, note that the relative sub-Hamiltonian
$H_{\mathrm{rel}}(g,c)$ can be re-expressed as the direct sum of three
Hamiltonians:

$\lim_{g\to\infty}H_{\mathrm{rel}}(g,c)=H_{I}\oplus H_{I\\!I}\oplus
H_{I\\!I\\!I}$ (23)

where $H_{R}$ is the restriction of the 1D harmonic oscillator Hamiltonian
restricted to region $R\in\\{I,I\\!I,I\\!I\\!I\\}$. These restricted
Hamiltonians commute, so instead of one time translation symmetry
parameterized by $t\in T_{t}\sim\mathbb{R}$, there are now three parameterized
by $(t_{I},t_{I\\!I},t_{I\\!I\\!I})\in T_{t}^{\times 3}\sim\mathbb{R}^{3}$ and
represented by products of the three unitary operators
$\hat{U}_{R}(t_{R})=\exp(-iH_{R}t_{R})$. Equivalently, the phase difference
between disjoint regions in the $g\to\infty$ limit is not an observable
quantity Harshman (2017a). Only when $g$ is finite, is there coupling between
adjacent regions that locks their relative phase.

However, this additional symmetry of tripled time evolution is still an
Abelian kinematic symmetry, and therefore not enough to explain the systematic
double and triple degeneracies that occur for the $g\to\infty$ limit of
$H_{\mathrm{rel}}(g,c)$. If regions $I$, $I\\!I$, and $I\\!I\\!I$ were all
intervals with inequivalent domains, then the spectrum would be singly-
degenerate except for so-called accidental degeneracies where two states (or
even more rarely, three states) would coincide in energy. For our system, the
regions $I$ and $I\\!I\\!I$ are equivalent intervals exchanged by relative
parity $\Sigma$. Therefore the spectrum of $H_{I}$ and $H_{I\\!I\\!I}$
coincide and the out states $\chi_{\nu}^{(out,I)}$ and
$\chi_{\nu}^{(out,I\\!I\\!I)}$ are double-degenerate pairs that can be
symmetrized into bosonic or antisymmetrized into fermionic states
$\phi_{n}^{(out)}$.

The kinematic symmetry group that incorporates the equivalence of the regions
$I$ and $I\\!I\\!I$ is formed from the time evolution operators
$\hat{U}_{I}(t_{I})$ and $\hat{U}_{I\\!I\\!I}(t_{I\\!I\\!I})$ (which commute)
and the relative parity operator $\hat{\Sigma}$, which satisfies

$\hat{\Sigma}\hat{U}_{I}(t_{I})=\hat{U}_{I\\!I\\!I}(t_{I\\!I\\!I})\hat{\Sigma}.$
(24)

Because $\hat{\Sigma}$ does not commute with the time translations, this
kinematic symmetry group is not the direct product $T_{t}^{\times 2}\times
S_{2}$ but must be expressed as the semidirect product $T_{t}^{\times
2}\rtimes S_{2}$, where $(t_{I},t_{I\\!I\\!I})\in T_{t}^{\times
2}\sim\mathbb{R}^{2}$ is the time translation group generated by $H_{I}$ and
$H_{I\\!I\\!I}$, $S_{2}$ is the group generated by $\Sigma$ that permutes the
two identical systems, and $\rtimes$ indicates that $S_{2}$ acts a non-trivial
automorphism on $T_{t}^{\times 2}$. This specific form of the semidirect
product is also known as the wreath product $T_{t}\wr S_{2}$ Bhattacharjee
_et al._ (1998); Harshman (2017a, b). The kinematic symmetry group $T_{t}\wr
S_{2}$ is non-Abelian and one can show that it has two-dimensional unitary
irreducible representations that explain the double-degeneracy of the out-
states for $g\to\infty$ and any $c$.

Further, triple degeneracies occur when states in the spectrum of region
$I\\!I$ have energies that coincide with states of the outer regions $I$ and
$I\\!I\\!I$. These spectral ‘accidents’ occur precisely when the displacements
$\pm c$ align with the $n$th Hermite polynomial root, $H_{n}(c/\sqrt{2})=0$.
When $n$ is even, the state in region $I\\!I$ is even, and the triply-
degenerate state is comprised of two bosonic states (one inside and two
outside the central region) and one fermionic state. When $n$ is odd, triple
degeneracy occurs between two fermionic states and one bosonic. Understanding
these accidental degeneracies as the result of dark states provides an
interesting perspective on an old subject Jauch and Hill (1940); McIntosh
(1959); Louck and Metropolis (1981); Moshinsky and Quesne (1983); Leyvraz  _et
al._ (1997), but not all accidental degeneracies can be understood this way.

##  V Summary and Conclusions

We have investigated the stationary properties of two harmonically trapped
particles interacting via contact potentials with a displacement $c$.
Depending on the value of the latter, the energy spectra are classified into
three regimes. The exclusion regime for $c$ smaller than the oscillator
length, where the energy levels of adjacent bosonic and fermionic states
cluster together as $g\to\infty$. The two particles are found at large
separations, and they are expelled from the interval dictated by the
interaction centers, $(-c,c)$. In the truncation regime occurring at $c$
larger than the oscillator length, all energy levels are singly-degenerate as
$c\to\infty$. The corresponding relative wavefunctions have a non-zero support
only at short interparticle distances, signalling a bunching effect. This is
understood in terms of the stationary properties at infinite interactions,
where the relative Hamiltonian in the interval $(-c,c)$ is equivalent to that
of a box potential superimposed with a harmonic oscillator. In the crossover
regime, as the name suggests, both singly- and doubly-degenerate eigenstates
coexist.

The boundary between these two kinds of eigenstates is set by dark states. The
latter occur whenever the interaction displacement $c$ lies close to a root of
an Hermite polynomial. Whenever such eigenstates are occupied, the two
particles do not experience any interaction. At finite interaction strengths,
there is a transition between singly- and doubly-degenerate eigenstates,
manifested as avoided crossings in the energy spectra when $c$ coincides with
Hermite polynomial roots. At $g=+\infty$, such a transition is prohibited
since the relative Hamiltonian partitions into three disjoint regions. Due to
the extra symmetries that this decomposition introduces, the avoided crossings
become exact. The doubly-degenerate adjacent bosonic and fermionic eigenstates
cluster with a dark state, i.e. triple degeneracy occurs.

Apart from the two-atom problem considered here, the few-body and many-body
aspects of such off-centered interactions are certainly intriguing. The
displacement $c$ introduces an additional length scale that competes with the
scales present in a many-body setup. Such a competition may lead to novel
phases, as in dipolar gases Chomaz  _et al._ (2022); Lahaye  _et al._ (2009)
for instance. Moreover, it is interesting to investigate which degeneracies
occur when considering external trapping potentials apart from the harmonic
oscillator. Dark states depend critically on trap shape and the number of
particles, leading to rearrangement of the energy levels.

Ultracold atoms trapped within optical tweezers Kaufman and Ni (2021);
Andersen (2022) may offer an experimental scheme for realizing such a model
Hamiltonian. In particular, the relative Hamiltonian is equivalent to the
Born-Oppenheimer description of a single trapped impurity interacting with two
non-interacting atoms, fixed at positions $\pm c$ from the impurity.
Capitalizing on optical tweezers, the displacement can be adjusted, realizing
the different regimes described above. In that regard, the density profile of
the impurity could be engineered, and induced interactions can occur between
the two fixed atoms.

The fact that the interactions among particles occur when they are apart at a
distance c could be considered as a decentered interaction. Indeed, such
decentered interactions are encountered for atoms in the presence of
perpendicular strong electric and magnetic fields Schmelcher and Cederbaum
(1993); Dippel  _et al._ (1994). The action at a given distance leads to so-
called giant dipole states in crossed fields in case of the attractive Coulomb
potential.

###### Acknowledgements.

This work is supported by the Cluster of Excellence ‘Advanced Imaging of
Matter’ of the Deutsche Forschungsgemeinschaft (DFG) - EXC2056 - project ID
390715994. NH additionally acknowledges the support of the Deutscher
Akademischer Austauschdienst and the U.S. Fulbright Specialist Program.

## References

  * Pethick and Smith (2008) C. J. Pethick and H. Smith,  _Bose–Einstein Condensation in Dilute Gases_ , 2nd ed. (Cambridge University Press, 2008).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Mistakidis  _et al._ (2023) S. Mistakidis, A. Volosniev, R. Barfknecht, T. Fogarty, T. Busch, A. Foerster, P. Schmelcher, and N. Zinner, Phys. Rep. 1042, 1 (2023).
  * Wenz  _et al._ (2013) A. N. Wenz, G. Zürn, S. Murmann, I. Brouzos, T. Lompe, and S. Jochim, Science 342, 457 (2013).
  * Serwane  _et al._ (2011) F. Serwane, G. Zürn, T. Lompe, T. Ottenstein, A. Wenz, and S. Jochim, Science 332, 336 (2011).
  * Sowiński and Ángel García-March (2019) T. Sowiński and M. Ángel García-March, Rep. Prog. Phys. 82, 104401 (2019).
  * Busch  _et al._ (1998) T. Busch, B.-G. Englert, K. Rzażewski, and M. Wilkens, Found. Phys. 28, 549 (1998).
  * Farrell and van Zyl (2009) A. Farrell and B. P. van Zyl, J. Phys. A: Math. Theor. 43, 015302 (2009).
  * Lieb and Liniger (1963) E. H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963).
  * Gaudin (2014) M. Gaudin,  _The Bethe Wavefunction_ (Cambridge University Press, 2014).
  * García-March  _et al._ (2014) M. A. García-March, B. Juliá-Díaz, G. E. Astrakharchik, T. Busch, J. Boronat, and A. Polls, New J. Phys. 16, 103004 (2014).
  * Harshman (2016a) N. L. Harshman, Few-Body Syst. 57, 11 (2016a).
  * Harshman (2016b) N. L. Harshman, Few-Body Syst. 57, 45 (2016b).
  * Sutherland (2004) B. Sutherland,  _Beautiful Models: 70 Years of Exactly Solved Quantum Many-Body Problems_ (World Scientific Publishing Company, River Edge, N.J, 2004).
  * Girardeau (1960) M. Girardeau, J. Math. Phys. 1, 516 (1960).
  * Albeverio  _et al._ (2012) S. Albeverio, F. Gesztesy, R. Hoegh-Krohn, and H. Holden,  _Solvable Models in Quantum Mechanics_ , softcover reprint of the original 1st ed. 1988 edition ed. (Springer, Place of publication not identified, 2012).
  * Cheon and Shigehara (1998) T. Cheon and T. Shigehara, Phys. Lett. A 243, 111 (1998).
  * Cheon and Shigehara (1999) T. Cheon and T. Shigehara, Phys. Rev. Lett. 82, 2536 (1999).
  * Loudon (1959) R. Loudon, Am. J. Phys. 27, 649 (1959).
  * Harshman (2017a) N. L. Harshman, Phys. Rev. A 95, 053616 (2017a).
  * Louck (1965) J. D. Louck, J. Math. Phys. 6, 1786 (1965).
  * Niederer (1973) U. Niederer, Helv. Phys. Acta 46, 191 (1973).
  * Abramowitz and Stegun (1948) M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Vol. 55 (US Government printing office, 1948).
  * Aouadi  _et al._ (2016) O. Aouadi, Y. Chargui, and M. S. Fayache, J. Math. Phys. 57, 023522 (2016).
  * Avakian  _et al._ (1987) M. Avakian, G. Pogosyan, A. Sissakian, and V. Ter-Antonyan, Phys. Lett. A 124, 233 (1987).
  * Belloni and Robinett (2014) M. Belloni and R. W. Robinett, Phys. Rep. 540, 25 (2014).
  * Budewig  _et al._ (2019) L. Budewig, S. Mistakidis, and P. Schmelcher, Mol. Phys. 117, 2043 (2019).
  * Salzer  _et al._ (1952) H. Salzer, R. Zucker, and R. Capuano, J. Res. Natl. Bur. Stand. (U. S.) 48, 111 (1952).
  * Tonks (1936) L. Tonks, Phys. Rev. 50, 955 (1936).
  * Kehrberger  _et al._ (2018) L. M. A. Kehrberger, V. J. Bolsinger, and P. Schmelcher, Phys. Rev. A 97, 013606 (2018).
  * Frost (1954) A. A. Frost, J. Chem. Phys. 22, 1613 (1954).
  * Frost (1956) A. A. Frost, J. Chem. Phys. 25, 1150 (1956).
  * Scott  _et al._ (1993) T. C. Scott, J. F. Babb, A. Dalgarno, and J. D. Morgan, III, J. Chem. Phys. 99, 2841 (1993).
  * Sakurai (1967) J. J. Sakurai,  _Advanced quantum mechanics_ (Pearson Education India, 1967).
  * Consortini and Frieden (1976) A. Consortini and B. R. Frieden, Nuov. Cim. B 35, 153 (1976).
  * Grosche (1993) C. Grosche, Annalen der Physik 505, 557 (1993).
  * Jafarov  _et al._ (2020) E. Jafarov, S. Nagiyev, and A. Jafarova, Rep. Math. Phys. 86, 25 (2020).
  * Aguilera-Navarro  _et al._ (1980) V. C. Aguilera-Navarro, E. L. Koo, and A. H. Zimerman, J. Phys. A: Math. and Gen. 13, 3585 (1980).
  * Aquino and Cruz (2017) N. Aquino and E. Cruz, Revista mexicana de física 63, 580 (2017).
  * Bhattacharjee  _et al._ (1998) M. Bhattacharjee, D. Macpherson, R. G. Möller, and P. M. Neumann, in  _Notes on Infinite Permutation Groups_, Lecture Notes in Mathematics (Springer, Berlin, Heidelberg, 1998) pp. 67–76.
  * Harshman (2017b) N. L. Harshman, Few-Body Syst. 58, 41 (2017b).
  * Jauch and Hill (1940) J. M. Jauch and E. L. Hill, Phys. Rev. 57, 641 (1940).
  * McIntosh (1959) H. V. McIntosh, Am. J. Phys. 27, 620 (1959).
  * Louck and Metropolis (1981) J. D. Louck and N. Metropolis, Adv. Appl. Math. 2, 138 (1981).
  * Moshinsky and Quesne (1983) M. Moshinsky and C. Quesne, Ann. Phys. 148, 462 (1983).
  * Leyvraz  _et al._ (1997) F. Leyvraz, A. Frank, R. Lemus, and M. V. Andrés, Am. J. Phys. 65, 1087 (1997).
  * Chomaz  _et al._ (2022) L. Chomaz, I. Ferrier-Barbut, F. Ferlaino, B. Laburthe-Tolra, B. L. Lev, and T. Pfau, Rep. Prog. Phys. (2022).
  * Lahaye  _et al._ (2009) T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, and T. Pfau, Rep. Prog. Phys. 72, 126401 (2009).
  * Kaufman and Ni (2021) A. M. Kaufman and K.-K. Ni, Nat. Phys. 17, 1324 (2021).
  * Andersen (2022) M. F. Andersen, Adv. Phys.: X 7, 2064231 (2022).
  * Schmelcher and Cederbaum (1993) P. Schmelcher and L. S. Cederbaum, Chem. Phys. Lett. 208, 548 (1993).
  * Dippel  _et al._ (1994) O. Dippel, P. Schmelcher, and L. Cederbaum, Phys. Rev. A 49, 4415 (1994).