Microscopic bosonization of band structures: X-ray processes beyond the Fermi edge

TL;DR

This paper extends bosonization to microscopic lattice models, enabling analysis of band structures and x-ray processes across all energy scales, including strong interactions, by representing fermionic kinetic energy in the energy domain.

Contribution

It introduces a microscopic bosonization formalism that overcomes low-energy limitations, applicable to entire band structures and capable of handling nonlinear bosonic interactions.

Findings

Successfully describes x-ray edge singularity at emission threshold.

Quantitative agreement with exact diagonalization results.

Reveals correlation signatures beyond the band bottom.

Abstract

Bosonization provides a powerful analytical framework to deal with one-dimensional strongly interacting fermion systems, which makes it a cornerstone in quantum many-body theory. Yet, this success comes at the expense of using effective infrared parameters, and restricting the description to low energy states near the Fermi level. We propose a radical extension of the bosonization technique that overcomes both limitations, allowing computations with microscopic lattice Hamiltonians, from the Fermi level down to the bottom of the band. The formalism rests on the simple idea of representing the fermion kinetic term in the energy domain, after which it can be expressed in terms of free bosonic degrees of freedom. As a result, one- and two-body fermionic scattering processes generate anharmonic boson-boson interactions, even in the forward channel. We show that up to moderate interaction strengths, these nonlinearities can be treated analytically at all energy scales, using the x-ray emission problem as a showcase. In the strong interaction regime, we employ a systematic variational solution of the bosonic theory, and obtain results that agree quantitatively with an exact diagonalization of the original one-particle fermionic model. This provides a proof of the fully microscopic character of bosonization on all energy scales for an arbitrary band structure. Besides recovering the known x-ray edge singularity at the emission threshold, we find strong signatures of correlations even at emission frequencies beyond the band bottom.

Equations

$$H_0={\int }_0^{\pi }dq\epsilon (q){\tilde{c}}_q^{†}{\tilde{c}}_q,$$

$$c_{\epsilon }^{†}=\sqrt{2\pi N(\epsilon )}{\tilde{c}}_{q(\epsilon )}^{†},$$

$$H_0={\int }_0^{\pi }dq\epsilon (q):{\tilde{c}}_q^{†}{\tilde{c}}_q:={\int }_{-D}^{D}d\epsilon \epsilon :c_{\epsilon }^{†}{c}_{\epsilon }:.$$

$$H_{0,\mathrm{enl}}$$

$$\psi (\tau )={\int }_{-\infty }^{\infty }\frac{d\epsilon }{\sqrt{2\pi }}{e}^{i\epsilon \tau }{c}_{\epsilon },$$

$$H_{0,\mathrm{enl}}=-i{\int }_{-\infty }^{+\infty }d\tau {\psi }^{†}(\tau ){\partial }_{\tau }\psi (\tau ).$$

$$b_{\epsilon }=\frac{1}{\sqrt{\epsilon }}{\int }_{-\infty }^{\infty }d\tau e^{-i\epsilon \tau }{\psi }^{†}(\tau )\psi (\tau ),\epsilon >0,$$

$$H_{0,\mathrm{enl}}={\int }_0^{\infty }d\epsilon \epsilon b_{\epsilon }^{†}{b}_{\epsilon },$$

$${\tilde{\psi }}_0^{†}={\int }_0^{\pi }\frac{dq}{\sqrt{\pi }}{\tilde{c}}_q^{†}={\int }_{-D}^{D}d\epsilon \sqrt{2N(\epsilon )}{c}_{\epsilon }^{†},$$

$${\tilde{\psi }}_0^{†}={\int }_{-D}^{D}d\epsilon \sqrt{2N(\epsilon )}{\int }_{-\infty }^{+\infty }\frac{d\tau }{\sqrt{2\pi }}{e}^{i\epsilon \tau }{\psi }^{†}(\tau )=2\sqrt{\pi }{\int }_{-\infty }^{+\infty }d\tau \Delta (\tau ){\psi }^{†}(\tau ),$$

$$\Delta (\tau )={\int }_{-D}^{+D}\frac{d\epsilon }{2\pi }{e}^{i\epsilon \tau }\sqrt{N(\epsilon )}.$$

$$\psi (\tau )=\frac{U}{\sqrt{2\pi a}}{e}^{i{\epsilon }_F\tau +{\int }_0^{\infty }\frac{d\epsilon }{\sqrt{\epsilon }}{e}^{-a\epsilon /2}\left(b_{\epsilon }{e}^{i\epsilon \tau }-b_{\epsilon }^{†}{e}^{-i\epsilon \tau }\right)},$$

$${\tilde{\psi }}_0^{†}=\sqrt{\frac{2}{a}}{U}^{†}{\int }_{-\infty }^{+\infty }d\tau \Delta (\tau {)}^{\ast }{e}^{-i{\epsilon }_F\tau -{\int }_0^{\infty }\frac{d\epsilon }{\sqrt{\epsilon }}{e}^{-a\epsilon /2}\left(b_{\epsilon }{e}^{i\epsilon \tau }-b_{\epsilon }^{†}{e}^{-i\epsilon \tau }\right)}.$$

$${\tilde{\psi }}_0^{†}{\tilde{\psi }}_0=\frac{1}{2}+2{\int }_{-\infty }^{\infty }d{\tau }_1{\int }_{-\infty }^{\infty }d{\tau }_2\mathcal{P}\frac{e^{i{\epsilon }_F({\tau }_2-{\tau }_1)}}{i({\tau }_2-{\tau }_1)}\Delta ({\tau }_1)\Delta ({\tau }_2{)}^{\ast }{e}^{i\left[{\phi }^{†}({\tau }_1)-{\phi }^{†}({\tau }_2)\right]}{e}^{i\left[\phi ({\tau }_1)-\phi ({\tau }_2)\right]},$$

$$\phi (\tau )=-i{\int }_0^{\infty }\frac{d\epsilon }{\sqrt{\epsilon }}{e}^{i\epsilon \tau }{b}_{\epsilon }.$$

$$H=\left(H_V+\frac{V}{2}-{\epsilon }_d\right)d{d}^{†}+H_0{d}^{†}d,$$

$$H_V={\int }_0^{\pi }dq\epsilon (q):{\tilde{c}}_q^{†}{\tilde{c}}_q:+V\left({\tilde{\psi }}_0^{†}{\tilde{\psi }}_0-\frac{1}{2}\right),$$

$${\tilde{\psi }}_0=\frac{1}{\sqrt{\pi }}{\int }_0^{\pi }dq{\tilde{c}}_q.$$

$$W={\gamma }^2\sum _{\nu }{\left|⟨{\Psi }_{\nu }^0|{\tilde{\psi }}_0|{\Psi }_0^V⟩\right|}^2\delta (E_{\nu }^0+{\epsilon }_d+\omega -E_0^V),$$

$$\epsilon =\omega +{\epsilon }_d-E_0^V,$$

$$\delta (\epsilon +E_{\nu }^0)=\frac{1}{\pi }\mathrm{Re}{\int }_0^{\infty }dt{e}^{-i(\epsilon +E_{\nu }^0)t},$$

$$W(\epsilon )$$

$$P(t)$$

$$W(\epsilon )=2{\gamma }^{2}N({\epsilon }_F+\epsilon )\theta (-\epsilon ),$$

$$N(\epsilon )=\frac{1}{2\pi }\frac{dq(\epsilon )}{d\epsilon }$$

$$W(\epsilon )\propto {\left|\frac{\epsilon }{{\epsilon }_F}\right|}^{[(f_0+1{)}^2-1]}\theta (-\epsilon ).$$

$$\varphi (\epsilon )=\arctan \left[\frac{2\pi N(\epsilon )V}{1+2V\mathcal{P}{\int }_{-\infty }^{\infty }d{\epsilon }^{\prime }\frac{N({\epsilon }^{\prime })}{{\epsilon }^{\prime }-\epsilon }}\right]$$

