Theory of tunneling spectra for a few-electron bilayer graphene quantum dot

The tuneability and control of quantum nanostructures in two-dimensional materials offer promising perspectives for their use in future electronics. It is hence necessary to analyze quantum transport in such nanostructures. Material properties such as a complex dispersion, topology, and charge carriers with multiple degrees of freedom, are appealing for novel device functionalities but complicate their theoretical description. Here, we study quantum tunnelling transport across a few-electron bilayer graphene quantum dot. We demonstrate how to uniquely identify single- and two-electron dot states' orbital, spin, and valley composition from differential conductance in a finite magnetic field. Furthermore, we show that the transport features manifest splittings in the dot's spin and valley multiplets induced by interactions and magnetic field (the latter splittings being a consequence of bilayer graphene's Berry curvature). Our results elucidate spin- and valley-dependent tunnelling mechanisms and will help to utilize bilayer graphene quantum dots, e.g., as spin and valley qubits.


I. INTRODUCTION
Carbon-based materials are considered promising candidates for spin-based quantum computation devices due to their low spin-orbit and hyperfine coupling entailing long spin coherence life times [1][2][3] . Any spin-qubit operation using a quantum dot will necessarily include the steps of controlled loading (transferring a charge carrier onto the dot) and storage (keeping the charge carrier on the dot). Such an operation hence requires understanding and control of the dot's few-electron states and tunnel transport processes.
In bilayer graphene, recent experiments achieve confinement of charge carriers in one-and zero-dimensional structures by electrostatic gating [4][5][6][7][8][9][10] . To electrostatically define a nanostructure in bilayer graphene multiple gates locally modulate the bilayer graphene band gap and charge carrier density, cf. Fig. 1a). Split gates can define a channel (pink strip in Fig. 1a)), while finger gates on top create a dot-like region within this channel (dark pink region), bounded by gapped regions acting as barriers (white regions). This confinement method offers immense gate-control of the nanostructure, e.g., the confinement width, depth, barriers, and bilayer graphene gap. It is now possible to operate such an electrostatically confined bilayer graphene dot controllably in the single and few-electron regime [11][12][13][14][15][16][17][18][19] . The rapid experimental progress in device design, quality, and control, calls for a theoretical investigation of single and few-electron tunnelling processes in such structures.
The two internal degrees of freedom, valley and spin, enrich the spectra of bilayer graphene-based devices. The result is highly degenerate multiplets split in various ways by a magnetic field and weak perturbations. In this work, we investigate tunnelling transport through a bilayer graphene quantum dot in the single and fewelectron regime as a tool to unravel some of the dot's a) Bilayer graphene lead-dot-lead setup. b) Single (N = 1) to two-particle (N = 2) tunnelling transitions allow characterising the dot's orbital, spin and valley states. The ground state transition defines the Coulomb diamond in the bias (VB) and gate voltage (VG) plane. Transitions contribute to transport depending on spin and valley selection rules. c) The slope of differential conductance lines with the magnetic field (dominated by the difference of the twoparticle and single-particle valley g-factors) teaches about the orbital and valley composition of the two-particle dot state. The splittings between transitions at B = 0 manifest the interaction-induced two-particle state gaps.
two-particle states' unusual characteristics. We demonstrate how the specifics of the dot's multiplets manifest in tunnelling current, cf. Fig. 1c), and how to link such experimentally observable transport features with interaction and field-induced gaps between different spin and valley configurations. We determine the particular tun- We consider transport through bilayer graphene quantum dots in regimes dominated by different types of interactions. Weak or strong long-range Coulomb interactions favour two-particle states with symmetric or antisymmetric orbital states and distinct degeneracies of the spin and valley multiplets. Possible short-range interaction mechanisms include couplings generated by inter-valley scattering or "current-current" interactions, and spin-orbit coupling. In brackets, we indicate the corresponding section number.
nelling sequences for spin and valley states of differently ordered multiplets and relate them to microscopic parameters, such as short-range interaction coupling constants, g ⊥ , g zz , g 0z , g z0 20,30-32 , and topological valley gfactors (the latter induced by Berry-curvature [22][23][24][25] ). Besides spin-and valley selection rules, these tunnelling sequences depend on the dot-lead coupling characteristics, such as asymmetric coupling to the source and drain and cotunnelling corrections. By combining the aspects of state multiplicity, electronic interactions, and dot-lead couplings, our results add to the understanding of tunnelling transport in complex few-electron systems.
The paper is structured as follows. In Sec. II, we introduce our theoretical model of the bilayer graphene quantum dot and the leads, discussing the dot's state structure in the single-and two-particle sector. Section III describes the rates for tunnelling between these states and the leads, and the calculation of tunnel current using rate equations. Section IV presents our calculations of tunnel transport through a bilayer graphene quantum dot. We provide maps of the differential conductance, dI/dV , in the plane spanned by the gate voltage and the magnetic field for representative cases of interaction parameters. This way, we characterise regimes in which different electronic interactions dominate, as tabulated in Fig. 2. The differential conductance in a proper bias interval reveals the transitions between the one-and two-particle states in the quantum dot. Levels energies are closely related to the symmetries of the corresponding orbital wave functions. The multiplicity and ordering of the two-particle levels depends on the orbital symmetry, the short-range part of interactions, and the external magnetic field. The latter allows one to affect the level ordering. The interpretation of such tunnelling data may depend on device characteristics, e.g., the lead-dot coupling strength or uniformity of source and drain coupling. Taking these device features into account, we show how to use the differential conductance maps to identify the dot's two-particle ground state and determine the dominant microscopic interaction parameters. Section V contains step-by-step instructions how to use our results to interpret differential conductance data for identifying the single-and two-particle state structure of a bilayer graphene quantum dot.

