Bosonization of the $\mathbf{Q}=0$ continuum of Dirac Fermions

We develop a bosonization formalism that captures non-perturbatively the interaction effects on the $\mathbf{Q}=0$ continuum of excitations of nodal fermions above one dimension. Our approach is a natural extension of the classic bosonization scheme for higher dimensional Fermi surfaces to include the $\mathbf{Q}=0$ neutral excitations that would be absent in a single-band system. The problem is reduced to solving a boson bilinear Hamiltonian. We establish a rigorous microscopic footing for this approach by showing that the solution of such boson bilinear Hamiltonian is exactly equivalent to performing the infinite sum of Feynman diagrams associated with the Kadanoff-Baym particle-hole propagator that arises from the self-consistent Hartree-Fock approximation to the single particle Green's function. We apply this machinery to compute the interaction corrections to the optical conductivity of 2D Dirac Fermions with Coulomb interactions reproducing the results of perturbative renormalization group at weak coupling and extending them to the strong coupling regime.


Introduction.
The remarkable success of bosonization in capturing the non-perturbative properties of interacting fermions in one-dimension [1] has long motivated the quest for extensions of this program to higher dimensions. One major such enterprise has been the development of higher dimensional bosonization of Fermi surfaces [2][3][4][5][6]. In this approach, particle-hole creation operators of a given total momentum Q, c † k+Q/2 c † k−Q/2 , are promoted to bosonic creation operators with a commutator that is approximated as a number. The resulting bosonized Hamiltonian only couples bosonic modes with momentum Q to bosons with either +Q or −Q. Namely, there is zero amplitude for a particle hole-pair with momentum Q to transition into two particle-hole pairs with momenta Q 1,2 and Q = Q 1 + Q 2 . The only allowed process are for the particle-hole pair with momentum Q to scatter into another one with the same Q, or to create pairs of particle-hole pairs with momentum +Q and −Q (for a succinct incarnation of this structure see, e.g., Eq.(7.1) in Ref. [7]). This assumption of separability of Hilbert spaces of particle hole pairs with different magnitudes of |Q|, lies at the heart of the higher dimensional bosonization approach to Fermi surfaces and it is believed to be an asymptotically correct description of particle-hole excitations of Landau fermi liquids at small |Q|.
Ordinary single-band Fermi liquids do not have low energy particle-hole excitations with total momentum Q = 0 and therefore this sector does not appear in the conventional problem of bosonization of Fermi surfaces. In contrast, nodal semimetals, in which the Fermi surface shrinks to a point, such as Weyl or massless Dirac semimetals, have a non-trivial set of gapless optical Q = 0 particle-hole excitations. The central purpose of the present study is to develop a systematic bosonization approach to this sector for gapless semimetals. For concreteness we will discuss only 2D massless Dirac fermions, such as those appearing in graphene and the surface of 3D topological insulators, but our ideas can be naturally extended to other cases and higher dimensions. To describe such excitations, we will borrow the central assumption of the bosonization approach of Fermi surfaces, ω + ν0 ω + ν1 ω + ν2 · · · ω + kn νn · · · ν2 ν1 ν0 k0s0 k1s1 k2s2 · · · knsn knsn · · · k2s2 k1s1 k0s0 (kx, ky) namely, that such optical particle-hole pairs are decoupled from the particle-hole pairs of finite momentum Q. We expect this simplification to be justified at low energies in phases which are adiabatically related to free fermions, in a similar sense to how such decoupling allows to describe Fermi liquids which are adiabatically related to free fermions in the higher dimensional bosonization of Fermi surfaces. We will, however, establish, a very explicit and solid connection between our bosonization approach and the conventional Feynman diagrammatic perturbation theory that demonstrates the validity of this central assumption of our approach. Specifically, we will prove that the solution of our effective bosonic arXiv:2002.05732v1 [cond-mat.str-el] 13 Feb 2020 Hamiltonian for the optical particle-hole pairs is exactly equivalent to the self-consistent Kadanoff-Baym resummation [8,9] of the particle-hole propagator at Q = 0, associated with the self-consistent Hartree-Fock approximation to the single particle-particle Green's function.
As an application of our approach we will compute the interaction corrections to the optical conductivity of 2D Dirac fermions with Coulomb interactions, whose strength is parametrized by the effective fine structure constant α = e 2 / v, where v is the velocity of the Dirac fermions and the dielectric constant of the surrounding medium. This optical conductivity at low energies is determined by fundamental constants of nature, and given by σ 0 = e 2 /16 per Dirac cone [10,11]. Its zero frequency limit is not expected to be renormalized by interactions, but, Coulomb interactions can produce a slow flow as a function of frequency to such value and a non-trivial non-analytic frequency dependence at low energies. Early perturbative calculations of such corrections where in mutual disagreement [12,13], but subsequent studies [14][15][16][17][18][19][20] validated the result of Ref. [12]. As we will see, our approach will recover the perturbative results of Ref. [12] at small interactions and extend them non-perturbatively to finite α. Non-perturbative attempts to understand the effects of Coulomb interactions in the optical conductivity of Dirac fermions have been scarce. A Quantum Monte Carlo effort [21] to compute the optical conductivity concluded that interaction corrections remain rather small even at α ∼ 2. Our analysis will also support this conclusion, which is broadly in agreement with experiments that have found values close to that for non-interacting fermions [22][23][24].
Effective Hamiltonian and Hilbert space. The microscopic Hamiltonian is ( = 1): where A is the system area, and V q is the Fourier transform of the interaction potential. It is convenient to imagine the Fermions moving in a 2D Torus so that its momentum is quantized on a lattice. In this momentum lattice the complete many-body Hilbert space is a tensor product of empty, singly and doubly occupied states: The kinetic term in Eq. (1) produces no fluctuations between the occupancy of the momentum sites, and favors a ground state with singly occupied states with a suitably oriented spin in the form of vortex around the Dirac point (see Fig. 1a). The interactions are a form of pair hopping terms that have a finite amplitude to induce transitions into states with doubly occupied sites and empty sites. Crucially, the subspace of the Hilbert space with singly occupied sites is equivalent to the space of particle-hole pairs with zero total momentum, Q = 0, while those states with doubly occupied and empty sites contain particle-hole excitations of finite momentum Q. Therefore, following the spirit of higher dimensional bosonization, we will project the Hamiltonian in Eq. (1) onto the Hilbert space of singly occupied sites in the momentum lattice, depicted in Fig. 1b. This Hilbert space contains a spin-1/2 at each momentum site: and the projection of the Hamiltonian from Eq. (1) leads to the following Heisenberg model: where s k = σ,σ ψ † k,σ σ σσ ψ † k,σ is a spin operator for the k site of the momentum lattice. The first term in Eq. (4) is a Zeeman vortex field and the second term is a long-range exchange coupling.
This Hamiltonian is not exactly solvable but the fluctuations around the non-interacting state can be described by a Holstein-Primakoff expansion [25]. To do so, we choose a spin basis that diagonalizes the kinetic energy at each momentum site s k = −s z kk + s x kẑ + s y kφ whereẑ is the out-of-plane direction andφ =ẑ ×k. The spin operators can be expanded as s z Up to boson bilinears the Hamiltonian becomes (see §A 2 of [26]): with E k = v|k| + Σ k and Σ k = k V k−k cos φ kk /2A is the Hartree-Fock self-energy and T kk is the interaction matrix in the band basis (for details see §A 1 of [26]) Connection to perturbation theory.
We will now demonstrate that the solution of the boson bilinear Hamiltonian in Eq. (5) is exactly equivalent to the calculation of the particle-hole propagator within the Kadanoff-Baym (KB) resummation of Feynman diagrams associated with the self-consistent Hartree-Fock (SCHF) approximation to the single-particle Green's function. In terms of electrons, the boson creation operator, b † k corresponds to the interband Q = 0 electronhole pair creation operator where the subindex s = ± denotes the conduction and valence bands. Therefore, our goal is to compute the electronhole pair propagator defined as: including all the terms of the KB SCHF resummation, which includes the entire Bethe-Salpeter ladder series with all internal single-particle Green's functions dressed with the SCHF self-energy [8,9]. These SCHF Green's functions are given by (see §A 3 in [26]): . (9) and the n-th order Feynman diagram of this series is shown in Fig. 1c. The zeroth order term of the series is: . (10) An important property of this series, that can be readily obtained by integrating over internal intermediate frequencies, is that the intermediate Green's functions are all constrained to satisfy s = s which physically means that the intermediate pairs always have one electron in the conduction and the other in the valence band. This allows to cast the series as a matrix geometric series involving χ 0 and T of eq. (7) of the form: and therefore the solution of the series has the form: where H k0k f is given in Eq. (6) and τ z is the diagonal Pauli matrix. The structure of this correlator is equal to the propagator of the HP bosons of the Hamiltonian (6).
We therefore see that the exciton propagator has an identical effective Hamiltonian to the one obtained from the HP bosonic Hamiltonian in Eq. (6), demonstrating that the bosonized Hamiltonian is equivalent to self-consistent KB resummation of the particle-hole propagator. Momentum space reparametrization. So far we have imagined our system to have a finite size so that momenta belongs to a discrete lattice. However, it is convenient to perform a reparametrization that manifestly displays the symmetries of the thermodynamic limit. If we parametrize momentum space by a new coordinate z(k), we can trade our boson Hamiltonian by one in a different lattice given by: As detailed in [26], in order to preserve the underlying microscopic normalization of the states, the boson operators and the Hamiltonian in the new lattice need to be rescaled as follows: where J(z) = D(z)(∆z 1 ∆z 2 )/(∆k 1 ∆k 2 ), ∆k i = 2π/L i , ∆z i is the discretization unit of the new coordinate system and D(z) is the Jacobian of the transformation. In particular, in order to exploit the emergent rotational invariance in the thermodynamic limit, we use the following polar parametrization z = (k, φ): where (k m , φ n ) are the polar coordinates of a given site in the polar momentum lattice depicted in Figs. 1d and 1e, K is the UV momentum scale, ∆θ = ( π /2)/(M + 1), ∆φ = 2π/(2L + 1), and n = 0, ..., 2L, m = 1, ..., M .
The radial discretization we are choosing is denser at small k and more dilute at large k. This is not crucial but allows faster numerical convergence at low energies. Notice also that we do not have a hard cutoff but the largest momentum apporaches infinity as M → ∞. We have verified that the results we will describe are independent of the specific choice of the radial discretization once the grids become sufficiently dense [26].
Applying the transformation from Eq. (13) to the boson Hamiltonian from Eq. (12) leads to the following decoupling into angular momentum channels: Therefore the problem reduces to a set of bosons moving in an effective one dimensional radial space for each angular momentum channel which in general needs to be solved numerically.
Optical Conductivity. As a concrete application of our formalism we study the Coulomb interaction corrections to the optical conductivity of Dirac fermions. We follow the Kubo approach to compute the conductivity from the current-current correlator χ µν (t) = iΘ(t)A [j µ (t), j ν (0)] . The total current operator carries Q = 0, so it can be represented exactly within the effective spin-1/2 Hilbert space of Eq. (4) as follows: Using the HP approximation for the spin operators, the current-current correlator then can be expressed as (see Eq. (C8) in [26])  Then if the Hamiltonian for the = 1 angular momentum channel is diagonalized by a transformation of the form (see also [27]): where Ω nn = diag ω 0 · · · ω M −ω 0 · · ·− ω M is the diagonal matrix of the eigenvalues of Eq. (5), the real part of the conductivity can be obtained from the following Lehmann-type representation (see Eq. (C12) in [26]): We will now describe the results for the optical conductivity obtained by numerically diagonalizing the = 1 angular momentum bilinear Hamiltonian of Eq. (15) for the Coulomb interaction V q = 2πe 2 / |q|. Further details on the numerics can be found in [26]. To isolate the interaction corrections to σ(ω) we define: where σ 0 = e 2 /16 is the non-interacting conductivity of Dirac fermions. The leading perturbative correction to this conductivity is expected to be of the form [12]: σ = Cα + O(α 2 ), with C = (19 − 6π)/12. We have been able to reproduce this perturbative correction numerically at small α as shown by solid horizontal lines in Fig.  2a along with the full numerical result from Eq. (18).
At larger values of α, clear deviations from the leading perturbative result are seen in Fig. 2b. One of the conspicuous deviations is a logarithmic decrease of the conductivivity at low frequencies (see [26]). This logarithmic decrease can be explained by the logarithmic running of the coupling constant at small frequencies expected from the perturbative renormalization group (RG) analysis: The predicted RG logarithmic correction is shown by a dotted line in Fig. 2a which is in good agreement with the numerical implementation of Eq. (20) at small α.
For larger values of α we see clear deviations from this leading RG perturbative result, as shown in Fig. 2b. Nevertheless, as shown in Fig. 2, even for a value of α as large as α = 5 the maximal deviation of the conductivity from the non-interacting value is only about 4%. This indicates a resilience of conductivity of Dirac fermions to interactions corrections even when non-pertubative effects are included, in agreement with experiments that have obtained values close to those of non-interacting fermions [22][23][24]. We even suspect that the interaction corrections in a full exact solution of the Hamiltonian in Eq. (1) would even be weaker than the corrections we obtained, because the RPA screening of the Coulomb interactions, roughly speaking, should lead to a reduction of the effective value of α → α RPA ≈ α/(1 + πN α/8), where N is the total number of Dirac cones (e.g. N = 4 for graphene [16]).
Discussion and Summary. We have developed a formalism that captures non-perturbatively the effects interactions on the continuum of Q = 0 particle-hole excitations of Dirac fermions. Our approach is constructed by projecting the full microscopic many-body Hamiltonian of Dirac fermions into the subspace of singly occupied momentum states, leading to an effective spin-1/2 Heisenberg-like model in a momentum lattice. This problem is subsequently reduced to a boson bilinear Hamiltonian by a standard Holstein-Primakoff transformation. We have provided a solid microscopic justification for this formalism by showing that it is equivalent to the Kadanoff-Baym resummation of the particle-hole propagator associated with the SCHF approximation to the single particle Green's function. This approximation is expected to capture the essential universal low energy properties of the semi-metallic phase that evolves adiabatically from Free fermions. We have applied this formalism to compute the Coulomb interaction corrections to the optical conductivity of Dirac fermions and found that it recovers the results of perturbative renormalization group at weak coupling [12] and extended them to strong coupling. Remarkably, we have found that the Coulomb interaction corrections remain very weak (∼ 4%) up to values of the effective fine structure constant α ∼ 5, in agreement with experiments in graphene that have measured a value of the optical conductivity that is consistent with the free electron theory [22][23][24]. Although our discussion has been restricted to 2D Dirac fermions, our approach can be naturally generalized to other multi-band semi-metals and higher dimensions, such as Weyl semimetals [28] and novel nodal fermions [29], providing an interesting tool to capture non-perturbative effects of interactions on the correlation functions of Q = 0 operators of these phases.

