Semiclassics in a system without classical limit: The few-body spectrum of two interacting bosons in one dimension

TL;DR

This paper develops a semiclassical approach to analyze the energy spectrum of two interacting bosons in one dimension, providing explicit analytical results that match numerical and exact solutions, including bound states.

Contribution

It extends semiclassical methods to a two-body problem with boundary conditions, offering analytical insights into the spectrum of interacting bosons in one dimension.

Findings

Analytical expressions for the smooth density of states

Excellent agreement with numerical calculations

Accurate reproduction of exact energy levels

Abstract

We present a semiclassical study of the spectrum of a few-body system consisting of two short-range interacting bosonic particles in one dimension, a particular case of a general class of integrable many-body systems where the energy spectrum is given by the solution of algebraic transcendental equations. By an exact mapping between $\delta$-potentials and boundary conditions on the few-body wave functions, we are able to extend previous semiclassical results for single-particle systems with mixed boundary conditions to the two-body problem. The semiclassical approach allows us to derive explicit analytical results for the smooth part of the two-body density of states that are in excellent agreement with numerical calculations. It further enables us to include the effect of bound states in the attractive case. Remarkably, for the particular case of two particles in one dimension, the discrete energy levels obtained through a requantization condition of the smooth density of states are essentially in perfect agreement with the exact ones.

Equations

$$\left(-\frac{{\hslash }^2}{2m}\frac{{\partial }^2}{\partial x^2}+\frac{{\hslash }^2\kappa }{m}\delta (x)\right)\psi (x)=E\psi (x)$$

$${\psi }_n^{(\mathrm{a})}(x)=\frac{1}{\sqrt{L}}\sin (k_n^{(\mathrm{a})}x),k_n^{(\mathrm{a})}=\frac{n\pi }{L},n\in \mathbb{N}.$$

$$\begin{split}\psi_{n}^{(\mathrm{s})}(x)&=\frac{A_{n}}{\sqrt{L}}\sin\bm{(}k_{n}^{(\mathrm{s})}(|x|-L)\bm{)},\\ k_{n}^{(\mathrm{s})}L&=n\pi-\arctan\!\bigg{(}\frac{k_{n}^{(\mathrm{s})}}{\kappa}\bigg{)},\end{split}$$

$${\psi }_0^{(\mathrm{s})}(x)=\sqrt{\frac{3}{2L^3}}(|x|-L)$$

$$\begin{array}{rl}{K}_{\lfloor \delta \rfloor }(x^{\prime },x,t)=& \sum _{n=1}^{\infty }{\mathrm{e}}^{-\frac{\mathrm{i}\hslash t}{2m}(k_n^{(\mathrm{a})}{)}^2}{\psi }_n^{(\mathrm{a})}(x^{\prime }){\psi }_n^{(\mathrm{a})}(x)\\ & +\sum _{n=0}^{\infty }{\mathrm{e}}^{-\frac{\mathrm{i}\hslash t}{2m}(k_n^{(\mathrm{s})}{)}^2}{\psi }_n^{(\mathrm{s})}(x^{\prime }){\psi }_n^{(\mathrm{s})}(x).\end{array}$$

$$\begin{array}{rl}{K}_{\delta }(x^{\prime },x,t)=& \frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}k{\mathrm{e}}^{-\frac{\mathrm{i}\hslash t}{2m}{k}^2}\cos (k(x^{\prime }-x))\\ & -\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}k{\mathrm{e}}^{-\frac{\mathrm{i}\hslash t}{2m}{k}^2}\mathrm{Re}\frac{{\mathrm{e}}^{ik(|x^{\prime }|+|x|)}}{1+i\frac{k}{\kappa }}\\ & -\Theta (-\kappa )\kappa {\mathrm{e}}^{\frac{\mathrm{i}\hslash t}{2m}{\kappa }^2}{\mathrm{e}}^{\kappa (|x^{\prime }|+|x|)},\end{array}$$

$$K_0(x^{\prime },x,t)=\sqrt{\frac{m}{2\pi \mathrm{i}\hslash t}}{\mathrm{e}}^{-\frac{m}{2\mathrm{i}\hslash t}(x^{\prime }-x{)}^2}.$$

$$\frac{1}{1+\mathrm{i}\frac{k}{\kappa }}={\int }_0^{\infty }\mathrm{d}ϵ{\mathrm{e}}^{-(1+\mathrm{i}\frac{k}{\kappa })ϵ}.$$

$$K_{\delta }(x^{\prime },x,t)=K_0(x^{\prime },x,t)+K_{\kappa }(x^{\prime },x,t),$$

