Effective model for a short Josephson junction with a phase discontinuity.

Effective model for a short Josephson junction with a phase discontinuity.

E. Goldobin gold@uni-tuebingen.de Physikalisches Institut and Center for Collective Quantum Phenomena in LISA, Universität Tübingen, Auf der Morgenstelle 14, D-72076 Tübingen, Germany    S. Mironov sermironov@rambler.ru [ LOMA, UMR-CNRS 5798, Université Bordeaux, 351, cours de la Liberation, F-33405 Talence Cedex, FRANCE    A. Buzdin a.bouzdine@loma.u-bordeaux1.fr LOMA, UMR-CNRS 5798, Université Bordeaux, 351, cours de la Liberation, F-33405 Talence Cedex, FRANCE    R.G. Mints mints@post.tau.ac.il The Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel    D. Koelle koelle@uni-tuebingen.de    R. Kleiner kleiner@uni-tuebingen.de Physikalisches Institut and Center for Collective Quantum Phenomena in LISA, Universität Tübingen, Auf der Morgenstelle 14, D-72076 Tübingen, Germany
July 14, 2019
Abstract

We consider a short Josephson junction with a phase discontinuity created, e.g., by a pair of tiny current injectors, at some point along the length of the junction. We derive the effective current-phase relation (CPR) for the system as a whole, i.e., reduce it to an effective point-like junction. From the effective CPR we obtain the ground state of the system and predict the dependence of its critical current on . We show that in a large range of values the effective junction behaves as a Josephson junction, i.e., has a unique ground state phase within each interval. For and near the middle of the junction one obtains a junction, i.e., the Josephson junction with degenerate ground state phase within each interval. Further, in view of possible escape experiments especially in the quantum domain, we investigate the scaling of the energy barrier and eigenfrequency close to the critical currents and predict the behavior of the escape histogram width in the regime of the macroscopic quantum tunneling.

current-phase relation
pacs:
74.50.+r, 85.25.Cp

Present address:] Moscow Institute of Physics and Technology, 147700, Dolgoprudny, Russia

I Introduction

Recently a lot of attention is attracted to Josephson junctions (JJs) with an unconventional current-phase relation (CPR)Golubov et al. (2004); Buzdin (2005). In particular, JJsKrive et al. (2005); Buzdin (2008); Reynoso et al. (2008); Alidoust and Linder (2013); Goldobin et al. (2015); Mironov et al. (2015) and JJsMints (1998); Mints et al. (2002); Buzdin and Koshelev (2003); Goldobin et al. (2007); Sickinger et al. (2012); Bakurskiy et al. (2013); Heim et al. (2013) and their combinationsGoldobin et al. (2015), proposed and/or demonstrated recently, show non-trivial physicsGoldobin et al. (2013a) and have potential for applications in the classicalOrtlepp et al. (2006); Khabipov et al. (2010); Goldobin et al. (2013b) and the quantum domainsFeofanov et al. (2010) similar to JJs. Here, JJs are defined as JJs having a unique ground state phase (single Josephson energy minimum situated at) within each phase interval, while JJs (sometimes denoted also JJ) have a doubly degenerate ground state phase (double-well Josephson energy with minima at) within each interval.

Currently, the classical properties of JJs made of a short 0- JJ are understood rather wellGoldobin et al. (2007); Sickinger et al. (2012); Goldobin et al. (2013b). For example, JJs have two critical currents and corresponding to the escape of the phase from and wells. In our group we are starting investigation of quantum properties of such JJs. The first step in this direction could be an observation of the macroscopic quantum tunneling (MQT) of the phaseVoss and Webb (1981); Bauch et al. (2005); Li et al. (2007) out of both and wells of the Josephson energy profile. For this purpose, one, usually, measures the phase escape statistics, by sweeping the bias current at a constant rate and measuring the exact value of the switching current many times. Assuming that in JJ at low temperatures (low damping) the initial state ( or ) is randomGoldobin et al. (2013a), the switching current histogram should have two peaks, each of them just below corresponding critical current . The widths of each histogram peak usually (when the damping is small) decreases with decreasing temperature . However, is expected to saturate at some value for temperature below some . Such behavior is usually interpreted as a transition from the regime of the thermal activation of the phase over the barrier to the regime of the MQT of the phase through the barrier. However, it is necessary to show that the observed is not related to the (noise in the) experimental setup and other trivial reasons. Usually, in such experiments one introduces some extra tuning parameter, e.g., a magnetic field, which allows to demonstrate that the setup is able to measure the histograms that are more narrow than . Simultaneously, for the MQT experiment with JJ, it would be advantageous to have a tuning parameter, which provides a continuous transition between (or ) JJ and a conventional 0 JJ, whose physics is well studied.

