Compressional mode softening and Euler buckling patterns in mesoscopic beams

We describe a sequence of Euler buckling instabilities associated with the transverse modes of a mesoscopic beam subjected to compressional strain. As the strain is increased, successively higher normal mode frequencies are driven to zero; each zero signals an instability in the corresponding normal mode that can be realized if all lower instabilities are suppressed by constraints. When expressed in terms of the critical buckling modes, the potential energy functional takes the form of a multimode Ginzburg–Landau system that describes static equilibria in the presence of symmetry breaking forces. This model is used to analyse the complex equilibrium shapes that have been observed experimentally in strained mesoscopic beams. Theoretically predicted critical strain values agree with the appearances of higher order mode structures as the length-to-width aspect ratio increases. The theory also predicts upper bounds on the individual mode amplitudes that are consistent with the data. Based on insights from the theory, we suggest possible origins of the buckling patterns.


Introduction
The ability to fabricate freely suspended structures of sub-micrometre dimensions has enabled the thermal and mechanical properties of precisely engineered doubly clamped beams to be studied [1]- [3]. Using materials of relatively low mass density and high Young's modulus, the geometry of the beams can be designed to produce high quality mechanical resonators with fundamental flexural motion in the gigahertz regime [4]. Doubly clamped structures have been discussed as the basis for highly sensitive electrometers [5,6], mass detectors [7,8], and memory schemes [9]- [12]. An intriguing aspect of sub-micrometre beams is that when cooled to the milli-Kelvin regime the thermal energy can become comparable to the phonon or mechanical energy scale set by one or more dimensions of the structure, thereby enabling access to quantum characteristics of the thermal (phonon) properties [13]- [16] and the study of macroscopically distinct quantum states [17,18]. Thermomechanical [19] and quantum squeezing have also been explored [3,20]. For a general discussion of quantum related aspects of mechanical structures, the reader is referred to a recent review [21].
Of importance to this paper are nonlinear mechanical effects associated with doubly clamped mesoscopic beams. Two nonlinear systems that have received attention are the Duffing nanomechanical resonator [22,23] and the Euler beam [24]. In the latter, compressive strain has been used to produce a tunable two-state system [25,26] with which it has been shown that the buckling properties have features that classical Euler buckling theory does not usually describe. In particular, the static buckling characteristics were found to depend on the length-to-width ratio of the beam. At low aspect ratios, the buckled shape is that expected for a macroscopic beam, whereas at large ratios the buckling assumes highly complex, stable shapes [27]. It was also demonstrated that external forces can induce transitions between the buckled shapes [28].
In this study, we provide a unified analysis covering both simple and complex buckling behaviour. The analysis is based upon the existence of Euler instabilities associated with higher transverse vibrational modes of the beam. The frequencies of all such modes are softened by longitudinal compression, giving rise to a succession of higher-order Euler instabilities as the frequency of the corresponding mode is driven to zero. In an ideal beam, a given higher instability is realized only if all lower ones are suppressed by constraints. The higher instabilities were observed in a recent experiment on the buckling of thin (but macroscopic) plates subject to 3 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 1. A scanning electron micrograph of a doubly clamped SiO 2 buckled beam from [28]. The scale bar represents 10 µm. The length-to-width ratio of the beam is 132.
transverse confinement [29]. Differential boundary conditions can also produce shapes in the form of a higher mode [30], but these are externally induced, not spontaneous. In the study reported in [27,28], the complex shapes are seen in the absence of known constraints, as shown in figure 1. While the microscopic origins of these shapes are unknown, in this study we provide an explanation of the evolution to increasingly more complex shapes as the aspect ratio is increased. We also suggest possible origins of the observed patterns.
We begin by describing the physical model of the ideal beam subject to fixed endpoint separation. In section 3, we describe the normal modes as a function of longitudinal strain, and show that their frequencies are driven to zero at a succession of critical strain values. We briefly discuss the interacting phonon interpretation [31] and the meaning of phonons as the first Euler instability is approached [24]. We then describe the basis set of critical buckling modes in section 4. Each of these is defined at its own critical strain value, where its frequency is driven to zero. We show that these static modes comprise an appropriate basis set for describing the static equilibrium observed experimentally [27]. We use these modes in section 5 to develop a multimode Ginzburg-Landau theory, adding symmetry breaking effects in the form of transverse forces. The Ginzburg-Landau equations are applied in section 6 to an analysis of recently observed complex beam shapes that exhibit higher modes, such as that shown in figure 1, leading to insights about their formation. The main conclusions are summarized in section 7.

