Perfect Transmission and Perfect Reflection of Bogoliubov Quasiparticles in a Dynamically Unstable Bose–Einstein Condensate

2.1 Model

We consider a binary BEC under a barrier potential at \(x=0\) (Fig. 1) and examine the tunneling properties of the quasiparticles through the barrier. The energy functional of the system is given by \begin{align} \mathcal{E}_{\text{binary}} &= \int d\boldsymbol{x}\sum_{m}\left[\frac{\hbar^{2}}{2M}\left|\frac{\partial}{\partial x}\Psi_{m}\right|^{2}{}+U(x)|\Psi_{m}|^{2}\right]\notag\\ &\quad + \int d\boldsymbol{x}\left[\sum_{m}\frac{g}{2}|\Psi_{m}|^{4}+ g'|\Psi_{1}|^{2}|\Psi_{2}|^{2}\right], \end{align} (1) where \(\Psi_{m}\) (\(m=1,2\)) is the condensate wave function, M is the atomic mass, \(U(x)\) is the barrier potential, and \(g>0\) and \(g'>0\) are the intra- and interspecies interaction strengths, respectively. For the sake of simplicity, we assume that the condensate wave function and the barrier potential depend only on x in a three-dimensional system so that the problem becomes essentially one-dimensional. We also assume that the strengths of the intraspecies interactions for the \(m=1\) and 2 components are the same. The two components are miscible (immiscible) for \(g>g'\) (\(g<g'\)).^25,26) The potential \(U(x)\) becomes zero and the condensate density \(n(x)=\sum_{m}|\Psi_{m}|^{2}\) converges to a constant value \(n_{0}\) at \(x\to\pm\infty\).

