{'fraction_non_alphanumeric': 0.04582994565449955, 'fraction_numerical': 0.010847704223544011, 'mean_word_length': 4.20189204229271, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Theoretical Analysis on Spectroscopy of Atomic Bose-Hubbard Systems

We provide a numerical method to calculate comprehensively the microwave and the laser spectra of ultracold bosonic atoms in optical lattices at finite temperatures. Our formulation is built up with the sum rules, up to the second order, derived from the general principle of spectroscopy. The sum rule approach allows us to discuss the physical origins of a spectral peak shift and also a peak broadening. We find that a spectral broadening of superfluid atoms can be determined from number fluctuations of atoms, while that of normal-state atoms is mainly attributed to quantum fluctuations resulting from hopping of atoms. To calculate spectra at finite temperatures, based on the sum rule approach, we provide a two-mode approximation assuming that spectra of the superfluid and normal state atoms can be calculated separately. Our method can properly deal with multi-peak structures of spectra resulting from thermal fluctuations and also coexisting of the superfluid and the normal states. By combining the two-mode approximation with a finite temperature Gutzwiller approximation, we calculate spectra at finite temperatures by considering realistic systems, and the calculated spectra show nice agreements with those in experiments.


I. INTRODUCTION
Ultracold atoms in an optical lattice allows us to simulate quantum phase transitions of lattice fermions and also bosons [1,2]. In fact, the superfluid (SF) to the Mott insulator (MI) transition of bosonic atoms has been demonstrated by using various measurement techniques [3][4][5][6][7][8][9][10][11][12][13][14][15]. The signature of the phase transitions can be observed in certain thermodynamic quantities [3][4][5][6][7][8][9]. One of examples is to characterize the transition by observing the disappearance of a coherent peak structure in the number distribution of atoms in the momentum space [3][4][5]. A spectroscopic measurement is another useful tool to detect phase transitions [10][11][12][13][14][15]. This is because much information is included in spectra that reflect the dynamical response of many-body systems after excitations caused by a certain external field. Furthermore, when the external field is very weak and perturbative, the dynamical response can be connected to thermodynamic quantities of thermal equilibrium states before excitations. In condensed matter physics, such a relationship, e.g., fluctuation-dissipation theorem, has been used to discuss quantum many-body phenomena. It is thus required to deeply discuss such spectroscopic relationships specific to cold-atom systems.
One of pioneering studies on spectroscopic measurement on atoms in a lattice is microwave spectroscopy experiments, where the Mott shell structure has been observed by spectroscopically distinguishing the different number states of atoms [11]. Theoretically, the corresponding spectra have been studied with an approximation satisfying the (first order) spectral sum rule [16][17][18], which is derived from the general principle of spectroscopy [19]. The first-order sum rule determines the relationship between the spectral peak position and the two-body correlation function of atoms [16]. This is a prominent example that connects thermodynamics to dynamics in cold atom systems. This calculation assumed that the system is at zero temperature, while the realistic experiments have been done at low but finite temperatures. In addition, such a first order approximation is insufficient to discuss important properties of spectra, such as, a standard deviation and a spectral broadening, which can be connected to fluctuations of atoms in thermal equilibrium. On the other hand, the laser spectroscopy is now being established [20][21][22]. The laser and the microwave spectroscopy are understood as a similar type spectroscopy based on electromagnet-field excitations. However, the laser spectroscopy cannot be straightforwardly described by the formulation of the previous studies [16]. This is mainly due to the difference in wavelengths of the external fields. A reliable theoretical method for comprehensively analyzing these spectroscopy at finite temperatures is now required.
In this paper, we theoretically discuss a common formulation for the microwave and the laser spectroscopy of ultracold bosonic atoms in a three-dimensional optical lattice. We start with analyzing the sum rules in the same way as the previous work [16], while we extend the approximation to the second order. This approach allows us to clarify that number fluctuations of atoms in thermal equilibrium can be connected to a broadening of spectra. Phenomenological discussions based on the sum rule approach allow us to establish a method for calculating spectra at finite temperatures. We propose a two-mode approximation assuming that the spectra of condensed SF atoms and uncondensed normal state (NS) atoms are separately dealt with. The multi-peak structures resulting from thermal fluctuations and also the coexisting of the SF and NS atoms can be appropriately taken account. Using this approximation combined with a finite temperature Gutzwiller approximation [23], we numerically calculate the microwave and the laser spectra by considering realistic experimental parameters [11,24]. We find that our approximations reproduce essential features of spectra seen in the microwave experiments [11], and we predict spectra of the realistic laser spectroscopy experiments [24].

II. THEORY OF SPECTROSCOPY
This section is devoted to the general theoretical framework of spectroscopy. We first explain the model Hamiltonian, and then, we show the common formulation to describe the microwave and the laser spectroscopy. To capture essence of the present spectroscopy, we discuss physical properties of spectra in simple model cases. For simplicity, we set = 1 and k B = 1.

A. Model Hamiltonian
Before spectroscopic excitations, thermal equilibrium properties of atoms in an optical lattice is well described by the following single-band Bose-Hubbard Hamiltonian [25,26] whereĉ † g,i (ĉ g,i ) is the creation (annihilation) operator of an unexcited atom at the i th site, andn g,i is the corresponding number operator. Here, contributions of higher orbitals can be neglected when we consider the low energy properties. We note that higher orbitals have a role in the spectroscopy as discussed later soon. The Hubbard parameters, i.e., the interaction strength U g,g and the hopping integral J g , are evaluated by the ab initio calculations based on the second quantization using experimental parameters: a lattice constant a L and a lattice depth V 0 , which are determined from a wavelength and an intensity of the lattice laser, respectively, and a g,g a scattering length between two unexcited atoms. The chemical potential µ is determined so as to fix the total number of atoms N tot = i n g,i and V g,i is the trapping potential. In the following, for simplicity, we omit to explicitly write down µ, which can be included in a global shift of V g,i .

B. Spectroscopy
In the microwave and the laser spectroscopy, excitation processes caused by an external electro-magnetic field are generally described by ΓÔ ex e iωt+iKex·r [16], where Γ, ω, and K ex are a non-dimensional normalized amplitude, an angular frequency, and a wavevector of the external field, respectively, and t and r are time and position. HereÔ ex is an excitation operator defined as follows. For convenience, we define α ρ αÔex,α as a second quantization ofÔ ex e iKex·r , and whereĉ eα,i is the annihilation operator of an excited atom in the α-th orbital at the i th site of the position r i . Note that k ex is a reduced wavevector in the first Brillouin zone defined by k ex ≡ K ex + G, where G represents any reciprocal vectors with e iG·ri = 1. An excitation matrix ρ α is defined by where W α (r − r i ) is the α-th Wannier orbital at the ith site for the excited atoms, and W 1 (r − r i ) is that for the unexcited atoms (i.e., α = 1). Here, |ρ α | 2 represents the probability that the orbital of atoms changes from the lowest to the α-th orbital during excitations. The orthogonality of the Wannier orbitals assures a condition α |ρ α | 2 = 1. In this paper, we neglect the probability that atoms are excited to the different lattice sites (i.e., inter-site excitation), because it is exponentially smaller than that of the onsite excitations. Namely, we assume that We focus on the weak excitation limit under the condition of |Γ| ≪ 1. The excitation spectra can be formally given by I(ω) = α |ρ α | 2 I α (ω), and where |n is the eigenstate of Hamiltonian H g in Eq. (1) with energy E n , and Ω(= −T ln n e −En/T ) is the grand potential. Note that the conservation law of the number of excited atoms allows us to decompose I(ω) into the sum of I α (ω). The excited state |n ′ is the eigenstate of HamiltonianĤ ≡Ĥ g +Ĥ e +Ĥ ge , and E ′ n is its energy. HereĤ e andĤ ge are given bŷ where J eα is the hopping integral of excited atoms in the α-th orbital, and U g,eα is the onsite interaction between the first orbital unexcited and the α th orbital excited atoms. Note that the interaction between two excited atoms U eα,eβ can be reasonably neglected in the limit of weak excitations. ∆ eα represents the energy difference between the unexcited atoms in the lowest orbital and the excited atoms in the α th orbital. We can always set ∆ e1 = 0 by appropriately choosing the origin of the spectral frequency. The spectral intensity proportional to |Γ| 2 is determined so as to satisfy the integral condition I(ω)dω = const., when we compare our analyses with the experimental observations. Thus, we can neglect the quantitative aspect of Γ by setting Γ = 1 without loss of generality.

C. Sum rules
We discuss the moment expansions of the spectral function by following the previous studies [16][17][18]. In general, spectra given by Eq. (4) should satisfy the following sum rules in terms of M (n) α the n-th order moment: where [[[Ô † ex,α ,Ĥ],Ĥ, · · · ] denotes the n-times commutator betweenÔ ex,α andĤ. These relations indicate that certain statistic quantities in thermal equilibrium have a relation with some properties of spectral functions reflecting a dynamical response of the system. For examples, by considering up to the second order moments, the spectral mean valueω α is given byω α = M Naively, we can stress that ω α and σ α determine a spectral peak position and its width, respectively, and the relations defined by Eq. (7) can describe the physical origin of the peak shift and also broadening caused by the many-body effects. In what follows, to discuss these important spectral properties, we analyze the sum rules up to the second order.
We can derive the following expressions with respect to the corresponding sum rules. The zeroth order is given by The first order moment is written as where d represents a summation over the adjacent sites, and the second order is The expressions in Eqs. (8)-(10) include the onsite multi-body correlation functions G ℓ,i ≡ (ĉ † g,i ) ℓ (ĉ g,i ) ℓ , where G 1,i is equivalent to the averaged number of atoms n i ≡ n g,i . The hopping Hamiltonian yields the intersite correlations such as G i,j ≡ ĉ † g,iĉ g,j , where G i,i = n i , and G 2,i,j ≡ ĉ † g,iĉ g,iĉ † g,iĉ g,j . In the second (or higher) order of moment, the hopping Hamiltonian yields also the onsite one-body correlation G 1,i (= n i ) [see the last term in Eq. (10)]. Note that this type of terms is caused by a round-trip hopping process. These terms are proportional to J g − J eα e ikex·d , where e ikex·d describes the momentum transfer from the external field to atoms. The momentum transfer can be regarded as a back action of the measurements, which plays an important role in spectra.

D. Physical meaning of the spectral deviations and spectral mean value
It is convenient to discuss the physical meaning of spectra, which can be figured out from the sum-rule approach. The spectral mean valueω α can be written as the sum of δU α , δV α , δ∆ α and δJ α , which are spectral energy shifts caused by the effects of interaction, trapping potential, band gap, and hopping, respectively. As discussed in the previous study [16], δU α ≡ (U g,eα − U g,g ) i G 2,i /N tot is a collisional energy shift. Two terms, δV α ≡ i (V eα,i − V g,i )n i /N tot and δ∆ α ≡ i ∆ eα n i /N tot , are related to the statistical average of the number of atoms. Here, δ∆ α reduces to a constant ∆ eα having no connection to any thermodynamic quantities, while δV α describes the effects of the inhomogeneity of the system. The hopping energy shift δJ α ≡ i d (J g − J eα e ikex·d )G i,i+d /N tot is related to the inter-site correlation. The Bloch band picture makes physical meanings of this term clear. We define the kinetic energy shift δK α ≡ δJ α + δ∆ α , which can be rewritten as k (ε eα,k+k ex −ε g,k )n g,k /N tot , where n g,k = i,j e ik·(ri−rj ) G i,j is the momentum space distribution, and ε g,k and ε eα,k are the dispersions of the unexcited atoms in the lowest orbital and of the excited atoms in the α-orbital, respectively. Because of the momentum conservation law, a wavevector k ex is transfered from an external field to atoms. It turns out that δK α becomes large for a large k ex . Note that the discrete translational symmetry imposes k ex to be a reduced wavevector in the first Brillouin zone.
We next discuss the standard deviation σ α related to the second order moment as σ 2 Roughly speaking, the standard deviation characterizes broadening of spectra. For clear vision, we now focus on a uniform system by setting V eα,i −V g,i = 0. Considering the physical origin of the deviation, we can rewrite σ α as σ 2 , and that caused by kinetic terms σ K,α is defined as and that originating from the cross terms σ √ UK,α is given by σ 2 In the following, we consider two simplified model cases and calculate a mean valueω α and a deviation σ α , so as to discuss what are the physical origins of spectral peak shifts and broadening. Here we focus on uniform systems at zero temperature for simplicity, and inhomogeneous systems at finite temperatures will be discussed in Sec. III. We first consider SF states at zero temperature and use the following simple mean-field approximation. We assume that ĉ g,i is finite and is given by a classical complex number c i , which leads to G ℓ,i ∼ |c i | 2ℓ , ĉ † g,iĉ g,j ∼ c * i c j , and n g,k ∼ δ k,k0 N SF , where N SF is the number of the condensed SF atoms, and k 0 is a wavevector at the bottom of a dispersion ε g,k (usually k 0 = 0 for a positive J g ). We here also assume that almost all atoms are in the condensed state N SF ∼ N tot . The spectral mean value is now given byω α = (U g,eα − U g,g )n i + (ε eα,k0+kex − ε g,k0 ). It is consistent with the previous study [16], while the additional kinetic energy shift is found. The deviation now reduces to σ 2 Interestingly, we find that the deviation has a contribution of an interaction term σ U,α only. Here because of the subtraction in M KU,α are canceled out under the pure-condensation condition n g,k ∼ δ k,k0 N tot . The number of the coherent SF atoms is indefinite, and number fluctuations ∆n i of such a coherent state are written as √ n i , where (∆n i ) 2 = n 2 g,i − n 2 i . Thus, we can conclude that the spectral deviation of coherent SF states is connected to number fluctuations in thermal equilibrium as We next consider MI states with m atoms in each site at zero temperature. For MI states, it is reasonable to set inter-site correlations G i,j for i = j to be zero. The spectral mean value is now given byω α = (U g,eα − U g,g )(m − 1) + ∆ eα , which is equivalent to the previous study [16]. The spectral deviation reduces to In the same manner as the above, the correlation-induced deviation σ U,α cancels out as σ 2 UK,α is also zero because of negligible inter-site correlations. For MI atoms at zero temperature, the number of atoms is definite, and there are no number fluctuations ∆n i = 0. On the other hand, the phase and the momentum is indefinite, and thus a spectral broadening of MI atoms is caused by kinetic fluctuations described by σ K,α . By executing the k-summation in σ K,α with a constant momentum distribution n g,k (= m) of uniform MI states, we obtain σ 2 This deviation can be connected to quantum fluctuations resulting from the round-trip hopping process given in the last term in Eq. (10).
As demonstrated by the previous experiments [11], the present spectroscopy has an ability to distinguish the different number states, when |U g,eα − U g,g | is large enough. This feature of spectra can be explained by the first order sum rule,ω α = (U g,eα − U g,g )(m − 1), as discussed in the previous study [16]. The above second-order sum rule approach further indicates that number fluctuations of atoms in thermal equilibrium ∆n i play an important role in such number-resolving spectroscopy. In fact, a spectral deviation σ α of coherent SF atoms can be deter-mined from number fluctuations: σ α = |U g,eα −U g,g |∆n i . Even for MI atoms with ∆n i = 0, a spectral deviation is finite owing to kinetic fluctuations, which is attributed to the indefinite phase and momentum in a reflection of the definite number and position of MI atoms. This deviation of the number definite states is not related to any thermodynamic quantities, and thus this constant σ α of d |J g − J eα e ikex·d | 2 is the intrinsic lower limit of a spectral linewidth (see Sec. III B). On the other hand, the deviations of general states with both phase and number fluctuations are given by the summation, We find that the spectral deviations of a specific ground state (SF and MI) with a definite quantity (phase and number) are characterized by fluctuations resulting from the conjugate indefinite quantity (number and phase, respectively). It should be noted that the second-order sum rule makes clear the fact that the spectral measurements are governed by the uncertainty principle.
These physical properties of spectra mentioned above are figured out from the general principle of weak excitation spectroscopy with given HamiltonianĤ and an operatorÔ ex,α . For example, an excitation operator O ex,α characterizes the intrinsic spectral broadening of where a momentum transfer of k ex , which is a back action of the measurements, determines a quantitative aspect of a spectral width. Bose-Hubbard Hamiltonian H g includes kinetic and interaction terms, and the competition between these two conjugate terms is the origin of the SF-MI transitions. The second-order sum rule approach clarifies that spectral deviations reflect the completely different properties of these two conjugate states, SF and MI. We can thus conclude that the present spectroscopy will be a sensitive tool for detecting the SF-MI transitions. Note that the first order sum rule approach is insufficient to clarify these important features of spectroscopy. However, the above simplified discussions cannot be straightforward applied to the finite temperature spectra. In the next section, we thus propose a two-mode approximation to numerically calculate spectra that satisfy the sum rules.

III. METHODS
In this section, we provide a numerical method for calculating spectra at finite temperatures. We first explain the finite temperature Gutzwiller approximation [23], which allows us to efficiently obtain the thermodynamic quantities in Eqs. (8)- (10). We next provide a two-mode approximation to numerically calculate finite temperature spectra in inhomogeneous systems. At the end of this section, we compare our method with the previous formulations [16,17].

A. Finite temperature Gutzwiller approximation
The Gutzwiller approximation allows us to efficiently analyze the thermal equilibrium properties described by the Bose-Hubbard Hamiltonian in Eq. (1). This is a mean-field approximation considering up to the firstorder collection in terms of J g and well describes the SF-MI transitions in high dimensional systems. Here, H g is then approximated by a set of the effective local HamiltonianĤ loc = iĤ loc,i witĥ Using exact diagonalization, we can numerically calculate statistical quantities such as c i ≡ ĉ g,i at finite temperatures [23].
As discussed in Sec. II D, a finite c i effectively describes the Bose-Einstein condensates (BEC) within the meanfield approximation. Here, the number of condensed SF atoms in each site can be defined by n SF,i = |c i | 2 . Both thermal fluctuations and interactions cause coexisting of condensed SF and uncondensed NS such as MI and normal fluid (NF). The annihilation operator of NS atoms at the ith site is effectively given byĉ NS,i =ĉ g,i − c i , where NS atoms satisfy always a condition ĉ NS,i = 0. The number of NS atoms n NS,i ≡ ĉ † NS,iĉ NS,i can be written as n NS,i = n i − n SF,i . This leads to the following reasonable condition: N tot = N SF + N NS , where N SF(NS) = i n SF(NS),i . The total number of atoms is given by the sum of the total number of SF and NS atoms.
It is useful to briefly explain how to calculate the thermal quantities in the moments Eqs. (9) and (10). The onsite correlation functions G n,i can be calculated straightforwardly by diagonalizing the effective local Hamilto-nianĤ loc,i . On the basis of this local approximation, the inter-site correlation is described by G i,j ∼ ĉ † g,i ĉ g,j (= c * i c j ) for i = j. The higher-order inter-site correlations G 2,i,i+d also reduce to ĉ † g,iĉ g,iĉ † g,i ĉ g,j . These expressions mean that the inter-site correlations of NS atoms ĉ † NS,iĉ NS,j for i = j are approximately set to be zero, and thus, two kinds of NS states, MI and NF, are dealt with approximately in the same way. We should note that the MI states appearing at lower temperatures can be characterized by focusing on the creation of the Mott shell structures and also the suppressed entropy per site [23]. We can thus effectively calculate that thermal fluctuations cause the MI-NF crossover within this local approximation.

B. Two-mode approximation
Next, we provide a two-mode approximation that helps us to calculate spectra at finite temperatures. We assume that I α (ω) in Eq. (4) can be decomposed into two com-ponents resulting from the contributions of SF and NS atoms: Two types of uncondensed NS (MI and NF) states appear at finite temperatures. As mentioned in Sec. III A, within the Gutzwiller approximation, NF states are approximately dealt with in the same way as MI states based on the local Hamiltonian picture. As discussed in Sec. II D, at zero temperature, spectra of coherent SF atoms show completely different properties by comparing with those of MI atoms. Note that the special characteristics of spectra of SF atoms result from the phase coherence caused by BEC. We thus deal with spectra of SF atoms in different way to two types of NS atoms. On the basis of the two-mode approximation, we reconsider the sum rules for the spectral moments: The sum rules up to the second order (i.e., up to n = 2) provide the following relations: whereω NS α andω SF α are the spectral mean value, and σ NS and σ SF are the spectral standard deviation for the NS and SF spectra, respectively. The zeroth order sum rule in Eq. (13) simply offers the condition associated with the total number of atoms, which is always satisfied within the Gutzwiller treatment as mentioned in Sec. III A. On the other hand, the first and second order sum rules require the balance conditions between I NS α (ω) and I SF α (ω), and these conditions allow us to properly calculate spectra.

Spectra of uncondensed normal state atoms
In what follows, we discuss the properties of I NS α (ω) and I SF α (ω) at finite temperatures, separately. Here we begin with I NS α (ω) by assuming that N SF = 0, and accordingly I SF α (ω) vanishes. We also assume that the intersite correlations are negligible G i,j for i = j by comparing to the on-site correlations n i . The density matrix of such a localized state is given by i ( m e −Em,i/T |m i m| i ) at finite temperatures, where E m,i = U g,g m(m − 1)/2 − V g,i m is the energy of the local number state |m i . The spectra can be obtained in a form of the exact representation: The spectral weight w α,i,m and the peak position p α,i,m are given by where e −(Em,i−Ωi)/T is the Boltzmann factor of the number state |m i , and Ω i (= −T ln m e −Em,i/T ) is the grand potential in the i th site. Note that, even for the uniform systems, spectra at finite temperatures have multi-peak structures depending on the thermal distributions of the number states |m i described by e −(Em,i−Ωi)/T . We now discuss that the zeroth and the first order sum rules are always satisfied in the above expression in Eq. (16) when G i,j ≪ n i . By using Eqs. (17) and (18), we obtain M These facts suggest that I NS α (ω) in Eq. (16) reproduces the zeroth and the first order moments in Eqs. (8) and (9) when we can neglect the term proportional to the inter-site correlations [the last term in Eq. (9)].
In contrast, the second order sum rule is not straightforward. When G i,j ≪ n i , almost all terms are reproduced in the same way as the above. Namely, M (2) α = i,m w α,i,m p 2 α,i,m is equivalent to the first four terms in Eq. (10). However, we cannot reproduce one of the terms in Eq. (10), which is the round-trip hopping term given by i d |J g − J eα e ikex·d | 2 n i . This means that, even though the inter-site correlations are negligible, quantum fluctuations resulting from the roundtrip hopping broaden the spectral width of each peak in Eq. (16). Namely, the sum rule requires that the delta function δ(ω) in Eq. (16) should be replaced with a certain function with a finite spectral width. We here use a Gaussian function, and I NS α (ω) is now given by where γ α is the spectral width defined by γ 2 α = d |J g − J eα e ikex·d | 2 . We comment that a Lorentzian function is not suitable for the substituting function, because the second order moment does not converge: ∞ −∞ dωω 2 γ/π(ω 2 + γ 2 ) → ∞. The extended representation in Eq. (19) with Eqs. (17) and (18) properly satisfies the sum rules up to the second order when G i,j ≪ n i . In the same way as the above, we can straightforwardly confirm that the zeroth and the first order sum rules are satisfied.
The second order M (2) α is extended as follows: i,m w α,i,m . The first term i,m w α,i,m p 2 α,i,m is equivalent to the second order moment obtained from the original representation in Eq. (16). The additional term γ 2 α i,m w α,i,m (= γ 2 α N NS ) properly describes the last term in Eq. (10).
We here estimate the magnitude of γ α that characterizes an intrinsic spectral broadening caused by hoppinginduced quantum fluctuations. For simplicity, we consider α = 1 and set J g ∼ J e1 , which leads to γ 2 1 = 2zJ 2 e1 [1− d cos(k ex ·d)/z], where z is the number of the neighboring lattice sites (z = 6 in the cubic lattice). For k ex ∼ 0, γ 1 reduces to zero. For k ex ∼ (π, π, π), γ 1 takes a maximum 2 √ z|J e1 | = W 1 / √ z, where W α = 2z|J eα | is a bandwidth of the α th orbital. We next consider higher orbitals α = 1 by assuming |J g | ≪ |J eα |, and then we obtain γ α = √ z|J eα | = W α /2 √ z. Simply put, the kinetic spectral broadening is proportional to the bandwidth, γ α ∝ W α , where the wavevector conservation law determines the proportionality coefficient ranging from 0 to 1/ √ z depending on k ex . Before closing the discussions on I NS α (ω), we consider the validity of the condition G i,j ≪ n i . For J g = 0, this condition is exactly satisfied: G i,j = 0 for i = j. For a finite but small J g (≪ U g,g ), where MI states will appear at low temperatures, the effects of interactions strongly suppress the inter-site correlations. Large potential differences strongly suppress the inter-site correlations (e.g., |V g,i −V g,j | ≫ J g ), and thermal fluctuations also decrease G i,j . We thus expect that the NS atoms in the realistic systems with interactions and trapping potential at finite temperatures will satisfy well the condition of small G i,j (≪ n i ). This condition is equivalent to flattened momentum distributions n g,k , which can be confirmed in experiments by using the time-of-flight measurements with the projection onto the first Brillouin zone [27].

Spectra of the superfluid atoms
Next, we consider the opposite limit, N SF ≫ N NS , where we neglect I NS α (ω). Taking account of the physical properties of BEC, we assume that I SF α (ω) can be described by the following single peak structure: where a spectral peak positionω SF α and a spectral width σ SF α are determined from the sum rules in Eqs. (14) and (15), respectively. For such a single peak structure, the deviation coincides with the spectral width. This single peak assumption may be oversimplification. It should be noted that we carefully take account of the spectral broadening caused by number fluctuations, which allows us to reasonably use this simple assumption.
To compare I SF α (ω) with I NS α (ω), here we mention again the properties of I NS α (ω) at finite temperatures. As shown in Eq. (19), I NS α (ω) has many peaks at positions p α,i,m (∝ U g,eα −U g,g ) with a spectral width of γ α (∝ W α ). Note that peak positions and a width are usually determined from the different energy scales. Thermal fluctuations change relative spectral weights and also increase the number of spectral peaks. This multi-peak structure is an essential feature of I NS α (ω) at finite temperatures, which can be properly dealt with in Eq. (19). In contrast, as discussed in Sec. II D, for SF atoms, a spectral mean value is proportional to the average number of atoms ω SF α ∝ (U g,eα − U g,g )n, and a spectral deviation is given by σ SF α ∼ |U g,eα − U g,g | √n . This fact suggests thatω SF α and σ SF α will be usually comparable. A main role of thermal fluctuations is a decrease in the number of condensed atoms N SF . We thus expect that the essence of I SF α (ω) at finite temperatures can be captured by a single peak broadened by large number fluctuations. The effects of decrease in N SF on the spectra are considered via calculations on the sum rules in Eq. (13), and also Eqs. (14) and (15), which properly describe decreases in a peak height, a peak shift and broadening, respectively.

Spectra for coexisting region
We next explain the formulation for the middle region, N SF = 0 and N NS = 0, where NS and SF atoms coexist. Here, I SF α (ω) and I NS α (ω) are separately calculated, and I NS α (ω) is the first. We assume that NS atoms are affected by a mean-field potential resulting from interactions with SF atoms. Thus, we here consider the following effective local Hamiltonian of NS atoms excluding SF atoms: where the potential V SF,i describes effectively mean-field interactions between SF and NS atoms given by V SF,i = 2U g,g n SF,i + δµ i . To impose a self-consistent condition n NS,i Ĥ NS = n g,i − | ĉ g,i | 2 , we further define a chemical potential shift δµ i , where · · · Ĥ NS is the statistical average at thermal equilibrium defined by the Hamilto-nianĤ NS , while · · · is that defined by the localized Hubbard Hamiltonian in Eq. (11). We note that δµ i ∼ 0 for n SF,i ≫ n NS,i or n SF,i ≪ n NS,i , because the meanfield treatment is appropriate for these dilute regions.
1. We first calculate thermal equilibrium states of Hubbard Hamiltonian based on the finite temper-ature Gutzwiller approximation. We use exact diagonalization to solve the localized Hamiltonian in Eq. (11) at finite temperatures. We obtain the moments M (n) α from Eqs. (8)- (10), and other statistical quantities such as n SF,i and n NS,i .
2. Next, we calculate I NS α (ω) in Eq. (19) by exactly diagonalizing the effective Hamiltonian in Eq. (21). After that, we can directly calculateω NS α and σ NS α from I NS α (ω). For the self-consistent condition mentioned above, we need n i and c i , which should be obtained in the previous process 1.
In this way, based on the sum rule approach, completely different features of I NS α (ω) and I SF α (ω) are properly dealt with, and multi-peak structures resulting from finite temperature effects and coexisting of SF and NS atoms will be taken account of precisely.

C. Role of the sum rules in the two-mode approximation
We finally discuss what a role the sum rules play in the present calculation procedure. The sum-rule approach combined with the two-mode approximation provides us with the reasonable relationship between I NS α (ω) and I SF α (ω). As discussed below,ω SF α and σ SF α can be determined reasonably within the level of the mean-field treatment.
We here consider the first order moment M (1) α (= N NSωNS,α + N SFωSF,α ) and focus on the effects of interactions [the first term in Eq. (9)]: (U g,eα − U g,g ) i G 2,i , and other terms are neglected for clarity. Using a definition ofĉ N S,i =ĉ g,i − c i , we can rewrite G 2,i as G NS 2,i + 4n NS,i n SF,i + n 2 SF,i , where G NS 2,i ≡ ĉ † NS,iĉ † NS,iĉ NS,iĉNS,i . Three terms, G NS 2,i , n NS,i n SF,i , and n 2 SF,i , represent the two-body correlations between two NS atoms, between NS and SF atoms, and between two SF atoms, respectively. As discussed in Sec. III B 1, the first order moment of I NS α (ω) is easily obtained: where the latter term results from the mean field potential V SF,i in Eq. (21). Consequently, the sum rule in Eq. (14) allows us to determine the correlation term in the first moment of I SF α (ω): i (n 2 SF,i + 2n SF,i n NS,i ).
The SF-NS correlation 4n NS,i n SF,i is shared equally between I NS α (ω) and I SF α (ω). We now again consider the uniform system to compare this expression with those in Sec. II D. The collisonal energy shift for the SF atoms is now given byω SF α = (U g,eα − U g,g )(n SF,i + 2n NS,i ). The first term (U g,eα − U g,g )n SF,i is equivalent to that of the pure SF atoms as discussed in Sec. II D and also in the previous study [16], while the second term is the additional contribution originating from the mean-field NS-SF interactions.
In the same way as the above, we can obtain the correlation terms in the second order moment of I NS α (ω) and For the uniform system, (σ SF α ) 2 is written as (U g,eα − U g,g ) 2 [n SF,i + n SF,i n NS,i + G NS 2,i + 4(∆n NS,i ) 2 ], where (∆n NS,i ) 2 = n 2 NS,i − n 2 NS,i . Note that ∆n NS,i ∼ 0 at low temperatures. The first term corresponds to number fluctuations of the SF atoms, and the other terms suggest that the NS-SF correlations enhance number fluctuations of the SF atoms and further broaden I SF α (ω).

D. Comparison with the previous studies
Here, to discuss difference between the present and the previous treatment, we summarize the previous formulation [16], which has been successfully applied to the microwave spectroscopy experiments [11]: This approximation satisfies the zeroth order sum rule M α = i G 2,i (U g,eα − U g,g ), while the second order sum rule is not satisfied at all. Note that, without a consideration of the deviation σ α , the spectral broadening is effectively described by the site-dependent G 2,i and n i resulting from the inhomogeneity of the system. When k ex = 0, |J g − J eα e ikex·d | = 0, V eα,i − V g,i = 0, T = 0 and the coexisting of NS and SF states is neglected (either N NS = 0 or N SF = 0), the spectral positionω SF α obtained from Eq. (9) or p i,α,m in Eq. (18) and the corresponding term (U g,e1 − U g,g )G 2,i /n i in Eq. (22) are equivalent with each other. Thus, within the first order approximation, our method is consistent with Eq. (22) in the following two limits; pure SF or MI phases at zero temperature. We note that, at least, the conditions k ex = 0, |J g −J eα e ikex·d | = 0, and V eα,i −V g,i = 0 are well satisfied in the microwave spectroscopy (see Sec. IV A).

IV. NUMERICAL SIMULATIONS
In this section, we show numerical results calculated by considering realistic parameters of the following two experiments; the microwave-spectroscopy of 87 Rb atoms and the laser spectroscopy of 174 Yb atoms. We first explain the parameters and then discuss obtained results.

A. Parameters
We first point out intrinsic differences of the microwave and the laser spectroscopy, and provide parameters used in the calculations. A length of microwave is much longer than a lattice constant a L (= λ L /2) and a lattice laser wavelength λ L (e.g., of 1064 nm), and as a result, a wavevector k ex can be set zero. Therefore, the parameter region of the microwave spectroscopy corresponds to the perfect Lamb-Dicke regime, where excitation matrices ρ α in Eq. (3) are given by ρ 1 = 1 and ρ α =1 = 0. Here the kinetic energy shift δK α and deviations σ K,α and σ √ UK,α are also negligible. In contrast, for the laser spectroscopy, a wavelength of the excitation laser λ ex (e.g., of 507 nm) is comparable to λ L (e.g., of 532 nm). Namely, k ex is the same order as a lattice wavevector 2π/a L . Thus, orbitalchanging excitations and kinetic contributions will be important in the laser spectroscopy.
The microwave spectroscopy of 87 Rb atoms uses an excitation between different hyperfine states [11], while the laser spectroscopy of 174 Yb atoms uses an excitation between different electron configurations, 1 S 0 and 3 P 2 states [20][21][22]. Rb atoms are trapped by the combination of optical and magnetic potential, and Yb atoms are trapped with optical potential. For the microwave spectroscopy, we use a harmonic trapping potential V g,i ∝ C x x 2 i + C y y 2 i + C z z 2 i , where the curvatures C x , C y , and C z are determined from the experimental parameters so as to reproduce the bottom of trapping potential [11], while for the laser spectroscopy, we use an anharmonic potential by carefully considering the laser configurations in experiments [24]. Note that the trapping potential of two hyperfine states was set to be the nearly same [11], while 1 S 0 and 3 P 2 states of Yb atoms are trapped in the different potential due to the greatly different polarizability. Thus, V eα,i − V g,i can be set zero for the microwave spectroscopy, while it is finite for the laser spectroscopy.
Differences in scattering lengths a g,e −a g,g are -0.13 nm and -30 nm for the microwave and laser spectroscopy, respectively. Both of them are negative, so that differences in the interaction strengths U g,eα −U g,g are also negative. In addition to the spectral broadening caused by quantum fluctuations (see Sec. II D), we consider linewidths of the excitation laser of about 1 kHz and of the microwave of about 5 Hz, including the Fourier width of the excitation pulse.

B. Comparison between two spectroscopy for the deep lattice
We first discuss spectra in a deep lattice and compare two kinds of spectroscopy. Figure 1 shows the spectra calculated at a temperature T of 100 nK. The other pa-  rameters in the microwave spectroscopy are V 0 = 35E r and N tot = 10 5 , and those in the laser spectroscopy are V 0 = 15E r and N tot = 2.2 × 10 4 , where E r is the recoil energy. A spectral peak appearing at ω = 0 always means that the m = 1 number state (|m = 1 ) is excited without orbital changing, because the origin of spectra is renormalized by setting ∆ e1 = 0. Since U g,eα −U g,g is negative for both spectroscopy, peaks of |m with m ≥ 2 appear orderly in the region of ω < 0. The orbital-changing excitation requires a large positive bandgap energy ∆ e2 . Thus, spectra in the laser spectroscopy show some peaks in ω > 0, which have the similar characteristics to those in ω < 0 but have the small intensities because of a small excitation probability |ρ 2 | 2 (∼ 0.1|ρ 1 | 2 ).
Next, we discuss a spectral peak width. For a deep lattice, spectral broadening is mainly attributed to the effects of inhomogeneity. Equation (9) shows that potential energy difference (V eα,i − V g,i )n i /N tot causes just a peak shift. However, due to the inhomogeneity, a variation of (V eα,i − V g,i )n i /N tot for different i effectively induces broadening of spectra. For the laser spectroscopy, an energy scale of this variation is estimated to be about 1 kHz, which is consistent with the obtained spectral features in Fig. 1. In the present parameter region, broadening effects resulting from the hopping terms are negligible by comparing with this inhomogeneous broadening. On the other hand, since |V eα,i − V g,i | ∼ 0 for the present microwave spectroscopy, a width of each peak nearly equals to a linewidth of the microwave pulse.
C. Lattice depth dependence of spectra in the microwave spectroscopy Next, by focusing on the microwave spectroscopy, we discuss how spectra change as the lattice depth varies. Figure 2 shows the spectra I(ω) calculated with V 0 = 5E r , 10E r , and 25E r for N tot = 10 5 at T = 25 nK. For V 0 = 25E r , the number of SF atoms N SF is zero, and a discrete peak structure resulting from I NS (ω) appears. In contrast, for V 0 = 5E r , almost all atoms are in condensed states N SF ∼ N tot and N NS ∼ 0. From discussions in Sec. II C, we can naively expect that the spectra of SF atoms I SF (ω) are centered at around (U g,eα −U g,g )n, and a peak width is determined from |U g,eα −U g,g | √n , wheren is the averaged number of atoms. From these features, we can estimaten ∼ 2. For a middle region V 0 = 10E r , we find an asymmetric broadened peak structure. We note that this characteristic asymmetric structure can be attributed to the coexisting of NS and SF atoms. Figure  2 also shows spectra of SF atoms I SF (ω) and those of NS atoms I NS (ω). The spectra of SF atoms I SF (ω) are centered at around 2.5(U g,eα − U g,g ). On the other hand, I NS (ω) has a large intensity at around ω ∼ 0 corresponding to the |m = 1 excitation. Here we find no discrete peak structures, because |U g,eα −U g,g | is smaller than the microwave linewidth. The sum of two spectra yields the characteristic asymmetric spectra.
The obtained spectra for all three parameters in Fig.  2 capture essential features of those observed in experiments [11]: The experimental spectra show the symmetric single peak structure for V 0 = 5E r , and the asymmetric broadened peak structure, or an overlapped doublepeak structure, for V 0 = 10E r , and the discrete peak structure for V 0 = 25E r . Our assumptions, e.g., the single peak assumption for the SF atoms, and the twomode approximation, appropriately reproduce characteristic structures of the spectra seen in experiments.
D. Temperature dependence of spectra in the microwave spectroscopy Next, we discuss the temperature dependence of spectra by focusing on the microwave spectroscopy again. In Fig. 3, we first show the spectra for V 0 = 35E r at T = 100  nK, 50 nK, and 10 nK. From Fig. 3, we find that spectral peak height of each |m state varies depending on temperature, which indicates the change in thermal distributions of the number states |m . At low temperatures (T =10 nK), the m = 2 peak is highest. On the other hand, at higher temperatures (T =100 nK and 50 nK), the m = 1 peak is dominant. To clearly discuss this behavior, in Fig. 4, we show the in-situ column-density distributions of the number states |m , which can be calculated by exact diagonalization for the localized Hamiltonian in Eq. (11). At higher temperatures (T =100 nK), the |m = 1 state spreads widely, while larger-m states prefer center of the potential. Thus, N m the total number of |m monotonically decreases with increasing m. This behavior is consistent with those seen in spectra in Fig. 3. In contrast, at low temperatures (T =10 nK), the Mott shell structures develop, and as a result, the number state distributions for m ≤ 2 show a dip structure as shown in Fig. 4. Here a three dimensional shell structure is projected onto the dip structure in the columndensity distributions. Thus, N m becomes nonmonotonic and N m=2 becomes maximum, leading to a large inten-sity of the peak of |m = 2 in spectra (see Fig. 3). These features suggest that the present spectroscopy can be a thermometer.
We finally compare these results with the experimental observations in Ref. [11], in which spectra and also in-situ number-state distributions have been measured. We find that both spectra and distributions agree well with the calculations at higher temperatures than those at lower temperatures. It thus suggests that temperature will be an important parameter to discuss the quantitative aspect of experiments. Although nearly pure BEC can be created before loading atoms into a lattice, the loading usually induces a certain amount of heating [23]. In the present calculations, a harmonic trapping potential is used by focusing on the bottom of the experimental trapping potential [11], whereas in Ref. [16], by carefully considering the anharmonicity of trapping potential, experimental spectra can be reproduced even at zero temperature. The calculations considering both anharmonic potentials and thermal fluctuations may be required for more quantitative agreements.

V. SUMMARY
In summary, we theoretically investigate the microwave and the laser spectroscopy on the bosonic Hubbard systems. We first discuss the sum rules of spectra up to the second order, which can be derived from the general principle of spectroscopy. This principle provides various useful information on the physical properties of spectra. The spectra of superfluid states with phase coherence is broadened by the many-body effects, and its broadened width can be characterized by number fluctuations in thermal equilibrium. In contrast, the spectra of the number definite Mott insulating states are broadened by quantum fluctuations caused by tunneling effects.
We next propose a two-mode approximation to calculate spectra at finite temperatures. This approximation assumes that spectra can be decomposed into two contributions originated from Bose-Einstein condensates and uncondensed normal states. Our method is built up by considering the spectral characteristics figured out from the sum rules, so that the multi-peak structures resulting from coexisting of superfluid and normal states at finite temperatures can be successfully dealt with.
Finally, by combining the two-mode approximation with the finite temperature Gutzwiller approximation, we numerically calculate spectra by considering realistic experimental parameters of the microwave and the laser spectroscopy. We find that our method can reproduce the essential features of spectra in experiments. We also discuss the lattice depth dependence and the temperature dependence of the microwave spectra. These results clarify that the present spectroscopy can be sensitive tools for investigating the quantum phase transitions at finite temperatures.