Manipulating Majorana zero modes in double quantum dots

TL;DR

This paper demonstrates how double quantum dot systems can be used to detect and manipulate Majorana zero modes through electrical control, advancing scalable quantum computing technologies.

Contribution

It introduces a method to manipulate MZMs in double quantum dots using spectral functions and parameter tuning, providing new control strategies.

Findings

Spectral functions reveal MZM presence and leakage into DQDs.

MZM signatures respond predictably to gate voltages and couplings.

Different coupling geometries enable switching MZMs between dots.

Abstract

Majorana zero modes (MZMs) emerging at the edges of topological superconducting wires have been proposed as the building blocks of novel, fault-tolerant quantum computation protocols. Coherent detection and manipulation of such states in scalable devices are, therefore, essential in these applications. Recent detection proposals include semiconductor quantum dots (QDs) coupled to the end of these wires, as changes in the QD electronic spectral density due to the MZM coupling could be detected in transport experiments. Here, we propose that multi-QD systems can also be used to manipulate MZMs through precise control over the QDs' parameters. The simplest case where Majorana manipulation is possible is in a double quantum dot (DQD) geometry. By using exact analytical methods and numerical renormalization-group calculations, we show that the QDs' spectral functions can be used to characterize the presence or not of MZMs "leaking" into the DQD. More importantly, we find that these signatures respond to changes in the DQD parameters such as gate-voltages and couplings in a consistent fashion. Additionally, we show that different MZM-DQD coupling geometries ("symmetric" , "in-series" and "T-shaped" junctions) offer distinct ways in which MZMs can be switched from dot to dot. These results highlight the interesting possibilities that DQDs offer for all-electrical MZM control in scalable devices.

Equations

$$H=H_{\mathrm{DQD}}+H_{\mathrm{lead}}+H_{\mathrm{DQD}-\mathrm{lead}}+H_{\mathrm{M}-\mathrm{DQD}}$$

$$H_{\mathrm{DQD}}=$$

$$+\sum _{\sigma }{t}_{dots}(d_{1\sigma }^{†}{d}_{2\sigma }+d_{2\sigma }^{†}{d}_{1\sigma }),$$

$$H_{\mathrm{lead}}=$$

$$H_{\mathrm{DQD}-\mathrm{lead}}=$$

$$H_{\mathrm{M}-\mathrm{DQD}}=\sum _{\sigma ,i=1}^2{t}_{i\sigma }\left(d_{i\sigma }^{†}{\gamma }_1+{\gamma }_1{d}_{i\sigma }\right).$$

$${\gamma }_1=\frac{1}{\sqrt{2}}\left(f_{\downarrow }^{†}+f_{\downarrow }\right),{\gamma }_2=\frac{i}{\sqrt{2}}\left(f_{\downarrow }^{†}-f_{\downarrow }\right).$$

$$H_{\mathrm{M}-\mathrm{DQD}}$$

$${\rho }_{i\sigma }(\omega )\equiv -\frac{1}{\pi }\text{Im}\left[G_{d_{i\sigma },d_{i\sigma }^{†}}(\omega )\right].$$

$$G_{d_{1\downarrow },d_{1\downarrow }^{†}}\left(\omega \right)=\frac{1}{\omega -{ϵ}_{DQD}^{+}-\frac{{\Vert T_{+}\Vert }^2}{\omega -{ϵ}_M-\frac{{\Vert T_{-}\Vert }^2}{\omega -{ϵ}_{DQD}^{-}}}},$$

