{'fraction_non_alphanumeric': 0.030566254199786938, 'fraction_numerical': 0.011445273020295556, 'mean_word_length': 4.508576587421005, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Strong vibration nonlinearity in semiconductor-based nanomechanical systems

We study the effect of the electron-phonon coupling on vibrational eigenmodes of nano- and micro-mechanical systems made of semiconductors with equivalent energy valleys. We show that the coupling can lead to a strong mode nonlinearity. The mechanism is the lifting of the valley degeneracy by the strain. The redistribution of the electrons between the valleys is controlled by a large ratio of the electron-phonon coupling constant to the electron chemical potential or temperature. We find the quartic in the strain terms in the electron free energy, which determine the amplitude dependence of the mode frequencies. This dependence is calculated for silicon micro-systems. It is significantly different for different modes and the crystal orientation, and can vary nonmonotonously with the electron density and temperature.


I. INTRODUCTION
The electron-phonon coupling strongly affects vibrational modes of nano-and micro-electro-mechanical systems. Much interest have attracted the effects of this coupling related to the reduced dimensionality of the electron system, as they make it possible to reveal interesting consequences of the electron correlations at the nanoscale, the Coulomb blockade being a simple example, cf. [1][2][3][4][5][6][7][8][9] and references therein.
Much less attention has been paid to the consequences of the electron-phonon coupling, which are related to the discreteness of the vibrational spectrum of a nanosystem, but emerge in the absence of size quantization of the electron motion. One of such consequences, which we study in this paper, is the coupling-induced change of the vibration nonlinearity. Strong nonlinearity is a generic feature of vibrations in small systems [10,11]. Its easily accessible manifestation is the dependence of the mode frequencies on the vibration amplitudes. This dependence corresponds to the "self-action" of the mode, and its familiar analog in bulk crystals are acoustic solitons [12,13]; however, the nonlinearity required for observing such solitons usually is sufficiently strong only for high-frequency phonons. Also, the change of the eigenfrequency with the mode amplitude is of interest for modes with a discrete frequency spectrum, such as standing waves in mesoscopic systems, but not for propagating waves with a quasi-continuous spectrum.
Much attention have been recently attracting Si-based nano-and micromechanical systems, see [14,15] and references therein. In such systems there was observed an unexpectedly large change of the amplitude dependence of the vibration frequency with the varying electron density [16,17]. When the doping level was increased from 2.8 × 10 18 cm −3 to 5.9 × 10 19 cm −3 , the nonlinearity parameter increased by more than an order of magnitude. Moreover, the nonlinearity change was different for the vibrational modes with different spatial structure.
In this paper we develop a theory of the nonlinearity of vibrational modes in semiconductor nano-and micromechanical systems with high electron density. We show that the electron-phonon coupling can lead to a strong self-action of the vibrational modes, which in turn significantly modifies the amplitude dependence of the mode frequencies. We find the dependence of the effect on the electron density and temperature.
For bulk semiconductors, the effect of the electronphonon coupling on the elastic properties, including the three-phonon coupling, was first analyzed by Keyes [18]. The analysis referred to n-Ge and was based on the deformation potential approximation. The idea was that deformation lifts the degeneracy of the equivalent electron valleys, which leads to a redistribution of the electrons over the valleys. In turn, such redistribution changes the speed of sound depending on the direction and polarization of the sound waves and also affects the sound speed in the presence of uniaxial stress. This theory was extended to silicon and the corresponding measurements were done by Hall [19]. However, Hall also observed the change of the speed of transverse sound waves and the effect of stress on sound propagation in the geometries, where these effects are due to shear deformation and do not arise in the deformation potential model. A theory of the change of the linear shear elastic constant in silicon due to the intervalley redistribution of the electrons was developed by Cerdeira and Cardona [20].
As we show, in mesoscopic systems the strain-induced redistribution of the electrons over the valleys of the conduction band leads to the previously unexplored strong fourth-order nonlinearity of the vibrational modes. This nonlinearity gives a major contribution to the amplitude dependence of the vibration frequency. The redistribution also leads to a temperature dependence of the frequencies. The magnitudes of the effects sensitively depends on the mode structure. We describe them for several types of modes, including those studied in the experiment [16,17] and qualitatively compare the results with the observations. The theoretical results refer to both degenerate and nondegenerate electron systems. Specific calculations are done for silicon resonators.
In Sec. II we give, for completeness, the expressions for the mode normalization and the amplitude-dependent frequency shift of coupled nonlinear modes in a nano-arXiv:1611.09912v1 [cond-mat.mes-hall] 29 Nov 2016 or micro-system. In Sec. III and Appendix A we provide expressions for the electron-phonon coupling induced change of the elasticity parameters, including the parameters of quartic nonlinearity. In Sec. IV we discuss the asymptotic behavior of the parameters of quartic nonlinearity for low and high electron density and give their explicit form for silicon. In Sec. V we calculate the nonlinear frequency shift for several frequently used vibrational modes in single-crystal silicon systems and show the dependence of this shift on the electron density and temperature. The explicit analytical expressions are given in Appendices C and D. Sec. VI contains concluding remarks.

II. NONLINEAR FREQUENCY SHIFT OF LOW-FREQUENCY EIGENMODES
Of primary interest for nano-and micro-mechanical systems are comparatively low-frequency modes with wavelength on the order of the maximal size of the system. Examples are provided by long-wavelength flexural modes of nanotubes, nanobeams, and nano/micromembranes, or acoustic-type modes in microplates or beams. These modes are easy to excite and detect. We will enumerate them by index ν. Their dynamics is described by the elasticity theory [21]. The spatial structure of the displacement field of a mode u (ν) (r) in the harmonic approximation is determined by the boundary conditions. We will choose u (ν) (r) dimensionless, so that in our finite-size system Here, V is the volume of the system. We assumed that the mode eigenfrequencies ω ν are nondegenerate; including degenerate modes is straightforward. For simplicity, we also assumed that the system is spatially uniform; an extension to spatially nonuniform systems is straightforward as well.
We emphasize the distinction of the normalization (1) from the conventional normalization for bulk crystals, where ν corresponds to the wave vector and the branch number, and the normalization integral is independent of the volume. The normalization (1) is convenient for the analysis of low-frequency modes with the discrete spectrum characteristic of mesoscopic systems. Such modes are standing waves, and therefore vectors u (ν) can be chosen real.
The low-frequency part of the displacement can be written as Functions Q ν (t) give the mode amplitudes. In the harmonic approximation the dynamics of the standing waves is described by the Hamiltonian where P ν is the momentum of mode ν and M is the mass of the system. The anharmonicity of the crystal leads to mode-mode coupling. Within the elasticity theory this coupling is described by the terms in the Hamiltonian, which are cubic and quartic in the strain tensor. We will not consider higher-order terms, which are small for the mode amplitudes of interest. From the expansion (2), we obtain the nonlinear part of the Hamiltonian in the form Equation (4) is essentially an expansion in the ratio of the mode amplitudes to their characteristic wavelength, which is of the order of the appropriate linear dimension of the system. This is why mesoscopic systems are of particular interest, as here vibrations of low-frequency eigenmodes become nonlinear for already small vibration amplitudes.
A familiar consequence of nonlinearity in nano-and micromechanical systems is the dependence of the vibration frequency of a mode on its own amplitude and on the amplitudes of other modes, see Ref. 11 for a review. In particular, the change δω ν of the mode frequency due to the vibrations of the mode itself, Q ν (t) = A ν cos ω ν t, is [22,23] where γ ν ≡ γ νννν and we kept the terms of the first order in γ and the second order in β.
The nonlinear mode coupling (4) leads also to the frequency shift due to thermal vibrations of the modes. The dominating contribution to this shift for low-frequency modes comes from their coupling to modes with frequencies ∼ k B T / , which have a much higher density of states. This shift is described by an expression that is similar to Eq.
and placed under the sum over ν , in the classical limit.

III. THE NONLINEARITY DUE TO THE ELECTRON-PHONON COUPLING
We will consider the vibration nonlinearity due to the electron-phonon coupling in multi-valley semiconductors with cubic symmetry, silicon and germanium being the best known examples. In such semiconductors, the energy valleys of the conduction band are located at highsymmetry axes of the Brillouin zone. Strain lifts the symmetry and thus the degeneracy of the valleys.
The simplest mechanism of the electron-phonon coupling is the deformation potential. Here, the energy shift δE α of valley α is determined by the deformation potential parameters Ξ u and Ξ d of the coupling to a uniaxial strain along the symmetry axis of the valley and to dilatation, respectively. In terms of the strain tensor ε ij we have , with e (α) being the unit vector along the symmetry axis of the valley. We use the hat symbol to indicate tensors and symbol "⊗" to indicate tensor products. The analysis below is not limited to the deformation potential approximation. An important extension will be discussed using silicon as an example.
We assume that the strain varies in time and space slowly compared to the reciprocal rate of intervalley electron scattering and the intervalley scattering length, respectively. Then the electron system follows the strain adiabatically. The electron density n (α) (r) in valley α is decreased or increased depending on whether the bottom of the valley goes up or down. In the single-electron approximation and for the deformation potential coupling, the electron free energy density for a given strain is is the free energy density for electrons with density n(r) in a valley in the absence of coupling to phonons.
The electro-neutrality requires that the total electron density summed over the valleys be constant. The free energy density F e has to be minimized over n (α) (r) to meet this constraint. This gives the change of the electron chemical potential δµ due to strainε. The resulting increment of the electron free energy density has the form of a series expansion in the strain tensor, Here Λ 1 , Λ 2 , Λ 3 , and Λ 4 are tensors of ranks 2, 4, 6, and 8, respectively. They are contracted with the tensor products of the strain tensorε. Respectively, Λ k are the electronic contributions to the linear (for k = 2) and nonlinear (for k > 2) elasticity parameters of the crystal. These contributions are isothermal, but since the change of the mode frequencies from the electron-phonon coupling is small and the nonlinearity is also small, the difference with the adiabatic expressions can be disregarded.
To the third order inε the expression for δF e in terms of the shift of the valleys was found by Keyes [18] in the analysis of sound wave propagation. However, to find the parameters of the quartic nonlinearity of resonant modes in small systems, which is of primary interest to us, we also need to keep quartic terms in Eq. (6).
As seen from the explicit form of the parameters of the expansion (6) given in Appendix A, where µ 0 is the electron chemical potential in the absence of strain; it is determined by the total (summed over the valleys) elec-tron density n. Of central importance for the analysis is that parameter Ξ u / max(µ 0 , k B T ) ∼ 10 3 for electron densities n ∼ 10 19 cm −3 and room temperatures, i.e.
As a consequence, the coefficients at the nonlinear inε terms in Eq. (6) quickly increase with the increasing order of the nonlinearity [the overall series (6) is converging fast because of the smallness of the strain tensor]. The increase of Λ k with k allows us to keep inε only the terms linear in the lattice displacement, i.e., to set ε ij = (1/2)(∂u i /∂x j + ∂u j /∂x i ), where u i and x i are the components of the displacement and the coordinates, respectively. Indeed, in this case a kth term of the series (6) is of order k in the displacement. If we included the quadratic in ∂u i /∂x j term into one of theε tensors in the kth term, this term would become of order k + 1 in the displacement. However, for linearε the (k + 1)th term in the series (6) is also of the (k + 1)th order in the displacement, but is larger by factor Ξ u / max(µ 0 , k B T ).
For linearε, the total strain is a sum of partial contributions of strain from individual modes. For mode ν, such partial contribution is expressed in terms of the We note that, in contrast to the dimensionless strain tensorε, ten- From Eq. (6) we find the electronic contributions to the nonlinearity parameters β whereε (ν) ≡ε (ν) (r); tensors Λ k are independent of r. Similarly, the electronic contribution to the eigenfrequency is Generally, the term ∝ Λ 2 leads to mode mixing; however, if the mode frequencies are nondegenerate, this mixing is weak and can be disregarded, to the leading order in the electron-phonon coupling. One can see that the effect of the static stress ∝ Λ 1 can be disregarded as well.
The frequency change (9) depends on temperature because of the temperature dependence of Λ 2 . The nonlinearity (8) also leads to a temperature dependence of the mode eigenfrequency. Together they modify the temperature dependence of the mode eigenfrequencies compared to that of undoped crystals. This modification often weakens the temperature dependence of the eigenfrequencies, which proves very important for applications of micro-mechanical systems in devices that work in a broad temperature range [24]. Equations (6) -(9) are generic and apply beyond the deformation potential approximation. This is of particular importance for silicon. Here, the electron band valleys lie on the 100 -axes close to the X-points on the zone boundaries where two electron energy bands cross. Lattice strain can lead to a band splitting at X-points and a shift of the valleys [25,26]. Importantly, this shift results from a shear strain, which does not lead to a linear in the strain shift in the deformation potential approximation. The valley shift is quadratic inε in this case, as explained in Appendix A, which corresponds to an effectively two-phonon coupling. The coupling parameter Ξ sh is quadratic in the strain-induced band splitting, see Eq. (A2). It is large, much larger than the constant Ξ u . Therefore the arguments given below Eq. (7) apply in this case as well. For purely shear strain in silicon, terms of odd order inε in δF e , Eq. (6), vanish.

IV. EXPLICIT FORM OF THE TENSORS OF NONLINEAR ELASTICITY
Tensors Λ n can be obtained by minimizing the free energy density of the electron system for a given strain and expanding the result in a series inε. A general procedure that allows one to find the components Λ n for n ≤ 4 is described in Appendix A. Using the symmetry arguments, the elasticity tensors are conveniently written in the contracted (Voigt) notation where the symmetric strain tensor is associated with a six-component vector. Then the nonlinear elasticity tensors Λ 3 and Λ 4 become tensors of rank three and four in the corresponding vector space. We use notation δĉ for tensors Λ in these notations to emphasize that we are calculating corrections to the nonlinear elasticity tensors due to the electron-phonon coupling.
The explicit expressions for the nonlinear elasticity tensors δĉ are given in Table I. They refer to silicon and include the contributions that come from both the deformation potential coupling and from the splitting of the electron bands due to shear strain. In the deformation potential approximation, the components of the third-rank tensor δĉ, which determine the cubic in the strain terms in the free energy, were found earlier [19]. Therefore we give only the components that contain a contribution from shear strain.
The fourth-rank tensor δĉ determines the quartic in the strain terms in the free energy and has not been discussed before, to the best of our knowledge. We give all independent components of this tensor. It is expressed in terms of the derivative of the electron density n over the chemical potential in the absence of strain µ 0 , which is a familiar thermodynamic characteristic. It is intuitively clear that the considered effect of the change of the electron density in different valleys in response to strain should be related to the derivative dn/dµ 0 . Interestingly, because we consider nonlinear response to strain, the expressions in Table I contain also higher-order derivatives of n over µ 0 . As we will see, this leads to a nontrivial behavior of the nonlinear frequency shift with varying temperature and density. The considered mechanism of the strain-induced inter-valley electron redistribution does not contribute to the components c 1123 and c 1456 , therefore δc 1123 = δc 1456 = 0.

A. Nonlinear elasticity in the limiting cases
The expressions for δĉ simplify in the case of low doping (or high temperature), where the electron gas is strongly nondegenerate, and in the opposite case of a strongly degenerate electron gas. For a nondegenerate gas, where the chemical potential in the absence of strain is µ 0 < 0, |µ 0 | k B T , we have in Table I Table I depend only on temperature, The decrease of the nonlinear elasticity parameters with increasing temperature in a nondegenerate electron gas is easy to understand. The effect we consider is determined by the competition between the energetically favorable unequal population of the electron energy valleys in a strained crystal and the entropically more favorable equal valley population. With increasing temperature the entropic factor becomes stronger, leading to a smaller population difference and thus smaller effect of the electron system on the vibrations.
The results on the asymptotic behavior of the corrections to nonlinear elasticity are not limited to silicon. Since parameters C 1,2,3,4 are given by the coefficients in the general expansion of the free energy in strain, (A1), these results can be applied to the nonlinear elasticity induced by the electron-phonon coupling in other multivalley semiconductors. To illustrate this point, in Appendix B we give δĉ tensor in germanium.
The difference between the asymptotic behavior of the tensors δĉ in the limits of nondegenerate and strongly degenerate electron gas can lead to a peculiar density and temperature dependence of the nonlinear frequency shift of the vibrational modes. It comes from the coefficients C 1,...,4 containing higher-order derivatives of n with respect to µ 0 . In the transition region µ 0 ∼ k B T , thinking of the competition between the entropic and energetic factors does not provide a simple insight into the behavior of δĉ, as both the energy and the entropy are TABLE I. The change of the components of the nonlinear elasticity tensors due to the strain-induced electron redistribution between equivalent energy valleys in doped silicon. The coordinate axes are chosen along the 100 axes. Parameter Ξ sh characterizes the effectively two-phonon coupling to shear strain. This parameter as well as function F 1/2 (x) are defined in Appendix A; x = µ0/kBT and n is the electron density.
complicated functions of density and temperature.

V. DOPING-INDUCED NONLINEARITY OF SIMPLE VIBRATIONAL MODES
The nonlinear elasticity tensors in Table I give the doping-induced contributions to the nonlinearity parameters of the eigenmodes of micro-and nanomechanical systems. These contributions are described by Eq. (8). As mentioned before, an important characteristic of the mode nonlinearity is the dependence of the mode frequency on the vibration amplitude. To the leading order, it is given by Eq. (5). This dependence has a contribution from the nonlinearity of an undoped crystal, which is quadratic in the parameters of the cubic nonlinearity; for example, if the latter is described by the Grüneisen constant, the corresponding contribution is quadratic in this constant. It is typically small. There is also a contribution from the quartic nonlinearity; the parameters of such nonlinearity are not known in undoped crystals and are not expected to be large. Respectively, the amplitude dependence of the vibration frequency for low-frequency modes in weakly doped single-crystal micro-mechanical systems is relatively weak [17].
A feature of the doping-induced nonlinearity described by Table I is that the quartic in the strain term in the free energy has a large coefficient compared to the cubic term, cf. Eq. (7) and the discussion below this equation. Therefore, in Eq. (5) for the amplitude dependence of vibration frequency one can keep only the Duffing nonlinearity constant γ ν . The contribution from the cubic nonlinearity terms ∝ β 2 ννν can be disregarded. For a mode ν, the doping-induced contribution to γ ν is equal to γ (e) νννν in Eq. (8).
To find the dependence of the mode frequency on the vibration amplitude we go through the following steps. First, we find the normal modes of interest for the given geometry of the system, with account taken of the boundary conditions, and normalize the displacements u (ν) (r) as indicated in Eq. (1). We use u (ν) (r) to find the strain tensorε (ν) (r). The result is substituted into Eq. (8) and is convoluted with tensor Λ 4 , giving the value of γ ν , which is then used in Eq. (5) to find the frequency dependence on the vibration amplitude δω ν .
Of particular interest is the relative frequency shift δω ν /ω ν . To find this shift to the leading order, one can disregard nonlinearity when calculating the eigenfrequency ω ν . Then, from Eq. (5), where Λ (f) 2 is the full tensor of linear elasticity, which includes the major term of the linear elasticity of the undoped crystal and the doping-induced correction Λ 2 .
An important feature of the relative shift δω ν /ω ν is its scaling with the size of the system. The vibration amplitude A ν in Eq. (10) can be scaled by the lateral dimension L, for example the length of a nanobeam or a nanowire for an extension mode, or the size of the square for a Lamé mode, or the diameter of a disk for a breathing mode in a disk. Respectively, we write A ν = η ν L. Then, if one takes into account the explicit form (8) of the parameter γ ν = γ (e) νννν , one finds from Eq. (10) that the ratio δω ν /(η 2 ν ω ν ) is independent of the system size for the aforementioned modes. In this estimate we used that the tensors Λ are material parameters and are independent of the geometry. We also used that the modes of interest have typical wavelength ∼ L, and thereforeε (ν) scales as L −1 .
Most of the experiments in nano-and micromechanics are done with nanobeams, nanowires, membranes, or thin plates. In such systems the thickness is much smaller than the length or, in the case of membranes or plates, the lateral dimensions. Then, from the boundary condition of the absence of tangential stress on free surfaces [21], it follows that the strain tensorε weakly depends on the coordinate normal to the surface. This simplifies the denominator in Eq. (10), making it proportional to the thickness. Similarly, from Eq. (8) γ ν is also proportional to the thickness, and the thickness drops out of Eq. (10).
The explicit expressions for M ω 2 ν and γ ν that determine the denominator and the numerator in Eq. (10), respectively, are given in Appendices C and D for Lamé and extension modes. These expressions are cumbersome, and it is convenient to use symbolic programming to obtain them. [27] A. Temperature and electron density dependence of the scaled nonlinear frequency shift The scaled ratio δω ν /(η 2 ν ω ν ) that characterizes the relative nonlinear frequency shift is shown in Fig. 1 for several modes that are often used in single-crystal silicon MEMS. This ratio depends on the type of the mode and the crystal orientation. Figure 1 refers to high-symmetry crystal orientations, in which case the modes have a comparatively simple spatial structure and the surfaces can be made smooth. We used the values Ξ u = 8.8 eV [28], Ξ sh = 300 eV, the effective mass for density of states m eff = 0.32m e [29], and the temperature-dependent linear elasticity parameters given in Ref. [30]. Figure 1 shows that the electron-redistribution induced nonlinearity of vibrational modes is very strong. For the ratio of the vibration amplitude to the system size η ∼ 10 −4 and the mode eigenfrequency ω ν /2π ∼ 10 MHz, the frequency change can be as a large as δω ν /2π ∼ 0.1 kHz. This explains, qualitatively, the observations [17]. A quantitative comparison with the experiment [17] is complicated, as the observations refer to different samples. Our preliminary results show an excellent quantitative agreement with the data obtained for the same sample at different temperatures and for different types of modes [31].
The nonlinear frequency shift displays several characteristic features, as seen from Fig. 1. One of them is the strong dependence of the shift on the type of the mode and the crystal orientation. For both the Lamé and the extension mode, the shift is much stronger for crystals cut out in 100 direction than in 110 direction. This is a consequence of the electron energy valleys lying along the 100 axes, making the system more "responsive" to the lattice displacement along these axes. Interestingly, in the both configurations the shifts for the Lamé modes are larger than for the extension modes.
A somewhat unexpected feature is the nonmonotonic dependence of the nonlinear frequency shift on the electron density and temperature. The nonmonotoncity occurs in the range where the electron system is close to degeneracy, µ 0 /k B T ∼ 1, and it strongly depends on the crystal orientation. It is much stronger for crystals cut in 100 than 110 directions. For a crystal cut in 110 direction, both the density and temperature dependence of the shift are monotonic in the case of the Lamé mode, whereas for the extension mode the nonmonotonicity is weak.
The nonmonotonicity of the frequency shift stems from the behavior of the parameters nC 2,3,4 in the range µ 0 ∼ k B T . As seen from Table I, parameter nC 2 exponentially increases with the increasing µ 0 /k B T for negative µ 0 /k B T , but for large positive µ 0 /k B T it falls off as (µ 0 /k B T ) −3/2 . It has a pronounced maximum for µ 0 /k B T ≈ 0.6. Parameter nC 3 also displays a maximum, which occurs for µ 0 /k B T ≈ 1.1. In contrast, parameters nC 1,4 depend on µ 0 /k B T monotonically.
The results of Appendices C and D show that, for the Lamé and extension modes in crystals cut in 100 direction, the relative shift δω ν /ω ν is determined by coefficient nC 2 , which explains the nonmonotonicity of the shift. For crystals cut in 110 , the shift of the Lamé mode is fully determined by coefficient nC 4 and is monotonic, whereas for the extension mode the expression for the shift has contributions from nC 2 , nC 3 , and nC 4 that partly compensate each other, leading to a comparatively small shift all together and its weak nonmonotonicity.

VI. CONCLUSIONS
The results of this paper show that the electronphonon coupling strongly affects the nonlinearity of vibrational modes in semiconductor-based nano-and micromechanical systems. The mechanism of the effect is the strain-induced redistribution of the electrons between the valleys of the conduction band. The redistribution results from lifting the degeneracy of the electron energy spectrum by the strain from a vibrational mode. The analysis refers to the range of temperatures where the rate of intervalley scattering strongly exceeds the frequencies of the considered modes. In this case the valley populations follow the strain adiabatically.
The change of the valley populations is a strongly nonlinear function of the strain tensor. The respective expansion of the free energy in the strain is an expansion in the strain multiplied by the ratio of the electron-phonon coupling energy (in particular, the deformation potential) to the chemical potential of the electron system or the temperature. This ratio is large, 10 3 . It is this parameter that makes the nonlinearity of the vibrational modes in doped semiconductor structures strong.
Of special interest in nano-and micromechanical systems is the amplitude dependence of the vibration frequency. To the leading order, it is determined by the quartic terms in the expansion of the free energy in strain. These terms are comparatively large in doped crystals.
We have calculated the nonlinear elasticity tensor that describes the electron contribution to the terms in the free energy, which are quartic in the strain. The explicit expressions for the tensor components refer to semiconductors with the valleys on 100 axes, in particular, to silicon. We have also found this tensor for germanium. In silicon, along with the deformation potential coupling, an important role is played by the coupling to shear strain. Such strain lifts the band degeneracy at the zone boundary and is effectively described by a two-phonon coupling. We show that this coupling also leads to strong nonlinearity of vibrational modes.
The parameter of the electron coupling to shear strain in silicon is not easy to access in the experiment [26,32]. Measurements of the nonlinear frequency shift provide a direct means for determining this parameter. In particular, the nonlinear frequency shift of the fundamen-tal Lamé mode in a silicon plate cut along 110 axes is determined by this parameter only, except for small corrections from the nonlinearity of the undoped crystal.
We found that the nonlinear frequency shift strongly depends on the type of a vibrational mode and the crystal orientation. We also found that the ratio of the frequency shift to the squared vibration amplitude can be profoundly nonmonotonic as a function of electron density and temperature. The results provide an insight into the experimentally observed strong mode nonlinearity in doped crystals [17]. In terms of applications, they enable choosing the appropriate range of doping and the temperature regime to optimize the operation of nano-and micromechanical resonators.
The major effect of a strain on the electron free energy comes from the shift of the energy valleys. We will assume that valley α is shifted in energy by δE α and the shift is small, |δE α | max(k B T, µ 0 ), where µ 0 is the chemical potential in the absence of strain. We further assume that the vibrations are slow compared to the time it takes the electron system to come, locally, to thermal equilibrium for given values of δE α , i.e., the temperature and the chemical potential are the same in all valleys. Since for high electron densities the thermal conductivity is high, the change of the temperature compared to the ambient temperature can be disregarded; also, as mentioned in the main text, the electron density n summed over all valleys is constant.
Expanding the electron free energy density to the 4th order in the strain-induced shifts δE α , we find that, in an N -valley semiconductor, the change δF e of the free energy density is Here, ∆ m ε = N −1 α (δE α /k B T ) m . We use the standard notation F 1/2 (x) = ∞ 0 dy y 1/2 /[1 + exp(y − x)]; primes indicate differentiation over x, for example, F 1/2 ≡ dF 1/2 /dx. Function F 1/2 and its derivatives are calculated for x = µ 0 /k B T . Equation (A1) immediately gives the tensors Λ n of the expansion of the free energy increment (6) if one expresses the shift δE α of the valleys in terms of the strain tensor. In the deformation potential approximation the relation between δE α andε is given in the main text, see also Eq. (A2) below.
In the case of Si crystals, which are often used in micromechanical resonators, an important contribution to δE α comes from the shear-strain induced splitting of the electron energy bands at the zone boundary. Shear strain does not lead to the valley shift in the deformation potential approximation. The overall shift of valley α, to the lowest order in the coupling that causes it (i.e., to the first order in the deformation potential where its contribution is nonzero and to the second order in the band splitting for shear strain) is [26]: Here we use that silicon has six valleys located at the 100 axes, and we chose the coordinate axes x, y, z along 100 . Respectively, the valley index α takes on three values that correspond to the x, y, z axes (the valleys lying on the same axis, but in the opposite directions, are equivalent). The strain ε α , which enters the second term in the right-hand side of Eq. (A2), is a component of the strain tensor ε ij with i, j such that i, j = α and i = j. The parameter 2Ξ u is the interband matrix element of the electron-phonon coupling calculated for the electron conduction bands ∆ 1 and ∆ 2 at the X point on the boundary of the Brillouin zone, where the bands cross; ∆E is the energy separation between the bands ∆ 1 and ∆ 2 at the value of the wave vector k that corresponds to the conduction band minimum. Parameter Ξ sh is the effective deformation potential of two-phonon coupling to shear strain. The numerical value of Ξ sh is not well known. The experimental data give Ξ u ≈ 7 − 8 eV [26,32] and the numerical data on the band splitting give ∆E ≈ 0.7 eV [33] so that Ξ sh is in the range of 280 − 360 eV; this is essentially an order of magnitude estimate.
In calculating δF e in Eq. (A1) we kept terms that are quartic inε. The components of the tensors Λ k in Eq. (6) are expressed in terms of δF e as Tensors Λ are symmetric with respect to the interchange of indices i k ↔ j k and the pairs (i k j k ) ↔ (i k j k ). For the considered long-wavelength strain, tensors Λ k are independent of coordinates. The corrections Λ 2 to the linear elasticity tensors were found previously [19,20] and are not discussed in this paper. In this section we provide the corrections to the nonlinear elastic constants of germanium, which are due to the redistribution of the electrons over the valleys. Germanium has four equivalent valleys in the conduction band, which are located on the boundary of the Brillouin zone along 111 axes. We use the Voigt notation and write the components of the corrections to the nonlinear elasticity tensor δĉ in the frame where the axes (x, y, z) are along the 100 directions of the crystal. Using the results of Appendix A, we obtain The notations are the same as in Appendix A and in Table I. The electron-phonon coupling does not contribute to the other third-and fourth-order elastic constants.
Corrections δc 44 and δc 456 for germanium were found by Keyes [18]; however, his final expression for δc 456 differs from Eq. (B1) by a factor of 4 (our expressions for δc 44 coincide with Ref. [18]). Parameters δc 4444 and δc 4455 have not been found before, to the best of our knowledge. In the limiting cases, corrections δc 4444 and δc 4455 have the same dependence on temperature and electron density as constant nC 2 discussed in Sec. IV A.
If we consider silicon and take into account only the contribution δĉ to the nonlinear elasticity tensorĉ, with the account taken of Table I, the expression for γ ν simplifies to γ ν = (27π 4 h/32L 2 )δc 1111 .
Expressions (D2) and (D4) were generated using a computer code to calculate the sums and integrals in Eq. (8).