Mechanism of the Resonant Enhancement of Electron Drift in Nanometre Semiconductor Superlattices Subjected to Electric and Inclined Magnetic Fields

We address the increase of electron drift velocity that arises in semiconductor superlattices (SLs) subjected to constant electric and magnetic fields. It occurs if the magnetic field possesses nonzero components both along and perpendicular to the SL axis and the Bloch oscillations along the SL axis become resonant with cyclotron rotation in the transverse plane. It is a phenomenon of considerable interest, so that it is important to understand the underlying mechanism. In an earlier Letter (Phys. Rev. Lett. 114, 166802 (2015)) we showed that, contrary to a general belief that drift enhancement occurs through chaotic diffusion along a stochastic web (SW) within semiclassical collisionless dynamics, the phenomenon actually arises through a non-chaotic mechanism. In fact, any chaos that occurs tends to reduce the drift. We now provide fuller details, elucidating the mechanism in physical terms, and extending the investigation. In particular, we: (i) demonstrate that pronounced drift enhancement can still occur even in the complete absence of an SW; (ii) show that, where an SW does exist and its characteristic slow dynamics comes into play, it suppresses the drift enhancement even before strong chaos is manifested; (iii) generalize our theory for non-small temperature, showing that heating does not affect the enhancement mechanism and accounting for some earlier numerical observations; (iv) demonstrate that certain analytic results reported previously are incorrect; (v) provide an extended critical review of the subject and closely related issues; and (vi) discuss some challenging problems for the future.


I. INTRODUCTION
Spatial periodicity plays a fundamental role in nature. In particular it governs quantum electron transport in crystals 1 . In a perfect crystal lattice, an electron in a constant electric field would undergo Bloch oscillations, moving forwards and backwards periodically so that its average drift speed would be zero 1 . But real lattices are imperfect and electrons may be scattered before reversing their motion, allowing them to acquire a steady drift. Typically, the Bloch oscillation period t B greatly exceeds the average scattering time t s , because t B is proportional to the reciprocal of the lattice period d l , which is very small. So Bloch oscillations are not observed in real crystals. Nanoscale superlattices 2 (SLs) impose on the crystal an additional periodicity with a period d greatly exceeding d l but still small enough for the quantum nature of the electron to be important: t B may then become comparable to or smaller than t s so that Bloch oscillations can manifest themselves, significantly suppressing the current, generating gigahertz/terahertz electric signals, and causing many other important effects [3][4][5][6][7] .
Studies of SLs now constitute a significant area within solid state physics: see e.g. the major reviews by Wacker 3 and by Bonilla and Grahn 4 , and the book by one of founders of the area Raphael Tsu 5 . Our present research falls into one of its sub-areas, studying how magnetic field affects electron transport in SLs. Most works in this subarea considered cases when the magnetic field was either perpendicular to the SL axis [8][9][10][11][12][13][14][15] or, more rarely, parallel to it 10,16 . It would be natural to expect that an intermediate tilt would give rise either to a mixture of phenomena characteristic of the perpendicular and parallel orientations, or to a less pronounced manifestation of one of them, rather than to distinctly new phenomena. Nevertheless, it was predicted as early as 1980 by Bass et al. 17 that, if a constant electric field is directed along the SL axis while the constant magnetic field possesses both parallel and perpendicular components, then a new phenomenon should occur: as the field parameters are varied, the electron drift velocity should undergo distinct "resonant" changes when the period of the cyclotron rotation in the transversal plane and the period of the Bloch oscillation or any its multiple approach each other. These results appeared to be clear, scientifically interesting, and promising for applications. Perhaps, it was initially expected that it would be relatively quick and straightforward to realise them experimentally. But this turned out not to be the case. Rather, it was just the beginning of a long and tortuous path towards the truth. Despite the publication of numerous papers (e.g. 9, ) and the involvement of many distinguished physicists, both theoreticians and experimentalists, there still remain fundamental unanswered questions. A detailed discussion of the intricate development of the subject is given in Sec. II below while here we restrict ourselves to a brief discussion of works of an immediate relevance to the new developments presented in this paper.
Before presenting the discussion and formulating the purposes and the outline of the present paper, we describe the foundational work on superlattices by Esaki and Tsu 2 in order to place the discussion in context. Using the semiclassical model for the motion of electrons between collisions, and introducing the relaxation-time approximation for the effect of collisions (i.e. scattering) on electrons, Esaki and Tsu 2 showed that, in the zerotemperature limit, the drift velocity v d vs. the constant electric field F along a one-dimensional SL possesses a peak with a maximum at F = F ET such that t B (being ∝ F −1 ) is equal to t s . In what follows, we will refer to it as the Esaki-Tsu (ET) peak. It has important consequences, e.g. causing a peak in the differential dc conductivity vs. voltage (corresponding to the pronounced maximum in dv d (F )/dF on the left side of the ET peak) and current oscillations at high voltages owing to the negative sign of dv d /dF on the right side of the peak 3,4 .
So, let us briefly overview works of an immediate relevance to our present research. Firstly, it is the aforementioned pioneering work by Bass et al. 17 . Its result had not been exploited for more than 30 years (until the paper 30 ) and the authors of the next key papers 19,21 were apparently unaware of it. In 2001, Fromhold et al. 19 generalized the approach of Esaki and Tsu 2 on the case with the tilted magnetic field and, for a given set of parameters, their numerical calculations revealed a large pronounced peak in v d (F ) near F = F 1 and few smaller ones near F n>1 , where F k corresponds to the equality between the Bloch frequency ω B ≡ 2π/t B and a product of the cyclotron frequency with any natural number k. To distinguish all these peaks from the ET peak, which is situated at much lower values of F , Fromhold et al. called them as "resonance peaks". Besides, it was shown in 19 that the intercollisional semiclassical dynamics reduced to the dynamics of an auxiliary classical harmonic oscillator subject to a travelling wave. It has been known since the end of the 80th [43][44][45][46] that, if the wave frequency is equal or almost equal to any of multiples of the oscillator frequency, then the phase plane of the oscillator is threaded by the so-called stochastic web (SW) through which the oscillator may diffuse chaotically very far from the origin. Fromhold et al. 19 suggested that the resonance peaks originated just in this chaotic diffusion. The peaks promised to give rise to interesting observable consequences, in particular to associated peaks in the differential conductivity of the sample 3 (analogous to that associated with the ET peak). The latter peaks were soon discovered indeed 21 as well as their numerical theory was developed, being in a reasonable qualitative agreement with the experiment. The next work of the immediate relevance was the paper by Sel'skii et al. 30 , who generalized the consideration to non-zero temperatures more consistently than it was done in 21 and numerically calculated for a given example the evolution of v d (F ) as temperature grew. Their calculations showed in particular that, as temperature grew, the ET peak quickly decayed while the largest (1st) resonance peak did it much slower so that it became the dominating feature of v d (F ) at moderate and moderately high temperatures. The authors gave no explanation of such an evolution.
In all papers on the subject since 19 until the beginning of 2015, it had been assumed that the mechanism of the resonant peaks was the chaotic diffusion through the SW. Moreover, many papers on other subjects  referred to the phenomenon as being a characteristic manifestation of dynamical chaos in quantum electron transport. In 2015, it was shown in our Letter 37 for the zerotemperature limit that the assumption was wrong. Based on regular approximations of the exact collisionless equations of motion, an analytic solution of the problem was obtained that agreed well with a solution based on numerical integration of the exact equations. Thus the resonant peak was successfully accounted for without any need to invoke chaos. The mechanism underlying the phenomenon was also explained in qualitative terms. It was demonstrated 70 that our work allows us very easily to find the values of the physical parameters needed to maximize the drift, and that they differ markedly from those required to optimise chaotic diffusion, if this were the operative mechanism 29 . The possible effect of the diffusion was also considered.
Finally, the theoretical paper by Bonilla et al. 40 was published in the end of 2017: it was a challenge to all previous theories as it aimed to show that there should be such redistribution of electrons in space which necessarily causes the appearance of the transversal component of the electric field with a large absolute value.
We embarked on the present work with two main motivations. One of them was to generalize our theory to encompass non-small temperatures, thus allowing us to explain the results of 30 and to predict all possible scenarios for the evolution of the resonant peak with changing temperature. The other motivation was to extend the earlier arguments in favour of the non-chaotic mechanism 37 . Need for the latter arose because some researchers seemed to remain unpersuaded by the arguments in the Letter 37 and had been explicitly 38,71,72 or implicitly 39,41 continuing to refer to the original conjecture 19,21 of chaotic diffusion as the mechanism underlying the resonant enhancement of electron transport. It also seems that a huge number of references to this exciting-sounding but incorrect idea done before the publication of the Letter 37 and its continued promulgation 38,39,41,71,72 after it keep misleading researchers in other areas [73][74][75][76][77][78][79] , who still mention the resonant enhancement of electron drift 19,21,26 as a manifestation of "chaotic dynamics in semiconductor SLs" 73,75,76,78,79 or as the ability of non-KAM chaos to "enhance electronic transport in semiconductor superlattices" 74,77 .
After the publication of the paper 40 , two more motivations appeared: (i) to analyse the subject in general and its major unsolved problems, (ii) to analyse a relevance of the conclusion of 40 about the necessity to introduce a strong transversal electric field to our immediate problem i.e. to the resonant enhancement of the dc drift.
In what follows, we will therefore: (i) present a detailed general analysis of the subject intricate development and list the major unsolved problems; (ii) show that, due to certain geometrical feature of samples used in experiments done to the date, the conclusion of 40 does not relate to the immediate subject of our paper; (iii) provide additional arguments demonstrating that the resonant enhancement of the dc drift velocity cannot be attributed to chaotic diffusion; (iv) present in greater detail the asymptotic theory 37 describing the enhancement, and elucidate the real drift enhancement mechanism in physically-motivated terms; and (v) generalize the asymptotic theory for arbitrary temperature, thereby demonstrating that possible heating of the electrons does not affect the non-chaotic nature of the phenomenon while, at the same time, accounting for some numerical results reported earlier 30 . In addition, we will demonstrate the incorrectness of the earlier analytic results 17,18,30 based on non-chaotic approaches differing from that introduced in 37 and discuss some challenging problems that remain to be addressed.
The rest of the paper is organized as follows. Sec. II gives the critical overview of the general subject devel-opment. Sec. III introduces the model and basic equations. In Sec. IV, our asymptotic theory for the zerotemperature limit is presented and compared with numerical results, including with those from earlier works. It unambiguously proves the validity of the regular mechanism and the inapplicability of the chaotic one. We emphasise that our numerical simulations use the exact equations of motion, so that they cannot be affected by the approximations underlying our analytic theory. In addition to the analytic theory, we present in Fig. 5 evidence of strong resonant enhancement of the drift velocity arising under conditions for which the SW does not exist, so that the chaotic diffusion does not exist either. Sec. V considers the chaos influence on the drift in cases when the SW does exist. Apart from the analytic proof, Fig. 9 and the animation 80 demonstrate that, if the resonant enhancement is pronounced, then the time-scale relevant to the resonant drift formation (scattering time, in the given case) is much smaller than that at which chaos would start to manifest itself. Besides, a comparison of Figs. 10 and 8(c) shows that, if chaos does manifest itself on the relevant time-scale, then it suppresses the resonant drift. Thus, taken together, Figs. 5, 8(c), 9 and 10 immediately demonstrate the invalidity of the chaotic concept of the resonant drift origin. The asymptotic theory is generalized for the case of arbitrary temperature in Sec. VI and the results are verified numerically (again using the exact equations of motion), thus confirming that the regular mechanism is valid at any temperature. Sec. VII presents a discussion, in particular showing the incorrectness of some earlier analytic results 17,18,30 . The validity of the model is also discussed. Conclusions are presented in Sec. VIII. The Supplemental Material 80 expands on some of the details as well as including the animation.

