Spin-orbit coupling is a fundamental quantum phenomenon that couples a particle’s intrinsic spin with its orbital angular momentum, and it finds applications in various fields, including spintronic devices, topological insulators, quantum computing, and high-precision measurements [1, 2]. Experimentally, SOC has been successfully implemented in Bose-Einstein condensates (BECs) through artificial gauge fields [3, 4, 5, 6, 7]. Given that the double-well potential is a standard model for studying tunneling and Josephson physics, SO-coupled BECs, where two internal atomic states are coupled to each other via a pair of counterpropagating laser beams, have received considerable attention when placed in such potentials. To date, research on this topic has been conducted from the viewpoints of both full-quantum and mean-field treatments, with particular emphasis on the dynamics of Josephson oscillations[8, 9, 10, 11, 12, 13, 14], self-trapping[15], flat bands in the energy spectrum[16], and dynamical suppression[17, 18, 19, 20, 21]. Here, two hyperfine atomic states are treated as spin (pseudospin), which serve as additional internal degrees of freedom, thereby opening up new possibilities for investigating the Josephson effects.

In recent years, the study of open systems has garnered significant interest among researchers. Experimentally, open Bose-Einstein condensate (BEC) systems have been effectively created, and the dissipation can be precisely engineered [22]. The theoretical frameworks for open systems encompass the Lindblad master equation and non-Hermitian Hamiltonians. Non-Hermitian Hamiltonians are introduced as a phenomenological approximation theory, and their implementation is achievable using the Feshbach projection operator method [23]. A rich variety of novel phenomena can emerge in open systems [24, 25, 26]. Among them, non-Hermitian Hamiltonians with parity-time ($\mathcal{PT}$) symmetry are particularly intriguing, as they exhibit a symmetry-breaking phase transition at which the spectrum changes from all real (the $\mathcal{PT}$ phase) to complex (the $\mathcal{PT}$-broken phase) when the non-Hermitian parameter exceeds a certain threshold [27]. Recently, the connection of the non-Hermitian physics to SO-coupled atomic systems has been explored, focusing primarily on ground states [28, 29], steady-state analytical solutions [30], phase transitions, and exceptional points (EPs) [24, 31, 32, 33]. The Floquet control of $\mathcal{PT}$-symmetric SO-coupled noninteracting BECs in a double-well potential has also been investigated [34]. A more encouraging advancement is that the role of dissipation in SO-coupled atomic systems has been studied experimentally [35], which demonstrates that non-Hermitian topological states can be controlled with SOC.

The quantum-classical correspondence (QCC) bridges the gap between macroscopic and microscopic perspectives, carrying significant theoretical and practical value. For Hermitian systems, the correspondence between the many-particle description and the mean-field approximation can be investigated in analogy with the usual quantum-classical correspondence for single-particle quantum mechanics. In particular, the mean-field dynamics of the Bose-Hubbard model and the QCC in this context have been extensively studied [8, 36, 37, 38, 39, 40]. Later, considerable attention has been paid to the open systems, and a series of mean-field results and beyond-mean-field approximations of the dissipative Bose-Hubbard model, based on a master equation description, have been reported in the literature [41, 42]. In parallel, there are numerous studies on the QCC of non-Hermitian Bose–Hubbard dimers [43, 44, 45, 46, 47, 48], where a complex Hamiltonian models a BEC in a double-well potential experiencing an effective decay from one of the modes. Such a purely decaying Hamiltonian can be mapped to a $\mathcal{PT}$-symmetric version including a source and sink of equal strength, and the corresponding classical (mean-field) dynamics can then be derived by investigating the quantum evolution of the SU(2) expectation values of the three angular momentum operators in the large $N$ limit [43, 44, 45, 46, 47]. However, due to the complexity of the non-Hermitian Bose-Hubbard model with SOC, we note that the non-Hermitian generalization of the discrete nonlinear Schrödinger equation resulting from the mean-field approximation has not been obtained, and the QCC in such systems has yet to be explored.

In this paper, we investigate the quantum-classical correspondence (QCC) of a non-Hermitian SO-coupled BEC in a double-well potential (a non-Hermitian SO-coupled bosonic junction), where non-Hermiticity arises from leakage in the right well [49]. Following the mean-field approximation methods from Refs. [41] and [46], we derive the equations of motion for the normalized two-point functions (also referred to as the single-particle reduced density matrix). We have conducted a rigorous mathematical proof confirming that the derived mean-field system embodies a classical Hamiltonian structure, and we have successfully derived the non-Hermitian version of the discrete nonlinear Schrödinger (Gross-Pitaevskii) equation. We demonstrate the effectiveness of the aforementioned calculation method for generalized classical Hamiltonian functions in a wider class of extended non-Hermitian, $N$-particle, $M$-mode Bose-Hubbard systems. By mapping the purely decaying model to a $\mathcal{PT}$-symmetric version, we find that when the effective SOC strength assumes half-integer values, $\mathcal{PT}$-symmetry-breaking appears even for arbitrarily small non-Hermitian parameters on both the mean-field and full-quantum sides. We subsequently examine the correspondence between full quantum dynamics and mean-field dynamics across different parameter regions. It is discovered that near the symmetry-breaking phase transition point, the quantum many-body fluctuation effects become significant, and the correspondence between the mean-field and quantum dynamics fails. While the system always resides in the symmetry-breaking regime when the SOC strength assumes half-integer values, the quantum dynamics are insensitive to changes in the number of particles. Additionally, we find that in both the mean-field and full-quantum systems, an extremely small SOC strength can induce synchronous periodic oscillations of the spin-up and spin-down components.

