Phase diagram of a pseudogap Anderson model with application to graphene

The Kondo effect, a correlation effect where the spin of an impurity is completely screened by conduction electrons of the host material at low temperature, is one of the most fundamental problems in the field of strongly correlated systems Hewson (1997). Since the theoretical work of Kondo studying the anomalous resistivity of metallic systems doped with magnetic impurities Kondo (1964) in the 1960s, the Kondo model has revealed interesting universal behaviors at temperatures below a Kondo energy scale $T_{K}$ related to the suppression of the impurity local moment Wilson (1975). The Kondo model and the more general Anderson model Anderson (1961) are considered as the most basic models to study many-body physics, which is rather well-understood especially with the invention of the numerical renormalization group (NRG) method in 1970s Wilson (1975); Bulla et al. (2008). The Kondo physics plays an important role in understanding more complicated strongly correlated problems such as the Kondo lattice model for heavy-fermion systems Doniach (1977); Löhneysen et al. (2007); Gegenwart et al. (2008), lattice correlated systems such as the Hubbard model or multi-orbital systems Georges and Kotliar (1992); Georges et al. (1996, 2013). From the methodology point of view, solving the correlated impurity problem allows scientists to further develop more powerful methods for correlated lattice systems, in particular the dynamical mean-field theory Georges et al. (1996); Gull et al. (2011).

Since 1990s, the topic of Kondo problems has been “revived” Kouwenhoven and Glazman (2001) thanks to new experimental realizations of Kondo systems in quantum dots Goldhaber-Gordon et al. (1998); Sasaki et al. (2000) and by controlling magnetic adatoms on metallic surfaces using scanning tunneling microscopy Crommie et al. (1993a, b); Li et al. (1998); Manoharan et al. (2000). Moreover, there are corners in the field of Kondo problems that are not thoroughly investigated and are still challenging to solve theoretically or to measure experimentally. One of them is the class of impurity problems where the host materials are semimetals, of which the low-energy density of states (DOS) vanishes at the Fermi level as a power law of frequency $\nu(\omega)\sim|\omega|^{r}$. It is called the “pseudogap” Kondo problem Withoff and Fradkin (1990); Gonzalez-Buxton and Ingersent (1998). The vanishing point at the Fermi level of the DOS of the host material may prevent the Kondo screening effect to occur even at nonzero coupling constant, thus there exists a quantum phase transition from the free local moment phase to the Kondo phase, where this local moment is fully screend, by varying the coupling constant or by increasing valence mixing Withoff and Fradkin (1990); Pixley et al. (2012). This is different from the conventional Kondo system (the DOS of the host is nonzero at the Fermi level), where the state of the free local moment only appears when the Kondo coupling vanishes. Previous works Withoff and Fradkin (1990); Bulla et al. (1997); Gonzalez-Buxton and Ingersent (1998) have shown that different quantum phases occur depending on the value of $r$, which is divided into three regions: small $r$ ($0<r<r^{*}$), intermediate $r$ ($r^{*}<r<r_{max}$) and large $r$ ($r>r_{max}$) with $r^{*}\approx 0.375$ and $r_{max}=\frac{1}{2}$. Small- and intermediate-$r$ regions are distinguished from the large-$r$ one by the quantum critical point (QCP) at the particle-hole symmetric point, which is absent in the latter case. The intermediate-$r$ region is differentiated from the small-$r$ region by an additional asymmetric fixed point when the particle-hole symmetry is broken Gonzalez-Buxton and Ingersent (1998). In this region, $r=1$ is the case of interest. Cassanello and Fradkin Cassanello and Fradkin (1996) has shown that, at $r=1$, there are logarithmic corrections included in the scaling laws of many kinds of physical observables, suggesting that the physics at $r=1$ be analogous to a system at an upper critical dimension at criticality. Further study Vojta and Fritz (2004); Fritz and Vojta (2004) provides more evidence for the characteristics of upper critical dimension at $r=1$ such as the nature of the phase transition is changed when crossing $r=1$ and perturbative approach can be feasible at $r>1$.