$$\begin{array}{rl}{K}_{\kappa }& (x^{\prime },x,t)=\frac{1\mp 1}{2}\kappa {\mathrm{e}}^{\frac{\mathrm{i}\hslash t}{2m}{\kappa }^2}{\mathrm{e}}^{-\kappa (|x^{\prime }|+|x|)}\\ & -\kappa \sqrt{\frac{m}{2\pi \mathrm{i}\hslash t}}{\int }_0^{\infty }\mathrm{d}ϵ{\mathrm{e}}^{-\kappa ϵ}{\mathrm{e}}^{-\frac{m}{2i\hslash t}(|x^{\prime }|+|x|\pm ϵ{)}^2}.\end{array}$$

$$-{\int }_0^{\infty }\mathrm{d}ϵ\kappa {\mathrm{e}}^{-\kappa ϵ}{K}_0(-|x^{\prime }|,|x|\pm ϵ,t),$$

$$\rho(E)=\mathcal{L}_{\beta}^{-1}\Big{[}\int\mathop{}\!\mathrm{d}^{d}{x}\,K(\mathbf{x},\mathbf{x},t=-\mathrm{i}\hbar\beta)\Big{]}(E).$$

$$K^{(d)}({\mathbf{x}}^{\prime },\mathbf{x},t)=$$

$$=$$

$$+K_0^{(d-1)}({\mathbf{x}}_{\parallel }^{\prime },{\mathbf{x}}_{\parallel },t)K_{\kappa }(x_{\perp }^{\prime },x_{\perp },t).$$

$$\mu =\frac{{\hslash }^2{\kappa }^2}{2m},$$

$$\begin{split}\int_{\Omega}\mathop{}\!\mathrm{d}^{d}{x}\,&K(\mathbf{x},\mathbf{x},t=-\mathrm{i}\hbar\beta)=V_{\Omega}\Big{(}\frac{m}{2\pi\hbar^{2}\beta}\Big{)}^{\frac{d}{2}}\\ &+\frac{S_{\delta}}{2}\Big{(}\frac{m}{2\pi\hbar^{2}\beta}\Big{)}^{\frac{d-1}{2}}\\ &\quad\times\big{(}-1\pm\mathrm{e}^{\mu\beta}\operatorname{erfc}(\sqrt{\mu\beta})+(1\mp 1)\mathrm{e}^{\mu\beta}\big{)}\end{split}$$

$$\begin{split}\rho(E)=&V_{\Omega}\Big{(}\frac{m}{2\pi\hbar^{2}}\Big{)}^{\frac{1}{2}}\frac{\Theta(E)}{\sqrt{\pi E}}-\frac{1}{2}\delta(E)\\ &\pm\frac{1}{2\pi}\sqrt{\frac{\mu}{E}}\frac{\Theta(E)}{E+\mu}+\frac{1\mp 1}{2}\delta(E+\mu).\end{split}$$

$$\begin{split}\rho(E)=&V_{\Omega}\left(\frac{m}{2\pi\hbar^{2}}\right)^{\frac{d}{2}}\frac{E^{\frac{d-2}{2}}}{\Gamma(\frac{d}{2})}\Theta(E)-\frac{S_{\delta}}{2}\left(\frac{m}{2\pi\hbar^{2}}\right)^{\frac{d-1}{2}}\frac{E^{\frac{d-3}{2}}}{\Gamma(\frac{d-1}{2})}\Theta(E)\\ &\pm\frac{S_{\delta}}{2}\left(\frac{m}{2\pi\hbar^{2}}\right)^{\frac{d-1}{2}}\left\{\mathcal{L}_{\beta}^{-1}\Big{[}\frac{1}{\beta^{\frac{d-1}{2}}}\mathrm{e}^{\mu\beta}\operatorname{erfc}(\sqrt{\mu\beta})\Big{]}(E)+(1\mp 1)\frac{(E+\mu)^{\frac{d-3}{2}}}{\Gamma(\frac{d-1}{2})}\Theta(E+\mu)\right\}.\end{split}$$

$$\mathcal{L}_{\beta}^{-1}\Big{[}\frac{1}{\beta^{\frac{d-1}{2}}}\mathrm{e}^{\mu\beta}\operatorname{erfc}(\sqrt{\mu\beta})\Big{]}(E)$$

$$\frac{\partial}{\partial x_{\perp}}\psi(\mathbf{x})\bigg{|}_{\mathbf{x}_{s}}=\kappa\psi(\mathbf{x}_{s}),\quad\mathbf{x}_{s}\in\partial\Omega.$$

