Long- and short-range interaction footprints in entanglement entropies of two-particle Wigner molecules in 2D quantum traps

TL;DR

This paper investigates how different interaction ranges affect entanglement entropies in a two-particle 2D quantum trap, revealing divergent behaviors and non-analytical features depending on trap anisotropy and interaction type.

Contribution

It introduces a method to exactly analyze the entanglement spectrum and entropies of two-particle systems with various interaction potentials in 2D traps, especially in the strong interaction limit.

Findings

Long-range interactions lead to finite entropies in anisotropic traps and logarithmic divergence in isotropic traps.

Short-range interactions cause divergence of entanglement measures for any anisotropy due to momentum uncertainty.

Non-analytical behavior of Rènyi entropies occurs when the reduced density matrix has finite support.

Abstract

The occupancies and entropic entanglement measures for the ground state of two particles in a two-dimensional harmonic anisotropic trap are studied. We implement a method to study the large interaction strength limit for different short- and long-range interaction potentials that allows to obtain the exact entanglement spectrum and several entropies. We show that for long-range interactions, the von Neumann, min-entropy and the family of R\'enyi entropies remain finite for the anisotropic traps and diverge logarithmically for the isotropic traps. In the short-range interaction case the entanglement measures diverge for any anisotropic parameter due to the divergence of uncertainty in the momentum since for short-range interactions the relative position width vanishes. We also show that when the reduced density matrix has finite support the R\'enyi entropies present a non-analytical behaviour.

Equations

$$H=-\frac{1}{2}\left({\nabla }_1^2+{\nabla }_2^2\right)+\frac{1}{2}\left\{(x_1^2+x_2^2)+{\epsilon }^2(y_1^2+y_2^2)\right\}+gV\left(r_{12};\left\{{\gamma }_i\right\}\right),$$

$$H^R=-\frac{1}{4}{\nabla }_R^2+\left(X^2+{\epsilon }^2{Y}^2\right),$$

$$H^r=-{\nabla }_r^2+V^{eff}(x,y;\epsilon ,\left\{{\gamma }_i\right\}),$$

$$V^{eff}(x,y;\epsilon ,\left\{{\gamma }_i\right\})=\frac{1}{4}\left(x^2+{\epsilon }^2{y}^2\right)+gV\left(\sqrt{x^2+y^2};\left\{{\gamma }_i\right\}\right).$$

$$\Psi (x,y,X,Y)={\psi }^R(X,Y){\psi }^r(x,y),$$

$$H^R{\psi }^R(\vec{R})=E^R{\psi }^R(\vec{R}),$$

$$H^r{\psi }^r(\vec{r})=E^r{\psi }^r(\vec{r}).$$

$$\vec{r}_{min}=\left(\pm x_{0},0\right)\;\;\;\;\mbox{with $x_{0}>0$ given by}\;\;\;\;\frac{1}{2g}=-\left.\left(\frac{1}{r}\frac{\partial V}{\partial r}\right)\right|_{x_{0}}\,.$$

$$H_{HA}^r=-{\nabla }_r^2+\frac{1}{2}\left\{{\omega }_x^2{\left(x-x_0\right)}^2+\frac{1}{2}\left({\epsilon }^2-1\right)y^2\right\},$$

$${\omega }_x^2=\frac{1}{2}\left(1+\frac{{\frac{{\partial }^{2}V}{\partial r^2}|}_{x_0}}{{-\frac{1}{r}\frac{\partial V}{\partial r}|}_{x_0}}\right),$$

$${\Lambda }_{l,\tilde{l}}={\Lambda }_l^x{\Lambda }_{\tilde{l}}^y,$$

$${\Lambda }_l^x=\frac{\left(1-\zeta \left({\omega }_x\right)\right)}{2\left(1+e^{-\frac{x_0^2{\omega }_x}{\sqrt{2}}}\right)}\zeta {\left({\omega }_x\right)}^l\text{mbox},\zeta \left({\omega }_x\right)={\left(\frac{{\left(2{\omega }_x^2\right)}^{\frac{1}{4}}-1}{{\left(2{\omega }_x^2\right)}^{\frac{1}{4}}+1}\right)}^2,$$

$${\Lambda }_{\tilde{l}}^y=\left(1-\xi (\epsilon )\right)\xi (\epsilon {)}^{\tilde{l}}\text{mbox},\xi (\epsilon )={\left(\frac{{\left({\epsilon }^2-1\right)}^{\frac{1}{4}}-\sqrt{\epsilon }}{{\left({\epsilon }^2-1\right)}^{\frac{1}{4}}+\sqrt{\epsilon }}\right)}^2.$$

$$S^{\alpha }=\frac{1}{1-\alpha }{\log }_2\text{mbox}Tr{\rho }^{\alpha }=\frac{1}{1-\alpha }{\log }_2\left(\sum _i{\Lambda }_i^{\alpha }\right),$$

$$S_{vN}=-\text{mbox}Tr\left(\rho {\log }_2\rho \right)=-\sum _i{\Lambda }_i{\log }_2{\Lambda }_i.$$

$$S_L=1-\text{mbox}Tr{\rho }^2=1-\sum _i{\Lambda }_i^2,$$

$$S^{\alpha }=S_x^{\alpha }({\omega }_x)+S_y^{\alpha }(\epsilon ),$$

$$S_x^{\alpha }\left({\omega }_x\right)=\frac{1}{1-\alpha }{\log }_2\left(\frac{(1-\zeta \left({\omega }_x\right){)}^{\alpha }}{(1-\zeta {\left({\omega }_x\right)}^{\alpha })}\right)+1,$$

$$S_y^{\alpha }(\epsilon )=\frac{1}{1-\alpha }{\log }_2\left(\frac{(1-\xi (\epsilon ){)}^{\alpha }}{(1-\xi (\epsilon {)}^{\alpha })}\right).$$

$$S_{vN}=S_x^1\left({\omega }_x\right)+S_y^1(\epsilon ),$$