In the experimental aspect, the pseudogap Kondo model has been used to interpret reasonably a new type of QCPs that occurs in heavy fermion systems Si et al. (2001); Gegenwart et al. (2008); Pixley et al. (2012) or to study impurity local moment in unconventional superconductors Vojta and Bulla (2001). However, the difficulty in realizing the model in experiments is the main issue that hinders the progress of research in this topic. Indeed there are not many materials with power-law DOS around the Fermi level, a few candidates are $d$-wave superconductors and graphene which are related to the case $r=1$. Studying impurities in $d$-wave superconductors is nontrivial because directly measuring properties at the impurity in the bulk system is not simple. Graphene on the other hand might be more prospective experimentally as it is a 2D material, allowing for surface measurement techniques such as scanning tunneling microscopy. Not so long after its discovery Novoselov et al. (2004), graphene has been used as the host material for magnetic adatoms, typically the cobalt adatoms Mattos (2009); Jacob et al. (2009); Brar et al. (2011); Wang et al. (2012). Much effort is devoted to the determination of the stable position of an adatom on graphene and the dependence of its magnetic state on each position Wehling et al. (2010a, b); Eelbo et al. (2013); Virgus et al. (2014). Furthermore, for most of the cases, graphene is deposited on a substrate, the substrate is thus another factor that strongly affects the magnetic behavior of the adatom on graphene Ren et al. (2014); Donati et al. (2014). The readers may also refer to Ref. Fritz and Vojta (2013) for a review of theories and experiments for the problem of an impurity on graphene lattice.

This paper is devoted to the understanding of the pseudogap Kondo problem where the host metal, a honeycomb lattice, simulating the low-energy physics of graphene, and the adatom is restricted to an $s$-wave single-orbital impurity. From our high-temperature simulation data generated by continuous-time quantum Monte Carlo method Gull et al. (2011), we employ Binder analysis Binder (1981); Pixley et al. (2012) to determine the phase transition at zero temperature and study physical properties as the temperature decreases. We focus mostly on the case of an impurity placed on top of a lattice site. First we build a full phase diagram for the $r=1$ pseudogap Anderson model, in which the Kondo screening phase can only be observed in a mixed-valent state by driving the impurity away from half-filling. We then discuss briefly the situations where the impurity is placed at different location on the lattice. Interpretation of the phase diagram and discussion of different impurity positions on the lattice give us a broader view about the possibility for the realization of the pseudogap Kondo model on graphene.

The rest of the paper is structured as follows. Section II describes the model and basic parameters of the system under investigation, as well as the simulation method in use (mainly the continuous-time quantum Monte Carlo method Gull et al. (2011)). Section III presents the full phase diagram of the pseudogap Anderson model when the impurity is placed on top of an atom in the honeycomb lattice. Section IV analyzes the behaviors of each phase boundary presented in Sec. III. Section V discusses the situations where the impurity is placed at different locations in the host material. Finally, Section VI concludes our study.

The main part of our study is the investigation for the case of an impurity on top of an atom in the honeycomb lattice (the $t$ site), which will be presented thoroughly in this and the next sections. We start by constructing the full phase diagram for this case. If considering the nearest neighbor hopping $t$ of conduction electrons as the basic energy scale and neglecting further hopping or hybridization, the phase diagram is a three-dimensional one with respect to the impurity energy level $\epsilon_{d}$, the Hubbard interaction strength $U$ and the nearest neighbor hybridization $v_{0}$.

Assuming the impurity is on top of an $A$ site in the honeycomb lattice, the hybridization function for the impurity is calculated using Eq. (8) with $V_{Ak}=v_{0}$ and $V_{Bk}=0$. It is indeed a replication of the corresponding honeycomb lattice Green function, differed only by a prefactor $\frac{v_{0}^{2}}{2}$. Thus the hybridization spectrum is related to the DOS of the lattice as

As we later examine in Fig. 8 the hybridization function for three impurity positions, the impurity at the $t$ site is the only case which exhibits symmetric hybridization spectrum. Thus it is the closest to the low-energy $r=1$-pseudogap Kondo model Gonzalez-Buxton and Ingersent (1998); Pixley et al. (2012). Therefore, in order to avoid unnecessary complication, the case of an impurity on top of a lattice site is under full consideration in this work.

