Nanomechanics driven by the superconducting proximity effect

Nanoelectromechanical systems (NEMS) provide a promising platform for investigations into the quantum mechanical interplay between mechanical and electronic subsystems [1, 2]. One of the most important phenomena providing the foundation of NEMS functionality is the generation of self-sustained mechanical oscillations by a dc flow [3, 4, 5, 7, 8, 6]. This effect is itself an interesting problem from a fundamental point of view, opening new possibilities for mass and force sensing [9, 10], while its underlying physical processes show potential applications for mechanical cooling [11]. Self-sustained mechanical oscillations were first observed in a carbon nanotube (CNT)-based transistor [12], with further studies later verifying their transport signatures [13, 14, 15]. Recently, the experimental observation of self-driven oscillations of a CNT-based quantum dot in the Coulomb blockade regime has been reported [16].

Superconducting elements incorporated into NEMS extend the horizon of this phenomenon, namely through the effects of superconducting phase coherence; see, for example, the following reviews [17, 18]. Interplay between electromechanical effects and phase coherence gives new and unusual properties to a number of normal metal/superconducting hybrid junctions [19, 20, 21]. In particular, it has recently been shown that in a normal metal–suspended CNT–superconductor transistor, Andreev reflection [22, 23] may give rise to a cooling of the mechanical subsystem [24] or generate a single-atom lasing effect [25] if certain conditions are fulfilled.

The mechanical functionality of NEMS is to a large extent determined by the physical principles underlying the interaction between the electronic and mechanical subsystems. In all studies mentioned above, this interaction was due to the localization of the charge [4, 5, 6] or spin [7, 8] carried by electrons in the movable part of the system. In the present paper, we consider a fundamentally new type of electromechanical coupling based on the quantum delocalization of Cooper pairs. We demonstrate that such coupling can promote a self-saturated mechanical instability resulting in the generation of self-sustained mechanical oscillations. It is also shown that these oscillations significantly increase the average current through the system, making it possible to carry out direct experimental detection.

A sketch of the NEMS investigated in this paper is presented in Fig. 1. A single-walled CNT is suspended above a trench in a bulk normal metal electrode biased by a constant voltage $V_{b}$. Two superconducting leads with the superconducting phase difference $\phi$ are positioned near the middle of the nanotube in such a way that the bending of the nanotube moves it closer to one electrode and further away from the other. The distance between the quantized electronic levels inside the nanotube is much greater than the other energy parameters, allowing one to consider the nanotube as a single-level quantum dot (QD). The bending dynamics of the CNT are reduced to the dynamics of the fundamental flexural mode. We suppose that the amplitude of this mode, $x$, is larger than the amplitude of zero-point oscillations, and will consider it as a classical mechanical oscillator with mass $m$ and frequency $\omega$. The dynamics is described by Newton’s equation,

$$\ddot{x}+\omega^{2}x=-\frac{1}{m}\text{Tr}\left\{\hat{\rho}\frac{\partial H(x)}{\partial x}\right\},$$

(1)

where

$$H=H_{d}+H_{l}+H_{t}$$

(2)

is the Hamiltonian of the electronic subsystem. The first term $H_{d}$ represents the single-level QD,

$$H_{d}=\sum_{\sigma}\varepsilon_{d}d^{\dagger}_{\sigma}d_{\sigma},$$

(3)

where $d_{\sigma}^{\dagger}(d_{\sigma})$ is the creation (annihilation) operator of an electron with spin projection $\sigma=\uparrow,\downarrow$ on the dot. The Hamiltonian $H_{l}=H_{l}^{n}+H_{l}^{s}$ describes the normal and superconducting leads, respectively, with

	$\displaystyle H_{l}^{n}=\sum_{k\sigma}(\varepsilon_{k}-eV_{b})a_{k\sigma}^{\dagger}a_{k\sigma},$		(4)
	$\displaystyle H_{l}^{s}=\sum_{kj\sigma}\left(\varepsilon_{k}c^{\dagger}_{kj\sigma}c_{kj\sigma}-\Delta_{s}(\text{e}^{\imath\phi_{j}}c^{\dagger}_{kj\uparrow}c^{\dagger}_{-kj\downarrow}+\text{H.c.})\right).$		(5)

