Lasing and antibunching of optical phonons in semiconductor double quantum dots

Figure 1(a) depicts our model of a DQD embedded in a semiconductor substrate, in which two single-level quantum dots, L and R, are connected by tunnel coupling V_C. The energy levels, ε_L and ε_R, are electrically tunable. We choose ε_R = −ε_L and denote the level spacing ε_L − ε_R by Δ. We assume that the total number of electrons in the DQD is restricted to one or zero due to Coulomb blockade. The electron couples to LO phonons of energy ℏω_ph in the substrate by the Fröhlich interaction. Our system Hamiltonian is $\mathcal {H} = \mathcal {H}_{\mathrm {e}} + \mathcal {H}_{\mathrm { ph}} + \mathcal {H}_{\mathrm { ep}}$,

$$\begin{equation} \mathcal{H}_{\mathrm{e}} = \frac{\Delta}{2} (n_L - n_R) + V_{\rm C} (d_L^\dagger d_R + d_R^\dagger d_L), \end{equation} \tag{ 1 }$$

$$\begin{equation} \mathcal{H}_{\rm ph} = \hbar \omega_{\rm ph} \sum_{\boldsymbol{q}} N_{\boldsymbol{q}}, \end{equation} \tag{ 2 }$$

$$\begin{equation} \mathcal{H}_{\rm ep} = \sum_{\alpha = L, R} \sum_{\boldsymbol{q}} M_{\alpha, \boldsymbol{q}} (a_{\boldsymbol{q}} + a_{-\boldsymbol{q}}^\dagger) n_\alpha \end{equation} \tag{ 3 }$$

using creation (annihilation) operators d^†_α (d_α) for an electron in dot α and a^†_q (a_q) for a phonon with wavevector q. n_α = d^†_αd_α and N_q = a^†_qa_q are the number operators. The spin index is omitted for electrons. The coupling constant is given by

$$\begin{equation} M_{\alpha, \boldsymbol{q}} = \sqrt{\frac{e^2 \hbar \omega_{\rm ph}}{2 \mathcal{V}} \left[ \frac{1}{\epsilon(\infty)} - \frac{1}{\epsilon(0)} \right]} \frac{1}{q} \int {\rm d}\boldsymbol{r} |\psi_{\alpha} (\boldsymbol{r})|^2 {\rm e}^{{\rm i} \boldsymbol{q} \cdot \boldsymbol{r}}, \end{equation} \tag{ 4 }$$

where (∞) ((0)) is the dielectric constant at high (low) frequency, $\mathcal {V}$ is the volume of substrate and ψ_α(r) is the electron wavefunction in dot α of radius $\mathcal {R}$. The LO phonons only around the Γ point, such as $|\bi {q}| {\lesssim } 1/\mathcal {R}$, are coupled to the DQD because of an oscillating factor in the integral over |ψ_α(r)|². This fact justifies the dispersionless phonons in $\mathcal {H}_{\mathrm { ph}}$. We assume equivalent quantum dots L and R, whence M_R,q = M_L,q e^iq·r_LR with r_LR being a vector joining their centers.

Zoom In Zoom Out Reset image size

In $\mathcal {H}_{\mathrm { ep}}$, an electron in dot α couples to a single mode of phonon described by

$$\begin{equation} a_\alpha = \frac{\sum_{\boldsymbol{q}} M_{\alpha, \boldsymbol{q}} a_{\boldsymbol{q}}} {(\sum_{\boldsymbol{q}} |M_{\alpha, \boldsymbol{q}}|^2)^{1/2}}. \end{equation} \tag{ 5 }$$

We perform a unitary transformation for phonons from a_q to S- and A-phonon modes:

$$\begin{equation} a_{\rm S} = \frac{a_{L} + a_{R}}{\sqrt{2 + (\mathcal{S} + \mathcal{S}^*)}}, \quad a_{\rm A} = \frac{a_{L} - a_{R}}{\sqrt{2 - (\mathcal{S} + \mathcal{S}^*)}} \end{equation} \tag{ 6 }$$

and others orthogonal to a_S and a_A, where $\mathcal {S}$ is the overlap integral between a_L and a_R phonons in equation (5).⁵ Since the phonon modes other than S- and A-phonons are decoupled from the DQD, they can be safely disregarded in the electron transport through the DQD. Hence, we obtain the effective Hamiltonian

$$\begin{equation}\fl H=\mathcal{H}_{\rm e} + \hbar \omega_{\rm ph} \left[ N_{\rm S} + \lambda_{\rm S} (a_{\rm S}+ a_{\rm S}^\dagger) (n_L + n_R) \right] + \hbar \omega_{\rm ph} \left[ N_{\rm A} + \lambda_{\rm A} (a_{\rm A}+ a_{\rm A}^\dagger) (n_L - n_R) \right], \end{equation} \tag{ 7 }$$

where N_S = a^†_Sa_S and N_A = a^†_Aa_A, with dimensionless coupling constants

$$\begin{equation} \lambda_{\rm S/A} = \frac{1}{2 \hbar \omega_{\rm ph}} \left( \sum_{\boldsymbol{q}} |M_{L, \boldsymbol{q}} \pm M_{R, \boldsymbol{q}}|^2 \right)^{1/2}. \end{equation} \tag{ 8 }$$

The mode functions for S- and A-phonons are shown in figure 1(b) along a line through the centers of the quantum dots. The definition and calculation of the mode functions are given in appendix A. Since the phonons are dispersionless, they do not diffuse and act as a cavity including the DQD.⁶ A-phonons play a crucial role in the phonon-assisted tunneling between the quantum dots and thus in the phonon lasing, as discussed below, whereas S-phonons do not since it couples to the total number of electrons in the DQD, n_L + n_R. Both phonons are relevant to the Franck–Condon effect.

Our Hamiltonian H in equation (7) is applicable to DQDs fabricated in a semiconductor substrate, where ℏω_ph = 36 meV and λ_S,A = 0.1–0.01 for $\mathcal {R} = 10\mbox {--}100\,{\mathrm { nm}}$ in GaAs. It also describes a DQD in a suspended CNT when an electron couples to a vibron, longitudinal stretching mode with ℏω_ph ∼ 1 meV, λ_A ≳ 1 and λ_S = 0 in experimental situations [14, 15], as shown in appendix B.