II. MILESTONES IN THE DEVELOPMENT OF THE SUBJECT
We now review in chronological order the key milestones in studies of the effect of an inclined magnetic field on electron transport in the SLs. Not surprisingly, some authors were apparently unaware of earlier key results. Apart from the inherent difficulty of the subject, this provided an additional reason for its development to be so intricate.
The first work in the area was the paper by Bass et al. 17 . Unlike Esaki and Tsu 2 , the authors took account of collisions by use of the kinetic equation for the distribution function of the electron quasi-momentum p in the so called ν-approximation, when only inelastic scattering is assumed to occur, characterized by a single scattering rate ν. In the absence of a magnetic field and for a temperature close to zero, such an approach leads to the result of Esaki and Tsu 2 . Bass et al. 17 proved that the implicit integral representation of the drift velocity via the quasi-clasical instantaneous velocity 2 can be generalized for the presence of a classical magnetic field and for non-zero temperature (details of the proof are given in their longer paper 9 ). Furthermore they noticed that, if the magnetic field is inclined. i.e. possessing both longitudinal (along the SL axis) and transverse components, then the longitudinal component of the quasi-momentum p x is affected by the Lorentz force, which oscillates with the frequency of the cyclotron rotation in the transverse plane ω c (being proportional to the longitudinal component of the magnetic field), while the amplitude of the force is proportional to the transverse component of the field. Given that it is the motion of p x which determines the Bloch oscillations, Bass et al. 17 suggested that the drift velocity might be expected to undergo sharp changes as parameters (e.g. the electric field F ) approached the vicinity of their "resonance"values corresponding to the resonance between the original (i.e. when only the electric field is present) Bloch oscillations and the cyclotron rotation i.e. ω B ≡ 2π/t B = nω c with n = 1, 2, 3, . . .. The authors did not verify their idea by numerical calculations. Rather they solved the problem analytically in the asymptotic limit of small inclination angle, and then extended the conclusions to the general case. The resultant expression for v d (F ) is a sum of the ET peak and of resonance contributions in the vicinity of values F n corresponding to the above equalities, and the shape of any of these contributions with a given n is an odd function of F − F n , while its magnitude decreases to zero as temperature goes to zero. The latter two features are in striking disagreement with numerical calculations 19,21,30,37 , our analytic formulas (see 37 and the present paper) and experiments 21 . As shown in Sec. VII.A below, the disagreement evidently originates in neglect of the feedback from the Bloch oscillations to the cyclotron oscillations, although the feedback plays a crucial role despite being weak 37 .
The next milestone in the subject was the theoretical work by Fromhold et al. 19 , who were evidently unaware of the pioneering paper by Bass et al. 17 . In one respect, their consideration was narrower as it related only to the zero-temperature limit but, in another respect, it was broader as they did not restrict consideration to small angles. Generalizing the Esaki-Tsu approach 2 , Fromhold et al. 19 calculated v d (F ) for a moderate angle numerically rather than analytically: contrary to the results of 17 for the zero-temperature limit, they did find strong resonant peaks and their shape clearly suggested that the resonant contribution for a given n was an even (and positive) function of F − F n . Furthermore, the authors had found a certain integral of motion for the intercollisional semiclassical dynamics which allowed them to reduce this rather complicated 3D dynamics to the relatively simple dynamics of a 1D classical harmonic oscillator subjected to a travelling wave, where the oscillator corresponds to the cyclotron rotation and the wave corresponds to the feedback from the Bloch oscillations to the cyclotron rotation, and the amplitude of the wave is proportional to the squared transverse component of the magnetic field. As is well known from the theory of dynamical systems, the phase plane of a harmonic oscillator subject to a travelling wave is threaded by a stochastic web (SW) when the ratio between the wave and oscillator frequencies is integer 43 or almost integer 44 . This web represents so-called weak (or non-KAM) Hamiltonian chaos and it plays an important role in many physical systems 45,46 . Fromhold et al. 19 suggested that the dynamical origin of the resonant peaks lay in a delocalization of electrons through chaotic diffusion along the web. Note that chaos in SLs had been predicted since 1995, starting from the theoretical work of Bulashenko and Bonilla 81 , followed by many experimental and theoretical works (references to those before 2005 are given by Bonilla and Grahn 4 , while some more recent ones are given by Li et al. 82 and by Alvaro et al. 83 ). However they related to a different kind of chaos -spatiotemporal chaos that occurs in a different kind of SL -with a very large number of weakly interacting quantum wells subjected to dc and ac electric fields 4,81 . In contrast, Fromhold et al. 19 referred for the first time to low-dimensional Hamiltonian chaos. Most strikingly, this chaos was claimed to be responsible for an enhancement of the unidirectional average shift of electron drift, thus playing a constructive role.
The next milestone was the initial experimental evidence for the main resonant peak. It was obtained in 2004 when Fromhold et al. 21 reported observation of a resonant peak in the differential dc conductivity, which can be shown 3 to be a direct consequence of the peak in v d (F ). This achievement had required the design of a superlattice with particular characteristics 20 . The authors also developed the theory in two important respects. First, there were both elastic and inelastic scatterings in their experimental sample so that calculation of v d immediately via the aforementioned integral representation 2,19 (sometimes known as the kinetic formula 19 ) was impossible. To overcome this difficulty the authors generalized the kinetic formula in such a way that the result of its application in the absence of a magnetic field coincided with the explicit result obtained by a different method by Ignatov et al. 84 . It was not rigorous, but the theoretical calculus based on the generalized kinetic formula exhibited reasonable qualitative agreement with the experimental results 21 , giving grounds for hope that, even when elastic scattering is dominant over inelastic, such a heuristic generalization could be used at least for qualitative predictions. Secondly, the authors assumed the existence of electron heating 3 and took account of it (their criterion for choice of electron temperature can be found in 85,86 ). Finally, the authors considered the collective dynamics by means of self-consistent drift-diffusion calculations 3 of dc current-voltage characteristics (I-V). They assumed that, as in the absence of the magnetic field 3 , the concentration of electrons changes only in the direction of the SL axis, and solved self-consistently the system of two equations: (i) the current continuity equation, in which the current is assumed proportional to the product of the single-electron drift velocity v d (F (x)) and the concentration of electrons n(x), and (ii) the Pois-son equation with the proper boundary conditions. Then they averaged the resulting current over time, thus obtaining the dc current-voltage characteristics I(V ), and calculated their derivatives with respect to V , yielding the differential dc conductivities G dc (V ). Their results are in reasonable qualitative agreement with the major experimental features.
These theoretical and experimental results 19,21 and, in particular, the dramatic conclusions drawn about the chaotic origin of the resonant peaks, further stimulated interest in the subject. It led to many new theoretical and experimental investigations some of which suggested the potential for applications.
Greenaway et al. 26 followed a similar theoretical approach to that used in 21 but, unlike the latter, they did not do the time-averaging and, moreover, focused on the current oscillations that arise for sufficiently large voltage. Such oscillations occur in the absence of the magnetic field too 3,4 , resulting from an instability associated with the range of F where dv d (F )/dF < 0. But the application of a strong, distinctly inclined, magnetic field can substantially change v d (F ), adding pronounced resonant peaks at values of F much higher than that at which the ET peak has its maximum. This suggests that the current oscillations at high voltages should have much larger frequency and intensity than in the absence of the magnetic field. Numerical results by Greenaway et al. 26 seemed to confirm the validity of these ideas, apparently implying the possibility of using SLs for the generation of electrical signals in the sub-THz range.
Because the resonant transport mechanism was believed to be chaotic diffusion through the SW, it was natural to search for easy ways of enhancing such diffusion. One such way was suggested by Soskin et al. 27 and further developed in 29 : (i) it was observed in numerical simulations that, regardless of the perturbation amplitude, the radius of the SW was necessarily limited; and (ii) it was shown analytically and verified in numerical simulations that both the radius of the SW and the diffusion rate may be greatly increased by the addition to the dc electric field of even a weak low frequency ac component. This appeared to promise a very convenient tool for the control of resonant transport.
The next important contribution was the theoretical paper by Sel'skii et al. 30 , studying the evolution of resonant transport with rising electron temperature. Their numerical results for a typical example show that the main resonant peak (the 1st one, corresponding to n = 1) in v d (F ) smears. This demonstrated once again the incorrectness of the analytic results by Bass et al. 17 which predicted the absence of a resonant contribution at zero temperature and its growth with rising temperature. Apart from that, it was numerically found in 30 that small higher-order resonances may appear, grow and/or become more distinct as the temperature increases in some range. The numerical results of Sel'skii et al. 30 showing the evolution of v d (F ) with temperature, especially those concerning the 1st resonance peak 87 , are valuable, but their suggested explanation is controversial: on the one hand, they base their arguments on the conjecture about the chaotic origin of the resonant transport but, on the other hand, they refer to the analytic result by Bass et al. 17 (which they re-derive in detail, making the same error as Bass et al. 17 ) although: (i) the latter assumes a purely regular mechanism; and (ii) the analytic result contradicts the numerical results, especially those related to the main peak. Furthermore the conjecture about its chaotic origin suggests that magnitude of the main peak grows as temperature increases whereas the numerical results demonstrate the opposite evolution.
The next milestone was the experimental work by Alexeeva et al. 32 , aiming to verify the promising predictions by Greenaway et al. 26 discussed above. Some predictions were confirmed qualitatively: the amplitude of the oscillations increases substantially if the inclination angle and the voltage enter ranges corresponding to the 1st (or even 2nd) resonance. At the same time, most comparisons between theory and experiment turned out to be disappointing: (i) the theory predicted a large oscillation amplitude about the dc value, while the experiment gave values ∼50 times smaller; (ii) the theory predicted that the frequency of the most powerful harmonic of the oscillations would be ∼10 GHz in the absence of the magnetic field and would rise to ∼100 GHz in the presence of the optimally inclined magnetic field (due to the doubling of the main frequency and the strong redistribution of the power in favor of high harmonics) while experiment showed that, regardless of the presence of the magnetic field, the most powerful harmonic was the first one and it was always at ∼1 GHz -thus differing from the theoretical prediction by a factor of up to 100×; (iii) the theory predicted that the oscillations would exist over a wide range of voltage V while the experiment demonstrated that oscillations were restricted to a rather narrow window in V . This disagreement appeared to indicate that the theoretical model had missed some important features, and the authors proposed possible explanations. They suggested that the discrepancy in the oscillation amplitude might be due to electron-electron scattering coming into play as the electron density increases, and that this type of scattering reduces the rise in electron density, while the theory does not take this process into account. Hence, within the theory 26 , the electron density can grow to very high values and it is these moving domains of high density that might give rise to the large amplitude current oscillations. To account for the discrepancies in the oscillation frequency and in the relevant voltage range, Alexeeva et al. 32 suggest the introduction into the model of an auxiliary resonant circuit consisting of a parasitic contact capacitance and an effective inductance and resistance. They managed to match the parameters of the circuit so that the calculated fundamental frequency remained at ∼1 GHz for all voltages and inclination angles while the relevant voltage range became limited from above. The hypothesis about parasitic impedance effects was, of course, rather speculative since it was not tested by independent measurements. The hypothesis concerning electron-electron scattering is similarly speculative. These issues require further study, both theoretical and experimental. Once the discrepancies between the experiments 32 and theory 26 have become well understood, it might be possible to improve the device so that its proposed role as a source of electrical signals in the sub-THz range becomes feasible.
Our Letter 37 in 2015 suggested a different direction for research on the subject. Until then, starting from the key work by Fromhold et al. 19 in 2001, the resonant drift enhancement had universally been assumed to originate in chaotic diffusion [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36] . Moreover, in papers on other subjects 47-69 , (references after 2015 will be discussed separately, below), authors referred to the phenomenon as being a characteristic manifestation of dynamical chaos in quantum electron transport. We realized 88 in 2011, however, that this widespread belief in the chaotic origin of the phenomenon (shared by ourselves up to then) was both surprising and unjustified, because there were many reasons for doubt.
First, the inference of the chaotic origin of enhanced transport 19,21 was based on arguments that were more intuitive than rigorous. Secondly, for typical parameter values used in the experiments and numerical simulations, the chaotic layers of the SW are narrow, so that the diffusion is bound to be slow. This is the case even at the centre of the web. So, how could such a slow diffusion cause a strong enhancement of the drift velocity? Thirdly, one would expect chaos to make all directions approximately equally probable. So, why should its onset lead to the observed unidirectional drift enhancement? A fourth reason for doubt is that chaos relates to collisionless electron motion, whereas drift is inherently impossible without collisions/scattering [2][3][4][5] . Moreover, the timescale on which this collisionless chaos would become pronounced for typical parameter values greatly exceeds the average scattering time. So, how could an electron's behaviour be strongly affected by something that typically cannot manifest itself prior to the scattering? Finally, if the resonant drift originated in chaos, then the strongest drift would take place when the inclination angle of the magnetic field θ approached π/2: chaos is maximized at this angle. But the experiments 21 demonstrate the opposite: the resonant drift vanishes as θ approaches π/2.
It was these questions about the dynamical mechanism underlying the phenomenon that triggered our own study of the issue, ultimately resulting in the publication of our Letter 37 in 2015. This resolved the problem in the zero-temperature limit 19 where the drift enhancement is maximal 30 . Based on regular approximations of the exact collisionless equations of motion, an analytic solution of the problem was obtained that agreed well with a solution based on numerical integration of the exact equations. Thus the resonant peak was successfully accounted for without any need to invoke chaos. The mechanism underlying the phenomenon was also explained in qualitative terms. It was demonstrated 70 that our work allows us very easily to find the values of the physical parameters needed to maximize the drift, and that they differ markedly from those required to optimise chaotic diffusion, if this were the operative mechanism 29 . The possible effect of the diffusion was also considered. It was demonstrated that, dependent on the parameter values chosen, chaos is either absent or too weak to be relevant; or, when it is strong on the relevant time-scales, it suppresses the electron drift rather than enhancing it.
The most recent milestone in the evolution of the subject was the theoretical work by Bonilla et al. 40 . Its primary aim was to describe the experimental results on current oscillations by Alexeeva et al. 32 without need to introduce an auxiliary parasitic resonant circuit. Moreover, it presents, in a sense, a challenge to all previous theoretical works on the subject. We therefore pay special attention to an analysis of this paper. The authors claim to have shown self-consistently that single-electron dynamics in the 3D space necessarily causes a redistribution of the electron density in space that gives rise to the appearance of a strong component of electric field in the z direction ( Fig. 1), so that the collective dynamics is 2-dimensional, rather than 1-dimensional, contrary to the earlier assumption 21,26,28,30,32 . Their numerical results for parameters similar to those used in 32 appear to capture qualitatively some of the characteristic features of the experimental current oscillations 32 , without artificial introduction of the resonant circuit: the frequency of the oscillations is independent of magnetic field while being ∼ 1GHz and the range of V for which the oscillations occur is limited from above. At the same time, some other features still strongly disagree, e.g. the amplitude of oscillations still exceeds that in the experiment by two orders of magnitude. We have carried out a thorough analysis of the paper 40 and summarise our conclusions below.
Because the resonant enhancement of the dc current studied previously 19,21,37 and in the present paper relies on the time-independence of the electric field, we are interested in whether the general claim by Bonilla et al. 40 about the generation of the strong electric field in the z direction is relevant to the resonant dc component of the electric current through the SL. Generally speaking, it might indeed be relevant, but all experiments to date were carried out on samples for which the longitudinal dimension was much smaller than the transverse dimension: L x ≪ L y ∼ L z ≡ L, where L x , L y and L z are the length scales of the sample in the direction of x, y and z respectively (note that Fig. 1 is schematic in this respect for, in reality, L/L x > 100 always 21,32 ). Let us show that this strong inequality provides for the y and z components of the dc electric field to have much smaller values than the x component almost over the whole SL. In experiments, the electric field in the SL is generated by the application of a voltage difference V between the contacts (collector and emitter) on the SL's right and left boundaries respectively. For the sake of simplicity, we assume below that the transverse crosssection of the SL is of rectangular form (|y| ≤ L y /2, |z| ≤ L z /2) and that the longitudinal dimensions of the emitter and collector are negligible compared to L x 89 . This means that the electric potential W at the contacts is fixed: Hence, the absolute value of the x-component of the electric field averaged over x is uniform within the y − z area of the SL and is equal to |F x |(y, z) = V /L x . The local electric field is determined by the gradient of the electric potential at a given point: F = ∇W , i.e. F x = ∂W/∂x, F y = ∂W/∂y, F z = ∂W/∂z. If we assume that the z component of the electric field does not change sign as z varies while its absolute value |F z | is ∼ |F x |, then, integrating the definition F z = ∂W/∂z and allowing for the strong inequality L z /L x ≫ 1, we would conclude that, typically, but the latter strong inequality contradicts the initial assumption |F z | ∼ |F x | thereby proving its invalidity. The huge values of |F x | in the above estimate could be avoided if we assume that the sign of F z alternates as z changes. However such an assumption does not make sense: it would inevitably lead to a non-stationarity (cf. 40 ) while we are interested in the stationary situation. F y may be considered analogously. So, we conclude that the dc electric current through the SL may be considered while neglecting small transverse components of the electric field generated by possible small inhomogeneities of the space electron density and the Hall effect 70 . It is worth mentioning also an important assumption used by Bonilla et al. 40 in solving their kinetic equation: it is that "Bloch, cyclotron and collision frequencies are of the same order" 40 , while the pronounced resonant increase of the drift velocity can occur only if the collision ("scattering") frequency is much smaller than two others 37 (see also below and cf. 19,21,30 ).
At the same time, when a non-stationarity comes into play, the assumption about the appearance of non-small transverse components of the electric field oscillating in space is no longer inconsistent. Therefore, the ideas and methods introduced in 40 may be very fruitful for the theoretical description of oscillations of the electric current. The results of 40 seem to be arguable since some of the assumptions used in their derivations raise questions. For example, the authors assume without explanation that F y ≡ 0, which is not obvious at all. Furthermore, one of the main assumptions which they use while solving the kinetic equation is, as mentioned above, the compa-rability of the Bloch, cyclotron and collision frequencies. But then they extrapolate the results to ranges of strong magnetic field and large voltage, where the cyclotron and Bloch frequencies greatly exceed the collision frequency and one of their main conclusions relates just to the latter limit. This appears to be inconsistent. However, the work by Bonilla et al. 40 is interesting, stimulating and serves to demonstrate that the subject is currently far from being fully understood. Moreover, the authors suggest a scheme for an experiment that could test their results and ideas in a straightforward manner. We note that two relatively recent theoretical works 41,42 studying the oscillations and synchronization in related systems still use the 1D model of the collective transport and do not even mention 40 . We hope that this brief review of 40 will draw the attention of other researchers to this stimulating work.
We conclude this brief review as follows. Despite 40 years of subject development, many important issues are still being debated and many fundamental problems remain unresolved. Special experiments should be carried out which would test in a straightforward way the assumptions of Alexeeva et al. 32 which allowed the authors to account for discrepancies between their experimental results on current oscillations at high electric fields and the theoretical predictions of Greenaway et al. 26 . The experiments suggested by Bonilla et al. 40 aiming to test their predictions in a straightforward way also would be interesting. The development of a rigorous framework for collective transport in the case when the cyclotron and Bloch frequencies greatly exceed scattering rates, especially when the elastic scattering dominates over the inelastic one, is the most challenging theoretical problem. An explanation of the decay of the main resonant peak as temperature increases 30 is required. The latter problem will be tackled in the present paper while other problems indicated above and those identified within the discussion in the end of the paper will remain open.

