Electrostatically Confined Monolayer Graphene Quantum Dots with Orbital and Valley Splittings.

The electrostatic confinement of massless charge carriers is hampered by Klein tunneling. Circumventing this problem in graphene mainly relies on carving out nanostructures or applying electric displacement fields to open a band gap in bilayer graphene. So far, these approaches suffer from edge disorder or insufficiently controlled localization of electrons. Here we realize an alternative strategy in monolayer graphene, by combining a homogeneous magnetic field and electrostatic confinement. Using the tip of a scanning tunneling microscope, we induce a confining potential in the Landau gaps of bulk graphene without the need for physical edges. Gating the localized states toward the Fermi energy leads to regular charging sequences with more than 40 Coulomb peaks exhibiting typical addition energies of 7-20 meV. Orbital splittings of 4-10 meV and a valley splitting of about 3 meV for the first orbital state can be deduced. These experimental observations are quantitatively reproduced by tight binding calculations, which include the interactions of the graphene with the aligned hexagonal boron nitride substrate. The demonstrated confinement approach appears suitable to create quantum dots with well-defined wave function properties beyond the reach of traditional techniques.

The charge carriers in graphene at low energies, described as massless Dirac quasiparticles,[1] are expected to feature long spin coherence times.[2−5] Exploiting this property requires precise manipulation of individual Dirac electrons. Quantum dots (QDs) present an essential building block, yet providing tailored confinement in graphene has remained challenging. So far, e-beam lithography[6] and various other techniques[7−12] have been used to design nanometer-sized devices. However, their performance lacks behind, for example, GaAs QDs,[13,14] as disordered sample edges of patterned graphene result in uncontrolled charge localization and scattering.[6,15−17] So far, no clear evidence for 4-fold degenerate charging sequences has been reported in transport measurements of tunable QDs. Moreover, failing to controllably lift graphene’s valley degeneracy renders spin qubits unfeasible.[2,18,19]

In principle, bilayer graphene could improve the situation because an electric displacement field opens a band gap at regular AB stacking.[20] Indeed, electrostatically confined QDs in bilayer graphene exhibit Coulomb blockade,[21−23] yet controlling the spin or valley degree of freedom of an individual state has also not been demonstrated. Moreover, confinement is still prone to parasitic conduction channels due to residual disorder in the band gap or conducting channels along domain walls of AB- and BA-stacked areas.[24] Another approach exploits whispering gallery modes in electrostatically confined QDs[25−27] but here the control of the wave functions by gates is difficult and dwell times are extremely short (<100 fs). On an even more intricate route, the tip of a scanning tunneling microscope (STM) is used to locally stretch a suspended monolayer graphene sheet.[28] The onset of charge quantization due to induced strain showcases confinement by pseudomagnetic fields. Adding a real magnetic field B leads to charging sequences with regular orbital but no valley splittings.[28] Creating multiple QDs in this fashion would require independent strain control for every QD on the suspended graphene. Thus, such an approach is barely scalable.

Landau quantization helps to overcome Klein tunneling by opening band gaps.[21−23] An elegant method to exploit this by combining a magnetic field and an electrostatic potential has been proposed theoretically.[29−31] Indeed, indications of such confinement have been found in metal contact-induced pnp junctions,[32] graphene on SiO2,[33,34] and a suspended graphene nanoribbon.[35] However, in these experiments the confinement potential was not tunable but was generated by electrostatic disorder.

Here, we demonstrate controlled confinement by a combination of magnetic and electrostatic fields. We use the tip-induced electrostatic potential of an STM[36,37] in a B field perpendicular to the graphene plane (Figure a). Scanning tunneling spectroscopy (STS) reveals sequences of charging peaks by means of Coulomb staircases that appear when these confined states cross the Fermi energy EF. The peaks systematically group in quadruplets for electrons and holes corresponding to the 4-fold (valley and spin) degeneracy in graphene (Figure c,d). Moreover, some quadruplets separate into doublets due to an additional valley splitting induced by the hexagonal boron nitride (BN) substrate. STS as a function of B reveals that the first confined states emerge from Landau levels (LLs) with indices ±1. A third-nearest neighbor tight binding (TB) calculation[38,39] reproduces the onset of charging events as a function of tip voltage Vtip and B as well as the magnitude of orbital and valley splittings.

Figure 1