By taking the functional derivative of the energy functional (1) with respect to \(\Psi^{*}_{m}\), we obtain the time-dependent GP equation: \begin{align} i\hbar \frac{\partial\Psi_{1}}{\partial t} &= \left[-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial x^{2}}+U(x)+g|\Psi_{1}|^{2}+g'|\Psi_{2}|^{2}\right]\Psi_{1}, \end{align} (2a) \begin{align} i\hbar \frac{\partial\Psi_{2}}{\partial t} &= \left[-\frac{\hbar^{2}}{2M}\frac{\partial^{2}}{\partial x^{2}}+U(x)+g|\Psi_{2}|^{2}+g'|\Psi_{1}|^{2}\right]\Psi_{2}. \end{align} (2b) Below, we consider an equal-population mixture of the two components and prepare a completely overlapped state: \(\Psi_{1}(x)=\Psi_{2}(x)\equiv\Phi(x)/\sqrt{2}\), where \(\Phi(x)\) satisfies \begin{equation} \mu \Phi = \left[-\frac{\hbar^{2}}{2M}\frac{d^{2}}{dx^{2}}+U(x)+\frac{g+g'}{2}|\Phi|^{2}\right]\Phi \end{equation} (3) with μ being the chemical potential. The initial state is a stationary solution of Eq. (2) regardless of whether the binary mixture is miscible or not. However, its stability markedly changes at \(g=g'\): When \(g>g'\) (miscible), the initial state is the ground state of the system; hence, the quasiparticle eigenfrequencies are all positive real numbers. When \(g<g'\) (immiscible), the initial state is dynamically unstable against phase separation, and complex eigenfrequencies appear.

Since \(U(x)\to 0\) and \(|\Phi|^{2}\to n_{0}\) (= const.) at \(x\to\pm\infty\), Eq. (3) leads to \(\mu=(g+g')n_{0}/2\). We describe the equations with rescaling the energy, length, and time scales by μ, \(\hbar/\sqrt{M\mu}\), and \(\mu/\hbar\), respectively. By introducing the interaction parameter β (\(-1\leq\beta\leq 1\)) as \begin{equation} g = (g+g')\frac{1+\beta}{2},\quad g' = (g+g')\frac{1-\beta}{2} \end{equation} (4) and rewriting the wave functions as \(\Psi_{1,2}(x)=\sqrt{n_{0}}\psi_{1,2}(x)\) and \(\Phi(x)=\sqrt{n_{0}}\phi(x)\), we obtain the dimensionless forms of Eqs. (2) and (3) as \begin{align} i\frac{\partial\psi_{1}}{\partial t} &= [\mathcal{L}_{0}+(1+\beta)|\psi_{1}|^{2}+(1-\beta)|\psi_{2}|^{2}]\psi_{1}, \end{align} (5a) \begin{align} i\frac{\partial\psi_{2}}{\partial t} &= [\mathcal{L}_{0}+(1+\beta)|\psi_{2}|^{2}+(1-\beta)|\psi_{1}|^{2}]\psi_{2}, \end{align} (5b) and \begin{equation} (\mathcal{L}_{0}-1+|\phi|^{2})\phi = 0, \end{equation} (6) where \begin{equation} \mathcal{L}_{0} = - \frac{1}{2}\frac{d^{2}}{dx^{2}}+U(x). \end{equation} (7)

The dynamics of quasiparticle excitations from a condensate at a low temperature is well described by the Bogoliubov theory. The BdG equation for a binary BEC is obtained by substituting \begin{equation} \begin{pmatrix} \psi_{1}\\ \psi_{2} \end{pmatrix}= e^{-it} \left[\frac{\phi}{\sqrt{2}} \begin{pmatrix} 1\\ 1 \end{pmatrix}+ \boldsymbol{u}e^{-iE t}+\boldsymbol{v}^{*}e^{+iE^{*} t}\right] \end{equation} (8) into Eq. (5) and linearizing the equation with respect to \(\boldsymbol{u}\) and \(\boldsymbol{v}\), where \(\boldsymbol{u}=(u_{1}(x),u_{2}(x))^{\text{T}}\) and \(\boldsymbol{v}=(v_{1}(x),v_{2}(x))^{\text{T}}\) are two-component spinors. The resulting eigenvalue equation is given by \begin{equation} \begin{pmatrix} (\mathcal{L}_{0}-1)\mathbf{1}+ H_{1} &H_{2}\\ -H^{*}_{2} &-[(\mathcal{L}_{0}-1)\mathbf{1}+ H^{*}_{1}] \end{pmatrix}\begin{pmatrix} \boldsymbol{u}\\ \boldsymbol{v}\end{pmatrix}= E \begin{pmatrix} \boldsymbol{u}\\ \boldsymbol{v}\end{pmatrix}, \end{equation} (9) where \(\mathbf{1}\) is the \(2\times 2\) identity matrix, and \begin{align} H_{1} &= \frac{|\phi|^{2}}{2} \begin{pmatrix} 3+\beta &1-\beta\\ 1-\beta &3+\beta \end{pmatrix}, \end{align} (10) \begin{align} H_{2} &= \frac{\phi^{2}}{2} \begin{pmatrix} 1+\beta &1-\beta\\ 1-\beta &1+\beta \end{pmatrix}. \end{align} (11)

When the two components completely overlap, the density- and spin-wave modes are decoupled in the BdG equation. Indeed, by defining the density-wave (phonon) mode \(u^{\text{d}}=u_{1}+u_{2}\), \(v^{\text{d}}=v_{1}+v_{2}\) and the spin-wave (magnon) mode \(u^{\text{s}}=u_{1}-u_{2}\), \(v^{\text{s}}=v_{1}-v_{2}\), BdG Eq. (9) is divided into the following two equations: \begin{align} &\begin{pmatrix} \mathcal{L}_{0} -1 + 2|\phi|^{2} &\phi^{2}\\ -(\phi^{*})^{2} &-(\mathcal{L}_{0}-1 +2|\phi|^{2}) \end{pmatrix}\begin{pmatrix} u^{\text{d}}\\ v^{\text{d}} \end{pmatrix}\notag\\ &\quad= E^{\text{d}} \begin{pmatrix} u^{\text{d}}\\ v^{\text{d}} \end{pmatrix}, \end{align} (12) \begin{align} &\begin{pmatrix} \mathcal{L}_{0} -1 + (1+\beta)|\phi|^{2} &\beta\phi^{2}\\ -\beta(\phi^{*})^{2} &-[\mathcal{L}_{0} -1 + (1+\beta)|\phi|^{2}] \end{pmatrix}\begin{pmatrix} u^{\text{s}}\\ v^{\text{s}} \end{pmatrix}\notag\\ &\quad= E^{\text{s}} \begin{pmatrix} u^{\text{s}}\\ v^{\text{s}} \end{pmatrix}. \end{align} (13) When \(U(x)=0\), the spectra of these equations exhibit gapless linear dispersions, corresponding to the NG phonon and NG magnon modes associated with the breaking of the U(1) gauge symmetry and the SO(2) spin rotation symmetry, respectively. (The spin here means the pseudo-spin-1/2 of the binary system.) In the case of \(U(x)\neq 0\), the existence of the NG modes is confirmed by the fact that \(u^{\text{d,s}}=-(v^{\text{d,s}})^{*}=\phi\) are the eigen solutions of Eqs. (12) and (13) with \(E^{\text{d,s}}=0\).

We note that GP Eq. (6) and BdG Eq. (12) for the phonon mode are identical to those for a scalar BEC. Hence, the tunneling problem is the same as that in the case of a scalar BEC discussed in Ref. 3. We therefore discuss the tunneling properties of the spin-wave mode in the rest of this section.

2.2 Bogoliubov spectrum in a uniform system

Before solving the tunneling problem, we analytically derive asymptotic forms of the quasiparticle wave function at \(x\to\pm\infty\). In the case of \(U(x)=0\), the stationary solution of Eq. (6) is \(\phi=1\). By substituting \((u^{\text{s}},v^{\text{s}})^{\text{T}}=(ue^{ikx},ve^{ikx})^{\text{T}}\) and \(\phi=1\) to Eq. (13), we obtain the BdG equation for spin-wave modes at \(U(x)=0\) as \begin{equation} \begin{pmatrix} \epsilon_{k}+ \beta &\beta\\ -\beta &-(\epsilon_{k} + \beta) \end{pmatrix}\begin{pmatrix} u\\ v \end{pmatrix}= E^{\text{s}} \begin{pmatrix} u\\ v \end{pmatrix}, \end{equation} (14) which has the eigenvalue \begin{equation} E^{\text{s}} = \sqrt{\epsilon_{k}(\epsilon_{k}+2\beta)}, \end{equation} (15) where \(\epsilon_{k}\equiv k^{2}/2\). Note that when \(\beta<0\), \(E^{\text{s}}\) becomes a pure imaginary value for small k and the system becomes dynamically unstable. This instability is caused by the immiscible components [\(g<g'\), which is equivalent to \(\beta<0\), see Eq. (4)] that are mixed in the initial state. Our interest in this paper is to clarify how such dynamically unstable modes are reflected or transmitted by the barrier potential. In the case of \(\beta\geq 0\), for which the system is dynamically stable, almost the same situation has been discussed in the previous works,^9,10) where the tunneling problem of spin waves in a spin-1 polar state has been investigated and the perfect transmission was observed in the low-energy limit. The BdG Eq. (12) for density-wave modes at \(U(x)=0\) has the eigenvalue \(E^{\text{d}}=\sqrt{\epsilon_{k}(\epsilon_{k}+2)}\), which is a positive real value for \(^{\forall} k\).

In the tunneling problem, we will solve the quasiparticle wave function injected from \(x=-\infty\) with an energy E. Here, we therefore calculate normalized eigenvectors of Eq. (14) as a function of an incident eigenvalue E. We should be careful in calculating the normalization constant for the BdG equation for a bosonic system since the BdG equation is a non-Hermitian matrix equation. For the case of Eq. (14), the eigenvector is normalized as \(|u|^{2}-|v|^{2}=1\) or −1 (\(u^{2}-v^{2}=1\) or −1) for real (pure imaginary) E. See Appendix for the details.

2.3 Tunneling properties of the spin-wave modes

We now solve BdG Eq. (13) and calculate the reflection and transmission probabilities of quasiparticles with energy E and momentum \(k_{\text{in}}\) injected from the left. As for the barrier potential \(U(x)\), we use the Gaussian potential given by \begin{equation} U(x) = U_{0}e^{-x^{2}/(2\sigma^{2})}. \end{equation} (26) We numerically solve BdG Eq. (13) by the finite element method with the asymptotic form of the wave function imposed at \(x\to\pm\infty\).

In the case of \(\beta\geq 0\), a quasiparticle with \(k_{\text{in}}=k_{1}\) injected from the left (\(e^{ik_{1} x}\)) is reflected to the left (\(e^{-ik_{1} x}\)) or transmitted to the right (\(e^{ik_{1} x}\)). In addition, localized modes at the potential barrier appear (\(e^{\pm q_{2} x}\)). The asymptotic form of the wave function at \(x\to\pm \infty\) is then given by

(i)

\(\beta\geq 0\), \(E\in\mathbb{R}\) with an incident wave \(e^{ik_{1} x}\): \begin{equation} \left\{ \begin{array}{l} \begin{pmatrix} u_{-\infty}\\ v_{-\infty} \end{pmatrix}= \begin{pmatrix} a_{r}\\ -b_{r} \end{pmatrix}(e^{ik_{1}x} + Re^{-ik_{1}x})+A \begin{pmatrix} b_{r}\\ a_{r} \end{pmatrix}e^{q_{2}x},\\ \begin{pmatrix} u_{+\infty}\\ v_{+\infty} \end{pmatrix}= T \begin{pmatrix} a_{r}\\ -b_{r} \end{pmatrix}e^{ik_{1}x} + B \begin{pmatrix} b_{r}\\ a_{r} \end{pmatrix}e^{-q_{2}x}. \end{array}\right. \end{equation} (27)

Here, R and T are the reflection and transmission coefficients, respectively, and A and B are the coefficients for the localized modes. The current conservation law requires \(|R|^{2}+|T|^{2}=1\).

In the case of \(\beta<0\), there are three possibilities: (ii) real positive \(E\in\mathbb{R}\) with \(k_{\text{in}}=k_{1}\), (iii) pure imaginary \(E=i\Delta\in i\mathbb{R}\) (\(\Delta>0\)) with \(k_{\text{in}}=-k_{4}\), and (iv) pure imaginary \(E=i\Delta\in i\mathbb{R}\) with \(k_{\text{in}}=k_{3}\). The results for \(\Delta<0\) are obtained by taking the complex conjugate of the wave functions obtained for \(\Delta>0\). The points corresponding to the incident wave are depicted in the dispersion relation in Fig. 2. Note that although the group velocity of a quasiparticle with a pure imaginary eigenvalue is zero, the causality is satisfied when we assume \(-d|\Delta|/dk\) to be a group velocity and choose the momentum of the incident, reflected, and transmitted waves. [The origin of the causality becomes clearer when the spectrum has a gap at \(k=0\), as in Eq. (43). We shall revisit this issue in Sect. 3.3.] Hence, the asymptotic forms of the wave function at \(x\to\pm \infty\) are given as follows:

(ii)	\(\beta<0\), \(E\in\mathbb{R}\) with \(k_{\text{in}}=k_{1}\): The asymptotic form is the same as Eq. (27).
(iii)	\(\beta<0\), \(E\in i\mathbb{R}\) with \(k_{\text{in}}=-k_{4}\): \begin{equation} \left\{ \begin{array}{l} \begin{pmatrix} u_{-\infty}\\ v_{-\infty} \end{pmatrix}= \begin{pmatrix} a_{c}\\ -b_{c} \end{pmatrix}(e^{-ik_{4}x} + Re^{ik_{4}x})+A \begin{pmatrix} b_{c}\\ a_{c} \end{pmatrix}e^{-ik_{3}x},\\ \begin{pmatrix} u_{+\infty}\\ v_{+\infty} \end{pmatrix}= T \begin{pmatrix} a_{c}\\ -b_{c} \end{pmatrix}e^{-ik_{4}x} + B \begin{pmatrix} b_{c}\\ a_{c} \end{pmatrix}e^{ik_{3}x}. \end{array}\right. \end{equation} (28)
(iv)	\(\beta<0\), \(E\in i\mathbb{R}\) with \(k_{\text{in}}=k_{3}\): \begin{equation} \left\{ \begin{array}{l} \begin{pmatrix} u_{-\infty}\\ v_{-\infty} \end{pmatrix}= \begin{pmatrix} b_{c}\\ a_{c} \end{pmatrix}(e^{ik_{3}x} + Re^{-ik_{3}x})+A \begin{pmatrix} a_{c}\\ -b_{c} \end{pmatrix}e^{ik_{4}x},\\ \begin{pmatrix} u_{+\infty}\\ v_{+\infty} \end{pmatrix}= T \begin{pmatrix} b_{c}\\ a_{c} \end{pmatrix}e^{ik_{3}x} + B \begin{pmatrix} a_{c}\\ -b_{c} \end{pmatrix}e^{-ik_{4}x}. \end{array}\right. \end{equation} (29)

Here, we further introduce the scaled coefficients \begin{equation} \tilde{A} = \sqrt{-\frac{k_{\text{sc}}}{k_{\text{in}}}} A,\quad \tilde{B} = \sqrt{-\frac{k_{\text{sc}}}{k_{\text{in}}}} B, \end{equation} (30) where \(k_{\text{sc}}\) is the real part of the momentum of the A term. For example, \(k_{\text{sc}}=0\), \(-k_{3}\), and \(k_{4}\) for the cases (ii), (iii), and (iv), respectively. Using \(\tilde{A}\) and \(\tilde{B}\), we write the current conservation law as \begin{equation} |T|^{2}+|R|^{2}+|\tilde{A}|^{2}+|\tilde{B}|^{2} = 1 \end{equation} (31) for all cases (i)–(iv).

In Fig. 3, we show the numerical results for the tunneling coefficients for the cases (i)–(iv). The results for the quasiparticle wave functions are depicted in Fig. 4. Figure 3(a) shows the transmission probability for real E (\(>0\)) with various values of β, obtained by imposing the asymptotic form of Eq. (27) [cases (i) and (ii)]. One can see that the perfect transmission occurs at the low-energy limit for \(\beta\geq 0\), i.e., \(|T|^{2}\) goes to unity as \(E\to 0\), which is consistent with the previous works.^9,10) On the other hand, for \(\beta<0\), the transmission probability becomes smaller than unity. This difference is understood as follows. The perfect transmission occurs when the quasiparticle wave function coincides with the condensate wave function: \(u(x)=-v^{*}(x)=\phi(x)\). In the present case, the incident momentum \(k_{1}\) goes to zero (nonzero) as \(E\to 0\) when \(\beta\geq 0\) (\(\beta<0\)); hence \(u(x)=-v^{*}(x)=\phi(x)\) is (is not) satisfied at \(E\to 0\). In Figs. 4(a) and 4(b), we show \(|u|^{2}-|v|^{2}\) at \(\beta=0.2\) and −0.2, respectively. [The integral \(\int (|u|^{2}-|v|^{2})\,dx\) gives the norm of the quasiparticle wave function for a real E, whereas it should vanish for a complex E. We therefore plot \(|u|^{2}+|v|^{2}\) for \(\beta<0\) as shown in Figs. 4(c) and 4(d). See also Appendix.] One can see from these figures that \(|u|^{2}-|v|^{2}\) for \(\beta=0.2\) becomes close to \(|\phi|^{2}\) as \(E\to 0\), whereas for \(\beta=-0.2\) at \(E\to 0\), it has a completely different x dependence from \(|\phi|^{2}\), which is consistent with the above discussion. The result in the previous works^4–6) that the phonon mode in a scalar BEC moving with the critical current does not exhibit the perfect transmission is due to the same reason, where u and \(v^{*}\) deviate from ϕ owing to local enhancement of density fluctuations around the potential barrier.⁶⁾

Note that \(k_{4}\) goes to zero and \(u(x)=-v^{*}(x)=\phi(x)\) is satisfied when E goes to zero along the imaginary axis. Hence, the perfect transmission at \(E\to 0\) occurs in this case. Figure 3(b) shows the E dependences of \(|T|^{2}\), \(|R|^{2}\), \(|\tilde{A}|^{2}\), and \(|\tilde{B}|^{2}\) at \(\beta=-0.2\). Here, we set the horizontal axis of Fig. 3(b) such that \(\epsilon_{k_{\text{in}}}\) is zero at the left end of the figure and \(\epsilon_{k_{\text{in}}}\) increases as one goes to the right. The region where E is pure imaginary is shaded in gray. We have numerically confirmed the current conservation law given by Eq. (31). In Fig. 3(b), the transmission probability \(|T|^{2}\) goes to unity as \(\mathop{\text{Im}}\nolimits E\to 0\) at the left end of the figure, indicating that the perfect transmission occurs even for the dynamically unstable modes. We also confirm that \(|u|^{2}+|v|^{2}\) shown in Fig. 4(c) becomes proportional to \(|\phi|^{2}\) as \(\mathop{\text{Im}}\nolimits E\to 0\), which is consistent with \(u(x)=-v^{*}(x)=\phi(x)\). On the other hand, at the point where E changes from pure imaginary values to real values [the point \(E=0\) and \(\epsilon_{k}>0\) in Figs. 2(c) and 2(d), which corresponds to the boundary between the gray shaded region and the white region in Fig. 3(b)], the incident momentum \(k_{3}\) remains nonzero in the limit of \(E\to 0\), which means \(u(x)=-v^{*}(x)=\phi(x)\) is not satisfied in this limit, resulting in the absence of the perfect transmission.

We also find from Fig. 3(b) that \(|\tilde{A}|^{2}\) goes to unity at \(E=i|\beta|\). At this point, \(|{\mathop{\text{Im}}\nolimits E}|\) takes its maximum value, i.e., \(E\in i\mathbb{R}\) and \(d\mathop{\text{Im}}\nolimits E/dk=0\). When E is pure imaginary, the A terms in Eqs. (28) and (29) represent the reflected wave with the momentum different from the incident one. Hence, \(|\tilde{A}|^{2}=1\) means the occurrence of perfect reflection. The origin of the perfect reflection is that the number of linearly independent solutions given in Eq. (22) decreases at \(k_{3}=k_{4}\), and the solution that satisfies the asymptotic forms of Eqs. (28) and (29) disappears. Thus, the incident wave and reflected wave (the A term) destructively interfere. The destructive interference can be confirmed in the quasiparticle wave function shown in Fig. 4(d), in which we plot \(|u|^{2}+|v|^{2}\) obtained for \(\beta=-0.2\) with the asymptotic form of Eq. (28). The wave function goes to zero as \(E\to i|\beta|\) even in the \(x<0\) region. Since we fix the amplitude of the incident wave, the reduction in \(|u|^{2}+|v|^{2}\) means the destructive interference. Perfect reflection was also predicted to occur for the long-wavelength limit of the spin-wave mode in a ferromagnetic BEC.^9,10) The perfect reflection we find in this study is different from that in the previous works in that the origin of the perfect reflection in the latter case is essentially the same as that of the ordinary quantum-mechanical tunneling of a free particle.

3.1 Model

Next, we consider the tunneling problem in a spin-1 polar BEC. The BdG equation for the spin-wave mode in this system is the same as that in the previous section except for the quadratic Zeeman energy \(q_{z}\) term. Since the quadratic Zeeman effect breaks the spin rotational symmetry, the spin-wave spectrum has a nonzero eigenvalue at \(k=0\). Here, the eigenvalue at \(k=0\) can be real or pure imaginary depending on the value of \(q_{z}\). We therefore focus on the \(q_{z}\) dependence of the tunneling properties. Below, we redefine the variables used in the previous section, so that the resulting BdG equation has the same form as that in the previous section except for the \(q_{z}\) term.

The energy functional of a spin-1 system is given by \begin{align} \mathcal{E}_{\text{spin-1}} &= \int d\boldsymbol{x}\sum_{m}\left\{\frac{\hbar^{2}}{2M}\left|\frac{\partial \Psi_{m}}{\partial x}\right|^{2}{}+[U(x)+q_{z} m^{2}]|\Psi_{m}|^{2}\right\}\notag\\ &\quad + \frac{1}{2}\int d\boldsymbol{x}[c_{0}n^{2}+c_{1}|\boldsymbol{F}|^{2}], \end{align} (32) where \(\Psi_{m}\) is the condensate wave function of the atoms in the magnetic sublevel \(m=1\), 0, and −1, M is the atomic mass, \(U(x)\) is the barrier potential, \(n(x)=\sum_{m}|\psi_{m}(x)|^{2}\) is the condensate density, and \(\boldsymbol{F}(x) =(F_{x}(x),F_{y}(x),F_{z}(x))\) is the spin density vector defined by \(\boldsymbol{F}(x)=\sum_{mm'}\Psi_{m}^{*}(x)\boldsymbol{S}_{mm'}\Psi_{m'}(x)\) with \(\boldsymbol{S}=(S_{x},S_{y},S_{z})\) being the spin-1 matrices. The interaction coefficients are given by \(c_{0}=4\pi\hbar^{2}(2a_{2}+a_{0})/(3M)\) and \(c_{1}=4\pi\hbar^{2}(a_{2}-a_{0})/(3M)\), where \(a_{\mathcal{F}}\) is the s-wave scattering length for the total spin \(\mathcal{F}=0,2\) channel.

The ground-state phase of this system is determined by the values of \(c_{1}\) and \(q_{z}\). (The phase diagram is given in, e.g., Ref. 27.) Below, we consider a condensate in the \(m=0\) state, i.e., a polar BEC. The polar BEC is the ground state of the system when \(q_{z}>\max(0,-2c_{1}n)\), and the system becomes dynamically unstable when \((c_{1},q_{z})\) is outside of this region.

The GP and BdG equations are obtained by the same manner as in the previous section. In the case of a polar BEC, the stationary solution \((\Psi_{1},\Psi_{0},\Psi_{-1})=(0,\Phi(x),0)\) satisfies \begin{equation} \mu \Phi = \left[-\frac{\hbar^{2}}{2M}\frac{d^{2}}{dx^{2}}+U(x)+c_{0}|\Phi|^{2}\right]\Phi, \end{equation} (33) where μ is the chemical potential. Assuming \(U(x)\rightarrow 0\) and \(|\Phi|^{2}\rightarrow n_{0}\) (= const.) at \(x\rightarrow\pm\infty\), we obtain \(\mu=c_{0}n_{0}\). We use this μ to scale the dimensionful variables, i.e., we rescale the energy, length, and time scales by μ, \(\hbar/\sqrt{M\mu}\), and \(\mu/\hbar\), respectively. By rewriting \(\Phi(x)=\sqrt{n_{0}}\phi(x)\), we reduce Eq. (33) to \begin{equation} (\mathcal{L}_{0}-1+|\phi|^{2})\phi = 0, \end{equation} (34) where \(\mathcal{L}_{0}\) is defined in Eq. (7). Equation (34) is identical to Eq. (6).

We introduce the interaction parameter β as \begin{equation} \beta \equiv \frac{c_{1}}{c_{0}}. \end{equation} (35) By rewriting the wave function as \(\Psi_{m}(x)=\sqrt{n_{0}}\psi_{m}(x)\), we obtain the time-dependent GP equation in the dimensionless form as \begin{align} i\frac{\partial\psi_{\pm 1}}{\partial t} &= [\mathcal{L}_{0}+1\pm\beta f_{z}+q_{z}]\psi_{\pm 1}+\frac{\beta}{\sqrt{2}}f_{\mp}\psi_{0}, \end{align} (36a) \begin{align} i\frac{\partial\psi_{0}}{\partial t} &= [\mathcal{L}_{0}+1]\psi_{0} +\frac{\beta}{\sqrt{2}}(f_{+}\psi_{+1}+f_{-}\psi_{-1}), \end{align} (36b) where \(f_{+}=f_{-}^{*}=\sqrt{2}(\psi_{1}^{*}\psi_{0}+\psi_{0}^{*}\psi_{-1})\) and \(f_{z}=|\psi_{1}|^{2}-|\psi_{-1}|^{2}\). By substituting \begin{equation} \begin{pmatrix} \psi_{+1}\\ \psi_{0}\\ \psi_{-1} \end{pmatrix}= e^{-it}\left[ \begin{pmatrix} 0\\ \phi\\ 0 \end{pmatrix}+\boldsymbol{u}e^{-iEt}+\boldsymbol{v}^{*}e^{+iE^{*}t}\right] \end{equation} (37) into Eq. (36) and linearizing the equation with respect to \(\boldsymbol{u}\) and \(\boldsymbol{v}\), with \(\boldsymbol{u}=(u_{+1},u_{0},u_{-1})^{\text{T}}\) and \(\boldsymbol{v}=(v_{+1},v_{0},v_{-1})^{\text{T}}\) being three-component spinors, we obtain the BdG equation for a polar BEC as \begin{equation} \begin{pmatrix} H_{0} + H_{1} &H_{2}\\ -H^{*}_{2} &-[H_{0} + H^{*}_{1}] \end{pmatrix}\begin{pmatrix} \boldsymbol{u}\\ \boldsymbol{v}\end{pmatrix}= E \begin{pmatrix} \boldsymbol{u}\\ \boldsymbol{v}\end{pmatrix}, \end{equation} (38) where \begin{align} H_{0} &= \begin{pmatrix} \mathcal{L}_{0}-1+q_{z} &0 &0\\ 0 &\mathcal{L}_{0}-1 &0\\ 0 &0 &\mathcal{L}_{0}-1+q_{z} \end{pmatrix}, \end{align} (39a) \begin{align} H_{1} &= |\phi|^{2} \begin{pmatrix} 1+\beta &0 &0\\ 0 &2 &0\\ 0 &0 &1+\beta \end{pmatrix}, \end{align} (39b) \begin{align} H_{2} &= \phi^{2} \begin{pmatrix} 0 &0 &\beta\\ 0 &1 &0\\ \beta &0 &0 \end{pmatrix}. \end{align} (39c) This \(6\times 6\) eigenvalue equation can be divided into three \(2\times 2\) equations. By defining the density-wave mode \(u^{\text{d}}=u_{0}\), \(v^{\text{d}}=v_{0}\) and spin-wave modes \(u_{\pm}^{\text{s}}=u_{\pm 1}\), \(v_{\pm}^{\text{s}}=v_{\mp 1}\), we reduce BdG Eq. (38) to \begin{align} \begin{pmatrix} \mathcal{L}_{0} -1 + 2|\phi|^{2} &\phi^{2}\\ -(\phi^{*})^{2} &-(\mathcal{L}_{0}-1 +2|\phi|^{2}) \end{pmatrix}\begin{pmatrix} u^{\text{d}}\\ v^{\text{d}} \end{pmatrix}&= E^{\text{d}} \begin{pmatrix} u^{\text{d}}\\ v^{\text{d}} \end{pmatrix}, \end{align} (40) \begin{align} \begin{pmatrix} \mathcal{L}_{0} -1 +q_{z}+ (1+\beta)|\phi|^{2} &\beta\phi^{2}\\ -\beta(\phi^{*})^{2} &-(\mathcal{L}_{0} -1 +q_{z}+(1+\beta)|\phi|^{2}) \end{pmatrix}\begin{pmatrix} u_{\pm}^{\text{s}}\\ v_{\pm}^{\text{s}} \end{pmatrix}&= E^{\text{s}} \begin{pmatrix} u_{\pm}^{\text{s}}\\ v_{\pm}^{\text{s}} \end{pmatrix}. \end{align} (41) Equation (40) has a zero-energy solution \(u^{\text{d}}=-(v^{\text{d}})^{*}=\phi\), indicating that the density-wave is the NG phonon associated with the spontaneous breaking of the U(1) gauge symmetry. On the other hand, the spin-wave mode is not the NG mode for \(q_{z}\neq 0\), since the SO(3) spin rotational symmetry is broken in the presence of the quadratic Zeeman energy. The previous works investigated the tunneling properties for BdG Eq. (41) with \(q_{z}=0\) and \(\beta>0\) and showed that the perfect transmission occurs in the low energy limit.^9,10)

Note that Eq. (41) at \(q_{z}=0\) is identical to Eq. (13). The role of \(q_{z}\) is that it effectively shifts the chemical potential from 1 to \(1-q_{z}\), opening an energy gap at \(k=0\) in a uniform system (see below). This is possible because the condensation occurs in a different internal state from the quasiparticles.

3.2 Bogoliubov spectrum in a uniform system

We analytically solve BdG Eq. (41) for \(U(x)=0\) and obtain propagating and damping/growing solutions. By substituting \((u^{\text{s}},v^{\text{s}})^{\text{T}}=(ue^{ikx},ve^{ikx})^{\text{T}}\) and \(\phi=1\) to Eq. (41), we obtain \begin{equation} \begin{pmatrix} \epsilon_{k} + q_{z} + \beta &\beta\\ -\beta &-(\epsilon_{k} + q_{z} + \beta) \end{pmatrix}\begin{pmatrix} u\\ v \end{pmatrix}= E^{\text{s}} \begin{pmatrix} u\\ v \end{pmatrix}, \end{equation} (42) which has the eigenvalue \begin{equation} E^{\text{s}} = \sqrt{(\epsilon_{k}+q_{z})(\epsilon_{k}+q_{z}+2\beta)}. \end{equation} (43) It follows that when \(q_{z}<\max(0,-2\beta)\), \(E^{\text{s}}\) becomes pure imaginary for a certain region of k and the system becomes dynamically unstable. The condition for the dynamical stability \(q_{z}>\max(0,-2\beta)\) agrees with the region for which the polar state is the ground state.

3.3 Tunneling properties of the spin-wave modes

We now investigate the tunneling properties of the spin-wave modes in a polar BEC. In the same way as in the previous section, for given energy E and momentum \(k_{\text{in}}\) of the incident wave, we construct the asymptotic form of the quasiparticle wave function at \(x\to\pm\infty\) and numerically solve BdG Eq. (41) by the finite element method. In the presence of the quadratic Zeeman effect, the energy dispersion is categorized into four types as shown in Fig. 5: (a) \(q_{z}>0\) and \(q_{z}+2\beta>0\), (b) \(q_{z}(q_{z}+2\beta)\leq 0\) and \(q_{z}+\beta>0\), (c) \(q_{z}(q_{z}+2\beta)\leq 0\) and \(q_{z}+\beta\geq 0\), and (d) \(q_{z}<0\) and \(q_{z}+2\beta<0\).

It is instructive to discuss the causality of the pure-imaginary-eigenvalue mode mentioned in Sect. 2.3 on the basis of the spectrum shown in Fig. 5(d). For \(0<E<\sqrt{q_{z}(q_{z}+2\beta)}\), \(k_{1}\) and \(k_{2}\) in Eqs. (45) and (46) are both real positive and satisfy \(k_{2}<k_{1}\). Since the group velocity at \(k=k_{1}\) (\(k_{2}\)) is positive (negative), \(e^{ik_{1}x}\) (\(e^{-ik_{2}x}\)) describes a right-going wave. When we continuously change E to pure imaginary, \(k_{1}\) and \(k_{2}\) are analytically continued to \(k_{3}\) and \(k_{4}\), respectively: Equations (48) and (49) are obtained from Eqs. (45) and (46), respectively, by extending the domain of E to a complex value. Hence, it is natural to regard \(e^{ik_{3}x}\) and \(e^{-ik_{4}x}\) having the same causality, obtaining Eqs. (28) and (29), which are the \(q_{z}\to 0\) limit of the present case.

We calculate the reflection and transmission probabilities for each case of (a)–(d) at various energies and momenta of the incident wave and find that the perfect transmission does not occur except for \(q_{z}=0\). Figures 6(a) and 6(b) show the behavior of the coefficients for the cases in Figs. 5(c) and 5(d), respectively, where we set the horizontal axis in the same manner as that in Fig. 3(b), i.e., \(\epsilon_{k_{\text{in}}}\) is zero at the left end of the figure and \(\epsilon_{k_{\text{in}}}\) increases as one goes right. We find that \(|T|^{2}<1\) for all regions of the figure. This result is consistent with the fact that the perfect transmission occurs when the quasiparticle wave function coincides with the condensate wave function: Since the quadratic Zeeman effect breaks the spin rotational symmetry, the spin-wave mode is not the NG mode of the system. On the other hand, one can see that the perfect reflection occurs at the maximum \(|{\mathop{\text{Im}}\nolimits E}|\) in both figures [at \(E=0.25i\) and \(0.15i\) for Figs. 6(a) and 6(b), respectively], indicating that this is a universal property of dynamically unstable modes.

Note that the transmission probability resonantly increases at \(E=0.2i\) in Fig. 6(a) and at \(E=0.2\) in Fig. 6(b). These points correspond to the incident energy \(E=\sqrt{q_{z}(q_{z}+2\beta)}\) with nonzero \(k_{\text{in}}\), which are depicted in Figs. 5(c) and 5(d) with filled circles. At these points, the momenta of the A and B terms in the asymptotic form become zero. Namely, the propagating modes for \(\mathop{\text{Im}}\nolimits E>0.2\) [\(\mathop{\text{Re}}\nolimits E<0.2\)] change into localized modes for \(\mathop{\text{Im}}\nolimits E<0.2\) [\(\mathop{\text{Re}}\nolimits E>0.2\)] in Fig. 5(c) [5(d)]. The increase in \(|T|^{2}\) at this point is understood as a resonance with these A and B terms. The peak value of the transmission probability depends on the barrier potential and decreases with increasing barrier potential.