III. BASIC CONCEPTS AND EQUATIONS
A. Miniband transport and resonant-time approximation.
We consider an electric field F applied along the SL axis together with a tilted magnetic field B, and we choose coordinate axes as illustrated: the x-axis coincides with the chosen SL axis, in turn, opposite to F so that the latter is F = (−F, 0, 0) where F > 0; the y-axis is perpendicular to the plane formed by the x-axis and B; the direction of the z-axis is chosen so that B can be written as B = (B cos(θ), 0, B sin(θ)) where θ is the angle between B and the x-axis.
In between the scattering events, an electron's motion is described by the semiclassical equations 1-3,5,19,24,28-30,37 : where e is the absolute value of the electronic charge.
The questions then arise as to how the scattering affects the drift velocity, and of how to describe this effect theoretically? When inelastic scattering (energy relaxation) dominates over elastic scattering (momentum relaxation), the relaxation-time approximation 2,3 is adequate. It assumes that each scattering event results in an immediate return of the electron's statistical momentum distribution to equilibrium, and that the scattering is characterized by just a single rate ν ≡ 1/t i , where t i is the average inelastic scattering time. Consider a given instant of time t (0) . The probability density for the last prior scattering to occur at t (0) − t with a given positive t is P (t) = ν exp(−νt) 2,3 . The semiclassical velocity of an electron that was scattered for the last time at the instant (2). The drift velocity results from averaging v x (t) (2) over all positive values of t with the weight P (t) 2,3 and over the equilibrium distribution of initial momenta in the dynamics (1) 3 . As will be shown in Sec. VI below (see also the corresponding numerical results in 30 ), the strongest drift enhancement occurs at temperatures approaching zero, when only zero initial momenta are relevant 3,28,30,37 , Thus the drift velocity in the zero-temperature limit is where v x (t) is given by Eq. (2) in which p x (t) is determined by the dynamical equations (1) with the initial conditions (3).
where ω B ≡ edF/h is the Bloch frequency, and an explicit integration in Eq. (4) gives the classical Esaki-Tsu (ET) result 2 : The functionṽ ET (ω B /ν) (5) has a maximum at ω B = ν. We shall call the functionṽ ET (x) the Esaki-Tsu peak.
If B = 0, the dynamics (1)- (3) is more complicated because the components of p are interwoven. It is this dynamics that leads to drift enhancement if the Bloch oscillations and the transverse cyclotron rotation are resonant with each other [19][20][21][22][23][24][25][26][28][29][30][31][32]37 and this will be a central topic in what follows. Before considering it, we comment on the ways in which elastic scattering may be taken into account. Generally speaking, this is a very complicated problem. It is simplified when both B = 0 and elastic scattering does not affect the transverse components of momentum. It can be shown 84 that the drift velocity is then described by an equation analogous to the Esaki-Tsu formula (5) in which the scattering constant and the maximum velocity are modified as follows: where t e is an elastic scattering time.
The equation for the drift velocity with the modified scattering constant and maximum velocity can be presented in a form similar to (4): Equation (7) (with ν (mod) and µ given in (6) and v x (t) given by Eq. (2) with (1) and (3)) was also used by Fromhold et al. 21 for the case B = 0. The same was true for most subsequent researchers using the semiclassical model 20,21,23,25,26,[28][29][30][31][32]37 . This generalization was not substantiated, however. A rigorous analytic representation of the drift velocity using the semiclassical trajectories and scattering constants is a very complicated problem which has yet to be solved. Surprisingly however, theoretical calculations for the current-voltage characteristics and for the differential conductivity in some particular SLs (where elastic scattering is significant) using Eq. (7) agree reasonably well with the experimental data, at least qualitatively 21 . Hence the use of (7) for the case when B = 0 and elastic scattering cannot be neglected may be considered as a useful heuristic approach. We will also adopt it in the present work, leaving the development of a rigorous treatment as a challenge for the future. In the case where elastic scattering is negligible, Eq (7) reduces to (4) and thus becomes rigorous.
B. Semiclassical dynamics: reduction to a classical oscillator driven by a travelling wave.
For a better understanding of the semiclassical dynamics (1)-(2), it is convenient to present Eq. (1) in a more explicit form -as a system of dynamical equations for the components of p 28,30,70 - , ω ≡ eB m * cos(θ). Formally, the quantities ω ⊥ and ω represent the frequencies of an autonomous classical cyclotron rotation under the action of magnetic fields equal to the magnetic field components perpendicular and parallel to the SL axis, respectively. If there was no motion along the SL, i.e. if v x was equal to zero, the 2nd and 3rd equations would merely correspond to autonomous cyclotron rotation in the transverse plane (i.e. y − z). On the other hand, such an autonomous rotation is possible only if the momentum in the transverse plane is non-zero, which is not in fact the case at the initial instant if temperature is zero: (3)). But even then, as time goes by, electron motion along the SL does occur which, given the presence of the z component of the magnetic field, results in a Lorentz force causing the electron to move in the y direction and, in turn, triggering cyclotron rotation in the transverse plane caused by the x-component of the magnetic field. At the same time, the onset of motion in the y direction plays another important role: together with the z component of the magnetic field, it generates a Lorentz force in the x direction which, in turn, changes p x (see the 2nd term on the r.h.s of the 1st equation of (8)). The latter alters the instantaneous velocity v x (t) by changing its angle (i.e. the argument of the sine in the definition ofṽ x in (8)). Thus, the complicated form of Bloch "oscillation" 90 dynamics results from the specific interaction between it and the cyclotron rotations.
To complete this physical picture of motion in the system, we note that those variations of p x and p z in time which relate to the magnetic field are mutually correlated because they are caused respectively by the x and z components of one and the same Lorentz force (generated by the magnetic field and electron motion in the y direction): these components therefore vary in time coherently (cf. the 2nd term on the r.h.s of the 1st equation with the r.h.s of the 3rd equation in (8)). That is why p x possesses, not only the conventional component proportional to t (which is what leads to the Bloch oscillations), but also a component proportional to p z . Substituting the expression for p x in terms of t and p z into the argument of the sine inṽ x , and using the second of the dynamical equations (8), we see that the dynamics of p z and p y reduces to cyclotron rotation which is, in addition, subject to a wave at the Bloch frequency 19,21 . Thus, remarkably, the complicated three-dimensional dynamics (8) reduces to the dynamics of p z in the form of a harmonic oscillator driven by a travelling wave, while p x and p y are expressed in terms of p z . It is convenient to present the reduced dynamics in terms of scaled quantities 30,37 where ǫ is the amplitude of the Lorentz force in dimensionless units. Its significance for the resonant drift will become clear in the theory presented below. Two other scaled components of the momentum are related top z (t) ≡p(t) as follows: It is convenient to present the problem of finding the drift velocity (given in the zero-temperature limit by Eqs. (1)- (7)) in a scaled form. It is this form that will be analysed in most of the rest of the present paper. For the sake of definiteness, we assume that θ < π/2, but we emphasize that it is not essential because it can readily be shown that the drift velocity is invariant to the replacement of while the invariance to a change of sign is obvious from a physical point of view. The scaled drift velocity in the zero-temperature limit is (the case of non-zero tempratures is given in wherep ≡p(t) is a solution of the differential equation for the initial conditions The parameterν is the scattering rate in terms of the dimensionless "time"t (9).
C. Stochastic web.
Consider a harmonic oscillator subject to an alternating force. The motion of the oscillator is 91 a superposition of natural oscillations (at the oscillator's natural frequency ω 0 ) and constrained vibrations (at the frequency ω of the alternating force). The closer ω is to ω 0 , the larger the amplitude of the constrained vibrations becomes. If the resonance is exact i.e. ω = ω 0 , then the amplitude grows linearly with time at large time-scales, and so that the oscillator energy diverges. Energy is pumped efficiently into the oscillator because the phase difference between the constraint vibration and the alternating force remains constant, equal to −π/2: on timescales exceeding the vibration period, this provides for a continuous increase in oscillator energy.
If, instead of the alternating force, the oscillator is perturbed by a travelling wave like that in (13), the difference between the angles can no longer remain constant because the oscillator coordinate (p in case of Eq. (13)), which enters the phase of the wave, changes with time. If the wave frequency ω is exactly equal to ω 0 ), then the amplitude of the oscillations grows especially fast -and it is natural to expect that, as soon as it becomes sufficiently large, the average shift between the wave and oscillator angles will change so that the change in oscillator energy alters from increasing to decreasing.
These intuitive arguments turn out to be true over the most of the phase plane, provided that the perturbation is weak. However, as discovered by Chernikov et al. 43 ,44 at the end of the 1980s, there is a layer in the phase space (in the Poincaré section) in which the trajectory can travel relatively far and, correspondingly, the variation of energy can be relatively large. Similar, but much narrower, layers are formed in the case of multiple wave frequencies i.e. when ω = nω 0 with n = 2, 3, 4, . . .. The properties of the layer are highly non-trivial. In particular, the dynamics within it manifests a distinct stochasticity at large time-scales. The shape of the skeleton of the layer is reminiscent of a cobweb, with the number of rays equal to 2n: see the example in Fig. 2. That is why such layers were named stochastic webs [43][44][45][46] .
It was postulated in 19 that it is chaotic diffusion along such a web that leads to the drift enhancement: "chaotic dynamics delocalizes the electron orbits and this increases the drift velocity and conductivity" 19 . Although this idea was not supported by any analytic estimates (thus being purely heuristic), it was assumed to be cor-rect in numerous further works on the subject (including our own) up until the Letter by Soskin et al. 37 .

IV. ASYMPTOTIC THEORY OF THE RESONANCE DRIFT FOR THE ZERO-TEMPERATURE LIMIT
In this section, we will studyṽ d (12)- (14) vs. ω which is equal to the ratio of the Bloch and cyclotron frequencies (see (9)) thus being proportional to the electric field F . We will show thatṽ d (ω) possesses the "resonance" peak at ω ≈ 1 and that it may be of magnitude ∼ 1 for arbitrarily small ǫ (while there is no distinct peak if ǫ is moderate or large). In contrast, the resonance contributions near multiple and rational frequencies vanish in the limit ǫ → 0. So, they are ignored in our asymptotic theory.
A. Necessary conditions and the key parameter α.
Necessary (but not sufficient) conditions for the distinct resonance peak are: If either of these conditions fails, the resonant component of v x (t) cannot accumulate for long. Besides, if the second condition fails, the peaks at multiple/rational frequencies are significant and/or the dynamics at the relevant time-scales is chaotic as shown in Sec. V below. We assume further that the conditions (15) hold true unless otherwise specified. i.e. the analytic theory presented below is asymptotic over the two small parameters in (15). We stress however that their ratio, may be arbitrary. Remarkably, this parameter alone is shown below to affect the magnitude and appropriately scaled shape of the resonance component ofṽ d (ω). The parameter α represents the ratio of the scattering timescale and the timescale of the strong modulation of the Bloch oscillation angle -which in terms of dimensionless time (9) are respectivelyt s =ν −1 andt SM = 4/ǫ 92 . To illustrate the latter timescale, consider the exact resonance ω B = ω c . The modulation amplitude A am then grows linearly with time, as A am = ǫt/2, until A am ∼ 1. The latter range is reached just byt ∼t SM , and so strong modulation changes essentially the dynamics (13). However, if α ≪ 1, then scattering occurs before the modulation can become strong so that the latter is irrelevant. Otherwise the strong modulation does come into play, and the resonant drift occurs differently.
B. Small α limit. Mechanism of the resonant enhancement of the drift.
We consider first the limit α ≪ 1. Not only does it allow us to obtainṽ d (ω) in explicit form but, even more importantly, it clearly illustrates the mechanism of drift enhancement.
In this case, the magnitude ofp at the scattering timescalet s ≡ν −1 is ∼ α ≪ 1, so that we can neglectp in sin(ωt −p) on the r.h.s of the equation of motion (13), which then reduces to the equation of the constrained vibration. For ω = 1, its solution with the initial conditions (14) can be presented in the following form (which can be checked by direct substitution into the reduced equation (13) and into Eq. (14)): For ω = 1, i.e. the case of exact resonance, the solution can be found either as the asymptotic limit of (17) for ω → 1, or independently. It reads as It follows from (12) that the relevant time-scale is the scattering timet s ≡ν −1 . It can be seen from (18) that, for this time-scale, |p| is less than or of the order of α ≪ 1.
Allowing for the smallness of |p|, we retain in the Taylor expansion ofṽ x overp only the terms of the zeroth and first orders: Substitutingp (17) intoṽ x (19) and then substituting the result into the integral in (12), performing the integration, doing some simple but cumbersome algebra, and omitting terms of the order of ǫν 2 , we derive: This is a superposition of the Esaki-Tsu (ET) peak (5) and the resonance peakṽ (ω). The approximation of the exact resonance peak byṽ (res) d1 (ω) is valid up to lowest order in α and up to first order inν. If we neglect the first-order corrections inν (vanishing in the asymptotic limitν → 0), then the peak reduces to a Lorentzian of half-widthν and maximum α, acquired at ω = 1: The physical mechanism of the peak formation is as follows.
If ω B is close to ω , then the instantaneous velocity is a sum of fast-oscillating terms (oscillating with frequencies close to integer values of ω ) and the slow term. Let us demonstrate it for the most distinct case -when the resonance is exact: ω B = ω (which is equivalent to ω = 1). Substituting (18) After the integration (12), each of the three first ("fastoscillating") terms inṽ x results in a contribution toṽ d which is proportional toν and to the amplitude of a fast-oscillating term. The term sin(t) has the largest amplitude and therefore, of the contributions resulting from the fast-oscillating terms inṽ x , we need retain only this one: it gives the first term on the r.h.s of (20) and can be interpreted as the Esaki-Tsu contribution. The resonant contribution toṽ d originates in the non-oscillating term inṽ x i.e. in the last term on the r.h.s of (22): it accumulates during the whole of the relevant (scattering) time. Moreover, the non-oscillating term grows linearly in time, thus making the accumulation particularly strong.
In other words, the angle of the instantaneous velocity v x (t) (i.e. angle of the Bloch "oscillation") is modulated by the cyclotron rotation in the transverse plane and, if the modulation frequency ω coincides with the Bloch frequency ω B , then the modulation-induced deviation of v x (t) possesses a component that retains its positive sign and, moreover, grows with time. The drift velocity v d is defined as the instantaneous velocity averaged over a statistical ensemble of electrons, while the statistics in the zero-temperature limit under consideration is gained only due to the random nature of the scat- the probability density for the interval between the given instant and the last preceding scattering to be equal 2,3 to t . The preservation of the (positive) sign of the slow component in the resonant modulation-induced deviation of v x (t), as well as the growth of its value with t, leads to a large contribution to the drift velocity as compared with fast oscillating components of a similar magnitude.
If ω B −ω = 0, the slow component in the deviation oscillates with the frequency |ω B − ω | (being proportional to sin((ω B −ω )t)/(ω B −ω )) and therefore, if |ω B −ω | ≫ ν, its sign changes many times during the scattering timescale t s ≡ ν −1 , so that the drift velocity averages almost to zero as compared to the exact resonance case.
We will now compare (20) with numerical simulations of the model (12)- (14) for parameters typical of the SLs used in most experiments 20,21,23 and for a typical magnetic field. Let us express the dimensional parameters ǫ, ν and α explicitly in terms of the SL physical parameters: As follows from the expression for α, it quickly decreases together with θ -approximately quadratically if θ ≪ π/2. For the case (25), α may be considered as small starting from θ ≈ 25 o , so that Eqs. (20) and (21) start working from this value and their accuracy quickly improves as θ further decreases. As θ increases further, the excess of the theoretical resonant peak (20)/(21) over that in the simulations grows: v d (ω = 1) in the simulations 30,31,70 for θ = 40 o is about half that given by (20)/ (21). The invalidity of (20)/ (21) here is unsurprising because α ≈ 0.77 is not small.

Transformation to slow variables.
To encompass arbitrary α, we develop an approach suggested earlier 43,44 in a different context. If ω ≃ 1 in (13), then, neglecting small fast oscillations, the dynamics reduces to the regular dynamics described by means of the "resonant"Hamiltonian 29,43-46 : where J 1 (x) is a Bessel function of the first order 94 . If |ω − 1| is sufficiently small, the Hamiltonian (26) possesses saddles generating separatrices ( Fig. 4(a),(b)). When ω = 1, the separatrices merge into a single infinite grid ( Fig. 4(a)). For the original system (13), the neglected fast-oscillating terms dress this grid with a chaotic layer, thus forming a stochastic web (SW). Formally, chaotic diffusion along the vertical filaments of the SW might transport the system to arbitrarily high values 95 of I , so that | p| might become arbitrarily large. In the earlier work on the subject 19-32,62,67,68 preceding 37 , it was this chaotic diffusion within the collisionless approximation of electron motion that was believed to be the origin of the resonant drift. But as explained and numerically demonstrated by Soskin et al. 37 (see also the next subsubsection and Sec. V below), this cannot be the case.
The true origin of the resonant drift is explained in the preceding subsection. It does not even relate to the wave-like form of the perturbation of the harmonic oscillator in the equation of motion (13): the possibility of approximating the wave by an ac force on the relevant time scale is in itself sufficient to give rise to resonant drift (in contrast, web formation necessarily requires the perturbation to be wave-like [43][44][45][46] ). At the same time, the wave-like form of the interaction between Bloch oscillations and the cyclotron rotation is essential for taking scattering into account, thus crucially affecting the drift.
We transform the equations of motion for the system (26) from I to ρ, and scale the time and frequency shift by the slow "time"t SM and its reciprocal, respectively: The initial conditions for (ρ,φ) corresponding to zero initial conditions (14) for (ṗ,p) are: The conditions (28) are derived as follows. From the definition of ρ in (26), its value corresponding to the initial conditions (14) is equal to zero. In order to avoid a singularity at the exact zero in the denominator on the r.h.s of the equation for dφ/dτ in (27), it is necessary to take an infinitesimal positive value instead, which we denote as +0. The initial angleφ is formally indefinite. However, if ρ = +0, then the equation for dφ/dτ shows that, for anyφ from the ranges ] − π/2, π/2[ and ]π/2, 3π/2[, the derivative dφ/dτ has a diverging absolute value while its sign is positive or negative respectively. Thus, the system (27) with an initial ρ equal to +0 and an initialφ lying beyond an infinitesimal vicinity of the value −π/2 is immediately transferred to an infinitesimal vicinity of the point (ρ = +0,φ = π/2), from which motion starts with finite derivatives of both variables. That is why the most convenient choice of the initial angle is π/2.

Pronounced resonant drift in the absence of the stochastic web.
Let us demonstrate that the resonant drift may be pronounced even when the SW does not exist at all (cf. Fig.  4(c)). As follows 44 from (27), the phase plane of the Hamiltonian system (26) lacks saddles if x (1) 1 is the first zero of the Bessel function of first order 94 . Therefore there are no separatrices, so that there is no chaotic diffusion associated with SWs. Despite the absence of diffusion at δ > δ (1) cr , resonant drift may still be pronounced for a broad range of such δ. This is manifested most clearly in the case of small α, where the resonant peak is described by Eq. (21). In this case the half-width of the peak of the resonant contribution is  (29), the definition of δ in (27), and the dependence (25) of ǫ on θ, that the SW for θ = 20 o exists only for |ω − 1| < ∆ω cr ≈ 0.0077. Fig. 5 relates to ω = 0.97, which lies well beyond this range. Fig. 5(a) shows several trajectories in Poincaré section: it is evident that there is no SW 96 . At the same time, it is clear from Fig. 5(b) that, for this same ω = 0.97, the resonant contribution to the overall drift is almost equal to that for ω = 1 (where the resonant contribution acquires its maximum) and, moreover, it exceeds the Esaki-Tsu contribution. Thus, it is pronounced, despite the absence of chaotic diffusion. This proves once again that the origin of the resonant drift does not lie in chaotic diffusion.
An SW does exist if ω lies within a narrow vicinity of the exact resonance ω = 1. However, for the case of α ≪ 1, the area of the Poincaré section relevant to resonant drift is situated close to the origin (see Sec. IV.B above), and thus is far from the SW saddles and circumference, i.e. far from the SW elements which can serve as sources of the characteristic chaotic diffusion [43][44][45][46] : hence the diffusion cannot affect the resonant drift in this case either. As α increases up to values ∼ 1, the area relevant to the drift widens so that the SW saddles and circumference come into play. However, rather than enhancing the drift, their involvement suppresses it as will be shown below. For α ≫ 1, the SW's influence on the drift becomes dominant and the drift then ceases.

Amplitude of the resonant peak.
It is intuitively obvious (and is proved in 70 ) that the maximum of the asymptotic resonant peak inṽ d (ω) occurs at the exact resonance i.e. for ω = 1. The purpose of the present sub-sub-section is to provide an analytic description and analysis of the maximum A.
So, let us put δ = 0. It then follows from the second of the equations of motion (27) thatφ retains its initial value π/2. Hence, the first of the equations of motion (27) is closed, and can be integrated in quadratures: Expressingp in terms of polar coordinates (26), i.e.p = ρ sin(ϕ) ≡ ρ sin(φ − π + ωt), and allowing forφ = π/2 and ω = 1, we obtain: where ρ(τ ) is the function implicitly defined by Eq. (30). Substituting the expression (31) into the equations (12) for the drift velocity with ω = 1, using the equality sin(t −p) = sin(t) cos(p) − cos(t) sin(p), expanding the functions cos(ρ cos(t)) and sin(ρ cos(t)) into Fourier series 94 , keeping only terms of lowest order in ǫ andν, and making the change of variablest → ρ at the integration, we derive (see details in Sec. 1 of 80 ): The small-α and large-α asymptotes of A(α) are: The function A(α) and its asymptotes are illustrated in Fig. 6(a). Fig. 6(b) showsṽ d (ω = 1) forν = 0.02 as a function of ǫ calculated (i) numerically using the exact equations (12)- (14), and (ii) by the asymptotic formula (32). The agreement is excellent up to ǫ ≈ 0.3 and good up to ǫ ≈ 0.7. The fluctuations and the faster average decay at larger ǫ are explained by the fact that, at such non-small ǫ, chaos comes into play 37,70 (see also the corresponding discussion in Sec. V below). The amplitude A of the resonant peak vs. α has two remarkable features. The first of these is its asymptotic universality i.e. independence of parameters. Another, even more striking, feature is its non-monotonicity: it first increases (approximately as A 0 (α) ≡ α at α ≪ 1), attains its maximum A max ≈ 0.38 at α max ≈ 1.16, and then gradually decays to zero approaching the asymptote A ∞ (α) ≈ 1.84/α at α ≫ 1. One reason for the nonmonotonicity is a suppression of the drift exerted by the SW saddle. The physical origin of this, together with another reason for the non-monotonicity, are outlined in 37 and a detailed explanation is given in Sec. 2 of 80 .

Shape of the resonant peak.
In this sub-section we provide an analytic description of the shape of the resonant peak for arbitrary α. Details can be found in 70 while we present here just results and a few illustrative figures.
We consider the case of inexact resonance, i.e. δ = 0 (but, if δ → 0, thenṽ (res) d obviously reduces to that for δ = 0). Similarly to the case of δ = 0 80 , one can show that, neglecting corrections of the order of small parameters (15),ṽ d may be presented as a superposition of the Esaki-Tsu peak and the resonant contribution: As a function of δ, the resonant contribution takes the form of a peak while, for a given value of δ, the latter depends only on α. The resonant contribution can be significant only in the vicinity of the resonance (i.e. where |ω − 1| ≡ |δαν| ≪ 1), where it is given by the following semi-explicit formula: where (ρ(τ ),φ(τ )) is the solution of the system of dynamical equations (27) with the initial conditions (28). Note that this solution depends only on the parameter δ.
In the cases of δ = 0 or δ ≫ 1, the solution can be found in quadratures or explicitly, respectively. In the general case of an arbitrary δ, it can easily be found numerically. Alternatively,ṽ where (ρ(τ ),φ(τ )) is the solution of the system of dynamical equations (27) with the initial conditions (28), while τ p is its period. In numerical calculations ofṽ the upper limit of integration (36) may be chosen as min(τ p , N α) with N ≫ 1 (in case of N α < τ p , the inaccuracy is ∼ exp(−N ), thus being practically independent of N provided it is large).
Eq. (36) (or, equivalently, Eq. (35)) gives a complete quantitative description of the resonant peak of the drift velocity in the asymptotic limitν → 0. It is worth noting also that it follows from Eq. (35) (or, equivalently, Eq. (36)) that a contribution to the drift velocity from segments of the trajectory with ρ close to x (1) 1 , i.e. to the value of ρ at the circumference of the SW skeleton, is close to zero, regardless of the value of α i.e. of the scattering. Thus, it is not only the saddle of the SW skeleton but also its circumference, that suppresses the drift.
In the asymptotic limits α → 0 and α → ∞, the expression (36) simplifies: 1. For α → 0, (36) can be shown to reduce toṽ (res) d0 2. For α → ∞, (36) can be shown to simplify as follows: Here, K(x) is a universal function which, to the best of our knowledge, had not been studied before 70 in any context. We will refer to it below as the K-form.
Its dependence on x is contained in the implicit dependence on x of the trajectory (27) maximum (apparently the modulus of the derivative diverges as |x| → 0) while its far wings decay as slowly as those of a Lorentzian i.e. ∝ x −2 (Fig. 7). It is exactly the same case as that considered by Sel'skii et al. 30 for T = 0 (the numerical curve in our Fig. 8(a) is identical to the curve for T = 0 in Fig. 4 of 30 ). The agreement between the theory and numerical simulations in the range of the resonant peak is excellent. At such values, (i) on the one hand, chaos comes into play resulting in an irregularity of the functionṽ d (ω) and its decrease on average, and (ii) on the other hand, the peaks at the multiple and rational frequencies become significant and wide, thus strongly affecting the peak at the main resonant frequency ω = 1 (for more details see 70 ). In the case of θ = 60 o , both ǫ andν become rather large, so that: the maximum of the main peak significantly decreases,ṽ d (ω) acquires a distinct irregularity, and the right wing starts to rise not far from the maximum (in contrast with the theory) since it is strongly affected by the peak at the rational frequency ω = 1.5.
3. Panel (c) shows the evolution ofṽ d (ω) in the vicinity of ω = 1, as α grows whileν = 0.02. In addition to perfect agreement for ǫ ≡ 4αν < ∼ 0.4, and reasonable agreement for higher ǫ up to 0.8, it illustrates the key features of the peak. One of them, namely the non-monotonic dependence of the peak amplitude on α, was discussed above in Sec. IV.C.3. Three other features are worth mentioning: (i) as α increases with fixedν, the shape of the resonant peak evolves from being Lorentzian at α ≪ 1 to the stretched K-form at α ≫ 1; (ii) as α increases with fixedν, the width of the peak increases; (iii) as ǫ becomes distinctly non-small, the resonant peak, apart from decaying if α becomes large, acquires an irregular structure. The features (i) and (ii) were considered above and in 70 . We discuss the feature (iii) in the next section.

V. ROLE OF CHAOS
As outlined above, even when an SW is present, chaotic diffusion along it cannot be the origin of the resonant drift. In the present section, we give further details for the most characteristic case -that of exact resonancefor which the resonant drift is maximal and the SW is at its most developed.
A. Small ǫ case: weak chaos, irrelevant to the drift.
We first analyse the relevant time-scale of chaotic diffusion for the case of ǫ ≪ 1. It is well known [43][44][45][46]97 that the SW is then exponentially narrow. Moreover, this statement is also relevant to moderately small values of ǫ: cf. Fig. 9. For the sake of definiteness, we assume that  (26) with H = 0). The motion is at first almost regular, being close to that along the SW skeleton. When the system reaches the vicinity of the skeleton saddle, however, it wanders within it in an irregular manner and, finally, moves in a random-like fashion onto one of the two filaments associated with the trajectories outgoing from the saddle. Such a dynamics is illustrated 98 in Fig. 9 and by the animation 80 . The time-scale t CH relevant to chaotic diffusion is the characteristic time between two sequential turns. In the context of the resonant drift, the time relevant to diffusion is the time of motion from the region close to the origin (p = 0,ṗ = 0) until the turn following the sojourn in the vicinity of the lower saddle: cf. Fig. 9(a) and Fig. 9(b)), or the animation 80 . This time may be estimated 45,46 as the time taken for movement along the corresponding part of the resonant trajectory (26) at the boundary of the SW chaotic layer. The energy (i.e. the value of the resonant Hamiltonian (26)) at this trajectory is equal to the half-width of the SW chaotic layer ∆E SW . This half-width has an exponentially small value 99 : The function D(ǫ) is rather complicated 97 but its most important property in the present context is that it diverges in the limit ǫ → 0. As follows from Eq. (27) at δ = 0, the exponent in the exponentially slow approach to the saddle along the filamentφ = π/2 is proportional to ǫ. Hence 45,46 the time-scale related to the chaotic diffusion t CH is proportional to ǫ −1 ln(1/∆E SW ) −→ ǫ −1 D(ǫ). In terms of the scaled time τ (27), it is τ CH ∝ D(ǫ). As for the characteristic time-scale τ SM (which provides an exponential cut-off for the contributions into the resonant drift), it is Moreover, this strong inequality between τ CH and τ SM also holds true for moderately small ǫ: it can be seen from Fig. 9, and from the animation 80 , which relate to ǫ = 0.407, thatt CH takes about 10 periods of the wave (i.e. 10 iterations of the section) whilet SM only takes about 2 periods (iterations). It follows from the strong inequality between τ CH and τ SM that the effect of chaotic diffusion on the resonant drift (regardless whether constructive or destructive) is exponentially small provided ǫ is small or moderately small. This is even more the case if the scattering timẽ t s ≡ν −1 is much smaller thant SM ≡ 4/ǫ, i.e. if α ≪ 1. Fig. 10 shows the evolution of the stroboscopic Poincaré section as ǫ increases. The sections are obtained by numerical integration of Eq. (13) with ω = 1 and the initial conditions (14), for 50000 periods (left panels) and 20 periods (right panels). The duration of the 20 periods is equal tot ≈ 2.5ν −1 withν = 0.02, exploited in Fig. 8(c) above. Therefore, for suchν, this duration is sufficient for a calculation of the drift velocity with an accuracy ∼ 10%. So, it is this short-time section which determines the main part of the resonant drift for the given value ofν.
The evolution exhibits some characteristic features. First, we see that the 20-periods section evolves quite differently from the 50000-periods one. At ǫ < ǫ cr ≈ 0.45, the trajectory seems to be regular even at very long times (this feature was observed earlier in 29 ). It is worth noting that the case exploited in Fig. 8(a) and in 30 corresponds to ǫ ≈ 0.407 which is below ǫ cr . Thus there is no way in which chaos could relate to the resonant peak: the longtime trajectory starting from rest is shown explicitly for this case in Fig. 11(a) (in magenta) and it obviously lies beyond the SW.
So, the long-time section practically does not change as ǫ varies within the range 0 < ǫ < 0.45 (panels (a)-(c)), unlike the short-time section. The pronounced change in resonant drift in this range is related immediately to the strong change in the short-time section. At ǫ = ǫ cr , a global bifurcation occurs in the long-time section: if ǫ < ǫ cr the trajectory is not absorbed by the SW while, otherwise, it is absorbed: cf. panels (c) and (d). However the short-time section is unaffected by this strong bifurcation which is why the resonant drift is not sensitive to it either: the drift is still well-described within the regular approximation. Chaos on the short time-scale manifested as a randomness in the direction of the turn of the "post-saddle" motion comes into play only at ǫ ≈ 0.8 (panel (e)). At ǫ ∼ 2 − 4 (panel (f)), chaos becomes very strong so that, starting from the latter range, the drift gradually vanishes as ǫ increases.