$$\begin{array}{rl}{K}_{\delta }^{+}(x^{\prime },x,t)=& K_{\delta }(x^{\prime },x,t)+K_{\delta }(-x^{\prime },x,t)\\ =& K_0(x^{\prime },x,t)+K_0(-x^{\prime },x,t)\\ & +2K_{\kappa }(x^{\prime },x,t)\end{array}$$

$$\frac{S_{\delta }}{4}{\left(\frac{m}{2\pi {\hslash }^2}\right)}^{\frac{d-1}{2}}\frac{E^{\frac{d-3}{2}}}{\Gamma (\frac{d-1}{2})}\Theta (E)$$

$$\beta =\frac{it}{\hslash }\text{and}\frac{{\hslash }^2}{2m}=1.$$

$$\hat{\mathcal{H}}(q_1,q_2)=-\frac{{\partial }^2}{\partial q_1^2}-\frac{{\partial }^2}{\partial q_2^2}+\sqrt{8}\tilde{\kappa }\delta (q_1-q_2).$$

$$x_1=\frac{1}{\sqrt{2}}(q_1-q_2),x_2=\frac{1}{\sqrt{2}}(q_1+q_2)$$

$$\hat{\mathcal{H}}(x_1,x_2)=-\frac{{\partial }^2}{\partial x_1^2}-\frac{{\partial }^2}{\partial x_2^2}+2\tilde{\kappa }\delta (x_1).$$

$$\frac{\partial}{\partial x_{1}}\psi(x_{1},x_{2})\bigg{|}_{x_{1}=0}=\tilde{\kappa}\psi(0,x_{2})$$

$$\begin{array}{rl}\rho (E)=& \frac{L^2}{8\pi }\Theta (E)-\frac{(2+\sqrt{2})L}{8\pi }\frac{\Theta (E)}{\sqrt{E}}\\ & \pm \frac{\sqrt{2}L}{4\pi }\frac{\Theta (E)}{\sqrt{E+\tilde{\mu }}}+{\theta }_1^0\frac{\sqrt{2}L}{2\pi }\frac{\Theta (E+\tilde{\mu })}{\sqrt{E+\tilde{\mu }}}\\ & +{\rho }_c(E),\end{array}$$

$$K_{\frac{\pi }{2}}({\mathbf{x}}^{\prime },\mathbf{x},t)=K_{\delta }^{+}(x_1^{\prime },x_1,t)K_{\delta }^{+}(x_2^{\prime },x_2,t).$$

$$\begin{split}K_{\frac{\pi}{4}}(\mathbf{x^{\prime}},\mathbf{x},t)=&\frac{1}{2}\Big{[}K_{\frac{\pi}{2}}\bm{(}(x_{1}^{\prime},x_{2}^{\prime}),(x_{1},x_{2}),t\bm{)}\\ &-K_{\frac{\pi}{2}}\bm{(}(x_{2}^{\prime},x_{1}^{\prime}),(x_{1},x_{2}),t\bm{)}\Big{]}.\end{split}$$

Full text

I Introduction

The discovery by Bethe of quantum many-body systems admitting analytical expressions for eigenstates and eigenenergies in terms of the implicit solutions of algebraic transcendental equations Bethe (1931); Lieb and Liniger (1963); Yang and Yang (1969) marked the birth of a whole branch of mathematical physics, namely, the theory of quantum integrable models (for an introduction see Korepin et al. (1997)). Since then, these models have served as the playground to study the role of symmetries and conservation laws for the stationary and dynamical properties of many-body systems Gaudin (1971); Batchelor et al. (2005); Sykes et al. (2007); Bolte and Garforth (2017).

Solvable models had been mainly restricted to the area of mathematical physics because they involve apparently unphysical $\delta$-type interparticle interactions and are commonly restricted to one-dimensional (1D) motion. This situation has drastically changed with the successful preparation of quantum states of interacting cold atoms Davis et al. (1995); Bloch (2008), especially in regimes where the system can be considered essentially 1D, as in elongated optical traps Morizot et al. (2006); Grießer and Ritsch (2013); Kim et al. (2016). As it turned out that interactions between neutral cold atoms are described with astonishing precision by $\delta$-type potentials Olshanii (1998); Chin et al. (2010); Cazalilla et al. (2011); Krutitsky (2016), the knowledge accumulated during almost one century of work on solvable models found its way into cold-atom physics during the last decade. Hence, the experimental study of solvable many-body systems has become a very active field that provides hints for other hitherto inaccessible regimes where external potentials destroy integrability.