$${ϵ}_{DQD}^{\pm }=\pm {ϵ}_1+\sum _{\mathbf{k}}\frac{V_1{V}_1^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}+\frac{{\Vert \pm t_{dots}+{\sum }_{\mathbf{k}}\frac{V_1{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}\Vert }^2}{\omega \mp {ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}},$$

$$T_{\pm }=\pm t_1\pm t_2\frac{\left(\pm t_{dots}+{\sum }_{\mathbf{k}}\frac{V_1{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}\right)}{\omega \mp {ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}},$$

$${ϵ}_M=\omega -\frac{{\Vert t_2\Vert }^2}{\omega -{ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}}-\frac{{\Vert t_2\Vert }^2}{\omega +{ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega +{ϵ}_{\mathbf{k}}}}.$$

$$G_{d_{1\uparrow },d_{1\uparrow }^{†}}\left(\omega \right)=\frac{1}{\omega -{ϵ}_{DQD}^{+}}.$$

$$-i{\Gamma }_i=\underset{s\to 0}{\lim }\sum _{\mathbf{k}}\frac{V_i^{\ast }{V}_i}{\omega +is-{ϵ}_{\mathbf{k}}}.$$

$${\rho }_{1\sigma }(\omega )=-\frac{1}{\pi }\text{Im}\left[G_{d_{1\sigma },d_{1\sigma }^{†}}(\omega )\right].$$

$$G_{A,B}(\omega )\equiv -i\int e^{i\omega t}\Theta (t)⟨\{A(t),B(0)\}⟩dt.$$

$$\omega G_{A,B}(\omega )={\delta }_{A^{†},B}+G_{\left[A,H\right],B}(\omega ).$$

$$\omega G_{f_{\downarrow },d_{1\downarrow }^{†}}(\omega )$$

$$=\sum _{i=1}^2\frac{t_i}{\sqrt{2}}\left(G_{d_{i\downarrow },d_{1\downarrow }^{†}}(\omega )-G_{d_{i\downarrow }^{†},d_{1\downarrow }^{†}}(\omega )\right).$$

$$\mathcal{T}{\vec{G}}_{d_1^{†}}={\hat{e}}_1,$$

$$\left[\begin{array}[]{ccccccc}\omega-\epsilon_{1}&-V_{1}^{*}&-t_{dots}&-t_{1}&0&0&0\\ -V_{1}&\omega-\epsilon_{k}&-V_{2}&0&0&0&0\\ -t_{dots}^{*}&-V_{2}^{*}&\omega-\epsilon_{2}&-t_{2}&0&0&0\\ -t_{1}^{*}&0&-t_{2}^{*}&\omega&t_{2}^{*}&0&t_{1}^{*}\\ 0&0&0&t_{2}&\omega+\epsilon_{2}&V_{2}^{*}&t_{dots}^{*}\\ 0&0&0&0&V_{2}&\omega+\epsilon_{k}&V_{1}\\ 0&0&0&t_{1}&t_{dots}&V_{1}^{*}&\omega+\epsilon_{1}\end{array}\right],$$

$$\displaystyle\Big{[}$$

$$\displaystyle G_{d_{2\downarrow}^{\dagger},d_{1\downarrow}^{\dagger}}(\omega),G_{c_{k\downarrow}^{\dagger},d_{1\downarrow}^{\dagger}}(\omega),G_{d_{1\downarrow}^{\dagger},d_{1\downarrow}^{\dagger}}(\omega)\Big{]}^{T}.$$

$$\left[\begin{array}[]{ccc}\omega-\epsilon_{1}&-V_{1}&-t_{dots}\\ -V_{1}^{*}&\omega-\epsilon_{k}&-V_{2}\\ -t_{dots}^{*}&-V_{2}^{*}&\omega-\epsilon_{2}\end{array}\right],$$

$$\left[\begin{array}[]{ccc}\frac{V_{1}^{*}V_{1}}{\omega-\epsilon_{k}}&-V_{1}^{*}&\frac{V_{2}V_{1}^{*}}{\omega-\epsilon_{k}}\\ -V_{1}&\omega-\epsilon_{k}&-V_{2}\\ \frac{V_{2}^{*}V_{1}}{\omega-\epsilon_{k}}&-V_{2}^{*}&\frac{V_{2}^{*}V_{2}}{\omega-\epsilon_{k}}\end{array}\right].$$

$$\left[\begin{array}[]{ccc}\omega-\epsilon_{1}-\frac{V_{1}^{*}V_{1}}{\omega-\epsilon_{k}}&0&-t_{dots}-\frac{V_{2}V_{1}^{*}}{\omega-\epsilon_{k}}\\ 0&0&0\\ -t_{dots}^{*}-\frac{V_{2}^{*}V_{1}}{\omega-\epsilon_{k}}&0&\omega-\epsilon_{2}-\frac{V_{2}V_{1}^{*}}{\omega-\epsilon_{k}}\end{array}\right]$$

$${ϵ}_{DQD}^{+}={ϵ}_1+\sum _{\mathbf{k}}\frac{V_1{V}_1^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}+\frac{{\Vert t_{dots}+{\sum }_{\mathbf{k}}\frac{V_1{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}\Vert }^2}{\omega -{ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}}$$

$${ϵ}_{DQD}^{\pm }=\pm {ϵ}_1+\sum _{\mathbf{k}}\frac{V_1{V}_1^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}+\frac{{\Vert \pm t_{dots}+{\sum }_{\mathbf{k}}\frac{V_1{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}\Vert }^2}{\omega \mp {ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}}$$

$$T_{\pm }=\pm t_1\pm t_2\frac{\left(\pm t_{dots}+{\sum }_{\mathbf{k}}\frac{V_1{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}\right)}{\omega \mp {ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}}.$$

$${ϵ}_M=\omega -\frac{{\Vert t_2\Vert }^2}{\omega -{ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega -{ϵ}_{\mathbf{k}}}}-\frac{{\Vert t_2\Vert }^2}{\omega +{ϵ}_2-{\sum }_{\mathbf{k}}\frac{V_2{V}_2^{\ast }}{\omega +{ϵ}_{\mathbf{k}}}}.$$

$$G_{d_{1\downarrow },d_{1\downarrow }^{†}}\left(\omega \right)=\frac{1}{\omega -{ϵ}_{DQD}^{+}-\frac{{\Vert T_{+}\Vert }^2}{\omega -{ϵ}_M-\frac{{\Vert T_{-}\Vert }^2}{\omega -{ϵ}_{DQD}^{-}}}}.$$

Full text

I Introduction

The search for Majorana quasiparticles in condensed matter systems has gained renewed attention in the last decade, motivated by the exciting prospects for achieving scalable, fault-tolerant topological quantum computation protocols in nanoscale devices Nayak et al. (2008); Alicea (2012); Beenakker (2013); Aguado (2017). From the landmark theoretical proposals for realizing Kitaev’s model of a 1D topological superconductor Kitaev (2001, 2003) using semiconductor quantum wires with strong-spin orbit coupling and proximity-induced superconductivity Lutchyn et al. (2010); Oreg et al. (2010), the field rapidly evolved towards the first experiments reporting data consistent with these predictions (Mourik et al., 2012; Das et al., 2012; Deng et al., 2012). As a result, the last few years have been full of excitement as continuing improvements in sample growth and characterization techniques have allowed for more consistent experimental evidence for Majorana bound states in semiconductor quantum wires Deng et al. (2016); Zhang et al. (2018); Deng et al. (2018).

In these set-ups, “Majorana signatures” are characterized by zero-bias signals in the conductance across the device due to the emergence of robust zero-energy modes localized at the edges of the quantum wire. An important requirement is to distinguish these so-called Majorana zero modes (MZM) from other zero-bias phenomena, such as the Kondo effect Hewson (1997), which have been found in similar systems Lee et al. (2012). A large effort in recent experimental proposals was put on ways to uniquely identify MZMs, including proposals for measuring the signatures of non-Abelian statistics Aasen et al. (2016); Sarma et al. (2015); Heck et al. (2012). Although this last property is crucial in the implementation of fault-tolerant quantum computers, its measurement has been elusive so far as it requires overcoming several experimental issues related to the need of “moving” Majorana quasiparticles in order to perform braiding operations Hoffman et al. (2016); Karzig et al. (2017).

A rather straightforward proposal to detect MZMs consists of attaching a quantum dot (QD) to the edge of a topological quantum wire and then measuring the electrical conductance through the QD Liu and Baranger (2011). In such an arrangement, the MZM at the end of the chain “leaks” into the attached QD Vernek et al. (2014), reducing the zero-bias conductance through the dot by the sizable amount of $\frac{e^{2}}{2h}$. This detection method offers two key advantages over other approaches: (i) no direct charge transfer between the MBS and the dot is necessary, thus preventing “quasiparticle poisoning” Rainis and Loss (2012), and (ii) a clear distinction with Kondo physics is warranted, even if the experiment is performed at temperatures below the Kondo temperature $T_{K}$ Lee et al. (2013); Ruiz-Tijerina et al. (2015); Gorski et al. (2018). Recently, topological quantum wire-QD junctions have been realized in experiments Deng et al. (2016, 2018), paving the way for further experiments involving the detection of MZMs using quantum dots.

Venturing beyond simple “detection” setups, the large degree of control over the QD parameters offers the unique possibility of manipulating MZMs inside multidot systems. The simplest case where Majorana manipulation is possible is in a double quantum dot (DQD). Tunneling Majorana modes in these basic structures have inspired theoretical studies Silva and Vernek (2016); Ivanov (2017); Rančić et al. (2019); Ricco et al. (2019) and experimental setups confirming the observations of Andreev molecules Su et al. (2017). However, despite the fact that DQDs offer several possibilities for manipulation of MZMs, there is still no complete analysis of the possible transitions of these Majorana signatures between the QDs even in a simple model.

In this paper, we explore the different possibilities for Majorana manipulation in a device consisting of a DQD coupled to a MZM and a metallic lead (see Fig. 1). The simplicity of this model allows us to analytically explore different geometries of QD’s from symmetric and “in-series” couplings to T-shaped junctions (Fig. 2). We considered both noninteracting and interacting regimes, observing major agreement between both approaches about the location of the Majorana signature.

We performed a detailed study of the noninteracting DQD limit, by using Zubarev’s procedure Zubarev (1960) to provide an exact formula to calculate the spectral functions. For the interacting case, we resort to numerical renormalization group (NRG)Bulla et al. (2008) calculations for this model. While the noninteracting regime is suitable to obtain exact expressions for the Green function, the interacting case shows how the Majorana signature co-exists with strongly correlated phenomena such as the Kondo effect Hewson (1997) and RKKY interactions Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957).

This paper is organized as follows. In Sec. II we describe the model of a DQD coupled to a MZM and to a metallic lead, as well as the methods used. The results are presented in Sec. III where we compare the noninteracting density of states (LDOS) (Sec. III.1) with the low-energy interacting results in Sec. III.3. Finally, our conclusions are given in Sec. IV.

II Model and methods

We consider the setup shown in Fig. 1, in which a single MZM $\gamma_{1}$ located at the edge of a 1D topological superconductor is coupled to a double quantum dot (DQD) attached to a single metallic lead. The Hamiltonian of the entire system can be expressed as:

[TABLE]

where the different terms describe, respectively, the (interacting) DQD, the (noninteracting) metallic lead, and the DQD-lead couplings, written as:

[TABLE]

while the DQD-MZM coupling is given by Liu and Baranger (2011); Lee et al. (2013); Ruiz-Tijerina et al. (2015); Hoffman et al. (2017); Prada et al. (2017):

[TABLE]

In the equations above, $\epsilon_{i}$ is the energy level of dot $i$, $U_{i}$ is the Coulomb repulsion and $t_{dots}$ is the coupling parameter between both QDs. The operator $d^{\dagger}_{i\sigma}$ creates a particle in dot $i$ with spin $\sigma$ and $\hat{n}_{i\sigma}=d_{i\sigma}^{\dagger}d_{i\sigma}$ is the particle number operator of state $i$, $c_{\mathbf{k}\sigma}^{\dagger}$ is the creation operator a particle with momentum $\mathbf{k}$ and spin $\sigma$ in the lead. Finally, $\epsilon_{\mathbf{k}l}$ is the corresponding energy and $V_{i}(\textbf{k})$ describes the tunneling coupling between the lead and dot $i$.

We take the length of the wire to be large so that we can safely neglect both the overalap between the two Majorana modes and the (much smaller) coupling between the DQD and the $\gamma_{2}$ MZM located at the other edge of the wire. We also note that the DQD-$\gamma_{1}$ coupling strength $t_{\sigma i}$ in Eq. (3) above is, in general, spin-dependent Hoffman et al. (2017); Prada et al. (2017); Deng et al. (2018) and can be written in terms of the $\gamma_{1}$ MZM “spin canting angle” $\theta_{1}$ as $\left(t_{\uparrow i},t_{\downarrow i}\right)\equiv t_{i}\left(\sin{\frac{\theta_{1}}{2}},-\cos{\frac{\theta_{1}}{2}}\right)$. For the purposes of this work, we take $\theta_{1}\!=\!\pi$ such that only spin down dot operators are coupled to the MZM, making $H_{\rm M-DQD}$ fully spin-conserving. This choice adds an extra symmetry (spin down parity) to the full Hamiltonian, which will turn out to be important in the NRG calculations presented in Sec. III.3.

It is also useful to recast the last term of Eq. (1) in terms of (Dirac) fermionic operators. Following Refs. Lee et al., 2013; Ruiz-Tijerina et al., 2015, we choose to write the MZMs $\gamma_{1}$ and $\gamma_{2}$ as a superposition of the creation ($f^{\dagger}_{\downarrow}$) and annihilation ($f_{\downarrow}$) operators of a (fully polarized) spin down fermion:

[TABLE]

In this representation, the effective coupling between the MZM and the DQD given by Eq. (3) becomes:

[TABLE]

where $t_{i}\equiv t_{\downarrow i}$ is the coupling strength between the MZM and QD $i$.

In order to identify the presence/absence of MZMs “leaking” from the edge of the TS into the dots Liu and Baranger (2011); Vernek et al. (2014); Lee et al. (2013); Ruiz-Tijerina et al. (2015), the quantities of interest are the spin-resolved spectral functions (or, equivalently, the local density of states) of the quantum dots. As usual, the spectral function for spin $\sigma$ in dot $i$ is defined as:

[TABLE]

where $G_{d_{i\sigma},d_{i\sigma}^{\dagger}}(\omega)\equiv\langle\langle d_{i\sigma},d_{i\sigma}^{\dagger}\rangle\rangle_{\omega}$ is the retarded (diagonal) Green’s function involving dot $i$ operators $d_{i\sigma}$ and $d_{i\sigma}^{\dagger}$. Next, we describe the procedures for calculating $\rho_{i\sigma}(\omega)$ in the regimes of weak ($U_{i}\ll V,t_{i}$) and strong ($U_{i}\!\gg\!V,t_{i}$) electron-electron interaction in the dots.

II.1 Noninteracting limit: Graph-Gauss-Jordan elimination

In the noninteracting limit ($U_{i}\!=\!0$), the Hamiltonian $H$ is quadratic in the fermionic operators and we can obtain analytic expressions for the spectral densities defined in Eq. (6). Using Zubarev’s equation of motion (EOM) approach Zubarev (1960), we can derive exact expressions for the Green functions associated to both quantum dot operators $(G_{d_{1}d^{\dagger}_{1}}(\omega),G_{d_{2}d^{\dagger}_{2}}(\omega))$.

The EOM equations define a $8\times 8$ linear system where the Hamiltonian parameters $(t_{1},t_{2},\epsilon_{1}\ldots)$ and the energy $\omega$ are taken as algebraic variables. The solution for these types of equations is a finite continued fraction of multivariate polynomials with maximum degree $8$, which makes it difficult to provide an exact solution using either analytic or numerical methods. To bypass this problem, we introduce a Graph-Gauss-Jordan elimination process Spielman (2010) to iteratively solve the coupled equations of motion. We briefly describe the procedure here.

We begin by representing the Majorana-DQD quantum dot system in a “flow graph”, where each spin-resolved fermion operator (e.g., $d^{\dagger}_{1\downarrow}$, $d_{1\downarrow}$, $f_{\downarrow}$, $f^{\dagger}_{\downarrow}$, etc.) is represented as a “vertex”. The coefficients of the quadratic terms (such as $d^{\dagger}_{1\downarrow}d_{1\downarrow}$ or $c^{\dagger}_{k\downarrow}c_{k\downarrow}$, etc.) are associated to each node as “self-energies” while the coupling terms involving two fermion operators (such as $d^{\dagger}_{1\downarrow}f_{\downarrow}$ or $d^{\dagger}_{1\downarrow}f^{\dagger}_{\downarrow}$, etc.) are associated to the “links” connecting the respective vertices (see Fig. 10 in the Appendix).

We then proceed to iteratively remove both vertices and links by rewriting the self-energies and couplings in terms of the eliminated variables, such that each vertex elimination depicts another step in the Gauss-Jordan process. In the end, the self-energy of the only remaining vertex will contain the full information needed in order to compute the target Green’s function.

This method proved to be efficient in solving complex systems of coupled Green’s functions as the graph elimination process provides a natural linear algorithm to compute the targeted continued fraction. Moreover, the graphic representation simplifies the procedure and allows one to readily identify minimal coupling points, which could reduce the complexity of the solution. A detailed description of the method is given in Appendix A.

After applying the Graph-Gauss-Jordan process, we obtain a closed form for the noninteracting Green’s functions. For instance, the GF for dot 1 (which is directly coupled to the MZM) will be given by:

[TABLE]

where the energies $\epsilon_{DQD}^{\pm}$ are given by

[TABLE]

[TABLE]

and

[TABLE]

The spin-up spectral density, which is not coupled to the MZM, can be obtained by taking $t_{1},t_{2}=0$ in Eqs. (7)-(10), hence giving

[TABLE]

The final results will depend on the broadening parameter of QD $i$ with the lead $(\Gamma_{i})$, given, in the broad-band limit, by:

[TABLE]

Finally, we compute the spin-resolved LDOS in dot 1 as:

[TABLE]

Similar results can be obtain for the LDOS of the second $\rho_{2\sigma}$ by exchanging the indexes $1$ and $2$ in Eq. (11).

II.2 Interacting limit: Wilson’s NRG

To address the case of interacting quantum dots, we employ the numerical renormalization group (NRG), one of the most successful methods to study interacting quantum impurity models (QIMs) Wilson (1975); Sindel (2005); Bulla et al. (2008). In general, a QIM describes a system spanning a finite and relatively small Hilbert space (the “impurity”) coupled to a much larger system (a “continuum”), spanning a large (typically infinite) Hilbert space. As it turns out, the Hamiltonian in Eq. (1) can be cast as a QIM where the impurity is the DQD coupled to the Majorana mode, which is then coupled to the continuum of electrons in the metallic leads.

We notice that the DQD-Majorana tunneling term given by Eq. (5) effectively breaks total spin $S_{z}$ and charge $Q$ conservation of the whole system, while it preserves spin-$\downarrow$ parity $P_{\downarrow}\!=\!\pm 1$ and spin up particle number $N_{\uparrow}$. To improve the efficiency of the method, we used these symmetries to maintain a block structure during NRG’s iterative diagonalization process Bulla et al. (2008); Lee et al. (2013); Ruiz-Tijerina et al. (2015). Both the states serving as a basis for the initial impurity Hamiltonian and the single-site Wilson chain states can be grouped in $(N_{\uparrow},P_{\downarrow})$ blocks. Thus, the $(N_{\uparrow},P_{\downarrow})$ block structure is preserved during the entire NRG iteration process Bulla et al. (2008). In order to compute the (interacting) spectral functions, we use the density matrix NRG (DM-NRG) procedure Hofstetter (2000) in combination with the z-trick method Oliveira and Oliveira (1994), which improves spectral resolution at high energies. We have checked the accuracy of the results by comparing the results with the complete Fock space method Peters et al. (2006) for some of the parameters used.

III Results

For the remainder of the paper, we will focus on the Majorana-DQD coupling geometries depicted in Fig. 2: a “symmetric coupling” arrangement [Fig. 2(a)], a “T-shaped” configuration [Fig. 2(b)] and the case where the Majorana and both dots are coupled “in-series” [Fig. 2(c)]. As we shall see, the intensity of the MZM signature in each dot can be controlled by external gate-voltages which change the position of the dot levels $\epsilon_{1,2}$ relative to the Fermi energy in the leads.

As mentioned previously, the spin-resolved spectral density (or local density of states LDOS) of each quantum dot provides significant information about the effective tunneling (or not) of a Majorana zero mode into the dot. By comparing the spectral densities for the cases with and without DQD-Majorana couplings, we could identify two generic types of signatures of the Majorana presence in quantum dot $i$, which we define as follows:

• •

**Type I: ** The spin-down LDOS is half of the spin-up LDOS at the Fermi energy $(\rho_{i\downarrow}(0)=\rho_{i\uparrow}(0)/2)$.

• •

**Type II: ** The spin-up spectral density of dot $i$ shows a zero mode of height $\rho_{i\downarrow}(0)=\frac{0.5}{\pi\Gamma_{1}}$ while no such signature appears in the spin-up spectral density.

We should point out that the identification of MZM signatures is much simpler in the case of single-quantum dots coupled to an MZM mode considered, e.g., in Refs. Liu and Baranger, 2011; Vernek et al., 2014; Ruiz-Tijerina et al., 2015. In that case, the generic MZM signature is essentially given by the ”type-II” condition defined above, which amounts to the spin-down spectral function value at the Fermi energy pinned at $0.5/(\pi\Gamma)$. If the QD energy level is tuned to the particle-hole symmetric case, the spin up spectral function pinned at $1/(\pi\Gamma)$ such that both “Type-I” and “Type-II” conditions apply.

In the double quantum dot set-up considered here, things are more complex as quantum interference effects give rise to situations where only one of these conditions is met. As we shall see in the following sections, either one of these two types of signatures appear over a wide range of parameters in our results. In practice, we find that the appearance of Type-I or Type-II signatures are related to the behavior of the spin-up spectral density near the Fermi energy $\rho_{i\uparrow}(\omega\!\sim\!0)$: Type I often appears when $\rho_{i\uparrow}(0)$ displays a peak, while Type II typically emerges in situations where $\rho_{i\uparrow}(0)\approx 0$.

Hereafter, we shall refer to “MZM manipulation” the changes in the Majorana signatures in the dot spectral functions induced by the tuning of the dot gate voltages $(\epsilon_{1},\epsilon_{2})$ in the three different setups depicted in Fig.2. In each case, we consider definite values of the couplings $\Gamma_{2}$, $t_{dots}$, $t_{1}$, and $t_{2}$, as follows. In the configuration shown in Fig.2(a), we coupled the QD symmetrically to the lead and the MZM by setting $t_{1}\!=\!t_{2}$. Within this setup, we expect the MZM signature to “split” due to quantum interference and identical signatures should appear in the spectral densities of both dots. We also considered setups in which only one of the dots is coupled directly the MZM or to the metallic lead. Hence, there are only two distinct coupling geometries: either both the MZM and the lead are coupled to the same dot, forming a “T-junction” or “side-dot” configuration ($t_{2(1)}\!=\!0$ and $\Gamma_{2(1)}\!=\!0$), as shown in Fig.2(b). Alternatively, the MZM can be coupled to one of the dots and the lead to the other, such that the MZM and dots are coupled in series [$t_{1(2)}\!=\!0$ and $\Gamma_{2(1)}\!=\!0$, see Fig.2(c)]. In the remainder of the paper, we take $\Gamma_{1}$ as the energy unit.

III.1 MZM manipulation in noninteracting quantum dots

The noninteracting results for setups (a),(b), and (c) of Fig. 2 are shown in Figs. 3, 4 and 5 respectively. In all cases, the left (right) panels depict the spectral density of dot $1$ (dot $2$). Each row represents a different gate voltage configuration in the dots, starting with $\epsilon_{1}\!=\!\epsilon_{2}\!=\!0$ (first row), $\epsilon_{1}\!=\!5\Gamma_{1}$, $\epsilon_{2}\!=\!0$ (second row), and finally $\epsilon_{1}\!=\!0$, $\epsilon_{2}\!=\!\!-5\Gamma_{1}$ (third row). The insets in each row show where the Majorana signature, represented by a red dashed circle inside the dot, is mainly located.

Figure 3 shows results for the symmetric coupling setup [Fig. 2(a)] in the noninteracting regime. For the particle-hole symmetric case (first row), both spin-down [$\rho_{\downarrow}(\omega)$, thin red line] and spin-up [$\rho_{i\uparrow}(\omega)$, bold blue line] spectral densities are identical in both dots, as expected. Notice, however, that the $\rho_{\downarrow}(\omega)$ shows a three-peak structure, a consequence of the coupling with the Majorana mode. Moreover, the spin-down LDOS value at the Fermi energy is half of the respective spin-up LDOS value $(\rho_{i\downarrow}(0)=\frac{1}{2}\rho_{i\uparrow}(0))$, which signals the MZM tunneling into the dots. This Majorana signature is similar to the one observed when a single dot is coupled to a Majorana mode Liu and Baranger (2011); Vernek et al. (2014) and falls in our “type-II” category mentioned above. We thus may conclude that the MZM is delocalizing into both dots in this symmetric configuration.

More interesting, we find that such delocalization can be reversed (and thus manipulated) by applying gate voltages in the dots. If a positive or negative gate voltage is induced in one of the dots, the spin-down LDOS at the Fermi energy can vanish at that dot while the type-I MZM signature $\rho_{i\downarrow}(0)=\frac{1}{2}\rho_{i\uparrow}(0)$ remains in the other dot. This is shown in panels (c)-(f) of Fig. 3 for the case of positive [Figs. 3(c) and 3(d)] and negative [Figs. 3(e) and 3(f)] gate voltages.

The location of the MZM signature can also be controlled by quantum interference, as illustrated in panels (a) and (b) of Fig. 4. Here, the MZM is coupled directly only to dot 1, which is then coupled to the lead, while dot 2 is coupled only to dot 1 via the interdot tunneling term, resulting in a “side-dot” configuration [see Fig. 2(b)]. Interestingly, if the energy level of dot 2 is fixed to be in resonance with the Fermi energy of the lead, quantum interference causes the spectral function in dot 1 to vanish at the Fermi level [Fig. 4(a)], while a type-I MZM signature appears in dot 2 only [Fig. 4(b)]. This interference-induced MZM signature in dot 2 is robust against shifts in dot 1’s gate voltage, as depicted in Figs. 4(c) and (d). While dot 1’s LDOS is pinned at zero at the Fermi energy, dot 2’s spin-down LDOS exhibits a robust zero-mode of height $\frac{0.5}{\pi\Gamma}$, which is a type-II MZM signature.

This qualitative picture is radically altered when dot 2’s gate voltage is shifted away from zero [Figs. 4(e) and 4(f)]. In this case, dot 2 is no longer in resonance with the leads, which changes the interference conditions such that dot 1 spectral function is no longer pinned at zero. The plots clearly show that the MZM signature, previously located in dot 2, now appears in dot 1. Moreover, the spin-up and spin-down LDOS in dot 1 become very similar to the spectral densities observed in the case of a single dot Liu and Baranger (2011); Vernek et al. (2014), which indicates that dot 2 is essentially decoupled from the MZM.

Finally, we consider the “in-series” configuration depicted in Fig. 2(c), in which is similar to the “side-dot” configuration [Fig. 2(b)] except for the fact that the (spin-down) MZM is coupled only to dot 2. Thus, results for the spin-up LDOS are identical to those shown in Fig. 4. However, the MZM signatures in the spin-down LDOS are quite distinct. As an example, when both dots are in resonance with the lead [Figs. 5(a) and 5(b)], the spin-down LDOS does not vanish at $\omega\!=\!0$ as in the previous case. Instead, both dots show $(\rho_{\downarrow}(0)=\frac{0.5}{\pi\Gamma})$, which leads to MZM signatures of type-I in dot 2 and type-II in dot 1.

An important feature of the “in-series” geometry is that dot 1 presents a robust, gate-voltage-independent MZM signature, despite the fact that it is not directly attached to the topological wire. As shown in Figs. 5(c) and 5(d), while a a shift in dot 1’s gate voltage erases the MZM signature in dot $2$, it does not affect the MZM signature in dot $1$. At the same time, the MZM signature in dot 1 is robust against changes in the dot 2’s gate voltage, as shown in Figs. 5(e) and 5(f). The only difference here is that MZM signature types are reversed: dot 1 now shows a type-I signature while dot 2 shows a type-II one.

III.2 Interacting dots: MZM-mediated indirect exchange

We now turn to the more realistic case of quantum dots in the Coulomb blockade regime where local electron-electron interaction terms are relevant. We consider the dots to be in an odd-$N$ Coulomb blockade valley where Kondo correlations are dominant at low-temperatures. The local Coulomb energy in the dots is accounted for by the terms $\frac{U_{i}}{2}(\sum_{\sigma}\hat{n}_{i\sigma}-1)^{2}$ in Eq. (2). For simplicity, we consider equal Coulomb repulsion energies $(U_{1}\!=\!U_{2}\equiv U)$ for both dots. For concreteness, the NRG calculations were performed with $U\!=\!17.3\Gamma_{1}$ in both dots and a half-bandwidth of the lead electrons set at $D=2U=34.6\Gamma_{1}$.

Let us review some of the main features of the spectral densities of the dots in the absence of the MZM coupling. For a single dot coupled to a metallic lead, the Kondo effect is characterized by the appearance of a sharp resonance in the spectral function near the Fermi energy with a width of order $k_{B}T_{K}\!\sim\!\sqrt{U\Gamma_{1}}\exp\left[-\pi\frac{|\epsilon_{1}||\epsilon_{1}+U|}{U\Gamma_{1}}\right]$. Here, $T_{K}\ll U$ is the Kondo temperature of the system Hewson (1997), which will be largest at the particle-hole symmetric point (phs) $\epsilon_{1}\!=\!-\frac{U}{2}$. In the case of two dots at phs $\left(\epsilon_{di}\!=\!-\frac{U_{i}}{2}\right)$, both symmetrically coupled to a single lead ($\Gamma_{1}\!=\!\Gamma_{2}$), there will be an additional effective exchange interaction between the dots mediated by the lead Žitko and Bonča (2007); Liao et al. (2015). Such exchange will compete with the antiferromagnetic Kondo coupling, producing a three-peak structure in the spectral density of both dots.

Figure 6(a) shows the spectral functions for both dots in this case. At large energies, the spectral density displays Hubbard peaks at $\omega\sim\epsilon_{di}\pm 8.6\Gamma_{1}=\pm\frac{U}{2}$, representing the single-particle hole and electron excitations whose width is of order $\sim 4\Gamma_{1}$. At low energies, the spin-independent spectral densities show a central Kondo peak accompanied by indirect-exchange-induced satellite peaks at $\omega\sim\pm 3.46\Gamma^{2}_{1}/U$, giving an energy separation that scales as $\sim\Gamma^{2}_{1}/U$ [see also insets in Fig. 6(a)].

Such exchange-driven three-peak structure remains when the MZM is coupled to the system in the symmetric coupling configuration, as shown in Figure 6(b). More striking is that the indirect-exchange splitting between the dots increases considerably with the MZM coupling up to $\sim\pm 8\Gamma_{1}^{2}/U$ : Our calculations show that the peak separation of the Majorana satellites increases quadratically with the MZM coupling $t_{1}\!=\!t_{2}$ as $4t_{1}^{2}/U$ and this effect enters in superposition with the indirect-exchange-induced satellites in Fig. 6(a). This indicates a MZM-mediated spin-spin correlation between the quantum dots. Thus, the coupling to a spin-down-polarized MZM (which is the case) affects the spin-up component of the spectral densities through this indirect spin-spin interaction. Additional details of these interesting features will be discussed elsewhere Cifuentes and Dias da Silva .

III.3 MZM manipulation in interacting dots

One of the key results of Fig. 6 is that, in the symmetric configuration, both quantum dots display type-I MZM signatures ($\rho_{i\downarrow}(0)=\frac{1}{2}\rho_{i\uparrow}(0)$) co-existing with a Kondo-related zero-energy peak. A similar result was reported in Refs. Lee et al. (2013); Ruiz-Tijerina et al. (2015) for the simpler case of an topological nanowire coupled to a single quantum dot. As in that case, here both Kondo and MZM signatures occur in the low-energy part of the spectral function $\omega\lesssim\Gamma_{1}$, as illustrated in the inset of Fig. 6(b). Within this scale, we can trace some interesting parallels with the noninteracting regime.

As an example, Fig. 7 shows the NRG results for the symmetric setup in Fig. 2(a). As in the noninteracting case [Fig. 3], type-I MZM signatures appear in both dots. These signatures can be manipulated by tuning the gate voltage of one of the dots to induce the MZM signature to appear only in the other dot. The LDOS at figures Fig. 7(d) shows a type-I MZM signature with $\rho_{\downarrow}(0)\approx\frac{1}{2}\rho_{\uparrow}(0)$. This MZM signature is stable against gate-voltage-induced energy shifts in dot 2 away from particle-hole symmetry ($\Delta\epsilon_{2}\equiv\epsilon_{2}+U/2$) in the range $\Delta\epsilon_{2}\lesssim 6\Gamma_{1}$ (see Fig. 7(e)). For larger values of $\Delta\epsilon_{2}$, dot 2 enters the mixed-valence regime and the Coulomb peak originally located at $\omega\sim\pm 8.7\Gamma_{1}$ for $\Delta\epsilon_{2}\!=\!0$ now overlaps with the Fermi energy and both Majorana and Kondo signals are lost.

Results for the interacting “side-dot” set-up [Fig. 2(b)] are shown in Fig. 8. As in the noninteracting case, the spin-up spectral density of dot 1 vanishes at the Fermi level due to single-particle quantum interference, as shown in Fig. 8(a). In dot 2, the spectral density is drastically reduced at the Fermi level, but it remains nonzero [Fig. 8(b) and inset], while still showing a type-I MZM signature, namely, $\rho_{2\downarrow}(0)=\frac{\rho_{2\uparrow}(0)}{2}$. This picture is qualitatively similar to the noninteracting case discussed previously, but it begs the question of what is the fate of the Kondo resonance in the dots in this configuration.

To try and answer this question, we note that a similar interplay between Kondo physics and single-particle interference on a T-shaped double dot geometry has been studied in earlier works by one of us Dias da Silva et al. (2006, 2008, 2017). It has been established that, for the case of the dot coupled to the lead (dot 1, in the present case) being noninteracting, its spectral density vanishes at the Fermi energy, while the spectral density in the second dot (dot 2) shows a “splitted” Kondo resonance for strong enough interdot coupling. The Kondo screening in this second dot, however, is still present. In fact, the Kondo temperature increases with the interdot coupling Dias da Silva et al. (2006, 2017). Here the situation is slightly different as dot 1 is also interacting but we believe the analogy still holds. This picture would explain why the up and down components of the spectral density in dot 2 do not vanish at the Fermi energy (although they are quite suppressed) while still showing the MZM type-I signature ($\rho_{2\downarrow}(0)=\frac{1}{2}\rho_{2\uparrow}(0)$).

When gate voltages are applied in either dot 1 or dot 2, a MZM signature appears in dot 1. This is shown in Figs. 8(c)-(f): a type-II MZM signature ($\rho_{1\downarrow}(0)=\frac{0.5}{\pi\Gamma_{1}}$, $\rho_{1\uparrow}(0)\approx 0$) appears in dot 1 while neither type-I or type-II signatures are evident in dot 2. This is clearly distinct from the noninteracting case, in which a shift in the gate voltage of dot 1 [Figs. 4(c) and 4(d)] leads to a type-II MZM signature in dot 2 and vice-versa. When interactions are present and the system is tuned out of the particle-hole symmetric point, no clear type-I or type-II MZM signatures appear in dot 2’s spectral density ([see Figs. 8(d) and 8(e)]. Instead, $\rho_{2\downarrow}(\omega)$ displays asymmetric resonances near the Fermi energy with comparable widths as the resonances in $\rho_{1\downarrow}(\omega)$. We attribute those to Fano-like single-particle interferences with dot 1 which are common in T-shaped structures Dias da Silva et al. (2008).

Finally, Fig. 9 depicts the NRG results for the “series” configuration shown in Fig. 2(c). In this configuration, the MZM is coupled directly to dot 2 only. As in the noninteracting case, the most consistent MZM signatures (type-II, in this case) occur in dot 1’s spectral properties. As an illustration, Figs. 9(a), 9(c) and 9(e) show robust zero-energy peaks in the spin down spectral densities of dot 1 obeying $\rho_{1\downarrow}(0)=\frac{0.5}{\pi\Gamma_{1}}$ while $\rho_{1\uparrow}(0)\approx 0$. The strong difference between spin up and down spectral densities clearly identifies this as a MZM signature rather than a Kondo peak.

The type-II MZM signature remains in dot 2 despite changes in either $\Delta\epsilon_{1}$ [Fig. 9(c)] or $\Delta\epsilon_{2}$ [Fig. 9(e)]. Moreover, a type-I MZM signature also appears in dot 2 in the particle-hole symmetric case, as depicted in Fig. 9(b). Away from particle-hole symmetry, the MZM traces in the spectral properties of dot 2 are less clear. While a shift in the dot 1 energy $\Delta\epsilon_{1}=+5\Gamma_{1}$ has little effect in the dot 2 spectral density [Fig. 9(d)], changing the energy of dot 2 by an amount $\Delta\epsilon_{2}=-5\Gamma_{1}$ gives a zero-energy peak in $\rho_{2\downarrow}(\omega)$ [Fig. 9(f)] which essentially meets the type-II MZM signature condition for these parameters.

Although this result seems to agree with the noninteracting case (see Fig. 5(f)), a closer inspection shows that this MZM signature in $\rho_{2\downarrow}(\omega\!\sim\!0)$ is parameter-dependent. In fact, the height of the zero-energy peak scales roughly as $\rho_{2\downarrow}(\omega\!\sim\!0)\sim(\Delta\epsilon_{2})^{2}$ for $\Delta\epsilon_{2}\lesssim 6\Gamma_{1}$. Thus, for interacting dots, the categorization of a MZM signature in dot 2 is clear only in the $4\Gamma_{1}\lesssim\Delta\epsilon_{2}\lesssim 6\Gamma_{1}$ range, in contrast with the noninteracting case where the MZM signature is robust and largely $\epsilon_{2}$-independent.

One way to understand these features is to use the “Majorana leaking” analogy of Ref. Vernek et al. (2014). In the series configuration of Fig. 2(c), both dots can be though as nontopological “extensions” of the Kitaev chain, with dot 1 being the “last site” or the “edge”. Thus, due to the leaking of the MZM to the neighboring sites (as it is the case of a MZM attached to a single quantum dot Vernek et al. (2014); Ruiz-Tijerina et al. (2015)), it would be expected that edge-mode signatures in dot 1 would be quite robust against changes in gate voltages.

IV Concluding remarks

In this paper, we have addressed the following question: Can one manipulate and detect Majorana zero modes (MZMs) in an all-electric set-up using semiconductor double quantum dots? To this end, we considered a minimal model of a MZM coupled to a double quantum dot (DQD) and metallic leads and calculated the spectral signatures in both strongly- and weakly-interacting regimes. By comparing exact analytical solutions in the noninteracting system and numerical renormalization-group results for interacting quantum dots, we were able to characterize the displacements of the MZM inside the double quantum dot for the three setups in Fig. 2.

Our results for both weakly- and strongly-interacting regime show that gate-voltage tuning in the dots allows for an effective manipulation of the tunneling of the MZM into the DQD system. By considering different MZM-DQD coupling geometries (“symmetric”, “T-shaped” and “in-series”) we found that the presence or not of the MZM in each dot can be monitored by two types of signatures in the spectral density (or local density of states) of the dots.

In the symmetric configuration, the MZM is equally coupled to both dots. As in a “double slit” set up, the MZM signature will appear in both dots if the gate voltages are tuned to the particle-hole symmetric (phs) point. By changing the gate voltage in one of the dots (the equivalent of “closing one of the slits”), the MZM signature will move to the other dot. In the “T-shaped” configuration, when the MZM is directly coupled only to dot 1, the MZM signature will appear only in one of the dots: either dot 2 (at phs or if a gate voltage is applied to dot 1) or in dot 1, when a gate voltage is applied to dot 2.

We also considered a configuration with the MZM coupled “in-series” with both dots, which is closely connected with recent design proposals for topological quantum computational circuits involving MZMs Hoffman et al. (2016); Karzig et al. (2017). In this case, there is a robust MZM signature in the “far dot”, (the one not directly coupled to the MZM) for all gate voltage configurations, while the MZM signature in the dot directly coupled to the MZM can be manipulated via gate-voltage tuning.

Electron-electron interactions will add some interesting effects to this picture. First, there will be the appearance of a Kondo resonance in the dots, which will split due to the indirect exchange between the dots mediated by the leads. More interestingly, we find that the coupling of the dots to the (spin-polarized) MZM will also contribute to the indirect exchange, thus creating a MZM-mediated spin exchange between the dots. These indirect exchange effects are more prominent in the symmetric configuration, where satellite peaks in the spectral density reflect the combined Kondo-Majorana physics at low energies.

An interesting question is how robust are these features with regards to single-particle effects such as Zeeman in the QDs or the coupling of the two MZMs at the ends in the quantum wire. In order to address this issue, we have performed additional calculations to include these effects in the noninteracting limit discussed in Sec. II.1. Our findings show that the overall effect of these terms is similar to the situation in the single QD coupled to a topological quantum wire Liu and Baranger (2011); Vernek et al. (2014); Ruiz-Tijerina et al. (2015). The energy shifts produced by Zeeman terms in the spin-resolved spectral densities are analogous to the changes in gate voltages shown in Figs. 3,4 and 5 such that the qualitative picture remains. As to the effect of the MZM-MZM coupling, we find a result similar to that of Ref. Liu and Baranger, 2011: that spin-up and spin-down spectral densities are similar for energies lower than the coupling energy scale such that the MZM signature disappears.

As both Zeeman splitting and the coupling between MZMs are single-particle in nature, we have no reason to believe that the same qualitative picture will not hold in the interacting regime. The main effect would be a suppression of the Kondo resonance, similar to the case discussed in Ref. Ruiz-Tijerina et al., 2015.

Finally, the results presented here provide a “recipe” for manipulating MZMs signatures present in double quantum dots by using gate voltages and different connection geometries. This opens the interesting prospect of using multi-QD set-ups as a high-control system displaying non-Abelian braiding Malciu et al. (2018), which can readily be integrated in recent architecture designs for scalable topological quantum computers Barkeshli and Sau (2015); Hoffman et al. (2016); Karzig et al. (2017).

Acknowledgements.

The authors thank Edson Vernek, Antônio Seridonio and Eric Andrade for enlightening discussions and suggestions. L. G. G. V. D. S. acknowledges financial support by CNPq (Grants Nos. 308351/2017-7 and 423137/2018-2), and FAPESP (Grant No. 2016/18495-4). J. D. C. acknowledges financial support by CNPq and CAPES.

Appendix A Non-interacting Green’s functions calculation

The spectral representation of the retarded Green function Zubarev (1960) associated to two fermion operators $A(t)$ , $B(t^{\prime})$ is

[TABLE]

Using the equations of motion technique we obtain the following relation Zubarev (1960)

[TABLE]

We apply this expression to Hamiltonian $H$ given by Eq. (1), with $B\equiv d^{\dagger}_{1\downarrow}$ and $A$ varying between the fermion operators $d^{\dagger}_{i\downarrow},f^{\dagger}_{\downarrow},c^{\dagger}_{k\downarrow},d_{i\downarrow},f_{\downarrow},c_{k\downarrow}$. Taking $(\omega,t_{1},t_{2},\epsilon_{1}\ldots)$ as fixed parameters, we obtain a closed linear system of eight equations with eight variables of the form $G_{A,d^{\dagger}_{1\downarrow}}(\omega)$. Hence, this system has a unique solution.

We are interested in computing an analytic expression for $G_{d_{1\downarrow},d^{\dagger}_{1\downarrow}}(\omega)$. The expected solution is a polynomial fraction of degree $8$, whose complexity depends on the number of couplings between the fermion operators. The method described in this paper borrows ideas from graph theory to simplify the Gauss-Jordan elimination process Spielman (2010). We use this method to deduce a simple algorithm to solve the equations of motion of Hamiltonian $H$ of Eq. (1).

Before describing the general procedure, we note that the equations of motion of Eqs. (15) for $A$ equal to $f_{\downarrow}$ and $f^{\dagger}_{\downarrow}$ are

[TABLE]

Since $G_{f_{\downarrow}^{\dagger},d_{1\downarrow}^{\dagger}}(\omega)=G_{f_{\downarrow},d_{1\downarrow}^{\dagger}}(\omega)$ it is possible to eliminate the variable $G_{f_{\downarrow}^{\dagger},d_{1\downarrow}^{\dagger}}(\omega)$ from the system even before starting the Gauss-Jordan elimination.

Writing the remaining EOMs in Eqs. (15) for $A$ varying between $d^{\dagger}_{i\downarrow},c^{\dagger}_{k\downarrow},d_{i\downarrow},f_{\downarrow},c_{k\downarrow}$, we obtain the following linear system

[TABLE]

where $\hat{e_{1}}$ is the vector with entries $\hat{e}_{1_{n}}=\delta_{1n}$, $\mathcal{T}$ is the matrix

[TABLE]

and $\vec{G}_{d^{\dagger}_{1}}$ is the column vector

[TABLE]

The graph associated to the matrix given by Eq. (19) is shown in Fig. 10. Each vertex depicts the first sub-index of the Green function. The values inside each node are obtained by subtracting the corresponding diagonal term from $\omega$. We usually refer to these terms as “self-energies”. The couplings are determined by the off-diagonal terms multiplied by $-1$.

A.1 Solution for a DQD attached to a metallic lead

Before attempting to solve the entire system, we will proceed to explain the Graph-Gauss-Jordan Spielman (2010) elimination process in a DQD-model without Majorana fermions $(t_{1}=t_{2}=0)$. This is equivalent to find the solution for the $3\times 3$ upper-left block matrix given in Eq. (19)

[TABLE]

which can be represented by the graph in Fig. 11(a). In order to eliminate the vertex $c_{k\downarrow}$ we just need to subtract from Eq. (20) the rank-$1$ matrix that cancels the row and the column corresponding to $c_{k\downarrow}$. This matrix is

[TABLE]

The result of Eqs. (20) through (21) is

[TABLE]

which is mapped to the graph in Fig. 11(b).

Note that it is possible to associate the correction to the energies and couplings in Fig. 11(b) to the “walks” passing through the vertex $c_{k\downarrow}$. For instance, $d_{1\downarrow}$’s energy $\epsilon_{1}$ gets an extra-term $\frac{V_{1}^{*}V_{1}}{\omega-\epsilon_{k}}$ representing an additional “walk” from $d_{1\downarrow}$ to $d_{1\downarrow}$ passing through $c_{k\downarrow}$. The terms $V_{1}^{*}$ and $V_{1}$ represent a movement from $d_{1\downarrow}$ to $c_{k\downarrow}$ and vice versa, while the division by $\omega-\epsilon_{k}$ can be thought of as a penalty for passing through $c_{k\downarrow}$. The same logic applies to the coupling terms. The correction to $t_{dots}$ is $\frac{V_{1}^{*}V_{2}}{\omega-\epsilon_{k}}$ which corresponds to a path from $d_{1\downarrow}$ to $d_{2\downarrow}$ passing through the removed vertex $c_{k\downarrow}$. Note that this term includes the multiplication of both couplings with the vertex divided by $\omega-\epsilon_{k}$. This correspondence between the energy correction and eliminated paths through the graph makes this process straightforward.

The next step is to remove the vertex $d_{2\downarrow}$ following the same procedure. At the end, the “self-energy” inside vertex $d_{1\downarrow}$ will be

[TABLE]

and the green function of $G_{d_{1\downarrow}d^{\dagger}_{1\downarrow}}(\omega)$ in a DQD is $\frac{1}{\omega-\epsilon^{+}_{DQD}}$ [see Fig. 11(c)].

A.2 The Graph-Gauss-Jordan algorithm

The previous method to compute the Green function $G_{d,d^{\dagger}}(\omega)$ of an operator $d$ can be summarized in the following steps:

Computing the equations of motion with the second term of the Green function fixed in the creation operator $d^{\dagger}$. The result is a linear system of the form $\mathcal{T}\vec{G}_{d^{\dagger}}=\hat{e}_{1}$ as in Eq. (18). 2. 2.

Mapping the linear system to the associated directed flow graph, labeling the vertices of the graph as $\nu_{n}$, with $\nu_{1}=d$. The self-energy $\epsilon_{n}$ of each vertex $\nu_{n}$ is initialized as $\omega$ minus the corresponding diagonal term $t_{nn}$ of $\mathcal{T}$ $(\epsilon_{n}=\omega-t_{nn})$. The coupling terms $c_{ij}$ connecting two vertices $\nu_{i}$ and $\nu_{j}$ are given by the $(i,j)$-off-diagonal terms $t_{ij}$ of the matrix $\mathcal{T}$ multiplied by $-1$ $(c_{ij}=-t_{ij})$. 3. 3.

Removing one-by-one the vertices of the graph, starting by the last vertex $\nu_{N}$. When a vertex $\nu_{n}$ is removed, an extra-term is added to each energy and coupling. These extra-terms are computed as follows:

1. (a)

Self-energy $\epsilon_{i}$: Let $c_{in}$, $c_{ni}$ be the coupling constants associated to the links from $\nu_{i}$ to $\nu_{n}$ and from $\nu_{n}$ to $\nu_{i}$ respectively. Note that $c_{ni}=c_{in}^{*}$ since the matrix $\mathcal{T}$ is hermitian. Then there is an indirect path from $\nu_{i}$ to itself passing through $\nu_{n}$. When $\nu_{n}$ is eliminated , the extra-term added to $\epsilon_{i}$ is $\frac{c_{in}c^{*}_{in}}{\omega-\epsilon_{n}}$. 2. (b)

Coupling $c_{ij}$ : Let $c_{in},\ c_{nj}\neq 0$ be the coupling constants associated to the links from $\nu_{i}$ to $\nu_{n}$ and from $\nu_{n}$ to $\nu_{j}$. Then there is an indirect path from $\nu_{i}$ to $\nu_{j}$ passing through $\nu_{n}$. When $\nu_{n}$ is eliminated , the extra-term added to $c_{ij}$ is $\frac{c_{in}c_{nj}}{\omega-\epsilon_{\nu}}$.

This process is iterated from $n=N$ until $n=1$. 4. 4.

The self-energy in the remaining vertex $\nu_{1}=d$ is related with the green-function as $\epsilon_{d}=\omega-\frac{1}{G_{d,d^{\dagger}}(\omega)}$.

This algorithm is equivalent to the Gauss-Jordan elimination process used in earlier real-space decimation methods Aoki (1982) in order to obtain noninteracting Green’s functions. Our approach has two additional insights: 1) The number of operations scales linearly with the number of vertices. 2) The graph structure allows one to identify minimal and maximal cutting points which simplifies the complexity of the solution. As pointed out in previous works Spielman (2010), selecting a good order of elimination of the vertices can improve the efficiency of the algorithm. In Fig. 10(a), for instance, it is preferable to start eliminating the vertices at the links, $c_{k\downarrow}$ and $c_{k\downarrow}^{\dagger}$, each one is coupled to just two nodes. Instead, the Majorana operator $f_{\downarrow}$ will be eliminated at last since it is the one with a higher number of couplings.

A.3 Solution for a DQD-Majorana system

From these ideas, we can execute the graph elimination process on the model in Fig. 10(a) . We start by removing the vertices $c_{k\downarrow},c^{\dagger}_{k},d_{2,\downarrow}$, and $d^{\dagger}_{2,\downarrow}$, in that order (See Fig. 10(b)). The energies associated to $d_{1,\downarrow}$ and $d^{\dagger}_{1,\downarrow}$ will be similar to those in Eq. (23), obtaining

[TABLE]

There is also a correction in the couplings between the Majorana mode and $d_{1,\downarrow}$, $d^{\dagger}_{1,\downarrow}$ given by

[TABLE]

In addition there appears a self-energy $\epsilon_{M}$ in the Majorana operator due to the coupling between $f_{\downarrow}$ and $d_{2\downarrow}$. This new term is

[TABLE]

With all the terms of the graph in Fig. 10.(b) computed, it only remains to remove the vertices $d^{\dagger}_{1\downarrow}$ and $f_{\downarrow}$, in that order. This will lead us to the final result (7).

[TABLE]

From this analytic expression we can compute rapidly dynamic quantities such as the density of states in the noninteracting regime. In this project, it allowed us to achieve a better understanding of the system in the different couplings, and also, to predict parameters that exhibit an interesting behavior. These parameters where simulated afterwards through NRG, which has a larger run-time.

We introduced the Graph-Gauss-Jordan algorithm as a simple, didactic and graphical method to solve the equations of motion of quadratic Hamiltonians. We hope for its extended use in condensed matter physics.