# Bose-Einstein condensate soliton qubit states for metrological applications

The Vinh Ngo  1 Dmitriy Tsarev  1 Ray-Kuang Lee  2,3,4 and Alexander
Alodjants1  1National Research University for Information Technology,
Mechanics and Optics (ITMO), St. Petersburg 197101, Russia
2Physics Division, National Center for Theoretical Sciences, Hsinchu 30013,
Taiwan
3Institute of Photonics Technologies, National Tsing Hua University, Hsinchu
30013, Taiwan
4Center for Center for Quantum Technology, Hsinchu 30013, Taiwan
<EMAIL_ADDRESS>
###### Abstract

By utilizing Bose-Einstein condensate solitons, optically manipulated and
trapped in a double-well potential, coupled through nonlinear Josephson
effect, we propose novel quantum metrology applications with two soliton qubit
states. In addition to steady-state solutions in different scenarios, phase
space analysis, in terms of population imbalance - phase difference variables,
is also performed to demonstrate macroscopic quantum self-trapping regimes.
Schrödinger-cat states, maximally path-entangled ($N00N$) states, and
macroscopic soliton qubits are predicted and exploited for the
distinguishability of obtained macroscopic states in the framework of binary
(non-orthogonal) state discrimination problem. For arbitrary phase estimation
in the framework of linear quantum metrology approach, these macroscopic
soliton states are revealed to have a scaling up to the Heisenberg limit (HL).
The examples are illustrated for HL estimation of angular frequency between
the ground and first excited macroscopic states of the condensate, which opens
new perspectives for current frequency standards technologies.

††journal: oe††articletype: Research Article

##  1 Introduction

Nowadays, nonlinear collective mode formation and interaction in Kerr-like
medium represent an indispensable platform for various practically important
applications in time and frequency metrology [1, 2], spectroscopy [3, 4],
absolute frequency synthesis [5], distance ranging [6]. In photonic settings
frequency combs are proposed for these purposes [7]. The combs occur due to
the nonlinear mode mixing in special (ring) microcavities, which possess some
certain eigenmodes. Notably, bright soliton formation emerges with vital
phenomena accompanying micro-comb generation [8]. Physically, such a soliton
arises due to the purely nonlinear effect of temporal self-organization
pattern occurring in an open (driven-dissipative) photonic system. However,
because of the high level of various noises in the system they can be hardly
explored for purely quantum metrological purposes.

Instead, atomic optics, which operates with Bose-Einstein condensates (BECs)
at low temperatures, provides a suitable platform for various quantum devices
that may be useful for metrology and sensing tasks [9]. In particular, so-
called Bosonic Josephson junction (BJJ) systems, established through two
weakly linked and trapped atomic condensates, are at the heart of the current
quantum technologies in atomtronics that considers atom condensates and aims
to design (on-chip) quantum devices. The condensates in this case represent
low dimensional systems and may be manipulated by magnetic and laser field
combinations. In this sense, they represent an alternative to optical
analogues.

The BJJs are intensively discussed and examined both in theory and experiment
[12, 11, 10, 13, 14]. The quantum properties of the BJJs are also widely
studied [15, 17, 21, 16, 18, 22, 19, 20], demonstrating spin-squeezing and
entanglement phenomena [23, 24, 25], with the capability in generating
$N00N$-states [26, 27] to go beyond the standard quantum limit [28].
Physically, the BJJs possess interesting features connected with the interplay
between quantum tunneling of the atoms and their nonlinear properties evoked
by atom-atom interaction [29, 30].

With the Kerr-like nonlinearities, solitons naturally emerge from atomic
condensates in low dimensions, see e.g. [31, 32, 33, 34, 35]. Especially, the
bright atomic solitons observed in lithium condensate possessing negative
scattering length [33, 34, 35] are worth noticing. Atomic gap solitons are
also observed in condensates with repulsive inter-particle interaction [32].
Based on soliton modes, we recently proposed the quantum soliton Josephson
junction (SJJ) device with the novel concept to improve the quantum properties
of effectively coupled two-mode system [37, 38, 39]. The SJJ-device consists
of two weakly-coupled condensates trapped in a double-well potential and
elongated in one dimension, cf. [36]. We demonstrated that quantum solitons
may be explored for the improvement of phase measurement and estimation up to
the Heisenberg limit (HL) and beyond [37]. In the framework of nonlinear
quantum metrology, we also showed that solitons permits Super-Heisenberg (SH)
scaling ($\propto N^{-3/2}$) even with coherent probes [38]. On the other
hand, steady-states of coupled solitons can be useful for effective formation
of Schrödinger-cat (SC) superposition state and maximally path-entangled
$N00N$-states, which can be applied for the phase estimation purposes[39]. It
is important that such superposition states arise only for soliton-shape
condensate wave functions and occur due to the existence of certain steady-
states in phase difference - population imbalance phase plane, cf. [38].

Remarkably, macroscopic states, like SC-states, play an essential role for
current information and metrology, see e.g. [40]. In quantum optics, various
strategies are proposed for creation photonic SC-states and relevant
(continuous variable) macroscopic qubits [41, 42, 43]. Special (projective)
measurement and detection techniques are also important here [44, 45, 46].
From this point of view, the fact that condensate environment, dealing with
mater waves, provides the minimally accessible thermal noises makes it
potentially promising for macroscopic qubits implementation, cf. [47, 48, 49].

In this work, we propose metrological applications for two soliton
superposition states as macroscopic qubits. The interaction between these
solitons comes from the nonlinear mode mixing in an atomic condensate trapped
in a double-well potential. In particular, we reveal the SC-states formation
and their implementation for arbitrary phase measurement prior HL and beyond.
Since SC-states are non-orthogonal states, a special measurement procedure is
applied by so-called sigma operators, as it enables us to estimate the unknown
phase parameter [50]. On the other hand, our approach can be also useful in
the framework of discrimination of binary coherent (non-orthogonal) states in
quantum information and communication, see e.g. [51, 52]. The non-
orthogonality of these states leads to so-called Helstrom bound for the
quantum error probability that simply indicates the impossibility for a
receiver to identify the transmitted state without some errors [53, 54]. In
quantum metrology, by means of various regimes of condensate soliton
interaction, we deal with a set of states for a quantum system, which may be
prepared before the measurement. Our results show that these SC-states
approach the soliton $N00N$-states to minimize the quantum error probability.

##  2 Two-soliton model

###  2.1 Coupled-mode theory approach

We start with the mean-field description of coupled mode-theory approach to
elongated (in $X$-dimension) BEC trapped in $V=V_{H}+V(x)$ potential, where
$V_{H}$ is a familiar 3D harmonic trapping potential, $V(x)$ is responsible
for condensate double-well confinement in one ($X$) dimension [39]. The
condensate (rescaled) wave function (mean field amplitude) $\Psi(x)$ obeys the
familiar 1D Gross-Pitaevskii Equation (GPE), cf. [36]:

$i\frac{\partial}{\partial t}\Psi=-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\Psi-u\left|\Psi\right|^{2}\Psi+V(x)\Psi,$ (1)

where $u=4\pi|a_{sc}|/a_{\perp}$ characterizes Kerr-like (focusing)
nonlinearity, $a_{sc}<0$ is the s-wave scattering length that appears due to
atom-atom scattering in Born-approximation,
$a_{\perp}=\sqrt{\hbar/m\omega_{\perp}}$ characterizes the trap scale, and $m$
is the particle mass. To be more specific, we only consider condensates
possessing a negative scattering length. In Eqs. (1) we also propose rescaled
(dimension-less) spatial and time variables, which are $x,y,z\rightarrow
x/a_{\perp},y/a_{\perp},z/a_{\perp}$, and $t\rightarrow\omega_{\perp}t$, cf.
[37, 38, 36].

The nonlinear coupled-mode theory admits solution of (1) that simply
represents a quantum-mechanical superposition

$\Psi(x,t)=\Psi_{1}(x,t)+\Psi_{2}(x,t),$ (2)

where wave functions $\Psi_{1}(x)$ and $\Psi_{2}(x)$ characterize the
condensate in two wells. For weakly interacting atoms one can assume that

$\Psi_{1,2}(x,t)=C_{1,2}(t)\Phi_{1,2}(x)e^{i\beta_{1,2}t},$ (3)

where $\Phi_{1}(x)$ and $\Phi_{2}(x)$ are ground- and first-order excited mode
state wave functions possessing energies $\beta_{1}$ and $\beta_{2}$,
respectively; $C_{1}(t)$ and $C_{2}(t)$ are time-dependent functions. If the
particle number is not too large, GPE (1) may be integrated in spatial ($X$ \-
dimension) leaving only two condensate variables $C_{1,2}(t)$, cf. [30]. In
particular, $\Phi_{1}(x)$ and $\Phi_{2}(x)$ may be time-independent Gaussian-
shape wave functions obeying different symmetry. Practically, this two-mode
approximation is valid for the condensates of several hundreds of particles
[55]. The condensate in this limit is effectively described by two
(macroscopically populated) modes as a result.

###  2.2 Quantization of coupled solitons

The sketch in Fig. 1 explains the two-soliton system described in our work. If
trapping potential $V(x)$ is weak enough and the condensate particles interact
not so weakly, the ansatz solution (3) is no more suitable. Especially we
would like to mention condensates with a negative scattering length which
admit a bright soliton solution for $\Psi_{1,2}(x,t)$ in (2). In fact, in this
case one can speak about two-soliton solution problem for GPE (1) without
trapping potential $V(x)$ known in classical theory of solitons [56].

