Research ArticlePHYSICS

Energy redistribution and spatiotemporal evolution of correlations after a sudden quench of the Bose-Hubbard model

See allHide authors and affiliations

Science Advances  30 Sep 2020:
Vol. 6, no. 40, eaba9255
DOI: 10.1126/sciadv.aba9255


An optical lattice quantum simulator is an ideal experimental platform to investigate nonequilibrium dynamics of a quantum many-body system, which is, in general, hard to simulate with classical computers. Here, we use our quantum simulator of the Bose-Hubbard model to study dynamics far from equilibrium after a quantum quench. We successfully confirm the energy conservation law in the one- and three-dimensional systems and extract the propagation velocity of the single-particle correlation in the one- and two-dimensional systems. We corroborate the validity of our quantum simulator through quantitative comparisons between the experiments and the exact numerical calculations in one dimension. In the computationally hard cases of two or three dimensions, by using the quantum-simulation results as references, we examine the performance of a numerical method, namely, the truncated Wigner approximation, revealing its usefulness and limitation. This work constitutes an exemplary case for the usage of analog quantum simulators.


Rapid advances in analog quantum simulation using highly controllable systems with long coherence time, such as ultracold gases in optical lattices [for example, see (1, 2)], Rydberg atoms in an optical tweezer array [for example, see (35)], and trapped ions [for example, see (6, 7)], have notably expanded possibilities for studying dynamics of quantum many-body systems. One of the recent targets of optical lattice quantum simulators has been the investigation of the nonequilibrium dynamics arising after a quantum quench (817), where a parameter of the system is varied rapidly and substantially. In the case of one dimension (1D) for a short time scale, quantum quench dynamics can be exactly computed with classical computers by means of the matrix product state (MPS) method [for example, see (18, 19)]. In pioneering works of quantum-simulation research, the outputs of experiments were directly compared with those of exact numerical simulations with classical computers to examine the performance of the quantum simulators (1012).

A two-point spatial correlation as a function of the distance of the two points has been the intense theoretical interest (2030), and in fact, in 1D systems, it has been shown that access to such a correlation function allows for exploring the dynamical spreading of quantum information, which is of great interest in connection with the Lieb-Robinson (LR) bound (11, 31). An exact computation of the spatiotemporal evolution of these two-point correlations is, however, generally intractable for a long time scale or in higher dimensions. While a more recent work has used outputs from a quantum simulator built with ultracold fermions in a Floquet-engineered optical lattice in 3D as a reference for examining the performance of an approximate numerical method, namely, the nonequilibrium dynamical mean-field theory (32), a direct comparison with quantitative theoretical approaches in the quench dynamics in higher dimensions is still lacking.

Here, we investigate the energy redistribution dynamics and the spatiotemporal evolution of the single-particle correlation function, which is one of the simplest two-point spatial correlations, in quantum quench dynamics starting with a Mott insulating state by using an optical lattice quantum simulator of the Bose-Hubbard model (BHM) in 2D and 3D as well as 1D. The observation of the redistribution of the kinetic and interaction energies turns out to be the confirmation of the energy conservation in the quench dynamics of a Bose-Hubbard quantum simulator. Furthermore, we successfully observe the correlation spreading after a rapid quench from a Mott insulating state toward the quantum critical region in 2D as well as toward the Mott region in 1D. We compare the measured propagation velocity of the correlation front, which is defined from the first peak in the time evolution of the correlation function at each distance, with the LR-like bound set by the maximum velocity of the quasiparticles. In the 2D case, we find that the former velocity exceeds the latter one. This happens because the single-particle correlation spreads with two typical velocities, namely, the group velocity and the phase velocity, as was pointed out in the recent theoretical work (30), and the measured velocity corresponds to the phase one. Since the first peak propagating with the phase velocity decays rapidly with the distance, our observation is not contradicting the existence of the LR bound implying that any correlation functions outside the LR light cone must be exponentially suppressed.

In addition to these experimental findings, to corroborate the quantitative performance of our quantum simulator, we present a thorough comparison between the quantum-simulation results and state-of-the-art theoretical calculations. We use the exact MPS method in the 1D case, finding excellent agreement with the observations. As for the case of the quench toward a deep superfluid region in 3D, the time evolution of the kinetic and interaction energies is directly compared with numerical results obtained using the truncated Wigner approximation (TWA) on the basis of the Gross-Pitaevskii mean-field theory (29). The good agreement between the experiment and the theory establishes the predictive power of the TWA for this type of quench. In contrast, in the case of the quench toward the quantum critical region in 2D, the TWA fails to capture quantitatively the experimental results, although it captures some qualitative features. This indicates that our quantum simulation goes beyond current classical computation and the data serve as a useful reference for pushing out its boundary.