SUPPLEMENTAL MATERIAL
Appendix A: Connection to perturbation theory

Band basis and interaction matrix
In this section, we provide details of the derivation of Eqs. (5) to (7) starting from Eq. (1). We begin describing the transfomation from pseudo-spin basis onto band basis. In the band basis s = {+, −} the kinetic term is (e ±iφ =k x ± ik y ): Band and pseudospin basis are related by: Fermion bilinears transform as where where φ i is the polar angle of k i , and φ 12 = φ 1 − φ 2 . The part of the Hamiltonian in Eq. (1) that produces only Q = 0 inter-band transitions in the band basis is Consequently, the Hamiltonian in Eq. (1) of the main text projected onto the subspace of singly occupied momentum sites, defined in Eq. (4), can be expressed as follows: with T s1s2 k1k2 given by: The Hamiltonian of Eq. (A5) can then be expressed in terms of spin operators and leads to the Heisenberg-like model introduced in Eq. (4) of the main text.

Holstein-Primakoff expansion
We select the following spin basis which diagonalizes the kinetic term. On this basis, the Hamiltonian can be expanded in a bosonic representation by means of the Holstein-Primakoff (HP) transformations (S = 1 /2): The term corresponding to the exchange coupling in Eq.
(4) can be transformed into pairing and hopping terms of bosons up to bilinears: The resulting bosonic Hamiltonian after applying the HP transformations is: The first line contains the kinetic and self-energy terms. The second line can be viewed as boson hopping terms in the momentum lattice. The third line can be viewed as pairing terms which change the number of bosons. Lastly, by using the Bogoliubov basis given by the Hamiltonian can be expressed as shown in Eq. (5) and (6) of the main text.

