Quantum breathing mode of trapped systems in one and two dimensions

We briefly recall the physical setting which has already been described in our previous works [22–24]. We aim at the quantum mechanical description of N identical particles in a d-dimensional space. The time-evolution of the corresponding wave function is governed by the time-dependent Schrödinger equation

$$\begin{equation} \mathrm{i} \frac{\mathrm{d}}{\mathrm{d}t} | \Psi(t) \rangle = \hat H(t) | \Psi(t) \rangle, \end{equation} \tag{ 1 }$$

where we set ℏ = 1. We assume that the Hamiltonian,

$$\begin{equation} \hat{H}(t) = \hat{H}_0 + \hat{H}_1(t), \end{equation} \tag{ 2 }$$

consists of the stationary part

$$\begin{equation} \hat{H}_0 = \hat{T} + \hat{V} + \hat{W} \end{equation} \tag{ 3 }$$

and an additional perturbation term, $\hat {H}_1(t)$. The explicit form of $\hat {H}_0$ in the spatial coordinates r ≡ (r₁,...,r_N) is given by the kinetic energy

$$\begin{equation} T(\mathbf{r}) = \sum_{i=1}^N - \frac{1}{2m} \frac{\partial^2}{\partial \mathbf{r}_i^2} \end{equation} \tag{ 4 }$$

the trap energy

$$\begin{equation} V(\mathbf{r}) = \sum_{i=1}^N \frac{1}{2}m\Omega^2 \mathbf{r}_i^2, \end{equation} \tag{ 5 }$$

and the interaction energy

$$\begin{equation} W(\mathbf{r}) = \sum_{i<j} \frac{K_\alpha}{|\mathbf{r}_i-\mathbf{r}_j|^\alpha}. \end{equation} \tag{ 6 }$$

If it is not further specified, we assume that the interaction is characterized by a repulsive power-law potential, w(r)∝1/r^α, with the proportionality constant K_α. Below we will concentrate on the important cases of Coulomb interaction (α = 1) and dipole interaction (α = 3). After rescaling all lengths, times and energies in terms of $l_0= \left (1 / (m\Omega ) \right )^{1/2}$, Ω⁻¹ and Ω, respectively, $\hat {H}_0$ takes the form

$$\begin{equation} {H}_0(\mathbf{r}) = \frac{1}{2}\sum_{i=1}^N \left \{-\frac{\partial^2}{\partial \mathbf{r}_i^2} + \mathbf{r}_i^2 \right \} +\sum_{i<j} \frac{\lambda}{|\mathbf{r}_i-\mathbf{r}_j|^\alpha}. \end{equation} \tag{ 7 }$$

The dimensionless coupling parameter λ > 0 determines the relative strength of the interaction energy [24]. It is given by

$$\begin{equation} \lambda = \frac{1}{\hbar\Omega} \frac{K_\alpha}{l_0^\alpha} \end{equation} \tag{ 8 }$$

with, e. g. K₁ = q²/(4π), for the Coulomb interaction of particles with charge q and the permittivity , and K₃ = C_dd/(4π) for the interaction of polarized magnetic or electric dipoles, with the corresponding proportionality constant C_dd [31]. While, in general, the interaction of two dipoles depends on their orientation, here we consider only the case of parallel dipoles [32, 33]. In experiments, parallel alignment perpendicular to the plane of the dipoles is typically realized with external fields.

To formally include a weak mode excitation by an arbitrary operator $\hat {Q}$, acting only at the time t = 0, we specify

$$\begin{equation} \hat{H}_1(t) = \eta \delta(t) \hat{Q}, \end{equation} \tag{ 9 }$$

where η is a small real parameter. In the following, we assume that the eigenstates of $\hat {H}_0$ are given by $\left \{ |0\rangle , |1\rangle , \ldots \right \}$ with the corresponding eigenvalues $\left \{ E_0, E_1, \ldots \right \}$. Furthermore, the system is supposed to be initially in the ground state |0〉. With the help of first-order time-dependent perturbation theory, it can be shown that the time-dependent expectation value of an operator $\hat {A}$ without explicit time-dependence takes the form [24]

$$\begin{equation} \langle \hat A \rangle(t) = \sum_{ij} c^*_i c_j \exp \left\{\mathrm{i} (E_i - E_j ) t \right\} \langle i| \hat A| j\rangle \end{equation} \tag{ 10 }$$

with the constant coefficients

