\section{Introduction}

The possibility of achieving of fault-tolerant quantum computing with qubits based on Majorana bound states (MBSs)~\cite{ReviewAguado,kitaev} started a new era in the domains of mesoscopic physics and quantum information. These exotic non-Abelian excitations~\cite{RevNonabelian} emerge as topologically protected mid-gap zero-energy modes in so-called topological superconductors~\cite{Alicea,RevMajorana}. The topological protection stems from the separation between individual MBSs, i.e, \textit{nonlocality}, which is also responsible for the immunity of a setup against local perturbations and consequent loss of the information due to the processes of decoherence~\cite{ReviewDasSarma,RevNonabelian}. However, it should be noticed that, for practical realizations of quantum computing systems, the MBSs qubit becomes vulnerable to the decoherence process caused by local perturbations when coupled to environment~\cite{Decaymemories,Failureprotection}, which can lead to unwelcome errors in the processing of quantum information.    

Topological superconductivity can be realized experimentally in hybrid superconductor-semiconductor nanowires with induced proximity effect in the presence of strong spin-orbit coupling and external magnetic field, favoring the formation of superconducting (SC) triplet states ~\cite{Nanowire1,Nanowire2}. In these hybrid devices, manifestation of a robust zero-bias conductance peak (ZBCP) has been considered as an experimental signature of the presence of highly nonlocal MBSs emerging at the opposite ends of a nanowire \cite{mourik2012,albrecht2016exponential,Deng2016,AguadoExperiment,Quantized}. However, it was argued later on that other physical mechanisms such as disorder~\cite{Disorder}, Kondo effect~\cite{Kondo1,Kondo2} and formation of Andreev bound states (ABSs)~\cite{ABS1,ABSLee,ABS2,psABS1,psABS2,Cayao2018,YSR,SNSjunction,Josepheson,SanJose2013} can be responsible for the appearance of analogs of ZBCP. In particular, there is ongoing controversy~\cite{nonHermitan} whether near zero-energy ABS, constituted by weakly overlapping MBSs, can mimic robust $2e^{2}/h$ ZBCP ~\cite{psABS1,psABS2}. Overall, there is currently consensus that observation of ZBCP only is not enough to guarantee the presence of topologically protected MBSs in the system.
\begin{figure}[t]
	\centerline{\includegraphics[width=3.5in,keepaspectratio]{hybridNW2.pdf}}
	\caption{Sketch of the system consisting of a hybrid superconducting (SC) nanowire (blue region) coupled to a quantum dot (QD) with energy level $\varepsilon_{d}$, which can be tuned by application of an external gate voltage $V_{\text{Dot}}$. The QD is coupled to both Majorana bound states (MBSs) $\gamma_{L}$ and $\gamma_{R}$ at the opposite ends of SC nanowire with strengths $\lambda_{L}$ and $\lambda_{R}$, respectively. The MBSs may be hybridized with each other by $\delta_{M}$ in the  presence of an external magnetic field applied longitudinally (purple arrow) due to the finite size effects. The QD levels are broadened due to the coupling $\Gamma$ with a normal metallic lead N.
		\label{Setup}}
\end{figure}

A possible way to clarify the origin of ZBCP is performing tunneling spectroscopy of a quantum dot (QD)~\cite{Hoffman,Denis} is assumed to be coupled to both ends of SC nanowire. In this type of experiment one can access the so-called degree of Majorana nonlocality~\cite{RAguado,DJClarcke,AguadoExperiment} characterizing \textit{``how topological''} are MBSs and distinguish between the cases of well-separated MBS and near zero-energy ABSs (overlapping MBSs)~\cite{Deng2016,ABS2,ABS3,ABS4}. 

Distinct from earlier works ~\cite{Deng2016,RAguado,DJClarcke,AguadoExperiment,MBSoscillations}, in the present paper we analyze how charging energy of the QD single-level coupled to a normal lead affects the energy spectrum of the device sketched in the Fig.~\ref{Setup}. To account for the correlation effects we go beyond the Hartree-Fock mean-field approximation used by Prada \textit{et al}~\cite{RAguado} for larger Zeeman fields, applying the method developed by Hubbard~\cite{HubbardI} to treat the charging energy which allows us to recover qualitatively the recent experimental profiles reported by Deng \textit{et al}~\cite{AguadoExperiment} for a set of tunable parameters. In this sense, we demonstrate that for highly nonlocal MBSs, a plateau at zero-energy is formed in the QD density of states in the wide range of the values of the dot level and charging energies. In the case of strongly overlapping MBSs forming an ABS, the spectrum is strongly modified and this plateau disappears. Moreover, changes in MBSs degree of nonlocality strongly affect spin resolved density of states of QD, which means that spectroscopic experiment with spin-polarized local probe~\cite{Bena,Bjornson,Tadeusz,Jeon,Andrei,Jelenaspinpolarized} will allow to identify whether ZBCP is induced by topological MBSs or ABSs (overlapping Majoranas).
 \section{The Model}
 
 The effective Hamiltonian describing the device depicted in Fig.~\ref{Setup} takes the following form \cite{RAguado,MBSoscillations}