II. MODEL
We consider a lead-dot-lead setup in which a bilayer graphene quantum dot is tunnel-coupled to bilayer graphene quantum point contacts as in Fig. 1a).
Single-particle states of the bilayer graphene quantum dot. We focus on the experimentally accessible regime of small and moderate displacement fields in the dot region. For a small gap, the bilayer graphene dispersion in the vicinity of the K-points is approximately quadratic, and a quantum dot's single-particle level structure resembles that of harmonic confinement, featuring an orbitally singly degenerate ground state 20 . These single-particle dot states are characterized by the orbital quantum number, n, and the electron's spin (s =↑ , ↓) and valley (t = +, −) degree of freedom. We denote a one-electron dot state by |n, is the electron creation operator and |0 is the empty dot state. The n-th spin and valley multiplet at zero magnetic field is characterized by energy, E n . Zero-point vibrations enhance Kane-Mele spin-orbit coupling 21 , ∆ SO , leading to reversed spin splitting in opposite two valleys. Each multiplet splits upon the application of a magnetic field, B, perpendicular to the bilayer graphene plane as, according to the free electron spin g-factor, g = 2, and valley g-factor, g n v (µ B being the Bohr magneton). The latter is a consequence of gapped bilayer graphene's nontrivial Bloch band Berry curvature entailing an topological orbital magnetic moment with opposite sign in the two different valleys [22][23][24][25] . As the orbital magnetic moment is a function of wave number in each valley, the topological valley g-factor depends on the gap and the states' momentum space distribution (and, consequently, on the orbital quantum number, n), determining how much orbital magnetic moment is picked up by the dot states 20,26-28 . The second term in Eq. (1) accounts for the presence of a gate with capacitance C G , which, at gate voltage V G , induces an effective charge on the dot, changing the dot's electrostatic potential by Here, N is the dot occupation number and C is the total capacitance of the dot. Two-particle states of the bilayer graphene quantum dot. The dot's two-particle sector is nontrivial, due to the large number of states arising from different combinations of the orbital, spin, and valley degrees of freedom. Moreover, these degrees of freedom are not independent since all three combined must form an overall antisymmetric two-particle wave function. As we showed in Ref. 20, Coulomb interaction between the two dot electrons further impacts the correlations between the different degrees of freedom.
The long-range Coulomb interaction on the scale of the dot state wave functions is given by, between the low-energy electronic fields, Ψ n (r), on the non-dimer sites A and B of the bilayer graphene lattice. We employ the 2D screened Coulomb interaction in a weakly gapped bilayer graphene 20,29 , with Fourier representation, V C (q) = e 2 4π 0 2π q(1+qR ) , where 0 is the vacuum permittivity, is the encapsulating substrate material's dielectric constant, R = √ 32 κ/ √ ∆, taking into account gapped bilayer graphene's polarisability 29 , κ 2 = 2me 2 /(4π 0 √ ∆) 2 , with m being the effective mass and ∆ the bilayer graphene gap. The Coulomb repulsion in Eq. (3) determines the spatial extent of the wave functions and the exchange energy. The competition of single-particle energies, direct-, and exchangeinteraction terms determines the mixing of single-particle orbitals forming orbitally symmetric or antisymmetric two-particle states 20 . For zero or weak Coulomb interaction (strong screening by the surrounding medium), two electrons on the dot form an orbitally symmetric wave function, both occupying the same single-particle orbital, n. If the Coulomb repulsion dominates (weak screening), the gain in exchange energy overcomes the cost of occupying higher single-particle orbitals, and the two-particle ground state forms an antisymmetric orbital wave function involving different single-particle orbitals, n and m.
In gapped bilayer graphene, where the gating needed to form the quantum dot lifts the layer symmetry, we take into account the lattice-scale symmetry breaking short-range interactions 20,30-32 , with ς AB i (ς +− i ) the Pauli matrices in sub-lattice (valley) space and (i, j) = (xx, xy, yx, yy, zz, z0, 0z). The interactions in Eq. (4) originate from symmetry breaking fluctuations and the relevant coupling constants g xx = g xy = g yx = g yy ≡ g ⊥ , g zz , g z0 , g 0z , (5) favour states with spontaneously broken symmetries 30,31 . Inter-valley scattering introduces the coupling g ⊥ . The couplings g 0z,z0 correspond to "current-current" interactions 33 , favouring states with spontaneously broken time-reversal invariance 31 . The case i = j = 0 is already included in Eq. The short-range interactions in Eq. (4) introduce anisotropies in the sublattice and valley space for twoparticle states with symmetric orbital wave function. For orbitally antisymmetric two-particle wave functions, contact interactions as in Eq. (4) are not relevant due to vanishing electronic density at small inter-particle distances. Short-range interaction induced splittings hence provide a way to distinguish orbitally symmetric and antisymmetric dot states.
Any theoretical estimation of the couplings' numerical values comes with inherent uncertainty since they depend on the relevant energy scale. The resulting renormalization and additional phonon-mediated effects can change the couplings g ij in absolute value and sign 31,32 . By studying tunnelling through two-particle multiplets for all possible combinations of values in Eq. (5) we demonstrate how to identify different parameters in transport. Our results will be relevant for unfolding experimental measurements using tunnelling spectroscopy of the bilayer graphene quantum dot's two-particle states as a tool to extract the microscopic short-range interaction parameters in Eq. (5).
Depending on the symmetry of the two-particle states' orbital part, any combination of spin and valley states is permissible that combines to an overall antisymmetric two-particle wave function. There are six combinations of spin/valley-singlet (σ −x /τ −x ) and -triplet (σ −z /τ −z , σ +x /τ +x , σ +z /τ +z ) states and an orbitally symmetric (s) two-particle state of orbital n: The energies of this two-particle multiplet are given by 20 , Here, E s nn , comprises the energy of the n-th singleparticle orbital and the screened electron-electron Coulomb interaction computed from Eq. (3). The factor, > 0 (for all combinations of t i corresponding to inter-and intra-valley scattering channels induced by Eq. (4)), captures specific dot state characteristics, i.e., dot shape, gap, and mode number. The short-range interaction constants, g ij , are a priori unknown and we discuss possible level orderings for different values of these couplings in Sec. IV. In a finite magnetic field, the two-particle levels split according to the g-factors in Eq. (7). The valley g-factor, 2g n v , of the two-particle states computes as the sum of the single-particle g-factors in the two valleys. For valley polarized states, 2g n v exceeds the single-particle valley and spin g-factors. Conversely, the g-factors from both valleys cancel for any valley coherent two-particle state.
The ten possible two-particle states with orbitally an-tisymmetric (a) wave function are, For brevity we consider the simplest case where the two-particle ground state consists of exactly two singleparticle orbitals n and m (substantial admixing of more than two orbitals is relevant only at higher energies 20 ).The energies of the states in Eq. (8) are, Here, E a nm is the energy of the orbitally antisymmetric states of two screened interacting electrons in singleparticle orbitals n and m (akin to the orbitally symmetric state described above), and (g n v + g m v ) is the valley g-factor of the two-particle multiplet.
Coupling to the leads. The point contacts in the bilayer graphene channel to the left and right of the quantum dot provide discrete lead modes due to the transverse confinement. These modes can couple to the quantum dot. Close to pinching off the lowest of their modes, we can treat the quantum point contacts as tunnel junctions with tunnelling amplitudes t L (t R ) for the left (right) quantum point contact. We describe these single-channel leads with a Hamiltonian, where c † lkst creates a lead electron with momentum k, energy l k , spin s, and valley quantum number t. The lead-dot tunnelling Hamiltonian is given by, In the following sections, we use this tunnelling Hamiltonian in Eq. (11) to compute the tunnelling current across the bilayer graphene quantum dot.