$$H_{V,\mathrm{enl}}={\int }_0^{\infty }d\epsilon \epsilon b_{\epsilon }^{†}{b}_{\epsilon }+2V{\int }_{-\infty }^{\infty }d{\tau }_1{\int }_{-\infty }^{\infty }d{\tau }_2\mathcal{P}\frac{e^{i{\epsilon }_F({\tau }_2-{\tau }_1)}}{i({\tau }_2-{\tau }_1)}\Delta ({\tau }_1)\Delta ({\tau }_2{)}^{\ast }{e}^{i\left[{\phi }^{†}({\tau }_1)-{\phi }^{†}({\tau }_2)\right]}{e}^{i\left[\phi ({\tau }_1)-\phi ({\tau }_2)\right]}..$$

$$\Phi (\epsilon )=-\frac{\sqrt{\epsilon }}{2}⟨b_{\epsilon }+b_{\epsilon }^{†}⟩$$

$$H_{V,\mathrm{enl}}\simeq {\int }_0^{\infty }d\epsilon \epsilon \left[b^{†}+\frac{f_1(\epsilon )}{\sqrt{\epsilon }}\right]\left[b^{†}+\frac{f_1(\epsilon )}{\sqrt{\epsilon }}\right]+V\left(⟨n⟩-\frac{1}{2}\right)-{\int }_0^{\infty }d\epsilon [f_1(\epsilon ){]}^2,$$

$$f_1(\epsilon )=\frac{2V}{\epsilon }{\int }_0^{\infty }d\omega \sqrt{N({\epsilon }_F-\omega )}\left[\sqrt{N({\epsilon }_F-\omega +\epsilon )}-\sqrt{N({\epsilon }_F-\omega -\epsilon )}\right].$$

Full text

1 Introduction

The general description of fermions at finite density constitutes a central problem in physics, requiring a microscopic understanding of a macroscopically large number of interacting particles. Standard ways to simplify the task often rely on a description of the low energy excitations above the many-body ground state in terms of weakly interacting quasi-particles with effective parameters. Predicting the behavior of excitations beyond this infrared regime remains challenging, and numerical techniques are typically used, with various limitations regarding the number of handled particles, the range of accessible temperatures, the strength of the Coulomb interaction, and the resolution of fine spectroscopic structures [1, 2, 3, 4]. However, modelling the response of strongly correlated solids up to large excitation energies has become a pressing issue, especially with advances in spectroscopic methods such as inelastic neutron scattering [5], angle-resolved photoemission [6], and resonant inelastic x-ray scattering [7]. In the field of cold atoms, various techniques such as Bragg [8] and momentum-resolved Raman [9] spectroscopies have also been developed to probe the spectrum of interacting Fermi gases.

At finite densities and low energies, one-dimensional fermion systems are elegantly described in terms of density fluctuations that behave like bosonic particles [10, 11, 12]. The mapping between fermions and bosons, known as bosonization, is usually applied in conjunction with a linearization of the fermion dispersion relation, which produces a low energy effective theory in which wavelengths of the order of the lattice constant are integrated out. Since information about the crystal lattice structure is lost in the process, it is difficult to describe phenomena resulting from an interplay between band stucture and many-body correlations. Incorporating band structure in the conventional bosonic picture introduces boson-boson interactions with a divergent perturbative expansion [13, 14, 15, 16, 17]. This has cast serious doubt on the usefulness of the bosonic picture beyond the infrared limit. In this paper, we take a new approach to bosonization by mapping free fermions with an arbitrary spectrum onto well-crafted free bosonic degrees of freedom. Unlike conventional bosonization, our approach is readily applicable to microscopic rather than low-energy effective models, either on tight-binding lattices or for continuous wave equations with a non-linear dispersion.

To demonstrate the utility of our approach, we investigate the stimulated emission of x-rays when an electron in a metal decays into a core-orbital inside one of the lattice ions, the so-called x-ray edge problem [18, 19] . For low frequency emitted photons (close to the Fermi edge), conventional bosonization provides an elegant answer [20]. There are however interesting effects, such as multiple-electron processes that produce a non-zero emission rate at frequencies beyond the band edge, that conventional bosonization cannot hope to capture. Using our new approach to bosoniztion, we calculate the full emission spectrum, from the low energy threshold (the edge) up to the high energy regime where band structure plays an important role. For weak interactions, we obtain for the first time an analytical solution that incorporates band structure effects. Beyond this regime, we implement a non-perturbative variational approach that yields near-perfect agreement with brute force numerical diagonalization (the current state of the art). Besides shedding light on the nature of many-body correlations and the natural degrees of freedom in the problem, our approach is more efficient than existing ones. The scaling of computation time with system size $\Omega$ is $\Omega\ln\Omega$ as opposed the the significantly more expensive $\Omega^{3}$ of the brute force approach [21, 22].

While we restrict our analysis to the specific physical problem of x-ray emission, in order to showcase our method, the ideas developed are general. A large class of low-dimensional fermion systems can be mapped onto new bosonic problems that faithfully reproduce physics from the infrared to the ultraviolet. While generically the bosonic models are interacting, they hold great advantages over the original fermionic description, because here the interactions do not invalidate a quasiparticle description in terms of bosons. In this work, we have developed a systematic variational treatment of the bosonized theory, since the x-ray response can be obtained solely from the knowledge of the underlying wavefunction. We anticipate that microscopic bosonization provides bosonic models that may be more efficiently treated by standard numerical methods than their original fermionic versions, since the interplay of band-structure and x-ray thresholds or Tomonaga-Luttinger physics is already incorporated at quadratic order in the bosonic theory.

The paper is organized as follows. Section 2 develops a general microscopic bosonization approach to cope with arbitrary electronic band structures, with a detailed reformulation of standards electronic operators in bosonic form. Section 3 discusses general aspects of the x-ray problem, presenting the microscopic Hamiltonian, the relevant optical response function, and the low-energy physics resulting from the orthogonality catastrophe. In Section 4, we use the microscopic bosonization method to derive an analytical solution of the x-ray edge problem for an arbitrary band structure, in the case where the interaction strength is not too strong. We show how conventional bosonization results are recovered for the threshold singularities in the case of a linear dispersion. In Section 5, the case of large interaction is then addressed by more advanced many-body wave function methods. This involves natural extensions of the analytical theory based on variational coherent states of the normal bosonic modes. In closing we provide a perspective that underlines the many possibilities that microscopic bosonization could open in the field of strongly correlated systems.

2 Microscopic bosonization formalism for lattice models

2.1 Arbitrary electronic spectrum in linear form

Our starting point is the free Hamiltonian for a single species of fermion in a one-dimensional crystal

[TABLE]

with $\tilde{c}^{\dagger}_{q}$ the fermion creation operator for the even band orbital that is a linear combination of Bloch states with crystal momenta $q$ and $-q$. We choose to work here with even modes (discarding odd modes) anticipating the fact that we will later consider a single-site impurity in a time-reversal invariant crystal, that scatters only even states. Alternatively, it is also possible to work with chiral species $\tilde{c}^{\dagger}_{q>0}$ and $\tilde{c}^{\dagger}_{q<0}$, if the problem at hand requires. Operators associated with the momentum and position representations are denoted with tildes, in order to distinguish them from the operators associated with the energy representation (and its conjugate time representation defined below). We will assume that the dispersion relation $\varepsilon(q)$ increases monotonically for $q\in(0,\pi)$, but that its precise form is arbitrary. For a non-monotonous dispersion relation, several electronic sub-bands would have to be introduced.

Instead of linearizing the dispersion relation around the Fermi energy, as one normally does when bosonizing, we invert it to express the momentum $q$ as a function of energy $\varepsilon$. We take the bottom and top of the band to lie at $-D$ and $D$ respectively. For $\varepsilon\in(-D,D)$, we thus define new fermionic operators

[TABLE]

where $N(\varepsilon)=\frac{1}{2\pi}\frac{dq(\varepsilon)}{d\varepsilon}$ is the density of even states per unit length, so that their anticommutator reads $\{c_{\varepsilon},c_{\varepsilon^{\prime}}^{\dagger}\}=\delta(\varepsilon-\varepsilon^{\prime})$ and