Investigating nonequilibrium dynamics of the BHM

We consider a system of ultracold bosonic atoms confined in an optical lattice. When an optical lattice potential is deep, the system is quantitatively described by the BHM (33, 34)Ĥ=JΣj,l(âjâl+h.c.)+U2Σjâjâjâjâj+Σj(Vjμ)âjâj(1)where âj and âj are the creation and annihilation operators at the site j, J is the tunneling-matrix element between nearest-neighbor sites, U is the on-site interaction energy, μ is the chemical potential, and Vj is the local potential offset at the site j, which originates from the trap potential and the Gaussian envelopes of optical lattice lasers. Σj, l represents the summation over all neighboring sites. The position of the site j is denoted by rj=Σα=1Dxjαeα, where x1, x2, and x3 mean x, y, and z, respectively. eα represents the unit vector in the xα-direction, and D is the spatial dimension.

When the atom number per site, namely, the filling factor n̄, is an integer and the ratio U/J is varied, the BHM exhibits a second-order quantum phase transition between the Mott insulator and the superfluid. The system favors the superfluid phase for a relatively small U/J, while it does the Mott insulator phase for a relatively large U/J. For the unit filling case (n̄=1), the quantum critical point has been determined with exact numerical methods as (U/J)c= 3.4 (1D), 16.7 (2D), and 29.3 (3D), respectively [for review, see (35)].

Our analog quantum simulator of the BHM is built with an ultracold Bose gas of 174Yb atoms confined in a 3D optical lattice. We use this 174Yb-atom-BHM quantum simulator to analyze dynamics after a quench of the ratio U/J starting with a Mott insulator state with unit filling. We convert a 174Yb Bose-Einstein condensate (BEC) in a weakly confining harmonic trap into the initial Mott insulator state by slowly ramping up the optical lattice depth up to sV0/ER = 15 for all the three directions, where V0 is the depth of the optical lattice and ER is the recoil energy of the optical lattice laser whose wavelength is 532 nm. See Methods for the preparation. The prepared state is deep in a Mott insulator regime (U/J = 100) and is well approximated as a product of local Fock statesΨMI=jâj0(2)

To realize a quench of U/J, we rapidly ramp down the lattice depth for some directions toward a final value. For instance, in the case of the 1D quench, we ramp down the lattice depth only for the x direction while in the 3D case, we do it for all the three directions. The ramp-down speed is set to be 100 ER/ms. By using the band-mapping techniques, we check that there is no discernible amount of the atoms in excited bands with this quench speed. We use the numerical values of U and J calculated as functions of lattice depth reported in (35) (see section S1).

After the quench process, we keep the lattice depth constant and let the system evolve. To obtain the single-particle correlation function at a certain distance Δ=Σα=1DeαΔxα, where Δxα ≥ 0, in the unit of a lattice spacing d (=266 nm) (36)KΔ=Σα=1DΣxjαxlα=dΔxαâjâl(3)after a certain hold time, we release the gas from the trapping and optical lattice potentials to measure the time-of-flight (TOF) image, from which we deduce the momentum distribution (see Methods for details). By performing a Fourier transform of the momentum distribution, we obtain KΔ. The kinetic energy of the BHM is equal to the sum of −JKΔ at Δ = 1, where Δ = ∣ Δ∣ (36).

Moreover, we measure the onsite-interaction energy of the BHM, U2Σjâjâjâjâj, by means of the atom number projection spectroscopy (36, 37) for the 3D case and photoassociation spectroscopy (38) for the 1D case where we confirm that there are almost no multiple occupancies larger than two (see Methods for details). Compared to the methods based on the quantum gas microscope techniques (11), our methods are rather efficient, especially in higher dimensions, for our current purposes of obtaining the ensemble average of the two-point correlation functions and the Hubbard energies, because less repetitions are needed thanks to the much larger number of atoms.

The experimental procedure and setup and typical high-resolution spectra are summarized in sections S1 and S2. It is worth noting that the dynamical evolution of the phase correlation, which is similar to the single-particle correlation, has been measured for weakly interacting Bose gases in 1D optical lattices by means of Talbot interferometry in (39). In contrast, the present work investigates the single-particle correlation in strongly correlated regimes in higher dimensions.