The model
The physical system we begin with [24] is an ideal elastic beam supported at both ends, having length L, width d, and intermediate dimension larger than d so that only transverse displacements in the d direction need be considered. The Lagrangian, as a functional of this transverse displacement, y(x, t), is where µ is the mass per unit length along the beam and V is potential energy functional. Assuming that the endpoint separation L is fixed, this functional is, to fourth order, where F is the linear elastic modulus of the beam. The first term is the bending energy (with b = d/ √ 12 the bending moment), and the second and third terms are the stretching (or compressional) energies, where = (L − L o )/L o is the longitudinal strain, with L o the equilibrium length of the bar. The third term accounts for the added strain due to the actual shape of the bar, as may be more apparent from the associated Euler-Lagrange equation of motion, The potential energy functional takes a different form depending upon whether one considers fixed endpoint separation or fixed applied tension (positive or negative). In the latter case, the quartic term in equation (2) is absent, and one must go beyond the quadratic contribution to the bending energy in order to stabilize the beam beyond its Euler instability. The full bending energy integral [32] is ds(y ) 2 (1 − (y ) 2 ) −1 , where s is the arclength along the bent beam, and the resulting quartic contribution is This contribution will not be included in our model because it is negligible compared to the quartic term in equation (2). We shall estimate their relative sizes in section 6 and show that equation (4) is orders of magnitude smaller for the beams discussed there. At this point it seems reasonable to assume that the physics of the fabricated beams conforms more closely to fixed endpoint separation than to fixed compressional stress.