[TABLE]

This is formally similar to the operator $\int_{-\infty}^{\infty}dq\,q\,:\tilde{c}_{q}^{\dagger}\tilde{c}_{q}:$ that is encoutered in the case of the linear dispersion [10], but with one important difference: The band of the $c_{\varepsilon}^{\dagger}$ fermions is bounded. However, one can always view a given system as a sub-system that is uncoupled from the rest of a larger system. In order to extend the energy spectrum infinitely downwards and upwards, we introduce spectator fermions created by operators $c_{\varepsilon}^{\dagger}$, with $|\varepsilon|>D$, that anti-commute with the physical fermion operators within the band, and whose Hamiltonian reads $H_{\rm spec}=\int_{|\varepsilon|>D}\!d\varepsilon\;\varepsilon:c_{\varepsilon}^{\dagger}c_{\varepsilon}:$. Here normal ordering is with respect to the state in which orbitals with $\varepsilon<-D$ are occupied and ones with $\varepsilon>D$ are empty. Rescaling transformations such as $\tilde{c}_{q}\to c_{\varepsilon}$ is standard in the analysis of infinite systems [3], but enlarging the Hilbert space with spectator orbitals is a new ingredient that is crucial for the microscopic bosonization approach. The Hamiltonian for the enlarged system reads

[TABLE]

Due to the normal ordering of the spectator fermion operators, $H_{0,{\rm enl}}$ has the same ground state energy as $H_{0}$. In the above expression (4), it is manifest that $H_{0,{\rm enl}}$ does not couple the spectator fermions and the band fermions. Thus, in the physical sector, the enlarged free Hamiltonian is fully equivalent to the original fermionic model.

2.2 Free bosonic representation of the electronic band

The second important step is to obtain a representation of the free electronic Hamiltonian $H_{0,{\rm enl}}$ (with an arbitrary non-linear dispersion) in terms of free bosonic degrees of freedom. With the energy representation (4) as our starting point, we define conjugate time representation operators through the Fourier transform:

[TABLE]

so that $\{\psi(\tau),\psi^{\dagger}(\tau^{\prime})\}=\delta(\tau-\tau^{\prime})$. Thus, in the time representation, the enlarged electronic band Hamiltonian (4) reads

[TABLE]

The presence of the first order differential operator $-i\partial_{\tau}$ in Eq. (6), reminiscent of a linearly dispersing wave equation, paves the way for expressing the exact kinetic energy of band electrons as a single-particle additive bosonic operator. In analogy with the standard bosonization identities [10] for a linearized spectrum in momentum representation, we construct bosonic annihilation operators

[TABLE]

that, in our case, are associated with the energy representation. This definition constitutes the essence of the microscopic bosonization method. In terms of these bosonic operators, the free electronic kinetic term (6) reads simply

[TABLE]

similarly to the case of the linearized theory.

2.3 Construction of the lattice fermionic operators

The final step of the microscopic bosonization is to express the local electronic fields of the original lattice in terms of the bosons associated with the energy representation (7). Without loss of generality, we consider the electronic field for the site at the origin:

[TABLE]

where we have used Eq. (2) to go from the momentum to the energy representation. From the Fourier transform (5), we can then obtain a faithful expression of the lattice field in terms of the time-local field:

[TABLE]

where we have defined

[TABLE]

Note that for a strictly infinite linear dispersion, the function $\Delta(\tau)$ becomes nearly local as the density of states $N(\varepsilon)$ is constant on all energy scales (up to the inverse of the short time cut-off given by the bosonization regulator $a$).

Finally, we can express the fermionic operators in the conjugate time representation using the usual bosonization identity [10]

[TABLE]

in terms of the bosonic modes in the energy domain. In this expression, $U$ is the usual unitary Klein factor, and $\varepsilon_{F}$ is the Fermi energy. It is crucial to note that the ultraviolet cut-off $1/a$ is assumed here to be much larger than the band width $2D$. This is in sharp contrast to the case of the linearized system, were $1/a$ represents the boundary between kept and discarded modes, which must lie sufficiently deep inside the band, for the linear approximation to the dispersion relation to hold. In the expressions that appear below, we will take the limit $a\to 0^{+}$, which here is well-behaved, because now the half-bandwidth $D$ plays the role of a natural physical cutoff. From Eqs. (10-12) we obtain a microscopic expression of the lattice fermion creation operator in terms of the collective bosonic fields:

[TABLE]

It is also useful to provide an explicit expression for the local electronic density operator on the lattice:

[TABLE]

in terms of the bosonic operators in the time domain:

[TABLE]

Note that in Eq. (14), boson operators are normal-ordered and limit $a\to 0^{+}$ has been taken. Eqs. (8), (13) and (14) are the central result of this section, allowing us to express faithfully standard lattice fermion operators in terms of collective boson modes that leave the kinetic energy as a purely harmonic term. The approach that we developed here is thus completely general, and could be applied to bosonize arbitrary one-dimensional models (including bulk interaction among fermions). We now turn in the next section to the simplest testbed for these ideas, namely the x-ray edge singularity in metals. This will allow us to extend previous analytical techniques [20] for an arbitrary density of states, and demonstrate the microscopic equivalence on all energy scales of the lattice fermionic model to our bosonic representation.

3 Formulation and review of the x-ray edge problem

3.1 Motivation

The x-ray edge problem concerns the stimulated transition of a fermion between a Fermi sea and a localized orbital  [18, 19]. In the solid state context, the electromagnetic radiation that stimulates the transition consists of soft x-rays, and we will use this terminology, regardless of the actual wavelength. Due to initial or final state interactions between the Fermi sea and the localized orbital, the x-ray transition rates involve the overlap between the calm Fermi sea, and one that is agitated. Such overlaps, and the associated Anderson orthogonality catastrophe [23], lead to a power-law singularity in x-ray emission and absorption spectra, close to the Fermi threshold. The low-energy physics of the Fermi edge singularity was elucidated theoretically in the Sixties [24, 25, 26, 27] by Nozières and others. Full spectroscopic calculations can be performed nowadays using brute force numerical diagonalization on large systems [21, 22], or using approximate diagrammatic methods [28]. One physical question that we wish to address here concerns the spectroscopic signatures of many-body physics away from the Fermi level, which cannot be described by Nozières’s low energy theory. Besides their obvious relevance in realistic aspects of x-ray spectra [29], it is worth mentioning that initial or final state interaction effects due to dynamical impurities can also occur in mesoscopic devices [30] and in cold atom gases [22, 31].

Historically, the x-ray problem represents an early success for the application of bosonization. The bosonic description [20] proves more succinct than the fermionic one [27], provided one accepts as starting point an effective low energy model with renormalized parameters, which are often difficult to express in terms of the original microscopic ones. In contrast, the parameters that appear in our bosonization approach, are the original microscopic ones. In principle, our microscopic bosonization approach could be used in conjunction with a variety of existing numerical techniques to compute density matrices, partition functions, or arbitrary dynamical response functions. Here we take a non-perturbative variational approach that is also physically intuitive. Since coherent states constitute a natural language for orthogonality catastrophe physics, we propose variational states based on multi-mode bosonic coherent states, which provide excellent x-ray emission spectra on all energy scales, and an appealing description of many-fermion states in the presence of dynamic impurities.

3.2 Modelling of a Fermi sea with a dynamical impurity