We present in Figure 3 the phase diagrams in the $(v_{0},\epsilon_{d})$ plane [panels (a) and (b)] and in the $(U,\epsilon_{d})$ plane [panel (c)] where the remaining parameter is kept constant. As noted previously, because these are quantum phase transitions, it requires low temperature for investigation. For nonzero-temperature methods such as CT-HYB, we can employ Binder analysis (see Sec. II.2) to extrapolate critical points from high-temperature calculations Binder (1981); Glossop et al. (2011); Pixley et al. (2012). All the phase boundaries in Fig. 3 are determined in this way. To justify the Binder analysis approach, we present in Appendix B additional measurements of the occupancy and the spin susceptibility for a specific case to confirm the existence of the critical points.

For the pseudogap model where the DOS of the host material is $\nu(\omega)\sim|\omega|^{r}$ with $r>r_{max}=1/2$, there is no phase transition at the particle-hole symmetric point $\epsilon_{d}=-\frac{U}{2}$, regardless of how large $U$ is Gonzalez-Buxton and Ingersent (1998). It means that at half-filling, the impurity local moment is always decoupled from the conduction electrons, i.e. the free local moment phase. Only when $\epsilon_{d}$ is shifted away from the particle-hole symmetry, a mixed-valence quantum critical point may occur Pixley et al. (2012) for the transition to the phase of screened local moment (the Kondo phase), as illustrated in Fig. 3. In all phase diagrams, the line $\epsilon_{d}=-\frac{U}{2}$ divides the phase diagram into two regions above and below half-filling, each of which contains only one phase boundary line separating the two phases. For the region of less than half-filling ($\epsilon_{d}<-\frac{U}{2}$), the left side of the phase boundary is the Kondo screening phase, while the right side is the free local moment phase. Because of the particle-hole symmetry at $\epsilon_{d}=-\frac{U}{2}$, the phase diagram is symmetric about this line. Thus, in Fig. 3(b) and (c), we plot the region of $\epsilon_{d}>-\frac{U}{2}$ as the mirror image of the other half-plane.

Focusing only on the region of less than half-filling $\epsilon_{d}<-\frac{U}{2}$, we find that the critical value $\epsilon_{d}^{c}$ is always in the range $[-U,-\frac{U}{2}]$ and it depends monotonically on $v_{0}$. In Fig. 3(a) and (b), $\epsilon_{d}^{c}$ approaches $-U$ as $v_{0}\to 0$ and tends to approach $-\frac{U}{2}$ as $v_{0}\to\infty$. This dependence has been proposed elsewhere Chowdhury and Ingersent (2015) but no quantitative calculation has been conducted until now. To understand this, one recalls that in conventional Kondo models, the Kondo scale depends strongly on the hybridization function around the Fermi level. Because the resonance width depends on the hybridization spectrum at the Fermi level $\Delta_{r}\sim-\frac{1}{\pi}\mathrm{Im}\Delta(\omega=0)$, the Kondo physics for symmetric hybridization spectrum can be observed if $-U-\Delta_{r}<\epsilon_{d}<\Delta_{r}$ Krishna-murthy et al. (1980a); Hewson (1997). For the pseudogap model, $\Delta_{r}=0$, the necessary condition to observe the Kondo physics is $-U<\epsilon_{d}<0$. However, because there is no phase transition at half-filling for $r=1$ pseudogap model, the line $\epsilon_{d}=-\frac{U}{2}$ is the upper limit for all phase boundary curves in this region. Besides, the local moment screening effect is enhanced by the hybridization strength $v_{0}$ and the valence mixing controlled by $\epsilon_{d}$. Therefore, at large $v_{0}$, a small amount of valence mixing is enough to screen the impurity local moment, which means $\epsilon_{d}^{c}$ is closer to $-\frac{U}{2}$ as $v_{0}$ increases.

Interestingly, we find that in Fig. 3(c) the phase boundary plotted in the $(U,\epsilon_{d})$ plane exhibits linear behavior

where $\alpha_{t}$ is defined as the slope of the phase boundary, which depends strongly on $v_{0}$. We will study several aspects of this linear feature. First, we examine the dependence of the $(U^{c},\epsilon_{d}^{c})$ line on $v_{0}$. Second, we compare it with that of low-energy pseudogap impurity models which are often used in literature, i.e. models where the DOS of the host is $\nu(\omega)\sim|\omega|^{r}$ Withoff and Fradkin (1990); Gonzalez-Buxton and Ingersent (1998).