For this purpose, we propose to use a short 1D conventional 0 JJ, equipped with a pair of tiny current injectors. By sending a current from one injector to the other, we can create a discontinuity of the Josephson phase at some point along the JJ, where injectors are attachedGoldobin et al. (2004); Gaber et al. (2005); Ustinov (2002); Malomed and Ustinov (2004). If , the system is similar to a superconductor-insulator-ferromagnet-superconductor (SIFS) 0- JJ, which becomes a JJ, if parameters are chosen correctlyGoldobin et al. (2011); Lipman et al. (2014). However, since is adjustable, one can in situ tune the junction from a 0 JJ to a JJ and also study all the states in between. Tuning one can also affect the widths of the hystograms.

The aim of this work is to develop a theoretical model for a short JJ with discontinuity of the phase and to predict or interpret the results of MQT experiment such as the one outlined above. Namely, we derive an effective (averaged) CPR for a short JJ with a phase discontinuity and obtain experimentally relevant quantities, such as the critical current or the escape histogram width as a function of .

The paper is organized as follows. In Sec. II we introduce the model and present the averaged CPR and the averaged Josephson energy derived in details in appendix A. In Sec. III we obtain several experimentally relevant dependences such as the ground state phase, critical current and escape-related characteristics as functions of . Sec. IV concludes the work.

Ii Model

We consider a short JJ of length ( in units of Josephson length ). The Josephson phase changes along the coordinate (). At () there is a discontinuity of the Josephson phase, created, e.g., by a pair of tiny (in theory infinitesimal) current injectorsGaber et al. (2005). The junction is biased by a uniform current density (given in the units of the critical current density). Our aim is to derive an effective (averaged over the JJ length) current-phase relation for this system, i.e., , where

(1)

is the average phase across the JJ. It is that is actually measured, if one considers the system described above as a black box with two electrodes.

In appendix A we derive the averaged CPR of the system under question by using the perturbation theory up to the second order in , treating (the half-length of the JJ) as a small parameter. It is convenient to write the resulting averaged CPR as function of the phase , which is related to the average phase across the JJ as

(2)

where . The averaged CPR can be written as

(3)

where

(4)

is the 0-th order result of the perturbation theory, the first order gives no correction, and

(5)

is the second order correction in the framework of the perturbation theory, and is introduced to make some formulas below more compact. The third order correction is zero and the terms and smaller will be neglected.

The effective Josephson energy of the system is an integral of the effective CPR (3) and is given by

(6)

where

(7)
(8)

Iii Results

iii.1 Comparison with the previous results

For we expect the result given by Eqs. (3), Eq. (4) and Eq. (5) to be similar to those previously obtained for a asymmetric 0- JJGoldobin et al. (2011). To compare both results, first, we have to convert the phases. In Ref. Goldobin et al., 2011 the phase is continuous, while in our case it is has a discontinuity. Instead of the discontinuous phase we can introduce the continuous phase , which behaves exactly like , but without a -jump, i.e.,

(9)

Then

(10)

It is that is used in Ref. Goldobin et al., 2011 (it is denoted as there). Rewriting our effective CPR in terms of and taking , we obtain

(11)

The quantities such as , and from Ref. Goldobin et al., 2011 can be expressed in terms of quantities used here as

(12)
(13)
(14)

By substituting this into expression (18) of Ref. Goldobin et al., 2011 and taking into account the definition of , see Eq. (17) of Ref. Goldobin et al., 2011, we arrive at the CPR (11) derived here. Thus, for the result of Ref. Goldobin et al., 2011 is reproduced exactly. Note however that in Ref. Goldobin et al., 2011 the small parameter is the deviation of the phase from its average value, while in our case the small parameter is . Although, they are related (one expects small deviations for small ), this relation is not straightforward.

In terms of variables used here

(15)

If then we have a JJ. This means that only for one obtains a JJ at .

iii.2 Ground state phase and the critical current