VI. GENERALIZATION TO ARBITRARY TEMPERATURE
A. Formal definitions and exact transformations.
In this section, we generalize our consideration of the zero-temperature limit to encompass the case of arbitrary temperature T . This generalization implies the possibility of a temperature dependence in ν and the replacement of the zero initial conditions (3) by their thermal statistical distribution, which may be approximated by a Gibbsian distribution 30 : where k B is Boltzmann's constant. Then the scaled drift velocity is the result of averaging over this statistical ensemble 30 : f ≡ f dp x (0)dp y (0)dp z (0) dp x0 dp y0 dp z0 = fh where u d may be called a scaled partial drift velocity related to a given set (p x0 , p y0 , p z0 ) while the momentum p ≡p(t) in the expression for the scaled partial instantaneous velocityṽ x (p x0 , p y0 , p z0 , ω,t) is a solution of Eq. (9) with initial conditions and a given value of p x0 ≡p x (0) (obviously, (41)- (43) and (9) reduce in case of T = 0 to (12)- (14)). Sel'skii et al. 30 numerically integrated equations of motion identical to those in Eq. (8) for the parameters of the superlattice and magnetic field used in Fig. 8(a) with a large number of initial conditions statistically weighted in the phase space in accordance with the distribution (41), then calculated v x (t) by use of an equation identical to Eq. (2) for each given set of initial conditions and afterwards averaged the result over the distribution of initial conditions. If properly done, such a procedure is equivalent to calculation of the drift velocity with Eq. (42). This study revealed a remarkable feature of the evolution: as the temperature increases, the ET peak quickly decreases, while the main (1st-order) resonant peak decays much more slowly so that, at moderate temperatures (∼ 100 − 400K), the resonant peak becomes a strongly dominating feature in the drift velocity vs. the electric field. However, no explanation was offered 30,31,36,38 , and no effort was made to study either how generic this feature is, or what characteristic types of resonant peak evolution may arise. We now address these important issues.
To facilitate the task, we transformṽ d (42) in the following way. We introduce the scaled temperaturẽ explicitly calculate the statistical integral Z, allow for the expression for ǫ, and transform from (p y0 , p z0 ) to the corresponding radial coordinates: Then the expression forṽ d (42)-(43) transforms tõ where I 0 (x) is the modified Bessel function of the zeroth order 94 .
In calculating the drift numerically, Sel'skii et al. 30 consider the case when the scattering is independent of temperature. Generally speaking, the scattering may vary with temperature, depending on the particular SL. However, since the main purpose of this section is just to demonstrate the regular-like mechanism of drift regardless of temperature, and for the sake of an immediate comparison with the results of 30 , we will assume like 30 thatν is independent of temperature: B. Asymptotic resonant theory

