The effect of elastic disorder on single electron transport through a buckled nanotube

We study transport properties of a single electron transistor based on elastic nanotube. Assuming that an external compressive force is applied to the nanotube, we focus on the vicinity of the Euler buckling instability. We demonstrate that in this regime the transport through the transistor is extremely sensitive to elastic disorder. In particular, built-in curvature (random or regular) leads to the ``elastic curvature blockade'': appearance of threshold bias voltage in the $I$-$V$ curve which can be larger than the Coulomb-blockade-induced one. In the case of a random curvature, an additional plateau in dependence of the average current on a bias voltage appears.


I. INTRODUCTION
A global trend of modern electronics is the design of nanodevices with ultra-low power consumption and a high level of integration. One of the most attractive candidates for this purpose is a single-electron transistor (SET) which is a sensitive electronic device based on the Coulomb blockade effect. Key points of fundamental physics of SET and its operation as a tunable nanodevice have been formulated about 20 years ago (for review see [1][2][3][4][5]). There has recently been an increased interest to this topic partially associated with the discovery of carbon systems with a transverse degree of freedom, such as suspended nanotubes and graphene. The interest is motivated by creation of nanoelectromechanical systems (NEMS) which are a class of devices integrating electrical and mechanical functionality on the nanoscale [6][7][8][9][10][11][12][13][14][15][16][17][18]. The simplest model of SET-based NEMS is provided by a harmonic mechanical oscillator coupled to an excess particle number on the SET island [16][17][18][19][20][21][22][23][24][25].
When suspended elastic nanotube is used as a SET island, coupling between mechanical and charge degrees of freedom provides an additional control via mechanical forces. This coupling can be strongly increased by applying a compressive force driving the nanotube towards the Euler buckling instability [26,27] (see more recent experimental [28][29][30][31][32][33][34][35][36] and theoretical [16][17][18][37][38][39][40][41][42][43][44] studies of this instability). Physics of such systems is captured by the model developed in Refs. [16][17][18]. This model was focused on the study of a SET made of a clean nanotube in which fluctuations arise only due to temperature T . Such fluctuations lead to "elastic blockade": appearance of a small threshold bias voltage in the SET I-V curve at the center of Coulomb blockade peak [16][17][18] (the term elastic blockade was introduced earlier for other materials [8,13]). However, the extreme sensitivity of the nanomechanical SET to an external environment implies that presence of disorder or some built-up deflection from ideal symmetry can strongly modify the elastic blockade.
Typically, nanotubes are not ideal and, therefore, are curved (in a random or regular way) in their equilibrium state. The effect of such built-in curvature is manifold. First of all, it can change the effective potential well confining electrons in the SET island and redistribute the electron density. This effect is weak provided that the curvature-induced electronic potential is small as compared to the Fermi energy. Secondly, built-in curvature can lead to an electron scattering. This effect is also weak, since a typical spacial scale of curvature is much larger than the Fermi wavelength. The most significant effect is related to the presence of elastic degree of freedom. In particular, as was shown in Refs. [41][42][43][44], builtin curvature can strongly affect buckling transition.
These studies have a great relevance to experiment. Experimentally, the SETs with elastic degree of freedom can be realized and tested in various systems, some of which are sketched in Fig. 1. The simplest realization is a nanotube located in a plane and coupled, both electrically and mechanically, with source and drain (see the left panel of Fig. 1). Mechanical coupling leads to buckling within the same plane. Physically, this system is fully equivalent to a graphene nanoribbon, which buckles in the direction perpendicular to plane (see the central panel in Fig. 1). Additionally, electro-mechanical coupling can be studied in an array of nanotubes shown in the right panel of Fig. 1.
Recently, a possibility of a complete control of nanotube-based SETs by means of an external mechanical forces has been clearly demonstrated experimentally [33][34][35][36]. In particular, in a very recent experimental work [35] a double-well elastic potential describing collective coordinate Y was created and controlled mechanically by using two mechanical forces: a compressive force along the nanotube which leads to buckling bistability and a force in the perpendicular direction, which makes the double-well potential asymmetrical thus "helping" the nanotube to choose the proper potential minimum. This experimental situation exactly corresponds to the theoretical model we shall focus in this paper. In Ref. [33] the lowest-lying flexural eigenmodes of nanotube were identified with an unprecedented precision by capacitively driv-FIG. 1. Sketch of a possible realization of SETs with an elastic degree of freedom. Left: A SET based on a suspended nanotube tunnel-coupled to a source and a drain and driven by an external force F. A gate electrode creates additional electric field E in the y direction. This transverse field bends the nanotube provided it is occupied by an excessive electron. Central: A SET with a graphene nanoribbon buckled by an external force. Right: A SET made of an array of nanotubes.
ing the nanotube-based mechanical resonator with external circuit. Excitation of such eigenmodes by placing the nanotube into the optical cavity shows surprisingly high quality factors (up to 10 4 !) [34]. The latter proves extreme sensitivity of suspended nanotubes to external forces. SET coupled to oscillating field was realized in Ref. [36] and field-driven Coulomb blockade peaks were used to make a single-electron "chopper".
In this paper, we study disordered nanotube-or nanoribbon-based SET under the action of a mechanical force (see Fig. 1). We model disorder by spontaneous local curvature (regular or random) [41][42][43][44]. We demonstrate that the results obtained in Refs. [16][17][18] for a clean nanomechanical SET are strongly modified in the presence of such curvature. This built-in curvature leads -due to electro-mechanical coupling -to two crucial effects: (i) existence of large threshold bias voltage in I-V curve at the center of Coulomb blockade peak for a fixed, built-in curvature (see Fig. 2), (ii) appearance of an intermediate plateau in bias-voltage dependence of the current averaged over realizations of random curvature (see bottom panel in Fig. 10).
Outline of the paper is as follows. We start with formulation of the model in Sec. II. In Sec. III the approximate treatment of the model is developed. The current for a fixed curvature of nanotube is calculated in Sec. IV. The case of a random curvature is considered in Sections V and VI. We end the paper with discussions and conclusions in Sec. VII. Some technical details are given in Appendices.

