Pump-probe Auger-electron spectroscopy of Mott insulators

TL;DR

This paper explores how pump-probe Auger-electron spectroscopy reveals the dynamics of two-hole bound states in Mott insulators, highlighting the method's sensitivity to electronic correlations and screening effects.

Contribution

It provides a theoretical analysis of nonequilibrium Auger spectra in Mott insulators, emphasizing the role of dynamical screening and orbital expansion effects.

Findings

Spectra of photo-doped systems resemble higher filling equilibrium spectra.

Bound two-hole states are highly stable over femtosecond timescales.

Oscillations in satellite barycenter reflect dynamical Coulomb screening.

Abstract

In high-resolution core-valence-valence (CVV) Auger electron spectroscopy from the surface of a solid at thermal equilibrium, the main correlation satellite, visible in the case of strong valence-electron correlations, corresponds to a bound state of the two holes in the final state of the CVV Auger process. We discuss the physical significance of this satellite in nonequilibrium pump-probe Auger spectroscopy by numerical analysis of a single-band Hubbard-type model system including core states and a continuum of high-energy scattering states. It turns out that the spectrum of the photo-doped system, due to the increased double occupancy, shares features with the equilibrium spectrum at higher fillings. The pumping of doublons can be watched when working with overlapping pulses at short $\Delta t$. For larger pump-probe delays $\Delta t$ and on the typical femtosecond time scale for electronic relaxation processes, spectra are hardly $\Delta t$-dependent, reflecting the high stability of bound two-hole states for strong Hubbard-$U$. We argue that taking into account the spatial expansion of single-particle orbitals when these are doubly occupied, as described by the dynamical Hubbard model, produces an oscillation of the barycenter of the satellite as a function of $\Delta t$. Pump-probe Auger-electron spectroscopy is thus highly sensitive to dynamical screening of the Coulomb interaction.

Equations

$$H_{\mathrm{Hub}}(t)=H_{\mathrm{hop}}(t)+U\sum _i{n}_{i\uparrow }{n}_{i\downarrow },$$

$$H_{\text{hop}}(t)=-T\sum _{i\sigma }\left(e^{i{A}_{\text{pmp}}\left(t\right)}{c}_{i\sigma }^{†}{c}_{i+1\sigma }+\text{h.c.}\right).$$

$$H_{\mathrm{core}}=-E_c\sum _{i\sigma }{n}_{i\sigma }^{(\mathrm{f})}+U_{cc}\sum _i{n}_{i\uparrow }^{(\mathrm{f})}{n}_{i\downarrow }^{(\mathrm{f})}+U_{cv}\sum _{i\sigma {\sigma }^{\prime }}{n}_{i\sigma }^{(\mathrm{f})}{n}_{i{\sigma }^{\prime }},$$

$$\mathcal{E}_{\gamma}(t)=\mathcal{E}_{\gamma}^{(0)}\cos\big{[}\Omega_{\gamma}\left(t-t_{\gamma}\right)\big{]}s_{\gamma}(t)\>,$$

$$s_{\gamma }(t)=\exp \left[-\frac{1}{2}{\left(\frac{t-t_{\gamma }}{{\tau }_{\gamma }}\right)}^2\right],$$

$$A_{\gamma }(t)=-{\int }_{-\infty }^{t}d{t}^{\prime }{\mathcal{E}}_{\gamma }\left(t^{\prime }\right).$$

$$H_{\mathrm{scatt}}=\int dϵ\rho \left(ϵ\right)ϵa_{\sigma }^{†}\left(ϵ\right)a_{\sigma }\left(ϵ\right).$$

$$V_{\text{prb}}\left(t\right)={\mathcal{E}}_{\text{prb}}\left(t\right)\sum _{i\sigma }\int dϵ\rho \left(ϵ\right)d_i\left(ϵ\right)a_{\sigma }^{†}\left(ϵ\right)f_{i\sigma },$$

$$H_A=\sum _{i\sigma }\int dϵ\rho \left(ϵ\right)\left(U_i^A\left(ϵ\right)a_{\sigma }^{†}\left(ϵ\right)f_{i-\sigma }^{†}{c}_{i-\sigma }{c}_{i\sigma }+h.c.\right),$$

$$H_{\text{el}}(t)=H_{\text{Hub}}\left(t\right)+H_{\mathrm{core}}+H_{\text{scatt}}+H_A+V_{\text{prb}}(t).$$

$$I_{\text{AES}}\left(\epsilon\right)=\>2\pi\rho_{\phi}\lim_{t_{\text{max}}\to\infty}\int\limits_{0}^{t_{\text{max}}}dt^{\prime}\>s^{2}_{\text{prb}}(t^{\prime})\cos^{2}[\Omega_{\rm prb}(t^{\prime}-t_{\rm prb})]\;\big{<}\Psi_{\phi}(t_{\text{max}},t^{\prime})\vphantom{n(\epsilon)\Psi_{\phi}(t_{\text{max}},t^{\prime})}\big{|}n(\epsilon)\big{|}\Psi_{\phi}(t_{\text{max}},t^{\prime})\vphantom{\Psi_{\phi}(t_{\text{max}},t^{\prime})n(\epsilon)}\big{>}\;,$$

$$\big{|}\Psi_{\phi}(t_{\text{max}},t^{\prime})\big{>}=U(t_{\text{max}},t^{\prime})V^{\left(1\right)}_{\text{prb}}\,U(t^{\prime},0)\big{|}\Psi(0)\big{>}\>,$$

$$V_{\text{prb}}^{(1)}\equiv d_0\sum _{i\sigma }{f}_{i\sigma }$$

$$\cos^{2}\big{[}\Omega_{\text{prb}}\left(t-t_{\text{prb}}\right)\big{]}\to\overline{\cos^{2}}\big{[}\Omega_{\text{prb}}\left(t-t_{\text{prb}}\right)\big{]}=\frac{1}{2}\;,$$

$$\int d{t}^{\prime }{s}_{\text{prb}}^2(t^{\prime })f(t^{\prime })\approx \sum _{i=1}^{N_t}{w}_{i}f(t_i).$$

$$V_{\text{prb}}\left(t\right)\to {\mathcal{E}}_{\text{prb}}\left(t\right)d_{i_0}\left(ϵ\right)\int \rho \left(ϵ\right)dϵa_{\uparrow }^{†}\left(ϵ\right)f_{i_0\uparrow }$$

