Laser spectroscopic probing of coexisting superfluid and insulating states of an atomic Bose–Hubbard system

A system of ultracold atoms in an optical lattice has been regarded as an ideal quantum simulator for a Hubbard model with extremely high controllability of the system parameters. While making use of the controllability, a comprehensive measurement across the weakly to strongly interacting regimes in the Hubbard model to discuss the quantum many-body state is still limited. Here we observe a great change in the excitation energy spectra across the two regimes in an atomic Bose–Hubbard system by using a spectroscopic technique, which can resolve the site occupancy in the lattice. By quantitatively comparing the observed spectra and numerical simulations based on sum rule relations and a binary fluid treatment under a finite temperature Gutzwiller approximation, we show that the spectra reflect the coexistence of a delocalized superfluid state and a localized insulating state across the two regimes.

U ltracold atoms confined in an optical lattice potential offer a novel way to study quantum many-body physics 1,2 . One of the most interesting problems is the quantum phase transition of ultracold bosonic atoms in a three-dimensional (3D) optical lattice from a superfluid (SF) state to a Mott insulating state. Various experimental techniques such as matterwave interference [3][4][5] , noise-correlation measurements 6,7 and quantum gas microscopes [8][9][10] are used to probe important statistical quantities such as phase coherence, the density-density correlation and the atom number distribution, which characterize the quantum phase transition. On the other hand, various spectroscopic measurement techniques such as Bragg spectroscopy [11][12][13] , lattice modulation spectroscopy [14][15][16] and radio-frequency (RF) spectroscopy 17,18 offer unprecedented potential for studying the dynamical response of interacting atoms as a result of an artificially induced perturbation. Excitation energy spectra measured with high precision allow us to obtain a deeper insight into many-body effects beyond the static and statistical observables. Furthermore, and impressively, high-resolution laser spectroscopy has recently revealed the novel SU(N) symmetry [19][20][21] .
In this paper, we report a laser spectroscopic study of the ultracold bosonic atoms in a 3D optical lattice across the weakly to strongly interacting regimes. By probing the very different onsite interactions between the different electronic states with the optical transition ( 1 S 0 -3 P 2 ) of an ytterbium (Yb) atom, our method is excellent for resolving different site occupancies, in both strongly [18][19][20] and weakly interacting regimes. To understand the excitation spectrum quantitatively, we develop a numerical method for calculating the spectrum, which is based on a finite temperature Gutzwiller approximation and formally derived sum rule relations 22 . The great and continuous change in the observed spectra as regards their shapes across the two regimes indicates that the spectroscopy captures the complicated change in the many-body state with high sensitivity, and a comparison of the spectra and the numerical simulations provides a comprehensive understanding of the observation.