$$S_x^1\left({\omega }_x\right)=-\frac{{\log }_2\left({\left(1-\zeta \left({\omega }_x\right)\right)}^{(1-\zeta \left({\omega }_x\right))}\zeta {\left({\omega }_x\right)}^{\zeta \left({\omega }_x\right)}\right)}{\left(1-\zeta \left({\omega }_x\right)\right)}+1,$$

$$S_y^1(\epsilon )=-\frac{{\log }_2\left({\left(1-\xi (\epsilon )\right)}^{(1-\xi (\epsilon ))}\xi (\epsilon {)}^{\xi (\epsilon )}\right)}{\left(1-\xi (\epsilon )\right)}.$$

$$S^{\infty }=\underset{\alpha \to \infty }{\lim }\left(S_x^{\alpha }\left({\omega }_x\right)+S_y^{\alpha }(\epsilon )\right)=\underset{\alpha \to \infty }{\lim }{S}_x^{\alpha }\left({\omega }_x\right)+\underset{\alpha \to \infty }{\lim }{S}_y^{\alpha }(\epsilon )=S_x^{\infty }\left({\omega }_x\right)+S_y^{\infty }(\epsilon ).$$

$$S_L=1-\frac{1}{2}\frac{1-\zeta \left({\omega }_x\right)}{1+\zeta \left({\omega }_x\right)}\frac{1-\xi (\epsilon )}{1+\xi (\epsilon )}.$$

$$S_{vN}\sim -\frac{\ln (\epsilon -1)}{\ln 16}\text{mbox}for\epsilon \sim 1^{+}.$$

$$V^{ip}\left(r;\beta \right)=\frac{1}{r^{2\beta }}.$$

$$x_0=(4g\beta {)}^{\frac{1}{2(\beta +1)}}\text{mbox}and{\omega }_x^2=\beta +1.$$

$$V^{il}\left(r\right)=\frac{1}{\ln (r+1)}.$$

$$2g=x_0\left(x_0+1\right){\ln }^2\left(x_0+1\right)\text{mbox}and{\omega }_x^2=\frac{1}{2}\left\{1+\left(\frac{\frac{2}{\ln \left(x_0+1\right)}+1}{\frac{1}{x_0}+1}\right)\right\}.$$

$$V^{sip}\left(r;\{\beta ,\gamma \}\right)=\frac{e^{-\gamma r}}{r^{2\beta }},$$

$$2g=\frac{e^{\gamma x_0}{x}_0^{2(1+\beta )}}{2\beta +\gamma x_0}\text{mbox}and{\omega }_x^2=\frac{1}{2}\left(1+\frac{2\beta }{2\beta +\gamma x_0}+2\beta +\gamma x_0\right).$$

Full text

I Introduction

The physics of systems of confined particles has attracted the interest of researchers from many different areas working on both, theoretical and experimental aspects wigner_1934 ; wineland_1987 . One feature of such quantum systems that has recently gained impulse is the study of entropic measures of entanglement osenda_2015 ; garagiola_2016 ; koscik_2015 ; lopezrosa_2015 ; simonovic_2015 . Among the many kind of systems that could be addressed using the physics of the confined systems, in the present work we focus on Wigner molecules, which are the finite-size analogue of Wigner crystals, named after the seminal work of E. Wigner wigner_1934 . Since the late 80’s, when the first confined linear chains of ions where reported by D. Wineland wineland_1987 , there has been an increasing capacity to confine, control and manipulate such entities, and has turned the Wigner molecules into a new platform to test the known, and look for new, traits of quantum mechanics drewsen_2015 .

Trapped ions are not the only physical systems that allow the formation of entities like the Wigner molecules (also named Coulomb crystals drewsen_2015 ). The experimental observation of strongly correlated states in quantum dots has attracted considerable interest cavaliere_2009 . Wigner molecules have also been observed in two-dimensional semiconductor heteroestructures andrei_1998 ; piot_2008 , semiconductor quantum dots kallikaos_2008 , one-dimensional quantum wires ellenberger_2006 ; singha_2010 ; kristinsdottir_2011 ; meyer_2009 , carbon nanotubes deshpande_2008 ; pecker_2013 , and in crystalline states for dusty plasma melzer_2003 . Several theoretical studies jauregui_1993 ; guclu_2008 ; mendl_2014 ; gambeta_2014 ; cavaliere_2014 ; cavaliere_2015 ; cavaliere_2015_2 ; kylanpaa_2016 have demonstrated that the physics of these systems with reduced dimensionality usually does not depend on the shape of the confinement but on its symmetries and strength (see Ref. pecak_2016 for an example where shape does have an influence). In the present work we use a harmonic confinement as a model potential, and represent different physical situations using different interaction potentials. Wigner molecules arise when the interparticle interaction strength is much larger than the kinetic energy. The latter can be related to the temperature of the system drewsen_2015 and also to the density, or confining energy, of the particles filinov_2001 . Hence Wigner localization is expected for low density systems or for large interaction strengths.

The Calogero and Moshinsky models are the most salient examples of analytically solvable models of confined particles, including the exact computation of the entanglement entropies osenda_2015 ; garagiola_2016 ; manzano_2010 ; dehesa_2012 ; benavides_2014 ; koscik_2010 ; koscik_2013 ; koscik_saha_2015 , in this sense it can also be mentioned the spherium model toranzo_2015 , and the quasi-solvable Hook model taut_1993 . In particular, the Calogero model has been widely studied in condensed matter physics and has experienced several revivals Polychronakos_1993 ; Polychronakos_2006 , such as the discovery of an explicit relation of the Calogero model with the fractionary quantum hall effect Azuma_1994 and fractional statistics Murthy_1994 . In a previous paper garagiola_2016 , we have studied the behavior of the von Neumann and Rényi entropies of the one- and two-dimensional Calogero model for two particles. By considering anisotropic confinement in the two-dimensional case we showed that the one-dimensional regime is reached when the anisotropy of the trap increases, and we also demonstrated that the Rényi entropies present a non analytical behavior in the neighborhood of those values of the interaction strength parameter for which the reduced density matrix has finite support.