We start with a description of our model, which consists of a single band of a one-dimensional crystal and a nearby localized orbital. (In the atomic gas context, the crystal would be engineered using an optical lattice.) The band is partially filled with non-interacting fermionic particles. The localized orbital can either be empty or filled. In the absence of x-ray stimulation (which we discuss in the next subsection), it is assumed that particles do not tunnel between the band and the localized orbital. When the localized orbital is empty, fermions in the crystal undergo potential scattering, and when the localized orbital is filled, the fermions in the crystal are not scattered. This form of interaction naturally arises in the solid state context if the localized orbital represents a state in a core shell of one of the lattice ions. In this case, the empty orbital state corresponds to a core hole that produces an attractive Coulomb potential. From a mathematical point of view however, the roles of the empty and filled states of the localized orbital can be reversed without affecting the applicability of our method. Since spin plays no role, we consider a single spin species. We assume that the static potential of the empty localized orbital is localized to site zero of the crystal. We assume inversion symmetry, so that spatial parity is a good quantum number, and only even parity band-orbitals couple to the localized orbital. We remove odd parity band-orbitals from our description at the outset. We define normal ordering $:\ldots:$ with respect to the clean Fermi sea associated with the band when it does not interact with the localized orbital. We measure crystal momentum $q$ in inverse units of the lattice constant. The system is thus described by the unperturbed Hamiltonian

[TABLE]

where $H_{0}=H_{V=0}$ with

[TABLE]

and

[TABLE]

The operator $d^{\dagger}$ creates a fermion in a localized orbital with energy $\varepsilon_{d}$. The operator $\tilde{c}_{q}^{\dagger}$ creates a fermion with positive momentum $q$ in the even band orbital, as defined in the previous section. The operator $\tilde{\psi}_{0}^{\dagger}$, also defined in the previous section, creates a fermion in the Wannier orbital associated with the site zero of the lattice. Again, the form of the electronic density of state is arbitrary. In the soft x-ray emission and absorption problem in metals, the orbital at energy $\varepsilon_{d}$ describes a core level, so that it lies below the bottom of the band in energy, and the interaction is attractive, i.e. $V<0$. We are however interested in the model for its own sake, and will not place any restrictions on $\varepsilon_{d}$ or $V$. The constant term $-V/2$ in the definition of $H_{V}$ is included for later convenience.

3.3 X-ray transition rate

We consider an initial state in which the localized orbital interacts with the particles in the crystal and the particles in the crystal have reached zero temperature equilibrium. The system is then subjected to incoherent electromagnetic radiation at frequency $|\omega|$, that stimulates a transition in which a fermion tunnels between the band and the localized orbital. We assume that the radiation only stimulates tunnelling between the localized orbital and the Wannier orbital associated with the site zero of the lattice, i.e. we neglect a possible momentum dependence of the optical matrix elements.

According to Fermi’s golden rule, the transition rate for this process is given by

[TABLE]

where $\gamma$ is a tunnelling amplitude (with dimensions of energy). Here, $\left|\Psi_{0}^{V}\right>$ is the ground state (Fermi sea) of $H_{V}$, and $E_{0}^{V}$ the associated energy, in the sector of Fock space that contains the same number of particles as when all single particle orbitals of $H_{0}$ up to the Fermi energy are filled. The states $\left|\Psi_{\nu}^{0}\right>$ are the complete set of eigenstates of $H_{0}$, in the sector of Fock space with one less particle than $\left|\Psi_{0}^{V}\right>$, and $E_{\nu}^{0}$ are their energies. When $\omega>0$ in (19), tunnelling is accompanied by the stimulated emission of a photon of frequency $\omega$, and when $\omega<0$, by the absorption of a photon of energy $-\omega$. We define a shifted frequency

[TABLE]

and consider the transition rate $W$ as a function of $\varepsilon$. Note that, due to normal ordering, the ground state energy of $H_{0}$ is zero. Thus, the transition rate $W(\varepsilon)$ vanishes for $\varepsilon>0$. At the threshold $\varepsilon=0$, the final state of the crystal is the clean Fermi sea.

By writing the $\delta$-function in (19) as

[TABLE]

and using the completeness of the states $\left|\Psi_{\nu}^{0}\right>$, we can rewrite formula (19) for the transition rate as

[TABLE]

For $V=0$, the transition rate evaluates to

[TABLE]

where

[TABLE]

is the density of even states per unit length, and $\varepsilon_{F}$ is the Fermi energy. At the threshold value of $\varepsilon=0$, the transition rate makes a discontinuous jump. Just below the threshold, the transition rate is finite. For the one-dimensional band that we consider, the transition rate diverges at $\varepsilon$ equal to the energy difference between the bottom of the band and the Fermi energy. This is due to the van Hove singularity in $N(\varepsilon)$ at the bottom of the band. For lower energies, the transition rate is zero (to order $\gamma^{2}$).

For non-zero $V$, the transition rate develops a prominent feature close to the $\varepsilon=0$ threshold. At energies $\varepsilon<0$ such that $|\varepsilon|$ is much less than both the energy differences between the top edge of the band and the Fermi energy, and between the Fermi energy and the bottom edge of the band, the transition rate acquires a power-law form (the so-called Fermi-edge singularity) [24, 25, 18]

[TABLE]

Here $f_{0}=\phi(\varepsilon_{F})/\pi$ with

[TABLE]

the exact scattering phase shift that the local potential $V\tilde{\psi}_{0}^{\dagger}\tilde{\psi}_{0}$ induces on a fermion in the band incident at energy $\varepsilon$, so that, according the the Friedel sum rule, $f_{0}$ is the average number of particles displaced by the potential. For negative $\phi$, the transition rate close to threshold is enhanced, while for positive $\phi$, it is a suppressed.

Further below the $\varepsilon=0$ threshold, the transition rate is also modified in important ways, but these effects are not accounted for by Nozières’s approach. For instance, transitions to final states containing multiple particle-hole excitations above the clean Fermi sea will cause a non-zero transition rate below the bottom of the band. In order to exactly calculate $W(\varepsilon)$ at values of $\varepsilon$ that are finite compared to the bandwidth, one has to employ brute force numerics [22]. Each evaluation of $P(t)$ in Eq. (23) involves the numerical evaluation of a Slater determinant, and other matrix operations, for which the computation time scales as the third power of the number of particles in the system. We provide more detail about this approach in the supplementary material that accompanies this Article.

Apart from the transition that we consider, namely one in which the initial state of the localized orbital interacts with the impurity and the final one does not, one can also consider the case where the interaction occurs in the final state. In this case, the initial state of the band is the clean Fermi sea $\left|\Psi_{0}^{0}\right>$, and after a particle is added to site zero, the system evolves with $H_{V}$. The wave function method we develop below is tailored to time evolution with $H_{0}$ rather than $H_{V}$, and we therefore consider here only the initial and not the final state interaction situation.

4 Weak interaction solution of the x-ray edge problem on all energy scales

4.1 Bosonized form of the x-ray edge model with arbitrary density of states

We express the electronic x-ray edge Hamiltonian (17) in terms of collective boson degrees of freedom, using the equations (8) and (14). We readily obtain

[TABLE]

We stress that the bosonic form of $H_{V,{\rm enl}}$ does not couple band fermions and spectator fermions, although this is no longer manifest, as these states are mixed in a complicated way within the bosonic fields. In its bosonized form, Hamiltonian (28) can in principle be analyzed in a variety of ways. For instance, numerical renormalization group [32] or quantum Monte Carlo [33, 34, 35, 36] calculations have previously been developed to deal with similar Hamiltonians. We choose to focus in what follows on wave function-based methods, because according to Eq. (23), we only need to find the ground state of $H_{V,{\rm enl}}$ in order to calculate the transition rate $W(\varepsilon)$. Since $H_{V,{\rm enl}}$ is not quadratic in the bosonic degrees of freedom, there is probably no simple expression for the exact ground state in the bosonic language. Our main task will be in what follows to provide controlled analytical and numerical computation of this bosonic many-body state.

4.2 Analytical solution at moderate interaction for an arbitrary band