We consider a non-Hermitian many-particle spin-orbit coupled bosonic junction described by the Hamiltonian[46, 18, 21]

$$\begin{split}\hat{H}&=-J(e^{-i\pi\gamma}\hat{a}_{L\uparrow}^{\dagger}\hat{a}_{R\uparrow}+e^{i\pi\gamma}\hat{a}_{L\downarrow}^{\dagger}\hat{a}_{R\downarrow}+h.c.)+\Omega(\hat{a}_{L\uparrow}^{\dagger}\hat{a}_{L\downarrow}\\ &+\hat{a}_{R\uparrow}^{\dagger}\hat{a}_{R\downarrow}+h.c.)+\sum_{j,\sigma}(\frac{g}{2N}\hat{a}_{j\sigma}^{\dagger}\hat{a}_{j\sigma}^{\dagger}\hat{a}_{j\sigma}\hat{a}_{j\sigma}\\ &+\frac{g}{N}\hat{a}_{j\uparrow}^{\dagger}\hat{a}_{j\downarrow}^{\dagger}\hat{a}_{j\uparrow}\hat{a}_{j\downarrow})-i\beta\sum_{\sigma}\hat{a}_{R\sigma}^{\dagger}\hat{a}_{R\sigma},\end{split}$$

(1)

where $\hat{a}_{j,\sigma}^{\dagger}$ ($\hat{a}_{j,\sigma}$) is the creation (annihilation) operator for the ground-state wave function with spin up ($\sigma=\uparrow$) or down ($\sigma=\downarrow$) in the left ($j=L$) or right ($j=R$) well, respectively, obeying the bosonic commutation relation $[\hat{a}_{i,\sigma},\hat{a}_{j,\sigma^{\prime}}^{\dagger}]=\delta_{i,j}\delta_{\sigma,\sigma^{\prime}}$. $\Omega$ represents the Raman coupling strength, $\gamma$ is the SOC strength, $J$ is the tunneling amplitude between the two wells, and $g$ is the strength of the on-site interaction. Note that the Hamiltonian (1) commutes with the total number operator $\hat{N}=\sum_{j,\sigma}\hat{a}_{j,\sigma}^{\dagger}\hat{a}_{j,\sigma}$, and the additional imaginary part of the mode energy, $\beta$, describes the decay over time of the probability to find the entire many-particle ensemble in the four modes[46].

Let us now introduce the mean-field description for the non-Hermitian Bose-Hubbard dimer, referred to as model (1), in the spirit of the classical approximation as we take the limit of infinite $N$. For mathematical conciseness, we use $\hat{a}_{1},\hat{a}_{2},\hat{a}_{3},\hat{a}_{4}$ to denote $\hat{a}_{L\uparrow},\hat{a}_{L\downarrow},\hat{a}_{R\uparrow},\hat{a}_{R\downarrow}$ in model (1), which are the annihilation operators for the four modes. Our mean-field treatment can be applied to the extended non-Hermitian Bose-Hubbard model with a number of particles $N$ and an arbitrary number of $M$ modes, including arbitrary inter-mode particle interactions. To begin with, we can write the Hamiltonian $\hat{H}$ of the extended non-Hermitian many-particle Bose-Hubbard model as follows:

$$\begin{split}\begin{aligned} \hat{H}=&\hat{H}_{h}+\hat{H}_{a}+\hat{H}_{n},\\ \hat{H}_{h}=&\sum_{i,j=1}^{M}h^{h}_{i,j}\hat{a}_{i}^{\dagger}\hat{a}_{j},\\ \hat{H}_{a}=&\sum_{i,j=1}^{M}h^{a}_{i,j}\hat{a}_{i}^{\dagger}\hat{a}_{j},\\ \hat{H}_{n}=&\frac{1}{N}\sum_{i,j=1}^{M}h_{i,j}\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}\hat{a}_{i}\hat{a}_{j},\end{aligned}\end{split}$$

(2)

where $\hat{H}_{h}$ represents the Hermitian single-particle term, $\hat{H}_{a}$ is the anti-Hermitian single-particle term, and $\hat{H}_{n}$ is the Hermitian interacting term. These parts combine to form the complete Hamiltonian of a many-body system and satisfy the following relation:

$$\begin{split}\begin{aligned} \hat{H}_{h}^{\dagger}=\hat{H}_{h},\hat{H}_{a}^{\dagger}=-\hat{H}_{a},\hat{H}_{n}^{\dagger}=\hat{H}_{n}.\end{aligned}\end{split}$$

(3)

It is important to note here that any non-Hermitian operator can be decomposed into the sum of a Hermitian part and an anti-Hermitian part, thus we do not impose additional requirements on the single-body terms. Using Eq. (3), we derive that $h^{h}$, $h^{a}$, and $h$ possess the following properties:

