Thermoelectric current in a graphene Cooper pair splitter

Generation of electric voltage in a conductor by applying a temperature gradient is a fundamental phenomenon called the Seebeck effect. This effect and its inverse is widely exploited in diverse applications ranging from thermoelectric power generators to temperature sensing. Recently, a possibility of thermoelectricity arising from the interplay of the non-local Cooper pair splitting and the elastic co-tunneling in the hybrid normal metal-superconductor-normal metal structures was predicted. Here, we report the observation of the non-local Seebeck effect in a graphene-based Cooper pair splitting device comprising two quantum dots connected to an aluminum superconductor and present a theoretical description of this phenomenon. The observed non-local Seebeck effect offers an efficient tool for producing entangled electrons.

This supplement provides additional details on experimental and theoretical analysis of our results. We begin with thermal modelling of our Cooper pair splitter device and discuss thermal transport in various parts of the sample, starting from the graphene heater and its electron-phonon coupling, down to small scale thermal gradients between the two quantum dots on SiO 2 (Supplementary Note 1). In Supplementary Note 2, we provide differential resistance data on the SGS superconducting junction thermometers and discuss how we can determine the thermal gradients in our samples using calibrations based on the SGS junctions. Additional thermopower data at finite bias is provided in Supplementary Note 3 in order to stress the importance of having a galvanically separated heater for inducing well-defined thermoelectric currents. Supplementary Note 4 contains details of our theoretical models. We present the theory of thermal transport in a system consisting of two quantum dots coupled to two separate normal leads and one common superconducting lead (cf. Supplementary Fig. 7 for the setting). In Supplementary Note 4 we describe the coherent transport model, which appears to agree well with the experimental results. In Supplementary Note 5, we discuss the incoherent transport regime. This model helps to account for additional sign changes in the local thermoelectric current at low heating power which are seen in the experiments, but not captured by the coherent model.  Fig. 1 displays a schematic heat flow diagram of our graphene Cooper pair splitter device. The heater, the graphene Josephson junctions thermometers (S GS L and S GS R ), and the quantum dots (QD L and QD R ) are all lying on a SiO 2 /Si substrate. We have left out the Al electrode which connects thermally QD L and QD R owing to quasiparticle resistance. According to theory, the important temperatures for non-local thermoelectric effects are those of the left and right normal reservoirs, T L and T R which are formed by the large graphene pieces before the narrow graphene constrictions making the connections to the dots. For the local thermopower also the temperature of the superconductor, T S , is relevant.
Thermal gradient between the two dots is produced by heat flow imposed across the substrate underneath the sample. Establishment of local substrate temperature requires sufficient amount of scattering of phonons, i.e. small enough mean free path for them (see Sect. D 1). In our analysis we assume that the local temperatures T L , T S , and T R can be approximately defined on sub-micron length scales, which is on the order of the separation of the two graphene reservoirs with the SGS thermometers.
The thermal analysis of our system can be divided into four parts: 1) thermal flow in the heater, 2) thermal flow from the heater to the substrate, 3) thermal flow on the SGS thermometer, and 4) thermal flow along the substrate.

A. Thermal flow in graphene heater
The simplest model for the temperature rise in the graphene ribbon heater is provided by the hot electron model. This allows us to estimate the electronic temperature of the graphene electrons using the formula: T e =  The scheme displays electron and phonon systems in graphene (heater, and superconducting junctions S GS L and S GS R ), as well as phonon system of the substrate. Heat is brought in by Joule heating P = V 2 h /R h with R h ∼ 20 kΩ at a typical charge carrier density of n = 4 × 10 11 cm −2 . The electron-phonon coupling in the graphene ribbon heater is expected to be dominated by impurity-assited acoustic phonon scattering (supercollisions) [1], while the phonon-phonon heat transport to the substrate is weakened by Kapitza mismatch at the interfaces. Compared with thermal conductivity of the strongly doped Si, the Kapitza resistance between Si and SiO 2 can be neglected (The Kapitza resistance corresponds to approximately 4 µm layer of our doped Si.) For electrically conducting parts, the electrical heat diffusion governed by the Wiedemann-Franz law is included as relevant. Warm (reddish) colors indicate higher temperatures (close to T e,max = 80 K), while bluish colors denote lower temperatures (close to the base temperature T = 90 mK). The total thickness of the sample is 525 microns (not to scale). temperature distribution over the wire is approximately quadratic with the end temperatures fixed according to the heat balance in the graphene reservoirs of the ribbon.
In order to evaluate the validity of the hot electron model, let us consider the behavior of electronic temperature at the ends of a wire coupled to two wide graphene leads as depicted in Supplementary Fig. 2. Model geometry of the graphene ribbon heater. The heater (size L × w) is connected to two wide, half-infinite-plane graphene reservoirs at the ends. This geometry was employed for the analysis of the relaxation of T e outside the heater section. The radius r 0 determines the border of the region where electronic heat diffusion can be regarded as small.
The temperature both in the ribbon and in the graphene leads can be found from the basic heat transport equation where σ is the conductivity per square, j is the current per unit width ( j = I/w) and the electron-phonon heat flow per unit area equals Σ sup = 7.5 × 10 −4 D 2ñ s k F according to the experimental supercollision results [1,2]; here D 70 eV denotes the deformation potential (in eV) andñ s is the charge carrier density (in units of 10 12 cm −2 ). Using the above numbers, we obtain approximately equal heat fluxes for electron-phonon coupling in the heater and the electronic heat flux at the heater ends. Thus, T max e becomes lowered to ∼ 50 K.
Outside the narrow section of the heater, j = I πr in the graphene lead at a distance r from the end of the ribbon. At large bias, electron-phonon coupling dominates and the temperature in the graphene near the heater ends equals to A characteristic length scale for the decay of the heat flow along graphene reservoir can be obtained by equating the local variation of the gradient term in Supplementary Equation 1 to the change in heat flux due to electron-phonon coupling. The characteristic thermal relaxation length scale L T (T S ) becomes At T S = 1 K, we obtain L T = 2.3 µm using a carrier concentration of n s = 4 × 10 11 cm −2 which is close to our experiments at a back gate voltage 5 V. Since this length decreases with temperature, it is evident that T e decreases quickly in the leads and that most of the heating power is deposited to the substrate within and near the heater. The Kapitza resistance between SiO 2 and graphene will enhance T of graphene phonons upto ∼ 2 K, when the electronic temperature is 10 K. This will change slightly the above relaxation length, but it is still the supercollision cooling with coupling Σ sup (T 3 e − T 3 ph ) that governs thermal relaxation of the electrons.
B. Thermal flow from the heater to the substrate On the basis of the smallness of the thermal relaxation length in Supplementary Equation 3, we assume that most of the heater power is deposited to the substrate within the area of the graphene ribbon heater and its immediate neighborhood. Using 2 µm 2 for the thermal contact area, A K = 500 Wm −2 K −4 for Kapitza conductance (P = 1 4 A K (T 4 h − T 4 s ) with temperatures T h and T s on the opposite sides of the interface) [3], we obtain a phonon temperature of T h ∼ 9 K in graphene at 40 mV heating voltage. This elevated phonon temperature will relax along the substrate via phononic thermal conductance.