Figure 1: (Color online) The and curves calculated using Eqs. (17) and (18) for (a) , (b) and (c) .

In the 0-th approximation the averaged CPR (4) can be rewritten as

(16)

where

(17)

is the maximum supercurrent. The critical current measured in experiment is . has maxima equal to 1 at ( is any integer) and minima equal to at , see Fig. 1.

In Eq. (16) the ground state phase , where

(18)

and the function returns the argument (phase angle) of a complex number . Obviously, the CPR given by Eq. (16) corresponds to a JJKrive et al. (2005); Buzdin (2008); Reynoso et al. (2008); Alidoust and Linder (2013); Goldobin et al. (2015); Mironov et al. (2015). Some examples of and dependences are shown in Fig. 1. For large asymmetry the modulation of is not as deep as for . The phase shift changes from 0 to as changes from 0 to . It is positive for and negative for .

When , the critical current given by the 0-th order formula (17) vanishes close to and one has to take into account the next (second) order corrections given by Eq. (5). This happens for , see the discussion after Eq. (15).

Figure 2: (Color online) The ground state phase . Comparison of the approximate dependence given by Eq. (20) (black line) with the exact dependence calculated by numerically solving , see Eq. (3), for each value of (symbols). (a) , (b) , (c) and the region close to zoomed. In all plots the regions with a positive slope (green) correspond to a stable solution (energy minimum, where ), while the regions with a negative slope (red) correspond to an unstable one (energy maximum, where ).

Next, we consider the second order approximation. The ground state phase is a solution of for this , i.e.,

(19)

This equation can be solved only numerically, see Fig. 2. It can be seen that multiple solutions appear in the vicinity of , see Fig. 2(b,c). To find the approximate analytical expression describing them we take (). Then we expand Eq. (19), up to the first order in and solve it for . Finally we obtain an approximate value of for any given ground state phase ,

(20)

i.e., the inverse of the ground state phase . This approximation is also shown in Fig. 2. One can see that approximation given by Eq. (20) is very good in the whole range of . Note, that the appearance of three solutions (two stable and one unstable) out of one near is a result of the competition of the -term with the -term in Eq. (20). From Eq. (20) one can figure out that the multiple solutions appear for , i.e., , which is in agreement with the discussion after Eq. (15). We note that if is so small it can be neglected in the definition of , so that for when the second order approximation becomes important, . If , one can omit the term in Eq. (20) and end up practically with expression (18) for from the 0-th order approximation. The difference is that (18) gives only the stable branch, cf. the ground state phase shown in Fig. 1 (only stable branch) and Fig. 2 (both branches).

From Eq. (20) one can find the range of where the double ground state exists, i.e., the points and in Fig. 2 where . We obtain that

(21)

which lay symmetrically with respect to . It follows from Eqs. (21) that the bifurcation point, where the curve switches from one stable to three solutions (two stable and one unstable), corresponds to , i.e., at , which is again in agreement with the result obtained directly from Eq. (20). The range of around where two stable solutions exist, is found by substituting Eq. (21) into Eq. (20).

Looking at the Fig. 2 one sees that the ground state phase has the values symmetrically placed around only at . From Eq. (20) they are given by

(22)

For the symmetry is brocken because the corresponding double-well potential Eq. (6) becomes asymmetric (one well is deeper than the other) relative to . The real (measurable) ground state phases are given by and at are symmetric with respect to the phase . Recalling that for doubly generate state to occur one need very small , the shift from is small. Thus, such a JJ can be called a JJ, where, at and .

Figure 3: (Color online) Examples of the dependence for . (a) , global behavior of ; (b) , zoom of the region of interest near ; (c) , zoom of the region near . Thick (red) lines/symbols show obtained by directly solving Eq. (23) numerically to find all and then calculating from Eq. (3). Thinner black lines/symbols correspond to the approximation given by Eq. (25). Gray dashed lines show , see Eq. (17), for the same parameters.

To find the critical current(s) in the second order approximation for each we search for an extremum of with respect to . It takes place at for which

(23)

Here and below the prime denotes . This equation can be solved for only numerically to find several (up to 4) for each value of . Then we substitute each of these into Eq. (3) to find . The result is presented in Fig. 3. The global behavior is defined mainly by , i.e., . However, near , where vanishes, results in a bistability and in the formation of a -like intersection of the branches. Such behavior is typical for a JJ made of 0 and partsGoldobin et al. (2011); Lipman et al. (2014); Goldobin et al. (2015).

