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

The Josephson junctions (JJs) are at the heart of modern quantum technologies and metrology. In this work we establish quantum features of an atomic soliton Josephson junction (SJJ) device, which consists of two weakly-coupled condensates with negative scattering length. The condensates are trapped in a double-well potential and elongated in one dimension. Starting with classical field theory we map for the first time a two-soliton problem onto the effective two-mode Hamiltonian and perform a second quantization procedure. Compared to the conventional Bosonic Josephson junction (BJJ) condensate system, we show that the SJJ-model in quantum domain exhibits unusual features due to its effective nonlinear strength proportional to the square of total particle number, $N^2$. A novel self-tuning effect for the effective tunneling parameter is also demonstrated in the SJJ-model, which depends on the particle number and rapidly vanishes as the JJ population imbalance increases. The formation of entangled Fock state superposition is predicted for the quantum SJJ-model, revealing dominant $N00N$-state components at the"edges"for $n=0, N$ particle number. We have shown that the obtained quantum state is more resistant to few particle losses from the condensates if tiny components of entangled Fock states are present in the vicinity of the major $N00N$-state component. This peculiarity of the quantum SJJ-model establishes an important difference from its semiclassical analogue obtained in the framework of Hartree approach.


I. INTRODUCTION
The Josephson junctions (JJs) represent an indispensable tool for current quantum technologies. The Josephson effect, that is at the core of various JJ-devices, presumes the macroscopic quantum effect of particle tunneling between two (or more) junctions; it was initially discovered in superconductors with Cooper pairs [1]. The JJs, which exhibit superfluidity, were demonstrated with 3 He, 4 He environments [2][3][4]. The atomic, so-called Bosonic Josephson junctions (BJJs), were shown for weakly-coupled atomic Bose-Einstein condensates (BECs) [5,6].
Since experimentally accessible temperatures of atomic BECs are very low (tens of mi-croKelvins and below), the dynamical regimes reveal an interplay between coherent Josephson tunneling and nonlinearity characterising weak atom-atom interaction in Born approximation [7,8]. The exciton-polariton BECs in narrow band semiconductor microstructures allow to obtain JJs at much higher temperatures like a few of Kelvins [9,10]. The photons possess tunneling between the contacts; however, in this case dissipation and temperature dependent phenomena must be taken into account [11]. Strictly speaking, such tunnelly-coupled waveguides represent purely photonic analogues of JJs far from thermal equilibrium [12,13].
From fundamental point of view, BJJ-devices may form macroscopic SC-states [26]. Some recent results associated with symmetry (rotational and inversion) are obtained for superposition states in [27,28]. Practically, BJJ quantum features are promising for realization of macroscopic quantum computation [29,30], as well as for quantum metrology and sensing purposes [31]. In particular, planar quantum squeezing, discussed in the framework of atomic JJs, allows to enhance the metrological sensitivity of measurements for arbitrary phase shifts [32]. Noteworthy, SC-states and especially maximally path-entangled N -particle N 00N -states represent a paramount tool to achieve the phase measurement at the Heisenberg sensitivity level [33].
The experimental generation of N 00N -states with a large number of particles, N , is one of the important and non-trivial problems, which significantly affects the state of art of current quantum metrology [34][35][36]. Investigating various aspects of appropriate N 00N -states is imperative in the presence of decoherence and losses, which can rapidly destroy N 00N -states [37,38]. For example, in the presence of losses balanced N 00N -states are non-optimal for the implementation in quantum metrology [39,40]. Designing macroscopic quantum states, which may be robust to losses, represents an essential and challenging task for current quantum technologies with atomic systems, cf. [41].
In this work we propose to create macroscopic entangled (Fock) states by means of tunnelycoupled bright atomic solitons, which demonstrate some useful features for quantum metrology purposes, cf. [42][43][44]. In particular, as reported in [43], quantum atomic solitons allow to surpass the Heisenberg limit and obtain N −3/2 accuracy in nonlinear phase shift estimation even with coherent probes in the framework of nonlinear quantum metrology. To overcome this limit we suggested SC-and N 00N -states formation in the scheme with two weakly-coupled condensate bright solitons.
However, in [42,43] we restricted ourselves by the Hartree approximation that cannot fully provide the advantages of the solitonic model in a quantum domain, cf. [44]. In this work we show that the system of coupled atomic solitons represents a new quantum soliton Josephson junction (SJJ) device with a nontrivial physical behavior in comparison with a conventional BJJ-device based on atomic condensates, which possess Gaussian-shape wave functions. In particular, paper aims to demonstrate that SJJs are suitable for the generation of entangled superposition Fock states, which may be resistant to some weak losses.
Remarkably, quantum solitons represent the Lieb-Liniger model, which was discussed in the literature for a long time ago, see e.g. [45][46][47]. Various aspects of solitons quantum fluctuations related to quadrature squeezing, entanglement, and photon number statistics were studied in a quantum optics domain [48][49][50][51]. However, optical solitons represent systems containing a large number of particles far from thermal equilibrium. In this sense, atomic solitons favorably differ from optical analogues and could be recognized as the best candidates for the experimental realization of macroscopic quantum superposition states.
The progress in theory and experiment achieved with matter-wave bright solitons is discussed and summarized in [52]. The bright solitons were experimentally demonstrated in several labs with atomic BECs possessing negative scattering length and allowing a moderate number of particles (from several tens to thousands) [53,54]. The collisions for such solitons were recently shown in [55]. Notably, the observation of bright matter-wave solitons, which admit collapsing phenomena due to attractive inter-particle interaction, represents a non-trivial task for experimenters, see e.g. Lev Khaykovich's article in Ref. [52]. Alternatively, bight solitons may be formed in periodic structures by manipulating the sign of effective particle mass [56]. However, exploiting a macroscopically large particle number (thousands and more) in solitons leads to various (manybody) losses [57]. In this sense, solitons with mesoscopic number of particles (few hundreds and less) represent a primary interest. Some of recent proposals with quantum atomic solitons are established in [58][59][60].
Finally, let us mention mesoscopic atomic Bose-systems, so-called quantum droplets, which occur by means of manipulation with scattering atomic length, and posess intriguing quantum features [61]. Ferrier-Barbut et al. in [62] demonstrated the important influence of quantum fluctuations on droplet stabilization beyond the mean-field level that occur for strong dipolar Bose gases containing several hundreds of particles. Physically, beyond-mean-field energy correction in low dimensions plays an essential role for quantum droplet formation and behavior [63]. In particular, as predicted in [64,65], quantum droplets represent soliton-like (bell-shape) objects and consist of two one-dimensional Bose-Bose mixtures with competing (repulsive and attractive) interactions.
In this work we establish a new approach to the two-soliton (SJJ) problem mapping it onto the effective two-mode Hamiltonian and performing a second quantization procedure.
The paper is arranged as follows. In Sec. 2 we discuss a classical field approach to coupled BECs possessing negative scattering length. We represent a general approach for tunneling BECs in elongated traps in Sec. 2.1. Important peculiarities and limitations for the JJ-models are also discussed. In Sec. 2.2 we represent a familiar BJJ-model problem valid for atomic condensates with a Gaussian wave function. We implement this model through the paper for validation of the results obtained with the SJJ-model. Then in Sec. 2.3 we establish the SJJ-model that operates with coupled 1D bright solitons. The Hamilton formalism is used to find mean-field equations.
The important features of these equations are revealed in a steady state. The quantum approach to the BJJ-and SJJ-models is given in Section 3. The macroscopic properties of the energy spectrum for superposition states are analyzed in the framework of the Hartree approximation. In Section 4 the full quantum theory of the SJJ-model is presented. The peculiarities of the quantum energy spectrum are discussed in detail. We elucidate the quantum features of entanglement and EPR steering in the crossover region where the Hartree approach is violated. The experimental conditions for quantum superposition state formation in the presence of many-body losses are discussed in Sec. 4.3. We also analyze the influence of few particle losses on predicted condensate states in the framework of the fictitious beam splitter approach, cf. [66]. Finally, in Conclusion we summarize the results obtained.