Experimental confirmation of our methods: Dynamics of the 1D BHM after a sudden quench

First, we investigate the behaviors of atoms after a sudden quench in 1D. The results for the dynamical redistribution of the Hubbard energies and the spatiotemporal evolution of the atom correlations are shown in Figs. 1 and 2, respectively. Specifically, we ramp the lattice depth in the x direction down to s = 5 implying U/J = 6.8, where the ground state is a Mott insulator state close to the quantum critical point. Figure 1 shows the time evolution of the kinetic energy, the onsite-interaction energy, and the sum of the two energies. On a short time scale, while the sum of the two remains almost constant, the kinetic energy decreases and the interaction energy increases. After making a small overshoot, each energy ends up with an almost steady value, i.e., the energies are redistributed. These behaviors are expected for an isolated system but have never been observed experimentally before.

Fig. 1 Energy redistribution after the quench in 1D.

The kinetic-energy term (red), the onsite interaction energy term (blue), and the sum of them (green) are shown as functions of the hold time t after a rapid quench into a Mott insulator region with U/J = 6.8 in 1D optical lattice tubes. The solid lines show the results of the numerical calculation at zero temperature with the MPS method using the time-dependent variational principle and local density approximation. The error bars for the kinetic-energy and onsite interaction energy terms denote the SE of 15 independent measurements.

Fig. 2 Spatiotemporal evolution of the single-particle correlation after the quench in 1D.

(A and B) The single-particle correlations for the distance in the unit of the lattice constant Δ up to 4 are shown as functions of the hold time t. Note that the displayed correlations are normalized by the maximum value of the correlation Cmax,Δ(1D) during 0 < t < 1.6ℏ/J for each distance Δ; (A) experiment and (B) numerical calculation. (C to F) Time evolution of the single-particle correlation KΔ after the quench. Solid blue lines show the results of numerical calculation. (C) Δ = 1, (D) Δ = 2, (E) Δ = 3, and (F) Δ = 4. The error bars denote the SE of five independent measurements. (G) Time of the first peak of the single-particle correlation is plotted as a function of the distance Δ. A fit with a linear function with a nonzero offset is shown as a solid line. The error bars denote the SE of five independent measurements.

Figure 2 (A to D) shows the time evolution of the single-particle correlation function for several values of Δ. As the time evolves, the correlations first grow and each of them has the first (local) maximum at a certain time. We extract the peak time for each Δ by numerically fitting to the experimental data, which is plotted against Δ in Fig. 2G. The peak time increases linearly with the distance, i.e., the correlation exhibits a light-cone-like propagation. From the peak time versus Δ, we extract the propagation velocity as v = 5.5(7)Jd/ℏ.

The maximum velocity of a particle-hole excitation is given as (11, 21)vmax6JdD[116J29U2](4)which can be interpreted as an LR-like bound. It is noted that vmax corresponds to the sum of the maximum velocities of the doublon and the holon, which are respectively given by 4JdD/ and 2JdD/ in the leading order with respect to J/U. As long as U/J ≫ 1, Eq. 4 is valid regardless of the spatial dimension. At U/J = 6.8, vmax = 5.8Jd/ℏ such that the condition v < vmax is satisfied, as expected. A similar propagation behavior has been also observed in the case of the density-density correlation (11). In contrast, we will see later that v > vmax in the 2D case. We will explain that this observation is still compatible with the LR-like bound.

While these observations reveal important features of the nonequilibrium dynamics of BHM, this 1D study is also important from another aspect. Since our BHM quantum simulator is analog, it is imperative to examine its accuracy through a direct comparison with exact numerical calculations in 1D before applying it to the cases of higher dimensions, in which exact computation on classical computers is currently unavailable. In Figs. 1 and 2, we compare the experimental results in 1D with the exact numerical ones at zero temperature obtained with the MPS method. For details of the MPS calculations, see section S4. We see that the experimental observations are in good agreement with the exact numerical calculations with no fitting parameters.

Dynamics of BHM after a sudden quench in higher dimensions