$$H_A\to U_A\int dϵ\rho \left(ϵ\right)\left(a_{\downarrow }^{†}\left(ϵ\right)f_{i_0,\uparrow }^{†}{c}_{i_0\uparrow }{c}_{i_0\downarrow }+h.c.\right)$$

$$\sigma\left(t_{\text{max}}\right)=2\pi\rho_{\phi}A_{s}\big{|}d_{0}\big{|}^{2}\big{<}n_{i\uparrow}n_{i\downarrow}\big{>}\left(1-e^{-2\Gamma t_{\text{max}}}\right),$$

$$\sigma\left(t_{\text{max}}\right)\approx 2\pi\rho_{\phi}A_{s}\big{|}d_{0}\big{|}^{2}\big{<}n_{i\uparrow}n_{i\downarrow}\big{>}\cdot 2\Gamma t_{\text{max}}\>,$$

$$2\Gamma=2\pi\rho_{A}\big{|}U_{A}\big{|}^{2}\;.$$

$$A_s=\frac{1}{2}\int d{t}^{\prime }{s}^2\left(t^{\prime }\right)=\frac{\sqrt{\pi }}{2}{\tau }_{\text{prb}}.$$

$${\tau }_c=\frac{1}{2\Gamma }.$$

$$\sigma \left(t_{\text{max}}\right)\sim \alpha \left(1-e^{-2\Gamma \left(t_{\text{max}}-t_{\text{prb}}\right)}\right).$$

$$\sigma \left(t_{\text{max}}\right)\sim 2\Gamma \left(t_{\text{max}}-t_{\text{prb}}\right).$$

$$d_{\text{tot}}\left(t\right)=\frac{1}{L}\sum_{i}\big{<}\Psi\left(t\right)\vphantom{n_{i\uparrow}n_{i\downarrow}\Psi\left(t\right)}\big{|}n_{i\uparrow}n_{i\downarrow}\big{|}\Psi\left(t\right)\vphantom{\Psi\left(t\right)n_{i\uparrow}n_{i\downarrow}}\big{>}$$

$$H(t)=H_{\text{el}}(t)+H_{\text{bos}}+H_{\text{int}}.$$

$$H_{\text{bos}}=\frac{{\omega }_0^2}{2}\sum _i{p}_i^2+\frac{1}{2}\sum _i{q}_i^2$$

$$H_{\text{int}}=k\sum _i{q}_i{O}_i$$

$$U_{i,\text{eff}}(t)=U+k{q}_i(t),$$

$$\overline{I}_{\text{AES}}=\int d\epsilon~{}\epsilon~{}I_{\text{AES}}\left(\epsilon\right)\Big{/}\int d\epsilon~{}I_{\text{AES}}\left(\epsilon\right)\>.$$

Full text

I Introduction

High-resolution core-valence-valence (CVV) Auger-electron spectroscopy (AES) Auger (1925); Lander (1953); Powell (1988) can be employed as a tool to examine the electronic structure of solid surfaces. It complements photoemission spectroscopy Feuerbacher et al. (1978); Cardona and Ley (1978); Courths and Hüfner (1984); Damascelli et al. (2003) as the standard technique to access the occupied states close to the Fermi energy. Compared to photoemission, AES has the advantage of being directly sensitive to local electron correlations: Verdozzi et al. (2001) Due to the core state involved in the transition, the CVV Auger process is predominantly local, i.e., the two final-state holes are mainly created at the same lattice site. For a “weakly correlated” simple metal, the Coulomb interaction among the valence electrons can be accounted for on a mean-field level. In this case the final-state holes quickly delocalize and propagate essentially independently through the lattice. Disregarding matrix-element effects, Potthoff (2001); Potthoff et al. (2001) the raw Auger spectrum is given by the occupied part of the self-convolution of the local density of the states. Lander (1953) The spectrum from a “strongly correlated” metal, on the other hand, may look substantially different. A sufficiently strong Hubbard-type local interaction $U$ causes a two-hole bound state which moves as a compound object through the lattice Rausch and Potthoff (2017) and is visible is the Auger spectrum as a correlation satellite. Rausch and Potthoff (2016) The important difference with respect to photoemission spectroscopy is that, for $U$ much larger than the valence bandwidth $W$, the satellite nearly takes the whole spectral weight.

A nice example is Ni, where the two-hole bound state, the famous “6-eV satellite”, is not easily seen in photoemission Guillot et al. (1977) while the Ni CVV Auger spectrum exhibits a strong correlation satellite. Mårtensson and Johansson (1980); Nolting (1990); Wegner et al. (2000) The Cini-Sawatzky theory Cini (1977); Sawatzky (1977) developed in the late 1970s for systems with a completely filled valence band predicts the position of the satellite as compared to the center of gravity of the two-hole scattering continuum to be approximately given by $U$. This has been exploited in the past in an attempt to measure $U$ for transition metals and compounds. de Boer et al. (1984); Brühwiler et al. (1993)

For the general case of a partially filled band, the theoretical description is challenging, even for simple models, such as the single-band Hubbard model. Gutzwiller (1963); Hubbard (1963); Kanamori (1963) Here, the Auger spectrum must be obtained from the local elements of the two-particle Green’s function, and this task requires to deal with the typical obstacles of a full many-body problem in the strong-coupling case. Various ideas have been used to generalize the Cini-Sawatzky theory, Cini (1979); Treglia et al. (1981); Drchal and Kudrnovsky (1984); Oleś et al. (1986); Cini and Verdozzi (1986); Nolting (1990) including perturbative, exact-diagonalization, decoupling and diagrammatic re-summation techniques. A numerically exact solution for the entire filling range of the one-dimensional Hubbard model could be provided recently. Rausch and Potthoff (2016)

In recent years, the idea of repulsively bound pairs of particles has also attracted the attention of researchers in the field of ultra-cold quantum-gases. Winkler et al. (2006) As a consequence of energy and momentum conservation, there is no phase space available for a dissociation of a pair of bosons in a Bose-Hubbard model. Bloch (2005); Jaksch and Zoller (2005) In case of a finite concentration of particles (bosons or fermions), a repulsively bound pair can only decay in a high-order scattering process involving several low-lying excitations of the continuum, and the lifetime is exponentially long in $U$, as verified experimentally SGJ+10 and theoretically. SGJ+10 ; Lenarčič and Prelovšek (2013); Chudnovskiy et al. (2012) On a short time scale, there is an initial “decay”, Hofmann and Potthoff (2012); Rausch and Potthoff (2017) which is related to the energy-time uncertainty principle.