General expressions.
The general analytic theory is developed in Sec. 3 of 80 while, here, we just present its main results and their qualitative analysis as well as illustrations in terms of specific examples, including the example exploited in 30 .
As in the T → 0 case, the drift velocityṽ d (ω) can be approximated by a sum of two terms which we will refer to as the generalized Esaki-Tsu (ET) and resonant contributions (peaks) respectively:ṽ  30 and by us in Fig. 8(a) above. (a) Stroboscopic Poincaré sections of a few characteristic collisionless trajectories starting from a variety of initial states on the horizontal axis and then numerically integrated for 50,000 periods of the wave using the exact equations of motion (13); the magenta trajectory starts from the origin ("zero initial conditions"), which is marked by the large magenta dot and corresponds to zero energy i.e. to the only relevant initial energy for temperature T = 0. Green indicates a chaotic trajectory within the stochastic web; and the olive, red and blue dashed circles show levels of the scaled cyclotron energy I ≡ (p 2 +ṗ 2 )/2 being equal to the accordingly scaled kBT for temperatures T of 4.2K, 100K and 400K respectively. (b) Scaled partial instantaneous velocity averaged over fast oscillations ṽx f = J1(ρ) sin(φ) with (ρ,φ) following the equations of motion (27) with a few initial conditions as function of the slow time τ . The dashed magenta line corresponds to the zero initial conditions, the solid lines correspond to states with ρ0 = 1.71 (which corresponds to the cyclotron energy being equal to kBT with T = 400K) and a few values ofφ0 ≡ ϕ0 − pz0 − px0 + π marked by different colors.