\begin{eqnarray}
H&=&\sum_{\mathbf{k}\sigma}\varepsilon_{\mathbf{k}\sigma}e_{\mathbf{k}\sigma}^{\dagger}e_{\mathbf{k}\sigma}+\sum_{\sigma}\varepsilon_{d\sigma}n_{d\sigma}+Un_{d\uparrow}n_{d\downarrow} \nonumber \\
& + & \sqrt{2}V\sum_{\mathbf{k}\sigma}(e_{\mathbf{k}\sigma}^{\dagger}d_{\sigma}+d_{\sigma}^{\dagger}e_{\mathbf{k}\sigma}) + H_{Nw},\label{eq:FullH}
\end{eqnarray}
where $n_{d\sigma}=d_{\sigma}^{\dagger}d_{\sigma}$ is the operator of the number of the electrons residing in the single-level QD, the operators  $e_{\mathbf{k}\sigma}^{\dagger}(e_{\mathbf{k}\sigma})$ correspond to the creation of the electrons (holes) in the normal (N) lead~\cite{MBSoscillations} with wave vector $\mathbf{k}$, spin $\sigma=\uparrow,\downarrow$ and energy $\varepsilon_{\mathbf{k}\sigma}$. The energy of an electron in the QD is spin dependent, $\varepsilon_{d\sigma}=\varepsilon_{d}-\sigma V_{Z}$, where
$V_{Z}$ is Zeeman energy splitting induced by an external magnetic field. $U$ is the charging energy of the QD. The dot is coupled to  the normal lead N with coupling strength $\sqrt{2}V$.

To describe the SC-nanowire hosting a pair of MBSs $\gamma_{i}$ at the opposite ends and coupled to the QD, we use low-energy effective model developed by Prada \textit{et. al}~\cite{RAguado}, and characterized by the following Hamiltonian:
\begin{eqnarray}
H_{Nw}&=&\imath\delta_{M}\gamma_{L}\gamma_{R}+(\lambda_{L}d_{\sigma}-\lambda_{L}^{*}d_{\sigma}^{\dagger})\gamma_{L}\nonumber \\ &+& (\lambda_{R}d_{\sigma} + \lambda_{R}^{*}d_{\sigma}^{\dagger})\gamma_{R}\label{eq:HNw}
\end{eqnarray}
where self-conjugated  operators $\gamma_{i}=\gamma_{i}^{\dagger}$ describe localized Majorana fermions and obey the algebraic relation $\{ \gamma_{i}, \gamma_{j} \}  = \delta_{ij}$~\cite{ReviewAguado,Alicea,ReviewDasSarma,RevMajorana}. $H_{Nw}$ can be rewritten in the regular fermionic basis, since Majorana operators can be decomposed into pairs of normal fermionic operators, $\gamma_{L}=(f_{\uparrow} + f_{\uparrow}^{\dagger})/\sqrt{2}$  and $\gamma_{R}=\imath(f_{\uparrow}^{\dagger}-f_{\uparrow})/\sqrt{2}$. The dot is coupled to the left and right MBSs, with coupling constants $\lambda_{L}$ and $\lambda_{R}$, respectively. The direct  hybridization $\delta_{M}$ between MBSs reads~\cite{Danon,Fleckenstein}
\begin{equation}
\delta_{M}=\frac{e^{-{l/2b}}}{\sqrt{b}}\cos(l\sqrt{b})E_{0}, \label{eq:deltaM}
\end{equation}
which is the function of both Zeeman energy splitting $b=V_{Z}/E_{0}$, $E_{0}=\left(2m^{*}\alpha^{2}\Delta_{\text{SC}}^{2}/\hbar^{2}\right)^{1/3}$ and $l=L\sqrt{2m^{*}E_{0}}/\hbar$ with $L$ being the length of the wire, $m^{*}$  being electrons effective mass, $\alpha$ the spin-orbit coupling constant and $\Delta_{\text{SC}}$ the induced SC gap~\cite{Danon}. 
The degree of MBSs nonlocality $\eta$ can be defined as ratio between QD-MBSs right/left coupling strengths~\cite{RAguado}:
\begin{equation}
\eta^{2} = \frac{|\lambda_{R}|}{|\lambda_{L}|}
\end{equation}
This parameter can be experimentally accessed through the measurement of the conductance as a function of the gate potential changing the energy of a QD and  drain-source voltage~\cite{AguadoExperiment} and estimated as the ratio between energy values in which the Majorana and QD states are on resonance (anticrossing points)~\cite{RAguado}, $\eta^2 \approx \epsilon_{\text{MBS}}^{\pm}/\epsilon_{\text{QD}}^{\pm}$ [See Fig~\ref{Fig2}(d-f)]. 
\subsection{Density of states calculations} 