The DQD is connected to external leads L and R in series, which enables electronic pumping by the electric current under a finite bias. The tunnel coupling between lead L and dot L is denoted by Γ_L and that between lead R and dot R by Γ_R. We also introduce the phonon decay rate Γ_ph to take into account natural decay of LO phonons into the so-called daughter phonons due to lattice anharmonicity [19]. We describe the dynamics of the DQD-phonon density matrix ρ using a Markovian master equation

$$\begin{equation} \dot \rho = - \frac{\rm i}{\hbar} [H, \rho] + \mathcal{L}_{\rm e} \rho + \mathcal{L}_{\rm ph} \rho, \end{equation} \tag{ 9 }$$

where $\mathcal {L}_{\mathrm { e}}$ and $\mathcal {L}_{\mathrm { ph}}$ describe the electron tunneling between the DQD and leads and the phonon decay, respectively. $\mathcal {L}_{\mathrm { e}}$ is written as

$$\begin{eqnarray}&&\fl \mathcal{L}_{\rm e} \rho = \sum_{\alpha = L, R; i, j} \frac{\Gamma_\alpha}{2} \Bigl[ f_\alpha(\epsilon_i - \epsilon_j) \Bigl( |i \rangle \langle i| d_\alpha^\dagger |j \rangle \langle j| \rho d_\alpha + d_\alpha^\dagger \rho |j \rangle \langle j| d_\alpha |i \rangle \langle i| \nonumber \\ &&- \rho |j \rangle \langle j| d_\alpha |i \rangle \langle i|d_\alpha^\dagger - d_\alpha |i \rangle \langle i| d_\alpha^\dagger |j \rangle \langle j| \rho \Bigr) \nonumber \\&&+ \skew4\bar{f}_\alpha(\epsilon_i - \epsilon_j) \Bigl( |j \rangle \langle j| d_\alpha |i \rangle \langle i| \rho d_\alpha^\dagger + d_\alpha \rho |i \rangle \langle i| d_\alpha^\dagger |j \rangle \langle j| \nonumber \\ &&- \rho |i \rangle \langle i| d_\alpha^\dagger |j \rangle \langle j| d_\alpha - d_\alpha^\dagger |j \rangle \langle j| d_\alpha |i \rangle \langle i| \rho \Bigr) \Bigr] \end{eqnarray} \tag{ 10 }$$

with |i〉 and _i being an eigenstate of H and the corresponding energy eigenvalue and f_α() $[\skew4\bar {f}_\alpha (\epsilon ) = 1 - f_\alpha (\epsilon )]$ being the Fermi distribution function for electrons (holes) in lead α [25]. The Fermi levels in leads L and R are given by μ_L = eV /2 and μ_R = −eV /2, respectively, with bias voltage V between the leads.

In the limit of large bias voltage, $\mathcal {L}_{\mathrm { e}}$ is reduced to

$$\begin{equation} \mathcal{L}_{\rm e} \rho = ( \Gamma_L \mathcal{D} [d_L^\dagger] + \Gamma_R \mathcal{D} [d_R] ) \rho, \end{equation} \tag{ 11 }$$

where $\mathcal {D}[A] \rho = A \rho A^\dagger - \frac {1}{2} \{ \rho , A^\dagger A \}$ is a Lindblad dissipator. In this case, an electron tunnels into dot L from lead L with tunneling rate Γ_L and tunnels out from dot R to lead R with Γ_R in one direction. We examine this situation in the main part of this paper. With finite bias voltages, equation (10) is evaluated in sections 3.3 and 4.2, where electron tunneling takes place in both directions unless eV is far beyond the temperature T.

The phonon dissipator $\mathcal {L}_{\mathrm { ph}}$ is given by

$$\begin{equation} \mathcal{L}_{\rm ph} \rho = \Gamma_{\rm ph} (\mathcal{D} [a_{\rm S}] + \mathcal{D} [a_{\rm A}]) \rho \end{equation} \tag{ 12 }$$

on the assumption that the temperature T in the substrate is much smaller than ℏω_ph, and daughter phonons immediately escape from the surroundings of DQD.

In the following, we adopt the BMS approximation [26] to equation (9). We diagonalize the Hamiltonian H in equation (7) and set up the rate equation in the energy eigenbasis,

$$\begin{equation} \dot{P_i} = \sum_{j} L_{ij} P_j \end{equation} \tag{ 13 }$$

for probability P_i to find the system in eigenstate |i〉. Here, $L_{ij} = \langle i| [ (\mathcal {L}_{\mathrm { e}} + \mathcal {L}_{\mathrm { ph}}) |j \rangle \langle j| ] |i \rangle .$ The solution of equation (13) with $\dot {P_i}=0$ determines the steady-state properties. The condition to justify the BMS approximation will be given in section 6.

First, we present our numerical results in the case of Γ_L,R ≫ Γ_ph. We consider the limit of large bias voltage. Figure 2(a) shows the current I through the DQD as a function of level spacing Δ, with λ_A = 0.1 (solid line) and 1 (dotted line). Beside the main peak at Δ = 0, we observe subpeaks at Δ = Δ_ν ≃ νℏω_ph (ν = 1,2,3,...) due to the phonon-assisted tunneling⁷. At the νth subpeak, electron transport through DQD is accompanied by the emission of ν phonons. As a result, the phonon number is markedly enhanced at the subpeaks, as shown in figure 2(b), in both cases of λ_A = 0.1 and 1. However, the physics is very different for the two cases, as we will show below.

Zoom In Zoom Out Reset image size

For λ_A = 0.1 and Δ = Δ_ν, the electronic state |L〉_e with n phonons is coherently coupled to |R〉_e with (n + ν) phonons [27], similarly to cavity QED systems, if the lattice distortion is neglected. To examine the amplification of A-phonons, we calculate the phonon autocorrelation function

$$\begin{equation} g^{(2)}_{\rm A} (\tau) = \frac{\langle :N_{\rm A} (0) N_{\rm A} (\tau): \rangle}{\langle N_{\rm A} \rangle^2}. \end{equation} \tag{ 14 }$$