It is tempting to ask if the decay of a repulsively bound pair of final-state holes can be watched on the time axis by means of time-resolved CVV Auger spectroscopy. While time-resolved photoemission Iwai et al. (2003); Perfetti et al. (2006); Wall et al. (2011); Braun et al. (2001) has become popular in the last decade to study the time evolution of electronic excitations of materials characterized by strong electron correlations well beyond the linear-response regime, there has not much attention been paid to Auger-electron spectroscopy – despite its high sensitivity to correlation effects. Atomic systems represent a notable exception. Drescher et al. (2002)

This is the starting point and motivation of the present paper. Here, we like to ask what be learnt from time-dependent AES from solid surfaces. To this end we propose an exemplary theoretical description of time-dependent AES covering essential features of strongly-correlated systems and focus on the time-dependence of the Auger correlation satellite in particular. To be specific, the single-band Hubbard model at half-filling and, as a variant, the dynamic Hubbard model Hirsch (2001, 2002); Werner and Eckstein (2016) will be considered. At strong Hubbard interaction $U$, these represent generic models for a correlation-induced Mott insulator. The present study profits from our previous work Rausch and Potthoff (2016) on the equilibrium two-hole spectrum of the Hubbard model at and off half-filling, where a comprehensive analysis of the various spectral features was provided.

The theory for the nonequilibrium case addresses a realistic experimental pump-probe setup, see Fig. 1. We consider an optical pump pulse with a mean frequency of the order of a few eV to bridge the Mott gap and with a duration in the femtosecond regime. The pump drives the system out of equilibrium and generates an initial nonthermal conduction-electron state with an enhanced fraction of doubly occupied sites. This state evolves in time. A relaxation of the electronic degrees of freedom generically takes place on the fs scale. After a time delay $\Delta t$, a hole in an inner core shell is created by an X-ray probe pulse which opens a channel for the subsequent CVV Auger decay. The dependence of the Auger spectrum on the time delay $\Delta t$ reflects the time evolution of the conduction-electron state. Our main focus is on the intense Auger correlation satellite and particularly on its time-dependent position and weight. Which kind of information on the electronic structure is encoded here? This fundamental question can be answered within simplified Hubbard-type models for the correlated electronic structure of a Mott insulator augmented by core states and a continuum of high-energy scattering states using exact-diagonalization techniques.

In the following section II, we introduce the model and the typical parameters involved in the calculations. Section III outlines our theoretical approach. We first present results on the core-hole decay in Sec. IV before discussing the pump-probe Auger spectra for the conventional Hubbard model in Sec. V and the dynamic Hubbard model in Sec. VI. Our findings are summarized in Sec. VII.

II Model

To be definite and to facilitate the numerical evaluation, the theory will be formulated for the one-dimensional single-band Hubbard model at half-filling as a prototypical model for the correlated electronic structure of a Mott insulator. Gebhard (1997) An extension to higher dimensions and multi-orbital models is straightforward in principle. The Hamiltonian is given by

[TABLE]

where $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$ is the occupation-number operator, and $c^{\dagger}_{i\sigma}$ creates an electron at lattice site $i=1,...,L$ with spin projection $\sigma=\uparrow,\downarrow$. The strength of the local Hubbard interaction is given by $U$. The tight-binding hopping term of the Hamiltonian $H_{\rm hop}(t)$ with the nearest-neighbor hopping-matrix element $T$ reads

[TABLE]

The energy scale (and with $\hbar\equiv 1$ the time scale as well) is fixed by setting $T\equiv 1$. The hopping Hamiltonian is time dependent as it includes the coupling of the light field to the electronic system, i.e., the pump pulse. This coupling is described in a standard manner in the velocity gauge by using the Peierls substitution. Thereby, the hopping matrix element between the valence orbitals acquires a complex phase which is given by the time-dependent vector potential $A_{\text{pmp}}(t)$.

In addition, we consider the core degrees of freedom involved in the Auger transition, i.e.,

[TABLE]

which describes fully localized core states with bare one-particle energy $-E_{c}$, the core-core Coulomb interaction $U_{cc}$ and the core-valence Coulomb interaction $U_{cv}$, assumed as local for simplicity. Accordingly, $n^{\rm(f)}_{i\sigma}=f^{\dagger}_{i\sigma}f_{i\sigma}$ is the occupation-number operator for the core states, and $f^{\dagger}_{i\sigma}$ creates a core electron at lattice site $i$ with spin projection $\sigma$.

For both the pump pulse and the probe pulse (index $\gamma=\text{pmp, prb}$), the electric field with frequency $\Omega_{\gamma}$ and amplitude $\mathcal{E}_{\gamma}^{(0)}$,

[TABLE]

is assumed to have a Gaussian shape with an envelope function,

[TABLE]

characterized by $t_{\gamma}$ and $\tau_{\gamma}$. The corresponding vector potential is obtained by integration:

[TABLE]

Note that for a harmonic pulse with Gaussian envelope, as assumed here, this results in the complex error function which can be numerically efficiently evaluated using the Faddeeva package. Johnson (2013) Furthermore, we define the pump-probe delay as the time difference between the respective pulse peaks: $\Delta t=t_{\text{prb}}-t_{\text{pmp}}$.

The continuum of high-energy scattering states above the vacuum level is represented by

[TABLE]

Here, $\rho\left(\epsilon\right)$ is the corresponding density of states. $a^{\dagger}_{\sigma}\left(\epsilon\right)$ creates an electron with spin projection $\sigma$ and energy $\epsilon$ in the scattering continuum, and $a_{\sigma}\left(\epsilon\right)$ is the corresponding annihilator. These operators obey fermionic anticommutation relations in the form $\{a_{\sigma}\left(\epsilon\right),a^{\dagger}_{\sigma^{\prime}}\left(\epsilon^{\prime}\right)\}=\rho\left(\epsilon\right)^{-1}\delta_{\sigma\sigma^{\prime}}\delta\left(\epsilon-\epsilon^{\prime}\right)$ and are thus dimensionless. Note that we assume two different continua, one for the scattering states occupied by the emitted photoelectron and another one for the Auger electron with one-particle energies $\epsilon_{\phi}$ and $\epsilon\equiv\epsilon_{A}$, respectively (usually suppressed for better readability). We assume that the two continua are well-separated in energy and that the photoelectron and the Auger electron can be treated as distinguishable particles. Effects that result from their non-distinguishability, such as post-collision interactions, Schütte et al. (2012) are beyond the scope of this work.