$$\begin{split}\begin{aligned} (h^{h})^{\dagger}=h^{h},(h^{a})^{\dagger}=-h^{a},(h)^{*}=(h)^{\mathrm{T}}=h.\end{aligned}\end{split}$$

(4)

We readily observe that the model (1) we are considering is a specific case of model (2), with $M=4$.

For a general quantum system with the Hamiltonian operator $\hat{H}$, which is not necessarily Hermitian, the evolution of a pure state $|\Psi\rangle$ is determined by the Schrödinger equation: $id{|\Psi\rangle}/dt=\hat{H}|\Psi\rangle$. From this, we can derive the non-Hermitian generalized equation of motion for the expectation value of an operator $\langle\hat{A}\rangle:={\langle\Psi|\hat{A}|\Psi\rangle}/{\langle\Psi|\Psi\rangle}$ as follows,

$$i\frac{d\langle\hat{A}\rangle}{dt}=\langle[\hat{A},\hat{H}]\rangle+2\Delta(\hat{A}),$$

(5)

where $\Delta(\hat{A})=\langle\hat{H}_{a}\hat{A}\rangle-\langle\hat{H}_{a}\rangle\langle\hat{A}\rangle$. To derive the mean-field approximation, we start with the two-point functions $\sigma_{ij}:=\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle/N$, which are also known as the reduced single-particle density matrix. In the limit of large $N$, we assume that all particles occupy the same state and can be described by a single macroscopic wave function. The fully condensed states can be expressed as SU(M) coherent states, namely, $|\vec{x},N\rangle_{c}=\frac{1}{\sqrt{N!}}\left(\sum_{i=1}^{M}x_{i}\hat{a}_{i}^{\dagger}\right)^{N}|0\rangle,$ where $\vec{x}:=(x_{1},x_{2},\ldots,x_{M})\in\mathbb{C}^{M}$ and we refer to $x_{i}$ as the coherent state parameters, which can be physically interpreted as the amplitude of a single particle occupying the $i$-th mode. The mean-field equations of motion can be obtained by replacing the expectation values of the relevant operators in the two-point functions with their values in SU(M) coherent states, namely, $\sigma_{ij}:={\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle_{c}}/N=\frac{{}_{c}\langle\vec{x},N|\hat{a}_{i}^{\dagger}\hat{a}_{j}|\vec{x},N\rangle_{c}}{N_{c}\langle\vec{x},N|\vec{x},N\rangle_{c}}={x_{i}^{*}x_{j}}/n$, where $n=|\vec{x}|^{2}$. In the following, unless otherwise specified, for the sake of simplicity in notation, we omit the subscript “$c$” and use “$\langle\rangle$” to refer to “$\langle\rangle_{c}$”, which denotes the expectation values in SU(M) coherent states. Accoding to Eq. (5), the equations of motion for the two-point functions can be formulated by

	$\displaystyle i\frac{d\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle/N}{dt}=$	$\displaystyle\frac{1}{N}\Big{(}\langle[\hat{a}_{i}^{\dagger}\hat{a}_{j},\hat{H}]\rangle+\Delta(\hat{a}_{i}^{\dagger}\hat{a}_{j})\Big{)}$
	$\displaystyle=$	$\displaystyle\frac{1}{N}\Big{(}\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{h}]\rangle+\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{a}]\rangle+\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{n}]\rangle$
	$\displaystyle+$	$\displaystyle\langle[\hat{a}_{i}^{\dagger},\hat{H}_{h}]\hat{a}_{j}\rangle+\langle[\hat{a}_{i}^{\dagger},\hat{H}_{a}]\hat{a}_{j}\rangle+\langle[\hat{a}_{i}^{\dagger},\hat{H}_{n}]\hat{a}_{j}\rangle$
	$\displaystyle+$	$\displaystyle 2\langle\hat{H}_{a}\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle-2\langle\hat{H}_{a}\rangle\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle\Big{)}.$		(6)

Substituting Eq. (2) into Eq. (II), the resulting terms related to $\langle[\hat{a}_{i}^{\dagger}\hat{a}_{j},\hat{H}]\rangle$ in Eq. (II) are calculated as follows:

		$\displaystyle\frac{1}{N}\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{h}]\rangle=\frac{1}{N}\sum_{l=1}^{M}h^{h}_{j,l}\langle a_{i}^{\dagger}\hat{a}_{l}\rangle=\sum_{l=1}^{M}h^{h}_{j,l}\frac{x_{i}^{*}x_{l}}{n},$
		$\displaystyle\frac{1}{N}\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{a}]\rangle=\frac{1}{N}\sum_{l=1}^{M}h^{a}_{j,l}\langle a_{i}^{\dagger}\hat{a}_{l}\rangle=\sum_{l=1}^{M}h^{a}_{j,l}\frac{x_{i}^{*}x_{l}}{n},$
		$\displaystyle\frac{1}{N}\langle[\hat{a}_{i}^{\dagger},\hat{H}_{h}]\hat{a}_{j}\rangle=-\frac{1}{N}\sum_{l=1}^{M}h^{h}_{l,i}\langle a_{l}^{\dagger}\hat{a}_{j}\rangle=-\sum_{l=1}^{M}h^{h}_{l,i}\frac{x_{l}^{*}x_{j}}{n},$
		$\displaystyle\frac{1}{N}\langle[\hat{a}_{i}^{\dagger},\hat{H}_{a}]\hat{a}_{j}\rangle=-\frac{1}{N}\sum_{l=1}^{M}h^{a}_{l,i}\langle a_{l}^{\dagger}\hat{a}_{j}\rangle=-\sum_{l=1}^{M}h^{a}_{l,i}\frac{x_{l}^{*}x_{j}}{n},$
		$\displaystyle\frac{1}{N}\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{n}]\rangle=2\frac{N-1}{N}\sum_{l=1}^{M}h_{l,j}\frac{|x_{l}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n},$
		$\displaystyle\frac{1}{N}\langle[\hat{a}_{i}^{\dagger},\hat{H}_{n}]\hat{a}_{j}\rangle=-2\frac{N-1}{N}\sum_{l=1}^{M}h_{i,l}\frac{|x_{l}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n},$		(7)

