Transport signatures of an Andreev molecule in a quantum dot–superconductor–quantum dot setup

Hybrid devices combining quantum dots with superconductors are important building blocks of conventional and topological quantum-information experiments. A requirement for the success of such experiments is to understand the various tunneling-induced non-local interaction mechanisms that are present in the devices, namely crossed Andreev reflection, elastic co-tunneling, and direct interdot tunneling. Here, we provide a theoretical study of a simple device that consists of two quantum dots and a superconductor tunnel-coupled to the dots, often called a Cooper-pair splitter. We study the three special cases where one of the three non-local mechanisms dominates, and calculate measurable ground-state properties, as well as the zero-bias and finite-bias differential conductance characterizing electron transport through this device. We describe how each non-local mechanism controls the measurable quantities, and thereby find experimental fingerprints that allow one to identify and quantify the dominant non-local mechanism using experimental data. Finally, we study the triplet blockade effect and the associated negative differential conductance in the Cooper-pair splitter, and show that they can arise regardless of the nature of the dominant non-local coupling mechanism. Our results should facilitate the characterization of hybrid devices, and their optimization for various quantum-information-related experiments and applications.

are the corresponding energy eigenvalues. To ensure the validity of the perturbative treatment, the relevant, quasiparticle-free, states must have lower energy than the virtual ones with a single quasi particle. For the on-site energy window which we study in this work, namely, −1.5U < ε L , ε R < 0.5U, this requirement translates to ∆ > 3.5U.
Tunneling gives rise to three terms in H eff . The first term is called local Andreev reflection (LAR). It represents processes where the electron number in one of the QDs is changed by 2 (e.g. hybridizes |0, 0 and |0, ↑↓ ), because two electrons forming a Cooper pair in the SC are transferred to the same QD. The second term is called crossed Andreev reflection (CAR). It represents processes where the number of electrons on both QDs change by 1 (e.g., |0, 0 ↔ | ↑, ↓ ), because a Cooper pair is split between the two QDs. The third term is called elastic cotunneling (EC). It represents processes where the total number of electrons on the QDs remains the same, but one electron tunnels between the QDs via the superconductor (e.g., |0, ↑ ↔ | ↑, 0 ).
Note that the EC processes couple the same pairs of states as the interdot tunneling (IT), still EC has different effects on the energy eigenstates and transport properties, because of the pronounced dependence of its matrix elements on the on-site energies, as it will be discussed below.

LAR matrix elements
Here, we derive an exemplary matrix element for LAR, coupling the |0, 0 and | ↑↓, 0 states, in detail. The first tunneling event splits a Cooper pair in SC, one of the electron tunnels to QD L , while the other electron will stay in the superconductor as a quasiparticle. The second tunneling event brings the system back to the low-energy subspace by annihilating the quasiparticle and creating a second electron on QD L . The corresponding matrix element is where (↑↓,0) = 2ε L +U L , and the energy of the intermediate virtual state is ε L + E k . The sum for the momentum k can be written as an energy integral, where ρ 0 is the normal density of states of the SC, assumed to be constant and D is the bandwidth of the SC, and the definitions of the Bogoliubov-amplitudes (see Eq. 4 in the main text) were used.
Generally, this integral cannot be performed analytically. However, in the large ∆ limit, the typical ε L and U L values are small compared to the quasiparticle energy, therefore they can be neglected in the denominators. Then, the integration can be carried out analytically, resulting in With the above approximation, all other matrix elements of LAR have the same form, i.e. Γ LAR,α = S3 −πt 2 Sα ρ 0 , which can be summarized as Note: To check the validity of the approximation applied between Eqs. (S3) and (S4), we have numerically integrated Eq. (S3), and found that the difference between the numerical and approximate values of the matrix element is less than 10% for ∆ > 5U L , and the variation of the matrix element is less than 3% in the −1.5U < ε L < 0.5U window (not shown here).

