Mesoscopic quantum superposition states of weakly-coupled matter-wave solitons

The atomic JJ system can be created basically by means of appropriate condensate trapping. In various experiments with JJ the BEC is confined in a potential V = V_H + V_W, where V_H is familiar 3D harmonic trapping potential, V_W is responsible for condensate double-well confinement [6, 67]. The aim of such a trapping is to split the condensate into two parts which possess tunneling through the barrier created by V_W potential.

The scheme of possible realization of JJ-device with two weakly coupled condensates is evident from figure 1, where for clarity we depicted the probability distribution for 2D condensates. In particular, the JJ platform consists of two highly asymmetric (cigar-shaped) condensates elongated in X direction and trapped in a double-well potential (green line in figure 1) allowing some distance d between trap centers. In experiments with dipole traps d is few micrometers [67]. The BEC particles possess an interaction that we suppose attractive in this work. This assumption allows to examine Gaussian (BJJ-model) and bright soliton (SJJ-model) solutions for the condensate wave functions in the framework of the same physical set-up illustrated in figure 1.

Zoom In Zoom Out Reset image size

We represent the condensate wave function Ψ(r_⊥, x, t) relevant to the scheme in figure 1 in a superposition state as

$$\begin{equation}{\Psi}\left({r}_{\perp },x,t\right)={{\Phi}}_{1}\left({r}_{\perp }-d/2\right){{\Psi}}_{1}\left(x,t\right)+{{\Phi}}_{2}\left({r}_{\perp }+d/2\right){{\Psi}}_{2}\left(x,t\right),\end{equation} \tag{ 1 }$$

where Φ_1,2(r_⊥ + d/2) and Φ₂(r_⊥ − d/2) are spatial (time independent) wave functions of condensates in transverse plane on either side of the barrier.

In equation (1) Ψ₁ ≡ Ψ₁(x, t) and Ψ₂ ≡ Ψ₂(x, t) are two condensate wave functions which characterize dynamical properties of JJ system and obey coupled Gross–Pitaevskii equations (cf [44])

$$\begin{equation}\mathrm{i}\frac{\partial }{\partial t}{{\Psi}}_{1}=-\frac{1}{2}\frac{{\partial }^{2}}{\partial {x}^{2}}{{\Psi}}_{1}-u{\left\vert {{\Psi}}_{1}\right\vert }^{2}{{\Psi}}_{1}+\frac{1}{2}{\nu }^{2}{x}^{2}{{\Psi}}_{1}-\kappa {{\Psi}}_{2};\end{equation} \tag{ 2a }$$

$$\begin{equation}\mathrm{i}\frac{\partial }{\partial t}{{\Psi}}_{2}=-\frac{1}{2}\frac{{\partial }^{2}}{\partial {x}^{2}}{{\Psi}}_{2}-u{\left\vert {{\Psi}}_{2}\right\vert }^{2}{{\Psi}}_{2}+\frac{1}{2}{\nu }^{2}{x}^{2}{{\Psi}}_{2}-\kappa {{\Psi}}_{1},\end{equation} \tag{ 2b }$$

where u = 2π|a_sc|/a_⊥ characterizes Kerr-like (focusing) nonlinearity, a_sc < 0 is the s-wave scattering length for attractive particles, ${a}_{\perp }=\sqrt{\hslash /m{\omega }_{\perp }}$ is a characteristic trap scale, and m is particle mass.

In equation (2) $\kappa \equiv \left\vert \mathcal{K}\right\vert /{\omega }_{\perp }{ >}0$ denotes an absolute value of the tunneling rate ($\mathcal{K}{< }0$), normalized on radial frequency, ω_⊥; $\mathcal{K}$ depends on the condensate wave functions Φ_1,2(r_⊥ ± d/2) overlapping in YZ-plane; ν = ω_x/ω_⊥ is a trap asymmetry parameter, ω_x and ω_⊥ are harmonic trap frequencies in X-axis and radial direction, respectively.

In the framework of field theory it is possible to consider a Hamilton function corresponding to equation (2), which looks like (cf [42, 43])

$$\begin{equation}H=\sum _{j=1}^{2}\left(\frac{1}{2}{\left\vert \frac{\partial {{\Psi}}_{j}}{\partial x}\right\vert }^{2}-\frac{u}{2}{\left\vert {{\Psi}}_{j}\right\vert }^{4}+{\left\vert {{\Psi}}_{j}\right\vert }^{2}{V}_{\mathrm{H}}\left(x\right)\right)-\kappa \left({{\Psi}}_{1}^{{\ast}}{{\Psi}}_{2}+{{\Psi}}_{1}{{\Psi}}_{2}^{{\ast}}\right),\end{equation} \tag{ 3 }$$

where ${V}_{\mathrm{H}}\left(x\right)=\frac{1}{2}{\nu }^{2}{x}^{2}$ is trapping potential in X direction. In (2) and (3), Ψ₁ and Ψ₂ obey a normalization condition

$$\begin{equation}{\int }_{-\infty }^{\infty }\left(\vert {{\Psi}}_{1}{\vert }^{2}+\vert {{\Psi}}_{2}{\vert }^{2}\right)\mathrm{d}x=N,\end{equation} \tag{ 4 }$$

where N is the average total number of particles.

In equations (2)–(4) we also propose rescaled (dimension-less) spatial and time variables, which are: x, y, z → x/a_⊥, y/a_⊥, z/a_⊥ and t → ω_⊥t, cf [42–44]. In other words, all lengths and time variables in the work are given in a_⊥ and ${\omega }_{\perp }^{-1}$ units respectively.

Let us briefly discuss limitations for equations (1)–(3) and conditions for bright atomic soliton formation that we are interested in this work. In particular, if the number of condensate particles, N, is not too large, spatial wave functions Φ_1,2(r_⊥ ± d/2) may be considered as Gaussian and still unchanging within appropriate time scales. Typically, this assumption is valid for weakly-coupled condensates possessing a moderate number of particles, N. In other words, equation (1) represents some ansatz for the wave function that provides so-called two-mode approximation used in this work. In particular, as it is shown in [68], characteristic trap scale, a_⊥, should be larger than nonlinear strength, N|a_sc|, in two-mode approximation. The results obtained in the experiments with condensate BJJs manifested validity of the two-mode approach up to few thousand particles [6, 8, 67].

However, condensates with attractive interaction between the particles possess instability and collapse if the number of atoms exceeds some critical value, N_c; references [69, 70] exhibit recent experimental efforts in this field. The critical particle number is N_c = 0.67a_⊥/|a_sc|, [52]. In particular, for condensate solitons with ⁷Li atoms possessing the negative scattering length, the effective nonlinear parameter is uN_c ≈ 4.2. In early experiment [54] authors demonstrated that N is limited by the number of atoms 5.2 × 10³ in the soliton providing N_c|a_sc| ≃ 1.105 μm. At the same time, it is shown that there exists some parameters window, where 1D bright soliton avoids collapsing.

Another way to obtain bright atomic solitons is based on formation of band gap solitons occurring due to the condensate confinement in a periodical optical lattice. Although the condensate possesses the positive scattering length, the negative effective mass m_eff < 0 of the particles that forms at the edge of Brillouin zone provides bright soliton occurrence in this case [71]. In [56] authors demonstrated band gap bright soliton formation with N ≃ 250 rubidium atoms interacting repulsively. Thus, the number of particles, N, in the soliton can be estimated from

$$\begin{equation}N=\frac{{a}_{\perp }^{2}m}{1.5\vert {m}_{\text{eff}}\vert {x}_{0}{a}_{\mathrm{s}\mathrm{c}}},\end{equation} \tag{ 5 }$$

where x₀ is the soliton width [56]. From (5) follows that the number of particles, N, may be reduced if we can enhance the value of scattering length a_sc by means of Feshbach resonance method, cf [72], or with tailoring the effective atom mass, m/|m_eff|, in the lattice, cf [73]. The ratio m/|m_eff| was equal to 10 in [56]. To be more specific, below we restrict ourselves by the matter-wave bright solitons with the negative scattering length.

Let us assume that the wave function shape in X-dimension is Gaussian and looks like

$$\begin{equation}{{\Psi}}_{j}=\frac{{\nu }^{1/4}\sqrt{{N}_{j}}}{{\pi }^{1/4}}{\text{e}}^{-\nu \enspace {x}^{2}/2}\enspace {\text{e}}^{\text{i}{\theta }_{j}},\end{equation} \tag{ 6 }$$

