Spectral properties of heterostructures containing half-metallic ferromagnets in the presence of local many-body correlations

In this work, we investigate models for bulk, bi- and multilayers containing half-metallic ferromagnets (HMFs), at zero and at finite temperature, in order to elucidate the effects of strong electronic correlations on the spectral properties (density of states). Our focus is on the evolution of the finite-temperature many-body induced tails in the half-metallic gap. To this end, the dynamical mean-field theory (DMFT) is employed. For the bulk, a Bethe lattice model is solved using a matrix product states based impurity solver at zero temperature and a continuous-time quantum Monte Carlo (CT-QMC) solver at finite temperature. We demonstrate numerically, in agreement with the analytical result, that the tails vanish at the Fermi level at zero temperature. In order to study multilayers, taken to be square lattices within the layers, we use the real-space DMFT extension with the CT-QMC impurity solver. For bilayers formed by the HMF with a band or correlated insulator, we find that charge fluctuations between the layers enhance the finite temperature tails. In addition, in the presence of inter-layer hopping, a coherent quasiparticle peak forms in the otherwise correlated insulator. In the multilayer heterostructure setup, we find that by suitably choosing the model parameters, the tails at the HMF/Mott insulator interface can be reduced significantly, and that a high spin polarization is conceivable, even in the presence of long-ranged electrostatic interactions.