Dynamic modes and compressional softening
The dynamic modes are the usual normal modes of small-amplitude oscillations, y(x, t) = f n (x) exp(−iω n t), which are eigenfunctions of the linearized equation (3), The precise form of the mode functions f n (x) depends on the boundary conditions: For the hinged case, with y(±l/2) = y (±l/2) = 0, the solutions are purely trigonometric, identical to those of a vibrating string under tension. For the clamped case, with y(±l/2) = y (±l/2) = 0, they involve hyperbolic functions as well, f n (x) = a n cos k n x + b n cosh q n x, (n = 1, 3, . . .), f n (x) = a n sin k n x + b n sinh q n x, (n = 2, 4, . . .).
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Table 1. Wavenumbers for the hinged and clamped boundary conditions, given in the latter case for two values of corresponding to the unstrained ( = 0) and critically strained ( = n ) cases. Other parameters, C n and γ n , defined in section 4, are given for the case of clamped boundary conditions at critical strain. Integer entries indicate exact values.
(k n L/π) (k n L/π) (k n L/π) C n γ n k n L/2 The wavenumbers in the hinged case are given by k n L = nπ (table 1), while in the clamped case there are two such parameters (k n and q n ), related by whose precise values are given by transcendental equations k n tan k n L/2 = q n tanh q n L/2, (n= 1, 3, . . .), k n cot k n L/2 = −q n coth q n L/2, (n= 2, 4, . . .).
The frequency depends explicitly on k n alone: where B is the bulk modulus and ρ is the mass per unit volume. Despite the appearance of transcendental equations in the clamped case, the wavenumbers nevertheless take simple approximate values in special cases. In the unstrained case ( = 0), equation (8) dictates that k n = q n , and equations (9) and (10) are then solved to close approximation by k n L = (n + 1 2 )π, with small corrections ≈2(−1) n−1 exp(−(n + 1 2 )π) that are negligible in all but the fundamental mode, as shown in column 2 of table 1. The other limiting cases of interest are the various critical strain values n at which the nth mode frequency is driven to zero by compression, or negative strain. These critical strain values are determined by equation (11): The frequencies (equation (11)) are conveniently re-expressed as where ω (o) n = bk 2 n √ B/ρ is formally the zero strain expression for the frequency, although with clamped endpoints k n itself has a weak strain dependence. As approaches the nth critical strain, n , both the wavenumber q n and frequency ω n are driven to zero (equations (8) and (11)), while the wavenumber k n remains finite. Equations (9) and (10) are now solved exactly by k n L = (n + 1)π for odd n, and approximately by k n L ≈ (n + 1)π − 4/(n + 1)π for even n, as reflected in column 3 of table 1.
The clamped mode eigenfunctions f n (x) reduce to simple forms at their respective critical strain values → n , because q n → 0. We refer to all such critical modes as static modes, since the corresponding eigenfrequencies vanish. As a comparison, one of the static mode functions, g 3 (x), is shown together with its dynamic counterpart, f 3 (x, = 0), in figure 2. The static modes will be discussed in detail in the next section.
The dynamic modes lead to an interacting transverse-phonon Hamiltonian, as introduced in [31]. Defining the usual Hamiltonian is the canonical momentum density, one can transform to normal mode amplitudes, Y n and P n , as defined by The mode functions f n (x) are normalized according to dx f n (x)f m (x) = δ nm to preserve canonical commutation relations, [Y n , P m ] = ihδ n,m . Substituting equations (14), (15) into (1) and (2), one obtains the Hamiltonian The single sum is the usual free phonon Hamiltonian. The motion in each mode corresponds to that of a simple harmonic oscillator with mass m * n and oscillation amplitude Y n . The effective mass parameter is arbitrary, and we define it so that Y n corresponds to the maximum displacement of the beam in the mode n. In the case of hinged endpoints, with purely trigonometric mode 7 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT functions, one thus finds that m * n /m = 1/2. The situation is less straightforward for clamped endpoints, where the effective mass varies with the strain, and also depends weakly on the mode index. For example, when = 0, the m * n /m values are close to 3/8, while for → n , they are close to 1/8.
The double sum represents phonon-phonon interactions. These interactions are not of the most general type; they arise from the anharmonicity inherent in the geometry of transverse oscillations of the beam. They do not account for scattering events that would be induced by microscopic anharmonicity of the material or from defects. Their physical effect is to enhance phonon frequencies due to the presence of other phonons, which increases the arclength of the beam. A dramatic example of this effect [24] is that of a beam compressed just to its first Euler instability. At this point, the frequency of the fundamental mode that appears in the freephonon part of the Hamiltonian is driven to zero. Thus, as phonons are added in that mode, their frequencies are determined entirely by the quartic potential energy, which arises from interactions among these phonons.
In the context of the present study, the anharmonic term represents mode-mode interactions whose effect is to induce competition between instabilities in different modes. Closer examination will show that, as mentioned earlier, allowing an instability in the fundamental mode prevents instabilities in higher modes. Conversely, inducing a displacement in a higher mode can suppress the instability in the fundamental.

Critical buckling (static) modes
In the previous section we showed that the fundamental and higher transverse vibrational modes of the bar are driven to zero sequentially, at a succession of Euler instabilities, as the compressional strain increases. For the case of clamped boundary conditions, where the mode functions f n (x) are strain dependent, it is necessary to define explicitly the critical buckling modes (for short, static modes), and to show that these comprise an appropriate basis set, even though each is defined at a different strain value, so that mutual orthogonality is not guaranteed.
We define each static mode function, g n (x), as the corresponding dynamic mode function evaluated at its critical strain value, f n (x, → n ). For the case of clamped boundary conditions, we also alter the normalization constant. The functions for this case are thus defined as with constants C n and γ n listed in table 1. While these functions are not mutually orthogonal, their first and second derivatives are (see equation (19)). To show this, note first that the linearized equation (5) is, at zero frequency, a Sturm-Liouville equation for y (x) itself, where now plays the role of the eigenvalue. This ensures that the g n are mutually orthogonal. It is straightforward to show, by direct evaluation, that this property also holds for the first derivatives.
The constants C n appearing in equation (18) are chosen so that all derivatives satisfy the orthonormality conditions, (g n , g m ) = k 2 n δ nm , (g n , g m ) = k 4 n δ nm , 8 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT the inner product being defined as the integral over the length of the bar. These conditions also hold, trivially, for the case of hinged boundary conditions. Fortunately, the potential energy functional depends only on derivatives with respect to x, so that its quadratic part is diagonal in the mode index, and the g n (x) indeed form an appropriate basis for time-independent problems.