Here, $a^{\dagger}_{k\sigma}(a_{k\sigma})$, and $c^{\dagger}_{kj\sigma}(c_{k\sigma})$ are the creation (annihilation) operators of an electron with quantum number $k$ and spin projection $\sigma$ in the normal and superconducting $j=1,2$ leads, respectively, and $\Delta_{s}\text{e}^{\imath\phi_{j}}$ is the superconducting order parameter (in the $j$ electrode). Note that the energies $\varepsilon_{d},\varepsilon_{k}$ are counted from the Fermi energy of the superconductors. In what follows, we set $\phi_{1}=-\phi_{2}=\phi/2$.

The Hamiltonian $H_{t}=H_{t}^{n}+H_{t}^{s}$ describes the tunneling of electrons between the dot and the leads, where

	$\displaystyle H_{t}^{n}=\sum_{k\sigma}t^{n}_{0}(a_{k\sigma}^{\dagger}d_{\sigma}+\text{H.c.}),$		(6)
	$\displaystyle H_{t}^{s}=\sum_{kj\sigma}t^{s}_{j}(x)(c^{\dagger}_{kj\sigma}d_{\sigma}+\text{H.c.}).$		(7)

The position-dependent superconducting tunneling amplitude $t^{s}_{1(2)}(x)=t^{s}_{0}\text{e}^{(-1)^{j}(x+a)/2\lambda}$, where $2\lambda$ is the characteristic tunneling length and $a$ is a parameter for asymmetry. For a typical CNT-based nanomechanical resonator, $2\lambda\sim 0.5$ nm [26]. We concentrate our attention on the symmetric case $a=0$ and leave the asymmetric one for discussion elsewhere.

The time evolution of the electronic density matrix $\hat{\rho}$ is described by the Liouville–von Neumann equation ($\hbar=1$),

$$\imath\partial_{t}\hat{\rho}=[H,\hat{\rho}],$$

(8)

which together with Eq. (1) forms a closed system of equations that describe the nanoelectromechanics of our system. In this paper, we restrict ourselves to the case $\Delta_{s}\gg|eV_{b}|\gg\Delta_{d}\sim\Gamma_{n}$, where $\Delta_{d}=(2\pi)\nu_{s}|t^{s}_{0}|^{2}$ and $\Gamma_{n}=(2\pi)\nu_{n}|t^{n}_{0}|^{2}$, with $\nu_{s(n)}$ the density of states in the superconducting (normal) electrode.

To describe the electronic dynamics of the QD, we use the reduced density matrix approximation in which the full density matrix of the system is factorized to the tensor product of the equilibrium density matrices of the normal and superconducting leads and the density matrix of the dot as $\hat{\rho}=\hat{\rho}_{n}\otimes\hat{\rho}_{s}\otimes\hat{\rho}_{d}$. Using the standard procedure, one can trace out the degrees of freedom of the leads and obtain the following equation for the reduced density matrix $\hat{\rho}_{d}$ [6] (in the deep subgap regime $\Delta_{s}\rightarrow\infty$),

$$\partial_{t}\hat{\rho}_{d}=-\imath\left[H_{d}^{eff},\hat{\rho}_{d}\right]+\mathcal{L}_{n}\{\hat{\rho}_{d}\},$$

(9)

where

$$H^{eff}_{d}=H_{d}+\Delta_{d}(x,\phi)d_{\downarrow}d_{\uparrow}+\Delta_{d}^{*}(x,\phi)d_{\uparrow}^{\dagger}d_{\downarrow}^{\dagger},$$

(10)

$$\Delta_{d}(x,\phi)=\Delta^{\prime}(x,\phi)+i\Delta^{\prime\prime}(x,\phi)=\Delta_{d}\cosh(x/\lambda+i\phi/2).$$

Above, $\Delta_{d}(x,\phi)$ is the off-diagonal order parameter induced by the superconducting proximity effect [27, 24], and $\Delta^{\prime,\prime\prime}(x,\phi)$ are real functions. The Lindbladian term in Eq. (9) reflects the incoherent electron exchange between the normal lead and QD. The latter in the high bias voltage regime, $|eV_{b}|\gg\varepsilon_{0},k_{B}T$, takes the form