Resonant peak evolution. Two characteristic types of the evolution.
Let us now turn to the resonant contributionṼ (res) d in (48)- (49). For the sake of brevity and simplicity, we discuss below mostly the range α > ∼ 1 of which α ∼ 1 is the most important part because resonant drift is stronger there than for other ranges. The evolution of the resonant peak in the case α ≪ 1 is somewhat different. It will be discussed in this section just briefly while further details are given in Sec. 3 of 80 .
Let α > ∼ 1 and, for clarity, consider only the top of the peak i.e. ω = 1. Fig. 11(a) presents the corresponding stroboscopic Poincaré sections of a few characteristic exact trajectories (13) whose parameters are the same as those used in Fig. 8(a) and by Sel'skii et al. 30 while initial states lie on the horizontal axis, namely: (i) the trajectory starting from the origin (ṗ = 0,p = 0), which is evidently non-chaotic (at least, there is no evidence of chaos during the huge integration time of 50000 periods); (ii) a chaotic trajectory within the SW, forming a layer which passes close to the origin while remaining distinctly separated from it; (iii) a few typical regular/quasi-regular trajectories inside the (two) meshes of the SW within its first circumference. If T → 0, then there is only one relevant trajectory: it starts from the origin of the section and, regardless of whether or not it is absorbed by the SW, it then goes approximately down the vertical axis (cf. the right-hand panels of Figs. 10(a-d)). In terms of the slow angle, this direction of motion corresponds tõ ϕ retaining the value π/2 and, as is clear from Eq. (48), trajectories with such values make the largest partial resonant contribution to the drift velocity. Asφ deviates, the contribution decreases and, ifφ lies in the interval ] − π, 0[, the contribution is negative. It can be seen from the Poincaré section that trajectories starting from the vicinity of the origin, namely when ρ 0 ≪ 2, stay close to the trajectory starting from the origin until they reach the vicinity of the saddle. The partial drift velocities for such initial states are therefore almost equal to that for the initial state at the origin of the section. In contrast, trajectories starting beyond the vicinity of the origin, i.e. with ρ 0 > ∼ 2, behave very differently. Many of them are close to the elliptic points, whereφ = 0, and thereforeφ remains close to zero along the entire trajectory, so that the magnitude of the partial instantaneous velocity averaged over fast oscillations ṽ x f ≡ J 1 (ρ) sin(φ) is small: cf. Fig. 11(b). Furthermore, if ρ 0 is close to x (1) 1 i.e. if the initial state lies in the vicinity of the circumference, then the partial drift velocity is close to zero because There is also a third relevant factor. The partial velocities ṽ x f undergo oscillations and, if ρ 0 > ∼ 2, their amplitudes, periods and relative shifts in angle vary strongly with the position of the initial state: cf. Fig. 11(b). These three factors give rise to a significant decrease in the overall resonant contribution to the drift velocity if a significant part of the thermal distribution of ρ 0 lies in the range ρ 0 > ∼ 2. As follows from the expression for the exponent in the integrand of the integral representingṼ than, or of the order of, the characteristic temperaturẽ while, ifT ≪T res , thenṼ res (T ) ≈Ṽ res (T = 0). A necessary condition for the existence of a pronounced resonant peak is the smallness of the parameter ǫ/4 (see Eq. (15) or the further details in Secs. IV.C.3 and IV). Therefore the scale ofT res strongly exceedsT ET ≡ 1: It is this strong separation of temperature scales that accounts for the seemingly more distinct manifestation of the resonant peak as temperature increases to moderately high values. The separation of scales has also allowed us to derive the relatively simple approximation (48). Fig. 12(a) demonstrates its effectiveness even for the example used by Sel'skii et al. 30 and in our Fig.  8(a) 100 , despite the parameter of smallness only being moderately small so that the separation of scales is itself only moderate strong:T res /T ET ≈ 5. Note that the real temperature corresponding toT ET ≡ 1 is T ET ≈ 70K, so that T res ≈ 350K. That is why, for T = 100K, the ET peak has already significantly decreased while the resonant contribution to the drift velocity is only a little lower than its zero-temperature limit.
The case of α ≪ 1 is rather different. The relevant timescale is then much smaller, being in terms of τ of the order of α rather than 1. For the angleφ to manage to get close toφ = π/2 (which is the value of the initial angle relevant to the caseT = 0) from most of the 2πrange of possible values of the initial angle, the condition ρ 0 ≪ 2 is insufficient: a stronger restriction ρ 0 ≪ 4α/π is relevant. Accordingly, ifT significantly exceeds (the meaning of the notation will become clear from further consideration), then the majority of electrons are scattered long before their angle reaches π/2 and therefore their collisionless trajectories differ strongly from that in the case ofT = 0. Intuition suggests that an averaging of the partial drift velocity over the whole 2πrange of possible values of the initial angle should lead to the resonant drift velocity averaging to zero, unlike the case of zero temperature. So, it is natural to expect that the scaleT 1/2 plays the same role asT res in case of α > ∼ 1 i.e. it marks the beginning of the drop in the resonant drift towards zero. Paradoxically, however, the drop to zero still occurs on the scaleT res while the scaleT 1/2 marks the drop ofṼ (res) d only twice, being then preserved over a wide range of temperature: This non-trivial evolution is explained in Sec. 3 of the Supplemental Material 80 . Ifν is of the same order in the cases α ≪ 1 and α > ∼ 1, the value of ǫ in the case α ≪ 1 is much smaller than that in the case α > ∼ 1, and thereforeṼ (res) d in the case of α ≪ 1 remains of the order of its zero-temperature limit till much higher temperatures as compared to the case of α > ∼ 1. This can be seen by comparing panels (a) and (b) of Fig. 12, which show the cases of θ = 40 o (for which α ≈ 0.77 while ǫ ≈ 0.407) and θ = 20 o (for which α ≈ 0.177 while ǫ ≈ 0.0766) respectively: the value ofT res differs in this cases by a factor of 5.3 and, as a consequence,Ṽ The first analytic results for the resonant peaks were obtained by Bass et al. 17 in the asymptotic limit of small tilt angle θ. Their result was reproduced in the review 18 and relatively recently re-derived in detail by Sel'skii et al. 30 . Two correct general conclusions may be drawn from this result: (a) an indication of possible resonant contributions to the drift; and (b) the fact of the multiple-resonance contributions in the limit θ → 0 being negligible compared to the contribution at the main resonance. Otherwise, however, the analytic solution 17,18,30 itself and the conclusions that one might draw from it are incorrect. The reason lies in their neglect of the feedback from the Bloch oscillation on the cyclotron rotation. The specific error in the derivation leading to the incorrect result will be identified below but, first, we demonstrate that the result is wrong.
The result is given by Eq. (7) in 17 , by Eq. (5.12) in 18 and by Eq. (10) in 30 . After proper scaling, all of these formulae reduce to one and the same form. For the sake of clarity and readers' convenience, we present it here using the notation of the present paper and, moreover, we consider only the asymptotic limit of small temperature. Allowing for the asymptotes of modified Bessel functions 94 I n (x) for x → 0 and x → ∞, taking into account that 94 I −n (x) = I n (x) andν ≪ 1, and skipping the contributions at multiple resonances (which vanish in the relevant asymptotic limits, in comparison to the 1st-order resonance), the approximationṽ where the normalized Bloch frequency ω and the small parameter ǫ are defined in Eq. (9) while the normalized scatteringν and temperatureT are defined in Eqs. (12) and (44) respectively. Eq. (55) suggests that: (i) forT = 0, there is no resonant drift at all; (ii) the drift magnitude increases asT grows; and (iii) even forT = 0, the resonant contribution to the drift at the exact resonance is zero while changing its sign as the sign of the deviation from the exact resonance changes, and attaining its maximum absolute value at |x| = 1. These major features ofṽ (B,S) res are all in striking disagreement with results based on numerical solution of the exact equations of motion. Indeed, not only is the resonant drift atT = 0 non-zero but it is at its most pronounced there because the drift decreases asT grows: see Sec. VI above and, in particular, Fig.12. Furthermore, the resonant contribution is necessarily positive regardless of the sign of the deviation and attains its maximum at the exact resonance: see 19,21,37 as well as the present paper.
To account for the discrepant result in the earlier papers 17,18,30 , we will refer to its derivation in 30 , which is where it is presented in its most detailed form. The exact equations of collisionless motion are those given in Eq.
(2) of that work, being identical to our Eq. (8). Considering the case of small angles θ -for which ω ⊥ ≪ ω -the authors neglect the term ∝ ω ⊥ on the r.h.s. of the equation forṗ y , assuming this to be justified because its amplitude is much smaller than that of the other term: see Eq. (A.2) in their Appendix. However, in doing so, they "throw out the baby with the bath-water". The omitted term represents the feedback between the Bloch oscillation and the cyclotron rotation, and it is just this feedback that gives rise in the resonant case to a mutual growth of the amplitude of the cyclotron rotation and the modulation of the Bloch oscillation. Because the authors neglect this feedback, the amplitude of cyclotron rotation in their theory never changes, thus remaining the same as at the initial instant. But, if T = 0, then the initial amplitude is equal to zero i.e. p y (t = 0) = p z (t = 0) = 0. Then, p y in the r.h.s. of their Eq. (A.1) (i.e. the equation of motion for p x , being the same as that in our Eq. (8)) is identically equal to zero i.e. there is no driving of p x related to the magnetic field and therefore there is no resonant drift at all. This scenario is simply wrong. The origins of the incorrect temperature evolution, and of the incorrect shape of the resonant contribution as a function of the shift from the exact resonance, are the same.  37 ) is subject to the numerous restrictions reviewed in detail by Soskin et al. 70 (see Sec. 6.2). We direct interested readers to this review and to the references therein while, here, we briefly discuss two issues that have not hitherto been properly considered.