In quantum theory we deal with a bosonic field operator
$\hat{a}(x,t)\propto\hat{a}_{1}+\hat{a}_{2}$ instead of (2), where
$\hat{a}_{1,2}\equiv\hat{a}_{{1,2}}(x,t)$ are field operators which correspond
to mean-field amplitudes $\Psi_{1,2}(x,t)$ defined in (2). We assume that
experimental conditions allow the formation of atomic bright solitons in each
well of harmonic double-well potential. In particular, these conditions may be
realized be means of manipulation with weak trapping potential $V(x)$. In an
experiment this manipulation may be performed with the help of a dipole trap
and laser field.

Figure 1: Sketch of probability density distribution $|\Psi|^{2}$ versus
spatial coordinates $X$ and $Y$, 2D projection of the 3D coupled condensates
trapped in a double-well (dashed green curve) and harmonic (dashed magenda
curve) potentials, respectively. Shadow regions display 1D condensate wave
packets projections; they represent a secant-shape in $X$-direction, and
Gaussian-shape in the transverse directions.

Then, considering linear superposition state like (4) we can write the total
Hamiltonian $\hat{H}$ for two BEC solitons in the form

$\hat{H}=\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{int},$ (4a) where $\hat{H}_{j}$
($j=1,2$) is the Hamiltonian for condensate particles in $j$-th well; while
$\hat{H}_{int}$ accounts the coupling between two wells due to the soliton
overlapping. In the second quantization form we explicitly have
$\displaystyle\hat{H}_{j}=\int_{-\infty}^{\infty}\hat{a}_{j}^{\dagger}\left(-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\right)\hat{a}_{j}dx;$ (4b)
$\displaystyle\hat{H}_{int}=-\frac{u}{2}\int_{-\infty}^{\infty}\Big{(}\hat{a}_{1}^{\dagger}+\hat{a}_{2}^{\dagger}\Big{)}^{2}\Big{(}\hat{a}_{1}+\hat{a}_{2}\Big{)}^{2}dx.$
(4c)

The annihilation (creation) operators of bosonic fields denoted as
$\hat{a}_{j}$ (${\hat{a}^{\dagger}_{j}}$) with $j=1,2$ obey the commutation
relations:

$[\hat{a}_{i}(x),{\hat{a}^{\dagger}_{j}(x^{\prime})}]=\delta(x-x^{\prime})\,\delta_{ij};\quad
i,j=1,2.$ (5)

In the Hartree approximation for a large particle number, $N>>1$, one can
assume that the quantum $N$-particle two-soliton state is the product of $N$
two-soliton states and can be written as [57, 58, 59]

$\left|\Psi_{N}\right\rangle=\frac{1}{\sqrt{N!}}\left[\int_{-\infty}^{\infty}\left(\Psi_{1}(x,t)\hat{a}_{1}^{\dagger}e^{-i\beta_{1}t}+\Psi_{2}(x,t)\hat{a}_{2}^{\dagger}e^{-i\beta_{2}t}\right)dx\right]^{N}\left|0\right\rangle,$
(6)

where $\Psi_{j}(x,t)$ are unknown wave functions,
$|0\rangle\equiv|0\rangle_{1}|0\rangle_{2}$ is a two-mode vacuum state. The
state (6) is normalized as
$\left\langle\Psi_{N}\big{|}\Psi_{N}\right\rangle=1$, and the bosonic field-
operators $\hat{a}_{j}$ act on it as

$\hat{a}_{j}\left|\Psi_{N}\right\rangle=\sqrt{N}\Psi_{j}(x,t)e^{-i\beta_{j}t}\left|\Psi_{N-1}\right\rangle.$
(7)

Applying variational field theory approach based on the ansatz $\Psi_{j}(x,t)$
we obtain the Lagrangian density in the form:

$L_{0}=\frac{1}{2}\sum_{j=1}^{2}\left(i\left[\Psi_{j}^{*}\dot{\Psi}_{j}-\dot{\Psi}_{j}^{*}\Psi_{j}\right]-\left|\frac{\partial\Psi_{j}}{\partial
x}\right|^{2}\right)+\frac{uN}{2}\left(\Psi_{1}^{*}e^{i\beta_{1}t}+\Psi_{2}^{*}e^{i\beta_{2}t}\right)^{2}\left(\Psi_{1}e^{-i\beta_{1}t}+\Psi_{2}e^{-2i\beta_{2}t}\right)^{2},$
(8)

where we suppose $N-1\approx N$ and omit the common term $N$.

Noteworthy, from (8) one can obtain the coupled GPEs for $\Psi_{j}$-functions
as

$\displaystyle i\frac{\partial}{\partial t}\Psi_{1}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\Psi_{1}-uN\left(\left|\Psi_{1}\right|^{2}+2\left|\Psi_{2}\right|^{2}\right)\Psi_{1}$
$\displaystyle-$ $\displaystyle
uN\left(\left|\Psi_{2}\right|^{2}+2\left|\Psi_{1}\right|^{2}\right)\Psi_{2}e^{-i\beta
t}-uN\Psi_{1}^{*}\Psi_{2}^{2}e^{-2i\beta
t}-uN\Psi_{2}^{*}\Psi_{1}^{2}e^{i\beta t};$ $\displaystyle
i\frac{\partial}{\partial t}\Psi_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\Psi_{2}-uN\left(\left|\Psi_{2}\right|^{2}+2\left|\Psi_{1}\right|^{2}\right)\Psi_{2}$
$\displaystyle-$ $\displaystyle
uN\left(\left|\Psi_{1}\right|^{2}+2\left|\Psi_{2}\right|^{2}\right)\Psi_{1}e^{i\beta
t}-uN\Psi_{2}^{*}\Psi_{1}^{2}e^{2i\beta
t}-uN\Psi_{1}^{*}\Psi_{2}^{2}e^{-2i\beta t},$

where $\beta=\beta_{2}-\beta_{1}$ is the energy (frequency) spacing.

Set of Eqs. (9) lead to the known problem for transitions between two lowest
self-trapped states of condensates in the nonlinear coupled mode approach if
we account (3) for $\Psi_{j}(x,t)$ condensate wave functions representation
[29, 30].

On the other hand, Eqs. (9) can be recognized in the framework of soliton
interaction problem that may be solved by means of perturbation theory for
solitons [56]. In particular, in accordance with Karpman’s approach we can
recognize in (9) terms proportional to
$\epsilon_{jk}=\Psi_{j}^{*}\Psi_{k}^{2}+2|\Psi_{j}|^{2}\Psi_{k}$, $j,k=1,2$,
$j\neq k$ as perturbations for two fundamental bright soliton solutions.
Physically, $\epsilon_{jk}$ implies the nonlinear Josephson coupling between
the solitons.

In this work we establish a variational approach for solution of Eqs. (9), cf.
[38]. For the weakly coupled condensate states, i.e. for $\epsilon_{jk}\simeq
0$, set of Eqs. (9) reduces to two independent GPEs:

$i\frac{\partial}{\partial t}\Psi_{j}=-\frac{1}{2}\frac{\partial^{2}}{\partial
x^{2}}\Psi_{j}-uN\left|\Psi_{j}\right|^{2}\Psi_{j},$ (10)

which possess bright (non-moving) soliton solutions

$\Psi_{j}(x,t)=\frac{\sqrt{uN_{j}}}{2}\operatorname{\textrm{sech}}\left[\frac{uN_{j}}{2}x\right]e^{i\frac{u^{2}N_{j}^{2}}{8}t}.$
(11)

In the case of $\epsilon_{jk}\neq 0$ and for non-zero inter-soliton distance
$\delta$, we examine ansatzes for $\Psi_{j}(x,t)$ in the form

$\displaystyle\Psi_{1}(x,t)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN_{1}}}{2}\operatorname{\textrm{sech}}\left[\frac{uN_{1}}{2}x-\delta\right]e^{i\theta_{1}};$
(12a) $\displaystyle\Psi_{2}(x,t)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN_{2}}}{2}\operatorname{\textrm{sech}}\left[\frac{uN_{2}}{2}x+\delta\right]e^{i\theta_{2}}.$
(12b)

In particular, our approach presumes the existence of two well distinguished
solitons (separated by the small distance $\delta$, with the shape preserved)
interacting through dynamical variation of the particle numbers ($N_{j}\equiv
N_{j}(t)$) and phases ($\theta_{j}\equiv\theta_{j}(t)$), which occurs in the
presence of weak coupling between the solitons. In other words, $N_{j}$ and
$\theta_{j}$ should be considered as time-dependent (variational) parameters.
Setting (12) in (8) we obtain

$L=\int_{-\infty}^{\infty}L_{0}dx\simeq-z\dot{\theta}+\Lambda
z^{2}+\frac{\Lambda}{2}\left(1-z^{2}\right)^{2}I(z,\Delta)\left(\cos[2\Theta]+2\right)+\Lambda\left(1-z^{2}\right)J(z,\Delta)\cos[\Theta],$
(13)

where $z=(N_{2}-N_{1})/N$ ($N_{1,2}=\frac{N}{2}(1\mp z)$) is the particle
number population imbalance;
$\Theta=\theta_{2}-\theta_{1}-(\beta_{2}-\beta_{1})t\equiv\theta-\Omega t$ is
an effective time-dependent phase-shift between the solitons.

Physically, $\Omega$ is an angular frequency spacing between the ground and
first excited macroscopic states of the condensate; it represents a vital
(measured) parameter for metrological purposes in this work. In (13) we
introduce the notation $\Lambda=N^{2}u^{2}/16$ and define the functionals

$\displaystyle I\equiv I(z,\Delta)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\operatorname{\textrm{sech}}^{2}\left[\left(1-z\right)\left(x-\Delta\right)\right]\operatorname{\textrm{sech}}^{2}\left[\left(1+z\right)\left(x+\Delta\right)\right]dx;$
(14a) $\displaystyle J\equiv J(z,\Delta)$ $\displaystyle=$
$\displaystyle\sum_{s=\pm
1}\left(\int_{-\infty}^{\infty}\left(1+sz\right)^{2}\operatorname{\textrm{sech}}^{3}\left[\left(1+sz\right)\left(x+s\Delta\right)\right]\operatorname{\textrm{sech}}\left[\left(1-sz\right)\left(x-s\Delta\right)\right]\right),$