Our main goal is to investigate how the spectral properties of the QD accessible in spin-resolved measurements are changing when the degree of MBSs nonlocality characterized by the parameter $\eta$~\cite{RAguado} is modified. Hence, it is suitable to evaluate the total density of states (DOS) in the QD, which reads:
\begin{equation}
\text{DOS}(\omega) = \pi\Gamma \sum_{\sigma}\rho_{\sigma}(\omega),\label{eq:DOS}
\end{equation}
where the constant $\Gamma = 2\pi V^{2} \rho_{0}$ is the QD-N lead effective coupling~\cite{Anderson}, with $\rho_{0}$ being the DOS of the lead. The quantity
\begin{equation}
 \rho_{\sigma}(\omega)=-\frac{1}{\pi}\text{Im}[G^{r,\sigma}_{d,d}(\omega)]   
\end{equation}
denotes the DOS corresponding to a given spin orientation, which is determined by the retarded Green's function of the QD $G^{r,\sigma}_{d,d}(\omega)$ in the spectral domain. The application of the equation of motion (EOM) method~\cite{Jauho} leads to the following equation \textcolor{blue}{(Appendix \ref{AppendixA}):}
\begin{eqnarray}
(\omega^{+}-\varepsilon_{d\sigma}-\Sigma_{\text{M},\sigma}^{\text{U=0}}+\imath\Gamma)G_{d,d}^{r,\sigma}(\omega)=1+UG_{d_{\sigma}n_{d\bar{\sigma}},d_{\sigma}}^{r}(\omega)\nonumber\\+U(|\lambda_{L}|^{2}-|\lambda_{R}|^{2})\bar{K}^{\sigma}G_{d_{\sigma}^{\dagger}n_{d\bar{\sigma}},d_{\sigma}}^{r}(\omega),\nonumber\\
\label{eq:E1}
\end{eqnarray}
where $\omega^+=\omega+i0^{+}$, $\Sigma_{\text{M},\sigma}^{\text{U=0}}=K_{1} +(|\lambda_{L}|^{2}-|\lambda_{R}|^{2})^{2}K\bar{K}^{\sigma}$ is the self-energy~\cite{MBSoscillations} due to QD-MBSs hybridization in the absence of the charging energy, and 
\begin{equation}
K=\frac{1}{2}\left(\frac{1}{\omega^{+}+\delta_{M}}+\frac{1}{\omega^{+}-\delta_{M}}\right)\label{eq:K},    
\end{equation}

\begin{equation}
\bar{K}^{\sigma} = \frac{1}{2} \left(\frac{K}{\omega^{+} + \varepsilon_{d\sigma}-K_{2} + \imath \Gamma}\right),\label{eq:Kbar}
\end{equation}
\begin{equation}
K_{1} = \frac{1}{2}\cdot \left[\frac{(|\lambda_{L}|-|\lambda_{R}|)^{2}}{\omega^{+}-\delta_{M}} + \frac{(|\lambda_{L}|+|\lambda_{R}|)^{2}}{\omega^{+}+\delta_{M}}\right] \label{eq:K1}
\end{equation}
and
\begin{equation}
K_{2} = \frac{1}{2}\cdot\left[\frac{(|\lambda_{L}|-|\lambda_{R}|)^{2}}{\omega^{+}+\delta_{M}} + \frac{(|\lambda_{L}|+|\lambda_{R}|)^{2}}{\omega^{+}-\delta_{M}}\right]. \label{eq:K2}
\end{equation}
The presence of the two-particle operator corresponding to the charging energy term in the  Hamiltonian [Eq.~(\ref{eq:FullH})] leads to the appearance of the two particle Green's functions in the Eq.~(\ref{eq:E1}). The iterative application  of the EOM procedure to such higher order functions will produce an infinite chain of the equations which should be truncated at some point~\cite{Jauho}. Distinct from the earlier work~\cite{RAguado} in which the charging energy of the dot was accounted for using a mean-field approximation, here we take a step further by following Hubbard-I truncation scheme~\cite{HubbardI}.  This allows us to account for the appearance of the so-called Hubbard peaks and thus describe better the physics of Coulomb blockade regime. Note, however, that in our approach Kondo-type correlations are fully neglected, and it is applicable only if $T_{K}/\Delta_{\text{SC}}\gtrapprox0.6$~\cite{RAguado} or $T\gg T_{K}$, wherein $T_{K}$ is the Kondo temperature~\cite{Kondo1,Kondo2}. Further details of the calculations can be found in the Appendix \ref{AppendixB}.  

