Emergence of damped-localized excitations of the Mott state due to disorder

A key aspect of ultracold bosonic quantum gases in deep optical lattice potential wells is the realization of the strongly interacting Mott insulating phase. Many characteristics of this phase are well understood, however little is known about the effects of a random external potential on its gapped quasiparticle and quasihole low-energy excitations. In the present study we investigate the effect of disorder upon the excitations of the Mott insulating state at zero temperature described by the Bose-Hubbard model. Using a field-theoretical approach we obtain a resummed expression for the disorder ensemble average of the spectral function. Its analysis shows that disorder leads to an increase of the effective mass of both quasiparticle and quasihole excitations. Furthermore, it yields the emergence of damped states, which exponentially decay during propagation in space and dominate the whole band when disorder becomes comparable to interactions. We argue that such damped-localized states correspond to single-particle excitations of the Bose-glass phase.


Introduction
Ultracold quantum gases in optical lattices offer a unique possibility to investigate strongly interacting many-body systems [1,2]. When loaded into a deep lattice potential at commensurate fillings, a cloud of ultracold bosons undergoes a quantum phase transition from a superfluid to a Mott insulating ground state [3,4,5]. The latter state is characterized by a gap for quasiparticle and quasihole excitations due to onsite repulsive interactions. In the additional presence of disorder, rare condensate regions emerge inside an insulating background characterizing a Bose-glass ground state [3,6,7,8,9]. Such an effect of spatially random fields can be mimicked by a timealternating potential where a nonequilibrium granular condensate appears [10]. The phase transitions between theses states was studied both numerically, with Monte-Carlo simulations [11,12,13,14,15,16,17,18,19] and stochastic [20,21] as well as local [22] mean-field techniques, and analytically, where mean-field theory [23,24,25] and field theoretical methods [26] were applied. Although the nature of the excitations of the Mott state has been the subject of extensive investigation [27,28,29,30], the effect of disorder on their energy spectrum remains so far unclear. Furthermore, a concrete characterization of the Bose-glass excitation spectrum is still lacking. Here we demonstrate how the properties of the excitation spectrum in the disordered case can be obtained from the spectral function A(k, ω), which generalizes the concept of dispersion relations [31,32,33].
In a perfect lattice, excitations of the Mott state with well-defined wavevector k are eigenstates of the underlying Hamiltonian with undamped propagation and finite effective mass inversely proportional to their tunneling energy [34,35]. Such states are associated with Dirac distribution peaks of the spectral function centered at the corresponding energy ℏω + 0 (k), represented by the dashed-black line in Fig. 1. A detailed analysis of the spectral function in view of characterizing low-energy excitations of the Mott phase can be found, for instance, in Refs. [29,36,37,38]. When disorder is introduced translational invariance is destroyed. However, if one is not interested in local properties of a single realization of the random potential, one can define a disorder ensemble average which ensures that translational symmetry is recovered [39,40]. This motivates the definition of the spectral function as the imaginary part of the averaged single-particle Green's function via [41] A(k, ω) = − 1 π Im⟨G(k, ω + i0 + )⟩, where we consider the analytic continuation from the Matsubara to the real frequency domain consistent with retarded response, i.e. iω m → ω + i0 + , and ⟨· · ·⟩ stands for the disorder ensemble average. The spectral function (1) satisfies both the bosonic positive-definite property and the summation rule [32] (signω)A(k, ω) ≥ 0, also in the presence of disorder. If disorder and tunneling energy are small compared to the interactions, fluctuations of the particle density become energetically costly and the Mott state still prevails [3,7,14]. In this strongly interacting regime, one can apply a perturbative approach based on the strong-coupling expansion. This method was first proposed in the context of lattice bosons by Ref. [42], where it was also shown that its first orders provide a remarkable comparison to Monte Carlo simulations. Subsequently, this method served as the basis for the worm algorithm in the Monte Carlo simulations of Refs. [11,14,15], where all orders are taken into account. Recently, the eighth order term was calculated in Ref. [43], demonstrating that the strong-coupling approximation works extremely well in predicting the quantum phase boundary. In this work, we use a strong-coupling approach based on field-theoretical methods, such as those developed in Refs. [44,45,26], Figure 1: Qualitative sketch of the spectral function A(k, ω) for fixed a|k| ≪ 1. The sharp dashed-black peak at ω + 0 (k) corresponds to stable excitations of the clean case. Disorder shifts this peak ω + (k) (red) towards lower energies while it also creates damped states with dispersion ω r (k) which appear as a broad peak (blue) of finite width Γ(k).
to construct an expansion of the Green's function for small tunneling energy values, considering only tree-level corrections. This approach leads to a mean-field phase boundary equivalent to the one found in [23,26]. Therefore, for sufficiently small tunneling energy, such a perturbative method is sufficient to describe the low-energy excitations in the disordered case.
In what follows we report that, similarly to the clean case, stable excitations are still present and are associated with sharp peaks in the spectral function. However, as the region occupied by the Mott states in the phase diagram shrinks due to increasing disorder [23,26], the gap for these stable states decreases leading to a shift of its dispersion ω + (k) towards lower energies. Additionally, we find that a broad peak of finite width Γ(k) emerges in the spectral function due to scattering effects with the random potential. Such a broad distribution corresponds to damped states with dispersion ω r (k). The qualitative features of our results for the quasiparticle branch of the spectrum are schematically depicted in Fig. 1. For the quasihole branch of the spectrum, which occurs at negative energies, the qualitative results are analogous. For strong disorder, the damped states occupy the whole spectrum. In this limit, the dispersive nature of the excitations ceases to exist implying a large distribution in momentum, which is typical of localized states. We argue that the new set of damped-localized states correspond to single-particle excitations of the Bose glass. Moreover, we demonstrate that in the case of bounded uniform randomness the effective mass of stable excitations increases with increasing disorder. At the same time, the lifetime of the damped states also increases.
We proceed by first deriving a hopping expansion to the single-particle Green's function in Sec. 2. In Sec. 3 we discuss the characteristics of the spectral function. We then consider in Sec. 4 the case of a uniform disorder distribution and obtain the band structure for the case where the average particle density at each site is n = 0 and n = 1. In Sec. 5 we discus the spatio-temporal profile of the Green's function. Finally, we summarize the implications of our results and conclude in Sec. 6.