$$\begin{equation} c_k = \delta_{k,0} - \mathrm{i} \eta \langle k| \hat Q | 0\rangle. \end{equation} \tag{ 11 }$$

The breathing mode is excited by the monopole operator $\hat {Q} = \hat {\mathbf {r}}^2 = \sum _i \hat {\mathbf {r}}_i^2$. For the typical observable $\hat {A}=\hat {\mathbf {r}}^2$, one can reduce equation (10) to the expression

$$\begin{equation} \langle \hat{\mathbf{r}}^2\rangle(t) = \langle 0 | \hat{\mathbf{r}}^2|0\rangle - {2\eta} \sum_i | \langle 0 | \hat{\mathbf{r}}^2|i\rangle|^2 \sin(\omega_{i,0} t),\end{equation} \tag{ 12 }$$

where ω_i,0 = E_i − E₀ are the mode frequencies. Equation (12) reveals that the quantum breathing mode is characterized by a superposition of sinusoids with different frequencies. For a full characterization of the breathing mode, one has to calculate the eigenvalues of $\hat {H}_0$ and the matrix elements $\langle 0 | \hat {\mathbf {r}}^2|i\rangle$.

Although the breathing motion comprises a variety of possible frequencies, it was shown in recent works [22–24] that the breathing mode is dominated by just two frequencies. One of these has the universal value 2 Ω. This value can be explained by a formal decoupling of the wave function into a center-of-mass (CM) and a relative part

$$\begin{equation} |\Psi (t)\rangle = |\Psi_{\mathrm{cm}}(t) \rangle \otimes |\Psi_{\mathrm{rel}}(t) \rangle. \end{equation} \tag{ 13 }$$

Such a decoupling is induced by the splitting of the Hamiltonian [34]

$$\begin{equation} \hat{H}(t) = \hat{H}_{\mathrm{cm}}(t) + \hat{H}_{\mathrm{rel}}(t), \end{equation} \tag{ 14 }$$

where the contributions read

$$\begin{equation} \hat{H}_{\mathrm{cm}}(t) = \frac{N}{2} \hat{\mathbf{P}}^2 + \frac{N}{2} \hat{\mathbf{R}}^2 + \eta \delta(t) N \hat{\mathbf{R}}^2, \end{equation} \tag{ 15 }$$

$$\begin{equation} \hat{H}_{\mathrm{rel}}(t) = \sum_{i<j} \left\{ \frac{1}{2N} \hat{\mathbf{p}}_{ij}^2 + \frac{1}{2N} \hat{\mathbf{r}}_{ij}^2 + \hat w(|\hat{\mathbf{r}}_{ij}|) + \eta \delta(t) \frac{1}{N}\hat{\mathbf{r}}_{ij}^2 \right\}. \end{equation} \tag{ 16 }$$

Here, we have introduced the CM and relative contributions, according to

$$\begin{equation} \hat{\mathbf{O}} = \frac{1}{N} \sum_{i=1}^{N} \hat{\mathbf{o}}_i, \end{equation} \tag{ 17 }$$

$$\begin{equation} \hspace{-4pt}\hat{\mathbf{o}}_{ij } = \hat{\mathbf{o}}_i - \hat{\mathbf{o}}_j. \end{equation} \tag{ 18 }$$

Since equation (15) describes a non-interacting d-dimensional oscillator problem for the observables $\hat {T}_{\mathrm {cm}} = N \, \hat {\mathbf {P}}^2 / 2$ and $\hat {V}_{\mathrm {cm}} = N \,\hat {\mathbf {R}}^2 / 2$, one can conclude that the quantities $\langle \hat {V} \rangle (t)$ and $\langle \hat {T}\rangle (t)$ always contain oscillations with the frequency ω_cm = 2 Ω. The contributions from the relative system, however, are non-trivial. Depending on the coupling parameter λ, the dominating frequency obtains the values $\sqrt {3} \, \Omega \leqslant \omega _{\mathrm {rel}} \leqslant 2 \, \Omega$, for Coulomb interaction, and $2 \, \Omega \leqslant \omega _{\mathrm {rel}} \leqslant \sqrt {5} \, \Omega$, for dipole interaction [22, 23, 35]. This frequency corresponds to the first monopole excitation in the relative system.

The weighted moments are defined by

$$\begin{equation} m_k = \sum_{i=1}^\infty (E_i - E_0)^k \, |\langle 0 |\hat{Q}|i\rangle|^2 \end{equation} \tag{ 19 }$$