Having corroborated the quantitative validity of our BHM quantum simulator by the comparison between the theory and experiment in 1D, we now discuss the main result of this work, i.e., the quench dynamics in higher dimensions. Figure 3 shows the energy redistribution dynamics for the 3D case after the ramp-down of the lattice depth to s = 5 (U/J = 3.4), where the ground state is deep in the superfluid phase. The general tendency of the time evolution is similar to the 1D case: The two energies are redistributed on a time scale smaller than ℏ/J and the sum of the two remains almost constant within the displayed time window t ≲ ℏ/J.

Fig. 3 Energy redistribution after the quench in 3D.

The kinetic-energy term (red), the onsite interaction energy term (blue), and the sum of them (green) are shown as functions of the hold time t after a rapid quench into a superfluid region with U/J = 3.4 in a 3D optical lattice. The solid lines show the results of the numerical calculation obtained using the TWA (29). The error bars for the kinetic-energy (onsite-interaction energy) terms denote the SE of 15 (3) independent measurements.

We next investigate the dynamical spreading of the single-particle correlation after a quantum quench in 2D. The final lattice depth in this case is s = 9 implying U/J = 19.6, where the ground state is a Mott insulator phase near the quantum critical point. Figure 4 (A to D) shows the spatial distribution of the single-particle correlation at several hold times after a quench. We clearly observe that the correlation first grows between nearest-neighbor sites and then it propagates for larger distances at later times. More directly, Fig. 4 (E to G) shows the time evolution of the single-particle correlation function KΔx, Δy for several values of (Δx, Δy), where Δxy) denotes the distance in the x (y) direction in units of the lattice spacing d. The delay in the growth of the correlation for longer distance is clearly observed along the directions of x, y, and x + y, in Fig. 4 (E, F, and G, respectively). In the same manner as the 1D case, we extract the position of the first peak in the time evolution of the correlation at each distance, which is plotted against the Euclidean distance Δx2+Δy2 in Fig. 4H. We further extract the propagation velocities from the linear fitting to Fig. 4 (H and I) as v = 13.7(2.1)Jd/ℏ (peak) and v = 10.2(1.4)Jd/ℏ (trough). According to Eq. 3, the maximum velocity of the particle-hole excitation is vmax = 8.4Jd/ℏ, which is slower than the observed propagation velocity in Fig. 4.

Fig. 4 Spatiotemporal evolution of the single-particle correlation after the quench in 2D.

(A to D) 2D plots of the single-particle correlation as functions of the distances Δx and Δy for several hold times of tJ/ℏ = 0 (A), 0.12 (B), 0.23 (C), and 0.35 (D). Data with Δx2+Δy24 are shown. Note that the displayed correlations are normalized by the maximum value of the correlation Cmax,Δ(2D) during 0 < t < 1.0ℏ/J for each distance (Δx, Δy). (E to G) Time evolution of the single-particle correlation KΔ after the quench for (Δx, Δy) = (1,0) [(E), red square], (2,0) [(E), green circle], (3,0) [(E), yellow diamond], (0,1) [(F), red square], (0,2) [(F), green circle], (0,3) [(F), yellow diamond], (1,1) [(G), green circle], and (2,2) [(G), yellow diamond]. The solid lines are the numerical results obtained using the TWA method. The error bars denote the SE of 15 independent measurements. (H and I) Time at the first peak (H) or the first trough (I) of the single-particle correlation as a function of the Euclidean distance Δ=Δx2+Δy2. A fit with a linear function with a nonzero offset is shown as a solid line both in (H) and (I). The first peak and trough are obtained by fitting the experimental data with the empirical function described in section S6. The error bars denote the fitting errors.


The observation that v > vmax in the 2D case shown in Fig. 4 can be interpreted along the line explained in (30). The correlation spatially propagates as a wave packet, whose width spreads in time. This means that the velocity of the first peak in the time evolution of the correlation function at each distance, namely, the phase velocity, is faster than that of the center of the wave packet, namely, the group velocity. Moreover, the first peak decays rather rapidly as the distance becomes larger. The velocity extracted from the experimental data in the way described above corresponds to the phase velocity, while the meaningful propagation velocity, which should be compared with the LR-like bound, does to the group velocity. In the case of the final lattice depth s = 9 in 2D, we cannot accurately extract the group velocity because of the unclear separation of the two velocities. Instead, In section S5, we show an example, in which the phase velocity is well separated from the group velocity in the 1D case with large U/J. There, we also see that v > vmax. Hence, the behavior that v > vmax is not unique to the 2D case but can emerge regardless of the spatial dimension as long as U/J after the quench is sufficiently large. Notice that we observed v < vmax in the case of the final lattice depth s = 5 (U/J = 6.8) in 1D because the phase velocity is approximately equal to the group velocity at U/J = 6.8 (30).