Results
Theoretical framework for the laser spectroscopy. Before excitation, the electronic ground 1 S 0 state atoms in an optical lattice are in thermal equilibrium, which is described by the many-body ground-state Hubbard Hamiltonian: whereĉ y g;iĉ g;i À Á is the creation (annihilation) operator of an Yb atom in the electronic ground 1 S 0 state at the i-th site andn g;i is its number operator, and the summation hi,ji is taken over the nearest-neighbour sites. J g , V g,i and U gg represent the hopping energy, the trapping potential at the i-th site determined from the Gaussian beam shape of the trapping lasers and the on-site interaction energy, respectively. The single orbital model well describes the low energy properties of atoms before excitation except for very shallow lattices, whereas the higher orbitals are inevitably involved after excitation. The excitation operatorÔ ex is written aŝ O ex ¼ P a Z aÔex;a with second quantization. HereÔ ex;a describes the orbital-changing excitations from the lowest to the ath bands, and is given bŷ where n and k ex are the excitation laser frequency and wave vector, andĉ y e;ia is the creation operator of the electronic excited 3 P 2 state atoms in the a orbital at the i-th site of a position r i . Z a is defined by R drC Ã e;ia r ð Þe ik ex Á r À r i ð Þ C g;i r ð Þ, where C g;i r ð Þ and C e;ia r ð Þ are the Wannier orbitals of the states corresponding tô c g;i andĉ e;ia , respectively. An excitation probability from the lowest to the a th orbitals is given by |Z a | 2 , where P a |Z a | 2 ¼ 1. After the excitation, the Hamiltonian is written aŝ H¼Ĥ g þĤ e þĤ ge , wherê where D a is the energy difference between the lowest and a orbitals, and the other parameters are labelled in the same manner as equation (1). Here U ge,a represents the on-site interaction of atoms in the electronic ground state and lowest orbital with atoms in the excited state and orbital a (Fig. 1a). V e,i denotes the trapping potential for the electronic excited state atoms, which is generally different from that for the ground-state atoms. Furthermore, the interaction between two electronic excited atoms U ee can be neglected under a weak excitation condition. The excitation energy spectrum I(n) is obtained by measuring the number of electronic excited atoms as a function of n. The number conservation law of electronic excited atoms allows us to decompose I(n) into the sum of I a (n): where I a (n) is formally given in a weak excitation limit by Fermi's golden rule as where |mi represents the many-body eigenstates of the HamiltonianĤ g in equation (1) with energy E m , and Oð¼ À k B T lat ln P m e À E m = k B T lat ð Þ Þ is the grand canonical potential. Here, k B is the Boltzmann constant, T lat is the temperature after lattice loading and h is Planck's constant. The many-body excited state |n 0 i with energy E 0 n is associated with the HamiltonianĤ. For simplicity, assuming the probe time t probe to be infinite, the delta function is used to describe spectrum. Note that spectral broadening caused by the finite probe time and also certain physical origins mentioned below will be taken into account in the numerical simulations.
As it is difficult to calculate the spectrum I a (n) by directly evaluating equation (5), we use the formally derived sum rule relations of spectral moments M ';a ¼ R dn hn ð Þ ' I a n ð Þ (ref. 23). Here, M ';a can be given by certain thermodynamic quantities, that is, M 0;a ¼ P i hn g;i i, and the explicit forms of M 1,a and M 2,a are found in the Methods. As pointed out in refs 22,24, the spectral peak position p a can be determined from the first and zeroth order moments: p a ¼ M 1,a /M 0,a . For example, at T lat ¼ 0, the numberdefinite states (NDS) with m-atom filling show p NDS,iam ¼ (m À 1)dU a þ dV i,a , where dU a ¼ U ge,a À U gg and dV i,a ¼ V e,i þ D a À V g,i . This feature clearly explains that the present laser spectroscopy can be used to resolve site occupancy 24 .
The second moment allows us to determine the spectral variance s 2 a ¼M 2;a =M 0;a À ðM 1;a =M 0;a Þ 2 . For simplicity, we consider the two limited cases at T lat ¼ 0 in the uniform system, which allows us to obtain a simple form of s a from equations (10) and (11) in the Methods. These simple calculations clarify the physical origins of spectral broadening caused by the many-body effects. First, for the local NDS, the long-range correlation can be negligible: hĉ y g;iĉ g;j i ( hn i i i= ¼j ð Þ. Therefore, we obtain: where d represents a vector to the adjacent site. The spectrum for NDS is broadened by the tunnelling effects including a momentum transfer of k ex . In contrast, for the coherent phasedefinite states (PDS), where we can use a mean field approximation hĉ g;i i¼c g;i , we get the '-body correlation function G ';i ¼hðĉ y g;i Þ ' ðĉ g;i Þ ' i¼ n ' i and hĉ y g;iĉ g;j i¼c Ã g;i c g;j . Then, we obtain: For PDS, the spectrum should be broadened by number fluctuations resulting from the delocalized nature of the coherent state. Equations (6) and (7) show that the definite quantity does not cause spectral broadening, and its conjugate should be an origin of the broadening because of large fluctuations resulting from uncertainty 22 .
In general, at finite temperatures in inhomogeneous systems, phase-definite and number-definite states can coexist, and so we use binary fluid approximation in ref. 22, in which I a (n) is decomposed into two contributions from incoherent normal fluid (NF) and coherent SF atoms I a (n) ¼ I NF,a (n) þ I SF,a (n). In addition, by considering the spectral broadening mentioned above, we use the following expressions for the spectral functions: For the NF, where p NF,iam (Bp NDS, iam ) and s NF,iam (Bs NDS, a ) are the spectral position and width of the m-atom-occupying state at ith site. The spectral weight is given by i=kBT lat . On the other hand, for the SF, where spectral properties W SF,a , p SF,a and s SF,a are determined from the sum rule relations mentioned in the Methods. Note that because of the finite temperature effects and the correlations between SF and NF, the spectral properties (for example, p NF,iam and s SF,a ) are not always equivalent to the simplified forms (for example, p NDS,iam and s PDS,a , respectively). By further considering the spectral broadening caused by the finite probe time t probe ; s NF SF ð Þ will be replaced as with h=t probe of 1 kHz, which effectively describes the competition between the probe time and the tunnelling time. Supplementary Note 9 and ref. 22 provide more details.
In the numerical calculation, we adopt the finite temperature Gutzwiller approximation to obtain the thermodynamic quantities in the spectral moments M ';a . This approximation is a mean-field approximation considering up to the first-order correction in terms of J g , and thenĤ g is approximated by a set of the effective local Hamiltonians (Supplementary Note 7). Therefore, on the basis of the local approximation, the onsite '-body correlation G ';i can be obtained by diagonalizing the effective local Hamiltonian, and the long-range correlation hĉ y g;i c g;j i is described by hĉ y g;i ihĉ g;j i.
Site-occupancy resolving laser spectroscopy. In the experiments, a small portion of the 1 S 0 state atoms was directly excited to the 3 P 2 state by a narrow-linewidth laser. Since the high detection efficiency of our fluorescence measurement method enables us to detect as few as 10 atoms, we can perform the spectroscopy in a sufficiently weak excitation regime. The spectroscopy method and the pulse sequence are described in detail in the Methods, Supplementary Fig. 1 and Supplementary Note 1. The novel feature of our spectroscopy technique is its extremely high siteoccupancy resolving power (Fig. 1). A Yb atom in the 3 P 2 state has a significantly different scattering property from a 1 S 0 state atom 25,26 . Therefore, we can tune the difference of the on-site interactions |dU 0 | ¼ |U ge,0 À U gg | sufficiently large compared with the other energy scales in the Hubbard Hamiltonian 27 , such as U gg , the hopping term J g and the site-dependent inhomogeneous trapping potential V g(e),i (Fig. 1b). As a result, the optical resonances associated with different site occupancies are well resolved as shown in Fig. 1c. Note that the peaks on the negative frequency side correspond to the doubly and triply occupied sites, and this is confirmed by controlling the total atom number or the site occupancies with the photo-association technique 28 . The coexistence of such peaks demonstrates the inhomogeneous nature of our system with a harmonic trap. In addition, the high site-occupancy resolving power is highlighted by the ability to separately induce Rabi oscillations for different fillings, which also demonstrates the novel controllability of our method, as shown in Fig. 1d. We should note that as shown in Fig. 1b, the difference |dU 0 | exceeds other energy scales included in the Hubbard Hamiltonian over a wide range from the weakly to the strongly interacting regimes. This is in good contrast to the RF spectroscopy of alkaline atoms 18,24,29 , in which high site-occupancy resolution is available only in a deep Mott insulating regime because of the relatively small difference between the on-site interaction energies of the hyperfine states.
Spectra in weakly to strongly interacting regimes. In Fig. 2, we show the spectra I(n) observed by widely varying the lattice depth V L from 1 to 15E r , where E r denotes the recoil energy associated with the lattice photon. The spectrum without the lattice potential (V L ¼ 0) can be found in Supplementary Fig. 2 and Supplementary Note 3. In Fig. 2, we also show the results of a numerical simulation in which the actual experimental parameters such as the atom number, temperature and inhomogeneous trapping potential are taken into account and there is no fitting parameter except for the normalization factor. The agreement between the numerical calculation and the observed spectra is satisfactory in Fig. 2d-o. In addition, we show the matter-wave interference data in the standard absorption images in the insets in Fig. 2. The shapes of the spectra change greatly compared with the interference peaks, which simply become blurred as the lattice depth is increased.
At lattice depths greater than 11E r (Fig. 2k-o), we observe discrete peaks in the spectra, similar to those in Fig. 1c because of the negative interaction energy difference dU 0 , and the small peaks on the positive side are resonances including orbital excitation. The discrete nature suggests that the number states are realized in each site, and this can be understood from the fact that the Hubbard parameter U gg /J g exceeds 30 at these lattice depths, and the system enters the strongly interacting regime 30 . Moreover, as shown in Fig. 2d-k, the spectra undergo a complicated change from a single broad peak to several discrete peaks. At a glance, the broad spectra observed in this regime seem to be inconsistent with our site occupancy resolving power as shown in Fig. 1b at relatively deep lattice depths. However, the broad spectrum is a manifestation of the superfluidity in the system, and its width can be explained by the number fluctuation of the SF components in each lattice site (see equation (7)). Furthermore, the coexistence of the discrete and broad peaks directly reflects the coexistence of the delocalization and the localization in the system caused by the effects of the inter-atomic interaction, finite temperature and inhomogeneity 29,31,32 . It should be noted that the spectral structure also depends on the time scales of the probe and tunnelling (Supplementary Note 10). In this study, the probe time t probe ( ¼ 1 ms) is longer than or comparable with the tunnelling time of h/(12J g ) up to the deep lattice V L ¼ 11E r . Therefore, the spectrum reflects the detail of the many-body state, such as the SF, NF or their coexistence. The competition between the finite probe time and the tunnelling time is effectively considered in the numerical simulation via the spectral width as discussed above, which allows us to properly calculate the spectra in the strongly correlated regime with a small amount of or no SF component, even though t probe ( h= 12J g À Á . To understand the observed spectra including the coexistence regime quantitatively, we discuss the numerical simulation results in further detail. Decomposition of SF and NF contributions. We examine the spectral features in detail by comparing experimental and calculated results. We show the numerically calculated contributions from the SF and NF components separately in Fig. 3. The full results can be found in Supplementary Figs 6 and 7 and Supplementary Note 11. As shown in Fig. 3a (V L ¼ 5E r ) and also Fig. 2d-f (4-6E r ), we find characteristic triangular spectra centred at n ¼ 0, and the shape is reasonably explained by the numerically obtained spectra as follows. In such a shallow lattice, the SF atoms are dominant, but the NF atoms with low filling can also be seen in Fig. 3b. The filling n SF,i is distributed up to three with an average of n ave so that the SF contribution I SF (n) has one broad peak around hn¼ À n ave dU 0 j j$ À 2 dU 0 j j (Fig. 3a, SF 0 ) and another broad peak resulting from the orbital excitations around hn ¼ À n ave |dU 0 | þ D 1 (Fig. 3a, SF 1 ). In addition, the NF atoms with low filling contribute to the additional spectral component around hn ¼ 0 (Fig. 3a, red curve). Thus, the triangular spectra can be understood as a result of the SF-NF coexistence.
The spectral widths for the NF and SF components are determined by equations (6) and (7), respectively, and we also consider the inhomogeneous broadening because of the trapping potential and the excitation laser linewidth. The inhomogeneous broadenings :o inh are almost constant and B1.7 kHz in this work (Fig. 1b) and the laser linewidth is B1 kHz (Supplementary Fig. 4 and Supplementary Note 5). On the other hand, s NF,a and s SF,a depend on the Hubbard parameters, hence, the lattice depth. s 2 NDS;a is 13.5|J g | 2 in our experimental setting and it contributes B1.0 kHz to the spectral width for the NF component in Fig. 3a. The spectral width for the SF component depends on n ave and is B8.0 kHz for the a ¼ 0 resonance (SF 0 ) in Fig. 3a.
For comparison, we also performed a simple calculation at zero temperature with the following expression 24 : I n ð Þ / P i;a Z a j j 2 n i d hn À dU a G 2;i =n i À dV i;a À Á , where we substitute d(n) for a Lorentzian with a 1-kHz linewidth. Figure 3a shows an apparent discrepancy between the calculated spectrum (green dashed line) and the experimental results. In particular, a peak at around n ¼ 0 is missing, which is the NF contribution within the binary fluid calculations. Note that the large intensity of I NF n $ 0 ð Þ results from the dilute and widely spread distributions of the NF atoms (see Fig. 3b for details).
As the lattice depth is increased further, the triangular shape is distorted and discrete peaks appear that reflect a noticeable    (8) and (9), respectively. The grey circles denote the same data shown in Fig. 2 for comparison.
increase in the incoherent component I NF as shown in both Fig. 4a (V L ¼ 7E r ) and Fig. 2g-k (7-11E r ). In addition, the discrepancy between the experimental and numerical results is more pronounced: the discrete peaks near n ¼ 0 are more contrasted in the experimental result ( Fig. 2h-k). In the numerical calculation presented in Figs 2 and 3, we assume the adiabatic atom loading into the lattice potential with the initial entropy of S ini calculated from the initial temperature T ini by neglecting the interaction effects for the Bose gas. Here, T ini is determined by a standard time-of-flight measurement and the laser spectroscopy without lattice potential, both of which allow us to calibrate the thermal fraction of Bose gases. To cover non-adiabatic heating effects and also an approximation for the estimation of S ini , we also performed the numerical simulation with higher initial temperatures T ini ¼ 110, 125 nK (Supplementary Note 8 and Supplementary Fig. 5). More comprehensive study including lower T ini can be found in Supplementary Fig. 8.
The numerical results at higher T ini show a distinct structure resulting from the incoherent NF contributions. As shown in Fig. 4a-c,d-f, thermal fluctuations destroy the SF and increase the NF, and the discrete structure are more visible. The experimental results with V L ¼ 7E r and 10E r are well explained with the numerical results with T ini ¼ 95 and 110 nK, respectively. In the deep lattice as shown in Fig. 4g-i (V L ¼ 13E r ), the higher temperatures modify the relative height of the discrete peaks, and the experimental result is well explained by the results calculated at T ini ¼ 125 nK. The heating behaviour is consistent with the visibility measurement in the matter-wave interference, which also indicates heating due to the atom loading sequence (Supplementary Fig. 3 and Supplementary Note 4). We finally discuss the noticeable deviation of the peak positions of multi-occupancies at deep lattices (for example, see Fig. 4g-i at around n ¼ À 30 kHz). As regards the m-atomoccupancy in our simulation, we simply assume that a peak position is determined from (m À 1)dU 0 . A more elaborate calculation that considers the interaction-broadening of the Wannier function, where the interaction strengths depend on m as dU 0 (m), is necessary to accurately predict the peak position.

Discussion
We have employed a laser spectroscopy technique for ultracold bosonic atoms in a 3D optical lattice and performed a numerical analysis based on the formally derived sum rule relations. By using the very different interaction energies in the different electronic states, our spectroscopy technique provides outstanding resolution from the weakly to the strongly interacting regimes, and the observed spectra exhibit great changes in their shapes. Our numerical method enables us to calculate the excitation spectrum in finite temperature as a steady-state solution in the weak excitation regime, and effectively deals with the effect of the finite probe time in the experiment as a laser linewidth. The combination of the experimental and numerical study gives us a quantitative understanding of the coexistence of the SF and insulating states. Our approach paves the way to an exploration of the finite temperature Bose-Hubbard phase diagram including the quantum critical point, and can be applied to other types of the quantum many-body system, such as fermionic or Bose-Fermi mixtures in the optical lattice.

Methods
Preparation of ultracold Yb atoms in an optical lattice. We first prepare lasercooled 174 Yb atoms in a crossed optical dipole trap 26 . The wavelength of the dipole trap laser is 532 nm. After evaporative cooling, we obtain a Bose-Einstein condensate (BEC) of 174 Yb, and load the BEC into a 3D optical lattice with a simple cubic geometry and a lattice constant of 266 nm by ramping up the optical lattice laser intensity. The lattice potential is linearly increased to a desired depth at a ramping rate of 15 À 1 E r /ms. The lattice potential is superimposed with a harmonic confinement with a trapping frequency of 140 Hz at 15E r . The intensity of the lattice laser is precisely controlled by the feedback of a monitored laser power. The  Fig. 3. The grey circles denote the same data shown in Fig. 2 for comparison. lattice depth is calibrated by observing the diffraction of the trapped atoms induced by a pulsed lattice potential. In all of the spectroscopy experiments shown in Fig. 2, we fix the number of atoms N at B2.2(3) Â 10 4 , the BEC transition temperature T c at B160 nK, and the initial atom temperature T ini at B95 nK.
Spectroscopy method. To achieve the weak excitation of the atoms, we measure the number of excited atoms by using the repumping method. The ground-state atoms are directly excited to the 3 P 2 (m ¼ 2) state by using a laser pulse with a wavelength of 507 nm and a pulse duration of 1 ms, where m denotes the magnetic quantum number. After the excitation pulse, we remove the ground-state atoms by using a blasting pulse that is resonant to the strong electric dipole 1 S 0 2 1 P 1 transition (399 nm). Then, we repump the atoms in the excited state by using two repumping transitions, namely the 3 P 0 2 3 S 1 transition (649 nm) and the 3 P 2 2 3 S 1 transition (770 nm). The pulse durations for the blasting and the repumping are 0.2 and 0.5 ms, respectively. To detect the repumped atoms with a high signal-tonoise ratio, we recapture the atoms by using a magneto-optical trap (MOT) with the electric dipole 1 S 0 2 1 P 1 transition, and measure the fluorescence from the MOT with an electron-multiplying charge-coupled-device camera with an exposure time of 100 ms. Although the atoms in the 3 P 2 (m ¼ 2) state suffer from a Zeeman-sublevel-changing collision (Supplementary Note 6) and acquire kinetic energy during the inelastic collision, the MOT can capture the repumped atoms thanks to the large capture volume and the strong trapping force.
Explicit formulae of sum rule relations. The sum rule imposes the condition that the first-order moment M 1;a ¼h½Ô y ex;a ;ĤÔ ex;a i should be expressed by a combination of various thermal quantities or correlation functions for the atoms in the ground state before excitation: It is noted that the sixth term represents correlated hopping. The first-order sum rule provides information about a spectral peak position, which is discussed in refs 24,29, and the second-order sum rule indicates the variance of the spectra. The comparison between the moments calculated from numerical and theoretical results is discussed in Supplementary Fig. 9 and Supplementary Note 12.
Hubbard parameters. The Hubbard parameters, namely the hopping matrix J g(e) , the trapping potential V g(e),i and the interaction strength U gg,(ge) and the excitation probability Z a , are determined in an ab initio manner by numerically calculating the Bloch and Wannier orbitals. We include the difference in the polarizability of the 1 S 0 and 3 P 2 states in the calculation. In addition, we consider the formation of two-body bound states consisting of the 1 S 0 and 3 P 2 atoms induced by the confinement of the lattice potential 33 , which causes the renormalization of the interaction strength U ge . By using the ab initio Bloch waves, we determine the two-body bound states 33 . This bound state formation only occurs when the energy scale of the onsite interaction becomes comparable to the Hubbard band gap 27,33,34 . We neglect other types of interaction renormalization, such as self-trapping 35 , because the number of atoms is small. The temperature in the lattice is determined on the assumption that the atom loading process is adiabatic. The thermal equilibrium temperature in the lattice, T lat , depends on V L .