In Fig. 3(a), the three phase boundaries are equally distanced as $U$ increases but the differences between the three $\epsilon_{d}^{c}$ decreases as $v_{0}$ increases, implying the dependence of $\alpha_{t}$ on $v_{0}$. Based on the data for Fig. 3(a) and the fact that the phase boundary line goes towards the origin of the $(U,\epsilon_{d})$ plane, we plot in Figure 4 the phase boundaries for various $v_{0}$ values. The phase boundaries approach $\epsilon_{d}=-U$ (the black dashed line) illustrating the lower limit described above. They are all below the upper limit $\epsilon_{d}=-\frac{U}{2}$ (the black dotted line). However it requires larger $v_{0}$ to observe the approach to the upper limit.

Because the Kondo effect is a low-energy physical phenomenon, we expect that even in the honeycomb lattice (graphene), the physics is only affected by changes in proximity to the Fermi level. In our case, the only change in low energy ($\omega\sim 0$) is the slope of the hybridization spectrum. Around the Fermi level, the tight-binding DOS of graphene is linear $\nu(\omega)=\frac{A_{c}}{\pi v_{F}^{2}}|\omega|$ (where $A_{c}=\frac{3\sqrt{3}}{2}$ is the area of the unit cell and $v_{F}=\frac{3t}{2}$ is the Fermi velocity) Castro Neto et al. (2009). Employing Eq. (10), the slope $\gamma_{t}$ of the hybridization spectrum (a dimensionless parameter) depends only on $v_{0}^{2}$

We verify our statement by comparing our results for honeycomb lattice with those from standard pseudogap problem, in which the hybridization spectrum is set to have the same slope

Here the step function $\theta(x)$ limits the spectra to be nonzero only in the range $(-D,D)$, defining the bandwidth $2D$. For simplicity, in this test we set $D=t$. We test with $v_{0}=0.5t,0.7t$ and $1.4t$. In Fig. 4 the phase boundary of the conventional model corresponding to $v_{0}=0.5t$ (the red dotted line) is plotted for comparison. This conventional pseudogap phase boundary stays rather closely to the $v_{0}=0.5t$ phase boundary of the honeycomb lattice model, the difference between the two phase boundary slopes is $4\%$. For other cases corresponding to $v_{0}=0.7t$ and $1.4t$ (not shown), the differences are $6\%$ and $5\%$, respectively. In addition to the error in CT-HYB measurements, the difference of the phase boundaries between the honeycomb lattice model and the simplifed pseudogap model might be due to the higher energy effect that may involve, resulting in the contribution of nonlinear higher order terms in the tight-binding honeycomb lattice DOS. Nevertheless, the difference is small, thus the similarity between the two models provides evidence that, because the Kondo effect is a low-energy effect, the physics near the Fermi level plays a crucial role. In this case, it is represented by the slope of the hybridization spectra, which depends mostly on $\frac{v_{0}}{t}$.

In the inset of Fig. 4, the dependence of the slope $\alpha_{t}$ of the phase boundary on $\left(\frac{v_{0}}{t}\right)^{2}$ is plotted. Interestingly, this slope depends linearly on $v_{0}^{2}$ for small $v_{0}$ ($v_{0}<t$) and approaches the lower limit value $1$ as $v_{0}\to 0$ (corresponding to $\epsilon_{d}^{c}=-U$). For $v_{0}>t$, $\alpha_{t}$ slowly approach the upper limit value $2$ (the upper limit at half-filling $\epsilon_{d}=-\frac{U}{2}$) as $v_{0}$ increases. Thus $v_{0}$ is the key factor that controls the slope $\alpha_{t}$ of the $(U^{c},\epsilon_{d}^{c})$ phase boundary. At different $v_{0}$ scales, $\alpha_{t}$ behaves differently.

We will discuss in details the behavior of each type of phase boundary in the phase diagram presented in Fig. 3. In these phase boundries, there are two interesting features for further analysis: (1) the linear relation between $U^{c}$ and $\epsilon_{d}^{c}$ for fixed $v_{0}$, and (2) the dependence of the slope $\alpha_{t}=-U^{c}/\epsilon_{d}^{c}$ as well as $\epsilon_{d}^{c}$ itself on $v_{0}$. The impurity self-energy is used intensively to understand these features.