where $\Delta\equiv~{}\frac{Nu}{4}\delta$ is a normalized distance between
solitons.

Finally, by using Eq. (13) for the population imbalance and phase-shift
difference, $z$ and $\Theta$, we obtain the set of equations

$\displaystyle\dot{z}$ $\displaystyle=$
$\displaystyle\left(1-z^{2}\right)\left\\{\left(1-z^{2}\right)I\sin[2\Theta]+J\sin[\Theta]\right\\};$
(15a) $\displaystyle\dot{\Theta}$ $\displaystyle=$
$\displaystyle-\frac{\Omega}{\Lambda}+2z+\frac{\rm{d}}{\rm{d}z}\left\\{\frac{1}{2}\left(1-z^{2}\right)^{2}I\left(\cos[2\Theta]+2\right)+\left(1-z^{2}\right)J\cos[\Theta]\right\\},$
(15b)

where dots denote derivatives with respect to renormalized time $\tau=\Lambda
t$.

In contrast to the problem with coupled Gaussian-shape condensates (cf. [29,
30]), the solutions of Eqs. (15) crucially depend on features of governing
functionals $I(z,\Delta)$ and $J(z,\Delta)$. In Appendix we represent some
analytical approximations for $I(z,\Delta)$ and $J(z,\Delta)$, which we
exploit further.

##  3 Steady-state (SS) solutions

###  3.1 SS solution for $z^{2}=1$

The SS solutions of (15) play a crucial role for metrological purposes with
coupled solitons cf. [37]. We start from the SS solution $z^{2}=1$ of Eq.
(15a) by setting the time-derivatives to zero. As seen from (14), in the limit
of maximal population imbalance, $z^{2}=1$, $I$ and $J$ do not depend on
$\Delta$ and approach

$\displaystyle I(z,\Delta)$ $\displaystyle=$ $\displaystyle 1;$ (16a)
$\displaystyle J(z,\Delta)$ $\displaystyle=$ $\displaystyle\pi.$ (16b)

Substituting $z^{2}=1$ and (16) into (15b) we obtain

$\displaystyle z^{2}$ $\displaystyle=$ $\displaystyle 1;$ (17a)
$\displaystyle\Theta$ $\displaystyle=$
$\displaystyle\arccos\left[\frac{2\Lambda-\operatorname{\textrm{sign}}[z]\Omega}{2\pi\Lambda}\right].$
(17b)

Notably, in quantum domain the SS solutions (17) admit the existence of
quantum states with maximal population imbalance $z=\pm 1$ and phase
difference. The latter depends on the frequency spacing $\Omega$, which is the
subject of precise measurement with maximally path-entangled $N00N$-states in
this paper.

Below we perform the analysis of the SS solutions of Eqs. (15) in two limiting
cases $\Omega\neq 0$, $\Delta\simeq 0$ and $\Omega\simeq 0$, $\Delta\neq 0$.

###  3.2 SS solutions for $\Theta=0,\pi$ and $\Delta\simeq 0$

To find the SS solutions we rewrite (15b) as

$\frac{\Omega}{\Lambda}=2z-6z\left(1-z^{2}\right)I+\frac{3}{2}\left(1-z^{2}\right)^{2}\frac{\partial
I}{\partial z}-2zJ+\left(1-z^{2}\right)\frac{\partial J}{\partial z}$ (18)

for $\Theta=0$ and

$\frac{\Omega}{\Lambda}=2z-6z\left(1-z^{2}\right)I+\frac{3}{2}\left(1-z^{2}\right)^{2}\frac{\partial
I}{\partial z}+2zJ-\left(1-z^{2}\right)\frac{\partial J}{\partial z}$ (19)

for $\Theta=\pi$, respectively.

In Appendix we represent a polynomial approximation for $I,J$ functionals
(14). Since the equations obtained from (18) and (19) are quite cumbersome,
here we just briefly analyze the results.

In the limit of closely spaced solitons and $\Theta=0$, the population
imbalance $z$ at equilibrium depends only on $\Omega$ and obeys

$\frac{\Omega}{\Lambda}=1.2z^{7}-8z^{5}+15z^{3}-12.5z.$ (20)

Similarly, for fixed soliton phase difference $\Theta=\pi$ we have

$\frac{\Omega}{\Lambda}=1.2z^{7}-3.2z^{5}+12.3z^{3}-2z.$ (21)

We plot the graphical solutions of Eqs. (20), (21) in Fig. 2; blue and red
curves characterize the right parts of Eqs. (20), (21), respectively. Straight
lines in Fig. 2 correspond to different values of $\Omega/\Lambda$ ratio.
These lines cross the curves in the points which indicate the solutions of
Eqs. (20), (21). Notice, that the solid blue and red curves denote the values
of $\Omega/\Lambda$ and $z$ corresponding to the stable SS solutions, while
the dotted ones describe parametric unstable solutions. As seen from Fig. 2 at
phase difference $\Theta=0$ there exists one stable SS solution for any
$z\in[-0.7;0.7]$ and only unstable solutions for $|z|>0.7$. At
$|\Omega/\Lambda|>1.55\pi$ no SS solutions exist.

Figure 2: Normalized frequency spacing $\Omega/\Lambda$ (dashed lines) versus
reduced population imbalance $z$ for Eq. (20) (blue line) and Eq. (21) (red
line), respectively. Dashed regions correspond to unstable solutions.

On the other hand, at $\Theta=\pi$ there exists a tiny region
$-0.1\pi\leq\Omega/\Lambda\leq 0.1\pi$ possessing two SS solutions
simultaneously. One stable SS solution exists within the domain
$0.1\pi<|\Omega/\Lambda|\leq 2.64\pi$.

###  3.3 SS solutions for $\Theta=0,\pi$ and $\Omega\simeq 0$

At $\Omega=0$ Eqs. (20), (21) admit the SS solutions, which look like:

$\displaystyle z=0,$ $\displaystyle\Theta=0;$ (22a) $\displaystyle z=0,$
$\displaystyle\Theta=\pi;$ (22b) $\displaystyle z^{2}\approx 0.17$
$\displaystyle\Theta=\pi.$ (22c)

As seen from Eq. (22b), at relative phase $\Theta=\pi$ (19) possesses three
solutions: a parametrically unstable solution occurs at $z=0$ and two
degenerate SS solutions appear for $z=\pm z_{0}$. $z_{0}$ varies from $0.41$
at $\Delta\approx 0$ to $0.64$ at $\Delta\approx 2.8$ for non-zero soliton
inter-distance, respectively. For $\Delta>2.8$ these SS solutions do not
exist.

In Fig. 3 we represent the more general analysis of SS solutions for
$\Theta=0$ as functions of inter-soliton distance $\Delta$ for different
$\Omega$. For that we exploit the sixth-order polynomial approximation, see
Appendix. In particular, at $\Omega\simeq 0$ there exists one solution at
$z=0$, stable at $\Delta\leq\Delta_{c}\approx 0.5867$. For $\Delta>\Delta_{c}$
this solution becomes parametrically unstable.

On the other hand, for $\Delta>\Delta_{c}$ Eqs. (18) possesses the degenerate
SS solutions similar to the ones at $\Theta=\pi$. The bifurcation for
population imbalance $z$ occurs at $\Delta=\Delta_{c}$; in Fig. 3 the $z_{+}$
(upper,positive) and $z_{-}$ (lower, negative) branches characterize this
bifurcation. In the vicinity of $\Delta_{c}$ we can consider $z_{\pm}=\pm
z_{0}$, where

$z_{0}=1.2\sqrt{\Delta-\Delta_{c}}.$ (23)

At $\Omega\neq 0$ the behavior of SS solutions depending on distance $\Delta$
complicates - see green curves in Fig. 3. The solid curves correspond to SS
solutions for different $\Delta$, while the dotted ones describe the unstable
solutions. From Fig. 3 it is clearly seen that for $\left|\Omega\right|>0$
there is no bifurcation for population imbalance $z$ and two stationary
solution branches $z_{\pm}$ occur with $|z_{-}|>|z_{+}|$.

At relatively large values of parameter $\Omega/\Lambda$ only one SS solution
exists - see red curve in Fig. 3.

Figure 3: Population imbalance $z$ versus $\Delta$ for $\Theta=0$ and
different $\Omega$. Solid curves denote SS solutions of (18) and the dotted
line represents (parametrically) unstable solutions.

##  4 Mean-field dynamics

###  4.1 Small amplitude oscillations

We start our analysis here from small amplitude oscillations close to the SS
solutions (22). For that we linearize Eqs. (15) in the vicinity of (22),
assuming $0\leq\Delta<0.6$ and $\Omega<<1$. The first assumption allows to use
the approximation of $I,J$-functionals by the fourth-degree polynomial, see
Appendix.

For zero-phase oscillations, i.e. for $\Theta\approx 0$
($\cos\left[\Theta\right]\approx 1$, $\sin\left[\Theta\right]\approx\Theta$),
from (15) we obtain

$\ddot{z}+\omega_{0}^{2}(\Delta)z=f_{0}(\Delta)\Omega$ (24)

with the solution

$z(\tau)=A\cos\left[\omega_{0}\tau\right]-\Omega\frac{f_{0}}{\omega_{0}^{2}},$
(25)