(a) Sketch of the experiment. Graphene covers a 30 nm thick hexagonal boron nitride flake on graphite. The magenta line represents the tip-induced confinement potential of graphene Φgrel for electrons, calculated as the numerical solution of Poisson’s equation (Supporting Information). (b) Energy diagram in real space: Fermi energy EF, black dashed line; local band bending Egr, magenta line; states belonging to electron (hole) LLs, blue (red); bulk LLs, 1, 0, −1. States embedded in the LL0–LL+1 gap (thin blue lines) are electrostatically confined. (c) Energy level diagram for the first two orbital states of a graphene QD exhibiting an orbital splitting Δ1o. Both orbitals are 4-fold degenerate, as indicated by black arrows representing physical spin. (d,e) Charging peak sequence in the differential conductance dI/dV corresponding to the level diagrams in panels c and f, respectively. Charging peaks are separated by the addition energy Eadd = EC + Δ where EC ≈ EC is the charging energy and Δ is comprised of Δo and/or the valley splittings Δτ. In panel d, quadruplet ordering showcases a dominant Δ1o, while Δτ become sizable in panel e, further separating quadruplets into doublets. (f) Same as panel c but including additional Δτ. The spin splitting Δσ is neglected, as Δσ < Δτ, Δo, EC in experiment.

We now sketch the principle of our experiment. A homogeneous, perpendicular B field condenses the electronic states of graphene into LLs at energieswhere νF is the Fermi velocity and is the LL index.[1] Consequently, energy gaps between the LLs emerge in the electronic spectrum. The smooth electrostatic potential Φgrel (magenta line in Figure a) induced by the STM tip locally shifts the eigenenergies ε(Φgrel) of charge carriers relative to the bulk LL energy (eq ). Shifting ε into the Landau gaps creates confined states (Figure b).[30] The shape of Φgrel determines the single-particle orbitals and energy levels, as in the case of artificial atoms.[14] Orbital splittings Δo separate the energy levels (Figure c), which we deduce experimentally to be Δo = 4–10 meV (see below) and thus Δo is small compared to the first LL gap E1 – E0 ≈ 100 meV at 7 T. While pristine graphene exhibits a 4-fold degeneracy, varying stacking orders of graphene on top of BN induce an additional valley splitting Δτ, which turns out to be smaller than Δo in our experiment. The finite B field creates a small Zeeman splitting estimated as Δσ = gμBB ≈ 800 μeV at 7 T (g-factor of 2, μB: Bohr’s magneton). Accordingly, the orbital splittings separate quadruplets of near-degenerate QD states, which exhibit a subtle spin-valley substructure (Figure f).

We use the STM tip not only as source of the electrostatic potential and thus as gate for the QD states but also to sequence the energy level spectrum of the QD as the states cross EF, that is, as the charge on the QD changes by ±e. This leads to a step in the tunneling current I(Vtip) and a corresponding charging peak in the differential conductance dI/dVtip. In addition to the single particle energy spacings, every additional electron on the dot needs to overcome the electrostatic repulsion to the electrons already inside the QD,[40] given by the charging energy EC. Thus, we probe the total energetic separation of charge states i and i + 1, given by the addition energy Eadd = EC + Δ, where Δ consists of Δo, Δτ, and/or Δσ. As we experimentally find EC ≈ EC ≈ 10 meV ≳ Δo (nearly independent of the charge state i, see below), the quadruplet near-degeneracy of the QD states translates to quadruplet ordering of the charging peaks (Figure d). Whenever either Δτ or Δσ significantly exceeds the other and temperature, quadruplets separate into doublets (Figure e).

We prepare our sample (see Figure a and Supporting Information) by dry-transferring[41,42] a graphene flake onto BN.[43−45] During this step, we align both crystal lattices with a precision better than 1° (Supporting Information). Then we place this graphene/BN stack on a large graphite flake to avoid insulating areas and simplify navigating the STM tip. Any disorder potential present in the sample will limit the confinement as long as it is larger than the Landau level gaps, thus larger gaps (e.g., the LL0–LL±1 gap) result in improved confinement. Moreover, the induced band bending will only be well-defined if the disorder potential is smaller than the maximum of Φgrel. By using the dry-transfer technique[41,42] and a graphite/BN substrate we reduce disorder in the graphene significantly.[46−48]