and the results of the last two terms corresponding to $\Delta(\hat{a}_{i}^{\dagger}\hat{a}_{j})$ in Eq. (II) are given by

		$\displaystyle\frac{1}{N}\big{(}2\langle\hat{H}^{a}\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle-2\langle\hat{H}^{a}\rangle\langle\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle\big{)}$
		$\displaystyle=2\sum_{k=1}^{M}h^{a}_{k,i}\frac{x^{*}_{k}x_{j}}{n}-2\sum_{k,l=1}^{M}h^{a}_{k,l}\frac{x_{k}^{*}x_{i}^{*}x_{l}x_{j}}{n^{2}}.$		(8)

In the derivation of Eq. (II), we have utilized the following relations: $\langle\hat{a}^{\dagger}_{k}\hat{a}^{\dagger}_{i}\hat{a}_{l}\hat{a}_{j}\rangle=N(N-1)\frac{x_{k}^{*}x_{i}^{*}x_{l}x_{j}}{n^{2}}$, $\langle\hat{a}^{\dagger}_{k}\hat{a}_{l}\rangle\langle\hat{a}^{\dagger}_{i}\hat{a}_{j}\rangle=N^{2}\frac{x_{k}^{*}x_{i}^{*}x_{l}x_{j}}{n^{2}}$, and $\langle\hat{a}_{k}^{\dagger}\hat{a}_{l}\hat{a}_{i}^{\dagger}\hat{a}_{j}\rangle=\delta_{i,l}\langle\hat{a}^{\dagger}_{k}\hat{a}_{j}\rangle+\langle\hat{a}^{\dagger}_{k}\hat{a}_{i}^{\dagger}\hat{a}_{l}\hat{a}_{j}\rangle$, where $\delta_{ij}$ is the Kronecker delta function.

To perform the mean-field approximation, we take the macroscopic limit $N\rightarrow\infty$ with $g$ fixed, and the last two terms of Eq. (II) are given by

$$\begin{split}\begin{aligned} &\lim_{N\to\infty}\frac{1}{N}\langle\hat{a}_{i}^{\dagger}[\hat{a}_{j},\hat{H}_{n}]\rangle=2\sum_{l=1}^{M}h_{l,j}\frac{|x_{l}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n},\\ &\lim_{N\to\infty}\frac{1}{N}\langle[\hat{a}_{i}^{\dagger},\hat{H}_{n}]\hat{a}_{j}\rangle=-2\sum_{l=1}^{M}h_{i,l}\frac{|x_{l}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n}.\end{aligned}\end{split}$$

(9)

Combining the above derivations, we obtain the classical coherent evolution of the two-point function:

$$\begin{split}\begin{aligned} i\frac{d\sigma_{ij}}{dt}=&\sum_{l=1}^{M}(h^{h}_{j,l}+h^{a}_{j,l})\sigma_{il}-\sum_{l=1}^{M}(h^{h}_{l,i}+h^{a}_{l,i}-2h^{a}_{l,i})\sigma_{lj}\\ -&2\sum_{k,l=1}^{M}h^{a}_{k,l}\sigma_{kl}\sigma_{ij}+2\sum_{l=1}^{M}(h_{l,j}-h_{i,l})\sigma_{ll}\sigma_{ij}.\end{aligned}\end{split}$$

(10)

Based on the equation (10), we can already describe the classical (mean-field) behavior of the non-Hermitian many-particle system.

As an alternative to performing the coherent state approximation in the generalized equations of motion for the two-point function, we find that the mean-field description can also be expressed in terms of the generalized canonical equations of motion. Mathematically, in terms of the coherent state parameters $x_{i}$, the non-Hermitian generalization of the discrete nonlinear Schrödinger equation, which is equivalent to the evolution equation (10), can be formulated as:

$$\begin{split}\begin{aligned} i\frac{dx_{i}}{dt}=\frac{\partial H}{\partial x_{i}^{*}}=\sum_{j=1}^{M}\widetilde{H}_{i,j}x_{j},\end{aligned}\end{split}$$

(11)

where the classical Hamiltonian function is given by

$$\begin{split}\begin{aligned} H=\lim_{N\to\infty}\frac{n\langle\hat{H}\rangle}{N}.\end{aligned}\end{split}$$

(12)