So far, we focus mostly on the situation where the impurity is placed on top of an atom on the honeycomb lattice. This is a special case in which the hybridization spectrum is symmetric around the Fermi level and there exists particle-hole symmetry at $\epsilon_{d}=-\frac{U}{2}$. In general, the hybridization spectrum depends on the impurity position in the lattice. In Appendix A, we show that when the impurity departs from this $t$ site, it hybridizes with atoms of both sublattices, the hybridization spectrum is no longer symetric around the Fermi level. When the hybridization spectrum changes, it obviously affects the physics. In this section, we discuss briefly the effects when the impurity is placed at the hollow $h$ and at the bridge $b$ sites [see Fig. 1].

Figure 8 shows the real part and the spectrum $-\frac{1}{\pi}\mathrm{Im}\Delta(\omega)$ of the hybridization function when the impurity is placed at $t$, $h$ and $b$ positions for the same $v_{0}$. Because the hybridization function is directly proportional to $v_{0}^{2}$, for simplicity, we fix $v_{0}/t=1$ in the figure. For different $v_{0}$, one simply multiplies the hybridization function with $v_{0}^{2}$ to obtain a new one. Fig. 8 shows that $\mathrm{Re}\Delta(\omega=0)$ increases as the impurity changes from $t\to b\to h$. Due to contributions from both sublattices (see Appendix A), the imaginary parts obviously show the asymmetry of the hybridization spectra when the impurity is not at the $t$ site. We have found that while the hybridization spectra for the $t$ and $b$ sites remain linear around the Fermi level, it behaves as $|\omega|^{3}$ for the $h$ site case. Expansion around the Fermi level shows the forms of $\Delta(\omega)$ around $\omega=0$, of which the results are summarized in Table 1.

Nonzero $\mathrm{Re}\Delta(\omega=0)$ values of impurity at the $b$ and $h$ sites result in the shifting of the critical impurity level $\epsilon_{d}^{c}$. By defining an effective energy level $\epsilon_{eff}=\epsilon_{d}+\mathrm{Re}\Delta(\omega=0)$, the critical value $\epsilon_{eff}^{c}$ is in the same energy range as $\epsilon_{d}^{c}$ for impurity at the $t$ site $(-U<\epsilon_{eff}^{c}<-\frac{U}{2})$, allowing to determine $\epsilon_{d}^{c}$ in these cases more easily. The difference in the real part of $\Delta(\omega)$, of which the main effect is to introduce a shifting in $\epsilon_{d}^{c}$, is thus not as important as the change in the imaginary part.

When placing the impurity at the $b$ site, the hybridization function exhibits the asymmetric spectrum, starting from the the second order of the frequency

Examining the $(U,\epsilon_{d})$ phase diagram at $v_{0}=0.5t$ for the tight-binding honeycomb lattice Eq. (8) and for the simplified model in Eq. (23), the phase boundaries (not shown) still exhibit linear relation. The difference of the boundary slope between the lattice model and the simplified model is about $2\%$, which we argue that it is due to the differences in the hybridization function at frequencies away from $\omega=0$. Apart from that, we do not observe significant changes when changing the location of the impurity from the $t$ site to the $b$ site. Thus, although the hybridization spectrum is asymmetric, as long as the spectrum still exhibits linear behavior around the Fermi level, there is no significant change in physics.

Impurity at the $h$ site is a different scenario. In this case, the impurity strongly hybridizes with conduction electrons of six nearest neighbor sites. We find that $|\omega|$ and $\omega^{2}$ terms in the hybridization spectrum vanish, leaving the $|\omega|^{3}$ term as the lowest order term in the spectrum (see Table 1). The impurity problem is transformed into an $r=3$ pseudogap Kondo model. Due to very small hybridization spectral weight around the Fermi level, it turns out to be not easy to observe the Kondo screening phase or determine the critical point with high accuracy in this case. Indeed, our numerical calculations contain large error bars, which prohibit us from reliable extrapolation for $\epsilon_{d}^{c}$.

Therefore, the physics of the pseudogap Kondo model of an impurity on the honeycomb lattice depends on the location of the impurity on the lattice. If it is placed at $b$ and $t$ sites, linear behavior of the host DOS is retained in the hybridization spectrum, allowing for the observation of the Kondo screening phase. However, around the $h$ site, the hybridization spectrum near the Fermi level is nearly depleted and becomes a different type of pseudogap model. At this impurity position, the conditions turn out to be more extreme in order to observe the Kondo physics.