Probing the sample in our custom-build UHV-STM system[49] at T = 8 K, we observe the superstructure with a = 13.8 nm periodicity, which develops due to the small lattice mismatch of 1.8% between graphene and BN.[47] An atomically resolved STM image of this superstructure is presented in Figure d. Prior to measuring dI/dV spectra, the tip–sample distance is adjusted at the stabilization voltage Vstab and current Istab and then the feedback loop is turned off (Supporting Information). Figure a shows exemplary dI/dV spectra, acquired at B = 7 T and adjusted to the same vertical scale by dividing dI/dV by the first value I0 of the respective I(V) curve (Supporting Information). We observe pronounced, regularly spaced peaks for Vtip < −170 mV and Vtip > 500 mV. A closer look at the sequences reveals the expected grouping in quadruplets, which can still be distinguished up to the 20th peak. This grouping becomes even more evident by directly comparing the voltage difference between adjacent peaks ΔV in Figure b,c; ΔV between quadruplets is up to twice as large as ΔV within the quadruplets indicating Δo ≲ EC while Δτ and Δσ are significantly smaller. To further elucidate grouping patterns, we measure 6400 dI/dV spectra at equidistant positions within a 60 nm × 60 nm area, thus probing all areas of the superstructure. The median ΔV values (orange circles in Figure b,c) portray the robust ordering into quadruplets on the hole side, implying Δo generally dominates over Δτ and Δσ. On the electron side of the spectra, the sequences are disturbed by a few additional charging peaks of defect states in the BN substrate[50] that are identified by their characteristic spatial development (Supporting Information). This limits the comparability of the electron and hole sector and hides possible smaller electron–hole asymmetries in the data. The dI/dV features in between the charging peaks most likely capture contributions from multiple orbital states of each LL, which are lifted in degeneracy by the tip-induced potential, but cannot be identified unambiguously (Supporting Information, Section 5).

Figure 2

(a) Representative differential conductance spectra dI/dV(Vtip), normalized by the first value I0 of the respective I(Vtip) curve (Supporting Information). Recording positions are X1, between AA and AB; X2, on AB; X3, between AB and BA (compare panel d). Spectra on other regions (e.g., AA, BA) look similar. Vstab = 1 V, Istab = 700 pA, Vmod = 4.2 mVrms and B = 7 T. Quadruplets of peaks are marked by “4” and the first charging peak on either Vtip side by an asterisk. Curves are offset for clarity, while horizontal gray lines mark dI/dV = 0 S. Inset shows a zoom with Gaussian fits (dashed lines) used to extract distances between adjacent peaks ΔV as marked. (b,c) ΔV as a function of consecutive peak index for spectrum X1 (blue, error bars smaller than symbol size) and the median values for 80 × 80 spectra recorded on 60 × 60 nm 2 (orange). (d) Atomically resolved STM image (raw data) of the aligned graphene on hexagonal boron nitride (BN). Vtip = 400 mV, I = 1 nA. Differently stacked areas AB, BA, and AA marked and sketched by ball models. Inset on the upper left shows a zoom into the AB stacked area, marked by the blue square, exhibiting an obvious sublattice symmetry breaking due to the underlying BN. Positions equivalent to those where spectra in panel a were recorded are marked by circles labeled X1, X2, X3.

To understand the origin of the charging peaks, we provide a detailed microscopic picture of the tip-induced gating of localized states. We will only discuss the case of positive Vtip, that is, electron confinement, because the arguments for negative Vtip are analogous. Increasing Vtip (orange arrow in Figure ) shifts the states underneath the tip energetically down. States originating from LLs with positive index are embedded in the LL0–LL+1 gap that provides electrostatic confinement (Figure a, see also Figure b). Within the bias window eVtip = μgr – μtip, electrons tunnel from the sample into unoccupied states of the tip. One current path (dashed green arrow Figure a) passes through states of the QD (blue lines). The other stronger current path (solid green arrow Figure a) originates from the quasi-continuous LDOS at lower energies where energetically overlapping LL states strongly couple to the graphene bulk. Though increasing Vtip gates QD states down (Figure b), the Coulomb gap around EF always separates the highest occupied from the lowest unoccupied state, prohibiting continuous charging of confined states. It is only when the next unoccupied level crosses μgr that the QD is charged by an additional electron. The electrostatic repulsion due to its charge abruptly increases the Hartree energy of all states, thereby shifting additional graphene states from below μtip into the bias window (Figure b, central transition). Consequently, the tunneling current I increases which translates to a charging peak in dI/dVtip (Figure c). This mechanism is called Coulomb staircase[40] and has been observed previously, for instance, for charging of clusters within an STM experiment.[51] In essence, charging peaks in dI/dV signal the coincidence of a charge level of the QD with μgr[52] and thus provide a clear signature of the addition energy spectrum of the QD.