In the canonical formulation of Eq. (11), $x_{i}$ serves as the generalized coordinates, with $x_{i}^{*}$ acting as the conjugate generalized momentum. To align it with the conventional canonical equations, we can formulate the equation conjugate to Eq. (11) as $i\frac{dx_{i}^{*}}{dt}=-\frac{\partial H^{*}}{\partial x_{i}}=-\frac{\partial H}{\partial x_{i}}+iQ_{i}^{(e)}$, where $iQ_{i}^{(e)}=2\frac{\partial H_{a}}{\partial x_{i}}$, and $H_{a}=\lim_{N\to\infty}\frac{n\langle\hat{H}_{a}\rangle}{N}$. Thus, $Q_{i}^{(e)}$ can be understood as an extra generalized force acting on the system, which is an active force not arising from potential energy.

By substituting Eq. (2) into Eq. (12), and taking the partial derivative of the classical Hamiltonian function $H$ with respect to $x_{i}^{*}$, we obtain the corresponding matrix appearing in Eq. (11), as follows:

$$\begin{split}\begin{aligned} \widetilde{H}=h^{h}+h^{a}+h^{n}(\vec{x}),\end{aligned}\end{split}$$

(13)

with

$$\begin{split}\begin{aligned} h_{i,j}^{n}(\vec{x})=2\sum_{l=1}^{M}\delta_{ij}h_{i,l}\frac{|x_{l}|^{2}}{n}-\sum_{k,l=1}^{M}\delta_{ij}h_{k,l}\frac{|x_{k}|^{2}|x_{l}|^{2}}{n^{2}}.\end{aligned}\end{split}$$

(14)

Utilizing the properties in Eq. (4), we can readily demonstrate that the nonlinear terms in Eq. (13) have the following properties: $(h^{n})^{*}(\vec{x})=(h^{n})^{\mathrm{T}}(\vec{x})=h^{n}(\vec{x})$.

Next, we will prove the equivalence of Eq. (10) and Eq. (11). To proceed, we need to express the dynamics of the two-point function in terms of the time derivatives of the coherent state parameters,

$$\begin{split}\begin{aligned} i\frac{d\sigma_{ij}}{dt}=&i\frac{d}{dt}\Big{(}\frac{x_{i}^{*}x_{j}}{n}\Big{)}\\ =&-\frac{x_{j}}{n}\Big{(}i\frac{dx_{i}}{dt}\Big{)}^{*}+\frac{x_{i}^{*}}{n}\Big{(}i\frac{dx_{j}}{dt}\Big{)}-\frac{x_{i}^{*}x_{j}}{n^{2}}\Big{(}i\frac{dn}{dt}\Big{)},\end{aligned}\end{split}$$

(15)

where

$$\begin{split}\begin{aligned} i\frac{dn}{dt}=\sum_{k=1}^{M}\Big{[}x_{k}^{*}\Big{(}i\frac{dx_{k}}{dt}\Big{)}-x_{k}\Big{(}i\frac{dx_{k}}{dt}\Big{)}^{*}\Big{]}.\end{aligned}\end{split}$$

(16)

By utilizing the nonlinear Schrödinger equation (11) and its complex conjugate, Eq. (15) is transformed as follows:

$$\begin{split}\begin{aligned} i\frac{d\sigma_{ij}}{dt}=&-\frac{x_{j}}{n}\Big{(}\sum_{l=1}^{M}\Big{[}h^{h}_{i,l}x_{l}+h^{a}_{i,l}x_{l}+h^{n}_{i,l}(\vec{x})x_{l}\Big{]}\Big{)}^{*}\\ +&\frac{x_{i}^{*}}{n}\Big{(}\sum_{l=1}^{M}\Big{[}h^{h}_{j,l}x_{l}+h^{a}_{j,l}x_{l}+h^{n}_{j,l}(\vec{x})x_{l}\Big{]}\Big{)}\\ -&\frac{x_{i}^{*}x_{j}}{n^{2}}\Big{(}i\frac{dn}{dt}\Big{)}.\end{aligned}\end{split}$$

(17)

Integrating Eqs. (11) and (16), and utilizing property (4), the total probability $n$ (the modulus of the coherent state parameters $\vec{x}$) decays as follows:

	$\displaystyle\frac{dn}{dt}=$	$\displaystyle-i\sum_{k=1}^{M}\Big{[}x_{k}^{*}\Big{(}\sum_{l=1}^{M}\Big{[}h^{h}_{k,l}x_{l}+h^{a}_{k,l}x_{l}+h^{n}_{k,l}(\vec{x})x_{l}\Big{]}\Big{)}$
	$\displaystyle-$	$\displaystyle x_{k}\Big{(}\sum_{l=1}^{M}\Big{[}h^{h}_{k,l}x_{l}+h^{a}_{k,l}x_{l}+h^{n}_{k,l}(\vec{x})x_{l}\Big{]}\Big{)}^{*}\Big{]}$
	$\displaystyle=$	$\displaystyle-i\sum_{k,l=1}^{M}\Big{(}2h^{a}_{k,l}x_{k}^{*}x_{l}\Big{)}.$		(18)