The numerator includes the normal product :N_A(0)N_A(τ): = a^†_A(0)a^†_A(τ)a_A(τ)a_A(0). g⁽²⁾_A(τ) is proportional to the probability of phonon emission at time τ on the condition that a phonon is emitted at time 0 [28, 29]. A value of g⁽²⁾_A(0) = 1 indicates a Poisson distribution of phonons which is a criterion of phonon lasing, whereas g⁽²⁾_A(0) < 1 (g⁽²⁾_A(0) > 1) represents the phonon antibunching (bunching). We thus find phonon lasing at the current subpeaks in figure 2(c) in the case of λ_A = 0.1 (solid line).

When λ_A = 1, the strength of the electron–phonon interaction is comparable with the phonon energy. In this case, the lattice distortion by the Franck–Condon effect seriously disturbs the above mentioned coherent coupling between an electron and phonons in the DQD and, as a result, suppresses the phonon lasing. Indeed, g⁽²⁾_A(0) > 1 at the current subpeaks, indicating phonon bunching.

To compare the two situations in detail, we present the number distribution of A-phonons in figures 3(a) and (b) at the current main peak and subpeaks. In the case of λ_A = 0.1, a Poisson-like distribution emerges at the subpeaks, whereas a Bose distribution with effective temperature T* is seen at the main peak. T* is determined from the number of phonons 〈N_A〉 in the stationary state as 1/[e^{ℏω_ph/(k_BT}*)  − 1] = 〈N_A〉. When λ_A = 1, on the other hand, the distribution shows an intermediate shape between Poisson and Bose distributions at the subpeaks and the Bose distribution at the main peak.

Zoom In Zoom Out Reset image size

In figures 4(a) and (b), we plot the autocorrelation function g⁽²⁾_A(τ) as a function of τ. In the case of λ_A = 0.1, g⁽²⁾_A(τ) ≃ 1, regardless of the time delay τ, which supports the phonon lasing at the current subpeaks. At the main peak, g⁽²⁾_A(τ) ≃ 1 + e^−Γ_phτ. This is a character of thermal phonons with temperature T*. When λ_A = 1, we find intermediate behavior, g⁽²⁾_A(τ) ≃ 1 + δ_ν e^−Γ_phτ (0 < δ_ν < 1), at the νth subpeak. This indicates that the phonons are partly thermalized by the Franck–Condon effect. For larger ν, the distribution is closer to the Poissonian with smaller δ_ν.

Zoom In Zoom Out Reset image size

To elucidate the competition between phonon lasing and thermalization by the Franck–Condon effect, we analyze the rate equation in equation (13), focusing on the current peaks in the large bias-voltage limit. We introduce polaron states |L(R),n〉_eA for an electron in dot L (R) and n phonons with lattice distortion:

$$\begin{equation} |L, n \rangle_{\rm eA} = |L \rangle_{\rm e} \otimes \mathcal{T}_{\rm A} |n \rangle_{\rm A}, \quad |R, n \rangle_{\rm eA} = |R \rangle_{\rm e} \otimes \mathcal{T}_{\rm A}^{\dagger} |n \rangle_{\rm A}, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} \mathcal{T}_{\rm A} = {\rm e}^{- \lambda_{\rm A} (a_{\rm A}^\dagger - a_{\rm A})} \end{equation} \tag{ 16 }$$

and its Hermitian conjugate $\mathcal {T}_{\mathrm { A}}^\dagger$ describes the shift of equilibrium position of the lattice when an electron stays in dots L and R. Note that the lattice distortion produces λ²_A extra phonons: _eA〈α,n|N_A|α,n〉_eA = n + λ²_A. When Δ = Δ_ν (≃ νℏω_ph), the eigenstates of Hamiltonian H are given by the zero-electron states |0,n〉_eA = |0〉_e⊗|n〉_A, bonding and anti-bonding states between the polarons,

$$\begin{equation} |\pm, n \rangle_{\rm eA} = \frac{1}{\sqrt2} (|L, n \rangle_{\rm eA} \pm |R, n + \nu \rangle_{\rm eA}) \end{equation} \tag{ 17 }$$

(n = 0,1,2,...), and polarons localized in dot R, |R,n〉 (n = 0,1,2,...,ν − 1). This is a good approximation provided that V_C ≪ ℏω_ph. The rate equations for these states are

$$\begin{eqnarray}&&\fl \dot{P}_{0, n} = -\Gamma_L P_{0, n} + \sum_{m = 0}^{\infty} \frac{\Gamma_R}{2} |_{\rm A} \langle n | \mathcal{T}_{\rm A}^\dagger | m + \nu \rangle_{\rm A} |^2 P_{{\rm mol}, m} + \sum_{m = 0}^{\nu - 1} \Gamma_R |_{\rm A} \langle n | \mathcal{T}_{\rm A}^\dagger | m \rangle_{\rm A} |^2 P_{R, m} \nonumber\\&&+ \Gamma_{\rm ph} \left[ (n + 1) P_{0, n + 1} - n P_{0, n} \right], \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\fl \dot{P}_{{\rm mol}, n} = -\frac{\Gamma_R}{2} P_{{\rm mol}, n} + \sum_{m = 0}^{\infty} \Gamma_L |_{\rm A} \langle n | \mathcal{T}_{\rm A} | m \rangle_{\rm A} |^2 P_{0, m} + \Gamma_{\rm ph} \left[ \left( n + 1 + \frac{\nu}{2} \right) P_{{\rm mol}, n + 1} - \left( n + \frac{\nu}{2} \right) P_{{\rm mol}, n} \right],\nonumber\\ \end{eqnarray} \tag{ 19 }$$

where P_mol,n = P_+,n + P_−,n (n = 0,1,2,...) and

$$\begin{equation} \dot{P}_{R, n} = - \Gamma_R P_{R, n} + \Gamma_{\rm ph} \bigl[ (n + 1) P_{R, n + 1} - n P_{R, n} \bigr], \end{equation} \tag{ 20 }$$