There are, however, at least two situations where the detailed level-by-level calculation of the many-body spectrum using methods of integrable quantum systems overshoots in the problem of calculating macroscopic properties like the microcanonical partition function. One example is the thermodynamic limit, where quantum fluctuations are suppressed and the spectrum behaves effectively smooth. In this regime the appropriate tool is the so-called thermodynamic Bethe ansatz Yang and Yang (1969). The second situation appears when studying interacting bosonic systems in the mesoscopic short-wavelength regime, where a separation of scales between the smooth and oscillatory part of the level density allows for the approximation of the spectrum as a smooth function. This is a well known procedure often referred to as the Thomas-Fermi approximation Brack and Bhaduri (2003).

The Thomas-Fermi density of states is of paramount importance as it fixes the energy-dependence of the mean-level spacing, the fundamental energy scale that heralds the appearance of quantum effects and determines the relative importance of external perturbations. In this context, cold atom systems pose an interesting problem: since the Thomas-Fermi approximation explicitly requires a classical limit for the quantum mechanical Hamiltonian, it is not well defined for systems with $\delta$-like interactions where the classical limit only exists away from the collisions.

All in all, in systems with few particles interacting through short-range interactions, the semiclassical limit appears to be difficult to handle already when we ask for a very fundamental quantity as the mean level spacing. Our objective in this paper is to show how, by using techniques imported from the calculation of Thomas-Fermi approximations in single-particle systems with mixed boundary conditions, one can obtain a rigorous definition of the smooth density of states in systems of few particles interacting through $\delta$-potentials, even when the latter do not have a classical limit.

II One particle in a $\delta$-potential

In this section we present the general formalism for the calculation of the smooth part of the density of states (DOS) for $d$-dimensional billiards with either $(d-1)$-dimensional $\delta$-barriers [i.e. a $\delta$-function potential along a ($d-1$)-dimensional manifold] inside the volume or with Robin (or mixed) boundary conditions on the surface. The latter case was already discussed by Balian and Bloch Balian and Bloch (1970, 1971) and by Sieber Sieber et al. (1995), but in all three cases the derivations were done via the energy dependent Greens function whereas we use a formalism based on 1D propagators Baltes and Hilf (1976); Hummel et al. (2013). For this we need the propagator for a $\delta$-potential in 1D, which is known exactly from a path integral calculation Bauch (1985) (see also Grosche (1994) for higher dimensions). Notably, quantum mechanical properties, like the appearance of a bound state for attractive interaction, are well hidden inside the results (see e.g. Goovaerts et al. (1973)). Therefore we take an alternative approach which is more instructive for our purposes.

II.1 Propagator for a particle in a $\delta$-potential

The 1D-propagator for the $\delta$-barrier is derived in a straightforward way by first calculating the exact propagator for a particle on a line with Dirichlet boundary conditions at $x=\pm L$ and a $\delta$-potential $V(x)=(\hbar^{2}\kappa/m)\,\delta(x)$ with $\kappa\in\mathbb{R}$ and then taking the limit $L\to\infty$. The solutions of the stationary Schrödinger equation

[TABLE]

for the confined system are well known (see e.g. Pedram and Vahabi (2010)) and can be separated into symmetric and antisymmetric solutions due to the symmetry of the problem. The latter are not affected by the $\delta$-potential and coincide with the antisymmetric solutions of a particle in a box:

[TABLE]

The symmetric solutions are given by

[TABLE]

where $A_{n}$ is a normalization constant that depends on $k_{n}^{(\mathrm{s})}$. The transcendental equation in Eq. (3) has exactly one solution for every $n\in\mathbb{N}$ and another nontrivial solution for $n=0$ if $\kappa$ is negative. In the case $\kappa L\in(-1,0)$ this solution is real whereas the case $\kappa L<-1$ yields a purely imaginary wave number which corresponds to a negative energy, referred to as a bound state in the following. The two types of states are continuously connected by the zero energy solution

[TABLE]

valid for $\kappa L=-1$. In the limit $L\to\infty$ the state $\psi_{0}^{(\mathrm{s})}$ will always be bound irrespective of the value of $\kappa$.

The exact propagator for the confined system can now be written as

[TABLE]

In the case $\kappa\geq 0$ both summations start with $n=1$. In the limit $L\to\infty$ this yields (see Appendix A)

[TABLE]

where the last term originates from the separate treatment of the (bound) state $\psi_{0}^{(\mathrm{s})}$ ($\Theta$ is the Heaviside step function). The first integral in Eq. (6) can be evaluated directly by means of a Gaussian (or Fresnel) integral and yields the well-known expression for the free propagator

[TABLE]

The second integral in Eq. (6) can be evaluated in the same way after replacing

[TABLE]

Finally, the propagator for the Hamiltonian in Eq. (1) reads

[TABLE]