Figure 3

Sketch of the Coulomb staircase. (a) The chemical potentials of graphene μgr (black dashed line) and tip μtip (black solid line) define the bias window eVtip within which graphene states tunnel into empty tip states. There are two current paths available: (i) a weak one (green dashed arrow) via quantum dot states (blue lines), (ii) a dominant one (solid green arrow) via states strongly coupled to the graphene bulk (marked LDOS). Left: bulk graphene LLs away from the tip-induced band bending. (b) Schematic diagram of change in QD energies (blue lines) and quasi-continuous LDOS underneath the tip (green and gray triangle) for increasing Vtip from left to right. Between the second and third frame, the QD changes its charge state shifting the energy of the QD states and the entire LDOS upward. (c) Tunneling current I displaying the staircase (green line) and differential conductance dI/dV (purple line) for increasing Vtip (aligned with panel b).

Because the measurement captures the QD level spacings as charging peak distances ΔV, they need to be converted to Eadd via the tip lever arm αtip. The latter relates a change of Vtip to its induced shift of the QD state energies. The lever arm is determined by the ratio of the capacitance between tip and dot Ctip, and the total capacitance of the dot CΣ, thus αtip = Ctip/CΣ. CΣ includes Ctip, the capacitance between dot and back-gate, and dot and surrounding graphene. We use a Poisson solver to estimate CΣ = 16.5 ± 3.2 aF and Ctip = 8 ± 1.5 aF for our QD (Supporting Information). Hence, we find EC = e2/CΣ ≈ 10 ± 2 meV and αtip = 0.51 ± 0.03 (close to values reported for a similar system by Jung et al.[33]). Consequently charging peaks dominantly separated by EC, that is, Eadd ≈ EC because Δ ≪ EC, should exhibit ΔV = EC/(eαtip) ≈ 20 mV, which is in close agreement with the values found within quadruplets at higher occupation numbers (Figure b,c). As expected, we also find significantly larger Eadd for every fourth charging peak. In the case of clear quadruplet ordering, the orbital splittings for our QD are deduced from Δo = Eadd4 – EC4 ≈ Eadd4 – Eadd4 and we find typical values of 4–10 meV for the first few orbitals (αtip = 0.51, Figure b,c). For this estimate, we neglect the additional Zeeman splitting or an even smaller valley splitting.

We next provide a theoretical framework to elucidate the details of the QD level spectrum. The eigenstates of bulk graphene LLs (eq ) feature different wave function amplitudes on sublattices[1] A and Bwhere K and K′ denote the two inequivalent K-points of the Brillouin zone associated with the two valleys. For N ≠ 0, the LL index differs by one for the two sublattices, while for N = 0 the part of the wave function with subscript |N| – 1 vanishes, resulting in polarized sublattices for each valley. The wave functions of bulk graphene (eq ) are modified by the tip-induced potential. Assuming a radially symmetric confinement potential, the eigenstates are described by radial and angular momentum quantum numbers (nr, m), with and . Adiabatically mapping a given LL with index N on to possible combinations of nr and m yields[53]with 0 ≤ nr ≤ |N| and m ≤ |N|.

We calculate eigenstates of a 120 nm × 100 nm commensurate graphene flake on BN using third-nearest neighbor TB,[38] where the substrate interaction enters via a periodic superstructure potential and local strain effects,[39] parametrized from DFT calculations.[54,55] We approximate the amplitude Φ0el and shape of Φgrel by a classic electrostatic solution of Poisson’s equation (Figure , Supporting Information) with the tip radius rtip as fit parameter. Comparing calculated charging energies to experiment yields a plausible value of rtip ≈ 120 nm, implying a full width at half-maximum (fwhm) of the QD confinement potential of 55 nm at 7 T. We independently determine the initially free parameter EF from the position of LL0 in STS as EF = −40 ± 5 meV (Supporting Information). Accordingly, the graphene is p-doped. We note that varying EF within the stated uncertainty range (see blue horizontal bar in Figure a) leads to no qualitative changes in the predictions of our model. We use open boundary conditions to simulate the coupling of the flake to the surrounding graphene. Consequently, eigenstates will feature complex eigenvalues E = ε + iΓ/2, where the real part ε represents the resonant energies and the imaginary part Γ the coupling to the delocalized bulk states.[56] Thus, we can readily distinguish states that are spread out over the flake (large Γ) from those localized near the tip (small Γ). We color code Γ in Figure a for a calculation with the tip-induced potential centered on an AB stacked area.