I. INTRODUCTION
A half-metal is a material that has a metallic density of states at the Fermi level for one spin channel and simultaneously a band gap for the other spin channel.This extreme asymmetry between the spin channels is the source of great promise for spintronic application 1,2 .Half-metallic electrodes could provide fully spin-polarized currents and large magnetoresistance in giant magnetoresistance and tunnel magnetoresistance devices 2 .
Density-functional theory (DFT) [3][4][5][6] studies have identified a number of half-metallic bulk compounds, including Heusler alloys 7 , double perovskites, transition metal oxides, chalcogenides, and pnictides.Some of these proposed materials have been realized in experiments.Within DFT the ground states of these materials are accessible from a single particle picture.However, whenever many-body effects are essential the band theory is expected to fail 8 .In particular, in metallic ferromagnets spin fluctuations play a crucial role 9 .Therefore, the scattering of charge carriers on such magnetic excitations is expected to influence the macroscopic properties of these materials including transport.
Contrary to the itinerant ferromagnets in which states near the Fermi level are quasiparticles for both spin projections, in half-metallic ferromagnets (HMFs) an important role is played by incoherent non-quasiparticle (NQP) states.These occur near the Fermi level in the energy gap [10][11][12][13] , and their tails crosses the Fermi level and produce significant depolarization effects 2 .The density of the NQP states vanishes at the Fermi level but increases strongly on an energy scale of the order of the characteristic magnon frequency.Note the clear distinction between the minority and majority half-metallic cases, corresponding to almost empty and almost full occupation of the insulating spin channel.While for the minority gap HMF the NQP states are located just above the Fermi level, for the majority gap HMF they appear just below the Fermi level.The NQP states are expected to contribute to several physical properties such as polarization, specific heat, and transport 2,14 .In the limit of very strong interactions and close to a completely polarized band, a significant logarithmic singularity exists in the imaginary part of the Green's function, which corresponds to a finite jump in the density of states 15 .Edwards and Herz investigated the stability of the saturated ferromagnetic state using Green's function methods, which agree in the limit of large interactions with variational calculations 16 .
It should be noted that dynamical mean-field theory (DMFT) yields qualitatively similar results in the limit of large U/t 17 : the saturated ferromagnetic state is stable, however for realistic values of U its stability is far from being obvious.(As usual, U and t denote the local interaction and the hopping amplitude, respectively.)DMFT [18][19][20][21] in combination with first principles 22,23 cal-culations have been performed for the prototype HMF, NiMnSb 24 , as well as for other Heusler alloys 25,26 , zincblende structure compounds 27,28 , and CrO 2 29 .Using cluster approaches beyond the local DMFT, the many-body features were found to be enhanced 30,31 .While these effects have been studied for bulk half-metallic systems, less is known for heterostructures containing half-metals except some specific cases of zinc-blende structures 32,33 .
In this paper, we describe the behavior of the manybody induced tails in the half-metallic gap for bulk systems as well as for multilayers using model studies.In our approach we consider bi-and multilayers consisting of a finite number of half-metallic layers in contact with different numbers of metallic or insulating layers.It is expected that away from the interface half-metallicity is preserved on the HMF side.At the interface region parameter optimization is important to preserve half-metallicity.This involves the control of electronic states in the vicinity of the gap to increase the spin-polarization (i.e., reduce the interaction-induced tails) at finite temperatures.In order to produce a half-metal in the interface, a band-gap at the Fermi level either in the minority-or in the majority-spin spectral function (density of states) needs to be created.
Using different single-band Hubbard Hamiltonians on distinct layers allows for optimization of their parameters such as the magnitude of the inter-layer hoppings, strength of local interactions, on-site energies, and Zeeman splittings.Therefore, we study such Hubbard Hamiltonians 34 for the multilayer using the DMFT and its real space extension (R-DMFT) 35,36 .R-DMFT considers a purely local self-energy for the strong electron correlation in the layers.To study electronic charge reconstruction, we extend the Hubbard Hamiltonian to include long-ranged Coulomb repulsion between the layers.We treat the latter on a mean-field level, calculating the electrostatic potential self-consistently from the Poisson equation.For a multilayer of five HMFs and the same number of Mott insulator layers, the effect of the long-ranged repulsion is found to lead to a slight redistribution of charges in the metallic channel.
The focus of our analysis is on a narrow energy range around the Fermi level.We show that by analytic continuation of the self-energy (instead of the Green's function) to real energies, robust numerical results can be obtained.Some preliminary results have been reported recently 37 .We demonstrate that many-body effects (described by DMFT) lead to a dynamical reduction of the Hartree part of the self-energy.Therefore, the splitting between majority and minority spin channels is reduced.Furthermore, a temperature dependent tail emerges in the half-metallic gap, reducing the polarization at high temperatures.The magnitude of these effects can be modified by the optimization of Hamiltonian parameters.We expect that our results will be useful for a systematic engineering of heterostructures containing half-metals with desired properties.
The paper is organized as follows.After the introductory section we provide in Sec.II the computational details, and discuss the relevant parameters and methods used to solve the bulk system as well as the multilayer setup.For completeness, we have included the derivation of the R-DMFT equations in Appendix B. The results section, Sec.III, starts with a discussion of the finite temperature behavior of the half-metallic gap, in particular, of the results for the spectral function and the susceptibility, and compares them with previous calculations.The analytic continuations of the self-energy and the Green's function for the finite temperature spectral functions are compared in Appendix A. This is followed by the results for the bilayer structure, Sec.III B 1, where a square lattice is considered within the layers.In Sec.III B 2 we consider a heterostructure of five half-metallic and five Mott insulator layers.Finally, Sec.IV presents the conclusions of our work.

II. COMPUTATIONAL METHOD
We use a single-band Hubbard model to describe correlation effects in bulk, bi-and multilayer systems.The system Hamiltonian reads: Here ĉ † iσ and ĉiσ are the fermionic creation and annihilation operators at site i with spin σ.We denote the number operator at site i with niσ = ĉ † iσ ĉiσ .The on-site contributions are ˜ iσ = i + σh i − µ, with the on-site energy i , the magnetic splitting h i , and the chemical potential µ.Furthermore, the parameters t ij are the hopping matrix elements, and U i is the Hubbard interaction.The hopping matrix is Hermitian, t ij = t * ji .While there exists a solution for a one-dimensional system 38 , we cannot solve this problem in general in higher dimensions.The dynamical mean-field theory (DMFT) [18][19][20][21] , however, provides a non-perturbative approach which is applicable for any range of parameters, and is exact in the limit of infinite coordination number.It is furthermore exact for both solvable limits, the non-interacting case, U i = 0, and the atomic limit, t ij = 0.
Continuous-time quantum Monte Carlo (CT-QMC) methods are a tool of choice to solve correlated electron problems 39 .In the context of DMFT, the Hubbard model in the limit of infinite coordination number maps onto that of the single-impurity Anderson model (SIAM), which leads to invaluable insight into the Mott transition 20 .Being an action-based method, it allows the simulation of effective low-energy models after integrating out high-energy degrees of freedom.Further applications of the CT-QMC include, e.g., formulations along the Keldysh contour, applications within the cluster extensions of DMFT to include spatial fluctuations 40 , the dual fermion approach 41 , or the dynamical vertex approximation 42 .The CT-QMC methods have different formulations: the interaction expansion (CT-INT) 43 , the auxiliary-field (CT-AUX) 44 , and the hybridization expansion (CT-HYB) 45 .We use the CT-HYB formulation for all finite temperature results presented here, since the efficiency of the segment picture was shown for single-band problems 39 .
CT-QMC operates on "imaginary time", therefore an analytical continuation is necessary to produce spectral functions on the real-frequency axis.This is an ill-conditioned problem and limits the precision of calculated spectral functions.This issue may be especially severe for multiorbital problems with complicated spectral lines.A solution to this problem is to do time evolution directly on the real time axis utilizing matrix product state (MPS) based impurity solvers [46][47][48][49] .These solvers allow for a precise discretization of the hybridization function (with several hundred bath sites per spin) and have been shown to yield excellent results even for sharp peaks at high energies in the spectral function 48 .They have been generalized to multi-orbital impurity solvers [50][51][52] by employing tensornetwork representations of the impurity model, while still keeping very good results at all energies, with moderate computational cost.For the present work, we have extended the MPS solver presented in Refs.49 and 50. to allow for magnetically polarized calculations at zero temperature, T = 0. We calculate the ground state using the density matrix renormalization group (DMRG) 53,54 and we perform the time evolution using time-evolving block decimation (TEBD) 55,56 .

III. RESULTS
Section III A addresses the correlation effects for the bulk setup for a semicircular density of states (DOS), which is realized by the Bethe lattice with infinite coordination number.The results presented here address the mechanism of gap closing as a function of temperature for the bulk HMF.Next, in Sec.III B, we investigate the spectral function of bilayers made of a metal, a bandinsulator, or a Mott insulator attached to the HMF.In particular, we study the changes induced by stacking a larger number of layers.Contrary to the Bethe DOS in the bulk case, for the bi-and multilayers we use the DOS of a square lattice within the layers.

A. Finite temperature behavior of minority spin gap in bulk
In the one-band model Eq.(1) a simple way to generate the half-metallic ferromagnetic state is to introduce a sufficiently strong spin splitting such that one spin subband is empty (or full) in the Hartree-Fock (Stoner) picture.Section III A discusses the results for a homogeneous Hubbard Hamiltonian Eq. (1) of a Bethe lattice with infinite coordination number with half-bandwidth D = 1 eV, spin splitting h = 0.5 eV, on-site energy − µ = 1.5 eV, and on-site interaction U = 2 eV.Difficulties in solving the Hubbard model for such a saturated ferromagnet are well known 10 .
For the real-frequency results at zero temperature, the hybridization of the Bethe lattice ∆ σ (E) = (D/2) 2 G σ (E) was discretized using 251 bath sites per spin.We find the ground-state |GS (T = 0) to be almost fully polarized (n ↓ ∼ 10 −6 ).The interacting Green's function and in turn the spectral function were calculated 48 from the time-evolved ĉ( †) ↓ |GS and ĉ( †) ↑ |GS .This was done using time steps of 0.05 eV −1 , up to a maximal time of t max = 150 eV −1 .A linear prediction 48,57 was performed for the Green's functions during the last 20 DMFT iterations up to a maximal time of t max = 1500 eV −1 , so that no dampening of the time series was required.For the singular-value decompositions, a truncated weight of t w = 10 −10 together with a maximal matrix dimension of 700 was chosen.This maximal dimension was reached during the time-evolution of ĉ( †) ↓ |GS at t = 100 eV −1 .The truncated weight always remained below 10 −8 .
For the QMC results at finite temperature we compute the self-energy via the ratio of the two-particle Green's function F σ and the one-particle Green's function G σ : 58 The brackets • S eff denote the average in the effective impurity model.This provides more accurate results than the Dyson equation, such that the Padé analytic continuation 59,60 of the self-energy is reasonably accurate.We calculate the spectral function from the analytically continued self-energy: where ρ(E) is the one-particle density of states of the non-interacting lattice, and Σ(E) is the Padé analytic continuation of the self-energy, Eq. ( 3). Figure 1 displays the results of the DMFT calculations for zero, low (T = 0.02 eV), and high (T = 0.25 eV) temperature.The dotted line shows the Hartree-Fock (HF) solution as a reference.We first discuss the T ≡ 0 spectrum.As the ↓-spin is completely depleted, the result for the ↑-spin are nearly identical to the HF result.The ↑-spin electrons are almost uncorrelated, the magnitude of the self-energy Σ ↑ (E) is negligibly small.For the ↓-spin we see two main effects of correlations.First, the size of the gap is reduced compared to the HF approximation.For low energies, there is a dynamical reduction of the (static) Hartree self-energy.Additionally, a many-body satellite appears at E ≈ 3.5 eV in A ↓ (E) as shown in Fig. 1c.At low temperature T = 0.02 eV the QMC result for the spectral function Eq. ( 4) is in good agreement with the real-frequency results for zero temperature.There is a deviation for the satellite, however analytic continuation is not expected to resolve features this high in energy well.At high temperature, T = 0.25 eV, we obtain a tail crossing the Fermi level E = 0 shown in Fig. 1b which depolarizes the HMF.Due to the tail the ↓-spin is now partially filled, resulting in correlation effects also in the ↑-spin.The many-body satellite is visible in both spin channels for the high temperature.
Previous calculations 24 used a simplified quantum Monte-Carlo scheme within the so-called exact enumeration technique 20 , therefore results for high temperature (T = 0.25 eV) only were accessible.Our high T results differ from the previous ones 24 which show additional peaks in the spectral function.In contrast to the previous calculations 24 , we determine the spectra from the ana-lytically continued self-energy using Eq. ( 4).In fact, we demonstrate in Appendix A that a Padé analytic continuation of the Green's function-instead of the self-energy Eq. ( 3)-causes the appearance of these spurious features in the spectral function.Figure 2 shows the temperature dependence of the spectral function for the minority spin, A ↓ (E), in particular its tail crossing the Fermi level.The first line in the front is at T = 0.16 eV, subsequent lines show always half the temperature.The disappearance of the spectral weight at the Fermi level with decreasing temperature is apparent.
A specific many-body feature in HMFs is attributed to spin-polaron processes 11 : the spin-down electron excitations forbidden in the one-electron description of HMFs arise due to the superposition of spin-up electron excitations and virtual magnons.In model calculations the existence of this feature has been shown by perturbationtheory arguments for the broad-band case 10 , and in the opposite, infinite-U limit 2,11 .An analytic approximation allows to explore the shape of the temperature dependence of the spectral function for the minority spins considering a contact electron-magnon interaction described by the exchange parameter 2,11,15 .According to this theory, a non-linear temperature dependence is obtained from the competing effects of the magnon contribution to the residue of the Green's function, ∼ T 3/2 , with the shift of the band edge states being proportional to T 5/2 .By a direct fit A ↓ (E) ∝ T α to the data in the inset of Fig. 1c an exponent α in the range of 3/2 to 2 is obtained.
Furthermore, we investigate the local spin-flip susceptibility which we calculate from the effective impurity model: where S eff is the same effective impurity model action from DMFT as in Eq. ( 2).At zero temperature, the spin-flip susceptibility was obtained directly on the real axis by time-evolving the matrix-product state |ψ = ĉ † ↓ ĉ↑ |GS using TEBD and then calculating the overlap χ +− (τ − τ ) = ψ(τ )|ψ(τ ) .Finite temperature results were sampled with worm-sampling in CT-HYB 61 ; the analytic continuation to real frequencies was performed using a sparse modeling approach 62,63 .Figure 3 shows the imaginary part of the susceptibilities χ +− (E) for different temperatures.For low and zero temperature, the imaginary part is gapped, i.e., it vanishes for a finite region around E = 0, in correspondence with the gapped spectral function shown in Fig. 2. For hightemperatures, on the other hand, we obtain a power-law behavior, lim E→0 (− Im χ +− (E)/E) > 0, as visible in the inset of Fig. 3; this is in agreement with the closing of the gap in Fig. 2. All curves have one peak; the peak position (in energy) seems to slightly increase with temperature.The real-frequency results show an additional shoulder around E ≈ 1.5 eV.In addition, a small satellite is found near E ≈ 3.5 eV, outside the area shown.

B. Bi-and multilayers
The starting point is the formulation of the Hamiltonian for the coupled layers.For a given number of layers l it has the form: The indices α, β denote sites within a given layer l.The first term in this Hamiltonian, containing Ĥl , describes isolated layers; analogous to Eq. ( 1), this is a sum of single-band Hubbard Hamiltonians.The second term, a double sum over nearest-neighbor layers ( Ĥll ), contains the hopping between adjacent layers as well as the interlayer Coulomb interaction.The latter is treated, for simplicity, within a mean-field approximation: This is equivalent to using Poisson's equation 64,65 to determine the potential self-consistently.Within the layers, we consider a two-dimensional square lattices as depicted in Fig. 4. The density of states in a single layers has a half-bandwidth of D = 1 eV, which corresponds to an in-plane hopping t l αβ = 0.5 eV for nearest-neighbors α, β.The inter-layer hoppings are chosen to have the same magnitude, t l,l+1 = 0.5 eV.Hence a stack of infinitely many square-lattice layers corresponds to a cubic lattice with half-bandwidth D = 1.5 eV.For the remainder of Sec.III B, we fix the temperature at T = 0.16 eV.

Bilayers
The systems studied next consist of two coupled layers; one of the layers (l = 1) is half-metallic and the other (l = 2) is either a metal, a band insulator, or a Mott insulator.The half-metallic layers have the same parameters as in Sec.III A: h l = 0.5 eV, l = −1.5 eV, and U l = 2 eV.We fix the filling of the bilayer to match the sum of the fillings of the isolated layers n iso l ; the HMF layer contributes a filling of n iso 1 = 0.355.The nearest-neighbor inter-layer hopping t 12 = t 21 = t couples the layers.
In the absence of interactions, U l = 0 and V l = 0, and in the presence of a splitting field h 1 acting on the HMF layer, the energy spectrum shows bonding (E − σ (k )) and anti-bonding (E + σ (k )) sub-bands:  For the Green's functions of the layers l = 1, 2 we get The magnetic field (h 1 ) splits the two spin channels.
For the bilayer heterostructure with one HMF layer coupled to a metallic layer (M), both layer spectral functions A lσ (E), l = 1, 2 are metallic; the gap in the minority channel of the HMF layer l = 1 closes.The essential physics is the charge transfer between the half-metallic and the metallic layer, which increases the filling in the minority spin channel of the half-metal that closes the gap.This effect also occurs in the absence of interactions.
Figure 5 shows the spectral function of a bilayer structure of a HMF layer interfaced with a band-insulating (BI) layer.The band-insulating layer l = 2 is non-interacting, U 2 = 0, non-magnetic, h 2 = 0, and completely empty, 2 = −2.25 eV, n iso 2 = 5 × 10 −5 ; these values imply a chemical potential of µ = −0.129eV.In the spin-up channel of Fig. 5a a clear bonding/anti-bonding gap is visible, while in the spin-down channel this gap is closed.The layer-resolved spectral functions show that the disappearance of the minority spin half-metallic gap is due to the interactions in the half-metallic layer.According to the HF solution of the bilayer, both layers show a gap for spin-down electrons, cf. the dotted lines in Figs.5b and 5c.The proximity to the correlated HMF layer causes the appearance of electronic states around the Fermi level of the band insulator.The many-body induced tail in the HMF is enhanced, decreasing the polarization of the HMF layer further.
Figure 6 shows the spectral functions of the bilayer formed by interfacing the HMF layer and a Mott insulating (MI) layer.Electrons in the MI layer are subject to a considerable interaction, U 2 = 5 eV, no magnetic splitting, h 2 = 0, and for the layer occupation the half-filled case ( 2 = 0, n iso 2 = 1) is considered; for these parameters, the chemical potential is µ = 0.013 eV.At the level of HF this corresponds to the interface between the half-metallic and the ordinary metallic layer as both spectral functions show states at and around the Fermi level.Within the insulating layer, Fig. 6b, the splitting into lower and upper Hubbard bands is visible (separated by ≈ U 2 ).The proximity to the HMF layer induces a slightly spin-polarized quasiparticle peak (QP) located at the Fermi level of the MI layer.In contrast, the isolated Mott layer, t ll ≡ 0, shows no QP peak for these parameters. 66,67In order to study the polarization of the QP peak we performed calculations increasing the magnitude of U 2 starting from U 2 = 1 eV.
In Fig. 7 we present the spectral function obtained for fixed parameters of the HMF layer (U 1 = 2 eV, 1 = −1.5 eV, h 1 = 0.5 eV), while increasing the strength of the Hubbard parameter U 2 = 1, 2, 3 eV towards a Mott insulator in the adjacent layer, l = 2.The quasiparticle peak and the lower and upper Hubbard bands are already seen for U 2 = 2 eV in Fig. 7b, their separation increases with increasing U 2 .The spectral function of the HMF layer shows, besides the expected satellite at about 3.5 eV, some additional spectral weight corresponding to the position of the lower Hubbard band of the Mott insulating layer.Likewise, at higher energies at the position of the upper Hubbard band a shoulder in the spectral function of the HMF layer is visible.Contrary to the homogeneous single layer, where increasing U 2 leads to a sharpening of the QP feature, the spectral weight induced by the chargetransfer seems to overlay the QP.While the spectral weight around the Fermi level decreases with increasing U 2 , it persists even for values as large as U = 10 eV.Accordingly, the double occupation of the MI layer is not completely suppressed in the bilayer case: while increasing the interaction U 2 reduces it, the double occupation is larger than in the isolated MI layer case.
We point out that we do not expect a strict Mott transition in the sense of a vanishing quasiparticle weight, respectively of a divergent effective mass.Instead, the mutual doping of Mott and HMF layer leads to metallic behavior of the whole bilayer.Thus the system favors a certain amount of charge fluctuations, and the hopping between the layers is never renormalized to zero.Such a behavior has been coined "electronic reconstruction." 68The common feature of these results indicates that the transfer of charge between the layers is a general phenomenon that produces metallic interfaces.

Half-metallic and Mott insulating multilayers
In the following, we scale up the system size and consider a heterostructure made up of five HMF layers coupled to five Mott insulator layers.We consider open boundary conditions.In order to preserve the halfmetallic gap, we scale the on-site parameters of the halfmetallic layers by a factor of two in comparison to the previous bilayer calculations, Sec.III B 1: U l = 4 eV, h l = 1 eV, and l − µ = −3 eV for all layers l ∈ {1, . . .5}.For the Mott insulating layers, we choose the same Hubbard interaction, U l = 5 eV, for all remaining layers l ∈ {6, . . .10}. Figure 8 shows the layer-resolved spectral function, A lσ (E), for this setup.The many-body effects in the half-metallic layers (Figs.8a and 8b) are qualitatively the same as in bulk: we observe a dynamical reduction of the Hartree part of self-energy, and a tail crossing the Fermi level.Within the HF approximation, the layers l = 6, . . . 10 (i.e., on the MI side, Figs.8c and 8d) are found to be metallic, as to be expected; in addition, the charge-transfer at the interface introduces small weight at the gap in the interface layer on the HMF side, l = 5 (Fig. 8b).
Within DMFT, we see that the spectral weight in the MI surface layer, l = 6, is strongly suppressed around the Fermi level, however, the layer remains metallic despite the strong interaction.On the other hand, a significant shift of spectral weight towards the Fermi level is apparent in the spin-down channel of the HMF interface layer, l = 5.Nevertheless, the polarization in this layer (l = 5) is close to the polarization obtained within HF.The closing of the Mott gap observed in the interface layer on the MI side is similar to the bilayer, cf.Fig. 6.This effect has the range of two layers, at l = 8 the gap is apparent again.The minimum of the spectral function A 6σ , Fig. 8c, of the MI interface layer shifts from zero energy to roughly E ≈ −1 eV for both spin channels.Contrary to the bilayer result, Fig. 6, there is no QP peak at E = 0, neither at the interface nor in the subsequent MI layers.Surprisingly, the spectral function shows a shoulder at the Fermi level for the down-spin only.
Next we include the long-ranged Coulomb repulsion in mean-field approximation for this multilayer setup.We apply the algorithm described in Ref. 64 and 65; the formula for the potential reads with the layer occupation n l = n l↓ + n l↑ ; the material parameter, e Schot , can be related to the screening length, as discussed previously. 64,65We use, however, a different update scheme to solve the Poisson equation, which avoids the thousands of iterations 64,65 necessary with a naive mixing scheme.Instead, after every DMFT iteration we temporarily fix the self-energy, Σ lσ (iω n ), to a selfconsistent potential V l .We start from the occupation numbers n l (e.g., given by the last DMFT iteration, or the non-interacting result).From the occupations, we calculate the potential using the above equation, where we introduce the vector notation V = {V l }, n = {n l }.Given the potential and the self-energy, we can calculate a new Green's function, with the vectors From the Matsubara sum of the Green's function we then calculate new occupations, giving us the self-consistency equation After every DMFT step, we solve for the self-consistent charge n l , Eq. ( 13), and therefore for a self-consistent potential V l for the given self-energy.This method significantly reduces the number of required DMFT iteration, however, it introduces costs for solving Eq. ( 13) after every iteration.The main cost of Eq. ( 13) is the evaluation of the Green's function matrix.Numerically, it is more efficient to solve the equivalent root-search problem, r(n * ) = 0, for the function A Newton-Krylov solver 69 , as implemented in Ref. 70, is found to be most suitable for this problem.Furthermore, we also include the search for the chemical potential, µ, necessary for fixing the total charge and therefore guaranteeing charge neutrality, l n l = l n bulk l , when performing the root search, Eq. ( 14).This is easily implemented using the modified equation where we append a row for the difference in total charge to the vector-valued function r(n). 71e fix the material parameter to the constant value e Schot l = 1 eV.The bulk occupations are n l = 0.205 for the HMF layers, l ∈ {1, . . .5}, and half-filling (n l = 1) for the MI layers, l ∈ {6, . . .10}.
. Change in spin-resolved occupations due to interlayer Coulomb mean-field potential, V l Eq. (10).The graph shows the deviation of the spin-resolved occupation of the layers from half the bulk occupation n bulk l /2.The blue marks show the occupations for V l ≡ 0, while for the yellow marks the self-consistent mean-field potential is included.The vertical bar indicates the interface between the HMF region (l ∈ {0, . . .5}) and the MI region (l ∈ {6, . . .10}).
Figure 9 shows the spin-resolved occupation for the multilayer, without (blue) and including (yellow) the longranged Coulomb repulsion.Apparently, for the present setup the long-ranged repulsion is of little relevance.The occupation of the metallic spin-channel σ =↑ is slightly affected, due to charge fluctuations in this channel, while the rest is nearly invariant with respect to inclusion of long-ranged effects.Likewise the spectral function is nearly identical to Fig. 8.

IV. CONCLUSION
In summary, we have presented detailed model studies for the spectral properties of bulk half-metallic ferromagnets (HMFs) as well as for bi-and multilayers containing half-metallic ferromagnets.Dynamical mean-field theory has been employed to describe the local correlations between charge carriers, while a mean-field approach was used to include the long-ranged Coulomb interactions.
Our numerical results show that the correlation-induced tails in the vicinity of the Fermi level in bulk HMFs are significantly reduced at zero temperature, in agreement with analytical predictions 2 .On the other hand, for biand multilayers we find an enhancement of the tail contribution at the half-metallic side, as well as coherent quasiparticle states on the Mott insulating side.In the multilayers these mobile carriers are confined to a relatively narrow region at the interface.Furthermore, the Fermi liquid states at the interface reduce the full spin polarization characteristic for bulk HMFs.Note that the formation of Fermi liquid states at such interfaces is similar to the LAO/STO interfaces, which have been theoretically [72][73][74][75] studied and experimentally observed 76,77 some time ago; however, to the best of our knowledge, such effects have not been studied for heterostructures containing HMFs.
On the technical side, we have demonstrated that the real-space DMFT allows, in a rather transparent way, the inclusion of long-ranged Coulomb interactions via the Poisson equation.In this approach, the charge distribution in the presence of strong short-range interactions and spatially inhomogeneous hoppings is determined selfconsistently.In contrast to previous implementations of the R-DMFT 36 , we use the Hubbard model and a stateof-the-art CT-QMC 39 implementation for the impurity solver.The Poisson equation is solved as an effective onedimensional problem in combination with the R-DMFT self-consistency condition as discussed above.For the biand the multilayer setup, we have considered the case where the layers can be modeled as square lattices.For our bilayer setup, we have considered a half-metallic monolayer in contact with either a metal, a band or a Mott insulator.We have seen that charge reconstruction at the interface causes the existence of metallicity, even in the presence of large Hubbard U parameters at the Mott insulator layer.In the R-DMFT analysis the HMF/MI bilayers are Fermi liquids with well defined quasiparticles, thus the present approach offers a way to access Fermi liquid quantities on the basis of a microscopic model.
On the experimental side, most studies have concentrated on the question of whether the half-metallic properties extend to the surface or interface.Using DFT calculations, a genuine half-metallic interface of NiMnSb with InP and CdS has been obtained only for the anion terminated (111) direction 78 .Interfaces of semi-Heuslers NiMnSb or NiMnSi with large gap insulators such as MgO have been also studied 79 .A high spin polarization has been obtained only under the prerequisite of structural optimization 79 .
However, the microscopic origin of the HMF/Mott insulator interface has never been addressed.In this context, we thus considered a minimal model in which half-metallic layers are in contact with correlated insulator layers.We solved the corresponding Hubbard Hamiltonian in Hartree-Fock (HF) approximation and beyond using dynamical mean-field theory (DMFT).Within the HF approximation, when crossing the interface from the HMF side into the metallic side, we find a sharp transition, i.e., the half-metallic layer is followed directly by a metallic layer.In contrast, by including dynamical correlations within DMFT, we find a continuous transition from the half-metallic region through a pseudo-gapped interface into an insulating region.Our simplified model thus indicates that a high spin-polarization within the interface region can be preserved in the presence of correlated Mott insulators.

ACKNOWLEDGMENTS
The zero-temperature calculations were performed at the Vienna Scientific Cluster (VSC).A. Weh   In the following, we show that a Padé analytic continuation of the Matsubara Green's function instead of the self-energy leads to artifacts in the spectrum.The self-energy, Eq. (3) in Sec.III A, for T = 0.25 eV can be fitted by the two pole function with the residues w σj and poles σj given in Table I.Below, we will use this analytic expression as a realistic test case for the quality of analytic continuation.First, we need to fit the parameters in Eq. (A1).They can be obtained using Padé analytic continuation of the self-energy Eq. ( 3), as it yields an analytic formula (with numerical coefficients) in form of a rational polynomial f (z) = p(z)/q(z), with polynomials p and q. 59 The poles σj can be calculated as the zeros of the denominator q.Instead of using Thiele's reciprocal difference method to determine the rational polynomial, we directly calculate the poles in Eq. (A1) employing the algorithm presented in Ref. 80.We write the linearized Padé approximation f (z)q(z) = p(z) in matrix form, where (F ) ij = f (z i )δ ij is the diagonal matrix of function values, V q and V p are the Vandermond matrices corresponding to q and p, and q and p are the polynomial coefficients.We rewrite the equation as The number of poles M is then determined such that the numerical null-dimension of the matrix C is one.To calculate the poles, we rewrite the rational polynomial f (z) = p(z)/q(z) by factorizing a pole m from q(z) = (z − m )q m (z), which leads to Again, we rewrite this set of equations in matrix form: This is a generalized eigenvalue problem, where the eigenvalues are the poles m , and the eigenvectors are the coefficients of the polynomials p and q.Using the knowledge of the poles σj , the residues w σj can be obtained solving the linear equation in w σj .To reproduce noisy input, Padé places artificial poles along the imaginary axis.We verify that the residues of these poles are small and neglect them subsequently, and recalculate the residues including physical poles only. 81sing Eq. (A1), we can evaluate the spectral function directly on the real axis A6) where the integral evaluates for the Bethe lattice with infinite coordination number to We compare it with the Padé analytic continuation In Fig. 10a, we use the fitted Σ σ (z) as input to Eq. (A8).The Figure shows that the analytic continuation of the Matsubara Green's function Eq. (A8) leads to spurious features similar to those in Ref. 24, in stark contrast to the analytic continuation of the self-energy used in Fig. 1. Figure 10b shows the Padé continuation using the noisy raw data Eq.(3) as input to Eq. (A8).The result qualitatively agrees with the previous results  The HF approximation is given as reference (dotted line).
We conclude that, independent of the presence of noise, Padé is unable to reproduce the branch-cut of the noninteracting Bethe DOS ρ(E) from Matsubara frequency data Eq.(A8), as it approximates the function with a finite number of poles.The sharp band-edges cannot be resolved, and oscillations similar to the Gibbs phenomenon in Fourier transform occur.We can further support this argument by looking into the analytic continuation of the hybridization function.From Eqs. (A1) and (A8) we can calculate the hybridization function on the Matsubara axis,  In the following, we derive the R-DMFT equations 35,36 used in Sec.III B. We use an expansion in the coupling between the layers in the action formalism.
We start from the action S of the multilayer heterostructure, which we split into two parts, S l is the action of the isolated layer l, and ∆S contains the hopping in between the layers.To ease the notation, we introduce the convention that indices with a bar are summed over: l is summed over layers, ᾱ, β over sites within a layer, and σ over spins.The contributions to the action Eq.(B1) reads Next we introduce auxiliary fields and expand the contributions of the isolated layers in cumulants.We truncate this expansion after the first order, keeping only the quadratic part in the auxiliary fields.This approximation results in a self-energy which is diagonal in the layers, Σ lα;l β = Σ l αβ δ ll .We perform the Grassmannian Hubbard-Stratonovich transformation following Ref.82.We rewrite the exponential e −∆S as Gaussian integral over auxiliary Grassmann fields ψ lασ (τ ) and ψ + lασ (τ ): Here the second equality is the Woodbury matrix identity 84 .Thus we can calculate the Green's function G of the full heterostructure from the Green's functions of the isolated layers G l .Written in terms of the self-energy, we get where the self-energies Σ l are determined by the Dyson equations of the isolated layers, We determine these self-energies Σ l by using DMFT for every distinct isolated layer.This approximation is evidently correct in the limit of isolated layers, T = 0, and in the limit of non-interacting layers, U l = 0.This approximation is supplemented by self-consistency.
Self-consistent R-DMFT equations.We are going to replace the self-energy of the isolated layers, Σ l , by a selfconsistently determined self-energy.We can also write Eq. (B18) in the form i.e., the inverse Green's function is written as sum over the diagonal and the off-diagonal elements.The Green's function then reads: This looks like a self-consistency equation, because it involves the same Green's function G on the left-and righthand-side of the equation.Finally, the Dyson equation for the self-energy reads We determine the self-energy Σ l by a DMFT scheme 35 , instead of using Σ l from Eq. (B19).