Motivated by this, we consider anisotropic harmonic confinement and compute the exact expression for the occupation numbers or occupancies of the two-dimensional ground state wave function in the large interaction strength limit for two particles which interact via different potentials depending on the distance between the particles. The exact natural orbitals are obtained from the Schmidt decomposition of the ground state wave function in the same limit and the occupancies are used to evaluate several quantum information measurements such as von Neumann and Rényi entropies in closed form. The method presented here is a generalization of the strategy developed in Refs. koscik_2015 ; koscik_2010 ; koscik_2015_2 ; glasser_2013 . The two particle one-dimensional systems with Coulomb and inverse powers interactions are addressed in koscik_2015 ; koscik_2015_2 ; glasser_2013 , while the natural orbitals and occupation numbers of elliptically deformed two-dimensional quantum dots are reported in koscik_2010 . Here we give the analytical expressions of the natural orbitals, occupation numbers, von Neumann and Rényi entropies in the strong interaction limit for any potential which depends only on the interparticle distance.

Our main purpose is to determine the influence of the anisotropy and the type of interparticle interaction by looking upon the linear, von Neumann, min-entropy, max-entropy and Rényi entropies as entanglement measures. We have a particular interest in which are the differences arising from a short-range interaction with respect to a long range one for which the emergence of Wigner Molecules has been widely described (see, for example, Refs. jauregui_1993 ; guclu_2008 ; cavaliere_2014 ). With this aim, we study two interacting potential cases for each interaction range, including one that can be exactly solved. In the long-range interaction case we consider the inverse power and the inverse logarithmic potential, and for the short-range interaction we solve the screened inverse power potential and a Gaussian repulsive interaction. It is important to emphasize that the inverse power interaction case is used to model quantum dots kallikaos_2008 or ion traps wineland_1987 where the large interaction regime can be achieved experimentally due to a strong interaction between the particles or a weak confinement energy scale, for inverse square power one gets the Calogero model, while the screened Coulomb interaction provides a simple model potential for ions and plasmas li_2012 .

The paper is organized as follows. The model is discussed in Section II. In Section III we show the derivation of the analytical occupancies of two interacting particles in a two-dimensional anisotropic harmonic trap, while in Section IV we calculate the entropic entanglement measures. We discuss the results for long- and short-range interaction potentials in Sections V and VI respectively. Finally, a summary and conclusions are presented in Section VII.

II Confined two-dimensional two-particle systems

The physics of confined particle systems is nowadays very relevant to understand the many recent experiments conducted in cold atom traps or in quantum dots, at least in a qualitative way andrei_1998 ; piot_2008 ; kallikaos_2008 ; ellenberger_2006 ; singha_2010 ; kristinsdottir_2011 ; meyer_2009 . The models for those systems contain two contributions to the potential energy: one is given by the trap potential and the other by the interaction between the particles. For small dots, containing few electrons the trap potential can be approximated by a harmonic one cavaliere_2009 , therefore we focus here on two interacting particles in a two-dimensional anisotropic harmonic traps, and implement a method to solve the entanglement spectrum in the large interaction limit for arbitrary interaction potentials. The Hamiltonian for two particles in an anisotropic trap, in atomic units, is

[TABLE]

where the frequency of the trap was taken equals to unity, $\varepsilon>1$ is the anisotropy parameter, $V\left(r_{12};\left\{\gamma_{i}\right\}\right)$ denote the interaction potential as a function of the interparticle distance $r_{12}$ and some parameters $\left\{\gamma_{i}\right\}$, and $g$ is the ratio between the interaction and the confinement energy scale. By introducing the center of mass $\vec{R}=\frac{1}{2}(\vec{r}_{1}+\vec{r}_{2})=(X,Y)$ and relative coordinates $\vec{r}=\vec{r}_{2}-\vec{r}_{1}=(x,y)$ the Hamiltonian (1) decouple as $H=H^{R}+H^{r}$, where

[TABLE]

and $V^{eff}$ is the effective potential of the relative Hamiltonian given by

[TABLE]

The total wave function is then the product of the center of mass wave function and the relative wave function

[TABLE]

and, consequently, the Schrödinger equation separates into two equations

[TABLE]

The solutions of the center of mass equation (Eq. (6)) are the eigenfunctions of the harmonic oscillator that are invariant under particle exchange.

The relative Hamiltonian, Eqs. (3) and (4), must be analysed on a case-specific basis. However, in the next section we present a method to obtain the large interaction strength limit of general potentials that fulfil simple requirements.

III Derivation of the analytical occupancies

The relative wave function may be obtained by solving the Schrödinger equation in the large interaction strength regime, $g\gg 1$, by using the harmonic approximation (HA) james_1998 ; balzer_2006 . In the framework of the harmonic approximation one has to find the minima of the effective potential Eq. (4) and then the potential is replaced by its Taylor expansion up to second order about its minima, which satisfy $\nabla V^{eff}(x,y;\varepsilon,\left\{\gamma_{i}\right\})=0$. If the potential is repulsive, decreases monotonously and $V\left(r;\left\{\gamma_{i}\right\}\right)\to 0$ for $r\to\infty$, with $\varepsilon>1$, the minima lie on the $x-$axis and can be written as

[TABLE]

It is important to notice that when the particles are confined in an isotropic trap, i.e $\varepsilon=1$, the minima degenerate into a circle of radius $x_{0}$.

Within the harmonic approximation, a Hamiltonian of uncoupled oscillators is obtained

[TABLE]

with a frequency associated to the $x$-coordinate given by

[TABLE]

where the dependence on the parameters $g$ and $\left\{\gamma_{i}\right\}$ is implicit in $x_{0}=x_{0}\left(g,\left\{\gamma_{i}\right\}\right)$.

The totally symmetric ground state wave function $\Psi^{GS}\left(\vec{r}_{1},\vec{r}_{2}\right)$ of the harmonic Hamiltonian $H^{R}+H^{r}_{HA}$ is a product of Gaussians. The eigenvalues of the one-particle reduced density matrix $\rho=\textrm{Tr}_{2}\left(\left|\Psi^{GS}\right\rangle\left\langle\Psi^{GS}\right|\right)$ are explicitly obtained (see supporting information), and are given by