II. MODEL
The sketch of the setup is shown in Fig. 1 (left panel). We consider a nanotube of length L suspended between left and right leads, which serve as source and drain, respectively. The tunneling coupling of the nanotube to the source (drain) is characterized by the tunneling rates Γ L (Γ R ). We assume that a finite voltage bias V is applied between source and drain shifting their chemical potentials to ±eV /2. The electrical potential V g of the nanotube can be tuned by a capacitively coupled electrode (gate). We consider regime of strong Coulomb blockade, T <E c , such that only a single excessive (in addition to a quasi-neutral Fermi gas) electron can occupy a nanotube. Here, E c is the charging energy (E c ∼10 meV for a typical nanotube length L∼0.1 µm). The Hamiltonian of a nanotube, including the elastic degrees of freedom parametrized by the bending angle θ as a function of the arc-length position s readŝ Here θ ≡dθ/ds, κ denotes the bending rigidity of the nanotube, F stands for the compressing force applied to the ends of the nanotube, E is an additional electric field created by the gate electrode, directed along y axis, and is the average displacement of the nanotube in y direction. Here, y(s) is transverse coordinate of the nanotube, related to bending angle as follows, We model disorder by built-in curvature described by function C(s), which can be regular or random. The electronic degrees of freedom enters the Hamiltonian (1) via the operator of the excess particle number,n d , and the wave function, ψ 0 (s), of an active energy level for an excessive electron on the nanotube. Provided the nanotube is occupied by an additional electron, the quasi-neutrality is broken, and the transverse electric field E bends the nanotube. Hence, there is the term eEn d Y in Eq. (1). The nanotube is assumed to be clamped on the left and right leads of the SET, such that Hamiltonian (1) has to be supplemented by the following curvature-dependent boundary conditions [41,42], In the Hamiltonian (1) we neglect the elastic kinetic energy. This is justified under assumption that a typical elastic frequency of transverse oscillations of the nanotube, ω∝ κ/(ρL 4 ), where ρ is the linear mass density of the nanotube, is small as compared to temperature, ω T . As a result, the only dynamical degree of freedom in the Hamiltonian (1) is the excess particle number,n d . In the strong Coulomb blockade regime, the operatorn d takes two values, n d =0 and n d =1. Solution of the coupled -electronic and elastic -dynamical problem can be essentially simplified in the adiabatic approximation provided the tunneling between the nanotube and the source and the drain is sufficiently fast, Γ L,R ω. In such adiabatic approximation and under assumption, Γ R,L max{T, |eV |, |eV g |}, the computation of the bending angle becomes fully classical and the bending angle profile is described by the following nonlinear equation (for C(s)≡0, this equation was derived in Ref. [17]) where E(s)=E L s ds |ψ 0 (s )| 2 and n d (Y ) is the average value of the operatorn d for a fixed configuration θ(s), Here f F (ε)=1/[exp(ε/T )+1] denotes the Fermi function, V ± = V g ±V /2, γ − =Γ L /(Γ L +Γ R ), and γ + =Γ R /(Γ L +Γ R ). Under the same assumptions, the SET current reads In order to compute the current by means of Eq. (7), one needs to solve the nonlinear second order differential Eq. (5). Although it can be done numerically, as we shall demonstrate below, one can construct an analytic solution near buckling instability.

