Coherent population trapping by dark state formation in a carbon nanotube quantum dot

Illumination of atoms by resonant lasers can pump electrons into a coherent superposition of hyperfine levels which can no longer absorb the light. Such superposition is known as a dark state, because fluorescent light emission is then suppressed. Here we report an all-electric analogue of this destructive interference effect in a carbon nanotube quantum dot. The dark states are a coherent superposition of valley (angular momentum) states which are decoupled from either the drain or the source leads. Their emergence is visible in asymmetric current−voltage characteristics, with missing current steps and current suppression which depend on the polarity of the applied source-drain bias. Our results demonstrate coherent-population trapping by all-electric means in an artificial atom.


SUPPLEMENTARY NOTE 1 EIGENSTATES AND DARK STATES OF A CARBON NANOTUBE
We discuss many-body configurations of a CNT quantum dot described by the Hamiltonian also given in Eq. (5) in the main manuscript. Here, m = 0, 1, 2 runs over three longitudinal shells, ε 0 is the shell spacing, U accounts for charging effects, and J is the exchange interaction. The occupation operatorn m z defines the number of electrons in shell m with angular momentum z , where a summation over the spin degree of freedom is performed. The eigenstates |N, E; S, S z , L z of the CNT Hamiltonian are classified according to their occupation N with respect to a reference state |0 given below, their energy E and total spin and angular momentum quantum numbers S, S z and L z . The reference state |0 has the shell m = 0 completely full, and the upper two shells m = 1, 2 completely empty and corresponds to the N = 0 groundstate: The left (right) states are for angular momentum values z = +(−) and the up/down arrows indicate opposite spin direction. The diagram in Eq. (2) fixes a set of occupations {n m z σ }in the representation |{n m z σ } = m z σ (d † m z σ ) nm z σ |∅ . Conventionally we fix the order of the creator operators such that electrons in the lower shells are created first. Moreover, within each shell we chose the ordering d We fix the energy of the state in Eq. (2) to be E 0 . Many-body excited states with N = 0 are obtained by creating electron-hole pairs starting from this ground state configuration for fixed electron number. In Table 1 we show both the groundstate as well as all the first excited states, which have energy ε 0 − J/2 above the groundstate. The N = 1 groundstate is trivially four-fold degenerate |1, E 10 = 0 ; with these states playing a crucial role in the coherent population trapping (CPT) effect discussed in this work. In Tabs. 1, 2 the lowest excited states with N = 1, corresponding to the energies ε 0 − J, ε 0 − J/2, and ε 0 are also shown. Finally, in Tabs. 2, 3 we depict the lowest energy states with N = 2, 3. Notice that in the numerics we have considered only states up to a cut-off of 1.5ε 0 above the respective groundstate. From the many-body energies, addition energies E add = U + 3J/2 and E add 3 = U − J/2 which in turn define the heights of the Coulomb diamonds.
A. Dark states for one electron Eq. (3) allows one to construct linear combinations |DS, σ; α of the single-particle groundstates |σ, z which are decoupled at given positions r α , and hence may act as dark states (DS) for a given bias polarity. Such states have the generic form |DS, σ; α = a(r α )|σ, + b(r α )|σ, − , where the coefficients satisfy the normalization condition |a(r α )| 2 + |b(r α )| 2 = 1 and are determined through the requirement whered ασ destroys a CNT electron of spin σ at position r α . We express such operator in the angular momentum basis,d ασ = m z r α |m z σ d m z σ , and observe that the orbital part φ m z σ (r α ) = r α |m z σ of the CNT wave function is complex, φ m z σ (r α ) = |φ m z σ (r α )|e iθα(m, z ) . Furthermore, |φ m z σ (r α )| = |φ m− z −σ (r α )| due to timereversal symmetry. Then, insertion in Eq. (4) yields for the coefficients the simple form a(r α ) = e −iθα(m, ) / √ 2, and a spin-triplet of second excited states with excitation energy E 22 − E 20 = J, The one-particle dark states in Eq. (5) can block transitions to the two-electron ground-state. However, whether the blocking is effective crucially depends on the exchange energy J. In fact, as soon as the two-particles first excited doublet enters the transport window, interference is destroyed since transport through the doublet can occur. In our simulations we took indeed J Γ = Γ L + Γ R , such that the splitting of ground-and excited state is large enough to destroy interference at least partially and small enough to not see an additional excitation line appearing. Notice that in this situation the secular approximation breaks down, and the dynamics of the reduced density matrix is governed by a more general set of equations accounting also for nonsecular terms. 1 Interestingly, for J = 0 (or at least J Γ α ), the two-particle ground state, which now is a sextuplet, can form a dark state itself which blocks transitions to the one-particle ground state at lead α since 1, E 10 ; 1 2 , ± 1 2 , ± |d ασ |2, DS = 0. Again we require that this dark state is not completely decoupled from the dynamics and can be reached from the left lead 1, E 10 ; 1 2 , ± 1 2 , ± |dᾱ σ |2, DS ∝ sin( (φ R − φ L )) = 0.