The system is probed by core-level photoemission, and the corresponding term in the Hamiltonian describing the probe pulse is given by

[TABLE]

Here, we employ the length gauge and assume a constant dipole matrix element $d_{i}\left(\epsilon\right)\approx d_{0}$. The subsequent Auger transition is mediated by a Coulomb-interaction term between core, valence and scattering states:

[TABLE]

where again the transition matrix element is approximated by a constant, i.e., $U^{A}_{i}\left(\epsilon\right)\approx U_{A}$. The full Hamiltonian then reads:

[TABLE]

III Theory

Assuming that the probe is weak and using low-order time-dependent perturbation theory in $V_{\text{prb}}\left(t\right)$, one straightforwardly arrives at the following one-step spectroscopic formula: Rausch (2017); Gunnarsson and Schönhammer (1980)

[TABLE]

where $\rho\left(\epsilon_{\phi}\right)\approx\rho_{\phi}=\text{const}$ is the (constant) density of states of the photoelectron continuum, $n\left(\epsilon\right)=\sum_{\sigma}a^{\dagger}_{\sigma}\left(\epsilon\right)a_{\sigma}\left(\epsilon\right)$ is the occupation-number operator for the Auger-electron continuum, and where

[TABLE]

is a state which is probed at time $t^{\prime}$ (note that $t^{\prime}$ is integrated over in Eq. (11)). Furthermore, $U(\cdot,\cdot)$ is the time-evolution operator corresponding to the Hamiltonian (10) including the pump pulse, but with $V_{\text{prb}}(t)\equiv 0$. At time $t=0$, the system is in its ground state \big{|}\Psi(0)\big{>}. Obviously, without a pump ($\mathcal{E}_{\text{pmp}}^{(0)}=0$), the first time evolution in Eq. (12) becomes trivial, and one would just create a core hole in the ground state at time $t^{\prime}$. For small $U_{A}$ this would essentially recover the Cini-Sawatzky theory of Auger emission, Cini (1977); Sawatzky (1977) with the sole difference of having a probe pulse rather than a continuous beam.

In Eq. (12) the perturbation describing the probe pulse reduces to

[TABLE]

when making use of the sudden approximation, i.e., assuming that a_{\sigma}(\epsilon)\big{|}\Psi(0)\big{>}=0. For an X-ray probe, calculations can be further simplified by exploiting that $\Omega_{\text{prb}}$ is very large on the energy scale of the conduction band set by the hopping $T$. Replacing the fast oscillating cosine square term stemming from Eq. (4) in the integral Eq. (11) by its time-average,

[TABLE]

only the envelope of the probe pulse remains. Since this a positive-definite function, it can be used as a weight function to create custom Gauss-Legendre integration weights using the Golub-Welsh algorithm, Golub and Welsch (1969) so that only few integration points $N_{t}$ are necessary:

[TABLE]

This has also the advantage that the pulse can be made arbitrarily short without numerical problems. In practice, having just $N_{t}=3$ integration points turns out as sufficient.

Another useful approximation relies on the locality of the photoemission and Auger process, i.e., the creation of the core hole and the subsequent Auger decay basically take place at the same lattice site (say, at $i=i_{0}$) opposed to a picture with a coherent superposition over all lattice sites. This approximation has also been justified a posteriori by a general calculation which shows that extra-atomic interferences cancel out when the inverse photoelectron momentum $k_{\phi}$ is much larger than the lattice spacing $a$, i.e., if $k_{\phi}^{-1}\gg a$. Abraham-Ibrahim et al. (1978); Rausch (2017) Furthermore, because of spin degeneracy, we can limit ourselves to just one spin direction (say, $\sigma=\uparrow$), and hence we may replace:

[TABLE]

and

[TABLE]

in Eqs. (8) and (9), respectively, and similarly in Eq. (13). Hence, the equilibrium spectral function, for instance, is of the type \big{<}c^{\dagger}_{i_{0}\downarrow}c^{\dagger}_{i_{0}\uparrow}\delta\left(\epsilon-H\right)c_{i_{0}\uparrow}c_{i_{0}\downarrow}\big{>} rather than \sum_{ij}\big{<}c^{\dagger}_{i\downarrow}c^{\dagger}_{i\uparrow}\delta\left(\epsilon-H\right)c_{j\uparrow}c_{j\downarrow}\big{>}.

The Auger-electron scattering states could be eliminated in the case of weak $U_{A}\ll T$ in very much the same way as for the photoelectron scattering states, which lead us to Eq. (13), namely truncating the Taylor expansion of the time evolution operator after the first order, the $a^{(\dagger)}_{\sigma}\left(\epsilon\right)$ disappear by commutation. We note, however, that the resulting equations are not necessarily easier to handle, since they involve several time propagations and subsequent integrations. It is therefore in fact easier to account for the Auger-electron continuum explicitly, which can be achieved without too much more computational effort and also allows for computations that are non-perturbative with respect to $U_{A}$. In practice, we adapt the mapping to a semi-infinite chain using orthogonal polynomials, as introduced by Chin et al. Chin et al. (2010) The procedure is described in detail in the Appendix A.

IV Core-hole lifetime and maximal propagation time

To test the above framework and to determine some of the parameters, we study the finite core-hole lifetime due to the Auger decay. The discussion of the complete pump-probe setup is shifted to the next section, and thus, for the time being, we assume the pump to be switched off. Using exact diagonalization, we calculate the total cross-section $\sigma=\int d\epsilon~{}\rho\left(\epsilon\right)I_{\text{AES}}(\epsilon)$ [see Eqs. (11) and (38)] as a function of $U_{A}$ with a finite probe pulse rather than a continuous beam. The density of states is assumed as constant, $\rho\left(\epsilon\right)=\rho_{A}=\mbox{const}$.

Assuming that the core-hole lifetime $\tau_{c}$ is long compared to the pulse duration, $\tau_{c}\gg\tau_{\text{prb}}$, there are two analytically accessible limits: (i) For the fully occupied band, $n\equiv\sum_{\sigma}\langle n_{i\sigma}\rangle=2$, the time propagation can be carried out exactly for any $U_{A}$ and $U_{cv}$. (ii) For weak $U_{A}$ and $U_{cv}=0$, perturbation theory can be used.