for any operator $\hat Q$ and any integer number $k \in \mathbb {Z}$. Containing the exact excitation energies, the moments can be used to define average excitation energies

$$\begin{equation} E_{k,l} = \left(\frac{m_k}{m_{k-l}} \right)^{1/l} \end{equation} \tag{ 20 }$$

for positive integer numbers l. In the literature [27], one often finds the quantities E_k,2 and E_k,1. The average excitation energies fulfill the relation

$$\begin{equation} \ldots \geqslant E_{k+2,1} \geqslant E_{k+2,2} \geqslant E_{k+1,1} \geqslant E_{k+1,2} \geqslant \cdots \end{equation} \tag{ 21 }$$

and, especially, [29]

$$\begin{equation} \lim_{k \to -\infty} E_{k,1} = \omega_{a,0}, \end{equation} \tag{ 22 }$$

where the index a corresponds to the lowest state excited by the operator $\hat {Q}$. Instead of directly evaluating the sum in equation (19), one can make use of the sum rules to simplify selected moments. In the following, we will concentrate on the calculation of the moments m₃, m₁ and m₋₁ for the monopole operator $\hat {Q} = \hat {\mathbf {r}}^2$. For m₁ and m₃, one finds the sum rules

$$\begin{equation} m_1 = \frac{1}{2} \langle 0 | [\hat{Q}, [\hat{H}_0, \hat{Q} ]] |0\rangle \end{equation} \tag{ 23 }$$

$$\begin{equation} \hphantom{m_1}= 2 \ \langle 0 | \hat{\mathbf{r}}^2 |0\rangle \end{equation} \tag{ 24 }$$

and

$$\begin{equation} m_3 = \frac{1}{2} \langle 0 | [ [\hat{Q},\hat{H}_0] , [\hat{H}_0, [\hat{H}_0, \hat{Q} ]]]|0 \rangle \end{equation} \tag{ 25 }$$

$$\begin{equation}\hphantom{m_1} = {8 \langle \hat{T} \rangle + 8 \langle \hat{V} \rangle + 2 \alpha^2 \langle\hat{W}\rangle}. \end{equation} \tag{ 26 }$$

The latter result can be obtained, reducing equation (25) to the evaluation of commutators with the types $[\hat {\mathbf {p}}_i, \hat {\mathbf {r}}_j]$ and $[\hat {\mathbf {p}}_i, 1 / | \hat {\mathbf {r}}_j - \hat {\mathbf {r}}_k|^\alpha ]$ for all occurring indices i, j, k. The moment m₋₁ can be calculated by

$$\begin{equation} m_{-1} = - \frac{\partial}{\partial \gamma}\langle \hat{\mathbf{r}}^2 \rangle_{\gamma=1}, \end{equation} \tag{ 27 }$$

where the expectation value refers to the ground state of the Hamiltonian

$$\begin{equation} \hat{H}_0(\gamma) := \hat{T} + \gamma \hat{V} + \hat{W}. \end{equation} \tag{ 28 }$$

This becomes clear if one writes the derivative in equation (27) as

$$\begin{equation} \frac{\partial}{\partial \gamma}\langle \hat{\mathbf{r}}^2 \rangle_{\gamma=1} = \lim_{\epsilon \to 0} \frac{\langle\hat{\mathbf{r}}^2 \rangle_{\gamma=1+\epsilon} - \langle\hat{\mathbf{r}}^2 \rangle_{\gamma=1}}{\epsilon} \end{equation} \tag{ 29 }$$

and evaluates $\langle \hat {\mathbf {r}}^2 \rangle _{\gamma =1+\epsilon }$ with stationary perturbation theory [29]. Henceforth, we will use the notation sr(k,k − l): = E_k,l. With the moments m₃ and m₁, one obtains the convenient sum rule formula

$$\begin{eqnarray} \mathrm{sr}(3,1)&= \left\{\left( 2+\alpha \right) + \left( 2-\alpha \right) \frac{\langle \hat{T} \rangle}{\langle \hat{V} \rangle} \right\}^{1/2}\nonumber\\ &= \left\{\left( 2+\alpha \right) + \left( 2-\alpha \right) \left( 1- \frac{\alpha\langle \hat{W} \rangle }{2\langle \hat{V} \rangle} \right) \right\}^{1/2}. \end{eqnarray} \tag{ 30 }$$

The two different representations of sr(3,1) in equation (30) are due to the virial theorem

$$\begin{equation} 2 \langle \hat{T}\rangle - 2 \langle \hat{V}\rangle + \alpha \langle \hat{W}\rangle = 0. \end{equation} \tag{ 31 }$$