C. Thermal flow on the SGS thermometer
The SGS thermometers are formed by two-terminal, Al/Ti-contacted graphene parts in the center of large pieces of graphene. They are made in such a way that there is an uninterrupted path of non-superconducting graphene connecting the thermometer to a quantum dot. This guarantees smooth thermal flow along the thermometer area without any strong suppression in the heat conductivity due to proximity-induced superconducting gap (note that the graphene electrode is fully crossed by a superconducting lead at the further end of the thermometer, see Fig. 1 of the main paper). The purpose of these SGS devices is to track the temperature of the electronic graphene reservoirs that govern the distribution of incoming electrons/holes on the graphene quantum dots. The large area of the SGS devices guarantees that their temperature will well track the temperature of the substrate, independent of a possible heat input coming along the graphene ribbons from the dots.
The operating range of our SGS thermometers is up to 0.8 K, which facilitates thermopower studies up to rather substantial thermal gradients. At a typical thermometer temperature of 0.5 K, we obtain a characteristic relaxation length of 3.3 µm using Supplementary Equation 3. This is larger than the extent of thermometer along the heat flow direction. Consequently, we will take a spatially independent temperature for the SGS thermometer. In order to determine how the SGS thermometer averages over the spatially dependent substrate temperature T s (x), we assume that the heat flow balance between SGS and the substrate remains zero. This yields the integral condition (T 3 S GS − T 3 s (x))dx = 0. At low temperature with small dT 3 s (x)/dx, one may obtain a simple estimate for T S GS for the temperature recorded by the thermometer; here T l and T r denote the temperature at the left and right edges of the thermometer region, respectively. D. Thermal flow along the substrate

Local phonon temperature and heat diffusion
For application of the regular heat diffusion equation, small mean free path of phonons is essential. In clean silicon, the mean free path of phonons can be on the order of 100 microns, which would lead to problems in defining a proper substrate temperature across our CPS device having a size smaller than 10 µm. However, the interfacial scattering at the SiO 2 /Si interface and the strong impurity doping of the Si ++ material decrease the mean free path of phonons substantially. Typically in case of thin layers and rough surfaces at an interface, the phonon mean free path is limited to a value of the order of layer thickness [4] (Casimir limit). Thus, in our case ph ∼ 300 nm in the SiO 2 that is in thermal contact with the CPS device. Owing to the relatively small ph , we argue that the heat diffusion equation can be employed to deduce approximate temperature distribution along the substrate across the area of the graphene reservoirs of our Cooper pair splitter.