with P_R,ν = P_mol,0/2 (n = 0,1,2,...,ν − 1). As shown in appendix C, these equations yield the current I and electron number in the DQD, 〈n_e〉 = 〈n_L + n_R〉, in terms of the number of polarons localized in dot R, $\langle \tilde n_R \rangle = \sum _{n = 0}^{\nu - 1} P_{R, n}$, as

$$\begin{equation} I = e \Gamma_R \frac{ 1 + \langle \tilde n_R \rangle }{2 + \gamma}, \quad \langle n_{\rm e} \rangle = \frac{2 - \gamma \langle \tilde n_R \rangle} {2 + \gamma} \end{equation} \tag{ 21 }$$

with γ = Γ_R/Γ_L. The number of A-phonons is given by

$$\begin{equation} \langle N_{\rm A} \rangle = (\nu + 2 \lambda_{\rm A}^2) \frac{I}{e \Gamma_{\rm ph}} + \lambda_{\rm A}^2 \langle n_{\rm e} \rangle. \end{equation} \tag{ 22 }$$

The first two terms in equation (22) indicate the emission of ν phonons by phonon-assisted tunneling (from dot L to dot R) and creation of 2λ²_A phonons by lattice distortion (with two tunnelings between the DQD and leads) per transfer of a single electron through the DQD. The last term describes the average number of polarons 〈n_e〉 in the stationary state.

When Γ_L,R ≫ Γ_ph, we obtain

$$\begin{equation} I = I_0 + \mathcal{O}(\Gamma_{\rm ph}/\Gamma_{L,R}), \end{equation} \tag{ 23 }$$

where I₀ = eΓ_R/(2 + γ) is the current at the main peak in the absence of electron–phonon interaction, and

$$\begin{equation} g^{(2)}_{\rm A} (0) = \frac{\nu + 4 \lambda_{\rm A}^2}{\nu + 2 \lambda_{\rm A}^2} + \mathcal{O} (\Gamma_{\rm ph}/\Gamma_{L,R}). \end{equation} \tag{ 24 }$$

These explain the numerical results in figure 2 at the current subpeaks. The formula in equation (24) indicates g⁽²⁾_A(0) ≃ 1 (phonon lasing) for λ²_A ≪ ν and g⁽²⁾_A(0) ≃ 2 (thermalized phonons by the lattice distortion) for λ²_A ≫ ν. In the latter case, the phonons follow Bose distribution with T* to deduce 〈N_A〉 in equation (22).

We comment on the peak width of the electric current in figure 2(a). The electron transfer around the νth current peak is dominated by the tunneling between polaron states |L,n〉_eA and |R,n + ν〉_eA with n ≃ 〈N_A〉. Thus, the peak width is determined by the effective tunnel coupling

$$\begin{eqnarray} W_\nu &=& |_{\rm eA} \langle R, n + \nu| H_{\mathrm{e}} | L, n \rangle_{\rm eA}| _{n \simeq \langle N_{\rm A} \rangle} \nonumber \\ & =& \left| \frac{n!}{(n + \nu)!} (-2 \lambda_{\rm A})^\nu {\rm e}^{-2 \lambda_{\rm A}^2} L^\nu_n (4 \lambda_{\rm A}^2) V_{\rm C} \right| _{n \simeq \langle N_{\rm A} \rangle} \end{eqnarray} \tag{ 25 }$$

(ν = 0,1,2,...), where L^ν_n(x) is the Laguerre polynomial⁸. The factor of e^{−2λ2 _A} in equation (25) stems from the electron localization by dressing the phonons in forming the polarons. This explains the narrower subpeaks in the case of λ_A = 1 than that of λ_A = 0.1.

When λ²_A ≪ 1, equation (25) yields

$$\begin{equation} W_\nu \simeq \sqrt{\langle N_{\rm A} \rangle + 1} \lambda_{\rm A}^\nu V_{\rm C}. \end{equation} \tag{ 26 }$$

This is in quantitative accordance with the peak widths in the case of λ_A = 0.1 (solid line in figure 2(a)).

So far we have considered the large bias-voltage limit. In this subsection, we examine the case of finite bias voltages to elucidate the Franck–Condon blockade [14] in our system. Figures 5(a) and (b) show the electric current as a function of bias voltage V when the level spacing Δ is tuned to the main and subpeaks in figure 2(a). The electron–phonon coupling is (a) λ_A = 0.1 and (b) 1. At the main peak (Δ = 0) in case (a), the current is almost identical to I₀ in the large bias-voltage limit when μ_L = eV /2 exceeds the interdot tunnel coupling V_C. (The current vanishes when eV /2 ≲ V_C, reflecting the formation of bonding and antibonding orbitals at energy level ±V_C, from two orbitals at ε_L = ε_R = 0 in the DQD.) The influence of electron–phonon coupling is hardly observable. At the subpeaks (Δ = Δ_ν ≃ νℏω_ph, ν = 1,2,3) in case (a), on the other hand, the current is suppressed at small V and it increases stepwise to the value in the large V limit. This is due to electron–phonon coupling, as explained below. The current suppression is much more prominent in the case of (b) with larger λ_A. We observe the suppression even at the main peak in this case.

Zoom In Zoom Out Reset image size

The reason for the current suppression is as follows. When an electron tunnels between the DQD and leads, the equilibrium position of the lattice is suddenly changed to form the polaron, |L,n〉_eA or |R,n〉_eA, in equation (15). While all phonon states participate in polaron formation in the large bias-voltage limit, the phonon states are limited under finite bias voltages due to energy conservation. This weakens the tunnel coupling between the DQD and leads and also between the quantum dots, which is known as the Franck–Condon blockade. In figures 5(a) and (b), the current increases stepwise as μ_L = eV /2 increases by ℏω_ph because higher-energy states become accessible (Franck–Condon steps) and converges to I = I₀ in the large bias-voltage limit. The larger voltage is required to lift off the Franck–Condon blockade for larger λ_A [16].