Heating in the context of the relaxation-time approximation.
It is obvious that the electric field may significantly accelerate electrons so that their average energy can exceed that in the absence of the field. This phenomenon is called "heating". One may introduce a corresponding "temperature" for the electron subsystem exceeding that of the lattice. The heating effect is widely discussed in the literature both in a general context 101 and specifically in the context of SLs 3 . However, in the latter case, the concept of heating is irrelevant, or only partly relevant, to the validity of the relaxation-time approximation -which is of course one of the two cornerstones of the theory of resonant drift. The irrelevance (or partial relevance) arises because conventional heating is already inherent within the relaxation-time approximation 2,3 . Let us illustrate this for the most characteristic case -when the lattice temperature is close to zero while the inelastic scattering dominates over elastic scattering (t i ≪ t e ). Starting its motion from the zero-energy state (3), the electron is drawn by the electric field as well as influenced by the magnetic field if present: see Eq.
(1). But as soon as the scattering occurs, the electron's energy is assumed to drop instantaneously to zero 2,3 . Thus, it is not primarily the average electron energy that matters; rather, the scattering should give rise to an immediate loss of the energy gained in order for the relaxation-time approximation to be valid.
From the physical point of view, the above situation may be relevant when the electron is scattered mainly by phonons 3 . However, if the scattering is mainly due to impurities and surface roughness, then most scattering events do not give rise to a loss of energy. As a result, the energy of the electron gradually increases despite the occurrence of scattering, i.e. the electron temperature grows. To the best of our knowledge this situation, when the elastic scattering dominates i.e. t e ≪ t i , applies to all experiments on resonant drift reported to date including, for example, those by Fromhold et al. 21 and Alexeeva et al. 32 . Intuition suggests that, if the energy of the electron subsystem exceeds that of the lattice, then there must be a flow of energy from the electrons to the lattice and, the larger the temperature difference is, the larger will be the flow. When this flow equals the field-related rate of energy gain, the growth in electron temperature saturates. This saturation temperature may appropriately be used as the temperature within the generalised relaxation-time approximation 3,84 .
It is very difficult to make a reliable theoretical estimate of the saturation temperature. However, it is apparently possible to estimate it from a combination of theory and experimental data. Thus, if we accept the simplified model of elastic scattering introduced by Ktitorov et al. 102 , and the theory of the drift velocity in the absence of magnetic field by Ignatov et al. 84 based on 102 (it is reviewed by Wacker 3 ), then the electron temperature in the relevant context might be estimated from a comparison of characteristic features of experimentally measured current-voltage curves with those obtained theoretically on the base of the relevant theoretical expression (it is reproduced in the present paper as the termṼ ET in Eq. (48)) provided the elastic and inelastic scattering times t e and t i have somehow been found independently. There is also a much more general and consistent quantum transport model based on non-equilibrium Green's functions 103 which may be used for the above aim in some ranges of parameters (see 3 for review of this and related papers) but it involves heavy numerical calculations.
To the best of our knowledge, neither of the aforementioned schemes for the estimate of the relevant electron temperature was used in any of the works on resonant drift. Let us demonstrate this taking as an example the first experimental report 21 about resonant drift where it was stated without justification that the electron temperature is 100K (regardless of the value of applied electric field). The origin of the estimate was explained in theses of some of the associated PhD students: see e.g. p.47 in the thesis by Greenaway 86 . It was stated that the authors roughly estimated the electron temperature in the conventional sense for the most characteristic 104 electric field. Thus, it has nothing to do with a temperature estimate that would be relevant to the relaxation-time approximation. Note that, in a more recent major paper 32 reporting related experiments, the theory developed for the comparison with the experiment did not assume any heating at all, though both the SL and electric/magnetic fields were similar to those used earlier in 21 .
An estimate of the relevant heating in any specific case was beyond the scope of 37 and also of the present paper. Rather our main purpose has been to reveal the true mechanism of the resonant enhancement of the drift. We have shown in particular (see Sec. VI) that, even when heating is present, the mechanism remains the same, and that it is regular in character. But the extent of the drift enhancement at the main resonance (which is much stronger and more important than others [21][22][23]28 ) decays as the temperature increases.