$$\mathcal{L}_{n}\{\hat{\rho}_{d}\}=\Gamma_{n}\sum\limits_{\sigma}\begin{cases}2d^{\dagger}_{\sigma}\hat{\rho}_{d}d_{\sigma}-\left\{d_{\sigma}d_{\sigma}^{\dagger},\hat{\rho}_{d}\right\},&\kappa=+1;\\ 2d_{\sigma}\hat{\rho}_{d}d^{\dagger}_{\sigma}-\left\{d^{\dagger}_{\sigma}d_{\sigma},\hat{\rho}_{d}\right\},&\kappa=-1;\end{cases}$$

(11)

where $\kappa=\text{sgn}(eV_{b})$.

Figure 2 represents the electronic dynamics on the dot for $\kappa=\pm 1$. From Fig. 2, one can see that not all electron processes are allowed due to the parameter scales in this work. In the subgap regime, single-electron transitions between the dot and the superconducting leads are prohibited, and thus only an exchange of Cooper pairs occurs. Moreover, because of the high bias voltage, single-electron tunneling between the dot and the normal leads is enabled exclusively in one direction (from the lead to the dot, see Fig. 2(a), or vice-versa, Fig. 2(b)), establishing that our model is electron-hole symmetric.

As a consequence, the QD density matrix $\hat{\rho}_{d}$ acts in the Hilbert space $\mathcal{H}_{4}$, which may be presented as a direct sum of two $\mathcal{H}_{2}$ spaces via $\mathcal{H}_{4}=\mathcal{H}_{e}\oplus\mathcal{H}_{CP}$ spanned over state vectors $|\uparrow\rangle=d^{\dagger}_{\uparrow}|0\rangle$, $|\downarrow\rangle=d^{\dagger}_{\downarrow}|0\rangle$, and $|0\rangle$, $|2\rangle=d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}|0\rangle$ (with $d_{\uparrow,\downarrow}|0\rangle=0$).

The superselection rule, which forbids the superposition of states with integer and half-integer spins, allows us to present the density matrix $\hat{\rho}_{d}$ as a direct sum of two density matrices $\hat{\rho}_{d}=\hat{\rho}_{e}\oplus\hat{\rho}_{CP}$ acting in the $\mathcal{H}_{2}$ Hilbert space spanned over state vectors $|\uparrow\rangle,|\downarrow\rangle$ and $|0\rangle,|2\rangle$, respectively. Moreover, taking into account spin-rotation symmetry, one can conclude that $\hat{\rho}_{e}$ should be proportional to the unit matrix, $\hat{\rho}_{e}=\rho_{e}\hat{I}$, while $\hat{\rho}_{CP}$ can be written in the form $\hat{\rho}_{CP}=\frac{1}{2}R_{0}\hat{I}+\frac{1}{2}\sum_{i}R_{i}\sigma_{i}$, where $\sigma_{i},(i=1,2,3)$ are the Pauli matrices.

Then by introducing the dimensionless time $\omega t\rightarrow t$ and displacement $x/\lambda\rightarrow x$, and taking into account the normalization condition $\text{Tr}\hat{\rho}_{d}=1$, we get the following closed system of equations for $x(t)$ and $R_{i}(t)$,

$$\ddot{x}+x=-\xi\left[\sinh{(x)}\cos{\left(\frac{\phi}{2}\right)}R_{1}-\cosh{(x)}\sin{\left(\frac{\phi}{2}\right)}R_{2}\right],$$

(12)

$$\alpha\dot{\vec{R}}=\hat{L}\vec{R}-\kappa\vec{e}_{3},$$

(13)

where $\vec{R}=(R_{1},R_{2},R_{3})^{T}$, $\vec{e_{3}}=(0,0,1)^{T}$, $\xi=\Delta_{d}/(m\lambda^{2}\omega^{2})$ is the nanoelectromechanical coupling parameter, and $\alpha=\omega/(2\Gamma_{n})$ is the adiabaticity parameter. For a typical CNT-based NEMS, one can estimate $\xi\sim 10^{-3}\ll 1$ [26, 28]. The matrix $\hat{L}$ is defined as follows,

$$\hat{L}(x)=\begin{pmatrix}-1&\tilde{\varepsilon}_{d}&-\tilde{\Delta}^{\prime\prime}(x,\phi)\\ -\tilde{\varepsilon}_{d}&-1&-\tilde{\Delta}^{\prime}(x,\phi)\\ \tilde{\Delta}^{\prime\prime}(x,\phi)&\tilde{\Delta}^{\prime}(x,\phi)&-1\end{pmatrix},$$

