Introduction

In recent years, the emergence of cloud-based experiment platforms has significantly expanded research methodologies in physics. This trend has been particularly prevalent in the area of quantum computing, with IBM Quantum (https://quantum.ibm.com), Amazon Braket (https://aws.amazon.com/braket), and Google Quantum AI (https://quantumai.google) serving as prominent examples. These services allow researchers to conduct experiments remotely without the need for direct access to expensive hardware. These cloud platforms facilitate the exploration of new phenomena and the rapid testing of theoretical concepts.

Oqtant is one such service, a cloud-accessible platform for a Bose-Einstein condensate (BEC) of ultracold neutral atomic gas, which is provided by a company in the USA, namely Infleqtion (https://oqtant-docs.infleqtion.com and https://www.infleqtion.com/oqtant). In the Oqtant platform, one can create a BEC of 87Rb atoms trapped in a cigar-shaped potential and utilize many useful functionalities. In particular, one can manipulate in real time repulsive potentials for atoms generated by blue-detuned lasers and observe the time evolution of the density distribution of atoms. This capability allows even theoretical researchers to experimentally explore the dynamics of a weakly interacting dilute Bose gas.

As an interesting phenomenon in a weakly interacting BEC, anomalous tunneling of Bogoliubov excitations, which are elementary excitations of the BEC, was theoretically predicted in the 2000s1,2. When one considers a one-dimensional (1D) scattering problem of Bogoliubov excitations across a repulsive potential barrier illustrated in Fig. 1a, the tunneling probability is higher for lower excitation energies. In the zero-energy limit, it even exhibits perfect transmission. This behavior is in stark contrast with the tunneling of a single particle, so it is called anomalous tunneling. This phenomenon is related to superfluidity3,4,5 and is conjectured to be a universal property of Nambu-Goldstone modes associated with the spontaneous breaking of a continuous symmetry6,7,8,9. Despite such extensive theoretical interest, experimental observation of anomalous tunneling has not been realized so far.

Fig. 1: Transmission coefficient of BECs in the tunneling problem.
Fig. 1: Transmission coefficient of BECs in the tunneling problem.
Full size image

a Schematic picture of the tunneling problem. The blue solid and green dashed curves represent the condensate density ∣ψ0∣2 and the Gaussian barrier. b Calculated transmission coefficient \(| {{{\mathcal{T}}}}{| }^{2}\) as a function of the energy of the incident wave ε for some barrier heights V0. We also show the halfwidths ΔεL and ΔεH for V0/μ = 0.7.

In this paper, using the Oqtant platform, we experimentally observe some signatures of anomalous tunneling in collective excitation modes of BECs in double-well potentials consisting of a trap potential and a repulsive potential barrier at the trap center. Previous theoretical work has shown that anomalous tunneling is closely related to how the lowest frequency of collective modes of BECs in the double-well potentials depends on the barrier height10. To provide theoretical references more specialized to the experimental setup of Oqtant, we theoretically extend this analysis to several higher-frequency modes. We show that anomalous tunneling is also relevant to the barrier-height dependence of the higher mode frequencies. In particular, in addition to the lowest frequency of the collective mode, the second and third lowest frequencies are suitable for measurement within the specifications of Oqtant. We also report experimental results obtained using Oqtant. In our experiment setup, we apply a repulsive potential to excite collective modes and then observe the time evolution of the density distribution of the atom cloud after removing the potential to evaluate the frequencies of collective modes. The above protocol is performed in the presence of a potential barrier for various barrier heights. We discuss the relation between the observed barrier-height dependence of collective mode frequencies and the anomalous tunneling effect.

Results

Theoretical analysis

In order to conduct numerical simulations with the physical parameters of BECs of the Oqtant platform, let us first introduce the specifications of Oqtant. Rubidium 87Rb atoms in the hyperfine state \(\left|F=2,{m}_{F}=2\right\rangle\) are trapped by the cigar-shaped harmonic potential \({V}_{{{{\rm{trap}}}}}=m[{\omega }_{\perp }^{2}({x}^{2}+{y}^{2})+{\omega }_{z}^{2}{z}^{2}]/2\), where m is the mass of Rubidium atoms. The transverse and axial trap frequencies are, respectively, given by ω⊥ = 2π × 400 Hz and ωz = 2π × 42 Hz. The value of ω⊥ is taken from the Oqtant manual, whereas that of ωz is determined from dipole-motion measurements. BECs are produced via evaporative cooling. The typical temperature of BECs is about 90 nK ( ≈ 0.4Tc, where Tc( ≈ 225 nK) is the critical temperature of a noninteracting Bose gas), and the condensate atom number Nc is about 8300. For more details, see Supplementary Note 1. We can measure the atom density distribution until 100 ms after BEC preparation.

Since ω⊥ ≫ ωz, the condensate is expected to show quasi-1D behavior. Although actual experiments are at relatively high temperatures ( ≈ 0.4Tc), we ignore finite temperature effects in our theoretical analysis. Thus, we use the effective quasi-1D Gross-Pitaevskii (GP) equation 11,12,13;

$${{{\rm{i}}}}\hslash \frac{\partial \psi (z,t)}{\partial t}=\left(-\frac{{\hslash }^{2}}{2m}\frac{{\partial }^{2}}{\partial {z}^{2}}+{V}_{{{{\rm{ext}}}}}(z)+\hslash {\omega }_{\perp }\sqrt{1+4a| \psi (z,t){| }^{2}}\right)\psi (z,t),$$
(1)

where Vext is an external potential, a ( ≈ 5.29 nm14) is the s-wave scattering length of 87Rb, and ψ(z, t) is the 1D condensate wave function, which is normalized as ∫ dz∣ψ(z, t)∣2 = Nc. Vext(z) represents the external potential term along the z direction, including the trap potential and the other potentials that we input.

To analyze anomalous tunneling of excitations through a barrier and collective excitations in a double-well potential, we assume that the condensate is weakly perturbed. We put the solution of Eq. (1) as \(\psi (z,t)=\left({\psi }_{0}(z)+u(z){{{{\rm{e}}}}}^{-{{{\rm{i}}}}\varepsilon t/\hslash }-{v}^{* }(z){{{{\rm{e}}}}}^{{{{\rm{i}}}}\varepsilon t/\hslash }\right){{{{\rm{e}}}}}^{-{{{\rm{i}}}}\mu t/\hslash }\), where ψ0(z) is the static solution of Eq. (1), (u(z), v(z)) the small fluctuations, μ the chemical potential, and ε the excitation energy. The static solution ψ0(z) satisfies the time-independent GP equation

$$\left({{{\mathcal{L}}}}+\hslash {\omega }_{\perp }\sqrt{1+4a| {\psi }_{0}(z){| }^{2}}\right){\psi }_{0}(z)=0,$$
(2)

where \({{{\mathcal{L}}}}=-\frac{{\hslash }^{2}}{2m}\frac{{\partial }^{2}}{\partial {z}^{2}}-\mu +{V}_{{{{\rm{ext}}}}}(z)\), and the fluctuations u(z) and v(z) are determined by the linearized GP equation (Bogoliubov equation)13

$$\left(\begin{array}{cc}{{{\mathcal{L}}}}+f(z) & -g(z)\\ -g(z) & {{{\mathcal{L}}}}+f(z)\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)=\varepsilon \left(\begin{array}{c}u\\ -v\end{array}\right),$$
(3)

where \(f(z)=| g(z)| +\hslash {\omega }_{\perp }\sqrt{1+4a| {\psi }_{0}(z){| }^{2}}\) and \(g(z)=2\hslash {\omega }_{\perp }a{\psi }_{0}^{2}(z)/\sqrt{1+4a| {\psi }_{0}(z){| }^{2}}\).

To check how anomalous tunneling occurs in the current quasi-1D situation, we consider a tunneling problem as schematically shown in Fig. 1a. We set a Gaussian potential; \({V}_{{{{\rm{ext}}}}}={V}_{0}\exp \left(-{z}^{2}/2{\sigma }^{2}\right)\). The width σ is set to be 0.5 μm, which is the smallest one available in Oqtant. Because we are considering the infinite system, we should specify the chemical potential μ rather than Nc. To be consistent with the chemical potential in experiments, we calculate the chemical potential by solving Eq. (2) with the trap potential \({V}_{{{{\rm{ext}}}}}=m{\omega }_{z}^{2}{z}^{2}/2\) and Nc = 8300, and use the resulting value μ/(ℏωz) ≈ 35.03. Using this chemical potential, the healing length \(\xi =\sqrt{{\hslash }^{2}/2m\mu }\) is estimated by ξ ≈ 0.2 μm. Thus, the width σ is expressed as σ ≈ 2.5ξ.

Under the above setup, we first determine the static condensate ψ0 from Eq. (2). Next, let us consider the situation where the plane waves of fluctuations u and v come from z → − ∞ and then are scattered by the Gaussian potential. The tunneling and reflection properties are encoded in the boundary condition as

$$\begin{array}{r}\left(\begin{array}{r}u(z)\\ v(z)\end{array}\right)=\left\{\begin{array}{ll}\left(\begin{array}{r}\widetilde{u}\\ \widetilde{v}\end{array}\right){{{{\rm{e}}}}}^{{{{\rm{i}}}}kz}+{{{\mathcal{R}}}}\left(\begin{array}{r}\widetilde{u}\\ \widetilde{v}\end{array}\right){{{{\rm{e}}}}}^{-{{{\rm{i}}}}kz} & z\to -\infty \\ {{{\mathcal{T}}}}\left(\begin{array}{r}\widetilde{u}\\ \widetilde{v}\end{array}\right){{{{\rm{e}}}}}^{{{{\rm{i}}}}kz}\hfill & z\to \infty \end{array}\right.,\end{array}$$
(4)

where the wave number k of the incident wave satisfies the dispersion relation of the Bogoliubov equation in the uniform regime; \({\varepsilon }^{2}=\frac{{\hslash }^{2}{k}^{2}}{2m}\left(\frac{{\hslash }^{2}{k}^{2}}{2m}+\frac{4\hslash {\omega }_{\perp }a{n}_{0}}{\sqrt{1+4a{n}_{0}}}\right)\) with the uniform condensate density \({n}_{0}={\lim}_{| z| \to \infty }| {\psi }_{0}(z){| }^{2}\). \(\widetilde{u}\) and \(\widetilde{v}\) are the solutions of the Bogoliubov equation in the uniform regime. \({{{\mathcal{T}}}}\) and \({{{\mathcal{R}}}}\) in Eq. (4) represent the transmission and reflection amplitudes satisfying \(| {{{\mathcal{T}}}}{| }^{2}+| {{{\mathcal{R}}}}{| }^{2}=1\). We numerically solve Eq. (3) under the boundary condition Eq. (4) to determine \({{{\mathcal{T}}}}\) for given ε. The method to numerically calculate \({{{\mathcal{T}}}}\) and \({{{\mathcal{R}}}}\), as well as ψ0(x), is summarized in the “Details of numerical analysis of tunneling coefficient” subsection in the Methods.

Before showing numerical results, let us briefly discuss the use of the term “tunneling” in our setup. In particular, the case where V0 < μ is confusing, because one may think that the effect of the potential only makes the condensate have a dent structure and therefore excitations are not significantly affected. To capture the effect of the potential on fluctuations u and v, we subtract the free part, which is the same as the Bogoliubov equation far away from the potential, from the original one:

$$\begin{array}{r}\left(\begin{array}{cc}{{{{\mathcal{L}}}}}_{0}+{f}_{0} & -{g}_{0}\\ -{g}_{0} & {{{{\mathcal{L}}}}}_{0}+{f}_{0}\end{array}\right)\left(\begin{array}{r}u\\ v\end{array}\right)+\left(\begin{array}{cc}{V}_{{{{\rm{ext}}}}}+\delta f(z) & -\delta g(z)\\ -\delta g(z) & {V}_{{{{\rm{ext}}}}}+\delta f(z)\end{array}\right)\left(\begin{array}{r}u\\ v\end{array}\right)=\varepsilon \left(\begin{array}{c}u\\ -v\end{array}\right),\end{array}$$
(5)

where \({{{{\mathcal{L}}}}}_{0}=-\frac{{\hslash }^{2}}{2m}\frac{{d}^{2}}{d{z}^{2}}-\mu\), \({g}_{0}=2\hslash {\omega }_{\perp }a{n}_{0}/\sqrt{1+4a{n}_{0}}\), \({f}_{0}={g}_{0}+\hslash {\omega }_{\perp } \sqrt{1+4a{n}_{0}}\), δf(z) = f(z) − f0, and δg(z) = g(z) − g0. The first term is the free part, and the second term can be regarded as the potential for fluctuations u and v. This expression suggests that, as for the scattering of fluctuations u and v, the relation of magnitude between μ and V0 is irrelevant. The bare potential Vext and the difference of condensate density together make the potential for fluctuations. In this sense, even when V0 < μ, fluctuations certainly feel the potential and then “tunnel” the potential even at low energies.

We numerically calculate the transmission coefficient \(| {{{\mathcal{T}}}}{| }^{2}\) for various barrier height V0 in the range 0 < V0/μ ≲ 2.28 (0 < V0/ℏωz ≤ 80) and show the energy dependence of \(| {{{\mathcal{T}}}}{| }^{2}\) in Fig. 1b for several barrier heights. As the energy is lowered from the high-energy side, the transmission coefficient first decreases but then increases in the low-energy region, finally approaching unity in the zero-energy limit. This is nothing but the anomalous tunneling effect.

A higher barrier potential V0 induces a more abrupt increase in the transmission coefficient, as seen in Fig. 1b. For convenience in discussing the relation between anomalous tunneling and collective excitation frequencies, we define the halfwidth Δε of the transmission coefficient, i.e., the energy at which the transmission coefficient equals 0.5. We note that there exist two halfwidths, ΔεL and ΔεH, as seen in Fig. 1b. The low-energy halfwidth ΔεL is due to anomalous tunneling, but the high-energy one ΔεH simply comes from the well-known fact that high-energy excitations are less sensitive to a potential. These halfwidths are plotted as functions of V0 in Fig. 2b, with V0 ranging from 0 to 80ℏωz ( ≈ 2.28μ). Excitations with energy ε < ΔεL or ε > ΔεH are not significantly affected by a barrier.

Fig. 2: Collective excitation frequencies of BECs in a double-well potential.
Fig. 2: Collective excitation frequencies of BECs in a double-well potential.
Full size image

a Schematic drawing of the situation in examining collective modes. The blue solid and red dash-dotted curves denote the condensate and the collective excitation, respectively. The green dashed line shows the potential energy, which is composed of the Gaussian barrier and the trap potential. b Calculated collective mode frequencies ε as functions of the barrier height V0. The solid and dashed lines show odd-parity and even-parity mode frequencies, respectively. The black dash-dotted line shows the halfwidth due to anomalous tunneling ΔεL(V0). The black dotted line is the halfwidth of high-energy excitations ΔεH(V0). We also put blue dots at (Vm, ε2i+1(Vm)) to indicate the characteristic barrier height Vm for the merging of two neighboring modes. c Magnified view of the low-energy region of (b). Excitations in this region are examined experimentally in the next section.

Next, we investigate the impact of anomalous tunneling on collective excitations under the Oqtant setup. We set the potential as \({V}_{{{{\rm{ext}}}}}=m{\omega }_{z}^{2}{z}^{2}/2+{V}_{0}\exp \left(-{z}^{2}/2{\sigma }^{2}\right)\) and solve Eqs. (2) and (3) to obtain the eigen frequencies ε/ℏ of collective excitations. The situation is schematically depicted in Fig. 2a.

We show calculated collective mode frequencies as functions of the barrier height in Fig. 2b. We see that the two neighboring frequencies ε2i(V0) and ε2i+1(V0) (i = 0, 1, 2, … ) tend to merge as V0 is increased. Notice that the lowest frequency ε1 tends to merge with the zero mode one ε0 = 0. This tendency is easily understood: When the barrier height is high enough, the condensate can be regarded as two independent BECs. In this case, the in-phase and out-of-phase oscillations of two BECs should have the same frequency. The even collective modes ε2i (odd ones ε2i+1) are indeed even (odd) parity modes and then converge to out-of-phase (in-phase) oscillations at V0 ≫ μ.

We introduce a characteristic barrier height Vm for the merging of two collective modes. We define Vm at which the difference of frequencies between two neighboring modes is half of that in the absence of the barrier; the values of Vm satisfy [ε2i+1(Vm) − ε2i(Vm)]/[ε2i+1(0) − ε2i(0)] = 1/2, which are shown by the dots in Fig. 2b. We also show the halfwidth Δε(V0) in Fig. 2b and find that Vm and Δε(V0) are well close to each other.

The coincidence between Δε(V0) and Vm suggests an interpretation of the merging behavior. In the region left of Δε(V0) in Fig. 2b, the transmission coefficient is greater than 1/2, so the effects of the barrier on the collective excitations are expected to be insignificant. Conversely, in the region right of Δε(V0), the barrier strongly affects the collective excitations, making two collective modes merge in the high barrier height limit. In this way, the merging behavior is also characterized by Δε(V0), which explains why Vm is close to Δε(V0).

Now, we discuss how to capture anomalous tunneling from measurements of collective mode frequencies. Let us first consider the case in the absence of the trap potential. As shown in Figs. 1b and 2b, for a given barrier height V0, two halfwidths appear. Of the two halfwidths, the signatures of anomalous tunneling manifest themselves in the low-energy branch: a larger barrier height corresponds to a more abrupt increase in the transmission coefficient, and hence to a smaller halfwidth, as shown by the dash-dotted line (ΔεL) in Fig. 2b and c. Therefore, the observation of the low-energy branch ΔεL would provide evidence for anomalous tunneling.

In this work, instead of directly measuring ΔεL in the trap-free system, we adopt the measurement of the merging barrier height Vm in the trapped system as a proxy, because Δε and Vm exhibit very similar behavior. In the low-energy region shown in Fig. 2c, Vm follows the same trend as ΔεL: for lower energy collective modes, Vm appears at the large barrier height V0. Thus, the observation of the merging barrier heights Vm for the lowest mode i = 0 and the second and third-lowest modes i = 1, together with the relation Vm,i=0 > Vm,i=1 would be indirect evidence of anomalous tunneling. The typical time scales of the first and second (third) lowest collective excitations are, respectively, T1 = 2π/ωz ≈ 24 ms and T2(3) ≈ 15 ms, which are feasible in Oqtant.

Experimental results

To observe low-energy collective excitations by using Oqtant, we perform experiments shown in Fig. 3 after preparing a BEC in the cigar-shaped trap mentioned above. The prepared BEC cloud is in dipole motion. To reduce this unwanted oscillation, we add a linear potential from 2 ms to 3 ms (0 ms corresponds to the end of the BEC preparation). We ramp up the barrier at the center of the trap linearly in time up to V0 over 8 ms to create a double-well potential. Notice that a BEC in a double-well potential has also been created in previous experiments15,16,17,18,19. While a few low-energy collective excitations have been observed in Refs. 16,19, their connection to anomalous tunneling has never been addressed. From 15 ms to 17 ms, we add a perturbation potential to excite collective motions. We use three types of potentials depending on the height of the barrier and which collective mode we intend to excite. After removing the perturbation potential, we measure the atom cloud distribution in time. The measurements were performed until 80 ms because the operation is limited to that time in Oqtant. (It was extended to 100 ms at some point. Some experiments were conducted up to 100 ms.) For more details of the choice of these parameters, see Supplementary Note 2.

Fig. 3: Experimental protocol after BEC preparation.
Fig. 3: Experimental protocol after BEC preparation.
Full size image

a Barrier schedule. We use three potentials: Vcanc, Vload, and Vpert. The horizontal axis represents the elapsed time after BEC preparation. The vertical axis shows the potential heights. b Potential shape. Vcanc is used for the cancellation of unwanted initial motion. Vload is the barrier at the center of the trap and linearly grows to V0 over 8 ms. Vpert is used to excite collective motions. We actually use three types of perturbation potentials and here show only one of them; for more details, see Supplementary Note 2.

Figure 4a shows an example of observed atom cloud distributions after removing background noise. To extract the frequency of collective excitations, we calculate the center of mass position (CoM) and the standard deviation (SD), as shown in Fig. 4b. We mainly use CoM (SD) for the first and third (second) mode frequency estimation because it is expected to be sensitive to odd (even) parity modes.

Fig. 4: Observed atom-cloud dynamics.
Fig. 4: Observed atom-cloud dynamics.
Full size image

a One example of an observed atom cloud distribution at an elapsed time 18 ms in the absence of a barrier. The color bar shows the optical depth corresponding to the atom density in an arbitrary unit. b Center of mass (CoM) and standard deviation (SD) of atom clouds along the z direction. The error bars represent the standard deviations of CoM and SD across measurements. We also show the corresponding power spectra (PS) in (c). In the power spectra, we remove the offsets of the original signals. The red squares with error bars show the estimated peaks of PS. The vertical dashed lines show the theoretical predictions for the lowest, second-lowest, and third-lowest collective mode frequencies in ascending order of frequency. The error bars represent the standard deviation of the peak positions obtained from noise-added data based on the measured standard deviation. For details, see the “Peak analysis of Fourier spectra of experimental data” subsection in the Methods.

We determine frequencies and their errors in two ways. The first method is the fitting of the observed time evolution of CoM and SD with the function

$$f(t)={A}_{1}{{{{\rm{e}}}}}^{-{B}_{1}t}\sin ({C}_{1}t+{D}_{1})+{A}_{2}{{{{\rm{e}}}}}^{-{B}_{2}t}\sin ({C}_{2}t+{D}_{2})+E.$$
(6)

We take the parameter C1(2) exhibiting larger oscillation amplitude Ai as an estimate of the lowest mode. The C1 or C2 value close to the theoretical prediction is taken as an estimated frequency of the second and third modes. The error is estimated from the standard deviation of C1(2).

The second method is a peak analysis in the Fourier spectra. Figure 4c shows the Fourier spectra of CoM and SD (blue curves), and their peaks (red squares with error bars). We also plot the theoretical predictions for the frequencies of the lowest three collective excitations as the vertical lines in Fig. 4c. We use the frequency and standard deviation of the dominant peak as estimates of the lowest mode. In addition, we use those of the peaks whose frequencies are close to the theoretical predictions as estimates of the second and third modes.

The estimated frequencies are summarized in Fig. 5. As for the second and third collective excitations, the measured frequencies agree well with the theoretical predictions, i.e., with increasing the barrier height, the two frequencies approach each other. We also find that the merging behavior becomes significant when the excitation energy ε is larger than the halfwidth of anomalous tunneling Δε.

Fig. 5: Estimated collective mode frequencies.
Fig. 5: Estimated collective mode frequencies.
Full size image

Red circles, blue triangles, and magenta squares denote the estimated lowest, second-lowest, and third-lowest frequencies, respectively. Unfilled (filled) symbols represent frequencies estimated from fitting (Fourier) analysis; their x-axis values are shifted by  − 1.5 ( + 1.5) for better readability. The error bars represent the uncertainties of the fitted parameters and the standard deviation of the peak positions, respectively. The data of the lowest mode from fitting at V0/ℏωz ≈ 70 and those of the third-lowest mode from fitting at V0/ℏωz ≈ 60 and 70 are absent because we cannot determine them. Solid, dashed, and black dash-dotted lines are the same as those in Fig. 2c. The purple cross and diamond show the collective mode frequency merging points for the lowest mode \({V}_{{{{\rm{m}}}},i=0}^{\exp }\) and the second-third lowest modes \({V}_{{{{\rm{m}}}},i=1}^{\exp }\), determined from the experimental data, respectively.

As for the lowest-frequency collective excitations, experimental results and theoretical predictions show the same qualitative tendency: As the barrier height increases, both decrease. Although they quantitatively show a large deviation compared to the second and third collective excitation case, we obtain that the lowest mode merging barrier height (\({V}_{{{{\rm{m}}}},i=0}^{\exp } \sim 3.58\,{{{\rm{kHz}}}}\)) is larger than the second-third mode one (\({V}_{{{{\rm{m}}}},i=1}^{\exp } \sim 1.01\,{{{\rm{kHz}}}}\)), which is also a signature of anomalous tunneling.

Let us briefly discuss the discrepancy between theory and experiment of the lowest-frequency mode. Finite-temperature effects may be responsible for this. However, finite-temperature analysis in a cigar-shaped and double-well trap20 shows that finite-temperature effects reduce the frequency of the lowest-frequency mode. The tendency is the opposite since our measured frequency data are larger than the zero-temperature theoretical data.

Another possible origin of the discrepancy is a three-dimensional (3D) effect. Since the radial confinement is not so strong (aspect ratio ω⊥/ωz ≈ 10), the lowest-mode discrepancy may arise from 3D leakage or transverse modes. To examine this possibility, we performed 3D GP simulations using an existing 3D GP simulation package21,22. The details of the simulation are summarized in the “Details of 3D GP simulation” subsection in the Methods.

In this simulation, the ground state Ψ0(x, y, z) is obtained by imaginary-time evolution in the Oqtant trap potential setup. To break the axial symmetry, we add small fluctuations with amplitude \(\alpha \sqrt{{N}_{{{{\rm{c}}}}}/{a}_{z}^{3}}\) to the condensate, where α is the parameter which controls the magnitude of fluctuations and \({a}_{z}=\sqrt{\hslash /m{\omega }_{z}}\) is the harmonic oscillator length along the z direction (for details, see the “Details of 3D GP simulation” subsection in the Methods). We take the average of estimated frequencies over 10 random fluctuation configurations.

The estimated frequencies from the 3D GP simulations are shown in Fig. 6 for the cases of V0 = 1.5 kHz and V0 = 2 kHz. We confirm that the frequencies predicted by the quasi-1D GP equation and those estimated from the 3D GP simulations agree well with each other without the initial condensate fluctuations. We also find that, even when 3D effects are included via initial fluctuations in the condensate, the estimated frequency is not sensitive to the fluctuation amplitude. In addition, as the magnitude of fluctuations increases, the estimated frequency tends to decrease, so that the discrepancy from the experimental results increases. These observations suggest that the 3D effects on the frequency are not strong and would not be the origin of the discrepancy.

Fig. 6: Estimated frequency of the lowest-frequency mode from 3D GP simulations.
Fig. 6: Estimated frequency of the lowest-frequency mode from 3D GP simulations.
Full size image

The horizontal axis denotes the magnitude of the fluctuations in the condensate. Red circles and blue triangles show the estimated frequencies for V0 = 1.5 kHz from fitting and Fourier analyses, respectively. Likewise, red squares and blue diamonds represent estimated frequencies for V0 = 2.0 kHz from fitting and Fourier analyses, respectively. The x-axis values from fitting (Fourier) analysis are shifted by  − 0.002( + 0.002) for better readability. We also show the standard deviations for random-fluctuation configurations as error bars, but they are smaller than the symbol sizes for the Fourier-analysis results. The green solid and magenta dashed lines represent the estimated frequencies from Fourier analyses of the experiment for V0 = 1.5 kHz and V0 = 2.0 kHz, respectively. The shaded regions with the same colors denote the corresponding standard deviations. We also show the predictions from the quasi-1D GP equation by black dotted (V0 = 2.0 kHz) and purple dash-dotted (V0 = 1.5 kHz) lines.

To address the discrepancy between theoretical and experimental results for the frequency of the lowest-frequency mode, more sophisticated frequency analysis, e.g., using machine learning23, and a careful setup suitable for the detection of the lowest-frequency mode would be helpful, but we leave these as future problems.

Conclusion

We theoretically and experimentally studied the collective excitation modes of a BEC in a double-well potential using the modified GP equation for a quasi-1D geometry and the cloud-accessible BEC platform named Oqtant. In the theory part, we showed that two neighboring collective mode frequencies merge as the barrier height increases, and that the merging is related to anomalous tunneling in the low-energy region. In the experiment part, we observed the merging behavior of collective excitation frequencies by using Oqtant. The observed frequencies of the second and third collective modes agree well with the theoretical predictions. As for the lowest collective mode, the measured frequency considerably deviates from the theoretical prediction, but they show the same qualitative tendency: as the barrier increases, both frequencies decrease. Moreover, we confirmed that the experimentally determined merging points satisfy \({V}_{{{{\rm{m}}}},i=1}^{\exp } < {V}_{{{{\rm{m}}}},i=0}^{\exp }\). From these observations, we conclude that we captured indirect evidence of anomalous tunneling in our Oqtant experiment.

Although we obtained indirect evidence of anomalous tunneling, it is beneficial to observe anomalous tunneling via the transmission coefficient in the setup of a standard 1D scattering problem depicted in Fig. 1a. Since direct measurement setups are within reach of modern experimental technology in cold-atom systems, we hope that our results will stimulate future experiments aiming at direct observation of anomalous tunneling. Performing other interesting BEC experiments using Oqtant would also be exciting future work. In doing so, since the Oqtant setup is at finite temperatures, simulations based on the Zaremba-Nikuni-Griffin equation 24 would be useful. Various cloud experiment platforms have recently become available, such as Aquila (a Rydberg atom system provided by QuEra)25. The use of these platforms is another interesting future direction (the use of Aquila has already been reported26,27,28).

Methods

Details of numerical analysis of tunneling coefficient

For the numerical evaluation of the tunneling coefficient, we first introduce the nondimensional variables: \(\bar{z}=z/\xi\), \(\bar{\varepsilon }=\varepsilon /\mu\), \({\bar{V}}_{{{{\rm{ext}}}}}={V}_{{{{\rm{ext}}}}}/\mu\), \(\bar{{a}}=a/\xi\), \({\bar{\omega }}_{\perp }=\hslash {\omega }_{\perp }/\mu\), \({\bar{\psi }}_{0}={\psi }_{0}/\sqrt{\xi }\), \(\bar{{u}}=u/\sqrt{\xi }\), and \(\bar{v}=v/\sqrt{\xi }\), where \(\xi =\sqrt{{\hslash }^{2}/2m\mu }\) is the healing length. Using these nondimensional variables, the time-independent GP equation and Bogoliubov equation are expressed as

$$\left(\bar{{{{\mathcal{L}}}}}+{\bar{\omega }}_{\perp }\sqrt{1+4\bar{{a}}| {\bar{\psi }}_{0}(\bar{z}){| }^{2}}\right){\bar{\psi }}_{0}(\bar{z})=0,$$
(7)
$$\begin{array}{l}\left(\begin{array}{cc}\bar{{{{\mathcal{L}}}}}+\bar{f}(\bar{z}) & - {\bar{g}} (\bar{z})\\ - {\bar{g}} (\bar{z}) & \bar{{{{\mathcal{L}}}}}+\bar{f}(\bar{z})\end{array}\right)\left(\begin{array}{l} \bar{u} \\ \bar{v}\end{array}\right)=\bar{\varepsilon }\left(\begin{array}{c}\bar{{u}}\\ -\bar{v}\end{array}\right),\end{array}$$
(8)

where \(\bar{{{{\mathcal{L}}}}}=-\frac{{\partial }^{2}}{\partial {\bar{\bar{z}}}^{2}}-1+{\bar{V}}_{{{{\rm{ext}}}}}(\bar{z})\), \(\bar{f}(\bar{z})=| {\bar{g}} (\bar{z})| +{\bar{\omega }}_{\perp }\sqrt{1+4\bar{{a}}| {\bar{\psi }}_{0}(\bar{z}){| }^{2}}\) and \({\bar{g}} (\bar{z})=2{\bar{\omega }}_{\perp }\bar{{a}}{\bar{\psi }}_{0}^{2}(z)/\sqrt{1+4\bar{{a}}| {\bar{\psi }}_{0}(\bar{z}){| }^{2}}\). Hereafter, we omit the bars of nondimensional variables for simplicity.

We discretize the space at the interval Δz, z → zi ≡ Δz × i − L/2 (i = 0, 1, …, N − 1), where L is the length of the simulation domain. Δz is set such that zN−1 = L/2. L is chosen to be large enough that the effects of the potential can be ignored at z = ± L/2. We use L = 30 and N = 2001. We denote ψ0(zi) ≡ ψ0,i, and ui, vi, Vext,i, fi, gi are defined similarly. The second derivative term is approximated by the finite difference: ∂2ψ0(z)/∂z2 = (ψ0,i+1 − 2ψ0,i + ψ0,i−1)/(Δz)2.

We determine ψ0,i using the imaginary-time evolution instead of directly solving the time-independent GP equation. The Bogoliubov equation (3) can be expressed by linear equations:

$$\begin{array}{rcl} \displaystyle -\frac{{u}_{i-1}-2{u}_{i}+{u}_{i+1}}{{(\Delta z)}^{2}}+({V}_{{{{\rm{ext}}}},i}-1+{f}_{i}-\varepsilon ){u}_{i}-{g}_{i}{v}_{i} & = & 0,\\ \displaystyle -{g}_{i}{u}_{i}-\frac{{v}_{i-1}-2{v}_{i}+{v}_{i+1}}{{(\Delta z)}^{2}}+({V}_{{{{\rm{ext}}}},i}-1+{f}_{i}+\varepsilon ){v}_{i} & = & 0,\end{array}$$
(9)

for i = 0, 1, …, N − 1. The number of equations is 2N, but the number of unknown variables is 2N + 4 because of the existence of u−1, v−1, uN, and vN. We should eliminate these four variables by imposing the boundary condition. The boundary condition (4) implies

$$\left(\begin{array}{r}{u}_{0}\\ {v}_{0}\end{array}\right)= (1+{{{\mathcal{R}}}})\left(\begin{array}{r}1\\ \widetilde{v}\end{array}\right),\,\left(\begin{array}{r}{u}_{-1}\\ {v}_{-1}\end{array}\right)=({{{{\rm{e}}}}}^{-{{{\rm{i}}}}k\Delta z}+{{{\mathcal{R}}}}{{{{\rm{e}}}}}^{{{{\rm{i}}}}k\Delta z})\left(\begin{array}{r}1\\ \widetilde{v}\end{array}\right),\\ \, \left(\begin{array}{r}{u}_{N-1}\\ {v}_{N-1}\end{array}\right)={{{\mathcal{T}}}}\left(\begin{array}{r}1\\ \widetilde{v}\end{array}\right),\,\left(\begin{array}{r}{u}_{N}\\ {v}_{N}\end{array}\right)={{{\mathcal{T}}}}{{{{\rm{e}}}}}^{{{{\rm{i}}}}k\Delta z}\left(\begin{array}{r}1\\ \widetilde{v}\end{array}\right),$$
(10)

where \(\widetilde{v}\) satisfies the Bogoliubov equation in the uniform regime with \(\widetilde{u}=1\):

$$\begin{array}{r}\left(\begin{array}{cc}{k}^{2}-1+{f}_{0} & -{g}_{0}\\ -{g}_{0} & {k}^{2}-1+{f}_{0}\end{array}\right)\left(\begin{array}{r}1\\ \widetilde{v}\end{array}\right)=\varepsilon \left(\begin{array}{c}1\\ -\widetilde{v}\end{array}\right).\end{array}$$
(11)

We can incorporate these conditions by substituting

$$\begin{array}{r}\left(\begin{array}{r}{u}_{-1}\\ {v}_{-1}\end{array}\right)={{{{\rm{e}}}}}^{{{{\rm{i}}}}k\Delta z}\left(\begin{array}{r}{u}_{0}\\ {v}_{0}\end{array}\right)-2{{{\rm{i}}}}\sin (k\Delta z)\left(\begin{array}{r}1\\ \widetilde{v}\end{array}\right),\,\left(\begin{array}{r}{u}_{N}\\ {v}_{N}\end{array}\right)={{{{\rm{e}}}}}^{{{{\rm{i}}}}k\Delta z}\left(\begin{array}{r}{u}_{N-1}\\ {v}_{N-1}\end{array}\right).\end{array}$$
(12)

Thus, we can eliminate four unknown variables, and then the problem is reduced to solving the linear equation. We can determine \({{{\mathcal{T}}}}\) and \({{{\mathcal{R}}}}\) simply as \({{{\mathcal{T}}}}={u}_{N-1}\) and \({{{\mathcal{R}}}}={u}_{0}-1\). A similar technique to numerically calculate the transmission coefficients (more generally, scattering wave functions under various boundary conditions) for the Schrödinger equation can be found in the literature29.

Peak analysis of Fourier spectra of experimental data

In analyzing peaks of the Fourier spectra of the center of mass (CoM) and standard deviation (SD) of the observed atom cloud distributions, statistical errors should be treated carefully. To account for the influence of statistical errors on peaks, we calculate the averages and standard deviations of CoM and SD for each time and randomly sample them according to Gaussian distributions. We generate 1024 sampled time-series datasets and calculate their Fourier spectra. We also apply zero-padding to the input signal before performing the Fourier transform to improve the frequency resolution. For each peak in the Fourier spectra, we evaluate the sample averages of peak frequencies and their standard deviations. Figure 4c in the main text shows the sample-averaged Fourier spectra of CoM and SD (blue curves), and the sample-averaged peaks (red squares with error bars). The error bars of the sample-averaged peaks show their standard deviations.

Extract merging points from experimental data

To determine the merging points from the experimental data, we apply linear fitting to the data near the merging points and evaluate Vm from the fits. For the lowest-frequency mode, we use the data for V0 = 2.5, 3.0, and 3.5 kHz. For the second–third mode merging point, we use the data for V0 = 0.7605, 1.207, and 1.5 kHz.

Details of 3D GP simulation

The 3D GP equation is given by

$${{{\rm{i}}}}\hslash \frac{\partial \Psi }{\partial t}=\left[-\frac{{\hslash }^{2}}{2m}{\nabla }^{2}+{V}_{{{{\rm{ext}}}}}+g| \Psi {| }^{2}\right]\Psi ,$$
(13)

where g = 4πℏ2a/m and Ψ is the three-dimensional condensate wave function which satisfies ∫ dxdydz∣Ψ∣2 = Nc. The external potential term is composed of the trap potential \({V}_{{{{\rm{trap}}}}}=m({\omega }_{\perp }^{2}{x}^{2}+{\omega }_{\perp }^{2}{y}^{2}+{\omega }_{z}^{2}{z}^{2})/2\) and potentials along the z-direction, such as the Gaussian potential and perturbation potentials. We introduce the nondimensional variables in terms of ℏωz and the corresponding harmonic oscillator length \({a}_{z}=\sqrt{\hslash /m{\omega }_{z}}\): \(\bar{t}={\omega }_{z}t\), \(\bar{x}=x/{a}_{z}\), \({\bar{V}}_{{{{\rm{ext}}}}}={V}_{{{{\rm{ext}}}}}/(\hslash {\omega }_{z})\), \(\bar{{a}}=a/{a}_{z}\), \({\bar{\omega }}_{\perp }={\omega }_{\perp }/{\omega }_{z}\), and \(\bar{\Psi }=\Psi /\sqrt{{N}_{{{{\rm{c}}}}}/{a}_{z}^{3}}\). \(\bar{y}\) and \(\bar{z}\) are defined similarly to \(\bar{x}\). Using these variables, the 3D GP equation becomes

$${{{\rm{i}}}}\frac{\partial }{\partial \bar{t}}\bar{\Psi }=\left[-\frac{1}{2}{\bar{\nabla }}^{2}+{\bar{V}}_{{{{\rm{ext}}}}}+ {\bar{g}} | \bar{\Psi }| ^{2}\right]\bar{\Psi },\,\int \,d{\bar{x}}d{\bar{y}}d{\bar{z}}| \bar{\Psi }| ^{2}=1,$$
(14)

where \({\bar{g}} =4\pi {N}_{{{{\rm{c}}}}}\bar{{a}}\) and the trap potential is expressed as \({\bar{V}}_{{{{\rm{trap}}}}}=[{\bar{\omega }}_{\perp }^{2}({\bar{x}}^{2}+{\bar{y}}^{2})+{\bar{z}}^{2}]/2\). In the following, we omit the bars of nondimensional variables for simplicity.

We perform the 3D GP simulation using an existing 3D GP simulation package21,22. This package implements the Crank-Nicolson scheme for the discretized version of the 3D GP equation. We discretize x, y, and z, and denote the resulting discrete values by xi, yj, and zk, respectively. We use the spatial grid dx = dy = dz = 0.1 and the time step dt = 2π × 0.00042, which corresponds to 0.01 ms. The ground state Ψ0(xi, yj, zk) is obtained by imaginary-time evolution in the Oqtant trap potential setup. To break the axial symmetry, we add small fluctuations to the condensate at each grid point as

$$\begin{array}{r}{\Psi }_{0}({x}_{i},{y}_{j},{z}_{k})+\alpha [{r}_{{{{\rm{re}}}}}({x}_{i},{y}_{j},{z}_{k})+{{{\rm{i}}}}{r}_{{{{\rm{im}}}}}({x}_{i},{y}_{j},{z}_{k})]w({x}_{i},{y}_{j},{z}_{k}),\end{array}$$
(15)

where \({r}_{{{{\rm{re}}}}}({x}_{i},{y}_{j},{z}_{k})\) and \({r}_{{{{\rm{im}}}}}({x}_{i},{y}_{j},{z}_{k})\) are random numbers uniformly distributed over [ − 1, 1], α is the parameter which controls the magnitude of fluctuations, and \(w(x,y,z)=[1-\tanh ((s-1)/\delta )]/2\) with \(s=\sqrt{{(x/{R}_{{{{\rm{F}}}},x})}^{2}+{(y/{R}_{{{{\rm{F}}}},y})}^{2}+{(z/{R}_{{{{\rm{F}}}},z})}^{2}}\) is the window function to suppress fluctuations in the low-density regions. RF,x = RF,y( ≈ 1.43 μm) and RF,z( ≈ 13.6 μm) are the Thomas-Fermi lengths along the x, y and z directions. δ is the width of the window function and is chosen as 5dx. After introducing fluctuations into the condensate, it is time-evolved under the time-dependent potential shown in Fig. 3 up to 80 ms. Notice that we do not apply Vcanc since there is no unwanted initial dipole motion.