II. CLASSICAL FIELD MODELS FOR COUPLED ATOMIC CONDENSATES
A. General approach to JJ-models 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 Fig. 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 Fig. 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 Fig. 1.
We represent the condensate wave function Ψ(r ⊥ , x, t) relevant to the scheme in Fig. 1 in a superposition state as where Φ 1,2 (r ⊥ + d/2) and Φ 2 (r ⊥ − d/2) are spatial (time independent) wave functions of condensates in transverse plane on either side of the barrier.
In the framework of field theory it is possible to consider a Hamilton function corresponding to Eqs. (2), which looks like (cf. [42,43]) where V H (x) = 1 2 ν 2 x 2 is trapping potential in X direction. In (2), (3), Ψ 1 and Ψ 2 obey a normalization condition where N is the average total number of particles.
In Eqs. (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][43][44]. In other words, all lengths and time variables in the work are given in a ⊥ and ω −1 ⊥ units respectively. Let us briefly discuss limitations for Eqs. (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, Eq. (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 ; Refs. [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 7 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 the number of atoms 5.2 × 10 3 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 ef f < 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 where x 0 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 ef f |, in the lattice, cf. [73]. The ratio m/|m ef f | 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.

B. Classical BJJ-model
Let us assume that the wave function shape in X-dimension is Gaussian and looks like where N j and θ j are the average particle number and phase of j-th 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 It is instructive to represent the Hamiltonian (7) in terms of new variables z = (N 2 − N 1 )/N and θ = θ 2 − θ 1 . Omitting the constant energy term for (7) we obtain where z and θ are normalized population imbalance and phase difference, respectively. In Eq. (8) we introduce the effective nonlinear parameter The Eq. (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 [70], physically correct tunneling rate K must be negative, cf. [42][43][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 Fig. 1. Thereafter, we completely neglect the trapping term in Eq. (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 Substituting (10) where The Hamiltonian (11) in terms of z and θ variables reads as and looks similar to Eq. (8) obtained for the BJJ-model. In (13) we introduced the effective nonlinear parameter is the vital parameter that determines various regimes of soliton dynamics.
Eq. (13) derived for the SJJ-model can be understood in the framework of the BJJ-model Hamil-

that depends on solitons population
imbalance, z. The nonlinear effects in the tunneling process vanish for z 2 << 1. In this limit the BJJ-and SJJ-models possess the same features.
However, when z 2 → 1, the effective tunneling rate, κ ef f , goes to zero and effective param- eter Λ ef f essentially increases. In other words, we deal here with self-tuning effect for effective tunneling parameter, κ ef f , 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.
These equations permit non-trivial steady-state solutions which are valid for 1.58 < Λ ≤ 2.42. For another specific value of phase θ we have 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.

A. Effective Hamiltonian quantization
The quantum theory approach proposes a second quantization procedure for the classical effective Hamiltonians (8), (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 BJJmodel. The quantization procedure that we propose here is based on the approach described in [76]. We construct a quantum Hamiltonian from Eq. (8) by introducing particle number op-eratorsN 1 =â †â andN 2 =b †b , respectively; one can define the reduced particle number dif- are usual Bosonic annihilation (creation) operators for condensate modes in two-mode approximation, cf. [76]. Then one can 1 +ẑe −iθ/2 for annihilation operatorsâ andb, respectively.
The quantized Hamiltonian obtained from (8) for the BJJ-model now has the form For the SJJ-model we start from the classical Hamiltonian (13) rewriting it as that is suitable for quantization procedure described above. The last term in (19) we treat by using (formal) Taylor expansion representation for the operator

Finally, the quantum version of Hamiltonian (19) reads aŝ
Thereafter, we refer to Eq. (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 Eq. (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 where we introduced parameter ℘ ≡ πκ/2ν that characterizes differences in the BJJ-and SJJ-models experimental realization.

B. Macroscopic superposition states in Hartree approximation
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 |ψ N = |ψ ⊗ |ψ ⊗ ... |ψ , where |ψ 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 where |0 ≡ |0 a |0 b is a two-mode vacuum state and α, β are unknown (variational) wave functions obeying normalization condition |α| 2 + |β| 2 = 1. The mean energy in respect of state (21) is obtained from (20) and looks like The energy E in Eq. (22) implies Lagrangian, L, given by By using Eq. (23) we derive the equations of motion for the field variables α, β in the form, where it is assumed that |α| = 0, or |β| = 0 that implies S = ±1.
Eqs. (25) lead to a single equation: In general, the functions α and β in Eq. (26) can be represented as for θ = π, and to for θ = 0, respectively. Notice, if we change the sign of tunneling rate, κ, the shift of phase θ by π is expected for (27), (28).
Eqs. (27) The energies corresponding to Eqs. (29) are given in the form It is important that E ± ≥ E 0 for any Λ, when S ± solution exists.
Substituting (29b) into (21) one can obtain "two halves" of the SC-state The states (31) can form macroscopic superpositions possessing the same energy (30b).
In other words, the N 00N -state appears for S = ±1 that implies |α| = 0 or |β| = 0 in (21)-(23) (see also Eqs. (17)). These solutions are where "±" correspond to the sign of S that defines two "halves" of the N 00N -state: In (35) The phase θ in (35) obey Eq. (17b) and defines the domain 0 < Λ ≤ 1.58 for Λ-parameter, where the N 00N -state can appear in the framework of the Hartree approach discussed above.
Noteworthy, the state (35) is fragile enough and rapidly collapses to states |N − 1 a |0 b or |0 a |N − 1 b even in the presence of one particle loss, cf. [38]. As we see below, the SJJ-model exhibits some new specifics for the N 00N -state formation in a purely quantum limit.

IV. MESOSCOPIC QUANTUM SUPERPOSITION STATES
A. Quantum energy spectrum 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 tunnelling 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 Fig. 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), (20) has to be fully quantum.
We represent the quantum two-mode state, |Ψ(τ ) , in the Fock basis as (cf. [16]): where |N −n, n ≡ |N −n a |n b denotes the two-mode atom-number state; the coefficients A n (τ ) obey normalization condition N n=0 |A n (τ )| 2 = 1. The A n (τ ) fulfills the Schrödinger equation In connection with quantum Hamiltonians (18), (20) the Eq. (38) can be represented as: where the following notations are introduced: for the BJJ-model, and for the SJJ one.
Then let us focus on stationary solution A n (τ ) = A n e −iEnτ of (39) that enables to find out energy spectrum E n of the quantum Hamiltonian (20) for the SJJ-model. Main peculiarities with the coupled soliton energies occur in the vicinity of point Λ ≡ Λ c ≈ 2.0009925. 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 selftrapping, cf. [8,42,44]. In particular, at Λ < Λ c energetically favorable state for the quantum SJJsystem in Fig. 2(a) is described by the blue dashed (or the red solid) line exhibiting the minimal energy E 0 . As seen from Fig. 2(a), at Λ ≥ Λ c the energies E ± , E (N 00N ) , which correspond to the superposition states in the Hartree approximation, are comparable with E 0 but placed a little bit above it.
Situation changes completely in a purely quantum limit when parameter Λ crosses critical its value Λ c . The red solid line in Fig. 2(a) clearly shows that the quantum N 00N -state is energetically favorable at Λ ≥ Λ c for the SJJ ground state (it is demonstrated below that this state is not exactly a N 00N -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: Similar results may be obtained for the BJJ-model. Diagonalization of the relevant Hamiltonian (18) leads to The energy eigenvalues (42) and (43) are plotted in Fig 2(b) (44) we obtain for the energy gaps ∆E BJJ ≈ 4κ 2 /u, and ∆E SJJ ≈ 5κ 2 /N u 2 ≈ 5∆E BJJ /4N u, respectively. The suppression of this gap with increasing atom-atom scattering length (u-parameter) is clearly seen in Fig 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 Fig. 3 we represent the ground state probability distributions, |A n | 2 , (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. Fig. 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.
The differences between distributions considered for the SJJ-and BJJ-models become evident and important in Fig. 3(b)-(d) due to the particle number dependence of the effective tunneling rate, κ ef f , see (14). Figures 3(b) and 3(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 Fig. 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 | 2 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 Fig. 3(c,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 Fig. 3

(d). The
N 00N -state, occurs for the BJJ-model in the limit of a very large λ-parameter, see below Fig. 4(b) and cf. [16,18].
In general, the SJJ-model demonstrates evident advantages to achieve a N 00N -state. Roughly speaking, state (37) for the SJJ-model after the crossover represents a superposition of entangled Fock states possessing a large N 00N -state component -see. Fig. 3(c). A physical explanation for this phenomenon is illustrated as follows.
First, for a given tunneling rate, κ, the Λ-parameter depends on N 2 , which allows to achieve a required crossover point for the moderate number of particles, N .
Second, as it follows from Fig. 3(b,c) the crossover for the SJJ-model happens abruptly with a slightly increasing Λ parameter. The dependence of the effective tunneling rate, κ ef f , on the atom number leads to a significant improvement of the "edge" states occupation. In particular, as clearly seen in Fig. 3(c), two large peaks of probability distribution occur at the "edges" with n = 0 and n = N , respectively. They correspond to the N 00N -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 Fig. 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" N 00N -state (45) at moderate values of Λ, cf. Fig. 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].

B. Entanglement, steering and planar spin-squeezing
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 prerequisits 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 E (m) The entanglement occurs when inequalities are satisfied. The value E (m) HZ = 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 Eq. (46) in terms of N -particle spin-operators as following: where δJ j ≡ is the variance of fluctuations for j-th component of spin, established by operatorŝ Combining Eq. (46) with (47) Fig. 4(b)) in this limit we can exploit the two-mode condensate state like (21) with The EPR steering entanglement is achieved in the region, where E (1) HZ < 0.5, cf. [21,25]. From Fig. 4(a) it is evident that this area occurs close to the point Λ c 2. Noteworthy, the sharp behavior of E HZ within the crossover region appear in the absence of losses, see the inset in Fig. 4(a). Evidently, broadening of this region is expected if one-and, especially, three-body losses are taken into account.
The minimal value E (1) HZ,min of parameter E (1) HZ , 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 In Fig. 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 .
The coefficient C J theoretically can be estimated by using the Gaussian distribution function possessing width, σ, (see Fig. 3(b)), cf. [17]. to the value Λ c 2 (or, λ c 1). However, as clearly seen from Fig. 3(b) the distribution for the SJJ-model is sufficiently narrower than for the BJJ one, i.e. σ SJJ < σ BJJ , and we expect to take place for curves exhibiting E HZ in Fig. 4(a,b). Notably, the suppression of variance δJ manifests planar squeezing defined as (cf. [32]) where is the magnitude of the spin component in XY -plane. In particular, for both SJJ-and BJJ-models we have Ĵ X ≈ J, Ĵ Y = Ĵ Z = 0 within the domain 0 < Λ ≤ 2 for Fig. 4(a) and 0 ≤ λ ≤ 1 for Fig. 4(b), respectively. Thus, in this area E (1) HZ 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 Fig. 4(a) (green dashed line) the SJJ-system undergoes the abrupt transition to N 00N -state (45) at the crossover region for Λ c 2 that is consistent with Fig. 3(c).
The E

C. Quantum state engineering in the presence of losses
We now outline the possibility of experimental observation of N 00N -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 32ms 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 Fig. 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∇ · 32ms. The Λ parameter achieves its critical value, Λ c 2, at tunneling rate |K|/2π 77 Hz or |K| 3.7 nK, 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 K or scattering length |a sc | allows to obtain the crossover to the N 00N -state with Λ ≥ Λ c .
For the BJJ-model the critical value λ c = 1 is achieved with |K|/2π 83 Hz. However, the N 00N -state requires essential enhancement of λ in the experiment (60 times and more, see inset in Fig. 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 Fig. 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 quan-tum metrology, where N 00N -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 4 atoms. However, the N 00N -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 N 00N -state formation. Since vital parameter, Λ, is proportional to N 2 , 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, κ ef f , 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 ∝ a 4 sc /m that approaches to L 3 2.6 × 10 −28 cm 6 /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 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 13 cm −3 for the effective rate of three-body losses γ 3 = 2L 3 ρ 2 we obtain γ 3 = 5.2 × 10 −2 s −1 .
The one-body losses rate may be estimated as γ 1 = 0.2 s −1 , see e.g. [41]. At times of soliton observation in [54,55] we effectively obtain γ 1,3 t << 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 N 00N -state applicability for quantum metrology tasks, see e.g. [38][39][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) N 00N -state may be non-optimal. It is shown that entangled Fock states, unbalanced N 00N -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 beam splitter (BS) approach. This approach represents a useful tool to study coupling of quantum macroscopic superposition states with environment -see Fig.6. The "input" two-mode Fock state (37) after two BSs transforms into (cf. [40]) where l a and l b are the numbers of particles lost from "a" and "b" wells, respectively. In (53) we introduce a coefficient 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 |Ψ out after particle losses reads as (see Fig.6) is a normalization constant for state |Ψ out . 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 |Ψ In Fig. 7 we represent the probabilities |C n 1,0 | 2 ( C n la,l b ≡ 1 √ p A n B n la,l b ) for the SJJ-and BJJmodels versus n in the presence of weak and equal losses in each channel, η a = η b = 0.999. To be more specific, in Fig.7(a,b), we use the same condensate parameters as in Fig.3(c,d), respectively.

The multiplier
√ N − n in (56) plays a crucial role in the behavior of relevant probability |C n 1,0 | 2 that characterizes the conditional state |Ψ (1,0) out . In the presence of weak losses, the tendency to occupation of "edge" states continues with increasing of parameter Λ, see Fig.7(a). However, in this case the state with n = N cannot be occupied due to one particle loss, cf. Fig. 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 |C n . On the contrary, the probability of occupation of new "edge" state with n = N − 1 is more than √ N times lower.
Physically, this case may by recognized as so-called partial "collapse" of N 00N -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 Fig. 7(b). Notice that the states with n = N − 1, N − 2, ... are poorly populated even without any losses -see insets to Fig. 3(c) and Fig. 3(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 Fig. 7(d) and cf. Fig. 3(d).
However, for l a = l b = 1, i.e., if two particles are detected simultaneously after BSs in Fig. 6, the behavior of the conditional state (55) becomes fundamentally different. In this case for (55) we obtain: 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 N 00N -state that possesses N − 2 total number of particles and has a form The unbalanced N 00N -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 Fig. 6. Tracing out the ancillary modes l a , l b in (55) for the resulting state we obtain (cf. [38]) In Fig. 8, we plot the probabilities |C n la,l b | 2 in N − n − l a and n − l b axis for the SJJ-model. Each column in Fig. 8 corresponds to the probability of state |N − n − l a a |n − l b b occupation. The parameter Λ = 4 in Fig. 8 is the same as for Fig. 3(d). Two big (yellow) columns in Fig. 8 correspond to the N 00N -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 Fig. 8. In particular, the columns depicted in Orange and Green-colors in Fig. 8  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 Fig. 6. However, relevant probabilities essentially vanish with increasing l a and l b .

V. CONCLUSION
In this work we have considered the problem of macroscopic and mesoscopic states formation with two weakly-coupled atomic condensate bright solitons. The condensates are trapped in a double-well potential in the Y Z-plane and are significantly elongated in the X-axis forming the specific Josephson junctions (JJs) geometry. We assume that condensates admit a weak interaction between the atoms possessing a negative scattering length of atom-atom interaction. This geometry allows to examine bright atomic soliton features in one (X) dimension. We start from the classical (mean) field approach to the problem that admits bright soliton solutions with preserving soliton shapes. In this limit these solutions permit to effectively map the system of weakly-coupled solitons onto the two-mode problem.
Transferring to the quantum picture of the problem we perform the second quantization proce- Comparing all the results obtained for the SJJs with the ones of Bosonic Josephson junction (BJJ) model consisting of two weakly-coupled (Gaussian-shape) condensates, we have shown that there exist two important differences between the SJJ-and BJJ-models. First, the SJJ-model permits a particle dependent tunneling rate. Second, the effective nonlinear parameter for the SJJmodel depends on the total particle number squared, N 2 , not N , that takes place for the BJJ-model.
To demonstrate the advantages of the SJJ-device we have examined the quantum analysis of the superposition states formation.
Initially, we have considered the quantum field variational (Hartree) approach that predicts Schrödinger-cat states formation within the vital parameter, Λ, domain. In this limit we explore the atomic coherent state as the basis for the superposition state analysis. We have shown that the cat states approach the N 00N -state in the limit when the cat becomes macroscopic.
We have performed the full quantum analysis of the energy spectrum for the SJJ-and BJJmodels to consider the N 00N -state formation problem for the JJs in detail. We have demonstrated that the crossover to the N 00N -state may be realized for the both models at some combination of atomic nonlinear strength and tunneling rate. However, due to the particle dependence, the effective tunneling rate, κ ef f , that governs this crossover, demonstrates self-tuning effect, revealing an abrupt transition to superposition of entangled Fock states. The N 00N -state appears at the "edges" with n = 0, N of this superposition and represents its major component. The existence of the satellite entangled Fock state components with non-zero probabilities in the vicinity of "edges" represent an important difference from the variational (Hartree) approach. It is worth noting that probabilities of these components are vanishing, and the SJJ-model approaches approximately "ideal" N 00N -state (45) for moderate values of Λ-parameter. The BJJ-model approaches the same N 00N -state only with a quite large nonlinear strength, λ ≥ 60.
To confirm the results obtained we have examined the m-order Hillery-Zubairy criterion and taken it with m = 1 and m = N , respectively. The numerical calculations reflect the existence of the strong entanglement and the planar quantum squeezing, which occur up to the crossing point defined by the vital parameter Λ (or λ). The EPR steering happens in the crossover point vicinity where Hartree approach is useless. Notably, the depth of the steering is much higher for the familiar atomic JJs, which may be explained through the additional particle number fluctuations occurring from the tunneling process for the SJJ-model. This process leads to a narrower distribution function appearing immediately at the crossover point. The results obtained for the HZ-criterion indicate the crossover from the quantum steering to the N 00N -state. This crossover exhibits a very sharp and large variation in the first order HZ-criterion that takes place for the SJJ-model. Finally, we examine how losses influence the formation of state (37). The qualitative estimations of influence from one-and three-body losses demonstrate feasibility of observation on predicted phenomena with SJJs. To take into account the losses of a few number of particles in condensate solitons, we have used the fictitious beam splitter approach that represents the essence of the state preparation procedure in Fig. 6. In particular, we have shown that the existence of tiny satellite entangled Fock state components provides some resistance of the N 00N -state to complete collapsing if we are able to detect (non-equal) number particles which are lost in each well (channel). We have also shown that the final state (56) preserves the N 00N -state structure if two or even number of particles are lost and may be equally detected in both channels. Especially, such an analysis seems to be important for practical applications of the soliton JJs in "real world" quantum metrology.
We are confident that the quantum state features of the SJJ-model described by Hamiltonian in Eq. (20) potentially may be significantly larger than we describe in this work. In particular, established quantum SJJ-model may be also useful in quantum communication and quantum information science, where quantum optical solitons represent an appropriate tool, cf. [97]. In this sense quantum properties of the SJJ-model proposed in this work represent a necessary step to further studies.