COMSOL simulations for temperature difference between the dots
The heat diffusion equation was solved using COMSOL multiphysics. The geometrical model used in the simulations (see Supplementary Fig. 3) is a simplified version of the real device presented in Fig. 1 of the main text. Our simulations using exact sample dimensions focused on the temperature profile on the surface along the line drawn through the heater and both quantum dots. Best available data on materials parameters were employed in calculating thermal gradients along the substrate. The graphene heater power was inputted only over the actual heater size, which slightly increased the maximum temperature of the SiO 2 but this was deemed irrelevant further away from the heater which was the main region of interest.
In order to achieve most realistic temperature profile, we employed temperature dependent heat conductivities both for the phonon thermal conductance as well as for the metallic heat diffusion. The mean free path of phonons influences strongly the thermal conductivity (see Sect. D 1). We employed κ S iO 2 = 0.01 × T 2 W/(mK 3 ) [5], κ gr = 0.02 × T 2 W/(mK 3 ) [6,7], and κ S i ++ = 0.01(T/K) 2 W/mK 3 [8] for thermal conductivities of SiO 2 , graphene, and strongly doped silicon, respectively. For dynamical calculations, we also specified heat capacities C S iO 2 = 3 × 10 −3 × T 2 J/(kgK 3 ), C gr = 1 × T J/(kgK 2 ), and C S i ++ = 10 −3 J/(kgK 3 ), respectively. Supplementary Fig. 3 displays simulation results for the spatially dependent temperature T (x) on the surface of the SiO 2 layer as a function of the distance x from the heater along the central symmetry axis intersecting the two quantum dots and the center of the graphene ribbon heater. The simulated trace at V h = 40 mV yields ∆T = 72 mK for the temperature difference between the two red lines marking the edges of the SGS graphene thermometers. The flatter regions before and after the steep section between the red lines is due to the large phononic conductance of graphene. Note that the regions would be even flatter if the electron-phonon coupling and the electronic heat conductance κ e in graphene would have been included. This omission of κ e also renders the modification of electrical thermal conductivity by the induced superconducting gap irrelevant for the simulation.
The COMSOL simulation was also used to estimate the thermal time constant of the system. This was done by looking at the phase difference between the heater temperature T h and T QD L as the heater modulation frequency was increased. The results are presented in Supplementary Fig. 3d at frequencies 0.1, 1, and 2 MHz. The simulation demonstrates that at 100 kHz the phase difference is still negligible, while at higher frequencies it starts to grow. Thus, we can conclude that the time constant τ < 0.01 ms and steady state approximations can safely be applied to our measurements at f < 10 Hz. Once calibrated, the heater voltage could be used for specifying the temperature difference across our Cooper pair splitter. This was useful, in particular because of the poor resolution of the S GS R thermometer, which imposed long measurement times for good accuracy. The poor resolution of S GS R was due to lack of measurement wires in the cryostat, which forced us to connect the S GS R thermometer in a 2-wire configuration. The temperature of S GS L w.r.t. the heater voltage was calibrated by measuring the IV characteristics of the S GS L as a function of the heater voltage V h and the equilibrium cryostat temperature T cryostat . The obtained differential resistance R d characteristics are presented in Supplementary Fig. 4, in which the decrease of the superconducting gap with increasing V h and T cryostat can be seen clearly. At the edge of the gap, we traced the maximum differential resistance, R peak d = max R d (I DC ), both as a function of V h and T cryo . The observed traces can be made to overlap each other by adjusting the x-scale, i.e. by mapping the voltage scale to a corresponding temperature scale. The obtained scaling relation for S GS L is given by: where V h is in volts. Note that there is an offset of 90 mK that we interpret as the effective electron temperature in graphene at the cryostat base temperature when no heating is applied. For S GS R , we obtain the relation that would in turn be consistent with the cooling P ∝ T 3 being dominated by acoustic phonons, for which the thermal conductivity κ ph ∝ T 3 .
As mentioned in the Methods section, the length of our SGS thermometers is intermediate between the short and long junction limits. Unfortunately, no analytic formulas exist for this regime. For modeling purposes, we have employed the low temperature expression for the critical current of long junctions I c = E th b(1-1.3 exp [−bE th /3.2k b T ])/eR n , using b as a fitting parameter; here R n is the normal state resistance and k b denotes Boltzmann's constant. Using b 2.2, this formula agrees approximately with the measured behavior of I c (T ) and the I c R n product at the lowest temperatures. The fit also corroborates with the observed optimum thermometer sensitivity around T 0.3 K, with tendency towards saturation below T < 0.2 K. However, the thermometer works even below 0.2 K as the maximum differential resistance in the sub gap region still grows with decreasing temperature below T = 0.2 K. The absence of clear hysteresis of the SGS junctions is assigned to the built-in shunts provided by the surrounding, non-proximitized graphene.
The calibration between heating-induced temperature and cryostat temperature is bound to the site where the SGS junctions are located. It is true that our thermometers are some distance apart from the quantum dots. Thus, there may be small thermal gradients along the graphene part in contact with the thermometer. The gradient depends on the ratio of electron-phonon/Kapitza conductance to the substrate (as discussed in relation to Supplementary Fig. 1) and the electronic conductance along graphene. At the lowest temperatures, the electronic conductance will dominate and we may neglect the thermal gradients along the graphene flake. At higher temperatures, thermal gradients will appear and the measured temperature difference will be an upper bound for the temperature difference. However, if we plot the local thermoelectric data for the peak 4 in Fig. 4a of the main text as a function of temperature difference ∆T = T L − T S , we find a linear dependence with ∆T up to the largest heating powers employed in the experiment, as demonstrated by the data in Supplementary Fig. 5. We consider this as evidence that thermal gradients are not disturbing our calibration significantly at those temperatures relevant for our analysis. The thermoelectric current measured in the main paper was purely driven by temperature gradient, without any bias applied on the Cooper pair splitter. We also performed experiments where the thermoelectric response was investigated without a separate heater. Bias current in the graphene ribbon connecting the QDs to the graphene reservoirs leads to heating which creates basically local thermopower phenomena with symmetry properties close to the non-local phenomena. Because of this, the central figures in the main paper (Figs. 3 and 4) contain only data with plain thermal drive. Supplementary Fig. 6 shows results on thermoelectric current in a Cooper pair splitter driven by a voltage bias on one of the quantum dots. In this case, thermoelectric effects and the bias-induced CPS and EC will mix together in the behavior. Supplementary Fig. 6a and 6b display the currents of the two dots I L and I R measured at ω as function of the dc bias V b and the side gate voltage V sg,L on QD L ; here dc and ac voltage bias is applied only on QD R with the middle superconcucting lead grounded. The current I L shows different sign with positive and negative dc bias V b on QD R . In addition, I L changes its sign when V sg,L crosses V sg,L = 0.72 V, at which the energy level of QD L matches the Fermi energy. Ac current I R , however, is insensitive to the V sg,L .
Figs. 6c and 6d display I L and I R as functions of dc bias V b on QD R and the side gate voltage V sg,R . Here I L changes its sign, when V b crosses zero, but it retains its sign when V sg,R crosses V sg,R = 0.535 V, at which the energy level of QD R matches the Fermi level. Note that the observed variation of I L as function of V sg,L , V sg,R and V b can be explained by considering heating due to V b and local thermoelectric currents induced by the ensuing thermal gradients between QD L and the superconductor. Unfortunately, there are no means to single out the non-local thermoelectric effects in these sets of data.
If QD R is ac biased at ω without any dc, we obtain Supplementary Fig. 6e for the currentĨ L in QD L at 2ω. The thermoelectric current in QD L changes its sign while the energy level tuned by V sg,L crosses the Fermi energy. On the contrary,Ĩ L keeps its sign as the energy level is tuned by V sg,R across the Fermi surface. These results are consistent with those in Figs. 6a-d.
Although the observed data in Supplementary Fig. 6 are similar to those in Figs. 3 and 4 in the main paper, and apparently non-local effects according to this similarity in symmetry and asymmetry, they can be explained using bias-induced heating and local thermoelectric effects induced by thermal gradients. In Figs. 6a and 6c, the dc and ac bias applied to QD R causes heating in the middle Al lead. When the dc bias is positive, the heating in the Al lead will be in the same phase with the ac bias of QD R : i.e. the heating in Al follows the maximum and minimum of the ac bias on QD R . When the dc bias is negative, however, the maximum heating will coincide with minimum of the ac, and the heating in the middle Al will be phase shifted by π from the ac bias on QD R . Hence, in Figs. 6a and 6c, the thermoelectric current I L has an opposite sign at positive or negative dc bias. In Supplementary Fig. 6a, I L changes its sign due to the same reason as in Fig. 2a of the main text when V sg,L crosses V sg,L = 0.675 V. Because the heating in the middle Al lead is equal at left and right fo the conductance peak of QD R , I L keeps the same sign while energy level of QD R crosses Fermi level in Supplementary Fig. 6c.
Even when using the 2ω method, driving QD R at ω, measuring I L at 2ω, the same heating argument can be employed to explain the observed bahavior in Supplementary Fig. 6e. In the absence of bias voltage, the heating is now at 2ω, at which the signal then appears. These data have strong resemblance with the result in Fig. 4d in the main paper. The current I R at ω is depicted in Supplementary Fig. 6f.
Ac current of QD R as a function of V sg,L while sweeping dc bias on QD R at V sg,R = 0.535 V. c Ac current of QD L as a function of V sg,R while sweeping dc bias on QD R at V sg,L = 0.72 V. d Ac current of QD R as a function of V sg,R while sweeping dc bias on QD R at V sg,L = 0.72 V. In all of these measurements a-d, QD R was also biased using 20 µV rms at 2.1 Hz for ac measurements. e Apparent non-local thermoelectric effects due to variation in Joule heating across the resonance. CurrentĨ L (at 2ω) of QD L measured over V sg,L and V sg,R plane using an ac bias of 50 µV on QD R . f The corresponding current I R of QD R measured at ω over V sg,L and V sg,R plane. Notice the resemblance of frame e with the non-local thermoelectric current patterns in Fig. 3 of the main paper.
Here, the large ac excitation leads to gate voltage dependent heating in QD R , which leads to variation of T with the same symmetry properties as the non-local thermoelectric effect with constant thermal gradient.  Fig. 7. Schematics of the N-dot-S-dot-N system. Two quantum dots with the energy levels ε L,n and ε R,n are coupled to the common superconducting lead with the gap ∆ and the temperature T S , and to the two different normal bulk leads with the temperatures T L (left lead) and T R (right lead). The tunneling rates between the dots and the superconductor are γ L,n and γ R,n , while the tunneling rates between the dots and the normal leads are Γ L,n and Γ R,n . Bias volatges V L and V R can be applied to the outer leads, and the side gate potentials V sg,L and V sg,R allow one to shift the energy leves ε L,n and ε R,n . The potential of the superconductor is supposed to be zero. The electric currents I L and I R flow through the quantum dots into the superconductor.
In this note we outline the theory of thermal transport in the system consisting of two quantum dots coupled to two separate normal leads and one common superconducting lead. The schematics of it is shown in Supplementary Fig. 7. Here, we describe the coherent transport regime, which rather well describes the main findings of our experiment. The additional sign changes in the local thermoelectric current at low heating power which are not captured by the present description are explained within the incoherent model discussed in Supplementary Note 5.
We assume that the conductances of the dots are high enough so that the dots enter the Fabry-Pérot regime, and one can ignore Coulomb blockade. We split the energy levels in both dots into two groups: well coupled to the leads levels ε j,n (here the index j = L, R indicates the left and right dot, respectively, and n enumerates the energy levels within one dot), and poorly coupled to the leads levels ε j,n . The hopping rates of electrons between the levels ε j,n and the normal leads, which in our experiment are big graphene flakes, are denoted by Γ j,n , while the hopping between the dots and the superconducting lead is set by rates γ j,n . The energy levels can be tuned by side gate potentials V sg, j as ε j,n (V sg, j ) = a j (V sg, j − V max, j,n ), ε j,n (V g j ) = a j (V sg, j − V j,n ).
Here the constants a j are determined by the ratios of the gate and the dot capacitances, and V max, j,n are the gate voltages at which the conductance of the corresponding dot exhibits a peak and reaches a local maximum, and V j,n are the gate voltages at which the energies of the dark states ε j,n become equal to zero. From the experiment we approximately find Bias voltages V j are applied to the metallic leads, while the potential of the superconductor equals zero.
A. Local electric current in the normal state In this section, we recall some well known results and consider fully normal system assuming ∆ = 0. For a nice general review on thermoelectric effects in normal three terminal systems we refer the reader to Ref. [9]. Since we neglect the Coulomb blockade, transport of charge through the quantum dots is well described by the Landauer formula, where the superscript loc indicates the local character of the current (i.e. I loc j can be evaluated assuming that the other quantum dot does not exist), h is Planck's constant, τ j (E, V sg, j ) is the transmission probability of the dot, are the electron distribution functions in the bulk normal lead j and in the superconductor, V j is the electric potential applied to the lead j, T j is the temperature of the lead j, and T S is the temperature of the central electrode. The factor 2 in front of Supplementary Equation (7) accounts for the spin degeneracy. We use two models for the transmission probability: (i) simple resonant tunneling model, in which τ j (E, V sg, j ) is given by a sum of Lorentzian peaks, each coming from one of the energy levels, and (ii) a model with Fano resonances, in which Here t j,n is the hopping amplitude between the energy level ε j,n , which is coupled to the leads, and to one of the closely lying uncoupled (dark) levels with the energy ε j,n . Of course, for some peaks one can have t j,n = 0, and in this case the transmission has the same Lorentzian shape as in Supplementary Equation (9). According to Supplementary Equation (7), the zero-bias dimensionless conductance of the dot in the zero temperature limit has the form In the limit k B T j , k B T S Γ j,n + γ j,n and for V j = 0, Supplementary Equation (7) reproduces the Mott formula for thermoelectric currents The derivatives of the transmission probabilities over the energy E can be converted into the derivatives over the gate potentials if one uses Supplementary Equations (9,10) in combination with the gate efficiency relations in Supplementary Equation (6). In this way, one finds The corresponding thermopower, or the Seebeck coefficient, is given by Supplementary Equations (13,14) allow us to estimate the expected values of thermoelectric currents and the thermopower once the dependence of the zero bias conductance on the gate voltage has been measured. Supplementary Equation (13) also determines the sign of the thermoelectric current and the values of the gate voltage at which it equals to zero.
B. Local electric current in the superconducting state In this section, we assume that the central lead of the system shown in Supplementary Fig. 7 is superconducting. In this case, the local transport properties of the individual dots are described by the theory of Andreev reflection at NS interfaces [10,11], which needs to be generalized to the case of energy dependent transmission probability. For the resonant tunneling model, such generalization has been presented, for example, in Ref. 12 which is convenient for the analysis of the experimental data. Here we have introduced the re-normalized energyẼ and the density of states in the superconductor ν S (E), We have also replaced the hopping rates of the individual levels Γ j,n and γ j,n by the gate-voltage-dependent expressions Γ j (V sg, j ) and γ j (V sg, j ). In fact, Supplementary Equation (15) provides a full description of the local transport through the quantum dots, and it can be used beyond the resonant tunneling model. For example, Supplementary Equation (15) accounts fully for the contribution of Andreev bound states, which are formed in the dots due to their coupling to a superconducting lead, and which have been observed in InAs [13] and in carbon nanotube quantum dots [14]. If the transmissions τ j do not depend on the energy and Γ j ∆, Supplementary Equation (15) reduces to well known result by Blonder Tinkham and Klapwijk [11]. Thermoelectric properties of a quantum point contact in this regime have been recently studied in detail in Ref. [15]. In the limit ∆ = 0, Supplementary Equation (15) takes the Landauer form (7). Finally, according to Supplementary Equation (15) in the low temperature limit, k B T j , k B T S ∆, Γ j,n + γ j,n , the zero bias conductance of the dot j acquires the form which, as expected, coincides with the prediction of the theory of Andreev reflection at low energies [10].
For an asymmetric quantum dot with low transmission, τ j (E) 1, and for eV, k B T j , k B T S , Γ j , γ j ∆, Supplementary Equation (15) reduces to a simple expression which contains the product of the transmission probabilities for an incoming electron, τ j (E, V sg, j ), and for a reflected hole, τ j (−E, V sg, j ). Supplementary Equation (17) in combination with the Lorentzian transmission probability of Supplementary Equation (9) has been demonstrated to reproduce accurately the shape of the conductance peaks in carbon nanotube quantum dots [14,16]. We find that this model is also working quite well for our graphene quantum dots. In Supplementary Fig. 8a,b we display the zero-bias conductance of the left and right quantum dot, measured at T L = T R = T S = 90 mK, as a function of the side gate voltages V sg,L , V sg,R . Inverting the expression (17) for zero bias conductance, we obtain normal state transmission probabilities in the form Applying this transformation to the experimental data sets depicted in Figs. 8 (a,b), we get the normal state transmission probabilities at zero energy τ j (0, V sg, j ). The latter are plotted in Figs. 8 (c,d) with blue lines. In the same figures, by red dotted lines we also show the fits of the transmissions τ j (0, V sg, j ) with multiple Lorentzian peaks in accordance with the predictions of the resonance tunneling model (Supplementary Equation (9)). The quality of these fits is rather good, especially for the left quantum dot. We conclude, therefore, that the resonance tunneling model describes the quantum dots well. A closer look at some of the conductance peaks reveals deviations from the form predicted by Supplementary Equation (17) with the Lorentzian transmission probabilities (9). In Supplementary Fig. 9 we plot two representative peaks, one for each dot, and indicate by the green arrows the features not captured using the resonant tunneling model. These features may be caused by various reasons. For example, they may result from Coulomb interaction, inelastic relaxation processes, or Kondo effect, etc. They may also arise from overlapping of the relevant energy levels. Here we fit these features using the Fano resonance model (10), in which a conductance peak splits into two closely lying peaks with different heights and widths. We believe that this is the most plausible explanation, although we cannot fully exclude other options. This uncertainty, however, does not affect the expression for the non-local current given below in Supplementary Equation (25) and in the main text. Indeed, the latter contains only the effective transmission probabilities τ j , which can be inferred from the conductance data, and, therefore, are not very sensitive to a specific model describing the features of τ j (E, V sg, j ) and its parameter dependence. With the Fano resonance model, we have fitted the conductance peak of the left quantum dot quite accurately, see Fig 9a. However, for the right dot we have adopted a more pronounced Fano resonance than the conductance data would suggest (Fig 9b) in order to reproduce extra sign changes of the non-local thermoelectric current in the vicinity of this conductance peak (cf. Fig. 4 in the main text).
Having determined the zero-energy transmission probabilities of the quantum dots, τ j (0, V sg, j ), as functions of the gate voltages, we make the replacement V g j → V g j + E/α j , where E is the energy of an electron, and in this way recover the full energy dependence of the transmission probabilities. Then, we use the obtained energy-dependent transmission probability τ L (E, V sg,L ) to calculate theoretical thermoelectric currents, based on Supplementary Equation (15), for different heating voltages. The result, displayed in Fig. 4a of the main text, agrees rather well with the experimental findings apart for the magnitude of the current. We would like to emphasize that, in this model, only quasiparticles with the energies |E| > ∆ contribute to the thermoelectric current. That is why the predicted values of the current at low heating powers, I j ∝ exp(−∆/T j ), are very small. In the experiment, however, we observe significant thermoelectric currents even at the lowest heating voltages 1 mV < V h < 5 mV. This observation may be explained, for instance, by significant heating of the superconducting lead even at small heating voltages. Alternatively, it may be a result of the induced thermal voltages, which can be explained by the incoherent tunneling model discussed in Sec. C.
For completeness, we provide the low temperature expressions for the local thermoelectric currents, which follow from the general expression of Supplementary Equation (15) for V L = V R = 0 and in the limit T L , T R , T S ∆/k B , C. Non-local contributions to the thermoelectric current The main focus of our study is Cooper pair splitting (CPS). It results from the crossed Andreev reflection process, in which a single Cooper pair splits into two electrons injected into different quantum dots. It has been shown (see e.g. Refs. [17][18][19][20] that the probability of such a process is proportional to the product of the transmission probabilities of the two quantum dots, τ L,R (E), and the effective transmission probability of the superconducting electrode τ S , Strictly speaking, Refs. 17-20 deal with metallic normal leads, in which the transmission probabilities τ L , τ R do not depend on the electron energy E. Correct energy dependence of these transmissions can be obtained by comparing the crossed Andreev reflection probability of Supplementary Equation (21) with the probability of ordinary Andreev reflection, in which both electrons are emitted in to the same quantum dot. The latter probability appears in Supplementary Equation (18) for the Andreev current through a low transmission quantum dot and contains the product τ j (E)τ j (−E), in which an electron and a reflected hole have opposite energies. Clearly, the same rule should apply to the crossed Andreev reflection, and in this way one arrives at the symmetry of Supplementary Equation (21). The effective transmission probability of the superconductor, τ S , does not depend on the energy E if the quantum dots are placed close to each other and the distance between them is less than the coherence length of the superconductor. This condition is approximately satisfied in our experimental configuration. The parameter τ S has the meaning of the probability for an electron emitted from one of the dots to reach the other dot instead of flying away into the bulk of the superconductor. From the theory of disordered superconductors based on the Usadel equation one formally finds [19] where R S is the resistance of the superconducting lead in the normal state, d is the distance between the two dots, ξ S = √ D/2∆ is the coherence length inside the superconductor and D = v F L e /3 is the diffusion constant. Supplementary Equation (22) is valid if the electron mean free path l e is much smaller than the distance between the dots. In practice, the value of τ S is sensitive to the sample geometry, granularity of the superconductor, quality of the contacts between the superconductor and the quantum dots, etc. Therefore, we treat it as a fit parameter. We find that τ S = 0.1 provides good fit of our data.
It is well known that in addition to the Cooper pair splitting process, elastic cotunneling also contributes to the non-local transport. The probability of this process is given by the product [17][18][19][20] which differs from τ CPS only by the replacement E → −E in the argument of τ R . The approximate expressions in Supplementary Equations (21,23) are valid provided τ CPS , τ EC 1. This condition is well satisfied in our experiment. The currents flowing through the quantum dots are given by the sum of local and non-local contributions, where the local currents are given by Supplementary Equation (15), and the non-local corrections in the low temperature limit k B T L , k B T R ∆ and zero voltage drops across the dots, V j = 0, are expressed as One can also define the linear combinations of the non-local currents, which are determined solely by Cooper pair splitting or elastic cotunneling, These expressions presented and discussed in the main text. In order to understand the overall behaviour of the non-local currents, we again consider the low temperature limit k B T L , k B T R ∆, γ j,n + Γ j,n . In this limit one can derive the expressions analogous to Mott's formula (13), Supplementary Equation (27), for example, shows that the non-local contribution to the thermoelectric current of the left quantum dot should change sign at the values of the gate voltage V sg,L corresponding to the positions of the conductance peaks and the minima of the conductance valleys of the left dot. At the same time, no sign change in ∆I nl L occurs if one varies the gate potential of the right dot V sg,R . The current ∆I nl L depends on the right gate potential V sg,R roughly in the same way as the conductance of the right dot g R (V sg,R ) does. Thus, in order to restore the sign pattern and relative magnitudes of the thermoelectric currents one can use a simple rule ∆I nl L ∝ (dg L /dV sg,L )g R , ∆I nl R ∝ g L (dg R /dV sg,R ). Such Mott-type behaviour of the non-local currents agrees well with gate voltage dependence observed in the experiment. Similar sign pattern is predicted by a slightly different model, in which a quantum dot is coupled to two normal and superconducting terminals [21]. Res RR are the electron-reservoirs, Dot L(R) is the left (right) quantum dot, R nr L and R nr R are the additional resistors, S is the superconducting region, δ L(R) is δ-barrier on the left (right) NS-interface. Res LL , Res RR and S are grounded. b Side view (schematic) of the graphene-aluminum (superconductor) junction. The patterned graphene piece has a characteristic size of ∼ 300 nm and constitutes a quantum dot. Here L and d are the lengths of the graphene parts which, respectively, are either overlapping or extending beyond the superconductor. c Schematics of the theoretical model accounting for the effect of the additional Fano resonance structure. The quantum dot is modelled by two δ-function potentials (the distance between barriers is d) connected to a fork with a stub of length L e ∼ L 2 /λ. The transmission and reflection amplitudes of the beam splitter (blank square) are denoted by t (2) and r (2) ; the amplitudes incorporating the internal features of the fork and the stub are denoted by t (1) and r (1) . We set the distance a between the double barrier and the beam splitter equal to zero.
The model described in Supplementary Note 4 is based on the assumption that the transport in the Dot-S-Dot system is fully coherent. In reality, the system may be subject to dephasing and relaxation. To demonstrate how the decoherence processes may affect the experiment, we present and analyze an alternative theoretical model based on the scattering matrix approach. This model is predicated on the assumption that the system may be split into the coherent subsystems which, in turn, are joined incoherently as connected circuit elements. For this reason, thermal gradients induce both electric currents and voltage drops on the dots. That is why, in contrast to the coherent transport model, the incoherent one predicts significant value of the local thermoelectric currents at low temperatures, where the quasiparticles in the superconducting lead disappear, and, thus, it may explain the experimental observations in this regime. However, our analysis shows that decoherence has little effect on the measured non-local phenomena.
To account for possible decoherence, we assume that the schematics includes additional reservoirs (Res LR and Res RL ) between the superconductor and both the left and the right dot as depicted in Supplementary Fig. 10a. Similar models of decoherence have been considered elsewhere, see e.g. Ref. [22]. We model the whole structure as a one-dimensional conducting structure with ballistic motion of electrons in the superconductor and in the quantum dots. The superconductor (S) is separated by δ-barriers (δ L and δ R ) from the adjacent reservoirs. The transmission probabilities of the quantum dots (Dot L and Dot R ) may exhibit Fano resonances. In order to model such resonances, the dots are assumed to be composed of two elements -Fabry-Pérot double barrier structure and a stub. Finally, we assume that the leftmost and rightmost reservoirs (Res LL and Res RR ) along with the superconductor are grounded.

A. Local electric current
Let us fist discuss the characteristic features of the experimental plots showing the dependence of the local thermoelectric current on the gate voltage in the quantum dot. For the present model we assume that the local thermoelectricity in the left (right) side of the structure is governed by the difference between the population distributions in Res LL and Res LR (Res RL and Res RR ) due to the temperature gradient. For simplicity, we will consider only the subgap regime.
As mentioned in the main paper, we can explain the extraordinary behavior of the local thermoelectric current (i.e., the appearance of a secondary extremum on one side of the main inflexion point) by the Fano resonant effect emerging in the connection between the dot and the superconductor. The arrangement of the experimental setting (as shown in Supplementary  Fig. 10b) suggests that it can be modeled by the scheme depicted in Supplementary Fig. 10c. Here the interface between the superconductor and the double delta barrier (the distance between barriers is d) is considered to have a fork structure with one of its contacts being a stub of effective length L e . L e is much larger than d which ensures the interaction between the discreet spectrum due to the stub and the semi-continuum in the double delta barrier (which encompasses a large level width). Drawing a parallel with the schematics in Supplementary Fig. 10b, the stub relates to the part of graphene overlapping with the superconductor, which has length L. The length of the stub L e , however, should not necessarily be equal to L: since the mean free path λ of the particles is smaller than L, in reality they may undergo multiple reflections inside this part of the structure. Then where D ∼ v F λ is the diffusion coefficient and t t is the time that the particle spends in the stub. Accordingly, the effective length of the stub L e ∼ v F t t ∼ L 2 /λ used in the model can exceed the typical size of the dot L 300 nm. We assume that the barrier in the end of the stub completely blocks the propagation of the particles, i.e., the probability density for finding the particle beyond the barrier is zero. For clarity, we separate the right delta barrier from the beam splitter, but imply that the distance a between them is zero. If L e is finite, the appearance of Fano resonances affects the transparency between left and right terminals changing the behaviour of the thermoelectric current.
is the wave vector of the particle inside the double barrier (here, V g is the voltage on the gate).
Fork.-Supposing that the wave function of the particles goes to zero on the boundary in the stub and the stub is grounded (i.e., it is not affected by the gate voltage on the dot), we may obtain the transmission and reflection amplitudes t (1) , r (1) on the contact 1-3 (incorporating the internal features, i.e., the stub, see Supplementary Fig. 10c): where index (2) denotes the transmission and reflection amplitudes of the beam splitter, and L e is the length of the stub; we assume that the wave vector k = k F + E v F in the stub does not depend on the gate voltage. For the numerical calculations we will parametrize the beam splitter's scattering matrix with only two parameters: the elements of the matrix can be expressed in terms of the transmission probabilities in contacts 1-3 and 2-3: where Whole structure.-The final transmission probability τ = |t| 2 of the dot can be calculated by substituting Supplementary Equations (30-36) into the formula The scattering on the outer (O) and inner (I) points of the structure is described by the reflection amplitudes In the experiment, the local conductance and current are measured in the network of the dot along with the graphene nanoribbon (represented by R nr L(R) in Supplementary Fig. 10a) and the N-δ-S boundary. For the theoretical model, we should assume that the nanoribbon, together with the N-δ-S boundary, have a constant electrical resistance R. This allows us to write the following equation for the local current through the left dot (the transparency of the dot τ L is calculated as discussed above): where V denotes the external voltage applied to the system (we will need V later to calculate the differential conductance), T LL(LR) is the temperature of Res LL(LR) and α (0 ≤ α ≤ 1) describes the voltage distribution across the dot. Since δT L = T LL − T LR T LR and eI loc L R k B T LR (as one will be able to see later), we obtain whence The parameters of our model are such that the dependence of τ L (E) on E is relatively weak. We can therefore write the following equations for the current I loc L (V = 0) and the differential conductance g loc L = ∂I loc L ∂V : 13 , the behaviour of the thermoelectric current predicted theoretically compares well with the experimental data. Moreover, we should note that the parameters of the fits are such that the length L e = 1810 nm comply with our prediction of ∼ L 2 /λ ∼ 2000 nm given that L 300 nm and λ ∼ 30 nm. Thus, the above model, accounting for quasi-particles in the superconductor at elevated temperatures, allows for a faithful theoretical description of the measured data. However, the model does not reproduce the experimental observation that the secondary features of the curves tend to disappear at large heating voltages.

B. Non-local electric current
The non-local currents in this model are contributed by two factors: • The voltage difference between Res LR and Res RL which, in turn, is induced by the local thermoelectricity in the left and right leads. For instance, the temperature difference between Res LL and Res LR generates a nonzero electric potential in Res LR .
• The temperature gradient between Res LR and Res RL .
Our numerical analysis shows that the latter factor does not play a significant role and can be neglected. Physically, this is due to the fact that the transmission through the δ L -S-δ R structure depends on the energy of the incident particles E rather weakly compared to the transmission through the dots. The left-to-right current in the left lead, which flows from Res LR to S and Res RL is given by [23] where T LR(RL) and V L(R) are, respectively, the temperature and electric potential of Res LR(RL) , R LL eh(ee) (E) is the probability of Andreev (normal) reflection on the N-δ L -S interface, T RL eh(ee) is the probability of the Cooper pair splitting (elastic co-tunneling) with regard to the particles incident from the right lead (here the upper subscript indicates the direction of the particle motion; e.g., RL means that the particle incident from the right lead is transmitted into the left one).
The probabilities T LL eh(ee) = |r eh(ee) | 2 and R RL eh(ee) = |t eh(ee) | 2 are determined by the transmission and reflection amplitudes given by [24]t eh = t L [t ee r R r eh + r eh r L t hh ] t R /D; (47) t ee = t L [t ee (1 − t 2 hh r L r R ) + r eh r L t hh r R r he ] t R /D; (48) r eh = t L r eh [1 + (t ee t hh − r eh r he ) r R r R ] t L /D; (49) r ee = r L + t L [r eh r L r he + t ee r R t ee − (t ee t hh − r eh r he ) 2 r L r R r R ] t L /D, where t L(R) and r L(R) are the transmission and reflection amplitudes of δ L(R) ; D is determined by multiple reflections inside the δ L -S-δ R structure: D = 1 − t 2 ee r L r R − t 2 hh r L r R − r eh r he (r L r L + r R r R ) + (t ee t hh − r eh r he ) 2 r L r R r L r R ; and t ee(hh) and r ee(hh) are the trasmission and reflection amplitudes of the superconducting part of the structure: t ee(hh) = e ±ipl S sin α sin (α − iql S ) ; (52) r eh(he) = sinh ql S i sin (α − iql S ) .
Here p 2 − q 2 = k 2 F and 2pq = (2m/ 2 )∆ sin α with α = arccos (E/∆); l S is the size of the superconducting region. Bearing in mind that I L (T, T, 0, 0) = 0 and that the temperature difference between Res LR and Res RL has a small contribution to the thermoelectric current, from Supplementary Equation (46) we obtain where the quantities defined as Using Kirchhoff's law for the currents at the intersection of Res LR , Res RL , and the superconductor, we obtain the following equations: where g L(R) is the conductance of Dot L(R) plus additional resistor (nanoribbon) R nr L(R) , and I δT L(R) is the locally induced thermoelectric current: We thus obtain V R = −I δT R + I δT L G NS N /(g L + G NS + G NS N ) g R + G NS + G NS N − (G NS N ) 2 /(g L + G NS + G NS N ) ; (65) The currents can be found by substituting V L and V R into Supplementary Equation (58). To select the non-local contribution from the I L(R) , one has to subtract I L(R) (averaging should be performed over different gate voltages on the right (left) dot) from I L(R) : As indicated by the calculations above, in the case of the incoherent system, the non-local electrical currents mostly originates from the electric potential difference between the intermediate reservoirs Res LR and Res RL . Therefore, in a strict sense, ∆I nl L(R) is not thermoelectric current but rather it is produced by the locally induced thermoelectric voltages V L and V R on Res LR and Res RL , which in this regard act as proxies. However, the plots for ∆I nl L(R) obtained from the present model are fundamentally different from the experimental data, which suggests that the observed non-local currents have a direct thermoelectric nature: Supplementary Fig. 12 displays the theoretically and experimentally obtained non-local current ∆I nl L on the left side as function of the gate voltages V sg,L and V sg,R in a simple situation where τ L (E) and τ R (E) have a Lorentzian form (i.e., the Fano resonance effects are absent). Both plots are characterized by similar patterns, but one can easily observe that they have different orientations of their symmetry axes. This indicates that while the experimental setup can still be subject to decoherence, which may enable the above discussed mechanism for the local thermoelectricity, the non-local current is mostly determined by the coherent electrical transport. Note that the exact values of ∆I nl L(R) in Supplementary Fig. 12b are not important since the purpose is rather to show the distinctive behaviour of the non-local electricity conditioned by the incoherent transport.