Hopping expansion
We start our analysis by defining the Bose-Hubbard Hamiltonian, which describes spinless bosons with short-range interactions on a lattice potential Hereâ † i andâ i satisfy standard bosonic commutation relations,n i =â † iâ i is the particle number operator, µ denotes the chemical potential. In addition, U represents the on-site interaction and J ij stands for the coupling between neighboring sites, i.e., it takes the value J only when i and j are first neighbors and vanishes otherwise. Furthermore, we assume that the local imperfections ϵ i are uncorrelated at different sites and randomly distributed over the lattice according to a distribution p(ϵ i ) bounded in the interval [−∆/2, ∆/2]. Therefore, ∆ represents the scale associated with the energy shift caused by the random potential on the lattice.
In the decoupled limit of J = 0, the HamiltonianĤ 0 = i [ U 2n i (n i − 1) + (ϵ i − µ)n i ] is diagonal in the particle number operator basis. Considering the hopping term V = − ij J ijâ † iâ j as a perturbation, we use the Dirac interaction picture representation to write the imaginary-time evolution operator as where we have set ℏ = 1. We employ the standard definition for the imaginary-time dependent operatorV (τ ) = e τĤ 0V e −τĤ 0 andT is the time-ordering operator. Thus, we define the single-particle Green's function as One can obtain a perturbative approximation to the above quantity by expanding the exponential of the operatorÛ in powers of the hoppingV and evaluating the corresponding traces with respect to the equilibrium ensemble of the unperturbed HamiltonianĤ 0 for fixed ϵ i . This approximation is analogous to the strong-coupling expansion used in Ref. [46] considering additionally frozen disorder. A similar method was applied in Ref. [42] which provided a significant comparison against Monte-Carlo calculations for the prediction of the Mott-lobes phase boundary both in the clean case and in the presence of box disorder. Each term in the expansion of (6) can be associated with a diagram corresponding to the path of an excitation which starts from site j at imaginary time zero and hops along the lattice reaching site i at imaginary time τ [44,45,26]. After combining the expansions of the numerator and the denominator of (6) order by order, one can write the remaining terms only using the connected parts of the Green's function. This result follows from the so-called linked cluster theorem [33,46,47]. Considering only the first-order correction in the hopping expansion we get where we define the unperturbed Green's function as By transforming into the Matsubara frequency domain we get that (7) can be rewritten as In the zero-temperature limit, the Matsubara frequencies become continuous and the unperturbed Green's function reads where n ∈ N 0 is the average particle density that minimizes the energy at each site. Note that the unperturbed Green's function is characterized by two simple poles that correspond to stable quasiparticle and quasihole states. The presence of disorder leads to randomly distributed shifts in the lattice potential. Such shifts act as scattering centers for the excitations. By averaging over all possible disorder configurations and considering independent scattering events for the propagation of the single-particle excitations, we obtain that the inverse of the Green's function can be written as where the disorder average is defined according to ⟨· · ·⟩ = i ∞ −∞ dϵ i · · · p(ϵ i ). By using the Fourier transform we invert (12) exactly which gives Note that this result has the form of a Dyson equation, where the self-energy is given by the dispersion J(k) = 2J D α=1 cos(ak α ) for a D-dimensional hypercubic lattice with spacing a between sites. Such a perturbative result is equivalent to the tree-level approximation used, for instance, in Ref. [26]. Thus, it amounts to considering an infinite amount of paths for an excitation created at a given site j with wavevector k to hop to a different site i and be annihilated with the same wavevector k by suffering independent scattering processes at each site. It should be noted that the partial summation result of (14) tends to underestimate fluctuations in the hopping expansion. Therefore, it is considered to be accurate only for sufficiently small tunneling energy, where second-order loop corrections can be disregarded. In the absence of disorder, it has been demonstrated that this result can effectively describe the collective excitations of superfluid ground states. Specifically, it predicts a gapless Goldstone mode and a gapped amplitude mode. This was achieved in Ref. [45] through an effective action approach. Here, however, we concentrate on the analysis on the strong interacting limit where the ground state is a Mott insulator.
We therefore define the retarded Green's function by taking the analytic continuation from the Matsubara to real frequency domain, iω → ω + i0 + , where the disorder average of the unperturbed Green's function at zero temperature is given by with P denoting the principal value of the integral of the first term. Note that the imaginary part of (16) is controlled by the disorder distribution. The clean case is recovered by setting p(ϵ i ) = δ(ϵ i ). Hence, we observe that the disorder averaging amounts to summing up infinitely many simple poles, each one coming from a single realization of the unperturbed Green's function, thus transforming the simple-pole structure of the excitations of the clean case into a branch cut. In the strong disorder limit ∆ > U , the particle density n never sticks to an integer value implying that no Mott insulating state exists [3,7]. Thus, we focus on the case of a Mott state where ∆ is finite and smaller than the on-site interaction strength, i.e., ∆ < U . In the next section we analyze the implications of our results (15) and (16) on the spectral function.

