Impurity dephasing in a Bose-Hubbard model

We study the dynamics of a two-level impurity embedded in a two-dimensional Bose-Hubbard model at zero temperature from an open quantum system perspective. Results for the decoherence across the whole phase diagram are presented, with a focus on the critical region close to the transition between superfluid and Mott insulator. In particular we show how the decoherence and the deviation from a Markovian behavior are sensitive to whether the transition is crossed at commensurate or incommensurate densities. The role of the spectrum of the Bose-Hubbard environment and its non-Gaussian statistics, beyond the standard independent boson model, is highlighted. Our analysis resorts on a recently developed method [Phys. Rev. Research 2, 033276 (2020)] - closely related to slave boson approaches - that enables us to capture the correlations across the whole phase diagram. This semi-analytical method provides us with a deep insight on the physics of the spin decoherence in the superfluid and Mott phases as well as close to the phase transitions.


Introduction
Understanding the dynamics of an open quantum system, i.e., a quantum system coupled to its environment, is relevant in a variety of domains including condensed matter physics, quantum computing, quantum optics and ultracold gases [1][2][3][4]. When the open system and its environment are weakly coupled, it is often a good approximation to describe the latter 1 INTRODUCTION 2 as a set of harmonic oscillators linearly coupled to the system. This class of problems is well described by the so-called Caldeira-Leggett model, when the open system is described in terms of continuous variables, or by the spin-boson model, when it is a discrete system. In any of these models, the influence of the environment on the system depends only on a single-particle spectral density, and this strongly simplifies the description of the system. The past few decades have seen the development of a large variety of methods to describe the open system dynamics in this context, including path integrals [5,6], stochastic Schrödinger equations [7,8], hierarchical systems of equations [9,10] or, when computing the full dynamics of both the system and its environment, chain mapping representations [11,12] or quantum Monte Carlo techniques [13,14].
However, when the environment is strongly correlated or non-harmonic, the above picture may no longer be accurate and more involved approaches are required to account for the resulting non-Gaussian environment statistics. The state-of-the-art methods to numerically study these systems are based on matrix-product states [15][16][17][18]; nevertheless, due to the rapid entanglement growth, these methods become highly inefficient beyond one-dimensional cases or when approaching to a critical regime.
The recent advances in locally manipulating ultra-cold gases in optical lattices has made such a platform ideal for the study of impurities coupled to a non-trivial bath [17,[19][20][21][22][23] either per se or as quantum simulators of toy models for less clean systems. In this paper, we analyze the pure dephasing dynamics of a two-level impurity whose environment is represented by a single-band Bose-Hubbard (BH) model. This problem has been recently analyzed for a one-dimensional BH environment away from its critical transition [17]. Here we take a leap forward by considering a 2D BH model and characterizing the impurity dynamics along the whole phase diagram, focusing on the critical regions. Our goal can be reached thanks to the use of a Gutzwiller technique that we recently developed [24]. The method allows us to include the relevant correlations of the bath -in particular the ones responsible for non-Gaussian effects -without being computationally demanding.
One of the main findings of our study is the strong dependence of the dephasing dynamics on the universality class of the Mott insulator-superfluid transition of the BH environment. In particular, we show that: (1) when the quantum phase transition is due to particle number change, also known as commensurate-incommensurate transition, the impurity dynamics is perfectly Markovian, being the environment dynamics dominated by single particle processes, despite the strong interactions; (2) on the other hand, when the transition occurs at fixed (integer) density, the spectrum of the bath contains multiple low-energy collective modes. Their presence leads to a non-Markovian dephasing dynamics, strongly affected by twoparticle processes in the environment, which make the standard Gaussian statistics fail. Most importantly -in close analogy with the findings of a related work on one-dimensional quantum spin baths [25] -we find that both the short and long-time behaviour of the dephasing dynamics are precise detectors of the type of universality class of the transition.
The paper is organized as follows. Section 2 is devoted to introducing the pure dephasing model, the quantum Gutzwiller approach used to access the relevant Bose-Hubbard correlations and the so-called BLP non-Markovianity measure of dephasing processes, which 2 MODEL AND THEORY 3 is taken as a reference for our analysis. In Section 3, we present our predictions for the dephasing dynamics across the phase diagram of the Bose-Hubbard environment, focusing on the intrinsic non-Markovian effects due to the lattice setting and the consequences of the spectral properties of the bath. Specifically, the role of the superfluid-Mott insulator transition is highlighted. We conclude in Section 4 including an outlook on future studies of experimental interest.

Quantum impurity in a Bose-Hubbard bath
We consider a two-level impurity coupled to a two-dimensional single-band Bose-Hubbard (BH) model [26,27] with HamiltonianĤ BH , hereafter referred to as the bath. The total Hamiltonian of the system can be written asĤ =Ĥ BH +Ĥ imp +Ĥ c witĥ where the operatorsâ r â † r annihilate (create) a boson on the lattice site r, J is the hopping energy, U the on-site bath interaction and µ the chemical potential, while r, s labels all pairs of nearest-neighboring sites. The impurity is governed by the HamiltonianĤ imp with a resonant frequency ω 0 and is coupled to the bath density at site zeron 0 via a local interaction HamiltonianĤ c with strength g.
We assume that initially the system's state is separable ρ(t = 0) = ρ 0 BH ⊗ ρ 0 imp , where ρ 0 BH is the zero-temperature ground state of the BH HamiltonianĤ BH and ρ 0 imp is the initial state of the impurity. As usual in the study of open quantum systems, we assume that the bath and the impurity are weakly coupled so that the bath's state is not too altered with respect to ρ 0 BH . Under such approximation, it is well-known that the impurity dynamics is fully characterized by the time correlation function of the environment coupling operator,n 0 . We estimate the latter by using a recently developed quantum Gutzwiller (QGW) approach [24], that has been proven to be very accurate to describe the quantum correlations of the BH model across the whole phase diagram. For completeness, we briefly review the method in the following.

The Quantum Gutzwiller method
The QGW approach combines the successful features of the Gutzwiller approximation [28] and the Bogoliubov theory of weakly-interacting gases [29] in order to develop a robust quantum many-body theory of a generic interacting lattice model. Building on the solution of the time-dependent Gutzwiller approximation [30], fluctuations on top of the mean-field ground state are quantized in terms of the elementary many-body excitations of the system 2 MODEL AND THEORY 4 and systematically included in the calculation of ground state expectation values. In spite of the local nature of the underlying Gutzwiller ansatz -see Eq. (A.1) in Appendix A -, the QGW approach accurately reproduce both local and non-local correlations across the different phases of the BH model with minimal numerical effort, showing a remarkable agreement with quantum Monte Carlo predictions concerning density correlations. Let us also mention that the QGW, when only quadratic fluctuations are considered, coincides essentially with including quantum fluctuations by slave boson approaches (see in particular [31], where the slave boson approach has been applied to the BH Hamiltonian to determine its entanglement entropy along its phase diagram).
Within the QGW, the BH environment -aside from a constant energy term -can be recast as the quadratic Hamiltonian where the operatorb α,k (b † α,k ) annihilates (creates) an excitation in the branch α with momentum k, whose energy is ω α,k . The quadratic nature of the bath Hamiltonian allows us to easily estimate its quantum correlations.
Before proceeding, we briefly review the structure of the BH excitation spectrum ω α,k along the phase diagram, since its knowledge gives important insights in the dephasing dynamics of the spin impurity, as we show in Section 3. The spectrum is well-known and can be obtained also from linear-response theory applied to the time-dependent Gutzwiller approximation [24,32]. For convenience, in Figure 1 a summary of the phase diagram and of the excitation spectra in different regimes is shown. The most relevant feature of the BH model is the existence of a quantum phase transition between a Mott insulator (MI) -which favours localized particles and occurs at integer fillings for U/J larger than a critical value -and a superfluid (SF) delocalized phase with broken U(1) symmetry. The quantum criticality is characterized by two different universality classes [26,27], depending on whether the transition point is crossed by tuning the density to a commensurate (i.e., integer) lattice filling -the so-called commensurate-incommensurate (CI) transition [at the edge of the Mott lobe: see point 2 on the blue dashed line in Figure 1 In the MI incompressible phase, the two lowest excitation branches are the gapped particle and hole excitations (not shown in Figure 1). As the SF phase is approached along a CI transition line one of the excitations becomes gapless and transforms into the superfluid gapless Goldstone mode. The low momentum dispersion relation of the Goldstone mode becomes quadratic at the transition point, while is linear in the SF phase (collisionless sound mode) [ Figure 1(b)]. Therefore, at the CI critical point the BH system, although strongly interacting, behaves as a free Bose gas of quasi-particles.
Instead, at the fixed-density O(2) critical point, both the lowest-energy modes are gapless [ Figure 1 [ Figure 1(c)]. The other gapped excitation is often referred to as the Higgs mode and it is related to the amplitude fluctuations of the order parameter [33,34].
The QGW approach provides a recipe to express operators and observables of the BH bath in terms of the excitations operatorsb α,k (see [24] and Appendix A). In particular, the impurity dynamics due to the weak coupling with the bath as described by Eq. (2) is fully characterised by the time dependent density correlation function at the impurity position. The expression for the density operatorn 0 can be written within the QGW approach aŝ where n 0 is the mean-field density and we separate the single quasi-particle contribution where V is the lattice volume. The coefficients N α,k , W αk,βp , U αk,βp and V αk,βp are given explicitly in Appendix A.
It is worth noticing that the inclusion of two-particle processes due to δ 2n into the bath description generalizes the independent boson model, where the impurity polarizationσ z couples only to linear contributions of the form (4) (see, e.g., [35]).
The two-particle contributions dominate the density correlation functions in the MI phase and close to the MI-SF transition [24]. In the following we show that this is the case also for the impurity dephasing, but not at the CI transition point.
Let us stress that, compared to other approaches like strong-interaction perturbative methods [36] and the standard Bogoliubov approximation [29], the QGW approach provides a unified description from the deep MI state to the weakly-interacting superfluids. Moreover, it yields an insight into the spectral composition of quantum expectation values.

Non-Markovianity measure of pure dephasing
Having reduced the BH environment to the effective quadratic model (2), the theoretical investigation of pure dephasing dynamics becomes tractable in the limit in which the presence of the impurity does not perturb significantly the behaviour of the environment, i.e. when the bath-impurity coupling g is small compared to all the other energy scales of the problem. For the purpose of this study, we choose to work in such weak coupling limit. Using the timeconvolutionless projection operator technique up to second order in the coupling constant g [37], the evolution of the density matrix of the impurity is proved to obey a time-local master equation [38] whereω 0 = ω 0 + g n 0 is the impurity energy splitting renormalized by the mean local density of the BH bath n 0 . As anticipated before, the dephasing rate γ(t) is completely determined by the time-dependent correlations of the bath operator coupled to the impurity -local density fluctuations in the present case - where we have defined · · · = Tr{· · · ρ 0 BH }. We recall here that the derivation of (6) does not require any assumption about the statistical properties of the environment, so that in principle the rate (7) can account also for weak-coupling effects of non-Gaussian correlations. Now, we highlight that the integrated rate is key to understanding the dephasing dynamics, as it establishes a direct connection between the decay rate γ(t) and the physical consequences of its non-Markovian features.
In the framework of the open quantum system formalism, Breuer, Laine and Piilo (BLP) have proposed a rigorous definition for non-Markovianity of a generic quantum channel [39]. Indeed, for the dephasing model studied in this work, the BLP non-Markovianity measure depends directly on the decoherence function Γ(t) via the so-called Loschmidt echo [40,41] driving the off-diagonal evolution of the impurity state ρ imp (t) ‡. In particular, the amount of non-Markovianity corresponds to the information back-flow [42][43][44] where the sum is taken over the set of time intervals [t i , t i+1 ] in which the echo increases, i.e. when γ(t) < 0. During these intervals, some of the previously lost information regarding the state of the impurity is temporarily recovered. Conversely, the Markovian character of the dynamics N + is quantified by summing L(t i+1 ) − L(t i ) over the time intervals in which quantum information is lost. It is worth underlining that, for the special open quantum system that we consider here, all non-Markovianity measures agree in distinguishing Markovian from non-Markovian evolution [45,46]. In the following sections, we will describe how a non-Markovian dephasing dynamics emerges due to strong correlations in the BH environment, focusing on the role of the universality classes of the MI-SF transition and on the importance of including non-Gaussian correlations beyond linear coupling between the bath excitations and the impurity (twoparticle contributions). In this regard, we start our analysis by illustrating how the QGW approach provides semi-analytical expressions for the dephasing rate γ(t) and the decoherence function Γ(t), with a clear distinction between single-particle and non-Gaussian correlations.

QGW expressions of γ(t) and Γ(t) and short-time behaviour of the Loschmidt echo L(t)
In this section we report for completeness the explicit expressions of the relevant quantities introduced above within the QGW formalism. Inserting the expression of the density operator (3) into the definition of the dephasing rate γ(t), we can distinguish two contributions γ(t) = γ 1 (t) + γ 2 (t). The first term is due to the linear-order part of the density operator (4), ‡ See Appendix B for a detailed definition of the BLP non-Markovianity measure and its calculation in the pure dephasing model considered in this paper.

NUMERICAL RESULTS
8 while the second contribution is generated by the two-particle density operator (5), in particular at zero temperature. Analogously, the decoherence function is given by The off-diagonal elements of the impurity density matrix will evolve according to the Loschmidt echo . From Eqs. (11)- (12) we see that the expected short-time the Gaussian behaviour e −λ t 2 [47] of the Loschmidt echo is recovered with In the following we show how both λ and the BLP non-Markovianity measure are not only extremely sensitive to the phase transition points, but behave differently depending on the universality class of the phase transition.

Numerical results
In the following we present the numerical results obtained by computing the dephasing rate functions (11)- (12) and the Loschmidt echo L(t). All the calculations have been performed on a 400 × 400 square lattice, which well approximates the thermodynamic limit and is made possible by the low numerical complexity of the QGW approach. For brevity, hereafter we will refer to the CI transition as edge transition, while the O(2) critical point will be indicated as tip transition.

Dephasing in the superfluid phase
We start our analysis about the dephasing dynamics starting from the weakly-interacting limit (deep SF phase) of the BH bath. In Figure 2(a) we report the behaviour of the dephasing rate function γ(t) [black solid line] for 2 d J/U = 1. As expected, in this regime the contribution from the single-particle gapless Goldstone mode [red dashed line] saturates the time evolution of γ(t). The dephasing rate exhibits broad oscillations around zero at short times, signalling the occurrence of non-Markovian effects, simply due to the finite bandwidth of the model. For the sake of clarity, we argue a little bit on such result, which can be better understood by expressing the dephasing γ(t) = ∞ 0 dω J(ω) sin (ω t)/ω [48] in terms of the single-particle spectral density This quantity for 2 d J/U = 1 is shown in Figure 2(b). Being the Goldstone spectrum gapless and linear at small momenta, the spectral density scales as J(ω) ∼ ω d at low frequencies §.
Nevertheless, in contrast with the non-Markovianity criterion generally adopted -obtained in [44] and fixing to d > 2 the necessary condition for memory effects in gapless baths -, we observe that γ(t) has negative values in our d = 2 model. The reason is that usually an environment with infinite-bandwidth modes is considered in the literature [44], resulting in a smooth cutoff of the spectral density. The finite bandwidth of the BH model excitations implies a sharp frequency cutoff of J(ω) corresponding to the Goldstone mode energy at the edge of the Brillouin zone, ω G,π . Correspondingly, we observe that the oscillations of γ(t) occur on a time scale τ G = 2 π/ω G,π [vertical dotted line in Figure 2(a)] set by the bandwidth of the Goldstone excitation .

Dephasing dynamics at the MI-SF transition
Moving away from the deep SF phase and approaching the MI-SF critical region, the fate of the SF non-Markovian dynamics turns out to strongly depend on the type of crossed critical point. In particular, crossing the edge transition [blue dashed line in Figure 1] the amplitude of memory effects decreases with increasing interaction U/J until the dynamics becomes purely Markovian on the Mott boundary. On the contrary, crossing the tip transition [red dashed line in Figure 1(a)], the non-Markovianity is even more enhanced by quantum fluctuations with respect to the deep SF phase.
In panel (c) of Figure 2 we display the evolution of γ(t) for different values of 2 d J/U upon approaching the edge transition. We observe that, close to the critical point (2 d J/U ) edge c = 0.08, γ(t) becomes strictly positive, the dynamics slows down significantly -as compared with the evolution for the deep SF bath in panel (a) -and the dephasing rate reaches a constant value γ(t) ∼ η at asymptotically large times. Hence, a transition from a non-Markovian to a Markovian regime occurs and, at the transition point, the Loschmidt echo acquires the typical exponential behaviour L(t) ∼ e −η t of a Lindbladian evolution. The origin of the Markovian behaviour is due to the peculiar spectral properties of the BH model on the edge of the Mott lobe. In particular, as illustrated in Subsection 2.2: (i) the Goldstone mode turns into an effective quasiparticle branch with quadratic energy dispersion; (ii) the Higgs mode keeps a finite energy gap. It follows that the strongly-correlated superfluid sitting close to the edge critical point can be described as a dilute free-boson gas with an effective mass renormalized by the vicinity of the Mott phase [24,27]. Indeed, it is easy to check that for a free Bose gas -with or without a lattice -the Loschmidt echo decays always exponentially as L(t) ∼ e −η t for d = 2 ¶. As in the deep SF case, the Goldstone single-particle contribution to γ(t) is the dominant one, but, in this case, the two-particle contributions to γ(t) are non-negligible in the edge critical region. However, we find that such a contribution integrates to zero See Appendix D for an extensive discussion on the difference between lattice and continuous models at the level of the spectral density J(ω) and the dephasing function γ(t). ¶ We refer again the reader to Appendix D for the explicit expressions of γ(t) and L(t) of a free boson gas on a lattice and on the continuum. See also Figure D1 for exact numerical results on the behaviour of γ(t) for lattice free bosons in one and two dimensions.

11
identically in the time integral of the decoherence function Γ(t) = t 0 dτ γ(τ ). In this respect, the irrelevance of non-Gaussian bath correlations can be seen as a natural consequence of the effective single-particle description when crossing the CI critical region.
The result is very different when approaching the commensurate transition at the tip of the Mott lobe, as shown in panel (d). The dynamics appears to be always non-Markovian and the memory effects are amplified with respect to the deep SF regime. The dephasing rate γ(t) gets a relevant contribution from the Higgs excitation and, most importantly, from the two-particle couplings [black dot-dashed line]. Specifically, the competition between the Goldstone and Higgs branches is evidently due to the closing of the Higgs gap at the tip critical point. For the same reason, one gets a sizable contribution to the dynamics from twoparticle correlations due to the coupling between the Goldstone and Higgs modes encoded in the structure factors W G,k;H,p and W H,p;G,k in the two-mode part of the density operator (5). Decreasing further 2 d J/U towards the critical point, non-Gaussian correlations eventually become the dominant contribution to γ(t), since the order of magnitude of the single-particle weights N α,k is totally suppressed on the brink of the MI-SF transition [24].
In this respect, we want to stress that two-particle processes become the only nonvanishing contributions to density correlations when the BH environment enters the MI phase [24]. Therefore, the dephasing dynamics undergoes a substantial change across the edge transition, where the single-particle picture is abruptly replaced by non-Gaussian correlations, while at the tip transition the single-to-two particle transfer of spectral weight appears to be a smoother crossover.

Short-time dephasing process and non-Markovianity measure
A concise way to visualize the previous results is provided by inspecting the dephasing dynamics from the point of view of the Loschmidt echo. Specifically, we focus our analysis on two complementary features of the decoherence process, namely (i) the short-time behaviour of the impurity decoherence L(t → 0) = exp (−λ t 2 ) and (ii) the estimation of the information back-flow N − . More precisely, we renormalize the information back-flow by the overall coherence loss as R = N − /N + , which provides a more effective measure of non-Markovianity while changing the bath parameters [49].
Our numerical results for the short-time decoherence rate λ, given by the expression (13), are reported in panel (a) of Figure 3. Reaching the MI-SF critical region from the deep SF phase, the decoherence rate λ decreases as a consequence of the stronger non-Markovianity driven by interactions in the BH bath. Reducing further the hopping energy, we observe that λ presents different behaviours depending on the type of approached transition. At the CI critical points, the decoherence rate quickly drops to a small value (decreasing by almost two orders of magnitude) entering the MI phase, where we find that λ ∝ (J/U ) 2 . The first derivative of λ with respect to J/U presents a discontinuity at the critical point. Conversely, when crossing the transition at the lobe's tip, λ is a smooth function of the hopping energy. We notice that our latter result nicely resembles what has been found for the impurity decoherence process in a d = 1 interacting quantum spin bath [25], which has a critical point of the same O(2) universality class. Therefore, as for the static properties [24], our method is able to capture the strong correlation also in this time dependent scenario, importantly beyond the one-dimensional case and without strong numerical requirements. The time-integrated dephasing dynamics, in the form of the non-Markovianity measure R, is even more affected by the type of critical correlations than the short-time decoherence. Our numerical results for R across the edge and tip transitions are reported in panel (b) of Figure 3 with the same color code of panel (a).
In the deep SF limit J/U 1, we find that both the information flows N ± tend to zero scaling as (J/U ) −1 , such that their ratio R is a constant. This indicates that, when embedded in a weakly-interacting gas, the impurity dephases according to a fixed fraction of information loss. When approaching the strongly-interacting regime, the renormalized backflow R reaches a maximal value well before the MI-SF transition. This suggests that, away from critical region, the primary effect of stronger interactions is to increase the amount of information recovered by the impurity during the dynamics. When approaching the critical point the non-Markovianity measure R starts decreasing and its behaviour depends on how the MI-SF is crossed.
Crossing the CI transition R rapidly vanishes, being zero within a small window in the SF region. This result perfectly mirrors the non-Markovian to Markovian transition displayed in Figure 2(c) and the effective free-particle description of the SF at the CI critical points. The quantity R show a discontinuous behaviour, when entering the insulating phase. This result finds a straightforward interpretation in terms of the particle-hole excitations of the Mott phase [17]. Due to their incoherent character, these modes excite doublon-holon 4 SUMMARY AND OUTLOOK 13 pairs with a finite correlation length, so that density fluctuations are localized in real space. Therefore, when particle-hole excitations couple to the impurity, the information flowing to the BH environment remains localized in a small neighbourhood of the impurity and is likely to be restored after a short time due to another particle-hole excitation. As the amplitude of density fluctuations in the Mott phase increase with 2 d J/U , the absolute value of both the information flows N ± increases accordingly; on the other hand, the renormalized back-flow R decrease as a consequence of the increasing BH correlation length, which prevent part of the lost information from flowing back to the impurity. However, since at the edge transition either the particle or the hole branch remains gapped, a finite correlation length still controls the dynamics exactly at the critical point [24,27], before diverging in the SF phase. This discontinuous behaviour of the correlation length is at the roots the finite jump in R across the non-Markovian to Markovian transition.
The behaviour is different at the tip transition. As shown before, in this regime critical fluctuations are mainly due to non-Gaussian correlations, whose main effect is to amplify the oscillation amplitude of the dephasing rate γ(t). Therefore, the amount of total information flowing both from and to the impurity grows accordingly. Nevertheless, the renormalized backflow R still converges to zero at the critical point, meaning that eventually the BH environment becomes effectively Markovian at the critical point. It follows that, in contrast with the edge case, R is found to be a continuous function of the hopping 2 d J/U across the tip transition, but with a very sharp non-monotonic profile [red line in Figure 3(b)].

Summary and outlook
In this paper, we present an exhaustive account of the non-Markovian effects characterizing the dephasing dynamics of an impurity embedded in a Bose-Hubbard environment undergoing the superfluid-Mott transition.
Our analysis addresses the impurity problem beyond the standard formalism of open quantum systems. The two main new features are the inclusion of the effects of the strong correlations and phase transitions in the environment and the extension beyond the onedimensional case in a flexible and numerically cheap way. Thereby, our method is, to the best of our knowledge, the first one that allows an efficient description of an open system that is coupled to an environment undergoing a critical transition.
Strong signatures of deviation from a Markovian behavior due to the spatial discreteness of the lattice setup, not explicitly discussed in previous works, have also been highlighted in the interacting superfluid phase and related to key features of the spectral density J(ω). This suggests the idea that the very same phenomenon could take place in different lattice models whose dynamics is governed by common spectral properties. Furthermore, we observed that the amount of non-Markovianity of the dephasing process is particularly large when approaching the O(2) critical region, where two-particle effects become more relevant in the physical picture and thus the environment differs more significantly from the standard spinboson description. This opens the path for further investigations into the role of strong non-Gaussian, i.e. two-particle, correlations in the presence of strong memory effects.

SUMMARY AND OUTLOOK 14
More importantly, we have found that, when the BH environment approaches the superfluid-Mott criticality, the dephasing dynamics is extremely sensitive to the universality class of the superfluid-Mott transition. In this regard, we have shown that not only the deviation from Markovianity, but also the short-time behaviour of the dephasing dynamics carries strong signatures of the type of criticality approached by the environment. This remarkable result agrees with similar findings for interacting quantum spin baths [25] with a complementary approach, suggesting a generality which goes beyond the precise nature and the dimensionality of the bath.
Finally, from an experimental perspective, the sharp difference between the dephasing processes at the different superfluid-Mott transitions discussed in this work identifies the study of the decoherence dynamics and, in particular, non-Markovianity measures of impurity dephasing as an unambiguous probe of the type of critical behaviour experienced by the environment.

Appendix A. Quantum Gutzwiller approach in a nutshell
Following the main derivation steps of [24], in this Appendix we briefly review the quantum Gutzwiller technique that we employ for a systematic evaluation of quantum correlations in the BH environment.
Our starting point is the Gutzwiller ansatz where the wave function is site-factorized and the complex amplitudes c n (r) of each local Fock state |n, r are variational parameters with normalization condition n |c n (r)| 2 = 1. In our specific case, we draw on the simple form of (A.1) to reformulate the BH model in terms of the following Lagrangian functional In the previous equation, the dot indicates the temporal derivative, are the matrix elements of the on-site terms of the BH HamiltonianĤ BH in Fock space and is the mean-field order parameter. In this formulation, the conjugate momenta of the parameters c n (r) are c * n (r) = ∂L/∂ċ n (r). The classical Euler-Lagrange equations associated to Lagrangian (A.2) are the so-called time-dependent Gutzwiller equations as derived, e.g., in [32,50]. In a uniform system, the stationary solutions are homogeneous: in particular, the system is in a Mott Insulator (MI) state if ψ(r) = 0 and in a superfluid (SF) state otherwise.
In order to go beyond the Gutzwiller approximation introduced above, it is natural to consider how quantum effects populate the excitation modes of the system and to investigate how they affect the observable quantities. We include quantum fluctuations by building a theory of the excitations starting from Lagrangian (A.2) via canonical quantization [51,52], namely promoting the coordinates of the theory and their conjugate momenta to operators and imposing equal-time canonical commutation relations ĉ n (r),ĉ † m (s) = δ r,s δ n,m . (A.5) In analogy with the Bogoliubov approximation for dilute Bose-Einstein condensates [53,54], we expand the operatorsĉ n around their ground state values c 0 n , obtained by minimizing the energy Ψ G Ĥ BH Ψ G , asĉ n (r) =Â(r) c 0 n + δĉ n (r) .
The normalization operatorÂ(r) is a function of δĉ n (r) and δĉ † n (r) and ensures the proper normalization nĉ † n (r)ĉ n (r) =1. By restricting to local fluctuations orthogonal to the ground state n δĉ † n (r) c 0 n = 0 one haŝ In a homogeneous system, it is convenient to work in momentum space by writing where V is the lattice volume. Inserting Eq. (A.8) in Ψ G |Ĥ BH Ψ G and keeping only terms up to the quadratic order in the fluctuations, we obtain where E 0 is the mean-field ground state energy, the vector δĈ(k) contains the components δĈ n (k), andL k is a pseudo-Hermitian matrix, for the explicit expression of which we refer the interested reader to [24]. A suitable Bogoliubov rotation of the Gutzwiller operators in terms of the fundamental excitation modes of the system where eachb α,k corresponds to a different many-body excitation mode with frequency ω α,k , labeled by its momentum k and branch index α. Bosonic commutation relations between the annihilation and creation operatorsb α,k andb † α,k , are enforced by choosing the usual Bogoliubov normalization condition where the vectors u α,k (v α,k ) contain the components u α,k,n (v α,k,n ). The effective, quadratic description of the BH environment in terms of its collective modes (A.11) provided the QGW not only allows for a direct reinterpretation of the pure dephasing model (1), but also opens a simple route to the calculation of any expectation value of the bath operators. Based on the quantization procedure outlined before, the evaluation of average value of any observable Ô â † r ,â r consists in applying a four-step procedure that we summarize as follows: (iv) Taking advantage of the quadratic character of the QGW Hamiltonian, invoke Wick theorem to compute the expectation value of products of operators on Gaussian states -such as ground or thermal states obtained from H QGW . The very same protocol determines the expansion of the BH local density operator (3) in terms of single (4) and two-particle (5) operator-valued expressions of the collective modeŝ b α,k , from which the bath correlation functions are systematically extracted. For the sake of completeness, we report the exact expressions of the one and two-particle structure factors of the density channel, whose derivation is extensively discussed in [24].

Appendix B. Pure dephasing and BLP non-Markovianity measure
The definition of the Breuer-Laine-Piilo (BLP) measure [39] derives from considering non-Markovian those systems in which a back-flow of information from the environment to the open system occurs during the dynamics. This information recovery is formally identified by an increase in the distinguishability of pairs of evolving quantum states of the system. In detail, a system is non-Markovian if there is a pair of system initial states ρ S (0) and ρ (2) S (0), such that for certain times t > 0 their distinguishability grows, namely σ ρ (1) S ; t is called the information flux at time t and D ρ (1) S (t)   Tr ρ (1) where λ i are the eigenvalues of the matrix ρ (1) S . The physical interpretation of the trace distance (B.16) is that it is related to the maximum probability of distinguishing between two quantum states. In an open quantum system, this probability in general tends to decrease in time, as the system information is lost to the environment, except when the dynamics is non-Markovian. In this case, the system regains part of the previously lost information. According to the BLP criterion, the amount of non-Markovianity of a quantum process Λ can be quantified through the measure which reflects the maximum amount of information that can flow back to the system for a given process. As proven in [55], for all finite-dimensional quantum systems the evaluation of (B.18) can be optimized by considering initial states ρ S (0) that are orthogonal and lie on the boundary of the subset of physical states. In the case of the two-level impurity undergoing pure dephasing studied in this paper, the open system dynamics is driven by the master equation (6), which allows for a simple rewriting in the vector representation of the density matrix, where we have defined ρ ij = Tr S {ρ S (t) |i j|}, with |i = |1 , |2 standing for the two possible states of the impurity, and neglected the unitary evolution terms set by the renormalized transition frequencyω 0 . The analytical integration of (B.19) yields where φ t is the dynamical map of the system density matrix associated to the pure dephasing dynamics. The function coincides with the so-called Loschmidt echo [41], defined as L(t) = | ψ(t)|ψ 0 (t) | 2 , where |ψ 0 (t) is the bath ground state evolved according to its own Hamiltonian, while |ψ(t) is the time-evolved bath state in presence of the open system. Indeed, the off-diagonal matrix elements of the system density matrix ρ S are given by L(t) exactly. Choosing two initial states that are orthogonal and lie on the Bloch sphere of the two-level system ρ (1) we find that the the trace distance (B.16) reads Therefore, we obtain that the distinguishability rate is given by and σ ρ S ; t > 0 for some t when the dephasing rate γ(t) is negative, leading to non-Markovian dynamics. Finally, it is straightforward to deduce that the non-Markovianity measure (B.18) is provided by the values of the Loschmidt echo L(t) at the boundaries of those time intervals [t i , t i+1 ] over which γ(t) < 0, namely     which takes into account time periods for which γ(t) > 0.

Appendix C. Dephasing dynamics at incommensurate filling
In this Appendix, we report and discuss the quantitative evolution of the dephasing rate γ(t) and of the Loschmidt echo L(t) as the BH bath becomes strongly-interacting without entering the Mott phase and, on the contrary, retaining a superfluid character. Specifically, this corresponds to reach the hard-core boson limit of the BH model by increasing the boson interaction U at fixed non-commensurate density. Typical constant-density contours in the strongly-interacting SF phase are shown in Figure C1. Figure C2(a) shows the change in the dephasing rate γ(t) for decreasing hopping energy 2 d J/U at fixed density n = 0.6 (see the corresponding solid black line in Figure C1). We observe that, upon approaching the hard-core limit 2 d J/U → 0 from the deep SF phase, the order of magnitude of γ(t) increases significantly, while the time scale of the dephasing dynamics slows down, in such a way that the profiles of γ(t) at different values of 2 d J/U are related by a simple scaling relation. On the other hand, the strongly-correlated SF regime still exhibits an evident non-Markovian character, as recognizable also in the oscillating behaviour of the Loschmidt echo L(t), see Figure C2(b). Here, we can appreciate how non-Markovianity and the overall magnitude of γ(t) compete in controlling the amount of dephasing of the impurity. However, at very small 2 d J/U , the strong enhancement of the amplitude of γ(t) wins over revival effects and induces almost complete dephasing in a small time interval. These results find an intuitive explanation in the physical properties of the hard-core SF state. For t 1/J, strong bath correlations prevent the density excitations induced by the presence of the impurity from leaving a neighbourhood of the impurity itself, therefore leading to the strong-positive density correlations observed in Figure C2(a). However, being the hard-core phase still coherent in character, hopping process are favoured at larger times and invert the sign of γ(t) in analogy with what we observe in the deep SF regime. Therefore, the total amount of dephasing depends on whether local density correlations are sufficiently strong to overcome non-Markovian effects due to long-range coherence.
The dependence of the dephasing rate on the lattice filling can be understood by looking at Figure C2(c)-(d), referring to a larger filling n = 0.8. In particular, we notice that the oscillation amplitude of γ(t) and the speed of the dephasing process decreases as the bath density is increased towards the integer value n = 1 required for crossing the MI-SF transition.
Finally, we report the remarkable fact that, upon reaching the hard-core SF regime, the Goldstone mode alone still provides the most important part of γ(t), which is essentially given by its Gaussian contribution γ 1 (t) (see the discussion of Section 2. that a single-particle description of the BH bath is a good approximation for the dephasing dynamics when the impurity is embedded in a strongly-interacting superfluid away from the MI-SF criticality.