Multi-mode Ginzburg-Landau potential
Suppose that the bar is subjected to a time-independent transverse force, distributed according to the profile F(x), and independent of the transverse displacement y(x). The equilibrium shape of the bar is determined by the equation (see equation (3)) The equilibrium shape,ỹ(x), that the bar would take if released at its endpoints, is the solution of this equation in the absence of the total strain (applied plus induced), as represented by the two second derivative terms. That is, This relation allows us to parameterize the applied transverse force in terms of the equilibrium shape functionỹ(x). That is, givenỹ(x), one can substitute equation (21) into equation (20) and solve the latter for the shape of the bar, y(x), in the presence tension or compression. This solution y(x) minimizes the potential energy functional Thus, the bending potential energy is a quadratic functional of the departure of the exact shape y(x) from the unstrained equilibrium shapeỹ(x). Replacing the transverse applied force profile F(x) by this shape functionỹ(x) allows us to treat externally applied forces and (possibly unknown) internal stresses on the same footing. At this point, it is useful to expand y(x) in terms of the static mode functions introduced in section 4, where the arbitrary prefactor is chosen to make Y n the physical peak to peak displacement, or more precisely, the caliper, 3 of the nth mode contribution. There is no distinction between the caliper and the usual (one-sided) displacement for the clamped static mode functions. This is clear from figure 2, and it also reflects the fact that the effective mass ratio of the dynamic modes 9 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT (equation (14)) approaches 1/8 in the limit as → n (exactly for odd n and approximately for even n). But clearly there is a factor of two difference for the hinged case, in which the mode functions are purely trigonometric for all values of the strain (recall that the effective mass ratio is 1/2 for this case). Adopting the caliper as the definition of Y n has direct physical appeal for characterizing static configurations of the beam. In addition, it has the formal advantage of making the potential energy functional look exactly the same for the two sets of boundary conditions. All numerical coefficients are the same; boundary conditions enter only through the values of the wavenumbers and critical strain parameters. Thus, expanding equation (22) in the mode functions of equations (17) and (18) for the clamped case, or their purely trigonometric counterparts for the hinged case, one obtains where a term ∼O(Ỹ 2 ) was dropped because it is a constant. This potential energy functional has the form of a general multi-mode Ginzburg-Landau potential, whereα nỸ n plays the role of an applied field that couples to the nth component of the order parameter, Y n . In the analogous magnetic system, these correspond to Fourier components of an applied magnetic field profile, H n , and the magnetization of the system, M n , respectively. 4 One could imagine the field to be applied externally, or internally (for example, arising from fixed microscopic magnetic moments that are not part of the system being considered). Because of the analogy with applied fields in general, we shall refer to the shape function componentsỸ n as biasing parameters. These could arise either from externally applied forces, or from internal stresses in the beam. In our model, theα n parameters depend only on the bending energies, while the α n depend in addition on the elastic strain energies: Thus, each α n can change its sign as | | increases beyond | n |, signalling a potential instability in the nth buckling mode. The applied strain, , is the single available tuning parameter in the system. A special aspect of the Euler beam problem is that theα and β parameters are all related by geometry. It follows from equations (24) and (26) that We are now in a position to compare the model of fixed endpoint separation, as represented by equation (27), with that of fixed applied stress. Equation (27) follows from the quartic term in equation (2). This contribution is absent in the fixed stress model, in which the quartic potential energy is provided by the bending energy, equation (4). It is easily verified that equation (4) Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT would lead to coefficients δβ that are smaller than those of equation (27) by (bk n ) 2 , which involves two inverse powers of the aspect ratio. As will be seen in the next section, this implies an approximate maximum value of δβ/β ∼ (d/L) 2 ∼ 10 −3 for those beams that show instabilities. This demonstrates that the quartic contribution in equation (2) is crucial in treating the case of fixed endpoints, while the quartic contribution to the bending energy may be ignored.
Relations (26) and (27) simplify V(Y n ) dramatically upon rewriting in dimensionless form: If we define the dimensionless displacement Y n ≡ √ 3Y n /2d, then (recalling the critical strain formula n = −b 2 k 2 n of equation (12)), the dimensionless potential energy V ≡ V/2FL is simply It is apparent that the relative strengths of all terms involve only the ratios of critical strains, which are just the squared ratios of wavenumbers. The relevant numbers from which to compute these ratios, for both sets of boundary conditions, are found in columns 1 and 3 of table 1. The extrema of this potential energy are given by the coupled Ginzburg-Landau-like equations One may conclude immediately from these equations that in ideal beams, where allỸ n = 0, solutions with nonzero Y n are confined to single modes. This is because the factor in brackets can vanish for only one mode at a time. These solutions are the nth nonzero solution appearing when | | > | n |. Furthermore, direct substitution into equation (28) shows that only the Y 1 solution is stable; all others correspond to saddle points of the potential and are unstable. The appearance of modes with n > 1, alone or in combination with others, implies that the corresponding biasing parametersỸ n are nonzero.