Similar to the case of the ground state phase, one can find an approximate expression for near . By substituting () into Eq. (23), Taylor-expanding it up to terms , and expressing , we obtain the critical value of corresponding to for given .

(24)

To calculate the critical current, we substitute Eq. (24) into Eq. (3), which was preliminary expanded near up to . We obtain

(25)

By sweeping in the range , we can now calculate and and make a parametric plot of vs. , see Fig. 3. The agreement with direct numerical calculations near is excellent, see Fig. 3(b,c). The deviations become noticeable as approaches 0 or , see Fig. 3(a).

This approximate analytical expression for allows us to calculate some key features in the plot. For example, one can find out the value of (and ), for which the branches meet each other, see points A and B in Fig. 3. The analysis of the dependence (24) shows that this happens when . Differentiating Eq. (24) we obtain the following equation for .

(26)

This cubic equation with respect to has only one suitable root, which (after some lengthy algebra) can be expressed as

(27)

where may be explicitely written as

(28)

It turns out that for typical parameters corresponding to a bistable case (small ), the value of changes from to . Then by expanding Eq. (28) as

(29)

we get . In fact, this limit of small corresponds to neglecting the term in Eq. (26), so that one obtains right away from Eq. (26). After finding , the value of can be found as from Eq. (24). Using the approximation (29), i.e., in the worst case neglecting in comparison with (accuracy ) in Eq. (26), we can write

(30)

iii.3 Energy barrier

We consider the thermal escape or the quantum tunneling of the phase out of the potential well, when the bias current . Since our model reduces the system to an effective point-like JJ, for calculation of the escape rate one can use standard thermal or quasi-classical quantum formulas. In these formulas, the key parameters are the barrier height and the eigenfrequency . The aim of this section is to obtain the expressions for them.

In general, we proceed as follows. Given the Josephson energy profile , the total potential energy of the biased JJ can be written as a tilted potential

(31)

The static solution(s) correspond(s) to

(32)

In essence this is a CPR. The critical current is reached for , when , i.e.,

(33)

From here one can, in principle, find (one or more) values of . Imagine that we have found all values of . Then, the value of the critical current is found from Eq. (32), i.e.,

(34)

Now we assume that the value of is slightly undercritical, i.e., , where for any sign of . We are interested to expand the potential in the vicinity of the bending point . We, therefore, write () and substitute this into Eq. (32) and Taylor-expand up to . We get

(35)

Here the first terms in the l.h.s. and the r.h.s. cancel because of Eq. (34), the second term in the l.h.s. vanishes because of Eq. (33). As a result we obtain new static solutions for an undercritical shifted from by

(36)

to the positive or negative direction. One of them is stable and corresponds to the minimum of , another, unstable one, corresponds to the maximum of . The energy barrier

(37)

After expanding up to , we see that the terms and cancel, and we obtain

(38)

Using Eq. (34), the definition of and the expression (36), we finally obtain

(39)

Now let us apply this general result to our system.

Figure 4: The energy barrier prefactor given by Eq. (43) vs. (black) and the one calculated using the 0-th order approximation Eq. (40) (gray). (a) shows global behavior in the interval , while (b) shows the zoom of the area close to , where multiple solutions appear. Parameters are: , .

If or if is far away from , then we can use only the 0-th order term in Eq. (3) and in Eq. (6). In this limit the CPR is sinusoidal, see Eq. (16). Although it is shifted by (), it is irrelevant for calculation of the escape barrier and eigenfrequency. Thus, the system behaves as a conventional JJ with sinusoidal CPR, critical current and Josephson energy . Thus, , and we obtain the usual approximation for energy barrier and eigenfrequency in the limit :

(40)

In the case when the second order correction is important, i.e., and , we use the same approach, but again, like in the section about critical current and ground state phase, we approximate for (). In this case the energy is given by

(41)

From here

(42)

We have to take in Eq. (42) and substitute this into Eq. (39). The dependence of on is obvious from Eq. (39), so our aim is to see how prefactor in Eq. (39) depends on (or ). Since, for each can be found only numerically, we, as in the previous sections, sweep from to and find the corresponding from Eq. (24) and then calculate from Eq. (42). Then we make a parametric plot of the energy barrier prefactor