[TABLE]

where

[TABLE]

and

[TABLE]

where $l,\tilde{l}=0,1,2,\ldots$. Each eigenvalue, or occupancy, is doubly degenerate due to the particle exchange symmetry.

The limiting values and behavior of $\zeta\left(\omega_{x}\right)$ and $\xi(\varepsilon)$ are needed to compute the entropic quantities. We note then that for $\omega_{x}>0$, $\zeta\left(\omega_{x}\right)$ is always below unity and $\zeta\left(\omega_{x}\right)\to 1$ when $\omega_{x}\to\infty$, while for $\varepsilon>1$, $\xi(\varepsilon)$ remains below one and $\xi(\varepsilon)\to 1$ for $\varepsilon\to 1^{+}$, then we must be specially careful in the isotropic confinement case (see Eq. (13)). For large anisotropy parameter $\varepsilon\gg 1$ one gets $\xi(\varepsilon)\to 0$ and the occupancies reach the asymptotic values of the one dimensional model $\Lambda_{l}^{x}$.

IV Entropies in the large interaction strength limit

The entanglement can be measured using different entropic quantities. If $\{\Lambda_{i}\}$ is the complete set of eigenvalues, then the Rényi entropies are a family of such entropies defined by

[TABLE]

which are widely used in many-body or extended systems Calabrese2011a ; Alba2010 . Special values of the parameter $\alpha$ allow to recover other entropies, being the min- and max-entropy good examples obtained by taking the limits $\alpha\rightarrow\infty$ and $\alpha\rightarrow 0$, respectively. The min-entropy serves as a lower bound to the entanglement measures obtained from the whole family of entropies. The Hartley or max-entropy, $S^{0}=\log_{2}R$, only depends on the Schmidt rank $R$ of the spectrum distribution and is a measure of bipartite entanglement which serves as a criterion for efficient classical representation of the state amico_2014 . The distribution of the entanglement spectrum can be better understood by computing the Rényi entropies for many different values of the parameter $\alpha$ garagiola_2016 . The von Neumann entropy is given by

[TABLE]

The von Neumann entropy has been used to study entanglement in continuous variables systems and spin models Alba2010 ; osenda_2007 ; pont_2010 . It can be recovered from the Rényi entropies in the limit $\alpha\to 1$. Finally, some authors use the linear entropy, defined by

[TABLE]

since for continuous variable systems the calculation of $\mbox{Tr}\,\rho^{2}$ is reduced to a single integral. Even though, the linear entropy has no relevant information for the systems studied in the present work, we compute it for the sake of completeness.

Once we have obtained the occupancies it is possible to calculate the quantum entropies. Since these calculations involve geometric series in $\zeta\left(\omega_{x}\right)$ and $\xi(\varepsilon)$, the limiting values must be carefully computed.

Let us start with the Rényi entropies defined by Eq. (14). It is straightforward to show that due to the separability of the wave function, the Rényi entropies in the large interaction strength limit are the sum of the entropy associated to $\psi_{x}(x_{1},x_{2})$ and $\psi_{y}(y_{1},y_{2})$, then

[TABLE]

where

[TABLE]

and

[TABLE]

Again, due to the separability of the wave function, we can write the two-dimensional von Neumann entropy of Eq. (15) as

[TABLE]

where each one of the terms in the sum has the form of a one-dimensional von Neumann entropy koscik_2015 , i.e.

[TABLE]

[TABLE]

The super-index points that the von Neumann entropy can be obtained as a limiting case of the Rényi entropies when $\alpha\to 1$.

It is worth to notice that from Eq. (18) and (19) it is straightforward to show that the min-entropy $S^{\infty}$, can also be written as a two-term sum:

[TABLE]

The Hartley or max-entropy in the large interaction strength limit can also be calculated as a limiting case with $\alpha\to 0$, $S^{0}=\log_{2}R$, and has finite value only when the one-particle reduced density matrix has finite support.

The two-dimensional linear entropy defined by Eq. (16) gives

[TABLE]

For the isotropic model $\varepsilon\to 1^{+}$, $\xi(\varepsilon)\to 1$ and the linear entropy goes to one, while for any other value of $\varepsilon$ the linear entropy remains below one.

A comment on the extension of the previous results to dimension $D$ is in place. They can be extended if one considers $D-1$ anisotropy parameters (see supporting information). The von Neumann and Rényi entropies are the sum of $D$ terms each one associated to one cartesian coordinate and, as we demonstrated for the two-dimensional case, the $x$-entropy term depends on the parameters of the interaction potential through $\omega_{x}$ and each one of the remaining terms depend on only one of the $D-1$ anisotropy parameters.

The two terms of the two-dimensional von Neumann, min-entropy and Rényi entropies ($\alpha=0.2,\,0.4,\,0.8,\,1.5,\,2$) are depicted in Fig. 1. Let us first discuss the behavior of the entropies with respect to the anisotropy of the trap, and afterwards the influence of the interparticle interaction.

As can be appreciated in Fig. 1(a), for the isotropic model ($\varepsilon\to 1^{+}$) the entropies $S^{\alpha}_{y}(\varepsilon)$ diverge logarithmically, whilst for any other values of $\varepsilon$ they remain finite. By calculating the first derivative of Eq. (22) it is straightforward to show that

[TABLE]

This asymptotic leading term is depicted in the figure as a yellow dashed line which makes the logarithmic divergence of the isotropic von Neumann entropy evident. Actually, for $\varepsilon\to 1^{+}$ the von Neumann, min-entropy and the family of Rényi entropies present this same behaviour. The figure also shows that for large anisotropy parameter the entropies $S^{\alpha}_{y}(\varepsilon)$ vanish. In other words, for $\varepsilon\gg 1$ the one dimensional problem is recovered and the von Neumann, Rényi and min-entropy reach the one dimensional values $S_{x}^{\alpha}(\omega_{x})$.