Details of connection to perturbation theory
The propagator of the Dirac fermions without interactions is diagonal in the band basis and is given by: Moreover, the Hartree-Fock self-energy of the fermions is: We can then sum the Dyson series to get the dressed fermionic propagator: where E k = v|k| + Σ k . Equivalently in matrix notation we can write: Our goal is to compute the specific resummation of Feynman diagrams for the particle-hole propagator associated with the Kadanoff-Baym (KB) conserving approximation that results from the self-consistent Hartree-Fock (SCHF) approximation to the single particle Green's function. This resummation consists of the sum of the infinite series of the Bethe-Salpeter ladder for the particlehole propagator with internal Green's functions dressed by the Hartree-Fock self energy from Eq. (A15). The series is depicted in Fig. 3. The particle-hole propagator of interest is defined in Eq. (8). The zeroth order non-interacting term of the series is given by: .
or expressed as a matrix in the band basis . (A17) We will now illustrate the leading terms of the series involving the interaction matrix from Eq. (A6). The diagrams involved at first order in the Bethe-Salpeter ladder are shown in the Fig. 3b, and given by (we assume summation on any repeated index): Similarly the n-th term of the series, shown in Fig. 3c, is given by The full summation can therefore be expressed as a geometric series: χ(ω) = χ 0 (ω) + χ 0 (ω)T χ 0 (ω) + χ 0 (ω)T χ 0 (ω)T χ 0 (ω) + · · · = χ 0 (ω) + χ 0 (ω)T χ 0 (ω) + χ 0 (ω)T χ 0 (ω) + · · · , k1 ω + ν0 ν0 k0s0 k0s0 k1 which correspond to a Dyson-like equation for the dressed particle-hole propagator χ(ω): whose solution is given by: Replacing the results from Eq. (A17) and (A6) we get: or, by using the definition of the HP boson Hamiltonian in Eq. (6) of the main text we get the final expression of the exciton propagator, given by The structure of this correlator is identical to the propagator of the HP bosons of the Hamiltonian (A10). From the above, we can assert that the full resummation of the KB conserving approximation associated with SCHF is equivalent to solving the HP bilinear boson problem. Therefore, in summary, the relation between operators and the Hamiltonian matrix in the new lattice defined by the discretization of the coordinates z(k), with the original operators and Hamiltonian of the square lattice is: The idea is that the Hamiltonian H HP in Eq. (B5) will produce the same physical results as the one in the square lattice in Eq. (5) of the main text in the thermodynamic limit.