where $A$ and $\omega_{0}(\Delta)=13.4\sqrt{0.37-\Delta^{2}-0.25\Delta}$ are
the amplitude and angular frequency of oscillations, respectively. The last
term in (25) with $f_{0}(\Delta)=5.36-0.8\Delta-4.22\Delta^{2}$ plays the role
of constant "external downward displacement force" that vanishes at
$\Omega\simeq 0$. Notably, at $\Delta>0.5$ oscillations become anharmonic and
$z(\tau)$ diverges at $\Delta>0.5867$. For $\Delta=0$ the frequency of
oscillations approaches $\omega\approx 8.15$, that agrees with numerical
solution of (15).

At $\Delta=\Delta_{c}\simeq 0.5867$ the SS solution (22) splits into two
degenerate solutions with $z=\pm z_{0}$ and $z_{0}$ determined by (23), see
Fig. 3. Near these points the equation, similar to (24), has a form

$\ddot{z}+\omega^{2}z=-18\Delta_{-}\sqrt{\Delta_{-}}-f(\Delta_{-})\Omega$ (26)

that implies a solution

$z(\tau)=\pm\left(1.2-\frac{18\Delta_{-}}{\omega^{2}}\right)\sqrt{\Delta_{-}}+A\cos[\omega\tau]-\Omega\frac{f(\Delta_{-})}{\omega^{2}},$
(27)

where $\Delta_{-}\equiv\Delta-\Delta_{c}$,
$\omega=14.53\sqrt{\Delta_{-}-4.48\Delta_{-}^{2}+17.8\Delta_{-}^{3}-53.5\Delta_{-}^{4}}$
is the angular frequency of oscillations,
$f=3.4-7.26\Delta_{-}+11\Delta_{-}^{2}$ is the "external" force. The relative
error for Eq. (27) is less than $5\%$.

In the vicinity of SS points determined by Eq. (22c), we obtain $\pi$-phase
oscillations characterized by

$z(\tau)=\pm
z_{0}+A\cos\left[\omega_{\pi}\tau\right]+\frac{f_{\pi}}{\omega_{\pi}^{2}}\Omega$
(28)

with $\omega_{\pi}=\sqrt{2-0.9\Delta^{2}-0.3\Delta}$,
$f_{\pi}=0.1\left(\Delta^{2}+0.38\Delta+5.5\right)$, and $z_{0}$ determined in
(22c). For $\Omega\simeq 0$ and $\Delta=0$ angular frequency is
$\omega_{\pi}\approx 1.42$ that is much smaller than in zero-phase regime.

The analysis of (15) in the vicinity of (22b) reveals that this solution is
parametric unstable, and highly nonlinear behavior is expected. Indeed, direct
numerical simulation demonstrates anharmonic dynamics plotted in Fig. 4. For
$0<|z|<0.5$ the nonlinear regime of self-trapping is observed, which turns
into nonlinear oscillations at $|z|>0.5$.

a)

b)

Figure 4: (a) the population imbalance $z(\tau)$ and (b) effective phase
difference $\Theta(\tau)$ versus reduced time $\tau$ for $\Theta(0)=\pi$.

The analysis of SS solution (17) reveals the strong sensitivity to
$z$-perturbation, when the condition $z^{2}=1$ is violated, the high-amplitude
nonlinear oscillations occur. On the other hand, the solution (17) is robust
to phase perturbations, which is an important property for metrology.

###  4.2 Large separation limit, $\Delta>>1$

For very large distance $\Delta$ between the solitons, $\Delta>>1$, the atom
tunneling between them vanishes and the solitons become independent. Strictly
speaking, in the limit $\Delta\rightarrow\infty$ the functionals
$I,J\rightarrow 0$ and Eqs. (15) look like

$\displaystyle\dot{z}$ $\displaystyle=$ $\displaystyle 0;$ (29a)
$\displaystyle\dot{\Theta}$ $\displaystyle=$ $\displaystyle-\Omega+2z,$ (29b)

i.e. the population imbalance is a constant in time and the running-phase
regime establishes.

For large but finite $\Delta$, SS solution $z=\pm z_{0}$ with
$z_{0}\rightarrow 1$ exists for the zero-phase regime, $\Theta=0$; for
example, for $\Delta=10$ the SS population imbalance is $z_{0}\approx 0.96$.

###  4.3 Phase-space analysis

The dynamical behavior of the coupled soliton system can be generalized in
terms of a phase portrait of two dynamical variables $z$ and $\Theta$, as
shown in Figs. 5 and 6.

In Fig. 5 we represent $z-\Theta$ phase-plane for $\Omega=0$ and for different
(increasing) values of distance $\Delta$. We distinguish three different
dynamic regimes. Solid curves correspond to oscillation regime when $z(\tau)$
and $\Theta(\tau)$ are some periodic functions of normalized time, see e.g.
(25) and red curve in Fig. 4. The dashed curves in Fig. 5 indicate the self-
trapping regime, when $z(\tau)$ is periodic and the sign of $z$ does not
change, c.f. (28) and blue curve in Fig. 4. Physically, this is the
macroscopic quantum self-trapping (MQST) regime characterized by a nonzero
average population imbalance, when the most of particles are "trapped" within
one of the solitons. At the same time, the behavior of phase $\Theta(\tau)$
may be quite complicated but periodic in time. On the other hand, for the
running-phase regime depicted by dashed-dotted curves $\Theta(\tau)$ grows
infinitely, see green curve in Fig. 5(b). Due to the symmetry that takes place
at $\Omega=0$, the running-phase can be achieved only with the MQST regime,
see Fig. 5.

As seen from Fig. 5, central area of nonlinear Rabi-like oscillations between
the ground and first excited macroscopic states happen for a relatively small
inter-soliton distance $\Delta$ and are inherent to zero-phase oscillations,
see Fig. 5(a). As we discussed before, at $\Delta=\Delta_{c}\approx 0.5867$
this area splits into two regions characterized by the MQST regimes, Fig.
5(b). This splitting occurs due to the bifurcation of population imbalance,
cf. black curve in Fig. 3. These regions are moving away from each other with
growing $\Delta$, see Fig. 5(c-f). It is worth noticing the bifurcation effect
and occurrence of MQST states at zero-phase regime for coupled solitons in
Fig. 1 disappear for the condensates described by Gaussian states, cf. [29,
30].

The phase trajectories inherent to $\pi$-phase region
$\frac{\pi}{2}<\Theta<\frac{3\pi}{2}$ stay weakly perturbed until the second
critical value $\Delta\approx 2$, when the MQST regime in Fig. 5(d) changes to
Rabi-like oscillations (Fig. 5(e)) and then approaches the running-phase at
$\Delta\approx 6$, see Fig. 5(f).

At large enough $\Delta$ the particle tunneling vanishes and the zero-phase
MQST domains arise in the vicinity of population imbalance $z=\pm 1$, Fig.
5(f). The phase dynamics corresponds to the running-phase regime with
$z=\operatorname{\textrm{const}}$, see Fig. 5(f) and (29).

a)

b)

c)

d)

e)

f)

Figure 5: Phase plane $z-\Theta$ at $\Omega=0$ and for (a) $\Delta=0$; (b)
$\Delta=0.75$; (c) $\Delta=1.2$; (d) $\Delta=1.5$; (e) $\Delta=3.0$; (f)
$\Delta=10$.

For non-zero $\Omega$, the phase portrait becomes asymmetric, Fig. 6. To
elucidate the role of $\Omega$ we study the soliton interaction for a given
inter-soliton distance $\Delta=0.75>\Delta_{c}$ that corresponds to the one
after the bifurcation. As seen from Fig. 6(a), the phase portrait does not
change significantly for small $\Omega$, cf. Fig. 5(b).

One of the SS solutions for zero and $\pi$\- phase regimes disappears with
increasing $\Omega$; the running-phase regime establishes, Fig. 6(b). Further
increasing of $\Omega$ leads to vanishing of the SS solution for zero-phase,
Fig. 6(c).

Thus, phase portraits in Figs. 5, 6 provide the existence of degenerate SSs
for coupled solitons by varying inter-soliton distance $\Delta$ and $\Omega$.
Such solutions, as we show below, may be exploited for the macroscopic
superposition soliton states formation in the quantum approach.

a)

b)

c)

Figure 6: Phase plane $z-\Theta$ at $\Delta=0.75$ and for (a)
$\Omega=0.05\pi$; (b) $\Omega=\pi$; (c) $\Omega=1.5\pi$.

##  5 Quantum metrology with two-soliton states.

###  5.1 Phase estimation with macroscopic qubit states.

Suppose that some quantum system is prepared in state $|\psi\rangle$, which
carries information about some parameter $\Gamma$ that we would like to
estimate; in this work we are interested in fundamental bound for POVM
measurements and consider pure states of the quantum system.

In quantum metrology the sensitivity of some parameter $\Gamma$ estimation is
described by the error propagation formula given as (cf. [28])

$\sigma_{\Gamma}=\frac{\sqrt{\left\langle\psi|(\Delta\hat{\Pi})^{2}|\psi\right\rangle}}{\left|\frac{\partial\left\langle\psi|\hat{\Pi}|\psi\right\rangle}{\partial\Gamma}\right|},$
(30)

where
$\left\langle\psi|(\Delta\hat{\Pi})^{2}|\psi\right\rangle={\left\langle\psi|\hat{\Pi}^{2}|\psi\right\rangle-\left\langle\psi|\hat{\Pi}|\psi\right\rangle}^{2}$
is the variance of fluctuations of some operator $\hat{\Pi}$ that corresponds
to the measurement procedure. Typically, such procedures are based on
appropriate interferometric schemes and use quantum superpositions, which
contain required information about estimated parameter $\Gamma$. In the case
of SC-states, which presume macroscopic (non-orthogonal in general) states,
the measurement procedure requires some specification. In particular, we
assume that quantum system may be prepared in the state $|\psi\rangle$ that we
represent as