For the special case of Coulomb interaction (α = 1), equation (30) has already been used by Sinha [39]. A similar result for contact interaction was derived by Stringari [38]. Furthermore, we mention that equation (30) is a quantum generalization of the formula by Olivetti et al [21] for classical systems. A revealing property of equation (30) is the fact that it allows one to interpret the behavior of the breathing frequency with the ratios of the contributions to the total energy. Especially the classical limit, with $\langle \hat {T}\rangle /\langle \hat {V}\rangle =0$, and the ideal quantum limit, with $\langle \hat {W}\rangle /\langle \hat {V}\rangle =0$, are included.

Another sum rule formula we will use in this paper is given by

$$\begin{equation} \mathrm{sr}(1,-1)= \left\{ -2 \frac{\langle \mathbf{r}^2 \rangle}{ \left[ \partial \langle \mathbf{r}^2 \rangle/\partial \gamma \right]_{\gamma=1}} \right\}^{1/2}. \end{equation} \tag{ 32 }$$

It has been presented as a rigorous upper bound of the breathing frequency by Menotti and Stringari [19].

So far, we expressed the moments with the eigenstates and the eigenvalues of the full Hamiltonian $\hat {H} = \hat {H}_{\mathrm {rel}}+\hat {H}_{\mathrm {cm}}$. This may be disadvantageous for the estimation of ω_rel if the weights of the CM terms are comparable to those of the relative terms. With the notations

$$\begin{equation} \hat{\mathbf{R}}_{\mathrm{cm}}^2 := N \hat{\mathbf{R}}^2 \end{equation} \tag{ 33 }$$

and

$$\begin{equation} \hat{\mathbf{r}}_{\mathrm{rel}}^2 := \frac{1}{N} \sum_{q<r} \hat{\mathbf{r}}^2_{qr}, \end{equation} \tag{ 34 }$$

this can be understood more formally, expressing the moments as

$$\begin{eqnarray}&&\fl m_l =\sum_{i=1}^\infty ( E_{\mathrm{rel},i} - E_{\mathrm{rel},0} )^l |\langle 0_{\mathrm{rel}} | \hat{\mathbf{r}}_{\mathrm{rel}}^2 | i_{\mathrm{rel}}\rangle |^2 + \sum_{k=1}^\infty ( E_{\mathrm{cm},k} - E_{\mathrm{cm},0} )^l |\langle 0_{\mathrm{cm}} | \hat{\mathbf{R}}_{\mathrm{cm}}^2 | k_{\mathrm{cm}}\rangle |^2. \end{eqnarray} \tag{ 35 }$$

Using the formulae sr(k,j), the contributions from the second sum in equation (35) may cause considerable inaccuracies of the estimator for ω_rel = E_rel,1 − E_rel,0. Furthermore, it cannot be guaranteed that the sum rule formulae are upper bounds of the frequency ω_rel if there exist non-vanishing contributions E_cm,k − E_cm,0 < ω_rel. For example, this is the case for dipole interaction, where ω_cm ⩽ ω_rel is valid for all coupling parameters. A simple solution to this problem is to eliminate the second sum in equation (35). As all terms are known analytically, this can simply be accomplished by introducing the corrected moments

$$\begin{equation} m^*_k := m_k - 2^{k-1}d. \end{equation} \tag{ 36 }$$

Using these moments, we introduce the following improved sum rule formulae:

$$\begin{equation} \mathrm{sr}^*(k,k-l) := \left( \frac{m^*_k}{m^*_{k-l}} \right)^{1/l}. \end{equation} \tag{ 37 }$$

These formulae are not only expected to yield more accurate results than the conventional formulae, they also restore the character of the approximation as an upper bound for all α. The special case sr*(3,1) is equivalent to the equilibrium formula presented in [26].

To demonstrate the difference between the formulae sr and sr*, we exactly determine the spectrum of the two-particle system. Expressing the relative vector r_1,2 in the hyperspherical coordinates $\left (\rho , \phi _1, \ldots , \phi _{d-1} \right )$ [40], one can reduce the relative problem to the equation

$$\begin{equation} \left\{-\frac{1}{2}\frac{{d^2}}{{d}\rho^2} + \frac{\rho^2}{2} + \frac{\left(l+(d-2)/2 \right)^2 - 1/4}{2\rho^2} + \frac{\lambda}{2^{\alpha/2}\rho^\alpha} - E_{\mathrm{rel}}\right\}u_l(\rho) = 0. \end{equation} \tag{ 38 }$$

