A charge self-consistent LDA+DMFT study of the spectral properties of hexagonal NiS

The electronic structure and spectral properties of hexagonal NiS have been studied in the high temperature paramagnetic phase and low temperature antiferromagnetic phase. The calculations have been performed using charge self-consistent density-functional theory in local density approximation combined with dynamical mean-field theory (LDA+DMFT). The photoemission spectra (PES) and optical properties have been computed and compared with the experimental data. Our results show that the dynamical correlation effects are important to understand the spectral and optical properties of NiS. These effects have been analyzed in detail by means of the computed real and imaginary part of the self-energy.


Introduction
Electron correlation effects in materials have over the past several decades emerged as one of the most important areas of current condensed matter research. Systems with localized electrons, usually associated with narrow 3d or 4f bands, are fundamental for magnetism, superconductivity and several other exotic phenomena. For weakly correlated materials where the band width (W) is much larger than the Coulomb interaction strength (U), the density functional theory within the local density approximation (LDA) or the generalized gradient approximation is usually successful in obtaining the ground state properties [1]. However, LDA/LSDA becomes inadequate in the presence of strong correlations [2,3] and the LSDA+U method [4][5][6] was developed to include the static Hartree-Fock like part of these effects. In the LSDA+U approach, an orbital dependent interaction term, describing the screened Coulomb interaction (U) between the correlated orbitals, is added to the LSDA functional. This method works surprisingly well in the limit ≫ U W and many interesting effects, such as orbital and charge ordering [7][8][9][10][11] in transition metal compounds, are successfully captured by the LSDA +U approach (for a review, see [12].) While there exist methods to study weakly as well as strongly correlated (Mott insulators) systems, addressing the physics of materials in the intermediate correlated regime ( ≃ W U) is still challenging. In this respect, dynamical mean field theory has proved to be a breakthrough for the realistic modeling of correlated materials due to its applicability in all the three regimes. The combination of conventional LDA based electronic structure calculations and dynamical mean field theory, the so called LDA+DMFT scheme, has emerged as a powerful method to understand the electronic structure and excitation spectrum of strongly correlated metals as well as insulators [13][14][15][16][17][18][19].
In view of the success of the LDA+DMFT method in computing the ground state and excitation spectrum of correlated metals, hexagonal NiS seems to be an ideal candidate to be studied in this framework. At a critical temperature T t = 265 K, this compound exhibits a first order magnetic phase transition [20], and transforms from a high-temperature Pauli paramagnetic metal to a low-temperature antiferromagnetic metal. The transition from the paramagnetic state to the antiferromagnetic state is accompanied by a volume enhancement of about 2%, where the hexagonal lattice parameters a and c increase by 0.2% and 1%, respectively [20]. A neutron diffraction study [21] in the low temperature phase revealed that all the Ni atoms are ferromagnetically coupled in the ab-plane and antiferromagnetically coupled along the c-axis (A-type antiferromagnet.) In the literature, there are several theoretical [4,[22][23][24][25][26][27][28] and experimental studies [21,[29][30][31][32][33][34][35][36][37] of the electronic structure of NiS for both the phases. It is possible to show that the ground-state magnetic properties of NiS are well described when considering correlation effects explicitly, even at the simplified level of LSDA +U. In fact, recent LSDA+U calculations have revealed that NiS is a novel example of a dilute self-doped antiferromagnetic p-d metal [28]. However, the inclusion of only static correlation effects as embodied in the LSDA+U approach may not be adequate to quantitatively reproduce all features of experimental photoemission spectra or optical conductivity. This was already exemplified by the fact that the optical conductivity of NiS in the low temperature phase calculated within LSDA+U required a contraction of the energy scale by about 30% for an agreement with the experimental data [28]. In addition, the calculated density of states in [28] does not have good quantitative agreement with experimental photoemission spectra as far as the main peak position is concerned. This points to the fact that a more sophisticated treatment of many-body effects is needed for understanding the electronic structure of NiS and this is indeed the main point of this investigation.
Recently, Amadon et al [38] performed LDA+DMFT calculations within the projector augmented wave and the mixed basis pseudo-potential approach for the high temperature phase of NiS. This study was restricted to the evaluation of the local spectral functions for selected values of U and did not focus on any direct comparison with experiments, which is instead the main aim of this work. In addition, the previous study emphasized the importance of the S 3p states and suggested further investigation of correlation effects in NiS. Investigation of the dynamical correlation effect for the low temperature metallic AFM phase of NiS was also beyond the scope of [38].
In the present paper, we have used the LDA+DMFT method to study the electronic structure of NiS in both the high temperature and low temperature phases. These results have then been employed to calculate the photoemission spectra, as well as optical properties. Comparison with experimental data clearly show that dynamical correlation effects are necessary to fully understand the spectroscopic and optical properties of NiS.

Computational details and crystal structure
All of the LDA/LSDA calculations reported in this work have been carried out using a full potential linear muffin-tin orbital method [39]. The Brillouin zone integration was performed using the thermal smearing method with 131 k points in the irreducible part of the Brillouin zone. The converged LDA/LSDA simulations were used as a starting point for the LDA +DMFT calculations, which we performed by means of the implementation presented in [39][40][41]. In the following we will summarize the most important technical aspects relevant for this study, and we direct the reader to the aforementioned references for a more extended description.
The first step of the LDA+DMFT method is to identify a set of local orbitals, which are not properly described by LDA/LSDA, and which for convenience we label with a quantum number ξ. Then the Kohn-Sham Hamiltonian Ĥ KS is corrected with an explicit Hubbard term describing the local Coulomb interaction U between those states The last term is the double counting correction, which removes from the Kohn-Sham Hamiltonian those contributions that are associated to the added local Coulomb interaction U in order to avoid counting them twice [15]. The new Hamiltonian is an effective Hubbard model, where the one particle states contain all the states included in the DFT problem (no downfolding is done). A crucial point in this scheme is the choice of the local orbitals ξ. Here we use the muffin-tin heads (see MT basis in [40] for a more comprehensive description) for the Ni-3d states. The Hubbard model can then be solved through DMFT, and convergence is considered for both the local self-energy and the full electron density. The effective impurity problem arising in the DMFT cycle was solved through the spin polarized T-matrix fluctuation-exchange (SPTF) solver [42]. The SPTF solver was chosen due to its efficiency and to that it presents many technical advantages, which are discussed in the next section. The accuracy of the SPTF solver for the moderately correlated systems ( ⩽ U W) is indirectly proven by several successful applications to various materials [40,[42][43][44][45][46]. In the present paper we use the SPTF version, whose perturbative contributions are calculated by means of the self-consistent Hartree-Fock Greenʼs function. The double counting correction ξ ξ H , DC 1 2 has been set as the orbitally averaged static part of the self-energy, as commonly done for calculations based on the SPTF solver [40,47]. The high and low temperature phases have been simulated for temperatures of T ≃ 300 K and T ≃ 150 K, respectively. A simultaneous convergence of the charge and self-energy has been achieved with 2048 Matsubara frequencies. The fully rotationally invariant U-matrix has been constructed from the Slater parameters F 0 , F 2 and F 4 . These have been reduced to only two parameters, the Hubbard U and the Hund exchange J, by using fixed atomic ratios [48]. Concerning the value of U, we have performed the calculations for two possible choices, i.e., U = 2.0 eV, and U = 3.0 eV, within the range suggested in a previous study [28]. As we shall see below, the calculations based on the value U = 3.0 eV compare best to the measured data. The value of J was fixed to 0.8 eV for all the calculations. To extract the density of states, the Padé approximation method has been used for the analytical continuation of the self-energy function from the Matsubara frequencies to the real energy axis.
The photoemission spectra have been calculated within the so-called single-scatterer finalstate approximation [49,50]. Here the photocurrent is a sum of local (atomic-like) and partial (llike) density of states weighted by the corresponding cross sections. The cross sections are calculated using the muffin-tin part of the potential over the energy range of the corresponding DOS for the fixed photon energy. As the cross section depends on the photon energy, the theoretically computed spectra will depend on the photon energy. The spectra has been broadened with a Gaussian to simulate the effect of spectrometer resolution.
Hexagonal NiS crystallizes in the NiAs structure. The lattice parameters [20] are a = 3.4456 Å, c = 5.405 Å for the low temperature phase, and a = 3.4395 Å, c = 5.3514 Å for the high temperature phase. Each Ni 2+ is in an octahedral environment of sulfur (S) atoms and these NiS 6 octahedra are edge shared within the ab plane and face shared along the c axis forming a three-dimensional network of Ni-S connectivity in this system. All of the calculations reported here have been performed with the above structural parameters for the respective phases.
3. Electronic structure of the high temperature paramagnetic phase As described below, a value of U = 3.0 eV seems to reproduce the valence band spectrum of both the PM and AFM phase. Hence, we discuss here the calculated spectral properties and quasi-particle dispersion for this value of U. The LDA and LDA+DMFT (U = 3.0 eV) partial densities of states (PDOS) of Ni 3d and S 3p are shown in figures 1(a) and (b), respectively. It follows from figure 1(a) that DMFT leads to a noticeable renormalization of the Ni 3d PDOS, by narrowing and shifting the main peak toward a lower binding energy. The DOS at the Fermi level does not change significantly, while a (weak) high energy satellite is formed at about −9 eV. A closer look at the Ni 3d spectrum reveals a small renormalization of the peak at −4 eV, whose position is pinned by the S 3p states. Notice that, due to the strong hybridization, such an effect is also visible in the S 3p partial DOS ( figure 1(b)), and these states are slightly shifted toward the Fermi level. Changes in the PDOS due to the inclusion of the dynamical correlations can be understood by analyzing the real and imaginary part of the energy dependent selfenergy, shown in figures 1(c) and (d), respectively. The quasiparticle energies are renormalized by the real part of the self-energy that shifts the peak positions of the LDA spectrum. As seen in figure 1(c), the amplitude of the real part of the self-energy is quite high and positive within the energy range from 0 eV to −7 eV which gives the shift in the DOS spectrum toward the Fermi level (Fermi liquid renormalization). The sharp negative peak at ∼ −9 eV leads instead to the formation of the high-energy satellite. The imaginary part of the self-energy, shown in figure 1 (d), causes broadening of the peaks, leading to an effective decrease of their intensities. Around the Fermi level the imaginary part of the self-energy is small, implying that the quasi-particles have long lifetimes, in accordance to Fermi liquid theory. For higher energies, we can observe the signature of a satellite formation in the imaginary part of the self-energy, i.e. the peak at −8 eV. This high energy peak, observed for both the real and imaginary part of the self-energy, is analogous to what is found in the late transition metals elements, especially for fcc Ni, and has been object of much investigation [40,47]. In DMFT the satellite peak is related [16] to the formation of Hubbard bands, i.e. noncoherent features in the spectrum corresponding to atomiclike transitions, although some authors have suggested [51] that the effective Hubbard model with an energy-independent U-matrix, i.e. the one used in the LDA+DMFT approach, could lead to the formation of a spurious plasmon excitation around those energies. In general, the behavior of the self-energy of NiS is similar to what is observed in Fe, Co and Ni [40] when using LDA+DMFT and the same SPTF solver. Since NiS, bcc Fe and fcc Ni have similar bandwidth and size of U, it seems reasonable that the influence of electron correlation and the overall shape of the self-energy are similar for the three materials, which is what we observe.
The features outlined in the previous paragraph can be clearly seen in the k-resolved spectral density (figure 2). The extent of energy-smearing, as an indication of width-broadening, is evident at the energies where the imaginary part of the self-energy is large, mainly at high excitation energies. In contrast, a well-defined, sharp dispersion relationship is seen close to the Fermi energy. It should be noted that the bands in figure 2 are hybridizing states, involving primarily correlated Ni d-states and uncorrelated S p-states, where the admixture of these two different orbitals depends on position in k-space. Hence, one observes that the correlation induced life-time broadening reflects this behavior, so that for a given binding energy states at some parts of the Brillouin zone have more broadened states (the Ni dominating ones) than other parts (states richer in S p-character.) Before discussing the calculation of the photoemission spectrum, it is relevant to compare the presented LDA+DMFT results for the high temperature phase with those of [38]. We must first emphasize that these two studies present several technical differences. The main computational scheme in the present work is somehow more complete, because the LDA +DMFT equations are solved directly in the full Hilbert space (without any truncation of the basis set) until convergence of the full self-consistent cycle in the electron density. This aspect is important to describe the feedback effects of the Ni-3d self-energy into the hybridization between Ni-3d and S-3p states, although these effects are not very large. More significant differences are related to the solution of the effective impurity model arising in DMFT. In [38] the authors have used the Hirsch-Fye quantum Monte-Carlo solver (HF-QMC) [16], which is fully general with respect to the strength of U, and is numerically exact, i.e. it leads to the exact solution in the limit of infinite computational effort. However, in practical applications the HF-QMC solver has some limitations. First of all the U matrix is restricted to only contain the density-density contributions. Moreover, the analytical continuation of the Greenʼs function from the Matsubara axis (where the the real calculation is done) to the real energy axis (used for the spectral properties) is performed by means of the maximum entropy method. The latter is very sensitive to the parameters chosen for the analytical continuation, and also introduces an additional smearing of the spectral features. The latter can become relevant for high energy excitations far from the Fermi level [52,53]. Keeping these factors in mind, and considering that the simulations of [38] were carried out at a temperature of around 1100 K and for U = 4 eV, the agreement between our study and [38] is rather satisfactory. The main difference concerns the position of the Hubbard satellite, which is located at slightly higher excitation energies in our simulations. This is a typical deficiency of the SPTF solver, and is observed also for other materials, most notably fcc Ni [40]. We must finally notice that when the Feynman diagrams defining SPTF [42] are constructed by means of the bare propagator (here coinciding with the bath Greenʼs function), then the effective Coulomb interaction U is slightly under screened [54]. This is the reason why in our computational scheme we need to use a slightly reduced value of U when comparing with HF-QMC results. To complete our analysis, in figure 3 we report the orbitally resolved PDOS for the Ni-3d states, and corresponding real and imaginary parts of the self-energies. We note that the m l projected spectral functions and selfenergies actually are quite different from each other. As figures 3(a), (b) and (c) show, both the spectral intensity and features of the self-energy are m l dependent. In particular, the m l = ±1 projected imaginary self-energy has a significantly weaker feature at −8 eV, which is due to that these orbitals have a stronger hybridization with S p-states. This hybridization leads to an effective broadening of the band-width for d-states with = ± m 1 l . Hence, for these states the influence of a Hubbard U is less significant, compared to more narrow partial DOS, and this is , for all energies. The importance of the correlation effects associated with the self-energy becomes clearer when comparing our theoretical spectroscopic data with the experimental ones. Computed photoemission spectra within LDA and LDA+DMFT are shown along with the experimental data in figures 4(a) and (b), respectively. From these figures, it follows that the inclusion of the dynamical correlation effects significantly improves on the LDA results. Photoemission spectra have been computed for two reasonable choices of U to find out the best choice of U which can correctly describe the spectral properties of NiS within the LDA+DMFT framework. Our results, as displayed in figure 4(b) reveal that both the choices of U improve the description of the experimental spectra quite significantly but the best description is obtained for U = 3.0 eV.

Electronic structure of the low temperature antiferromagnetic phase
The low temperature antiferromagnetic phase has been studied in detail theoretically and experimentally by various groups. It is well established that LSDA cannot stabilize the A-type antiferromagnetic ground state and the inclusion of U in the framework of LSDA+U is essential to obtain the observed magnetic structure [27,28,34]. It has been difficult to study the magnetic phase of NiS by means of the LSDA+DMFT scheme with numerically exact QMC solvers. In fact the need of including the entire d-shell at a relatively low temperature would lead to a prohibitive computational effort. The SPTF solver does not suffer from this limitation, and therefore we can investigate the importance of the dynamical many-body effects in the low temperature phase. In figure 5, we compare the Ni 3d and S 3p PDOS for both spin channels, obtained using the LSDA+U and LSDA+DMFT methods with U = 3.0 eV. Figure 5 (a) shows that the PDOS for both the spin channels of Ni 3d states calculated with LSDA+DMFT deviate substantially from those obtained within LSDA+U. Here it is also clear that with the inclusion of dynamical correlations the prominent peak of Ni 3d states in both the spin channels becomes narrow and shifts toward the lower binding energy by about 0.8 eV and the DOS at the Fermi level increases. Due to the strong hybridization, the Ni self-energy also affects the S 3p states at around −4 eV, as shown in figure 5(b), shifting them toward the Fermi level. The sizable modification of the peak is related to the strong amplitude of the real part of the self-energy, shown in figure 5 (c). Sharp features, visible in LSDA+U PDOS, are smoothened out in LSDA +DMFT results, due to the quasiparticle life-time effects associated with the strong imaginary part of the self-energy ( figure 5(d)). In figure 6, we show the k-dependent, spin-integrated spectral density along different high symmetry directions. As expected, the quasi particle features are still well defined near the Fermi level, where there are electron pockets and hole pockets, as was observed in the LSDA+U band dispersion [28]. These spectral features support the identification of NiS as a dilute self-doped nearly-compensated metal, as suggested in [28]. Overall the spectral densities reported in figure 6 show several differences with that shown in figure 2 due to the AFM ordering. However, a significant difference between low and high temperature phases is the absolute value of the self-energy. The values in figures 5(c) and (d) seem to be a factor of 2 smaller than those reported in figures 1(c) and (d). These differences are too large to appear physical, given that the values of the bandwidth and the U in the two phases are rather similar. The origin of this discrepancy is most likely due to that an unordered local moment arises in the high temperature phase. In fact, as discussed in [56], SPTF performs well when starting from a magnetic solution, since the exchange splitting includes already a large part of exchange and correlation effects. On the other hand, when starting from a non-magnetic phase, the impurity solver tries to reproduce the strong, and essentially mean-field, effect associated to the formation of a local magnetic moment through dynamical fluctuations around a non-spin polarized state. This description is adequate for a solver with no formal restriction on the amplitude of the dynamical fluctuations, such as QMC, but is clearly inappropriate for a method based on perturbations around the non-interacting limit, such as SPTF. As a result, one observes an overestimation of the self-energy correction, as was previously analyzed for γ-Mn [56].
As we discussed above, the general spectral features observed in LSDA+DMFT are also captured by LSDA+U. However, the lack of dynamical correlation effects has a clear impact on the calculated photoemission spectrum. It should be noted here that the increase of inelastic scattering of the photoelectrons at lower photon energies reduces the applicability of the singlescatterer approximation [18]. However, at moderate inelastic scattering the main effect is an asymmetry in the peak shapes, while the positions remain relatively unchanged. Since the main focus of our comparison between theory and experiment is the positions of the peaks, the same single-scatterer approximation as mentioned earlier has been used to calculate the PES of the low temperature phase for a photon energy of 21.2 eV. The theoretical PES spectra for LSDA +U and LSDA+DMFT with U = 3.0 eV are shown in figure 7, where they are compared with the experimental data. Figure 7(a) clearly shows that the position of the theoretical peaks obtained from the LSDA+U calculation deviate substantially from the measured ones. The inclusion of the dynamic correlations, however, significantly improves the computed spectra, reproducing the main peak at the correct position. Since the shifts in the peak positions are driven by the real part of the self-energy, as already discussed with reference to figures 5(a) and (c), the substantially improved agreement of the LSDA+DMFT results with the experimental spectrum is a clear evidence for the importance of the dynamical correlation effects in describing the low temperature phase of NiS.
Finally, we have computed the optically allowed joint density of states (JDOS) [57] by convoluting the occupied and unoccupied parts of the partial density of states, and imposing the dipole selection rules. We have used the following expression of the JDOS.
where, ρ l and ρ l + 1 are the PDOS corresponding to l and + l ( 1) angular quantum numbers, respectively. All of the possible spin conserved → + (JDOS) l l 1 have been added to get the total JDOS. While we do not have the possibility of calculating the matrix elements relevant for the optical conductivity data within the LSDA+DMFT approach, we note that the JDOS, when weighted by the matrix elements, yields the optical data. Thus, it is expected that the qualitative features of the optical data to be reasonably well captured by the energy distribution of the JDOS. The calculation has been restricted to transition from l to + l ( 1) angular quantum numbers due to the fact that the transition probability is much larger for l to + l ( 1) transition than the l to − l ( 1) transition [58]. Such a calculation of the optically allowed JDOS without matrix elements within LSDA+U and LSDA+DMFT will provide direct evidence of the importance of dynamical correlations for the optical property of NiS. The absorptive part of the optical conductivity, calculated within LSDA+U and LSDA+DMFT ( = U 3.0 eV), are displayed in figure 8 along with the experimental data. In figure 8 it may be noted that the DMFT spectrum has a main peak shifted by about 0.50 eV toward lower energies in comparison to the LSDA+U result. This downward shift brings the calculation in closer agreement with the experimental spectrum, although the agreement is not perfect. This discrepancy is possibly due to the neglect of the matrix element effects in our calculations or may be due to the approximations introduced by the SPTF solver. The computed Drude peak height (see the inset of figure 8) also shows a good qualitative agreement with experimental data. We note here that the sharp peak in the optical conductivity, observed around 0.04 eV photon energy is likely to result from excitation of an optical phonon [59]. Our calculation cannot capture this sharp peak since it is beyond the scope of the present calculation. An improved description of the optical data within the present LSDA+DMFT approach further shows the importance of dynamical correlation effects for the low temperature phase of NiS.

Summary
A detailed analysis of the low and high temperature phase of NiS reveals that LDA and LDA+U are less suited to describe in detail the electronic structure of this system. Dynamical correlations are important to reproduce the main features of the photoemission spectra, and we find that DMFT in combination with a solver based on spin polarized T-matrix fluctuationexchange is sufficient for this. Peak positions are well reproduced by the incorporation of Figure 8. The absorptive part of the optical conductivity of the low temperature phase, both from theory and experiment. The experimental data has been taken from [59]. Inset shows the same for narrow energy range close to 0 eV to clearly display the Drude part of the optical conductivity. dynamic correlations within the LSDA+DMFT approach, if U = 3.0 eV is used. With these parameters both spectra based on electron spectroscopy as well as from the optical conductivity are reproduced. We also find that, due to hybridization between Ni d-states and S p-states, the importance of electron correlations and life-time broadening is anisotropic over the Brillouin zone.