We take a finite cutoff time $t_{\text{max}}$ [see Eq. (11)] and find the same result for both cases, namely:

[TABLE]

(see Refs. Rausch, 2017; Gunnarsson and Schönhammer, 1980 for details) or, for $\Gamma t_{\text{max}}\ll 1$:

[TABLE]

with the width Gunnarsson and Schönhammer (1980)

[TABLE]

$A_{s}$ is the area under the probe envelope squared, which for the given Gaussian pulse is given by

[TABLE]

The core-hole lifetime can be defined as

[TABLE]

We have checked the validity of the above formula by comparing with the fully numerical evaluation, see Fig. 2. Some deviations are expected as the core hole starts to decay already during probe pulse, but these get smaller with shorter pulse duration $\tau_{\rm prb}$. For $\tau_{\rm prb}=0.1$, the agreement is almost perfect in the two mentioned limits.

Calculations are done for a finite cutoff time $t_{\text{max}}$. Note that the cross-section simply grows in the cutoff time (linearly at first, not directly visible on the log scale). The convergence of the Auger spectrum with increasing $t_{\text{max}}$ is seen in the inset of Fig. 2. Obviously, finer spectral features are resolved at larger $t_{\text{max}}$. In fact, using a finite $t_{\text{max}}$ even before the spectrum is converged is not necessarily unphysical, but simply smoothens out the result further and can be regarded as finite numerical resolution. A system of $L=10$ sites (in a linear chain geometry with periodic boundary conditions) is already enough to observe the basic spectral features, i.e., the band-like part (“B”) and the doublon satellite (“D”) at $U=6$; see also the discussion in the following section and Ref. Rausch and Potthoff, 2016.

To get the core-hole lifetime for the non-trivial case of $n=1$, $U_{cv}\neq 0$ and large or moderate $U_{A}$ to be studied here, we calculate the spectrum for several $t_{\text{max}}$ and employ a two-parameter fit assuming the same functional dependence (and shifting the time-zero to the probe peak):

[TABLE]

Since in the linear regime ($U_{A}\lesssim 0.01$ and $t_{\text{max}}-t_{\rm prb}\lesssim 10$), an exponential fit is ill-defined, we rather use a linear one-parameter fit:

[TABLE]

The results are shown in Fig. 3. Since there are less conduction electrons available to fill the core hole in the half-filled case $n=1$, its lifetime is at least an order of magnitude higher as compared to $n=2$. However, increasing $U_{cv}$ reduces the lifetime, as the core-hole potential is quickly screened by surrounding electrons. As a function of $U_{A}$ (for $U_{A}\gtrsim 0.01$), we find a power-law behavior of the lifetime, $\tau_{c}\propto U_{A}^{\alpha}$, where the exponent increases from $\alpha\approx 1.6$ for $U_{cv}=0$ to $\alpha=2$ for strong $U_{cv}$, i.e., for strong core-hole screening we recover the $n=2$ result.

V Pump-probe spectra for the pure Hubbard model

We proceed to the calculation of pump-probe spectra and at first discuss the pump-pulse parameters in Eq. (4). The pump frequency and the intensity of the pulse will be chosen as to maximize the resulting double occupancy, see below. The pump duration is fixed at $\tau_{\text{pmp}}=1$, i.e., at the time scale set by the inverse nearest-neighbor hopping $1/T$. A typical value of $T\sim 100$ meV would correspond to about $\tau_{\rm pmp}\sim 40$ fs.

For the calculations, i.e., for the time-propagation during the pump pulse in particular, see Eq. (12), we use the CFET algorithm Alvermann and Fehske (2011); Alvermann et al. (2012) based on the Magnus series expansion. Magnus (1954) Results for the double occupancy

[TABLE]

are shown in Fig. 4.

The optimal $\Omega_{\text{pmp}}$ is expected to be given by the energy of an excitation from the middle of lower Hubbard band to the middle of upper Hubbard band. For $U=6$ and half-filling this amounts to $\Omega_{\text{pmp}}\approx W/2+\Delta+W/2\approx 5.45$ where $W=4$ is the noninteracting band width and where the charge gap $\Delta\approx 1.446$ at $U=6$ is taken from the Bethe ansatz solution. Ovchinnikov (1969) Based on the data for $d_{\text{tot}}\left(t\right)$ for frequencies $\Omega_{\text{pmp}}$ in the vicinity of this value, it turns out that the best choice is achieved with a slightly smaller value, $\Omega_{\text{pmp}}\approx 4$, see Fig. 4. With higher frequencies, the resulting double occupancy goes down again, since the electronic system becomes unable to react to such fast driving. The driving field strength $\mathcal{E}^{(0)}_{\text{pmp}}$ has been varied simultaneously and maximizes the double occupancy for $\mathcal{E}^{(0)}_{\text{pmp}}\approx 8$. These optimal values could be somewhat affected by finite-size effects.

Fig. 5 shows the pump-probe Auger spectrum for several delay times $\Delta t$ at $U=6$ and $n=1$, compared with the equilibrium case, i.e., for the pump pulse switched off. We first discuss the equilibrium case for $n=2$ (solid black line). The broad peak at highest kinetic energies of the Auger electron is the band-like part (“B”), where the doublon dissociates into two final-state hole which propagate more or less independently and delocalize. Note that there are some spurious oscillations due to the finite system size. The maximal kinetic energy that the Auger electron can obtain is given by the difference of $E_{0}\left(N-1\right)$, the ground-state energy with one core hole, and $E_{0}\left(N-2\right)$, the ground-state energy with two holes in the valence band. For the full band and strong $U$, this is given by $\epsilon_{\text{max}}=2U+W+E_{c}-U_{cc}$, whereby $\epsilon=2U+E_{c}-U_{cc}$ marks the barycenter of the band-like part which has a width of $2W$. Separated by an energy of about $U$, the spectrum features a satellite corresponding to the case that the doublon does not decay and forms a repulsively bound object (“D”). At fillings $n<2$ (not shown), an additional quadruplon satellite (“Q”) with very low spectral weight appears, separated yet again by about $U$ from the doublon one. Rausch and Potthoff (2016) This corresponds to a bound state of four holes in the system. At half-filling $n=1$, all these features have collapsed into a single peak, see the dashed line for $n=1$, where again some spurious finite-size effects are visible.