Figure 1 .
Figure 1.Spin-resolved spectral function Aσ(E) for the bulk half-metal.Part (a) shows the overall behavior, while (b) and (c) focus on the region around the Fermi level and the high energy satellite, respectively.Black dotted lines correspond to the Hartree-Fock (HF) approximation.At high temperatures (T = 0.25 eV) the tail of the spectral function A ↓ (E) crosses the Fermi energy, while for low enough temperatures (T ≤ 0.02 eV) the half-metallic gap is preserved.The inset in (c) shows the T dependence of the spectral weight of the ↓-spin.Crosses indicate finite temperature results at T = 0.2, 0.4, 0.8, 0.16, 0.25 eV; the circle corresponds to T = 0.

TFigure 2 .
Figure 2. Evolution of the tail in the minority spin spectral function A σ=↓ with temperature.The temperature decreases on a logarithmic scale from the front to the back.The last line in the back corresponds to the low T result of Fig. 1.

Figure 3 .
Figure 3. Imaginary part of the local spin-flip susceptibility at zero temperature (dashed green line), and for a selection of finite temperatures (T = 0.02 and 0.08 eV: purple and blue solid lines).The inset shows the susceptibility divided by energy: − Im χ +− (E)/E. t

Figure 4 .
Figure 4. Illustration of coupled monolayers of square lattices.The in-and inter-layer hopping integrals are indicated.