where N_j and θ_j are the average particle number and phase of jth condensate, respectively. The approach discussed here is valid at zero temperature in thermodynamic limit when the number of particles and occupied volume are extremely (infinitely) large but their ratio is finite. The losses and non-equilibrium phenomena are neglected.

Substituting (6) into (3) and omitting constant energy terms proportional to N for effective classical Hamiltonian ${H}_{\text{BJJ}}={\int }_{-\infty }^{\infty }H\enspace \mathrm{d}x$ one can obtain

$$\begin{equation}{H}_{\text{BJJ}}=-\frac{u\sqrt{\nu }}{2\sqrt{2\pi }}\left({N}_{1}^{2}+{N}_{2}^{2}\right)-2\kappa \sqrt{{N}_{1}{N}_{2}}\enspace \mathrm{cos}\left[\theta \right].\end{equation} \tag{ 7 }$$

It is instructive to represent the Hamiltonian (7) in terms of new variables z = (N₂ − N₁)/N and θ = θ₂ − θ₁. Omitting the constant energy term for (7) we obtain

$$\begin{equation}{H}_{\text{BJJ}}=\kappa N\left(-\frac{\lambda }{2}{z}^{2}-\sqrt{1-{z}^{2}}\enspace \mathrm{cos}\left[\theta \right]\right),\end{equation} \tag{ 8 }$$

where z and θ are normalized population imbalance and phase difference, respectively. In equation (8) we introduce the effective nonlinear parameter

$$\begin{equation}\lambda =\frac{\sqrt{\nu }\enspace uN}{2\sqrt{2\pi }\enspace \kappa }\equiv \frac{\sqrt{\nu }\enspace uN\enspace {\omega }_{\perp }}{2\sqrt{2\pi }\enspace \left\vert \mathcal{K}\right\vert }.\end{equation} \tag{ 9 }$$

The equation (8) with the key parameter λ represents a conventional JJ Hamiltonian and characterizes an atomic condensate BJJ-model intensively studied in the framework of various quantum and nonlinear features occurring as a result of tunneling and inter-particle interaction interplay [8, 14, 16, 17, 74].

It is worth noticing that the Hamiltonian (8) is invariant under the transformation κ → −κ and θ → θ + π. However, as it is argued in [75], physically correct tunneling rate $\mathcal{K}$ must be negative, cf [42–44].

Now let us consider the case when two condensates in X-dimension may be described by bright soliton solutions. We suppose that condensate trapping in X-dimension is weak enough that can be realized for highly elongated condensates obtained by means of asymmetric trapping potential with ν ≪ 1, see figure 1. Thereafter, we completely neglect the trapping term in equation (3), cf [44].

If coupling parameter κ is switched-off, the two condensates behave independently possessing non-moving bright soliton solutions in X-dimension, which are

$$\begin{equation}{{\Psi}}_{j}=\frac{{N}_{j}\sqrt{u}}{2}\enspace \mathrm{s}\mathrm{e}\mathrm{c}\mathrm{h}\left[\frac{u{N}_{j}}{2}x\right]{\text{e}}^{\text{i}{\theta }_{j}}.\end{equation} \tag{ 10 }$$

Substituting (10) into (3) for effective classical Hamiltonian ${H}_{\text{SJJ}}={\int }_{-\infty }^{\infty }H\enspace \mathrm{d}x$ one can obtain

$$\begin{equation}{H}_{\text{SJJ}}=-\frac{{u}^{2}}{24}\left({N}_{1}^{3}+{N}_{2}^{3}\right)-\frac{4\kappa {N}_{1}{N}_{2}}{N}I\left(z\right)\mathrm{cos}\left[\theta \right],\end{equation} \tag{ 11 }$$

where

$$\begin{equation}I\left(z\right)={\int }_{0}^{\infty }\frac{\mathrm{d}x}{{\mathrm{cosh}}^{2}\left[x\right]+{\mathrm{sinh}}^{2}\left[zx\right]}\approx 1-0.21{z}^{2}.\end{equation} \tag{ 12 }$$

The Hamiltonian (11) in terms of z and θ variables reads as

$$\begin{equation}{H}_{\text{SJJ}}={\kappa }_{\text{eff}}N\left(-\frac{{{\Lambda}}_{\text{eff}}}{2}{z}^{2}-\sqrt{1-{z}^{2}}\enspace \mathrm{cos}\left[\theta \right]\right),\end{equation} \tag{ 13 }$$

and looks similar to equation (8) obtained for the BJJ-model. In (13) we introduced the effective nonlinear parameter

$$\begin{equation}{{\Lambda}}_{\text{eff}}=\frac{{u}^{2}{N}^{2}}{16{\kappa }_{\text{eff}}}\equiv \frac{{\Lambda}}{\left(1-0.21{z}^{2}\right)\sqrt{1-{z}^{2}}},\end{equation} \tag{ 14a }$$

where

$$\begin{equation}{\Lambda}\equiv \frac{{u}^{2}{N}^{2}}{16\kappa }=\frac{{u}^{2}{N}^{2}\enspace {\omega }_{\perp }}{16\left\vert \mathcal{K}\right\vert }\end{equation} \tag{ 14b }$$

is the vital parameter that determines various regimes of soliton dynamics.

Equation (13) derived for the SJJ-model can be understood in the framework of the BJJ-model Hamiltonian (8) with tunneling rate, ${\kappa }_{\text{eff}}\equiv \kappa \left(1-0.21{z}^{2}\right)\sqrt{1-{z}^{2}}$, that depends on solitons population imbalance, z. The nonlinear effects in the tunneling process vanish for z² ≪ 1. In this limit the BJJ- and SJJ-models possess the same features.

However, when z² → 1, the effective tunneling rate, κ_eff, goes to zero and effective parameter Λ_eff essentially increases. In other words, we deal here with self-tuning effect for effective tunneling parameter, κ_eff, that establishes various regimes for interacting solitons. As we can see below, this limit is of particular interest for the quantum (superposition) states formation, which we specify in the paper.

The equations $\dot {\theta }=-\partial {H}_{\text{SJJ}}/\partial z$ and $\dot {z}=\partial {H}_{\text{SJJ}}/\partial \theta$ for canonical variables, z and θ, can be obtained with the Hamiltonian (13) as [42, 44]:

$$\begin{equation}\dot {z}=\left(1-{z}^{2}\right)\left(1-0.21{z}^{2}\right)\mathrm{sin}\left[\theta \right];\end{equation} \tag{ 15a }$$

$$\begin{equation}\dot {\theta }={\Lambda}z-2z\left(1.21-0.42{z}^{2}\right)\mathrm{cos}\left[\theta \right].\end{equation} \tag{ 15b }$$

In equation (15) dots denote derivatives with respect to effective (dimensionless) time τ = κNt. These equations permit non-trivial steady-state solutions

$$\begin{equation}{z}_{{\pm}}={\pm}\sqrt{\frac{1}{0.42}\left(1.21-\frac{{\Lambda}}{2}\right)};\end{equation} \tag{ 16a }$$

$$\begin{equation}\mathrm{cos}\left[\theta \right]=1,\end{equation} \tag{ 16b }$$

which are valid for 1.58 < Λ ⩽ 2.42. For another specific value of phase θ we have

$$\begin{equation}{z}_{{\pm}}={\pm}1;\end{equation} \tag{ 17a }$$

$$\begin{equation}\mathrm{cos}\left[\theta \right]=\frac{{\Lambda}}{1.58},\end{equation} \tag{ 17b }$$

that requires 0 < Λ ⩽ 1.58.

It is important that classical mean-field theory admits realization of one of solutions: either with z₊ or with z₋. On the contrary, the quantum approach, that we are going to discuss below, permits the simultaneous realization of these states, being a clearly superposition state.

The quantum theory approach proposes a second quantization procedure for the classical effective Hamiltonians (8) and (13), and relevant field amplitudes. Notice that especially for nonlinear Hamiltonians containing phase-dependent terms this procedure can be non-unique because of the absence of rigorous phase quantization procedure and non-commutativity of relevant operators. However, as known, various approaches to Hamiltonian quantization permit difference in terms by the factor proportional to 1/N and less, which is negligibly small for large N. This difference occurs due to non-commutativity of annihilation and creation operators forming the Hamiltonian.