We provide here a controlled analytical solution of the x-ray problem in an arbitrary band and on all energy scales (even beyond the band edge), provided that the impurity interaction is weak enough. The orthogonality catastrophe induced by the impurity will be treated non-perturbatively, thanks to the bosonic representation (28) of the problem. This solution will also serve as a starting point for further variational refinements that we will develop in the following section in order to address the regime of strong coupling. The starting point is the observation that for a weak interaction strength $V\ll D$, the bosonic modes $b_{\varepsilon}^{\dagger}$ in (28) fluctuate very mildly around their undisplaced configuration $b_{\varepsilon}^{\dagger}=0$. Thus, it is legitimate to expand the exponential term in (28) to first order in the bosonic modes. This produces a quadratic bosonic Hamiltonian, reminiscent of the case with linear dispersion [20, 10], while still encoding non-trivial features of the full band structure, in contrast to the usual bosonization solution. In this approximation, the exact Hamiltonian matrix element for creating a single particle-hole pair of energy $\varepsilon$ is included in the description, for all energies $\varepsilon$ up the extreme ultraviolet limit $2D$. The x-ray problem is thus accurately captured at an arbitrary energy, as long as $|V|$ is sufficiently small.

In order to visualize the configuration of the bosonic degrees of freedom it is useful to define a rescaled average bosonic displacement

[TABLE]

with respect to the exact bosonic ground state, a quantity which vanishes when the potential $V$ is turned off. In the case of a linear dispersion relation, $\Phi(\varepsilon=v_{F}k)$ corresponds to the $k$-component of the Fourier transform of the fermion density ($v_{F}$ is the Fermi velocity). In general, $\Phi(0)$ equals the charge displaced by the impurity potential, which according to the Friedel sum rule, equals the phase shift at the Fermi energy, divided by $\pi$. The bosonic displacement can be calculated by exact diagonalization, in the original fermionic representation of $H_{V}$, by using the definition (7) for $b_{\varepsilon}$. This reveals that its maximum value as a function of $\varepsilon$ is always of the order of $\Phi(0)$. Thus, $\Phi(\varepsilon)$ is small when the impurity induces a small phase shift at the Fermi energy, namely when $V$ is small enough. This confirms our initial argument that expanding the exponentials $e^{-i\left[\varphi^{\dagger}(\tau_{1})-\varphi^{\dagger}(\tau_{2})\right]}e^{-i\left[\varphi(\tau_{1})-\varphi(\tau_{2})\right]}$ in (28) to first order in boson operators will yield an accurate approximation. The approximation should work particularly well for a particle-hole symmetric band at half-filling, where the dropped terms in the expansion contain at least three normal-ordered boson operators. Whereas the linear dispersion approximation [20, 10] can only predict $\Phi(0)$ and the power-law behavior of the transition rate $W(\varepsilon)$ close to the Fermi edge threshold, here we expect to obtain quantitatively correct results for $\Phi(\varepsilon)$ and $W(\varepsilon)$ for all $\varepsilon$ as long as the potential strength $V$ is sufficiently small.

Expanding the potential term in (28) to first order in $b_{\varepsilon}$ and $b_{\varepsilon}^{\dagger}$, we obtain the lowest order Hamiltonian

[TABLE]

where $\left<n\right>=2\int_{-D}^{\varepsilon_{F}}d\varepsilon\,N(\varepsilon)$ is the number of particles per lattice cell of the clean band, and the displacement

[TABLE]

In this approximation, the normalized ground state is a coherent state

[TABLE]

parametrized by the displacement (31). While this wavefunction is similar in form to that of the linearly dispersive model [20, 10], the average oscillator displacements $f_{1}(\varepsilon)$ now clearly encode information about the whole band structure. In particular, from (31) follows that $f_{1}(\varepsilon)$ vanishes for $\varepsilon>2D$ and that, owing to the van Hove singularities in $N(\varepsilon)$ at $\varepsilon=\pm D$, $f_{1}(\varepsilon)$ has cusps at $\varepsilon=D\pm\varepsilon_{F}$. As we will see below, these are properties shared with the exact solution of the problem.

Within this first order approximation, the rescaled bosonic displacement (29) is simply given by $\Phi(\varepsilon)\simeq f_{1}(\varepsilon)$. Taking the limit $\varepsilon\to 0$ in (31), we find

[TABLE]

As already mentioned, the exact answer is $\Phi(0)=\phi(\varepsilon_{F})/\pi$, with the phase shift given by (27). We see that the approximate answer (33) is correct to first order in $V$, which is expected due to our first order expansion of the potential in the bosonic fields.

The advantage in our method becomes more evident when investigating the full profile of the bosonic displacement $\Phi(\varepsilon)$ as a function of $\varepsilon$. In Figure 1 we compare approximate and exact results for $\Phi(\varepsilon)/\Phi(0)$ in the case of a cosine dispersion relation (namely for a nearest neighbor tight-binding band), for which

[TABLE]

We consider here half-filling, i.e. $\varepsilon_{F}=0$. The exact results were obtained by exact diagonalization of the fermionic Hamiltonian, for respectively $V=D/4,\,D$ and $4D$. For $V<D/4$, the exact result is nearly indistinguishable from the approximate result, for all energies $\varepsilon$. We see that even when $V$ is so large that the approximate result is no longer quantitatively accurate, there are still strong qualitative similarities between the exact and approximate result.

Now we discuss the calculation of the transition rate $W(\varepsilon)$ for the approximate ground state $\left|f_{1}\right>$ given in Eq. (32). According to Eq. (23), we obtain $W(\varepsilon)$ from the Fourier transform of the function

[TABLE]

In the time representation, the above expression reads

[TABLE]

We have encountered the operator $\psi(\tau_{1})^{\dagger}\psi(\tau_{2})$ before in the potential energy term of $H_{V}$, and we use its bosonized form. After some algebra, we obtain

[TABLE]

The expectation value is evaluated using standard coherent state technology. The result is an analytical expression for $P(t)$ in terms of $f_{1}(\varepsilon)$. It reads

[TABLE]

where

[TABLE]

The above integrals can be implemented as fast Fourier transforms so that the computation time for $P(t)$ scales like $\Omega\ln\Omega$ where $\Omega$ is the size of the energy grid used to discretize $f_{1}(\varepsilon)$.

Before addressing our result for a realistic band structure, we show that our approach reproduces the standard bosonization results [20, 10] for a linearized spectrum. For this purpose, we replace the microscopic density of states with an effective one in which modes far from the Fermi energy are suppressed by a soft cut-off:

[TABLE]

with $\Lambda\ll D$. The effective density of states should be thought of as the result of integrating out high energy modes, and is therefore accompanied by a renormalization of the system parameters $V$ and $\gamma$. This leads to a displacement (31), $f_{1}(\varepsilon)=2N(\varepsilon_{F})Ve^{-\varepsilon/2\Lambda}$. For $t\gg 1/\Lambda$, this gives

[TABLE]

For $|\varepsilon|\ll\Lambda$, the transition rate $W(\varepsilon)$ is proportional to the Fourier transform of (41). By making a change of integration variable $t\to\varepsilon t$ in the Fourier transform and noting the analiticity in the lower half of the complex $t$-plane of (41), one then readily obtains

[TABLE]

with some cut-off dependent constant $W_{\Lambda}$. As discussed previously, the microscopic power-law behavior at the threshold can be predicted by standard bosonization provided the exact phase shift is used in place of linearized expression $\pi f_{1}(0)=2\pi VN(\varepsilon_{F})$. On the other hand, it is usually a daunting task to relate the prefactor $W_{\Lambda}$ to the bare system parameters within the standard bosonization method, as this requires explicitly working out how the system parameters are renormalized when ultraviolet modes are eliminated. The microscopic bosonization method can however get access to this prefactor, and in fact to the whole energy dependence of the transition rate, as we demonstrate now.

For the cosine band (34), we have calculated the x-ray transiton rate $W(\varepsilon)$ in the first order approximation (38) of the bosonized theory, and compared to exact results obtained by direct diagonalization of the original fermion Hamiltonian. (See Fig. 2.) We find excellent quantitative agreement for $|V|=\pm 0.16D$ on all energy scales, not only near the threshold, but also within and outside the band. The van Hove singularity at the band bottom is also taken into account properly by this simple analytical treatment. Because the first order expansion only captures the phase shift correctly for $V\lesssim D/4$, the method cannot be trusted regarding the x-ray emission spectrum for larger $V$ values. A more accurate variational approach for, valid even for $V>D$ will now be developed in Sec. 5, and we defer to this section a thorough discussion of the x-ray edge spectra.