We are now prepared to make a comparison with the time-resolved (pump-probe) case. Here, the spectroscopy is performed on a photo-doped system at half-filling. Due to the increased double occupancy, the spectrum shares certain features with the equilibrium spectrum at higher fillings. In particular, we observe an intense peak at the energy position of the doublon satellite. Its flat upper shoulder reaches into the band-like part “B”, while the lower shoulder exhibits a peak at the quadruplon energy (“Q”), albeit an order or magnitude smaller, which will be difficult to observe experimentally. The doublon satellite “D”, on the other hand, is the dominating spectral structure.

However, with regard to the delay time $\Delta t$, we note that there is almost no dependency whatsoever; the spectra calculated for $\Delta t=5$ up to $\Delta t=60$ basically fall on top of each other. The physical reason is that the doublons created during the pump pulse are extremely stable at $U=6$. Their decay time is expected to increase exponentially with increasing Coulomb interaction $U$ in the strong-coupling regime. SGJ+10 ; Eckstein and Werner (2011) Thus, the double occupancy, to which the Auger process is sensitive, more or less stays constant after the pump is over (see also the right inset of Fig. 4), so that it does not matter when, on the typical (femtosecond) timescale of electron dynamics, the probe is applied.

This result has important implications: Apparently, the Auger spectrum is highly robust even when recorded from a system in an excited nonequilibrium state, its overall shape being simply determined by the value of the double occupancy $d_{\text{tot}}$ and its position by the value of $U$. Hence, pump-probe spectra will only yield time-resolved information if either of these quantities fluctuates in time.

The former, however, is a quasi-conserved quantity in the interesting strong-coupling regime, except, of course, during the pump pulse. Hence, there is the possibility to work with overlapping pulses, i.e., probe the system while pumping of doublons is still active. Therewith, the $\Delta t$-dependence would reflect the time dependence of the doublon-pumping process.

As an even more interesting alternative, one may look into cases with a time-dependent (effective) interaction $U(t)$. A time-dependent Hubbard-$U$ can be realized in a number of ways, in particular within the dynamic Hubbard model. This possibility is discussed in the following.

VI Pump-probe spectra for the dynamic Hubbard model

The dynamic Hubbard model Hirsch (2001, 2002); Werner and Eckstein (2016) was introduced by Hirsch to account for the effect that whenever two electrons are on the same site, they will not simply occupy the same single-electron orbital. Rather, the orbital will relax and expand due to the Coulomb interaction between them. To account for this mechanism on the level of the Hubbard model, one can add a bosonic degree of freedom that couples to the double occupancy. Here we adopt the simple approach where this is given by a classical harmonic oscillator. Hence, our extended model now reads as:

[TABLE]

Here, $H_{\text{el}}(t)$ comprises the electronic system, Eq. (10), discussed before. Furthermore,

[TABLE]

is the Hamilton function of the classical oscillators with frequency $\omega_{0}$, and

[TABLE]

is the interaction term with coupling strength $k$. We use rescaled coordinates and momenta as described in the Appendix B.

The choice of the operator $O_{i}$ leads to different physical situations: For $O_{i}=n_{i}$ we obtain the semiclassical Hubbard-Holstein model, which can be of interest for a study of doublon decay due to the presence of phonons, while for $O_{i}=n_{i\uparrow}n_{i\downarrow}$ we obtain the dynamic Hubbard model. This is characterized by a modified, effective Coulomb interaction,

[TABLE]

which in general is time-dependent.

To time-propagate the quantum-classical hybrid system, we exactly map the quantum degrees of freedom onto an additional “virtual” set of coupled classical harmonic oscillators. The dynamics of this large nonlinear classical system is then treated using the standard splitting integrator approach.Ruth (1983); McLachlan (1995); Grochowski and Lesyng (2003); Blanes et al. (2008); Lubich (2008) The detailed procedure is outlined in Appendix B.

In Fig. 6 we show the Auger spectra of the dynamic Hubbard model for the adiabatic case. This can be realized with a frequency $\omega_{0}=0.1$ much smaller than the bare electronic time scale which is given by the inverse nearest-neighbor hopping $1/T=1$. Calculations are done for a moderately strong coupling $k=2$. With this choice of the coupling constant the system is deep in the homogeneous phase as can be seen from the equilibrium phase diagram of the model discussed in the Appendix C. During the pump pulse ($\tau_{\text{pmp}}=1$), electrons are excited into the upper Hubbard band, i.e., double occupations are created. The doubly occupied orbitals then start to oscillate, basically with the frequency $\omega_{0}$, after some initial reaction to the pulse. As shown in the upper inset of the figure, this means that there is relatively long timescale $\tau_{\text{osc}}=2\pi/\omega_{0}\approx 62.8$, on which the effective interaction [see Eq. (29)] varies.

It now matters at which point of this cycle the Auger decay is initiated via the probe pulse. Namely the resulting reduction or enhancement of $U_{i,\text{eff}}\left(t\right)=U+kq_{i}\left(t\right)$ is responsible for a corresponding energy shift of the correlation satellite of the Auger spectrum to lower or higher binding energies, respectively, see the main part of Fig. 6. This observation can be made even quantitative. At $U=6$ and for strong $U$ in general, the correlation satellite basically takes almost the whole spectral weight, so that its position is to a very good approximation given by the barycenter of the complete Auger spectrum:

[TABLE]

Since the core hole takes a certain time to decay, we compare the barycenter with $U_{i,\text{eff}}\left(t\right)$, averaged over some time interval $\rm t_{\rm eff}$,

[TABLE]

see the lower inset of Fig. 6. Taking $t_{\text{eff}}=10$ appears to produce the best results, so that the barycenter even matches $\overline{U}_{i,\text{eff}}$ almost perfectly.

Finally, in Fig. 7 we show results for the antiadiabatic case, with $\omega_{0}=1$. The effective interaction now oscillates on the much shorter timescale of $\tau_{\text{osc}}\approx 6.28$ (upper inset). This oscillation is in fact fast enough to complete at least one cycle within the lifetime of the core hole. Consequently, $\overline{U}_{i,\text{eff}}$ almost averages out to $U_{i,\text{eff}}\left(0\right)$, with a weak remaining oscillation. Again, the barycenter of the spectrum follows $\overline{U}_{i,\text{eff}}$ almost perfectly (see lower inset of Fig. 7).