After Hubbard-I truncation, the two-particle Green's functions take the following form:
\begin{equation}
G_{d_{\sigma}n_{d\bar{\sigma}};d_{\sigma}}^{r}(\omega)=\frac{\langle n_{d\bar{\sigma}}\rangle}{\omega^{+}-\varepsilon_{d\sigma}-U-\Sigma_{\text{M},\sigma}^{\text{U\ensuremath{\neq}0}}+\imath \Gamma},\label{eq:E2}
\end{equation}
and
\begin{equation}
G_{d_{\sigma}^{\dagger}n_{d\bar{\sigma}}d_{\sigma}}^{r}(\omega)=-(|\lambda_{L}|^{2}-|\lambda_{R}|^{2})\bar{K}_{\text{U}}^{\sigma}G_{d_{\sigma}n_{d\bar{\sigma}},d_{\sigma}}^{r}(\omega),\label{eq:E3}
\end{equation}
wherein
\begin{equation}
\langle n_{d\bar{\sigma}}\rangle = \int_{-\infty}^{0} d\omega\rho_{\bar{\sigma}}(\omega) \label{eq:occupation}
\end{equation}
gives the occupation number of the dot per spin $\bar{\sigma}$ (opposite to $\sigma$) at $T=0$. The self-energy term provided by the presence of MBSs and charging energy $U$ is given by
\begin{equation}
\Sigma_{\text{M},\sigma}^{\text{U\ensuremath{\neq}0}}=K_{1}+(|\lambda_{L}|^{2}-|\lambda_{R}|^{2})^{2}K\bar{K}_{\text{U}}^{\sigma}    
\end{equation}
where
\begin{equation}
\bar{K}_{\text{U}}^{\sigma}=\frac{1}{2}\cdot\frac{K}{\omega^{+}+\varepsilon_{d\sigma}+U+\-K_{2} + \imath \Gamma}
\end{equation}

After some algebra we get from Eqs.~(\ref{eq:E1}), (\ref{eq:E2}) and (\ref{eq:E3}) the following expression for the retarded Green's function of the dot:
\begin{equation}
G_{d,d}^{r,\sigma}(\omega)=\frac{\lambda(\omega,\sigma\bar{\sigma})-U(|\lambda_{L}|^{2}-|\lambda_{R}|^{2})^{2}\mathcal{M}(\omega,\sigma\bar{\sigma})}{\omega^{+}-\varepsilon_{d\sigma}-\Sigma_{\text{M},\sigma}^{\text{U=0}} + \imath \Gamma}\label{eq:Gddfinal}
\end{equation}
with
\begin{equation}
\lambda(\omega,\sigma\bar{\sigma})=1+\frac{U\langle n_{d\bar{\sigma}}\rangle}{\omega^{+}-\varepsilon_{d\sigma}-U-\Sigma_{\text{M},\sigma}^{\text{U\ensuremath{\neq}0}} + \imath \Gamma}\label{eq:lambda}
\end{equation}
and
\begin{equation}
\mathcal{M}(\omega,\sigma\bar{\sigma})= \frac{\langle n_{d\bar{\sigma}}\rangle\bar{K}^{\sigma}\bar{K}_{\text{U}}^{\sigma}}{\omega^{+}-\varepsilon_{d\sigma}-U-\Sigma_{\text{M},\sigma}^{\text{U\ensuremath{\neq}0}}+\imath \Gamma}. \label{eq:M}
\end{equation}
 \section{Results and Discussion}
 
 We investigate the energy spectrum of the device depicted in Fig.~\ref{Setup} analyzing the DOS of the QD [Eq.~(\ref{eq:DOS})] as a function of spectral frequency $\omega$ and dot energy level $\varepsilon_{d}$ for several regimes corresponding to the different ratios between the parameters of the system. The relevant parameters of our model, in units of $E_{0}$, are the charging energy $U$, hybridization $\lambda_{L}(\lambda_{R})$ between the dot and MBS(left/right) and Zeeman energy splitting $V_{Z}$, which modulates the direct overlap $\delta_{M}$ between the MBSs at the opposite nanowire ends. The length of the SC was chosen as $L\approx0.1l\si{\micro\metre}$, in accordance with the results presented in the Ref.~\onlinecite{Danon}. The occupation numbers for each spin [Eq.~(\ref{eq:occupation})] were self-consistently computed. In all the situations, the QD-left MBS coupling strength is kept fixed ($\lambda_{L}=1.0E_{0}$). Concerning the charging energy strength $U$, we follow Prada \textit{et al.}~\cite{RAguado} effective Hamiltonian, assuming that $U>\Delta$.