Figures 5(c)–(f) show the phonon number and autocorrelation function g⁽²⁾_A(0) as a function of V . The phonon number shows the Franck–Condon steps in both cases of (c) λ_A = 0.1 and (d) 1. The autocorrelation function, on the other hand, is qualitatively different for the two cases. In figure 5(e) with λ_A = 0.1, g⁽²⁾_A(0) ≃ 1 even at the first Franck–Condon step except for anomalous behavior around the beginning of the step. This indicates that phonon lasing is robust against the current suppression by the Franck–Condon blockade and hence it is observable under finite bias. In figure 5(f) with λ_A = 1, g⁽²⁾_A(0) changes slowly with V , reflecting V -dependence of the thermalization due to the Franck–Condon effect.

In sections 3.1–3.3, we have restricted ourselves to the case of Γ_L,R ≫ Γ_ph to examine phonon lasing. If tunnel coupling is tuned to be Γ_L,R ≲ Γ_ph, we observe another phenomenon, antibunching of A-phonons [30]. Figure 6(a) presents a color-scale plot of g⁽²⁾_A(0) in the λ_A–(Γ_L,R/Γ_ph) plane when Δ is tuned to be at the first current subpeak (Δ = Δ₁ ≃ ℏω_ph). We assume that Γ_L = Γ_R ≡ Γ_L,R, λ_S = 0, and large limit of bias voltage. At λ_A = 0.05 and Γ_L,R/Γ_ph = 0.1, for example, g⁽²⁾_A(0) ≪ 1, representing a strong antibunching of phonons. This is because the phonon emission is regularized by the electron transport through the DQD. In figure 6(b), we plot the autocorrelation function of the electric current

$$\begin{equation} g^{(2)}_{\rm current} (\tau) = \frac{\langle :n_R(0) n_R(\tau): \rangle}{\langle n_R \rangle^2}, \end{equation} \tag{ 27 }$$

where n_R is the electron number in dot R. It fulfill g⁽²⁾_current(0) = 0, indicating the antibunching of electron transport, since dot R is empty just after the electron tunnels out [29]. Remarkably, g⁽²⁾_A(τ) almost coincides with g⁽²⁾_current(τ). When Γ_L,R ≪ Γ_ph, the emitted phonon escapes from the natural cavity soon after the electron tunneling between the quantum dots. Thus, the stimulated emission for the lasing does not take place.

Zoom In Zoom Out Reset image size

At strong couplings of λ_A ≳ 1, neither phonon antibunching nor phonon lasing can be observed because of effective phonon thermalization due to the Franck–Condon effect. More than one phonon is created by the polaron formation, which spoils the regularized phonon emission by single electron tunneling and results in phonon bunching.

Even with small λ_A, bunched phonons are emitted if Γ_L,R/Γ_ph is too small. Then the number of phonons created by the tunneling is exceeded by that accompanied by the polaron staying in dot R (the first two terms are much smaller than the last term in equation (22)), as discussed in appendix (C.4.). The analytical expression of g⁽²⁾_A(0) is also given for Γ_L,R ≪ Γ_ph in the appendix.

Figure 6(c) shows a color-scale plot of g⁽²⁾_A(0) when Δ is tuned to be at the second current subpeak (Δ = Δ₂ ≃ 2ℏω_ph). The antibunching does not occur even when Γ_L,R ≪ Γ_ph because two phonons are emitted simultaneously by electron tunneling, which are bunched to each other.

We begin with the large bias-voltage limit. The electric current has a single-peaked structure as a function of Δ (Lorentzian with center at Δ = 0 and width of $V_{\mathrm { C}}\sqrt {2+\gamma }$, as will be seen in equation (30)). We do not observe subpeaks at Δ ≃ νℏω_ph since S-phonons are not relevant to the phonon-assisted tunneling between the quantum dots because they couple to the total number of electrons, n_L + n_R in the DQD. The polaron states involving S-phonons are given by

$$\begin{equation} |L, n \rangle_{\rm eS} = |L \rangle_{\rm e} \otimes \mathcal{T}_{\rm S} |n \rangle_{\rm S}, \quad |R, n \rangle_{\rm eS} = |R \rangle_{\rm e} \otimes \mathcal{T}_{\rm S} |n \rangle_{\rm S} \end{equation} \tag{ 28 }$$

for an electron in dot L or R, with n phonons, where the lattice distortion

$$\begin{equation} \mathcal{T}_{\rm S} = {\rm e}^{- \lambda_{\rm S} (a_{\rm S}^\dagger - a_{\rm S})} \end{equation} \tag{ 29 }$$

is common for |L,n〉_eS and |R,n〉_eS. S-phonons show neither phonon lasing nor antibunching.

We derive the rate equation for arbitrary level spacing Δ in appendix D. By tracing out S-phonon degrees of freedom, we obtain the reduced rate equation for electrons, which is the same as that in the absence of electron–phonon coupling. We obtain the electric current and electron number in the DQD,

$$\begin{equation} I = \frac{e \Gamma_R}{(\Delta / V_{\rm C})^2 + 2 + \gamma}, \quad \langle n_{\rm e} \rangle = \frac{(\Delta / V_{\rm C})^2 + 2} {(\Delta / V_{\rm C})^2 + 2 + \gamma}. \end{equation} \tag{ 30 }$$

The number of S-phonons is given by

$$\begin{equation} \langle N_{\rm S} \rangle = 2 \lambda_{\rm S}^2 \frac{I}{e \Gamma_{\rm ph}} + \lambda_{\rm S}^2 \langle n_{\rm e} \rangle. \end{equation} \tag{ 31 }$$

The first term in equation (31) indicates the creation of 2λ²_S phonons by the lattice distortion with two tunnelings between the DQD and leads per single electron transfer through the DQD. The second term describes the average number of polarons. In contrast to equation (22) for A-phonons, S-phonons are not created by the interdot tunneling.

We also examine the autocorrelation function

$$\begin{equation} g^{(2)}_{\rm S} (\tau) = \frac{\langle : N_{\rm S}(0) N_{\rm S}(\tau):\rangle}{\langle N_{\rm S} \rangle^2}. \end{equation} \tag{ 32 }$$

g⁽²⁾_S(0) is independent of λ_S, for arbitrary Δ. When Γ_L,R/Γ_ph ≫ 1, we find