From Eqs. (17) and (II), using the property $(h^{n})^{*}(\vec{x})=(h^{n})^{\mathrm{T}}(\vec{x})=h^{n}(\vec{x})$, and the equality

	$\displaystyle\sum_{l=1}^{M}\Big{[}h^{n}_{j,l}(\vec{x})\frac{x_{i}^{*}x_{l}}{n}-h^{n}_{i,l}(\vec{x})\frac{x_{l}^{*}x_{j}}{n}\Big{]}$
	$\displaystyle=2\sum_{k=1}^{M}h_{j,k}\frac{|x_{k}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n}-2\sum_{k=1}^{M}h_{i,k}\frac{|x_{k}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n},$		(19)

we arrive at

	$\displaystyle i\frac{d\sigma_{ij}}{dt}=$	$\displaystyle-\Big{(}\sum_{l=1}^{M}\Big{[}h^{h}_{l,i}\frac{x_{l}^{*}x_{j}}{n}-h^{a}_{l,i}\frac{x_{l}^{*}x_{j}}{n}\Big{]}\Big{)}+\Big{(}\sum_{l=1}^{M}\Big{[}h^{h}_{j,l}\frac{x_{i}^{*}x_{l}}{n}+h^{a}_{j,l}\frac{x_{i}^{*}x_{l}}{n}\Big{]}\Big{)}$
	$\displaystyle+$	$\displaystyle 2\sum_{k=1}^{M}h_{j,k}\frac{|x_{k}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n}-2\sum_{k=1}^{M}h_{i,k}\frac{|x_{k}|^{2}}{n}\frac{x_{i}^{*}x_{j}}{n}$
	$\displaystyle-$	$\displaystyle\frac{x_{i}^{*}x_{j}}{n}\sum_{k,l=1}^{M}\Big{(}2h^{a}_{k,l}\frac{x_{k}^{*}x_{l}}{n}\Big{)}.$		(20)

Given the definition $\sigma_{ij}={x_{i}^{*}x_{j}}/{n}$, a comparison of Eq. (II) with Eq. (10) readily reveals their equivalence. Since Eq. (II) is derived from the nonlinear Schrödinger equation (11) solely with the aid of the definition of the two-point function, this completes our proof of the equivalence of Eq. (10) and Eq. (11).

For our specified model (1), the coefficient matrix that appears in the extended model (2) is explicitly given by

$\displaystyle\begin{split}&h^{h}=\begin{pmatrix}0&\Omega&-Je^{-i\pi\gamma}&0\\ \Omega&0&0&-Je^{i\pi\gamma}\\ -Je^{i\pi\gamma}&0&0&\Omega\\ 0&-Je^{-i\pi\gamma}&\Omega&0\end{pmatrix},\\ &h^{a}=\begin{pmatrix}0&0&0&0\\ 0&0&0&0\\ 0&0&-i\beta&0\\ 0&0&0&-i\beta\end{pmatrix},h=\frac{g}{2}\begin{pmatrix}1&1&0&0\\ 1&1&0&0\\ 0&0&1&1\\ 0&0&1&1\end{pmatrix}.\end{split}$

(21)

It is straightforward to verify that the three explicit matrices in Eq. (21) satisfy property (4). According to Eq. (12), the classical Hamiltonian function corresponding to our model (1) is given by

	$\displaystyle H$	$\displaystyle=\lim_{N\to\infty}\frac{n\langle\hat{H}\rangle}{N}$
		$\displaystyle=-J(e^{-i\pi x}x_{1}^{*}x_{3}+e^{i\pi x}x_{2}^{*}x_{4}+h.c.)+\Omega(x_{1}^{*}x_{2}+x_{3}^{*}x_{4}+h.c.)$
		$\displaystyle+\frac{g}{2n}(|x_{1}|^{4}+|x_{2}|^{4}+|x_{3}|^{4}+|x_{4}|^{4})+\frac{g}{n}(|x_{1}|^{2}|x_{2}|^{2}+|x_{3}|^{2}|x_{4}|^{2})$
		$\displaystyle-i\beta(|x_{3}|^{2}+|x_{4}|^{2}),$		(22)

with $n=|x_{1}|^{2}+|x_{2}|^{2}+|x_{3}|^{2}+|x_{4}|^{2}$. Accoding to Eq. (II), the total probability $n$ decays as

$\displaystyle\frac{dn}{dt}=-2\beta(|x_{3}|^{2}+|x_{4}|^{2}).$

(23)

Using Eq. (11) in conjunction with the classical Hamiltonian function specified in Eq. (II), the components of the unnormalized wave function, represented by the vector $|\varphi\rangle=(x_{1},x_{2},x_{3},x_{4})^{\mathrm{T}}$, adhere to the discrete non-Hermitian nonlinear Schrödinger equation

$$\begin{split}\begin{aligned} i\frac{d}{dt}|\varphi\rangle=\widetilde{H}|\varphi\rangle,\end{aligned}\end{split}$$

(24)

where

$$\begin{split}\begin{aligned} &\widetilde{H}=\begin{pmatrix}\xi_{1}&\Omega&-Je^{-i\pi\gamma}&0\\ \Omega&\xi_{2}&0&-Je^{i\pi\gamma}\\ -Je^{i\pi\gamma}&0&\xi_{3}&\Omega\\ 0&-Je^{-i\pi\gamma}&\Omega&\xi_{4}\end{pmatrix},\end{aligned}\end{split}$$

(25)