The behavior of the $x$-entropies (denoted by $S^{\alpha}_{x}(\omega_{x})$) as a function of the frequency is shown in Fig. 1(b). The figure shows that the entropies are decreasing functions of the frequency for $0<\omega_{x}^{2}<1/2$, and increasing functions for $\omega_{x}^{2}>1/2$. Actually, the entropies diverge logarithmically for large frequencies and also for $\omega_{x}\to 0$, because in these limits one gets that $\zeta(\omega_{x})\to 1$.

The entropy of a given system is computed using the frequency $\omega_{x}$ obtained by the harmonic approximation Eq. (10). If it remains finite for large interactions parameters $g\gg 1$, then the von Neumann, min-entropy and the family of Rényi entropies are finite for the anisotropic model and diverge logarithmically for the isotropic model. In the deformed or anisotropic case the particles crystallize around the two classical minima of the relative Hamiltonian giving rise to a Wigner molecule, while for the isotropic model those minima degenerate into a circle, the particles are no longer localized around discrete minima and this lack of information is reflected in the divergence of the entanglement entropies. If the obtained frequency increases monotonously for large interactions, the von Neumann, min-entropy and the family of Rényi entropies diverge logarithmically for any anisotropy parameter. In this sense, the behaviour of the system is defined by the one-dimensional entropy $S_{x}^{\alpha}(\omega_{x})$.

The previous analysis can be understood more qualitatively by using the Heisenberg uncertainty principle. The width of the Gaussian wave packet in the relative coordinate $\psi^{r}(\vec{r})$ (ground state of the Hamiltonian Eq. (9)) goes to zero when the frequency increases. Actually, the relative position and momentum uncertainty are $\Delta x^{r}_{HA}=\sqrt{\langle\left(x_{2}-x_{1}\right)^{2}\rangle-\langle x_{2}-x_{1}\rangle^{2}}=2^{\frac{1}{4}}/\sqrt{\omega_{x}}$ and $\Delta p^{r}_{HA}=\sqrt{\omega_{x}}/2^{\frac{5}{4}}$, then if $\omega_{x}\to 0$ we obtain that $\Delta x^{r}_{HA}\to\infty$ and, conversely, when $\omega_{x}\to\infty$ it is straightforward to show that $\Delta p^{r}_{HA}\to\infty$. Thus, we see that for $\omega_{x}\to\infty$ the entropy of the Wigner molecule diverge because the position is completely determined and consequently the momentum uncertainty diverges, we refer to this limit as strong crystallization. The opposite case, $\omega_{x}\to 0$ leads to a well defined momentum state and hence we have no knowledge of the position. In both cases the divergence in the position or momentum width leads to the divergence of the entanglement entropies. Furthermore, the divergence of the $y$-entropies could also be explained in a similar way: for the isotropic model the minima degenerate into a circle and consequently the particles are no longer localized around any definite angular positions, but the state has definite angular momentum.

For $\omega_{x}^{2}=1/2$ the entropies have their minimum value equal to unity. Around this point, the von Neumann, min-entropy and Rényi entropies with $\alpha>1$ present an analytical behaviour while the Rényi entropies with $\alpha<1$ have a non-analytical behaviour. The von Neumann and Rényi entropies with $\alpha=0.4,\,0.5,\,0.6$ and their first derivatives around the point $\omega_{x}^{2}=1/2$ are shown in Fig. 2 (a) and (b) respectively. It shows that the Rényi entropies present an infinite derivative for $\alpha=0.4$, discontinuous derivative for $\alpha=0.5$ and a continuous derivative with infinite second derivative for $\alpha=0.6$, while the von Neumann entropy ($\alpha\rightarrow{}^{+}1$) is an analytical function of the frequency.

Recent studies by Amico and co-workers in $1/2$-spin chains show the physical implications of non-monotonous properties of the Rényi entropies in many-body systems with topological order due to a truncation of the support of the reduced density matrix amico_2013 ; amico_2013_2 ; amico_2014_2 . In Ref. osenda_2015 some of the present authors found that the Calogero model in one dimension has a finite number of non-zero occupancies for a discrete set of values of the interaction parameter, and in Ref. garagiola_2016 we demonstrated that in those particular values of the interaction parameter the Rényi entropies present a non-analytical behaviour.

Summarizing, non-analytical behaviour of the Rényi entropies exposes the finite support of the reduced density matrix. In the present case, taking $\omega_{x}^{2}=1/2$ in Eq. (12) it is straightforward to see that for this particular frequency there is only two non vanishing occupancies $\Lambda^{x}_{0}$ associated to the two lowest natural orbitals in the $x$-coordinate.

In the following sections we apply our findings to study the behavior of the occupancies and entropic entanglement measures in the large interaction strength limit for different cases divided as long- or short-range potentials. From now on we calculate only the one dimensional entropy $S_{x}^{\alpha}(\omega_{x})$, since the behaviour of the entropy terms $S^{\alpha}_{y}(\varepsilon)$ were already analysed.

V Long-range interaction potentials

In the present section we consider two long-range interactions to exemplify our results: the inverse power interaction and inverse logarithmic interaction.

V.1 Inverse power interaction

The inverse power potential is

[TABLE]

For this potential $x_{0}$ and $\omega_{x}$, Eqs. (8) and (10), can be obtained exactly and give

[TABLE]

Thus, $x_{0}$ increases when increasing the interaction strength parameter $g$, but the frequency remains invariant. For $\beta=\frac{1}{2},1$ one gets the Hook and the Calogero model respectively.