We conclude that the frequency $\omega_{0}$ manifests itself directly in the delay-time dependence of the main correlation satellite in pump-probe Auger spectroscopy. A significant effect that should be observable is expected if the period is much longer than the core-hole lifetime, i.e., if $\tau_{\text{osc}}=2\pi/\omega_{0}\gg\tau_{c}$.

VII Discussion and summary

To summarize, let us compare CVV Auger-electron spectroscopy from an equilibrium electronic state with nonequilibrium pump-probe AES. We consider a strongly correlated system with a valence electronic structure, described by a Hubbard model with interaction $U$ larger than the bare bandwidth $W$. In the equilibrium case, the local sensitivity of AES can be employed to determine the Hubbard-$U$ as the energy distance of the correlation satellite from the band-like part of the spectrum, since the satellite is shifted in energy by the two-hole binding energy of about $U$. The satellite almost takes the whole spectral weight if $U\gg W$ which makes AES much more sensitive to correlations as compared to single-electron photoemission spectroscopy. This argument builds on the simple single-band Hubbard model and must be refined for real materials (i) in case of systems with partially filled bands where an explicit many-body calculation is necessary to determine the precise satellite position and (ii) by taking into account all relevant low-energy electronic degrees around the Fermi energy, i.e., in by considering a multi-orbital Hubbard-type model.

For the nonequilibrium case, we have seen that the Auger process starts from a photo-doped nonequilibrium state with enhanced double occupancy. In a mode with overlapping pump and probe pulses, the system is probed while still pumping doublons, and the dependence on the pump-probe delay time $\Delta t$ would reflect the time dependence of the doublon-pumping process. After the pump, and in the interesting strong-coupling regime $U\gg W$, however, the double occupancy is a quasi-conserved quantity due to missing phase space for decay in a single scattering event. The important consequence is that the pump-probe Auger spectrum is basically independent of $\Delta t$. Actually, one would expect at least a weak $\Delta t$ dependence due to a relaxation of doublons to the bottom of the upper Hubbard band. Eckstein and Werner (2011); Werner and Eckstein (2016) To see this effect, however, one would have to consider much larger systems or use mean-field-type methods working directly in the thermodynamical limit.

Interestingly, the nonequilibrium spectrum exhibits features that are characteristic of the equilibrium spectrum from a system with higher band filling and thus higher equilibrium double occupancy. This has been verified for a Mott insulator at half-filling where the calculation shows a quadruplon peak in the nonequilibrium spectrum at much higher binding energies compared to the doublon satellite, which is impossible in the equilibrium case at half-filling but shows up at higher fillings. Rausch and Potthoff (2016) As the location of the main doublon satellite and additional features, such as a quadruplon peak, sensitively depend on the average number of doubly occupied sites, pump-probe AES could be used as a tool to measure the doublon concentration in the initial transient state as a function of the pump parameters.

It is well known Werner and Casula (2016) that, in the context of quantitative studies for real materials, the single-band Hubbard model is, at best, an effective low-energy theory whose parameters must be determined in a downfolding process by integrating out higher-energy degrees of freedom, eventually by starting from an all-electron Hamiltonian in the continuum. This downfolding generates a dynamically screened Coulomb interaction, i.e., in the simplest case a frequency- or time-dependent Hubbard-$U$, which leads to important observable effects such as a reduction of the effective bandwidth and of the spectral weight of low-energy excitations as has been pointed out recently. Casula et al. (2012) A first step to capture the effects of dynamical screening can already be done with the dynamic Hubbard model Hirsch (2001) involving a single bosonic mode of frequency $\omega_{0}$ as has been considered here. The screening frequency $\omega_{0}$ is the essential ingredient of the effective low-energy theory as it determines the observable bandwidth renormalizations. Casula et al. (2012) We have demonstrated that pump-probe AES is in fact directly sensitive to the dynamical screening process, subject to the condition that $2\pi/\omega_{0}$ is larger than the core-hole lifetime. The essential observation is that the position of the intense correlation satellite, as a function of the pump-probe delay $\Delta t$, just oscillates with the screening frequency $\omega_{0}$. This is easily understood in a simple picture, where the Hubbard interaction, after integrating out the bosonic mode, becomes time-dependent and oscillates with frequency $\omega_{0}$, such that the satellite position $\sim U_{\rm eff}(t)$ must oscillate with basically the same frequency. Simple estimates as well as state-of-the-art many-body calculations Hirsch (2001); Casula et al. (2012) show that $2\pi/\omega_{0}$ is in the relevant femtosecond regime for strongly correlated materials, such as Mott insulators, and is thus well accessible to modern pump-probe setups.

Acknowledgements.

We would like to thank M. Martins and W. Wurth (University of Hamburg) for instructive discussions. Financial support of this work through the Deutsche Forschungsgemeinschaft within the Sonderforschungsbereich 925 “Light induced dynamics and control of correlated quantum systems” (project B5) is gratefully acknowledged.

Appendix A Orthogonal polynomial mapping

A set of orthogonal polynomials $\{p_{n}\}$ is defined either via a recurrence relation or via a positive-definite weight function, see Ref. Chin et al., 2010. We can treat the density of states of a continuum $\rho\left(\epsilon\right)$ (having a compact support) as such a weight function and define the scalar product between two polynomials $p_{n}$ in the following fashion:

[TABLE]

One can then define new discrete creation and annihilation operators in the following way:

[TABLE]

Note that they are still fermionic: $\{b_{n\sigma},b^{\dagger}_{n^{\prime}\sigma^{\prime}}\}=\delta_{nn^{\prime}}\delta_{\sigma\sigma^{\prime}}$. Using the orthogonality of the polynomials, the inverse transformation is found to be

[TABLE]

The scattering continuum now becomes a semi-infinite chain,

[TABLE]

while the Auger term couples all sites to the first one of the chain:

[TABLE]

Here, $\alpha_{n}$ and $\beta_{n}$ are the recurrence coefficients of the corresponding monic polynomials $\pi_{n}$, and \left|\pi_{0}\right|=\sqrt{\left<\pi_{0},\pi_{0}\right>}=\big{|}U_{A}\big{|}\sqrt{\int\rho\left(\epsilon\right)d\epsilon} is the norm of the first one. In the case of a constant density of states, $\rho\left(\epsilon\right)=\text{const}=\rho_{A}$, and a support limited between $\epsilon_{\text{min}}$ and $\epsilon_{\text{max}}$, the polynomials $p_{n}$ are given as shifted Legendre polynomials whose recurrence relations are explicitly given by

[TABLE]