Polar re-discretization
We choose z = (k, φ) where k is the radius of the momentum vector and φ its polar angle. We will discretize the radial direction in a non-uniform way, to make it denser at small momenta and more dilute at large momenta. We have checked numerically that the precise form of the discretization is not crucial, but the choice we are making produces faster convergence to the thermodynamic limit. Therefore we choose the radius to be: where K is a UV momentum scale, and θ ∈ (0, π/2) is another parameter labeling the radial coordinate that we will choose to be uniformly discretized. The corresponding Jacobian for this parametrization is: We choose θ and φ to be uniformly discretized as follows: where ∆θ = π /2 M + 1 , ∆φ = 2π 2L + 1 .
After replacing (B7) and (B8) into (B6) we get the expression for B k and H kk in the polar lattice where t(θ m ) = tan(θ m ) sec(θ m ) and k mn = k(θ m , φ n ). Finally, the whole Hamiltonian is Because the Hamiltonian matrix H kk that enters into the Hamiltonina H HP in Eq. (5) of the main text only depends on the difference between the polar angles φ − φ we have conservation of the angular momentum l of the bosons. Consequently, we perform Fourier transforms on the polar angles for the fields B mn and the matrix H mn,m n such that the total Bogoliubov Hamiltonian decomposes into a direct sum for different angular mommentum channels, as follows: The corresponding susceptibility is given by Inversion symmetry guarantees that = 1 contributes the same as = −1, so that (C9)