Let us start with the Calogero model. The occupancies defined in Eq. (11) for large interaction strength parameter $g\gg 1$ are shown as a function of $\ln(\varepsilon-1)$ in Fig. 3, where the grey dashed line is the dominant one-dimensional occupancy in the large interaction strength limit, obtained from Eq. (12), that is in agreement with the value reported in Ref. koscik_2015 . The figure shows that for $\epsilon\to 1^{+}$ (isotropic model) the occupancies go to zero, but note that their sum is always equal to $1/2$ due to the mentioned double degeneracy koscik_2010 . When the anisotropy increases all the occupancies $\Lambda_{l,\tilde{l}}$ with $\tilde{l}\neq 0$ present a local maximum. For fixed $l$ the value of the anisotropy parameter at which the maximum occurs decreases when $\tilde{l}$ increases, while for fixed $\tilde{l}$ this value is the same for each $l$. For $\varepsilon\gg 1$ the Hamiltonian reduces to a one dimensional oscillator and the occupancies $\Lambda_{l,0}$ reach the asymptotic values of the one dimensional model. For values of $\varepsilon$ near $\varepsilon_{c}=\sqrt{5}$ the occupancies with $\tilde{l}=0$ stabilize on the one dimensional values and those with $\tilde{l}\neq 0$ saturate at vanishingly small values. This feature can be explained if one takes into account that for $\varepsilon=\varepsilon_{c}$ the relative Hamiltonian Eq. (9) reduces to a harmonic oscillator in polar coordinates around each minimum. More generally, for arbitrary $\beta$, the one dimensional regime is reached at the value $\varepsilon_{c}=\sqrt{2\left(\beta+1\right)+1}$. In this case the effective potential of the relative Hamiltonian Eq. (4) is isotropic in a small neighborhood around its minima. For $\varepsilon>\varepsilon_{c}$, the largest occupancy $\Lambda_{00}$ reaches the value $\sim 0.4853$, and the sum of all the remaining occupancies is only $\sim 0.0147$; this means that the two natural orbitals associated to this eigenvalue are the only two that are occupied while all the others natural orbitals contribution are negligible, and consequently, the spatial wave functions are quite similar to those two natural orbitals koscik_2010 .

As we mentioned above, the dependence with $\beta$ is present only through the entropy term $S^{\alpha}_{x}(\omega_{x})$. The width of the Gaussian wave packet in the relative coordinate is finite, $\Delta x^{r}_{HA}=2^{\frac{1}{4}}/(\beta+1)^{\frac{1}{4}}$ and consequently the von Neumann, min-entropy and Rényi entropies are finite. However, the max-entropy diverges due to the infinite support of the one-particle density matrix. Notice that in the limit $\beta\to\infty$ the entropies diverge logarithmically due to the divergence in the momentum uncertainty. This behavior can be seen in Fig. 4 where the von Neumann, the min-entropy and Rényi entropies are depicted as a function of the parameter $\beta$. The Rényi entropies increases for decreasing $\alpha$, and the von Neumann entropy is a limiting case with $\alpha\to 1$. It is important to emphasize that taking the limit $\beta\to 0$ in the entropies does not result in the same entropies obtained for a system with harmonic confinement and a constant interaction (Eq. (26) with $\beta=0$), since this limit does not commute with the large interaction limit.

V.2 Inverse logarithmic interaction

The potential for inverse logarithmic interparticle interaction is

[TABLE]

In this case $x_{0}$ and $\omega_{x}$ satisfy the following equations

[TABLE]

For large interaction strength parameter the value of $x_{0}$ increases when $g$ increases, and consequently, the frequency goes to unity. Therefore, for $g\gg 1$, the one-dimensional von Neumann and Rényi entropies with $\alpha>0$ remain finite, but once more the max-entropy diverges. We included two figures in the supporting information showing the qualitative behavior of $x_{0}$ and the Rényi entropy as a function of $\alpha$. The same analysis performed in the inverse power interaction case can be done for the inverse logarithmic potential.

VI Short-range interaction potentials

In the present section we consider two particles in a two-dimensional anisotropic harmonic trap with two different short-range interactions: the screened inverse power interaction and Gaussian repulsive interaction.

VI.1 The screened inverse power interaction

For the screened inverse power interaction the potential is

[TABLE]

where $1/\gamma$ is the cut-off distance. In this case $x_{0}$ and $\omega_{x}$ are given by

[TABLE]

Notice that taking $\gamma=0$ the minima and the frequency of the inverse power interaction is recovered, and for $\beta=0$ the interaction has exponential decay. For large interaction strength parameter $g$ the minima and the frequency increase monotonously with $g$, and consequently, the one-dimensional von Neumann and Rényi entropies diverge logarithmically. As we mentioned in section IV the divergence of the entropies can be explained as arising from the momentum uncertainty $\Delta p^{r}_{HA}=\sqrt{\omega_{x}}/2^{\frac{5}{4}}$, which diverges when $\omega_{x}\to\infty$. Actually, the larger the $\gamma$ parameter is, the larger the frequency is and the higher the entanglement entropies are, this behaviour is shown in Fig. 5 where the one-dimensional von Neumann entropy is depicted as a function of the interaction strength for $\beta=1$ and $\gamma=0,\,1/2,\,1,\;2$.

VI.2 The Gaussian repulsive interaction

In this subsection we consider the following interaction potential

[TABLE]

where $\sigma$ is the half width of the potential. In this case $x_{0}$ and $\omega_{x}$ can be found exactly

[TABLE]

They are increasing functions of the interaction strength parameter $g$, thus, the one-dimensional von Neumann and Rényi entropies in the large interaction strength limit diverge logarithmically. We interpret this divergence in the same way as for the screened inverse power interaction. It is worth to mention that the limit $\sigma\to 0$ does not reproduce the results of a delta interaction, which have a finite von Neumann entropy avakian_1987 , since this limit does not commute with the large interaction strength limit.

From Eq. (33) it is straightforward to show that for $g=g_{c}$ with

[TABLE]

$x_{0}=\sigma$ and $\omega_{x}^{2}=\frac{1}{2}$, therefore, as was explained in section III all the occupancies vanish except two of them. The reduced density matrix has then a finite support and the Rényi entropies with $\alpha<1$ have a non-analytical behaviour, while the von Neumann entropy and Rényi entropies with $\alpha>1$ present a minimum at $g=g_{c}$. The behaviour of the first and second degenerate occupancies (Eq. (12)) for $\sigma=10$ is depicted in Fig. 6 together with the value $g_{c}$ as a gray dashed line. The first two occupancies $\Lambda_{0}^{x}$ (Fig. 6(a)) reach the maximum value $1/2$ for $g=g_{c}$, value at which all the others occupancies vanish as can be appreciated for $\Lambda_{1}^{x}$ in Fig. 6(b).