Figure 5 .
Figure 5. (a) Spin-resolved total spectral function l A lσ (E) for one HMF layer interfaced with one band insulating (BI) layer.The solid lines are DMFT (CT-HYB), and the dotted lines the HF results.The green lines show the spectral function for isolated layers (t ≡ 0).(b) Spin-resolved spectral function A lσ (E) for the HMF layer of the bilayer, and (c) for the BI layer.

Figure 6 .
Figure 6.Spin-resolved spectral function A lσ (E) for one HMF layer (a) interfaced with one Mott insulator (MI) layer (b).The solid lines are DMFT (CT-HYB), and the dotted lines the HF results.The green lines show the spectral function for isolated layers (t ≡ 0).(c) Spectral weight at the Fermi level A lσ (E = 0) as function of the hopping t between the layers.

Figure 7 .
Figure 7. Spin-resolved spectral function A lσ (E) for one HMF layer interfaced with one layer of different interacting strengths U2.The green lines show the spectral function for isolated layers t ≡ 0.

Figure 8 .
Figure 8. Spin-resolved spectral function, A lσ (E), for five HMF layers interfaced with five MI layers.(a) shows the results for the first four layers (HMF); the spectra are shifted, their respective baselines are plotted in the same color.Panels (b) and (c) represent the HMF interface layer, l = 5, and the MI interface layer, l = 6, respectively; in both cases, the HF results are included as dotted lines for easy reference.The results for the last four MI layers are shown in (d); again the spectra are shifted.The spectral functions are evaluated by analytic continuation, E → E + 0.04i eV.

Figure 10 .
Figure 10.(a) Comparison of the Padé analytic continuation of the Matsubara Green's function, Eq. (A8), with the direct evaluation of Eq. (A6), using the fit Eq. (A1) as input to both.(b) Padé analytic continuation of the raw data of the Matsubara Green's function not employing the fit Eq. (A1).The HF approximation is given as reference (dotted line).

) 1 E
This function encapsulates the effect of the DOS ρ(E); if we perform the Padé analytic continuation ∆ σ (E), we get the spectral function− σ − ∆ σ (E) − Σ σ (E)(A10) which shows the same features as A Padé σ (E) = − Im G σ (E)/π in Fig.10.An analytic continuation of the self-energy avoids this problem; the analytic expression for ρ(E) is used directly, and the self-energy lacks such sharp features.