(43)

as a function of , see Fig. 4. The global behavior is given by the 0-th order approximation, see prefactor in Eq. (40). The second order approximation, where we expanded all expressions near , works well near , but deviates substantially from the real solution given by the 0-th order approximation when or . In Fig. 4(b), as increases, one sees two branches, given by Eq. (43). One of them corresponds to the negative critical current branch, cf. Fig. 3(b), another to the positive one. At slightly larger than the positive branch crosses zero, see Fig. 3(b), so that we see that the prefactor also vanishes at this point, see Fig. 4(b). Then, at somewhat larger both mentioned branches join, see Fig. 3(b). At this point the prefactor diverges, see Fig. 4(b). The other two branches in Fig. 4(b) show similar behavior.

iii.4 Eigenfrequency

Figure 5: The eigenfrequency prefactor given by Eq. (47) vs. given by Eq. (24) (black) and the one calculated using the 0-th order approximation Eq. (46) (gray). (a) shows global behavior in the interval , while (b) shows the zoom of the area close to , where multiple solutions appear. Parameters are: , .

In general, the eigenfrequency of phase oscillations around one of the static solution , see Eq. (36), is given by

(44)

The first term vanishes because of Eq. (33). So, using Eq. (36) we get

(45)

In the 0-th approximation , and we arrive to the well-known result

(46)

In the second order approximation we again sweep and make an implicit plot of the prefactor

(47)

as a function of given by Eq. (24). The behavior of the eigenfrequency prefactor is shown in Fig. 5. Similar to the energy barrier prefactor, the eigenfrequency prefactor given by Eq. (47) describes the multiple solutions near well. However, the 0-th approximation is better outside this vicinity.

iii.5 Escape histogram width in the MQT regime

Figure 6: The width of the escape current histogram in MQT regime calculated using the 0-th order approximation Eq. (50) (gray line) and the second order approximation linearized near , see Eq. (49). (a) shows the global behavior in the interval , while (b) shows the zoom of the area close to , where multiple solutions appear. Parameters are: , .

The dependences and allow to directly calculate not only the escape rate, but also the width of the escape histogram as a function of . This dependence can be directly compared with the experimentally measured one. For the sake of simplicity we limit ourselves to the case of MQT, so that the temperature is excluded. The approximate, but rather precise, formula for the histogram width (dispersion) was derived by GargGarg (1995) in the general case of a particle in a tilted potential. For MQT regime the GargGarg (1995) expression reduces to

(48)

where we have omitted -terms that are much weaker than power-terms. We use a sign as we are interested not in the width itself but in its scaling as a function of . This is a dispersion of defined above, i.e., it assumes that the critical current is equal to 1. If the critical current is equal to , the sigma (measured in the same units as ) is

(49)

When , or is not very close to , we can use the 0-th order approximation. In this case, by substituting the prefactors from Eq. (40) and Eq. (46) into Eq. (49), we get

(50)

This dependence is shown by the gray line in Fig. 6.

For small and in the vicinity of , we have to use the second order formulas. Again by substituting prefactors from Eq. (43) and Eq. (47) into Eq. (49) we obtain , which we plot vs. as a parametric plot, see Fig. 6. One can see that vanishes at the bifurcation point where the two branches join. Note also that at the points where vanishes (crosses zero), i.e. is linear, both and vanish, however does not have zero or any other peculiarity at these points, as can be seen from Eq. (49).

iii.6 Experimental relevance

On one hand the range of (), required to create a JJ near , is very tiny. This was already pointed out in the previous worksBuzdin and Koshelev (2003); Goldobin et al. (2011); Sickinger et al. (2012); Lipman et al. (2014); Goldobin et al. (2015). This makes it very difficult to controllably fabricate the desired — a small technological shift can drastically change the junction.

However even a nominally symmetric junction (), due to a tiny technological misalignment can get . As a result an experimental curve exhibits asymmetric minima for positive and negative bias current, as in Fig. 3. Also, in experiment it should be easy to conclude whether the asymmetry is smaller than (and we deal with a JJ) or larger (and we deal with single state JJ). In the former case the dependence should have a cusp-like minimum with branch crossing, while in the latter case the minimum will be smooth.

Iv Conclusions