Fig.~\ref{Fig2}(a) corresponds to the highest nonlocal situation ($\eta=0$), where SC nanowire is long enough ($L\approx\num{2.0}\si{\micro\meter}$) to ensure formation of isolated MBSs at the ends ($\delta_{M}=0$). In this regime, the dot couples only with the left MBS ($\lambda_{R}=0$).  As predicted by the earlier works~\cite{RAguado,DJClarcke,AguadoExperiment}, the Majorana states remain unperturbed under variations of the QD energy level, since the latter can not cross the topologically protected zero-energy MBSs. 

For shorter wires ($L\approx\num{0.4}\si{\micro\meter}$), Fig.~\ref{Fig2}(b), DOS of the dot reveals the so-called \textit{``bowtie''} profile, characteristic to the regime when the overlap between MBSs is finite and the dot is only weakly hybridized with rightmost Majorana ($\lambda_{L},\delta_{M} \gg \lambda_{R}$)~\cite{RAguado}. In this situation the topological protection is absent and the energies of the overlapping MBSs are strongly perturbed in the vicinity of the resonance with the QD state. The splitting of near-zero states is ruled by direct hybridization between MBSs $2\delta_{M}$ (See yellow bar in the panel (b)). 

Fig.~\ref{Fig2}(c) demonstrates the spectra for the case of the  local fermionic zero-mode ($\delta_{M}=0$ and $\eta = 1$), corresponding to the highest localization of MBSs (lowest topological quality factor~\cite{DJClarcke}), for which any pronounced structure at $\omega=0$ is absent. 

\begin{figure}[t]
	\centerline{\includegraphics[width=3.6in,keepaspectratio]{Result2.pdf}}
	\caption{Color scale plots of the DOS of a QD as a function of QD energy level $\varepsilon_{d}$ and spectral frequency $\omega$ for distinct regimes. In all situations we chose $U=5E_{0}$ and $\lambda_{L}=1.0E_{0}$. (a) highly nonlocal isolated MBSs ($\delta_{M}=\lambda_{R}=0$) for $V_{Z}=0.8E_{0}$ and $L=\SI{2.0}{\micro\meter}$; (b) bowtie profile, wherein $V_{Z}=1.72E_{0}$, $L=\SI{0.4}{\micro\meter}$, $\delta_{M}=0.12E_{0}$ and $\lambda_{R}=0.003E_{0}$. (c) regular fermionic zero mode ($\lambda_{L}=\lambda_{R}=E_{0}$, $\delta_{M}=0$), with $L=\SI{0.4}{\micro\meter}$ and $V_{Z}=1.38E_{0}$; (d)-(f) diamond profiles ($\delta_{M} \ll \lambda_{R},\lambda_{L}$), with actual nanowire length $L=\SI{0.4}{\micro\meter}$. Panel (d) exhibits the situation for  $V_{Z}=1.4E_{0}$, $\delta_{M}=0.004E_{0}$ and $\lambda_{R}=0.03E_{0}$. In (e) we have set $V_{Z}=1.5E_{0}$, $\delta_{M}=0.04E_{0}$ and $\lambda_{R}=0.25E_{0}$, while in panel (f) $V_{Z}=1.6E_{0}$, $\delta_{M}=0.08E_{0}$ and $\lambda_{R}=0.5E_{0}$. 
		\label{Fig2}}
\end{figure}

The panels (d)-(f) of the Fig.~{\ref{Fig2}} correspond to the case of the shorter SC nanowires ($L\approx\num{0.4}\si{\micro\meter}$), for which $\delta_M\neq0$ but $\delta_{M} \ll \lambda_{R}$.  This regime corresponds to the situation wherein the wave function describing the right MBS moves towards the QD due to the application of the magnetic field \cite{Quantized,psABS2}.
One can notice the presence of the previously reported \textit{``diamond''} profiles~\cite{RAguado,AguadoExperiment}.  Fig.~{\ref{Fig2}}(d) shows the diamond lineshape for a quasi-ideal case of the isolated MBSs ($\eta = 0.17$), while panels (e) and (f) illustrate the situations where the nonlocal feature was suppressed by enhancing $\lambda_{R}$ and, consequently $\eta$. The loss of the nonlocality ($\eta \rightarrow1$) is related to the displacement of the right MBS ($\gamma_{R}$) wave function towards the left MBS, increasing the overlap between such states and enhancing the hybridization $\lambda_{R}$ of right Majorana mode with the dot state. By comparing the Figs.~{\ref{Fig2}}(d)-(f), it can be noticed that ABS formation due to the strong localization of the right MBS near the QD gives rise to the disintegration of the diamond shapes. In other words, the closer $\epsilon_{\text{MBS}}^{\pm}$ is to $\epsilon_{\text{QD}}^{\pm}$ (See panels (d) and (f)), the higher is the local nature of MBSs and the lower is the topological quality of the device.