with the deviation from free propagation given by

[TABLE]

Here, we identified $\kappa=|\kappa|$ and the upper (lower) signs stand for a repulsive (attractive) potential. This result generalizes the propagator found in Manoukian (1989); Gaveau and Schulman (1986), which is restricted to the repulsive case, and is equivalent to the result from path integral approaches Bauch (1985). Note that the second term in Eq. (10) can be written as

[TABLE]

which is closely related to the correction $-K_{0}(-x^{\prime},x,t)$ obtained for a Dirichlet boundary condition ($\kappa\to\infty$) at $x=0$ Baltes and Hilf (1976). It thus can be interpreted as the propagation from $x$ to $x^{\prime}$ via the $\delta$-potential, taking a detour or shortcut of length $\epsilon$ weighted with the density $\kappa\mathrm{e}^{-\kappa\epsilon}$.

II.2 The DOS for a billiard with a $\delta$-barrier

The DOS of a $d$-dimensional system can be written as the inverse Laplace transform of the trace of the propagator Hummel et al. (2013):

[TABLE]

Let $\Omega$ be the configuration space of a $d$-dimensional billiard with a classically impenetrable thin barrier inside the volume which can be approximated by a $\delta$-potential along a $(d-1)$-dimensional smooth manifold. By only taking into account short-time propagation inside the billiard we can locally approximate the barrier by $(d-1)$-dimensional planes and thus treat the coordinate perpendicular to the barrier as independent from the remaining $d-1$ tangential coordinates. In this case the local approximation for the propagator can be written as

[TABLE]

Here, $K_{0}^{(d)}(\mathbf{x^{\prime}},\mathbf{x},t)$ is the $d$-dimensional free propagator and $x_{\perp}$ and $\mathbf{x}_{\parallel}$ are the coordinates perpendicular and tangential to the barrier. The trace is now calculated assuming that the integration of the perpendicular direction converges rapidly and is thus independent of the position at the barrier. By introducing the interaction strength in units of energy,

[TABLE]

this yields

[TABLE]

with the $d$-dimensional volume $V_{\Omega}$ of the configuration space $\Omega$ and the surface $S_{\delta}$ of the barrier. The DOS is then given by the inverse Laplace transform of Eq. (15). For $d=1$ we have to set $S_{\delta}=1$ and the DOS is given by

[TABLE]

In all other cases the result is

[TABLE]

A closed formula for the inverse Laplace transform

[TABLE]

for arbitrary dimensions $d>1$ is given in Appendix B.

For varying strength of the $\delta$-potential along the surface of the barrier, i.e., $\kappa=\kappa(\mathbf{x})$, the surface $S_{\delta}$ has to be replaced by the integral operator $\int_{S_{\delta}}\mathop{}\!\mathrm{d}^{d-1}{x}\,$. Note that the boundary conditions at the boundary $\partial\Omega$ of the billiard are not yet included and the approximation of a flat barrier may fail near the boundary. Furthermore, the above approximation does not include curvature corrections and contains only information on the smooth part of the DOS. The result could, in principle, be improved by using periodic orbit theory following Sieber et al. (1995), but this would be at the expense of generality.

Now consider a different setup without a $\delta$-barrier inside the billiard but with Robin- (or mixed) boundary conditions

[TABLE]

In 1D this is equivalent to a $\delta$-potential $(\hbar^{2}\kappa/m)\,\delta(x_{s})$ at the surface (i.e., the end points of the line segment) while only allowing for solutions symmetric to the endpoints in a coordinate space extended beyond the latter. This means that in the approximation of a locally flat surface of the boundary we only have to replace the propagator $K_{\delta}$ in the above derivation by its symmetry-projected equivalent

[TABLE]

while taking the trace perpendicular to the boundary to one side only. One can easily see that this only adds an additional surface term

[TABLE]

to Eq. (17), which corresponds to a Neumann boundary condition, while the other terms remain unchanged. This is due to the fact that the Robin boundary condition is equivalent to a $\delta$-potential on the surface combined with a Neumann boundary condition in the sense of reflection symmetry in the extended space. The above result reproduces the first term in the expansion derived in Balian and Bloch (1971). Note that our approach also comprises the two-dimensional (2D) case, which was treated separately in Balian and Bloch (1971). This will be essential in the next section, where we will use the results to derive the Weyl expansion for an interacting system.

III Two particles on a line segment

III.1 Configuration space

In this section we consider two identical particles on a line with Dirichlet boundary conditions at $q=0$ and $q=L$ [two particles in a box, see Fig. 1(a)] as an idealized model for confined particles. The particles shall interact only when they are at the same point, which is realized by a $\delta$-potential. Furthermore we restrict ourselves to either bosons with zero spin or fermions with spin $1/2$. In order to shorten notation, we use imaginary time and choose scaled units, i.e.,