We have derived an effective model, which describes a short JJ with a phase discontinuity at an arbitrary point along its length. This model reduces the system considered here to a point-like JJ with an unconventional CPR. One can relatively easy obtain all desired characteristics of such a point-like JJ. For example, we analyzed the ground state and found that close to one obtains a JJ, while far from this point it is JJ. We also calculated the dependence of the critical current of such a JJ as a function of and found multiple branches close to , corresponding to states. Further, we have calculated the behavior of the energy barrier and eigenfrequency close to the critical current, which allow to make estimations of the width of the switching current histogram in the regime of macroscopic quantum tunneling.

Appendix A Derivation of the averaged CPR

We use the following phase ansatz

(51)

which corresponds to the phase discontinuity at and a Taylor expansion of the phase in the left and right region (subscripts and ) relative to the discontinuity .

In statics, the phase should satisfy the Ferrel-Prange equation

(52)

subject to boundary conditions at the edges , corresponding to zero applied magnetic field, and the field continuity at :

(53a)
(53b)
(53c)

The prime denotes . Below, we use the junction half-length as a small parameter and develop a perturbation theory with respect to .

a.1 0-th approximation in

We substitute the ansatz Eq. (51) into the Ferrel-Prange Eq. (52). After calculating , we would like to expand with respect to the small parameter . The key point is to make this expansion correctly. For this, we transform the argument of the sine function to explicitly pull out from all terms. Namely, we define that and from ansatz Eq. (51) depend on as

(54)
(55)
(56)

Here and below the subscripts separated by a comma mean that it is actually two equations: one is obtained by taking the first subscript in the whole equation, the second equation is obtained by taking the second subscript in the whole equation. The higher order and -terms from Eq. (51) are not relevant in -th approximation. Initially, the scaling of and with is actually not obvious, but later we will see that the scaling given by Eqs. (55) and (56) is consistent. After the above substitution we expand in Eq. (52) relative to , keeping only constant terms (neglecting and smaller). We arrive at the following expression(s).

(57)

where in our -th approximation. From the Eqs. (57) it is obvious that scales as written in Eq. (56). From the boundary conditions Eq. (53) we have

(58a)
(58b)
(58c)

where . It is Eqs. (58a) and (58b) where it becomes obvious that , as it was correctly written in Eq. (55), otherwise the l.h.s. and the r.h.s. would have different orders in . By substituting from Eqs. (57) into Eqs. (58a) and (58b) and then and from Eqs. (58a) and (58b) into Eq. (58c), we finally get the current-phase relation

(59)

a.2 2-nd order approximation

In the next order () approximation we use all the terms in ansatz (51) to substitute into the Eq. (52). After calculating , we explicitly extract from all terms using the following substitutions

(60a)
(60b)
(60c)
(60d)

Here , and are from the 0-th order approximation and and are the next order corrections. Other powers, e.g., . After the above substitution we expand in Eq. (52) relative to , keeping only terms and larger (neglecting and smaller). We arrive at the second order polynomial in equal to zero. Obviously, it can be equal to zero for any only if each coefficient in front of , and constant are all equal to zero. We, thus, obtain

(61a)
(61b)
(61c)

As a next step we substitute the ansatz (51) into the boundary conditions (53). Then we substitute the definitions (60), cancel the terms from 0-th order approximation (if any), substitute , and from Eqs. (61) and obtain

(62a)
(62b)
(62c)

By substituting Eqs. (62a) and (62b) into Eq. (62c) we finally obtain

(63)

Finally, we substitute from Eq. (57) (0-th approximation) and obtain the final expression for the second order correction to the current

(64)

a.3 From to the average phase

Up to now both the 0-th order CPR Eq. (59) and the second order CPR Eq. (59) are given as functions of the , while our aim is to express those as a function of the average phase . To find , we substitute ansatz Eq. (51) into Eq. (1), integrate on each interval, substitute the definitions of , , , from Eqs. (60), then substitute expressions for , , , and from Eqs. (57), (61), (62) and (59). Then we keep only the terms and larger, after some simplifications arrive at

(65)

Our aim is to invert this expression, i.e., to express to substitute to CPR and obtain . We again act following the perturbation theory with respect to the small parameter . In the 0-th approximation

(66)

Here we introduced the angle , which makes expressions more compact.

In the next (second) approximation . By substituting this into Eq. (65) and expanding up to , we obtain and therefore