CAR matrix elements
The CAR mechanism is similar to the LAR, with a difference that the second electron is transferred to the opposite QD as the first one. Therefore, the matrix elements have a similar form, the difference is that both tunneling amplitudes appear in the result. Following a derivation similar to the previous one, we find that the CAR effective Hamiltonian can be written as where Γ CAR = −πt SL t SR ρ 0 .
Within the framework our model, where the spatial separation of the QDs is neglected Γ CAR is linked to the strength of the LAR couplings, i.e. Γ CAR = Γ LAR,L · Γ LAR,R . However taking into account the finite separation of the QDs Γ CAR become suppressed by the distance (see e.g. [2]), hence we treat them as an independent parameters.

EC matrix elements
The EC mechanism transfers one electron from one of the QDs to the other through a quasiparticle state in the SC. Here we show the derivation of the |0, ↑ ↔ | ↑, 0 matrix element, noting that further terms can be obtained similarly. There are two different sequences of underlying processes that build up the EC mechanism. In the first sequence, the electron tunnels from QD R to the SC, where S4 it occupies a quasiparticle state, and then tunnels out to QD L . In the second sequence, a Cooper pair is split, the up-spin part tunneling out to QD L , while the down-spin remains in the SC as a quasiparticle, and then in the second tunneling event the quasiparticle in annihilated by the electron from QD R as they together form a Cooper-pair in the SC. Having two processes leads to the doubling of terms in the matrix element: Using the same approximations as previously, i.e. neglecting ε L/R next to E k , the matrix element vanishes: As this approximation led to a vanishing result, we do the calculation without the approximation.
In fact, the sum in Eq. (S7) can be converted into an integral which can be performed analytically, Taking the infinite-bandwidth limit, D → ∞, this result simplifies to One can further simplify the formula by neglecting ε α next to ∆ in the denominators, yielding (S11)

S5
The main difference compared to the LAR, CAR and IT mechanisms is the dependence of the matrix element on the on-site energies. On the ε L = −ε R line, this EC matrix element vanishes; it is possible that methods using less approximations (e.g., higher-order perturbation theory) yield non-vanishing contributions.
Eq. S11 implies that the strenghts of CAR and EC mechanisms are not independent. We introduce the dimensionless quantity to characterize the strength of the EC mechanism. Even though in our model, and most probably also in an actual experiment, the strengths of CAR and EC are linked, in the main text we treat Γ CAR and γ EC as independent parameters to analyze their effect separately.
With derivations similar to the previous one, the list of the non-zero matrix elements of the EC term of the effective Hamiltonian are as follows: Particle-hole transformations and particle-hole symmetry In the main text, we discussed that the different terms of the Hamiltonian might have particle-hole symmetry (PHS), which can be generated by various particle-hole transformations (PHTs). Here we will show how the individual terms of the Hamiltonian transform under different PHTs.
To analyze the presence and consequences of the PHS of our effective Hamiltonian, it is instructive to rewrite the on-site terms using the normal ordered particle number operators: where : n ασ := d † ασ d ασ − 1 2 is the normal ordered particle particle number operator.H QD differs in two ways from the original H QD . First, the reference of the on-site energy is shifted to 'the center of the odd charge state', or 'particle-hole symmetric point', i.e.,ε α = ε α −U α /2.
All PHTs listed in Table. 1 preserves H INT , but flips the sign of H 0 . Therefore all PHTs introduced here connects the properties at (ε L ,ε R ) to (−ε L , −ε R ), i.e. corresponding to an inversion with respect to (ε L = 0,ε R = 0). Hence the (ε L = 0,ε R = 0) point is called the particle-hole symmetric point.  Our study focuses on a two-site model. One work in which the role of PHTs was elaborated on is Ref. [3], where Sznajd and Becker showed that the generalization of PHT (iii) in Table 1 is a PHS of the usual non-superconducting Hubbard-model on a bipartite lattice.