5 Variational bosonic solution of the microscopic x-ray edge problem

5.1 Single coherent state variational theory

For moderate interaction strengths $V<D/4$, the analytical approximation that we presented in the previous section was found to capture the physics of the x-ray problem quantitatively on all energy scales. For larger interaction strengths, the approximation remained only qualitatively predictive, as seen in Fig. 1. This is already a remarkable success for a microscopic bosonization approach of electronic lattice models. In the first part of this section we investigate a straightforward variational generalization of the analytical theory that builds on the natural structure of the wave function in terms of bosonic coherent states. We will see that the variational Ansatz presented in this section captures important aspects of the x-ray edge physics at large $V$, especially beyond the band edge, but still leaves room for improvement. In the second part of this section, we will formulate an improved variational Ansatz which produces a rate $W(\varepsilon)$ that, though not exact, is very accurate up to phase shifts nearly equal to the maximal value $\pi/2$. This will serve to illustrate that a fully microscopic bosonization calculation can be performed on tight-binding models even beyond the weak interaction regime, an aspect that clearly clashes with the common wisdom about bosonization.

Our first variational Ansatz is simply a generic coherent state

[TABLE]

which is obtained by promoting $f_{1}(\varepsilon)$ in the approximate ground state (32) to a function $f(\varepsilon)$ that is determined by minimizing the energy $E_{\rm var}=\left<f\right|H_{V,{\rm enl}}\left|f\right>$. Although the computations are performed entirely within the bosonized theory, we provide in Appendix: Physical meaning of the coherent state wave function a physical interpretation of this coherent state in terms of the original fermions. This coherent state Ansatz is guaranteed to produce a Fermi edge singularity of the form (26), with a phase shift $\phi(\varepsilon_{F})=\pi f(0)$. Our main goal in using this variational state is to account for the non-linear behavior of the phase shift as a function of $V$, Eq. (27), which lies beyond the leading-order expression (33). One can verify that $\pi f(0)$ corresponds to the phase shift at the Fermi level, by studying the overlap with the vacuum state $|{\rm vac}\rangle$:

[TABLE]

which vanishes because the integral in the exponent is logarithmically divergent at small $\varepsilon$. In a finite system, the divergence is cut off at an energy $\sim v_{F}/L$, where $L$ is the system size. The overlap $\left<\rm vac\right.\left|f\right>$ thus vanishes like $L^{-f(0)^{2}/2}$. In view of Anderson’s orthogonality theorem [23], we again identify $\pi f(0)$ as the phase shift at the Fermi energy.

For the single coherent state Ansatz, the energy functional takes the following explicit form:

[TABLE]

where

[TABLE]

is real. The functional derivative with respect to $f(\varepsilon)$ is

[TABLE]

where

[TABLE]

The final lines of (47) and (49) are convenient for numerical calculations, because all integral transforms that are involved can be implemented as fast Fourier transforms. The computation time for calculating $E_{\rm var}$ and its gradient again scales like $\Omega\ln\Omega$, where $\Omega$ is the size of the energy grid used to discretize $f(\varepsilon)$. Note that (47) can be solved for $f(\varepsilon)$ to linear order in $V$, by setting $f(\varepsilon)=0$ in $A(\omega)$. By doing so, we recover $f(\varepsilon)=f_{1}(\varepsilon)$, where $f_{1}(\varepsilon)$ is given by the first order analytical formula (31). Thus, at small $V$, the single coherent state variational Ansatz reduces to our previous approximation where the potential is expanded to first order in boson operators. Once the variational state is optimized, we need to calculate the correlation function $P(t)$ in Eq. (23), and from there, the transition rate $W(\varepsilon)$. Since the variational state has the same coherent state form as the analytical approximation of the previous section, we can do so simply by replacing $f_{1}(\varepsilon)$ by $f(\varepsilon)$ in the equations (38) and (39) of the previous section. The details of our numerical implementation of the variational calculation can be found in the supplementary material. We were comfortably able to perform the variational calculation at an energy resolution of $10^{-3}D$.

The results we present are for a cosine band (34), leading to the explicit form of

[TABLE]

where $\Gamma(z)$ is the Gamma function, and $J_{\nu}(z)$ is the Bessel function of order $\nu$. Because of particle-hole symmetry, i.e. $N(\varepsilon)=N(-\varepsilon)$, $\Delta(\tau)$ is real and even. The Fermi energy is set in the middle of the band throughout, and the phase shift is given exactly by:

[TABLE]

In the left panel of Figure 3, we compare the minimized variational energy $E_{\rm var}$ for various $V$ to the true ground state energy of the infinite system, which is given by

[TABLE]

At the smallest potential strength that we considered, namely $V=0.16D$, the coherent state Ansatz yields an energy that is accurate to within the discretization error $\sim 10^{-3}D$. However, the error $E_{\rm var}-E_{0}^{V}$ grows to a significant fraction of the band width $2D$ when $V$ becomes large, indicating that the Ansatz does not provide a quantitatively accurate description of the ultraviolet modes that are affected by the potential (this problem will be cured by an improved Ansatz in what follows).

In the right panel of Figure 3, we plot the variational parameter $f(0)$ as a function of $V$, which gives $\pi^{-1}$ times the phase shift of fermions at the Fermi energy, for the variational state. We compare it to the exact phase shift (also divided by $\pi$), given by $\arctan(V/D)/\pi$, when the Fermi energy is in the middle of the band. From (31) it follows that to lowest order in $V$, the variational state reproduces the correct phase shift $2\pi N(\varepsilon_{F})V$, as we can check on the plot. At large positive $V$, the single coherent state Ansatz nicely corrects for the unbounded growth found at large $V$ in Schotte and Schotte’s solution, although it significantly underestimates the exact phase shift. The variational Ansatz respects the Friedel sum rule, so that the average number of particles displaced by the potential energy term is $f(0)$. The underestimation of $f(0)$ at large positive $V$ therefore implies that the single coherent state Ansatz displaces too few particles at large $V$.

We now use the optimized coherent state trial wave function to compute the x-ray transition rate $W(\varepsilon)$ and compare the variationally obtained rate to numerically exact diagonalization results. (See Supplementary Material for details about the implementation.) Before discussing the variational results, it is useful to highlight the following features of the numerically exact results. As $\varepsilon$ approaches the threshold $\varepsilon=0$ from below, $W(\varepsilon)$ displays power-law behavior (the Fermi edge singularity) with an exponent in agreement with the analytical prediction $[\phi(0)/\pi+1]^{2}-1$, where $\phi(0)$ is given by (51). In the $V=0$ limit, we know that $W(\varepsilon)$ diverges as $1/\sqrt{1+\varepsilon/D}$ when $\varepsilon$ approaches $-D$ from above, due to the van Hove singularity in $N(\varepsilon)$. When $V$ is varied, this peak does not seem to significantly broaden in our exact diagonalization data, and we conclude that the van Hove divergence remains present when $V\not=0$. For $V=0$, the transition rate is strictly zero below the band bottom for $\varepsilon<-D$. However, the rate develops a tail in the region $\varepsilon<-D$ for non-zero $V$. Since $\varepsilon<-D$ corresponds to more excitation energy than a single particle-hole pair can carry, the tail is necessarily associated with multiple-particle excitations. It seems from our numerically exact data that $W(\varepsilon)$ tends to a finite value as $\varepsilon$ approaches $-D$ from below.