with

	$\displaystyle\xi_{1}=\xi_{2}=-\frac{g}{2}\Big{[}\frac{(|x_{1}|^{2}+|x_{2}|^{2})^{2}}{n^{2}}+\frac{(|x_{3}|^{2}+|x_{4}|^{2})^{2}}{n^{2}}\Big{]}+\frac{g}{n}\Big{(}|x_{1}|^{2}+|x_{2}|^{2}\Big{)},$
	$\displaystyle\xi_{3}=\xi_{4}=-\frac{g}{2}\Big{[}\frac{(|x_{1}|^{2}+|x_{2}|^{2})^{2}}{n^{2}}+\frac{(|x_{3}|^{2}+|x_{4}|^{2})^{2}}{n^{2}}\Big{]}+\frac{g}{n}\Big{(}|x_{3}|^{2}+|x_{4}|^{2}\Big{)}$
	$\displaystyle-i\beta.$		(26)

Apparently, the mean-field Hamiltonian $\widetilde{H}$ in Eq. (25) aligns with the expressions in (13) and (14) with the coefficient matrix (21) corresponding to our model (1).

By gauging away the (irrelevant) global phase via the following transformation:

$$|\varphi^{\prime}\rangle=e^{i\theta}|\varphi\rangle,\frac{d\theta}{dt}=-\frac{g}{2}\Big{[}\frac{(|x_{1}|^{2}+|x_{2}|^{2})^{2}}{n^{2}}+\frac{(|x_{3}|^{2}+|x_{4}|^{2})^{2}}{n^{2}}\Big{]},$$

(27)

Eq. (24) becomes

$$\begin{split}\begin{aligned} i\frac{d}{dt}|\varphi^{\prime}\rangle=H_{m}^{\prime}|\varphi^{\prime}\rangle,\end{aligned}\end{split}$$

(28)

where

$$\begin{split}\begin{aligned} &H_{m}^{\prime}=\begin{pmatrix}\xi_{1}^{\prime}&\Omega&-Je^{-i\pi\gamma}&0\\ \Omega&\xi_{2}^{\prime}&0&-Je^{i\pi\gamma}&\\ -Je^{i\pi\gamma}&0&\xi_{3}^{\prime}&\Omega\\ 0&-Je^{-i\pi\gamma}&\Omega&\xi_{4}^{\prime}\end{pmatrix},\end{aligned}\end{split}$$

(29)

with

$$\begin{split}\begin{aligned} &\xi_{1}^{\prime}=\xi_{2}^{\prime}=\frac{g}{n}\Big{(}|x_{1}|^{2}+|x_{2}|^{2}\Big{)},\\ &\xi_{3}^{\prime}=\xi_{4}^{\prime}=\frac{g}{n}\Big{(}|x_{3}|^{2}+|x_{4}|^{2}\Big{)}-i\beta,\end{aligned}\end{split}$$

(30)

which forms the basis for our subsequent numerical simulations of the mean-field dynamics. For vanishing dissipation, $\beta=0$, and the conserved norm $n=1$, Eq. (28) reduces to the Hermitian mean-field model of SO-coupled bosonic junction, as described in Refs. [50, 51, 52].

We first investigate the many-particle spectrum of our model (1), which is crucial for the dissipative dynamics. Because the system (1) is purely lossy, the imaginary parts of the eigenenergies are all negative, i.e., Im($E$) < 0, which signifies that the eigenstates are going to decay with time. Due to the normalized physical quantities of interest, we map the system to a $\mathcal{PT}$-symmetric version, as the explicit matrix structure in the Fock basis is essentially equivalent to the original model (1). By applying an imaginary energy shift, $\hat{H}\rightarrow\hat{H}+\frac{i\beta\hat{N}}{2}$, the $\mathcal{PT}$-symmetric model is given by

$$\begin{split}\hat{H}_{PT}&=\hat{H}+\frac{i\beta\hat{N}}{2}\\ &=-J(e^{-i\pi\gamma}\hat{a}_{1}^{\dagger}\hat{a}_{3}+e^{i\pi\gamma}\hat{a}_{2}^{\dagger}\hat{a}_{4}+h.c.)+\Omega(\hat{a}_{1}^{\dagger}\hat{a}_{2}\\ &+\hat{a}_{3}^{\dagger}\hat{a}_{4}+h.c.)+\sum_{i,j=1}^{4}(\frac{g}{2N}\hat{a}_{i}^{\dagger}\hat{a}_{j}^{\dagger}\hat{a}_{i}\hat{a}_{j})\\ &-\frac{i\beta}{2}[(\hat{a}_{3}^{\dagger}\hat{a}_{3}+\hat{a}_{4}^{\dagger}\hat{a}_{4})-(\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2})].\end{split}$$

(31)