(14)

where $\tilde{\varepsilon}_{d}\equiv\varepsilon_{d}/\Gamma_{n},\tilde{\Delta}_{d}\equiv\Delta_{d}/\Gamma_{n}$.

The system of Eqs. (12) and (13) has an obvious static solution $x_{st}=0+\mathcal{O}(\xi)$, $\vec{R}_{st}=\kappa L^{-1}(0)\vec{e}_{3}+\mathcal{O}(\xi)\vec{R}^{(1)}$, here $\|\vec{R}^{(1)}\|=1$. The stability of this solution can then be investigated in a standard way, see for example Ref. [29]. However, to simplify this procedure, we will consider the adiabatic case when $\alpha\ll 1$, which corresponds to a typical experimental situation [16] and reduces the problem to one that allows the use of Poincare analysis. More specifically, this inequality allows one to find a solution of Eq. (13) to the accuracy $\alpha$,

$$\vec{R}(x,t)=\kappa L^{-1}(x(t))(1+\alpha\dot{x}\partial_{x}L^{-1}(x(t))+\mathcal{O}(\alpha^{2}))\vec{e}_{3},$$

(15)

and then substituting this solution into Eq. (12) gives (to accuracy $\alpha$) the following nonlinear differential equation for $x(t)$,

$$\ddot{x}-\eta(x,\phi)\dot{x}+x=F(x,\phi),$$

(16)

the solution of which may be analyzed via Poincare’s theory. Here, the nonlinear force $F(x,\phi)$ and friction coefficient $\eta(x,\phi)$, which in what follows we refer to as a pumping coefficient, generated by interaction with the nonequilibrium electronic environment take the form,

	$\displaystyle F(x,\phi)=\kappa\xi\frac{\tilde{\Delta}_{d}}{2\tilde{\mathcal{D}}^{2}}\left[\sin{\phi}-\tilde{\varepsilon}_{d}\sinh{\left(2x\right)}\right],$		(17)
	$\displaystyle\eta(x,\phi)=\kappa\alpha\xi\frac{\tilde{\Delta}_{d}}{2\tilde{\mathcal{D}}^{6}}\left[\tilde{\Delta}_{d}^{2}\sinh{(2x)}\left(\tilde{\mathcal{D}}^{2}+4\right)\left\{\sin{\phi}-\tilde{\varepsilon}_{d}\sinh{(2x)}\right\}+8\tilde{\varepsilon}_{d}\tilde{\mathcal{D}}^{2}\left\{\sin^{2}{(\phi/2)}+\sinh^{2}{x}\right\}\right].$		(18)

Here $\tilde{\mathcal{D}}^{2}\equiv\mathcal{D}^{2}(x,\phi)/\Gamma_{n}^{2}=\tilde{\Delta}_{d}^{2}\left[\sinh^{2}{x}+\cos^{2}{(\phi/2)}\right]+\tilde{\varepsilon}^{2}_{d}+1$.

In order to find the stationary solutions $x_{c}(t)$ in Eq. (16), it is natural to use the smallness of the parameter $\xi$ and look for such solutions [30] in the form $x_{c}(t)=x_{st}+\sqrt{A}\sin\left(t+\varphi(t)\right)+\mathcal{O}(\xi)$, where $x_{st},\dot{A}(t),\dot{\varphi}(t)\sim\xi$. Then, with the accepted accuracy $\xi$, one can get the following equations:

	$\displaystyle\dot{A}=A\bar{\eta}(A,\phi),$		(19)
	$\displaystyle\dot{\varphi}=-A^{-1/2}\bar{F}(A,\phi),$		(20)

where

	$\displaystyle\bar{\eta}(A,\phi)=\kappa(\pi)^{-1}\int_{0}^{2\pi}d\psi\cos^{2}{(\psi)}\eta\left(\sqrt{A}\sin\psi,\phi\right)\equiv\kappa\xi\alpha W(A,\phi),$		(21)
	$\displaystyle\bar{F}(A,\phi)=\kappa(2\pi)^{-1}\int_{0}^{2\pi}d\psi\sin{(\psi)}F\left(\sqrt{A}\sin\psi,\phi\right).$		(22)