III. FUNDAMENTAL MODE APPROXIMATION
In order to satisfy the boundary conditions (4), it is convenient to expand the functions θ(s) and C(s) in the Fourier series, where q n =πn/L. Euler buckling instability is clearly seen from Eq. (5) at E=0. Indeed, linearization of this equation yields divergence of θ 1 , for F =F c : for the critical force of the instability (see Refs. [26,27]).
Hereafter we assume that | | 1. In fact, the growth of θ 1 for | |→0 is limited by nonlinear terms in the Eq. (5). The modes θ n with n 2 are finite at F =F c and can be found within the linear approximation: θ n ∼C n . Such crucial difference between the behavior of the mode with n=1 and modes with n 2 allows us to project Eq. (5) onto the fundamental mode cos(q 1 s) as it was done in the clean case [16,17]. In the most of our paper, we throw out all but the first mode (weak effects coming from modes with n 2 are briefly discussed at the Appendix A). The shape of nanotube, within fundamental mode approximation can be found from Eq. (3) and is given by We shall also make several simplifying assumptions throughout the paper. We assume the symmetric SET, Γ R =Γ L =Γ (i.e., γ + =γ − =1/2). Also we set T =0. We shall estimate the effect of thermal fluctuations at the end of the paper and demonstrate that for realistic values of parameters this effect is small. The latter assumption allows one to replace Fermi functions entering Eqs. (6) and (7) with step-functions. We shall also assume that disorder is weak so that |θ 1 | 1. Then, expanding sinus in Eq. (5) as sin θ≈θ−θ 3 /6, keeping only the first harmonic, θ≈θ 1 cos(q 1 s), and projecting thus obtained equation onto the first mode, we obtain the closed nonlinear balance equation for the amplitude of the first harmonic θ 1 : where Here Θ(x) denotes the Heaviside step function. Within the same approximation the current is given by In Eqs. (12), (13), and (14) we introduced dimensionless strength of the electron-phonon interaction and the dimensionless voltages [45] where a 1 = L 0 ds |ψ 0 (s)| 2 sin(q 1 s) is the numerical coefficient which fully encodes the dependence on the wave function of the active electron quantum level: it changes from a 1 =8/3π for ψ 0 (s)= 2/L sin(q 1 s) to a 1 ≈2/π for ψ 0 (s)= 2/L sin(q N s) with N ≈k F L/π 1. Here k F stands for the Fermi momentum.
Importantly, the typical magnitude of α is very small, α∼10 −4 ÷ 10 −2 [17], that will allow us to study the effect of the electro-mechanical coupling perturbatively.
We are interested in the solution θ 1 (C 1 ) of Eq. (11) that provides the minimum of the dimensionless energy related to f (θ 1 ) as follows: f (θ 1 )=∂W (θ 1 )/∂θ 1 . Here we introduce We note that in the single mode approximation the energy, E(θ 1 ), corresponding to Hamiltonian (1) is connected with the dimensional energy W (θ 1 ) as follows For α=0 the energy W (θ 1 ) has the form of the Landau expansion of the free energy in series of the order parameter θ 1 near the second order phase transition at =0. The curvature C 1 plays the role of the symmetry-breaking field that breaks the symmetry between two states of minimal energy at >0. Hence, the system can show the hysteretic behavior. However, since we focus here in dc current, we assume that the system always has time to reach the absolute energy minimum.
Most importantly, the single-mode approximation allows us to find the I-V curve for a fixed curvature (i.e. for given C 1 ). For random curvature which is characterized by a certain distribution function P(C 1 ), this approximation allows us to calculate analytically the distribution function for the current and analyze the competition of the random curvature with the term αw d (θ 1 ), which suppresses the bistablity. Below, we discuss the cases of fixed and random curvature in Sections IV and V, respectively.