Next, we discuss the usefulness and limitation of some numerical methods based on the quantum simulation results. Since there is no exact computation method applicable to the 2D and 3D cases, it is meaningful to examine the accuracy of some approximate methods by using the quantum simulation results as a quantitative reference. In (32), the time-dependent dynamical mean-field theory has been examined by comparison with quantum simulation results for real-time dynamics of the Fermi-Hubbard model. This method is not suited for computing the nonlocal spatial correlations analyzed in the present work because it ignores the momentum dependence of the correlation functions. Instead, we choose the TWA approximation, which is supposed to accurately capture semiclassical dynamics of the BHM at least on a short time scale [see (29) and references therein]. In Fig. 3, where the energy redistribution dynamics in 3D is depicted, we also show the numerical calculations as solid lines obtained with the TWA (29). In the TWA calculations, we take the Mott insulator state of Eq. 2 as the initial state and set the system size to be 303 sites. We ignore the trapping potential because it is irrelevant to the dynamics within the time window t ≲ ℏ/J as was discussed in the 1D case. The TWA results are in good agreement with the experimental observations. More details of the TWA calculations are described in (29).

Let us turn our attention to the correlation spreading in 2D shown in Fig. 4. The solid lines in Fig. 4 (E to G) represent the results obtained by using the TWA. The TWA agrees with the experiment on a very short time scale (t < 0.1ℏ/J). Moreover, the peak positions and the values at a relatively long time (t > 1ℏ/J) for a short distance, say Δ = 1, are reasonably captured. However, it fails to capture some important properties of the correlation dynamics, such as the locations of the correlation troughs and the almost converged value of the correlation for Δ > 1. This disagreement is consistent with the general fact that the TWA is less accurate when U/(Dn̄J) or tJ/ℏ is larger. This failure of the TWA indicates that one needs to push out the boundary of currently available numerical techniques for quantitative description of the quantum simulation results. One possible candidate is to extend the SU(N) TWA (40, 41) for analyzing the BHM with unit filling.

In both of the 1D and 2D quench cases, we observed that the peaks propagated linearly with a constant velocity (see Figs. 2G and 4, H and I). However, extrapolations to t = 0 have nonzero offsets. In addition, our numerical results also support the existence of the offsets. The offsets reflect the difference between the speed for the creation of a particle-hole pair and that for its propagation. The former speed determines the time giving the first peak at Δ = 1, while the latter does those at Δ > 1. It is noted that the dependence of the propagation velocity on distance in the case of the 1D quench was already numerically discussed in (11).

Our quantum simulation platform for studying nonequilibrium dynamics can be straightforwardly applied to other quantum many-body systems such as the Fermi-Hubbard model [with SU(N) symmetry (42, 43)], the Bose-Fermi Hubbard model, and the spinful BHM. In addition, it is interesting to extend our work to a study of quench dynamics on a quantum system with controlled dissipation, which has recently attracted much interest (44).


Preparation of the initial Fock state

Details of our experimental setup are described in (36). We first prepare a BEC of 174Yb atoms confined in an optical far-off resonant trap (FORT) whose wavelength is 532 nm. The trap frequencies of the FORT are given by (ωx, ωy, ωz) = 2π × (28,130,160) Hz, where the x and y axes were tilted from the x and y axes, to which two of the optical lattices are directed, by 45. Then, we slowly ramp up the optical lattice depth for all the three directions from s = 0 to 5 in 100 ms and from 5 to 15 in another 100 ms. A typical number of atoms is chosen to be N = 1.3 × 104 such that the filling factor is unity.

Lattice quench

We perform the quench by sudden decrease of the optical lattice with depth of s ER in 0.01(15 − s) ms. See also section S1. The excitation of the atoms into higher bands is negligible with this procedure. For the cases of the 1D and 2D quench, we decrease the lattice depth along the one direction of x and two directions of x and y, respectively. It is to be noted that when lattice depth is 10.6 ER, U/J is equal to 29.34, which is the critical lattice depth for the superfluid-Mott transition at n̄=1.

Measurement of the ensemble average of the nonlocal atom correlation