Additional data on finite bias transport simulation
In the main text we discussed that in general the non-local coupling terms hybridize the Andreev states on the separate QDs, leading to the appearance of anticrossings in the transport spectrum.
Here we will show the similarities and the differences of the anticrossings originating form different coupling terms. Generally two dominating anticrossings open, one on diagonal, i.e. on the ε R = ε L line and one on the skew-diagonal, ε R = −U − ε L , their relative size and their evolution as ε α is tuned differ for each coupling mechanisms.
On Fig. S1 we show the differential conductance of the two side, G L and G R , in the absence of non-  Fig. S2 shows the conductance for the three latter case along the diagonal and the skew-diagonal, to illustrate the evolution of the size of the anticrossings.
In case of CAR coupling (Γ CAR = 0.1U) the more pronounced anticrossing is located on the skewdiagonal at ε L = 0.2U, constituting of 3 conductance lines, appearing along the dashed line on   At ε L = ε R = −1.2U the ground state is dominantly |gs = | ↑↓, ↑↓ + ... and the excited states are |es 1 = |σ , ↑↓ + ... and |es 2 = | ↑↓, σ + .... These two components cannot directly by coupled by CAR, since they have the same number of electrons. However they can be coupled in second order, using the fact that LAR mix these state with |σ , 0 and |0, σ respectively. The CAR matrix element S10 is finite between the |σ , 0 (|0, σ ) and | ↑↓, σ (|σ , ↑↓ ) states. The absence of the direct coupling between the excited states results in the suppression of the size of the anticrossing compared to the one on the skew-diagonal.
The splitting of further conductance lines is comparable to the line width, which is further broadened in an experiment due to the coupling to the normal leads -which is neglected in the present model -, therefore we neglect the detailed analysis of the fine structure here.
In case of EC the anticrossings are positioned similarly as for CAR (see Fig. S1c with γ EC = 0.15), a larger one constituting of 3 conductance line positioned on the skew-diagonal (dashed line) and a smaller one on the diagonal (dotted line). The main difference between EC and CAR lies in behavior of the size of the anticrossing as ε R is tuned. While the pronounced one have the same size for CAR (Fig. S2d), for EC is maximal in the particle-hole symmetric point and it decreases (see Fig. S2e) as ε R is tuned away from this point. For the anticrossing positioned on the diagonal the behavior is strictly different from the case of CAR: On one hand the higher energy anticrossing line has a negative differential conductance (see the black lines on Fig. S2b) for ε L > 0 (negative bias) and ε L < −U (positive bias). Furthermore the two excitation lines move apart as ε L increased from 0 (or decreased below −U) for EC, while the distance between the two lines decreases for CAR (Fig. S2a). Note that second excitation is not visible for opposite bias directions as discussed previously on Fig. S2b.
The third conductance line in the anticrossing positioned on the skew-diagonal is originating from the excitation to the third excited state similarly to the case of CAR, but positioned at higher energy than the first two excitations for EC.
Note that due to the different definition of Γ CAR and γ EC one cannot directly compare the magnitude of the splittings (see Eq. 11 Supplementary Information File 1.).

S11
Finally Fig. S1d shows the effect of the IT, with t LR = 0.1U. On contrary to the CAR and EC, in this case only one anticrossing can be found -along the diagonal (dotted line) -and there is no anticrossing on the skew-diagonal (dashed line). The size of the anticrossing does not depend on the value of ε L , except the ε L ≈ −U/2 region (see Fig. S2c). On the skew diagonal the first two excitations are exactly degenerate, independently of ε R (see Fig. S2f).
As we have shown above the main, common signature of the non-local couplings is the appearance of anticrossings in the excitation spectrum of the QD-SC-QD system, and the main difference of the three couplings lie in the behavior of these anticrossings. For CAR and EC they are positioned similarly, but their size change differently when ε R is changed, while for IT only one anticrossing is present. This difference in the anticrossings allows one to determine the dominant coupling term in an experiment.
A not triplet related NDC line is shown on Fig. S1b at ε L ≈ −0.2U, µ N ≈ −0.3U (marked with a white + sign), where the blocking occurs when the |es 1 → |es 4 transition becomes available, where |es 1 is a singlet, and |es 4 is a doublet state.