Presentation of the drift in the form similar to that in
the absence of the magnetic field.
The approach of Esaki and Tsu 2 to the problem of the electron drift in SLs is phenomenological. A more consistent approach is to use the Boltzmann kinetic equation (BKE) for the distribution function 3 . Nevertheless it can be shown that the result based on the BKE reduces to the form (4) for inelastic scattering at low temperatures (cf. Ignatov and Romanov 105 ).
In one of the key theoretical papers on resonant drift, Fromhold et al. 19 generalized the phenomenological Esaki-Tsu approach to the case when a magnetic field is present. Such a generalization is quite reasonable when only inelastic scattering is present while the temperature is low. But when elastic scattering cannot be neglected, the problem of drift in the presence of a magnetic field becomes much more complicated. Palmier et al. 11 , Miller and Laikhtman 13 , and Cannon et al. 14 , applied different approaches to the problem for the case when the magnetic field is perpendicular to the SL axis. Resonant drift enhancement is then irrelevant. In the present context it is important to note, however, that they give no indication that the scattering may be taken into account in a simple integral form like that of the Esaki and Tsu 2 type. For a tilted magnetic field, one more dimension is involved: the motion becomes even more sophisticated and it is yet more doubtful whether the scattering can be taken into account in a simple integral form. Nevertheless, in the first experimental work on resonant drift 21 , it was decided to adapt the Esaki-Tsu formula to the case when the elastic scattering dominates by use of the modified parameters introduced by Ignatov et al. 84 for the case when only an electric field was present. Such an approach can be substantiated theoretically for the latter case (on the assumption that the elastic scattering is described within the model introduced by Ktitorov et al. 102 ). But when a magnetic field is present too, it may be considered as heuristic at best. Surprisingly, de-spite its heuristic nature, it provided good qualitative agreement with the experiments 21 . Probably because of this, and due to its simplicity, the same integral formula was used in most of the subsequent work, including the present paper (see Eqs. (6)-(7) above). It is a challenging problem for the future either to substantiate it rigorously or to find a different but rigorous way of taking the scattering into account in a simpler form than that suggested by the kinetic equation approach.
On the other hand, intuition suggests that the most general qualitative theoretical conclusions about resonant drift in the case when elastic scattering cannot be neglected should not depend on the validity of the simplified representation of the scattering, since they relate first of all to the specific collisionless dynamics and to the fact that scattering does occur somehow. Our present work suggests that these most general conclusions relate to: (i) a non-chaotic mechanism underlying resonant drift; (ii) a non-monotonic dependence of the resonant drift velocity on the magnetic field and on some other relevant parameters; (iii) a suppression of the drift by chaos within the collisionless dynamics; and (iv) non-chaotic mechanisms of a decay in the drift at the main resonance as temperature of electrons increases. To the best of our knowledge, all items except the last one conform with the experiments done to date, whereas experiments that could test the last item have yet to be undertaken (it may, however, be non-trivial to interpret their results).

VIII. CONCLUSIONS
We have unambiguously confirmed our earlier findings 37 concerning the non-chaotic nature of the resonant electron drift along the 1D nanometre-scale superlattice in the constant longitudinal electric and tilted magnetic fields, contrary to the widely accepted opinion about the chaotic diffusion along the stochastic web (SW) in the phase plane of one of the electron quasimomenta. Not only have we presented more details of our asymptotic regular theory 37 , but we have also explicitly demonstrated both via numerical simulations of the exact equations of motion and by means of analytic estimates that the pronounced resonant drift may exist in a broad range parameters where the SW does not exist at all, while, for the ranges where the SW does exist, chaos either comes into play at the time-scale greatly exceeding that relevant to the pronounced drift formation or suppresses the resonant drift otherwise.
An interaction between Bloch oscillations and cyclotron rotations in the transverse plane and in the plane perpendicular to the transverse component of the mag-netic field determines the intercollisional dynamics of the electron. It can be reduced to a transverse rotation driven by the "rotation-Bloch" wave: the latter represents a kind of feedback from the oscillations to this rotation via the second type of rotation. If the oscillations are resonant with the transverse rotation while the feedback is weak, then a positive contribution to the drift velocity caused by the feedback is accumulated on a long timescale. This is the true mechanism, and it was neglect of the weak feedback that led to the failure of the earlier attempts to describe the resonant dc drift analytically 17,30 .
Generalizing our asymptotic theory for the case of nonzero temperature, we have accounted for the earlier numerical observation 30 that an increase of electron temperature causes a decrease in resonant drift at the main resonance. Electrons with a variety of initial momenta (immediately after scattering) are involved and, because the initial momenta are non-zero, the relevant component in v x (t) caused by the resonant modulation of the angle of v x (t) typically oscillates even when the resonance is exact, unlike the case of zero initial conditions where it is sign-preserving. The amplitudes, periods and initial angles of such oscillations vary, leading to a decorrelation between them. The larger the temperature is, the larger the extent of the decorrelation and therefore the smaller the sum of such components becomes. Thus, heating of whatever nature does not change the mechanism of resonant drift but leads to decay of the drift due to the increasing decorrelation between contributions from electrons with differing initial conditions within the thermal ensemble. Our explicit and semi-explicit analysis allows us to classify characteristic scenarios of evolution of the resonant peak as temperature increases.
We have provided an extended analysis of the developments since the pioneering work by Bass et al. 17 up until the latest key work by Bonilla et al. 40 and we have indicated some of the unsolved problems in what remains a challenging subject.
For the lower saddle, the allowed directions of the turns are "to the left" and "to the right" (see Fig. 9(a) and Fig. 9(c) or the animation 80 ). For the upper saddle, the allowed directions are "down" and "up". However, the filament along which the "up" turn may pass is so extremely narrow that this turn does not occur during a reasonable integration time (Fig. 9). It becomes possible only if the amplitude of the wave is much larger (cf. Fig. 10(f)). 99 The intuitive, non-rigorous, derivations in 43,45,46 give D ∝ ǫ −1 . A mathematically rigorous calculation 97 suggests a different function D(ǫ) but it still diverges as ǫ → 0. 100 After being properly scaled, the numerically calculated curves for T = 0 in Fig. 4 of 30 and in our Fig. 12(a) are identical. However, as T increases, the results start to deviate from each other, so that the deviation in the range of the resonant peak reaches ∼ 50% at T = 400K. Our computation of the drift velocity is based on the immediate definition (42). In contrast, Sel'skii et al. 30 used a kind of Monte-Carlo method instead, motivating this by the necessity of faster computation. If properly done, the method does give correct results, but its convergence was checked in 30 only for the simple case of zero magnetic field (as was also the case in relation to verification of the analytic results of 17,18 ). The dynamics in the presence of a magnetic field is far more complicated, however, necessitating much larger statistics to ensure validity of the Monte-Carlo method. We have checked the result by application of the Monte-Carlo method for the magnetic field used in 30 and the results converge reasonably well to the result obtained by straightforward integration (42), but only if the number of initial states exceeds 10 6 i.e. is at least four times larger than that used in 30  The most characteristic electric field is apparently the field corresponding to the maximum of the first resonant peak for the case when the tilt of the magnetic field is such that this maximum of the peak attains its largest value (a tilt of about 45 − 50 o for the SLs and B considered in 21,86 ). However, taking temperature in the conventional sense, the figure of 100K is valid only for the specified range of electric fields, whereas it is being assumed 21,86 to be relevant for any electric field within what is a huge range, including in particular the resonant ranges for magnetic field tilt angles within the full range 0 − 90 o ). This approach is therefore inconsistent even for heating in the conventional sense.