Spectral function
We obtain the spectral function by applying the definition (1) to (15) The states of the excitations are associated with peaks of the spectral function.
Assuming that Im⟨g i (ω)⟩ is a smooth function of ω, the peaks of A(k, ω) coincide with the frequencies that satisfy Thus, the solutions of (18) for fixed n correspond to the dispersion relations of the respective excitations inside the Mott state.
If Im⟨g i (ω)⟩ approaches zero, the spectral function will only be non-zero if (18) is satisfied. When this condition is met, the spectral function can be expressed as where the plus (minus) sign relates to stable quasiparticle (quasihole) excitations. The vanishing width of the δ-distribution characterizes a stable propagation of these excitations. From the dispersion relations ω ± (k) we define the effective-mass tensor with components [48] ( On the other hand, when Im⟨g i (ω)⟩ is finite, resonances occur whenever (18) is fulfilled. Near such resonances ω ≈ ω r (k) we expand (17) according to This consists of approximating the spectral function by a Cauchy-Lorentz distribution in that region. Such an expansion yields an estimate for the renormalization factor and for the width of the distribution whose inverseτ is interpreted as the lifetime of these states [41]. As the region where Im⟨g i (ω)⟩ is finite is controlled by the disorder distribution p(ϵ i ), the emergence of these damped states is directly linked to the presence of the disorder potential. Notice that these resonances appear as poles in the lower half plane of the retarded Green's function ⟨G(k, ω + i0 + )⟩. Therefore, for each Mott lobe with fixed n, the energy band splits into two regions. The first one is associated with stable states. The second one emerges due to the disorder. It corresponds to damped states with a characteristic lifetimeτ (k). Thus, we separate the corresponding terms of the spectral function as Note that we have made no assumption so far on the specific value of n nor on the concrete shape of the bounded distribution p(ϵ i ). Therefore, all the conclusions drawn so far can be considered to be generally valid within each Mott lobe, subject to the limitations of our approximations. We remark that, even though our method is better suited for high number of dimensions, it is still applicable in D = 1 for sufficiently small J [44].
We now turn our attention to the case of a uniform distribution.

Uniform disorder distribution
To better demonstrate the above described results, we investigate further the case of uniform disorder distribution p( is the Heaviside function. Such a distribution can be generated, for instance, by customizing the intensity in speckle laser experiments [49]. In this case, (18) reduces to Note that not all solutions of (26) are solutions of (18). For n ≥ 1 we expect from the above arguments that only four solutions of (26) correspond to real excitations of each Mott lobe. Two of them are associated with stable quasiparticle and quasihole dispersions and the remaining two correspond to damped states. Thus, in order to proceed with the analysis, we must analyze (26) in each Mott lobe with fixed integer particle density n separately.

Mott lobe n = 0
For the sake of simplicity, we consider first the case of n = 0, i.e. for µ < −∆/2. This case can be interpreted as a Mott lobe where the energy spectrum only contains the quasiparticle branch, and the energy gap depends on the potential barrier, i.e. it depends both on the tunneling energy and the disorder strength. Thus, we find that in this configuration, Eq. (26) admits two solutions. The first one corresponds to stable states with diagonal components of the effective-mass tensor given according to (20) By taking the limit of ∆ → 0 in (27), we recover the clean-case dispersion The second solution corresponds to damped states with a lifetime obtained from (24) Using (25), we express the spectral function as It can straightforwardly be checked that (32) satisfies the general properties (2). Furthermore, the density of states follows from integrating (32) over the first Brillouin zone In Fig. 2 we show a plot of the resulting band structure for D = 1, D = 2 and D = 3. Although the results (27)−(33) for the Mott lobe n = 0 do not depend on the interaction energy U we use it as a measure for the energy scales in Fig. 2. This allows to localize the equilibrium points in the phase diagram defined in the J × µ plane. In Fig. 2(a), (d) and (g) the special case of the spectral function for k = 0 resembles the qualitative sketch depicted in Fig. 1. We observe a sharp peak shifted due to disorder towards lower energies for the stable states (red) along with broad distribution for the damped states (blue). In Fig. 2(b), (e) and (h) we plot (27), (29), and (30), where, for D = 2, the critical points of the Brillouin zone are defined as Γ = (0, 0), X = (0, π), M = (π, π), while for D = 3 such critical points are defined as Γ = (0, 0, 0), X = (0, π, 0), M = (π, π, 0), and R = (π, π, π). There we observe the damped states inside the lightblue-shaded region of width ∆. Note that the dispersions of both stable and damped states have jumps (dotted red lines). We can understand these jumps considering the random potential as a superposition of many spatial Fourier components. As the disorder is uncorrelated at different sites, no Fourier component with frequency larger than π/a can exist. Therefore, a matter wave excitation, which propagates, for instance, in the k x direction and interacts separately with each component of the disorder potential, would meet the condition for Bragg scattering exactly at k x = ±π/2a. Thus, we interpret these jumps as the scattering experienced by each wavevector component of the excitations interacting with the corresponding spatial frequency component of the disordered potential [50]. We point out that the group velocity described as the gradient of the dispersion relations, v g = ∇ k ω(k), vanishes at those points, which implies a localization of the wave packets. In accordance with the present results, previous studies [51] demonstrated the emergence of sub bands due to the disorder in the Bose-Hubbard model. In Fig. 2(c), (f) and (i) we plot the density of states (33). In the D = 1, Fig. 2(c), we recognize that between the Van Hove singularities corresponding to the stable excitations (red), damped states (blue) emerge, thus increasing the band width to ∆ coth(∆/4J). In D = 2 the strong peak that appears in the clean-case band center, gets divided into two peaks exactly at the frequencies where the dispersion of the stable states meets the dispersion of the damped states, namely at ω = −µ ± ∆/2. An analogous situation occurs in D = 3.
In order to understand what happens for increasing disorder, we plot the effective mass of the stable states and the lifetime of the damped states in Fig. 3. In Fig. 3(a) we find that the effect of disorder is to increase the effective mass of the stable states, thus indicating that their dispersion is becoming more flat. In Fig. 3(b), we plot the lifetime of the damped states. We deduce from such a plot that the lifetime (blue) increases for large ∆, which corresponds to the resonance becoming sharply peaked. Thus, these states become more stable when the disorder strength ∆ is of the order of the interaction energy U .
We make the connection to the quantum phase transitions by analyzing the excitations gap. For damped states the gap is given by the lower bound of the broad distribution in (32), namely E r 0 = −µ−∆/2, while the gap for stable states is obtained by expanding (27) For sufficiently small tunneling, we get that E r 0 − E + 0 ∼ ∆e −∆/2JD , which is always positive. Thus, increasing ∆ has the effect that the stable states gap closes before the gap for damped states. However, in this limit, the dispersive nature of these excitations disappears indicating a broad distribution in momentum, which is characteristic of localized states. For sufficiently strong disorder, the damped sates occupy the whole band. This corresponds to a transition from Mott to Bose glass. Therefore, we come to the fundamental conclusion that the damped states correspond to single-particle excitations of the Boseglass state. In the presence of disorder no direct Mott-superfluid transition is possible [3,52,7,14]. In order to get information on the complete quantum phase diagram, one would have to consider higher number of scattering processes [26]. Quantum and thermal fluctuations could also be included using an effective-action approach [45]. It has been proposed that one could define fluctuations in the disorder average of the mean particle density as an order parameter to identify the Bose-glass phase, in analogy to the Edwards-Anderson order parameter in the spin glass theory [53,54]. The incorporation of such an order parameter could lead to precise results inside the Bose-glass phase.
Next, we analyze the case of n = 1, where interactions become important.

Mott lobe n = 1
In the case of n = 1, assuming that the arguments of the absolute values are real, (26) simplifies to where ξ = ±1. Due to the different conditions imposed by a uniform disorder distribution on the frequencies that satisfy (34), we distinguish the intervals −µ − ∆/2 < ω < −µ + ∆/2, Real solutions of (34) inside these intervals correspond to damped states. On the other hand, frequencies that satisfy (34) outside the intervals (35) correspond to stable states. We note that (34) can be rewritten as a cubic equation of the form which admits three solutions for each value of ξ. Using the reduced form for this cubic equation and applying Cardano's formula, we get the solutions in the form where α = −1+i √ 3 2 is the primitive cubic root of unity and we have introduced the abbreviations . (41) Note that the on-site interaction energy U appears explicitly for n = 1 according to Eqs. (39)−(42) in contrast to the n = 0 case treated in Sec. 4.1. The dispersion relation of each excitation is obtained from (37) considering the different exponents of the primitive cubic root defined as l ∈ {0, 1, 2}. For the case of ξ = −1 all three solutions are real. However, for ξ = 1 only the solution for l = 1 is real. Therefore, these are the solutions of (36) which correspond to excitations of the n = 1 Mott lobe.
In order to illustrate the results, we plot in Fig. 4 the D = 1 band structure for the n = 1 Mott lobe. We observe that such a plot resembles qualitatively the n = 0 case in Fig. 2, i.e., the gap for stable states decreases and damped states emerge in the middle of the band. Additionally, as a result of the scattering with the random potential, the dispersions have jumps exactly at k x = π/2a analogously to the previous n = 0 case. However, instead of having only the quasiparticle branch, in this case the band structure shows for negative energies a quasihole branch as well.
Applying equations (20) and (24), we plot in Fig. 5 the effective mass and the lifetime corresponding to the stable and damped states in the quasiparticle and quasihole branches, respectively. We observe that, analogously to the previous case of n = 0 in Fig. 3, the effective mass of the stable states as well as the lifetime of the damped states increase in the strong disorder limit. We read off from Fig. 5 that both the effective mass of the stable states and the lifetime of the damped states in the quasihole branch (doted red line in (a) and dotted blue line in (b)) depend more sensibly upon the disorder strength.
It is important to note what happens when the gap for creating these excitations closes. In our present case, the excitation spectrum comprises both the quasiparticle and quasihole branches. During a generic phase transition, the gap of one of these branches will close while the gap for the other remains open. Restricting the analysis to the quasiparticle branch of the spectrum, the gap for stable excitations can be computed by expanding (37) near k = 0 for l = 1 and ξ = −1. In the asymptotic limit of vanishing tunneling energy, such a gap can be expressed as E + 1 ∼ U − µ − ∆(1/2 + e −∆/2JD ). The gap for damped states, however, can be computed from the left-hand side of the second line of (35), which yields E r 1 = U − µ − ∆/2. Analogously to the case analyzed in Section 4.1, the difference between the two energy gaps reads E r 1 − E + 1 ∼ ∆e −∆/2JD , which is again always positive. By applying the same reasoning as previously discussed, we observe that increasing ∆ causes the gap for stable states to close before the gap for damped states. However, in this scenario, as depicted in Fig. 5, the effective mass of stable excitations increases, and their dispersive nature vanishes, indicating a broad distribution in momentum, which is a characteristic of localized states. As the disorder strength becomes sufficiently high, the damped states fill the entire band, leading to a transition from Mott to Bose glass. This confirms our fundamental conclusion that the damped states correspond to single-particle excitations of the Bose-glass state. In the clean case, it was shown that in the first Mott lobe, when the quasiparticle gap closes during a generic phase transition, it transforms continuously into a Goldstone mode of the superfluid phase, while the quasihole gap, which remains open, transforms continuously into a gapped amplitude mode [35]. We, therefore, expect that this scenario will also hold when disorder is present, such that the stable and damped excitations will transform continuously into excitations of the superfluid phase. However, further analysis is required to draw definitive conclusions regarding the effects of disorder on these superfluid excitations. Our next area of focus is to examine the spatio-temporal propagation of stable and damped excitations.

Spatio-temporal profile of the Green's function
Thus far, we have focused on the properties of the available states for an excitation in Fourier space. We now turn our attention to the implications of the findings demonstrated on the previous section in real space and time. To this end, we investigate the full Green's function within our approximation.
First, we obtain the long-wavelength behavior of the Green's function in space by integrating (15) over the first Brillouin zone in the limit of large space separations, which yields where the length scale associated with the exponential decay of ⟨G ij (ω)⟩ turns out to be This corresponds to the mean free path of excitations between each scattering event [41]. Note that only the states at k = 0 contribute to the asymptotic behavior of ⟨G ij (ω)⟩ for large separations. The exponential decay of the single-particle Green's function is a general consequence of the mass gap present in the Mott insulating phase [47]. However, the mean free path diverges as Im⟨g i (ω)⟩ → 0, so the amplitude of propagation for stable states decays algebraically as |r i − r j | −(D−1)/2 . For finite Im⟨g i (ω)⟩, i.e., in the energy range of damped states, there is an additional exponential decay with the characteristic length ℓ defined in (43). We remark that, within our approach, the asymptotic decay of ⟨G ij (ω)⟩ given in (42) is generally valid for all Mott lobes and any form of bounded disorder distribution p(ϵ i ).
In order to get the complete picture, we write the retarded Green's function in terms of the spectral function by using the following representation [32] ⟨G Hence, the spatio-temporal profile can be written as For t = 0, the integration over the frequency domain becomes the sum rule (2). In this case, solving the integral in k yields where r (q) ij is the q-ht component of the vector r i − r j . Note that this result is the same for the clean and the disordered cases. Furthermore, using (32) for the case of n = 0 we distinguish for t > 0 the contributions corresponding to the stable and damped states, respectively. The absolute squared value of these quantities is plotted in Fig. 6. In Fig. 6(a) we observe two distinct regimes for the time evolution of |⟨G ij (t)⟩| 2 . At short time scales its amplitude is localized around r i − r j = 0. However, at time t = 16/U the amplitude starts to become extended. For long time scales only the extended states contribute to |⟨G ij (t)⟩| 2 . In Fig. 6(b) we read off that the contribution ⟨G + ij (t)⟩ of the stable states spreads through space as time increases. Therefore, in the long-time limit, there exists a finite probability of finding the excitation on a site arbitrarily distant from the site where it was created. We conclude from (42) that the algebraic decay is not enough to localize theses excitations in space. In Fig. 6(c) we notice no diffusion of the contribution ⟨G r ij (t)⟩ together with a rapid decay in time. Thus, the additional exponential decay contributes to the localization of the damped states. Therefore, the behavior of the full Green's function in space is dominated by the damped states at short time scales, while for long time scales it is dominated by the stable states.

Summary and conclusions
In conclusion, we investigated the effect of disorder on the low-energy excitations of the Bose-Hubbard model in the strongly interacting limit at zero temperature applying a perturbative field-theoretical approach to obtain a resummed expression for the spectral function. By analyzing the peaks of the spectral function we demonstrated that two different kinds of excitation states are present, namely stable states which are extended in space with algebraically decaying amplitude and damped states which are localized with exponentially decaying amplitude. By considering the limit of strong disorder, where the damped states dominate the spectrum with lifetime increased by the disorder, we argued that they correspond to low-energy single-particle excitations of the Boseglass phase. Our results inside each Mott lobe for small values of the tunneling energy are general, and therefore independent of the exact form of the bounded disorder distributions. Furthermore, by analyzing the case of uniform distribution we showed that disorder increases the effective mass of the stable states. We point out that the spectral function can experimentally be determined, for instance by Bragg spectroscopy [55,56] or by the radio frequency transfer method [57]. Moreover, it has been recently proposed that the spectral function could be used to probe key aspects of the excitation spectrum in the disordered case in quench spectroscopy experiments [58,59]. Further insight into the transition to the superfluid phase could be obtained by including loop corrections to the resummation method developed here. Future research is required to determine the role of the damped states on the multiple matter wave interference pattern which can be measured in time-of-flight experiments. It was demonstrated that the effect of spatially random fields in bosonic ultracold quantum gases can be mimicked by a time-alternating external potential resulting in a Bose-glass like nonequilibrium granular condensate [10]. Since it has been reported that a free expanding out-of-equilibrium condensate resembles a propagating optical speckle [60], exploring the connection between the effects of spatially random and time-varying potentials could be an interesting direction for future work.