Analysis of observed beam shapes
In the experiments described in [27,28], one starts with a Si wafer coated with a 500 nm thermal oxide. SiO 2 beams are then delineated using electron-beam lithography and released from the Si substrate by plasma etching. There is a lattice mismatch between the Si and SiO 2 which leaves a compressive strain in the suspended SiO 2 beam. After fabricating and analysing more than 100 beams of varying aspect ratio in this manner, it appears that when a beam buckles under these conditions, the endpoints (which are fixed to the substrate) remain fixed to a very good approximation and the stress is at least partially relieved. The model of fixed endpoint separation is applicable on physical grounds for these beams, while a model of fixed longitudinal stress seems inappropriate. As shown in the paragraph following equation (27), the latter model would seriously underestimate the quartic terms in the Ginzburg-Landau potential.
References [25] and [28] show that the compressive strain due to the lattice mismatch is independent of the aspect ratio. The critical strains, on the other hand, decrease in magnitude as Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT the aspect ratio increases: the exact numerical factors given by column 3 of table 1. As the aspect ratio increases, first | 1 | drops below | | and the fundamental mode amplitude increases rapidly. As the aspect ratio increases further, successively higher critical strains drop below | |. For ideal beams, one would only expect the fundamental amplitude to continue to increase, with no further transitions involving higher modes. However, the data reported in [27] and reproduced here in figure 3 show something quite different. Transitions involving higher modes are actually observed, but within a probabilistic distribution allowing for coexistence with lower modes. This is compelling evidence for biasing. The solid lines are the theoretical upper bounds on individual mode coefficients. With the appropriate choice of biasing coefficients, the theory could reproduce the actual mode coefficients as represented by the points in figure 3, which are obtained from best fits to the individual beam shapes. The upper bounds could only be realized if all other biasing coefficients vanished, which is highly unlikely. Biasing was first identified through more subtle effects involving only the fundamental mode, where it was recognized as intrinsic asymmetry in the beams [25]. A brief review of this study on shorter aspect ratios will provide useful background.