\begin{figure}[t]
	\centerline{\includegraphics[width=3.6in,keepaspectratio]{ColormapFullFev.pdf}}
	\caption{Color scale plots of the DOS of a QD as a function of QD energy level $\varepsilon_{d}$ and spectral frequency $\omega$ for distinct regimes. In all situations we chose $\lambda_{L}=1.0E_{0}$. Left panels show (a)isolated MBSs, (c)bowtie and (e)diamond profiles, corresponding to the cases illusrtated by Fig.~\ref{Fig2}(a), (b) and (e), respectively, but for the distinct values of Zeeman splitting $V_{Z}$ and charging energy $U$. Right panels reproduce the corresponding left ones in the reduced gray scale. 
		\label{Newcolormap}}
\end{figure}

It is worth noting that the results presented in the Fig.~\ref{Fig2} differ from analogous results of the earlier works~\cite{RAguado,AguadoExperiment} due to the presence of extra crossing points appearing in the middle region of all the panels. This feature is direct outcome of the theoretical treatment afforded by Hubbard-I approximation~\cite{HubbardI} to charging energy term of the system Hamiltonian [Eq.~(\ref{eq:FullH})]. It is known that such approximation makes the condition of the transition to ferromagnetic state more restrictive compared to Hartree-Fock mean-field approximation, since it accounts for the higher order correlations thus reducing the energy of non-magnetic states with respect to ferromagnetic ones~\cite{HubbardI}. For this reason, Hartree-Fock mean-field approximation works well for larger Zeeman fields, as previously noticed by Prada \textit{et al}~\cite{RAguado}. However, within such a mean-field approach for any value of $V_{Z}$, the information about correlated motion of electrons is only taken into account with the mean occupation~\cite{Bruus}. Hubbard-I decoupling scheme accounts for such correlated motion, which gives rise to the appearance of the so-called Hubbard peaks at $\varepsilon_{d\sigma}$ and $\varepsilon_{d\sigma} + U$, describing the regime of the Coulomb blockade. From the experimental perspective, the work of Deng \textit{et al}~\cite{AguadoExperiment} shows the validity of the Hartree-Fock approach, since such an experiment was performed under relatively larger Zeeman fields, which are enough to resolve the QD spin-degrees of freedom. In this work, by using the Hubbard I decoupling scheme, we predict the system behavior in the scenario of weaker Zeeman splitting, which could be addressed in future experiments. 

Detuning of the parameters $V_{Z}$ and $U$ from the value corresponding to the Figs.~\ref{Fig2}(a), (b) and (e) allows us to recover low-energy spectrum theoretically predicted within the Hartree-Fock mean-field approximation~\cite{RAguado}, which was experimentally verified by Deng \textit{et. al}~\cite{AguadoExperiment}, thus ensuring the comprehensiveness of the approximation employed here. Figs.~\ref{Newcolormap}(a)-(b) illustrate the case of isolated MBSs, but with $V_{Z}$ lower than that of Fig.~\ref{Fig2}(a). Such a Zeeman splitting is unable to resolve spin up and down states of the QD. Consequently, instead of a four-peak structure, profiles with only two peaks resembling those presented in the Fig.~4(d) of the Ref.~\onlinecite{AguadoExperiment} appear. Figs.~\ref{Newcolormap}(c)-(d) show bowtie profiles corresponding to $V_{Z}=2.5E_{0}$ and $U=20E_{0}$ and comparable to those presented in the Fig.~3(b) of Ref.~\onlinecite{AguadoExperiment}. In the corresponding scenario only two peaks are present as well, since higher charging energy sets other peaks outside the considered range of $\varepsilon_{d}$. Increase of the value of $U$ results in the profiles shown in the Figs.~\ref{Newcolormap}(e)-(f) corresponding to the diamond situation shown in the Fig.~\ref{Fig2}(e) and comparable to those presented in the Fig.~3(c) of the Ref.~\onlinecite{AguadoExperiment}.

\begin{figure}[t]
	\centerline{\includegraphics[width=3.5in,keepaspectratio]{Result3.pdf}}
	\caption{ Density of states as a function of the QD level $\varepsilon_{d}$ at $\omega=0$ for (a) isolated MBSs, (b) bowtie and (c) diamond configurations. Plots for several  values of charging energy [$U=2.5E_{0}$ (blue squares), $5E_{0}$ (magenta triangles) and $7.5E_{0}$ (green crosses)] are presented.
		\label{Fig3}}
\end{figure}

\begin{figure}[t]
	\centerline{\includegraphics[width=3.5in,keepaspectratio]{IsolatedMBSs.pdf}}
	\caption{(a)-(d): (a)-(d): total DOS as a function of $\omega$ for MBSs (corresponding to the Fig.2(a)) for various values of $\varepsilon_{d}$ corresponding to the vertical dashed white lines in Fig.~\ref{Fig2}(a). (e)-(h) spin resolved DOS $\rho_{\sigma}(\omega)$ for same parameters as in the panels (a)-(d)
		\label{Isolated}}