We obtain the spectrum by solving the corresponding eigenvalue problem in the matrix representation that arises from the expansion of u_l in terms of finite-element discrete variable representation (FEDVR) basis functions (see section 3.1). As the spectrum of the Hamiltonian $\hat {H}_{\mathrm {cm}}$ is known analytically, we can directly reconstruct the moments. In figure 1, we compare the exact breathing frequencies with the sum rule estimators for Coulomb and dipole interaction. In both cases, the improvement of the sum rules leads to a higher accuracy. For dipole interaction, it is shown that, in fact, the frequencies from the estimators sr(3,1) and sr(1,−1) are below the exact values. At the same time, it can be seen that this problem is solved by the improved formulae.

Zoom In Zoom Out Reset image size

To conclude, we remark that the subtractive correction of the moments in equation (36) does not depend on the particle number. Therefore, it can be expected that the improvement of the sum rules is only important for small systems, where the eigenvalues of $\hat {H}_{\mathrm {rel}}$ and $\hat {H}_{\mathrm {cm}}$ are of the same order. Nevertheless, we stress the difference between the results, because the influence of the CM subsystem was not mentioned in some other works, e.g. [37] for the case of dipole interaction.

The HF approximation [41, 42] reduces the N-body problem to an effective one-body problem, where the interactions are taken into account on the mean-field level. The central idea of the HF method is to assume that the solution of the Schrödinger equation is a single Slater determinant

$$\begin{equation} |\Psi\rangle = |\phi_1 \ldots \phi_N \rangle. \end{equation} \tag{ 39 }$$

Requiring that the set of single-particle spin orbitals |ϕ₁〉,...|ϕ_N〉 minimizes the total energy $E=\langle \Psi |\hat {H}_0|\Psi \rangle$, one can derive the HF equations

$$\begin{equation} \hat{F} |\phi_k\rangle = \epsilon | \phi_k\rangle, \end{equation} \tag{ 40 }$$

where

$$\begin{equation} \hat{F} = \hat{h} + \sum_{i=1}^N \hat{J}_i - \hat{K}_i \end{equation} \tag{ 41 }$$

is the Fock operator. For the trapped system, $\hat {h}$ takes the form

$$\begin{equation} h(\mathbf{r}) = -\frac{1}{2}\frac{\partial^2}{\partial \mathbf{r}^2} + \frac{1}{2}\mathbf{r}^2. \end{equation} \tag{ 42 }$$

The interaction is incorporated by the operators $\hat {J}_i$ and $\hat {K}_i$, which are defined by their actions

$$\begin{equation} \hat{J}_i |\phi_k\rangle = \lambda \langle \phi_i |\hat{w}|\phi_i\rangle|\phi_k\rangle, \end{equation} \tag{ 43 }$$

$$\begin{equation} \hat{K}_i |\phi_k\rangle = \lambda \langle \phi_i | \hat{w} |\phi_k\rangle | \phi_i \rangle. \end{equation} \tag{ 44 }$$

Expanding the spin orbitals in terms of a single-particle basis, we transfer the HF equations to a matrix equation. For 1D systems, the utilized basis is an FEDVR  [43, 44]. To avoid divergences, the Coulomb potential is regularized, according to the standard formula [22, 23]

$$\begin{equation} w(|r-r'|)=\frac{1}{\left[ (r-r')^2+ \kappa^2 \right]^{1/2}} \end{equation} \tag{ 45 }$$

with the screening parameter κ. In our recent work [24], we checked that converged results are obtained with κ = 0.1. In 2D systems, by contrast, such a screening is not necessary, as the basis is formed by the eigenfunctions of the harmonic oscillator in spherical coordinates. Our implementation for the calculation of the matrix elements is based on the code OpenFCI by Kvaal [45].

For the extension of the HF results, we use the well-known TF approximation [46–49]. It is expected to yield the correct trend for large particle numbers, because the oscillations of the single-particle density—which are not shown by the density in TF approximation—become negligible [24].

In the 1D case, the TF equation for spin-polarized particles reads

$$\begin{equation} \frac{\pi^2}{2} n(r)^2 + \frac{1}{2}r^2 + \lambda \int \frac{n(r')}{\left[ (r-r')^2+ \kappa^2 \right]^{1/2}}\, \mathrm{d}r' = \mu, \end{equation} \tag{ 46 }$$