IMPACT OF PRECESSION, TEMPERATURE AND RELAXATION ON COHERENT POPULATION TRAPPING
In this section we give a closer look at the role played by the precession frequencies ω L/R , Eq. (11) of the Methods, on the shape of the current as a function of the bias voltage. Furthermore, we investigate the role of temperature and inelastic relaxation. We focus on the vicinity of the 0 ↔ 1 resonance where, as seen from the comparison in Fig. 2b, the analytical expression for the current, Eq. (4) of the main text, well reproduces the numerics.
We first start with Figs. 2a, which shows the bias dependence of the numerically evaluated current for different temperatures. While the traditional step-like behavior of the current at positive bias for is temperature broadened via the Fermi function, the interference peak at negative bias is not. To understand this feature, we have depicted in Fig. 2c the precession frequencies ω L/R in the same bias voltage range of panel (a) and for the same temperatures. It is clear that the current changes occur in the correspondence of changes in the peaks in the precession frequencies. In particular, the temperature basically only changes the height of the resonance peaks and leaves the tails invariant.
We observe that at fixed temperature the broadening of the precession frequencies peaks is largely dominated by the charging energy U , since ω L/R vanish when the bias becomes of the order of U/e, as shown in Fig. 2d. Eq. (4) of the main text proofs that the current is dominated by a single precession frequency in the numerator, where the bias direction defines which one. If this precession frequency becomes zero, interference perfectly blocks the current. In the experiment U is so large that this behavior cannot be seen; excited states enter the bias window before this value of the bias voltage is reached.
In the main manuscript we have discussed two analytic expressions for the current, valid near the 0 ↔ 1 resonance and corresponding to the case of vanishing relaxation rate, and vanishing phase difference, respectively. In principle Eq. (10) from the Methods can be solved fully analytically; however the resulting expression is very lengthy and not useful here. Therefore we display the full analytical current in Fig. 3 and compare it to the limiting cases. We show the dependence of the current on the phase difference ∆φ for various values of the relaxation rate. As an example we choose a large negative bias with precession frequencies ω L = 3µeV and ω R = 2µeV. For low relaxation rates (Γ rel Γ L/R ) the simple limit of Eq. (4) is recovered at ∆φ = 0. This agreement is expected since in the shown bias range excited states are far-off in energy. At ∆φ = 0 the current is always given by the incoherent limit from Eq. (13) of the Methods section.

SUPPLEMENTARY NOTE 3 TUNNELING RATE MATRIX
We observe that in some nanotube based devices the quantum dot region is defined electrostatically, and not, as it is assumed in this work, by the presence of tunneling barriers between the metallic contact and the CNT. Electrostatic confinement is characterized by a strong variation of the tunneling coupling and charging energy with the gate voltage, leading to different behaviors ranging from Fabry-Perot interference to Kondo physics and sequential tunneling as the gate is swept from negative to positive voltages. In our experiment we did not observe strong variations of the current as the gate voltage was swept across the band gap down to the deep n-doped regime considered here (11 v to 13 V), cfg. Fig. 4. This supports our assumption that the bottleneck for transport is given by the contact region between the metallic leads and the CNT.
In this section we want to show that, in general, the rate matrix defined in Eq. (7) of the Methods, is non diagonal in the angular momentum basis. We then exemplify our considerations to the special case of a ring coupled to a metal, where some close analytical expressions for the rate matrix can be obtained.