Figure 4

(a) Tight binding eigenenergies of a 120 × 100 nm 2 graphene sample with open boundaries as a function of tip-induced potential amplitude Φ0el at B = 7 T with the tip-induced potential centered on an AB area (BA and AA yield very similar behavior, not shown). Line color encodes coupling to the boundary (imaginary part Γ of eigenenergies); black (red) indicates strong (weak) localization underneath the tip. States from LL±1 and the split LL0 are labeled by ±1 and 0, respectively. The LL0 splitting reduces the confining gap to E0 – E–1 ≈ 50 meV. First states crossing EF from LL±1 are highlighted in orange. Uncertainty in EF indicated as blue horizontal bar (Supporting Information). The green rectangle marks the zoom shown in panel e. (b–d) Color plot of the wave function amplitude of states marked by orange crosses in panel e. Φ0el at the crossing point ε(Φ0el) = EF is marked. Solid (dashed) white lines are line cuts along the dotted white line in panel b for contributions from sublattice A (B), as marked. All scale bars identical. (e) Zoom into area marked by a green box in panel a. Colored lines identify valley K (cyan) and K′ (purple). Orange crosses mark crossing of EF (blue dashed line) of selected states, which are displayed in panels b–d. First two orbital Δo and valley Δτ splittings marked by arrows. (f) Comparison of length scales: tip-induced potential, magenta; calculated wave function amplitude |Ψ| of first state crossing EF (same as panel b) for sublattice A (gray line) and B (dashed line); superstructure lattice constant a = 13.8 nm; magnetic length lB (7 T) = 9.7 nm.

At B = 7 T and vanishing band bending (Φ0el = 0), we find only delocalized states whose eigenenergies cluster around the bulk LL energies (eq , Figure a). As we increase Φ0el, states begin to localize at the tip and shift in energy, with smaller Γ (darker curves) pointing to stronger localization (see Figure a). Comparing hole states originating from LL–1 for negative and positive Φ0el, we find, as expected, stronger localization in case of negative Φ0el. The potential is always attractive to one kind of charge carriers that will localize underneath the tip. The other kind is repelled by the induced potential (see also ref (31)) which results in stronger coupling to the bulk. In order to classify our TB wave functions in terms of the quantum numbers N, nr, and m, we consider sublattice A and B separately. Tracing the states back to their LL of origin reveals N, constraining possible nr ≤ |N|. The value of nr is then determined by counting radial minima in the line cuts of the wave function amplitude for each sublattice (Figure b–d). The distance of the first radial maximum from the center of the wave function is finally sufficient to assign the possible m quantum numbers of the LL (eq ). Additionally, the (nr,m) combinations need to be consistent with N differing by one on the two sublattices (eq ). For instance, the line cuts in Figure b portray (0,0) and (0,1) on sublattice A and B, respectively. As expected, small angular momentum states are the first ones to localize with increasing Φ0el, which is in line with calculations by Giavaras et al.[30] Notice that the applied B naturally lifts the orbital degeneracy in QDs.[57] Delocalized states remain at bulk LL energies (red horizontal lines in Figure a).