$|\psi\rangle=\frac{1}{\sqrt{2}}(|\pi_{0}\rangle+e^{-i\phi_{N}}|\pi_{1}\rangle).$
(31)

In (31) $\phi_{N}$ is a relative (estimated) phase between states
$|\pi_{0}\rangle$ and $|\pi_{1}\rangle$, which are defined as

$\displaystyle|\pi_{0}\rangle=c_{1}|\Phi_{1}\rangle-c_{2}|\Phi_{2}\rangle;$
(32a) $\displaystyle|\pi_{1}\rangle=c_{2}|\Phi_{1}\rangle-
c_{1}|\Phi_{2}\rangle,$ (32b)

$|\Phi_{1}\rangle$ and $|\Phi_{2}\rangle$ are two macroscopic states
representing two "halves" of the SC-states. In particular, operators
$\hat{\Pi}_{i}=|\pi_{i}\rangle\langle\pi_{i}|$ realize a projection onto the
superposition of states $|\Phi_{1,2}\rangle$, which generally are not
orthogonal to each other obeying the condition

$\left\langle\Phi_{1}\Big{|}\Phi_{2}\right\rangle=\eta.$ (33)

Simultaneously, we require the states in Eqs.(32) to fulfill the normalization
condition

$\left\langle\pi_{i}\Big{|}\pi_{j}\right\rangle=\delta_{ij},i,j=0,1.$ (34)

Now, we are able to determine the coefficients $c_{1,2}$, which fulfill to
Eqs.(33), (34) and look like

$c_{1,2}=\sqrt{\frac{1\pm\sqrt{1-\eta^{2}}}{2(1-\eta^{2})}}.$ (35)

In Eqs.(33), (35) parameter $\eta$ defines the distinguishability of the
states $|\Phi_{1,2}\rangle$. The case of $\eta=1$ in (33), (35) corresponds to
completely indistinguishable states $|\Phi_{1,2}\rangle$. In this case one can
assume that $|\Phi_{1}\rangle$ and $|\Phi_{2}\rangle$ represent the same
state.

On the other hand, the situation with $\eta=0$ characterizes in (33), (35)
completely orthogonal states $|\Phi_{1,2}\rangle$ that become possible if
$|\Phi_{1,2}\rangle$ approaches two-mode Fock states. In other words, this is
a limit of the $N00N$-state for which coupled solitons are examined.

Then, we define a complete set of operators $\hat{\Sigma}_{j}$, $j=1,2,3$ (cf.
[50])

$\displaystyle\hat{\Sigma}_{0}$ $\displaystyle=$
$\displaystyle\left|\pi_{1}\right\rangle\left\langle\pi_{1}\right|+\left|\pi_{0}\right\rangle\left\langle\pi_{0}\right|,$
(36a) $\displaystyle\hat{\Sigma}_{1}$ $\displaystyle=$
$\displaystyle\left|\pi_{1}\right\rangle\left\langle\pi_{1}\right|-\left|\pi_{0}\right\rangle\left\langle\pi_{0}\right|,$
(36b) $\displaystyle\hat{\Sigma}_{2}$ $\displaystyle=$
$\displaystyle\left|\pi_{0}\right\rangle\left\langle\pi_{1}\right|+\left|\pi_{1}\right\rangle\left\langle\pi_{0}\right|,$
(36c) $\displaystyle\hat{\Sigma}_{3}$ $\displaystyle=$ $\displaystyle
i(\left|\pi_{0}\right\rangle\left\langle\pi_{1}\right|-\left|\pi_{1}\right\rangle\left\langle\pi_{0}\right|),$
(36d)

which obey the SU(2) algebra commutation relation.

The meaning of sigma-operators is evident from their definitions (36). Due to
properties (34) the states $|\pi_{i}\rangle$ are suitable candidates for
macroscopic qubit states, which we can define by mapping
$|\pi_{0}\rangle\rightarrow|\mathbf{0}\rangle$ and
$|\pi_{1}\rangle\rightarrow|\mathbf{1}\rangle$, respectively, cf. [43, 60]. In
this form we can use them for POVM measurements defined with operators [61]

$\displaystyle E_{1}$ $\displaystyle\equiv$
$\displaystyle\frac{\sqrt{2}}{1+\sqrt{2}}|\mathbf{1}\rangle\left\langle\mathbf{1}\right|=\frac{1}{\sqrt{2}(1+\sqrt{2})}(\hat{\Sigma}_{0}+\hat{\Sigma}_{1});$
(37a) $\displaystyle E_{2}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\sqrt{2}(1+\sqrt{2})}(|\mathbf{0}\rangle-|\mathbf{1}\rangle)(\langle\mathbf{0}|+\langle\mathbf{1}|)=-\frac{1}{\sqrt{2}(1+\sqrt{2})}(\hat{\Sigma}_{1}+i\hat{\Sigma}_{3});$
(37b) $\displaystyle E_{3}$ $\displaystyle\equiv$ $\displaystyle
I-E_{1}-E_{2}.$ (37c)

Importantly, current quantum (photonic) technologies permit POVM tomography,
cf. [44].

Average values of sigma-operators in (36) may be obtained by means of (31),
(34) and look like

$\displaystyle\left\langle\hat{\Sigma}_{1}\right\rangle$ $\displaystyle=$
$\displaystyle 0,$ (38a)
$\displaystyle\left\langle\hat{\Sigma}_{2}\right\rangle$ $\displaystyle=$
$\displaystyle\cos[\phi_{N}],$ (38b)
$\displaystyle\left\langle\hat{\Sigma}_{3}\right\rangle$ $\displaystyle=$
$\displaystyle\sin[\phi_{N}].$ (38c)

From (38) it follows that only $\left\langle\hat{\Sigma}_{2,3}\right\rangle$
contain the information about desired phase $\phi_{N}$.

To estimate the sensitivity of phase measurement it is possible to assume that
$\phi_{N}=N\Gamma$ and use Eq. (30) with measured operator
$\hat{\Pi}\equiv\hat{\Sigma_{2}}$. Taking into account
$\left\langle\hat{\Sigma}_{2}^{2}\right\rangle=1$ for the variance of
fluctuations $\left\langle(\Delta\hat{\Sigma}_{2})^{2}\right\rangle$ we obtain

$\left\langle(\Delta\hat{\Sigma}_{2})^{2}\right\rangle=\sin^{2}[N\Gamma].$
(39)

Finally, from Eqs. (30), (39) for the phase error propagation we obtain

$\sigma_{\Gamma}=\frac{1}{N},$ (40)

that clearly corresponds to HL of arbitrary ($N$-linearly dependent) phase
estimation and explores the sigma-operator measurement procedure. Notice this
procedure can be mapped onto the parity measurement, cf. [37, 50].

###  5.2 Soliton SC-qubit states.

The phase estimation procedure described above enables to use two-soliton
quantum states (6) for frequency quantum metrology purposes. It is instructive
to represent soliton wave functions (12) in the form

$\displaystyle\Psi_{1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN}}{4}(1-z)\operatorname{\textrm{sech}}\left[\left(1-z\right)\left(\frac{uN}{4}x-\Delta\right)\right]e^{-i\frac{\theta}{2}};$
(41a) $\displaystyle\Psi_{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN}}{4}(1+z)\operatorname{\textrm{sech}}\left[\left(1+z\right)\left(\frac{uN}{4}x+\Delta\right)\right]e^{i\frac{\theta}{2}}.$
(41b)

Degenerate SS solutions of population imbalance obtained before (see e.g. (17)
and Fig. 3) enable to prepare various superposition soliton states for quantum
metrology purposes. In particular, for $\Theta=0$ from (6) we obtain

$\displaystyle\left|\Phi_{1}\right\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\left[\int_{-\infty}^{\infty}\left(\Psi_{1}^{(+)}\hat{a}_{1}^{\dagger}+\Psi_{2}^{(+)}\hat{a}_{2}^{\dagger}\right)dx\right]^{N}\left|0\right\rangle,$
(42a) $\displaystyle\left|\Phi_{2}\right\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\left[\int_{-\infty}^{\infty}\left(\Psi_{1}^{(-)}\hat{a}_{1}^{\dagger}+\Psi_{2}^{(-)}\hat{a}_{2}^{\dagger}\right)dx\right]^{N}\left|0\right\rangle,$
(42b)

for two "halves" of the SC-state, where

$\displaystyle\Psi_{1}^{(\pm)}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN}}{4}(1-z_{\pm})\operatorname{\textrm{sech}}\left[\left(1-z_{\pm}\right)\left(\frac{uN}{4}x-\Delta\right)\right];$
(43a) $\displaystyle\Psi_{2}^{({\pm})}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN}}{4}(1+z_{\pm})\operatorname{\textrm{sech}}\left[\left(1+z_{\pm}\right)\left(\frac{uN}{4}x+\Delta\right)\right].$
(43b)

In Eqs. (43) $z_{+}$ and $z_{-}$ are two SS solutions corresponding to upper
and lower branches in Fig. 3, respectively. In (43) we omit the common
unimportant term $e^{-iN(\theta/2+\beta_{1}\tau)}$. In particular, for
$\Omega\approx 0$, we have $z_{\pm}\rightarrow\pm z_{0}$.

The scalar product for state (42) is

$\eta\equiv\left\langle\Phi_{1}\Big{|}\Phi_{2}\right\rangle=\left[\int_{-\infty}^{\infty}\left(\Psi_{1}^{(+)}\Psi_{1}^{(-)}+\Psi_{2}^{(+)}\Psi_{2}^{(-)}\right)dx\right]^{N}\equiv\epsilon^{N},$
(44)