A. Tunneling amplitude
We start from the tunneling Hamiltonian, whose form is given in Eq. (6) of the Methodŝ numerics analytics  The tunneling amplitude is proportional to the overlap of a wave function of the CNT φ m z σ (r) = r|m z σ and one of the lead ψ αkσ (r) = r|αkσ . Explicitly, αkσ|ĥ|m z σ = t αkm z σ δ σσ , whereĥ =p 2 2m el + v(r) is the single-particle Hamiltonian of the CNT-leads complex. By decomposing the electrostatic potential into a contribution from the CNT and one from the leads, v(r) = v CNT (r) + v leads (r), see the schematics in the Supplementary Figure 5, the tunneling amplitude can be written as whereĥ CNT is the single particle part of the CNT Hamiltonian from Eq. (5) in the Methods section. Since the wave functions of the CNT are much more localized than the lead ones, the contribution containing the overlap of lead and CNT wave function in the lead region (where the potential v leads is finite) can be neglected, yielding the simple expression Hence, the evaluation of the tunneling amplitude requires to take a closer look at the CNT wave functions as well as the lead wave functions in the tunneling region. Electrostatic potential along the tunneling direction, chosen to be along the x-axis. In the lead, the electrons are considered to be free electrons in the direction parallel to the surface while they experience a confinement potential in the x-direction. E b is the energy at the band bottom, EF the Fermi energy and φ0 the work function. Notice that the zero of the energy has been set to the vacuum. The CNT is located at a distance d from the lead and features localized bound states.
We start from the latter. We assume an adiabatically smooth variation of the lead surface in the contact region, such that the lead wave functions locally factorize in a contribution parallel to the surface and in an exponentially decaying part perpendicular to it:

Conservation of energy in the lead's potential well and in the tunneling region yields
2m el , where E el is the energy of the lead electron with respect to the vacuum. Moreover, the energy at the band bottom is Thus, according to Eq. (16), the smallest values of κ x , and hence the largest penetration in the CNT, are obtained when k x ≈ k F yielding κ x ≈ 2m el φ 0 / 2 . Since the total energy is bound to be E α F , this simultaneously implies that the longitudinal components k y , k z should be vanishingly small (i.e. in the vicinity of the surface Γ -point).
Regarding the CNT wave functions, we assume that they are well described as a linear combination of atomic orbitals (LCAO) localized at the atomic positions R j = (X j , Y j , Z j ). In particular, the low energy properties are already well captured by considering a single p-orbital for each atomic position. 2 We denote |jσ such atomic state and the associated wave function as p σ (r − R j ) = r|jσ . Hence, |m z σ = jσ |jσ jσ|m z σ = jσ |jσ c j (m z σ), where the LCAO coefficients c j (m z σ) ≡ jσ|m z σ have been introduced. Notice that they are chosen in such a way that the CNT wavefunctions obey proper boundary conditions at the ends of the tube. 3 Furthermore, due to time reversal symmetry, it holds c j (m z σ) = c * j (m − z − σ). It follows In the absence of spin-orbit coupling, as in our case, the spatial and spin parts factorize both for the leads as well as the CNT wave functions, yielding spin-independent coefficients c j (m z σ) = c j (m z ). Similarly, the scalar product αkσ|jσ becomes spin independent, yielding the final form for the tunneling amplitude where in the last step we approximated the localized p-orbitals to Dirac-delta functions, p(r) = aδ(r), centered at the atomic position R j . Here a is a normalization factor. The last approximation neglects the nodal plane of the p z orbitals, but it is justified by i) the selection of j given by the lead wave function and ii) the negligible contribution to the integral given by the CNT wave function inside the tube.

B. Rate matrix of a CNT-metal complex
We calculate the rate matrix according to Eqs. (11) and (18), i.e., in the absence of spin-orbit coupling. We then obtain Whether the rate matrix is diagonal in the angular momentum basis, crucially depends on the geometry of the contact region. The exponential e −κx(Xj +X j ) in fact selects in the sums over the atomic positions those CNT atoms closest to the leads. Furthermore, in the summation over the momenta k, it selects the smallest values of κ x compatible with the requirement that the energy of the tunneling lead electron is resonant with the CNT chemical potential ∆E. As discussed in the previous subsection, this yields k y , k z ≈ 0 and κ x ≈ 2m el φ 0 / 2 := κ min , such that Thus, this so called surface Γ-point approximation 4 enables us to decouple the sums over j and j into two independent sums. We introduce the density of states at the Fermi level D = k δ(ε k − E F ) and the tunneling coefficients yielding In the surface Γ-point approximation, the wave function ψ (Y j , Z j ) is independent of k y and k z and hence real. Furthermore, the LCAO coefficients are related by time-reversal symmetry, c j (m, z ) = c * j (m, − z ), yielding the result τ α (m z ) = τ * α (m − z ). Accounting for this symmetry we finally obtain the final form for the rate matrix used in the main text The result in Equation (23) strongly relies on the surface Γ-point approximation, which allows one to decouple the double sum over the atomic positions j and j . In the next subsection we have explored the consequences of keeping a finite contribution for the parallel momenta k y and k z on the example of a ring of N carbon atoms coupled in three different ways to a metal. As we shall see, if the ring is lying flat on the substrate, such that all atoms are equally distant from the lead, the rate matrix becomes diagonal, see Supplementary Figure 6a. The result in Eq. (23) is in contrast recovered when the ring is orthogonal to the substrate, in a way that tunneling is dominated by only one closest atom, Supplementary Figure 6b. When two atoms are equally close to the surface, as for the case in figure 6c, the rate matrix is off-diagonal, but the modulus of the diagonal elements is smaller than that of the diagonal ones.
From this we conclude that in the CNT case, where only few atoms are close to the leads, the rate matrix is not diagonal. How good the simple form Eq. (23) describes the experiment, depends on various factors, among which tube's chirality. In the case of our experiment, we consider CNTs of the zig-zag class, which at the tube's end have non vanishing weights only for atoms of a given sublattice. 3 Thus, if at the left end only A atoms have non vanishing LCAO coefficients, this implies that the neighboring B atoms are not tunneling coupled, hence effectively achieving the situation described in the simple example in Fig. (6)b. At the right lead, the same considerations apply upon exchange of the role of atoms A and B.

