Radiation emission due to fluxon scattering on inhomogeneity in large two-dimensional Josephson junction
aa r X i v : . [ c ond - m a t . s up r- c on ] D ec Radiation emission due to fluxon scattering on an inhomogeneity in a largetwo-dimensional Josephson junction
Ivan O. Starodub ∗ and Yaroslav Zolotaryuk † Bogolyubov Institute for Theoretical Physics, National Academy ofSciences of Ukraine, vul. Metrologichna 14B, 03680 Kyiv, Ukraine (Dated: July 12, 2018)Interaction of a fluxon in the two-dimensional large Josephson junction with the finite-area in-homogeneity is studied within the sine-Gordon theory. The spectral density of the emitted planewaves is computed exactly for the rectangular and rhombic inhomogeneities. The total emittedenergy as a function of the fluxon velocity exhibits at least one local maximum. Connection to thepreviously studied limiting cases including the point impurity and the one-dimensional limit hasbeen performed. An important feature of the emitted energy as a function of the fluxon velocityis a clear maximum (or maxima). The dependence of these maxima on the geometric properties ofthe impurity has been studied in detail.
PACS numbers: 03.75.Lm,74.50.+r,74.62.En
I. INTRODUCTION
Studies of the fluxon (Josephson vortex) dynamics inlarge Josephson junctions (LJJs) is an important prob-lem in modern superconductivity. The LJJs can be spa-tially inhomogeneous either due to the production defectsor can be manufactured in such a way on purpose. Thus,the problem of the fluxon interaction with the spatial in-homogeneity (microshort,microresistor, Abrikosov vortexetc.) is of remarkable importance . As a result of thefluxon-impurity interaction the radiation of the small-amplitude linear waves (Josephson plasmons) occurs .The issue of the linear wave radiation due to the fluxoncollision with the spatial inhomogeneity has been studiedin detail for the one-dimensional case (1D). Most of these(both theoretical and experimental) studies have focusedon the scattering on the point-like inhomogeneity beingeither a microshort or a microresistor , or a mag-netic impurity . An extended inhomogeneity has beeninvestigated as well as the interface separating two dif-ferent junctions .An important thing to note is that a 1D Joseph-son junction is only a 1D approximation of the two-dimensional (2D) LJJ of the finite width. Thus, a naturalquestion is to take the transverse direction into accountand to study the fluxon scattering on an impurity in thissituation. Moreover, fluxon dynamics in the large areaJJ is an interesting and important problem in its ownright. It has been studied in different contexts such asdynamical properties , pinning on impurities andapplications . However, up to now the radiation emis-sion due to the 2D fluxon scattering on the impurity hasbeen studied in detail only for the special case of thepoint-like impurity described by the Dirac δ - function .Thus, the aim of this paper is to study the propertiesof the small-amplitude wave radiation that appears asa result of the fluxon transmission through the inhomo-geneity of the general shape.The paper is organized as follows. In the next sec-tion, the model is described. Section III is devoted to the studies of the radiation emitted due to the fluxon-impurity interaction. In the last section, the discussionand conclusions are presented. II. THE MODEL
We consider fluxon dynamics in the LJJ with spatialinhomogeneities. The main dynamical variable is thedifference between the phases θ ( x, y ; t ) − θ ( x, y ; t ) . = φ ( x, y ; t ) of the macroscopic wave functions of the the su-perconducting layers of the junction, also known as theJosephson phase. In the bulk of the junction this variablesatisfies the equation ∂ x H y − ∂ y H x = j c [1+ f I ( x, y )] sin φ + ~ C ( x, y )2 e ∂ t φ , (1)where the function f I ( x, y ) describes the critical currentchange on the spatial inhomogeneity and the magneticfield components H x,y are related to the Josephson phaseas H x = − ~ eµ l ( x, y ) ∂ y φ , H y = ~ eµ l ( x, y ) ∂ x φ . (2)The junction capacitance C ( x, y ) is spatially inhomoge-neous due to the impurity. Among other parameters j c isthe critical current density away from the impurity, e isthe electron charge, µ is the vacuum permeability and ~ is Planck’s constant. The value l ( x, y ) describes thethickness of the layer that allows magnetic field penetra-tion. It varies in space due to the presence of the impu-rity and can be written as l ( x, y ) = 2 λ L + d i ( x, y ), where λ L is the superconductor London penetration depth and d i ( x, y ) is the insulating layer thickness. Away from theimpurity d i ( x, y ) = d = const while d i ( x, y ) = d + d = const inside the impurity. For the impurity of the gen-eral shape that covers a certain segment Ω ∈ R of thejunction one can write f I ( x, y ) = (cid:26) µ I if ( x, y ) ∈ Ω , x, y ) / ∈ Ω . (3)Similarly, the spatial change of the magnetic length andcapacitance is given by l ( x, y ) = (cid:26) d + 2 λ L + d if ( x, y ) ∈ Ω ,l = d + 2 λ L if ( x, y ) / ∈ Ω . (4)and C ( x, y ) = C d d i ( x, y ) = C [1 + f C ( x, y )] ,f C ( x, y ) = (cid:26) µ C = − d d + d if ( x, y ) ∈ Ω , x, y ) / ∈ Ω . (5)where C is the junction capacitance per unit area awayfrom the impurity. For the sake of convenience the fol-lowing function can be introduced l l ( x, y ) = 1 + f H ( x, y ) = 1 + (cid:26) x, y ) / ∈ Ω µ H if ( x, y ) ∈ Ω ,µ H = l l + d − − d d + d + 2 λ L . (6)Equation (1) can be rewritten in the dimensionless formby normalizing the spatial variables x and y to theJosephson penetration depth λ J and the time t to theinverse Josephson plasma frequency ω − J . As a result,the two-dimensional perturbed sine-Gordon (SG) equa-tion is obtained: {− ∂ x [1 + f H ( x, y )] ∂ x − ∂ y [1 + f H ( x, y )] ∂ y + (7)+[1 + f C ( x, y )] ∂ t (cid:9) φ + [1 + f I ( x, y )] sin φ = 0 . For details one might consult the textbooks . The im-purity is a microshort if µ I > d < µ I < d >
0. Hence µ H /µ I > µ C /µ I > the insu-lating layer thickness d ∼ A , while the London pene-tration depth λ L is of the order of several tens of ˚ A , theinequality | µ H | < | µ C | holds. III. RADIATION EMISSION
Fluxon interaction with the spatial inhomogeneity isnormally accompanied with the radiation of the small-amplitude electromagnetic waves (Josephson plasmons).Below we present the general scheme for the calcula-tion of the radiation created by the fluxon-impurityinteraction which is based on the method developedfor the delta-like impurity or for the respective 1Dproblems . Only the main points of the derivationprocedure are presented. For the details the interestedreader can consult the above-mentioned papers. A. General framework
Both sides of the SG equation (7) can be divided by[1 + f C ( x, y )], and- as a result it can be rewritten as ∂ t φ − ∆ φ + [1 + ¯ f I ( x, y )] sin φ = ¯ f H ( x, y )∆ φ ++ 11 + f C ( x, y ) [ ∂ x f H ( x, y ) ∂ x φ + ∂ y f H ( x, y ) ∂ y φ ] , (8)where ∆ = ∂ x + ∂ y and¯ f I ( x, y ) = 1 + f I ( x, y )1 + f C ( x, y ) − (cid:26) x, y ) / ∈ Ω¯ µ I if ( x, y ) ∈ Ω , ¯ f H ( x, y ) = (cid:26) x, y ) / ∈ Ω¯ µ H if ( x, y ) ∈ Ω , ¯ µ I,H = µ I,H − µ C µ C . (9)We seek the solution of the SG equation (7) as a su-perposition of the exact soliton solution and the plas-mon radiation on its background: φ ( x, y, t ) = φ ( x, t ) + ψ ( x, y, t ). The spatial inhomogeneity is considered asa small ( | µ I,H,C | ≪
1) perturbation. Here φ ( x, t ) =4 arctan h exp (cid:16) x − vt √ − v (cid:17)i is the exact soliton solution ofthe unperturbed 1D SG equation and ψ ( x, y, t ) is the ra-diative correction, | ψ | ≪ φ . It is convenient to work inthe reference frame that moves with the fluxon velocity v : ξ = x − vt √ − v , τ = t − vx √ − v . In these new variables wehave φ ( x, t ) = φ ( ξ ) = 4 arctan (exp ξ ).In the moving reference frame the equation that de-scribes the emitted radiation reads (cid:8) ∂ τ − ( ∂ ξ + ∂ y ) + cos[ φ ( ξ )] (cid:9) ψ = R ( ξ, y ; τ ) , (10)where the right-hand side of Eq. (10) is completely de-fined by the impurity: R ( ξ, y ; τ ) = 2 (cid:20)(cid:18) − − v ¯ µ H ¯ µ I (cid:19) tanh ξ cosh ξ × (11) × ¯ f I (cid:18) ξ + vτ √ − v , y (cid:19) + h H (cid:16) ξ + vτ √ − v , y (cid:17) √ − v cosh ξ ,h H ( x, y ) = ∂ x f H ( x, y ) . In this expression it has been taken into account thatsin[ φ ( ξ )] = ∂ ξ φ ( ξ ) = − ξ/ cosh ξ . Also, forany two functions of the type f α ( x, y ) or ¯ f α ( x, y ) ( α = I, C, H ) the following equality is true: f α ( x, y ) = µ α f β ( x, y ) /µ β . Here the last term of R ( ξ, y ; τ ) that con-tains the function h H ( x, y ) is associated with the fluxoninteraction with the borders of the impurity because h H ( x, y ) = 0 only there, i.e., if ( x, y ) / ∈ ∂ Ω. The firstterm corresponds to the radiation produced when thefluxon passes the bulk of the impurity.The solution of Eq. (10) can be represented as ψ ( ξ, y, τ ) = Z + ∞−∞ Z + ∞−∞ a ( q ξ , q y ; τ ) ϕ ( ξ, y ; q ξ , q y ) dq ξ dq y , (12)where ϕ ( ξ, y ; q ξ , q y ) is the eigenfunction of the homo-geneous part of this equation: ϕ ( ξ, y ; q ξ , q y ) = e i ( q ξ ξ + q y y ) (2 π ) / q ξ + i tanh ξ (1 + q ξ ) / , (13) Z + ∞−∞ Z + ∞−∞ ϕ ∗ ( ξ, y ; q ξ , q y ) ϕ ( ξ, y ; q ′ ξ , q ′ y ) dξdy == 12 π δ ( q ξ − q ′ ξ ) δ ( q y − q ′ y ) . (14)Here δ is the Dirac delta function, q ξ and q y are thecomponents of the plasmon wave vector in the movingframe, and ¯ ω = q q ξ + q y , (15)is the plasmon dispersion law in that frame. The function a ( q ξ , q y ) is the radiation amplitude. It is convenient tointroduce another function which also describes the emit-ted radiation, namely b ( q ξ , q y ; τ ) . = ( a τ − i ¯ ωa ) exp( i ¯ ωτ ).As a result, the following equality holds: ∂ τ b = e i ¯ ωτ (cid:0) ∂ τ a + ¯ ω a (cid:1) . (16)Multiplying both sides of Eq. (10) by ϕ ∗ ( ξ, y ; q ′ ξ , q ′ y ) andintegrating simultaneously over y ∈ R and ξ ∈ R weobtain δ ( q ξ − q ′ ξ ) and δ ( q y − q ′ y ) on the left-hand side [theorthogonality condition (14) is used] of Eq. (10). Afterremoving the integration over q ξ and q y one arrives atthe following expression ∂ τ b = 2 π e i ¯ ωτ Z + ∞−∞ Z + ∞−∞ R ( ξ, y ; τ ) ×× ϕ ∗ ( ξ, y ; q ξ , q y ) dξ dy . (17)The total radiation over the whole time is defined by thefunction B ( q ξ , q y ) = Z + ∞−∞ ∂ τ b ( q ξ , q y ; τ ) dτ . (18)Thus, with the pair of equations (17) and (18) one has thecomplete formula for the energy calculation. From thispoint it is possible to proceed to the emitted radiationstudies for the particular shapes of Ω. The return tothe laboratory frame is performed with the help of thefollowing Lorentz transformation: q x = q ξ + v ¯ ω √ − v , ω = vq ξ + ¯ ω √ − v , (19) q ξ = q x − vω √ − v , ¯ ω = ω − vq x √ − v . (20)The q y component remains unchanged. Taking intoaccount that the emitted energy density equals E ( q x , q y ) ≃ | B ( q x , q y ) | / (4 π ), the total energy is givenby the integral E = Z + ∞−∞ Z + ∞−∞ E ( q x , q y ) dq x dq y . (21) The following simplification can be achieved if Ω has theproperties defined below. Suppose the impurity coversthe area that is limited by the lines x = x and x = x along the y axis and by the continuous and single-valuedfunctions y = g ± ( x ) along the x axis, as is shown in Fig.1. In this case f I,H,C ( x, y ) = µ I,H,C [ θ ( x − x ) − θ ( x − x )] ×× { θ [ y − g − ( x )] − θ [ y − g + ( x )] } , (22)and the integral over y can always be taken. As a result,the computation of the radiation function b is reducedconsiderably. Here θ ( x ) is the Heaviside function. g (x) − x x g (x) + Ω fluxon propagationy x FIG. 1: Schematic top view of the impurity area Ω.
Below we consider the concrete examples when the im-purity area Ω is limited by the piecewise functions.
B. Rectangular impurity
In this subsection the rectangular impurity of finite sizein both x and y directions, f I,H,C ( x, y ) = µ I,H,C (cid:20) θ (cid:18) x + d x (cid:19) − θ (cid:18) x − d x (cid:19)(cid:21) ×× (cid:20) θ (cid:18) y + d y (cid:19) − θ (cid:18) y − d y (cid:19)(cid:21) , (23) h H ( x, y ) = µ H (cid:20) δ (cid:18) x + d x (cid:19) − δ (cid:18) x − d x (cid:19)(cid:21) ×× (cid:20) θ (cid:18) y + d y (cid:19) − θ (cid:18) y − d y (cid:19)(cid:21) , (24)is considered. The parameters d x and d y are the impuritylength and width, respectively.
1. Spectral density of the emitted waves
At this point we can substitute the actual expressions(23) and (24) that corresponds to the rectangular impu-rity into Eqs. (17) and (18). Then the radiation function(18) in the moving frame is obtained after the consecutive integration over the y , τ , and ξ variables: B ( q ξ , q y ) = i √ πµ I q y q q ξ (1 − v ) / v sin (cid:18) q y d y (cid:19) sin ¯ ω √ − v v d x ! sech h π v ( q ξ v + ¯ ω ) i ×× ((cid:18) − v − µ H µ I + v µ C µ I (cid:19) [¯ ω − (1 + q ξ ) v ]1 + µ C + 2 µ H µ I (1 − v )¯ ω ) . (25)The first term in the curly brackets in Eq. (25) appearsdue to the first term in R [see Eq. (11)] and can betreated as a result of the fluxon interaction with the bulkof the impurity. The second term in the curly brack-ets appears due to the second term (associated with the function h H ) in Eq. (11) and can be considered as theradiation that appears due to the fluxon interaction withthe border of the impurity. After returning to the labo-ratory frame of reference with the help of Eqs. (19)-(20)the final formula for the spectral density reads: E ( q x , q y ) = 2 µ I v (cid:20) sin ( q y d y / q y (cid:21) (cid:26) sin [ d x ( ω − vq x ) / v ] ω − q x v (cid:27) sech (cid:16) πω v p − v (cid:17) ×× (cid:26) − v − µHµI + v µCµI µ C (cid:2) ( ω − q x v ) + q y v (cid:3) + 2 µ H µ I ( ω − vq x ) (cid:27) ( ω − q x v ) + ( v − q y , (26) ω = q q x + q y . (27)This function is symmetric with respect to the mirrorsymmetry q y → − q y and to the transform q x → − q x , v → − v . Therefore, it is sufficient to restrict the plots of E ( v ) to the interval 0 ≤ v ≤
1. In order to compute thetotal emitted energy E ( v ) [see Eq. (21)] it is necessary touse numerical methods because it is not possible to takethe respective double integral explicitly.
2. 1D limit
Before embarking on the investigation of the full 2Dproblem it is instructive to recall the corresponding one-dimensional (1D) case of the fluxon scattering on theimpurity with the length d x . Formally this limit can beachieved if d y → ∞ . The energy density in this case isalready known from the previous work : E ( q ) = πv µ I − v − µ H µ I + v µ C µ I µ C + 2 µ H µ I ! ×× sin (cid:20) d x v (cid:16)p q − qv (cid:17)(cid:21) × × sech π √ − v v p q ! . (28)However, in the paper, cited above, the spatial inhomo-geneity of the capacitance was not taken into account.We note that Eq. (28) can be obtained in the limit q y → µ I should be renormalized as µ I d y → µ I ).This means that the impurity width d y tends to infinity,and, as a result, the scattering does not create any radi-ation in the y direction, leaving the problem completelyinvariant in that direction, i.e., one-dimensional.Typical dependencies of the spectral density E = E ( q )for the different values of the fluxon velocity are given inFig. 2. It is easy to see that the energy density E ( q ) [Eq.(28)] has an infinite countable set of global minima forwhich E ( q min ) = 0. They are the roots of the equation d x ( p q min − q min v )2 v = πn, n = n , n + 1 , . . . ,n = ⌈ d x (1 − v ) / / (2 vπ ) ⌉ > , (29)where ⌈ x ⌉ is the ceiling function of x . Similarly, thereare maxima that are placed between those minima at thevalues of q that are the roots of the equations d x ( p q max − q max v )2 v ≈ π (2 n − , (30) n = n , n + 1 , . . . . The minima and maxima are associated with the con-structive and destructive interference of the plasmons,emitted when the fluxon enters and exits the impurity.Depending on the length of the impurity and the fluxonvelocity, the radiated plasmons can either cancel eachother if their phases differ by ± π or can enhance eachother if their phases coincide. The radiation consist ofthe forward ( q >
0) and backward ( q <
0) emitted plas-mons, and the energy of these plasmons is distributednon-homogeneously with respect to q . First of all, themost of the energy is concentrated in the long-wavelengthmodes due to the presence of the sech ( · · · ) term in Eq.(28). Secondly, as can be seen from Fig. 2, the distribu-tion of the backward radiation is defined by the extrema(29) and (30) that lie on the negative half-axis ( q < πv/ [ d x (1 + v )]; therefore, the smallchange of v will lead to the small change in the area underthe E ( q ) curve. On the contrary, the forward radiationdepends strongly on v , especially if v is not small ( v < v ≪ q ’s the extrema are dis-tributed with the almost fixed step 2 πv/ [ d x (1 − v )]. Theminima of E ( q ) given by Eq. (29) come in pairs, num-bered by the index n . These pairs are placed on the dif-ferent sides from the value q = v/ √ − v , which is theminimum of the left-hand side of Eqs. (29) and (30). Thepair with n = n is the pair of the minima, that are theclosest to each other. There always should be a maximumbetween these minima. If the above-mentioned minimaare very close to each other (2 πn v/d x > ∼ √ − v ), themaximum between them cannot be associated with Eq.(30), as seen in Figs. 2(a) and 2(c); thus, the respectivevalue of E lies not on the sech ( · · · ) envelope function,but significantly below it. As a result, for these valuesof v the forward radiation can be insignificant, as can beobserved from the area below the curve E ( q ) at q >
0. Inanother case, the pair of minima that correspond to n are significantly separated, and the maximum betweenthem belongs to the set (30). It is again the first max-imum at the positive axis, and it attains the value of E which is quite large comparing to the previous case, ascan be seen in Figs. 2(b) and 2(d).The dashed lines 6 and 7 in Fig. 3 show the depen-dence of the total emitted energy on the fluxon velocity(the solid lines correspond to the 2D case which will bediscussed later). The values of v which correspond tothe minima of the E ( v ) in line 6 in Fig. 3, have theminimal forward emission, and the respective spectralenergy distributions are shown in Fig. 2(a), 2(c). Thevalues of v that are the maxima of E ( v ) correspond tothe maximal forward emission and the respective spectraldistributions are given in Figs. 2(b) and 2(d). Thus, themaxima of the total energy coincide with the maximal FIG. 2: (Color online). Energy density [see Eq. (28)] for the1D junction with d x = 8, µ H = µ C = 0 at the fluxon velocity v = 0 .
398 (a), v = 0 .
488 (b), v = 0 .
552 (c) and v = 0 .
676 (d).The red dashed line depicts the sech “envelope” term in Eq.(28). forward emission while the minima of E ( v ) correspondto the minimal forward emission. It should be notedthat the minima [Eq. (29)] and maxima [Eq. (30)] ofthe energy density are distributed approximately equidis-tantly for the short-wavelength ( | q | ≫
1) modes but withthe different step for q > q <
0. In the limit | v | ≪ πv/d x . Consequently, in the limit | v | → E ( v )dependence, and this can be noticed from the inset.Finally, we remark that in the relativistic limit v → E ( v ) → µ H = µ C = 0 and E ( v ) →∞ if µ H,C = 0. The details of this limit will be discussedbelow together with the 2D case.
3. Total emitted energy in the 2D case
First of all, we discuss the dependence of the totalemitted energy E ( v ) on the impurity parameters µ I , µ H and µ C . We remind the reader that µ I is associatedwith the change of the critical current, while µ C [see Eq.(5)] and µ H [see Eq. (6)] appear due to the narrowing ordistension of the insulating area. If µ H = µ C = 0 the im-purity corresponds only to the local change of the criticalcurrent without any changes in the insulating layer thick-ness. The total emitted energy for the different values of µ H and µ C is given in Fig. 3. We note the principal dif-ference in the behavior of the E ( v ) function in the limit v → µ C,H = 0 compared to the case µ C = µ H = 0.In the latter case E ( v ) tends to zero while in the formercase it diverges: E ( v ) v → → + ∞ . The same is observedin the 1D case (shown by the dashed lines). It is quite FIG. 3: (Color online). Total emitted energy (normalizedto µ I ) as a function of the soliton velocity for the impuritywith d x = 8, d y = 8 and µ H = 0, µ C = 0 . µ C /µ I = 1(curve 1), µ H = 0, µ C = − . µ C /µ I = 0 . µ H /µ I = 0 . µ C = 0 . µ C /µ I = 0 . µ H =0, µ C = 0 . µ C /µ I = 0 . µ H = µ C = 0(curve 5). The dashed lines 6 and 7 correspond to the samedependence but for the 1D problem [see Eq. (28) for thespectral energy density] with d x = 8 and µ H = µ C = 0 (curve6) and µ H = 0, µ C /µ I = 0 . µ C = 0 .
05 (curve 7). Thesedependencies are multiplied by a factor 10 for the sake ofconvenience. The inset shows the details of the curves 4 and7. obvious from the lines 1 and 4 that for the larger valuesof µ C the value of the emitted energy is larger. If onetakes two opposite values of µ C , the case of a microresis-tor ( µ C <
0, line 2) yields slightly larger energy emissioncompared to the case of a microshort ( µ C >
0, line 4)due to the presence of the (1 + µ C ) − coefficient in theenergy density (26). The effect of the spatial variationof the magnetic field, governed by the coefficient µ H isnegligible, as one can observe from the comparison of thelines 3 and 4. Therefore, we will assume µ H = 0 furtheron throughout the paper.The divergence at v → − v ) − and (1 − v ) − / . In the former term the function ¯ f I contains both the parameters µ C and µ H and is alwaysfinite, thus the divergence appears only due to the divi-sor. In the latter term, in addition, there is a function h H which is non-zero only on the edges of the inhomogene-ity, where it is proportional to the Dirac’s δ − function.This term generates the sharp growth of radiation whenthe fluxon interacts with the edges of the impurity. Inthe 1D case it produces such growth only at the entrance( x = − d x /
2) and exit ( x = d x /
2) points of the impurity.We would like to mention that the divergence at v → E ( v ) dependence such as the mul-tiple extrema will be discussed below. At this point weonly note that as µ C decreases, the positions of the ex-trema do not shift significantly, but the absolute valuesof E at the extrema decrease. This happens because thecontribution to the emitted radiation due to the narrow-ing/expansion of the insulating layer, decreases. Depend-ing on the value of µ C some extrema can disappear dueto the growth of E ( v ) as v → v → d x = 8 while its width d y is varied. The 1D result for the same length is plottedwith the dashed line as a reference. Naturally, the valueof the emitted energy decreases as d y decreases. Moreinterestingly, the extrema become less pronounced, and,finally no extrema are seen in curve 4 that corresponds tothe case d y = 2. In the case µ H = µ C = 0 we obtain thesame picture: compare curve 5 of Fig. 3 ( d y = 8), curve5 of Fig. 4 ( d y = 6), and curve 6 of Fig. 4 ( d y = 2). Themaxima become more shallow and gradually disappear.The following interpretation of the obtained results canbe made. The shape of the energy density distributionis given in Fig. 5. The absolute minima of the energydensity satisfy E ( q x , q y ) = 0 and these minimal valuesare attained at the following set on the ( q x , q y ) plane: q y = 2 πnd y , n = ± , ± , . . . for any q x , (31)(1 − v ) q x + q y = (cid:18) πmvd x (cid:19) − πmv d x q x , (32) m = n , n + 1 , . . . , where n is given by Eq. (29). Thus, the minima arelocated on the set of parallel lines (31) as well as on theset of embedded ellipses given by Eq. (32). The ridgesof the maximal E lie between the curves, defined by theroots of Eq. (31). For large d y these ridges are stronglylocalized in the q y direction [see Figs. 5(a) and 5(b)],while the decreasing of d y makes them concentric andcrescent-like as shown in Figs. 5(c) and 5(d).For the large values of d y the problem can be treated asan almost 1D, so that most of the emitted radiation trav-els in the x direction while the y - component of the radi-ation remains insignificant. This can be clearly observedin Figs. 5(a) and 5(b) where the spectral density E ( q x , q y )(26) is plotted for the values of velocity close to the min-imum [panel a] and maximum [panel b] of the curve 1 inFig. 4. Since the decay of the function [sin( q y d y / /q y ] with the growth of q y is quite fast for the large valuesof d y , the energy density function remains strongly local-ized along the q x axis in the neighbourhood of the q y = 0 FIG. 4: (Color online) Total emitted energy (normalized to µ I ) as a function of the fluxon velocity for µ H = 0, µ C /µ I =0 . µ C = 0 . d x = 8 and d y = 8 (curve 1), d y = 6 (curve 2), d y = 4 (curve 3), and d y = 2 (curve 4). The case µ H = µ C = 0is represented by the red curves 5 ( d y = 6) and 6 ( d y = 2).The dashed line corresponds to the respective 1D problemwith d x = 8 (for the sake of convenience it is multiplied by afactor of 10). line. Its behaviour along the q x axis is reminiscent of therespective 1D problem, see Eq. (28) and Fig. 2. Indeed,the minimum of the total emitted energy corresponds tothe minimal forward emission. It can be easily observedin Fig. 5(a) that the global maximum is placed on the q x axis at q x < q x > q x axis, and this happens at v = 0 .
73 which is quite closeto the maximum of the E ( v ) function (curve 1) in Fig.4. The further decreasing of d y smears out maxima inthe E ( v ) dependence (compare the curves 1-4 in Fig. 5)up to the point when only one local maximum can bespotted. The scattering problem cannot be consideredas a quasi-1D any more. The radiation distribution be-comes rather different as shown in Figs. 5(c) and 5(d).The maxima of E ( q x , q y ) still lie on the q x axis, but thecurves (32) that define the minimal values become dis-tinctly arc shaped. The y - component of the radiationbecomes more delocalized and the analogy with the 1Dpicture breaks down.It is interesting to note how the shape of the energydensity function varies in the extreme limits of the ve-locity value: v → v →
1. In the small velocitylimit | v | ≪ E are almost circles and the density func-tion is close to being radially symmetric, see Fig. 5(e).The increasing of v makes the ellipses Eq. (32) moreelongated in the x direction, as has been demonstratedpreviously [see Figs. 5(a)-(d)]. An interesting situationemerges in the opposite limit, namely if 1 − | v | ≪
1. Theglobal maximum that was positioned on the q x axis splitsup into two maxima that are now located off the q x axissymmetrically with respect to each other, as shown in FIG. 5: (Color online). Emitted energy density E for the datain Fig. 4, curve 1 at v = 0 .
55 (a) and v = 0 .
73 (b); curve 6at v = 0 . v = 0 .
75 (d), v = 0 . v = 0 .
99 (f).
Fig. 5(f). Physically this means the following. The slowfluxon “feels” the impurity as a wall and the emergentradiation moves mostly along the fluxon propagation di-rection. The fast (relativistic) fluxon interacts with theimpurity in such a way that the impurity acts like a groin(a wave- breaker) and the emitted radiation is split by theimpurity into two halves that have both x and y compo-nents. A the same time the x component of the radiationbecomes insignificant.The number of local extrema of E ( v ) depends on theimpurity length d x . This is easily demonstrated by Fig.6, where the number of the maxima decreases with thedecreasing of d x . This result is similar to the same situ-ation in the 1D model .
4. Limiting cases
It is of interest to check the limiting cases when in oneof the directions ( x or y ) the impurity becomes infinitelynarrow. In the first case the limit d x →
0, while µ ∗ = µ I d x remains constant, corresponds to the situation whenthe impurity becomes infinitely thin in the x direction.In Eqs. (23) and (24) the difference of the θ -functionsthat form the first factor in the functions f I,H,C ( x, y )(23) becomes the Dirac δ function while h H ≡
0. This
FIG. 6: Total emitted energy (normalized to µ I ) as a functionof the fluxon velocity for µ H = 0, µ C /µ I = 0 . µ C = 0 . d y = 8 and d x = 2 (curve 1), d x = 4 (curve 2), d x = 6 (curve3), and d x = 10 (curve 4). case with µ C = µ H = 0 has been studied previously .Yet another interesting limit can be considered if d y → µ ∗ = µ I d y . In other words, the impurity remainselongated in the x direction, but becomes infinitely thinin the y direction. The spectral density of the emittedplasmons in the cases mentioned above reads E ( q x , q y ) µ ∗ → sin (cid:16) q y d y (cid:17) v q y (cid:26) − v − µHµI + v µCµI µ C (cid:2) ( ω − q x v ) + q y v (cid:3) + 2 µ H µ I ( ω − vq x ) (cid:27) ( ω − q x v ) + ( v − q y ×× sech (cid:16) πω √ − v v (cid:17) , if d x → , v (cid:20) sin[ d x ( ω − q x v ) / v ] ω − q x v (cid:21) (cid:26) − v − µHµI + v µCµI µ C (cid:2) ( ω − q x v ) + q y v (cid:3) + 2 µ H µ I ( ω − vq x ) (cid:27) ( ω − q x v ) + ( v − q y ×× sech (cid:16) πω √ − v v (cid:17) , if d y → . (33)In these limits the modulation in q - space, caused by theinterference, disappears along the q x direction in the firstformula because the impurity length becomes infinitelysmall. For the same reason there is no interference alongthe q y component when d y → E ( v ) dependence disappearleaving only one local maximum. The limit of the point[ f I ( x, y ) = µ I δ ( x ) δ ( y )] impurity can be achieved eas-ily from the stripe impurity by taking in Eqs. (33) thelimits µ C,H → d x → d y → µ C = µ H = 0. The obtained formula coincides with theprevious result . C. Rhombic impurity
Now we consider the rhombus(diamond)-shaped impu-rity with d x and d y being its length and width respec- tively:Ω : | x | ≤ d x / \ | y | ≤ g ( x ) = d y (1 / − | x | /d x ) . (34)The tip of the rhombus is perpendicular to the fluxonline. Then h H ( x, y ) = µ H { [ δ ( x + d x / − δ ( x − d x / ×× [ θ ( y + g ( x )) − θ ( y − g ( x ))] ++ [ θ ( x + d x / − θ ( x − d x / ×× { δ [ y − g ( x )] + δ [ y + g ( x )] } g ′ ( x ) } , (35) g ′ ( x ) = − d y d x sign( x ) = − d y d x sign( ξ + vτ ) . (36)Substituting the formulas (35) and (36) into Eqs. (17)and (18) we obtain the radiation function in the movingframe: B ( q ξ , q y ) = i √ πµ I q y q q ξ (1 − v ) / v d x d y ((cid:18) − v − µ H µ I + v µ C µ I (cid:19) [¯ ω − (1 + q ξ ) v ]1 + µ C + 2 µ H µ I (1 − v )¯ ω ) ×× cos (cid:16) q y d y (cid:17) − cos (cid:16) ¯ ω √ − v v d x (cid:17) ( d x d y ¯ ω ) − v q y v − h π v ( q ξ v + ¯ ω ) i , (37)where the dispersion law ¯ ω = ¯ ω ( q ξ , q y ) in the movingframe is given by Eq. (15). The transition to the labo-ratory frame is performed in the standard way, and- as a result, the spectral energy density in the laboratoryframe is expressed by the following formula: E ( q x , q y ) = 2 µ I v (cid:18) d x d y (cid:19) cos( q y d y / − cos[ d x ( ω − q x v ) / v ] (cid:16) d x d y (cid:17) ( ω − q x v ) − v q y sech (cid:16) πω v p − v (cid:17) ×× (cid:20) − v − µHµI + v µCµI (1+ µ C ) (cid:0) ( ω − q x v ) + q y v (cid:1) + 2 µ H µ I ( ω − vq x ) (cid:21) ( ω − q x v ) + ( v − q y , (38)where the dispersion law ω = ω ( q x , q y ) in the labora-tory frame is given by Eq. (27). It may seem thatthis dependence has a singularity where the equation (cid:16) d x d y (cid:17) ( ω − q x v ) = v q y is satisfied. However, withthe help of the trigonometric formula cos a − cos b =2 sin[( a + b ) /
2] sin[( b − a ) /
2] it is straightforward to showthat the respective divergences cancel out.The total emitted energy as a function of the fluxonvelocity v is shown in Figs. 7-9. The first figure (Fig.7) focuses on the situation when the ratio d y /d x is fixedwhile the area covered by the impurity is varied. Themain figure correspond to the impurity with its narrowedge pointing towards the fluxon direction ( d x /d y = 4).The inset (a) describes the opposite situation: d y /d x = 4.In general, the dependence E ( v ) grows with v in the limit v ≪ v → µ C and µ H terms [otherwise, if µ C = µ H = 0, we have E ( v ) v → → µ C = µ H = 0 there is one well-established maximum of the E ( v ) dependence which is positioned very close to the value v = 1. As the impu-rity area is increased, the peak of the energy dependencesharpens, while the position of the maximum shifts to-wards the point v = 1. If µ C = 0 the main maximumdisappears due to the unbounded growth of the energydependence. There are other maxima of the E ( v ) depen-dence, however they are very weak and can be noticedonly if the respective region is zoomed [see the inset (b)].When the impurity area decreases, some of these maximadisappear [compare the curves 3 and 2 in the the inset(b)].Inset (a) of Fig. 7 corresponds to the situation whenthe impurity is elongated in the y - direction with the ra-tio d x /d y = 1 / d x d y / but the limit (33) is not restoredmathematically.It is possible to consider the limiting cases of the in-finitely narrow stripes: d x → d y →
0. If the impu-rity amplitude is redefined as µ ∗ = µ I d x (or µ ∗ = µ I d y ),the spectral density in these limits reads0 FIG. 7: Total emitted energy (normalized to µ I ) as a functionof the fluxon velocity for the rhombic impurity with the fixedratio d x /d y = 4. The solid lines correspond to µ C = µ H = 0, d x = 4 (curve 1), d x = 8 (curve 2) and d x = 12 (curve 3).The dashed line corresponds to d x = 12 and µ C = 0 .
01 and µ C /µ I = 0 . µ C = 0 . µ H = 0. The inset (a) correspondsto the case d x /d y = 1 / µ C = µ H = 0 d x = 1 (curve 1), d x = 2 (curve 2) and d x = 3 (curve 3) and µ C /µ I = 0 . µ H = 0, d x = 3 (dashed curve). The inset (b) shows thedetails of the main figure. E ( q x , q y ) µ ∗ → (cid:16) q y d y (cid:17) v d y q y (cid:26) − v − µHµI + v µCµI µ C (cid:2) ( ω − q x v ) + q y v (cid:3) + 2 µ H µ I ( ω − vq x ) (cid:27) ( ω − q x v ) + ( v − q y ×× sech (cid:16) πω √ − v v (cid:17) , if d x → , d x v (cid:26) sin[ d x ( ω − q x v ) / v ] ω − q x v (cid:27) (cid:26) − v − µHµI + v µCµI µ C (cid:2) ( ω − q x v ) + q y v (cid:3) + 2 µ H µ I ( ω − vq x ) (cid:27) ( ω − q x v ) + ( v − q y ×× sech (cid:16) πω √ − v v (cid:17) if d y → . (39)These limiting values of E are very similar to the analo-gous limits for the rectangular impurity (33). The onlyprincipal difference is the interference terms that are re-sponsible for the oscillations in the q x or q y directioncome with the power 4 and not 2 as in Eq. (33).Next we focus on the situation when the impuritywidth d y is fixed and its length d x is varied. In Fig.8(a) the dependence of the local maximum value (de-fined within the interval 0 ≤ v <
1) of the emitted en-ergy as a function of the rhombus angle arctan( d y /d x )is plotted. If µ C = µ H = 0 the max v ∈ [0 , E ( v ) de-pendence on the rhombus angle is a decaying functionalmost everywhere in the interval [0 , π/ d x → ∞ the maximum of E ( v ) grows as the amount ofthe emitted energy increases. Only in the neighborhoodof the angle π/ µ C = 0 such a dependence cannot be defined for the whole interval [0 , π/
2] and it starts from somecritical value of the angle (see the dependencies, markedby squares and inverted triangles) and continues till thevalue π/
2. Below this critical angle there is no local max-imum of E ( v ) because it becomes strictly monotonic. If d y is decreased, the dependence becomes a strictly decay-ing function (shown by the circles in Fig. 8) that coverthe whole interval [0 , π/
2] even if µ C = 0. In Fig. 8(b)the E ( v ) dependence is demonstrated in the limit of theextremely narrow rhombic impurity. If µ C > E ( v ) function is a monotonicallyincreasing function. If µ C = 0 there is a sharp maxi-mum very close to v = 1 and everywhere else the func-tion behaves almost identically to the case µ C >
0. Onecan notice a fine structure of multiple inflection points.These points are the remnants of the local maxima thatare clearly seen in the inset (b) of Fig. 7. The number1
FIG. 8: (a) The value of the local maximum max v ∈ [0 , E ( v )of the emitted energy (normalized to µ I ) as a function of theangle arctan( d y /d x ) for the parameters d y = 5, µ C = 0 . µ C /µ I = 0 . ∇ ), µ C = 0 . µ C /µ I = 0 .
05 ( (cid:3) ), µ C = 0 ( ◦ )and d y = 1, µ C = 0 ( ⋄ ). The solid line is used as a guide foran eye.(b) Emitted energy dependence (normalized to µ I ) as a func-tion of the fluxon velocity for d x = 100, d y = 5, µ C = 0 . µ C /µ I = 0 . µ C = 0 (curve 2). µ H = 0 every-where. of these inflection points increases as the length of therhombus d x increases. Here we observe a weak link withthe case of the rectangular impurity, studied in III B. Inthat case we reported the increasing of the number ofmaxima of E ( v ) when d x increased. For the rhombuswe see the maxima degenerate into the inflection points.The limit d x → ∞ means that the impurity acts as an ex-tremely narrow groin that does not cause much radiationdue to its narrowness for small and intermediate veloci-ties. Significant growth of the emitted radiation can bespotted only in the relativistic regime (1 − v ≪ d y → ∞ .In the limits d x → x direction. When this limit is approached the local maxi-mum of the E ( v ) dependence becomes less and less pro-nounced. The energy density is proportional to d x ; thus, it is not surprising that the total energy tends to zero inthis limit. The renormalization of the impurity ampli-tude µ ∗ = µ I d x and d x → d x = d y ) the localmaximum of the radiation becomes more pronounced ifthe area of the impurity increases, as shown in Fig. 9.Also, the decreasing of the impurity area makes the localmaximum less pronounced. The main maximum is dom-inant, although there exist secondary local maxima, tothe left from the main maximum, although they are verysmall. The position of the main maximum shifts to theleft as the impurity size is decreased; however this shift isinsignificant even if the area d x d y / FIG. 9: Total emitted energy (normalized to µ I ) as a func-tion of the fluxon velocity for the square rhombic ( d x = d y )impurity at µ H = 0, µ C = 0 . µ C /µ I = 0 . d x = d y = 5(curve 1), d x = d y = 10 (curve 2), d x = d y = 20 (curve 3).The dashed lines corresponds to the case µ H = µ C = 0 and d x = d y = 5 (curve 4), d x = d y = 10 (curve 5), d x = d y = 20(curve 6). The inset gives the details of curve 3 on the largerscale. ( d x , d y →
0, and µ C,H →
0) brings the spectral energydensity function (38) to the already known limit of thepoint-like impurity . The same limit can be obtainedfrom any of the Eqs. (39) by setting d y → µ C,H → d x → µ C,H → µ ∗ = µ I d y or µ ∗ = µ I d x , respectively.The energy density profiles E ( q x , q y ) that correspondto the rhombic impurity are presented in Fig. 10. Asa particular example, we consider an impurity that cor-responds to the curve 3 from Fig. 7, i.e., for d x = 12, d y = 3. This energy density distribution bears manyqualitative similarities with the energy density functionfor the rectangular impurity shown in Fig. 5. Theglobal minima of the energy density satisfy the condition2 FIG. 10: (Color online). Emitted energy density E for therhombic impurity with d x = 12, d y = 3, µ C = µ H = 0 (curve3 in Fig. 7) at v = 0 . v = 0 .
47 (b), v = 0 .
63 (c), v = 0 . v = 0 .
85 (e), v = 0 .
95 (f), and v = 0 .
993 (g). The panel(h) corresponds to v = 0 .
993 and µ C = 0 . µ C /µ I = 0 . E ( q x , q y ) = 0 and are given by the set of equations d x q q x + q y v − q x ± q y d y = 4 πn ∓ . (40)This set of equations describes the sequence of pairs ofellipses that are numbered by the integers n ± , n ∓ = n , n + 1 , . . . , n = d x π s v − d x + d y d x . (41) if | v | < d x q d x + d y . (42)Otherwise, the Eqs. (40) yield the set of hyperbolas thatare numbered with n ± = ± , ± , . . . . The two curves (el-lipses or hyperbolas) given by Eq. (40) that correspondto the opposite signs but with n + = n − are mapped intoeach other with the mirror symmetry with respect to the q x axis. If we consider the set of curves with the samesign, say +, they are embedded into each other and theyexpand with the growth of the index n + . Between thesecurves lie the ridges of the E ( q x , q y ) function, and thelocal maxima of the energy density lie on these ridges.The signatures of these curves can be spotted in allpanels of Fig. 10. For small and intermediate values ofthe fluxon velocity the emitted radiation is localized pre-dominantly in one peak in the q − space, as shown in Figs.10(a-d). This peak lies on the q x axis; thus, most of theradiation does not propagate in the perpendicular direc-tion. In the panel (a) one can observe the distribution forthe rather small value of the fluxon velocity ( v = 0 .
2) andthis distribution is close to being radial. At such smallvelocities the pair of ellipses (40) with n − = n + are veryclose to being circles and almost coincide with each other.For larger values of v these pairs start to separate, as il-lustrated in Figs. 10(b-e). The panel (b) corresponds tothe local minimum of E ( v ) (curve 3 of Fig. 7) at v = 0 . v = 0 .
63. The structure of both these functions is similarand the only difference is that the maximal peak in panel(c) lies in the area of backward radiation ( q x ≈ − . q x axis at q x ≈ .
3. Thus, for the intermediate ve-locities the situation is similar to the case of rectangularimpurity, where the minimum of E ( v ) corresponded tothe minimal forward radiation. Panel (d) corresponds tothe next local minimum of the E ( v ) curve at v = 0 . v leads to the appearance of the pairof equivalent local maxima off the q y = 0 axis [see panel(e)]. These maxima become global as v approaches thevalue v = 1 [see panels (f) and (g)]. Thus, we observethe increasing of the perpendicular radiation that reachesits climax in the relativistic limit v →
1. Panel (f) cor-responds to the maximum of the E ( v ) function (curve3 of Fig. 7) at v = 0 . q x axis. They appear to be strongly localized in the q y di-rection while their localization in the q x is significantlyweaker. In this limit the interaction time with the tip ofthe rhombus is too small to generate significant longitu-dinal radiation, and the shape of the obstacle breaks theincident fluxon as a groin and generates predominantlytransverse radiation.3Finally, we mention the dependence of the emitted en-ergy on the parameter µ C . Panel (h) corresponds to thesame parameters of the model as in panel (g) but with µ C >
0. Comparing panels (g) and (h) we see that thestructure of these functions is very similar while the ab-solute values of E are significantly smaller in the µ C = 0case. If µ C = 0, but for the same value of the fluxon ve-locity, the values of the maxima actually decrease with v .Thus, the total emitted energy tends to zero, in the sameway as shown by the dashed lines in Fig. 9. This hasbeen confirmed for the values of v even closer to unity aswell as for the different values of d x,y . The qualitativebehavior of E ( v ) in the limit | v | → IV. DISCUSSION AND CONCLUSION
The radiation emitted as a result of the fluxon interac-tion with the impurity of a general geometrical shape inthe large two-dimensional Josephson junction has beenstudied. The emitted energy distribution in the q − spacehas been computed as well as the total emitted energy.This energy distribution can always be represented as atriple integral. In principle, any geometrical shape canbe taken into consideration; however, the explicit integra-tion is not always possible, but if the inhomogeneity areacan be represented by the piecewise-linear functions, thisintegration can be done. In this article the rectangularand rhombic impurities have been studied.The main result of this work has been formulated in thedependence E ( v ) of the total emitted energy as a functionof the incident fluxon velocity. It appears that this de-pendence has local maxima that depend strongly on thegeometric properties of the impurity. These local max-ima do not exist if the impurity is treated as a point .Controlling the shape of the impurity one can remove theextrema or make them more pronounced. The limit ofthe 1D problem with the finite-size inhomogeneity canbe restored.First of all we would like to mention the differencesbetween the 1D and 2D cases. The 1D case appears tobe the limit of the 2D rectangular impurity case when d y → ∞ . While moving away from the 1D limit by de-creasing d y we observe gradual lowering and disappear-ance of the extrema of the E ( v ) dependence. Next, the2D model allows to take into account the impurity shapesthat are different from the rectangle. For the rhombicimpurity we have demonstrated that the emitted energydependence on the fluxon velocity is rather different fromthe rectangular case and does not possess the 1D limit.In principle, other geometrical impurity shapes can bestudied, including the asymmetric ones.In this article the junction thickness change due to thehomogeneity is taken into account. Its role is measured by the parameters µ H and µ C [see Eqs. (3) and (6)]. Theparameter µ C is responsible for the capacitance changeand plays the dominant role. In some papers these pa-rameters are ignored (especially they are always ignoredif the point impurities are considered), and, in general,are considered to be weak . However, the junction thick-ness change influences significantly the asymptotic be-havior of the total emitted energy in the “relativistic”(i.e., v →
1) limit of the fluxon velocity. If the thick-ness change is ignored, the total emitted energy goes tozero, while it exhibits unbounded growth if the thick-ness change is taken into account. This is true for theboth 1D and 2D junctions. The emitted energy has beencomputed under the assumption that it is a small pertur-bation on the fluxon background. Consideration of thehigher order corrections may block the infinite radiationgrowth. Also, the dissipative effects, which have beenignored in this work, should contribute to the decreasingof the emitted energy.Although the real large-area Josephson junctions havefinite dimensions, in this article the infinitely-sized junc-tion has been considered. This approximation is suffi-cient if the physical dimensions of the LJJ exceed by theorder of magnitude the Josephson penetration depth and,consequently, the fluxon length in the x direction. Theboundary conditions are also important, however , ifthe junction width is large enough (exceeds the Joseph-son length at least by the order of magnitude) the fluxondistortion from the linear shape is insignificant. In anycase, before focusing on the more concrete setup an ide-alized, but more easily solvable model should be studied.Finally, we discuss the possible application of the ob-tained results. Recently, a number of papers have fo-cused on the different application of the fluxon dynam-ics in the 2D LJJ, such as fluxon splitting on the T-shaped junctions , excitation of the different modes thatmove along the fluxon front , and the fluxon logicgates where the interaction with the spatial inhomo-geneity takes place. If the incident fluxon velocity islarge enough, the emitted radiation becomes sufficientand it should influence the fluxon motion. In particular,the non-monotonicity of the E ( v ) dependence may pro-duce the hysteresis-like branches on the current-voltagecharacteristics (IVCs) of the LJJ. Studies of these IVCsfor the different shapes of the inhomogeneity in the gen-uinely 2D case are in progress and will be published else-where. Acknowledgemets
Y.Z. acknowledges financial support fromUkrainian State Grant for Fundamental ResearchNo. 0112U000056.4 ∗ Electronic address: [email protected] † Electronic address: [email protected] A. Barone and G. Paterno,
Physics and Applications of theJosephson Effect (Wiley, New York, 1982). K. K. Likharev,
Dynamics of Josephson Junctions and Cir-cuits (Gordon and Breach, New York, 1986). D. W. McLaughlin and A. C. Scott, Phys. Rev. A , 1652(1978). Y. S. Gal’pern and A. T. Filippov, Sov. Phys. JETP ,1527 (1984). L. G. Aslamazov and E. V. Gurovich, JETP Lett. , 746(1984). R. Fehrenbacher, V. B. Geshkenbein, and G. Blatter, Phys.Rev. B , 5450 (1992). L. Balents and S. H. Simon, Phys. Rev. B , 6515 (1995). M. V. Fistul and G. F. Giuliani, Phys. Rev. B , 9348(1998). G. Mkrtchyan and V. V. Shmidt, Solit State Commun. ,791 (1979). Y. S. Kivshar, B. A. Malomed, and A. A. Nepomnyashchy,JETP , 356 (1988). B. A. Malomed, I. L. Serpuchenko, M. I. Tribelsky, andA. V. Ustinov, JETP Lett. , 591 (1988). B. A. Malomed and A. V. Ustinov, J. Appl. Phys. , 3791(1990). Y. S. Kivshar and O. A. Chubykalo, Phys. Rev. B , 5419(1991). Y. S. Kivshar, A. M. Kosevich, and O. A. Chubykalo, Phys.Lett. A , 449 (1988). Y. S. Kivshar, A. M. Kosevich, and O. A. Chubykalo, Fiz.Nizk. Temp , 800 (1987). P. S. Lomdahl, O. H. Olsen, J. C. Eilbeck, and M. R.Samuelsen, J. Appl. Phys. , 997 (1985). S. G. Lachenmann, G. Filatrella, T. Doderer, J. C. Fernan-dez, and R. P. Huebener, Phys. Rev. B , 16623 (1993). I. O. Starodub and Y. Zolotaryuk, Phys. Lett. A , 3101(2012). D. R. Gulevich, F. V. Kusmartsev, S. Savel’ev, V. A. Yam-pol’skii, and F. Nori, Phys. Rev. B , 094509(13) (2009). H. Nacak and F. Kusmartsev, Physica C , 827 (2010). B. A. Malomed, Physica D , 157 (1991). Y. S. Kivshar and B. A. Malomed, Phys. Lett. A , 443(1988). J. Rubinstein, J. Math. Phys. , 258 (1970). M. Salerno, E. Joergensen, and M. Samuelsen, Phys. Rev.B , 2635 (1984). R. L. Graham, D. E. Knuth, and O. Patashnik,
ConcreteMathematics (Addison-Wesley, Reading Ma., 1994). I. O. Starodub and Y. Zolotaryuk, Ukr. J. Phys. , 687(2013). J. C. Eilbeck, P. S. Lomdahl, O. H. Olsen, and M. R.Samuelsen, J. Appl. Phys. , 861 (1985). D. R. Gulevich and F. V. Kusmartsev, Phys. Rev. Lett. , 017004 (2006). D. R. Gulevich, F. V. Kusmartsev, S. Savel’ev, V. A. Yam-pol’skii, and F. Nori, Phys. Rev. Lett.101