Analytical model for the approximation of hysteresis loop and its application to the scanning tunneling microscope
AAnalytical model for the approximation of hysteresis loop and its application to the scanning tunneling microscope
Rostislav V. Lapshin “Delta” Microelectronics and Nanotechnology Research Institute, 2 Schelkovskoye Shosse, Moscow 105122, Russia
A new model description and type classification carried out on its base of a wide variety of practical hysteresis loops are suggested. An analysis of the loop approximating function was carried out; the parameters and character-istics of the model were defined – coersitivity, remanent polarization, value of hysteresis, spontaneous polarization, induced piezocoefficients, value of saturation, hysteresis losses of energy per cycle. It was shown that with piezo-manipulators of certain hysteresis loop types, there is no difference in heat production. The harmonic linearization coefficients were calculated, and the harmonically linearized transfer function of a nonlinear hysteresis element was deduced. The hysteresis loop type was defined that possesses minimum phase shift. The average relative ap-proximation error of the model has been evaluated as 1.5%–6% for real hysteresis loops. A procedure for definition of the model parameters by experimental data is introduced. Examples of using the results in a scan unit of a scan-ning tunneling microscope for compensation of raster distortion are given.
PACS: 77.80.Dj, 75.60.-d, 07.79.Lh, 07.79.Cz, 07.05.Tp Keywords: hysteresis, hysteresis loop, parametric form, double hysteresis loop, area of hysteresis loop, linearization, piezoscanner, piezomanipulator, scanning tunneling microscope, STM, metrology, simulation, piecewise-linear hysteresis loop I. INTRODUCTION
The law of a closed hysteresis curve would emerge in many physical phenomena such as dielectric hysteresis [polarization curve P = f ( E )], magnetic hysteresis [magnetization curve B = f ( H )], elastic hysteresis [deformation curve ε = f ( F )], and some others. This phenomenon is widespread and important, however a simple analytical equation that could approximate it with a sufficient degree of precision has not existed. Therefore, rather often when analyzing various processes and systems with a hysteresis element, the solution was being searched for either graphically – by using experimental data, or with the help of straight-line approximation of the curve.
1, 2
The first method is inconvenient because of the need for table function representation and low precision of graphics. Low precision of the second method is caused by roughness of piecewise-linear approximation (certainly, if the number of the segments is not too great) and thereto it implies searching for the solutions at several intervals followed by “gluing” them with each other. Other methods of hysteresis loop approximation are known in both polynomial models and integral opera-tors
3, 4 classes. But their usage is limited because of the complexity of hardware support or great time of calcula-tions. This work pursued two objectives. First, to give a description of the suggested model and its properties and characteristics. The part of the paper devoted to this matter is of a wide, general scope inasmuch as the use of the model is considered for analyzing static nonlinearity of hysteresis loops of various types and physical nature en-countered in many scientific instruments. Second, to highlight the ways and manners of practical application of the conclusions obtained from the model description to particular purposes, namely, for description and compensation of nonidealities of piezoceramic manipulator
5, 6, 7, 8, 9, 10, 11, 12 of a scanning tunneling microscope (STM) – nonlinear- nalytical model for the approximation of hysteresis loop ity and ambiguity of static characteristic, piezoceramic creep, thermal drift. II.
DESCRIPTION OF THE MODEL A.
Analytical expression for hysteresis curve family. Formation of the types and their classification
The family of hysteresis loops can be described by a generalized transcendental equation in parametric form as follows ( )( ) ,sin ,sincos aa aaa y nxm by bax = ––= (1) where a is the split point coordinate; b x , b y are the saturation point coordinates; m , n are integer numbers (see Table I); a is a real parameter (- ¥£ a £ + ¥ ). Plus signs in Eq. (1) correspond to a hysteresis curve with saturation points in I, III quadrants (this is the case to be considered below) and minus signs in II, IV. Thus, the curves built with plus and minus signs would be sym-metrical to each other relative to oy axis. Classification of the hysteresis loop types supported by the model is presented in Table I. Each type is ascribed with a status of either main or derivative. The main types are “leaf,” “crescent” (see Ref. 13), “classical” [Figs. 1(a)-(c)]. They are determined by the n coefficient. The m and n ( n >3) powers in Eq. (1) define the steepness Table I. Classification of hysteresis loops. No. Type Type status n m b Equation Figure 1 Leaf Main 1 1 3 5 xy bb arctan ( )( ) aa aaa sin ,sincos y nxm by bax = += p or 1(b)-I 1(b)-II 1(b)-III 3 Classical Main 3 1 3 5 p ( )( ) ,cos ,sincos aa aaa y mnx by abx = -= £ a £ p а а
0… 2 p ( ) ( ) ( )( ) ( ) ( ) ,cossin ,sincos qaqaa qaqaa yxy yxx +-= += q =( p /2)- b а а а ( ) ( )( ) ( ) yx byy bxx –= –= aa aa , 1(e) 6 Bat Derivative from types 1/3/4 а а а ( ) ( )( ) ( ) aa aa yy xx == ~ ,~ 1(f) а Depending on the initial type. . V. Lapshin (I, II, III curves). The derivative types are “tilted classical,” “double loop,” “bat” [Figs. 1(d)-(f)]. They derive from main types or from other derivative types with implementing some extra operations. So, a curve of tilted classical derivative type [Fig. 1(d)], as the tangent in inflection point of the unsplit loop makes an angle b „p /2 with ox axis, can be built by rotating the coordinate system clockwise through the angle q = p /2- b . Thus, using the foregone formula for transformation of Cartesian coordinates with the axes rotated, the following equation is obtained: ( ) ( ) ( )( ) ( ) ( ) ,cossin ,sincos qaqaa qaqaa yxy yxx +-= += (2) where ( ) a x , ( ) a y are the coordinates of the rotated system. When rotating, the split point – a and the saturation point – b change their positions relative to the origin sys-tem. Therefore a preliminary distortion of their coordinates is needed to get these points to coincide with the origi-nals after rotation. Here, the following transformation formulas will be used: q cos aa x = (3) and .cossin ,sincos qq qq yxy yxx bbb bbb += -= (4) (a) (b) (c) (d) (e) (f) Fig. 1. Hysteresis loop types supported by the model: (a) Leaf ( a =0.2; b x =0.6; b y =0.8; m =1, 3, 5; n =1; q =0 (cid:176) ); (b) crescent ( a =0.2; b x =0.6; b y =0.8; m =1, 3, 5; n =2; q =0 (cid:176) ); (c) classical ( a =0.2; b x =0.6; b y =0.8; m =1, 3, 5; n =3; q =0 (cid:176) ); (d) tilted classical ( a =0.2; b x =0.6; b y =0.8; m =3; n =3; q =15 (cid:176) ); (e) double loop ( a =0.1; b x =0.4; b y =0.4; m =3; n =3; q =15 (cid:176) ); (f) bat ( a =0.2; b x =0.6; b y =0.8; m =3; n =3; q =15 (cid:176) ). nalytical model for the approximation of hysteresis loop As a , b x , b y in Eq. (1) are substituted with the corrected value of split constant x a and saturation constants x b , y b , equations for x ( a ) and y ( a ) can be obtained, which are needed in calculations by formula (2). To build the models of more complicated hysteresis loops [e. g., double loops, Fig. 1(e), which are a composi-tion of l/3/4 type] or to simulate thermal drift processes or creep effect, the rectangular coordinate system should be shifted by the value of x , y . Herein, the value of x is a bias voltage, and for creep and thermal drift, the value of y is a certain time function y = f ( t ).
6, 7, 8, 9
Thus, the shifted curve could be described as follows: ( ) ( )( ) ( ) ,, yyy xxx –= –= aa aa (5) where the equations for x ( a ), y ( a ) are determined by the composition base type. To form a model of bat type [Fig. 1(f)], it is sufficient to take the module of y ( a ). The direction of passing along a hysteresis loop is supposed to be determined in the following way.
1, 2
For a coordinate lag system, the movement along the lower half of the curve is associated with the inequality dx / dt >0, and along the upper half with dx / dt <0, the direction obtained being counterclockwise. For a coordinate lead system, the signs in those inequalities would be swapped (clockwise direction). When a parametrically set model is being con-sidered, it is suitable to associate the positive (counter-clockwise) direction with increasing a parameter, and the negative one with decreasing. In the case of an unsplit hysteresis loop ( a =0), it is easy to pass from the parametric record over to an explicit function definition ( ) ., xxnn xy bxbxbbxy +££-= (6) By using expression (6), loops embedded into the limit cycle can be built. The procedure consists of the following (see Fig. 2): for * x b chosen within the (- b x , b x ) interval, the value of * y b is being found by Eq. (6), which lies at the unsplit curve. Then, the embedded loop being assumed, at first approximation, as similar to the limit cycle, the value of split a * = a * x b / b x is being set and the curve is being built using the parameters m , n , q of the limit cycle. In the general case, instead of Eq. (6), any suitable function that passes through the points – b and the origin of coordi-nates can be used. To trace the history of motion, it should be done in the following way. First, the moment is defined when the sign of derivative of the input signal is changing (this moment corresponds to a saturation point of a particular cy-cle). Next, as the value of the a parameter is known at that moment, the parameters – * x b = x ( a ), – * y b = y ( a ) can be calculated. Then, as the law of disposal of particular cycles inside limit one is known, it becomes possible to define Fig. 2. Embedded loops. Determination of spontane-ous polarization value ( a =0.2; b x =0.6; b y =0.8; m =1; n =3; * x b =0.3). . V. Lapshin a * and q * parameters. Finally, the a parameter is assigned the value of p k /2, where k =1, 3, 5, … Here are all the parameters needed to proceed with the motion along the new cycle. Beginning with m =3 and n =1, function (1) cannot be resolved in explicit form y = f ( x ) because it would result in an equation of more than fourth degree. However, an explicit record can be obtained for the function, inverse to Eq. (1), which is yielded by simply swapping x ( a ) and y ( a ) ( )( ) .sincos ,sin aaa aa nxmy bay bx +== (7) That inverse function will serve as the base for building a hysteresis compensation system. An explicit record for Eq. (7) can be given by ( ) ., xxmymynnyx bxbxbbaxbby +££--–= (8) Taking into consideration, on one hand, that the operations of division and square root are not facile in hard-ware implementation, require much time of calculation, and are not typical in digital signal processor (DSP) algo-rithms and, on the other hand, the convenience of definition of the inverse function and the contour pass direction, the parametric form (1) must be admitted to be the most suitable to be used further. B. Analysis of the curve: Coersitivity, remanent polarization, dielectric spontaneous polarization, induced piezocoefficients, value of saturation
Function (1), which was taken as the model of hysteresis, is a continuous, nonlinear, ambiguous, limited one defined on the [- b x , b x ] segment. Besides, this function is periodical with the period 2 p , since the following equali-ties are true: ( ) ( )( ) ( ) ,2 ,2 kyy kxx paa paa += += (9) where k =1, 2, 3, … Owning to the property of periodicity of function (1), use of the model will prove to be most effective with cycle processes, e. g., in scanning systems (see Sec. III.A). Let us find the zeros of function (1). Here, the value of the a parameter should be defined such that y ( a ) be-comes equal to 0, i. e., b y sin a =0, whence a = p k . Substituting this value of a in x ( a ), it will yield x = – a . Note that the physical sense of a zero of function (1) is coersitivity. Now, the coordinates of the – c point are to be defined as the intersection points of the loop and oy axis. The point – c defines the value of remanent polarization. Writing .sin ,0sincos cb ba y nxm –= =+ a aa (10) In the case of m = n , this will result in the following solution: ( ) .1 mx y abbc + –=– (11) The value of hysteresis H y is defined as c / b y (100%); by using expression (11), it can be written as ( ) .1 %100 mxy abH += (12) Note, that the ratio b x / a in expression (12) is none other than 100%/ H x , where H x is “hysteresis” along ox axis. nalytical model for the approximation of hysteresis loop Thus, a correlation is obtained that links the values H y and H x with each other ( m = n ) .1%100%100 = - mxy HH (13) Let us find the coordinate of the c ´ point (see Fig. 2). It is the intersection point of the oy axis and bd tangent drawn through the saturation point of the unsplit curve (6). Note that the point c ´ defines the value of dielectric spontaneous polarization. The equation for the bd straight line is given by y = kx + c ´, where k =tan g . The derivative k = dy / dx of the unsplit curve (6) in the point x = b x is defined by b y /( nb x ). From the equation for the bd line with x = b x , y = b y , taking into account the found k value, the value of c ´ is defined as .11 -=¢ nbc y (14) It will be shown (not rigorously) that function (1) is odd (except for crescent type). By definition, a function y = f ( x ) would be odd provided f (- x )=- f ( x ). In regard to function (1) ( ) ( ) ( ) .sin ,sincos aa aaa y nxm by bax -=- +-=- (15) Using the identities -sin a =sin(- a ), cos a =cos(- a ), making a formal replacement of the - a parameter with z , and considering the oddness of m and n as well as the fact that the sign of the split parameter a in formula (1) defines only the beginning point and the curve pass direction and, therefore, can be chosen as negative here, the following statements can be formulated: ( )( ) .sin ,sincos zz zzz y nxm by bax =- +=- (16) Thus, the right-hand parts of the equations in (16) have no changes beside those in Eq. (1); consequently the func-tion is odd. As the graph of an odd function is symmetrical relative to the origin of coordinates, the information of a half is sufficient to have the whole function set. The property of symmetry of function (1) gives the opportunity, for example at estimating the approximation error (see Sec. II.E), to carry out the calculations at a half of the domain. Let us find first and second derivatives of function (1), which in case of a parametrically set function are given by ( )( ) aa xydxdy && = (17) and ( ) ( ) ( ) ( )( ) , a aaaa x yxyxdxyd & &&&&&& -= (18) where ( ) „= aa ddxx & , ( ) aa dxdx = && , ( ) aa ddyy = & , ( ) aa dydy = && are first and second derivatives of the functions x ( a ) and y ( a ) by a parameter, respectively. For the function being considered ( ) ,cotsintancos aaaaa nxm nbamx +-= & (19) ( ) ( ) [ ] ( ) [ ] ,1cot1sin1tan1cos --+--= aaaaa nnbmamx nxm && (20) ( ) ,cos aa y by = & (21) ( ) .sin aa y by -= && (22) . V. Lapshin In accordance with Eq. (17), the first derivative is defined as aaaa a sincoscossin cos -- +-= nxm y nbam bdxdy (23) and, in accordance with Eq. (18) after simple transformations, the second derivative is defined as ( ) ( ) ( ) .sincoscossin cotsin1coscossin2
311 3111222 aaaa aaaaa -- ++- +- --+-= nxm nyxmymy nbam nnbbmabmmabdxyd (24) Expression (23) allows for calculating the induced piezocoefficients of a hysteresis curve in the coordinates of displacement versus electric field strength as well as the differential magnetic permeance at any point of the curve B = f ( H ). By analyzing expressions (23) and (24), it can be shown that for m =1 function (1) reaches its maximum value in the point (+ b x ,+ b y ) and minimum – in the point (- b x ,- b y ). For all m „
1, the lower and the upper parts of hysteresis loop (1) are monotonously increasing functions, therefore at the edges of the domain, i. e., in the points – b x , func-tion (1) has the greatest + b y and the least - b y values. Thus, the ( – b x , – b y ) points of curve (1) are the saturation points of a hysteresis loop. C. Square of a loop
Physically, the square of a hysteresis loop characterizes the heat losses that cause heating of the material and, therefore, define its efficiency coefficient. To find the square S of a loop, the following integral should be calcu-lated: ( ) ( ) ,21 ∫ -= aaaaa dddxyddyxS (25) then taking into account expressions (1), (19), and (21), writing ( ) ( ) [ ] .sincoscossinsincossincos21
20 11 ∫ -- +--+= p aaaaaaaaa dnbambbbaS nxmyynxm (26) Opening the parenthesis, grouping up the terms, then using De Moivre’s expansion of cos m +1 a and cos m -1 a , and solving the integral for odd n , the hysteresis loop square (26) will be given by ( ) ( ) ( ) ( ) ,2121 ymmmmmmm abCCmCS p -+= ----++ (27) where ( ) [ ] !!! klklC kl -= are binomial coefficients. Thus, the square of leaf and classical hysteresis loops would not apparently depend on the saturation value by x coordinate and is only determined by the split constant and the satu-ration value by y coordinate. Formula (27) is applicable to odd n , but the coefficient n itself does not participate in it explicitly. Therefore, the following theorem can be formulated: the quantity of the heat produced by the piezomanipulator for a cycle is the same for hysteresis loops of both leaf and classical types, provided the values of m , a , b y involved are the same. Since the square of a geometrical figure would be invariant to rotation, by substituting a and b y in formula (27) with their corrected values x a and y b [see Eqs. (3), (4)], a formula can be written to define the square of tilted classical hysteresis loop ( ) ( ) ( ) ( ) ( ) [ ] .12cos2sin2 121 ++ -+= +----++ qqp yxmmmmmmm bbaCCmCS (28) nalytical model for the approximation of hysteresis loop As the terms containing the m coefficient are considered at various values of m in formulas (27), (28), it can be concluded that with fixed values of the rest model parameters the inequality S m =5 < S m =3 < S m =1 is true, which is seen well in Figs. 1(a), (c) (curves I, II, III). Thus, the greater the m parameter is, the less heat the piezomanipulator produces. For example, with m = n =3 ( q =0) the losses of heat will be determined by the hysteresis loop square of S =3/4 p ab y . Note that with the product of a by b y equal to 4/(3 p ), the hysteresis loop square becomes equal to unit. D. Harmonic linearization coefficients. The harmonically linearized transfer function of a hysteresis element
Let the linear part of a nonlinear hysteresis system have the property of a low-pass filter, i. e., transmit the first harmonic with practically no losses and considerably attenuate all the others, which have been produced by a nonlinear element under x ( w t )= A sin( w t ) input harmonic action. When analyzing such a system, the method of harmonic linearization can be applied, in accordance with which the hysteresis function y = f ( x ) would be given in the following form:
1, 2 ( ) ( ) ,ˆ xpAqAqy += w (29) where ( ) ( )( ) ( ) aaap aaap pp dAfAAq dAfAAq cossin1ˆ ,sinsin1 ∫∫ == (30) are the coefficients of harmonic linearization ( a = w t ); p denotes the differentiation operator ( p = d / dt ). Let us define the harmonic linearization coefficients q ( A ), ( ) Aq ˆ for the inverse function (7). This function ob-viously matches the conditions of the task under consideration. At that, the input harmonic signal x ( a ) is given by b y sin a , and f ( b y sin a ) – by a cos m a + b x sin n a (here, instead of A amplitude, the denomination b y , is involved), then in accordance with formula (30) writing ( ) ( ) ( ) ( ) .cossincos1ˆ ,sinsincos1 aaaap aaaap pp dbabbq dbabbq nxmyy nxmyy ∫∫ += += (31) By solving the integrals in Eqs. (31) similarly to Eq. (26) for odd n , the coefficients of harmonic linearization of the inverse function (7) are given by ( ) ( ) ( ) ( ) .2ˆ ,2
211 211 ymmmy yn xnny b aCbq b bCbq ++++ == (32) In accordance with formula (29), the linearized equation for the hysteresis loop inverse function becomes the following: ( ) ( ) ,22 xpb aCb bCy ymmmyn xnn += ++++ w (33) . V. Lapshin and the transfer function W -1 ( b y , s )= Y ( s )/ X ( s ) becomes ( ) ( ) ( ) ( ) ( ) .22ˆ, sb aCb bCsbqbqsbW ymmmyn xnnyyy ww ++++- +=+= (34) Now, the required harmonically linearized transfer function of a hysteresis element (1) is to be found as an in-verse function to Eq. (34), i. e., W ( b x , s )=1/ W -1 ( b y , s ). Note that at that, the input amplitude will not be b y but b x instead. Finally ( ) ( ) ( ) .22 2, ww xnnmmmn ynmx bCasC bsbW ++++ + += (35) By substituting s = j w in Eq. (35) (where j is imaginary unit), the amplitude-phase characteristic of a hysteresis element will be given by ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ,22 22, xnnmmmn ymmmnyxnnnmxx bCaC abCjbbCbWjbW ++++ ++++++ +-== w (36) which only depends upon b x amplitude and does not upon w frequency. Since ( ) ( ) ( ) xxx bqjbqbW ˆ += , from Eq. (36) comes the next ( ) ( )( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ,22 2ˆ ,22 2 xnnmmmn ymmmnx xnnmmmn yxnnnmx bCaC abCbq bCaC bbCbq ++++ +++ ++++ +++ +-= += (37) which are the required harmonic linearization coefficients of function (1). For the case of m = n =3, coefficients (37) are given by ( ) ( ) [ ] xyxx babbbq += and ( ) ( ) [ ] xyx baabbq +-= ; their graphs are shown in Fig. 3. The coefficient ( ) x bq determines the steepness of the averaging line. From formula (37), it is easy to see that with the amplitude of input harmonic signal b x increasing, ( ) x bq approaches zero since in that case nonlinearity (1) gets saturated. The amplitude of first harmonic A ( b x )=| W ( b x )| can be defined by the transfer function (36) as ( ) ( ) ( ) ( ) ( ) ( ) ( )
22 2ˆ xnnmmmn ynmxxx bCaC bbqbqbA ++++ + +=+= (38) and its phase j ( b x )=arg W ( b x ) – as ( ) ( )( ) ( )( ) .22arctanˆarctan
211 211 xnnm mmnxxx bC aCbq bqb ++ ++ -== j (39) For example, the phase shift of first harmonic at the output of classical nonlinear element shown in Fig. 1(c), curve II will make -18.4 (cid:176) . The minus sign in expression (37) for ( ) x bq ˆ [this coefficient is placed by the derivative in formula (29)] as well as the minus sign in Eq. (39) mean that the presence of a hysteresis element results in a phase lag of the output signal beside the input one. As is derived from Eq. (39), the wider the hysteresis loop is (the greater the split pa-rameter a ), the greater the phase shift becomes. By analyzing expression (39), it can be shown that the following theorem would take place: with the same val-ues of m , a , and b x , minimum phase shift would occur at the leaf-type loops, and at that | j m=5 |<| j m=3 |<| j m=1 |. Note nalytical model for the approximation of hysteresis loop here that it is the very type encountered in STM piezoma-nipulators [see Fig. 4(a) and Refs. 5, 9-12]. The harmonically linearized transfer function (35) of a hysteresis element having been defined, the further analysis of an automatic control system should be carried out by employing the existing theory of linear systems in order to define its stability, precision, and quality of tran-sient process. E. Estimates of approximation error
The error of approximating the experimental charac-teristics of hysteresis loops will be estimated by calculat-ing the following quantities (1) Maximal absolute approximation error ( ) ( ) ,max xyxy embx x -=D ££ (40) where y m ( x ) are the model data; y e ( x ) are the experimental data. (2) Maximal relative approximation error ( ) ( ) %,100%100max maxmax0 ee embx yy xyxy x D=-= ££ d (41) where ye by = max is the y coordinate of the saturation point of experimental curve. (3) Average relative approximation error ( ) ( ) .%100 ∫ -= x b emex dxxyxyyb d (42) (4) Quadratic mean approximation error ( s ) ( ) ( )( ) .1 ∫ -= x b emx dxxyxyb s (43) The results of definition of the model approximation error of experimental hysteresis characteristics are presented in Table II, where the type and the parameters of a model are given. F. The procedure of definition of model parameters by limit cycle
For building the model, the six parameters must be extracted from the experimental dependence: a , b x , b y , m , n , q by the following algorithm. (1) By the shape of the curve, its type is being identified (see Table I). If the type proves to be derivative, then the main type is being identified, from which the derivative type has derived. When the type is known, the value of parameter n and the kind of equation become known, also. (2) If the type is identified as tilted classical or as its derivatives double loop or bat, then by the graph of experi-mental dependence, the angle q is being defined. Fig. 3. Graphs of ( ) x bq , ( ) x bq ˆ harmonic linearization coefficients vs b x for the nonlinear hysteresis element ( a =0.2; b y =0.8; m = n =3). . V. Lapshin (3) Provided that in a vicinity of the saturation point the curve behaves the way shown in Figs. 1(a)-(c) graph I, it means that m =1; and if like in graphs II, III (the same figure), then m >1. (4) By the experimental dependence, the coersitivity value a is being defined. (5) By the experimental dependence, the b x and b y coordinates of the saturation point are being defined. If the model is of leaf type, then the angle b can be calculated. (6) If the model type is double loop, then the displacement values b x , b y are being defined by the experimental de-pendencies. (7) By varying the coefficients m (if only m „ n (if only n „ q angle (if only the type identified is not leaf or crescent) and using the estimates of Sec. II.E, the minimal value of approximation error is being achieved. III.
APPLICATION OF THE MODEL
The hysteresis compensation elements (HCE), as they are inserted into the STM control system, are generally shown in Fig. 5. The compensation system described below relating to the class of open-loop systems, its HCE elements are serial nonlinear correcting units which could be realized in either hardware (in analog or digital em-bodiment) or software. (a) (b) (c) (d) (e) Fig. 4. Accuracy of approximation: I – model data; II – experimental data. (a) Displacement of STM PZT X, Y piezomanipulators vs applied voltage (leaf; a =32.6; b x =300; b y =955; m =3; n =1; b =72.6 (cid:176) ; H y » 〈 d 〉 =1.5%); (b) ceramic polarization vs electric field strength (Ref. ) (tilted classical; a =54; b x =130; b y =36.4; m =1; n =3; q =27.7 (cid:176) ; 〈 d 〉 =6%); (c) ceramic polarization vs electric field strength (Ref. ); (tilted classical; a =90; b x =340; b y =39.6; m =5; n =3; q =30 (cid:176) ; 〈 d 〉 =4.1%); (d) ceramic polarization vs electric field strength (Ref. 16) (leaf; a =180; b x =195; b y =63.8; m =1; n =1; b =18.1 (cid:176) ; 〈 d 〉 =3.9%); (e) ceramic polarization vs electric field strength (Ref. 16) (tilted classical; a =122.5; b x =341; b y =23.2; m =1; n =3; q =14 (cid:176) ; 〈 d 〉 =2.7%). nalytical model for the approximation of hysteresis loop In order to compensate the hysteresis effect in all the three ways of realization, the function inverse to Eq. (1) will be employed so that to have the identity f ( f -1 ( X )) ” X work (here, -1 designates an inverse function). In terms of the system being described, the source function f is the piezomanipulator static characteristic approximated by Eq. (1), and the inverse function f -1 is function (7) synthesized by the HCE element. The inverse function argument X is a control signal, for example “XIn” scan signal, which undergoes preliminary distortion in the HCE and then is applied, through a high voltage amplifier (HVA), to the piezomanipulator. Thus, the resulting displacement of the manipulator will correspond to the input function of the scan. A. Compensation of raster distortion in the STM scan unit
The flow chart of the hysteresis compensation unit is shown in Fig. 6(a). The scheme is built on analog ele-ments and consists of the following units: GEN sinusoidal generator; a channel of MUL multipliers that carries out the operation of raising the input oscillation sin( w t ) to m and n powers; PS phase-shifting element that converts the sin m ( w t ) oscillation into cos m ( w t ) by shifting it by a quarter of the period; AMP2, AMP3 operational amplifiers with the gains K = a and K = b x , respectively; SUM summing amplifier; AMP1 amplifier with the gain K = b y ; COMP comparator, and AS analog storage unit. The last three components of the scheme are destined to set the accor-dance between the “XIn” input scan signal and the pa-rameter a = w t by successive approximation [such em-bodiment is simpler than an explicit calculation of the function arcsin( X ( t )/ b y )]. In fact, the scheme implements the calculation of the inverse function (7). Its operation principle is clear with Table II. Results of calculation of the approximation error. Source of experimental data and figure Type of model and its parameters: a ; b x ; b y ; m ; n ; q (cid:176) or b (cid:176) D d (%) 〈 d 〉 (%) s Fig. 4(a) Leaf 32.6; 300; 955; 3; 1; 72.6 49.5 5.2 1.5 20.5 Ref. 15, Fig. 4(b) Tilted classical 54; 130; 36.4; 1; 3; 27.7 5.4 14.8 6.0 2.8 Ref. 16, Fig. 4(c) Tilted classical 90; 340; 39.6; 5; 3; 30 3.3 8.3 4.1 1.9 Ref. 16, Fig. 4(d) Leaf 180; 195; 63.8; 1; 1; 18.1 9.0 14.1 3.9 3.2 Ref. 16, Fig. 4(e) Tilted classical 122.5; 341; 23.2; 1; 3; 14 1.7 8.1 2.7 0.9 Fig. 5. HCE hysteresis compensation elements in-cluded in STM control system. . V. Lapshin no extra elucidation, to only point out that, to ensure the condition of in phaseness, the SCAN unit forming the scan works out an “XStart” signal at the beginning of each half-period of the scan, which starts the waiting generator GEN so that it would generate the proper half of the sine function accordingly to the sign of the scan signal deriva-tive [the probing frequency of GEN is much greater than the frequency of the scan signal X ( t )]. Pay attention to the triangle (not ramp) impulses of the scan signal shown in Fig. 6(a), the circumstance pointing out to operation with no idle stroke. Let us show now what the scheme of the HCE element will look like in case of a circular scan microscope. Note that a circular scan is more preferable beside a raster scan for the following reasons: the design of X and Y manipulators proves to be symmetrical, an even employment of the manipulators occurs during the operation, the scan grows smoother and to those the possibility appears of getting the additional features of the microscope which have been described in paper 17. As a disadvantage of a circular scan, there can be admitted the necessity of con-verting the points of the circle into a conventional rectangular display window when visualizing. To have a circular scan, it is necessary to apply a sin( w t ) voltage to the X manipulator and cos( w t ) to Y. Let us, instead of the triangle voltage [see Fig. 6(a)], mentally apply the X ( t )= b y sin a voltage, where a = w t . It can be easy to see that the chain of units – AMP1 amplifier, COMP comparator, and AS analog storage unit – becomes out of use and, therefore, can be excluded from the scheme. Now, instead of the triangle voltage Y ( t ), let the voltage (a) (b) Fig. 6. Flow chart of analog device for hysteresis compensation for the cases of: (a) raster scan; (b) circular scan. nalytical model for the approximation of hysteresis loop Y ( t )= b y cos a be applied (structurally, X and Y channels are identical). Since a cosine function would outstrip a sine at p /2, formula (7) can be given by y ( a )= a cos m ( a + p /2)+ b x sin n ( a + p /2), whence the inverse function can be written in another record ( )( ) ,sincos ,cos aaa aa mnxy aby bx -== (44) which can be used abreast with Eq. (7). From the last reasonings, it can be concluded that the circular scan genera-tion scheme with simultaneous hysteresis compensation will take the outline presented in Fig. 6(b), where SUB is an operational amplifier connected in a differential mode. If comparing Figs. 6(a) and (b), the decrease of hardware expenditures can be revealed with the circular scan implementation; note that the parameter b y becomes unneces-sary. The circular scan scheme described here would serve as the base for building cycloid and spiral scans. The units of hysteresis compensation described above are the simplest ones and, generally speaking, do not permit us to work with loops of derivative types or various image sizes either to shift the image window along the scanning field or to implement vector access to a surface point. A hardware realization of the features mentioned above would result in a substantial complication of the equipment. The way out could be to build a digital system or to compute the model by a program. A digital hysteresis compensation system will be designed by using the structures and the conclusions yielded while synthesizing the analog system. So, for the schemes presented in Figs. 6(a), (b), the analog multipliers MUL and the operational amplifiers AMP1/2/3/4 should be replaced with digital multipliers; the summing operational amplifier SUM and the differentiating amplifier SUB – with an arithmetic-logic unit; the analog storage unit AS – with a strobed register with a digital-to-analog converter (DAC) connected up to its output; the sinusoidal generator GEN – with a read-only memory (ROM) scheme with a sine table written down in it; the analog comparator COMP – with a digital one. A ring counter should be connected to the ROM address inputs. The function of the phase-shifting element PS consists of shifting the value of the address worked out by the ring counter so as to skip exactly a quarter of the period of sine in the ROM table. At program realization of the model, the functions cos m a , sin n a are stored in the computer data memory as some table structures. The scan voltage can be transformed into the corresponding a value by either successive approximation (similar to the operation of the unit chain AMPl-COMP-AS) or immediately by an arcsin table also kept in the data memory, or by directly calculating the arcsin by the foregone identity ( ) ( )( ) ,arcsin ∫ -= t yy tfb tdfb tf where f ( t ) is a scan function (| f ( t )|< b y ). At a half-period of the triangle scan, f ( t ) is given by kt . So, by passing from the integral over to a sum, the following expression for a is obtained ,1
10 2222 ∑ -= D-D=
Ni y itkbtk a (45) where N is the number of the samples taken; D t = t / N is the time discrete. The calculation of the model could be fulfilled at a universal processor but most effectively this task could be resolved at a DSP because the operation AB i + C (where A , B , C are real numbers) is encountered in the calculations which is typical of DSP algorithms. Here, A = a , B i =cos m a i , C = DE j + F , and D = b x , E j =sin n a j , F =0. The indexes i and j . V. Lapshin refer to the addresses of the memory cells where the tabu-lated values of cos m a i and sin n a j are contained, respec-tively. The flow chart of a hysteresis compensation digital unit is shown in Fig. 7. This unit is realized on a soft-hardware base and ensures working with all the hysteresis types and performing a scan along an arbitrary trajectory so as to support any scan type, to carry out rotation of the scan window around oz axis (in order, e. g., to reduce moire distortion) as well as to vary the window size and shift it within the scan field. A ring programmable counter PCT2 is involved in the scheme which, in the rate defined by the “Clock” sig-nal and the division coefficient “DC,” generates the address “Addr2” and the read signal “Read2” for the Dual-Port RAM (random-access memory) scheme with the model data written in it. From the RAM output, the “MData2” code, which corresponds to the preliminary distorted current value of raster voltage, comes into the DAC from where the signal, after having passed through a high voltage amplifier, is applied to the manipulator. Another RAM port is intended for writing the model data “MData1” calculated by the microcomputer processor. The use of a dual-port RAM permits us to get the processes of calculating the model and generating the control signal for the manipulator coincided in time and, therefore, to increase the unit fast acting. Built on a dual-port RAM base, the scan subsystem is capable of being transformed from synchronous into asynchronous. To do that, it is sufficient, instead of the “Clock” signal, to apply a signal pointing out to readiness of the data in the tunnel junction stabilization system. Besides, if that system also uses a dual-port RAM as a Z-Buffer, then for the next memory cell sampling, the “Addr2” signal generated by the PCT2 counter can be used. Thus, there can be seen a good mutual coordination in work between the tunnel junction stabilization system and the scan signal generation system. It is appropriate to note that when using the scan modes that require frequent changing the model parameters a , b x , b y (e. g., in order to change the image size), in a certain operation time, some error will have been accumulated in the model. At the moment when it reaches a certain threshold defined by the admissible approximation error, a compulsory correction must be done: for the model it is assigning p /2 to the a parameter and for the piezoelement – applying the saturation voltage to it so as to ensure it set at a fixed point of the limit characteristic – the saturation point. IV.
EXPERIMENTAL RESULTS
The surface image of a test pattern is shown in Figs. 8(a), (c), which was taken by a STM when scanning with idle stroke and without, accordingly. As the test pattern, a diffraction grating with 0.3 µm period coated with gold was used. Since the grating structure changes along one direction only and scan piezomanipulators, because of hysteresis, introduce distortions along two directions, then for reflecting these distortions on the image taken, the test object was fitted up so that its stripes would make an angle of some 45 (cid:176) with the manipulator X axis. On the images, the distortions of the test pattern are seen well, which consist in curving the stripes and chang- Fig. 7. Flow chart of digital device for hysteresis compensation (soft-hardware realization). nalytical model for the approximation of hysteresis loop ing their width from one stripe to another [Fig. 8(a), cf. Ref. 11], splitting the stripes and formation of a double-sided comb-like structure [Fig. 8(c)], as well as parallel stripes looking like divergent ones [Fig. 8(e), cf. Ref. 12]. Figs. 8(b), (d) present the corrected images of the same area that were taken using the hysteresis compensation system shown in Fig. 7 [the model parameters were extracted from the data presented in Fig. 4(a)]. A visual com-parison of Figs. 8(a) and (b) [(c) and (d)] will show that nonlinear distortions caused by hysteresis of manipulator piezoceramic are practically removed. Hysteresis loops with parameter m equaling to unit would have a section of negative derivative [see Fig. 4(d)]. The behavior of piezomanipulator at that part of the curve is similar to creep: on achieving e point the voltage be-gins to decrease, the displacement still increasing for a certain time until the point b is reached. [The parameter a that corresponds to e point can be found by setting the denominator in Eq. (23) equal to zero.] From works 7, 8, it is (a) (b) (c) (d) (e) Fig. 8. 128 ·
128 STM scan of the same 1 · area of test pattern surface (diffraction grating with 0.3 µm period coated with gold) (a), (c) Surface image data obtained by scanning with and without idle stroke, accordingly (the test object has been rotated by 45 (cid:176) counterclockwise), the hysteresis compensation system turned off. (b), (d) Sur-face image data corresponding to (a), (c) obtained with the hysteresis compensation system turned on. (e) Surface image data obtained by scanning with idle stroke (the test object has been rotated by 45 (cid:176) clockwise). After compen-sation of hysteresis, the image will look like the one shown in (b). . V. Lapshin known that creep would cause distortions at the image boundaries where the tip motion is reversed (e. g., straight lines would get hooks at the ends). Distortions of the same kind would appear when a piezomanipulator with the hysteresis loop said above is applied. Thus, the resulting picture will contain distortions of true creep and hysteresis mixed up together. V. DISCUSSION
The approximating model suggested belongs to the real time model class. It allows for compensation of the STM piezomanipulator nonidealities such as nonlinearity and ambiguity, which makes it possible to get rid of the image distortions (see Fig. 8 and Refs. 10-12), which is especially important with great scans, to refuse the idle stroke in the raster scan, and consequently to reduce the scanning time; in addition, the model supports the mode of direct access to a desired point. Therein, the reconstruction of the true image is available for either preprocessing (preliminary distortion of the scan trajectory) or post-processing of the distorted image taken. As was shown in Sec. II.C, the Q quantity of the heat produced by the piezomanipulator for a hysteresis cycle t determined by the system timer or by the scan frequency can be easily calculated, for instance, by formula (28). The values Q and t might serve as the input data for some model considering the heat exchange proceeding among the STM construction elements. After having defined the difference in the manipulator temperature D T by using such model, the difference in the manipulator length D l can be found by the formula , Tll
D=D b (46) where b is the heat expansion coefficient (1/K); l is the initial manipulator length (m). Calculated in such a way, the thermal drift D l can be compensated in the scan unit by properly shifting the scan window. Thermal drift compen-sation is especially important in a scan device since there is not any active servo system mounted in it. Note that the STM heat processes, depending much upon construction features of the microscope and the materials used, are of a rather complicated nature and deserve a special investigation which comes out of the limits of this paper. The main advantage of the model, its simplicity, is called forth by no need of calculating the parameters – they are taken right from the experimental curve, neither the inverse function is to be calculated. To that, the model pa-rameters make clear physical sense and have a simple geometrical interpretation. A distinctive feature of the model is the possibility of comparing different hysteresis loop types with each other (see theorems in Secs. II.C, II.D). As a drawback of the model, there could be mentioned a restriction on the approximation accuracy to be achieved, i. e., its dependence upon the particular shape of the hysteresis curve. Though, the reverse task could be apparently resolved, the piezoceramic hysteresis curve being adopted to the desired model by changing the techno-logical parameters, namely, the chemical composition, the baking and cooling conditions, the mechanical influ-ences, etc. Owing to application of the model, it becomes possible to linearize the STM tunnel junction stabilization sys-tem contour by methods similar to those described in Sec. III.A, which allows us to prevent the appearance of auto-oscillations and to get rid of distortions in Z direction, although the method that was used in Ref. 18 must be ad-mitted as the most acceptable solution here. When measuring the frequency characteristic of Z manipulators, an auxiliary piezoelement is often employed, which is supposed to modulate the tunnel junction at a harmonic law. With this method, if a sine voltage is applied to the auxiliary piezoelement, the law of its mechanical displacement would not be exactly the sine because of the presence of hysteresis. The model for approximation of hysteresis loops would allow to increase the precision of nalytical model for the approximation of hysteresis loop this method. To do that, the voltage applied to the auxiliary piezoelement must match x ( a ) by formula (1). In that case, the displacement y ( a ) will be a harmonica1 function. The model described may prove useful with the tasks of imitation modeling as well as in engineer calculations of nonlinear control systems containing hysteresis elements. ACKNOWLEDGMENTS
I want to thank Oleg E. Lyapin, Valery V. Efremov, Oleg D. Cnab, Vladimir N. Yakovlev, and Oleg V. Obyedkov for their helpful advice and discussions.
APPENDIX
Beside smooth loops, the model suggested can be implemented for description of piecewise-linear loops as well (see Refs. 1, 2). To obtain a set of piecewise-linear hysteresis loop primitives, in formula (1) there must be used some piecewise-linear functions instead of sin a and cos a . They could be, for instance, trapeziumlike pulses with unit amplitude ( ) ( ) ( ) ( ) ( )( ) ,4trptrp ,,rect1,rect124trp += -+- --= ∑ ¥= T iiTidD sc i iis aa aaaa (A1) where the subscripts s and c refer to sine and cosine, respectively; d and D are the upper and lower bases of the trapezium, respectively; T = d + D is the period of pulses; rect ( a , i )=1( a +( D - d )/4- iT /2)-1( a -( D - d )/4- iT /2) and rect ( a , i )=1( a -( D - d )/4- iT /2)-1( a -( D - d )/4- d - iT /2) are the i th rectangular pulses; 1( a , i ) is the i th unit step function. Triangle functions [tri s ( a )= ( ) a sd trplim fi , tri c ( a )=tri s ( a + T /4)] or rectangle functions [rect s ( a )= ( ) a sDd trplim fi , rect c ( a )=rect s ( a + T /4)] could be also used to that purpose. References E. P. Popov,
The theory of nonlinear automation control systems (Nauka, Moscow, 1988) (in Russian). E. P. Popov and I. P. Paltov,
Approximation methods for investigation of nonlinear automatic systems (Pizmatgiz, Moscow, 1960) (in Russian). B. Jankovi ć , Proceedings of the 5th International Conference on Nonlinear Oscillations, Kiev, 1969, vol. 4, p. 503. R. Bouc, Proceedings of the 5th International Conference on Nonlinear Oscillations, Kiev, 1969, vol. 4, p. 100. L. E. C. van de Leemput, P. H. H. Rongen, B. H. Timmerman, and H. van Kempen, Calibration and characteriza-tion of piezoelectric elements as used in scanning tunneling microscopy, Rev. Sci. Instrum. , 989 (1991). S. Vieira, The behavior and calibration of some piezoelectric ceramics used in the STM, IBM J. Res. Dev. , 553 (1986). E. P. Stoll, Restoration of STM images distorted by time-dependent piezo driver aftereffects, Ultramicroscopy - , 1585 (1992). E. P. Stoll, Correction of geometrical distortions in scanning tunneling and atomic force microscopes caused by piezo hysteresis and nonlinear feedback, Rev. Sci. Instrum. , 2864 (1994). D. W. Pohl, Some design criteria in scanning tunneling microscopy, IBM J. Res. Dev. , 417 (1986). . V. Lapshin O. Nishikawa, M. Tomitori, and A. Minakuchi, Piezoelectric and electrostrictive ceramics for STM, Surf. Sci. , 210 (1987). R. C. Barrett and C. E. Quate, Optical scan-correction system applied to atomic force microscopy, Rev. Sci. In-strum. , 1393 (1991). L. Libioulle, A. Ronda, M. Taborelli, and J. M. Gilles, Deformation and nonlinearity in scanning tunneling mi-croscope images, J. Vac. Sci. Technol. B , 655 (1991). S.-Y. Liu and I-W. Chen, Fatigue deformation mechanisms of zirconia ceramics, J. Am. Ceram. Soc. , 1191 (1992). G. A. Korn and T. M. Korn,
Mathematical handbook (McGraw-Hill, New York, 1968). S. Hirano, T. Yogo, K. Kikuta, Y. Araki, M. Saitoh, and S. Ogasahara, Synthesis of highly oriented lead zircon-ate-lead titanate film using metallo-organics, J. Am. Ceram. Soc. , 2785 (1992). C. D. E. Lakeman and D. A. Payne, Processing effects in the sol-gel preparation of PZT dried gels, powders, and ferroelectric thin layers, J. Am. Ceram. Soc. , 3091 (1992). D. W. Pohl and R. Möller, “Tracking” tunneling microscopy, Rev. Sci. Instrum. , 840 (1988). R. V. Lapshin and O. V. Obyedkov, Fast-acting piezoactuator and digital feedback loop for scanning tunneling microscopes, Rev. Sci. Instrum.64