where $\epsilon$ characterizes solitons wave functions overlapping. Assuming
non-zero and positive $\Omega$ for $\epsilon$ one can obtain

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigg{(}\left(1-\frac{z_{+}+z_{-}}{2}\right)\left(1-z_{+}\right)\left(1-z_{-}\right)\left(1-0.21\left[\frac{z_{+}-z_{-}}{2-(z_{+}+z_{-})}\right]^{2}\right)$
(45) $\displaystyle+$
$\displaystyle\left(1+\frac{z_{+}+z_{-}}{2}\right)\left(1+z_{+}\right)\left(1+z_{-}\right)\left(1-0.21\left[\frac{z_{+}-z_{-}}{2+(z_{+}+z_{-})}\right]^{2}\right)\Bigg{)}.$

In Fig. 7 we establish the principal features of coefficients (35) and
parameter $\epsilon$ (see the inset in Fig. 7) as functions of $\Delta$. The
value of $\Omega$ plays a significant role in the distinguishability problem
for states $|\Phi_{1}\rangle$ and $|\Phi_{2}\rangle$. In particular, at
$\Omega=0$, in the bifurcation point $\Delta=\Delta_{c}=0.5867$ we have
$\epsilon=1$ that implies indistinguishable states $|\Phi_{1}\rangle$ and
$|\Phi_{2}\rangle$, see red curve in the inset of Fig. 7. The coefficients
$c_{1,2}\rightarrow\infty$ in this limit.

However, even for the small (positive) $\Omega$ it follows from (45) that
$\epsilon\neq 1$ for any $\Delta>\Delta_{c}$, and states $|\Phi_{1}\rangle$,
$|\Phi_{2}\rangle$ are always distinguishable. In particular, it follows from
zero-phase solution (27) that
$z_{\pm}=\pm\left(1.2-\frac{18\Delta_{-}}{\omega^{2}}\right)\sqrt{\Delta_{-}}-\Omega\frac{f(\Delta_{-})}{\omega^{2}}$
and $|z_{+}|\neq|z_{-}|$. This situation is displayed by green curves in Fig.
7; the SS solutions possess $c_{1}=1.057$, $c_{2}=0.203$ for $c_{1,2}$ that
correspond to $\Delta\approx 0.647$, $\epsilon\approx 0.9056$ for
$\Omega=0.05\pi$.

From Fig. 7 it is evident that the coefficients $c_{1,2}$ rapidly approach
(due to the factor $N$) the levels $c_{1}=1$, $c_{2}=0$ (completely
distinguishable macroscopic SC soliton states), when $\Delta$ increases. In
this limit, as seen from Fig. 3, $z_{\pm}$ approaches $\pm z_{0}$, and from
(45) we obtain

$\epsilon\approx(1-z_{0}^{2})(1-0.21z_{0}^{2}).$ (46)

Practically, in this limit red and green curves coincide in Fig. 7.

Thus, we can exploit states (42) for metrological measurement purposes for
arbitrary phase $\phi_{N}$ estimation that we describe in Sec. 5.1. The phase
$\phi_{N}$ may be created after soliton SC-state formation by means of some
additional soliton interaction or collisions.

###  5.3 Frequency measurement, $\Gamma\equiv\Omega$.

Now we represent a particularly important case of measurement of angular
frequency $\Omega$ that characterizes energy spacing between the ground and
first excited macroscopic states. The SS solutions (17), which correspond to
the maximal population imbalance $z^{2}=1$, allow to prepare the maximally
path-entangled superposition state, a.k.a. $N00N$-state. As seen from (17b),
the solution with $z=1$ exist, when $-2(\pi-1)\leq\Omega/\Lambda\leq
2(\pi+1)$. Similarly, the domain of solution $z=-1$ is
$-2(\pi+1)\leq\Omega/\Lambda\leq 2(\pi-1)$. To achieve the superposition
$N00N$-state formation we require both solutions to exist simultaneously. This
restricts the domain of allowed $\Omega$ as $-2(\pi-1)\leq\Omega/\Lambda\leq
2(\pi-1)$.

Figure 7: Coefficients $c_{1,2}$ versus normalized inter-soliton distance
$\Delta$ for $N=10$ and different $\Omega$. The inset demonstrates behavior of
$\epsilon$ for different $\Delta$.

Substituting $z=\pm 1$ into (41) we obtain

$\displaystyle\Psi_{1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN}}{2}\operatorname{\textrm{sech}}\left[\left(\frac{uN}{2}x-2\Delta\right)\right]e^{-i\frac{\theta}{2}},$
(47a) $\displaystyle\Psi_{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{uN}}{2}\operatorname{\textrm{sech}}\left[\left(\frac{uN}{2}x+2\Delta\right)\right]e^{i\frac{\theta}{2}},$
(47b)

which are relevant to the $N00N$-state’s two "halves" defined as

$\displaystyle\left|0N\right\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\left[\int_{-\infty}^{\infty}\Psi_{1}\hat{a}_{1}^{\dagger}dx\right]^{N}\left|0\right\rangle,$
(48a) $\displaystyle\left|N0\right\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N!}}\left[\int_{-\infty}^{\infty}\Psi_{2}\hat{a}_{2}^{\dagger}dx\right]^{N}\left|0\right\rangle.$
(48b)

Considering the superposition of states (48) and omitting unimportant common
phase $e^{iN\left(0.5\theta^{(+)}-\beta_{2}t\right)}$ we arrive to (cf. (31))

$\left|N00N\right\rangle=\frac{1}{\sqrt{2}}\left(\left|N0\right\rangle+e^{-iN\Theta^{\prime}}\left|0N\right\rangle\right),$
(49)

that represents the $N00N$-state of coupled BEC solitons for our problem. In
(49)

$\Theta^{\prime}=\frac{\Theta^{(+)}+\Theta^{(-)}}{2}=\frac{1}{2}\left(\arccos\left[\frac{2\Lambda-\Omega}{2\pi\Lambda}\right]+\arccos\left[\frac{2\Lambda+\Omega}{2\pi\Lambda}\right]\right)$
(50)

is the phase shift that contains $\Omega$-parameter required for estimation.
Comparing Eq. (50) with (31) we can conclude that the $N00N$-state’s "halves"
$\left|N0\right\rangle$ and $\left|0N\right\rangle$ in (50) may be associated
with states $|\pi_{0}\rangle$ and $|\pi_{1}\rangle$, respectively. To estimate
the sensitivity of $\Omega$ measurement we use (30) with measured operator
$\hat{\Pi}\equiv\hat{\Sigma}$ defined as (cf.(36)(c))

$\hat{\Sigma}=\left|N0\right\rangle\left\langle
0N\right|+\left|0N\right\rangle\left\langle N0\right|.$ (51)

Since the states (48) are orthogonal, the mean value of (51) is

$\left\langle\hat{\Sigma}\right\rangle=\cos[N\Theta^{\prime}].$ (52)

Fig. 8 demonstrates $\left\langle\hat{\Sigma}\right\rangle$ as a function of
$\Omega/\Lambda$. Notice, the interference pattern in Fig. 8 exhibits
essentially nonlinear behavior for measured
$\left\langle\hat{\Sigma}\right\rangle$.

Figure 8: Mean value $\left\langle\hat{\Sigma}\right\rangle$ vs. normalized
angular frequency $\Omega/\Lambda$ for $N=200$.

The variance of fluctuations
$\left\langle(\Delta\hat{\Sigma})^{2}\right\rangle$ for the measured sigma-
operator reads as

$\left\langle(\Delta\hat{\Sigma})^{2}\right\rangle=\sin^{2}[N\Theta^{\prime}].$
(53)

Now, by using (30), (53) we can easily find the propagation error for
frequency $\Omega$ estimation as

$\sigma_{\Omega}=\frac{2\Lambda}{N}\left|\frac{\sqrt{4\pi^{2}-(2+\Omega/\Lambda)^{2}}\sqrt{4\pi^{2}-(2-\Omega/\Lambda)^{2}}}{\sqrt{4\pi^{2}-(2+\Omega/\Lambda)^{2}}-\sqrt{4\pi^{2}-(2-\Omega/\Lambda)^{2}}}\right|,$
(54)

Equation (54) is non-applicable for $\Omega=0$ since the denominator in (54)
turns to zero. We choose the optimal estimation area for the frequency
$\Omega$, where the best sensitivity is reached, in the vicinity of the domain
border at $\Omega/\Lambda\rightarrow 2(\pi-1)$. In this limit (54) approaches

$\sigma_{\Omega}\simeq\frac{10\Lambda}{N}\frac{1.65\sqrt{4.28-\Omega/\Lambda}}{1.65-\sqrt{4.28-\Omega/\Lambda}}.$
(55)

Equation (55) exhibits one of the important results of this work: for a given
$\Lambda$ that characterizes atomic condensate peculiarities Eq. (55)
demonstrates Heisenberg scaling for frequency measurement sensitivity, cf.
(40).

## Conclusion

In summary, we have considered the problem of two-soliton formation for 1D
BECs trapped effectively in a double-well potential. These soliton Josephson
junctions analytical solutions and corresponding phase portraits exhibit the
occurrence of novel macroscopic quantum selft-trapping (MQST) phases in
contrast to the condensates with only Gaussian wave functions. With these
soliton states, we have also explored the formation of the Schrödinger-cat
(SC) state in the framework of the Hartree approximation. In particular, we
have analyzed the distinguishabiltiy problem for binary (non-orthogonal)
macroscopic states. Compared to the results known in the literature, see e.g.
[37], finite frequency spacing $\Omega$ leads to the distinguishable
macroscopic states for condensate solitons. This circumstance may be important
for the experimental design of the SC-states.