We distinguish two regimes in the sequence of spin degenerate states crossing EF for negative Φ0el. The first regime (Figure e) exhibits Δτ ≲ Δo ≲ EC, while the second at higher Φ0el is characterized by densely spaced states, thus Δo ≈ Δτ ≪ EC. The sequence within the first regime corresponds to about five orbital pairs from valley K and K′, which is in line with about five quadruplets in our experimental spectra (see labels “4” in Figure a and ΔV sequences in Figure b,c). The quite uniform spacing of peaks for larger Vtip (Figure a) agrees with the second regime. In order to extract Δo and Δτ within the first regime, we carefully assign the valley index to the states. Using the previously determined nr and m in eq , the first state crossing EF (Figure b) features LL index NA = 0 + 1/2(0 + |0|) = 0 on sublattice A and NB= 1 + 1/2(0 + |0|) = 1 on sublattice B, as predicted by eq for a LL|1| state in valley K. The role of the sublattices interchanges for the second state crossing EF (Figure c), placing it in valley K′. Consequently, states with NA = NB – 1 and NB = NA – 1 are assigned to valleys K and K′, respectively. The calculation therefore predicts a valley splitting of about Δ1τ = 3 meV on the AB and BA areas (see Figure b,c,e). Δ2τ is comparatively large (about 12 meV) and the respective orbital splitting Δ2o is only larger by 1–2 meV (see Figure e). Consequently, additional electrons may occupy the next orbital state of one valley prior to the same orbital state of the other valley at higher occupation numbers. Hence, we limit further comparison to experiment to Δ1τ. In our TB model, the strength of the valley splitting is dominated by the sublattice symmetry breaking term due to the BN substrate.[39] The calculations also show that the radial extent of the wave functions grows for the first couple of states crossing EF, as expected for increasing m (compare Figure panels b and c to panel d), explaining the decrease of EC toward higher peak indices at fixed B (see Figure b,c).

Theory and experiment can be directly compared for the B dependence of the onset voltage of charging peaks V*. Experimentally, V* shifts toward higher |Vtip| for increasing B (Figure a), thus gating the first state to EF requires stronger band bending for higher B. Because the curves for B > 0 T are offset proportional to √B, the straight line connecting the first charging peaks reveals that the energy distance of the first state to EF scales with √B. This corresponds to the increase in bulk LL energies for N ≠ 0 (eq ), strongly suggesting those LLs as source of the confined states. This analysis is confirmed by our TB calculations, as the first crossing points of LL±1 states with the Fermi level Φ0* also shift toward higher |Φ0el | with increasing B (Figure b). While the evolution of states with Φgrel in Figure a is (approximately) symmetric with respect to Φgrel →−Φgrel, the previously discussed p-doping induces an asymmetry in Φ0* for electrons and holes (see the lines highlighted in orange in Figure a) and thus accounts for the observed asymmetry in V*. In Figure c, we compare V* and Φ0* by using the Φ0el(Vtip) dependence from the Poisson solver (see inset Figure c, Supporting Information). Care must be taken to correctly account for the work function difference between the tip and the sample: the tip’s work function (4.5–4.8 eV[36,58]) exceeds that of graphene (4.5 eV), placing electric field neutrality in the positive Vtip sector. Moreover, it definitely has to lie in between the two charging peak regimes because the QD vanishes without band bending. Using a plausible work function difference of +50 meV in Figure c leads to satisfactory agreement between the theoretical predictions for the first state crossings and the experimental V*.

Figure 5

(a) dI/dV spectroscopy in the vicinity of an AA stacked area at varying B, marked on the right. Four spatially adjacent spectra are averaged and the ones for B > 0 T are offset by a value proportional to √B. Vstab = 1000 mV, Istab = 700 pA, Vmod = 4.2 mVrms. Green lines are guides to the eye, marking the onset voltage of charging peaks V*. At 7 T, an asterisk marks the first charging peak on either side. Inset shows zoom onto marked peaks. (b) Energy of first confined hole state as a function of induced potential amplitude Φ0el for different B as marked. At larger B, states cross EF at larger Φ0el, shifting V* to larger negative Vtip. Color codes imaginary part of the eigenenergy as in Figure a. (c) Comparison between measured and calculated V*. Inset shows the required Φ0el(Vtip) for conversion, taken from a Poisson-solver (Supporting Information) and including a reasonable work function difference of ΔΦ = 50 meV between graphene and tip. Error bars for measured V* reflect typical variation of V* on AA areas across a few superstructure unit cells. Error bars for calculation arise from the uncertainty in EF. (d) Plot of the B dependence of Eadd1 ≈ EC1 of 20 spectra (semitransparent dots) in the vicinity of an AA area. Data points are recorded at integer valued B fields (in Tesla) but displayed slightly shifted to the left (electrons, blue) and to the right (holes, red) for clarity. Median values are encircled in black. (e) Histogram of Δ1τ ≈ Eadd2 – Eadd3 (experimental error below 0.2 meV) at B = 7 T for the same AA area used in panel d. Electron (blue bars) and hole (red bars) contributions are colored. (f) Calculated |Ψ| of the first confined hole state (see panel b) crossing EF at different B as marked. The state originates from LL–1, i.e., Vtip < 0 V when crossing. All scale bars identical.