In this paper, we have studied the pseudogap Anderson model for an $s$-wave impurity placed on a graphene-like honeycomb lattice. We used the tight binding model to describe conduction electrons of the host material, and employed the hybridization-expansion continuous-time quantum Monte Carlo method (CT-HYB) to fully treat the correlation effect. We focused on the case the impurity is placed on top of an atom in the lattice (the $t$ position) and constructed the phase diagram with three parameters $U,v_{0}$ and $\epsilon_{d}$ for the quantum phase transition from the free local moment phase to the Kondo screening phase. We analyzed the impacts of the hybridization and the correlation strength on the phase diagram and provided insightful discussion for understanding of the phase diagram based on the self-energy of the impurity. We also discussed qualitatively the physics and the possibility to observe the Kondo screening phase when the impurity position is moved away from $t$ position.

The work presents calculations from high-temperature viewpoint and contains several interesting results. First, we have found that the full phase diagram exhibits interesting relations between each pair of $(U^{c},\epsilon_{d}^{c})$ and $(v_{0},\epsilon_{d}^{c})$. The $(U^{c},\epsilon_{d}^{c})$ phase boundary is linear for the range of $v_{0}$ under investigation; while in the $(v_{0},\epsilon_{d}^{c})$ plane, $\epsilon_{d}^{c}$ exhibits $v_{0}^{2}$ (or $\frac{1}{v_{0}^{2}}$) dependence in the limit of small (large) $v_{0}$. We have discussed the results by relating the critical value $\epsilon_{d}^{c}$ to the impurity self-energy $\Sigma(i\omega_{n})$ and obtained an empirical expression for the self-energy in order to gain more insights into the physics of the phase diagram. Lastly, we have found that changing the impurity from position on top of an atom to the bridge ($b$) or hollow ($h$) site might affect the physics. Specically, while the system remains in the $r=1$ pseudogap Anderson model if the impurity is at the $t$ or $b$ sites, it becomes an $r=3$ pseudogap model if the impurity is placed at the $h$ site.

In overall, the realization of the $r=1$ pseudogap model in graphene is still a challenging task for both experimental and theoretical sides. Indeed, the transition from the free local moment to Kondo screening phase depends on the position of the impurity on the honeycomb lattice. The prominent case to observe the Kondo physics is when the impurity is located at the $t$ site, where the hybridization spectrum is symmetric and directly proportional to the DOS of the honeycomb lattice. If it is located at the $h$ site, it becomes a different pseudogap model and is not easy to observe the Kondo screening effect. Although the hybridization at the hollow site for a more realistic $d$-wave impurity might be slightly different Wehling et al. (2010a, b), the small slope of the hybridization spectrum around the Fermi level remains and is an obstacle for any measurement of the Kondo screening phase.

Our work is restricted to the CT-HYB impurity solver, which is a high-temperature method, with an assumption of a single $s$-wave impurity for the sake of simplicity. Therefore, certain limitations in the work are inevitable. The major one is that, due to the CT-HYB impurity solver, we cannot access low temperature region in order to enter the quantum critical regime and measure critical behaviors more accurately. There is also a lack of analytic calculations for asymtotic limits, such as at $v_{0}\to 0$ and $v_{0}\to\infty$, in order to verify numerical results. Moreover, the logarithmic corrections to scaling Cassanello and Fradkin (1996) for $r=1$-pseudogap model is not considered in this work. However, these limitations open several interesting directions for future investigation. First, one might consider low-temperature approach, in particular the numerical renormalization group method Bulla et al. (2008), to treat the problem completely. Second, one might employ previous knowledge of perturbation theories Dai et al. (2005); Yosida and Yamada (1975); Horvatić and Zlatić (1980); Tong (2015) for analytic results in the asymtotic limits, which may justify our numerical results. Other questions of how the logarithmic corrections to scaling enter the numerical results, how the results change when using a more realistic $d$-orbital magnetic impurity, or about the effect of substrates to the graphene layer and the impurity Ren et al. (2014); Donati et al. (2014) are also interesting open challenges.

We thank T. A. Costi for helpful discussion. This research is funded by Vietnam National Foundation for Science and Technology Development (NAFOSTED) under grant number 103.01-2018.12. We also acknowledge the support for the allocation of computing time at Jülich Supercomputing Centre. Portion of our calculation as well as data post-processing tasks has been performed using the computer cluster of Phenikaa Institute for Advanced Study.