III. TUNNELING RATES AND RATE EQUATIONS
For our transport calculations, we consider the hightemperature regime where in the last term we compare to the energy difference of dot states with different particle number and Γ nst is the tunnel-coupling induced level broadening with lead density of states ν l . For a level broadening Γ nst much smaller than the thermal energy k B T , we can compute transport perturbatively in the tunnel Hamiltonian 34-36 , H T , in Eq. (11).
The lowest (first) order in the lead-dot tunnel coupling describes the single-electron processes involved in sequential tunnelling: an electron tunnels either from the leads to the dot or from the dot to the leads thereby changing the occupation number of the dot by one. Expanding to first order in H T and applying Fermi's golden rule, transition rates for a one-electron tunnelling, which induces a transition of dot from a single-particle state, |N = 1 : χ , to a two-particle dot state, |N = 2 : χ , read, Here, N indicates the dot particle number and χ identifies the state of the corresponding multiplet. Hence, a prefix N = 1 implies χ = (n, s, t) and for N = 2, χ indexes the orbital, spin, and valley combinations from the family of states in Eqs. (6) or (8), respectively. In Eq. (13), (nst) are the indices of the electron tunnelling into the dot, forming the two-particle state |2 : χ with the single electron previously on the dot (the latter having quantum numbers χ ). The initial and final states of the leads are |i = |i L |i R and |f = |f L |f R , the former weighted by a thermal distribution ρ i . Further, f denotes the Fermi function and µ l is the chemical potential of lead l, which depends on the bias voltage, V B . We consider the case where the dot is biased symmet-rically, µ L/R = ±eV B , with respect to the equilibrium chemical potential. The rates for the reverse transitions, We provide the explicit rates for each transition in Appendix A.
Going to second order in H T describes correlated twoelectron cotunnelling: an electron tunnels from one lead to the other (or the same lead) via the quantum dot, leaving the occupation number of the dot invariant. Within each particle number sector (N = 1 or N = 2), the dot's state may change (inelastic cotunnelling) or remain the same (elastic cotunnelling). The corresponding cotunnelling rates read, These rates involve the intermediate states of higher or lower dot occupation number, N ± 1, if they are allowed by spin and valley selection rules. In Eq. (14), we take into account transitions via the single-particle ground state multiplet and the two-particle ground state multiplets of Eq. (6) and (8). Projection onto these singleparticle and two-particle state spaces is valid for quantum dots where all other states are separated sufficiently in energy to exclude any virtual transitions to them. It is not straightforward to evaluate the cotunneling rates in Eq. (14) due to the second-order poles causing the integrals to diverge. These divergences are related to the intermediate state's zero width and hence infinite lifetime within this perturbative approach. We follow the standard regularization procedure to extract the correct cotunneling rates from Eq. (14) 36-40 : First, a level width γ ∼ Γ nst is introduced as imaginary parts in the denominators (accounting for the intermediate states' tunnelcoupling induced level broadening). These imaginary parts shift the poles away from the real axis, and the integrals can be carried out. Next, the resulting expression is expanded in powers of γ. The leading order term is a sequential-tunnelling contribution (reflecting that, at finite temperature, the final state of any cotunnelinginduced transition can also be reached via two successive single-electron tunnelings). This term is disregarded to avoid double-counting sequential tunnelling processes. The next-to-leading-order term in the γ expansion gives the regularized expression for the cotunneling rate, where the limit γ → 0 can be taken. We provide the regularization calculations and resulting expressions for the cotunnelling rates W N:χ←N:χ in appendix B. Given the rates for transitions between different dot states, we write a master equation describing the dynamics of the probabilities, P N:χ , for the state, |N : χ , to be occupied at a given time, where the terms with changing particle number, 1 2, describe current flow whereas cotunnelling terms introduce relaxation within the multiplets at fixed particle number. We solve these rate equations, Eq. (15), in the stationary limit,Ṗ N:χ = 0, using the normalization condition N:χ P N:χ = 1. From the probabilities we compute the total particle current I = I seq + I cot , with the sequential tunnel currents flowing from the dot to lead l, It depends on the tunnelling strength compared to the isolated dot's level splitting whether second-order cotunnelling processes contribute significantly to transport. We define the regime of purely sequential tunnelling for weak dot-lead tunnel coupling, and the regime of sequential + cotunnelling for stronger dot-lead tunnel coupling, where second order effects contribute. Numerically, we find that the first regime is realized for |t l nst | ∼ ∆E N ,N ±1 /1000 while reaching the latter regime requires approximately |t l nst | ∼ ∆E N ,N ±1 /100.

IV. RESOLVING THE TWO-PARTICLE DOT STATES
A. Spectroscopy of an orbitally symmetric two-particle ground state This section considers dots with orbitally symmetric two-particle ground state wave functions. We discuss the possible level orderings which can result from Eq. (7) and at zero and finite magnetic field and how to distinguish the spin and valley states in tunnelling transport.
Possible level orderings of orbitally symmetric two-particle dot states. We illustrate the various level orderings of the states in Eq. (6) for different signs and relative magnitudes of the short-range couplings g 0z , g z0 , and g ⊥ in Fig. 3. Generally, there are three levels at zero magnetic field, being singly, doubly, and three-fold The ordering of the dot's two-particle lowest-state multiplet with symmetric orbital wave function at zero bias voltage, Eq. (7), depend on the relative magnitude and sign of the short-range interaction coupling constants, g ⊥ , g0z, and gz0, and on the magnetic field.
degenerate, respectively. These degeneracies are lifted by a finite magnetic field, splitting different valley and spin states. According to Eq. (7), the coupling constant g zz shifts all energies equally. The mutual splitting between the two inter-valley coherent states, τ ±x , is proportional to the coupling g ⊥ , while these states are split from the valley polarized states, τ ±z , proportionally to the sum g 0z + g z0 .
Tunnelling transitions in the single-and two-particle sector allow identifying the spin and valley states and determining the short-range couplings by combining the two following considerations: Firstly, the single-and twoparticle states split in a magnetic field. Transition energies hence depend on the difference in single-and twoparticle valley g-factors. Besides, any single-particle-totwo-particle tunnelling transition is subject to spin and valley selection rules. Therefore, we can identify the twoparticle states that can be reached, e.g., from the singleparticle ground state. With the two-particle levels being identified, we can relate the level splittings to the shortrange interaction couplings g 0z , g z0 , and g ⊥ as in Fig. 3. Hence, classifying the dot's two-particle states and their mutual gaps is a way to quantify bilayer graphene's microscopic short-range interaction parameters.

Two-particle states with broken time-inversion symmetry
Single-to-two-particle transitions to the levels in Fig. 3 yield differential conductance features as in Fig. 4. Here, we consider sequential tunnelling and symmetric coupling to the leads. Differential conductance maps as the ones in Fig. 4 are cuts at finite bias voltage (we chose V B = 0.45 meV) through the Coulomb diamonds for different values of magnetic field (cf. Fig. 1). Each allowed single-totwo-particle transition manifests as an increase/decrease in conductance (red/blue lines) once this transition enters/leaves the bias window. The differential conductance features at zero magnetic field reflect the splittings of the two-particle multiplets in Fig. 3. At finite magnetic field, the conductance lines disperse according to the two-particle and single-particle g-factors. Hence, for similar zero-field splittings and similar g-factors, the conductance maps can coincide even for distinct two-particle level orderings.
To facilitate the electron transport through a dot, the bias window must allow single-particle-to-two-particle transitions between the ground state of the dot with one and two electrons, respectively. When the two-particle ground state is valley coherent, the corresponding lines in the differential conductance maps have positive slope in a magnetic field. While these coherent two-particle states do not disperse with B, the single-particle ground . Differential conductance maps for the two-particle multiplets in Fig. 3 for the different possible regimes of shortrange interaction constants g ⊥ , g0z, and gz0. (sequential tunnelling and symmetric coupling to both leads). The conductance increases/decreases when a single-particle-to-two-particle transition enters/leaves the bias window at fixed bias voltage VB = 0.45 mV (red/blue lines). Permissible single-to-two-particle transitions depend on spin and valley selection rules and whether the two-particle ground state is spin and valley coherent or polarized. The difference of the two-particle and single-particle valley g-factors dominates the slope of the lines with magnetic field. Here, kBT = 0.003 meV.
state, |n, ↓, − is pushed down and the energy required for this transition increases. Conversely, a K − valley polarized two-particle ground states is pushed down even faster with B (since 2g n v > g n v ), causing the transitions energy to decrease. This leads to lines with negative slopes limiting the bias window range in Fig. 4 for these cases. Within the bias window, whether energetically allowed single-particle-to-two-particle transitions contribute to transport is determined by spin and valley selection rules. For example, the K + excited singleparticle states can be populated via transitions to valley coherent two-particle states. However, if there are no such transitions available at equal or lower gate voltage, the K + single-particle states are depopulated, causing the corresponding lines to terminate in the differential conductance maps in Fig. 4.

Two-particle states preserving time-inversion invariance
In the following sections IV A 2, IV B and IV C, we exemplify the quantum dot's tunnelling characteristics for one specific level arrangement of the orbitally symmetric two-particle states and study different regimes of lead couplings as well as the impact of a finite spin-orbit coupling gap. Numerical values we have estimated previously in one specific dot model 20 , yielded g 0z = g z0 = 0, (preserving time-reversal invariance), g zz > 0, g ⊥ < 0 (favouring the spin and valley coherent ground state |nn, σ −x , τ −z ), and Jg zz 4|Jg ⊥ |. For this choice of short-range couplings, the two-particle triplet is equally spaced at B = 0 (top left panel of Fig. 5). A finite magnetic field splits the levels according to the spin and valley configuration (top row of Fig. 5).
The contrasting magnetic field coupling of valley polarized and valley coherent two-particle states leads to level crossings at finite B. For zero and small magnetic field the state |nn, σ −x , τ +x is the two-particle ground state. This spin-and valley-coherent state does not couple to the magnetic field. At sufficiently large B, the valley polarized state, |nn, σ −x , τ −z , being pushed down by the magnetic field, becomes the two-particle ground state. Being able to identify the differential conductance characteristics in Fig. 5 with the possible single-particleto-two-particle transitions allows extracting information about a symmetric two-particle dot state.
In the regime of sequential tunnelling, transport is possible, once the gate voltage sufficed to induce the ground state-to-ground state transition. For zero or weak magnetic field, this is the transition, The involved two-particle state occupies all four different spin and valley states. Hence, when one electron leaves the dot in the subsequent tunnelling process, the remaining electron can be in any of the single-particle 5. Transport characteristics of a bilayer graphene quantum dot with an orbitally symmetric two-particle ground state preserving time-inversion invariance (in which case, gz0 = g0z = 0, compared to Fig. 4). Top: Single-particle and two-particle dot levels for different magnetic fields. Bottom left: Differential conductance across the dot in different tunnelling regimes and symmetric or asymmetric lead coupling (suppressing coupling to the right lead = source or the left lead = drain, respectively) depending on a magnetic field, B. For the conductance maps, we fix VB = 0.45 mV and kBT = 0.003 meV. The sequential tunnelling regime is realized for |t l nst | ∼ ∆EN ,N ±1/1000 while significant cotunnelling contributions require approximately |t l nst | ∼ ∆EN ,N ±1/100. Cotunnelling induces relaxation processes within each fixed particle number multiplet and hence opens additional transport channels compared to purely sequential tunnelling. The panels on the right consider the potential influence of a finite spin-orbit coupling gap, ∆SO, possibly smaller (here |∆SO| = 0.02 meV) or larger (here |∆SO| = 0.1 meV) than the splitting induced by g ⊥ and of different sign. The labels A ○, B ○, C ○, D ○ indicate transitions which we discuss in detail in the main text.
states. Consequently, with increasing gate voltage, all single-particle-to-two-particle transitions become possible and manifest in the differential conductance maps within the bias window. At higher magnetic fields, the ground state-to-ground state transition changes to, where the valley K − polarized two-particle state entails that after the next tunnelling process, the remaining elec-tron occupies one of the K − single-particle states. As a consequence, transitions from the K + single-particle states do not contribute to conductance in this regime if there is no transition to a valley-coherent two-particle state lower in gate voltage. The corresponding lines C ○ terminate in the differential conductance maps in Fig. 5. We note that coupling stronger to the source (left lead in our convention) and suppressing the drain coupling (right lead) suppresses transport features from P 1: P 1: B=0.05T B=0.25T Regime of sequential tunneling Single-particle states Two-particle states Probabilities for the single-particle and orbitally symmetric two-particle states to be occupied for fixed magnetic field cuts along the gate voltage axis through the symmetrically coupled differential conductance maps with zero spin-orbit coupling in Fig. 5 (leftmost differential conductance maps) in the sequential tunnelling regime. the transitions involving valley polarized two-particle states. In comparison, stronger coupling to the drain decreases the amplitudes of all transport channels.A finite spin-orbit coupling gap, ∆ SO , further splits the states and corresponding transitions depending on its sign and magnitude relative to the short-range splittings. When the spin-orbit gap overcomes the short-range couplings, |∆ SO | > 4|g ⊥ J|, transitions may occur in a different order, exemplified by the transition D ○ in Fig. 5. We depict representative differential conductance maps in the regime of sequential tunnelling and symmetric leadcoupling in the bottom right panels of Fig. 5.
Cotunnelling leads to relaxation processes within the multiplets of each seperate particle number sector and can hence make additional transport channels available. In the coupling regime where cotunnelling processes play a significant role, the transitions from the K + singleparticle states to the τ −x two-particle states ( C ○ in Fig. 5) reappear compared to the regime of purely sequential tunnelling as a result of population of these singleparticle states via the |n, ↓, − → |n, ↑, + , |n, ↓, + cotunnelling transitions (cotunnelling assisted sequential tunnelling 34 ). Besides, we observe features outside the Coulomb diamond, where cotunnelling events populate states that do not yet fall into the bias window for a certain gate voltage value. Increasing magnetic field and any asymmetry in the lead couplings suppress cotunnelinginduced transport features. The former is due to en- ergy differences between states growing with B, suppressing inelastic events. The latter suppression comes from the fact that at finite bias, the relevant contributions to cotunnelling scattering rates involve tunnelling at both leads (cf. Eq. (14)). The occupation probabilities, P N:χ shown in Fig. 6 for different values of magnetic field support the conclusions above about states contributing to transport in different regimes. Allowed transitions manifest as steps where state occupation numbers change. Cotunnelling processes alter these steps by introducing alternative transitions between dot states (see appendix C).
Each transition in Eq. (20) contributes a line to the differential conductance maps in Fig. 8, the slopes of which are given by the g-factor difference of the involved single-particle and two-particle states. The first transition listed in Eq. (20) is the ground state-to ground state transition. The single-particle excited state |n, ↑, − is populated from the spin coherent two-particle state |nm, σ +x , τ −z via the tunnelling sequence |n, ↓, − → |nm, σ +x , τ −z → |n, ↑, − . The transitions F ○ : |n, ↓, + , |n, ↑, + → |nm, σ −z , τ +x , |nm, σ +x , τ +x , |nm, σ +z , τ +x , are absent in the sequential tunnelling differential conductance maps as the K + single-particle states are not populated at the values of gate voltage needed for these transitions. Electrons cannot reach the K + singleparticle states because all transitions lower in gate volt- Regime of sequential tunneling Single-particle states Two-particle states Occupation probabilities for the single-particle and orbitally antisymmetric two-particle states for cuts at B = 0.5 T along the gate voltage axis through the symmetrically coupled differential conductance maps in Fig. 8 in the sequential tunnelling regime.
age, including the ground state-to-ground state transition |n, ↓, − → |nm, σ −z , τ −z , occur between valley K − polarized states. Cotunnelling transitions, when relevant, enable the transitions in Eq. (21), by populating the K + single-particle states via inelastic cotunnelling |n, ↓, − → |n, ↑, + , |n, ↓, + . This cotunnelling-induced repopulation makes sequential tunnelling from the K + single-particle states possible leading to weak features in the differential conductance maps at the gate voltages required for the transitions in Eq. (21) (bottom row of Fig. 8). Additionally, we observe cotunnellinginduced transport features outside the Coulomb diamond similar to the case of the orbitally symmetric multiplet, Sec. IV A. Similarly, all cotunnelling features are suppressed by magnetic field and asymmetric coupling to the leads. Figure 9 demonstrates the cotunnelling-mediated redistribution of electrons among the states by comparing the occupation probabilities in the purely sequential tunnelling and sequential tunnelling + cotunnelling regimes. We note that a finite spin-orbit coupling gap ∆ SO leads to two split copies of fanning lines in the differential conductance maps as those in Fig. 8.
C. Interplay of ground-and excited two-particle state multiplets The dot's two-particle ground and first excited state can be sufficiently close in energy for both to contribute transport signatures within the bias window 20 . We consider the cases in which the two-particle ground state is either orbitally symmetric or antisymmetric, while the first excited state's orbital wave function is of the opposite symmetry. These scenarios yield distinct cases compared to the isolated two-particle ground states discussed in the previous sections. The ground state-to-ground state transitions originating from the two-particle states of opposite symmetry can enable different transitions in  10. Differential conductance maps for tunnelling transport through a bilayer graphene quantum dot in which the twoparticle ground-and first excited state are close enough in energy for both to be reached within the bias window are not merely superpositions of the two maps in Figs. 5 and 8 for the two states due to interplay of the different multiplet states' occupation numbers. a) Orbitally symmetric two-particle ground state and antisymmetric excited state, b) Orbitally antisymmetric twoparticle ground state and symmetric excited state. Parameters as in Fig. 5.
the excited state multiplet compared to the isolated case. Also, we can clearly distinguish the orbitally symmetric and antisymmetric two-particle states by their zero-field splittings or absence thereof. Hence, investigating both simultaneously reveals changes in the orbital composition when comparing the ground and excited two-particle states. Figure 10a) depicts the differential conductance across a dot with an orbitally symmetric two-particle ground state and orbitally antisymmetric first excited state. Transitions to both two-particle multiplets manifest in the differential conductance maps. Notably, tunnelling channels involving valley-coherent two-particle states in the orbitally symmetric two-particle manifold lead to a population of the K + valley at sufficiently high gate voltages. These populations enable transitions to all the orbitally antisymmetric two-particle states by purely sequential tunnelling. The differential conductance lines originating from either multiplet in Fig. 10a) have distinct slopes with B due to the different orbital composition of the symmetric and antisymmetric orbital twoparticle wave functions yielding different valley g-factors.
Since the orbital composition is unequal also for ground and excited states, the valley g-factors differ for the orbitally antisymmetric states in Fig. 10a) and Fig. 8.
Similar statements apply to the case of an orbitally antisymmetric two-particle ground state and orbitally symmetric first excited state, Fig. 10b). Also here, an orbitally symmetric state occupies different orbitals, n, when being an excited state compared to a ground state, leading to different valley g-factors and different magnetic field splittings compared to Fig. 5.