The important part of this work is devoted to the applicability of predicted
states for quantum metrology. By utilizing the macroscopic qubits problem with
interacting BEC solitons, one can apply the sigma-operators to elucidate the
measurement and subsequent estimation of arbitrary phase, that linearly
depends on the particle number, up to the HL. It is worth mentioning the
sigma-operators relate to the POVM detection tomography procedure. On the
other hand, the phase estimation procedure for the phase-dependent sigma-
operator can be realized by means of the parity measurement technique that
produces the same accuracy for phase estimation. We have shown that in the
limit of soliton state solution with the population imbalance $|z|=\pm 1$ the
coupled soliton system admits the maximally path-entangled $N00N$-state
formation. The feasibility of frequency $\Omega$ estimation at the Heisenberg
level is also demonstrated.

In this paper we have not examined the loses and decoherence effects for the
quantum soliton system depicted in Fig. 1. Recently in [39] we examined this
problem for quantum solitons possessing simple Josephson coupling. From the
experimental point of view, the recent BEC soliton experiments with lithium
condensates demonstrated, that the collisions may be recognized as one of the
most detrimental effects [35]. In particular, as we established in [39], the
three-body and one-body losses may be unimportant at the time scales of few
tens of milliseconds that relevant to experimental conditions in [35].
Moreover, the purely quantum analysis of the problem demonstrated that the
superposition of Fock states, occurring at some specific parameters of the
system, behaves robust to few particle losses. The Detailed analysis of this
problem for interaction of solitons depicted in Fig. 1 we will publish
elsewhere.

## Funding

This work was partially supported by the Ministry of Science and Technology of
Taiwan under Grant Nos. MOST 108-2923-M-007-001-MY3 and 109-2112-M-007-019-MY3
and financially supported by the Grant of RFBR, No 19-52-52012 MHT_a.

## Disclosures

The authors declare no conflicts of interest.

## Appendix: Approximation of functionals $I$ and $J$.

Solutions of Eqs. (15) strictly depend on functionals $I$ and $J$ and their
derivatives $I^{\prime}\equiv dI/dz$ and $J^{\prime}\equiv dJ/dz$ defined in
(14). In Fig. 9 we represented them as two-dimension surfaces given in
$z-\Delta$ plane. From Fig. 9 (a,b) it is seen that $I,J$ approaches zero for
large $\Delta$ excluding edge domains where $|z|\simeq 1$. The behavior of
$I,J$, as it is follows from Fig. 9 (a,b), is not so evident for small
$\Delta$ inherent to $0\leq\Delta<1.5$. Thus for numerical estimations we use
the polynomial approximations of $I$ and $J$.

a)

b)

c)

d)

Figure 9: (a-b) The $I,J$ functionals and (c-d) their derivatives with respect
to $z$,$I^{\prime},J^{\prime}$, versus $z$ and $\Delta$.

In particular, within the domain $0\leq\Delta<0.6$, $I,J(z)$ can be
effectively approximated by the forth-order polynomials as follows

$\displaystyle I(z,\Delta)$ $\displaystyle\approx$ $\displaystyle
a_{I}(\Delta)z^{4}+b_{I}(\Delta)z^{2}+c_{I}(\Delta);$ (56a) $\displaystyle
J(z,\Delta)$ $\displaystyle\approx$ $\displaystyle
a_{J}(\Delta)z^{4}+b_{J}(\Delta)z^{2}+c_{J}(\Delta),$ (56b)

where the coefficients are the polynomials themselves:

$\displaystyle a_{I}$ $\displaystyle=$
$\displaystyle-\Delta^{2}-0.52\Delta+0.1;$ (57a) $\displaystyle b_{I}$
$\displaystyle=$ $\displaystyle 2\Delta^{2}+0.76\Delta-0.42;$ (57b)
$\displaystyle c_{I}$ $\displaystyle=$
$\displaystyle-1.16\Delta^{2}-0.24\Delta+1.33;$ (57c) $\displaystyle a_{J}$
$\displaystyle=$ $\displaystyle-2\Delta^{2}-0.72\Delta+0.4;$ (57d)
$\displaystyle b_{J}$ $\displaystyle=$ $\displaystyle
3.9\Delta^{2}+1.03\Delta+0.07;$ (57e) $\displaystyle c_{J}$ $\displaystyle=$
$\displaystyle-1.9\Delta^{2}-0.32\Delta+2.7.$ (57f)

At $0.6<\Delta\leq 1.5$ in (56) the sixth-order polynomial approximation is
required, that is

$\displaystyle I(z,\Delta)$ $\displaystyle\approx$ $\displaystyle
a_{I}(\Delta)z^{6}+b_{I}(\Delta)z^{4}+c_{I}(\Delta)z^{2}+d_{I}(\Delta);$ (58a)
$\displaystyle J(z,\Delta)$ $\displaystyle\approx$ $\displaystyle
a_{J}(\Delta)z^{6}+b_{J}(\Delta)z^{4}+c_{J}(\Delta)z^{2}+d_{J}(\Delta),$ (58b)

where

$\displaystyle a_{I}$ $\displaystyle=$ $\displaystyle
0.31\Delta^{2}-2.57\Delta+1.43;$ (59a) $\displaystyle b_{I}$ $\displaystyle=$
$\displaystyle 0.9\Delta^{2}+1.24\Delta-1.6;$ (59b) $\displaystyle c_{I}$
$\displaystyle=$ $\displaystyle-1.9\Delta^{2}+3.5\Delta-0.67;$ (59c)
$\displaystyle d_{I}$ $\displaystyle=$ $\displaystyle
0.69\Delta^{2}-2.21\Delta+1.85;$ (59d) $\displaystyle a_{J}$ $\displaystyle=$
$\displaystyle-1.5\Delta^{2}-0.13\Delta+0.89;$ (59e) $\displaystyle b_{J}$
$\displaystyle=$ $\displaystyle 4.62\Delta^{2}-4.78\Delta+0.15;$ (59f)
$\displaystyle c_{J}$ $\displaystyle=$
$\displaystyle-4\Delta^{2}+8.4\Delta-1.45;$ (59g) $\displaystyle d_{J}$
$\displaystyle=$ $\displaystyle 0.94\Delta^{2}-3.52\Delta+3.56.$ (59h)

An approximation that we use here provides error less than $4\%$ for any
$-1\leq z\leq 1$ and $\Delta$ in the mentioned domains.