The pumping coefficient $\bar{\eta}(A,\phi)$ has an obvious physical meaning: it gives the ratio between the energy supplied into the mechanical degree of freedom for one period of mechanical oscillation with amplitude $\sqrt{A}$ and the total mechanical energy.

It is evident from Eq. (19) that stationary regimes $\dot{A}=0$ are given by equations $A=0$ ($x(t)=x_{st}$) and $\bar{\eta}(A,\phi)=0$. The first one is a static state of the nanotube, and the second one corresponds to periodic oscillations with the amplitude $\sqrt{A_{c}}$, where $W(A_{c},\phi)=0$. The static regime is stable when $\bar{\eta}(0,\phi)<0$ and unstable otherwise. The stability of the periodic solution is defined by the sign of the derivative $\partial_{A}{\bar{\eta}(A,\phi)}|_{A=A_{c}}$: if it is negative (positive), then the periodic regime is stable (unstable). Analyzing Eqs. (16) and (21), one can conclude that the pumping coefficient $\bar{\eta}(A,\phi)\propto\kappa$ is an odd function of $\varepsilon_{d}$ and takes the following limit values,

	$\displaystyle\bar{\eta}(0,\phi)=\kappa\alpha\xi W(0,\phi)=+\kappa\alpha\xi\frac{4\tilde{\varepsilon}_{d}\tilde{\Delta}_{d}}{\tilde{\mathcal{D}}^{4}(0,\phi)}\sin^{2}{\left(\frac{\phi}{2}\right)},$		(23)
	$\displaystyle\eta(A\rightarrow\infty,\phi)=\kappa\alpha\xi W(A\rightarrow\infty,\phi)=-\kappa\alpha\xi\frac{\tilde{\varepsilon}_{d}}{2\tilde{\Delta}_{d}},$		(24)

from which follows that at $\phi\neq 0$, the solution of the equation $W(A_{c},\phi)=0$, corresponding to the stationary periodic regime, exists at any values of the other parameters. The case when $\phi=0$ is very unstable with respect to the small asymmetry parameter $|a|\ll 1$ [see below Eq. (7)] and will be analyzed elsewhere. The function $W(A)$ and $A_{c}$ at different $\phi\geq 1$ and $\tilde{\varepsilon}_{d}>0$ are presented in Figs. 3 and 4.

From Eq. (23), it follows that if $\kappa\varepsilon_{d}>0$ and $\phi\neq 0$, the static mechanical state $x=x_{st}\ll 1$ is unstable with respect to the appearance of bending oscillations with amplitudes that exponentially increase in time with the increment $\gamma=\kappa\alpha\xi W(0,\phi)$. The latter takes its maximum at $\phi=\pi$ for the fixed values of all other parameters (notice that $x_{st}(\pi)=0$). However, the increase saturates at the amplitude $\sqrt{A_{c}}$, resulting in self-sustained oscillations at this amplitude. It should be noted that the amplitude saturation in the system under consideration is a completely internal effect and still takes place when the friction caused by interaction with a thermodynamic environment is zero. A ”self-saturation” effect was also reported in [8] where a special magnetic NEM system was considered.

The self-sustained oscillations considered above have a very specific transport signature. This raises the possibility of detecting the mechanical instability through electric current measurement. To explore such a possibility, let us consider the electric current through the system $I_{n}$, determined in a standard way, $I_{n}=e\kappa\text{Tr}\left\{\hat{\dot{N}}\hat{\rho}\right\}$, where $\hat{\dot{N}}=i[\hat{H},\hat{N}]$ and $\hat{N}$ is the operator of the number of electrons in the normal electrode. In the high bias voltage regime at $\kappa=+1$, where electron tunneling from the QD to the normal leads is forbidden, an expression for $I_{n}$ can be easily obtained by analyzing Fig. 2. From those diagrams, one can see that a decrease in the number of electrons in the normal electrode is defined by two different processes. The first one is the tunneling of an electron with spin up or down into the empty dot. The rate of this process is $2\Gamma_{n}\rho_{0}$, where $\rho_{0}=\left(R_{0}+R_{3}\right)/2$ is the probability that the dot is empty. The second one is the tunneling of an electron into the dot occupied by a single electron with spin up or down. The rate of this process is $2\Gamma_{n}\rho_{e}$. Taking into account the normalization condition $2\rho_{e}+R_{0}=1$, and using a similar speculation for $eV_{b}<0$, one gets from Eq. (15) the following equation for $I_{n}$,