IV. CURRENT FOR A GIVEN CURVATURE
In order to calculate θ 1 for fixed value of C 1 one needs to find all solutions of the balance Eq. (11), calculate corresponding energies with the use of Eq. (16) and find global energy minimum by choosing solution with the lowest energy.
A. Clean case in the absence of electro-mechanical coupling (C1=0, α=0) In order to set notations we start from the analysis of a clean nanotube without electromechanical coupling. Then, Eq. (11) simplifies drastically, 3. The effective energy for the clean nanotube at α=0 below (a) and above (b) the buckling instability threshold. Above the threshold the instability leads to a mechanical bistability in the clean nanotube. Disorder breaks the symmetry between two bistable states and the global energy minimum appears (c), which depends on a specific disorder realization.
This equation describes conventional buckling instability [27]. It has a single stable solution, θ 1 =0 for <0, and two stable solutions and an unstable one for >0: The effective energy for θ 1 mode looks: Function W (θ 1 ) is plotted in Fig. 3a for <0 and in Fig. 3b for >0. As one can see, the system shows mechanical bistability above the instability threshold, >0.
B. Disordered case in the absence of electro-mechanical coupling (C1 =0, α=0) Disorder modifies the balance equation and the effective energy as follows The effective energy for C 1 = 0 is shown in Fig. 3c. Solutions of Eq. (22) depend on C 1 , so that for >0 and at sufficiently weak disorder there are two stable solutions, θ ± (C 1 ) (see blue thick curves in Fig. 4), and one unstable solution (see dashed curve in Fig. 4), just as in the clean case. It is easy to check that corresponding to the absolute energy minimum of W (θ 1 ) is shown in Fig. 4 by red thick curve for > 0. As seen, there is a jump at C 1 = 0, corresponding to transition from θ − to θ + . Hence, in the case > 0 the dependence of θ 1 on C 1 can be written as C. Effect of electro-mechanical coupling For α =0, the step functions entering both f (θ 1 ) and W (θ 1 ) depend on v ± [see Eqs. (12), (13), (16), and (17)]. As follows from Eqs. (13) and (17), the quantities n d and w d can take different values depending on relation between θ 1 and v ± : It is convenient to introduce three functions, that correspond to different possible values of n d . Here, we denote These functions are shown in Fig. 5 by dashed lines. As seen from this figure, solutions of equations f s (θ 1 )=0 give (for not too large magnitudes of C 1 and α) six values of θ 1 (we do not consider here three unstable states with θ 1 close to zero): where functions θ ± (C 1 ) are shown in Fig. 4. We notice that only two solutions, θ ± b , belong to the interval v − <θ 1 <v + and, therefore, correspond to nonzero current through the nanotube [see Eq. (14)]. Using Eqs. (16) and (25), one can easily find energies corresponding to the solutions θ ± s : Importantly, ω ± s are functions of C 1 and α only and do not depend on v g and v.
Next step is to find regions in the plane (v g , v), where one or several solutions (28) are realized and, in the case of several solutions, chose solution with the lowest energy. This procedure is illustrated in Fig. 5. Function f (θ 1 ) defined by Eq. (12) is non-monotonous (see red curve in Fig. 5a) and have jumps at the points θ 1 =v + and θ 1 =v − . Equation f (θ 1 )=0 yields several solutions for θ 1 belonging to manifold (28). For example, in Fig. 5a the voltages v ± are such that there are three solutions: θ − b , θ − c and θ + c . These solutions give the three local minima of the energy W (see Fig. 5b). As one can see, the global minimum correspond to θ + c , so that current through nanotube is zero. This example illustrates physical origin of the elastic blockade: local minimum θ − b , corresponding to non-zero current, does not give absolute minimum. Hence, an electron placed in this local minimum should decrease its elastic energy by re-buckling from with simultaneous "jumping out" from the current-carrying window.
Let us now discuss the general case. First we notice, that (v g , v) plane can be divided onto several regions, where different solutions (28) can be realized. These regions corresponds to different regions in Eq. (25).
Since θ ± s do not depend on v g and v, these regions are limited by lines: v/2=±v g +const. They are shown by dashed lines in Fig. 6, where we plotted only half plane v>0 (the diagram is symmetric with respect to change v→ − v and change I→ − I). Dashed lines connected by arrowed arcs indicates regions where different solutions (28) exist, e.g. regions a ± corresponds to solutions θ ± a etc. As seen, these regions overlap, so that several solutions can coexist in agreement with (c − , b − , c + ). Therefore, one has to calculate energies of coexisting states (W − c , W − b , and W + c in the above example) and to find the global minimum. Within the pink triangle in Fig. 6 the global minimum is given by one of the energies W − b or W + b , corresponding to non-zero current, I = 1/2. Boundary of this triangular region, shown by two red thick lines, represents the threshold voltage, v Th =v Th (C 1 , v g ). This voltage separates the region with zero current (v<v Th ) from the region (v>v Th ), where I = 1/2. Position of the triangle depends on C 1 in a nontrivial way. Analyzing Eqs. (29) and (30) one can demonstrate that the left (L) and the right (R) boundaries (red lines) are determined by the following conditions: Using these relations and Eq. (29), we find that the threshold voltage can be written as follows: Dependence of v Th =v Th (v g , C 1 ) on v g for fixed C 1 is shown by red lines (boundaries of the current-carrying triangle) in Fig. 6. The dependence on C 1 is fully en- coded in functions v 0 =v 0 (C 1 ) and v g0 =v g0 (C 1 ): Equations (32), (33), and (34) provide a full solution of the problem for a fixed disorder. Functions ω ± s entering these equations are given by Eq. (30) with θ ± s found from Eq. (28).
For a fixed disorder the current-voltage dependence becomes Equations (33) and (34) allows one to find some general properties of the current-carrying triangle. In particular, one can easily demonstrate that the curve v 0 (C 1 ) has a maximum at C 1 =α/16 and is symmetric with respect to this point: v 0 (C 1 )=v 0 (α/8−C 1 ). Also, one finds v g0 =0 at C 1 =α/16 and is anti-symmetric with respect to this point v g0 (C 1 )=−v g0 (α/8−C 1 ).