\end{figure}

\begin{figure}[t]
	\centerline{\includegraphics[width=3.5in,keepaspectratio]{Bowtie.pdf}}
	\caption{(a)-(d): total DOS as a function of $\omega$ for bowtie profile (corresponding to the Fig.2(b)) for various values of $\varepsilon_{d}$ corresponding to the vertical dashed white lines in Fig.~\ref{Fig2}(b). (e)-(h) spin resolved DOS $\rho_{\sigma}(\omega)$ for same parameters as in the panels (a)-(d)
		\label{Bowtie}}
\end{figure}

Fig~\ref{Fig3} shows the DOS of the dot as a function of QD level $\varepsilon_{d}$ at $\omega=0$ for (a) isolated MBSs, (b) bowtie and (c) diamond situations, for distinct strengths of the dot charging energy [$U=2.5E_{0}$ (blue squares), $5E_{0}$ (magenta triangles) and $7.5E_{0}$ (green crosses) ]. In the highest nonlocal case [Fig~\ref{Fig3}(a)], the insensitivity of the zero frequency peak to the tuning of the QD level and variations of the dot charging energy is verified~\cite{RAguado,DJClarcke,AguadoExperiment}. There is a plateau in the total DOS characteristic to ZBCP. This scenario breaks down for the  possible situation of the formation of ABS due to the overlap between MBSs [Fig~\ref{Fig3}(b-c)]. In this case the plateau in the DOS is destroyed and positions of the peaks change with variations of both position of the dot level $\varepsilon_{d}$ and charging energy $U$. Fig~\ref{Fig3}(b) describes a linecut at $\omega=0$ for the bowtie configuration, wherein the four resonances for $U=5E_{0}$ (magenta triangles) correspond to the anticrossings between the dot level and near zero-energy states appearing in the Fig~\ref{Fig2}(b) (indicated by magenta rectangles). These anticrossing points (resonance positions) are strongly dependent on the charging energy, since near zero-energy ABS, which can be a trivial non-protected state, is affected by the QD energy levels. Similar behavior is found for a diamond profile with degree of nonlocality $\eta =0.5$, as it is shown in the Fig.~\ref{Fig3}(c). Moreover, the plateau depicted in the Fig.~\ref{Fig3}(a) also allows to distinguish the Majorana ZBCP from that induced by usual Kondo effect, since Kondo resonance only appears when a QD is single occupied, $\varepsilon_{d}<\varepsilon_{F}$, as verified by one of us in the Ref.~\onlinecite{Vernek} (see Fig.~1(g)-(i) of that paper).

\begin{figure}[t]
	\centerline{\includegraphics[width=3.50in,keepaspectratio]{Diamond.pdf}}
	\caption{(a)-(d): total DOS as a function of $\omega$ for diamond profile for various values of $\varepsilon_{d}$ corresponding to the vertical dashed white lines in Fig.~\ref{Fig2}(e). (e)-(h) spin resolved DOS $\rho_{\sigma}(\omega)$ for same parameters as in the panels (a)-(d)
		\label{Diamond}}
\end{figure}

The impact of each spin component to the total DOS is presented in the Figs.~\ref{Isolated}, \ref{Bowtie} and \ref{Diamond} for isolated MBSs, bowtie [Fig.\ref{Fig2}(b)] and diamond [Fig.\ref{Fig2}(e)] configurations, respectively. Panels (a)-(d) of the  Fig.~\ref{Isolated} show the DOS of the dot for isolated MBSs configuration [Fig.\ref{Fig2}(a)] as a function of $\omega$ for several values of $\varepsilon_{d}$ indicated by the white dashed lines in the Fig.~\ref{Fig2}(a). As can be seen, a peak at $\omega=0$ emerges for all values of the dot energy as observed in the panels (a)-(d). Panels (e)-(h) of the same figure reveal that zero-peak structure in the total DOS is spin degenerated: at $\omega=0$,  $\rho_{\uparrow}=\rho_{\downarrow}$. 

This degeneracy is broken when MBSs overlap ($\delta_{M}\neq0$) and the dot hybridizes with the rightmost Majorana ($\lambda_{R} \neq 0$) as well, as it is demonstrated for the bowtie configuration in the Fig.~\ref{Bowtie}. In this case a zero-peak structure in the  total DOS emerges only when $\varepsilon_{d}$ crosses zero-energy as it happens in the panel (b).  In this case DOS around $\omega=0$ becomes spin sensitive, $\rho_{\uparrow}\neq\rho_{\downarrow}$, as one can verify from the panel (f), thus suggesting the situation of the  spin-dependent transport. 