$$I_{n}(t)=\kappa I_{0}\left(1+\kappa R_{3}\right)=\kappa I_{0}\left[\frac{|\Delta_{d}(x,\phi)|^{2}}{\mathcal{D}^{2}(x,\phi)}+\alpha\dot{x}f(x)+\mathcal{O}(\alpha^{2})\right],$$

(25)

where $I_{0}=e\Gamma_{n}$. In the stationary regime corresponding to the generation of self-sustained oscillations with amplitude $\sqrt{A_{c}}$, the average electric current $\bar{I}_{n}$ is defined by the equation

$${\bar{I}}_{n}(\kappa,\varepsilon_{d})=\kappa I_{0}\left[\frac{\Delta_{d}^{2}\cos^{2}{(\phi/2})}{\Delta_{d}^{2}\cos^{2}(\phi/2)+\Gamma_{n}^{2}+\varepsilon_{d}^{2}}+\theta(\kappa\varepsilon_{d})\delta\bar{I}(A_{c})\right],$$

(26)

where the first term corresponds to the static dc current, which crucially depends on the superconducting phase difference $\phi$. In particular, the first term is equal to zero at $\phi=\pi$, in contrast to the second term,

$$\delta\bar{I}(A_{c})=\dfrac{1}{2\pi}\dfrac{\Delta_{d}^{2}(\Gamma_{n}^{2}+\varepsilon_{d}^{2})}{\mathcal{D}^{2}(0,\phi)}\int_{0}^{2\pi}d\psi\dfrac{\sinh^{2}{\left(\sqrt{A_{c}}\sin{\psi}\right)}}{\mathcal{D}^{2}\left(\sqrt{A_{c}}\sin{\psi},\phi\right)}>0,$$

(27)

which emerges exclusively due to the self-sustained oscillations and equals zero if the static state is stable, as indicated by the Heaviside step function $\theta(\kappa\varepsilon_{d})$. Plots of $\bar{I}_{n}$ as a function of $\tilde{\varepsilon}_{d}$ at $\tilde{\Delta}_{d}=1$ and $\phi=\pi/3,2\pi/3,\pi$ are presented in Fig. 5. These graphs show that the nanomechanical instability discussed in this article leads to the emergence of significant diode and transistor effects. The effects are most pronounced at $\phi=\pi$ when in the static regime $A_{c}=0$ where the Cooper pair exchange between the dot and the superconducting leads is completely blocked. In such a situation, a jump in the average current from zero to a finite value $\sim I_{0}$ (or vice-versa) occurs if the direction of the bias voltage changes (diode effect) or if the position of the level $\varepsilon_{d}$ controlled by the gate voltage passes zero (transistor effect). Note that the discontinuity of the average current as a function of $\varepsilon_{d}$ must be treated to the accuracy $\xi$ accepted in this paper.

In summary, we considered a nanoelectromechanical system comprising a carbon nanotube suspended above a trench in a normal metal electrode that undergoes bending vibrations in the gap between two superconductors. The nanotube was treated as a movable single-level quantum dot in which the superconducting order parameter is induced as a result of Cooper pair exchange with the superconductors. The latter essentially depends on the bending of the nanotube and the phase difference between the superconductors. We have shown that in such a system, the static, straight configuration of the nanotube is unstable regarding the occurrence of self-sustained bending vibrations in a wide range of parameters if a bias voltage is applied between the normal and superconducting leads. It was demonstrated that the occurrence of this instability crucially depends on the direction of the bias voltage and the relative position of the QD level. This makes it possible to govern the operating mode of the system by changing the bias and gate voltages. We have also shown that the appearance of self-sustained mechanical vibrations strongly affects the dc current through the system, leading to transistor and diode effects. The latter can be used for the direct experimental observation of the predicted phenomena.

O.M.B. thanks A.V. Parafilo and O.A. Ilinskaya for helpful discussions. S.I.K. acknowledges the financial support from the NAS of Ukraine (grant F 26-4). This work was supported by IBS-R024-D1.