To be more specific, we start with quantization of the classical Hamiltonian (8) for the BJJ-model. The quantization procedure that we propose here is based on the approach described in [76]. We construct a quantum Hamiltonian from equation (8) by introducing particle number operators ${\hat{N}}_{1}={\hat{a}}^{{\dagger}}\hat{a}$ and ${\hat{N}}_{2}={\hat{b}}^{{\dagger}}\hat{b}$, respectively; one can define the reduced particle number difference operator $\hat{z}=\frac{1}{N}\left({\hat{b}}^{{\dagger}}\hat{b}-{\hat{a}}^{{\dagger}}\hat{a}\right)$, where $\hat{a}$ (${\hat{a}}^{{\dagger}}$) and $\hat{b}$ (${\hat{b}}^{{\dagger}}$) are usual Bosonic annihilation (creation) operators for condensate modes in two-mode approximation, cf [76]. Then one can suppose in (8) that $\sqrt{1-{z}^{2}}\enspace \mathrm{cos}\left[\theta \right]$ corresponds to $\frac{1}{N}\left({\hat{a}}^{{\dagger}}\hat{b}+{\hat{b}}^{{\dagger}}\hat{a}\right)$ if we assume decomposition $\hat{a}=\sqrt{\frac{N}{2}}\sqrt{1-\hat{z}}\enspace {\text{e}}^{\text{i}\hat{\theta }/2}$ and $\hat{b}=\sqrt{\frac{N}{2}}\sqrt{1+\hat{z}}\enspace {\text{e}}^{-\text{i}\hat{\theta }/2}$ for annihilation operators $\hat{a}$ and $\hat{b}$, respectively. The quantized Hamiltonian obtained from (8) for the BJJ-model now has the form

$$\begin{equation}{\hat{H}}_{\text{BJJ}}=\kappa N\left(-\frac{\lambda }{2}{\hat{z}}^{2}-\frac{1}{N}\left({\hat{a}}^{{\dagger}}\hat{b}+{\hat{b}}^{{\dagger}}\hat{a}\right)\right).\end{equation} \tag{ 18 }$$

For the SJJ-model we start from the classical Hamiltonian (13) rewriting it as

$$\begin{equation}{H}_{\text{SJJ}}=\kappa N\left(-\frac{{\Lambda}}{2}{z}^{2}-\left(1-0.21{z}^{2}\right)\enspace \left(\sqrt{1-{z}^{2}}\enspace \mathrm{cos}\left[\theta \right]\right)\enspace \enspace \sqrt{1-{z}^{2}}\right),\end{equation} \tag{ 19 }$$

that is suitable for quantization procedure described above. The last term in (19) we treat by using (formal) Taylor expansion representation for the operator $\sqrt{1-{\hat{z}}^{2}}={\sum }_{k=0}^{\infty }{C}_{1/2}^{k}{\left(-1\right)}^{k}{\hat{z}}^{2k}$. Finally, the quantum version of Hamiltonian (19) reads as

$$\begin{equation}{\hat{H}}_{\text{SJJ}}=\kappa N\left(-\frac{{\Lambda}}{2}{\hat{z}}^{2}-\frac{1}{2N}\left(\sum _{k=0}^{\infty }{C}_{1/2}^{k}{\left(-1\right)}^{k}\left(1-0.21{\hat{z}}^{2}\right)\left({\hat{a}}^{{\dagger}}\hat{b}+{\hat{b}}^{{\dagger}}\hat{a}\right){\hat{z}}^{2k}+\mathrm{h}.\mathrm{c}.\right)\right).\end{equation} \tag{ 20 }$$

Thereafter, we refer to equation (20) as quantum SJJ-model that characterizes the coupled quantum bright soliton problem in a two-mode approximation. Unlike familiar BJJ-model, atom number difference dependent tunneling is also at the core of the quantum SJJ-model, cf (18).

Formally, from equation (20) it is clear that Λ-parameter for the SJJ-model plays the same role as parameter λ in the conventional two-mode BJJ-model, see (18). However, it is possible to see that Λ = ℘λ² ∼ N², where we introduced parameter ℘ ≡ πκ/2ν that characterizes differences in the BJJ- and SJJ-models experimental realization.

Now we consider the features of macroscopic states for coupled solitons possessing a quite large number of atoms in two-mode approximation. These features can be revealed in the framework of the Hartree approximation applied to solitons [46, 77] and quantum dimers [78]. In particular, the Hartree approximation presumes factorization of N-particle state, which allows to represent it as $\left\vert {\psi }_{N}\right\rangle =\left\vert \psi \right\rangle \otimes \left\vert \psi \right\rangle \otimes \cdots \left\vert \psi \right\rangle$, where $\left\vert \psi \right\rangle$ is a single particle (qubit) state [14, 46].

To be more specific we are interested in the quantum field variational approach to the model described by (20) exploiting an atom-coherent state ansatz in the form

$$\begin{equation}\left\vert {\psi }_{N}\right\rangle =\frac{1}{\sqrt{N!}}{\left[\alpha {\hat{a}}^{{\dagger}}+\beta {\hat{b}}^{{\dagger}}\right]}^{N}\left\vert 0\right\rangle ,\end{equation} \tag{ 21 }$$

where $\left\vert 0\right\rangle \equiv {\left\vert 0\right\rangle }_{a}{\left\vert 0\right\rangle }_{b}$ is a two-mode vacuum state and α, β are unknown (variational) wave functions obeying normalization condition |α|² + |β|² = 1. The mean energy in respect of state (21) is obtained from (20) and looks like

$$\begin{equation}E\left(\alpha ,{\alpha }^{{\ast}},\beta ,{\beta }^{{\ast}}\right)=\left\langle {\hat{H}}_{\text{SJJ}}\right\rangle =\kappa N\left(-\frac{{\Lambda}}{2}{S}^{2}-2\vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)\left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)\right),\end{equation} \tag{ 22 }$$

where $S\equiv \left\langle {\psi }_{N}\vert \hat{z}\vert {\psi }_{N}\right\rangle =\vert \beta {\vert }^{2}-\vert \alpha {\vert }^{2}$.

The energy E in equation (22) implies Lagrangian, L, given by

$$\begin{equation}L=\frac{\mathrm{i}}{2}\left[{\alpha }^{{\ast}}\dot {\alpha }-{\dot {\alpha }}^{{\ast}}\alpha +{\beta }^{{\ast}}\dot {\beta }-{\dot {\beta }}^{{\ast}}\beta \right]-\kappa N\left(-\frac{{\Lambda}}{2}{S}^{2}-2\vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)\left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)\right).\end{equation} \tag{ 23 }$$

By using equation (23) we derive the equations of motion for the field variables α, β in the form,

$$\begin{align}\hfill \mathrm{i}\dot {\alpha }& =\kappa N\left(\alpha S{\Lambda}-\alpha \frac{\vert \beta \vert }{\vert \alpha \vert }\left(1-0.21{S}^{2}\right)\left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)\right.\hfill \\ \hfill & \quad \left.-0.84\alpha S\vert \alpha {\Vert}\beta \vert \left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)-2\beta \vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)\right);\hfill \end{align} \tag{ 24a }$$

$$\begin{align}\hfill \mathrm{i}\dot {\beta }& =\kappa N\left(-\beta S{\Lambda}-\beta \frac{\vert \alpha \vert }{\vert \beta \vert }\left(1-0.21{S}^{2}\right)\left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)\right.\hfill \\ \hfill & \quad \left.+0.84\beta S\vert \alpha {\Vert}\beta \vert \left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)-2\alpha \vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)\right),\hfill \end{align} \tag{ 24b }$$

where it is assumed that |α| ≠ 0, or |β| ≠ 0 that implies S ≠ ±1.

Setting in (24) stationary solutions α(t) ≡ α e^−iκNμt; β(t) ≡ β e^−iκNμt we get

$$\begin{align}\hfill \mu & =S{\Lambda}-\frac{\vert \beta \vert }{\vert \alpha \vert }\left(1-0.21{S}^{2}\right)\left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)\hfill \\ \hfill & \quad \left.-0.84S\vert \alpha {\Vert}\beta \vert \left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)-2\frac{\beta }{\alpha }\vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)\right);\hfill \end{align} \tag{ 25a }$$

$$\begin{align}\hfill \mu & =-S{\Lambda}-\frac{\vert \alpha \vert }{\vert \beta \vert }\left(1-0.21{S}^{2}\right)\left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)\hfill \\ \hfill & \quad \left.+0.84S\vert \alpha {\Vert}\beta \vert \left({\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha \right)-2\frac{\alpha }{\beta }\vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)\right),\hfill \end{align} \tag{ 25b }$$

equation (25) lead to a single equation:

$$\begin{equation}S{\Lambda}-S\left(1.42-0.63{S}^{2}\right)\frac{{\alpha }^{{\ast}}\beta +{\beta }^{{\ast}}\alpha }{2\vert \alpha {\Vert}\beta \vert }-\left(\frac{\beta }{\alpha }-\frac{\alpha }{\beta }\right)\vert \alpha {\Vert}\beta \vert \left(1-0.21{S}^{2}\right)=0.\end{equation} \tag{ 26 }$$