[TABLE]

Inside the configuration space $\Omega=\{\mathbf{q}\in\mathbb{R}^{2}\mid 0<q_{i}<L\}$ the Hamiltonian is given by

[TABLE]

After introducing relative and center of mass coordinates

[TABLE]

the Hamiltonian reads

[TABLE]

For bosons, the wave functions must be symmetric with respect to particle exchange and, in the case of fermions, an interaction can only occur if the particles have different spin, i.e., the wave functions must be symmetric, too. We can therefore restrict ourselves to bosons. For $\tilde{\kappa}=0$ the system of two indistinguishable particles is equivalent to a system of one quasi-particle of the same mass in the fundamental domain $\mathcal{F}=\{\mathbf{q}\in\Omega\mid q_{1}\geq q_{2}\}$ while requiring a Neumann boundary condition on the line $q_{1}=q_{2}$ Hummel et al. (2013). In the interacting case the same arguments yield the very same equivalence but with a Robin boundary condition

[TABLE]

on the symmetry line instead, as illustrated in Fig. 1(b),(c).

III.2 Smooth part of the DOS

Using the above equivalence of two particles on a line to one quasi-particle in the fundamental domain we can calculate the DOS directly from Eqs. (17) and (21):

[TABLE]

with $\tilde{\mu}=\tilde{\kappa}^{2}$ in scaled units (22). We introduced the abbreviation $\theta^{0}_{1}=(1\mp 1)/2$ here in order to shorten notation. The last term $\rho_{c}(E)$ represents the contributions coming from the corners in the fundamental domain. Here we need the contribution from a $\pi/4$ corner with Dirichlet and Robin boundary conditions along the rays. The exact propagator for such a corner can be derived from the propagator for a $\pi/2$ corner with Robin boundary conditions on the axes,

[TABLE]

The Dirichlet boundary condition at $x_{1}=x_{2}$ can be satisfied by antisymmetrizing $K_{\frac{\pi}{2}}$ with respect to this line. This yields the expression

[TABLE]

The factor $1/2$ was chosen such that the trace can be taken in the first quadrant instead of only integrating the inside of the corner. This is possible due to the symmetry of $K_{\frac{\pi}{2}}$. It is now convenient to write the symmetric propagator (20) as

[TABLE]

with

[TABLE]

Here, $K_{R}(x^{\prime},x,t)$ can be interpreted as the correction to the free propagator representing the propagation from $x$ to the point $x^{\prime}$ via the boundary with mixed boundary condition. The propagator $K_{\frac{\pi}{4}}$ now takes the form

[TABLE]

The interpretation of the different terms for $\mathbf{x}^{\prime}=\mathbf{x}$ is shown in Figs. 2 and 3.

Tracing the above propagator will lead to different contributions to the DOS. The first term can be identified as a volume term, and the third and fourth can be combined to a surface term corresponding to mixed boundary conditions. The surface terms for the Dirichlet boundary are given by the fifth and part of the sixth term. All these contributions have to be dropped to get the parts that arise only from the corner:

[TABLE]

The inverse Laplace transform (12) of the trace of this propagator yields the corner contribution. For the first term in Eq. (33) the trace can be calculated separately for each coordinate and the inverse Laplace transform can be calculated as a convolution of the density

[TABLE]

with itself. The remaining four terms of the propagator (33) are traced as a whole because most of the integrals do not have to be evaluated as they cancel mutually after some elementary manipulations. Altogether, the $\pi/4$ corner contribution is

[TABLE]

Note that all the $\tilde{\mu}$-dependent expressions give multiples of $\delta(E)$ in the limit $\tilde{\mu}\to 0$ (Neumann case), as expected. Combining this result with the contribution from a $\pi/2$ Dirichlet-Dirichlet corner (see, e.g., Eckhardt (1988)) we finally obtain the entire DOS for the system

[TABLE]

The above result can be represented in a shorter form if we consider positive energies only. Then the cases of an attractive and a repulsive interaction only differ in the overall sign of the $\tilde{\mu}$-dependent corner corrections:

[TABLE]

Note that the corner corrections are monotonous and nonzero for $E>0$ but vanish for $E\to\infty$. This means that the corner corrections are most important for energies near the ground state or, in the case of an attractive interaction, close to $E=0$.

In the regime $-\tilde{\mu}<E<0$ the DOS is conveniently expressed in terms of the excitation energy $E^{*}=E+\tilde{\mu}$:

[TABLE]

The first two terms correspond exactly to the DOS for one particle of mass $2m$ on a line with Dirichlet boundary conditions at the end points: for very strong interaction the two particles behave as one single particle. The third term represents the interplay of boundary reflections and the interaction. If integrated from $-\mu$ to [math] it reduces the level counting function by exactly $1/2$ bound state.

III.3 Comparison with numerical calculations

The comparison with the exact quantum mechanical DOS will be done using the level counting function

[TABLE]

In all plots $E$ and $\tilde{\mu}$ will be given in units of $1/L^{2}$ ($\tilde{\kappa}$ in units of $1/L$) which is equivalent to setting $L=1$ in $\mathcal{N}$.

In Fig. 4 the semiclassical level counting function is plotted for $-25\leq\tilde{\kappa}\leq 25$ in steps of $\Delta\tilde{\kappa}=10$. For $E<0$ the three curves representing the attractive cases $\kappa=-25,-15,-5$ resemble 1D single-particle counting functions.

The exact solutions of the Schrödinger equation for the Hamiltonian (25) can be found by using that the system is symmetric with respect to the line $x_{2}=L/\sqrt{2}$. This means that we can determine a full set of energy eigenstates that are either symmetric or antisymmetric to that line [see Fig. 5(a)]. Except for normalization, this is equivalent to considering only the lower part of the fundamental domain while requiring either Neumann or Dirichlet boundary conditions at the symmetry line [see Fig. 5(b)]. The Dirichlet boundary condition at the line $x_{1}=x_{2}$ corresponds to the antisymmetric solutions of the extended system shown in Fig. 5(c).

Moreover, the normalization constants for the original and the extended system are the same. The solutions can be written straightforwardly as antisymmetrized products of 1D wave functions.

[TABLE]

Here, the 1D wave functions and the sets $\text{Spec}^{\mathrm{N/D}}$ are defined by

[TABLE]

where $k_{0}\in\text{Spec}^{\text{D}}$ is chosen either positive or purely imaginary with positive imaginary part depending on the value of $\kappa d$ with $d=L/\sqrt{2}$.

The energy eigenvalues of the system are now given as $E_{mn}=k_{m}^{2}+k_{n}^{2}$ with $k_{m},k_{n}$ either both in $\text{Spec}^{\text{D}}$ or both in $\text{Spec}^{\text{N}}$. The energies have been calculated numerically and the comparison to the smooth level counting function is shown in Fig. 6 for the repulsive case ($\tilde{\kappa}=10$) and in Fig. 7 for the attractive case ($\tilde{\kappa}=-25$).

The value $\tilde{\kappa}=10$ for the repulsive case has been chosen such that the resulting level counting function lies well in between the two limits $\tilde{\kappa}=0$ and $\tilde{\kappa}\to\infty$ corresponding to non-interacting bosons and fermions (denoted as $\mathcal{N}_{\text{N}}$ and $\mathcal{N}_{\text{D}}$ in Fig. 6). For $\tilde{\kappa}\geq 0$ numerical calculations show that the ratio of corner corrections in the DOS to the full semiclassical DOS at the ground state has a maximum of about $8\%$ for $\tilde{\kappa}\approx 5.814$ which is very close to the value at which the equation $E_{0}=\tilde{\kappa}^{2}$ for the ground state energy holds ($E_{0}$ varies smoothly from $E_{0}=2\pi^{2}$ for $\tilde{\kappa}=0$ to $E_{0}=5\pi^{2}$ for $\tilde{\kappa}=\infty$). As this ratio decreases with the energy the corner corrections can be neglected in the DOS for high energies. This holds also true for the level counting function, where they decrease rapidly from $3/8$ to $-1/8$ as the energy increases.

Figure 7 (attractive case) shows that the corner corrections are very important for negative energies and result in a shift by one level in the level counting function at $E=0$. In fact the corner corrections give a contribution that nearly exactly reproduces the quantum mechanical energy eigenstates. By integrating Eq. (38) to get $\mathcal{N}_{-}(E^{*})$ and substituting $k^{*}=\sqrt{E^{*}}$ the equation $\mathcal{N}_{-}(k_{n}^{*})+1/2=n$ for $n\in\mathbb{N}$ yields

[TABLE]