Our TB simulations predict a strong reduction of Γ with increasing magnetic field, corresponding to the suppression of the radial tail of the wave function in Figure f and indicating the onset of localization between 1 and 3 T (Figure b). The first appearance of charging peaks in the experiment at around 2 T (Figure a) fits nicely. This finding is further corroborated by comparing the diameter of the LL state , being d1 = 89 nm (63 nm) for LL1 at 1 T (2 T) with the fwhm of the band bending region of 55 nm, providing an independent confirmation of the estimated Φgrel. At higher B, the diameter of the first QD state wave function is dominated by lB rather than by the width of Φgrel (Figure f). The compression of the wave function for increasing B (Figure f) also manifests itself as increase in addition energy, for instance, for Eadd1 = EC1 + Δσ in Figure d, where the increase in Eadd1 with B by about 4 meV cannot be explained by that of Δσ, being 460 μeV between 3 and 7 T. Consequently, increased Coulomb repulsion between electrons due to stronger compression and thus larger EC1 dominates Eadd1(B) . We observe a similar monotonic increase for the other Eadd with odd index i, independent of the position of the QD.

Experiment and theory also provide detailed insight into the valley splitting Δτ of the first confined states. The peaks of the first quadruplets in Figures a and 5a (see, e.g., inset) often group in doublets, suggesting sizable values of either Δτ or Δσ (Figure e,f). While Δσ is expected to be spatially homogeneous and only weakly varying between different quadruplets, the TB calculations predict strongly varying Δτ for different quadruplets (Figure e), which is in accordance with our observations in the experimental spectra. For a quantitative comparison, we focus on Eadd2, which separates the two doublets within the first quadruplet. In view of the small value of the Zeeman splitting (Δσ ≈ 800 μeV at 7 T), we approximate EC2 by Eadd3 to extract the valley splitting Δ1τ ≈ Eadd2 – Eadd3. We record 20 spectra in the vicinity of an AA stacked area at B = 7 T to obtain a histogram of Δ1τ for electrons and holes (Figure e), where Δ1τ could be determined with an experimental error smaller than 0.2 meV. The values strikingly group around the predicted Δ1τ ≈ 3 meV found in the TB calculations (Figure e), with a probable offset in the QD position relative to the tunneling tip (Supporting Information, Section 5) explaining the QD probing an area adjacent to the tunneling tip. We conclude that sizable Δτ separate quadruplets into doublets, while the smaller Δσ contributes to the odd addition energies within the doublets. Realizing such a controlled lifting of one of the two degeneracies in graphene QDs is a key requirement for 2-qubit gate operation.[2] It enables Pauli blockade in exchange driven qubits as required for scalable quantum computation approaches using graphene.[2] Our observation of valley splittings, so far elusive, provides a stepping stone toward the exploitation of the presumably large coherence time of electron spins in graphene QDs.[2−5]

In summary, we have realized graphene quantum dots without physical edges via electrostatic confinement in magnetic field using low disorder graphene crystallographically aligned to a hexagonal boron nitride substrate. We observe more than 40 charging peaks in the hole and electron sector arranged in quadruplets due to orbital splittings. The first few peaks on the hole and electron side show an additional doublet structure traced back to lifting of the valley degeneracy. Note that such a lifting is key for the use of graphene quantum dots as spin qubits.[2] Tight binding calculations quantitatively reproduce the orbital splitting energy of 4–10 meV as well as the first orbital’s valley splitting energy of about 3 meV by assuming a tip potential deduced from an electrostatic Poisson calculation. Also the onset of confinement at about 2 T is well reproduced by the calculation. Our results demonstrate a much better controlled confinement by combining magnetic and electrostatic fields than previously found in graphene. Exploiting the present approach in transport merely requires replacing the tip by a conventional electrostatic gate with a diameter of about 100 nm. Moreover, the approach allows for straightforward tuning of (i) orbital splittings by changing the gate geometry and thus the confinement potential, (ii) valley splittings based on substrate interaction, (iii) the Zeeman splitting by altering the magnetic field, and (iv) the coupling of dot states to leads or to other quantum dots by changing the magnetic field or selecting a different quantum dot state. Finally, our novel mobile quantum dot enables a detailed investigation of structural details of graphene stacked on various substrates, by spatially mapping the quantum dot energies.