V. DISCUSSION AND CONCLUSION
In summary, we have analysed quantum tunnelling across an electrostatically induced bilayer graphene quantum dot as a spectroscopic tool to resolve the dot's single and highly degenerate two-electron multiplets. Here, we summarise how to use tunnelling transport maps as a function of gate voltage and magnetic field to distinguish the interaction regimes specified in Fig. 2 and identify two-particle states with different orbital, spin, and valley compositions: • The number and the splittings of peaks in the differential conductance at zero magnetic field tell about the orbital symmetry of the two-particle wave function. An orbitally antisymmetric two-particle state (as for dots with weak screening and strong long-range Coulomb interaction, cf. Fig. 2) hosts a tenfold degenerate multiplet of spin and valley states at B = 0 (cf. Fig. 7), manifesting in one single transition. Conversely, the six possible spin and valley states of an orbitally symmetric two-particle state (which forms for strongly screened long-range interactions) are slightly split by short-range lattice-scale interactions (cf. Fig. 3). Such splittings manifest in multiple possible transitions and corresponding tun-nelling transport features within the bias window at zero magnetic field (Figs. 4, 5, and 8).
• The various spin and valley states couple differently to a perpendicular magnetic field. Hence a magnetic field allows us to identify them and infer their g-factors. Spin-and valley-polarized states split with B, while spin-and valley-coherent two-particle states do not couple to a magnetic field. In combination with spin and valley selection rules, this contrasting magnetic field dependence helps identify the dispersing lines in the magnetic field-dependent differential conductance maps with the corresponding single-particleto-two-particle transitions. The slope of these lines is proportional to the difference of the singleparticle and two-particle states' g-factors. The orbital magnetic moment induced valley g-factor being much larger than the free particle spin g-factor allows distinguishing spin and valley splittings.
• If multiple two-particle states can be reached within the bias window, their distinct valley g-factors will help identify them. The valley gfactor depends on the orbital wave function and its distribution in momentum space. Hence transitions from the same single-particle state to distinct twoparticle states show as lines with different slopes in the differential conductance maps, as in Fig. 10.
Our results will help to explain tunnelling transport experiments in bilayer graphene quantum dots in the oneand two-particle sectors. Identifying and controlling fewelectron states is a crucial step towards using their degrees of freedom for quantum information storage and processing in future devices.

VI. ACKNOWLEDGEMENTS
in terms of the Bose function n B and the polygamma function ψ.
Appendix C: Occupation probabilities Similar to Figs. 6 and 9 in the main text, here we discuss how the occupation probabilities of the singleparticle and the orbitally symmetric or the orbitally antisymmetric two-particle states change when varying the gate voltage at fixed values of the bias voltage and the magnetic field. Figures 11 and 12 compare the occupation probabilities of the dot states for sequential tunnelling to the regime of sequential + cotunnelling for the two different two-particle multiplets. We observe how the additional inter-multiplet transitions induced by cotunnelling processes alter the states' occupations and enable different tunnelling sequences compared to purely sequential tunnelling.  11. Probabilities for the single-particle and orbitally symmetric two-particle states to be occupied for fixed magnetic field cuts along the gate voltage axis through the symmetrically coupled differential conductance maps with zero spin-orbit coupling in Fig. 5 (leftmost differential conductance maps).

Regime of sequential tunneling
Regime of sequential + cotunneling Single-particle states Two-particle states FIG. 12. Probabilities for the single-particle and orbitally antisymmetric two-particle states to be occupied cuts at B = 0.5 T along the gate voltage axis through the symmetrically coupled differential conductance maps in Fig. 8.