where μ is the chemical potential. Fixing the particle number N and requiring the normalization

$$\begin{equation} \int n(r)\,\ \mathrm{d}r = N, \end{equation} \tag{ 47 }$$

one has to find the density n(r) and the corresponding chemical potential which solve equation (46). We obtained our results on a grid, according to the following scheme: for each considered particle number, we start our calculation with a very small coupling parameter λ. At this, we choose the initial trial density

$$\begin{equation} n(r) = \frac{1}{\pi} \sqrt{2N-r^2}, \end{equation} \tag{ 48 }$$

which is exact for a non-interacting system [50]. Updating the chemical potential and the density in a self-consistent procedure [51], we finally obtain the density that solves equation (46). After that, we slightly increase λ and start a new calculation, where the initial density is the final density from the previous calculation. This procedure is repeated until some final value of λ is reached.

For 2D systems, we follow a simpler approach by Sinha [39]. Making the ansatz

$$\begin{equation} n(\mathbf{r}) = n(r) = \frac{1}{2\pi \gamma} (r_0^2 - r^2) \end{equation} \tag{ 49 }$$

for the density, we determine the variational parameter γ that minimizes the total energy

$$\begin{equation} E = \frac{1}{3}N^{3/2} \frac{1}{\gamma^{1/2}} +\frac{1}{3}N^{3/2} \gamma^{1/2} + \lambda \frac{512}{315} \frac{\sqrt{2}}{\pi \gamma^{1/4}} N^{7/4}. \end{equation} \tag{ 50 }$$

In this equation, the single terms correspond to the kinetic energy, the trap energy and the interaction energy in this order.

The method for the 2D case has the advantage that the parabolic ansatz for the density allows us to express the interaction energy analytically. In the 1D case, however, the calculations are more complex, because the integral

$$\begin{equation} \lambda \int \frac{n(r')}{\left[ (r-r')^2+ \kappa^2 \right]^{1/2}}\, \mathrm{d}r' \end{equation} \tag{ 51 }$$

has to be calculated for each grid point r in each iteration step.

4.1.1. One-dimensional systems.

We start with an analysis of small 1D systems. In figure 2, we compare the breathing frequencies from different methods for the coupling parameters λ = 0.3 and 1. On the one hand, we show the time-dependent HF and CI results from our previous work [24]. On the other hand, we show the results obtained with static HF calculations in combination with the conventional and the improved sum rule formulae, respectively. One can draw the following conclusions: comparing with the exact CI results, the results from the improved sum rule formulae are more accurate than the TDHF results. The improved accuracy of the sum rules can be explained by the fact that the static HF calculations could be performed with a larger basis than the time-dependent HF calculations. Furthermore, the error induced by the sum rules appears to be less important than the error induced by the HF approximation. Finally, one notices that the improvement of the sum rule formulae does not only lead to a higher accuracy for small particles, it also enables one to find the frequency minimum—which was already discovered in our previous work [24]—at the correct position N = 6, instead of N = 7. However, the figure already reveals that the difference between both formulae tends to vanish for large particle numbers.

Zoom In Zoom Out Reset image size

4.1.2. Two-dimensional systems.

The step from 1D to 2D systems is usually numerically demanding, because—roughly estimated—the number of single-particle basis functions has to be squared. This is especially challenging for time-dependent calculations. However, using the sum rules and the HF method, we are able to investigate the finite-size effects of the breathing mode in 2D Coulomb systems for the first time. In figure 3, we show the frequencies for small (N ⩽ 22) systems with an intermediate coupling parameter λ = 0.5. Compared to the 1D systems, we observe that the non-monotonic behavior of the frequency has increased. More precisely, the frequency has local minima for 2, 6, 12 and 20 particles. These 'magic' particle numbers are well-known from various experimental and theoretical studies of quantum dots [8, 52, 53]. The corresponding configurations are very stable, because they are characterized by closed energy shells. Typically, experimental evidence for the occurrence of the magic configurations is given by the measurement of N-dependent addition energies. Apparently, the breathing mode provides an alternative tool for the diagnostics of these properties. Furthermore, it has been shown recently that important system properties, such as the kinetic and potential energy, can be obtained from the results of the breathing frequency [26].

Zoom In Zoom Out Reset image size

Using the HF approximation, we observe the same non-monotonic behavior for all coupling parameters λ ⩽ 1. Another important result is that—especially for two and three particles—the improvement of the sum rules is crucial to reproduce the correct behavior.

4.2.1. One-dimensional systems.