with $a=\left(\epsilon_{\text{max}}-\epsilon_{\text{min}}\right)/2$ and $b=\left(\epsilon_{\text{max}}+\epsilon_{\text{min}}\right)/2$, so that $\left|\pi_{0}\right|=U_{A}\sqrt{2a\rho_{A}}$. We also note that the total cross-section $\sigma$ can be written as

[TABLE]

where

[TABLE]

is the total particle number operator within the chain.

For the given problem, the initial state is a product state of the system ground state and the chain vacuum. Hence, cutting the chain to a finite length $L_{\text{ch}}$ does not introduce any error. For the time evolution, however, one has to ensure that the occupation of the last site \sum_{\sigma}\big{<}b^{\dagger}_{L_{\text{ch}}-1,\sigma}b_{L_{\text{ch}}-1,\sigma}\big{>} remains small to have a negligible error. Otherwise a reflection from the open end will cause unphysical results.

Note the relation to energy-time uncertainty: With increasing time an increasing number of higher-order polynomials is necessary to describe the dynamics of particle traveling through the chain. This corresponds to a successively better resolution of spectral features.

It is important to realize that $\epsilon_{\text{max}}$ and $\epsilon_{\text{min}}$ do not have to encompass the full support of the spectrum. We can choose to concentrate on a small spectral window only and measure the spectrum there. This implies a smaller effective hopping amplitude, both into the chain $\left|\pi_{0}\right|$ and along the chain $\beta_{n}$ (due to the smaller width $a$). Hence, one needs less chain sites and a smaller Hilbert space, but in turn, more spectral windows are needed to reconstruct the full spectrum. Because of the perturbative treatment of $V_{\text{prb}}$, there can only be a single electron occupying the scattering states, i.e., in the chain. Hence, the Hilbert-space dimension after the probe is given by $\dim\mathcal{H}=\dim\mathcal{H}^{\left(N-1\right)}_{\text{el}}+L_{\text{ch}}\times\dim\mathcal{H}^{\left(N-2\right)}_{\text{el}}$. In the extreme case we can actually have $L_{\text{ch}}=1$ (no increase of Hilbert space at all, apart from the sum), which only allows for a polynomial of degree zero (a constant) and quite a lot of spectral windows are needed to get the full spectrum. In practice, a compromise between memory and operations can be struck, setting $L_{\text{ch}}$ as high as possible for the available memory resources.

Appendix B Time propagation of the dynamic Hubbard model

We couple our quantum system to classical oscillators described by the Hamiltonian function $H_{\rm bos}+H_{\rm int}$ with

[TABLE]

Using the standard rescaling of the oscillator coordinates, i.e., $k=\lambda/\sqrt{m\omega_{0}^{2}}$, $q_{i}=\sqrt{m\omega_{0}^{2}}~{}\tilde{q}_{i}$, $\tilde{p}_{i}=\sqrt{m\omega_{0}^{2}}~{}p_{i}$, allows us to rewrite the terms as

[TABLE]

where $O_{i}=n_{i\uparrow}n_{i\downarrow}$. The superscripts $q,p$ indicate the respective dependencies.

Following Ref. Grochowski and Lesyng, 2003, we choose an orthonormal basis \{\big{|}m\big{>}\} of the many-body Hilbert space to expand the state vector in the Schrödinger picture,

[TABLE]

which can be used to define a classical Hamiltonian function for all the occurring degrees of freedom $\left\{\vec{q},\vec{p},\vec{c}\right\}$ where $\vec{q}\equiv(q_{1},...,q_{L})^{T}$ and $\vec{p}$ analogously and $\vec{c}\equiv(c_{1},c_{2},...c_{M})^{T}$, and $M$ is the Hilbert-space dimension:

[TABLE]

Here $T_{\rm el}$ and $V^{q}_{\rm int}$ are the matrix representations of $H_{\rm el}$ and $H_{\rm int}$ with respect to the chosen basis.

With the Liouville-Operator, given by

[TABLE]

we get the equations of motion for the classical degrees of freedom,

[TABLE]

as well as for the $c$ degrees of freedom, which reproduces the Schrödinger equation:

[TABLE]

The so-called PICKABACK splitting, Grochowski and Lesyng (2003) $\mathcal{L}=\mathcal{L}_{qT}+\mathcal{L}_{pV}$ with

[TABLE]

allows us to perform alternate propagations with these two Liouville operators:

[TABLE]

Here, we employ the fourth-order method ($r=4$), characterized by the choice $a_{1}=a_{4}=1/\left(2x\right)$, $a_{2}=a_{3}=\left(x-1\right)/\left(2x\right)$, $b_{1}=b_{3}=1/x$, $b_{2}=b_{4}=\left(x-2\right)/x$, with $x=2-2^{1/3}$.Ruth (1983); McLachlan (1995); Grochowski and Lesyng (2003); Blanes et al. (2008); Lubich (2008) Note that $H_{\rm int}$ [Eq. (41)] only contains the interaction term with the oscillators and not the Hubbard-$U$, which ensures that no additional error will be introduced in the case $l=0$, which otherwise would result from splitting non-commuting terms of the Hamiltonian.

For our current problem, we have $T_{\text{el}}=T_{\text{el}}(t)$ due to the pump pulse, which must be seen as a complication. We choose a pragmatic solution by taking small time steps (typically, $\delta t=0.01$) for the time propagation during the pump pulse. As this regards a short part of the propagation, it is not very expensive computationally. Once the pulse is over, longer steps are chosen ($\delta t=0.05$). As an error measure one can take the total energy, which becomes constant after the pump is over. The relative error |E-\big{<}H(t)\big{>}|/E is typically around $10^{-7}$ for the parameters above.

Appendix C Ground-state phase diagram of the dynamic Hubbard model

In equilibrium, we have $\dot{p}_{i}=0$, and from Eq. (45) we can see that the displacement of the oscillators is given by q_{i}^{(0)}=-k\big{<}n_{i\uparrow}n_{i\downarrow}\big{>}, i.e., the ground state has to be determined self-consistently. If the coupling strength $k$ becomes comparable to $U$, the system undergoes a phase transition to a symmetry-broken charge-density wave (CDW) phase, where every second site is doubly occupied. This maximizes the double occupancy and reduces the total energy compared to that of a homogeneous state. The $k$-$U$ phase diagram for a system of $L=10$ sites is shown in Fig. 8. One observes that the phase boundary is very sharp despite the small system size.