For $|V|\lesssim D/4$, the variational results are nearly identical to to the analytical calculation of Sec. 4, and thus are quantitatively accurate on all energy scales. Figure 4 demonstrates more-than-qualitative agreement of the transition rate at moderately large values of the potential $V=\pm D$ for the single coherent state Ansatz (43). At both positive and negative $V$, we find that the variational results capture both the infrared physics of the Fermi edge singularity close to $\varepsilon=0$, and the band structure physics of the van Hove singularity at $\varepsilon=-D$. Remarkably, the tail beyond the band edge at $\varepsilon<-D$ is also reproduced. At large negative $V$, deviations from the exact result become more pronounced, but the general shape and scale of the variational curve is still similar to the exact result. For instance, the deviation $\left|1-W_{\rm var}(\varepsilon)/W_{\rm exact}(\varepsilon)\right|$, averaged over the interval $\varepsilon\in(-D,0)$, peaks at $\sim 15\%$ as $|V|$ is increased, and the largest contribution to the error comes from the vicinity of the singularities at $\varepsilon=0$ and $\varepsilon=-D$. At large positive $V$, we find a more important mismatch to the exact results, with the variationally calculated rate significantly larger than the true one. Curiously, the shape of the variationally determined transition rate remains in excellent agreement with the exact answer. For instance, at $V=D$, scaling the variational rate by $0.55$ produces a result nearly identical to the exact result for all $\varepsilon$ (not shown).

According to the definition of the transition rate (19), one finds the sum rule:

[TABLE]

i.e. the area under the curve of $W(\varepsilon)$ is proportional to the average number of particles at the impurity site. The fact that the single coherent state Ansatz significantly overestimates the transition rate $W(\varepsilon)$ at large positive $V$, implies that it predicts the wrong average number of particles on the site zero of the lattice. Yet, it is quite surprising that the overall line-shape of the emission spectrum is so well described. This useful piece of information will lead to a drastic improvement of the variational Ansatz, that we consider next.

5.2 Improved variational Ansatz as superposed coherent states

We now propose an improved Ansatz, based on the previous considerations. We found that the single coherent state Ansatz (43) produces an x-ray transition rate $W(\varepsilon)$ that is qualitatively correct up to large negative $V$, and that has almost the perfect line shape at large positive $V$ (but not the correct scale overall). Owing to the form of the correlation function (23) involved in the transition rate, the component of the wave function in which unit cell zero is empty will not affect the shape of the transition rate, only its magnitude. For $V>0$, we therefore consider an Ansatz of the form

[TABLE]

Here $\left|f\right>$ and $\left|g\right>$ are single coherent states (32), with real functions $f(\varepsilon)$, $g(\varepsilon)$ and a real relative weight $c$, all to be optimized variationally. Since the trial state consists of two terms with distinct configurations of collective degrees of freedom, we refer to it as the superposed Ansatz. It is interesting to note that the second term in the Ansatz can be written in the form $\int_{-\infty}^{\infty}d\tau\,K(\tau)\left|g_{\tau}\right>$ where $\left|g_{\tau}\right>$ is a coherent state with displacement $g_{\tau}(\varepsilon)=g(\varepsilon)+\exp[-(i\tau+a/2)\varepsilon]$. Of course, an arbitrary state can be written as a superposition of coherent states. However, in the general case, the exact decomposition consists of an infinite-dimensional integral, with two dimensions for every bosonic mode, because each mode displacement can vary independently over the whole complex plane. Here in contrast, we are dealing with a single one-dimensional integral, that limits the the degree of entanglement in the trial state. A similar sparseness was encountered previously in the systematic coherent state expansion that the authors developed to deal with Kondo-type impurity models in an infinite flat band [37, 38]. There, in fact, a handful of terms sufficed to account very accurately for the proper many-body ground state.

For $c=0$, the superposed Ansatz reduces to the single coherent state trial wave function of the previous section, but additional control over the occupation of the impurity site is provided by the second term, where the operator $U^{\dagger}\tilde{\psi}_{0}$ depletes unit cell zero of particles, without removing any particles from the band as a whole. For negative $V$, $U^{\dagger}\tilde{\psi}_{0}$ in the second term of (54) is replaced by $U\tilde{\psi}_{0}^{\dagger}$. It is also convenient to reverse the signs of $f(\varepsilon)$ and $g(\varepsilon)$, so that the negative $V$ Ansatz reads

[TABLE]

Under particle-hole interchange, $b_{\varepsilon}$ maps onto $-b_{\varepsilon}$. If the band possesses particle-hole symmetry, i.e. if $N(\varepsilon)=N(-\varepsilon)$, the same optimal values of $c$, $f(\varepsilon)$ and $g(\varepsilon)$ then minimize the energy at $(\varepsilon_{F},V)$ and $(-\varepsilon_{F},-V)$.

Using the same arguments as before, one readily deduces that the second term in the superposed Ansatz (54) is associated with a phase shift $\phi(\varepsilon_{F})=\pi[g(0)+1]$. When $g(0)\not=f(0)-1$, all cross-terms in the expectation value of the Hamiltonian vanish, owing to the Anderson orthogonality theorem. As a result, one must have

[TABLE]

in the optimized Ansatz for the energy actually to be lowered. Note that unlike the single coherent state Ansatz, the superposed Ansatz (54) has to be normalized by hand. Expressions for the associated energy functional and its gradient, in terms of the variational parameters, can be found in the supplementary material. Although lengthier than the expressions for the single coherent state Ansatz of Sec. 5, they hold the same advantage over the exact diagonalization approach, namely that they can be implemented using a sequence of fast Fourier transforms, and the scaling of execution time with energy grid size is again linear (up to logarithmic corrections), in contrast to the cubic cost for the exact diagonalization method.

For the superposed Ansatz (54), it is easier to calculate $P(t)$ and hence the transition rate $W(\varepsilon)$ for positive $V$ than for negative $V$. For positive $V$, the second term on the right hand side of (54) only contributes to $P(t)$ through a normalization constant, because $U^{\dagger}$ commutes with $\tilde{\psi}_{0}$ and $\tilde{\psi}_{0}^{2}=0$, so that

[TABLE]

The numerator equals the right hand side of the expression (38) that we derived for $P(t)$ in the context of the analytical approximation of Sec. 4.2. The denominator, which we also have to calculate when we minimize $E_{\rm var}$, can be constructed from the expressions presented in the supplementary material. In order to calculate $P(t)$ at negative $V$ on the other hand, we must evaluate the expectation value

[TABLE]

which contains non-vanishing cross-terms between the first and second terms of the Ansatz, because compared to the $V>0$ case, products like $\tilde{\psi}_{0}^{2}=0$ are replaced by $\tilde{\psi}_{0}^{\dagger}\tilde{\psi}_{0}\not=0$. We have not been able to write the resulting expressions in such a way that integral transforms can be implemented as fast Fourier transforms. For this reason, the results for $W(\varepsilon)$ that we present in the next section are restricted to positive $V$.

We discuss now our results for a cosine dispersion relation with the Fermi energy in the middle of the band. In the previous Figure 3 we compare both the minimized energy $E_{\rm var}$ and the phase shift $f(0)=\phi(\varepsilon_{F}=0)/\pi$ of the superposed Ansatz to the exact results for an infinite system and to the single coherent state Ansatz. For the superposed Ansatz, in contrast to the single coherent state Ansatz, the error $E_{\rm var}-E_{0}^{V}$ saturates to $\sim 1\%$ of $D$ at large $V$, suggesting that ultraviolet modes are now accurately accounted for. For the rescaled phase shift $f(0)$, the superposed Ansatz (54) also yields a significant improvement over the single coherent state solution. In view of this success of the superposed Ansatz (54), it is worth studying the lowest energy configuration of the bosonic degrees of freedom further. To this end, we consider again the rescaled bosonic displacement (29), defined in Sec. 2. For the single coherent state Ansatz, we obviously have $\Phi(\varepsilon)=f(\varepsilon)$. For the superposed Ansatz, the relationship between $\Phi(\varepsilon)$ and the variational parameters is more complicated, and reads for a particle-hole symmetric band at half-filling:

[TABLE]

where

[TABLE]

The function $B$ seems to depend on an as yet undefined parameter $\Lambda$ with dimensions of energy. However, this is not the case. The $\Lambda$-dependence of the pre-factor is exactly cancelled by the $\Lambda$-dependence of the argument of the exponential. Note that for the superposed Ansatz, it still holds that $\Phi(0)=f(0)$, i.e. $\Phi(0)$ corresponds to the average number of fermions displaced by the impurity potential.

