Table Of Content

Head-to-Head Domain Wall Structures in Wide Permalloy Strips Virginia Estévez∗ and Lasse Laurson COMP Centre of Excellence and Helsinki Institute of Physics, Department of Applied Physics, Aalto University, P.O.Box 11100, FI-00076 Aalto, Espoo, Finland. WeanalyzetheequilibriummicromagneticdomainwallstructuresencounteredinPermalloystrips of a wide range of thicknesses and widths, with strip widths up to several micrometers. By per- forminganextensivesetofmicromagneticsimulations,weshowthattheequilibriumphasediagram of the domain wall structures exhibits in addition to the previously found structures (symmetric 5 and asymmetric transverse walls, vortex wall) also double vortex and triple vortex domain walls 1 for large enough strip widths and thicknesses. Also several metastable domain wall structures are 0 found for wide and/or thick strips. We discuss the details of the relaxation process from random 2 magnetizationinitialstatestowardsthestabledomainwallstructure,andshowthatourresultsare robust with respect to changes of e.g. the magnitude of the Gilbert damping constant and details n of the initial conditions. a J PACSnumbers: 75.60.Ch,75.78.Cd 0 3 I. INTRODUCTION ergy may be found18,28–31. For H >HW or J >JW, the ] DW structures exhibit dynamical evolution: for TWs, l l repeated nucleation and propagation of an antivortex a During the last decade, a lot of effort has been de- h across the strip width takes place19. Similarly, in VWs voted to understand static and dynamic properties of - the vortex core performs oscillatory back and forth per- s magnetic domain walls (DWs) in ferromagnetic nanos- e tructures such as nanowires and -strips. These studies pendicular motion19. m have been largely driven by promising technological ap- In Permalloy strips with even larger widths and/or . plications based on domain walls and their dynamics, in thicknesses, one might expect also other, possibly more t a particular memory1,2 and logic devices3–5. In typical ex- complicatedequilibriumDWstructures. Forwiderstrips m periments DWs are driven by either applied magnetic shape anisotropy is less important, implying that energy - fields6,7 or spin-polarized electric currents8–11. The re- minima with more complex spin structures closing the d sulting DW dynamics depends crucially on the micro- flux more efficiently than TWs, ATWs or VWs may ap- n magnetic DW structure, typically involving various in- pear. Indeed, e.g. double and triple vortex DWs have o c ternal degrees of freedom. These are essential e.g. for been observed in experiments on wide strips28, but they [ the emergence of the Walker breakdown12, an instability have been attributed to current-induced vortex nucle- occurring when the DW internal degrees of freedom get ation resulting in metastable DW structures. Conse- 1 v excited by a strong external drive (a magnetic field H quently, a pertinent and fundamental question is what 1 or a spin-polarized current J exceeding the Walker field are the possible intermediate equilibrium DW structures 3 H or current J , respectively), limiting the propaga- observable when the lateral Permalloy strip dimensions W W 7 tion velocity of the DWs. increase from those corresponding to the typical nanos- 7 tripgeometry(withTW,ATWorVWasthestableDW Twomainclassesofferromagneticmaterialshavebeen 0 structure) to strip widths of micrometers and beyond. extensively studied within the strip geometry. Ma- . 1 terials with a high perpendicular magnetic anisotropy In this paper we present an extensive numerical study 0 (PMA13–16) exhibit simple and narrow DWs of the of the equilibrium and metastable micromagnetic DW 5 Bloch and/or Néel type. For H > H or J > J , structuresinPermalloystrips,withthestripwidthsupto 1 W W : repeated transitions between these two structures are an order of magnitude larger than before19–24. Contrary v observed16. The second class of systems includes soft to previous studies focusing on comparing the energies i X (lowanisotropy)magneticmaterials13,14 suchasPermal- of different a priori known DW structures19–24, we per- loy, where in-plane domain magnetization along the long form micromagnetic simulations of relaxation dynamics r a axis of the strip is induced by shape anisotropy. By us- from random initial states towards the stable DW struc- ing various experimental techniques17,18 and micromag- tures. In addition to the previously observed TW, ATW netic simulations, it has been established that the equi- andVWDWs,wefindalsoDWswithequilibriumdouble libriumDWstructuresseparatingthesein-planedomains andtriplevortexstructuresforwideand/orthickenough are more complex, and depend crucially on the sample strips. Thelaststructureisencounteredonlyinthevery geometry19–22. Transverse DWs (TWs) and asymmetric largest system sizes we were able to simulate. Moreover, transverse DWs (ATWs) are observed for narrow and for wide strips we find a rich variety of metastable DWs thin strips22–27, while in wider and thicker strips one with even more complex micromagnetic structures. We encounters the vortex DW (VW)22,24–26,28,29. In addi- demonstrate that our results are robust with respect to tion, various metastable DW structures with higher en- changes of the magnitude of the Gilbert damping con- 2 (a) thickness Δ 0 ns 0.2 ns z width w length l 1.0 ns 2.0 ns l (b) r 3.0 ns 4.0 ns FIG. 1. (color online) (a) Geometry of the Permalloy strip. (b)Atopviewofthemagnetizationintheinitialstate. Mag- 6.0 ns 12.0 ns netization points along the long axis of the strip within the two domains (as indicated by the arrows) forming a head- FIG. 2. (color online) An example of the temporal evolution to-head configuration. In between them, a region of random of the relaxation, with w=420 nm, ∆ =10 nm, α=3 and magnetization (of length l ) has been included. z r l = 2w. Relaxation towards the equilibrum DW structure r (here, a VW) takes place via coarsening dynamics of the de- fectstructureinthemagnetictexture. Thecolorwheelinthe stant or using different initial conditions for the relax- middle shows the mapping between magnetization directions ation process. Our results underline the crucial role of and colors. topological defects for physics of DWs in soft strips, and that of micromagnetic simulations for finding the true equilibrium DW structure. contributions due to exchange, Zeeman, and demagne- tizing energies. The size of the discretization cell used depends on the system size, but is always bounded by II. MICROMAGNETIC SIMULATIONS the exchange length, Λ=(2A/µ0Ms2)1/2 ≈5 nm, in the in-planedirections,andequals∆ inthetheout-of-plane z direction. ThesystemstudiedisaPermalloystripofwidthwand thickness ∆ , satisfying ∆ (cid:28) w, see Fig. 1 (a). In the z z micromagnetic simulations, magnetic charges are com- III. RESULTS pensatedontheleftandrightendsofthestrip,tomimic an infinitely long strip; the actual simulated length sat- isfies l ≥ 4w for all cases considered. The initial state We start by considering the effect of varying α and lr from which the relaxation towards a stable DW struc- on the relaxation process. Fig. 2 shows an example of ture starts is an in-plane head-to-head domain struc- the time evolution of m(r,t) for w = 420 nm, ∆z = 10 ture, with a region of random magnetization of length nm, α = 3 and lr = 2w. The initially random magne- l in the middle of the sample, see Fig. 1 (b). If not tization evolves via coarsening of the defect structure of r specified otherwise, we consider l = 2w. Material pa- the magnetization texture towards the stable DW (here, r rameters of Permalloy are used, i.e. saturation magne- a VW). During the relaxation, the total energy E of the tization M = 860 × 103 A/m and exchange constant system decreases in a manner that for a given geometry s Aex = 13 × 10−12 J/m. The typical Gilbert damping (w and ∆z) depends on both α and lr, see Fig. 3 (a) and constant for Permalloy is α = 0.01, but here we analyze (b) where w =5120 nm and ∆z =20 nm, is considered. alsotheinfluenceofαontherelaxationprocess,andthus For instance, E decreases faster for an intermediate α consideralsoothervalues. Forsimplicity,wesetthetem- [Fig. 3 (a)]. We attribute this behavior to the balance perature T to zero, and focus on the ideal case of strips between inertial effects related to precession favored by free of any structural disorder or impurities. a small α, helping to overcome energy barriers, and the The simulations are performed using the GPU- higher rate of energy dissipation due to a large α. Thus, accelerated micromagnetic code MuMax332–34, offering the relaxation time to reach a (meta)stable DW struc- a significant speedup as compared to CPU codes for the ture depends on α. Fig. 3 (b) illustrates that for a fixed largesystemsizesweconsiderhere. Tocalculatethemag- α, systems with a larger lr relax more slowly. Fig. 3 netization dynamics of the system, the Landau-Lifshitz- (c) shows that on average, the early-time relaxation of Gilbert equation35,36, E towards its final value Ef exhibits temporal power-law decay, (cid:104)E−E(cid:105)∝t−β with β ≈1.3 for the α=0.3 case f ∂m/∂t=γH ×m+αm×∂m/∂t, (1) shown, possibly related to collective effects due to inter- eff actions between several topological defects during early is solved numerically. Here, m is the magnetization, γ stages of relaxation (Fig. 2). the gyromagnetic ratio, and H the effective field, with In general, the final (meta)stable DW structure may eff 3 9×10-13 (a) α = 3 (b) l = w 6×10-13 α = 0.3 lr = 2w 25 E [J]63××1100--1133 lr =α 2 =w 0.03 α l=rr =0 .33w 3×10-13E [J] 20 Vortex wall DTroipulbel ev ovrotretxe xw walalll ] m n15 0 0 Δt [z Asymmetric transverse wall 0 t [s]2×10-9 0 t [s] 2×10-9 10 (0)> 100 (c) 5 E Transverse wall E> / <f1100--21 lllrr === 23www 100 w [nm10]00 - r α = 0.3 Et() 10-3 t-1.3 FIG. 4. (color online) Phase diagram of the equilibrium < 10-10 10-9 10-8 DWstructureinPermalloystripsofvariousthicknesses(from t [s] ∆z =5to25nm)andwidthsrangingfromw=120nmupto 5120nm. Thesymbolscorrespondtoobservationsofthevari- ousequilibriumDWstructures,withphaseboundariesshown FIG. 3. (color online) The energy E(t) as a function of time as solid lines. Examples of the DW structures corresponding t for w=5120 nm and ∆ =20 nm. (a) For different values z to the 5 different phases are shown in Fig. 5. of α and l = 2w. (b) For different values of l , and α = 0.3 r r [resulting in the fastest relaxation in (a)]. (c) shows that on average, the early time decay of E(t) towards its final value Transverse Wall E obeys(cid:104)E(t)−E(cid:105)∝t−β. Fortheα=0.3caseshownhere, f f β ≈1.3. Empty(filled)symbolsin(c)correspondtow=420, ∆ =10 nm (w=860, ∆ =20 nm). z z (a) depend on the realization of the random initial state. Asymmetric Transverse Wall Vortex Wall Thus, we consider 21 realizations of the initial random magnetization for each w and ∆ , and compare the en- z ergies of the resulting relaxed configurations. The struc- (b) (c) ture with the lowest energy is chosen as the equilibrium structure,whileotherswithhigherenergyaremetastable Double Vortex Wall Triple Vortex Wall states. Although, as discussed above, the relaxation times depend on α and l , the equilibrium DW struc- r ture is found to be independent of α and l in the range r (d) (e) considered, i.e. α ∈ [0.01,3] and l ∈ [w,3w]. Thus, in r what follows, we will use α=3 and l =2w. r The main results of this paper are summarized in FIG. 5. (color online) Examples of the different equilibrium Figs. 4 and 5, showing the phase diagram of the equi- micromagnetic DW structures: (a) TW for w =120 nm and librium DW structures for w ranging from 120 to 5120 ∆z =5 nm, (b) ATW for w =160 nm and ∆z =10 nm, (c) VWforw=640nmand∆ =15nm,(d)DVWforw=2560 nm, and ∆ from 5 to 25 nm, and examples of these z z nm and ∆ = 20 nm, and (e) TVW for w = 5120 nm and structures, respectively. Forsmallw, werecoverthepre- z ∆ = 25 nm. The colorwheel (top left) shows the mapping viousresults19–22,i.e. phasescorrespodingtoTW,ATW z between magnetization directions and colors. and VW, shown in Fig. 5 (a), (b) and (c), respectively. For larger strip widths (w approaching or exceeding 1µm,dependingon∆z,seeFig.4),anewequilibriummi- is found for the very largest system sizes we have been cromagneticDWstructure,adoublevortexwall(DVW), able to simulate. The middle vortex of the TVW has an is observed. This structure consists of two vortices with oppositesenseofrotationtotheothertwo. Forw =5120 opposite sense of rotation of the magnetization around and∆ =25nm, DVWandTVWhavethesameenergy z the vortex core, see Fig. 5 (d). At the phase boundary (thecyansquaresymbolinthetoprightcornerofFig.4), (blue triangle symbols pointing left in Fig. 4), VW and suggestingthepresenceofaphaseboundarybetweenthe DVW have the same energy. The DVW phase spans a two structures. Indeed, by performing a set of 10 ad- relatively large area within the (w,∆z) space, highlight- ditional simulations with w = 6144 and ∆z = 25 nm ing the robustness of our results. (i.e. outside the phase diagram in Fig. 4), suggests that In addition, a second new phase, with a triple-vortex TVWistheequilibriumDWstructureforverylargestrip wall (TVW) as the equilibrium structure [see Fig. 5 (e)], widths. Thisstructurehasbeenobservedinexperiments 4 2V+AV 3V+AV vortices(withthesamesenseofrotation)andanantivor- (a) (b) tex [2V+AV, see Fig. 6 (a)], there are two vortices and only two edge defects. Thus, in order to compensate the topological defects, also an antivortex appears. In general,wehaveobservedthatinaDWwithN vorticeswith the same sense of rotation, there must be N−1 antivor- 3V+AV 3V+2AV tices to get a zero total winding number, see Fig. 6 (a) (c) (d) and (d) for examples with 2V+AV and 3V+2AV configurations, respectively. When some of the vortices have oppposite sense of rotation, more complex scenarios are encountered,withexamplesshowninFigs.6(b),(c),(e), 4V+AV 4V+2AV (f),(g)and(h). NoticealsothattwoDWstructureswith the same elements can look very different, see e.g. the two 3V+AV DWs shown in Figs. 6 (b) and (c). All the (e) (f) DWstructuresfound,boththeequilibriumonesinFig.5 and the metastable states in Fig. 6 obey the principle of 4V+3AV 5V+2AV compensationoftopologicaldefectstoyieldatotalwind- ingnumberofzero. Therichnessoftheequilibriumphase diagram and the large collection of metastable states in- (g) (h) dicatethatforwide/thickstripsinparticular,themicro- magneticenergylandscapeisquitecomplex,withalarge number of local minima. This is also in agreement with FIG.6. (coloronline)ExamplesofmetastableDWstructures our observations of power-law energy relaxation. observed for a system with w=5120 nm and different thicknesses ∆ : (a) Two vortices and an antivortex (2V+AV), z ∆ = 25 nm, (b) and (c) three vortices and an antivortex z (3V+AV),∆z =5nm,(d)threevorticesandtwoantivortices IV. SUMMARY AND CONCLUSIONS (3V+2AV), ∆ =10 nm, (e) four vortices and an antivortex z (4V+AV), ∆ = 5 nm, (g) four vortices and three antivor- z To summarize, we have performed an extensive set tices (4V+3AV), ∆ = 5 nm, and (h) five vortices and two z of micromagnetic simulations to study the equilibrium antivortices (5V+2AV), ∆ =5 nm. z and metastable DW structures in Permalloy strips of a wide range of widths and thicknesses, as well as the relaxation dynamics starting from random magnetization asametastablestateforsmallersystems28,29. Noticealso initial states. The general trend of our results is that that the middle part of the TVW [Fig. 5 (e)], exhibiting both the equilibrium and metastable DW configurations four line-like 90◦ DWs meeting at a vortex core in the become increasingly complex (i.e. they consist of an in- middle of the TVW, resembles the typical Landau flux- creasingnumberoftopologicaldefects)asthelateralstrip closure magnetization patterns observed for rectangular dimensions increase. We note that somewhat analogous Permalloy thin films37–39. behaviour - i.e. existence of equilibrium magnetization Followingtherelaxationfromarandommagnetization configurations with increasing complexity as the system initial state, the system may in general end up into var- size increases - is observed also in some other systems ious metastable states with higher energy than that of such as three-dimensional cylindrical elements with per- the equilibrium DW. Sometimes these metastable states pendicular anisotropy41,42. have even a higher probability than the equilibrium one, Several remarks are in order: first, for strips with estimated here from the sample of 21 relaxed configura- even larger lateral dimensions one may in principle ex- tions. Fig. 6 shows some of the metastable states found pect more complex DW patterns - possibly with four or foralargestripwithw =5120nmanddifferentvaluesof more vortices with alternating sense of rotation. These, ∆ ; for strips with smaller lateral dimensions, different however, are currently beyond the reach of our available z metastable states tend to be less numerous and have a computing resources. Second, our phase diagram allows simpler structure. Despite their apparent complexity, all one to check if experimental observations of the various themetastableDWstructuresshowninFig.6respectthe DW structures in wider strips are equilibrium configura- basic principles of topology of DWs. Each of the DWs tions or metastable states. According to our review of are composed of topological defects, with an associated the experimental literature, most observations of DVWs windingnumber: +1forvortices, -1forantivortices, and and TVWs appear to be metastable states28,29. Third, ±1/2 for edge defects40. In a DW all the topological whiletheequilibriumstructureswefindarecertainlysta- defects have to be compensated, i.e. the total winding ble in the absence of external perturbations such as ap- number is equal to zero. In the case of the DVW, the pliedmagneticfields,itremainstobeseenhowtheirfield two topological vortex defects are compensated by four driven dynamics is like, and whether wide strips with a edge defects [Fig. 5 (d)]. For the metastable state of two relatively weak shape ansitropy are able to support the 5 DWs as compact objects also when external perturba- ACKNOWLEDGMENTS tions are being applied43. We thank Mikko J. Alava for a critical reading of the manuscript. This work has been supported by the Academy of Finland through its Centres of Excel- lence Programme (2012-2017) under project no. 251748, and an Academy Research Fellowship (LL, project no. 268302). We acknowledge the computational resources provided by the Aalto University School of Science “Science-IT” project, as well as those provided by CSC (Finland). ∗ virginia.esteveznuno@aalto.fi 33, 4167 (1997). 1 S. S. P. Parkin, M. Hayashi and L. Thomas, Science 320, 21 Y. Nakatani, A. Thiaville and J. Miltat, J. Magn. Magn. 190 (2008). Mater. 290-291 750 (2005). 2 S. E. Barnes, J. Ieda and S. Maekawa, Appl. Phys. Lett. 22 M. Klaüi, C. A. F. Vaz, J. A. C. Bland, L. J. Heyderman, 89, 122507 (2006). F. Nolting, A. Paviovska, E. Bauer, S. Cherifi, S. Heun, 3 R. P. Cowburn and M. E. Welland, Science 287, 1466 and A. Locatelli, Appl. Phys. Lett. 85, 5637 (2004). (2000). 23 S. Jamet, N. Rougemaille, J. C. Toussaint, and O. 4 D. A. Allwood, G. Xiong, M. D. Cooke, C. C. Faulkner, Fruchart, arXiv preprint arXiv:1412.0679 (2014). D. Atkinson, N. Vernier and R. P. Cowburn, Science 296, 24 Y.Nozaki,K.Matsuyama,T.Ono,andH.Miyajima,Jpn. 2003 (2002) J. Appl. Phys. Vol. 38 (1999) 6282. 5 D.A.Allwood,G.Xiong,C.C.Faulkner,D.Atkinson,D. 25 M. Kläui, C. A. F. Vaz, J. A. C. Bland, T. L. Monchesky, Petit and R. P. Cowburn, Science 309, 1688 (2005). J.Unguris,E.Bauer,S.Cherifi,S.Heun,A.Locatelli,L.J. 6 T. Ono, H. Miyajima, K. Shigeto, K. Mibu, N. Hosoito Heyderman, and Z. Cui, Phys. Rev. B 68, 134426 (2003). and T. Shinjo, Science 284, 468 (1999). 26 M. Kläui, U. Ru¨diger, C. A. F. Vaz, J. A. C. Bland, S. 7 G. S. D. Beach, C. Nistor, C. Knutson, M. Tsoi and J. L. Cherifi,A.Locatelli,S.Heun,A.Pavlovska,E.Bauer,and Erskine, Nature Mater. 4, 741 (2005). L. J. Heyderman, Appl. Phys. Lett. 99, 08G308 (2006). 8 M. Klaüi, C. A. F. Vaz, J. A. C. Bland, W. Wernsdorfer, 27 D.Backes,C.Schieback,M.Klaüi,F.Junginger,H.Ehrke, G. Faini, E. Cambril and L. J. Heyderman, Appl. Phys. P. Nielaba, U. Ru¨diger, L. J. Heyderman, C. S. Chen, T. Lett. 83, 105 (2003). Kasama, R. E. Dunin-Borkowski, C. A. F. Vaz, and J. A. 9 A. Thiaville, Y. Nakatani, J. Miltat and N. Vernier, J. C. Bland, Appl. Phys. Lett. 91, 112502 (2007). Appl. Phys. 95, 7049 (2004). 28 M.Klaüi,M.Laufenberg,L.Heyne,D.Backes,U.Ru¨diger, 10 M. Klaüi, P.-O. Jubert, R. Allenspach, A. Bischof, J. A. C.A.F.Vaz,J.A.C.Bland,L.J.Heyderman,S.Cherifi, C. Bland, G. Faini, U. Ru¨diger, C. A. F. Vaz, L. Vila and A. Locatelli, T. O. Mentes, and L. Aballe, Appl. Phys. C. Vouille, Phys. Rev. Lett. 95, 026601 (2005). Lett. 88 232507 (2006). 11 N.Vernier,D.A.Allwood,D.Atkinson,M.D.Cookeand 29 E.-M.Hempe,M.Klaüi,T.Kasama,D.Backes,F.Jungin- R. P. Cowburn, EPL 65, 526 (2004). ger,S.Krzyk,L.J.Heyderman,R.Dunin-Borkowski,and 12 N.L.SchryerandL.R.Walker.JournalofAppliedPhysics U. Ru¨diger, phys. stat. sol.(a) 204, 3922 (2007). 45,5406 (1074). 30 C. A. F. Vaz, M. Klaüi, L. J. Heyderman, C. David, F. 13 K.H.J. Buschow and F.R. de Boer, Physics of Mag- Nolting, and J. A. C. Bland, Phys. Rev. B. 72, 224426 netism and Magnetic Materials. Kluwer Academic Pub- (2005). lishers (2003). 31 M.H.Park,Y.K.Hong,B.C.Choi,M.J.Donahue,H.Han, 14 J.M.D. Coey, Magnetism and Magnetig Materials, Cam- and S.H. Gee, Phys. Rev. B 73, 094424 (2006). bridge University Press. 32 http://mumax.github.io 15 O. Boulle, J. Kimling,P. Warnicke, M. Klaüi, U. R. 33 A. Vansteenkiste and B. Van de Wiele, J. Magn. Magn. Ru¨diger, G. Malinowski, H. J. M. Swagten, B. Koop- Mater. 323, 2585 (2011). mans,C.UlysseandG.Faini.Phys.Rev.Lett.101216601 34 A. Vansteenkiste et al., AIP Advances 4, 107133 (2014). (2008). 35 T. L. Gilbert. IEEE Trans. Magn. 40, 3443 (2004). 16 E.Martinez.J.Phys.: Condens.Matter24(2012)024206. 36 W. F. Brown, Jr., Micromagnetics. (New York, Wiley, 17 H. Kronmu¨ller and S. S. P. Parkin (eds), Handbook of 1963). Magnetism and Advanced Materials 3 (Chichester, Wiley, 37 T.Shinjo,T.Okuno,R.Hassdorf,K.Shigeto,andT.Ono, 2007). Science 289, 930 (2000). 18 M. Klaüi, J. Phys.: Condens. Matter 20, 313001 (2008). 38 A Wachowiak, J. Wiebe, M. Bode, O. Pietzsch, M. Mor- 19 A. Thiaville and Y. Nakatani, in Spin Dynamics in Con- genstern, and R. Wiesendanger, Science 298, 577 (2002). finedmagneticStructuresIII,editedbyB.Hillebrandsand 39 P. Fischer, M.-Y. Im, S. Kasai, K. Yamada, T. Ono, and A. Thiaville, Topics in Applied Physics 101 (Springer, A. Thiaville, Phys. Rev. B 83, 212402 (2011). Berlin, 2006). 40 O. Tchernyshyov and G.-W. Chern, Phys. Rev. Lett., 95, 20 R.D.McMichaelandM.J.Donahue,IEEETrans.Magn. 197204 (2005). 6 41 C. Moutafis, S. Komineas, C. A. F. Vaz, J. A. C. Bland, (2008). T. Shima, T. Seki, and K. Takanashi, Phys. Rev. B 76, 43 C. Zinoni, A. Vanhaverbeke, P. Eib, G. Salis, and R. Al- 104426 (2007). lenspach, Phys. Rev. Lett. 107, 207204 (2011). 42 N. Vukadinovic and F. Boust, Phys. Rev. B 78, 184411