Beams with shorter aspect ratios
In the earlier experimental results [25], it was found that the fundamental displacement Y 1 underwent a rapid, but not perfectly sharp transition as a function of aspect ratio. The slight rounding of the transition was attributed to inherent asymmetry, as modelled by a constant biasing parameterỸ 1 . This model successfully accounted for the rounding, including the small displacements observed at subcritical aspect ratios. In later experiments [28] with longer etching times, a shift in the transition to lower aspect ratios was observed. This shift could be well accounted for by an additional biasing contribution proportional to the observed Y 1 . This dependence suggests that asymmetry is somehow built in to the beam during or shortly after the deformation process. It could be that the etching itself is asymmetric, or it could be that the deformation process that occurs upon release of the beam during etching is partially plastic, although the mechanical response of the beam after the completion of this process has been shown to be elastic [28]. An analysis of the Y 1 data at low aspect ratios (L/d < 80) (see the first panel of figure 3) was performed in [28]. While the transition point depends on both and the biasing parameter Y 1 /Y 1 , the asymptotic slope of Y 1 as a function of aspect ratio well beyond the transition depends only on . This asymptotic slope is found from equations (30) and (31) and, in the units of figure 3, it is It was found that a good fit was possible, both for the position of the transition and the slope above the transition, withỸ 1 /Y 1 = 0.7 and | | = 3.8 × 10 −3 . Although the biasing parameter is different in [25], the values of are consistent, providing evidence for the sample independence of in the cited experiments. We refer to this regime of shorter aspect ratios as deterministic, because all of the samples appear to fall on a single curve of Y 1 as a function of aspect ratio. Although the biasing parameter 12 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 3. Experimentally observed coefficients Y n for modes n = 1-6 as a function of the aspect ratio. For a given mode, each point represents a different beam, as reported in [27]. The solid curves are the predicted equilibrium displacements (equation (30)) assuming all other mode coefficients are zero. Y 1 may vary with some randomness from sample to sample, the curve is only sensitive to this parameter close to the transition. For this reason the intrinsic strain | | is accurately determined in the low aspect ratio regime, somewhat beyond the transition.

Longer aspect ratios
As noted above, the data of [28] (figure 3) become probabilistic at larger aspect ratios, say L/d > 80, in that the Y 1 values no longer fall on a single curve. This onset of randomness is associated with a rapid rise of Y 2 and corresponding turnover in Y 1 , in conformity with a sum rule introduced in [27], which states that the sum over all modes appearing in equation (29) is restricted by It was shown that the data obey this inequality, and theoretical arguments were given based on the reasonable assumption that the biasing parameters are not too large. Its meaning is, in effect, that the actual arclength of the beam may not exceed its equilibrium length. It was argued further that the sum approaches | | asymptotically for large aspect ratios, which means that the arclength of the beam approaches its equilibrium length in this limit. Large L/d corresponds to | 1 | | |, as is satisfied well beyond the first Euler instability. The data are suggestive of this sum rule, but not conclusive, because we were unable to resolve amplitudes beyond Y 6 , and such amplitudes would be expected to contribute to the sum for L/d > 120.
The sum rule is consistent with the crossover from deterministic to probabilistic behaviour with increasing aspect ratio. That is, assuming that some randomness always exists in the biasing parameters, the sum rule leads to deterministic outcomes when only a single mode is admitted, while random outcomes are expected when two or more are significant.
Beneath this globally probabilistic behaviour, the individual mode amplitudes exhibit some striking features: (i) every higher mode that could be analysed (n = 2-6) undergoes a transition at the predicted ideal Euler instability. These are indicated by the smooth curves in figure 3, which are the monomodal solutions of equation (30). (ii) The amplitudes coexist with those of lower modes, so that the envelopes of the distributions lie below the ideal upper bounds of the smooth curves. In order for Y n to reach the smooth curve, all other amplitudes would have to vanish, an event that is increasingly improbable as n increases. (iii) The onsets are all more gradual than that of the fundamental, and the amplitudes show random distributions starting at their onsets. Note that the initial randomness of Y 2 is not reflected immediately in Y 1 because Y 2 is small and the sum rule involves the squares of amplitudes.
These properties lead to the following conclusions about both the existence and the nature of biasing, based on equations (28) and (29). Firstly, the existence of any higher mode (n 2) in a static equilibrium configuration requires the presence of biasing in that mode. To show this, first note that formal solutions of equation (29) may be written as where the mode-independent factor is determined self-consistently, The sum is equal to the fractional extension of the beam due to transverse distortion [27], and this partially relieves the compressive strain. We therefore refer to as the 'unrelieved strain.' An unbiased (spontaneous) displacement in any higher mode n would require that = n . But this in turn would require that all lower Y m 's take signs opposing those of theỸ m 's (or vanish, ifỸ m = 0). Clearly, such a solution is unstable. Stable solutions fall into two categories, depending on whether or notỸ 1 vanishes. If it does, then the fundamental mode may exhibit a spontaneous displacement, which is bistable because it may take either sign. Its magnitude must be such that = 1 , and this implies that higher amplitudes are Y n =Ỹ n /(1 − 1 / n ). On the other hand, ifỸ 1 = 0, then equation (35) has a unique stable solution with | | < | 1 |, provided that the biasing parameters are not so large as to violate the sum rule of equation (34).
Equation (35) leads to a second conclusion: A random distribution of Y n 's implies a random distribution ofỸ n 's, albeit a slightly different distribution. The only caveat regards the fundamental mode, for whichỸ 1 may be zero when Y 1 is nonzero. However, if the linear biasing mechanism found at lower aspect ratios also exists at higher aspect ratios, then the Y 1 distribution reflects that ofỸ 1 as well.
Equation (35) leads equally well to a third conclusion: the abruptness of the onsets observed in the higher Y n 's implies similarly abrupt onsets in theỸ n 's. In other words, the absence of a measurable Y n below its critical aspect ratio implies a similar absence of biasing. This means that theỸ n 's themselves must arise spontaneously when their ideal critical aspect ratios are exceeded.
The last conclusion suggests that the formation mechanism of higher mode biasing may be analogous to the linear biasing mechanism operating at lower aspect ratios [28], in which a spontaneous displacement is partially built in and becomes permanent. However, it also raises the question of how a higher mode amplitude Y n , being unstable, could survive long enough to be built in. One possibility is that the etching process releases different parts of the beam at different times. Once any critical segment for the first Euler instability is released, a spontaneous displacement occurs and leaves a permanent bias before much more of the beam is released. The details of the underlying scenario are unknown: the initial displacement might have plastic character, so that a fraction of it would become permanent. Alternatively, the initial displacement could be elastic, but continued etching might effectively anneal the beam, with the same end result. Whatever the microscopic mechanism, the combination of a number of such first Euler instabilities in different locations along the beam would result in a complex shape describable by a superposition of static mode functions. Random distributions would naturally arise from randomness in the exact locations and magnitudes of the separate Euler instabilities. The abruptness of the onsets of higher mode amplitudes is also consistent with the differential release mechanism, because there is a minimum wavelength (2πb/ √ ) in the biasing shape function. Further study on the kinetics of the etching process is needed to resolve the question of the formation of these patterns.