which gives exactly the allowed real wave numbers in $\text{Spec}^{\text{D}}$ and $\text{Spec}^{\text{N}}$. Since the negative quantum-mechanical energy levels are always given as $E=k_{0}^{2}+k_{n}^{2}$ with $k_{0}$ purely imaginary and $k_{n}$ real we can see that the error in the approximation (41) of the eigenenergies lies only in the assumption that the imaginary wave numbers $k_{0}$ in both $\text{Spec}^{\text{D}}$ and $\text{Spec}^{\text{N}}$ are equal to $-\mathrm{i}\tilde{\kappa}$. This approximation is very good for large absolute values of $\tilde{\kappa}$ and is still reasonable as long as the requantization with $\rho_{-}$ makes sense, i.e., the ground state energy is negative.

The error $\Delta E=|k_{0}^{2}-\tilde{\kappa}^{2}|$ of the semiclassically requantized energies is plotted in Fig. 8 for $-20<\tilde{\kappa}<-2$. The ground state energy changes sign at $\tilde{\kappa}=-(3\pi)/(2\sqrt{2})\approx-3.332$ (semiclassical prediction, while the quantum mechanical result yields $\tilde{\kappa}\approx-3.286$). At this point the error is less than $0.45/L^{2}$ which is small compared to the mean level spacing $[\rho_{-}(\tilde{\mu})]^{-1}\approx 19/L^{2}$. This shows that the corner corrections are essential for $E<0$ and can be used to requantize the system in this regime. Moreover one can use the level counting function directly to find the number of bound states in the system.

IV Conclusion

We presented an alternative derivation for the propagator for a $\delta$-potential in Sec. II.1. The resulting expression has a natural interpretation by means of free propagation and hard-wall reflection, taking exponentially weighted detours or shortcuts in the cases of repulsive or attractive interaction, respectively. The propagator can be used to calculate the single-particle DOS for arbitrary shaped billiards of any dimension with mixed boundary conditions and/or $\delta$-barriers inside the volume. Although the calculations do not include curvature terms the formalism is straightforward and allowed us to calculate the contributions from a Dirichlet-Robin $(\pi/4)$-corner to the DOS which, to our knowledge, have not been calculated before. The corrections given by these corners are, as expected, of lowest order in the energy but give very accurate corrections at low positive energies, where the repulsive and the attractive case differ only in an overall sign. In the attractive case the corner corrections allow for the requantization of the system for negative energies, which shows the accuracy of the approximations used in the formalism.

By using the equivalence of the 1D two-particle problem with $\delta$-interactions to the 2D problem in the fundamental domain with mixed boundary conditions we presented a powerful tool for the treatment of few-body systems. Its natural application lies in the approximation of the DOS for an arbitrary number of $\delta$-interacting particles. It is important to note that the scattering and bound states of such a system with infinite volume are known exactly McGuire (1964). This means that the short-time propagators required for the treatment of higher particle numbers can, in principle, be calculated directly from them. Our approach presented here can be extended to include smooth external potentials, which allowed us to examine the thermodynamics of a system of $\delta$-interacting bosons in a harmonic confinement Hummel et al. (2016).

Acknowledgements.

B.G. thanks the Studienstiftung des deutschen Volkes for support. We further acknowledge support through the Deutsche Forschungsgemeinschaft.

Appendix A Continuum limit of the propagator

The antisymmetric solutions of the confined system have equidistant $k$’s, i.e., $\Delta k=k_{n+1}-k_{n}=\pi/L$ for arbitrary $n$. So the first line of Eq. (5) is easily verified to have the continuum limit

[TABLE]

whereas the second line needs a more careful treatment. Omitting the indices we can write

[TABLE]

Using trigonometric identities and identifying

[TABLE]

this yields

[TABLE]

where the absolute values in the first two arguments can be dropped due to symmetry. Now, observing that $A^{2}=1+O(L^{-1})$ and $\Delta k_{n}=k_{n+1}-k_{n}=\pi/L+O(L^{-2})$ for all $n$ as $L$ goes to infinity, the continuum limit can be performed exactly in the same way as for the antisymmetric solutions. The bound state solution of the confined system has an exact limit for the free space, i.e., an exponentially decaying wave function localized at the potential. Summing up all the contributions we get Eq. (6).

Appendix B Inverse Laplace transformation

In Eq. (36) we need the inverse Laplace transform $f_{n}(E)$ of functions of the form

[TABLE]

with $2n\in\mathbb{N}_{0}$. To this end one can use the property of the two-sided Laplace transformation

[TABLE]

to get $f_{n}$ recursively from $f_{0}$ or $f_{\frac{1}{2}}$. We will take a different approach that yields the full expression. It is straightforward to show

[TABLE]

This can be used to verify by induction the formula

[TABLE]

where

[TABLE]

Now one can use the linearity of the Laplace transformation and the identity

[TABLE]

to get

[TABLE]

The functions $f_{\gamma}$ are given by Abramowitz and Stegun (1965)

[TABLE]