$$\begin{equation} g^{(2)}_{\rm S} (0) = 2 + \mathcal{O}(\Gamma_{\rm ph} / \Gamma_{L, R}), \end{equation} \tag{ 33 }$$

which indicates the thermalization induced by the Franck–Condon effect.

Next, we examine the Franck–Condon blockade under finite bias voltages. Figures 7(a) and (b) show the current as a function of V when Δ = 0. The dimensionless coupling constant is (a) λ_S = 0.1 and (b) 1. While the Franck–Condon blockade suppresses the current under small bias voltages in the case of λ_S = 1, the current suppression is negligible in the case of λ_S = 0.1. In the former case, a larger voltage is needed to lift the Franck–Condon blockade for larger Γ_L,R/Γ_ph.

Zoom In Zoom Out Reset image size

When Γ_L,R/Γ_ph = 100, the V dependence of the current is almost the same as in figure 5 for A-phonons with Δ = 0. The phonon number and its autocorrelation function also change with the bias voltage V in a similar manner to those at the current main peak for A-phonons.

First, we calculate the current I = eΓ_L〈n₀〉. For the purpose, we sum up both sides of equation (C.1) over n. Using

$$\begin{eqnarray*} &&\sum_{n = 0}^\infty |_{\rm A} \langle n | \mathcal{T}_{\rm A}^\dagger | m \rangle_{\rm A} |^2 = \left._{\rm A} \langle m | \mathcal{T}_{\rm A} \left( \sum_n |n \rangle_{\rm AA} \langle n| \right) \mathcal{T}_{\rm A}^\dagger | m \rangle_{\rm A} \right. = 1, \end{eqnarray*}$$

we obtain

$$\begin{eqnarray*} &&0 = -\Gamma_L \langle n_0 \rangle + \frac{\Gamma_R}{2} \langle n_{\rm mol} \rangle + \Gamma_R \langle \tilde n_R \rangle. \end{eqnarray*}$$

Since $\langle n_0 \rangle + \langle n_{\mathrm { mol}} \rangle + \langle \tilde n_R \rangle =1$, we obtain

$$\begin{eqnarray*} &&\langle n_0 \rangle = \frac{\gamma}{2 + \gamma} (1 + \langle \tilde n_R \rangle), \quad \langle n_{\rm mol} \rangle = \frac{2}{2 + \gamma} \left[ 1 - (1 + \gamma) \langle \tilde n_R \rangle \right], \end{eqnarray*}$$

where γ = Γ_R/Γ_L. These equations result in equation (21), i.e.

$$\begin{eqnarray*} &&I = e \Gamma_R \frac{1 + \langle \tilde n_R \rangle}{2 + \gamma}, \quad \langle n_{\rm e} \rangle = \frac{2 - \gamma \langle \tilde n_R \rangle} {2 + \gamma}. \end{eqnarray*}$$

The summation of equation (C.3) over n yields

$$\begin{equation} \langle \tilde n_R \rangle = \frac{\nu \Gamma_{\rm ph}}{2 \Gamma_R} P_{{\rm mol}, 0}. \end{equation} \tag{ C.4 }$$

Next, we derive the phonon number

$$\begin{eqnarray*}\fl \langle N_{\rm A} \rangle\! &=&\!\! \sum_{n = 0}^\infty \! P_{0, n~{\rm eA}}\langle 0, n| N_{\rm A} |0, n \rangle_{\rm eA}\! +\! \sum_{\sigma = \pm} \sum_{n = 0}^\infty\! P_{\sigma, n~{\rm eA}}\langle \sigma, n | N_{\rm A} |\sigma, n \rangle_{\rm eA}\! +\! \sum_{n = 0}^{\nu - 1}\! P_{R, n~{\rm eA}}\langle R, n| N_{\rm A} |R, n \rangle_{\rm eA} \\ \fl &=& \sum_{n = 0}^\infty n P_{0, n} + \sum_{n = 0}^\infty \left( n + \frac{\nu}{2} + \lambda_{\rm A}^2 \right) P_{{\rm mol}, n} + \sum_{n = 0}^{\nu - 1} (n + \lambda_{\rm A}^2) P_{R, n} \\ \fl &\equiv &\langle N_{\rm A} n_0 \rangle + \langle N_{\rm A} n_{\rm mol} \rangle + \langle N_{\rm A} \tilde n_R \rangle. \end{eqnarray*}$$

We have used the relation, $\mathcal {T}_{\mathrm { A}}^\dagger N_{\mathrm { A}} \mathcal {T}_{\mathrm { A}} = ( \mathcal {T}_{\mathrm { A}}^\dagger a_{\mathrm { A}}^\dagger \mathcal {T}_{\mathrm { A}} ) ( \mathcal {T}_{\mathrm { A}}^\dagger a_{\mathrm { A}} \mathcal {T}_{\mathrm { A}} ) = (a_{\mathrm { A}}^\dagger - \lambda _{\mathrm { A}}) (a_{\mathrm { A}} - \lambda _{\mathrm { A}})$. We multiply both sides of equations (C.1)–(C.3) by n and sum up over n. Then we find

$$\begin{equation}\fl 0 = - (\Gamma_L + \Gamma_{\rm ph}) \langle N_{\rm A} n_0 \rangle + \frac{\Gamma_R}{2} \langle N_{\rm A} n_{\rm mol} \rangle + \Gamma_R \langle N_{\rm A} \tilde n_R \rangle + \frac{\nu \Gamma_R}{4} \langle n_{\rm mol} \rangle, \end{equation} \tag{ C.5 }$$

$$\begin{eqnarray}&&\fl 0 = \Gamma_L \langle N_{\rm A} n_0 \rangle - \left( \frac{\Gamma_R}{2} +\Gamma_{\rm ph} \right) \langle N_{\rm A} n_{\rm mol} \rangle + \lambda_{\rm A}^2 \Gamma_L \langle n_0 \rangle + \left( \frac{\nu + 2 \lambda_{\rm A}^2}{4} \Gamma_R + \lambda_{\rm A}^2 \Gamma_{\rm ph} \right) \langle n_{\rm mol} \rangle + \Gamma_R \langle \tilde n_R \rangle,\nonumber\\ \end{eqnarray} \tag{ C.6 }$$