C. Rate matrix of a ring-metal complex
In this subsection we calculate the rate matrix for a ring of N carbon atoms with radius R. To study the effects of the orientation of the ring with respect to the surface, we study three different configurations. In the first one, the ring is lying parallel to the surface of the lead, at a distance d from it, as shown in the Supplementary Figure 6a; in the second one it is standing on the x − y plane perpendicular to the lead plane like shown in Figure 6b; in the third one it is also standing but now in a way that two atoms are equally distant from the lead, as shown in Figure 6c. The rate matrix for the ring follows easily from the general expression Eq. (19) upon dropping the shell index m. Similar to the CNT, the ring has a C N symmetry and its single particle eigenstates can be classified in terms of angular momentum, z , and spin, σ, degrees of freedom. The LCAO coefficients follow from the diagonalization of the ring Hamiltonian, and have the form c j ( z ) = j| z = 1 √ N e i 2π N j z . Eq. (19) yields then for the ring the rate matrix To further simplify the calculation we consider in the following a plane wave behavior for the parallel wave function, ψ αk (R) = 1 √ S e ik ·R , with S a normalization constant.
Case a Let us consider the first case where the ring is lying planar on top of the lead at a distance X j = d. The rate matrix then reads where C =: |a| 2 /(SL x ). It is convenient to express j = j −j +j := −∆j +j and to observe that |R j −R j | := R ∆j only depends on the relative distance ∆j but not on the position j. This suggests to transform the sum over momentum into an integral, and to express this integral in cylindrical coordinates k y , k z → ϕ, k . This results in The integration over the angle ϕ results in a real function of k ; similarly, in the sum over ∆j for each finite positive ∆j there is a negative counterpart. This leaves us with the diagonal and real rate matrix which therefore does not support dark states. The absolute value of the diagonal parts is not important in this consideration. For simplicity and homogeneity with the other cases we have kept the δ approximation for the p z functions. The description of a CNT "slice" would rather require to consider the orbital structure of the radially distributed p orbitals. Eq. (28) is obtained, though, also out of more fundamental symmetry arguments: whereĈ N is the rotation of 2π/N around the x axis and the isotropy of the leads is assumed. Eq. (29) implies Γ α is diagonal. The form in Eq. (28) follows by requiring time reversal symmetry.

Case b
In the second case the result is quite different. The rate matrix for the standing ring is The ring is standing perpendicular to the lead with only one carbon atom closest to the surface. c Similar to the previous case, the ring is standing perpendicular to the surface but rotated such that two carbon atoms are equidistant from the lead.
where C = |a| 2 /(SL x ). One can see immediately that the trick used in the previous case a does not work here, due to the dependence on j (j ) of the variables X j (X j ). An estimation of κ x = 2m el for typical work functions φ 0 = O(eV) tells that the contribution to the rate matrix shrinks by one order of magnitude for a distance of 1Å of the atom to the lead surface. This suggests that perfect local tunneling to the atom closest to the lead j = j = J at distance X J = d is a good approximation. We then obtain Thus the single atom contact yields a rate matrix with maximal coherence, like in the surface Γ-point approximation discussed in the previous subsection.

Case c
In the third case the ring is rotated in a way that tunneling can occur through two atoms J and J which are both in contact with the lead (X J = X J = d). The rate matrix reads The diagonal and off-diagonal elements of the rate matrix can be simplified to with ∆Y = Y J − Y J . We used the fact that the sum over k is isotropic and therefore sin(k y ∆Y ) → 0. One can see directly that for ∆Y = 0 the result of case b is recovered. For ∆Y = 0 the amplitude of the off-diagonal terms in Γ α is smaller than the diagonal values since in general 1 + cos x cos y − cos x − cos y = (cos x − 1)(cos y − 1) ≥ 0. We obtain a rate matrix intermediate to cases a and b