In Figure 5, we plot the exact displacement $\Phi(\varepsilon)$ together with its estimated value according to the two variational Ansätze, for two large positive $V$ values (for negative $V$, the vertical axis is simply inverted). The displacement $\Phi(\varepsilon)$ shows several interesting and general features. At $\varepsilon=0$ it equals $\phi/\pi$, the average number of particles displaced by the potential. In addition, $\Phi(\varepsilon)$ has cusps at $\varepsilon=D$ and at $\varepsilon=2D$ which are related to the Friedel oscillations of the fermion density — despite dispersive effects, the Fourier transform of $\Phi(\varepsilon)$ roughly corresponds to the average particle density profile. The cusp value $\Phi(D)$ reaches a maximum at $V\sim 1.4D$ before decreasing again. This mirrors the amplitude of Friedel oscillations: Obviously, the amplitude is zero at $V=0$. The amplitude is also zero at $V\to\infty$, where the impurity cuts the crystal into two uncoupled semi-infinite sections. In between these two limits, the amplitude first increases and then decreases. Finally, the displacement $\Phi(\varepsilon)$ strictly vanishes for $\varepsilon>2D$, because it is associated to particle-hole excitations within a finite band. At larger $V$, the single coherent state Ansatz underestimates $\Phi(0)$, and overestimates $\Phi(D)$, while it is not well suppressed at $\varepsilon>2D$. Indeed the single coherent state does not allow enough freedom to entirely prevent the transfer of particles between band and spectator orbitals. In contrast, the spurious tail for $\varepsilon>2D$ is very small in the superposed Ansatz, because the spectator excitations above the vacuum are energetically expensive and the superposed Ansatz allows for enough freedom to adjust ultraviolet modes optimally.

Globally, the superposed Ansatz agrees very well with the exact result for all $V$ and at all energies. The largest disagreement is found at small $\varepsilon$ where, as we have seen in Figure 3, the superposed Ansatz slightly overestimates the total displaced charge. At very large $V$, the superposed Ansatz does underestimate the strength of the kink in $\Phi(\varepsilon)$ at $\varepsilon=D$. Nonetheless, the non-monotone behaviour of $\Phi(D)$ as a function of $V$ is well reproduced. Apart from elucidating how the superposed Ansatz (54) improves on the single coherent trial state, the above analysis also demonstrates that the bosonic description is less of a black box than the exact diagonalization method. It is possible to form an intuitive understanding of physical properties of the ground state, not only in the infrared, but at all scales, by inspecting the configuration of bosonic degrees of freedom.

A final result, and a highlight of our study, is the transition rate $W(\varepsilon)$ calculated at $V>0$ with the superposed Ansatz. As can be seen in Figure 6, there is excellent agreement with the exact result, up to $V$ as large as $6.3D$, which corresponds to a phase shift $\phi(0)=0.9\,\pi/2$. Not only do we capture the Fermi-edge singularity for $\varepsilon\to 0$, but multiple-particle excitations corresponding to $\varepsilon<-D$ are also well-accounted for. This implies that our variational treatment accurately describes all excitations that are produced when the Fermi sea is shaken up during an inelastic tunnelling transition, even at strong impurity interactions and large energy transfers.

6 Conclusion and perspectives

While bosonization has proven invaluable in the study of impurity and bulk one-dimensional systems [10], its practical applications has until now largely been confined to the study of linearly dispersing fermion models. More realistic calculations in the fermionic language are typically based on an infinite size extension [39, 40] of the density matrix renormalization group [4, 2], leading for instance to detailed studies of quantum spin chains [41]. Our work provides a proof of principle that the bosonization technique can be successfully applied to lattice models, incorporating a non-trivial band structure into a description valid on all energy scales.

As a testing ground for our ideas, we studied the x-ray edge singularity for a microscopic model of a one-dimensional electronic band, investigating the emission rate on all energy scales. We found, in addition to the known Fermi-edge singularity, a strong reorganization of many-electron states even at energy scales below the band edge, that can be accounted for simply in a bosonic language. At weak coupling, we found an analytic expression for the x-ray spectrum, that is quantitatively correct at all scales, and that, to our best knowledge, has not appeared in the literature before. At strong coupling we provide a non-perturbative variational solution, that quantitatively accounts for x-ray processes even at frequencies outside the conduction band, where the signal is exclusively produced by multi-electron processes. Besides building intuition by efficiently parameterizing Slater determinants in terms of bosonic coherent states, our scheme turned out to be numerically very efficient, with a scaling $\Omega\ln\Omega$ of the computation time, where $\Omega$ is the size of the energy grid used to discretize the system, instead of $\Omega^{3}$ for the exact diagonalization. Though not exact, the variational approach yields very accurate results (see our main result in Figure 6).

On a conceptual level, this leads to a drastic change of viewpoint in the many-body problem, in which a technique usually associated with effective low-energy theories, is used to predict microscopically the behavior of excitations across the whole band structure for lattice electronic models with arbitrary band dispersion. In future work we plan to exploit this to study bulk-interacting fermions on a lattice. Another avenue for future work is to extend the coherent state variational technique employed here to deal with time-dependent problems, as has already been done in the context of the spin-boson model [42].

Our analysis was based on trial wave functions built on coherent states. In future applications of our bosonization approach, the many-boson systems that results could in principle be studied by other numerical means as well (quantum Monte Carlo [33, 34, 35, 36] or the Numerical Renormalization Group [32]). An interesting open question concerns whether the microscopic bosonic description also holds advantages over the original fermionic one, as it did in the case we studied, when applied in conjunction with these techniques.

Appendix: Physical meaning of the coherent state wave function

We address here the meaning of the bosonic coherent states in terms of the original fermions for the problem of the dispersive band, in first quantization language. The bosonic coherent state $\left|f\right>$ (43) is the exact ground state of the bosonic parent Hamiltonian

[TABLE]

Using the relation (7) between the bosonic $b_{\varepsilon}$ operators and the fermionic $\psi(\tau)$ operators, the parent Hamiltonian can be refermionized

[TABLE]

where $\bar{f}(\tau)=\int_{-\infty}^{\infty}d\omega\,e^{i\omega t}f(|\omega|)$. The exact x-ray edge Hamiltonian can be expressed in the $\tau$ basis as:

[TABLE]

Thus we see that approximating the ground state of the true Hamiltonian (64) with a coherent state, amounts to replacing the non-local in $\tau$ scattering term, with an effective potential $\int_{-\infty}^{\infty}d\tau\,\bar{f}(\tau)\psi^{\dagger}(\tau)\psi(\tau)$ that is local in $\tau$. If we denote the single-particle eigenstate with energy $\varepsilon$ of this non-interacting parent Hamiltonian by $\left|\psi_{\varepsilon}\right>$, then the time-representation single-particle wave functions read

[TABLE]

where we may choose

[TABLE]

Let us return to the energy representation and denote the single-particle eigenstates of the clean system as $c_{\varepsilon}^{\dagger}\left|0\right>=\left|\varepsilon\right>$. In the energy representation the single-particle wave functions of the parent Hamiltonian read

[TABLE]

Here $e^{-iF}$ is an operator that acts on the single-particle Hilbert space such that $\left|\psi_{\varepsilon}\right>=e^{-iF}\left|\varepsilon\right>$. The operator $F$ has matrix elements

[TABLE]

A Slater determinant, constructed from the single-particle wave functions $\left<\varepsilon^{\prime}\right|\left.\psi_{\varepsilon}\right>$ with $\varepsilon<\varepsilon_{F}$, is fully equivalent to the single bosonic coherent state (43). Up to moderate interaction strength $V$, there is little mixing of physical and spectator degrees of freedom in the coherent state, implying that $\left<\varepsilon^{\prime}\right|\left.\psi_{\varepsilon}\right>$ is close to a delta-function for $\varepsilon^{\prime}>|D|$. In this regime a Slater determinant restricted to band orbitals only, and excluding spectator orbitals, is thus nearly equivalent to a bosonic coherent state.

This work is based on research supported in part by the National Research Foundation of South Africa (Grant Number 90657) and CNRS PICS contract “FermiCats”.