In our recent work [24], we found several hints that an increase of the particle number in 1D systems leads to a steady increase of quantum effects. This means that the breathing frequency slowly transitions to the quantum limit ω_rel = 2 Ω. However, the calculations of the frequencies were restricted to 20 particles with coupling parameters λ ⩽ 1. In this work, we extend the results, reporting the frequencies for up to 10 000 particles with λ ≲ 2.

To start the analysis for large systems, we compare the results from 1D HF calculations and the corresponding TF calculations for several hundred particles in the top plot of figure 4. As one can see, the difference between both approximations tends to vanish for large N. Hence, both approximations confirm the trend ω_rel → 2 Ω in the limit N → ∞. For an overview, the contour plot in figure 5 summarizes the behavior of the breathing frequency in the (λ,N)-plane. For all coupling parameters, the contours have a small step around N = 100, where the results from HF and TF calculations have been joined. Nevertheless, the plot is suitable to trace the finite size effects, with the frequency minima for N = 6, as well as monotonic behavior of large systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the frequencies approach the ideal limit only very slowly, we conclude with a remark that one can give a numerical proof for the large-N behavior in the TF approximation. For that purpose, one uses the ideal density to calculate the kinetic energy, the trap energy and the interaction energy (with arbitrary fixed λ) as a function of the particle number. Such a calculation reveals that, for each coupling parameter λ, there exists a large particle number for which the interaction energy becomes negligible compared to the other two contributions. Consequently, the ideal density becomes a correct solution of equation (46) in the limit N → ∞. In practice, that implies that an increase of N leads to a faster convergence of the density if the initial guess in the self-consistent procedure is the ideal density.

4.2.2. Two-dimensional systems.

In 2D systems, the N-dependent behavior of the breathing frequency is contrary to that of 1D systems. In the bottom plot of figure 4, one can see that the frequency always reaches the value $\sqrt {3}\,\Omega$, in the limit N → ∞. This behavior is confirmed for a broad range of coupling parameters in figure 6. Hence, 2D systems always approach the classical limit if either the particle number or the coupling parameter is increased.

Zoom In Zoom Out Reset image size

Compared to the results for the 1D systems, the ranges of numerically accessible λ and N are much broader for 2D systems. This can be explained by the fact that the parabolic ansatz for the density profile in the TF approximation, equation (49), is a drastic simplification of the problem. Hence, we can treat all coupling parameters with the same computational effort of a simple minimization procedure. Despite the restrictiveness of the parabolic ansatz, we can show that the results are fairly accurate. On the one hand, the good agreement of the TF and the HF results shown in figure 4 justifies the parabolic ansatz, at least for weakly coupled systems. On the other hand, to check that the ansatz can also be used to reproduce the correct energies of strongly coupled systems, we show the total energies for a broad range of coupling parameters and different particle numbers in figure 7. It turns out that—for sufficiently large particle numbers—the TF results converge to the results from classical molecular dynamics (MD) simulations in the limit of large λ.

Zoom In Zoom Out Reset image size

In the following, we provide some supportive explanation for the different behaviors of 1D and 2D systems. Firstly, we trace the characteristics with the help of a localization parameter that measures the extension of the system. Secondly, we introduce a simple estimator for the distinction between quantum and classical systems.

4.3.1. Localization parameter.

Being inspired by the degeneracy parameter of a macroscopic homogeneous electron gas, one can define a localization parameter for the trap. This parameter is meant to express the non-ideality of the system with geometric quantities. An estimator for the mean extension of the system is given by

$$\begin{equation} \sigma = \{ 2 \langle \hat{V} \rangle \}^{1/2}. \end{equation} \tag{ 52 }$$

For the non-interacting system, it has the value

$$\begin{equation} \sigma_{\mathrm{ideal}} =\frac{1}{\sqrt{2}} N, \quad \mathrm{(1D)},\end{equation} \tag{ 53 }$$

$$\begin{equation} \sigma_{\mathrm{ideal}} = \sqrt{2/3} N^{3/4}, \quad \mathrm{(2D)}. \end{equation} \tag{ 54 }$$

With this, one can define the localization parameter

$$\begin{equation} \chi = \frac{\sigma_{\mathrm{ideal}}}{\sigma}, \end{equation} \tag{ 55 }$$