D. Limiting cases
Next, we consider some limiting cases. General equations derived in the previous section are dramatically simplified when C 1 and α are small: Then, keeping in θ ± s terms up to linear order with respect to C 1 and α, and in ω ± s up to quadratic order we get Substituting Eq. (38) into Eqs. (33) and (34) we find that the latter equations can be written in a compact way (it is convenient to express v 0 and v g0 in terms of Dependencies of v 0 and v g0 on C 1 are schematically shown in Fig. 7. Most interestingly, there is a very sharp increase of v 0 (and, consequently, v Th ) within a narrow region of disorder strength: 0<C 1 <α/8. Hence, we predict the suppression of the current by the built-in curvature-the phenomenon, which we call elastic curvature blockade. Physically, curvature-induced threshold in the current arises due to the need to choose the bending amplitude θ 1 from minimization of W (θ 1 ). For |v|<v Th , the bending angle found from energy minimization (with account of the electro-mechanical coupling) turns out to be beyond the current-carrying window v − <θ 1 <v + .
One can also find dependence of v 0 on . The results of corresponding numerical solution is shown in Fig. 8. As one can see, this dependence is very sensitive to a built-in curvature. In particular, the change of C 1 from negative to positive values results in qualitative change of the behavior of v 0 with .
Equations (39) and (40) can be further simplify in the absence of disorder, C 1 =0. In this case, v 0 =α/4 v g0 =− √ −α/4 , so we obtain v Th = α 4 This result is in agreement with Refs. [16][17][18], where it was found that electro-mechanical coupling leads to elastic blockade. Indeed, we see that there is a finite minimal threshold voltage, v min Th = α/4 , proportional to electromechanical coupling. This means that at low temperatures and small driving voltage, the transport through SET is blocked.
Physically, approximation (43) neglects relatively small effects caused by elastic blockade in the absence of curvature, but captures much larger curvature-induced blockade. Within this approximation, position of the currentcarrying triangle is fully determined by parameter ξ = 16C b /α=16C 1 /α−1. With increasing ξ from −∞ to ∞ the lower tip of the triangle moves in the following way: it remains stationary in the position (v g0 =− √ , v 0 =0) for −∞<ξ<−1, then "climbs a hill" with a height 2 √ in the interval −1<ξ<0, moves down in the interval 0<ξ<1, and finally arrives at the position (v g0 = √ , v 0 =0) and remains stationary for 1<ξ<∞. This approximation neglects very slow motion of the triangle at large |ξ|>1 (see the bottom panel in Fig. 7). The change of v g0 from − √ to √ with increasing ξ is due to curvature-induced re-buckling of the nanotube.

V. DISORDER-AVERAGED I − V CURVE
Next we assume that a random C 1 is described by a distribution function P(C 1 ). Changing C 1 in the interval, −∞<C 1 <∞, one can find the curve (v g (C 1 ), v 0 (C 1 )) in (v g , v) plane parameterized by C 1 . This curve yields limiting value v Th (v g ), below which the current is equal to zero for any C 1 and consequently, for disorder-averaged I − V curve. Using general properties of functions v 0 (C 1 ) and v g0 (C 1 ) [see discussion after Eq. (35)], one can easily find that v Th (v g ) has a maximum at v g =0 and is symmetric with respect to this point: v Th (v g )=v Th (−v g ). Dependence v Th (v g ) is shown schematically in Fig. 9. For small α obeying inequality (36), this dependence reads where v g = √ [1+α/(4 3/2 )]. One can also show that for very large gate voltages, v g √ , the threshold voltage decays as v Th (v g )≈α/(6v 2 g ). We stress that v Th (v g ) does not depend on the distribution function of the random curvature, P(C 1 ), provided that this function is non-zero within the whole interval −∞<C 1 <∞. To be specific, we assume that the random curvature has Gaussian distribution with the zero mean and where C (s)=dC(s)/ds. Then the eigenmodes C n are independently correlated with the dispersion determined by ∆, C n C m = ∆ n δ nm , ∆ n = 2∆L π 2 n 2 , n = 1, 2, . . . . (46) Hence, C 1 has the normal distribution, characterized by the variance ∆ 1 . The point {v g0 (C 1 ), v 0 (C 1 )} is the position of the lower tip of triangle region of non-zero current. This point moves as C 1 is increasing from negative to positive values.
As shown in Fig. 2, for v>v Th (v g ), there are two values of curvature C R and C L (C L >C R ) for which, respectively, right and left boundaries of the current-carrying triangle cross a point (v g , v). Since at any point belonging to the triangle the current equals I = 1/2, and equals to zero outside the triangle, we get the following expressions for averaged current and the current distribution function As follows from Eqs. (32), C R,L are implicitly defined by the following equations: Here functions v 0 (C) and v g0 (C) are determined by Eqs. (28), (30), (33), and (34). Equation (48) suggests that I cannot exceed the value 1/2.
The density plot of the average current calculated numerically with the use of Eq. (48) as the function of v and v g is shown in the top panel of Fig. 10. Numerically calculated current-voltage curve is presented in the bottom panel of Fig. 10 for two different values of variance ∆ 1 . One can also obtain analytical solution for I − V by using approximation (43). Corresponding calculations are delegated to Appendix B.
Analysis of the numerical and analytical solutions demonstrate that at fixed v g , the value 1/2 for the current is reached at large values of v when C L →∞, C R →−∞. One of the most interesting properties of the averaged I − V curve is appearance of an intermediate plateau at I=1/4 with increasing disorder variance. This property is illustrated in the bottom panel of Fig. 10 and in Fig. 16. In order to understand physics behind this phenomenon, we will discuss in the next section the distribution function of the angle θ 1 .

VI. THE DISTRIBUTION FUNCTION OF THE BENDING ANGLE, θ1
The probability distribution function (PDF) for the bending angle θ 1 is given by where function θ 1 (C 1 ) gives solution of equation f (θ 1 ) = 0, corresponding to the global minimum of W and angular brackets stand for the averaging with function P(C 1 ) [see Eq. (47)]. Average current can be found (instead of integration over C 1 , according to Eq. (48)) with the use of P (θ 1 ) by integration within the current-carrying window:

A. Distribution of the bending angle in the absence of electro-mechanical coupling
Let us first consider the case of zero electro-mechanical coupling. For α=0, θ 1 (C 1 ) is shown in Fig. 4 and for >0 is expressed in terms of θ ± (C) according to Eq. (24).
In the limit of a very weak disorder, the distribution function consists of two delta peaks, Here ∆ = √ Θ( ) is the buckling-induced gap. The disorder broadens the delta functions in Eq. (55).
The function P (θ 1 ), described by Eq. (53), is plotted in Fig. 11. for a fixed weak disorder and two values of , below and above instability threshold. Above the threshold, >0, there is a gap |θ 1 |<∆ inside which P (θ 1 )≡0. We also notice that slightly above the gap the distribution function increases despite of the presence of the damping exponent in Eq. (53). This increase is due to the Jacobian, J(θ 1 ).
Although the terms involving higher curvature harmonics C n>1 are small, they can modify the result (53) in the region where P (θ 1 ) is very small. In particular, they lead to appearance of tails inside the gap, |θ 1 |<∆ , for >0 and to some modification of P (θ 1 ) for <0. We discuss in detail the effects of the higher harmonics in Appendix A and demonstrate that these modifications do not lead to significant change of the distribution function. In particular, as one can see from Fig. 14, even in the case of a moderate disorder, when delta peaks in the distribution function at θ 1 =∆ are essentially broadened, the smearing of the gap edges is quite small. Therefore, we shall neglect the corrections due to the high-order harmonics to P (θ 1 ) in the rest of the paper.

B. Effect of the electro-mechanical coupling
As we demonstrated above, equation f (θ 1 )=0 has several solutions belonging to the manifold θ ± s (C 1 ) (see Fig. 5). These solutions and corresponding energies depend on C 1 . Since function f (θ 1 ) has jumps at points θ 1 =v − and θ 1 =v + , the function θ 1 (C 1 ), which, by definition, represents the solution with the minimal energy, might have mini-jumps in addition to a large jump existing in the case α=0. Existence of these mini-jumps as well as their positions and position of the large jump depend on v g , v, and on the value of α. For illustration, in Fig. 12 we plotted θ 1 (C 1 ) and corresponding distribution function, (here C 1 (θ 1 ) is the inverse function for θ 1 (C 1 )), for the case v>v Th (v g ) and relatively large v g such that C R >α/8. As seen from this figure, there are three major effects of the electromechanical coupling: (i) the curve θ 1 (C 1 ) becomes asymmetric with respect to the inversion: θ 1 → −θ 1 ; (ii) several vertical mini-jumps form on this dependence in addition to the large jump; (iii) large jump between values of the bending angle close to − √ and √ shifts from C 1 =0 to the point C 1 =α/8 (this specific value depends on a position of point in the (v g , v) plane).
Let us discuss these properties in more detail. First, we notice that W + s −W − s change sign when C s =0 (in particular, for |C s | 3/2 , we get W + s −W − s ≈−16C s √ ). Since, by assumption, C R >α/8, large jump occurs with increasing C 1 at the point C a =0, i.e. at C 1 =α/8, where W + a =W − a and, consequently, there happens a jump from θ − a to θ + a . With further increase of C 1 it reach the value C 1 =C R , where solution jump into current-carrying triangle θ + a →θ + b , and further jumps out the triangle, The values of C L,R can be found from Eqs. (39), (40) and (49): The averaged current is given by Eq. (48), while the amplitudes of jumps read: . For small α, obeying Eq. (36), we obtain As seen from Fig. 12, probability that the bending angle is close to the value − √ is larger than the probability that θ 1 ≈ √ . The probability imbalance increases with decreasing ∆ 1 . For √ ∆ 1 α, the bending angle with exponential precision is concentrated near − √ . Let us now discuss the physics behind the additional random-curvature-induced plateau in the I − V curve. The key point is the competition between two symmetry breaking mechanisms: electro-mechanical coupling and random curvature. Remarkably, the limits α→0 and ∆ 1 →0 do not commute. Indeed, for α≡0 and variance tending to zero, ∆ 1 →0, energy W (θ 1 ) has two degenerate minima for >0. Thus, distribution function is given by the sum of two delta-functions: Although a weak disorder slightly broadens the delta functions, the symmetry P ∆ (θ 1 )=P ∆ (−θ 1 ) remains intact. In the considered limit, ∆ 1 →0 and α=0, expression for average current reads One can easily check that this expression contains an intermediate plateau with I = 1/4 in a wide interval of v in addition to the plateau with I = 1/2 for v → ∞. By contrast, setting ∆ 1 =0 first and, then, sending α→0, we find asymmetrical distribution function, Corresponding expression for the average current, has only a single plateau with I = 1/2 for v → ∞. The competition between two mechanisms is controlled by the parameter √ ∆ 1 /α (see also Appendix B). An additional step in the I − V curve appears when this parameter becomes large: √ ∆ 1 /α 1. Remarkably, the current can be enhanced by disorder contrary to naive expectations (see the bottom panels in Fig. 10 and Fig. 16b).

VII. DISCUSSIONS AND CONCLUSIONS
In this paper, we predict curvature-induced enhancement of elastic blockade. The maximal position of the tip of the current-carrying triangle is given by (15) and (39)]. Hence, For E 10 kV/cm, L 0.1 µm and ∼1 it gives eV max 0 ∼0.1 eV that is larger than the charging energy E c ∼10 meV. We note that enhanced bias-voltage threshold of the order of eV max 0 exists for a sufficiently weak built-in curvature, in a narrow window, 0<C 1 <α/8. Away from this window the bias-voltage threshold is parametrically reduced (for α 3/2 ) and becomes of the order of the threshold in the clean case [16][17][18]. From Eq. (41), we find v cl 0 =α/4 , so that Eqs. (59) and (60) imply that bias threshold voltage very sharply depends on disorder: it changes by a large factor 3/2 /α in a very narrow region of small C 1 . For the experiments of Ref. [46] we can roughly estimate absolute value of curvature |C 1 | 0.1÷1. For this estimate we used Eq. (9) with F =0 (i.e. =−8) and values of θ built−in 1 estimated from images of bent nanotubes in Ref. [46]. We also assumed that in equilibrium C 1 ∼θ built−in 1 (this estimate follows from Eq. (5) for F =0 and small α). Although the microscopic model of built-in curvature is absent so far, one may expect that disorder-induced bending angle is not very sensitive to length of the nanotube. Due to smallness of α, the above estimates imply that built-in curvature might dominate electron transport through a nanotube-based SET for typical values of experimental parameters. It is worth noting that the threshold voltage given by Eq. (59) does not depend on the curvature but increases with length. Therefore, for long nanotubes used in Ref. [46], the effect of elastic blockade, which we discuss in the current work, is even higher and value eV max 0 is even larger compared to the estimate presented above. A more detailed comparison with experiment requires development of the microscopical theory of the elastic disorder, and therefore is out of scope of the current work devoted to development of the phenomenological approach.
In the previous sections we assumed that T =0. Let us now discuss the effect of nonzero temperature. For a clean SET these effects were analyzed in Refs. [16][17][18]. The current depends on T due to at least two reasons. At first, temperature enters the Fermi distribution functions in the expressions (6) and (7) for the average excess electron number and the current, respectively. Secondly, there are thermal fluctuations of the generalized coordinate Y , i.e. θ 1 . In the fundamental mode approximation these thermal fluctuations of the bending angle are described by the Gibbs factor exp[−E eff (θ 1 )/T ]= exp(−W eff (θ 1 )/T ) where W eff (θ 1 ) and E eff (θ 1 ) are obtained from Eqs. (16) and (18) by replacement [see Eqs. (19) and (25) of Ref. [18]] The effective dimensionless temperature is defined as As usual, a non-zero temperature results in smearing of the bias-voltage threshold in the current. For a fixed curvature C 1 , our results for the current are valid provided temperature is much smaller than the threshold voltage, Here voltage V 0 (C 1 )=(a 1 /π)ELv 0 (C 1 ), cf. Eq. (32), is a very sharp function of curvature within the window 0<C 1 <α/8, changing within the interval As follows from the estimates presented above, T is always smaller than V max 0 up to the room temperatures but can be much larger than V cl 0 . Hence, the effect of temperature is negligible for 0<C 1 <α/8, but it can wash out the bias-voltage threshold for C 1 outside this interval.
In the case of a randomly distributed curvature the enhanced bias-voltage threshold (59) occurs near zero gate voltage, V g 0.
Away from the elastic blockade peak the bias-voltage threshold is suppressed down to V cl 0 , cf. Eq.(60). The thermal smearing is not important at temperatures T eV Th (V g ) wherē V Th (V g )=(a 1 /π)ELv Th (v g ), cf. Eq. (44). For a random curvature distributed continuously, the average current is given as the average over P (C 1 ) and the Gibbs weight exp(−W eff (θ 1 )/T ). The result of the corresponding nu- and different values of disorder variance ∆1=8 · 10 −4 and ∆1=3 · 10 −6 (red and blue correspondingly) forT =0, 2.5 · 10 −3 , 5 · 10 −3 (solid, dashed, dotted curves, respectively). merical computation is shown on the Fig. 13. It illustrates blurring of the disorder-induced step in the I − V curve at intermediate voltages by a finite temperature.
One of the most interesting problems to be solved in future is to study effect of disorder on high-quality resonances observed in Refs. [33][34][35][36].
To summarize, we considered transport properties of a SET based on elastic nanotube with a regular or random built-in curvature in the strong Coulomb blockade regime. We demonstrated that close to the buckling transition, the I-V curve of the transistor is extremely sensitive to such curvature.
Most importantly, we predict the following • for fixed built-in curvature: a giant curvatureinduced enhancement of the bias-voltage threshold (59) below which the current is exactly zero; • for random curvature: the existence of an additional intermediate curvature-induced plateau in average I − V curve, see Fig. 10.
Here we calculated Jacobian by using Eq. (A1). We also took into account that corrections due to high harmonics are small and can be fully neglected in the argument of the the exponent. The averaging in Eq. (A7) is taken over C n>1 . Further consideration depends on relation between ∆ 1 and . Below, we consider limits of both weak and strong disorder.
In this case, one can neglect disorder everywhere except the arguments of the step function, where we can also neglect the difference between eff and . Then disorder averaging is reduced to averaging of the step functions: (A8) Here we used ∆ 3 =∆ 1 /9 (cf. Eq. (46)). The distribution function becomes Hence, for weak disorder, the only effect of the higher harmonics is blurring the step-like edge of the gap (for >0) over a very narrow window δθ 1 ∼ √ ∆ 1 √ .