Conclusions
We have described the effects of longitudinal strain on the normal modes of transverse oscillations of a rectangular beam within the usual continuum elastic model of Euler buckling. All of these modes are softened by compression, resulting in a succession of potential Euler instabilities at successively higher aspect ratios subjected to fixed compressional strain. In an ideal symmetric beam, the nth mode instability is realized only if all lower instabilities are suppressed by constraints. We have shown that the critical buckling modes nonetheless provide an appropriate basis for the analysis of static equilibria in the presence of transverse forces or intrinsic asymmetry. We derived a multi-mode Ginzburg-Landau model and used it to interpret recently observed mesoscopic beam shapes.
This model successfully describes the regime of lower aspect ratios, where at most the fundamental mode exhibits large amplitudes and the behaviour is deterministic. It also successfully describes the behaviour at larger aspect ratios, where there is a high degree of randomness in the sample shapes, attributed to similar randomness in the biasing or intrinsic shape parameters. We showed that the observed transitions to beam shapes with higher modes occur at the critical aspect ratios corresponding to the predicted higher Euler instabilities of ideal beams. We argue that these transitions must occur in the biasing parameters themselves. This in turn suggests that the biasing shape function is formed spontaneously but nonuniformly, as the result of separate fundamental Euler instabilities in different segments of the beam.
The model presented here is phenomenological in the sense that it assumes the presence of intrinsic asymmetry or biasing, but it does not account for the origin of such biasing in terms of an underlying microscopic or kinetic theory. Nevertheless, the analysis of experimental data on fabricated mesoscopic beams enables us to characterize biasing quantitatively and provides important clues as to its origins. Future theoretical study should elucidate the formation mechanism and provide a quantitative estimate of the distribution of magnitudes of the biasing parameters. This and future study along these lines will have important implications for the design and operation of mechanical devices down to the nanoscale.