Diamond configuration also displays spin-dependent behavior, as it is shown in the Fig.~\ref{Diamond}. However, there is a remarkable difference from the bowtie case: the near zero-energy two-peaks having $\sigma=\uparrow$ are no longer pinned and merge with one another at $\omega=0$ [panels (e) and (g)] for certain values of $\varepsilon_{d}$, giving rise to zero-peak structure in the DOS, as displayed in the panels (a) and (c). Such a feature is consistent with the picture of coalescing ABSs (overlapped MBSs), which can mimic the behavior MBSs under certain conditions~\cite{ABS2,ABS3,ABS4,psABS1,psABS2}. It also can be verified for the situations where the QD level does not shift the near-zero peaks describing ABSs towards $\omega=0$ [panels(f) and (h)] and consequently there is no peak structure at zero-energy in DOS [panels (b) and (d)]. 

It is worth noting that the presence of ABSs in the total DOS of a QD [Eq.~(\ref{eq:DOS})] in the results presented in the Figs.~\ref{Fig3}, \ref{Bowtie} and \ref{Diamond} can be confirmed by the computation of the anomalous Green's function $G_{f_{\uparrow}^{\dagger},d_{\sigma}}^{r}(\omega)$. The latter enters into the Green's function of the QD according to the Eq.~(\ref{eq:A4}),  describes the correlation between  QD ($d^{\dagger}_{\sigma}$) and SC nanowire appearing due Andreev reflection (AR) ~\cite{ABSLee,Baraski2013,RevModPhysimpuritiesSC} and can be defined in the time domain as:~\cite{Bruus}

\begin{equation}
G_{f_{\uparrow}^{\dagger},d_{\sigma}}^{r}(t-t')=-\imath \theta (t-t')\langle \lbrace f_{\uparrow}^{\dagger}(t),d_{\sigma}^{\dagger} (t') \rbrace  \rangle, \label{eq:fdagd}  
\end{equation}wherein $\theta(t-t')$ is Heaviside function, $\lbrace...,...\rbrace$ denotes anticommutator and $f_{\uparrow}^{\dagger}$ corresponds to the creation of a non-local fermion in SC nanowire, which is formed by the linear combination of the left ($\gamma_{L}$) and right ($\gamma_{R}$) MBSs, \textit{i.e,} $f_{\uparrow}^{\dagger}=(\gamma_{L}-\imath\gamma_{R})/\sqrt{2}$ and $f_{\uparrow}=(\gamma_{L} + \imath\gamma_{R})/\sqrt{2}$.

AR can take place through different transport channels, once QD is coupled to the both ends of a nanowire ($\lambda_{L}$ and $\lambda_{R}\neq0$). The presence of the coupling asymmetry can give rise to the appearence of Fano antiresonances ~\cite{RevFano,CPS,Recher,MBSoscillations} as it can be seen in the Fig.~\ref{ARFanodips}, where we illustrate AR process by plotting $\rho_{AR}^{\sigma}(\omega) = -\text{Im} \lbrace G_{f_{\uparrow}^{\dagger},d_{\sigma}}^{r}(\omega)\rbrace/\pi$ for two representative \textit{``bowtie''} and \textit{``diamond''} configurations (see Figs.~\ref{Bowtie}(a)-(e) and~\ref{Diamond}(c)-(g), respectively).

\begin{figure}[t]
	\centerline{\includegraphics[width=3.5in,keepaspectratio]{ARFanodips.pdf}}
	\caption{DOS related to AR processes for each spin component, wherein two representative situations were considered. Panel (a) depicts the same case of Figs.~\ref{Bowtie}(a)-(e), while panel (b) is related to Figs.~\ref{Diamond}(c)-(g).} 
		\label{ARFanodips}
\end{figure}



\section{Conclusions}

We presented the theoretical study of MBSs nonlocality in the hybrid device sketched in the Fig.~\ref{Setup}. It was analyzed in detail how the charging energy affects spin resolved DOS of a QD coupled with MBSs. For highly nonlocal MBSs there is a plateau at zero-energy in the QD density of states for any values of the dot level and charging energy. For overlapping MBSs the spectrum of the dot is strongly modified.  It was shown that the zero-peak structure in the DOS reveals pronounced spin dependence if MBSs become hybridized. Our findings suggest that a spin-dependent local probe may be used as a tool to resolve an outstanding problem in experimental Majorana physics: discriminating between the cases when ZBCP is due to the isolated MBSs or ABSs.
\section{Acknowledgments} 

We thank the funding Brazilian agencies CNPq (Grants No.~307573/2015-0 and No.~305668/2018-8) and S\~ao Paulo Research Foundation (FAPESP) Grant No.~2015/23539-8. IAS acknowledges support by Horizon2020 RISE project CoExAN. MdS acknowledges the Austrian Academy of Science \"OAW for the JESH fellowship and Serdar Sariciftci for the hospitality.