With the definition of the parity operator $\hat{P}$ that interchanges the left and right wells: $\hat{a}_{1}\rightarrow\hat{a}_{3}$, $\hat{a}_{2}\rightarrow\hat{a}_{4}$, $\hat{a}_{3}\rightarrow\hat{a}_{1}$, $\hat{a}_{4}\rightarrow\hat{a}_{2}$, and the time-reversal operator $\hat{T}$ as $i\rightarrow-i$, it is evident that the Hamiltonian $\hat{H}_{PT}$ possesses $\mathcal{PT}$ symmetry: $[\hat{P}\hat{T},\hat{H}_{PT}]=0$. For $\hat{H}_{PT}$, the energy spectrum $\widetilde{E}$ remains fully real for $\beta$ below the critical value $\beta_{c}$, and becomes complex (appears in complex-conjugate pairs) for $\beta>\beta_{c}$. Thus, $\beta_{c}$ represents a phase transition between the broken and unbroken $\mathcal{PT}$ symmetry regions. Given that the energy spectrum of $\hat{H}_{PT}$ is essentially a shift in the imaginary part of the energy spectrum of $\hat{H}$, i.e., $\widetilde{E}=E+i\beta N/2$, this means that the imaginary parts of the eigenvalues of $\hat{H}$ are symmetrically distributed around $-\beta N/2$.

In Fig. 1, we have numerically investigated the dependence of the symmetry-breaking threshold $\beta_{c}$ on the interaction strength $g$, the particle number $N$, and the SOC strength $\gamma$. A key finding is that when the SOC strength $\gamma=0.5$, the $\mathcal{PT}$symmetry is fragile to nonzero non-Hermiticity, indicating that the system is in the $\mathcal{PT}$-broken phase for any arbitrarily small value of $\beta$. We also observe that the zero symmetry-breaking threshold $\beta_{c}$ for $\gamma=0.5$ is a universal characteristic applicable to any particle number and interaction strength. In Fig. 1 (c), we find that the symmetry-breaking threshold $\beta_{c}$ is very close to zero when $\gamma$ is approximately 0.2 and 0.8 for $N=20$, but this is accidental, as zero $\beta_{c}$ can be removed by changing the particle number.

To provide some initial physical insights into the key observations presented in Fig. 1, we first consider the limiting case of zero atomic interaction. We can analytically treat the noninteracting problem by employing the representation of bosonic creation and annihilation operators in generalized Bargmann space [53, 54, 55], with the following substitution:

$$\begin{split}&\hat{a}_{i}^{\dagger}\rightarrow Z_{i},\hat{a}_{i}\rightarrow\frac{\partial}{\partial Z_{i}},\\ &|\psi\rangle=\sum_{i,j,k,l}f(i,j,k,l)(\hat{a}_{1}^{\dagger})^{i}(\hat{a}_{2}^{\dagger})^{j}(\hat{a}_{3})^{k}(\hat{a}_{4})^{l}|0\rangle\\ &\rightarrow\psi(Z_{1},Z_{2},Z_{3},Z_{4})=\sum_{i,j,k,l}f(i,j,k,l)Z_{1}^{i}Z_{2}^{j}Z_{3}^{k}Z_{4}^{l}.\end{split}$$

(32)

In the Bargmann space, the eigenvalue equation $\hat{H}_{PT}|x_{1},x_{2},x_{3},x_{4},N\rangle=\widetilde{E}|x_{1},x_{2},x_{3},x_{4},N\rangle$ can be expressed as follows:

$$\begin{split}\begin{aligned} \hat{H}_{PT}^{B}f_{x,N}(Z)=&\widetilde{E}f_{x,N}(Z)\\ =&N\frac{1}{\sqrt{N!}}\Big{[}\sum_{i,j=1}^{4}\Big{(}h^{h}_{i,j}+h^{a}_{i,j}+\frac{i\beta}{2}\delta_{i,j}\Big{)}x_{j}Z_{i}\Big{]}\\ &\times\Big{(}x_{1}Z_{1}+x_{2}Z_{2}+x_{3}Z_{3}+x_{4}Z_{4}\Big{)}^{N-1},\end{aligned}\end{split}$$

(33)

where the $N$-particle coherent state $|x_{1},x_{2},x_{3},x_{4},N\rangle$ is represented by the analytical function $f_{x,N}(Z)=\frac{1}{\sqrt{N!}}(x_{1}Z_{1}+x_{2}Z_{2}+x_{3}Z_{3}+x_{4}Z_{4})^{N}$, and the operator corresponding to $\hat{H}_{PT}$ is represented by $\hat{H}_{PT}^{B}$ in the Bargmann space. From Eq. (33), it follows that

$$\begin{split}&(h^{h}+h^{a}+\frac{i\beta}{2})(x_{1},x_{2},x_{3},x_{4})^{\mathrm{T}}=\varepsilon(x_{1},x_{2},x_{3},x_{4})^{\mathrm{T}},\\ &\hat{H}_{PT}|x_{1},x_{2},x_{3},x_{4},N\rangle=N\varepsilon|x_{1},x_{2},x_{3},x_{4},N\rangle.\end{split}$$

(34)

This is a straightforward physical consequence: in the noninteracting case, if the state vector $|\psi\rangle=(x_{1},x_{2},x_{3},x_{4})^{\mathrm{T}}$ is an eigenstate of the single-particle Hamiltonian $h^{h}+h^{a}+\frac{i\beta}{2}$, then the corresponding $N$-particle coherent state $|x_{1},x_{2},x_{3},x_{4},N\rangle$ is also an eigenstate of the many-particle Hamiltonian $\hat{H}_{PT}$, with its eigenvalues being $N$ times the eigenvalues of the single-particle Hamiltonian. Consequently, insights into the many-particle spectrum can be derived directly from the single-particle Hamiltonian when $g=0$.