## References

  * [1] T. Udem, R. Holzwarth, and T. W. Hänsch, “Optical frequency metrology,” Nature 416,233–237 (2002).
  * [2] S. B. Papp, K. Beha, P. Del’Haye, F. Quinlan, H. Lee, K. J. Vahala, and S. A. Diddams, “Microresonator frequency comb optical clock,” Optica 1, 10-14 (2014).
  * [3] M. G. Suh, Q. F. Yang, K. Y. Yang, X. Yi, and K. J. Vahala, “Microresonator soliton dual-comb spectroscopy,” Science 354, 600-603 (2016).
  * [4] A. Dutt, X. Ji, J. Cardenas, Y. Okawachi, K. Luke, A. L. Gaeta, and M. Lipson, “On-chip dual-comb source for spectroscopy,” Science advances 4, e1701858 (2018).
  * [5] D. T. Spencer, T. Drake, T. C. Briles, J. Stone, L. C. Sinclair, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, A. Bluestone, N. Volet, T. Komljenovic, L. Chang, S. H. Lee, D. Y. Oh, M. G. Suh, K. Y. Yang, M. H. P. Pfeiffer, T. J. Kippenberg, E. Norberg, L. Theogarajan, K. Vahala, N. R. Newbury, K. Srinivasan, J. E. Bowers, S. A. Diddams, and S. B. Papp, “An optical-frequency synthesizer using integrated photonics,” Nature 557, 81-85 (2018).
  * [6] P. Trocha, M. Karpov, D. Ganin, M. H. P. Pfeiffer, A. Kordts, S. Wolf, J. Krockenberger, P. Marin-Palomo, C. Weimann, S. Randel, W. Freude, T. J. Kippenberg, and C. Koos, “Ultrafast optical ranging using microresonator soliton frequency combs,” Science 359, 887-891 (2018).
  * [7] A. L. Gaeta, M. Lipson, and T. J. Kippenberg, “Photonic-chip-based frequency combs,” Nature Photon 13, 158-169 (2019).
  * [8] T. J. Kippenberg, A. L. Gaeta, M. Lipson, and M. L. Gorodetsky, “Dissipative Kerr solitons in optical microresonators,” Science 361, 567 (2018).
  * [9] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, “Quantum metrology with nonclassical states of atomic ensembles,” Rev. Mod. Phys. 90, 035005 (2018).
  * [10] D. Ananikian and T. Bergeman, “Gross-Pitaevskii equation for Bose particles in a double-well potential: Two-mode models and beyond,” Phys. Rev. A 73, 013604 (2006).
  * [11] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy, “Coherent oscillations between two weakly-coupled Bose-Einstein condensates: Josephson effects, $\pi$ oscillations, and macroscopic quantum self-trapping,” Phys. Rev. A 59, 620 (1999).
  * [12] A. J. Leggett, “Bose-Einstein condensation in the alkali gases: Some fundamental concepts,” Rev. Mod. Phys. 73, 307 (2001).
  * [13] S. Kohler and F. Sols, “Oscillatory decay of a two-component Bose-Einstein condensate,” Phys. Rev. Lett. 89, 060403 (2002).
  * [14] R. Gati and M. K. Oberthaler, “A bosonic Josephson junction,” J. Phys. B: Atom. Mole. Opt. Phys. 40, R61 (2007).
  * [15] J. I. Cirac, M. Lewenstein, K. Mølmer, and P. Zoller, “Quantum superposition states of Bose-Einstein condensates,” Phys. Rev. A 57, 1208 (1998).
  * [16] T. J. Haigh, A. J. Ferris, and M. K. Olsen, “Demonstrating mesoscopic superpositions in double-well Bose–Einstein condensates,” Opt. Comm. 283 3540-3547 (2010).
  * [17] Q. Y. He, S. G. Peng, P. D. Drummond, and M. D. Reid, “Planar quantum squeezing and atom interferometry,” Phys. Rev. A 84, 022107 (2011).
  * [18] G. Mazzarella, L. Salasnich, A. Parola, and F. Toigo, “Coherence and entanglement in the ground state of a bosonic Josephson junction: From macroscopic Schrödinger cat states to separable Fock states,” Phys. Rev. A 83, 053607 (2011).
  * [19] A. S. Sørensen and K. Mølmer, “Entanglement and Extreme Spin Squeezing,” Phys. Rev. Lett. 86, 4431 (2001).
  * [20] A. S. Sørensen, L. M. Duan, J. I. Cirac, and P. Zoller, “Many-particle entanglement with Bose-Einstein condensates,” Nature 409, 63-66 (2001).
  * [21] Q. Y. He, P. D. Drummond, M. K. Olsen, and M. D. Reid, “Einstein-Podolsky-Rosen entanglement and steering in two-well Bose-Einstein-condensate ground states,” Phys. Rev. A 86, 023626 (2012).
  * [22] L. Rosales-Zárate, B. J. Dalton, and M. D. Reid, “Einstein-Podolsky-Rosen steering, depth of steering, and planar spin squeezing in two-mode Bose-Einstein condensates,” Phys. Rev. A 98, 022120 (2018).
  * [23] Q. Y. He, P. D. Drummond, M. K. Olsen, and M. D. Reid, “Einstein-Podolsky-Rosen entanglement and steering in two-well Bose-Einstein-condensate ground states,” Phys. Rev. A 86, 023626 (2012).
  * [24] G. Puentes, G. Colangelo, R. J Sewell, and M. W. Mitchell, “Planar squeezing by quantum non-demolition measurement in cold atomic ensembles,” New J. Phys. 15, 103031 (2013).
  * [25] G. Vitagliano, G. Colangelo, F. M. Ciurana, M. W. Mitchell, R. J. Sewell, G. and Toth, “Entanglement and extreme planar spin squeezing,” Phys. Rev. A 97, 020301(R) (2018).
  * [26] G. Mazzarella, L. Salasnich, A. Parola, and F. Toigo, “Coherence and entanglement in the ground state of a bosonic Josephson junction: From macroscopic Schrödinger cat states to separable Fock states,” Phys. Rev. A 83, 053607 (2011).
  * [27] T. J. Haigh, A. J. Ferris, and M. K. Olsen, “Demonstrating mesoscopic superpositions in double-well Bose–Einstein condensates,” Opt. Comm. 283, 3540–3547 (2010).
  * [28] G. Tóth and I. Apellaniz, “Quantum metrology from a quantum information science perspective,” J. Phys. A: Math. Theor. 47, 424006 (2014).
  * [29] E. A. Ostrovskaya, Y. S. Kivshar, M. Lisak, B. Hall, F. Cattani, and D. Anderson, “Coupled-mode theory for Bose-Einstein condensates,” Phys. Rev. A 61, 031601(R) (2000).
  * [30] P. V. Elyutin and A. N. Rogovenko, “Stimulated transitions between the self-trapped states of the nonlinear Schrödinger equation,” Phys. Rev. E 63, 026610 (2001).
  * [31] P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-González, Emergent Nonlinear Phenomena in Bose–Einstein Condensates (Springer-Verlag Berlin Heidelberg, 2008).
  * [32] B. Eiermann, T. Anker, M. Albiez, M. Taglieber, P. Treutlein, K. P. Marzlin, and M. K. Oberthaler, “Bright Bose-Einstein gap solitons of atoms with repulsive interaction,” Phys. Rev. Lett. 92, 230401 (2004).
  * [33] K. E. Strecker, G. B. Partridge, A. G. Truscott, and R. G. Hulet, “Formation and propagation of matter-wave soliton trains,” Nature 417, 150-153 (2002)
  * [34] L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and C. Salomon, “Formation of a Matter-Wave Bright Soliton,” Science 296, 1290-1293 (2002).
  * [35] J. H. Nguyen, P. Dyke, D. Luo, B. A. Malomed, and R. G. Hulet, “Collisions of matter-wave solitons,” Nature Phys. 10, 918–922 (2014).
  * [36] S. Raghavan and G. P. Agrawal, “Switching and self-trapping dynamics of Bose-Einstein solitons,” J. Mod. Opt. 47, 1155-1169 (2000).
  * [37] D. V. Tsarev, S. M. Arakelian, Y. L. Chuang, R. K. Lee, and A. P. Alodjants, “Quantum metrology beyond Heisenberg limit with entangled matter wave solitons,” Opt. Express 26, 19583-19595 (2018).
  * [38] D. V. Tsarev, T. V. Ngo, R. K. Lee, and A. P. Alodjants, “Nonlinear quantum metrology with moving matter-wave solitons,” New J. Phys. 21, 083041 (2019).
  * [39] D. Tsarev, A. Alodjants, T. V. Ngo, and R.-K. Lee, “Mesoscopic quantum superposition states of weakly-coupled matter-wave solitons,” 22, 113016 (2020).
  * [40] A. Gilchrist, K. Nemoto, W. J. Munro, T. C. Ralph, S. Glancy, S. L. Braunstein, and G. J. Milburn, “Schrödinger cats and their power for quantum information processing,” J. Opt. B: Quantum Semiclass. Opt. 6, S828-S833 (2004).
  * [41] A. Ourjoumtsev, R. Tualle-Brouri, J. Laurat, and P. Grangier, “Generating optical Schrödinger kittens for quantum information processingm,“ Science 312, 83-86 (2006).
  * [42] J. S. Neergaard-Nielsen, M. Takeuchi, K. Wakui, H. Takahashi, K. Hayasaka, M. Takeoka, and M. Sasak, “Optical Continuous-Variable Qubit,” Phys. Rev. Lett. 105, 053602 (2010).
  * [43] S. Izumi, M. Takeoka, K. Wakui, M. Fujiwara, Kazuhiro Ema, and M. Sasaki, “Projective measurement onto arbitrary superposition of weak coherent state bases,” Sci. Rep. 8, 1-8 (2018).
  * [44] J. S. Lundeen, A. Feito, H. Coldenstrodt-Ronge, K. L. Pregnell, C. Silberhorn, T. C. Ralph, J. Eisert, M. B. Plenio, and I. A. Walmsley, “Tomography of quantum detectors,” Nature Phys. 5, 27-30 (2009).
  * [45] M. Takeoka, “Discrimination of Quantum States with Linear Optics and Continuous Photon Counting,” Optics and Spectroscopy 99, 418-424 (2005).
  * [46] G. S. Thekkadath, B. A. Bell, I. A. Walmsley, and A. I. Lvovsky, “Engineering Schrödinger cat states with a photonic even-parity detector,” Quantum 4, 239 (2020).
  * [47] T. Byrnes, K. Wen, and Y. Yamamoto, “Macroscopic quantum computation using Bose-Einstein condensates,” Phys. Rev. A 85, 040306(R) (2012).
  * [48] S. S. Demirchyan, I. Yu. Chestnov, A. P. Alodjants, M. M. Glazov, and A. V. Kavokin, “Qubits Based on Polariton Rabi Oscillator,” Phys. Rev. Lett. 112, 196403 (2014).
  * [49] K. Pawłowski, M. Fadel, P. Treutlein, Y. Castin, and A. Sinatra, Phys. Rev. A 95, 063609 (2017).
  * [50] C. C. Gerry, A. Benmoussa, and R. A. Campos “Parity measurements, Heisenberg-limited phase estimation, and beyond,” J. Mod. Opt. 54, 2177-2184 (2007).
  * [51] M. Sasaki and O. Hirota, “Optimum decision scheme with a unitary control process for binary quantum-state signals,” Phys. Rev. A 54, 2728 (1996).
  * [52] M. Takeoka, M. Sasaki, P. van Loock, and N. Lütkenhaus, “Implementation of projective measurements with linear optics and continuous photon counting,” Phys. Rev. A 71, 022318 (2005).
  * [53] C. W. Helstrom, “Quantum Detection and Estimation Theory,” Mathematics in Science and Engineering 123 (Academic Press, New York, 1976).
  * [54] J.M. Geremia, “Distinguishing between optical coherent states with imperfect detection,” Phys. Rev. A 70, 062303 (2004).
  * [55] J. R. Anglin and A. Vardi, “Dynamics of a two-mode Bose-Einstein condensate beyond mean-field theory,” Phys. Rev. A 64, 013605 (2001).
  * [56] V. I. Karpman and V. V. Solov’ev, “A perturbational approach to the two-soliton systems,” Phys. D 3, 487–502 (1981).
  * [57] Y. Lai and H. A. Haus, “Quantum theory of solitons in optical fibers. I. Time-dependent Hartree approximation," Phys. Rev. A 40, 844 (1989).
  * [58] Y. Lai and H. A. Haus, “Quantum theory of solitons in optical fibers. II. Exact solution," Phys. Rev. A 40, 854 (1989).
  * [59] A. P. Alodjants and S. M. Arakelian, “Quantum chaos and its observation in coupled optical solitons," Zh. Eksp. i Teor. Fiz., 107, 1792 (1995).
  * [60] P. T. Cochrane, G. J. Milburn, and W. J. Munro, “Macroscopically distinct quantum-superposition states as a bosonic code for amplitude damping,” Phys. Rev. A 59, 2631 (1999).
  * [61] M. A. Nielsen, I. L. Chuang, Quantum Computation and Quantum Information (Cambridge University Press, 2010).