$$\begin{equation}\fl 0 = - (\Gamma_R + \Gamma_{\rm ph}) \langle N_{\rm A} \tilde n_R \rangle + \left[ (\nu - 1 + \lambda_{\rm A}^2) \Gamma_R + \lambda_{\rm A}^2 \Gamma_{\rm ph} \right] \langle \tilde n_R \rangle. \end{equation} \tag{ C.7 }$$

Here, we have used

$$\begin{equation}\fl \sum_{n = 0}^\infty n |_{\rm A} \langle n | \mathcal{T}_{\rm A}^\dagger | m \rangle_{\rm A} |^2 =_{\rm A} \langle m | \mathcal{T}_{\rm A} N_{\rm A} \left( \sum_n |n \rangle_{\rm AA} \langle n| \right) \mathcal{T}_{\rm A}^\dagger | m \rangle_{\rm A}=\sum_{m}{_{\rm A} \langle} m | \mathcal{T}_{\rm A} N_{\rm A} \mathcal{T}_{\rm A}^\dagger | m \rangle_{\rm A}. \end{equation} \tag{ C.8 }$$

From equations (C.5)–(C.7), we obtain equation (22), i.e.

$$\begin{eqnarray*} &&\langle N_{\rm A} \rangle = (\nu + 2 \lambda_{\rm A}^2) \frac{I}{e \Gamma_{\rm ph}} + \lambda_{\rm A}^2 \langle n_{\rm e} \rangle. \end{eqnarray*}$$

Finally, we derive the phonon autocorrelation function at τ = 0,

$$\begin{eqnarray*} &&g^{(2)}_{\rm A} (0) = \frac{\langle :N_{\rm A}^2: \rangle}{\langle N_{\rm A} \rangle^2} = \frac{\langle N_{\rm A}^2 \rangle - \langle N_{\rm A} \rangle}{\langle N_{\rm A} \rangle^2}, \end{eqnarray*}$$

where

$$\begin{eqnarray*}&&\fl \langle N_{\rm A}^2 \rangle\! =\! \sum_{n = 0}^\infty \!P_{0, n~{\rm eA}}\!\langle 0, n| N_{\rm A}^2 |0, n \rangle_{\rm eA} +\! \sum_{\sigma = \pm} \sum_{n = 0}^\infty\! P_{\sigma, n~{\rm eA}} \langle \sigma, n | N_{\rm A}^2 |\sigma, n \rangle_{\rm eA}\! +\! \sum_{n = 0}^{\nu - 1}\! P_{R, n~{\rm eA}}\langle R, n| N_{\rm A}^2 |R, n \rangle_{\rm eA} \\ &&= \sum_{n = 0}^\infty n^2 P_{0, n} + \sum_{n = 0}^\infty \left[ n^2 + \lambda_{\rm A}^2 (4n + 1 + \lambda_{\rm A}^2) + \nu \left( n + \frac{\nu}{2} + 2 \lambda_{\rm A}^2 \right) \right] P_{{\rm mol}, n} \\ &&\quad + \sum_{n = 0}^{\nu - 1} \left[ n^2 + \lambda_{\rm A}^2 (4 n + 1 + \lambda_{\rm A}^2) \right] P_{R, n} \\ &&\equiv \langle N_{\rm A}^2 n_0 \rangle + \langle N_{\rm A}^2 n_{\rm mol} \rangle + \langle N_{\rm A}^2 \tilde n_R \rangle.\ \end{eqnarray*}$$

We multiply both sides of equations (C.1)–(C.3) by n² and sum up over n. A similar technique to equation (C.8) leads to

$$\begin{eqnarray}&&\fl 0 = -(\Gamma_L + 2 \Gamma_{\rm ph}) \langle N_{\rm A}^2 n_0 \rangle + \frac{\Gamma_R}{2} \langle N_{\rm A}^2 n_{\rm mol} \rangle + \Gamma_R \langle N_{\rm A}^2 \tilde n_R \rangle + \Gamma_{\rm ph} \langle N_{\rm A} n_0 \rangle \nonumber\\ &&+ \frac{\nu \Gamma_R}{2} \langle N_{\rm A} n_{\rm mol} \rangle + \frac{\nu \lambda_{\rm A}^2 \Gamma_R}{2} \langle n_{\rm mol} \rangle, \end{eqnarray} \tag{ C.9 }$$

$$\begin{eqnarray}&&\fl 0 = \Gamma_L \langle N_{\rm A}^2 n_0 \rangle - \left( \frac{\Gamma_R}{2} + 2 \Gamma_{\rm ph} \right) \langle N_{\rm A}^2 n_{\rm mol} \rangle + 4 \lambda_{\rm A}^2 \Gamma_L \langle N_{\rm A} n_0 \rangle \nonumber \\ &&+ \left[ \frac{\nu + 4 \lambda_{\rm A}^2}{2} \Gamma_R + ( \nu + 1 + 8 \lambda_{\rm A}^2) \Gamma_{\rm ph} \right] \langle N_{\rm A} n_{\rm mol} \rangle \lambda_{\rm A}^2 (1 + \lambda_{\rm A}^2) \Gamma_L \langle n_0 \rangle \nonumber \\ &&+ \lambda_{\rm A}^2 \left[ \frac{1 - \nu - 3 \lambda_{\rm A}^2}{2} \Gamma_R + (1 - \nu - 6 \lambda_{\rm A}^2) \Gamma_{\rm ph} \right] \langle n_{\rm mol} \rangle - \Gamma_R \langle \tilde n_R \rangle, \end{eqnarray} \tag{ C.10 }$$

$$\begin{eqnarray}&&\fl 0 = -(\Gamma_R + 2 \Gamma_{\rm ph}) \langle N_{\rm A}^2 \tilde n_R \rangle + \left[ 4 \lambda_{\rm A}^2 \Gamma_R + (1 + 8 \lambda_{\rm A}^2) \Gamma_{\rm ph} \right] \langle N_{\rm A} \tilde n_R \rangle \nonumber \\ &&+ \left\{ \left[ (\nu - 1)^2 + \lambda_{\rm A}^2 (1 - 3\lambda_{\rm A}^2) \right] \Gamma_R + \lambda_{\rm A}^2 ( 1 - 6 \lambda_{\rm A}^2) \Gamma_{\rm ph} \right\} \langle \tilde n_R \rangle . \end{eqnarray} \tag{ C.11 }$$