In general, the functions α and β in equation (26) can be represented as $\alpha \equiv \vert \alpha \vert {\text{e}}^{\text{i}{\theta }_{a}}$, $\beta \equiv \vert \beta \vert {\text{e}}^{\text{i}{\theta }_{b}}$, where θ = θ_b − θ_a is phase difference. In particular, equation (26) reduces to

$$\begin{equation}S\left(\frac{{\Lambda}}{2}+1.21-0.42{S}^{2}\right)=0\end{equation} \tag{ 27 }$$

for θ = π, and to

$$\begin{equation}S\left(\frac{{\Lambda}}{2}-1.21+0.42{S}^{2}\right)=0,\end{equation} \tag{ 28 }$$

for θ = 0, respectively. Notice, if we change the sign of tunneling rate, κ, the shift of phase by π is expected for (27) and (28).

Equations (27) and (28) possess a trivial solution with S ≡ S₀ = |β₀|² − |α₀|² = 0 that is valid for any Λ. Equation (28) permits two additional solutions, namely solutions ${S}_{{\pm}}={\pm}\sqrt{\frac{1}{0.84}\left(2.42-{\Lambda}\right)}$ exist for 1.58 ⩽ Λ < 2.42, cf (16a). Variables α and β which correspond to these solutions are

$$\begin{equation}{\alpha }_{0}={\beta }_{0}=\frac{1}{\sqrt{2}};\end{equation} \tag{ 29a }$$

$$\begin{equation}{\beta }_{{\pm}}={\alpha }_{\mp }=\sqrt{\frac{1}{2}\left(1{\pm}\sqrt{1-{X}^{2}}\right)},\end{equation} \tag{ 29b }$$

where $X\equiv \sqrt{1-{S}_{{\pm}}^{2}}=\sqrt{\frac{1}{0.84}\left({\Lambda}-1.58\right)}$.

The energies corresponding to equation (29) are given in the form

$$\begin{equation}{E}_{0}=-\kappa N;\end{equation} \tag{ 30a }$$

$$\begin{equation}{E}_{{\pm}}=\kappa N\left(0.30{{\Lambda}}^{2}-1.44{\Lambda}+0.74\right).\end{equation} \tag{ 30b }$$

It is important that E_± ⩾ E₀ for any Λ, when S_± solution exists.

Substituting (29b) into (21) one can obtain 'two halves' of the SC-state