Here, we briefly describe a method for obtaining the ensemble average of the nonlocal atom correlation KΔ of Eq. 3. Details are described in (36). The atomic-density distribution n(r) after the TOF t is given byn(r)=(mt)3w˜0(k)2S(k)(5)where w0(k) is the Fourier transformation of the Wannier function in the lowest Bloch band w0(r) and k = mr/ℏt.

When t is long enough, the structure factor S(k) is expressed asS(k)=Σj,leik·(rjrl)âjâl(6)where 〈 · 〉 represents the ensemble average. Therefore, the ensemble average of the nonlocal atom correlation KΔ can be easily obtained by Fourier transformation. For real experiments, two factors should be taken into account: an interaction effect and the finite-TOF effect.

A careful estimation of our experimental conditions (36) shows that the ratio of the interaction energy Un(n − 1)/2 to the kinetic energy ℏωL is mostly far lower than 1, justifying our ignorance of the interaction effect during TOF. The finite-TOF effect is small but not negligible so that we determined the nonlocal atom correlation by extrapolation based on the theoretical model described in (36).

Measurement of ensemble average of the interaction energy

To measure the ensemble average of the interaction energy (1/2)UΣiâiâiâiâi = (1/2)UΣin̂i(n̂i1), a method for projecting the distribution of the atom number per site on an observable, namely, the atom number projection method, is required. Details are described in (36). First, we increase the optical lattice depth quickly to freeze the hopping of atoms. The ramp-up time is smaller than the hopping time but large enough to prevent the atoms from being excited into the higher band of the optical lattice. For example, the ramp-up time is 0.1 ms from 5ER to 15ER.

Subsequently, we perform a site-occupancy-resolved spectroscopy. We use two methods: the high-resolution spectroscopy using the optical transition between the 1S0 and 3P2 (mJ = 0) electronic states of Yb atoms and the photoassociation spectroscopy. The excellent resolution of the spectroscopy using the 1S0 and 3P2 (mJ = 0) transition allows us to distinguish different site occupancies, owing to quite different two-body interactions of Ueg/h = −8.5 kHz and Ugg/h = 3.2 kHz at 15ER. From the area of the spectra, we obtain the total number Nn of n-occupied sites. The interaction energy is obtained as U = (1/2)ΣnNnn(n − 1). The correlation factors induced by occupancy-dependent Rabi frequencies and a loss of atoms in the 3P2 (mJ = 0) state during the spectroscopy were studied in our previous work (36).

Another method that we use is the photoassociation, which is a process to create one molecule from two atoms by light. The created molecule rapidly escapes from the trap so that we can measure the total number of doubly occupied sites as the loss of the atoms. It is noted that the method is invalid in the case of triple and higher occupancies. For example, photoassociation in a triply occupied site induces only two-atom loss and one atom remains, which is the same result as the case of a doubly occupied site, concerning the loss of atoms.

In our experiment, we use the high-resolution spectroscopy using the optical transition between the 1S0 and 3P2 (mJ = 0) for the 3D quench experiment. In contrast, we use the photoassociation for the 1D quench experiment, where we additionally check the absence of triply occupied sites or higher by means of the high-resolution spectroscopy.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J. Sakamoto for experimental assistance and T. Kuno, Y. Watanabe, and T. Sagawa for comments and discussions. The MPS-based calculations in Figs. 1 and 2 and fig. S6 were performed using the ITensor Library ( The TWA simulation in Fig. 4 was carried out at the Yukawa Institute Computer Facility. Funding: This work was supported by the Grant-in-Aid for Scientific Research of the Ministry of Education, Culture Sports, Science, and Technology / Japan Society for the Promotion of Science (MEXT/JSPS KAKENHI) nos. JP25220711, JP16H00801, JP17H06138, JP18H05405, JP18K03492, and JP18H05228; the Impulsing Paradigm Change through Disruptive Technologies (ImPACT) program; Japan Science and Technology Agency CREST (no. JPMJCR1673); and MEXT Quantum Leap Flagship Program (Q-LEAP) grant no. JPMXS0118069021. Author contributions: Y.Takasu and Y.Takahashi conceived and planned the experiments. Y.Takasu, T.Y., H.A., and Y.F. carried out the experiments. Y.Takasu analyzed the experimental data. K.N., S.G., and I.D. developed the theoretical framework and planned and carried out the numerical simulations. Y.Takasu, I.D., and Y.Takahashi wrote the paper with input from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the corresponding author.
View Abstract

Stay Connected to Science Advances

Navigate This Article