From equations (C.9)–(C.11), we find

$$\begin{eqnarray*}&&\fl \langle N_{\rm A}^2 \rangle - \langle N_{\rm A} \rangle = 2 \lambda_{\rm A}^2 \frac{\Gamma_L}{\Gamma_{\rm ph}} \langle N_{\rm A} n_0 \rangle + \left( \frac{\nu + 2 \lambda_{\rm A}^2}{2} \frac{\Gamma_R}{\Gamma_{\rm ph}} + \frac{\nu + 4 \lambda_{\rm A}^2}{2} \right) \langle N_{\rm A} n_{\rm mol} \rangle+ 2 \lambda_{\rm A}^2 \left( \frac{\Gamma_R}{\Gamma_{\rm ph}} + 2 \right) \langle N_{\rm A} \tilde n_R \rangle\\ &&- \frac{\nu + 2 \lambda_{\rm A}^4}{2} \frac{\Gamma_L}{\Gamma_{\rm ph}} \langle n_0 \rangle - \frac{\lambda_{\rm A}^2 (\nu + 6 \lambda_{\rm A}^2)}{2} \langle n_{\rm mol} \rangle - \left[ \frac{\nu (2 - \nu)}{2} \frac{\Gamma_R}{\Gamma_{\rm ph}} + 3 \lambda_{\rm A}^2 \right] \langle \tilde n_R \rangle. \end{eqnarray*}$$

Now we evaluate g⁽²⁾_A(0) when Γ_L,R ≫ Γ_ph. In this case, $\langle \tilde n_R \rangle = \mathcal {O} (\Gamma _{\mathrm { ph}} / \Gamma _{L, R})$ from equation (C.4). Then

$$\begin{eqnarray*} &&I = \frac{e \Gamma_R}{2 + \gamma} + \mathcal{O}(\Gamma_{\rm ph} / \Gamma_{L, R}), \quad \langle N_{\rm A} \rangle = \frac{\nu + 2 \lambda_{\rm A}^2}{2 + \gamma} \frac{\Gamma_R}{\Gamma_{\rm ph}} + \mathcal{O}(1). \end{eqnarray*}$$

Equations (C.5) and (C.7) yield

$$\begin{eqnarray*} &&2 \langle N_{\rm A} n_0 \rangle = \gamma \langle N_{\rm A} n_{\rm mol} \rangle + \mathcal{O} (1), \quad \langle N_{\rm A} \tilde n_R \rangle = \mathcal{O}(\Gamma_{\rm ph} / \Gamma_{L, R}). \end{eqnarray*}$$

Using $\langle N_{\mathrm { A}} \rangle = \langle N_{\mathrm { A}} n_0 \rangle + \langle N_{\mathrm { A}} n_{\mathrm { mol}} \rangle + \langle N_{\mathrm { A}} \tilde n_R \rangle$, we have

$$\begin{eqnarray*} \hphantom{0\,}\langle N_{\rm A} n_0 \rangle &=& (\nu + 2 \lambda_{\rm A}^2) \frac{\gamma}{(2 + \gamma)^2} \frac{\Gamma_R}{\Gamma_{\rm ph}} + \mathcal{O} (1), \\ \langle N_{\rm A} n_{\rm mol} \rangle &=& (\nu + 2 \lambda_{\rm A}^2) \frac{2}{(2 + \gamma)^2} \frac{\Gamma_R}{\Gamma_{\rm ph}} + \mathcal{O}(1). \end{eqnarray*}$$

Using these relations, we obtain equation (24), i.e.

$$\begin{eqnarray*} &&g^{(2)}_{\rm A} (0) = \frac{\nu + 4 \lambda_{\rm A}^2}{\nu + 2 \lambda_{\rm A}} + \mathcal{O} (\Gamma_{\rm ph} / \Gamma_{L, R}). \end{eqnarray*}$$

Here, we comment on the opposite limit of Γ_L,R ≪ Γ_ph. The analytical expressions for the current, g⁽²⁾_A(0), etc can be obtained in a similar way to the case of Γ_L,R ≫ Γ_ph. At the main peak of the current (Δ = 0), the electron number and current are written as

$$\begin{eqnarray*} &&\langle n_{\rm e} \rangle = \langle n_{\rm mol} \rangle \simeq \frac{1}{2 + \gamma}, \quad I \simeq \frac{e \Gamma_R}{2 + \gamma}, \end{eqnarray*}$$

whereas the phonon number and its autocorrelation function are

$$\begin{eqnarray*} &&\langle N_{\rm A} \rangle \simeq \frac{\lambda_{\rm A}^2}{2 + \gamma}, \quad g^{(2)}_{\rm A}(0) \simeq 2 + \gamma. \end{eqnarray*}$$

At the νth subpeak of the current (Δ = Δ_ν ≃ νℏω_ph), we obtain

$$\begin{eqnarray*} &&\langle n_{\rm e} \rangle \simeq \langle \tilde n_R \rangle \simeq \frac{1}{1 + \gamma}, \quad \langle n_{\rm mol} \rangle \simeq 0, \quad I \simeq \frac{e \Gamma_R}{1 + \gamma} \end{eqnarray*}$$

and

$$\begin{equation} \langle N_{\rm A} \rangle \simeq \frac{\lambda_{\rm A}^2}{1 + \gamma}, \quad g^{(2)}_{\rm A}(0) \simeq 1 + \gamma. \end{equation} \tag{ C.12 }$$

As discussed in section 3.4, phonon bunching is observed even for small λ_A in the case of Γ_L,R ≪ Γ_ph. In this case, an electron is localized in dot R for a long time, forming a polaron |R,0〉_eA, after a phonon is immediately decayed. Thus, the number of phonons created by the interdot tunneling is much smaller than that accompanied by the polaron staying in dot R. This situation results in the bunched phonons.