Interpolation formula and gap formation
One can easily find a cross-over function that interpolates between weak and strong disorder limits. We note that both ∆ 1 and are still assumed to be smaller than unity. The cross-over function reads Equation (A13) describes formation of the gap in the distribution function. Due to coupling with higher harmonics, the edge of the gap is not infinitely sharp. Formally this is due to presence of ercf−functions in Eq. (A13) instead of Θ[θ 2 1 − ] in Eq. (53). Deep inside the gap, say for θ 1 = 0, perturbative expansion over harmonics fails and Eq. (A13) becomes invalid. Correct calculation of exponentially small tails of the distribution function in the middle of the gap can be done by using optimal fluctuation method. Such calculation is out of scope of the current work. We also notice that comparison of Eqs. (A13) and (53) shown in Fig. 14 shows that difference between distribution functions calculated with and without account of high harmonics is quite small.

Appendix B: Simplified model
Here, we discuss in more detail averaged current obtained within simplified model, when both C 1 and α tend to zero in such a way that the ratio C 1 /α remains finite, cf. Eq. (43) and >0. Approximation (43) allows us to find exactly values of C L and C R entering Eq. (48), and, consequently, to find simple analytical equation for the curvature-averaged I-V curve.
Using Eq. (43), one can easily find how left and right boundaries of the current-carrying triangle move with changing ξ in the interval −1<ξ<1: For ξ<−1, the triangle does not move and is limited by lines v L (v g , −1) and v R (v g , 0). For ξ>1, the triangle is also stationary and is bounded by lines v L (v g , 0) and v R (v g , 1).
We notice that in the interval −1<ξ<0, the left boundary changes according to Eq. (B1), while right boundary remains stationary and is still given by v R (v g , 0). For 0<ξ<1, left boundary stops and is given by v L (v g , 0), while right boundary starts to move according to Eq. (B2). This is very similar to Fig. 2 with two minor distinctions: in the simplified model the boundaries of the current-carrying triangle are parallel to blue lines and we neglect slow motion of triangle for C 1 < 0 and C 1 > α/8. region 1 are not covered by the current-carrying triangle, so that C L = C R = ∞ and I ≡ 0. By contrast, all points of the region 4 belong to the current-carrying cone for any C 1 , hence, I ≡ 1/2.
For region 2 we find C R = −∞, and Here C L is found from the condition v=v L (v g , 16C L /α−1). The averaged current in the region 2 reads In the region 3 we get C L = ∞ and This result for C R is derived from the condition v=v R (v g , 16C R /α−1). Averaged current in the region 3 reads Voltage-current curves calculated within simplified model for two values of the gate voltage v g =v g1 and v g =v g2 are shown in Fig. 16. The parameter α is kept fixed whereas the disorder strength varies. In both cases, for small driving voltage, in region 1, the current is equal to zero, while at very large v, in region 4, the current is equal to 1/2. In the region 2 (see Fig. 16a) the current is given by Eq. (B4). In the region 3 (see Fig. 16b) the current is given by Eq. (B6). The main difference of the simplified model as compared to the exact one is the presence of sharp jumps between different regions of the I − V curve. These jump are smoothed in a more accurate model (compare bottom panel of Fig. 10 and Fig. 16).