Summarizing, we found that two trapped particles with a Gaussian repulsive interaction between them, have a reduced density matrix with infinite support (infinite non-vanishing occupancies) for all the interaction strengths but for $g=g_{c}$, value at which all the occupancies vanish except two, and the support is finite. Nevertheless, it is important to notice that the largest occupancy $\Lambda_{0}^{x}\gtrsim 0.49$ throughout the considered range of $g$ values, and the sum of all the remaining occupancies is $\lesssim 0.01$. Therefore, in the neighborhood of $g_{c}$, the two natural orbitals associated to $\Lambda_{0}^{x}$ are the main contributions to the spatial wave function expansion.

VII Summary and Conclusions

In this work we present analytical expressions in the large interaction strength limit for the occupancies and quantum entropies for the ground state of a two-particle Wigner molecule in a two-dimensional anisotropic harmonic trap. Our main result is that one is able to determine the influence of the anisotropy and the range of the interparticle interaction looking upon the entropic entanglement measures.

The wave function is obtained within the framework of the harmonic approximation for large interaction strength values, and once we have the ground state wave function, we calculate the occupancies from the Schmidt decomposition of the reduced density matrix. We obtain doubly degenerate occupancies, and relate this to the equivalence between particle exchange in the wave function and the exchange between the two minima of the effective potential of the relative Hamiltonian. The linear, von Neumann, min-entropy, max-entropy and Rényi entropies are calculated exactly in terms of the occupancies as a function of the anisotropy parameter and the parameters of the interaction potential.

We found that, due to the coordinate separability of the wave function, the von Neumann, min-entropy, max-entropy and Rényi entropies are a sum of terms associated to each coordinate, and that only one of these terms depends on the anisotropy parameter and the other term is associated to the interaction potential. Consequently, the behavior of the entropies with respect to the anisotropy parameter can be analyzed without regard of the interaction potential, and the dependence on the interaction potential is entirely defined by the frequency obtained by the harmonic approximation of the one-dimensional problem. Moreover, we generalize these results to dimensions higher than two, see details in the supporting information.

We show that when the frequency remains finite for large interactions, then the von Neumann, min-entropy and the family of Rényi entropies remain finite for the anisotropic model and diverge logarithmically for the isotropic model. The divergence of the entanglement measure entropies of the isotropic model can be understood as follows: in the deformed or anisotropic case the particles locate around the two classical minima of the relative Hamiltonian forming a Wigner molecule, while for the isotropic model those minima degenerate into a circle, the particles are no longer localized and this lack of information is reflected in the divergence of the entanglement entropies. If the frequency increases monotonously for large interactions then, the von Neumann, min-entropy and the family of Rényi entropies diverge logarithmically for any anisotropy parameter. In this sense, the influence of the interaction potential is present only in the one-dimensional entropies.

The previous interaction-independent analysis allows us to apply them to different interactions straightforwardly. We group the interactions into short and long-range potentials and show the differences of the obtained results between the groups. For long-range interaction potentials, the frequencies remain finite in the large interaction strength limit, and the von Neumann, min-entropy and Rényi entropies are finite. In contradistinction, for short-range interaction potentials, the frequencies increase monotonously as a function of the interaction strength and, consequently, the one-dimensional von Neumann and Rényi entropies diverge in the large interaction strength limit. The divergence of the entanglement entropies can be explained as arising from the momentum uncertainty divergence at large frequencies. It is important to mention that the one-dimensional von Neumann, min-entropy and Rényi entropies of the inverse power interaction model diverge logarithmically when the power of the inverse interaction increases (see Eqs. (18) and (21)), since in this limit the interaction between the particles goes to a short range one.

We also demonstrate that when the frequency associated to the interaction potential, satisfy $\omega_{x}^{2}=1/2$, the entropies have their minimum value equal to unity. Actually, the von Neumann, min-entropy and Rényi entropies with $\alpha>1$ present an analytical behaviour around this point, while the Rényi entropies with $\alpha<1$ have a non-analytical behaviour, which exposes the finite support of the reduced density matrix. For this particular frequency there is only two non vanishing occupancies $\Lambda^{x}_{0}$. Similar features were also recently reported by Amico and co-workers for $1/2$-spin chains amico_2013 ; amico_2013_2 ; amico_2014_2 , and for the Calogero model garagiola_2016 by the present authors. We further illustrate these features showing that two trapped particles with a Gaussian repulsive interaction have a reduced density matrix with infinite support (infinite non-vanishing occupancies) for all Hamiltonian parameters except for those values that allow $\omega_{x}^{2}=1/2$, where all the occupancies vanish except two, and the reduced density matrix has finite support.

As a final remark, there is a very recent work concerning a system of two Coulombically interacting particles confined to a $D-1$ sphere, where the dependence of the entanglement measures on the radius of the system and the spatial dimensionality has been investigated toranzo_2015 . Thus, as future perspectives we would like to study the effects of the dimensionality and the interaction strength on the entanglement of two confined particles which interact via a general potential, taking as a starting point the results obtained in the generalization to dimensions higher than two presented in the second section of the supporting information.

Acknowledgments

We acknowledge SECYT-UNC and CONICET for partial financial support of this project. E.C. would like to thank Alvaro Cuestas for an exhaustive reading of the manuscript.

I Derivation of the analytical occupancies and natural orbitals

We obtain the relative wave function by solving the Schrödinger equation in the large interaction strength regime, within the harmonic approximation (HA) sjames_1998 ; sbalzer_2006 . As is pointed in the introduction, the method presented here is a generalization of the strategy developed in Ref. skoscik_2015 ; skoscik_2010 ; skoscik_2015_2 ; sglasser_2013 for some particular interaction potentials (Coulomb and inverse powers) to any interaction which depends only on the interparticle interaction. As is also mentioned in the main text, if the potential is repulsive, decreases monotonously and $V\left(r;\left\{\gamma_{i}\right\}\right)\to 0$ for $r\to\infty$, with $\varepsilon>1$, the minima lie on the $x-$axis and can be written as