$$\begin{equation}\left\vert {\psi }_{N,{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle =\frac{1}{\sqrt{N!}}{\left[{\alpha }_{{\pm}}{\hat{a}}^{{\dagger}}+{\beta }_{{\pm}}{\hat{b}}^{{\dagger}}\right]}^{N}\left\vert 0\right\rangle .\end{equation} \tag{ 31 }$$

The states (31) can form macroscopic superpositions

$$\begin{equation}\left\vert {\psi }_{{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle =\frac{1}{\sqrt{2\left(1{\pm}{\epsilon}\right)}}\left(\left\vert {\psi }_{N,+}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle {\pm}\left\vert {\psi }_{N,-}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle \right),\end{equation} \tag{ 32 }$$

possessing the same energy (30b).

The states $\left\vert {\psi }_{{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle$ in (32) are mutually orthogonal SC-states defined within the domain 1.58 < Λ < 2.42, cf (16a). The -parameter in (32) characterizes states $\left\vert {\psi }_{N,{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle$ overlapping (distinguishability) and is defined as

$$\begin{equation}{\epsilon}=\left\langle {\psi }_{N,+}^{\left(\mathrm{S}\mathrm{C}\right)}\vert {\psi }_{N,-}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle ={X}^{N}.\end{equation} \tag{ 33 }$$

At = 1 (X = 1, Λ ≃ 2.42) the states $\left\vert {\psi }_{N,{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle$ cannot be distinguished at all. In this case, we have ${\beta }_{{\pm}}={\alpha }_{\mp }=1/\sqrt{2}$ in (29b) and SC-states definition (32) becomes useless.

In another limit, setting X → 0 in (33) we obtain → 0 that corresponds to macroscopically distinguishable states $\left\vert {\psi }_{N,{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle$.

As it follows from (33), at = 0 (X = 0) the states $\left\vert {\psi }_{N,{\pm}}^{\left(\mathrm{S}\mathrm{C}\right)}\right\rangle$ are orthogonal, i.e. maximally distinguishable. In this limit the SC-states approaches the N00N-state that occurs for X = 0 (Λ = 1.58), see (29b), (31), and (32).

In other words, the N00N-state appears for S = ±1 that implies |α| = 0 or |β| = 0 in (21)–(23) (see also equation (17)). These solutions are

$$\begin{equation}\left\vert {\psi }_{N,-}^{\left(N00N\right)}\right\rangle ={\mathrm{e}}^{\mathrm{i}N{\theta }_{a}}{\left\vert N\right\rangle }_{a}{\left\vert 0\right\rangle }_{b}=\frac{{\left({\hat{a}}^{{\dagger}}\right)}^{N}\enspace {\mathrm{e}}^{\mathrm{i}N{\theta }_{a}}}{\sqrt{N!}}\left\vert 0\right\rangle ;\end{equation} \tag{ 34a }$$

$$\begin{equation}\left\vert {\psi }_{N,+}^{\left(N00N\right)}\right\rangle ={\mathrm{e}}^{\mathrm{i}N{\theta }_{b}}{\left\vert 0\right\rangle }_{a}{\left\vert N\right\rangle }_{b}=\frac{{\left({\hat{b}}^{{\dagger}}\right)}^{N}\enspace {\mathrm{e}}^{\mathrm{i}N{\theta }_{b}}}{\sqrt{N!}}\left\vert 0\right\rangle ,\end{equation} \tag{ 34b }$$

where '±' correspond to the sign of S that defines two 'halves' of the N00N-state:

$$\begin{equation}\left\vert N00N\right\rangle =\frac{1}{\sqrt{2}}\left(\left\vert {\psi }_{N,-}^{\left(N00N\right)}\right\rangle +{\mathrm{e}}^{\mathrm{i}N\theta }\left\vert {\psi }_{N,+}^{\left(N00N\right)}\right\rangle \right).\end{equation} \tag{ 35 }$$

In (35) we omit the common phase factor ${\mathrm{e}}^{\mathrm{i}N{\theta }_{a}}$. State (35) possesses phase independent mean energy

$$\begin{equation}{E}^{\left(N00N\right)}=-\frac{1}{2}\kappa N{\Lambda}.\end{equation} \tag{ 36 }$$

The phase θ in (35) obey equation (17b) and defines the domain 0 < Λ ⩽ 1.58 for Λ-parameter, where the N00N-state can appear in the framework of the Hartree approach discussed above.

Noteworthy, the state (35) is fragile enough and rapidly collapses to states ${\left\vert N-1\right\rangle }_{a}{\left\vert 0\right\rangle }_{b}$ or ${\left\vert 0\right\rangle }_{a}{\left\vert N-1\right\rangle }_{b}$ even in the presence of one particle loss, cf [38]. As we see below, the SJJ-model exhibits some new specifics for the N00N-state formation in a purely quantum limit.

Now let us consider purely quantum features of the superposition states for the SJJ-model beyond the Hartree approach. This approximation presumes a quite large number of particles, that can be thousands or more in experiments with atomic condensates, see e.g. [79]. However, for the condensates possessing the negative scattering length the number of atoms is limited due to the collapsing effect. For the BECs confined in optical lattices and representing familiar Bose–Hubbard system the number of atoms varies from several tens to several hundreds per site, see e.g. [80], representing a mesoscopic system with a moderate number of particles. In respect of the number of particles, mesoscopic (in Greek 'meso' means 'middle') regime is located between microscopic (up to tens of particles) and macroscopic (thousand particles and more) regimes and represents a great interest for current studies in quantum physics [81]. The formation of efficient quantum mesoscopic superposition states represents one of challenging problems both in theory and experiment [41]. It is worth mentioning that the crossover from a coherent to Fock (or number squeezed) regime possessing some important quantum features plays an important role in this case [85]. In particular, variation of the tunneling rate between the lattice sites allows to achieve the crossover from a coherent to Fock regime (Mott-insulator state). However, as it is shown below the (self)tuning of this crossover may be much more efficient for the SJJ-model that represents the simplest (two-site) Bose–Hubbard system.

Here we are interested in mesoscopic limit for the SJJ-model illustrated in figure 1. We assume that the total number of atoms, N, is not too large, and the characteristic temperatures, T, T_c, relevant to the condensates are sufficiently low in comparison with the characteristic energy space, δE, between the neighboring low-energy quantum levels of the system [86]; T_c is the critical temperature of the phase transition to the BEC state. Physically it implies the importance how a small number of atoms behave in the presence of a large total number, N. In such a limit the analysis of the two-mode system described by Hamiltonians (18) and (20) has to be fully quantum.

We represent the quantum two-mode state, |Ψ(τ)⟩, in the Fock basis as (cf [16]):

$$\begin{equation}\vert {\Psi}\left(\tau \right)\rangle =\sum _{n=0}^{N}{A}_{n}\left(\tau \right)\vert N-n,n\rangle ,\end{equation} \tag{ 37 }$$

where |N − n, n⟩ ≡ |N − n⟩_a|n⟩_b denotes the two-mode atom-number state; the coefficients A_n(τ) obey normalization condition ${\sum }_{n=0}^{N}\vert {A}_{n}\left(\tau \right){\vert }^{2}=1$. The A_n(τ) fulfills the Schrödinger equation

$$\begin{equation}\mathrm{i}\frac{\mathrm{d}{A}_{n}\left(\tau \right)}{\mathrm{d}\tau }=\langle N-n,n\vert \hat{H}\vert {\Psi}\left(\tau \right)\rangle .\end{equation} \tag{ 38 }$$

In connection with quantum Hamiltonians (18) and (20) the equation (38) can be represented as:

$$\begin{equation}\mathrm{i}{\dot {A}}_{n}={a}_{n}{A}_{n+1}+{\beta }_{n-1}{A}_{n-1},\end{equation} \tag{ 39 }$$

where the following notations are introduced:

$$\begin{equation}{\alpha }_{n}=-\frac{\lambda }{2}{\left(\frac{2n}{N}-1\right)}^{2},\end{equation} \tag{ 40a }$$

$$\begin{equation}{\beta }_{n}=-\frac{1}{N}\sqrt{\left(n+1\right)\left(N-n\right)}\end{equation} \tag{ 40b }$$

for the BJJ-model, and

$$\begin{align}\hfill {\alpha }_{n}& =-\frac{{\Lambda}}{2}{\left(\frac{2n}{N}-1\right)}^{2},\hfill \\ \hfill {\beta }_{n}& =-\frac{1}{{N}^{2}}\left(\left[1-0.21{\left(\frac{2n}{N}-1\right)}^{2}\right]\left(n+1\right)\sqrt{\left(N-n\right)\left(N-n-1\right)}\right.\hfill \\ \hfill & \quad \left.+\left[1-0.21{\left(\frac{2\left(n+1\right)}{N}-1\right)}^{2}\right]\left(N-n\right)\sqrt{n\left(n+1\right)}\right)\hfill \end{align} \tag{ 41a }$$

for the SJJ one.

Then let us focus on stationary solution ${A}_{n}\left(\tau \right)={A}_{n}\enspace {\mathrm{e}}^{-\mathrm{i}{E}_{n}\tau }$ of (39) that enables to find out energy spectrum E_n of the quantum Hamiltonian (20) for the SJJ-model.

Figure 2(a) demonstrates numerical calculations of the spectrum for N = 300 particles as a function of Λ-parameter. Practically, the dependence on Λ can be examined with the help of a Feshbach resonance through controlling atom–atom scattering length (u-parameter) [59, 72].

Zoom In Zoom Out Reset image size

The black dashed curve in figure 2(a) indicates the steady-state solutions for the SJJ-model in the energy dense region, where the Hartree approach is satisfied, cf [78]; this curve characterizes the energies (36) and (30b), and is valid for 0 < Λ ⩽ 1.58 and 1.58 < Λ ⩽ 2.42, respectively. The blue dashed line corresponds to the energy E₀ in accordance with equation (30a).

Main peculiarities with the coupled soliton energies occur in the vicinity of point Λ ≡ Λ_c ≈ 2.000 9925. In particular, at Λ = Λ_c the SJJ-system loses its coherence features entering Fock (Mott-insulator) regime modifying particle number fluctuations, cf [5, 18, 67]. Notice, that in this case, the SJJ-system treated classically changes its regime from Rabi-like oscillations to self-trapping, cf [8, 42, 44]. In particular, at Λ < Λ_c energetically favorable state for the quantum SJJ-system in figure 2(a) is described by the blue dashed (or the red solid) line exhibiting the minimal energy E₀. As seen from figure 2(a), at Λ ⩾ Λ_c the energies E_±, E^(N00N), which correspond to the superposition states in the Hartree approximation, are comparable with E₀ but placed a little bit above it.

Situation changes completely in a purely quantum limit when parameter Λ crosses its critical value Λ_c. The red solid line in figure 2(a) clearly shows that the quantum N00N-state is energetically favorable at Λ ⩾ Λ_c for the SJJ ground state (it is demonstrated below that this state is not exactly a N00N-state as it happens in the semiclassical limit).

To understand the main properties of the SJJ-model at low energies (n = 0) we examine the simplest case with N = 2 particles. The diagonalization of (20) gives three energy eigenvalues obtained analytically:

$$\begin{equation}{E}_{0}=-\kappa N\frac{{\Lambda}}{2};\end{equation} \tag{ 42a }$$

$$\begin{equation}{E}_{{\pm}}=\frac{1}{4}\kappa N\left[-{\Lambda}{\pm}\sqrt{{{\Lambda}}^{2}+2.5}\right].\end{equation} \tag{ 42b }$$

Similar results may be obtained for the BJJ-model. Diagonalization of the relevant Hamiltonian (18) leads to

$$\begin{equation}{E}_{0}=-\kappa N\frac{\lambda }{2};\end{equation} \tag{ 43a }$$

$$\begin{equation}{E}_{{\pm}}=\frac{1}{4}\kappa N\left[-\lambda {\pm}\sqrt{{\lambda }^{2}+16}\right].\end{equation} \tag{ 43b }$$

The energy eigenvalues (42) and (43) are plotted in figure 2(b) as functions of λ that is relevant to the key parameter combination uN/κ and inherent to both JJ-models. The gap ΔE = E₀ − E₋ between two lowest energy eigenstates obtained from equations (42) and (43) looks like

$$\begin{equation}{\Delta}{E}_{\text{SJJ}}=\frac{1}{4}\kappa N\left[-{\Lambda}+\sqrt{{{\Lambda}}^{2}+2.5}\right];\end{equation} \tag{ 44a }$$

$$\begin{equation}{\Delta}{E}_{\text{BJJ}}=\frac{1}{4}\kappa N\left[-\lambda +\sqrt{{\lambda }^{2}+16}\right].\end{equation} \tag{ 44b }$$

From equation (44) it follows ΔE_SJJ ≃ 0.4κN and ΔE_BJJ = 4κN for vanishing nonlinearity u ≈ 0. Notably, the limiting case with u = 0 is not valid for the SJJ-model since bright solitons do not exist in the absence of nonlinearity. On the contrary, for the relatively large nonlinear parameter, such as λ, Λ ≫ 1, from (44) we obtain for the energy gaps ΔE_BJJ ≈ 4κ²/u, and ΔE_SJJ ≈ 5κ²/Nu² ≈ 5ΔE_BJJ/4Nu, respectively. The suppression of this gap with increasing atom-atom scattering length (u-parameter) is clearly seen in figure 2(b).

For the moderate values of λ and Λ parameters ΔE_SJJ is N times larger than ΔE_BJJ. The lifetime ${T}_{{E}_{-}}$ for transition ΔE is proportional to ℏ/ΔE₋, see e.g. [87]. It is evident that in general even for a small number of particles the lifetime of soliton junctions may be sufficiently larger than for the conventional BJJ-model, cf [78].

In figure 3 we represent the ground state probability distributions, |A_n|², (eigenfunctions of (20) with the minimal energies) for different values of the key parameters Λ and λ for the SJJ- and BJJ-models, respectively. The plots for the SJJ-model are established by red and can be easily compared with the ones for the BJJ-model indicated by blue. Figure 3(a) demonstrates a small difference in distributions for the SJJ- and BJJ-models in an ideal (non-interacting) atomic gas limit. Since the total number of particles is large enough, the initial distribution approaches Gaussian for the BJJ-model.

Zoom In Zoom Out Reset image size

The differences between distributions considered for the SJJ- and BJJ-models become evident and important in figures 3(b)–(d) due to the particle number dependence of the effective tunneling rate, κ_eff, see (14). Figures 3(b) and (c) exhibit features of examined JJ-systems immediately before and after the crossover, respectively. In the vicinity of the critical value Λ_c ≃ 2 the distribution for the SJJ-model becomes significantly non-Gaussian, see figure 3(b).

The SJJ- and BJJ-models behave differently at the crossover region for parameters Λ ⩾ Λ_c and λ ⩾ λ_c, respectively. In particular, deepening in ground state probability |A_n|² takes place for the BJJ-system with increasing λ.

Behavioral scenarios after crossover points Λ_c ≃ 2 and λ_c ≃ 1 for the SJJ- and BJJ-models, respectively, are represented in figures 3(c) and (d). In particular, two probability peaks corresponding to the macroscopic SC-state arise for BJJs. Further increasing of λ results in the growth of the 'distance' between the cats 'moving' to the 'edges', n = 0 and n = N, see figure 3(d). The N00N-state,

$$\begin{equation}\left\vert N00N\right\rangle =\frac{1}{\sqrt{2}}\left({\left\vert N\right\rangle }_{a}{\left\vert 0\right\rangle }_{b}+{\text{e}}^{\text{i}{\theta }_{N}}{\left\vert 0\right\rangle }_{a}{\left\vert N\right\rangle }_{b}\right),\end{equation} \tag{ 45 }$$

occurs for the BJJ-model in the limit of a very large λ-parameter, see below figure 4(b) and cf [16, 18].

Zoom In Zoom Out Reset image size

In general, the SJJ-model demonstrates evident advantages to achieve a N00N-state. Roughly speaking, state (37) for the SJJ-model after the crossover represents a superposition of entangled Fock states possessing a large N00N-state component—see figure 3(c). A physical explanation for this phenomenon is illustrated as follows.

First, for a given tunneling rate, κ, the Λ-parameter depends on N², which allows to achieve a required crossover point for the moderate number of particles, N.

Second, as it follows from figures 3(b) and (c) the crossover for the SJJ-model happens abruptly with a slightly increasing Λ parameter. The dependence of the effective tunneling rate, κ_eff, on the atom number leads to a significant improvement of the 'edge' states occupation. In particular, as clearly seen in figure 3(c), two large peaks of probability distribution occur at the 'edges' with n = 0 and n = N, respectively. They correspond to the N00N-state. Simultaneously, there exist small (non-vanishing symmetric) peaks for n = 1, N − 1, n = 2, N − 2, etc, in the vicinity of 'edge' states—see the inset in figure 3(c).

Furthermore, the increment in Λ leads to the suppression of small peaks probabilities, but to the enhancement of the 'edge' states n = 0, N population. The distribution approaches approximately 'ideal' N00N-state (45) at moderate values of Λ, cf figure 3(d). However, as we show below, the presence of small components in the satellite entangled Fock state plays an important role for the resistance of state (37) to weak particle losses, cf [39].

The complete analysis of the highly nonclassical states formation requires investigating high (two- and much higher) order correlation functions. Entanglement and spin squeezing are important prerequisites of non-classical behavior of multiparticle atomic system. To be more specific we examine here so-called m-order Hillery–Zubairy (HZ) criterion [17, 23] defining as

$$\begin{equation}{E}_{\mathrm{H}\mathrm{Z}}^{\left(m\right)}=1+\frac{\left\langle {\hat{a}}^{{\dagger}m}{\hat{a}}^{m}{\hat{b}}^{{\dagger}m}{\hat{b}}^{m}\right\rangle -{\left\vert \left\langle {\hat{a}}^{m}{\hat{b}}^{{\dagger}m}\right\rangle \right\vert }^{2}}{\left\langle {\hat{a}}^{{\dagger}m}{\hat{a}}^{m}\left({\hat{b}}^{m}{\hat{b}}^{{\dagger}m}-{\hat{b}}^{{\dagger}m}{\hat{b}}^{m}\right)\right\rangle }.\end{equation} \tag{ 46 }$$

The entanglement occurs when inequalities

$$\begin{equation}0{\leqslant}{E}_{\mathrm{H}\mathrm{Z}}^{\left(m\right)}{< }1\end{equation} \tag{ 47 }$$

are satisfied. The value ${E}_{\mathrm{H}\mathrm{Z}}^{\left(m\right)}=1$ should be discussed separately for some specific states. We examine two limiting cases for definition (46). In particular, for m = 1 we can recast equation (46) in terms of N-particle spin-operators as following:

$$\begin{equation}{E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}=\frac{\delta {J}_{X}+\delta {J}_{Y}}{\langle \hat{J}\rangle },\end{equation} \tag{ 48 }$$

where $\delta {J}_{j}\equiv \left\langle {\left({\Delta}{\hat{J}}_{j}\right)}^{2}\right\rangle =\left\langle {\hat{J}}_{j}^{2}\right\rangle -{\left\langle {\hat{J}}_{j}\right\rangle }^{2}$ is the variance of fluctuations for jth component of spin, established by operators

$$\begin{align}\hfill {\hat{J}}_{X}& =\frac{1}{2}\left({\hat{a}}^{{\dagger}}\hat{b}+{\hat{b}}^{{\dagger}}\hat{a}\right); {\hat{J}}_{Y}=\frac{1}{2i}\left({\hat{b}}^{{\dagger}}\hat{a}-{\hat{a}}^{{\dagger}}\hat{b}\right); {\hat{J}}_{Z}=\frac{1}{2}\left({\hat{a}}^{{\dagger}}\hat{a}-{\hat{b}}^{{\dagger}}\hat{b}\right);\hfill \\ \hfill \hat{J}& =\frac{1}{2}\left({\hat{a}}^{{\dagger}}\hat{a}+{\hat{b}}^{{\dagger}}\hat{b}\right)=\frac{1}{2}\hat{N}.\hfill \end{align} \tag{ 49 }$$

Combining equation (46) with (47) one can see that strong entanglement can be obtained if variances δJ_X and δJ_Y are minimal for a given atom number, N. Strictly speaking, δJ_∥ = δJ_X + δJ_Y represents the variance of fluctuations of atomic spin in XY-plane. Red curve in figure 4(a) demonstrates behavior of ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}$ versus vital parameter Λ. Since the SJJ-model is invalid for ideal atomic gases, we represent a relevant curve dashed in the vicinity of u = 0. For the BJJ-model (see blue curve in figure 4(b)) in this limit we can exploit the two-mode condensate state like (21) with |α|² = |β|² = 0.5 that gives ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}=0.5$.

The EPR steering entanglement is achieved in the region, where ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}{< }0.5$, cf [21, 25]. From figure 4(a) it is evident that this area occurs close to the point Λ_c ≃ 2. Noteworthy, the sharp behavior of ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}$ within the crossover region appear in the absence of losses, see the inset in figure 4(a). Evidently, broadening of this region is expected if one- and, especially, three-body losses are taken into account.

The minimal value ${E}_{\mathrm{H}\mathrm{Z},\mathrm{min}}^{\left(1\right)}$ of parameter ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}$, that we denote as C_J, is determined through uncertainties in J_X and J_Y components which do not commute. The C_J-coefficient is estimated from inequality

$$\begin{equation}{E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}{\geqslant}{C}_{J}.\end{equation} \tag{ 50 }$$

In figure 5 we plot the parametric dependence of the normalized C_J-coefficient on a spin vector 'length', J, that is simply a half of the particle number, N/2. In particular, the BJJ-model exhibits the largest depth of steering determined by C_J.

Zoom In Zoom Out Reset image size

The coefficient C_J theoretically can be estimated by using the Gaussian distribution function possessing width, σ, (see figure 3(b)), cf [17]. In particular, for the BJJ-model σ_BJJ = (2J)^−2/3 and ${C}_{J}^{\left(\mathrm{B}\mathrm{J}\mathrm{J}\right)}\propto 1/{\sigma }_{\text{BJJ}}=0.6{J}^{-1/3}$ giving the minimum of the blue curve in figure 4(b) about 0.163. The magnitude ${C}_{J}^{\left(\mathrm{S}\mathrm{J}\mathrm{J}\right)}$ (that is about 0.319 24) for the SJJ-model in figure 4(a) (red curve) is not so simple to estimate analytically because of essentially nonlinear behavior happening close to the value Λ_c ≃ 2 (or, λ_c ≃ 1). However, as clearly seen from figure 3(b) the distribution for the SJJ-model is sufficiently narrower than for the BJJ one, i.e. σ_SJJ < σ_BJJ, and we expect ${C}_{J}^{\left(\mathrm{S}\mathrm{J}\mathrm{J}\right)}{ >}{C}_{J}^{\left(\mathrm{B}\mathrm{J}\mathrm{J}\right)}$ to take place for curves exhibiting ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}$ in figures 4(a) and (b).

Notably, the suppression of variance δJ_∥ manifests planar squeezing defined as (cf [32])

$$\begin{equation}\delta {J}_{{\Vert}}{< }{J}_{{\Vert}},\end{equation} \tag{ 51 }$$

where ${J}_{{\Vert}}\equiv \sqrt{{\left\langle {\hat{J}}_{X}\right\rangle }^{2}+{\left\langle {\hat{J}}_{Y}\right\rangle }^{2}}$ is the magnitude of the spin component in XY-plane. In particular, for both SJJ- and BJJ-models we have $\left\langle {\hat{J}}_{X}\right\rangle \approx J$, $\left\langle {\hat{J}}_{Y}\right\rangle =\left\langle {\hat{J}}_{Z}\right\rangle =0$ within the domain 0 < Λ ⩽ 2 for figure 4(a) and 0 ⩽ λ ⩽ 1 for figure 4(b), respectively. Thus, in this area ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}$ reflects the planar spin squeezing that occurs due to significant enhancement of fluctuations in the Z-component of spin [17].

It is important to stress that in figure 4(a) (green dashed line) the SJJ-system undergoes the abrupt transition to N00N-state (45) at the crossover region for Λ_c ≃ 2 that is consistent with figure 3(c). The ${E}_{\mathrm{H}\mathrm{Z}}^{\left(m\right)}$-parameter calculated for m = 1 and m = N by using N00N-state (45) gives ${E}_{\mathrm{H}\mathrm{Z}}^{\left(1\right)}=1$, ${E}_{\mathrm{H}\mathrm{Z}}^{\left(m\right)}=0.5$, respectively.

From figure 4(b) (green dashed line) it is evident that the transition to the N00N-state for the BJJ-model occurs upon the large scales of λ ≫ 1.

We now outline the possibility of experimental observation of N00N-states with solitons with feasible condensate parameters. To be more specific we discuss lithium condensate solitons [54]. In addition, we would like to mention the work [55] where authors recently studied lithium atomic soliton collisions on relatively long times, i.e. on 32 ms and more. In particular, we focus on two weakly coupled condensates of N = 300 atoms (|a_sc|≃ 1.4 nm, uN ∼ 1.88), confined in an asymmetric trap as it is shown in figure 1, with ω_x/2π = 70 Hz and radial trapping frequency ω_⊥/2π = 700 Hz that imposes a_⊥≃ 1.4 μm. The solitons in [54, 55] observed within characteristic time 8 ÷ 32 ms. The Λ parameter achieves its critical value, Λ_c ≃ 2, at tunneling rate $\vert \mathcal{K}\vert /2\pi \simeq 77\enspace \mathrm{H}\mathrm{z}$ or $\vert \mathcal{K}\vert \simeq 3.7\enspace \mathrm{n}\mathrm{K}$, which correlates the current experiments with atomic JJs, see e.g. [67]. Notice, such a tunneling rate promotes effective coupling of solitons within their observation time. Some moderate tuning of $\mathcal{K}$ or scattering length |a_sc| allows to obtain the crossover to the N00N-state with Λ ⩾ Λ_c.

For the BJJ-model the critical value λ_c = 1 is achieved with $\vert \mathcal{K}\vert /2\pi \simeq 83\enspace \mathrm{H}\mathrm{z}$. However, the N00N-state requires essential enhancement of λ in the experiment (60 times and more, see inset in figure 4(b)) that may require a non-trivial experimental task. On the other hand, SC-states can be observed for the BJJ-model with λ ⩾ 1, see figure 3(d) and [41].

Let us briefly discuss the loss issues for the SJJ-model described in this work. This problem is important for applications of the SJJ-model in current quantum technologies, especially in quantum metrology, where N00N-states play an essential role, cf [33, 38, 82].

Typically, it is instructive to consider one- and three-body losses occurring in condensates [83]. As shown in [84], the role of various losses becomes critical at a relatively large number of particles in the condensate. For example, in [84] it is shown that spin squeezing degrades for a bimodal condensate with N ⩾ 10⁴ atoms. However, the N00N-state formation in the presence of weak losses requires some special attention in this work.

First, let us estimate possible losses for the SJJ-model as a prerequisite for N00N-state formation. Since vital parameter, Λ, is proportional to N², many-body, especially three-body, inelastic collisions bringing to losses may be important [88].

Physically three-body losses occur as a result of three-body recombination process appearing with condensate atoms. Such processes primarily remove condensate particles possessing low kinetic energy and those located at the trap center. The three-body losses play a very important role in the vicinity of Feshbach resonance where scattering length essentially enhances [89].

As with one-body losses we expect here to have dynamical adiabatic changing of regimes from Fock to Josephson or Rabi, if losses are weak [90]. On the other hand, since effective tunneling parameter, κ_eff, depends on the particle number, the SJJ-model in the presence of three-body losses may be also relevant to some kind of dissipative tunneling problem [91].

In zero temperature limit the three-body recombination rate is ${L}_{3}\propto \hslash {a}_{\mathrm{s}\mathrm{c}}^{4}/m$ that approaches to L₃ ≃ 2.6 × 10⁻²⁸ cm⁶ s⁻¹ for lithium atoms, see e.g. [92]. A simple rate equation for particle number leads to decay of condensate atoms with the law as

$$\begin{equation}N\left(t\right)\simeq \frac{N\left(0\right)}{\sqrt{1+2{L}_{3}{\rho }^{2}t}},\end{equation} \tag{ 52 }$$

where N(0) is initial atom number in the condensate, ρ is atom number density that we assume here as a constant. Taking it as ρ ≃ 10¹³ cm⁻³ for the effective rate of three-body losses γ₃ = 2L₃ρ² we obtain γ₃ = 5.2 × 10⁻² s⁻¹.

The one-body losses rate may be estimated as γ₁ = 0.2 s⁻¹, see e.g. [41]. At times of soliton observation in [54, 55] we effectively obtain γ_1,3t ≪ 1 that allows to examine a non-perturbed SJJ-system. Notice, that in [54] observation time was less than 10 ms.

This conclusion is supported by recent experiments performed in a more realistic case at finite temperatures [93]. In particular, as it is shown in [94], three-body losses for lithium condensate atoms with scattering length |a_sc|≃ 10.6 nm and tuned by magnetic field at the temperature T = 5.2 μK are important at the time scales t ≃ 100 ms and more. Especially at times t ⩽ 32 ms relevant to matter-wave soliton observations in [55] the particle number reduction is practically negligible.

The analysis represented above is valid in the mean-field limit. In this sense, it is instructive to examine few particle losses in a purely quantum way for state (37). We can consider an approach to the losses, which is typically used in the framework of N00N-state applicability for quantum metrology tasks, see e.g. [38–40, 66]. Various strategies are proposed to achieve maximally accessible scaling for phase estimation; they manifest that in the presence of losses the ideal (balanced) N00N-state may be non-optimal. It is shown that entangled Fock states, unbalanced N00N-states, and some specific two-component states may be more useful in the presence of different level of losses [37, 39, 40]. Here we model the losses of condensate particles by means of the fictitious BS approach. This approach represents a useful tool to study coupling of quantum macroscopic superposition states with environment—see figure 6. The 'input' two-mode Fock state (37) after two BSs transforms into (cf [40])

$$\begin{equation}{\left\vert N-n\right\rangle }_{a}{\left\vert n\right\rangle }_{b}\to \sum _{{l}_{b}=0}^{N}\sum _{{l}_{a}=0}^{N-{l}_{b}}\sqrt{{B}_{{l}_{a},{l}_{b}}^{n}}{\left\vert N-n-{l}_{a}\right\rangle }_{a}{\left\vert n-{l}_{b}\right\rangle }_{b}\left\vert {l}_{a}\right\rangle \left\vert {l}_{b}\right\rangle ,\end{equation} \tag{ 53 }$$

where l_a and l_b are the numbers of particles lost from 'a' and 'b' wells, respectively. In (53) we introduce a coefficient

$$\begin{equation}{B}_{{l}_{a},{l}_{b}}^{n}=\left(\begin{matrix}\hfill N-n\hfill \\ \hfill {l}_{a}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill n\hfill \\ \hfill {l}_{b}\hfill \end{matrix}\right){\eta }_{a}^{N-n}{\left({\eta }_{a}^{-1}-1\right)}^{{l}_{a}}{\eta }_{b}^{n}{\left({\eta }_{b}^{-1}-1\right)}^{{l}_{b}}.\end{equation} \tag{ 54 }$$

that characterizes the relevant probabilities in the presence of particle losses. Here η_a and η_b (η_a,b ⩽ 1) are the transmissivities of BSs in the channels 'a' and 'b', respectively. Then state $\left\vert {{\Psi}}_{\text{out}}\right\rangle$ after particle losses reads as (see figure 6)

$$\begin{equation}\left\vert {{\Psi}}_{\text{out}}\right\rangle =\frac{1}{\sqrt{p}}\sum _{{l}_{b}=0}^{N}\sum _{{l}_{a}=0}^{N-{l}_{b}}\sum _{n={l}_{b}}^{N-{l}_{a}}{A}_{n}\sqrt{{B}_{{l}_{a},{l}_{b}}^{n}}{\left\vert N-n-{l}_{a}\right\rangle }_{a}{\left\vert n-{l}_{b}\right\rangle }_{b}\left\vert {l}_{a}\right\rangle \left\vert {l}_{b}\right\rangle ,\end{equation} \tag{ 55 }$$

where $p={\sum }_{{l}_{b}=0}^{N}{\sum }_{{l}_{a}=0}^{N-{l}_{b}}{\sum }_{n={l}_{b}}^{N-{l}_{a}}{\left\vert {A}_{n}\right\vert }^{2}{B}_{{l}_{a},{l}_{b}}^{n}$ is a normalization constant for state $\left\vert {{\Psi}}_{\text{out}}\right\rangle$.

Zoom In Zoom Out Reset image size

Let us suppose that we can somehow efficiently detect the numbers, l_a,b, of particles leaving the condensate state. In particular, this situation is similar to the conditional state preparation by means of subtraction of a small number of particles. For example, in quantum optical experiments real BS devices are used instead of fictitious ones, cf [95, 96]. Consider l_a = 1 and l_b = 0, i.e. if only one particle is detected in channel 'a' after BS_a, the conditional state $\vert {{\Psi}}_{\text{out}}^{\left(1,0\right)}\rangle$ obtained from (55) at η_a = η_b = η takes the form

$$\begin{equation}\left\vert {{\Psi}}_{\text{out}}^{\left(1,0\right)}\right\rangle =\frac{1}{\sqrt{p}}\sum _{n=0}^{N-1}{A}_{n}\sqrt{N-n}\sqrt{{\eta }^{N-1}\left(1-\eta \right)}{\left\vert N-n-1\right\rangle }_{a}{\left\vert n\right\rangle }_{b}.\end{equation} \tag{ 56 }$$

In figure 7 we represent the probabilities $\vert {C}_{1,0}^{n}{\vert }^{2}$ (${C}_{{l}_{a},{l}_{b}}^{n}\equiv \frac{1}{\sqrt{p}}{A}_{n}\sqrt{{B}_{{l}_{a},{l}_{b}}^{n}}$) for the SJJ- and BJJ-models versus n in the presence of weak and equal losses in each channel, η_a = η_b = 0.999. To be more specific, in figures 7(a) and (b), we use the same condensate parameters as in figures 3(c) and (d), respectively.

Zoom In Zoom Out Reset image size

The multiplier $\sqrt{N-n}$ in (56) plays a crucial role in the behavior of relevant probability $\vert {C}_{1,0}^{n}{\vert }^{2}$ that characterizes the conditional state $\vert {{\Psi}}_{\text{out}}^{\left(1,0\right)}\rangle$. In the presence of weak losses, the tendency to occupation of 'edge' states continues with increasing of parameter Λ, see figure 7(a). However, in this case the state with n = N cannot be occupied due to one particle loss, cf figure 3(c). Strictly speaking, the state with n = 0, which is state |N − 1⟩_a|0⟩_b, becomes macroscopically populated for the SJJ-model with probability $\vert {C}_{1,0}^{n}{\vert }^{2}\propto \sqrt{N}{\left\vert {A}_{N-1}\right\vert }^{2}{\eta }^{N-1}\left(1-\eta \right)$. On the contrary, the probability of occupation of new 'edge' state with n = N − 1 is more than $\sqrt{N}$ times lower. Physically, this case may by recognized as so-called partial 'collapse' of N00N-state (45) to state |N − 1⟩_a|0⟩_b, when one particle is detected in 'a' channel; some small occupation of the state with n = N − 1 still exists, see the inset in figure 7(b). Notice that the states with n = N − 1, N − 2, ... are poorly populated even without any losses—see insets to figures 3(c) and (d).

It is worth noticing that since the particles possess some (Gaussian-like) distribution in the vicinity of 'edges', the SC-states for the BJJ-model seem to be more robust to the particle losses, see the inset in figure 7(b) and cf figure 3(d).

However, for l_a = l_b = 1, i.e., if two particles are detected simultaneously after BSs in figure 6, the behavior of the conditional state (55) becomes fundamentally different. In this case for (55) we obtain:

$$\begin{equation}\left\vert {{\Psi}}_{\text{out}}^{\left(1,1\right)}\right\rangle =\frac{1}{\sqrt{p}}\sum _{n=1}^{N-1}{A}_{n}\sqrt{n\left(N-n\right){\left(\frac{{\eta }_{b}}{{\eta }_{a}}\right)}^{n}{\eta }_{a}^{N}\left({\eta }_{a}^{-1}-1\right)\left({\eta }_{b}^{-1}-1\right)}{\left\vert N-n-1\right\rangle }_{a}{\left\vert n-1\right\rangle }_{b}.\end{equation} \tag{ 57 }$$

In this limit, new 'edge' states with n = 1 and n = N − 1 are occupied with comparable probabilities. In particular, for η_a ≠ η_b superposition (57) turns to the unbalanced N00N-state that possesses N − 2 total number of particles and has a form

$$\begin{equation}\left\vert {{\Psi}}_{N00N}^{\left(1,1\right)}\right\rangle \propto \sqrt{\left(N-1\right)\left(1-{\eta }_{a}\right)\left(1-{\eta }_{b}\right)}\left(\sqrt{{\eta }_{a}^{N-2}}{\left\vert N-2\right\rangle }_{a}{\left\vert 0\right\rangle }_{b}+\sqrt{{\eta }_{b}^{N-2}}{\left\vert 0\right\rangle }_{a}{\left\vert N-2\right\rangle }_{b}\right).\end{equation} \tag{ 58 }$$

The unbalanced N00N-state (58) represents a useful tool for the quantum metrology purposes in the presence of weak losses if the corresponding transmissivities exceed a certain threshold value, cf [40].

To generalize these results, let us examine the state (55) in the case, when the numbers of lost particles, l_a and l_b, are undetected (and unknown), as it is shown in figure 6. Tracing out the ancillary modes l_a, l_b in (55) for the resulting state we obtain (cf [38])

$$\begin{equation}\left\vert {{\Psi}}_{\text{out}}\right\rangle =\sum _{{l}_{b}=0}^{N}\sum _{{l}_{a}=0}^{N-{l}_{b}}\sum _{n={l}_{b}}^{N-{l}_{a}}{C}_{{l}_{a},{l}_{b}}^{n}{\left\vert N-n-{l}_{a}\right\rangle }_{a}{\left\vert n-{l}_{b}\right\rangle }_{b}.\end{equation} \tag{ 59 }$$

In figure 8, we plot the probabilities $\vert {C}_{{l}_{a},{l}_{b}}^{n}{\vert }^{2}$ in N − n − l_a and n − l_b axis for the SJJ-model. Each column in figure 8 corresponds to the probability of state ${\left\vert N-n-{l}_{a}\right\rangle }_{a}{\left\vert n-{l}_{b}\right\rangle }_{b}$ occupation. The parameter Λ = 4 in figure 8 is the same as for figure 3(d). Two big (yellow) columns in figure 8 correspond to the N00N-state if no condensate particles are lost. In general, in the absence of losses the occupied modes are located on the main diagonal between these columns. In the presence of losses, the terms with l_a ≠ 0 and/or l_b ≠ 0 appear; they contain N' = N − l_a − l_b particles and are characterized by the columns on the lines parallel to the main diagonal in figure 8. In particular, the columns depicted in orange and green-colors in figure 8 form the set of the N00N-states with N' < N number of particles, see (58) and cf [38]. These superposition states still exist if equal numbers of particle l_a and l_b are measured after BSs in figure 6. However, relevant probabilities essentially vanish with increasing l_a and l_b.

Zoom In Zoom Out Reset image size