Representation of the optical conductivity
Let us assume that we diagonalize the Hamiltonian from Eq.(15) for the l=1 angular momentum channel via a Bogoliubov transformation, expressed as follows: where Ω nn = diag ω 0 · · · ω M − ω 0 · · · − ω M is the diagonal matrix of the eigenvalues of the = 1 block of the HP Hamiltonian [27]. Replacing such transformations in Eq. (C9) we get which, because of D 1 † n , D 1 † n = I nn , yields Then, we take the Fourier transform of χ(t) to get the frequency-dependent susceptibility Finally, we take the imaginary part of Eq. (C12) to get the optical conductivity depending on frequency as the following Lehmann-type representation: (C12)

Dependence on the IR cutoff
In Fig. 5 we plot the behavior of the conductivity extrapolated to M → ∞ for different values of the IR cutoff (each panel is for a different fixed value of α). We see in Fig. 5, that at extremely low frequencies the conductivity has a bump followed by drop that changes with the value of the IR cutoff. Therefore this bump and the drop which occur at very low frequencies is a consequence of the IR cutoff which models screening of the Coulomb at long distances as captured by Eq. (D2). Therefore, although this low frequency behavior is not completely unphysical as it models the behavior of the optical conductivity in the presence of a perfect metallic screening gate, it is not part of the universal behavior of the ideal unscreened Coulomb interaction that we are interested in. However, Fig. 5 demonstrates clearly that the conductivities follow a universal curve because they agree perfectly at higher frequencies. Eventually the conductivity escapes from this universal curve at low frequencies when the IR cutoff becomes important. Therefore, we conclude from Fig. 5, that if we see that two curves with different IR cutoff overlap until some low frequency, and below this frequency they start to deviate from each other, the behavior for frequencies above this frequency at which they deviate is universal and represents the behavior for the ideal ideal Coulomb interaction problem without any IR cutoff. The results presented in the main text correspond to frequencies ranges where we observe independence of the IR cutoff, and where we are confident that we are simulating the ideal behavior of the unscreened Coulomb interaction.
3. Non-perturbative effects of large α on the optical conductivity According to the perturbative RG result described in Ref. [12], the optical conductivity is expected to have the following behavior at small frequencies and small α: where C = 19−6π 12 . In the second line of this equation we have expanded the denominator in α to get the leading logarithmic correction to the conductivity. Other logarithmic corrections are expected to contain higher powers of α. We have indeed observed such a weak logarithmic drift of the conductivity with frequency at small values of α, as depicted in Fig. 6. By fitting the conductivity with a logarithmic dependence: we obtained the coefficients which are listed in Table I. As we see there is excellent agreement between the value ofσ 0 at weak coupling and also a reasonable agreement for the value ofσ 1 with those of the perturbative analysis of Ref. [12]. Therefore we have been able capture the logarithmic running of the coupling constant expected from the RG analysis at weak coupling.
FIG. 6. Numerical calculation of the conductivity (color lines) and the the expected value from the leading order perturbative RG (dotted lines). The logarithmic running of the coupling constant leads to a visible linear logarithmic drift of the conductivity at weak coupling.