[TABLE]

Within the harmonic approximation, a Hamiltonian of uncoupled oscillators is obtained

[TABLE]

with a frequency associated to the $x$-coordinate given by

[TABLE]

By using Eq. (1), the frequency can be rewritten as

[TABLE]

where the dependence on the parameters $g$ and $\left\{\gamma_{i}\right\}$ is implicit in $x_{0}=x_{0}\left(g,\left\{\gamma_{i}\right\}\right)$.

The solutions to the corresponding Schrödinger equation are

[TABLE]

with energies

[TABLE]

By using the solutions of the center of mass equation, given by

[TABLE]

with energies

[TABLE]

and (5) taking $n=m=\tilde{n}=\tilde{m}=0$, the totally symmetric ground state wave function is

[TABLE]

where $C$ is the normalization constant

[TABLE]

The total wave function, Eq. (9), is separable in the $x$ and $y$ coordinates as $\Psi^{GS}\left(\vec{r}_{1},\vec{r}_{2}\right)=\psi_{x}(x_{1},x_{2})\psi_{y}(y_{1},y_{2})$ where,

[TABLE]

with

[TABLE]

and

[TABLE]

where $C_{x}$ and $C_{y}$ are the first and second factors of Eq. (10) respectively.

Since we are interested in the occupancies of the natural orbitals, we solve the integral equation for the eigenvalues of the one-particle reduced density matrix obtained from the totally symmetric ground state $\Psi^{GS}\left(\vec{r}_{1},\vec{r}_{2}\right)$, Eq. (9).

[TABLE]

It is possible to reduce the computation effort by noting that the iterated kernel ($\rho$) of a symmetric kernel ($\Psi^{GS}$) has the same eigenfunctions as the kernel, while the iterated eigenvalues are the squared eigenvalues of the kernel stricomi_1957 . This means that instead of directly solving Eq. (14) one can solve the following equation

[TABLE]

with $\Lambda_{i}=\lambda_{i}^{2}$ (see Eq. (14)). Solving the eigenvalue problem Eq. (15) is equivalent to finding the Schmidt decomposition of the functions $\psi_{x}(x_{1},x_{2})$ and $\psi_{y}(y_{1},y_{2})$.

To this end we use the Mehler’s formula serdelyi_1957 ,

[TABLE]

to write the Schmidt decomposition of Eqs. (12) and (13),

[TABLE]

Then $\psi_{x}(x_{1},x_{2})$ can be written as

[TABLE]

where $\varphi_{l}$ are the harmonic oscillator states,

[TABLE]

Taking

[TABLE]

Eq. (18) can be re-written as

[TABLE]

It is worth to mention that $\left\langle\left.\phi_{l}^{+}\left(u\right)\right|\phi_{\tilde{l}}^{-}\left(u\right)\right\rangle=0$, but $\left\langle\left.\phi_{l}^{+}\left(u\right)\right|\phi_{\tilde{l}}^{+}\left(u\right)\right\rangle=\delta_{l,\tilde{l}}$ and $\left\langle\left.\phi_{l}^{-}\left(u\right)\right|\phi_{\tilde{l}}^{-}\left(u\right)\right\rangle=\delta_{l,\tilde{l}}$ only if the overlap $\left\langle\left.\varphi_{l}\left(u-\frac{x_{0}}{2}\right)\right|\varphi_{\tilde{l}}\left(u+\frac{x_{0}}{2}\right)\right\rangle=0$. This overlap decreases when $x_{0}$ increases and it vanish for $x_{0}\to\infty$. Therefore, Eq. (21) is the Schmidt decomposition of the wave function of Eq. (11) up from some interaction strength value large enough to guarantee that the minima $\pm x_{0}$ are sufficiently far from each other. Moreover, it is the Schmidt decomposition of the reduced density matrix if $\lambda_{l}$ is replaced by $\Lambda_{l}=\lambda_{l}^{2}$. From Eq. (21) it is clear that each occupancy $\Lambda_{l}$ is doubly degenerate with natural orbitals $\phi_{l}^{+}\left(u\right)$ and $\phi_{l}^{-}\left(u\right)$, then, the normalization is given by $\sum\limits_{l}^{\infty}\Lambda_{l}=1/2$.

Due to the separability of the wave function, the natural orbitals are the product of a natural orbital associated to $\psi_{x}(x_{1},x_{2})$ i.e. $\phi_{l}^{+}\left(u\right)$ and $\phi_{l}^{-}\left(u\right)$, and a natural orbital associated to $\psi_{y}(y_{1},y_{2})$ with the following form,

[TABLE]

After performing some algebra, we get the occupancies in the limit of large interaction strength parameter $g$ as

[TABLE]

where

[TABLE]

and

[TABLE]

II Generalization to dimension higher than two

Although all the calculations were heretofore carried out for two dimensions, the generalization to three and more dimensions is straightforward, we present here the generalization for dimension $D$. A system consisting of two particles in a $D$-dimensional anisotropic harmonic trap interacting via some potential which depends on the distance between particles, is given by the following Hamiltonian,

[TABLE]

with $\varepsilon_{D-1}>\varepsilon_{D-2}>...>\varepsilon_{1}>1$, and $x_{ij}$ denoting the $i$-th coordinate of the particle $j$. The minima satisfy

[TABLE]

It is worth to notice that in dimension $D$ the minima of the potential in the isotropic case lie on the $D$-dimensional shell of radius $x_{0}$ (see Eq. (27)).

The occupancies are

[TABLE]

where

[TABLE]

and

[TABLE]

and the replacement $\varepsilon\mapsto\varepsilon_{i}$.

The linear entropy is

[TABLE]

with $\zeta(\omega_{x})$ and $\xi(\varepsilon)$ as in Eqs. (24) and (25) respectively.

Finally, the Rényi entropies are given by

[TABLE]

where

[TABLE]

and

[TABLE]

The von Neumann entropy can be obtained as a limiting case of the Rényi entropies when $\alpha\to 1$, with

[TABLE]

[TABLE]

III Inverse logarithmic interaction