which measures how much the extension of the system deviates from that of an ideal non-interacting quantum system. Starting with the value χ = 1 in the ideal system, the localization parameter decreases to zero in the course of the transition to the strongly coupled regime. In figure 8, we compare the localization parameter with the breathing frequency for 1D and 2D systems in the (λ,N)-plane. The data were produced with the TF approximation and the frequencies were estimated with the sum rule formula sr*(3,1). We concentrate on particle numbers N ⩾ 100 to make sure that finite-size effects are reduced. Furthermore, since the 1D TF method is reliable for small coupling parameters, we restrict the illustration to the regime λ ≲ 1.

Zoom In Zoom Out Reset image size

The figure shows that the localization parameter is well-suited to track the qualitative behavior of the breathing frequency. Equal values of χ correspond to equal values of the breathing frequency. To illustrate this, the plot shows exemplary dotted lines which indicate equal values of each quantity. In the logarithmic plot, these are straight lines. Remarkably, the lines for 1D systems have a positive slope, while the lines for 2D systems have a negative slope. The corresponding best fits in the ranges 100 ⩽ N ⩽ 10 000 and 0.1 ⩽ λ ⩽ 1 read

$$\begin{equation} N = (1.3 \times 10^4) \lambda^{2.45}, \quad \mathrm{(1D)}, \end{equation} \tag{ 56 }$$

$$\begin{equation} N = 100 \lambda^{-4}, \quad \mathrm{(2D)}. \end{equation} \tag{ 57 }$$

4.3.2. Estimator for intermediate couplings.

Finally, we provide another rough estimator for a supportive understanding of the observed behavior. For the 2D system, the straight lines in figure 7 demonstrate that the total energies can be approximated by those from non-interacting Fermi gases for weak couplings,

$$\begin{equation} E_{\mathrm{ideal}}^{\mathrm{2D}}= \frac{2}{3} N^{3/2} \end{equation} \tag{ 58 }$$

and those from purely classical systems for strong couplings,

$$\begin{equation} E_{\mathrm{classical}}^{\mathrm{2D}}= K \lambda^{2/3} N^{5/3}. \end{equation} \tag{ 59 }$$

The last equation is derived by setting the kinetic energy to zero and minimizing the remaining terms in equation (50). The constant K has the value

$$\begin{equation} K= \left( \frac{256\sqrt{2}}{315\pi} \right)^{2/3} + \frac{512\sqrt{2}}{315\pi} \left( \frac{256\sqrt{2}}{105\pi} \right)^{-1/3}. \end{equation} \tag{ 60 }$$

As has been shown by Astrakharchik and Girardeau [50], one can also derive the corresponding terms for 1D systems, where one obtains

$$\begin{equation} E_{\mathrm{ideal}}^{\mathrm{1D}}= \frac{1}{2} N^2 \end{equation} \tag{ 61 }$$

for non-interacting systems and

$$\begin{equation} E_{\mathrm{classical}}^{\mathrm{1D}}= \frac{3}{10} (3 \lambda N \ln N)^{2/3}N \end{equation} \tag{ 62 }$$

for classical systems. In both 1D and 2D, we use the quantities E_ideal and E_classical to define a regime of intermediate coupling. Taking account of the λ-dependence of the total energy shown in figure 7, we mark this region by the coupling parameter $\tilde {\lambda }$ for which the classical estimator E_classical is equal to the ideal estimator E_ideal. One obtains

$$\begin{equation} \tilde{\lambda} = \left( \frac{2}{3K} \right)^{3/2} N^{-1/4} \end{equation} \tag{ 63 }$$

for 2D systems, which is in agreement with the functional form of equation (57). For the 1D systems, one obtains

$$\begin{equation} \tilde{\lambda} =\frac{1}{3} \left( \frac{5}{3} \right)^{3/2} \frac{N^{1/2}}{\ln N}. \end{equation} \tag{ 64 }$$

In this equation, the relation between N and $\tilde \lambda$ slightly differs from that in equation (56). However, we expect that both definitions come to agreement if one covers a larger area of the (λ,N)-plane shown in figure 8.

The coupling parameter $\tilde {\lambda }$ roughly divides the systems into the ones with dominating quantum-like behavior ($\lambda < \tilde {\lambda }$) and the other ones with dominating classical behavior ($\lambda > \tilde {\lambda }$). Apparently, the trend for large systems is as follows: an increase of N leads to a growing classical regime in 2D and a growing quantum regime in 1D. In figure 9, equations (63) and (64) are plotted. Remarkably, for 1D systems, $\tilde {\lambda }$ has a minimum for N = 7 and thus even reproduces the non-monotonic behavior.

Zoom In Zoom Out Reset image size