Particle-Based Computer Simulations of Magnetic Skyrmions Dissertation zur Erlangung des Grades „Doktor der Naturwissenschaften“ am Fachbereich Physik, Mathematik und Informatik der Johannes Gutenberg-Universität in Mainz Jan Rothörl geb. in Wiesbaden Mainz, den 21.09.2023 Abstract Magnetic skyrmions are topologically stabilized magnetic structures found in multilayer thin film systems. They are often described as two-dimensional quasi- particles evolving according to the Thiele equation. This description allows for computer simulations of skyrmion dynamics using Brownian dynamics simulations with an additional term that takes care of the skyrmion Hall effect. Such simulations are significantly more efficient than conventional methods such as micromagnetic or atomistic spin dynamics simulations and therefore allow for simulating larger skyrmion systems for longer. Modeling specific skyrmion systems with this equation requires knowledge of interaction potentials of skyrmions with each other and with magnetic material boundaries. This thesis presents an approach to determine these potentials directly from experiment without prior assumptions on the potential shape by using the iterative Boltzmann inversion method, which has been previously established in soft-matter systems. The interaction potentials are found to be purely repulsive for the micrometer sized skyrmions studied here. Based on these potentials, further simulations of skyrmion systems are performed showing a good agreement with experiments and allowing for a comparison of simulated and experimental skyrmion lattices. Besides the analysis of phase transitions and ordering, the analysis of large skyrmion systems also comprises the effect of the gyroscopic Magnus force on diffusion of skyrmions. Contrary to other systems containing diffusive particles, an increased skyrmion density can lead to an increase in diffusion if a sufficiently strong Magnus force is present. Moreover, pinning effects and their influence on free and current-induced skyrmion diffusion are investigated. This investigation motivates the development of two different methods to increase skyrmion diffusion at constant temperature by reducing the impact of pinning on skyrmions. The first method uses periodic perpendicular magnetic field excitations to change skyrmion sizes and thereby the effective pinning landscape in which skyrmions move. The second one uses periodic current excitations to move skyrmions out of pinning sites. Finally, the effect of tight confinements on skyrmion dynamics is investigated. Here, the skyrmion dynamics is not only affected by the skyrmion density but also by the commensurability of the skyrmion number with the confinement. When the number of skyrmions is commensurate with the confinement geometry, diffusion is significantly decreased. Building upon these results, the changes in ordering and dynamics of confined skyrmions due to applied currents are analyzed in the context of employing such confined skyrmion systems for reservoir computing. iii Zusammenfassung Magnetische Skyrmionen sind topologisch stabilisierte magnetische Objekte, welche in Mehrschichtsystemen auftreten. Sie werden oft als zweidimensionale Quasiteil- chen beschrieben, deren Bewegung der Thielegleichung folgt. Diese Beschreibung ermöglicht Computersimulationen der Dynamik von Skyrmionen mit Hilfe von Brownschen Bewegungsgleichungen mit einem zusätzlichen Term, welcher den Skyrmion-Hall-Effekt berücksichtigt. Diese Simulationsmethode ist deutlich effizien- ter als konventionelle Simulationsmethoden wie Mikromagnetik oder atomistische Spindynamiksimulationen und ermöglicht daher deutlich längere Simulationen von größeren Skyrmionensystemen. Für die Modellierung spezifischer Systeme mit dieser Gleichung werden Wechselwirkungspotentiale von Skyrmionen miteinander sowie mit dem Rand des magnetischen Materials benötigt. Diese Arbeit stellt einen Ansatz zur Bestimmung dieser Potentiale direkt aus experimentellen Daten ohne weitere Annahmen an die Form des Potentials vor, welcher die Methode der iterativen Boltzmann-Inversion nutzt, die im Bereich der weichen Materie entwickelt wurde. Dabei ergeben sich für die hier verwendeten mikrometergroßen Skyrmionen vollständig repulsive Wechselwirkungspotentiale. Auf Basis dieser Potentiale wer- den weitere Simulationen von Skyrmionensystemen durchgeführt, welche eine gute Übereinstimmung mit Experimenten zeigen und daher einen Vergleich von expe- rimentellen und simulierten Skyrmionengittern ermöglichen. Neben der Analyse von Phasenübergängen und dem Ordnungsverhalten umfasst die Analyse großer Skyrmionensysteme auch eine Betrachtung der Auswirkung der gyroskopischen Magnuskraft auf die Diffusion von Skyrmionen. Im Unterschied zu anderen Syste- men diffusiver Teilchen kann eine erhöhte Dichte zu einer erhöhten Diffusion der Skyrmionen führen, falls die Magnuskraft stark genug ist. Darüber hinaus werden Pinningeffekte und deren Einfluss auf freie und strominduzierte Skyrmiondiffusion betrachtet. Dies motiviert die Entwicklung von zwei Methoden zur Erhöhung der Diffusivität bei konstanter Temperatur durch Reduktion des Einflusses von Pinning. Die erste Methode nutzt periodische Anregungen des angelegten magnetischen Feldes, um die Größe der Skyrmionen und damit die effektiv wirkende Pinningland- schaft zu verändern. Die zweite Methode nutzt periodische Stromanregungen, um Skyrmionen aus Regionen mit starkem Pinning herauszubewegen. Zuletzt wird die Wirkung von engem Confinement auf Skyrmionendynamik untersucht. Diese hängt nicht nur von der Skyrmionendichte, sondern auch von der Kommensurabilität der Skyrmionenzahl mit dem Confinement ab. Wenn die Skyrmionenzahl kom- mensurabel mit der Geometrie des Confinements ist, führt dies zu einer deutlich gesenkten Diffusion. Aufbauend auf diesen Ergebnissen werden auch die Verände- rungen der Ordnung und Dynamik von Skyrmionen in eingeschränkter Geometrie durch angelegte Ströme im Kontext der Nutzung eines solchen Skyrmionensystems als Reservoir-Computer analysiert. iv Contents Contents 1. Introduction 1 2. Theoretical Background 6 2.1. Magnetism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.1. Magnetic Energy Terms . . . . . . . . . . . . . . . . . . . . 6 2.1.1.1. Exchange Interaction . . . . . . . . . . . . . . . . . 6 2.1.1.2. Dzyaloshinskii-Moriya Interaction . . . . . . . . . . 9 2.1.1.3. Magnetocrystalline Anisotropy . . . . . . . . . . . . 10 2.1.1.4. External Magnetic Field . . . . . . . . . . . . . . . . 11 2.1.1.5. Stray Field . . . . . . . . . . . . . . . . . . . . . . . 11 2.1.2. Magnetic Structures . . . . . . . . . . . . . . . . . . . . . . 13 2.1.2.1. Domain Walls . . . . . . . . . . . . . . . . . . . . . 13 2.1.2.2. Skyrmions . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.3. Dynamics Described by the Landau-Lifshitz-Gilbert Equation 16 2.1.4. Spin Hall Effect . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.1.5. Spin-Transfer and Spin-Orbit Torques . . . . . . . . . . . . 18 2.2. Molecular and Brownian Dynamics Simulations . . . . . . . . . . . 19 2.2.1. Thermal White Noise . . . . . . . . . . . . . . . . . . . . . 20 2.2.2. Two-Body Interactions and Neighbor Lists . . . . . . . . . . 21 2.2.3. Periodic Boundary Conditions . . . . . . . . . . . . . . . . . 23 2.3. Skyrmion Dynamics Using the Thiele Equation . . . . . . . . . . . 24 2.3.1. Skyrmion Hall Angle . . . . . . . . . . . . . . . . . . . . . . 25 2.3.2. Dissipation Tensor . . . . . . . . . . . . . . . . . . . . . . . 25 2.3.3. Gyrocoupling Constant . . . . . . . . . . . . . . . . . . . . 26 2.3.4. Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3.5. Mean Squared Displacement (MSD) and Diffusion Constant 28 2.4. Skyrmion Pinning . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5. 2D Ordering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5.1. Phases in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.5.2. Radial Distribution Function . . . . . . . . . . . . . . . . . 31 2.5.2.1. Boundary Distribution Function . . . . . . . . . . . 33 2.5.2.2. Two-Dimensional Radial Distribution Function . . . 33 2.5.3. Voronoi Tesselation . . . . . . . . . . . . . . . . . . . . . . . 34 2.5.4. Hexatic Order Parameter . . . . . . . . . . . . . . . . . . . 35 2.5.5. Translational Order . . . . . . . . . . . . . . . . . . . . . . 37 2.6. Experimental Methods . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.6.1. Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 v Contents 2.6.2. Magneto-Optic Kerr Effect (MOKE) Microscopy . . . . . . 38 2.6.3. Skyrmion Tracking with Trackpy . . . . . . . . . . . . . . . 39 2.6.4. Linking of Skyrmion Trajectories . . . . . . . . . . . . . . . 40 3. Skyrmion Interaction Potentials 41 3.1. Derivation of Interaction Potentials via Iterative Boltzmann Inversion 41 3.1.1. Experimental Setup and Data . . . . . . . . . . . . . . . . . 41 3.1.2. Computer Simulation . . . . . . . . . . . . . . . . . . . . . 44 3.1.3. Iterative Boltzmann Inversion . . . . . . . . . . . . . . . . . 45 3.1.4. Analysis of Results . . . . . . . . . . . . . . . . . . . . . . . 47 3.1.5. Potentials for Larger Skyrmions . . . . . . . . . . . . . . . . 50 3.1.6. Applicability to Other Systems . . . . . . . . . . . . . . . . 51 3.2. Pinning of Skyrmions at Material Defects . . . . . . . . . . . . . . 51 3.2.1. IBI with Pinning Effects . . . . . . . . . . . . . . . . . . . . 54 3.2.2. Modeling of Ideal Pinning Sites . . . . . . . . . . . . . . . . 55 4. Dynamics and Phase Behavior of Skyrmion Lattices 57 4.1. Pinned Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.2. Skyrmion Hall Angle in Different Regimes . . . . . . . . . . . . . . 59 4.3. Influence of the Magnus Force on Skyrmion Diffusion . . . . . . . . 61 4.4. Ordering of Lattices . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.5. Analysis of Experimental Lattices . . . . . . . . . . . . . . . . . . . 66 5. Increasing Skyrmion Diffusion by External Stimuli 69 5.1. Oscillating Magnetic Field . . . . . . . . . . . . . . . . . . . . . . . 69 5.1.1. Experimental Results . . . . . . . . . . . . . . . . . . . . . . 69 5.1.2. Theoretical Calculations . . . . . . . . . . . . . . . . . . . . 72 5.1.3. Effective Reduction of Pinning . . . . . . . . . . . . . . . . 75 5.2. Alternating Current . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.3. Comparison of the Methods to Increase Diffusion . . . . . . . . . . 84 6. Skyrmion Dynamics in Confinements 86 6.1. Confinements Without Excitations . . . . . . . . . . . . . . . . . . 86 6.1.1. Circular Confinement . . . . . . . . . . . . . . . . . . . . . 87 6.1.2. Triangular Confinement . . . . . . . . . . . . . . . . . . . . 90 6.1.3. Square Confinement . . . . . . . . . . . . . . . . . . . . . . 92 6.2. Confined Dynamics with Applied Currents . . . . . . . . . . . . . . 94 6.2.1. Single Skyrmion . . . . . . . . . . . . . . . . . . . . . . . . 94 6.2.2. Usage of Multiple Skyrmions . . . . . . . . . . . . . . . . . 97 vi Contents 7. Conclusion and Outlook 99 7.1. Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 7.2. Limitations of the Particle-Based Model for Skyrmions . . . . . . . 100 7.3. Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 A. Appendix 104 A.1. Additional Functions . . . . . . . . . . . . . . . . . . . . . . . . . . 104 A.1.1. Running Average . . . . . . . . . . . . . . . . . . . . . . . . 104 A.1.2. Standard Error of the Mean . . . . . . . . . . . . . . . . . . 104 A.2. Technical Explanations . . . . . . . . . . . . . . . . . . . . . . . . . 104 A.2.1. Determination of the Material Boundary Position . . . . . . 104 A.2.2. Derivation of the Size Dependence of Skyrmion Diffusion . . 106 A.3. Additional Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 A.4. Polymer Physics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 B. Bibliography 132 vii 1. Introduction 1. Introduction Magnetic skyrmions [1–9] are topologically non-trivial nano- to micrometer sized chiral spin structures stabilized typically by the Dzyaloshinskii-Moriya interaction (DMI) [10–13], which are found in thin film multilayer systems [14–18]. Due to their stability [2, 19, 20], soliton property, and approximately 2D nature as their in-plane diameter is usually much larger than the thickness of the multilayer stack, skyrmions can be described as two-dimensional quasi-particles [21] exhibiting ther- mal Brownian motion [22–24]. Given their quasi-particle property and the easy manipulation of skyrmions by spin-transfer [25–27] and spin-orbit torques [28, 29], they are of interest for both fundamental research and novel applications. In fundamental research this interest is fueled in particular by the 2D phase behav- ior, which is distinctly different from the 3D case due to the possible occurrence of an hexatic phase. The hexatic phase is found in between the liquid and solid phase [30, 31]. This leads to a two-stage melting scenario from solid to liquid with partially unknown phase transitions, which have been theoretically described by the KTHNY theory [30–34]. The development of this theory was awarded with the 2016 Nobel Prize in Physics. Subsequent work has led to additional insight into this phase transition, including the discovery of systems in which the liquid- hexatic transition is first order in contrast to a continuous transition predicted by the KTHNY theory [35–38]. While experiments on these phase transitions have already been presented for a highly optimized system of colloids [39], the targeted manipulation of magnetic skyrmions is much easier as one can e.g. change their size by changing an external magnetic field [40] or employ different mechanisms to enhance their diffusion at constant temperature [41–43]. Therefore, high-density states of skyrmions have been subject of extensive recent research [44–50]. As the study of skyrmion lattices requires the analysis of large numbers of skyrmions, simulating lattices of micrometer-sized skyrmions involves systems sized hundreds of micrometers in each lateral dimension. For such systems, conventional approaches such as micromagnetic or atomistic simulations are not feasible due to their immense computational cost. This leads to the use of particle-based skyrmion models [21, 51–57] that exploit the quasi-particle properties of skyrmions enabling a description based on the Thiele equation [58]. The use of these models to describe skyrmion systems requires knowledge of the respective interaction potentials of skyrmions with other skyrmions as well as with the boundaries of the magnetic material. These interactions have been a subject of research in theory and simulations [21, 55, 59–61] resulting in analytic formulae [21, 59] as well as fits of these formulae for specific systems [54, 55, 60, 61]. A common disadvantage of these approaches is that they do not use experimental data for the determination of the potential shape 1 1. Introduction but purely rely on theoretical considerations. In addition, concrete micromagnetic interaction potentials were only determined for nanometer-sized skyrmions [54, 55, 60, 61] and are therefore mostly unsuitable for commonly studied micrometer-sized skyrmions. In this thesis, the above-mentioned limitations are addressed by ascertaining skyrmion-skyrmion and skyrmion-boundary interaction potentials directly from ex- periments. This is possible by exploiting the quasi-particle properties of skyrmions, which enables the application of methods known from soft matter physics. Such methods, which do not require prior knowledge of the shape of the interaction potentials, include iterative Boltzmann inversion [62, 63], which is used here. While this method is typically applied to coarse-grain atomistic potentials in molecular systems, its use for the analysis of experimental distributions has already been discussed [64]. One key advantage of this method is that it uses experimental data directly and can consequently predict experimental behavior more reliably than purely theoretical models. In addition to the simulation of extended skyrmion lattices, the same potentials and techniques can also be applied to describe the behavior of skyrmions in confinement. While tight confinement leads to commensurate ordering of skyrmions [54, 65], the shape of the confinements can also strongly affect the diffusive dynamics [66]. This thesis includes analysis of experiments and simulations on skyrmions in con- finements with different shapes, studying the connection between geometrically commensurate ordering and diffusion. Beyond skyrmion-skyrmion interactions, skyrmions in real experimental systems are affected by pinning effects [57, 67]. These effects limit skyrmion motion by attractively pinning skyrmions to specific positions in the sample. Pinning can be caused by several types of material inhomogeneities in the magnetic sample hosting skyrmions such as variations in different magnetic materials constants [68–71]. While pinning effects are usually mitigated by strong current-induced motion [56], they have a large impact on thermal diffusive motion [57], motion under weak applied currents [72–75] and skyrmion lattice formation [44, 57]. Pinning energy scales in systems exhibiting skyrmion diffusion are usually comparable to the ther- mal energy as the systems only show significant thermal dynamics when thermal fluctuations can overcome pinning effects [22, 43]. Furthermore, recent research has shown that pinning effects are strongly dependent on the skyrmion size [67, 76] and are therefore not universal for all skyrmions in a specific sample. In this thesis, different skyrmion systems affected by pinning effects are simu- lated. This includes different regimes of pinned skyrmion motion depending on the strength of the applied forces [74, 77–82] as well as a quantitative analysis of pinning-impeded diffusion. Furthermore, by applying knowledge of the properties of pinning interactions [67], this thesis presents two different methods to mitigate 2 1. Introduction pinning effects and enhance skyrmion diffusion achieving a type of motion close to free diffusion. Besides skyrmion lattices and interactions, another reason for a fundamental inter- est in skyrmions is the gyroscopic nature of their motion. Due to the skyrmion Hall effect [14, 58, 72, 73, 83], when a force is applied to a skyrmion, it moves with a certain angle respective to the applied force. This effect is caused by the chiral internal structure of the skyrmions and opens up new dynamic effects known as odd diffusion [84–86]. In this context, this thesis investigates how pinning affects this type of motion and how increasing the density of a skyrmion system can lead to an increased diffusion of skyrmions also in bulk systems without commensurability effects. With the knowledge of the motion and interactions of skyrmions, several novel computing and memory concepts can be realized. A popular concept of using skyrmions in the field of computing is the skyrmion racetrack, which is a memory device where information is encoded by the presence or absence of a skyrmion [87– 90]. A major advantage of this kind of memory is the very low current required to move the skyrmions through the racetrack when writing and reading informa- tion [87–89]. Here, the skyrmions do not perform thermal motion and therefore allow for a deterministic system as is required for current CMOS-based computing architectures. As this development led to a growing interest in magnetic skyrmions, additional concepts for skyrmion-based computing were developed [91]. In contrast to the racetrack memory, many of these proposed devices share the idea to exploit (room temperature) thermal Brownian motion of skyrmions, which was observed in specific thin film multilayer systems [22]. Thermal skyrmion dynamics has enabled a variety of non-deterministic computing concepts including a skyrmion reshuffler device for probabilistic computing [22, 92], random number generators [93], cellular automata [94] and token-based logic devices [43, 95, 96] such as a half-adder [43, 97, 98]. Another concept to harness thermal motion of skyrmions for computing is skyrmion-based reservoir computing [91, 99]. While this approach was initially proposed without the use of thermal diffusion [100–106], including thermal effects leads to a lower energy consumption [91, 99]. In reservoir computing, a magnetic skyrmion system acts as a part of a neural network to solve tasks such as audio recognition [104] or handwritten digit recognition [102]. An experimental realiza- tion of such a reservoir designed for performing logical operations by exploiting thermal diffusion is described in this thesis. The key advantage of these Brownian computing approaches, which exploit thermal random motion for computation, is their high energy efficiency as the random dynamics is fueled by the thermal energy of the surrounding [43, 99]. Skyrmions are a particularly promising system in this field due to their room-temperature thermal activity and electrical manipulability and readability [91, 107]. 3 1. Introduction The results presented in this thesis are organized into chapters describing published and unpublished projects. Chapter 2 describes the theoretical background necessary to perform and under- stand the computer simulations performed in this thesis. This includes insight into magnetism, computer simulations, analysis methods for two-dimensional systems, and a short part describing experimental methods required to obtain the data used here. The main part of the thesis contains both published and unpublished results and is presented in Chapters 3 to 6. Chapter 3 describes the quantitative determination of skyrmion interaction poten- tials both for the interaction with each other as well as with material boundaries and pinning sites. Most of the contents of this chapter are published in • Y. Ge, J. Rothörl, M. A. Brems, N. Kerber, R. Gruber, T. Dohi, M. Kläui, and P. Virnau, “Constructing coarse-grained skyrmion potentials from ex- perimental data with iterative Boltzmann inversion,” Commun. Phys. 6, 30 (2023). Chapter 4 presents various results on skyrmion dynamics and motion including pinned diffusion and skyrmion lattices. Predominantly, the results in these chapter are purely simulation results, but it also includes an analysis of experimental skyrmion lattices. Some of the results are published in • D. Schick, M. Weißenhofer, L. Rózsa, J. Rothörl, P. Virnau, and U. Nowak, “Two levels of topology in skyrmion lattice dynamics,” in press, 2023. Chapter 5 contains results on the enhancement of skyrmion diffusion by applying external stimuli. These include periodically alternating magnetic fields and electric currents. The part on magnetic fields is already published while results on alter- nating currents are subject of ongoing research. The first half of this chapter is published in • R. Gruber, M. A. Brems, J. Rothörl, T. Sparmann, M. Schmitt, I. Kononenko, F. Kammerbauer, M.-A. Syskaki, O. Farago, P. Virnau, and M. Kläui, “300- times-increased diffusive skyrmion dynamics and effective pinning reduction by periodic field excitation,” Adv. Mater. 35, 2208922 (2023). Chapter 6 describes skyrmions in small confinements where commensurability effects are key. Further parts of this chapter describe simulations performed to understand viability of systems for skyrmion-based computing. Parts of these results are published. The major part of this chapter is published. Relevant publications for this chapter are: 4 1. Introduction • C. Song, N. Kerber, J. Rothörl, Y. Ge, K. Raab, B. Seng, M. A. Brems, F. Dittrich, R. M. Reeve, J. Wang, Q. Liu, P. Virnau, and M. Kläui, “Commen- surability between element symmetry and the number of skyrmions governing skyrmion diffusion in confined geometries,” Adv. Funct. Mater. 31, 2010739 (2021), • K. Raab, M. A. Brems, G. Beneke, T. Dohi, J. Rothörl, F. Kammerbauer, J. H. Mentink, and M. Kläui, “Brownian reservoir computing realized using geometrically confined skyrmion dynamics,” Nat. Commun. 13, 6982 (2022), • T. B. Winkler, J. Rothörl, M. A. Brems, H. Fangohr, and M. Kläui, “Coarse- graining collective skyrmion dynamics in confined geometries,” 10.48550/ arXiv.2303.16472 (2023). Chapter 7 discusses the results presented in this thesis including an outlook on other interesting projects which could be performed using the Thiele model. Additionally, limitations of the model and methods used in this thesis are discussed. The appendix contains additional information relevant to understand the thesis including some formulae, a mathematical derivation of skyrmion diffusion with an alternating external magnetic field and additional figures. Beyond this, the appendix also discusses projects on which I worked during the time of my PhD studies which are not related to magnetic skyrmions. The relevant publications in this field are • J. Rothörl, S. Wettermann, P. Virnau, and A. Bhattacharya, “Knot formation of dsDNA pushed inside a nanochannel,” Sci. Rep. 12, 5342 (2022), • Y. Zhao, J. Rothörl, P. Besenius, P. Virnau, and K. C. Daoulas, “Can polymer helicity affect topological chirality of polymer knots?” ACS Macro Lett. 12, 234–240 (2023). 5 2. Theoretical Background 2. Theoretical Background 2.1. Magnetism As this thesis studies magnetic skyrmions, I start with an introduction on the mag- netic interactions, which occur in metallic thin-film systems hosting skyrmions. Af- terwards I continue with an explanation of magnetic structures including skyrmions as well as a description of how magnetic systems evolve dynamically. 2.1.1. Magnetic Energy Terms A magnetic system consists of electrically charged particles with an angular momen- tum such as protons and electrons, which have a magnetic moment ®̀ originating from their spin as well as their orbital angular momentum. Summing up the magnetic moments in a volume and dividing by the volume gives a magnetization "® (A®) which has an absolute value between 0 and the saturation magnetization "S, which is the magnetization the system reaches in case all magnetic moments are aligned in parallel. As describing every individual magnetic moment is not feasible for a system of a size of up to hundreds of micrometers, one can describe it using the micromagnetic approximation [114]. In this approximation, the interactions of individual spins are coarse-grained to obtain a continuous field describing the local magnetization direction <® (A®). In this approximation, the Hamiltonian of a simple macroscopic magnetic system can be calculated by adding up the following energy terms [115, 116] which are explained in the paragraphs below  = Exchange + DMI + Anisotropy + Zeeman + Stray. (2.1) While there are further terms which could affect the energy of a magnetic system such as external stress or magnetostriction [116], I limit the description to the terms mentioned above as these are the most relevant ones for describing material systems hosting skyrmions. 2.1.1.1. Exchange Interaction The dominant interaction in many magnetic systems is the exchange interaction. It is a purely quantum-mechanical effect [117] which does not occur in classical systems but which can be nonetheless described in a classical manner. To describe this effect microscopically, one needs to model the spatial wave functions of two 6 2. Theoretical Background electrons sharing an orbital. This total wave function Ψ can be either symmetric ΨS(A®1, A®2) 1 = √ [Φa(A®1)Φb(A®2) +Φb(A®1)Φa(A®2)] (2.2) 2 or antisymmetric ΨA(A®1, A® 1 2) = √ [Φa(A®1)Φb(A®2) −Φb(A®1)Φa(A®2)] (2.3) 2 under spatial permutation with Φa and Φb being orthonormal atomistic wave functions. Due to the fermionic nature of electrons, they are suspect to the Pauli exclusion principle. The total wave function of a two-electron orbital has to be antisymmetric, allowing for either a symmetric spatial and an antisymmetric spin wave function or vice versa. This constraint to the spatial wave functions depending on the spins yields a different energy depending on the relative sign of the two spins. The total spin of two electrons that each have a spin quantum number B = 12 and a squared spin angular momentum 〈B®2〉 2= ℏ2B(B + 1) = 3ℏ 8 4 is [118] 3ℏ2〈(®2〉 = 〈(B®1 + B® 22) 〉 = 〈B®21〉 + 〈B®22〉 + 2〈B®8 · B®9 〉 = 2 + 2〈B®8 · B®9 〉. (2.4) The total spin quantum number ( is either 0 or 1 leading to 〈(®2〉 = ℏ2((( + 1) being either 0 or 2ℏ2 and therefore 〈B®8 · B®9 〉 ℏ 2 = 4 for parallel spins (triplet state, ( = 1) and 〈B®8 · B®9 〉 = −3ℏ 2 4 for antiparallel spins (singlet state, ( = 0). The energy per spin is different for the triplet state (with energy t)) and for the singlet state (with energy s) and can be written as [118] 1  = (s + 3t4 ) + ( t − s) 12 〈B®8 · B®9 〉, (2.5)ℏ As the first term is constant, it can be neglected as constant terms in a Hamiltonian do not affect the physical system. Therefo∑re, one can write [118]# Exchange = − 8 9 B®8 · B®9 , (2.6) 8> 9 where the sum goes over all pairs of particles in the system. Here, 8 9 is the exchange constant for a specific pair of spins and B8 the spin of particle 8. To prevent overcounting, one only needs to address each pair once, as indicated by the 8 > 9 in the sum. In a simple model of spins on a lattice one can assume 8 9 =  for next neighbors 7 2. Theoretical Background (a) (b) Figure 2.1.: Illustration of a ferromagnetic (a) and antiferromagnetic (b) state of a 2D lattice of spins. These states are ground states of the lattice at zero temperature. Colors indicate the spin directions. and 8 9 = 0 else. For a 2D system with spins placed on a square lattice, this leads to four neighbors per spin in an infinite system. In this case  > 0 leads to ferromagnetic ordering with all spins being aligned in the same direction in the ground state as shown in Figure 2.1(a), while for  < 0 the ground state is antiferromagnetic as in Figure 2.1(b). Situations where 8 9 depends on 8 and 9 can also lead to e.g. a ferrimagnetic or a paramagnetic ground state. To obtain a continuum approximation of this equation, assuming the above simpli- fication one can consider particles 8 and 9 to be placed at positions (A®8 − X4̂U) and (A®8 + X4̂U) with X being∑sma# ∑ll: Exchange ∼ <® (A®8 − X4̂U) · <® (A®8 + X4̂U) (2.7) ∑8 U=∑{G,H}# ≈ (<® (A®8) − XmU<® 8 (A8)) · (<® (A®8) + XmU<® 8 (A8)) (2.8) ∑8 U=∑{G,H}# = 2<® (A® )28 − (XmU<® 8 (A®8)) . (2.9) 8 U={G,H} The same calculation can also be performed analogously for systems with three dimensions by summation over U = {G, H, I}. Due to the total magnetization 8 2. Theoretical Background ∑# 8 <® (A®8)2 being constant,∫one can neglect the∫first term∑and arrive at[116] Exchange = 3 3 2A (∇<®)2 = 3 32A mU 0, thermal fluctuations of positions of individual particles are observed. The strength of the noise in a system with given temperature and damping is given by the fluctuation-dissipation theorem. Thermal white noise 5® has to obey [146, 151] 〈 58 (C)〉 = 0, (2.56) 〈 58 (C) 5 ′ ′9 (C )〉 = 2:B)WX8 9X(C − C ), (2.57) with :B being the Boltzmann constant, X8 9 the Kronecker delta for spatial indices, X(C − C′) the delta function for times C and C′ and angle brackets indicating the average. The noise acting on different particles as well as the noise acting on the same particle in different points in time is independent. Because damping and white noise are changing the energy of the system, it is able to relax into energy minima. Due to the noise obeying the fluctuation-dissipation theorem with given 20 2. Theoretical Background damping W, the average temperature 〈)〉 of the system is kept constant over time achieving a system in the canonical ensemble. Given the stochastic nature of the noise, the temperature ) at a specific time step is however still a stochastically fluctuating quantity. To add this kind of noise to a discretized computation scheme at a specific time step, one needs to consider the time step ΔC [148] 〈 58 5 9 〉 2:B)W= X8 9 . (2.58)ΔC While typically, the random force applied to a particle is Gaussian distributed, implementing this would take a large computational effort. To save this effort, instead a uniformly distributed random noise, which fulfills Equations 2.56 and 2.58, suffices to describe the motion of Brownian particles and gain consistent results [152]. This requires the same mean and variance as the Gaussian noise would have. With the minimum and maximum of a uniform distribution being 0 and 1, mean and variance are given by 1 ` = 2 (0 + 1) (2.59) 2 1f = 12 (1 − 0) 2, (2.60) respectively. From the mean being ` = 0 follows clearly that 0 = −1 and th√erefore f2 = 1 2 √ 31 . This leads to boundaries of the uniform distribution of ± 3 = ± 6:B)Wf ΔC for each coordinate. 2.2.2. Two-Body Interactions and Neighbor Lists Molecular Dynamics simulations usually contain interactions between different particles. The most simple kind of these interactions are two-body interactions where one particle interacts with another leading to a force applied to one or both interaction partners. In a system with # particles, this means that there are up to # (# − 1) sets of interactions to be computed. For 1000 particles, this would require nearly a million calculations to get the interaction force which would be the most time-consuming part of most MD simulations. For symmetrical interactions between two particles, this number could be halved by just considering the interaction of 8 with 9 and not 9 with 8, but the calculation of interactions would still have a computation time of O(#2) [153]. This problem can be handled by using a cutoff distance Acut on the two-body interactions beyond which no interaction force is calculated. This is usually justified by the decay of most interactions with distance and therefore their negligible impact 21 2. Theoretical Background rcut Figure 2.7.: Graphical example for a cell neighbor list. The simulated area is separated into squares of side length Acut. A particle depicted by the red dot interacts with all other particles in a circle of radius Acut around it. at large spatial separations. This works for e.g. Lennard-Jones interactions or exponentially decaying interactions while for Coulomb interactions, which are not considered in this thesis, this simplification would not work. When just using a cutoff, one must of course still calculate distances between all pairs of particles and therefore only saves the computation of interaction forces. To only consider the relevant interactions one can use a neighbor list instead of calculating all individual distances. The most straightforward approach of constructing such a neighbor list is a cell neighbor list which is used in most simulations performed for this thesis (simulations using HOOMD [148] use a more elaborate neighbor list while simulations with only a few skyrmions do not use a neighbor list for efficiency reasons). Using this kind of neighbor list, the system is separated in quadratic non-overlapping cells with a side length of Acut as depicted in Figure 2.7. Every particle in the system is now assigned to one of the boxes. Due to the cutoff, a particle can only interact with other particles in its own cell or in one of the eight neighboring cells (in a 2D system) instead of the whole simulated area. With a constant density of particles, this leads to a constant number of interactions per particle and therefore the system scales as O(#) instead of O(#2) allowing for much larger simulations [153]. 22 2. Theoretical Background Figure 2.8.: Graphical example for periodic boundary conditions. The particles inside the center box are mirrored into neighboring boxes and interacting with particles in the center box. The circle around the red particle contains all interaction partners of this particle including the gray particle in the left box which is a periodic copy of the original gray particle in the central box. 2.2.3. Periodic Boundary Conditions When implementing a simulation system, one needs to define the behavior of particles reaching the boundary of the simulation box. In this thesis, two different ways of dealing with this are used: A repulsive boundary and periodic boundary conditions. A repulsive boundary applies a repulsive force to each particle which is large enough that particles cannot reach the boundary of the simulation box. The other way of implementing a system limit is using periodic boundary conditions. Here, particles can freely cross the system’s boundary leading to them behaving as if they were in an infinite system. To achieve this, particles outside the initial cell are projected back into the initial cell and interact with other particles as if they were in this position. This interaction is shown schematically in Figure 2.8. When calculating interactions, the projected position is calculated for each coordinate via [153] A′U = AU − bAU/3Uc · 3U, (2.61) 23 2. Theoretical Background with A′U being an element of the projected position, AU an element of the real position, b..c indicating the floor operation and 3U being the size of the system in U-direction. The projected positions are all part of the initial simulation box. For consistency, this also requires particles e.g. at the right border of the box to interact with particles at the left border. This is achieved by modifying the relative distance A®8 9 of particles to each other with a correction similar to Equation 2.61 [153] A′8 9 ,U = A8 9 ,U − round(A8 9 ,U/3U) · 3U, (2.62) where round indicates rounding to the nearest integer number. To prevent double counting of interactions, this type of boundary conditions is however only viable if there is a cutoff to interactions and if the initial cell dimensions are at least twice the cutoff radius: 38 ≥ 2Acut. Additionally, when using a cell list as described in the section above, this list also needs to take into account periodic boundary conditions for determining neighboring cells. In case the dimension of the system is not a multiple of Acut and has a periodic boundary, one also needs to account for this by using smaller boxes at the boundary of the system. With this, every box except for the ones touching the boundary is normally sized and interactions between two normally sized boxes can be calculated as described above. If a particle is located in a box neighboring one of the smaller boxes, the box behind the neighboring box is then also considered for possible neighbors. With this, one can again guarantee that all particle pairs within a distance of Acut are considered in the calculations. 2.3. Skyrmion Dynamics Using the Thiele Equation Due to their high stability, circular shape and soliton property a skyrmion can be approximated as a point-like quasi-particle in many situations. This allows for using methods from standard Molecular Dynamics simulations to simulate the behavior of magnetic skyrmions. These simulations are performed using the Thiele equation which was developed explicitly for∑the movement of magnetic quasi-particles [21,58] WE®8 − Î × E® = ®8 SkSk(A®8 − A®9 ) + ® ®4GC + 5 , (2.63) 9 where W is the damping constant, Î×E®8 is the Magnus force with  the gyrocoupling constant, ®SkSk the skyrmion-skyrmion interaction force, ®4GC contains external forces such as interactions with material defects and sample boundaries or forces due to applied torques and 5® is thermal white noise following the fluctuation-dissipation theorem. Note that the Magnus force acts perpendicular to the skyrmion velocity and therefore does not contribute to the energy of the system. For this reason, it 24 2. Theoretical Background affects only dynamics but not statics [52]. To perform computer simulations using this equation, one needs to know the values of W and  as well as functional forms for the forces applied to the skyrmions. The Thiele equation presented here does not contain an inertia term and therefore uses Brownian dynamics as in Equation 2.53. This is a typical model used to represent skyrmions [21] as inertia effects only play a role for frequencies at least of the same order as the intrinsic frequencies of the magnetic system which are in the gigahertz regime [51] while the frequencies of excitations used here do not exceed 10 kHz. 2.3.1. Skyrmion Hall Angle Due to the gyroscopic part of the Thiele equation, skyrmions do not move parallel to an applied force but instead with some angle relative to the force. This angle is called the skyrmion Hall angle and given by ([58)] \Sk = arctan  . (2.64) W A Hall angle close to zero indicates movement approximately parallel to an applied force while an angle close to 90°or -90° indicates very strong gyroscopic motion and therefore motion distinctly different from e.g. colloids. 2.3.2. Dissipation Tensor The damping constant used in simulations W is defined as W = UT (2.65) with the Gilbert damping U and 2T the trace of the dissipation tensor 8 9 . The dissipation tensor is defined by an inte∫gral over the magnetic surface [58] "S 8 9 = 2 ∗ 3 Am8<® · m9<® , (2.66)3W with "S the saturation magnetization, 3 the thickness of the magnetic layer. For isotropic skyrmions the dissipation tensor is also isotropic such that T = GG = HH. In most cases the trace of the dissipation tensor is proportional to the skyrmion radius. 25 2. Theoretical Background 2.3.3. Gyrocoupling Constant The gyrocoupling tensor  is also∫calculated( via an in)tegral over the magneticsurface as [58] "S  28 9 = ∗ 3 A<® · m8<® × m9<® (2.67)3W and has zero arguments in the main diagonal due to the symmetry of the cross product. Assuming skyrmions being cylindrically symmetric along the I-axis this can be simplified to ¬E® = Î × E® with  = GH. For a skyrmion this can be further simplified to " = ( ∗4c= (2.68)3W due to the radial symmetry of the skyrmion. = is the winding number defined in Equation 2.39 which is equal to one for the skyrmions discussed in this thesis. This result does not depend on the size of a skyrmion for given material properties. With the gyrocoupling being constant and the dissipation tensor usually being larger for larger skyrmions, a strong gyroscopic motion is especially found for small skyrmions. 2.3.4. Diffusion Knowing the damping and gyrocoupling constants, one can calculate the average thermal diffusion of a skyrmion. Here, I derive the well-known result for the diffusion constant of skyrmions [154]. To do so, only a random force 5® is considered simplifying Equation 2.63 to WE®8 − Î × E®8 = 5®. (2.69) In components one obtains WEG + EH = 5G (2.70) WEH − EG = 5H . (2.71) Solving these equations for EG and EH gives: W 5G −  5H EG = (2.72) W2 + 2 W 5H +  5G EH = . (2.73) W2 + 2 26 2. Theoretical Background This can be used to calculate the v〈elocity aut〉ocorrela2 [t〈ion 〉 〈 〉] ′ W 〈 5G (C) 5G (C ′)〉 + 2 5H (C) 5H (C′) − W 5 (C) 5 ′ ′〈 ( ) ( )〉 G H (C ) + 5H (C) 5G (C ) EG C EG C = ( 2 + 2)2 2 〈 ′ 〉 W  [〈 〉 〈 (2.7〉4]) ′ W 5H (C) 5H (C ) +  2 〈 5G (C) 5 ′G (C )〉 + W 5G (C) 5H (C′) + 5H (C) 5G (C′)〈EH (C)EH (C )〉 = 2 〈 〉 〈 (W2〉+ 2) [ 〈 (2.75) W2 〉] ′ 5G (C) 5 ′ H (C ) − 2 5 ′G (C ) 5H (C) + W 〈 5G (C) 5 (C′G )〉 − 5H (C) 5 ′H (C )〈EG (C)EH (C )〉 = . ( 2W2 + 2) (2.76) U〈 sing〉 the properties of white noise given in Equations 2.56 and 2.58 and especially 58 5 9 ∼ X8 9X(C − C′), this simplifies to: W2 · 2: )W + 2〈 ( ) ( ′)〉 B · 2:= B)W ( − ′) 2:= B)WEG C EG C 2 X C C 2 + 2 X(C − C ′) (2.77) (W2 + 2) W  W2 · 2: )W + 2〈 ( ) ( ′)〉 B · 2:B)W ( − ′) 2:B)WE ′H C EH C = 2 X C C = 2 + 2 X(C − C ) (2.78)(W2 + 2) W  〈 ( ) ( ′)〉 W · [2:B)W − 2:B)W]EG C EH C = 2 X(C − C ′) = 0, (2.79) (W2 + 2) which results in the same average velocity square in both spatial directions 〈E ∫(C)E (C′G G ∫)〉 = 〈EH (C)E (C′H )〉. Inte∫gratin∫g this autocorrelation in time givesC C C C ′〈 ( ) ( ′)〉 〈∫ ∫′ 2:B)W ( 〉− ′) 2:B)W3g 3g E8 g E8 g = 3g 3g X g g = C (2.80)0 0 0 0 W2 + 2 W2 + 2C C = 3gE (g) 3gE ′ 28 8 (g ) = 〈ΔG8 (C)〉, (2.81) 0 0 where ΔG8 (C) is the displacement of a particle in direction 8 over a time C. This can be identified with the diffusion constant defined by [154] 〈ΔG2(C)〉 = 8 2 , (2.82)C 27 2. Theoretical Background which is calculated in simulations via the mean squared displacement (MSD) explained below. This leads to [154] W  = :B) 2 + 2 , (2.83)W  which shows that the gyroscopic force reduces diffusion compared to a system which does not contain a gyroscopic component. In the case of  = 0, this diffusion equation simplifies to the Einstein-Smoluchowski equation which is the standard form for dissipative motion of Brownian particles [145, 155] 1  = :B) (2.84) W 2.3.5. Mean Squared Displacement (MSD) and Diffusion Constant To calculate the diffusion constant from a trajectory, the mean squared displacement (MSD) is used. It is given by [153] 〈 〉 MSD(g, C) = ( 2A®(C + g) − A®(C)) (2.85) # with A®(C) the position of a particle at a time C and the average taken over all particles in the system. Usually, this average over all particles is combined with a running average over different times C with the lag time g being kept constant. The lag time is the time difference between the two frames where skyrmion positions are compared [153]: 〈 〉 MSD(g) = (A®(C + 2g) − A®(C)) . (2.86) #,C This can be done due to equilibrated systems being Markovian and therefore their dynamics depending only on time-differences and not on the absolute value of the time C. When describing equilibration processes, this method is not appropriate and the formula in Equation 2.85 has to be used instead. The MSD of free Brownian diffusion is connected to the diffusion constant  of a system via [153] MSD(g) = 4:B)g, (2.87) for a 2D system with the Boltzmann constant :B and temperature ) . 28 2. Theoretical Background 2.4. Skyrmion Pinning While the above sections describe models for perfectly homogeneous materials and skyrmions, in reality properties of magnetic materials are inhomogeneous and therefore skyrmion behavior does not exactly follow the equations given above. In particular, material inhomogeneities usually lead to a dependence of the energy of a skyrmion on its position. For example, a spatially dependent exchange constant would lead to a spatial dependence of the energy of a domain wall. Due to such effects, skyrmions are more likely to be located in energy-minimizing regions, i.e. pinning sites, and therefore do not isotropically explore the whole experimental sample. This leads to changes in static ordering of skyrmions [110] as well as their dynamics [22, 57, 72] compared to the case that does not include pinning effects. While in systems with strong applied currents pinning effects are often mitigated leading to an increased diffusion similar to the unpinned case [56, 72], pinning effects are found to be of a similar magnitude as thermal effects in many cases and therefore need to be considered when describing applications using skyrmion- based devices [67]. The resulting mode of diffusion with such pinning effects is a hopping-like diffusion in most cases where a skyrmion moves between different pinning sites [41]. To further understand the dynamics of pinned skyrmions, one needs to investigate their internal structure beyond the simplifications of the Thiele model. While the skyrmion core as well as the area outside a skyrmion is homogeneously magnetized, the domain wall at the skyrmion boundary contains non-parallel spins and therefore interacts differently with changes in material parameters. Colleagues have shown that for the kind of skyrmions used in the experiments described in this thesis, most pinning effects apply to the domain wall and not to the skyrmion core and that therefore the pinning affects different parts of the skyrmion inhomogeneously [67] leading to e.g. a strong size dependence of pinning effects [41]. 2.5. 2D Ordering To understand the lattice ordering of two-dimensional particles such as skyrmions, one must consider both general properties of lattices and phase transitions as well as specific properties arising from the dimensionality. 2.5.1. Phases in 2D The phase behavior of two-dimensional systems is distinctly different from 3D systems. In three dimensions, a first-order phase transition between the liquid and the solid phase occurs [156]. In two dimensions, depending on the shape of 29 2. Theoretical Background interactions [36, 157], an additional phase is observed between these two phases [30, 31]. The properties of phase transitions between the three phases are still not fully understood and subject of theoretical [35, 36, 38, 158] and experimental research [39, 44, 47, 159, 160]. Liquid The (isotropic) liquid phase is a phase without any long-range order. Therefore, spatial order decays exponentially in this phase. Short-range interactions however still lead to some degree of short-range order as e.g. an excluded volume interaction prevents too close proximity of different particles. Hexatic The hexatic phase is a special property of 2D systems in which a decoupling of orientational and translational order can occur. It is also called anisotropic liquid due to the exponentially decaying translational ordering. On the other hand, the orientational order decays algebraically leading to a quasi-long range order of the local orientation of a lattice but not of the particle positions. Usually, this phase only occurs in a quite narrow range of thermodynamic parameters [36]. This can change in systems with strong quenched disorder such as disorder induced by pinning which causes a broader hexatic phase as pinning inhibits long range translational order [161–163]. Solid The solid phase is characterized by quasi-long range translational order decaying algebraically with distance. This differentiates 2D solids from 3D solids which may display a perfect long range order. This is caused by logarithmically diverging displacements in a 2D crystal at finite temperature as described in the Mermin-Wagner theorem. This theorem says that a spontaneous symmetry breaking, which would be required for a perfect long-range order, cannot occur [164– 167]. The orientational order in a solid is truly long range and does no longer decay algebraically as in the hexatic phase [36]. Phase Transitions While the KTHNY theory [30–34] proposes that two contin- uous phase transitions happen when transitioning from a crystallized solid to a liquid stat, other publications propose a first-order liquid-hexatic and an continuous hexatic-solid transition based on computer simulation results [35, 36, 158, 168]. To differentiate these two kinds of phase transitions, one has to check for the coexistence of domains with different ordering which would indicate a first-order transition. The solid to hexatic phase transition is considered to be first order both in the KTHNY theory and in later works reproducing the phase transitions in computer simulations [35, 36, 158]. 30 2. Theoretical Background 2.5.2. Radial Distribution Function The radial distribution (or pair correlation) function 6(A) describes the dimensionless density of particles at a certain distance to another particle. For a completely disordered system like an ideal gas the radial distribution function is 6(A) = 1, indicating that there is no preferred distance between different particles. Systems with particle-particle interactions display a different shape of the radial distribution with peaks indicating ideal distances and therefore some degree of ordering. In a perfect solid crystal at ) = 0 the distribution function therefore has distinct delta-peaks at ideal distances and the structure of the distribution function is conserved for long range. For particles with no long-range order (e.g. liquids) 6(A) converges to 1 for large values of A. Some examples of the radial distribution function for systems of repulsively interacting particles at finite temperature are given in Figure 3.6. In the general case of 3 dimensions, one can define the radial distribution function as [169] 〈∑ 〉+3 6(A®) = 2 X (3) (A® − A®8 9 ) (2.88) # 8≠ 9 where +3 is the 3-dimensional volume, # is the number of particles in this volume, A®8 9 is the vector connecting the positions of particles 8 and 9 and X(3) is the 3-dimensional Kronecker-delta function. Note that this definition of the radial distribution function only applies to periodic systems and does not account for closed boundaries of any kind. A method to handle boundaries is explained below. Also note that this function depends on the vector A® and thereby also explicitly on angles instead of just the radial distance A. To simplify the above given equation for 2D systems such as skyrmion systems described in this thesis, one can define the two-dimensional density d = #/+2. Additionally, the delta function can be split into individual coordinates and then transformed to polar coordinates〈 〉 ( 1 ∑ 6 G, H) = 〈 X(G − G8 9 )X(H − H8 9 ) (2.89)d# 1 ∑8≠ 9 〉 6(A, i) 1= | X(A − A8 9 )X(i − i8 9 ) , (2.90)d#  | 8≠ 9 31 2. Theoretical Background 4 (a) 4 (b) 4 (c) 1 (d) 3 3 3 2 2 2 1 1 1 0 0 0 0 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 r r r r Figure 2.9.: Radial distribution functions of a two-dimensional system of 40,000 repulsive particles (soft discs as described in Ref. [36] with a pairwise interaction potential + (A) = A−8) at different densities d with periodic boundary conditions: (a) Low-density liquid (d = 0.5). (b) High-density liquid (d = 1.0). (c) Solid at finite temperature (d = 1.3). (d) Solid at ) = 0 (d = 1.3). As the radial distribution function becomes a sum of delta distributions for a solid at zero temperature, 6(A) is given in arbitrary units for this case. where  is the Jacobian of the coordinate transformation which is  = A for 2D polar coordin∫ates. By integrati∫ng over the po〈l∑ar coordinate, one now〉arrives at [44]2c 2c 1 1 3i6(A, i) = 3〈i X〉(A − A8 9 )X(i − i8 9 ) (2.91)0 0 ∑d# A 8≠ 9 2 1 1c6(A) = d# A 8≠〈 X(A − A8 9 ) 〉 (2.92)9 1 1 ∑ 6(A) = 2 X(A − A8 9 ) , (2.93)d# cA 8≠ 9 which is a typical form used for 2D systems [44]. This function however still uses a delta-function which is the correct form for infinite statistics but does not describe binning of data which is required to get results from real systems with limited statistics. To introduce binning to the description, one can discretize the delta-function by using the equality to the Heaviside-step-function X(G) = 3 \ (G) 3G and discretizing the derivative 3 \ (A8 9 − (A + ΔA)) − \ (A8 9 − A) X(A − A8 9 ) = X(A8 9 − A) = − \ (A8 9 − A) = − . (2.94) 3A ΔA 32 g(r) g(r) g(r) g(r) [arbitrary units] 2. Theoretical Background Using this, one arrives at 1 1 ∑# [ ] 6(A) = \ (A8 9 − A) − \ (A8 9 − A − ΔA) (2.95) d# 2cAΔA 8≠ 9 where ΔA is the bin width of the histogram. The sum goes over all pairs of particles in the system. The first peak of the radial distribution function refers to the distance between nearest neighbors. Now, to apply this function to an experiment, one also needs to deal with the boundary conditions. To achieve this, instead of using the whole experimental picture to choose particles 8 and 9 with equal weight, one has to either limit the region where particles are chosen from or apply a correction. For simplicity, in this thesis, particles 8 are chosen from a region with distance to the system boundaries of at least the maximum distance at which 6(A) is computed. Particles 9 are then chosen from a larger region so that all possible neighbors of 8 particles are accounted for. Inserting this into Equation 2.95, # is the number of 8 particles while d can be calculated either from the inner or the outer part assuming the density is not dependent on the location. Here, the outer part was chosen to gain more statistics. 2.5.2.1. Boundary Distribution Function The logic behind the radial distribution function can be extended to the distribution of particles relative to the boundary of a sample containing 2D particles. For a homogeneous straight boundary, t∑he boundary distribution function is defined as ( ) 16Bnd A = [\ (A8 − A) − \ (A8 − A − ΔA)] , (2.96) d!ΔA 8 where ! is the length of the boundary with which the interaction is calculated and A8 is the distance of particle 8 perpendicular to the boundary. For large distances 6Bnd(A) converges to 1 analogously to 6(A). 2.5.2.2. Two-Dimensional Radial Distribution Function Expanding the concept of a radial distribution function to two dimensions, one can define a measure for the structu∑re of 2D lattices [35, 36]1 # [ ] 6(A) = [ ( )2 \ (ΔG8 9 − G) − \ (ΔG8]9 − G − ΔA)d# ΔA 8≠ 9 (2.97) · \ (ΔH8 9 − H) − \ (ΔH8 9 − H − ΔA) , 33 2. Theoretical Background Figure 2.10.: Example Voronoi tesselation of a region of a high-density two- dimensional liquid. Particles are indicated by black dots with their Voronoi cell colored according to their number of neighbors. Cells with six neighbors are drawn in gray, cyan indicates seven neighbors, purple five neighbors and blue eight neighbors. with ΔG8 9 being the distance between particles 8 and 9 in x-direction and ΔH8 9 being defined analogously. As the one-dimension√al function, the 2D radial distribution function 6(G, H) converges to one for A = G2 + H2 → ∞ in a structure without long-range order. In short distances, this function does not only contain favored distances as the one-dimensional 6(A) but also shows the lattice orientation with six peaks in the first layer of a hexatic structure. 2.5.3. Voronoi Tesselation To find neighbors of a particle in an arbitrary system, one needs a measure to define what a neighbor is. Here, neighbors are defined as particles with a common edge in a Voronoi tesselation or Voronoi diagram. A Voronoi diagram is created by adding all points in a plane which are closest to one specific particle 8 to a cell which also contains the respective particle 8 and no other particle. When using the standard Euclidean distance, this leads to polygons with straight edges where two neighboring polygons are connected by exactly one common edge [170]. Such a Voronoi tesselation is shown in Figure 2.10. Voronoi diagrams can also be used to analyze phases of particle ensembles because 34 2. Theoretical Background in a perfectly ordered mono-crystalline phase, each particle has exactly the same number of neighbors. When a solid order decays, different kinds of defects lead to different numbers of neighbors of individual particles in the Voronoi tesselation. In this thesis, Voronoi tessellations are computed using the freud library for python [170, 171]. Probability of Sixfold Orientation To quantify the order based on a Voronoi tesselation, one can count the number of particles whose neighbor number is the same as in the ideal lattice of some kind. For a hexagonal lattice, which is the ideal structure for skyrmions, this number is 6. Therefore, one can use the probability of sixfold orientation %6 as an indicator of lattice ordering. This parameter is especially relevant to quantify liquid-hexatic and hexatic-solid phase transitions where a strong increase in %6 is observed. 2.5.4. Hexatic Order Parameter Hexatic order in a 2D system can be specified by using the hexatic order parameter k6 which calculated for every individual p∑article 9 [36, 160, 172]1 = k6( 9) = 468\ 9: , (2.98) = :=1 where \ 9 : is the angle between the connecting vector A®9 : and a fixed axis in the system (usually the G-axis, choosing a different axis only changes the complex angle and not the absolute value of k6). The sum goes over all = neighbors of a particle found via a Voronoi tesselation. Particles with a number of neighbors different from six usually have a k6 close to zero while a perfect hexagonal order is indicated by |k6 | = 1 for every particle. For instance, this local parameter can be used to detect e.g. domain boundaries as particles at the boundary have different numbers of neighbors and therefore a low k6. A simple indicator for the total order of a system is the average absolute value of k6 [44] 〈| |〉 1 ∑# k6 = |k6( 9) |. (2.99) # 9=1 A value close to zero indicates either a disordered phase or a phase with an order strongly different from hexagonal while a value close to one indicates a hexagonally ordered solid. As this is just an average local order parameter, 〈|k6 |〉 does not reliably describe phase transitions as it does not account for long-ranged correlations. However, it is a valuable estimator for the degree of ordering in a small system or 35 2. Theoretical Background a system which is only partially accessible. Local Orientation Angle To describe local orientation in a 2D system, one can use the polar angle \6 of the hexatic order parameter k6. This angle describes the local orientation of a particle with respect to its neighbors. Due to the sixfold symmetry of this lattice, only angles between 0 and 60 degree indicate a different lattice so that \ ± c6 3 =̂ \6. In a perfect lattice, this angle is uniform while a disordered state shows randomly distributed angles. The local orientational angle can be used to visualize domains of different orientations. Angular Susceptibility Besides the average absolute value of k6, one can also calculate the absolute average of the hexatic order parameter which is known as angular susceptibility [173] 1 ∑# j6 = k6( ) 9 . (2.100)# 9=1 As long as the lattice is not oriented homogeneously, this quantity is going to be very close to zero. This changes when the system transitions to hexatic order and long-range orientational order becomes present. In the solid phase, this value further increases as the orientational order improves. The angular susceptibility is primarily used to detect the phase transition from liquid to hexatic as it becomes finite at this transition. When using this parameter one needs to be aware that increasing the lattice size leads to reduced parallel alignment even for a solid and therefore to a reduced value of j6. For this reason, it can only be used reliably to compare lattices of the same size [173]. Orientational Correlation Function The spatial correlation of the hexatic order parameter 6 can be used to describe the long-range order of a 2D system and especially the phase transition between the liquid and the hexatic phase [33, 36]. When including binning,〈 it is calcula[ted via ]〉 6(A) = k6(8)k6( 9) \ (A8 9 − A) − \ (A8 9 − A − ΔA) , (2.101) where the variables are defined as in Equation 2.95 and k6 as in Equation 2.98. In the liquid phase, this ordering is expected to decay exponentially while in the hexatic phase the function has an algebraic decay  (A) ∼ A−[6 6 . Following from the KTHNY theory, the liquid-hexatic phase transition is at [6 = 0.25 [36, 39]. In the solid phase 6(A) is constant for large distances which is however difficult to determine in experiments or simulations due to finite size of the analyzed samples. 36 2. Theoretical Background Therefore the hexatic-solid phase transition has to be detected via other means. 2.5.5. Translational Order Solid order is defined by long-range translational correlations. These correlations can be described based on the transla tional order pa rameter [174] 1 ∑# $Tr = exp (8@® · A®8)# 1 , (2.102)9= where @® is the wave vector where the structure factor of a perfect hexagonal lattice has its main peaks. It can be calculated via @® = (cos(\), sin(\)) · @ (2.103) with \ the orientation angle of the lattice and @ = √4c with the lattice spacing 0. 30 This order parameter is zero in both, the liquid and the hexatic phase, and only reaches finite values in the solid phase. For this reason, it is used to determine the hexatic to solid phase transition. Similar to the angular susceptibility, this parameter depends on the lattice size and is smaller for larger lattices. Translational Correlation Function In addition to the order parameter given above, one can also define〈 the translational correlation function〉 @ as [175]  (A) = 48@®·(A®8−® [ ] A 9 ) @ \ (A8 9 − A) − \ (A8 9 − A − ΔA) . (2.104) While the translational order parameter is most useful for smaller systems, this correlation function is best suited for describing large systems. A 2D system is assumed to be in the solid phase when @ decays algebraically  (A) ∼ A−[@@ with an exponent [ 1@ < 3 . In the hexatic as well as the liquid phase, @ is expected to decay exponentially. This correlation function is used to determine the hexatic to solid phase transition. 2.6. Experimental Methods While this thesis is a fully theoretical work without any experiments performed by myself, I shortly mention experimental methods which are required to generate the experimental results used in this thesis. This includes the metallic thin-film systems used in this thesis as well as the microscopy method used for imaging of 37 2. Theoretical Background skyrmions. Afterwards I describe the particle tracking method which I also have been using myself. 2.6.1. Materials Skyrmions presented in this thesis were observed in different multilayer stacks consisting of nm-scale layers deposited via sputtering. All stacks used to generate results which are part of this thesis have a similar structure of Ta/Co20Fe60B20/Ta/ MgO/X where X is either Ta or HfO2. Thicknesses of different layers vary between different stacks [22, 41, 44, 67, 99, 108]. The Ta layer at bottom is a heavy metal layer, which is required to get a significant strength of the DMI to allow for formation of chiral skyrmions. Additionally, this layer is responsible for a strong spin Hall effect allowing for large spin-orbit torques. The Co20Fe60B20 layer is a ferromagnetic layer creating the ferromagnetic exchange interaction in this material stack which is required for imaging magnetic structures. The MgO layer is used to enhance the perpendicular magnetic anisotropy which is further tuned by the thin Tantalum interlayer separating CoFeB and MgO. The final layer of Ta or HfO2 is required to prevent oxidation of the sample and therefore guarantee longevity of the sample. Using HfO2 additionally leads to a hugely increased magnetic contrast in the below-discussed magnetic microscopy methods as it is transparent to visible light which simplifies tracking of skyrmions. 2.6.2. Magneto-Optic Kerr Effect (MOKE) Microscopy Videos of skyrmions are recorded via the out-of-plane magnetic field component measured by a Kerr microscope by evico magnetics GmbH which uses the Magneto- Optic Kerr Effect. This effect describes the change of intensity and polarization of light when it is reflected from a magnetic surface [176]. The polarization of light is changed by the Kerr angle \K which is proportional to the magnetization on a point of the sample. The effect is similar to the Faraday effect which describes transmission of light instead of reflection. MOKE microscopy does not give information about the absolute magnetization at some point but only about the relative values which is however sufficient to locate skyrmions and material defects on a sample. More specifically, in the experiments described in this thesis, a polar MOKE microscope is used. This device captures the out-of-plane magnetic field by using linear polarized light which is changing its polarization to elliptical when being reflected. An example picture from a MOKE video of skyrmions is shown in Figure 2.11. Due to the large size of the skyrmion core and the thin domain wall in comparison to the picture shown in Figure 2.6, the size of the dots on the MOKE image is close to the size of the skyrmion core. The radius in which a skyrmion interacts with other skyrmions however is larger than the size visible in the video. 38 2. Theoretical Background Figure 2.11.: Background-subtracted image of a skyrmion lattice between two magnetic material boundaries recorded by MOKE. The colors indicate the magne- tization at individual pixels relative to the background image. White circles are skyrmions. The black and white lines at top and bottom arise due to a slight shift of the sample compared to the background measurement. 2.6.3. Skyrmion Tracking with Trackpy Videos generated by MOKE microscopy are gray scale movies of the sample magnetization in which skyrmion cores are depicted by a different color than the magnetic background. Therefore, the task of tracking skyrmions in MOKE videos is basically the task of tracking white points on a black surface or vice- versa. This includes complications caused by thermal noise and by material defects which also lead to differently colored points. To perform this analysis, the python package trackpy [177, 178] is used. Trackpy fits 2D Gaussians to local maxima (or minima) in intensity depending on thresholds for the peak size, width, total peak intensity and distance between distinct peaks. Usually, one uses tracking masks slightly larger than the object one wants to track. To increase contrast between skyrmions and the magnetic background, at first the magnetic background is measured without skyrmions and then this picture is subtracted from the pictures containing skyrmions. Further preprocessing is done by using a blurring kernel to average out noise on short distances. Other objects like material boundaries or defects are tracked manually. Major defects are recognizable by a much stronger signal than the signal of a skyrmion and their incapacity to move. Boundaries are fitted by fitting the intensity on both sides of a boundary with a straight line and using the point where the intensity is the average of both sides as the position of the boundary. This process of finding material boundaries is done using the background image. A more thorough explanation of the boundary fitting process is given in Appendix A.2.1. 39 2. Theoretical Background Further improvements of the detection technique are required for long trajectories which include movement of the sample. Here, the position of easily detectable fixed points like material defects is detected. Then, their movement is subtracted from the positions in the videos to ensure successful background subtraction and a correct sampling of other position-dependent features like an energy landscape of the sample. 2.6.4. Linking of Skyrmion Trajectories After tracking positions of individual skyrmions in experimental videos, these positions can be linked to obtain skyrmion trajectories which are required to analyze dynamic effects such as diffusion. Linking is done using the trackpy package [177, 178] which has already been used for tracking. The linking algorithm is trying to minimize the error in linking using an assumption based on the type of motion present in the data. In the case of standard Markovian Brownian motion, the velocity of particles at a specific point in time is independent from the velocity before and therefore cannot be used for the prediction of particle positions. For this reason, the optimal predictor for a position of particle 8 at C = C1 based on a position at C = C0 is %(C1, C0, G®8 (C0)) = G®8 (C0), (2.105) i.e. the position in the previous frame. To limit the computation time required for linking, a maximum displacement of a particle from the prediction Asearch is defined so that particles are only considered for linking if they fulfill | |%(C1, C0, G®8 (C0)) − G®8 (C1) | | < Asearch. (2.106) This limitation reduces the computation time of O(# !) which would be required for this algorithm to a polynomial one. The ideal linking is now found by minimizing the quadratic d∑istance between the prediction∑and the detected positions# # | |%(C , C , G® (C )) − G® (C ) | |21 0 8 0 8 1 = | |G®8 (C0) − G®8 (C1) | |2. (2.107) 8=1 8=1 In case the number of particles differs between C0 and C1 or distances are too large to be assigned, particles are considered missing and no connection is made. Due to possible errors in tracking, the a variable can be set to allow for a number of skipped frames within a connected trajectory. In this case, a particle missing at C1 uses the position at C0 to be linked with a particle at C2. 40 3. Skyrmion Interaction Potentials 3. Skyrmion Interaction Potentials The description of skyrmions in a particle-based model, such as the Thiele model, requires description of interactions between different skyrmions, as well as between skyrmions and parts of the magnetic material such as material boundaries and defects. While there are several publications calculating the skyrmion-skyrmion interaction potential using theoretical considerations [21, 59] and micromagnetic simulations [21, 54, 55, 59], the method presented here uses experimental data to calculate skyrmion interaction potentials. 3.1. Derivation of Interaction Potentials via Iterative Boltzmann Inversion The following section presents results for the determination of skyrmion-skyrmion and skyrmion-boundary interaction potentials. The bulk of these results is published in [108]. Experiments presented here were predominantly performed by Yuqing Ge without my participation. Tracking of experimental data was done primarily by Yuqing Ge with my help. Analysis and computer simulations were performed by me. 3.1.1. Experimental Setup and Data Skyrmion-skyrmion and skyrmion-boundary interaction potentials are determined in a system where the skyrmions are confined between two material boundaries in the y-direction. This setup allows for simultaneously determining both kinds of interaction potentials in one experiment. The length of the system in x-direction is much larger than the camera frame and can therefore be considered infinite in this analysis. Positions of skyrmions are detected using the trackpy [177] package while the material boundary is determined directly from the background image. These positions are used to calculate the distribution functions for skyrmions with respect to each other using Equation 2.95 and with respect to the boundary using Equation 2.96. The geometrical constraints in which distribution functions are calculated are displayed in a snapshot from the experimental video shown in Figure 3.1. These constrains are chosen such that the skyrmion-skyrmion radial distribution is mostly independent from effects caused by the material boundary and can therefore also be transferred to bulk systems. To calculate the boundary distribution function, an average over the top and bottom boundary is taken as the two boundaries have the same physical 41 3. Skyrmion Interaction Potentials Figure 3.1.: Background-subtracted MOKE image of skyrmions in a band between two material boundaries. The grayscale indicates the local magnetization relative to the background image. The stripe is significantly longer than the field of view and approximately 80.6 µm wide. Skyrmions are the light-gray circular objects on the dark-gray background. They are tracked using Trackpy and their tracked locations are denoted by red circles with constant radius. This radius shows the size of the tracking mask which is not equal to the radius of skyrmions. The material boundaries are marked by the yellow dashed lines at top and bottom. The cyan and green lines define the areas used for the calculation of 6(A). Reference particles for the calculation of 6(A) are located in the cyan area, particles used as interaction partners for this calculation are found in the green area. The distance between the two lines is 20 µm which is also the maximum distance for which 6(A) is computed. This figure is adapted from our publication [108]. properties and therefore are expected to interact identically with skyrmions. One however has to note that initially the calculated boundary distribution functions exhibited a shift of 0.6 µm with respect to each other which can most probably be explained by drifting of the sample compared to the background image taken before the experiment and used for determination of the boundaries. Because the exact drift amplitude cannot be determined directly from the Kerr videos, both sides are moved by 0.3 µm each so that the first maxima of both sides are at the same distance from the corrected boundary. The top and bottom boundary distribution functions as well as their sum are shown in Figure 3.2. Due to pinning effects in the material sample as described in Section 2.4, before calculating an interaction potential, one needs to analyze the radial distribution function and make sure that the experimental results are actually reproducible for other samples and not dominated by pinning-induced disorder. While for the skyrmion-skyrmion interaction the distribution has a shape that does not show significant anomalies, the boundary distribution function indeed contains an anomaly affecting the first peak of 6Bnd(A) at whose falling side another small 42 3. Skyrmion Interaction Potentials Figure 3.2.: Experimental 6Bnd(A) plotted individually for the top and bottom boundary as well as combined. The plotted curves already contain the correction of wall positions described in the text. The red circle indicates an additional peak which is probably caused by pinning and therefore not used for potential determination. increase in skyrmion numbers is found at around 7.5 µm as indicated by the red circle in Figure 3.2. This can be explained by multiple pinning sites with a similar distance to the boundary. In a system without pinning, there would be no physical reason why there should be a further peak of 6Bnd(A) at this distance. This error is corrected by applying a running average to this function before running any simulations. The running average as defined in Appendix A.1 is performed over 8 neighbors from the first peak to the first minimum. To better understand why such a correction is necessary for 6Bnd(A) and not for 6(A), one needs to take a look at the different statistics which can be used to determine the two functions. The skyrmion-skyrmion interaction is an interaction between pairs of skyrmions and therefore every skyrmion contributes to its statistics around cA2cutd/2 times assuming an equal distribution of skyrmions in a circle of radius Acut around every skyrmion and not counting pairs twice. On the other hand, the skyrmion-boundary interaction can only count each skyrmion at most once per frame when analyzing distances smaller than half the channel width. Additionally, pinning has a stronger influence on the boundary interaction because a strong pinning site with a given distance to the boundary always leads to a skyrmion at fixed distance to the boundary, while a fixed point for the skyrmion-skyrmion interaction requires two strong pinning sites pinning both skyrmions at a given distance to each other. 43 3. Skyrmion Interaction Potentials 3.1.2. Computer Simulation To calculate skyrmion interaction potentials from these distribution functions, one has to reproduce the experimental system in computer simulations. Therein, simulations in the Thiele model are performed with standard parameters of :B) = 1, W = 1,  = 0.25 (accounting for a small skyrmion Hall angle of around 14° in this system) and ΔC = 10−5. The simulated system contains repulsive boundaries in y-direction (described by a skyrmion-boundary interaction potential) and periodic boundaries in x-direction due the large extension of the system in this direction. The skyrmion density of 0.0228 skyrmions per µm2 as well as the distance between the top and bottom boundary of 80.6 µm are chosen to be the same as in the experiment. With this, one is only missing interaction potentials for a simple particle-based model of skyrmions. To obtain information about these potentials, one can use the distribution functions described above. In a previous work by my colleagues, the problem of determining interaction potentials was tackled by assuming a spe(cific potential shape - in this case asoft-sphere potential f )−= +power(A) = (3.1) A - and testing different parameters to obtain the parameter values with best agree- ment to the experimental distribution [44]. This, however, comes with the major limitation that such testing relies on assumptions on the potential shape limiting the agreement between experiment and simulation. Additionally, when also simu- lating interactions with the material’s boundary, the number of fitting parameters doubles leading to a four-dimensional parameter space and thereby complicating the parameter fitting procedure. For these reasons, a different approach for deter- mining interaction potentials without these limitations was chosen here. The most straightforward method for calculating interaction potentials from distribution functions without such assumptions is using the potential of mean force [169] + (A) = −:B) ln 6(A), (3.2) which has been used in several works to e.g. determine interaction potentials for colloids [178] or superconducting vortices [179]. Here, interaction potentials are fitted up to the first minimum after the first peak of the interaction function because especially for the skyrmion-boundary interaction this region is sampled best. Further peaks of this function have worse statistics and are also more influenced by effects of skyrmion-skyrmion interactions and therefore not able to reproduce a realistic skyrmion-boundary interaction. The potential of mean force, however, is only an exact solution at zero density as it neglects collective effects of particle lattices [62, 180]. For finite density, it is 44 3. Skyrmion Interaction Potentials (a) Experiment 2.52.5 (b) ExperimentSimulation Corrected 2.0 2.0 Simulation 1.5 1.5 1.0 1.0 0.5 0.5 0.0 0.0 0 5 10 15 20 0 5 10 15 20 r [ m] r [ m] Figure 3.3.: Experimental 6(A) and simulated 6(A) using potentials of mean force for both, skyrmion-skyrmion interaction (a) and skyrmion-boundary interaction (b). The bin width is ΔA = 0.1. only a first estimate and usually not sufficient to reproduce the spatial distribution of particles as in the experiment described above. Figure 3.3 displays the radial distribution function of simulations performed using the potential of mean force both for skyrmion-skyrmion and skyrmion-boundary interactions in comparison to the experimental radial distribution. One can see large differences already in the first peak indicating that indeed the potential of mean force does not reproduce the experimental ordering of skyrmions. The model of skyrmions described here uses an isotropic interaction potential. This comes with some limitations as skyrmions are assumed to be monodisperse, isotropic, and non-deformable. Also effects such as skyrmion annihilation and material inhomogeneities are not accounted for in this minimal model. A throughout discussion of these limitations is found in the last chapter of this thesis. 3.1.3. Iterative Boltzmann Inversion While the potential of mean force alone is not enough to describe skyrmion in- teractions, potentials determined this way can be improved using the Iterative Boltzmann Inversion (IBI) method [62]. This method matches a target radial distribution function with the radial distribution function of a simulation based on the Boltzmann weight of some state 8 exp(−8/:B)) ∼ %8, (3.3) where 8 is the energy of this state. The Boltzmann weight is proportional to the equilibrium occurrence probability %8 of the state. Inverting this relation, one 45 g(r) gBnd(r) 3. Skyrmion Interaction Potentials (a) Experiment (b) 2.0 Simulation 2.0 1.5 1.5 1.0 1.0 0.5 0.5 Experiment Corrected 0.0 0.0 Simulation 0 5 10 15 20 0 5 10 15 20 r [ m] r [ m] Figure 3.4.: Experimental 6(A) and simulated 6(A) using potentials determined by the iterative Boltzmann inversion method for skyrmion-skyrmion interaction (a) and skyrmion-boundary interaction (b). The simulated distribution function is generated by a particle-based Brownian dynamics simulation after performing 500 update steps. The bin width is ΔA = 0.1. This figure is taken from our publication [108]. arrives at  ∼ :B) ln(%) which motivates the derivation of the IBI method. IBI has first been used in soft matter physics to calculate coarse-grained potentials for much more detailed atomistic systems [62]. Like the potential of mean force, the IBI method also does not use any assumptions on the shape of the potentials and can therefore in general describe both, purely repulsive and partly attractive interactions. While this method is broadly used in soft matter theory [181] and the possibility to be used for experimental systems has also been discussed [64], it has so far rarely been used for experimental analysis. Using this method, the interaction potential is updated iteratively depending on the ratio of the initial radial distribution function 6(A) one wants to reproduce and on the simulated distribution function 68 (A) calculated from simulations in iteration step 8 as 68 (A) +8+1(A) = +8 (A) + U:B) ln ( ) , (3.4)6 A with a parameter U ∈ (0, 1] determining the speed with which the potential is updated [62]. Very large values of U lead to faster convergence but also to oscillations around the correct potential while too low U leads to slower convergence. For the given system, simulations where performed with U = 0.2, achieving stable results within reasonable time. Using this method, after 500 iteration steps of the IBI method, one arrives at a potential which can very well reproduce the ordering of skyrmions with respect to each other and the system boundary as it is shown in Figure 3.4. Not only 46 g(r) gBnd(r) 3. Skyrmion Interaction Potentials 10 1 Skyrmion-Skyrmion Skyrmion-Wall 10 2 10 3 10 4 10 5 0 100 200 300 400 500 Run number F( igure 3.5.: Mean squared difference between simulated and experimental∑dis-tribution function)s over 500 IBI steps with a logarithmic y-scale: 1 ## 8=12 6exp(A8) − 6sim(A8) . This figure is taken from our publication [108]. the first peak, which is used for fitting the potential, but also the second and third peak of the skyrmion-skyrmion distribution are reproduced considerably well. This shows that also longer-range ordering of the skyrmions is reproduced by the short-range two-particle interaction calculated by IBI. For the skyrmion-boundary interaction the fitted region of the first peak is also reproduced well. Beyond this peak, however, one sees significant deviations between the experimental and simulated functions. This is due to the pinning-induced peaks in the experimental system which cannot be described in the pinning-free simulation as well as the smaller amount of statistics compare to 6(A). The disorder of the peaks in the experimental boundary distribution function is also visible when comparing the top and bottom functions in Figure 3.2 which do not agree to each other very well. To analyze the convergence of this method, the mean squared difference between the experimental and the simulated results can be used. This difference is showing a clear convergence after around 300 steps both for the skyrmion-skyrmion and the skyrmion-boundary interaction as shown in Figure 3.5. 3.1.4. Analysis of Results To further test the capability of the performed computer simulations and the IBI method to reproduce the experiment, a few more quantifiers are considered. The hexagonal order of a 2D system is often described by the local order parameter k6. For the inner part of the system not affected by boundaries (cyan region in Figure 3.1), the average order parameter 〈|k6 |〉 is 0.486 ± 0.008 for the experiment 47 Mean-Squared Difference 3. Skyrmion Interaction Potentials Figure 3.6.: Two-dimensional radial distribution function 6(G, H) for (a) experiment and (b) simulation. The functions are normalized to a value of one for the point with most skyrmion occurrences as in our published work [108]. A difference between the two results is visible e.g. in the second and third layer which are more pronounced in the simulation. The same figure with the standard normalization from Equation 2.97 is given in Appendix A.3. This figure is taken from our publication [108]. and 0.496 ± 0.007 for the simulation (errors are the standard error of the mean given in Appendix A.2) indicating a good agreement of the local structure. Another indicator for the local structure is the percentage of skyrmions with exactly six neighbors, %6, whose value in the same region is %6 = (59.2±0.8)% in the experiment and %6 = (62.0 ± 1.1)% in the simulation. To visualize skyrmion ordering beyond the first shell, the two-dimensional radial distribution function can be used which is displayed in Figure 3.6. This function is calculated within the same constraints as the 1D skyrmion-skyrmion distribution. While the first layer displays a good sixfold order, the ordering is quickly decaying beyond this in experiment and simulation. This structure of the two-dimensional distribution clearly shows a system without long-range order and therefore a liquid [44, 48]. These different quantifiers show a remarkable agreement between experiment and simulation considering that pinning effects have been neglected in the simulation. This indicates that the influence of pinning is reduced for skyrmion lattices compared to the effect of pinning on skyrmions at low density. Additionally, this agreement demonstrates that the quantitative description of skyrmion lattices as the one observed in this experiment is possible in a simple Thiele model not accounting for deformation, size polydispersity or multi-skyrmion interactions. Good agreement between experimental and simulated distribution functions is a 48 3. Skyrmion Interaction Potentials 17.5 Skyrmion-Skyrmion 15.0 Fit Skyrmion-Skyrmion 12.5 Skyrmion-Boundary Fit Skyrmion-Boundary 10.0 7.5 5.0 2.5 0.0 4 5 6 7 8 9 r [ m] Figure 3.7.: Skyrmion-skyrmion and skyrmion-boundary potentials determined by iterative Boltzmann inversion and fits to the potentials with an exponential potential + = 0 exp(−A/1). Fit parameters are 0SkSk = 735.1 :B) and 1SkSk = 1.079 µm for skyrmion-skyrmion interaction and 0SkBnd = 176.7 :B) and 1SkBnd = 1.673 µm for skyrmion-boundary interaction. Errors of the fit are smaller than the last digit. This figure is taken from our publication [108]. requirement for meaningful interaction potentials to be calculated using the IBI method. The potentials are given in Figure 3.7 and are both fully repulsive as expected for chiral skyrmions [21, 59]. Due to this shape, the potentials can be fitted to a falling exponential function + (A) = 0 · 4−A/1exponential (3.5) with quite high accuracy. The mean squared difference between the result and the exponential fit is 0.07(:B))2 for the skyrmion-skyrmion interaction and 0.43(:B))2 for the skyrmion-boundary interaction. When instead using a Modified Bessel function of second kind 0(A) as proposed by Lin et al. [21], the mean squared difference is 0.08(: 2B)) and 0.51(:B))2, respectively. These values are very close to each other indicating that based on this study one cannot distinguish between these two shapes. On the other hand, a fit to a power-law potential as in Equation 3.1 leads to mean squared differences of 0.29(: 2B)) and 1.14(:B))2, respectively. This indicates that this description is less suitable to describe interactions of skyrmions than using an exponential. While the analysis above leads to potentials, which can well reproduce skyrmion interactions and ordering, one must still make sure that the resulting potential is indeed the only correct potential describing skyrmion interactions. From a mathematical perspective, it has been proven that there is a unique relation 49 V [kBT] 3. Skyrmion Interaction Potentials between a radial distribution function and a pairwise potential of particles leading to this function [182]. However, this proof is only valid in the limit of infinite data and therefore perfectly smooth distribution functions. For a real system with limited data, similar potentials often lead to indistinguishable experimental distribution functions. For this reason, one cannot state that the exact shape of the resulting potential is the only correct one, but as simulations with the calculated potential reproduce the experimental structure, this potential is a good approximation for describing skyrmion interactions. 3.1.5. Potentials for Larger Skyrmions As the experimental conditions such as temperature, magnetic field and exact composition of the multilayer stack vary, the size of skyrmions is also varying. This variation in size leads to different interaction potentials and therefore a limited applicability of the exact results from the above simulations. One can, however, rerun the algorithm described above on a different system to gain potentials de- scribing it. Here, a system where skyrmions are larger by a factor of around 1.5 and where no material boundaries are present in the image of the Kerr microscope is used. For this reason, only the skyrmion-skyrmion interaction is analyzed here. As before, the goal is to reproduce the experimental 6(A) which is shown together with the simulated value in the left panel of Figure 3.8. The right panel shows the skyrmion-skyrmion interaction potentials determined here. The potential has an 25 3.0 Experiment Skyrmion-Skyrmion Simulation 2.5 20 Fit Skyrmion-SkyrmionPaper Skyrmion-Skyrmion 2.0 15 1.5 10 1.0 0.5 5 0.0 0 5 10 15 20 05 6 7 8 9 r [ m] r [ m] Figure 3.8.: Left: Experimental 6(A) and simulated 6(A) using potentials deter- mined by the iterative Boltzmann inversion method for larger skyrmions. Right: Interaction Potentials determined by IBI for these skyrmions. The blue dotted line shows an exponential fit to the results with fit parameters 0SkSk = 2141.8 :B) and 1SkSk = 1.0826 µm. Errors of the fit are smaller than the last digit. The red dotted line shows the fitted potential for smaller skyrmions for reference. 50 g(r) V [kBT] 3. Skyrmion Interaction Potentials exponential form similar to the one determined for smaller skyrmions but a much larger distance at which a significant repulsion is observed. 3.1.6. Applicability to Other Systems The analysis and discussion provided so far focuses on the direct determination of skyrmion interaction positions from an experimental video. The IBI method in general is, however, able to determine potentials describing a much larger class of different systems. This includes e.g. nanometer-scale magnetic systems in which a pair distribution function 6(A) of skyrmions cannot be directly measured. In these systems, a scattering pattern can be recorded using neutron diffraction [46, 124]. This scattering pattern can be used to calculate the pair distribution function via a simple Fourier transform. In these systems one could use this model to validate interaction potentials calculated from micro-magnetic or atomistic calculations which are only feasible for nanometer-sized skyrmions. Thereby, one cannot only describe repulsive skyrmion-skyrmion interactions of DMI-stabilized chiral skyrmions but also attractive skyrmion-skyrmion interactions predicted [183–185] and experimentally determined for skyrmions in the conical phase [186, 187], as well as for skyrmions in frustrated magnets [125, 126]. Additionally, while the IBI method is only used for skyrmions here, it can in general be used for all quasi-particles or particles with isotropic interactions whose distribution functions can be accessed either directly or indirectly by analyzing scattering patterns. These systems include e.g. colloids [178] or other spin structures like skyrmionium [188]. Different from these cases, the model is explicitly not applicable for systems with strongly non-isotropic interactions. This includes e.g. skyrmions under strong applied currents leading to deformations of the skyrmion shape [189]. Also worm domains or asymmetric skyrmion-like objects such as antiskyrmions are not described by this model. 3.2. Pinning of Skyrmions at Material Defects Beyond the description of skyrmion-skyrmion and skyrmion-boundary interactions, a comprehensive description of skyrmion motion additionally requires realistic pinning profiles. The approach presented here describes pinning as skyrmion- center pinning which comes with limitations as discussed in Section 2.4. While this leads to a description of pinning, which is only applicable to skyrmions of a specific size, this simple model can still describe pinning effects such as reduced diffusion or obstructed ordering. The determination of such profiles is based on the analysis of skyrmion positions on a given sample as shown in the occurrence 51 3. Skyrmion Interaction Potentials Figure 3.9.: (a) Logarithmically plotted histogram of skyrmion occurrence used to determine the pinning strength. (b) Energy landscape of the same sample calculated from the occurrence distribution. The potential at bins with zero occurrences is set to  = 0 preventing nonphysical divergences and allowing for a smooth energy landscape. The data is taken from a paper by Raphael Gruber et al. [67]. map in Figure 3.9(a). In this sample, the observed region is limited by a material boundary a few micrometers away from the edges shown in the figure. To create a potential landscape for simulations incorporating pinning effects, one can take the negative logarithm of the skyrmion occurrence distribution [67] + (G, H) = −:B) ln # (G, H), (3.6) where # (G, H) is the number of skyrmions found in one bin over the course of the experiment. The calculation of a reliable energy landscape is however impeded by poor statistics. In several bins no skyrmion has been measured during the whole 52 3. Skyrmion Interaction Potentials Figure 3.10.: (a) Logarithmically plotted histogram of skyrmion occurrence recorded by Tobias Sparmann in a different experiment. (b) Energy landscape for the same experiment. measurement as shown by the white pixels in Figure 3.9(a). This statistics also cannot be reliably improved by increasing the skyrmion density because for a high density, skyrmion positions are primarily governed by the ordering of a skyrmion lattice and not by pinning effects. This would lead to a skyrmion occurrence map, which just shows lattice positions and a few pinning sites together and therefore cannot be used to model pinning distinct from skyrmion-skyrmion interactions. This effect is also shown in the profile given in Appendix A.4, which shows occurrence peaks placed on lattice positions of a hexagonal skyrmion lattice. To address the 53 3. Skyrmion Interaction Potentials issue of empty bins, the energy of the former is set to + = 0 (corresponding to # = 1). This way the skyrmions can still diffuse to these locations and are not repelled by an unphysical infinite potential. The resulting energy landscape is given in Figure 3.9(b). This approach, however, comes with the limitation that possible repulsive positions on the sample and under-sampled areas cannot be distinguished and thereby attractive skyrmion pinning effects are likely overestimated. Improving this would require much more statistics and thereby no remaining empty points on the sample. Another big issue when determining pinning landscapes is the limited time resolution of the MOKE setup, which usually allows for only 16 pictures per second. This leads to recording of pictures averaged over 62.5 ms. Therefore positions where a skyrmion stays for only e.g. 10 ms are not detected as their signal is too weak to be tracked in one frame. This especially applies to regions close to a strong attractive potential where the skyrmion usually crosses with high velocity. This effect leads to an overestimation of pinning depths as flanks of pinning sites are underestimated in the pinning landscape. Without further investigation including a higher time resolution of the experimental setup it is unclear how large these two effects are and whether they can mitigate each other. Even with these issues, the pinning landscape is already a useful tool for analyzing the impact of realistic pinning on skyrmion motion and ordering. Due to the limited size and statistics of this experimental recording, further analysis is performed on different experimental data and the landscape shown in Figure 3.10. While the maximum depth of the pinning potential is similar to the potential determined above, it is more homogeneous and therefore apart from the largest pinning sites, skyrmion motion is less impeded here. Computer simulations of single skyrmions using the potential determined via Equation 3.6 lead to an occurrence map which looks very similar to the experimental occurrence map as shown in Appendix A.5. This indicates that the described method is appropriate to at least reproduce experimental pinning of individual skyrmions. The only difference between the experimental and the simulated landscape besides statistical effects are the points with 0 occurrences in the experiment. Due to the assignment of a zero potential to these points, they are occasionally filled with skyrmions in the simulation. This is also required for ergodicity as without setting this potential to a finite value some of the deepest pinning sites would be surrounded by an impassible wall. 3.2.1. IBI with Pinning Effects As skyrmions which are used to determine skyrmion interactions are also subject to pinning effects, one should test how this affects the determination of interaction potentials. This is realized by adding the pinning landscape to the simulations used to determine potentials with the IBI method. As the pinning landscape in 54 3. Skyrmion Interaction Potentials Without Pinning 10 2 With Pinning 10 3 10 40 100 200 300 400 500 Run Number Figure 3.11.: Mean squared difference (MSD) between simulated and experimental radial distribution functions 6(A) over 500 IBI steps for simulations with and without pinning. this sample has been determined more accurately, the following simulations are performed for the data from Section 3.1.5. The resulting potential parameters after 500 steps are 0SkSk = 2139.2 :B) and 1SkSk = 1.1150 µm which is less then a 5% deviation in the fit parameters from potentials determined without pinning given in Figure 3.8. A similarity between the results can also be observed in the mean squared difference between experimental and simulated 6(A) shown in Figure 3.11 which decreases slightly slower when adding pinning but reaching a very similar final value. These results show that equilibration is slowed down by pinning effects but pinning does not lead to a systematic breakdown of the IBI potential determination method. This demonstrates that one can indeed neglect pinning when determining skyrmion-skyrmion interaction potentials. This result was expected as skyrmion pinning is suppressed by high-density skyrmion lattices and the IBI method typically requires such lattices. 3.2.2. Modeling of Ideal Pinning Sites As an alternative to the realistic pinning landscape, the general impact of pinning on skyrmions can also be analyzed using a simpler model. For this, a periodic square lattice of pinning sites is used. Due to the hexagonal ordering of unpinned skyrmions, the chosen square pinning lattice is incommensurate with the typical ordering of the system and therefore impedes skyrmion lattice formation similar to a realistic pinning background. To model pinning sites in such a lattice, one can 55 Mean-Squared Difference 3. Skyrmion Interaction Potentials assume harmonic pinning potentials [21] + (A) = − :2 (A − 0) 2, (3.7) 0 where 0 is the radius of the pinning site and A is the distance to the center of the pinning site. This potential shape is used in later chapters to describe homogeneous pinning with known properties. 56 4. Dynamics and Phase Behavior of Skyrmion Lattices 4. Dynamics and Phase Behavior of Skyrmion Lattices Skyrmions are a very promising system for analyzing 2D phase behavior. They are readily accessible experimentally by direct measurements of positions in a Kerr microscope. The skyrmion size and therefore density can be changed easily by changing the applied magnetic field. Additionally, skyrmions can be readily moved by applying spin-orbit torques and allow for controlled creation and annihilation by magnetic fields while exhibiting a high stability when not excited. The tunability of skyrmions can also be used to speed up and enhance ordering of a skyrmion lattice. Experiments performed on skyrmions can easily be accompanied by computer simulations providing an even better insight into details of phase transitions. In this chapter, I show several computer simulation results, which improve the understanding of skyrmion dynamics and skyrmion lattices. Some of these results have already been published in a similar manner and have only been reproduced by me, most of them are not related to any publication. I start by analyzing pinned diffusion by tuning pinning strength and skyrmion density. Afterwards I show simulations of pinned skyrmions with an applied pushing force. The next section contains analysis of the interplay between Magnus force and skyrmion density and contains results published in a publication to which I contributed [109]. After this I present different parameters describing ordering of a 2D lattice and their evolution close to the liquid-hexatic and hexatic-solid phase transitions. These results are then used to describe a well-ordered experimental skyrmion lattice. 4.1. Pinned Diffusion Pinning effects lead to a significant decrease in diffusion. To quantify the impact of such pinning effects, simulations using the pinning potential given in Figure 3.10 are performed. Here, skyrmions are simulated using the realistic potentials determined in Section 3.1.5. For these simulations, the number and therefore density of skyrmions is varied. Simulations are performed with and without a pinning background generating the results presented in Figure 4.1. In the two cases, an increased skyrmion density leads to a reduced diffusion of individual skyrmions. When comparing the diffusion constants in the pinned and the unpinned case, the pinned diffusion is significantly lower than the unpinned one except for very high densities. The ratio between pinned and unpinned diffusion is conserved over a relatively large range of densities from zero up to around 0.005 µm−2. In this region, 57 4. Dynamics and Phase Behavior of Skyrmion Lattices (a) (b) Not Pinned 0.8 Pinned 3 0.6 2 0.4 0.2 1 00.0.0000 0.0025 0.0050 0.0075 0.0100 00.0000 0.0025 0.0050 0.0075 0.0100 Skyrmion Density [ m 2] Skyrmion Density [ m 2] Figure 4.1.: (a) Skyrmion diffusion constants for different densities of skyrmions confined in a rectangular area with and without the pinning energy background displayed in Figure 3.10. Errors are smaller than the points in the figure. Lines are an aid for visualization. (b) Ratio of unpinned to pinned diffusion constants for different skyrmion densities. skyrmion motion is dominated by pinning and not by lattice effects. For higher densities, the relative difference declines as a consequence of two factors. On the one hand, the limited amount of deep pinning sites available in the sample leads to some skyrmions being located at positions with weaker to no pinning. This allows for them to diffuse in a way similar to the pinning-free case. The other effect is the formation of a skyrmion lattice which leads to skyrmion diffusivity being more limited by other skyrmions than by pinning sites. This is due to lattices stipulating preferred skyrmion positions different from preferred positions due to pinning. At high densities, the lattice energy becomes larger than the pinning energy and therefore the lattice effects dominate [50, 57]. As this effect only depends on density and not on pinning, it reduces the difference between the two cases. This effect also explains the even further reduced difference for very high density because a high-density skyrmion lattice approaching the hexatic or solid phase restricts skyrmion positions irrespective of pinning effects by preventing movement beyond a lattice position. Another property to look at in the context of pinned diffusion is how the strength of pinning affects diffusion. This calculation is performed here by multiplying the pinning potential given in Figure 3.9 by a constant factor. With this potential, simulations of skyrmions without skyrmion-skyrmion interactions are performed depicting the case of skyrmion density d = 0. The dependence of the diffusion constant on pinning strength is given in Figure 4.2. The plot shows clearly that 58 Diffusion Constant D [ m2/ t] Dunpinned/Dpinned 4. Dynamics and Phase Behavior of Skyrmion Lattices 100 No Pinning (Theory) 10 1 Simulations 10 2 10 3 10 4 0 1 2 3 4 5 Pinning Strength Factor Figure 4.2.: Dependence of the diffusion constant  of individual skyrmions on pinning strength for the energy landscape given in Figure 3.9 multiplied with some factor. The red line indicates the theory value without pinning from the Thiele equation as given in Equation 2.83. the diffusion is monotonically decreasing with increasing pinning strength. Beyond this, no simple functional dependence between diffusion and pinning strength in a realistic pinning landscape could be identified. Very weak pinning below half the original pinning strength leads to only a small reduction of diffusion as in this case a lot of pinning sites are of a similar amplitude as thermal fluctuations. From this point on the slope falls faster as skyrmions are more and more strongly pinned. Beyond approximately twice the original value, this effect becomes weaker again as occasional movement between the deep pinning sites, which is responsible for most diffusion in this region, is suppressed less. 4.2. Skyrmion Hall Angle in Different Regimes One important property differentiating chiral skyrmion systems from other 2D systems such as colloids is the skyrmion Hall angle which was defined in Equa- tion 2.64. While single skyrmions without pinning effects move in an angle equal to the intrinsic skyrmion Hall angle \Sk,int with respect to the force, skyrmion lattices on a non-flat energy landscape react differently to applied forces. This behavior can in general be described by three different regimes, which are depending on the relation between the applied force and the pinning force [78, 81, 82]. Results similar to the ones presented here have already been published by Reichhardt et al. for similar systems [78, 81, 82] and part of the behavior shown here has already been demonstrated experimentally [72, 74]. 59 Diffusion Constant D [ m2/s] 4. Dynamics and Phase Behavior of Skyrmion Lattices • Pinned regime: At very low applied force in relation to pinning, the force is not able to push skyrmions out of their pinning sites and therefore does not generate any movement of skyrmions. As there is neither parallel nor perpendicular movement, the skyrmion Hall angle cannot be determined reliably. If pinning forces in a system are smaller than thermal fluctuations, this regime is not present. • Creep regime: When increasing the force, some skyrmions are moved out of their pinning sites and along the applied force. Movement perpendicular to the applied force is however still strongly suppressed as skyrmions are heavily pinned. The skyrmion Hall angle in this regime is close to zero. • Flow regime: When the force goes beyond forces in the creep regime, individ- ual skyrmions start moving with some angle respective to the force. However, the forces are typically not strong enough to induce movement of the whole lattice and therefore leads to an angle below the ideal skyrmion Hall angle of an unconfined single skyrmion. In this regime, \Sk increases monotonically with the applied force up to \Sk,int. Beyond some critical force, this movement becomes strong enough to displace all skyrmions and overcome most pinning in the system. In this case, the skyrmions move with an angle \Sk,int with respect to the force. When further increasing the applied force, this leads to skyrmion deformation which increases the skyrmion Hall angle beyond \Sk,int [74]. As the model used here is not able to describe deformation, this effect in omitted in the analysis presented here. To demonstrate this behavior, a simulation of skyrmions with a large Hall angle of \Sk,int = 80.06° according to the choice by Reichhardt [78] was performed ( = 0.985 and W = 0.176). Skyrmions are placed in a periodic system with a density of 0.033 µm−2 and harmonic pinning sites as defined in Equation 3.7 with : = 10:B) and radius 0 = 1 µm arranged in a square lattice with distance 6 µm. Skyrmion- skyrmion interactions are described by the exponential function fitted in Figure 3.7. The results for computer simulations with different applied forces  given in Figure 4.3 indicate the occurrence of the regimes mentioned above as the system’s reaction to applied forces. For very low forces   1, there is at first no significant movement in any direction as expected for the pinned regime. For somewhat larger forces still no movement perpendicular to the applied force and only very small movement in direction of the force is present. This results in a skyrmion Hall angle of \Sk ≈ 0 as predicted for the creep regime. In a large intermediate region of forces up to around  = 20, the velocity of skyrmions both perpendicular and parallel to the applied force increases monotonically as one would expect for the flow regime. Here, the perpendicular velocity increases significantly faster than 60 4. Dynamics and Phase Behavior of Skyrmion Lattices (a) 80 (b) 102 v|| 60 v 100 40 6 4 10 2 20 2 0 40 5 10 15 10 0 F 100 101 102 100 101 102 Force F Force F Figure 4.3.: (a) Skyrmion Hall angle \Sk as a function of the applied force . The applied force is plotted logarithmically. The red line indicates the intrinsic skyrmion Hall angle \Sk,int of the system. The inset shows tan(\Sk) with the red line indicating tan(\Sk,int). (b) Velocities parallel and perpendicular to the direction of the applied force plotted against the force on a logarithmic scale. the parallel velocity corresponding to an increasing skyrmion Hall angle. For the smaller forces in this range up to  ≈ 10, the tangent of \Sk is linearly dependent on the applied force with the increase slowing down for higher forces. For forces beyond  ≈ 20 the lattice of skyrmions is flowing freely without hindrance caused by pinning sites. In this state the observed skyrmion Hall angle is equal to the intrinsic angle \Sk = \Sk,int. These theoretical results as well as similar experimental results [72, 74] indicate that in many experiments pinned skyrmions are not going to diffuse according to their intrinsic skyrmion Hall angle but with a smaller one. This is relevant when e.g. using currents to push skyrmions through confinements. Additionally, the dependence of skyrmion motion on the relation between pinning and applied forces can help in determining a scale for applied forces in a simulated system. 4.3. Influence of the Magnus Force on Skyrmion Diffusion Beyond the impact of the Magnus force term on the direction of skyrmion motion, it also impacts the absolute value of diffusion. According to Equation 2.83, a higher Magnus force amplitude  leads to a smaller diffusion constant. This simple formula however only holds for non-interacting skyrmions in a flat energy landscape 61 Skyrmion Hall Angle Sk tan( Sk) Velocity v 4. Dynamics and Phase Behavior of Skyrmion Lattices 1.0 G = 0.00 0.8 G = 0.50 G = 1.00 G = 2.00 0.6 G = 4.00 G = 8.00 0.4 G = 16.00 G = 32.00 0.2 0.00.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Skyrmion Density [ x 2] Figure 4.4.: Average diffusion constant of skyrmions for different densities and different values of the Magnus force amplitude . Connecting lines are added for readability. while skyrmion-skyrmion interactions lead to a non-monotonous impact of the Magnus force. This effect has been analyzed in a shared paper with the group of Ulrich Nowak [109]. As most of the results in this paper were calculated by Daniel Schick, I present similar results from my own simulations that were not included in the paper. To analyze the relationship between Magnus force and skyrmion-skyrmion inter- actions, a minimal model which consists of soft disks with a simple power-law potential + (A) = A−8, Acut = 2.5 and W = 1 is used. While this simple model does not describe skyrmions in any specific system, it is a simple description of skyrmion-like particles. This system exhibits a well-known phase transition from liquid via hexatic to solid between d = 1.193 and d = 1.206 [36] and acts as a good model system similar to the one used in Zázvorka et al. [44]. Simulations are performed for different values of  and different densities to analyze the interplay of skyrmion density and Magnus force. As the blue line in Figure 4.4 shows, without Magnus force, the diffusion is monotonically decreasing with increasing density. At densities beyond the liquid phase, the diffusion constant decreases by orders of magnitude as the system transitions into a solid where diffusion is limited to movement around a particle’s lattice position. For small  (as e.g.  = 0.5 which leads to a realistic Hall angle for micrometer-sized skyrmions of 26.57◦ with W = 1.0) the biggest difference to the Magnus-free case is found for very low values of the skyrmion density. At these densities, the diffusion is however still close to free diffusion as described in Equation 2.83. For higher densities, the difference to the 62 Diffusion Constant D [ x2/ t] 4. Dynamics and Phase Behavior of Skyrmion Lattices Magnus-free case is diminishing. This indicates that the topological suppression of diffusion is limited by a skyrmion lattice. To get a better understanding of this effect, higher values of  are considered. For values of  & 1 the diffusion constant is no longer monotonically decreasing but has a maximum at some finite density. This density at which this maximum occurs becomes larger for larger values of . However, even for very large  the maximal diffusion density is still far below the liquid-hexatic transition density. The increase in diffusion for higher densities is caused by a rotational motion of skyrmions around each other resulting from the interplay of Magnus force and pair repulsion. Instead of moving away from each other in a straight trajectory, skyrmions spiral around themselves. Additionally, while the diffusion constant is increasing, it never goes beyond the diffusion constant found for the same value of W and  = 0. This indicates that while the topological suppression due to the Magnus force is lifted, the damping constant always limits the diffusion and this limit cannot be overcome by the interplay of skyrmion-skyrmion interaction and Magnus force. When increasing the density beyond the maximum of diffusion, all simulated systems show a roughly linearly decaying diffusion constant up to the transition from liquid to liquid-hexatic coexistence at a density of 1.193 [36]. The skyrmion density at the phase transition is independent of the Magnus force amplitude  because the Magnus force does not contribute to the energy of the system [52]. Therefore is alters only dynamic properties of the system such as the diffusion constant but not static properties such as the ordering at a given temperature. The results discussed here fully agree with Daniel Schick’s calculations in our joint publication [109]. 4.4. Ordering of Lattices In computer simulations, systems without any material defects, perfect point-like particles and statistics only limited by the available computational resources can be depicted. Experimental systems on the other hand have several limitations in these regards. While pinning effects and deviations of skyrmions from a model of isotropic quasi-particles impair lattice ordering, the limited resolution of the camera and changes in sample properties over time limit the potential amount of statistics which can be recorded in experiments. For these reasons, some long-range specifiers of lattice ordering described in Section 2.5 such as 6(A) or  (A), which require a large amount of statistics, are difficult to determine. Therefore, pinpointing the location of phase transitions may require analysis of different properties. One local parameter is the number of lattice defects which can be described by Voronoi diagrams and by %6 as given in Section 2.5.3. To get a better understanding of the properties of these diagrams for different phases, simulations are performed using the well-understood = = 8 power-law potential from reference [36] as in the section 63 4. Dynamics and Phase Behavior of Skyrmion Lattices Figure 4.5.: Voronoi tesselation of equilibrated lattices of soft disks at densities of 1.193 (just below the liquid-hexatic transition) (left) and 1.206 (just above the hexatic-solid transition) (right). The color bar is the same for both plots. The plots only show a part containing around 2,000 particles out of a simulation of 40,000 particles in total. above. As Figure 4.5 shows, defects in the liquid phase form several types of different clusters including dislocations consisting of one fivefold-defect and one sevenfold- defect, disclinations as well as larger disordered regions. In the solid phase, defects occur with a much smaller probability and mostly as dislocation pairs consisting of two sevenfold-defects and two fivefold-defects. This different abundance and type of defects in different phases is in accordance with the KTHNY theory [30–34] which predicts that solids primarily feature dislocation pairs [36]. This difference in defect occurrence can be used to differentiate different phases also for systems where long-range order is not clearly distinguishable from the field of view. The effect of strongly increased order around the phase transitions is also clearly visible in the plot in Figure 4.6 showing an increased slope of several measures for lattice ordering. While there is also an increase of these measures inside the liquid or the solid phase, the changes occur much faster close to the phase transitions. Local order of the system is described by %6 as well as 〈|k6 |〉. While a perfect lattice without defects would exhibit %6 = 1, already high-density liquids contain mostly particles with six neighbors as indicated by a %6 of up to around 85% while in solids, %6 is at least around 95% indicating even smaller amounts of defects. When increasing the density beyond the hexatic-solid phase transition density, the value of %6 gets very close to 100% as the energy penalty for topological defects 64 4. Dynamics and Phase Behavior of Skyrmion Lattices 1.0 liquid solid 0.8 0.6 0.4 P6 0.2 | 6| 6 0.0 OTr 1.100 1.125 1.150 1.175 1.200 1.225 1.250 1.275 1.300 Density Figure 4.6.: Different measures for the ordering of a lattice (percentage of sixfold orientation %6, mean orientational order parameter 〈|k6 |〉, angular susceptibility j6 and translational order parameter $Tr) as defined in Section 2.5 for a lattice of 40,000 particles as a function of density. The dashed lines indicate the boundary of the liquid and solid regions [36] with both phase transitions happening in between these lines. Errors of these parameters (calculated as an error of the mean over multiple frames) are smaller than the points in the plot. increases with higher density and therefore the lattice gets closer to perfect lattice order. The density dependence of 〈|k6 |〉 has a similar shape with a value of around 0.6 for a liquid increasing to around 0.8 for a solid. This is in agreement with a value of 〈|k6 |〉 = 0.69 found for the liquid to liquid-hexatic coexistence transition for different potentials including the one used here [44]. In contrast to %6, this parameter does not reach a value of 1 for finite temperatures as there are always local fluctuations at lattice sites even for solids. Beyond these local parameters, the global orientational order of the system is analyzed using the angular susceptibility j6. This parameter reaches zero for an isotropic liquid phase as there is no global orientation in a liquid and reaches finite values only in proximity to the phase transition. As the local parameters, j6 increases strongly along the hexatic and solid phase and then increases much slower for high density. The other global property, the translational order parameter $Tr, shows a similar density dependence. As it describes translational and not orientational order, it reaches finite values only in the solid phase and is still close to zero for a hexatic structure. The value of $Tr for a solid is still smaller than j6 because translational order is fluctuating more strongly than orientational order. Due to the decaying long range order as a consequence of the Mermin-Wagner theorem [164–166], j6 and $Tr are smaller for 65 Measures of Ordering 4. Dynamics and Phase Behavior of Skyrmion Lattices larger numbers of particles and never reach 1 for finite temperature but increase with increasing density. 4.5. Analysis of Experimental Lattices One aim in performing experiments on skyrmion lattices is to create lattices with good ordering including skyrmion crystals with long-range translational order. While in previous publications the onset of a hexatic phase was already found [44], this structure still contained phase boundaries and was likely not fully equilibrated. To improve the understanding of skyrmion ordering, further experiments have been performed which build upon these results. A snapshot from an experiment with particularly good ordering is shown in Fig- ure 4.7. Here, one can already see good local and also reasonable global ordering in the field of view. Analysis of this lattice based on tracking results by Raphael Gruber can be performed with both, local and global parameters. As presented in sections 3.1.4 and 4.4, local order is described by %6 and 〈|k6 |〉 as well as radial and two-dimensional distribution functions. When excluding a boundary area of 20 µm at every side from the video shown in Figure 4.7, one obtains %6 = (87.1 ± 0.9)% and 〈|k6 |〉 = 0.819 ± 0.007. While this value of %6 is within the liquid-hexatic Figure 4.7.: Snapshot from a Kerr video of an experimental skyrmion lattice with hexatic ordering. White circles indicate skyrmions while the darker gray is the magnetic background. This experiment has been performed by Raphael Gruber on a sample created by Fabian Kammerbauer. 66 4. Dynamics and Phase Behavior of Skyrmion Lattices 4 3 2 1 0 0 5 10 15 20 Skyrmion-Skyrmion distance r [ m] Figure 4.8.: Left: Radial distribution function averaged over the full measurement of 5 minutes for the lattice shown in Figure 4.7. Right: Two dimensional radial distribution function of the same system with arbitrary scaling. coexistence region according to Figure 4.6, the value of 〈|k6 |〉 is in the range of values one would expect from simulations in the solid region. This difference in values to the ideal case in Section 4.4 is probably induced by pinning which leads to some local displacements of individual skyrmions. In addition, while patches of homogeneous orientation are larger than before, grain boundaries are still visible and influence results. Further analysis of the local structure can be performed using the radial distribution functions given in Figure 4.8. 6(A) shows a clear indication of at least hexatic order as the second peak is clearly split into two parts as it is an established feature for a two-dimensional solid with hexagonal order but not for a liquid [190] as already shown in Figure 2.9. This ordering is worse than for an ideal solid presented there as the splitting of the second peak is less pronounced. The two-dimensional distribution function 6(G, H) gives further insight into this. Compared to the lattice in Figure 3.6, this structure has a much stronger local ordering which is also preserved over a few layers beyond the first. This indicates a local order which is better than in an isotropic liquid. To gain more reliable results on the phase of the observed system, the global ordering of the system has to be considered. One quantity to do so is the angular susceptibility j6. Its value for the given system is 0.385 ± 0.023 which is a clear indication of at least hexatic order as this value is close to zero in liquids. Another quantity used to check if a system is in the hexatic phase is the radial distribution function 6(A) defined in Equation 2.101. As Figure 4.9 shows, this function decays slowly with distance A as  (A) ∼ A−0.056 . This indicates that this system has long-range orientational correlation and therefore at least hexatic order. As j6 and 6 only describe the orientation, they do not describe reliably if the system is 67 g(r) 4. Dynamics and Phase Behavior of Skyrmion Lattices 100 10 1 = 0.05 10 2 5 7 10 20 30 50 r [ m] Figure 4.9.: Orientational correlation function 6(A) averaged over the full mea- surement of 5 minutes for the lattice shown in Figure 4.7. The dashed line is a fit of 6(A) = 0A−[ against the upper envelope of this function. in the hexatic or solid phase. To check if the system is a solid, one needs to analyze the translational orientation of the system described by the translational order parameter $Tr and the translational correlation function @ (A). Using the same constraints as for the calculations above, one gets a translational order parameter of $Tr = 0.045 ± 0.016. As a value of exactly zero is improbable for finite systems, this is still small enough to assume a lack of translational order. This result can also be verified by the translational correlation function shown in Appendix A.6. As this function does not show a clear curve, performing a fit is difficult. A fit using only the highest peaks would give an exponent of 0.36 which indicates a non-solid structure as it is slightly above 1/3. Using different peaks would, however, yield a different number which could be larger or smaller than this. This ambiguity originates from the structure consisting of different domains without common orientation. Based on these results, the system is not in a homogeneous solid state spanning the whole sample. Looking at individual grains, it is, however, difficult to judge if these grains themselves are solid. Despite this ambiguity in assessing the solid order, it is clear that the system is at least in the hexatic phase as 6(A) indicates long-range orientational order. Additionally, local parameters also indicate at least hexatic order of the system. 68 G6(r) 5. Increasing Skyrmion Diffusion by External Stimuli 5. Increasing Skyrmion Diffusion by External Stimuli Skyrmion motion and especially thermal skyrmion diffusion is impeded by non-flat energy landscapes pinning skyrmions to specific positions in the sample. To reach faster and ideally free diffusion, one must overcome these pinning effects. To achieve this, different methods which lead to a significant increase in diffusion are described here. 5.1. Oscillating Magnetic Field One approach to increase skyrmion diffusion is applying an oscillating out-of-plane (OOP) magnetic field to the sample. This can lead to a diffusion increase by a factor of more than 300 as our published experimental results show [41]. The experiments and particle tracking have been done predominantly by Raphael Gruber without my contribution. I have performed analysis of occurrence maps, energy landscapes and displacement distributions as well as the analytical calculations presented in this section. 5.1.1. Experimental Results While the effect of static fields on skyrmion diffusion and ordering is already quite well known [22, 44, 189], effects of oscillating fields are so far less researched. Here, it is shown that using an oscillating field can lead to a huge increase of the diffusion constant by a factor of around 330. This increase was found for a field variation of 2.1 mT at a frequency of 2 kHz as indicated in Figure 5.1. A higher time resolution of the MOKE could even improve this diffusion increase as the currently used setup and trajectory linking limit the diffusion that can be recorded. Beyond this, an even higher increase in diffusion is prevented by skyrmion annihilation at too high fields and lower frequencies. In general, a higher amplitude of the alternating field leads to a larger increase in diffusion. Besides this effect one can see a strong frequency-dependence in diffusion for all choices of the magnetic field. For very high frequencies above 5 kHz, there is next to no effect of any kind of oscillation and also for very small frequencies the effect is very weak. In between, the diffusion peaks at 5 = 10 ± 2 Hz. These effects go beyond what a simple size-dependence of the diffusion constant as described in Figure 5.2 could explain. This figure shows the dependence of skyrmion size as well as the diffusion constant on the magnetic 69 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.1.: Experimentally determined skyrmion diffusion constant  depending on the frequency 5 of an applied sinus-shaped magnetic field variation around a field of 0 = 27 µT for different values of the peak-to-peak amplitude of the field variation . For excitations with low frequencies and high amplitudes skyrmions start to annihilate leading to the gray shaded area being inaccessible experimentally. This figure is taken from our publication [41]. field for the static case without an alternating field. While an increasing magnetic field leads to a reduced skyrmion size, the diffusion constant increases with an increasing field. Combining these two results indicates that smaller skyrmions diffuse faster but the difference in diffusion between small and large skyrmions is only by a factor of around 1.5. Besides the different order of magnitude in diffusion changes for the static and the excited case one also notices that skyrmions are only stable in a relatively narrow field region of around 25 µT without excitations as given in Figure 5.2. In the excited case, however, a field region nearly 100 times as wide in which skyrmions are stable with fast enough excitation is found. This result also shows that just using the static description does not suffice to explain the increase in diffusion for periodically excited skyrmions. These differences between static and alternating fields as well as the peak of the diffusion constant around 10 Hz can be explained by skyrmion pinning where an inhomogeneous energy landscape leads to some positions being preferably occupied by skyrmions and impeding their diffusion. This leads to an MSD of skyrmions which is several orders of magnitude lower than it would be on a flat energy landscape. As described in Section 2.4, Raphael Gruber et al. showed that pinning effects mostly interact with the domain wall at the skyrmion border and not with 70 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.2.: (a) The skyrmion radius depending on the applied magnetic field . (b) The experimentally determined diffusion coefficient of skyrmions  depending on the magnetic field . This general result is not true for the brown areas of the figure. For very high field, the lifetime of free skyrmions is in the scale of several seconds leading to an overrepresentation of strongly pinned skyrmions whose average lifetime is longer and therefore to lower diffusion. For very low field, worm domains begin to form and therefore skyrmions are no longer spherical. In this case, the tracking algorithm which is focused on spherical objects cannot reliably track the real positions of skyrmions anymore and therefore introduces artificial noise to the real diffusion leading to an overestimation of . The diffusion constant is fitted to = 2 0− giving the dotted line with 20 = 8.91 µm 2/(µT s) and 2 21  1 = 114.33 µT. The derivation of this formula is explained below. This figure is taken from our publication [41]. the skyrmion core [67]. Therefore, a pinning site which is favorable for a skyrmion of given size by pinning several parts of its domain wall is often no longer favorable for a smaller or larger skyrmion. For this reason, during the process of changing the external field and thereby the skyrmion size, a skyrmion is subject to many different energy landscapes instead of one constant landscape. A position which minimizes the energy at a given size and therefore prevents diffusion likely does not have the same effect on a smaller or larger skyrmion and therefore the diffusion of size-oscillating skyrmions is much less impacted by pinning. This effect is largest for very large changes in the magnetic field because this also leads to the largest changes of skyrmion sizes and in the effective pinning landscape. Also the limit for high frequencies can be explained by pinning effects, because the time for which skyrmions change their size is not long enough to significantly diffuse and thereby leave a given pinning site. For very low frequencies, the timescale on which a skyrmion could leave a pinning site due to size changes is similar to the timescale on which it leaves pinning sites due to random thermal excitations and therefore the diffusion is not significantly increased here. In between these limits, the skyrmions 71 5. Increasing Skyrmion Diffusion by External Stimuli can leave a pinning site aided by the oscillating field and move much closer to free diffusion than before. As an example, a skyrmion increasing its size leaves a pinning site which is strong for small skyrmions to move to a strong pinning site for large skyrmions. In the next period when the skyrmion becomes smaller again, it can now easily leave the pinning site for large skyrmions and can diffuse further until it finds another pinning site. This increases the diffusion coefficient way beyond the value determined with fixed pinning sites. 5.1.2. Theoretical Calculations Before analyzing the increase in skyrmion diffusion further using e.g. trajectories and displacement distributions, I perform some calculations to analyze the size dependence of diffusion. These calculations demonstrate that given the employed approximations the size dependence as in Figure 5.2 is insufficient to explain the large increase in diffusion found for this system. As this is a simplified model, effects of pinning are not considered here. I start by calculating the dissipation tensor of a skyrmion which is an integral part of the calculation of the skyrmion diffusion constant. Here I assume Néel skyrmions with a tanh-shaped domain wall [55] as this shape minimizes the e√nergy of the domain wall:©cos(i)√1 − tanh2( A−0 )1 ª <® = ­­­ ®«sin(i) 1 − tanh2( A−0 )®® , (5.1)1tanh( A−0 )1 ¬ where 0 is the radius of the skyrmion (i.e. the distance from the skyrmion center to the point in the domain wall where the I-magnetization is zero) and 1 is a constant proportional to the domain wall width given in Equation 2.37: X = c1. Using this spin structure and the definition of the dissipation tensor in Equation 2.66, one can calculate the diagonal element of the tensor with few approximations. Due to the length of this calculation, it is given in Appendix A.2.2. Doing the calculation and using some simple approximation∫s, one arrives at GG = HH ∼ (m (®G )2 2c0 32A ≈ . (5.2) 1 This result for the dissipation tensor can now be used to calculate the diffusion constant. In the case of large skyrmions with narrow domain walls, one can assume UT   and therefore neglect the Magnus force for a rough estimate of the 72 5. Increasing Skyrmion Diffusion by External Stimuli diffusion. Thereby one obtains by plugging this into Equation 2.83 U = GG :B) 1  :B) (U )2 + 2 ≈ ∼ . (5.3)GG  UGG 0 Mathematically, this relation only holds true in a case of 0  1, i.e. where the skyrmion radius is much larger than the width of the domain wall. However, for experimentally relevant skyrmions this is not an issue as micromagnetic simulations show that 0  1 for all fields under which skyrmions are stable in the experiments presented in this thesis. Assuming that the skyrmion size depends linearly on the applied field and that the domain wall width is constant, which is in good agreement with Figure 5.2(a), one can write the dependence of the diffusion constant on the field as 2 = 0 21 − , (5.4)  with two numerical constants 20 and 21 which have to be determined by a fit to the data given in Figure 5.2(b). The fit is able to reproduce the data in the stable region quite well, indicating that the theoretically determined relation holds true in this case with a reduced j2 of 0.67. In the experimentally accessible region with constant parameters one can therefore modify both size and diffusion by up to 40%. A higher increase in diffusion requires additional effects besides the described size effect. Going beyond these static effects, one can also take a look at dynamic effects of alternating perpendicular magnetic fields. At first, the effect of an alternating magnetic field compared to a static field is analyzed. Here, the field dependency of the diffusion constant established in Equation 5.4 and a sinusoidal variation around the field 0 are used: 1 ∫ ) 20 2= 0 3C 1 ( ) = √ , (5.5)) 0 2 −  −  sin 2cC1 0 ( )2 ) (21 − 2 )2 − 10 2 where  is defined as the peak-to-peak variation of the field leading to the factor 1/2 in the equation. The result shows that an alternating field always increases the diffusion compared to the constant field case  = 0. Additionally, in this ideal case the diffusion is independent from the frequency of the oscillation. This assumption only holds in regions where the frequency is low enough to allow for the skyrmion size to follow the excitation. Outside the stable region the applicability of this model is limited as Equation 5.5 would indicate an imaginary and thus unphysical value of the diffusion for a too large field variation. This happens because when assuming a linear dependence of skyrmion size on the magnetic 73 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.3.: Example trajectories of skyrmions with their initial position placed at the coordinate origin without excitation (a) and in the case of highest observed diffusion with an applied field change of 2.1 mT at a frequency of 2 kHz (b). The background shows a Kerr image of one skyrmion placed at the coordinate origin and its surrounding area. Note the larger scale on the axes for excited skyrmions. This figure is taken from our publication [41]. field, a low enough field would lead to skyrmions with negative diameter which is physically not possible. Due to this limitation, the case of a field variation of  = 30 µT around 0 = 27 µT is chosen to compare the experimental results to the theoretical calculation. This case is close enough to the experimentally determined sizes and diffusion constants to rely on the model derived above. For this excitation one obtains the highest diffusion speedup of around 300% for 5 = 10 Hz in experiments, which is two orders of magnitude lower than the maximum speedup found for higher fields. When using the theoretical formula for this values, the expected increase due to size oscillations is however only 1.6%, more than two orders of magnitude smaller than the measured increase. Also when just assuming that the skyrmion would always diffuse at the speed it has at the smallest possible field of  − 2 = 12 µT, this would lead to an increase of only 31%. These results show that the increase in skyrmion diffusion by applying a field is orders of magnitude larger than what could be explained by size-dependent free diffusion. 74 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.4.: Logarithmically plotted histograms of skyrmion occurrences for the unexcited case (a) and an excitation with  = 2.1 mT and 5 = 2 kHz (b). Both histograms are normalized in a way that the same total number of positions is used to determine the histogram. White spaces indicate zero occurrence in the analyzed frames. The bin width of the histogram is 1 µm. This figure is taken from our publication [41]. 5.1.3. Effective Reduction of Pinning To further analyze how pinning affects the diffusion of skyrmions with and without periodic excitations, one can check different other properties of the system such as individual trajectories, skyrmion occurrence profiles or distributions of displace- ments to gain more information on the kinds of diffusion the skyrmion exhibits. The trajectories in Figure 5.3 for the unexcited case indicate a hopping-like motion on longer length scales while between the jumps a skyrmion stays close to a specific position and is only displaced slightly. The diffusion in the excited case is closer to normal Brownian diffusion and spans a much larger area. Pinning regions are still visible in this case but the time for which a skyrmion stays in one is significantly reduced. One can conclude that with excitations the diffusion is less hopping-like and approaches free diffusion but does not reach the ideal free diffusion limit. The favoring of specific positions in which skyrmions are pinned can also be ob- served in a histogram of skyrmion occurrences given in Figure 5.4. Without an excitation skyrmions tend to stay at specific positions and are not able to explore the whole sample as the large white regions with zero skyrmion occurrence in the picture indicate. When applying an alternating field, these unexplored regions are significantly less frequent as skyrmions are able to explore the whole sample and are much less bound to specific positions or pinning sites. This difference in the occurrence maps is another clear indication that applying field excitations leads to a strongly reduced effective pinning. The occurrence maps can also be used to calculate effective energy profiles of the 75 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.5.: Logarithmically plotted histogram of skyrmion occurrence for the unexcited case using all available data from 225 minutes of measurements. This figure is taken from our publication [41]. Figure 5.6.: Effective energy landscapes of skyrmion pinning calculated using Equation 3.6 at a constant field (a) and with a field excitation of  = 2.1 mT and 5 = 2 kHz (b). The potential at bins with zero occurrences is set to  = 0. This figure is taken from our publication [41]. samples using Equation 3.6. To also describe regions with low skyrmion occurrence reliably, all available statistics in the unexcited case, which is more than 20 times more data as used for the excited case in Figure 5.4, is used here as displayed in Figure 5.5. The energy landscapes calculated from this data in Figure 5.6 show a result very similar to the occurrence maps. Pinning sites in the unexcited case are significantly deeper than in the excited case indicating that the effect of pinning on skyrmions with size oscillations is diminished. Due to the Arrhenius-like dynamics of particles in pinning sites, the probability to leave a site is proportional to the negative exponential of the depth leading to slight changes in the potential depth already having a strong effect on dynamics. Adding to these qualitative results one can also directly analyze the displacements of skyrmions after a fixed time to understand their modes of movement. For free 76 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.7.: Histograms of skyrmion displacement after ΔC = 2 s determined as a running average over longer experiments for the excited case (a, b) with 5 = 2 kHz and  = 2100 `T and for the unexcited case (c, d). (a) and (c) show the displacement in x-direction which is an arbitrarily chosen direction and has roughly the same shape as the displacement e.g. in y-direction. (b) and (d) show the absolute displacement. Displacements in x-direction are fitted to a sum of two Gaussians indicated with A and B in (a). Absolute displacements are fitted to a sum of two Rayleigh distributions. Note the significantly smaller displacements for the unexcited case. This figure is taken from our publication [41]. diffusion, displacement distributions in one dimension are Gaussian distributed. Summing the diffusion in two dimensions this leads to a Rayleigh distribution of displacements when adding up the parts in (x- and y)-direction:2 ?(ΔG) ∼ exp −(ΔG)ΔG 2 , (5.6)f2 77 5. Increasing Skyrmion Diffusion by External Stimuli where f is the same as f of the Gaussian distributions in x and y. The derivation of this distribution is equivalent to the Maxwell-Boltzmann distribution commonly known for three-dimensional systems. Displacements as used here follow the same probability distribution as instantaneous velocities, which are commonly used for this kind of analysis, but experimentally inaccessible here due to the overdamped motion of skyrmions and too low time-resolution of the MOKE setup. In contrast to the theory for free diffusion, the results in Figure 5.7 show a very good separation into two different Gaussians for the excited x-distribution. The inner Gaussian (A) describes the movement of pinned skyrmions inside the pinning site while the outer Gaussian (B), which contains significantly more trajectories, describes freely diffusing skyrmions. The same can also be observed in the Rayleigh distributions where the left peak is again the movement of pinned skyrmions and the right one is the movement of freely diffusing skyrmions. Using two Rayleigh distributions here builds on the assumption that a pinned skyrmion is pinned in both, G and H while an unpinned skyrmion is found in the unpinned distribution for G and H. Besides these two peaks, another smaller peak is found for an absolute dis- placement of around 4 µm. Different from the two other peaks, the position of this peak is also conserved when changing the lag time ΔC indicating an origin beyond dynamic effects. The most plausible cause for this effect is the size of skyrmions in the Kerr videos being around 4 µm (the size determined by the tracking routine and shown in Figure 5.2 is slightly smaller than the size one would determine by eye) and the peak therefore describing a displacement by one skyrmion diameter. Displacements of exactly one diameter are caused by the domain wall pinning which can pin either side of a skyrmion so that the skyrmion center moves by one diameter when different sides of the skyrmion are pinned. The distribution without an alternating field shows much smaller displacements on the same timescale and also a much smaller separation between the two different Gaussians used in the fit. The double Rayleigh distribution cannot be fitted here well indicating that there is no clear separation into two kinds of diffusion as before. Instead, all skyrmions exhibit pinning-impeded diffusion which can cover different distances based on the size and depth of individual pinning sites. Beyond the general shape of the displacement, one notices a small shift of both x-displacements to negative values. This is probably caused by a small field gradient present in the sample. Furthermore, one needs to consider that for ΔC → ∞ both distributions are expected to exhibit a single-Gaussian shape due to the central limit theorem applying in this case where ΔC is much larger than the time a skyrmion stays in one pinning site. This effect is displayed for the excited case in Appendix A.7. Here, the diffusion would however still be much larger for the excited case. On the other hand, for ΔC → 0, diffusion inside and outside pinning sites happens on very similar length 78 5. Increasing Skyrmion Diffusion by External Stimuli scales which are shorter than a pinning-site diameter. In this case, the two different kinds of motion cannot be clearly distinguished anymore as shown in Appendix A.8. For this reason, a timescale in which one sees a clear difference between the two different kinds of diffusion is chosen. In the case of free diffusion, diffusion would be described by only one Gaussian displacement distribution in one dimension or one Rayleigh distribution in two dimensions. This peak is expected to be similar to the wider Gaussian in Fig- ure 5.7(a) while the narrow peak around a displacement of zero disappears for free diffusion. Additionally, the shape of displacement distributions for free diffusion does not depend on ΔC. In conclusion, one can state that pinning is not completely removed by applying an alternating field but instead it is effectively reduced so that unpinned skyrmion dynamics dominates over pinned dynamics. The good agreement of the fit functions with the data in the excited case indicates that the model of two different regimes is indeed appropriate to describe the motion. 5.2. Alternating Current Besides changing the magnetic field and thereby changing the effective pinning landscape for skyrmions due to a change of the skyrmion radius, diffusion of skyrmions can also be accelerated by using currents to push skyrmions out of pinning sites. As the goal of this study is mitigating pinning effects and not directed motion, a current excitation with a mean current of zero is used. This can be achieved by using a periodic function, which contains negative and positive currents of the same amplitude. The investigations presented here are motivated by preliminary experiments conducted by Tobias Sparmann which are not published, yet. As these experiments were using a sinus-shaped current excitation, computer simulations were also performed with a sinus excitation  (C) = 0 sin(2c 5 C). The main results of these experiments are given in Appendix A.9. As expected an applied oscillating current leads to increased diffusion of skyrmions. While this effect is quite small for high frequencies, diffusion is increased significantly for smaller frequencies of the sinus excitation. This relation between frequency and diffusion can be fitted to a falling power law  ∼ 5 −U with U > 0. To investigate this behavior, simulations were performed for a toy system with regularly distributed pinning sites and no skyrmion-skyrmion interactions as well as for interacting skyrmions on the pinning landscape shown in Figure 3.9. In the most simple model, the only relevant interactions are the interactions between skyrmions and pinning sites and between skyrmions and the external force. As the acceleration of skyrmions only works with pinning effects, these two interactions are required for a minimal model. Periodic boundary conditions were applied. 79 5. Increasing Skyrmion Diffusion by External Stimuli 1.0 No Excitation 0.7 Excited 0.4 0.2 0.1 100 101 102 103 Frequency [1/ t] Figure 5.8.: Diffusion constant of skyrmions with an applied sinus force with fixed amplitude and varying frequency  (C) = 0 sin(2c 5 C). The red line shows the diffusion without an applied force. Errors are smaller than the size of the points. Skyrmion-skyrmion interactions were not considered so the simulation effectively emulates the behavior of single skyrmions in the presence of pinning sites. Here, the latter are placed on a periodic square lattice and modeled with a harmonic pinning potential as given in Equation 3.7 with : = 10:B) and 0 = 1 µm. The distance between pinning sites is 6 µm. This pinning strength and distance is similar to experimental samples used in different parts of this thesis. The applied force is sinusoidal leading to a net force of zero, has a constant amplitude and a varying frequency. The amplitude of the sinus is set so that skyrmions are moved by an absolute value of 25 µm per second, which is similar to the distance moved by the force in the preliminary experiments. The simulated diffusion coefficient for different excitation frequencies in this system is given in Figure 5.8. It can be described in three different frequency regimes. For very high frequencies beyond 100ΔC−1 the current excitation has no visible effect on the diffusion. In this regime the current flowing in a given direction does not suffice to allow the skyrmion to leave a pinning site during one period. Here, the diffusion is therefore not altered compared to a case with no applied current and skyrmions stay fully pinned. The second region is found between 100 and around 1ΔC−1. In this region, lowering the excitation frequency leads to a significant increase of the diffusion constant of the system. Here, skyrmions are able to sometimes leave their pinning sites with aid of the applied current. However, the force is not applied long enough to allow for this movement every time and therefore the diffusion is still lower than free diffusion and pinning effects are visible. At around 1ΔC−1 the system reaches a diffusion constant close to the unpinned diffusion of 80 Diffusion Constant [ m2/ t] 5. Increasing Skyrmion Diffusion by External Stimuli = W :B) 2+ 2 = 0.94µm 2/ΔC [154] after which it does not increase further. In W  agreement with the experiment, one observes no increased skyrmion diffusion for high frequencies. In the intermediate area, simulations do not show a linear drop in the logarithmic plot as it was found in the experiment. This could be the case either due to the experiment only containing a small part of the rising slope which can be approximated by a linear function or due to systematic differences between the simulated and the experimental system such as very deep or large pinning sites keeping skyrmions for even longer. To learn which of these explanations is the correct one, one would have to perform experiments for lower frequencies to see if the slope decreases and the diffusion reaches a plateau as it does in the simulation. This, however, is difficult due to the limited size of the experimental setup, as a lower frequency would correspond to a larger distance being traversed by one skyrmion during one excitation. Still, some additional simulations can be performed to show that the increase in diffusion and also the plateau for low frequencies are robust. As the explanation above describes the increase in diffusion when applying a current purely as a consequence of mitigating pinning effects, one would expect that the same current excitation does not impact diffusion in a system without pinning. Simulations to investigate this effect are presented in Appendix A.10. This result shows that there is indeed no effect caused by current in this case as the calculated diffusion constants are within the theory value considering their error bars. While the results above describe the behavior of individual skyrmions, many experiments use skyrmions at a finite density. To simulate this effect, simulations with skyrmion-skyrmion interactions as determined in Section 3.1 at a density of 0.00833 µm−2 are performed. At this density skyrmions do not form a lattice but skyrmion-skyrmion interactions already play a significant role for dynamics. While this addition leads to overall reduced diffusion due to lattice effects, it is assumed to not change the shape of the curve. The frequency dependence of the diffusion constant in this case is given in Appendix A.11. As expected, the shape is very similar to the case without skyrmion-skyrmion interactions. The amplitude is smaller over all values of the frequency by a factor of around 1.6 as a consequence of the skyrmion-skyrmion interaction limiting the mobility of individual skyrmions. Further analysis of the effect of applied currents can be performed by using a realistic experimental pinning lattice for the same kind of simulation. To achieve this, the lattice described in Figure 3.9 and the same interactions and density as above are used. The overall shape of the curve given in Figure 5.9 is very similar to previous results. However, the less uniform pinning present in the landscape leads to some differences. As pinning sites are larger compared to the model with ideal sites, longer excitations are required for skyrmions to leave their individual positions. This leads to a smaller frequency of around 50 ΔC−1 above which current excitations 81 5. Increasing Skyrmion Diffusion by External Stimuli 0.6 No Excitation 0.4 Excited 0.3 0.2 0.1 100 101 102 103 Frequency [1/ t] Figure 5.9.: Diffusion constant of skyrmions with an applied sinus force with fixed amplitude and varying frequency  (C) = 0 sin(2c 5 C) at density 0.00833 µm−2 on the energy landscape given in Figure 3.9. The red line shows the diffusion without an applied force. Errors are smaller than the size of the points. do not lead to increased diffusion. Additionally, the minimum amount of diffusion found here is significantly larger than in the case with a lattice of pinning sites. This is caused by the pinning landscape used here not only containing deep pinning sites but also areas in which these sites are shallower. Therefore skyrmion diffusion in these regional is less impeded by pinning effects. The maximum diffusion constant at low frequency is, however, similar to the case with ideal pinning sites as lattice effects still reduce diffusion in the same way as before. To further analyze the effect of applied currents on skyrmion dynamics methods similar to what was done for applied oscillating fields are employed. As occurrence maps are especially suitable for analysis of experiments, this analysis is focused at displacement distributions which can be analyzed in a similar fashion as in Figure 5.7. As the diffusion within full periods of the sinus does not depend on the direction, the y-diffusion has the same distribution as the x-diffusion for all three plots and is therefore not presented here. To exclude lattice effects, the case without skyrmion-skyrmion interaction and two different excitations as well as the unexcited case are considered. In the unexcited case in Figure 5.10(a) one can see that the diffusion happens on very small length scales but has a mostly Gaussian shape in 1D (and Rayleigh shape in 2D). This is caused by pinning effects, which lead to nearly all skyrmions remaining in their individual pinning sites and diffusing inside those. Adding a periodic excitation of 5 = 10ΔC−1 leads to a qualitative change in diffusion as Figure 5.10(b) shows. Instead of a shape which could be fitted to a single Gaussian this displacement distribution is clearly 82 Diffusion Constant [ m2/ t] 5. Increasing Skyrmion Diffusion by External Stimuli Figure 5.10.: Displacements of skyrmions after 1ΔC in x-direction fitted to a Gaussian (left) and absolute displacements fitted to a Rayleigh function (right). Results for no excitation (a), 5 = 10ΔC−1 (b) and 5 = 1ΔC−1 (c). (b) is fitted to a sum of two Gaussians or Rayleigh functions instead of one, respectively. 83 5. Increasing Skyrmion Diffusion by External Stimuli consisting of two parts - pinned and non-pinned diffusion. This splitting is easily visible as the diffusion in x-direction can be fitted reliably by two Gaussians - one for pinned and one for unpinned diffusion - while the absolute diffusion is fitted to a sum of two Rayleigh functions. Increasing the diffusion further by reducing the frequency of the sinus to 5 = 1ΔC−1 leads to a simpler shape again as given in Figure 5.10(c). In this regime, nearly all skyrmions are able to move in an unpinned fashion and therefore the wider Gaussian from the previous case is now describing the behavior of skyrmions. However, some pinning effects are still present as the slightly increased height of the peak in x-diffusion indicates. In summary, one sees a very similar effect as shown for applied currents in Figure 5.7. A strong enough excitation leads to a kind of diffusion, which is dominated by free-like diffusion and for this reason reaches a Gaussian-like shape. In this situation, pinning effects only play a minor role. 5.3. Comparison of the Methods to Increase Diffusion Above, I presented two different approaches which reliably allow for skyrmions to diffuse in an unpinned fashion even on a strongly pinned sample. Besides this similarity, they however differ strongly in how they affect the skyrmions. This also affects how they can be implemented for e.g. skyrmion-based computing approaches or experiments on formation of skyrmion lattices. The major advantage of using out-of-plane magnetic fields to increase diffusion is that skyrmions are not subjected to directed movement during the excitation. This makes it possible to use field excitation in tight confinements. Current excitations would not be suitable under such circumstances as the current applied during half a period could push the skyrmions into a material boundary and therefore lead to skyrmion annihilation. Field variations with a spatially uniform out-of plane field on the other hand do only change the skyrmion size within a certain range. Skyrmions can keep their individual positions and are not pushed into structures that would annihilate them such as material boundaries. Another advantage of magnetic field excitations is that they can possibly lead to diffusion beyond the free diffusion. Changing the relevant energy landscapes could lead to pinning sites pulling skyrmions and making them move faster than they would normally. However, estimating the exact strength of this effect is quite difficult. Using current excitations, on the other hand, also comes with some advantages. When using a nearly 1D channel in which a skyrmion moves, increasing the size could become obstructed by the confinement while a current along the channel would not be impacted. In this case, current excitations could lead to a larger increase in diffusion. Additionally, modeling of current effects via computer simulations is much simpler as the effective energy landscape is calculated by superposition of the measured 84 5. Increasing Skyrmion Diffusion by External Stimuli pinning landscape and a pushing force. For field excitations one would on the other hand need landscapes for all different skyrmions sizes to be able to reliably model the system. Additionally, as long as currents are not too large, skyrmion shapes are not changed by current excitations. This could be relevant e.g. for skyrmion ordering and lattice formation. Both methods share the property that by increasing the diffusion, they are able to both speed up ordering of skyrmions as well as destroy an ordered skyrmion lattice. This is comparable to increasing the temperature in other kinds of systems. For this use-case however, it is not obvious, which of the two effects would be more reliable as long as the system is large enough that annihilation at boundaries does not play a huge role. Obtaining insights into this effect would be an interesting subject for new research. 85 6. Skyrmion Dynamics in Confinements 6. Skyrmion Dynamics in Confinements Many skyrmion-based applications in computing require skyrmions constrained to some confinement. To improve the understanding of such confinement effects, experiments and computer simulations of skyrmions in specific confinements have been performed. Here, static structures of skyrmions in confinement as well as the dynamics are analyzed. The bulk of these results has already been published [110]. Additionally, the effects of applied currents on this ordering and dynamics are analyzed to improve the understanding of manipulation of confined skyrmions. These results are part of two publications [99, 111] which are not presented in this thesis in their entirety as I contributed to specific parts of these works only. 6.1. Confinements Without Excitations Skyrmion motion is strongly influenced by confinement effects. This has already been shown in different micromagnetic simulations of e.g. triangular confine- ments [54, 65]. Extending on these results, experiments in micrometer-scale experimental setups with micrometer-sized skyrmions were performed mostly by Chengkun Song without my participation [110]. To better understand the experi- mental results, particle-based Brownian dynamics simulations as described in the theory chapter are performed using an exponential potential with constants deter- mined in micromagnetic simulations by Schäffer et al. [54] for skyrmion-skyrmion and skyrmion-wall interactions of nanometer-size skyrmions: + (A) = 04−A/1 (6.1) Parameters of this potential are 0SkSk = 1.246 eV and 1SkSk = 1.176 nm for skyrmion- skyrmion and 0SkBnd = 0.211 eV and 1SkBnd = 1.343 for skyrmion-boundary interactions [54, 110]. These potentials describe skyrmions in the nanometer scale which is a distinctly different length scale than in the experiment. This is, however, not a big issue as these simulations are not focused on reproducing quantitative results but on reproducing commensurability effects qualitatively in experiment and simulation. The most important result of this investigation is how the dynamics of skyrmions are not only affected by the density but also by commensurability effects. While an increase in density leads to less diffusive motion of skyrmions, commensurability 86 6. Skyrmion Dynamics in Confinements Figure 6.1.: (a) MOKE image of two, five and eight skyrmions in a circle of radius A = 8.4 µm averaged over 10 minutes. Bright areas are regions where skyrmions are likely to be located while dark regions are on average less populated by skyrmions. (b) Skyrmion trajectories corresponding to the MOKE images shown above. Individual skyrmions are encoded by different colors. (c) Normalized magnetic field for different distances A from the disk center. The signal is averaged over 10 minutes and all angular positions. This figure is taken from our publication [110]. leads to different dynamics for skyrmion numbers commensurate or incommensurate with the sample geometry. 6.1.1. Circular Confinement A circle is the geometrically simplest two-dimensional confinement, because it has maximum symmetry with infinite mirror axes and no preferred directions. This leads to a rotationally invariant ordering of skyrmions both in simulation and experiment. Also, there is no specific number of skyrmions which is either commensurate or incommensurate with the system as circles can form containing any number of particles. To analyze the effects of a circular confinement, experiments were performed for a circle of 8.4 µm containing two, five and eight micrometer-sized skyrmions. For 87 6. Skyrmion Dynamics in Confinements 1 2 3 4 5 6 7 8 9 10 Figure 6.2.: Histograms of one to ten skyrmions in circular confinement of radius 10 nm. Skyrmions are displayed as white circles with a radius equal to the skyrmion- skyrmion interaction width 1. Brightness is normalized in the same way for every picture. All structures are radially symmetric. two and five skyrmions a circle without a center skyrmion forms while for eight skyrmions one skyrmion is placed in the center as shown in the histograms and trajectories in Figure 6.1(a,b). The configuration with two skyrmions allows for significantly more diffusion of individual skyrmions as it is indicated by both trajectories spanning the all 360°of the circle during the recorded measurement time of 10 minutes while for eight skyrmions the same measurement time only allows for each skyrmion moving in a small part of the full circle. Also, no position swaps between the center skyrmion and one of the outer skyrmions happened during the measurement time. The non-isotropic occupation of different positions at the same radial distance from the center is caused by pinning effects in the sample which lead to some positions preferably occupied by skyrmions. Figure 6.1(c) shows the averaged brightness as a function of radial position. The peaks in this plot are caused by the ring of skyrmions and in the eight-skyrmion configuration by the center skyrmion. The constant relative value of 1 close to the center for eight skyrmions indicates that the center skyrmion never leaves the central position during the experiment while there is a negatively magnetized area between the skyrmions within the ring. Beyond the experimental results presented above, computer simulations have been performed for a circle of radius 10 nm containing one to ten skyrmions. This radius of the circle leads to a similar ratio between skyrmion size and confinement radius as found in the experiment as simulated skyrmions have a diameter of a few nanometers. As shown in Figure 6.2, a circle without a skyrmion in the center 88 6. Skyrmion Dynamics in Confinements 1 6 2 7 1 3 810 4 9 5 10 100 10 1 100 101 102 t Figure 6.3.: Short-time MSD for different numbers of skyrmions in a circle in a computer simulation calculated using the windows-sliding method. The MSDs are averaged over all skyrmions. ΔC is given in simulation units. forms for two to five skyrmions while for six to nine skyrmions one skyrmion is placed at the center and the remaining ones form a circle again. For ten skyrmions, the state with minimal energy contains two skyrmions in an inner circle and eight in an outer circle. As another way to gain a deeper understanding of skyrmion dynamics in confinement, the short-time mean-square displacement (MSD) of different numbers of skyrmions is investigated. Without geometrical effects it is expected that a higher skyrmion density leads to less diffusion and therefore a lower MSD. Figure 6.3 however shows a significantly lower saturation MSD for one skyrmion compared to larger numbers after some time. This effect is caused by the energy-minimizing position of one skyrmion in the center region of the disk. For larger numbers of skyrmions there is no single ideal position which leads to collective movement of skyrmions in a circle. The same effect is visible in the case of ten skyrmions compared to eight or nine skyrmions where two skyrmions are placed in the center instead of one and therefore there is no skyrmion with a fixed position anymore. Due to the high symmetry of the circle, these effects are very small and the density dependence of diffusion is dominant. To support the hypothesis that the crossing of MSDs for ten skyrmions with MSDs for lower numbers is indeed a consequence of the different behavior of the skyrmions in the center, the MSD of only the outer skyrmions can be considered. To achieve this, all skyrmions positions which are placed less than 3 nm from the center are 89 MSD [nm2] 6. Skyrmion Dynamics in Confinements omitted and only parts of trajectories outside this area are used. The result in Appendix A.12 shows that now all MSDs are in the inverse order of the skyrmion numbers and the MSD of ten skyrmions is no longer an exception. Beyond this, one sees that the MSD for five and six skyrmion as well as for nine and ten skyrmions, respectively, is nearly the same. This is a consequence of these states containing the same number of skyrmions in the outer ring and only differentiating themselves by the skyrmion in the center. 6.1.2. Triangular Confinement In comparison to a circle, equilateral triangles exhibit much lower symmetry with only three symmetry axes. As a consequence, a triangular confinement does not allow for a collective rotation of all skyrmions. This leads to skyrmions having strongly preferred positions in two dimensions instead of just a preferred radial position. To investigate this behavior, experiments are performed in a triangle with side length 26.8 µm while simulations are performed using a triangle with side length 30 nm with skyrmion sizes as in the circle in both cases. This leads to a comparable ratio between confinement size and skyrmion radius in experiment and simulation. Due to the symmetry of the confinements, skyrmions show a lattice-like ordering as in Figure 6.4. The skyrmions form a regular structure inside the triangle which is mirror-symmetric along the three symmetry axes of the triangle in simulations and nearly symmetric in experiments containing pinning. A triangular ordering is found for three, six and ten skyrmions in agreement with the confinement both in experiment and simulation. Structures of nine skyrmions occupy similar positions as ten skyrmions, but because of the missing skyrmion one lattice site is left empty and skyrmions can move between different lattice sites. This effect also occurs for five skyrmions in the simulation while structures with four, seven or eight skyrmions form somewhat different configurations. The largest difference between simulation and experiment is visible for two, four and five skyrmions where pinning effects have a large impact on the experimentally determined structure, leading to an asymmetric placement of skyrmions compared to the triangular symmetric case observed in simulations. Pinning-induced dif- ferences are also visible for seven and eight skyrmions where some skyrmions are placed at the same positions as in the experiment with ten skyrmions while the simulated structures are structurally different. In general these results show that a commensurate order of skyrmions is able to overcome pinning in this sample while the simulated ideal ordering of non-commensurate numbers is overshadowed by pinning effects. This finding is in agreement with the finding that skyrmion lattices are less affected by pinning than individual skyrmions. Figure 6.5 contains the MSD of the different skyrmion configurations in experiment 90 6. Skyrmion Dynamics in Confinements Figure 6.4.: (a) MOKE pictures of one to ten skyrmions in a triangular confinement of side length 26.8 µm. The top row displays initial configurations just after nucleation. The bottom row displays configurations time-averaged over 10 minutes. (b) Time-averaged histograms of skyrmion positions generated from simulations of a 30 nm triangle using positions from 1 million frames. Skyrmions are displayed as white circles with a radius equal to the interaction width 1 of the exponential potential. This figure is taken from our publication [110]. and simulation. Commensurate states of three, six, and ten skyrmions have a significantly lower peak value of the MSD compared to non-commensurate states which is a distinct difference to the results for the higher-symmetry circle. The one-skyrmion state also reaches a plateau value due to its confinement to center of the triangle similar to the one-skyrmion state in the circle. Analysis of the simu- lated result for the four-skyrmion MSD is more difficult, because the four-skyrmion state with one skyrmion in the middle is meta-stable but still allows for switching between positions after some time leading to the MSD incorporating two timescales of motion and a lower MSD in between. This effect is not visible in the experiment where only the square-like state occurs frequently. Other skyrmion numbers are ordered in a way indicating that the MSD is lower for systems with higher skyrmion density both in experiment and theory as already shown for the circle. This effect is visible when comparing non-commensurate states to each other as well as when comparing commensurate states for different skyrmion numbers. For this reason the highest MSD is achieved for two skyrmions both in experiment and simulation. Using these results one can choose different numbers of skyrmions for applications depending on the strength of diffusion and the amount of possible states required for a specific application. 91 6. Skyrmion Dynamics in Confinements Figure 6.5.: Short-time MSDs for one to ten skyrmions in a triangle in (a) experiment with radius 26.8 µm and (b) simulation with radius 30 nm calculated using the windows-sliding method. The MSDs are averaged over all skyrmions. This figure is taken with minor adaptations from our publication [110]. 6.1.3. Square Confinement The confinement with the next higher number of corners which is also often used in applications is a square. Similar to a triangle a square is also expected to have commensurate occupation numbers which in this case are the square numbers. 92 6. Skyrmion Dynamics in Confinements 1 2 3 4 5 6 7 8 9 10 Figure 6.6.: Histograms of 1 to 10 skyrmions in square confinement. Skyrmions are displayed as white circles with a radius equal to the skyrmion-skyrmion interaction width 1. This figure is taken and rearranged from our publication [110]. For this case no experiments but only simulations have been performed aiming at expanding the results for commensurability to another geometry. The simulated square has a side length of 21.21 nm leading to an area similar to the area in the triangle. Figures 6.6 and 6.7 show that a commensurate state is indeed found for four and nine skyrmions. This commensurate state is a square lattice with the same symmetry axes as the confinement which leads to a significantly lower MSD. As already seen for the triangular confinement, non-commensurate numbers of skyrmions often have similar favored positions as the commensurate state with the next-higher number of skyrmions. However, because there are less skyrmions than lattice positions available this leads to skyrmions moving between lattice sites and therefore a higher MSD compared to a fully filled lattice. The structures are also generally somewhat less stable than for a triangle and allow for collective rotation with a lower energy barrier as the slight increase in the long-time MSD for four skyrmions indicates. The method described here could also be applied to various other geometries and would also show commensurability effects for different polygons as a hexagon or for lower-symmetry confinements as a rectangle with different sides or a star geometry. However, because the physical principle behind this kind of commensurability effects is the same for all geometries and experiments and computing applications mainly use the three geometries mentioned above [99, 104], analysis of further geometries is omitted. 93 6. Skyrmion Dynamics in Confinements Figure 6.7.: Short-time MSD for different numbers of skyrmions in a square con- finement. The MSDs are averaged over all skyrmions. ΔC is given in simulation units. This figure is taken and slightly adapted from our publication [110]. 6.2. Confined Dynamics with Applied Currents 6.2.1. Single Skyrmion Expanding on the results for skyrmion ordering in confined geometries, one can also analyze skyrmion diffusion in such geometries with an applied current in some direction. A very simple case is using a triangle with three input voltages placed at each of the triangle’s corners. While in general all three inputs can have a continuous spectrum of applied voltages, one can further simplify the state space to levels + (=̂ true) and 0 (=̂ false) similar to the two possible values of a bit in a computer. Before doing experiments with applied currents, computer simulations are performed to gain insight in possible states that can be reached this way. In this system, simulations are performed using the repulsive exponential function determined by fitting the iterative Boltzmann inversion result given in Figure 3.7. Thereby, simulations and experiments are performed on the same length scale, which allows for direct matching of the results. In this case, a triangle of side length 40 µm is used. Beyond the interactions included in previous simulations, the simulations are extended by using a current-profile for the given material properties and system geometry generated by Maarten A. Brems using COMSOL [191]. Here, the current is transmitted into the sample via gold contacts on rectangles which 94 6. Skyrmion Dynamics in Confinements Figure 6.8.: Plot of the current density distribution for a triangle with a positive voltage applied to the bottom right corner and ground applied to the other two corners. The left and center plot show the directional currents in x- and y-direction while the right plot shows the absolute current. Current densities are given in units of A/m2. Current density distributions were calculated with COMSOL [191] by Maarten A. Brems. Figure 6.9.: Histograms of one skyrmion simulated in a triangular confinement of side length 40 µm without applied voltage (left), with plus voltage applied to the bottom right corner (center) and with plus voltage applied to both bottom corners (right). Skyrmions are displayed as white circles with a radius equal to the skyrmion-skyrmion interaction width 1. expand from the corners of the triangle. The applied currents induce spin-orbit torques which move the skyrmions. Due to the strength of spin-orbit torques being proportional to the applied current, one can directly input the current profiles with a heuristic proportionality factor to determine skyrmion motion in the sample. An example current profile is given in Figure 6.8 for applying a positive voltage at the bottom right position and ground at top and left. The strongest current is present at the corners of the triangle where small notches separate the triangle from the rectangle at which the current is applied. This can potentially lead to skyrmion annihilation at these position when the current through the sample is too large. For smaller currents, however, this is not problematic because the edge repulsion is strong enough to prevent skyrmions from moving into the corner. Besides these notches, the current is flowing into the bottom right corner and is highest close to this corner and lowest on the opposing triangle edge. The voltage distribution described above is one of the three different current distributions possible for the given system when not considering rotated configu- 95 6. Skyrmion Dynamics in Confinements Figure 6.10.: MOKE pictures of preferred skyrmion positions for five different combinations of +, - and 0 inputs at the three corners of a 36 µm triangle. The combinations are shown in the pictures and below in the brackets. Possible positions of four magnetic tunnel junctions which would be capable to read out skyrmion positions for different input combinations are displayed in (e). These positions can be used for reliable readout of Boolean logic operations. The maximum current density in these pictures is around 108 Am−2 which is around twice the required current for visible skyrmion motion but below the current leading to skyrmion annihilation in the corners. This figure is taken from our publication [99]. rations separately. The other two relevant distributions are three ground-inputs or equivalently three plus-inputs in which case no current flows at all and two plus-inputs where the current flows symmetrically into the two corners as displayed in Appendix A.13. Using these current profile, one can create the histograms of skyrmion occurrence in Figure 6.9. These histograms show that with only one positive voltage, the skyrmion is pushed into a corner. With two corners set to a positive voltage, the skyrmion is pushed onto the triangle edge between the two corners with positive voltage. The range of possible positions on this side of the triangle also includes both corners, thereby making distinction from the one-positive case difficult when only observing individual frames. This result indicates that experiments performed with this system probably require some time averaging to be able to reliably distinguish e.g. the ++0 state from the 0+0 and +00 states in which a skyrmion is located in one corner. Transferring this system to experiments, my collaborators recorded not only videos of the cases described above but also different structures by allowing for an ad- ditional negative input which has the same potential difference with respect to the ground but the opposite sign. As displayed in Figure 6.10, skyrmions without applied voltage move to the triangle center while when applying voltages, skyrmions tend to move to the corner with positive voltage and away from the negative voltage. Building on these investigations on skyrmion motion with applied currents, a simple reservoir computing device capable of reliably calculating Boolean logic operations has been constructed [99]. It uses plus and ground voltages as Boolean input signals and calculates results out of them by weighting the readout of skyrmion occurrences 96 6. Skyrmion Dynamics in Confinements at specific positions defined in Figure 6.10e. This computer is able to calculate all Boolean two-input calculations as well as several three-input calculations reliably. Further development of this system also aims to perform more difficult tasks such as gesture recognition. 6.2.2. Usage of Multiple Skyrmions The number of available states in a skyrmion system consisting of only one skyrmion is usually quite small due to the lack of skyrmion-skyrmion interaction. When increasing the number of skyrmions in a system, this interaction allows for more distinguishable states. However, when using a larger number of skyrmions, one needs to consider commensurability effects, which may lead to a fixed lattice state and therefore a reduced number of available states as discussed in Section 6.1.2. Examples of different states with an applied current and varying skyrmion numbers are given in Appendix A.14. This also includes simulations incorporating wires at each corner of the triangle in Appendix A.15 showing that in this more realistic setup, wall repulsion still suffices to keep skyrmions inside the triangle. As modern computing concepts such as reservoir computing usually require a large number of states to accomplish complicated tasks [91, 99], a more elaborate study of skyrmion states with applied current was performed for the non-commensurate number of four skyrmions in a triangular confinement. To obtain a clear direction of skyrmion movement, this study is focused on the 0+- voltage combination in particular, in which skyrmions are asymmetrically pushed into one corner. The resulting current density distribution calculated for this state is shown in Appendix A.16. To quantify the states in this system, the simulation temperature :B) is varied to obtain fixed lattice states at a low temperature as well as more possible states which are reached at higher temperatures. These structures as well as a low-temperature structure without applied current are shown in Figure 6.11. At low temperature, the four skyrmions have fixed positions which are stabilized by high energy barriers. The bottom left skyrmion (blue) has some freedom of movement while all other positions are relatively fixed. When increasing the temperature, some collective rotation of the three skyrmions outside the corner becomes possible while the skyrmion in the bottom right corner still more or less keeps its position. The histogram shows that while this switching allows for a certain degree of rotation, the most common skyrmion positions are still very similar to the lower temperature case before. At even higher temperature the separation of the three and one skyrmions also ends and all skyrmions are able to switch their position in the given time. This increase in energy also opens up more shapes of the three skyrmions which no longer align in the same structure in all frames but can perform some degree of collective rotation. While the temperature dependence of the results is clearly visible, due to the pinning-dominated diffusion of skyrmions as well as the temperature dependence 97 6. Skyrmion Dynamics in Confinements Figure 6.11.: Top: Positions of the skyrmion centers in a simulation over two million simulation frames recorded every ΔC = 0.1 (a) at a low temperature of :B) = 0.1 without applied current, (b) at the same temperature with applied current, (c) at a higher temperature of :B) = 0.25 with the same current, (d) at a temperature of :B) = 1.0 with the same current. The black dots indicate cluster centers of a k-means algorithm used on this data which is explained in our publication [111]. Bottom: The same data plotted using the same style as in Figure 6.4(b). The upper part of this figure is taken from our publication [111]. of the spin structure, this temperature variation in simulations cannot be directly mapped to real temperatures in the experiment. Instead different temperatures describe systems with more or less thermal diffusion which is a tunable parameter in experimental realizations. Expanding on this qualitative analysis described above, one can also use a variety of techniques for cluster recognition such as k-means or Generalized Perron Cluster Cluster Analysis (GPCCA) to get a more quantitative insight in the prediction power of such a skyrmion reservoir. While these methods are not described in this thesis, one can find them in the related paper [111]. The GPCCA method suggests a separation of the given data into two different clusters as shown in Appendix A.17 where one is an arrow-like state pointing into the right corner and the other one is a rhombus like shape. 98 7. Conclusion and Outlook 7. Conclusion and Outlook 7.1. Conclusion In this thesis I have presented results obtained from computer simulations of mag- netic skyrmions as well as from experimental analysis of skyrmions. The simulations presented here are based on the Thiele model which describes skyrmions as rigid quasi-particles diffusing according to Brownian equations of motion. This descrip- tion of skyrmions allows for very efficient simulations compared to conventional atomistic or micromagnetic models and thereby enables simulations of systems containing large numbers of skyrmions. A substantial constituent of this description are interaction potentials which describe the interactions of skyrmions with each other and their interactions with the mag- netic material including skyrmion-boundary interactions and pinning landscapes. Skyrmion-skyrmion and skyrmion-boundary interactions have been determined using iterative Boltzmann inversion, a method originally known from soft matter physics. This way, the potentials are directly ascertained based on experimental data without making any assumptions about their functional shape. The deter- mined potentials were found to be purely repulsive and have a shape similar to a falling exponential which is in agreement with theoretical predictions [21]. This is an excellent example of how experimental data can be employed to improve simulation models such as the Thiele model used here. The model reproduces and predicts skyrmion behavior reliably and can match the properties of a skyrmion lattice recorded in an experiment in simulations. With this, the model greatly improves the predictive power and accuracy of computer simulations of large skyrmion lat- tices. The empirical potentials are first step towards comprehensive and predictive quantitative modeling of the dynamics of micrometer sized skyrmions in computer simulations. Such modeling is particularly valuable as state of the art atomistic and micromagnetic models are incapable to accurately describe such large systems due to the exceeding computational cost. Furthermore, the effect of skyrmion chirality is investigated, which introduces an altered diffusion mechanism. In the case of strong gyroscopic force and therefore a large skyrmion Hall angle, this results in increased diffusion at higher skyrmion densities. Given this, skyrmion systems are an excellent experimentally accessible example for a novel diffusion mechanism called odd diffusion [84]. Furthermore, this thesis presents two different methods to increase skyrmion diffusion by miti- gating pinning effects present in the sample. This increase in diffusion provides a large speedup for several skyrmion-based devices and is an important step to make such devices feasible. Another ingredient in skyrmion-based computing is 99 7. Conclusion and Outlook the understanding of skyrmion ordering and diffusion in confinements which is investigated in this thesis as well. The analysis provides insight into promising choices of confinements and skyrmion numbers for e.g. reservoir computing. In summary, all of these different applications described in this thesis show how Thiele model simulations of skyrmions and experiments can be used together to thoroughly predict statics of skyrmion lattices on experimental length scales. Combined with better insight into matching of experimental and simulated time scales, this experimentally parametrized Thiele model could also predict skyrmion dynamics at the same time and length scales as in experimental systems. In the next section, I discuss the limitations of this model and after that I present further points which could be analyzed in the context of the model and the results presented here. 7.2. Limitations of the Particle-Based Model for Skyrmions All simulations performed in the context of this thesis approximate skyrmions as circular isotropic, non-deformable quasi-particles and neglect processes of skyrmion creation and annihilation. These simplifications lead to a robust, low-complexity particle-based model accessible by standard methods from soft matter physics. Such particle-based models of skyrmions require significantly lower computation times compared to more detailed models like micromagnetic or atomistic simula- tions. This allows for simulations of lattices and dynamics of micrometer-scale skyrmions as demonstrated here. Micromagnetic simulations would only be able to capture static effects at this scale and atomistic simulations are not feasible for systems of hundreds of micrometers. The limitation of micromagnetic simulations is particularly relevant here due to the requirement of a certain number of cells to correctly describe the domain wall of a skyrmion. In the systems described in this thesis the domain wall only comprises a very small part of the skyrmion spanning some nanometers for a micrometer-sized skyrmion leading to a very large number of cells required for each skyrmion. Besides the computational reasons for using a particle-model, the model is limited to parameters which can be measured with high reliability in experiments. Since the MOKE setup only gives access to the magnetization with a spatial resolution of slightly less than 1 µm and a time resolution of up to tens of milliseconds, skyrmions can only be recorded as homogeneous disks without resolution of their domain walls or internal excitations. Theoretical models and calculations for the skyrmions used in the experiments also only grant limited insight. For instance, they usually do not include all relevant effects in the samples such as defects or temperature. For 100 7. Conclusion and Outlook these reasons, the model used here, while having some limitations, is well capable of predicting static and dynamic properties of skyrmions seen in experiments. This includes diffusion, lattice ordering and motion under applied current. Effects such as deformations of skyrmions due to close packing or external cur- rents [189] and internal excitations like the breathing mode [192–194] are not described by the model. However, such excitation modes are typically most relevant in systems with large applied currents [72, 189] or fluctuating in-plane magnetic fields [194]. In this thesis, all applied magnetic fields are directed out-of-plane and therefore impact the whole skyrmion isotropically and do not lead to deformation. Also, skyrmions are assumed to be mono-disperse which is in good agreement with our experiments where the standard deviation of sizes is below 20% of the average skyrmion size in all observed cases. Skyrmion creation and annihilation processes are negligible in skyrmion lattices and for the diffusion of individual skyrmions due to the generally high stability of magnetic skyrmions. These processes would become relevant for large applied external forces [99] or for confinements which are smaller than one skyrmion diameter. Such systems are therefore not well-described by this model. Also, transitions e.g. between skyrmions and worm domains are not modeled, which limits the applicability of this model to systems which host only skyrmions. Another limitation of the particle-based model are the types of interactions used to describe skyrmion dynamics. The interactions used here are two-body interac- tions between skyrmions and interactions between individual skyrmions and their surroundings such as material boundaries, non-flat energy landscapes or external currents. This neglects possible multi-particle interactions and by applying a potential cutoff long-range two-body interactions are neglected as well. This is justified by how the interaction potentials calculated for a skyrmion lattice are capable to reproduce local skyrmion ordering without higher-order or long-range interactions. Similar limitations also apply to pinning effects which are simulated by pinning of the skyrmion center while in reality the domain wall at the skyrmion boundary is pinned [67]. While this different mechanism does not have a large effect on diffusion and lattice formation which can be modeled well by pinning of point-like particles, effects like the quantitative size dependence of pinning [41, 67] cannot be reliably described by this quasi-particle model. While some of the limitations discussed here are intrinsic to the particle model, others could generally be overcome by conducting further research. For example, one could use the distribution of skyrmion sizes and assume some size-dependence of the skyrmion interaction potential to describe lattice formation of a polydisperse skyrmion system. One could also add a probabilistic annihilation of skyrmions that move too close to a sample edge or a probabilistic creation at specific locations. 101 7. Conclusion and Outlook 7.3. Outlook While computer simulations of skyrmions in the Thiele model are already providing insights into many phenomena involving skyrmions, there is still a lot more to be explored. At the fundamental level, the phases and phase transitions of skyrmion systems are still not fully understood. This is, on the one hand, a consequence of pinning which promotes formation of grain boundaries and therefore often pre- vents reaching a solid phase. On the other hand, experimental probing of phase transitions is difficult as changes in temperature or skyrmion density impact the skyrmion size and thus the skyrmion-skyrmion interactions. Alternatively, one could induce phase transitions by applying field variations which induce a change in skyrmion size and thus the effective density. Further analysis of this topic would include improving the understanding of pinning and of skyrmion deformations and size variations. In addition, one could analyze topological defects in skyrmion lattices to understand how defects are related to phase transitions and how pinning and domain boundaries affect topological defects. This study of topological defects in a lattice of topological objects would introduce a second level of topology to skyrmion research. While so far the analysis of phases has been performed in bulk systems, another promising path would be experiments and simulations in large confinements containing hundreds to thousands of skyrmions. As - due to the sixfold symmetry of a skyrmion lattice - a hexagon is assumed to improve ordering, other shapes such as a square are incommensurable with a skyrmion lattice. This incommensurability is expected to lead to grain boundaries in the skyrmion lattice and thus to worse overall ordering. A hexagon, on the other hand, would induce the correct ordering at the edges and thereby simplify reaching a state of solid ordering compared to a bulk system. While this state is not the same as a solid phase in the continuum, it can still provide insight into ordering processes and possibly also phase transitions. Besides the analysis of high-density skyrmion structures and phase transitions, the analysis of the diffusion of individual skyrmions also promises further insights. One point is to match time scales between experiments and simulations more reliably. Knowledge of the time scales of skyrmion motion in the Thiele model would also open up the possibility of determining the scale of spin-orbit torque-induced forces in the system. Knowing these two scales as well as interaction potentials deter- mined with the method described in this thesis would allow for the creation of a comprehensive quantitative computational model of skyrmions. This model would be able to predict experimental behavior realistically on experimental time and length scales. Another very important part of understanding skyrmion diffusion is understand- ing the mechanisms leading to pinning. While it has already been found that micrometer-sized skyrmions are pinned at their domain walls [67] and methods to 102 7. Conclusion and Outlook overcome pinning have been described in this thesis, the exact physical mechanism behind pinning is still unknown. Learning about this would improve the model- ing of pinning sites beyond the simple models that are currently currently used. Additionally, it would likely make it easier to tune pinning effects in a material and thereby to increase or reduce skyrmion diffusion. Further analysis of pinning sites could also include an interaction between the gyroscopic force and pinning which has already been explored in theory [57, 86, 195, 196] but not explicitly in experiments. Beyond these directions for fundamental research on skyrmions, improved and novel skyrmion-based computing devices are also of great interest. Current approaches of reservoir computing [99] are limited by the time resolution of the current ex- perimental setup as well as pinning effects. By solving these hindrances, reservoir computing schemes could be extended to allow for experimental realization of more complicated tasks such as audio recognition [104] or handwritten digit reading [102]. Token-based Brownian computing, on the other hand, is mostly limited by pinning effects as these prevent fast and reliable skyrmion motion [43]. As mechanisms of pinning as well as depinning via external excitations are researched, this improved understanding can also be used to speed up token-based computing. Beyond these already known schemes, skyrmions are also a promising candidate for the realization of other unconventional computing schemes due to the variety of ways to manipulate them as well as their high stability. In conclusion, magnetic skyrmions are a very interesting physical system relevant for fundamental science as well as promising for various applications. Further work will hopefully be able to improve the understanding of this interesting 2D system and enable the realization of competitive skyrmion-based computing devices. 103 A. Appendix A. Appendix A.1. Additional Functions A.1.1. Running Average A running average 〈G〉# over a specified numb∑er of neighbors # on each side iscalculated as :+# 〈G〉# (:) 1 = 2 + 1 G8 . (A.1)# 8=:−# A.1.2. Standard Error of the Mean The standard error of the mean describes the statistical error of an estimated mean value calculated from normally distributed variables. It is given as √ Δ〈G〉 = f/ #, (A.2) where f is the sample standard deviation and # is the number of sampled values. A.2. Technical Explanations A.2.1. Determination of the Material Boundary Position To determine the position of the material boundary, a background image shown in Figure A.1 is used. This image is taken while the sample is in a homogeneous magnetic state without skyrmions. Due to the homogeneous magnetization during this measurement, different grayscale values indicate the magnetic material, the non-magnetic part of the sample and defects of either the material or the lens. Using the grayscale difference between the magnetized and non-magnetized part, one can fit the intensity distribution in y-direction to a linear function determining the point where the intensity is exactly the average between the two parts. This position can then be used as the material boundary position. It is, however, still probably not the exact position of the boundary and therefore subject to further corrections. The linear fit is shown in Figure A.2. The procedure described here was done by Yuqing Ge [108] and is part of this thesis only for completeness. 104 A. Appendix Figure A.1.: Background image of the magnetic material used for determining skyrmion interaction potentials. The material boundaries are clearly visible as well as some material defects and camera lens or CCD sensor defects. 0.8 averaged intensity 0.7 0.6 10 m 0.5 0.4 0.3 10.0 12.5 15.0 17.5 20.0 22.5 25.0 27.5 y [ m] Figure A.2.: The boundary position is determined by using the average intensity at each y-position. The red and blue lines show the intensity at two different x-positions indicating the uneven illumination of the sample. The crosses show the average intensity for both curves and therefore the estimated boundary position for these x-positions. The average over all x-pixels is the boundary position used here. The inset is a Kerr image showing the x-positions of the red and blue curve as well as the determined boundary position as the cyan line. This figure is taken from our publication [108]. 105 Normalized intensity A. Appendix A.2.2. Derivation of the Size Dependence of Skyrmion Diffusion Deriving the diagonal components of the skyrmion diffusion tensor requires some calculations, which are presented in this part of the appendix. Due to the rotational isotropy of the skyrmion, the result of this calculation also applies to HH. As stated in the main part of this thesis in Section 5.1.2, a Néel skyrmion with a tanh-shaped domain wall is assumed. Starting with the standard definition of the diffusion tensor and inputting the spin structure, one obtains the following calculat∫ion: ∫ 2c ∫ ∞ [ ] GG ∼∫ (mG(®) 2 2 ∫3 A = ( 3i √3AA (m 2 2 2 G(G)( + (mG)() H) + (mG(I)0 02c ∞  2A − 0 = 3i 3AA  m cos i 1 − tanh2G0( 0√ (A − 0 )) ( 1 2 (  + sin 1 − tanh2 + tanh A − 0 ))2 ∫ mG i m  G 1 1  2c ∫ ∞ [ ( ) ( (A − 0 )) = ( √3i 3AA ( (m 2 2 G co)s)i) (+ (mG sin(i) 1))−tanh 2 0 0  12 2 + mG 1 − tanh2 A − 0 + mG tanh A − 0  ∫ 2c ∫ ­©( 1 ) ( √ )1 ∞ 2 G 2« + 1 − G 2 ®ª ( (¬ 1 − tanh2 A − 0 )) = 3i 3AA mG mG 0 A A2 1( ) 0©©  ( ) ª2 ( )2ª2 tanh A−0 + mA ­­ √ 1 mG ­««­ 2 ( − ) ( ) ®® + 1 2 ( A− )0 ®® ∫ ∫ 1 cosh A 0 ( 1)− tanh 2 A−0 1 cosh  1  2c ∞  1 1 ¬ ¬ =  ( ) 3i 3AA (  © «­­ 1 G2 2 © GA− 3 + ­) «­ √ 2 1 − G2 A A3 ®¬® ª ®ª ( ( ))¬® 1 − tanh 2 A − 0 0( ) 0 ( A) A 1 − G2 ] 1A22 tanh2 A−0 + G (1 ) + 1 · 1 ( ) ∫ A ∫1 − tanh 2 A−0 12 cosh4 A−0 2 ( 1 ) 1c ∞ ©­­ 1 − cos2 2i= 3i 3AA « + ­ © ( )­co« √ s 1−cos 2 2i i A ª® ª® ( (1 − tanh2 A − 0 )) 0 0 A 1 − cos2 ® ®i ¬ ¬ 1 106 A. Appendix ∫+ cos 2 1i 2[( ] ((− ) · 1 1 − tanh A 0 12 ( ) cosh4 A−0 1 2 ) ( 1 ( ) ))c ∫ ∞ 1 − cos2 2i 1 + cos 2i 1= 3i 3AA 0 0 ( ) A ] 1 − cos2 ( )i cosh2 A−01 + cos2 cosh2 A − 0 1∫ i ∫ [ · ( )1 1(2 cosh4 A−01 ) ]2c ∞ ( )sin4 i 1 + cos 2i cos2 i 1=∫ 3i ∫ 3AA0 0 [ A2 sin] + ( )i 12 cosh2 A−012c ∞ sin2 i cos2 i 1 =∫ 3i[ 3AA0 0∞ c cA ] 2 + A 12 2 ( ) cosh A−0 1 =∫ 3A + 1( ) 0 A 12 cosh2 A ∞ c c [−01 (A − 0 ) (A − 0 )]∞ =∫ 3A 2 ( − ) + A tanh − 1 log cosh0 A cosh A 0 1 1 1 01∞ [ ( A−0 )c ( ) c 4 (−0 )]1=∫ 3A 2 − + lim A − 1 log + 1 log cosh0 A cosh A 0 A→∞ 1 2 11∞ c =∫ 3A 2 ( − ) + lim c [ (A − 0 ) ( )] A − 1 − log 2 + 1 log cosh 0 0 A cosh A 0 A→[∞ 1 1 ( 1( 1∞ c c 0 )]= 3A 0 A cosh2 A− ) + 0 + 1 log 2 + 1 log cosh 0 1 1 1 ∫Using the assumption that the skyrmion size 0 is much larger than the domainwall width 1( , o)ne can further simplify this calculation using cosh( 0 ) ≈ 1 0/1∞ 1 24 and3A c0 cosh2 A−0 ≈ 0 for 0 [ 1:A 1 c (0 )] 2c0 GG ∼ 0 + 1 log 2 + 1 − log(2) = (A.3) 1 1 1 107 A. Appendix A.3. Additional Figures Figure A.3.: Two-dimensional radial distribution function 6(G, H) for (a) experiment and (b) simulation. The functions are normalized as defined in Equation 2.97. Figure A.4.: Experimental occurrence landscape of skyrmions at a high skyrmion density. 108 A. Appendix Figure A.5.: Logarithmically plotted histogram of skyrmion occurrence in a com- puter simulation of single skyrmions using the experimental pinning profile in Figure 3.10. 100 = 0.36 10 1 10 2 5 7 10 20 30 50 r [ m] Figure A.6.: Translational correlation function @ (A) averaged over the full mea- surement of 5 minutes for the lattice shown in Figure 4.7. The dashed line is a fit of @ (A) = 0A−[ against the upper envelope of this function. 109 Gq(r) A. Appendix 0.015 Rayleigh 0.010 Data Gauss 0.010 Data 0.005 0.005 0.000 100 0 100 0.0000 50 100 150 x-Displacement [ m] Absolute Displacement [ m] Figure A.7.: Histograms of skyrmion displacement after ΔC = 20 s determined as a running average over longer experiments for the excited case with 5 = 2 kHz and  = 2100 `T. (a) shows the displacement in x-direction, (b) shows the absolute displacement. Fits are performed as in Figure 5.7. Gauss Rayleigh 0.4 Data Data0.4 0.2 0.2 0.0 10 0 10 0.00 5 10 x-Displacement [ m] Absolute Displacement [ m] Figure A.8.: Histograms of skyrmion displacement after ΔC = 62.5 ms determined as a running average over longer experiments for the excited case with 5 = 2 kHz and  = 2100 `T. (a) shows the displacement in x-direction, (b) shows the absolute displacement. Fits are performed as in Figure 5.7. 110 Probability Density Probability Density Probability Density Probability Density A. Appendix Figure A.9.: Diffusion coefficient for skyrmions excited by a sinus-shaped alternating current of different frequency and constant amplitude  (C) = 0 sin(2c 5 C) plotted on a logarithmic scale. The red line indicates the diffusion without an applied current. The falling slope of the data is fitted against  (C) ∼ 5 −U with U = 0.51± 0.02. This plot is taken with permission from the bachelor thesis of Tobias Sparmann. 111 A. Appendix 0.96 No Excitation (Theory) Excited 0.95 0.94 0.93 100 101 102 Frequency [1/ t] Figure A.10.: Diffusion constant of skyrmions with an applied sinus force with fixed amplitude and varying frequency  (C) = 0 sin(2c 5 C) without pinning or skyrmion- skyrmion interactions. The red line shows the diffusion without an applied force. Error bars are the standard error of the mean. 0.6 No Excitation 0.4 Excited 0.3 0.2 0.1 100 101 102 103 Frequency [1/ t] Figure A.11.: Diffusion constant of skyrmions with an applied sinus force with fixed amplitude and varying frequency  (C) = 0 sin(lC) at density 0.00833 µm−2. The red line shows the diffusion without an applied force. Errors are smaller than the size of the points. 112 Diffusion Constant [ m2/ t] Diffusion Constant [ m2/ t] A. Appendix 2 7 3 8 1 4 910 5 10 6 100 10 1 100 101 102 t Figure A.12.: Short-time MSD for different numbers of skyrmions in a circle in a computer simulation calculated using the windows-sliding method. All parts of trajectories closer than 3 nm to the center of the circle are omitted for the calculation. The MSDs are averaged over all skyrmions. The one-skyrmion case is omitted as the only skyrmion in this case is located in the center. ΔC is given in simulation units. Figure A.13.: Plot of the current density distribution for a triangle with positive voltage applied to the two bottom corners and ground applied to the top corner. The left and center plot show the directional currents in x- and y-direction while the right plot shows the absolute current. Current densities are given in units of A/m2. Current density distributions were calculated with COMSOL [191] by Maarten A. Brems. 113 MSD [nm2] A. Appendix Figure A.14.: Histograms of 1 to 5 skyrmions (from top to bottom) placed in a triangle drawn in the same style as in Figure 6.2. In the left column a positive voltage is applied to the two lower corners and ground to top. In the center column a positive voltage is applied only at the bottom right corner with zero voltage on the other two corners. In the right column a positive voltage is applied to the bottom right corner and a negative voltage of the same absolute value to the bottom left corner with zero voltage at top. Current density distributions were calculated with COMSOL [191] by Maarten A. Brems. 114 A. Appendix Figure A.15.: Histograms of 1 to 5 skyrmions placed in a triangle with wires on each corner drawn in the same style as in Figure 6.2. A positive voltage is applied to the bottom right corner while zero voltage is applied to the two remaining corners. For all numbers of skyrmions, the probability that a skyrmion moves into a wire is zero. The current density distribution was calculated with COMSOL [191] by Maarten A. Brems. Figure A.16.: Plot of the current density distribution for a triangle with positive voltage applied to the bottom right, a negative voltage of the same absolute value applied to the bottom left and ground applied to the top corner. The left and center plot show the directional currents in x- and y-direction while the right plot shows the absolute current. Current densities are given in units of A/m2. Current density distributions were calculated with COMSOL [191] by Maarten A. Brems. 115 A. Appendix Figure A.17.: Occurrence distribution of skyrmions in the two different states deter- mined with the GPCCA method for different temperatures. The blue probability distribution is the rhombus-like state, the red one is the arrow-like state. The blue ring indicates the rotation center of the three central skyrmions. This figure is taken from our publication [111]. 116 A. Appendix A.4. Polymer Physics This section contains works on polymer physics in which I participated during my PhD studies. As these are side projects of my studies, I do not include a throughout explanation of these works and the underlying theory but instead just attach the papers as a part of this thesis. The first work in this chapter is J. Rothörl, S. Wettermann, P. Virnau, and A. Bhattacharya, “Knot formation of dsDNA pushed inside a nanochannel,” Sci. Rep. 12, 5342 (2022). Here, Langevin dynamics as well as Monte Carlo simulations of a polymer chain in a confined nanochannel were performed to analyze how the channel as well as a force applied on one side affect the structure and knotting behavior of the polymer chain. The main result of this work is that knotting of polymers is increased when pushing the polymer through a channel. Simulations of this system were performed by Aniket Bhattacharya and Sarah Wettermann without my participation. I performed most analysis of the simulation results. Additionally, I participated in writing the paper. The second work presented in this chapter is Y. Zhao, J. Rothörl, P. Besenius, P. Virnau, and K. C. Daoulas, “Can polymer helicity affect topological chirality of polymer knots?” ACS Macro Lett. 12, 234–240 (2023). In this work, knotting behavior of chiral polymers was investigated by Monte Carlo simulations. The main result of this paper is that the chirality of knots found in a polymer chain depends on the local chirality of the polymer introduced by a local chiral inter- action. Additionally, a general reduction of knotting probability as well as an increased abundance of specific knots in a system with chiral interactions was found. Simulations in this work were performed by Yani Zhao and Kostas Ch. Daoulas without my contribution. I performed the knot analysis and participated in writing the paper. The third polymer physics paper I worked on is J. Rothörl, M. A. Brems, T. J. Stevens, and P. Virnau, “Reconstructing diploid 3D chromatin structures from single cell Hi-C data with a polymer-based approach,” in press, 2023. It describes a simulation protocol for creating 3D chromatin structures for diploid cells. It includes a method to assign ambiguous contacts between chromosome positions based on 3D structures and allows for efficient simulations at resolutions up to 5,000 base pairs per bead. The simulations presented in this work were performed by myself mostly during the time of my master thesis. Analysis was mostly performed by myself with help of Maarten A. Brems. As this paper is not yet published at the date of submission of this thesis, it is not included here. 117 www.nature.com/scientificreports OPEN Knot formation of dsDNA pushed inside a nanochannel Jan Rothörl1, Sarah Wettermann1, Peter Virnau1* & Aniket Bhattacharya2* Recent experiments demonstrated that knots in single molecule dsDNA can be formed by compression in a nanochannel. In this manuscript, we further elucidate the underlying molecular mechanisms by carrying out a compression experiment in silico, where an equilibrated coarse-grained double-stranded DNA confined in a square channel is pushed by a piston. The probability of forming knots is a non-monotonic function of the persistence length and can be enhanced significantly by increasing the piston speed. Under compression knots are abundant and delocalized due to a backfolding mechanism from which chain-spanning loops emerge, while knots are less frequent and only weakly localized in equilibrium. Our in silico study thus provides insights into the formation, origin and control of DNA knots in nanopores. In many biological processes a double-stranded DNA (dsDNA) is confined in a geometry much shorter than its contour length in a highly organized and compact state and often under high p ressure1,2. A classic example is an organized state of a dsDNA strand in a viral c apsid3–8. The viral DNA uses the stored elastic energy for its inva- sion process. Intriguingly, this DNA was found to be highly knotted particularly in a mutant variant for which both sticky ends are allowed to reside within the capsid3,4. It is in general difficult to develop an experimental protocol to study an actual system in vitro, although there have been studies to measure the force and the organ- ized topology of the dsDNA inside a c apsid9,10. During the last decade advancements in nanotechnology have enabled us to prepare nanochannels of sub-persistence length dimensions11. DNA pushed inside nanofluidic devices12–16 is now used for mapping genomes, sequence motifs, structural v ariations17,18. Recently, nanochannels were even used for the detection of knots in DNA as demonstrated in experiment19 and simulation20. Due to a controllable and simpler geometry, nanochannels offer immense promise to under- stand universal aspects of biological phenomena using well established concepts from polymer p hysics21,22. Besides problems of biological significance and of human health, nanochannel based experiments claimed the occurrence of j amming9,10- which indicates that confined bio-polymers offer yet another platform to study slow relaxation and glassy dynamics. Thus studies of chain compression in nanochannels appeal to broad areas of science. Many numerical studies of knots have established numerous results on generic23–38 and biopolymers1,6,8,39–45. While knots are, e.g., abundant in single ideal c hains24,25, the addition of excluded volume typically reduces the knotting probability significantly25,29, while spherical confinement or globular states enhance k notting26,29. For the latter, knots tend to be delocalized, i.e. average sizes scale linearly with the chain length29, while knots are weakly localized and scale sub-linearly for under ideal or good solvent c onditions28,29. Previously, we have studied compression of semi-flexible polymers in nanochannels using a Langevin dynam- ics (LD) scheme46,47. These LD simulation studies and another recent s tudy48 have provided substantial insights about many details at smaller length scales unattainable experimentally, but are essential for microscopic under- standing and interpretation of nanochannel experiments using fundamental laws of physics. One of the advan- tages of these simulation studies is that, unlike an actual experiment, one can vary confining dimensions and chain stiffness easily and is thus capable of extending the simulation studies for a broader parameter space which is often quite expensive to design experimentally. In our recent LD simulation study47 we mimicked a recent experiment in silico where dsDNA—modeled as a semi-flexible polymer—was pushed inside a rectangular open-ended nanochannel much longer than required to attain a steady state. By varying the bending rigidity of the chain we showed how the structure evolves from a disordered state to a highly organized spooled state. Furthermore, the LD simulation revealed a detailed picture of how the fold nucleation originates at the piston end and expands during the compression process47. An important and relevant question in these compression studies in the biological context is to study how the formation of knots are initiated and once formed how they spatially evolve under confinement. Theoretical and 1Institut für Physik, Johannes Gutenberg-Universität, Staudinger Weg 9, 55099 Mainz, Germany. 2Department of Physics, University of Central Florida, Orlando, FL 32816-2385, USA. *email: virnau@uni-mainz.de; Aniket.Bhattacharya@ucf.edu Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 1 Vol.:(0123456789) www.nature.com/scientificreports/ Figure 1. Knots and radii of gyration of steady-state configurations for different input parameters. (a) Root- √ mean-squared radius of gyration 〈R2g (t)〉 for chains of different stiffness. Error bars are omitted as their size is smaller than the size of the points. (b) Knotting probability and occurrence probability of trefoil knots for different values of the bond stiffness κ . Error bars are determined by taking the standard deviation over the square-root of the number of runs. Lines are added for readability. Knotting probability for v = 0 is strictly zero √ for κ ≥ 50 . (c) Trefoil knot lengths found for different κ with a moving piston and in equilibrium. (d) 〈R2g (t)〉 for different piston velocities at κ = 4 . Again the size of error bars would be smaller than the size of the points. The red dashes in (d–f) indicate the average result in MC simulations without a moving piston for κ = 4 which is equivalent to the leftmost equilibrium values in plots (a–c). (e) Knotting probabilities for different piston velocities at κ = 4 . (f) Trefoil knot lengths for different piston velocities at κ = 4. simulation studies have been further fueled by a recent experiment that demonstrates that knots indeed occur in compression experiments16. In computer simulations studies Orlandini and Micheletti have already investigated equilibrium knot formation of coarse-grained DNA models in nanonchannels31,32. Of particular relevance to our work is a recent study38 in which the non-equilibrium formation of knots and so-called geometrical entangle- ments as measured by counting crossings under projections are investigated in closed nanochannels exposed to compression and relaxation cycles. It was found that the two types of entanglements evolve with different dynamics and are for the most part uncoupled. Here, we investigate knot formation when a confined dsDNA is being pushed by a nano-dozer in a nanochan- nel whose width is much smaller than the contour length of the dsDNA. An important difference from previous studies is that our system is open ended in one direction and that we study the evolution of knots in a constantly moving steady state. We also vary the bond stiffnesses to investigate the influence of changes in salt conditions. The key result is that the confined chain in the nanochannel pushed by a nano-dozer will progressively become highly knotted with delocalized knots. The knotting probability is greatly enhanced compared to corresponding equilibrium simulations, which in addition to compactification can be traced back to a backfolding mechanism for semi-flexible chains. Next, we describe the model, some essential facts about the LD simulation scheme, how our coarse-grained chains can be mapped onto DNA and the method that we use to analyze knots. Results Emergence of knots in nanochannels. Figure 1 summarizes the main findings of our study. Applying a pushing force leads to a compactification of the polymers (Fig. 1a), which in turn dramatically increases the occurrence of knots in the steady state in comparison to equilibrium values (Fig. 1b). Likewise, the amount of trefoil knots is reduced for configurations with a higher total knotting probability because the high density induces the formation of multiple or more complex knots. In this compact state, knots are delocalized and span over the whole chain as indicated for the example of trefoil knots in Fig. 1c where for a bending stiffness κ > 20 the average length of the knot is approximately 80% of the contour length or higher, which implicates that knots are formed preferentially near each end (please see Fig. 2e), while knots in equilibrium conformations are sig- nificantly smaller. These findings in a sense mirror previous observations, e.g. in Ref.29, which demonstrated that a θ-transition from a swollen coil to a globular state is not only accompanied by an increase in knotting but also by a delocalization of the latter. Figure 1d investigates the influence of the piston velocity for the experimentally relevant case of κ = 4 (DNA in a nanochannel, see “Methods” section). Again, compactification with increasing velocity is directly related to an increase of overall knotting. These results suggest that the occurrences of knots can be tuned by the speed Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 2 Vol:.(1234567890) www.nature.com/scientificreports/ Figure 2. Simulated structures for different parameters. (a, b) Structures for κ = 4 (a) (length parallel to the channel 73σ ) and κ = 20 (b) (length 91σ ) visualized using VMD49. Beads are colored from blue to red according to their monomer number. (c, d) The plots show the average position along the tube for all beads averaged over a simulation time of 50, 000τ for simulations at κ = 4 (c) and κ = 20 (d) with the piston at position 0. The insets show the relative density of beads along the structure averaged over the same frames. The structure for κ = 20 is significantly larger. The highest density for both structures is found to be close to the piston on the right. (e) Equilibrium structure without piston for κ = 20. of the piston and converge towards the equilibrium values for small piston velocities (Fig. 1e). This outcome is similar to results of Michieletto et al.38 indicating that an increased piston force will lead to an increase in the overall knotting probability and knot complexity as shown by a decrease in the occurrence of simple trefoil knots. As indicated above the decrease of knotting towards the equilibrium state at slow piston velocities is again accompanied by a trend towards a weak localization of trefoil knots (Fig. 1f)29. Figure 2 sheds light on these findings from a molecular basis. For κ = 4 the structure is disordered but the position of the monomers is still correlated with their sequence as indicated by the color scheme in Fig. 2a and the bead positions in Fig. 2c. For κ = 20 , the persistence length already exceeds the width of the tube which in conjunction with compactification leads to backfolding (Fig. 2b,d). The backfolding on the other hand creates loops which are a prerequisite for knots and in turn explains the initial rise in knotting with κ as well as their delocalization. In this context it would be interesting to study if backfolding creates a prevalence of torus-type knots as e.g. observed in the DNA located in viral capsids4,7. Unfortunately, the statistics of our simulations do not allow for a meaningful comparison. For large persistence lengths, backfolding becomes more difficult resulting in a lower knotting probability (Fig. 1b). In the equilibrium case, the compactification from the piston is no longer present and for κ = 20 , the chain can already spread throughout the channel which leads to a low knotting probability and weakly localized knots (Fig. 2e). Discussion In this manuscript we investigate velocity induced knot “production” in a nanochannel in comparison to those under equilibrium conformations. Both knotting probability and knot sizes depend strongly on piston velocity and resulting compactification as well as chain stiffness which can be, e.g., mitigated by adjusting ionic conditions and screening of charges in DNA. We observe that if the chain’s persistence length is greater than the width of the nanochannel, knots form by a backfolding mechanism. Since backfolding becomes harder for larger stiffness, the probability of knot formation decreases which explains the observed non-monotonic characteristic of knot formation in a nanochannel. We also study relative occurrences of complex knots as a function of the piston velocity and the chain stiffness. Our study thus sheds some new light on recent experiments in which DNA knots were created in a flow channel16 and provides insight on the molecular origin and control of self-entanglements under these conditions. Finally, we would like to point out that the coarse-grained simulation does not include hydrodynamic effects. It is worth considering how the results will change if we had incorporated it in the simulation. Dorfman has argued that for flexible chains hydrodynamic interactions are important. But for the semi-flexible chains with persistence length larger than the pore width, the chain is fully extended and is described by the free-draining limit50. Thus, for most channel sizes which result in a significant extension of the DNA compared to its bulk conformation, the hydrodynamic interactions between segments of the chain are almost completely screened. For the parameters used here the chain conformation lies in the transition region between deGennes blobs and Odjik limit and there is no theoretical argument for the effects of the hydrodynamic interaction in this regime. However, from Fig. 2 we observe that the chain conformations are mostly extended and hydrodynamic effects are likely to be small. For large velocities folded conformations are very different from deGennes blobs and therefore, one would expect the hydrodynamic effects to be small also, and the conclusions of this manuscript will essentially remain the same. Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 3 Vol.:(0123456789) www.nature.com/scientificreports/ Figure 3. Schematics of the simulation model. (a) A semi-flexible bead-spring chain confined inside a rectangular nano-channel is pushed by the green piston from right to left at velocity v = vpiston . The confinement potentials are imposed on four (two xy and two xz) planes, and by the moving piston in the yz plane in the negative x direction. The chain is free to move on the side opposite the piston. This figure was created using Xfig (version 3.2.8a, URL: https:// sourc eforge. net/p roje cts/ mcj/). (b) Demonstration of the knot closure. The termini of the polymer are connected with the center of mass (black dot) as indicated by the dashed red lines. The solid red lines are then appended to the polymer and connected by a closing arc drawn with red √ dots. (c, d) The root-mean-squared radius of gyration 〈R2g (τ )〉 of the polymer as a function of LD time τ for 10 simulations each and two different values of the chain stiffness, κ = 4 and 20, respectively. The polymer is compressed until it reaches a steady state with a constant radius of gyration. At this point it moves like a blob of √ a fixed shape. The steady state 〈R2g 〉 for the polymer with lower stiffness is significantly smaller in its final state. (e) The knotting probability as a function of time τ . During compression, the knotting probability increases while it is constant in the steady state. A higher stiffness leads to a higher knotting probability in the steady state. The results are averaged over 20 independent runs each. Methods Coarse-grained polymer model. The coarse-grained polymer model for LD simulations used here is exactly the same as in our previous publication47 where a bead-spring model polymer chain is confined to an open-ended rectangular channel and pushed from the right with a piston in the negative x-direction (Fig. 3a). The semi-flexible chain (Fig. 3b) is represented by a generalized bead-spring model51 where the beads (mono- mers) interact via excluded volume (EV), a Finite Extension Nonlinear Elastic (FENE) spring potential and a bond-bending potential enabling variation of ℓp as implemented p reviously46,47. The excluded volume interac- tion between any two monomers i and j of diameter σ is given by a short range truncated and shifted Lennard- Jones (LJ) potential ULJ (Eq. 1) of strength ε with a cutoff distance r = 21/6c σ is given by: [ ] N ( )12 ( )6 ∑ σ σ U = ε − + ε r ≤ 1/6LJ 4 for ij 2 σ rij ri ij 1/60 for ij 2 σ , where rij = |�ri − �rj| is the distance between any pair of beads. Successive monomers are connected by a FENE spring potential N−1 1 ∑ ( ) U 2FENE = − kFENER0 ln 1− r 2 2 2 i,i+1 /R0 , (2) i where kFENE is the spring constant and R0 is the maximum allowed bond length. The parameters kFENE and R0 along with ε and σ determine the bond-length. The chain stiffness is controlled by a bond-bending potential N−1 ∑ Ubend = κ (1− cos θi). (3) i=2 ( ) Here −1 b � i−1·b�θi = cos i is the angle between two successive bond vectors b�i−1 = �ri − �ri−1 and b�i = �ri+1 − �ri , |b�i−1||b�i | respectively, as shown in Fig. 3b. In three dimensions, for κ = 0 , the persistence length ℓp of the chain is related to κ via39,52,53 κσ ℓp≈ , (4) kBT for the values of κ considered in this work where kB is the Boltzmann constant and T is the temperature. Langevin dynamics simulation. We use the following Langevin dynamics equation of motion to advance the position of the ith monomer Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 4 Vol:.(1234567890) www.nature.com/scientificreports/ mr̈i = −∇(ULJ + UFENE + Ubend + Uwall + Upiston)+ γ vi +Wi , (5) where γ is the monomer friction coefficient, and Wi is a Gaussian random force with zero mean at temperature T which satisfies the fluctuation-dissipation relation. The numerical integration is implemented using the algo- rithm introduced by Gunsteren and Berendsen54. Our previous experience with LD simulations suggests that √ appropriate parameter specifications are γ = 0.7 mε/σ 2 , kFENE = 30ε/σ , R0 = 1.5σ , and a temperature kBT/ε = 1.2 . For a time step �t = 0.01τ (with τ being the standard Lennard-Jones time) these parameter values produce stable trajectories over a very long period of time and do not lead to an unphysical crossing of a bond by a m onomer55,56. The average bond length stabilizes at bl = 0.970± 0.002 with negligible fluctuation regardless of chain size and r igidity55. The piston is moved with a constant velocity of v0 = 0.005 σ if not noted otherwise τ after an initial equilibration of the chain. We ensure that the MD time for the pushing phase is long enough for the chain to attain a steady state shown in Fig. 3c–e that displays the connection between chain extension (Fig. 3c,d) and knot formation (Fig. 3e) in approach to the steady state. Times to reach the latter depend on bond √ stiffness κ as seen from the behavior of 〈 〈R2g (t)〉 in Fig. 3c,d. While for κ = 4 reaching a steady state takes less than 50,000τ , it takes around 160,000τ for κ = 20 . For each κ and v, physical quantities are averaged over at least ten independent runs. Our analysis indicates that knots can form and dissolve both during the initial compres- sion and in the steady state. Even compact structures can still change their knot type in accordance with other results on knot formation under applied force38. Of course, the probability of forming a knot, however, is signifi- cantly larger in the compressed (steady) state. Reptation Monte Carlo simulation to study the equilibrium limit. Note that piston speeds in coarse-grained implicit solvent simulations are typically orders of magnitude faster when compared to experi- ments. Therefore, we have also undertaken reptation Monte Carlo simulations of a slightly simplified model system of a single semi-flexible bead-spring chain confined inside a rectangular nanochannel of fixed size. In the reptation move which resembles the movement of a slithering snake, one segment is deleted from a ran- domly chosen chain end and attached to the other e nd57. Moves are accepted based on the Metropolis criterion. The repulsive Lennard-Jones and bond-bending potentials were matched with those of the LD simulation as described above. However, contrary to the LD simulation model’s confinement potential imposed onto the tube walls we use non-interacting walls and fixed bond length bl = 0.967 . These MC simulations allow for a compari- son of our dynamical investigations with equilibrium values (corresponding to piston velocity v → 0). Knot analysis. Knots in a closed chain are typically characterized by the minimum number of crossings observed when projecting a 3D chain onto a plane and can be considered as a fine gauge for the overall structure. Apart from the unknotted ring, the so-called unknot, the simplest knot is the trefoil (31) knot, which contains three crossings. There is one knot type with four crossings (41) and two with five crossings, and from there on the number of different knots with the same number of crossings increases e xponentially58,59. In our setup the polymer chain is open, and therefore, a closure connecting both ends of the chain has to be defined. First, we connect the end-points of each polymer with its center of mass. Along these lines we define a closure which emerges from one terminus follows the first line connects to the second one far away from the polymer and ends at the second t erminus40. After closure, the Alexander polynomial can be determined as described in detail i n60 (compare Fig. 3b). Knot sizes are determined by successively removing monomers from the ends of a polymer until the knot type changes29. Mapping onto DNA and comparison with experiments. Mapping our semi-flexible chain onto DNA is based on Eq. (4). For simplicity, we assume a solvent-independent persistence length of 50 nm or 150 base pairs. Furthermore, we assume that our beads describe the locus of a double-stranded DNA strand. In high salt conditions (1M NaCl), charges of DNA are completely screened and σ ≈ 2.5 nm. In physiological conditions charges are only partially screened and σ ≈ 5 nm, and for low salt conditions σ increases even further to about 15 nm at 0.01 M NaCl33,39,61. With a simulation temperature of T = 1.2 used throughout we obtain (in simulation units) κ = 24 for high salt, κ = 12 for physiological and κ = 4 for low salt conditions. This allows us to put our simulations in the context of recent experiments by Amin et al.16 undertaken at an estimated ionic strength of 8 mM which corresponds to our low salt scenario. Our chain has a contour length of L = Nσ = 1024σ = 15,360 nm or 46,080 base pairs, while our confining tube has a width of 16σ ≈ 240 nm. This compares to 168,903 base pairs and tube dimensions of 325× 415 nm used in Ref.16. Note that the mapping changes drastically with ionic conditions. The time scale of the simulated quantities can be translated into an experimental time scale via √ m tMD = σ . (6) kBT For low salt conditions we use σ = 15  nm as stated above, a mass m of 618u per base pair, a persistence length of κσk T = 3.33 beads and therefore 45 base pairs per bead. We assume room temperature of T = 300 . This results in aB time scale of 1 simulation time which equals to approximately 1.6ns. Our simulation time of 400,000τ for κ = 4 is therefore equivalent to 6.4× 10−4 s. Hence without explicit solvent one obtains a much faster time scale compared to experiments which often take place at a scale of seconds16. A typical experimental piston velocity for these experiments is 0.1–1 µm/s14,15. In the simulation we advance the piston with a velocity v0 = 0.005 σ ≈ 0.05 τ m/s, several order of magnitude faster than the experimental piston speed. Thus dynamics in our coarse-grained Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 5 Vol.:(0123456789) www.nature.com/scientificreports/ simulation are accelerated by several orders of magnitude in comparison to experiments and cannot be com- pared directly. Received: 13 December 2021; Accepted: 17 March 2022 References 1. Marenduzzo, D., Micheletti, C. & Orlandini, E. Biopolymer organization upon confinement. J. Phys. Condens. Matt. 22, 283102 (2010). 2. Reisner, W., Pedersen, J. N. & Austin, R. H. DNA confinement in nanochannels: Physics and biological applications. Rep. Prog. Phys. 75, 106601 (2012). 3. Arsuaga, J., Vazquez, M., Trigueros, S., Sumners, D. W. & Roca, J. Knotting probability of DNA molecules confined in restricted volumes: DNA knotting in phage capsids. PNAS 99, 5373–5377 (2002). 4. Arsuaga, J. et al. DNA knots reveal a chiral organization of DNA in phage capsids. PNAS 102, 9165–9169 (2005). 5. Petrov, A. S. & Harvey, S. C. Packaging double-helical DNA into viral capsids: Structures, forces, and energetics. Biophys. J. 95, 497 (2008). 6. Marenduzzo, D. et al. DNA–DNA interactions in bacteriophage capsids are responsible for the observed DNA knotting. PNAS 106, 22269–22274 (2009). 7. Reith, D., Cifra, P., Stasiak, A. & Virnau, P. Effective stiffening of DNA due to nematic ordering causes DNA molecules packed in phage capsids to preferentially form torus knots. Nucleic Acids Res. 40, 5129–5137 (2012). 8. Marenduzzo, D., Micheletti, C., Orlandini, E. & Sumners, D. W. Topological friction strongly affects viral DNA ejection. PNAS 110, 20081–20086 (2013). 9. Berndsen, Z. T., Keller, N., Grimes, S., Jardine, P. J. & Smith, D. E. Nonequilibrium dynamics and ultraslow relaxation of confined DNA during viral packaging. PNAS 111, 8345 (2014). 10. Keller, N., Grimes, S., Jardine, P. J. & Smith, D. E. Single DNA molecule jamming and history-dependent dynamics during motor- driven viral packaging. Nat. Phys. 12, 757 (2016). 11. Kim, Y. et al. Nanochannel confinement: DNA stretch approaching full contour length. Lab Chip 11, 1721 (2011). 12. Stein, D., van der Heyden, H. J., Koopmas, W. J. A. & Dekker, C. Pressure-driven transport of confined DNA polymers in fluidic channels. PNAS 3, 15853 (2006). 13. Mikkelsen, M. B., Reisner, W., Flyvbjerg, H. & Kristensen, A. Pressure-driven DNA in nanogroove arrays: Complex dynamics leads to length- and topology-dependent separation. Nano Lett. 11, 1598 (2011). 14. Khorshid, A. et al. Dynamic compression of single nanochannel confined DNA via a nanodozer assay. Phys. Rev. Lett. 113, 268104 (2014). 15. Khorshid, A., Amin, S., Zhang, Z., Sakaue, T. & Reisner, W. Non-equilibrium dynamics of nanochannel confined dna. Macromol- ecules 49, 1933 (2016). 16. Amin, S., Khorshid, A., Zeng, L., Zimny, P. & Reisner, W. A nanofluidic knot factory based on compression of single DNA in nanochannels. Nat. Commun. 9, 1506–1506 (2018). 17. Lam, T. E. et al. Genome mapping on nanochannel arrays for structural variation analysis and sequence assembly. Nat. Biotechnol. 30, 771 (2012). 18. Zhou, J. et al. Enhanced nanochannel translocation and localization of genomic DNA molecules using three-dimensional nano- funnels. Nat. Commun. 8, 807 (2017). 19. Plesa, C. et al. Direct observation of DNA knots using a solid-state nanopore. Nat. Nanotechnol. 11, 1093–1097 (2016). 20. Suma, A. & Micheletti, C. Pore translocation of knotted DNA rings. Proc. Natl. Acad. Sci. 114, E2991–E2997 (2017). 21. Sakaue, T. Semiflexible polymer confined in closed spaces. Macromolecules 40, 5206 (2007). 22. Jun, S., Thirumalai, D. & Ha, B.-Y. Compression and stretching of a self-avoiding chain in cylindrical nanopores. Phys. Rev. Lett. 101, 138101 (2008). 23. Micheletti, C., Marenduzzo, D. & Orlandini, E. Polymers with spatial or topological constraints: Theoretical and computational results. Phys. Rep. 504, 1–73 (2011). 24. Vologodskii, A. V., Lukashin, A. V. & Frank-Kamenetskii, M. D. Topological interaction between polymer chains. Sov. Phys. JETP 40, 932–936 (1974). 25. Koniaris, K. & Muthukumar, M. Knottedness in ring polymers. Phys. Rev. Lett. 66, 2211–2214 (1991). 26. Mansfield, M. L. Knots in Hamilton cycles. Macromolecules 27, 5924–5926 (1994). 27. Grosberg, A. Y. Critical exponents for random knots. Phys. Rev. Lett. 85, 3858–3861 (2000). 28. Marcone, B., Orlandini, E., Stella, A. L. & Zonta, F. What is the length of a knot in a polymer?. J. Phys. Math. Gen. 38, L15–L21 (2005). 29. Virnau, P., Kantor, Y. & Kardar, M. Knots in globule and coil phases of a model polyethylene. JACS 127, 15102 (2005). 30. Foteinopoulou, K., Karayiannis, N. C., Laso, M., Kröger, M. & Mansfield, M. L. Universal scaling, entanglements, and knots of model chain molecules. Phys. Rev. Lett. 101, 265702 (2008). 31. Orlandini, E. & Micheletti, C. Knotting of linear DNA in nano-slits and nano-channels: A numerical study. J. Biol. Phys. 39, 267–275 (2013). 32. Micheletti, C. & Orlandini, E. Knotting and unknotting dynamics of DNA strands in nanochannels. ACS Macro Lett. 3, 876–880 (2014). 33. Trefz, B., Siebert, J. & Virnau, P. How molecular knots can pass through each other. PNAS 111, 7948 (2014). 34. Coronel, L., Orlandini, E. & Micheletti, C. Non-monotonic knotting probability and knot length of semiflexible rings: The compet- ing roles of entropy and bending energy. Soft Matter 13, 4260 (2017). 35. Meyer, H., Horwath, E. & Virnau, P. Mapping onto ideal chains overestimates self-entanglements in polymer melts. ACS Macro Lett. 7, 757–761 (2018). 36. Zhang, J., Meyer, H., Virnau, P. & Daoulas, K. C. Can soft models describe polymer knots?. Macromolecules 53, 10475–10486 (2020). 37. Tubiana, L. et al. Comparing equilibration schemes of high-molecular-weight polymer melts with topological indicators. J. Phys. Condens. Matter. 33, 204003 (2021). 38. Michieletto, D., Orlandini, E., Turner, M. S. & Micheletti, C. Separation of geometrical and topological entanglement in confined polymers driven out of equilibrium. ACS Macro Lett. 9, 1081–1085 (2020). 39. Rieger, F. C. & Virnau, P. A Monte Carlo study of knots in long double-stranded DNA chains. PLoS Comput. Biol 12, e1005029 (2016). 40. Virnau, P., Mirny, L. & Kardar, M. Intricate knots in proteins: Function and evolution. PLoS Comput. Biol. 2, e122 (2006). 41. Sulkowska, J. I., Sulkowski, P., Szymczak, P. & Cieplak, M. Stabilizing effect of knots on proteins. PNAS 105, 19714–19719 (2008). 42. Sulkowska, J. I., Sulkowski, P. & Onuchic, J. Dodging the crisis of folding proteins with knots. PNAS 106, 3119–3124 (2009). 43. Boelinger, D. et al. A stevedore’s protein knot. PLoS Comput. Biol. 6, e1000731 (2010). Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 6 Vol:.(1234567890) www.nature.com/scientificreports/ 44. Jamroz, M. et al. Knotprot: A database of proteins with knots and slipknots. Nucleic Acids Res. 43, D306–D314 (2015). 45. Wüst, T., Reith, D. & Virnau, P. Sequence determines degree of knottedness in a coarse-grained protein model. Phys. Rev. Lett. 114, 028102 (2015). 46. Huang, A., Reisner, W. & Bhattacharya, A. Dynamics of DNA squeezed inside a nanochannel via a sliding gasket. Polymers 8, 352 (2016). 47. Bernier, S., Huang, A., Reisner, W. & Bhattacharya, A. Evolution of nested folding states in compression of a strongly confined semiflexible chain. Macromoecules 51, 4012 (2018). 48. Hayase, Y., Sakaue, T. & Nakanishi, H. Compressive response and helix formation of a semiflexible polymer confined in a nano- channel. Phys. Rev. E. 95, 052502 (2017). 49. Humphrey, W., Dalke, A. & Schulten, K. VMD—Visual molecular dynamics. J. Mol. Graph 14, 33 (1996). 50. Dorfman, K. D., Gupta, D., Jain, A., Muralidhar, A. & Tree, D. R. Hydrodynamics of DNA confined in nanoslits and nanochannels. Eur. Phys. J. Spec. Top. 223, 3179–3200 (2014). 51. Grest, G. S. & Kremer, K. Molecular dynamics simulation for polymers in the presence of a heat bath. Phys. Rev. A 33, 3628–3631 (1986). 52. Hsu, H.-P. & Binder, K. Stretching semiflexible polymer chains: Evidence for the importance of excluded volume effects from Monte Carlo simulation. J. Chem. Phys. 136, 024901 (2012). 53. Fixman, M. & Kovac, J. Polymer conformational statistics. III. Modified Gaussian models of stiff chains. J. Chem. Phys. 58, 1564–1568 (1973). 54. van Gunsteren, W. F. & Berendsen, H. J. C. Algorithms for Brownian dynamics. Mol. Phys. 45, 637 (1982). 55. Huang, A., Bhattacharya, A. & Binder, K. Conformations, transverse fluctuations, and crossover dynamics of a semi-flexible chain in two dimensions. J. Chem. Phys. 140, 214902 (2014). 56. Huang, A., Hsu, H.-P., Bhattacharya, A. & Binder, K. Semiflexible macromolecules in quasi-one-dimensional confinement: Discrete versus continuous bond angles. J. Chem. Phys. 143, 243102 (2015). 57. Wall, F. T. & Mandel, F. Macromolecular dimensions obtained by an efficient Monte Carlo method without sample attrition. J. Chem. Phys. 63, 4592–4595 (1975). 58. Rolfsen, D. Knots and Links. Mathematics Lecture Series (Publish or Perish, 1976). 59. Adams, C. The Knot Book (W.H. Freeman, 1994). 60. Virnau, P. Detection and visualization of physical knots in macromolecules. Phys. Proc. 6, 114 (2010). 61. Rybenkov, V. V., Cozzarelli, N. R. & Vologodskii, A. V. Probability of DNA knotting and the effective diameter of the DNA double helix. PNAS 90, 5307–5311 (1993). Acknowledgements The authors acknowledge partial funding from TopDyn. We are grateful to the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for funding this research: Project Number 233630050-TRR 146. The LD simulation were carried out at the UCF’s high performance computing cluster STOKES. The authors gratefully acknowledge the computing time granted on the supercomputer Mogon offered by the Johannes Gutenberg Uni- versity Mainz (hpc.uni-mainz.de), which is a member of the AHRP (Alliance for High Performance Computing in Rhineland Palatinate, http:// www. ahrp. info) and the Gauss Alliance e.V.. A.B. acknowledges travel support from Institut für Physik, Mainz. A.B. thanks Aiqun Huang for the preliminary runs. Author contributions A.B. performed the LD simulations. S.W. performed the Monte Carlo simulations. J.R and P.V. analyzed the results. All authors reviewed the manuscript. Funding Open Access funding enabled and organized by Projekt DEAL. Competing interests The authors declare no competing interests. Additional information Correspondence and requests for materials should be addressed to P.V. or A.B. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://c reati veco mmons.o rg/ licens es/ by/4. 0/. © The Author(s) 2022 Scientific Reports | (2022) 12:5342 | https://doi.org/10.1038/s41598-022-09242-5 7 Vol.:(0123456789) pubs.acs.org/macroletters Letter Can Polymer Helicity Affect Topological Chirality of Polymer Knots? Yani Zhao, Jan Rothörl, Pol Besenius, Peter Virnau,* and Kostas Ch. Daoulas* Cite This: ACS Macro Lett. 2023, 12, 234−240 Read Online ACCESS Metrics & More Article Recommendations *sı Supporting Information ABSTRACT: We investigate the effect of helicity in isolated polymers on the topological chirality of their knots with computer simulations. Polymers are described by generic worm-like chains (WLC), where helical conformations are promoted by chiral coupling between segments that are neighbors along the chain contour. The sign and magnitude of the coupling coefficient u determine the sense and strength of helicity. The corrugation of the helix is adjusted via the radius R of a spherical, hard excluded volume around each WLC segment. Open and compact helices are, respectively, obtained for R that is either zero or smaller than the length of the WLC bond, and R that is a few times larger than the bond length. We use a Monte Carlo algorithm to sample polymer conformations for different values of u, spanning the range from achiral polymers to chains with well-developed helices. Monitoring the average helix torsion and fluctuations of chiral order as a function of u, for two very different chain lengths, demonstrates that the coil−helix transition in this model is not a phase transition but a crossover. Statistical analysis of conformations forming the simplest chiral knots, 31, 51, and 52, demonstrates that topological mirror symmetry is broken�knots formed by helices with a given sense prefer one handedness over the other. For the 31 and 51 knots, positive helical sense favors positive handedness. Intriguingly, an opposite trend is observed for 52 knots, where positive helical sense promotes negative handedness. We argue that this special coupling between helicity and topological chirality stems from a generic mechanism: conformations where some of the knot crossings are found in “braids” formed by two tightly interwoven sections of the polymer. Molecular chirality can be classified into1−3 geometrical, achiral 41 knot). Since the closure is imaginary, knots evolve,chemical, and topological chirality. Geometrical chir- disappear, and re-emerge, as linear polymers sample their ality applies to molecular objects that cannot1−3 be super- conformational space.6 For a polymer with chemically achiral imposed with their mirror image by translation and/or rotation conformations (such as the bead−spring polymer in Figure 1), operations. Chemical chirality incorporates the influence of one expects that knots and their mirror images are equally actual dynamics: a molecular configuration is chemically chiral represented in the accessible conformational space. In other when it cannot3 be deformed into its mirror image by words, the mixture of knotted conformations is racemic. intramolecular transformations that are physically feasible Conversely, the presence of chemical chirality might favor left- under the given conditions, e.g., temperature. Topological or right-handed knots, breaking the mirror symmetry of the chirality applies to molecular configurations that cannot1−3 be knot mixture. transformed into their mirror image by a continuous Although this hypothesis sounds reasonable, the actual deformation. Here, actual dynamics is irrelevant: bond lengths, knowledge on connections between chemical and topological valence angles, and dihedral angles can change arbitrarily, as chirality in polymers is very limited. Among others, the long as the deformed molecule does not intersect itself. influence of chirality on interactions between knots 7 and their 8 Intriguingly, polymer chains form structures that are most capability to meander through helical channels have been natural to study topological chirality: knots. Figure 1a presents studied, and recently, a 819 knot with controlled handedness the simplest possible knot, the 3 (known as trefoil) knot, has been created artificially. 9 Many studies have investigated 1 created, for illustration, on a generic bead−spring chain. knotted structures in biopolymers such as proteins 10−15 and Formally, knots are defined only in closed loops, but the concept is applicable to linear polymers after introducing an Received: October 13, 2022 imaginary closure, indicated in Figure 1a by dashed lines. The Accepted: January 23, 2023 closed 31 knot is topologically chiral because it cannot 1,4,5 be continuously transformed into its mirror image, Figure 1b. Generic cartoons in Figure 1c present the next two simplest 51 and 52 (known as fivefold) knots, that are chiral (we omit the © XXXX The Authors. Published by American Chemical Society https://doi.org/10.1021/acsmacrolett.2c00600 234 ACS Macro Lett. 2023, 12, 234−240 Downloaded via UNIV MAINZ on January 30, 2023 at 09:48:14 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. ACS Macro Letters pubs.acs.org/macroletters Letter In this work, we consider a special, but very basic, example of chemical chirality: polymers that form right- or left-handed helices (Figure 1a and b explain the definition of a left- and right-handed helix). Our goal is to find generic relationships between helical sense and topological chirality of polymer knots. We use a generic coarse-grained model to describe isolated helical polymer chains found in Θ solvent (ideal chains) and good solvent conditions. Polymers are represented by worm- like chains (WLC) with N segments (bonds). The interactions are defined through a Hamiltonian expressed in units of thermal energy kBT as N 1 N d i+d H = t ·t + u r k T b i i+1 ij ·[ti × tj] B i=1 i=1 j=i+1 N [4R/b] N + U(2R |rij|) i=1 j=i+[4R/b] (1) We arbitrarily choose one of the two directions along the contour of the WLC to number the segments. Accordingly, ti is a unit vector oriented along the i-th segment of the WLC and rij = rij/|rij|, where rij is the distance vector between the center- of-mass (COM) of the i-th and the j-th segment (rij is oriented Figure 1. Panels (a) and (b) show, respectively, a conformation of an from the i-th toward the j-th COM). achiral bead−spring chain forming a left-handed 3 knot and its right- The first term in eq 1 manipulates polymer stiffness by1 handed mirror image. Panel (a) explains the closure: two lines (green adjusting the parameter ϵb. The length of the segments is fixed dashed lines) are determined by extrapolating to infinity two line to b. The second term is a chiral potential that couples only segments (dashed purple lines), connecting the center-of-mass (green those segments that are separated, along the WLC contour, by circle) of the conformation with the two ends of the chain. The less than d + 1 segments and favors conformations with helical closure is accomplished by connecting the two extrapolated lines far twist of prescribed sense. The magnitude of u controls the from the chain (green dashed line). The two-dimensional projection strength of the chiral coupling and its sign determines the of the knots on the xy-plane (gray) is sketched below. Red, green, and helical sense. Namely, u > 0 and u < 0 result, respectively, in gray mark the three arcs of each projected knot. The ± signs indicate the handedness of each crossing. Panel (a) illustrates the unit vectors right-handed and left-handed helices, whereas u = 0 along z-axis ez, overpassing bond o, and underpassing bond u, used to corresponds to an achiral polymer. Similar chiral interactions determine the handedness of one representative crossing. The insets are used in studies of cholesteric phase, 24,25 chiral block of panels (a) and (b) explain, respectively, the definition of a left- and copolymers, 26,27 simulations of chiral liquid crystals,28 chiral 29 right-handed helix. They show a top view of a helix which rotates into aggregates, coil−helix transitions,30 and Go-like models of the page along the direction of the red arrow. Panel (c) shows generic proteins.31 The third term assigns to each segment a hard sketches of the next two simplest chiral knots, 51 and 52 (both excluded volume with radius R, centered at its COM. examples show left-handed knots). Specifically, U(2R − |rij|) = +∞ when 2R − |rij| ≥ 0 and U(2R − |rij|) = 0 otherwise. There are no excluded volume DNA.16−22 While confining viral DNA in a capsid increases the interactions between segments that are separated by less than knotting probability (≈95% for the 10 kilobase pairs (kbp) [4R/b] segments along the WLC contour (angular brackets DNA strand of the P4 phage vir1 del22 tailless mutants16,17), define the integer-part function). We sample the conforma- knotting can be substantial even without confinement, tional space using a Monte Carlo (MC) reptation algorithm.32 provided that chains are long: a 166 kbp phage T4 GT7 Details are provided in the Supporting Information. DNA contains a knot in 70% of all cases.21 Although Exploring the behavior of the generic model across the four- biopolymers are a canonical example of chemically chiral dimensional space of parameters ϵb, u, d, and R is outside the macromolecules, research effort concentrated on detection and scope of this study. We choose three subsets of parameter classification of knots, without correlating their topology with space that are interesting for studying the behavior of knots. molecular-level features.18 In particular, protein knots with They have the same ϵb = 4 and d = 10 but differ regarding the both positive and negative handedness have been ob- size of the excluded volume: R = 0, 0.5, and 2.5, respectively (R served.14,15 Intriguingly, chiral DNA knots have been is given in units of b). u is a free parameter. constructed as early as 1995 by an enzymatic closure reaction For R = 0 and 0.5 the chosen parameters lead to open using a left-handed Z-DNA to craft trefoils with positive helices where the pitch p is substantially larger33 than the handedness and right-handed B-DNA for trefoils with negative excluded volume of the segments, as illustrated in the main handedness.23 However, because the knotted state was made panel of Figure 2. This snapshot stems from MC simulations permanent by the cyclic (closed) molecular architecture, the with u = 0.5. Studying knots in open helical polymers for both question regarding which knot handedness is thermodynami- R = 0 and 0.5 is important because R = 0.5 retains some cally favored by a right- or left-handed polymer chain was out excluded volume. This situation might be more straightforward of the scope of that study. to realize in experiments (we provide further discussion later 235 https://doi.org/10.1021/acsmacrolett.2c00600 ACS Macro Lett. 2023, 12, 234−240 ACS Macro Letters pubs.acs.org/macroletters Letter Figure 2.Main panel: Snapshot of an open WLC helix with N = 2000, ϵb = 4, d = 10, and u = 0.5 obtained for small radius of excluded volume R = 0.5. The inset presents a section of a compact WLC helix obtained for the same N, ϵb, d, and u but with a large radius of excluded volume R = 2.5. on). For R = 2.5, we obtain compact helices. The inset of Figure 2 presents a small part of a compact helix generated during MC simulations at u = 0.5. Before exploring the behavior of knots, we must understand how helicity changes as a function of u. We quantify the helicity of chains using the ensemble-averaged torsion ⟨τ⟩ of chains as an order parameter, defined as34 1 N 1 (ri × ri )·r= i | × |2 Figure 3. (a) Average torsion ⟨τ⟩ shown as a function of u for chainN 3 i=3 ri ri (2) lengths N = 1000 and 2000 (dashed and solid lines). For each N, data corresponding to R = 0, 0.5, and 2.5 are presented, as indicated by the Here ri is the position vector of the i-th “monomer”; a WLC legends. (b) Variance of chiral energy Cc (normalized by u) presentedwith N segments has N + 1 “monomers”, i.e., N − 1 internal for polymers from panel (a) as a function of u. The color code is the junctions and two ends. The derivatives are discretized34 as r′i same as in panel (a). The insets show snapshots of a WLC with N = = 1 (ri + 1 − ri − 1), r″ 1i = ri + 1 − 2ri + ri − 1, ri‴ = (ri + 2 − 1000, at u = 0.1 (left) and 0.5 (right). To visualize the helical2 2 − − structure more clearly, parts of the polymer are shown enlarged in the2(ri + 1 ri − 1) ri − 2). Figure 3a presents ⟨τ⟩ as a function of dashed frames. u, for open, R = 0 and 0.5, and compact, R = 2.5, helices. For each R we present plots for two representative chain lengths, N = 1000 and 2000. Error bars are estimated from the standard transition.44 Some studies35,45 have emphasized that explaining error of the mean and are smaller than symbol size. the phenomenology of mirror-symmetry breaking in single Interestingly, ⟨τ⟩ grows smoothly as u becomes larger, polymers through the thermodynamics of 1D systems might be irrespective of R. Moreover, the plots of ⟨τ⟩ for the two an oversimplification because, typically, there are long-range different chain lengths are, practically, on top of each other; interactions between repeat units. In principle, the excluded that is, we do not observe finite system-size effects. These volume interaction in eq 1 correlates segments that are far observations suggest that in our model the transition from apart along the WLC contour.46 However, the trends in Figure achiral to chiral state is not a phase transition but a continuous 3a suggest that these correlations are not sufficient to cause a crossover.35 Formally, mirror symmetry is broken even for very phase transition in our case. small u, although in this limit chirality is weak. We quantify the strength of fluctuations in chiral order One expects36 that phase transition is absent when the through the variance of chiral energy, a heat capacity-like breaking of mirror symmetry in isolated chains is caused by property, normalized by u: chiral intramolecular interactions between a limited number of N d i+d consecutive monomers. In eq 1 the chiral potential couples d C = 1 [ E2 E 2], E = r ·[t × t ] consecutive segments only. Such potentials resemble spin- c N c c c ij i j couplings in Ising-like models of helical polymers.37−42 They i=1 j=i+1 (3) create equivalence36,43 to an effective one-dimensional (1D) Figure 3b presents Cc calculated for the same systems as in system with local interactions, which cannot exhibit a phase Figure 3a. Error bars are again small; we estimate them from 236 https://doi.org/10.1021/acsmacrolett.2c00600 ACS Macro Lett. 2023, 12, 234−240 ACS Macro Letters pubs.acs.org/macroletters Letter the expression of the standard error of the mean for variances (“fluctuation averages”).47 Figure 3b demonstrates that the fluctuations in chirality are stronger for small u. Cc exhibits a peak near u = 0.1, which is, however, weak. Importantly, the magnitude of the peak and the overall shape of Cc do not depend on N (the plots for N = 1000 and 2000 are practically on top of each other). This phenomenology of Cc is consistent with an absence of a phase transition. Visual inspection reveals that the stronger fluctuations at low u are associated with more “fluffy” conformations, exhibiting weakly helical and molten nonhelical regions. In contrast, helices are well formed for large u and dominate chain conformations. Figure 3b presents two snapshots of a chain with N = 1000 and R = 0.5, taken at u = 0.1 and 0.5. We can now analyze knots. Generally, one can define the handedness of a knot by considering the minimal projection of the knot onto a plane.48 Then, the handedness of a single crossing is defined as h = ez · ([o × u])/|[o × u]|, where ez, o, and u are, respectively, the unit vector along the z-axis, overpassing bond (at the projected crossing), and under- passing bond. The sign of the sum of all h in a minimal projection determines the handedness of the knot.48,49 For example, the minimal projection of a trefoil has only three essential crossings (see Figure 1). The sum of h is positive for right-handed and negative for left-handed knots. In practice, we calculate HOMFLYPT polynomials.50,51 First, chain conformations, generated by MC, are “closed” by using the closure49 explained in Figure 1. This closure enables calculations of HOMFLYPT polynomials with the Topoly package52 to determine knot type and handedness. Figure 4a presents the probability P Figure 4. (a) The main panel shows the total knotting and occurrencek (black lines) to find a probability, both indicated by P , of 3 knots as a function of u, for knotted conformation in a N = 2000 chain, as a function of u k 1chains with R = 0 and N = 2000. Error bars are smaller than the width for R = 0 (main panel), 0.5 (left inset), and 2.5 (right inset). of the line. The left and right insets show the total knotting and The plots demonstrate that Pk is very small for chains with an occurrence probability for 31 knots for R = 0.5 and 2.5, respectively. excluded volume and highlight the difficulties in collecting data (b) Occurrence probability Pk,norm of 41, 51, and 52 knots for chains for analyzing knot handedness in this case. For R = 0 and 0.5, with R = 0 and N = 2000 normalized by the total knotting probability. the Pk has a clear maximum at u = 0, whereas for R = 2.5, the (c) Probability Prk that a knot formed at given u is right handed. Prk = plot (despite the large error bars) suggests a nonmonotonous 50% (horizontal black dashed line) indicates no preferred handedness. dependence, i.e., P has a maximum at u ≠ 0. Increasing u Data for 31 knots obtained with R = 0, 0.5, and 2.5, and for fivefoldk makes the chains stiffer, so both trends in Figure 4a are knots obtained with R = 0 are presented, as indicated near each plot. consistent with dependencies of Pk on chain stiffness that have been reported for achiral ideal53 (monotonous decay) and regarding how Prk changes as a function of u for chains without achiral self-avoiding54,55 (nonmonotonous decay) chains. In and with excluded volume. Overall, we observe that for R = 0 Figure 4a, for R = 0, we separately show the probability of the effect of helicity on handedness of knots is much stronger observing a 31 (red line) knot, which is by far the most than for R > 0, that is, deviations from Prk = 50% are more common knot type at this length scale. For R = 0.5 and 2.5, the pronounced. Furthermore, Prk increases monotonously for R = preponderance of 31 knots is even stronger. 0 to saturate (at least for the considered range of u) to a Figure 4b displays the probability Pk,norm of observing 41 constant value. In contrast, for chains with excluded volume, (orange line), 51 (blue line), and 52 (green line) knots, the effect of helicity on knots is stronger for small u. normalized by the total knotting probability. We use this Specifically, Prk exhibits a peak near u = 0.1 and decays after normalization to emphasize the fraction of individual knots. that. Importantly, for the largest excluded volume R = 2.5, the The share of torus knots 31 (not shown here) and 51 increases excess (depletion) of right-handed knots is smaller than for R = with increasing |u|, while the portion of nontorus knots such as 0.5. 41 or 52 is approximately constant for small |u| and decreases Prk shows broken mirror symmetry also for 51 and 52 knots. for larger |u|. This trend resembles results obtained for DNA Because it is challenging to accumulate reliable statistics on 51 confined to bacteriophage capsids16,17 and simulations and 52 knots for chains with an excluded volume (see insets in describing such systems.18,19 Figure 4a), Figure 4c presents results only for R = 0. The Figure 4c is central for our work. It presents the probability dependence of Prk on u for 51 knots (blue line) qualitatively Prk that 31, 51, and 52 knots, when formed at given u, are right- reproduces the trends observed for 31 knots: we see handed. For 31 knots (red line), the mirror symmetry is clearly pronounced excess (depletion) of right-handed knots for u > broken for all three R considered in our study, because there is 0 (u < 0), which saturates at high u. This behavior is consistent an excess and depletion of right-handed knots for u > 0 and u < with 51 knot being a torus knot which (for open chains) can be 0, respectively. However, there is a qualitative difference obtained from a 31 knot by one extra winding around the knot 237 https://doi.org/10.1021/acsmacrolett.2c00600 ACS Macro Lett. 2023, 12, 234−240 ACS Macro Letters pubs.acs.org/macroletters Letter contour. Intriguingly, however, 52 knots (green line) show an synthesis of molecular knots. 58 Second, we take into account opposite trend; there is a surplus of right-handed knots for u < that in interwoven helices segments come close to each other. 0 and a depletion for u > 0. This, at first glance, unexpected Therefore, for R = 0, we eliminate from the sample knotted behavior of 52 knots demonstrates that the handedness of conformations where at least two segments (separated by more knots may not coincide with the helical sense of the molecule. than d segments along the WLC contour) are found closer We suggest that the coupling between helicity and than a cutoff distance rmin; a typical choice is r = 1.5. Fortopological chirality found in our simulations stems (to large minknots that survive this screening (and therefore have no tightly extent, at least) from a generic mechanism. Namely, it is caused by conformations where some of the knot crossings are packed braids), the deviations of Prk from 50% are significantly encapsulated in a “braid” formed by two interwoven helical reduced. Plots are available in the Supporting Information. subchains. Figure 5a illustrates three representative conforma- Consistent with the effect of screening for R = 0, we observe a reduced deviation of Prk from 50% in systems with excluded volume, especially R = 2.5 (Figure 4c). In the latter case, compact helices do not allow for molecular interdigitation sufficient for forming braids. Finally, we note that for R = 0.5 and 2.5 there is also a similarity between the nonmonotonous dependence of Prk on u and the behavior of Cc in Figure 3b. This observation suggests that strong fluctuations in local chain conformations and helicity promote interdigitation. Various studies59−63 have revealed that special packing of molecules with helical surface affects mesoscopic chiral order in multichain systems. In this respect, our observations regarding the relationship between local packing of helices and topological chirality of knots are not surprising. Still, it is rather unexpected that the fraction of conformations found in this particular knotted state is sufficiently large to cause perceptible mirror-symmetry breaking for the entire set of 31, Figure 5. (a) Snapshots of parts of an N = 2000 chain forming 3 , 5 , 51, and 52 knots. Detailed analysis of knot handedness versus1 1 and 5 knots (from left to right) for u = 0.5 and R = 0, illustrating the amount of monomers contained in a knot demonstrates2 typical shapes of knots with braids (for these parameters). The that broken mirror symmetry is more pronounced for smaller braided part of the knots is shown enlarged, in a separate frame. (b) knots. Hence, it is plausible to expect that the effects of helicity Illustrations of the same knots based on an idealized physical (wire) on topological chirality of the entire population of knots in model. The chain contour is traced from end A toward end B. Red longer chains (than those that have been considered here) will and blue arrows, respectively, mark the direction of motion along the first and the second polymer strand forming the braid. The topology be reduced. of 3 and 5 knots is such that the two strands are traveled in the same Our findings are based on a generic molecular model but can1 1 direction, whereas for the 52 knot, they are traveled in opposite be extrapolated to actual helical polymers. We expect that an directions. This difference affects the handedness of the crossings, excess of knots with one sense of handedness might be which is indicated by the ± signs. observed in chiral polymers where helices have well-separated “ridges” and “valleys”, the analog of open helices formed in our tions of 31, 51, and 52 knots with such a braid (the braided part model at a small excluded volume. Polyisocyanates64,65 might is enlarged, in a separate frame) for R = 0 (for clarity, we are be one example, taking into account that the formation of their showing only the knotted part of an N = 2000 chain). The sense (direction) of winding of the subchains around each lyotropic cholesteric phases can be explained 62 assuming a other is the same as the sense of the polymer helix (positive in strongly corrugated, screw-like, helical molecular surface. the examples of Figure 5a). Another candidate are biopolymers with polyproline helices Initially, our explanation is motivated by visual inspection of of type PP-II. In contrast, we do not expect strong preferred knotted conformations. Of course, visual analysis cannot be handedness for knots in polymers with compact helices. Here, systematic because of the significant amount of knotted representative examples are biopolymers with polyproline conformations and their variability. However, there are also helices of type PP-I or α-helices.66 The topology of single several more quantitative arguments that favor our conjecture knotted polymer conformations can be analyzed by modern based on indirect evidence. First of all, helices, say, with imaging techniques such AFM.67 positive helical sense, interwoven into a braid with positive twist can simultaneously explain the preference for positive ASSOCIATED CONTENT handedness in 31 and 51 and the negative handedness in 5 ■2 knots, respectively. Figure 5b provides explanatory illustrations *sı Supporting Information based on an idealized physical (wire) model of helical knotted The Supporting Information is available free of charge at chains. These illustrations may also indicate that for large |u| https://pubs.acs.org/doi/10.1021/acsmacrolett.2c00600. torus knots, 31 and 51 are easier to form and have a simpler braiding pattern than, for example, a 5 knot, providing an Details on the applied reptation algorithm; The effect of2 interpretation for results observed in Figure 4b. A similar exclusion of conformations with closely packed segments analysis of knots in terms of “braids” was performed56,57 to on handedness of trefoil knots; The effect of polymer explain the occurrence of certain knot types in template helicity on the prevalent type of fivefold knots (PDF) 238 https://doi.org/10.1021/acsmacrolett.2c00600 ACS Macro Lett. 2023, 12, 234−240 ACS Macro Letters pubs.acs.org/macroletters Letter ■ AUTHOR INFORMATION (11) Virnau, P.; Mirny, L. A.; Kardar, M. Intricate Knots in Proteins: Corresponding Authors Function and Evolution. PLoS Comput. Biol. 2006, 2, No. e122.(12) Lua, R. C.; Grosberg, A. Y. Statistics of Knots, Geometry of Peter Virnau − Department of Physics, Johannes Gutenberg Conformations, and Evolution of Proteins. PLoS Comp. Biol. 2006, 2, University Mainz, 55128 Mainz, Germany; Email: virnau@ 350−357. uni-mainz.de (13) Boelinger, D.; Sulkowska, J.; Hsu, H.-P.; Mirny, L.; Kardar, M.; Kostas Ch. Daoulas − Max Planck Institute for Polymer Onuchic, J.; Virnau, P. A Stevedore’s Protein Knot. PLOS Comp. Biol. Research, 55128 Mainz, Germany; orcid.org/0000-0001- 2010, 6, e1000731. 9278-6036; Phone: +49 (0)6131 379218; (14) Virnau, P.; Mallam, A.; Jackson, S. Structures and Folding Email: daoulas@mpip-mainz.mpg.de Pathways of Topologically Knotted Proteins. J. Phys: Cond. Matt. 2011, 23, 033101. Authors (15) Jarmolinska, A. I.; Perlinska, A. P.; Runkel, R.; Trefz, B.; Ginn, Yani Zhao − Max Planck Institute for Polymer Research, H. M.; Virnau, P.; Sulkowska, J. I. Proteins’ Knotty Problems. J. Mol. 55128 Mainz, Germany; orcid.org/0000-0003-1430- Biol. 2019, 431, 244−257. 4518 (16) Arsuaga, J.; Vazquez, M.; Trigueros, S.; Sumners, D. W.; Roca, Jan Rothörl − Department of Physics, Johannes Gutenberg J. Knotting Probability of DNA Molecules Confined in Restricted University Mainz, 55128 Mainz, Germany Volumes: DNA Knotting in Phage Capsids. Proc. Natl. Acad. Sci. Pol Besenius − Department of Chemistry, Johannes Gutenberg U.S.A. 2002, 99, 5373−5377. University Mainz, 55128 Mainz, Germany; orcid.org/ (17) Arsuaga, J.; Vazquez, M.; McGuirk, P.; Trigueros, S.; Sumners,D. W.; Roca, J. DNA Knots Reveal a Chiral Organization of DNA in 0000-0001-7478-4459 Phage Capsids. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 9165−9169. Complete contact information is available at: (18) Marenduzzo, D.; Orlandini, E.; Stasiak, A.; Sumners, D. W.; https://pubs.acs.org/10.1021/acsmacrolett.2c00600 Tubiana, L.; Micheletti, C. DNA−DNA Interactions in Bacteriophage Capsids are Responsible for the Observed DNA Knotting. Proc. Natl. Funding Acad. Sci. U.S.A. 2009, 106, 22269−22274. Open access funded by Max Planck Society. (19) Reith, D.; Cifra, P.; Stasiak, A.; Virnau, P. Effective Stiffening ofDNA due to Nematic Ordering Causes DNA Molecules Packed in Notes Phage Capsids to Preferentially Form Torus Knots. Nucleic Acids Res. The authors declare no competing financial interest. 2012, 40, 5129−5137. ■ (20) Rieger, F. C.; Virnau, P. A Monte Carlo Study of Knots in LongACKNOWLEDGMENTS Double-Stranded DNA Chains. PLoS Comp. Biol. 2016, 12, e1005029.(21) Plesa, C.; Verschueren, D.; Pud, S.; van der Torre, J.; We are grateful to Kurt Kremer for helpful discussions related Ruitenberg, J. W.; Witteveen, M. J.; Jonsson, M. P.; Grosberg, A. Y.; to this work and useful comments made after reading our Rabin, Y.; Dekker, C. Direct Observation of DNA Knots using a manuscript. We thank Wanda Niemyska and Pawel Rubach for Solid-State Nanopore. Nat. Nanotechnol. 2016, 11, 1093−1097. helping us with using the Topoly software for knot analysis. (22) Kumar Sharma, R.; Agrawal, I.; Dai, L.; Doyle, P. S.; Garaj, S. P.V. and K.C.D. acknowledge funding from the Deutsche Complex DNA Knots Detected with a Nanopore Sensor. Nat. Forschungsgemeinschaft (DFG, German Research Founda- Commun. 2019, 10, 4473. tion): Project Number 233630050 - TRR 146. (23) Du, S.; Stollar, B.; Seeman, N. A Synthetic DNA Molecule in 3 ■ Knotted Topologies. J. Am. Chem. Soc. 1995, 117, 1194−1200.REFERENCES (24) van der Meer, B.; Vertogen, G.; Dekker, A. J.; Ypma, J. G. J. AMolecular-Statistical Theory of the Temperature-Dependent Pitch in (1) Mislow, K. A Commentary on the Topological Chirality and − Cholesteric Liquid Crystals. J. Chem. Phys. 1976, 65, 3935−3943.Achirality of Molecules. Croat. Chem. Acta 1996, 69, 485 511. (25) Osipov, M. In Liquid Crystalline and Mesomorphic Polymers; (2) Chambron, J. C.; Dietrich-Buchecker, C.; Sauvage, J. P. From Classical Chirality to Topologically Chiral Catenands and Knots. Top. Shibaev, V. P., Lam, L., Eds.; Springer, 1994; pp 1−25. Curr. Chem. 1993, 165, 131−162. (26) Zhao, W.; Russell, T. P.; Grason, G. M. Chirality in Block (3) Bonchev, D.; Rouvray, D. H. Chemical Topology: Applications and Copolymer Melts: Mesoscopic Helicity from Intersegment Twist. Techniques; Gordon and Breach Science Publishers: Amsterdam, Phys. Rev. Lett. 2013, 110, 058301. 2000. (27) Grason, G. M. Chirality Transfer in Block Copolymer Melts: (4) Dehn, M. Die beiden Kleeblattschlingen. Math. Ann. 1914, 75, Emerging Concepts. ACS Macro Lett. 2015, 4, 526−532. 402−413. (28) Memmer, R.; Kuball, H.-G.; Schönhofer, A. Computer (5) Fielden, S. D.; Leigh, D. A.; Woltering, S. L. Molecular Knots. Simulation of Chiral Liquid Crystal Phases I. The Polymorphism of Angew. Chem., Int. Ed. 2017, 56, 11166−11194. the Chiral Gay-Berne Fluid. Liq. Cryst. 1993, 15, 345−360. (6) Tubiana, L.; Rosa, A.; Fragiacomo, F.; Micheletti, C. (29) Sutherland, B. J.; Olesen, S. W.; Kusumaatmaja, H.; Morgan, J. Spontaneous Knotting and Unknotting of Flexible Linear Polymers: W. R.; Wales, D. J. Morphological Analysis of Chiral Rod Clusters Equilibrium and Kinetic Aspects. Macromolecules 2013, 46, 3669− from a Coarse-Grained Single-Site Chiral Potential. Soft Matter 2019, 3678. 15, 8147−8155. (7) Najafi, S.; Tubiana, L.; Podgornik, R.; Potestio, R. Chirality (30) Kemp, J.; Chen, Z. Formation of Helical States in Wormlike Modifies the Interaction between Knots. EPL 2016, 114, 50007. Polymer Chains. Phys. Rev. Lett. 1998, 81, 3880−3883. (8) Ruskova, R.; Racko, D. Channels with Helical Modulation (31) Wolék, K.; Gómez-Sicilia, Á.; Cieplak, M. Determination of Display Stereospecific Sensitivity for Chiral Superstructures. Polymers Contact Maps in Proteins: a Combination of Structural and Chemical 2021, 13, 3726. Approaches. J. Chem. Phys. 2015, 143, 243105. (9) Carpenter, J. P.; McTernan, C. T.; Greenfield, J. L.; (32) Wall, F.; Mandel, F. Macromolecular Dimensions Obtained by Lavendomme, R.; Ronson, T. K.; Nitschke, J. R. Controlling the an Efficient Monte Carlo Method without Sample Attrition. J. Chem. Shape and Chirality of an Eight-Crossing Molecular Knot. Chem. Phys. 1975, 63, 4592−4595. 2021, 7, 1534−1543. (33) Glagolev, M. K.; Vasilevskaya, V. V.; Khokhlov, A. R. (10) Taylor, W. A. Deeply Knotted Protein Structure and How it Compactization of Rigid-Chain Amphiphilic Macromolecules with Might Fold. Nature 2000, 406, 916−919. Local Helical Structure. Polymer Science, Ser. A 2010, 52, 761−774. 239 https://doi.org/10.1021/acsmacrolett.2c00600 ACS Macro Lett. 2023, 12, 234−240 ACS Macro Letters pubs.acs.org/macroletters Letter (34) Magee, J. E.; Song, Z.; Curtis, R. A.; Lue, L. Structure and (58) Ayme, J.-F.; Beves, J. E.; Campbell, C. J.; Leigh, D. A. Template Aggregation of a Helix-Forming Polymer. J. Chem. Phys. 2007, 126, synthesis of molecular knots. Chem. Soc. Rev. 2013, 42, 1700−1712. 144911. (59) Ferrarini, A.; Moro, G.; Nordio, P. Shape Model for Ordering (35) Boehm, C. R.; Terentjev, E. M. Minimal Model of Intrinsic Properties of Molecular Dopants Inducing Chiral Mesophases. Mol. Chirality to Study the Folding Behavior of Helical Polymers. Phys. 1996, 87, 485−499. Macromolecules 2014, 47, 6086−6094. (60) De Michele, C.; Zanchetta, G.; Bellini, T.; Frezza, E.; Ferrarini, (36) Grosberg, A. Y.; Khokhlov, A. R. Statistical Physics of A. Hierarchical Propagation of Chirality through Reversible Polymer- Macromolecules; AIP Press: New York, 1994. ization: The Cholesteric Phase of DNA Oligomers. ACS Macro Lett. (37) Selinger, J. V.; Selinger, R. L. B. Theory of Chiral Order in 2016, 5, 208−212. random Copolymers. Phys. Rev. Lett. 1996, 76, 58−61. (61) Straley, J. Theory of Piezoelectricity in Nematic Liquid (38) Green, M. M.; Park, J.-W.; Sato, T.; Teramoto, A.; Lifson, S.; Crystals, and of the Cholesteric Order. Phys. Rev. A 1976, 14, 1835− Selinger, R. L. B.; Selinger, J. V. The Macromolecular Route to Chiral 1841. Amplification. Angew. Chem., Int. Ed. 1999, 38, 3138−3154. (62) Green, M. M.; Peterson, N. C.; Sato, T.; Teramoto, A.; Cook, (39) van Gestel, J.; van der Schoot, P.; Michels, M. A. J. R.; Lifson, S. A Helical Polymer with a Cooperative Response to Amplification of Chirality in Helical Supramolecular Polymers Chiral Information. Science 1995, 268, 1860−1866. Beyond the Long-Chain Limit. J. Chem. Phys. 2004, 120, 8253−8261. (63) Earl, D. J.; Wilson, M. R. Predictions of Molecular Chirality and (40) van Gestel, J. Amplification of Chirality in Helical Supra- Helical Twisting Powers: A Theoretical Study. J. Chem. Phys. 2003, molecular Polymers: The Majority-Rules Principle. Macromolecules 119, 10280−10288. 2004, 37, 3894−3898. (64) Green, M. M.; Reidy, M. P.; Johnson, R. D.; Darling, G.; (41) Jouvelet, B.; Isare, B.; Bouteiller, L.; van der Schoot, P. Direct O’Leary, D. J.; Willson, G. Macromolecular Stereochemistry: the Out- Probing of the Free-Energy Penalty for Helix Reversals and Chiral of-Proportion Influence of Optically Active Comonomers on the Mismatches in Chiral Supramolecular Polymers. Langmuir 2014, 30, Conformational Characteristics of Polyisocyanates. The Sergeants and 4570−4575. Soldiers Experiment. J. Am. Chem. Soc. 1989, 111, 6452. (42) van Gestel, J.; Palmans, A. R. A.; Titulaer, B.; Vekemans, J. A. J. (65) Sato, T.; Sato, Y.; Umemura, Y.; Teramoto, A.; Nagamura, Y.;Wagner, J.; Weng, D.; Okamoto, Y.; Hatada, K.; Green, M. M. M.; Meijer, E. W. ”Majority-Rules” Operative in Chiral Columnar Polyisocyanates and the Interplay of Experiment and Theory in the Stacks of C3-Symmetrical Molecules. J. Am. Chem. Soc. 2005, 127, Formation of Lyotropic Cholesteric States. Macromolecules 1993, 26, 5490−5494. 4551−4559. (43) Zimm, B. H.; Bragg, J. K. Theory of the Phase Transition (66) Pauling, L.; Corey, R. B.; Branson, H. R. The Structure of between Helix and Random Coil in Polypeptide Chains. J. Chem. Proteins: two Hydrogen-Bonded Helical Configurations of the Phys. 1959, 31, 526−535. Polypeptide Chain. Proc. Natl. Acad. Sci. U.S.A. 1951, 37, 205−211. (44) Landau, L. D.; Lifshitz, E. M. Statistical Physics; Pergamon: (67) Schappacher, M.; Deffieux, A. Imaging of Catenated, Figure-of- London, 1959; Vol. 5. Eight, and Trefoil Knot Polymer Rings. Angew. Chem. Int. Ed. 2009, (45) Hansmann, U. H.; Okamoto, Y. Finite-Size Scaling of Helix- 48, 5930−5933. Coil Transitions in Poly-Alanine Studied by Multicanonical Simulations. J. Chem. Phys. 1999, 110, 1267−1276. (46) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Recommended by ACS Claredon Press: Oxford, England, 1986. (47) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1989. Fast and Integral Nano-Surface-Coating of Various Fiber (48) Livingston, C. Knot Theory; Mathematical Association of Materials via Interfacial Polymerization America: WA, 1993. Aqiang Wang, Jian Jin, et al. (49) Virnau, P. Detection and Visualization of Physical Knots in JANUARY 03, 2023 Macromolecules. Phys. Procedia 2010, 6, 117−125. ACS MACRO LETTERS READ (50) Freyd, P.; Yetter, D.; Hoste, J.; Lickorish, W. B. R.; Millett, K.; Ocneanu, A. A New Polynomial Invariant of Knots and Links. Bull. High-Temperature Wetting and Dewetting Dynamics of Am. Math. Soc. 1985, 12, 239−246. Silver Droplets on Molybdenum Surfaces (51) Przytycki, J.; Traczyk, P. Invariants of Links of Conway Type. Kobe J. Math 1987, 4, 115−139. Lin Lin, Duu-Jong Lee, et al. (52) Dabrowski-Tumanski, P.; Rubach, P.; Niemyska, W.; Gren, B. JANUARY 09, 2023LANGMUIR A.; Sulkowska, J. I. Topoly: Python package to analyze topology of READ polymers. Briefings in Bioinformatics 2021, 22, bbaa196. (53) Virnau, P.; Rieger, F. C.; Reith, D. Influence of Chain Stiffness Influence of the Fabric Topology on the Performance of a on Knottedness in Single Polymers. Biochem. Soc. Trans. 2013, 41, Textile-Based Triboelectric Nanogenerator for Self-Powered 528−532. Monitoring (54) Coronel, L.; Orlandini, E.; Micheletti, C. Non-Monotonic Viraj U. Somkuwar and Bipin Kumar Knotting Probability and Knot Length of Semiflexible Rings: the JANUARY 09, 2023 Competing Roles of Entropy and Bending Energy. Soft Matter 2017, ACS APPLIED POLYMER MATERIALS READ 13, 4260−4267. (55) Uehara, E.; Coronel, L.; Micheletti, C.; Deguchi, T. Bimodality Photo-Thermal Superhydrophobic Sponge for Highly in the Knotting Probability of Semiflexible Rings Suggested by Efficient Anti-Icing and De-Icing Mapping with Self-Avoiding Polygons. React. Funct. Polym. 2019, 134, Bo Yu, Feng Zhou, et al. 141−149. JANUARY 15, 2023 (56) Polles, G.; Marenduzzo, D.; Orlandini, E.; Micheletti, C. Self- LANGMUIR READ assembling knots of controlled topology by designing the geometry of patchy templates. Nat. Commun. 2015, 6, 6423. (57) Marenda, M.; Orlandini, E.; Micheletti, C. Discovering privileged topologies of molecular knots with self-assembling models. Get More Suggestions > Nat. Commun. 2018, 9, 3051. 240 https://doi.org/10.1021/acsmacrolett.2c00600 ACS Macro Lett. 2023, 12, 234−240 B. Bibliography B. Bibliography 1A. N. Bogdanov and D. A. Yablonskiǐ, “Thermodynamically stable ”vortices” in magnetically ordered crystals. The mixed state of magnets,” Sov. Phys. JETP 68, 101–103 (1989). 2A. N. Bogdanov and A. Hubert, “Thermodynamically stable magnetic vortex states in magnetic crystals,” J. Magn. Magn. Mater. 138, 255–269 (1994). 3U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, “Spontaneous skyrmion ground states in magnetic metals,” Nature 442, 797–801 (2006). 4S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, “Skyrmion lattice in a chiral magnet,” Science 323, 915–919 (2009). 5X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, “Real-space observation of a two-dimensional skyrmion crystal,” Nature 465, 901–904 (2010). 6X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura, “Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe,” Nat. Mater. 10, 106–109 (2011). 7S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blügel, “Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions,” Nat. Phys. 7, 713–718 (2011). 8G. Finocchio, F. Büttner, R. Tomasello, M. Carpentieri, and M. Kläui, “Magnetic skyrmions: from fundamental to applications,” J. Phys. D: Appl. Phys. 49, 423001 (2016). 9K. Everschor-Sitte, J. Masell, R. M. Reeve, and M. Kläui, “Perspective: magnetic skyrmions—overview of recent progress in an active research field,” J. Appl. Phys. 124, 240901 (2018). 10I. Dzyaloshinsky, “A thermodynamic theory of ”weak” ferromagnetism of anti- ferromagnetics,” J. Phys. Chem. Solids 4, 241–255 (1958). 11T. Moriya, “Anisotropic superexchange interaction and weak ferromagnetism,” Phys. Rev. 120, 91–98 (1960). 12A. Fert and P. M. Levy, “Role of anisotropic exchange interactions in determining the properties of spin-glasses,” Phys. Rev. Lett. 44, 1538–1541 (1980). 132 B. Bibliography 13S. Rohart and A. Thiaville, “Skyrmion confinement in ultrathin film nanostruc- tures in the presence of Dzyaloshinskii-Moriya interaction,” Phys. Rev. B 88, 184422 (2013). 14N. Nagaosa and Y. Tokura, “Topological properties and dynamics of magnetic skyrmions,” Nat. Nanotechnol. 8, 899–911 (2013). 15R. Lo Conte, A. Hrabec, A. P. Mihai, T. Schulz, S.-J. Noh, C. H. Marrows, T. A. Moore, and M. Kläui, “Spin-orbit torque-driven magnetization switching and thermal effects studied in Ta\CoFeB\MgO nanowires,” Appl. Phys. Lett. 105, 122404 (2014). 16R. Lo Conte, E. Martinez, A. Hrabec, A. Lamperti, T. Schulz, L. Nasi, L. Lazzarini, R. Mantovan, F. Maccherozzi, S. S. Dhesi, B. Ocker, C. H. Marrows, T. A. Moore, and M. Kläui, “Role of B diffusion in the interfacial Dzyaloshinskii- Moriya interaction in Ta/Co20Fe60B20/MgO nanowires,” Phys. Rev. B 91, 014433 (2015). 17S. Jaiswal, K. Litzius, I. Lemesh, F. Büttner, S. Finizio, J. Raabe, M. Weigand, K. Lee, J. Langer, B. Ocker, G. Jakob, G. S. D. Beach, and M. Kläui, “Investigation of the Dzyaloshinskii-Moriya interaction and room temperature skyrmions in W/CoFeB/MgO thin films and microwires,” Appl. Phys. Lett. 111, 022409 (2017). 18W. Jiang, G. Chen, K. Liu, J. Zang, S. te Velthuis, and A. Hoffmann, “Skyrmions in magnetic multilayers,” Phys. Rep. 704, 1–49 (2017). 19A. N. Bogdanov and C. Panagopoulos, “Physical foundations and basic properties of magnetic skyrmions,” Nat. Rev. Phys. 2, 492–498 (2020). 20F. Zheng, H. Li, S. Wang, D. Song, C. Jin, W. Wei, A. Kovács, J. Zang, M. Tian, Y. Zhang, H. Du, and R. E. Dunin-Borkowski, “Direct imaging of a zero-field target skyrmion and its polarity switch in a chiral magnetic nanodisk,” Phys. Rev. Lett. 119, 197205 (2017). 21S.-Z. Lin, C. Reichhardt, C. D. Batista, and A. Saxena, “Particle model for skyrmions in metallic chiral magnets: dynamics, pinning, and creep,” Phys. Rev. B 87, 214419 (2013). 22J. Zázvorka, F. Jakobs, D. Heinze, N. Keil, S. Kromin, S. Jaiswal, K. Litzius, G. Jakob, P. Virnau, D. Pinna, K. Everschor-Sitte, L. Rózsa, A. Donges, U. Nowak, and M. Kläui, “Thermal skyrmion diffusion used in a reshuffler device,” Nat. Nanotechnol. 14, 658–661 (2019). 23R. E. Troncoso and Á. S. Núñez, “Brownian motion of massive skyrmions in magnetic thin films,” Ann. Phys. 351, 850–856 (2014). 133 B. Bibliography 24J. Miltat, S. Rohart, and A. Thiaville, “Brownian motion of magnetic domain walls and skyrmions, and their diffusion constants,” Phys. Rev. B 97, 214426 (2018). 25F. Jonietz, S. Mühlbauer, C. Pfleiderer, A. Neubauer, W. Münzer, A. Bauer, T. Adams, R. Georgii, P. Böni, R. A. Duine, K. Everschor, M. Garst, and A. Rosch, “Spin transfer torques in MnSi at ultralow current densities,” Science 330, 1648–1651 (2010). 26X. Z. Yu, N. Kanazawa, W. Z. Zhang, T. Nagai, T. Hara, K. Kimoto, Y. Matsui, Y. Onose, and Y. Tokura, “Skyrmion flow near room temperature in an ultralow current density,” Nat. Commun. 3, 988 (2012). 27K. Everschor, M. Garst, B. Binz, F. Jonietz, S. Mühlbauer, C. Pfleiderer, and A. Rosch, “Rotating skyrmion lattices by spin torques and field or temperature gradients,” Phys. Rev. B 86, 054432 (2012). 28W. Jiang, P. Upadhyaya, W. Zhang, G. Yu, M. B. Jungfleisch, F. Y. Fradin, J. E. Pearson, Y. Tserkovnyak, K. L. Wang, O. Heinonen, S. G. E. te Velthuis, and A. Hoffmann, “Blowing magnetic skyrmion bubbles,” Science 349, 283–286 (2015). 29S. Woo, K. Litzius, B. Krüger, M.-Y. Im, L. Caretta, K. Richter, M. Mann, A. Krone, R. M. Reeve, M. Weigand, P. Agrawal, I. Lemesh, M.-A. Mawass, P. Fischer, M. Kläui, and G. S. D. Beach, “Observation of room-temperature magnetic skyrmions and their current-driven dynamics in ultrathin metallic ferromagnets,” Nat. Mater. 15, 501–506 (2016). 30J. M. Kosterlitz and D. J. Thouless, “Long range order and metastability in two dimensional solids and superfluids. (Application of dislocation theory),” J. Phys. C: Solid State Phys. 5, L124–L126 (1972). 31J. M. Kosterlitz and D. J. Thouless, “Ordering, metastability and phase transi- tions in two-dimensional systems,” J. Phys. C: Solid State Phys. 6, 1181–1203 (1973). 32B. I. Halperin and D. R. Nelson, “Theory of two-dimensional melting,” Phys. Rev. Lett. 41, 121–124 (1978). 33D. R. Nelson and B. I. Halperin, “Dislocation-mediated melting in two dimen- sions,” Phys. Rev. B 19, 2457–2484 (1979). 34A. P. Young, “Melting and the vector Coulomb gas in two dimensions,” Phys. Rev. B 19, 1855–1866 (1979). 35E. P. Bernard and W. Krauth, “Two-step melting in two dimensions: first-order liquid-hexatic transition,” Phys. Rev. Lett. 107, 155704 (2011). 134 B. Bibliography 36S. C. Kapfer and W. Krauth, “Two-dimensional melting: from liquid-hexatic coexistence to continuous transitions,” Phys. Rev. Lett. 114, 035702 (2015). 37A. L. Thorneywork, J. L. Abbott, D. G. A. L. Aarts, and R. P. A. Dullens, “Two-dimensional melting of colloidal hard spheres,” Phys. Rev. Lett. 118, 158001 (2017). 38Y.-W. Li and M. P. Ciamarra, “Accurate determination of the translational correlation function of two-dimensional solids,” Phys. Rev. E 100, 062606 (2019). 39K. Zahn, R. Lenke, and G. Maret, “Two-stage melting of paramagnetic colloidal crystals in two dimensions,” Phys. Rev. Lett. 82, 2721–2724 (1999). 40X. S. Wang, H. Y. Yuan, and X. R. Wang, “A theory on skyrmion size,” Commun. Phys. 1, 31 (2018). 41R. Gruber, M. A. Brems, J. Rothörl, T. Sparmann, M. Schmitt, I. Kononenko, F. Kammerbauer, M.-A. Syskaki, O. Farago, P. Virnau, and M. Kläui, “300- times-increased diffusive skyrmion dynamics and effective pinning reduction by periodic field excitation,” Adv. Mater. 35, 2208922 (2023). 42G. Qin, R. Zhang, C. Yang, X. Lv, K. Pei, L. Yang, X. Liu, X. Zhang, and R. Che, “Magnetic-field-assisted diffusion motion of magnetic skyrmions,” ACS Nano 16, 15927–15934 (2022). 43M. A. Brems, M. Kläui, and P. Virnau, “Circuits and excitations to enable Brownian token-based computing with skyrmions,” Appl. Phys. Lett. 119, 132405 (2021). 44J. Zázvorka, F. Dittrich, Y. Ge, N. Kerber, K. Raab, T. Winkler, K. Litzius, M. Veis, P. Virnau, and M. Kläui, “Skyrmion lattice phases in thin film multilayer,” Adv. Funct. Mater. 30, 2004037 (2020). 45A. Fert, N. Reyren, and V. Cros, “Magnetic skyrmions: advances in physics and potential applications,” Nat. Rev. Mater. 2, 17031 (2017). 46T. Nakajima, H. Oike, A. Kikkawa, E. P. Gilbert, N. Booth, K. Kakurai, Y. Taguchi, Y. Tokura, F. Kagawa, and T.-h. Arima, “Skyrmion lattice structural transition in MnSi,” Sci. Adv. 3, e1602562 (2017). 47P. Huang, T. Schönenberger, M. Cantoni, L. Heinen, A. Magrez, A. Rosch, F. Carbone, and H. M. Rønnow, “Melting of a skyrmion lattice to a skyrmion liquid via a hexatic phase,” Nat. Nanotechnol. 15, 761–767 (2020). 48A. V. Ognev, A. G. Kolesnikov, Y. J. Kim, I. H. Cha, A. V. Sadovnikov, S. A. Nikitov, I. V. Soldatov, A. Talapatra, J. Mohanty, M. Mruczkiewicz, Y. Ge, N. Kerber, F. Dittrich, P. Virnau, M. Kläui, Y. K. Kim, and A. S. Samardak, “Magnetic direct-write skyrmion nanolithography,” ACS Nano 14, 14960–14970 (2020). 135 B. Bibliography 49P. Baláž, M. Paściak, and J. Hlinka, “Melting of Néel skyrmion lattice,” Phys. Rev. B 103, 174411 (2021). 50P. Meisenheimer, H. Zhang, D. Raftrey, X. Chen, Y.-T. Shao, Y.-T. Chan, R. Yalisove, R. Chen, J. Yao, M. C. Scott, W. Wu, D. A. Muller, P. Fischer, R. J. Birgeneau, and R. Ramesh, “Ordering of room-temperature magnetic skyrmions in a polar van der Waals magnet,” Nat. Commun. 14, 3744 (2023). 51F. Büttner, C. Moutafis, M. Schneider, B. Krüger, C. M. Günther, J. Geilhufe, C. v. K. Schmising, J. Mohanty, B. Pfau, S. Schaffert, A. Bisig, M. Foerster, T. Schulz, C. A. F. Vaz, J. H. Franken, H. J. M. Swagten, M. Kläui, and S. Eisebitt, “Dynamics and inertia of skyrmionic spin structures,” Nat. Phys. 11, 225–228 (2015). 52B. L. Brown, U. C. Täuber, and M. Pleimling, “Effect of the Magnus force on skyrmion relaxation dynamics,” Phys. Rev. B 97, 020405 (2018). 53B. L. Brown, U. C. Täuber, and M. Pleimling, “Skyrmion relaxation dynamics in the presence of quenched disorder,” Phys. Rev. B 100, 024410 (2019). 54A. F. Schäffer, L. Rózsa, J. Berakdar, E. Y. Vedmedenko, and R. Wiesen- danger, “Stochastic dynamics and pattern formation of geometrically confined skyrmions,” Commun. Phys. 2, 72 (2019). 55D. Foster, C. Kind, P. J. Ackerman, J.-S. B. Tai, M. R. Dennis, and I. I. Smalyukh, “Two-dimensional skyrmion bags in liquid crystals and ferromagnets,” Nat. Phys. 15, 655–659 (2019). 56C. Reichhardt and C. J. O. Reichhardt, “Depinning and nonequilibrium dynamic phases of particle assemblies driven over random and ordered substrates: a review,” Rep. Prog. Phys 80, 026501 (2016). 57C. Reichhardt, C. J. O. Reichhardt, and M. V. Milošević, “Statics and dynamics of skyrmions interacting with disorder and nanostructures,” Rev. Mod. Phys. 94, 035005 (2022). 58A. A. Thiele, “Steady-state motion of magnetic domains,” Phys. Rev. Lett. 30, 230–233 (1973). 59R. Brearton, G. van der Laan, and T. Hesjedal, “Magnetic skyrmion interactions in the micromagnetic framework,” Phys. Rev. B 101, 134422 (2020). 60D. Capic, D. A. Garanin, and E. M. Chudnovsky, “Skyrmion–skyrmion interac- tion in a magnetic film,” J. Phys. Condens. Matter 32, 415803 (2020). 61Y. Wang, J. Wang, T. Kitamura, H. Hirakata, and T. Shimada, “Exponential temperature effects on skyrmion-skyrmion interaction,” Phys. Rev. Appl. 18, 044024 (2022). 136 B. Bibliography 62D. Reith, M. Pütz, and F. Müller-Plathe, “Deriving effective mesoscale potentials from atomistic simulations,” J. Comput. Chem. 24, 1624–1636 (2003). 63T. C. Moore, C. R. Iacovella, and C. McCabe, “Derivation of coarse-grained potentials via multistate iterative Boltzmann inversion,” J. Chem. Phys. 140, 224104 (2014). 64G. Milano and F. Müller-Plathe, “Mapping atomistic simulations to mesoscopic models: a systematic coarse-graining procedure for vinyl polymer chains,” J. Phys. Chem. B 109, 18609–18619 (2005). 65R. A. Pepper, M. Beg, D. Cortés-Ortuño, T. Kluyver, M.-A. Bisotti, R. Carey, M. Vousden, M. Albert, W. Wang, O. Hovorka, and H. Fangohr, “Skyrmion states in thin confined polygonal nanostructures,” J. Appl. Phys. 123, 093903 (2018). 66J. C. B. Souza, N. P. Vizarim, C. J. O. Reichhardt, C. Reichhardt, and P. A. Venegas, “Skyrmion ratchet in funnel geometries,” Phys. Rev. B 104, 054434 (2021). 67R. Gruber, J. Zázvorka, M. A. Brems, D. R. Rodrigues, T. Dohi, N. Kerber, B. Seng, M. Vafaee, K. Everschor-Sitte, P. Virnau, and M. Kläui, “Skyrmion pinning energetics in thin film systems,” Nat. Commun. 13, 3144 (2022). 68J. Iwasaki, M. Mochizuki, and N. Nagaosa, “Universal current-velocity relation of skyrmion motion in chiral magnets,” Nat. Commun. 4, 1463 (2013). 69Y.-H. Liu and Y.-Q. Li, “A mechanism to pin skyrmions in chiral magnets,” J. Phys. Condens. Matter 25, 076005 (2013). 70A. Derras-Chouk and E. M. Chudnovsky, “Skyrmions near defects,” J. Phys. Condens. Matter 33, 195802 (2021). 71C. Navau, N. Del-Valle, and A. Sanchez, “Interaction of isolated skyrmions with point and linear defects,” J. Magn. Magn. Mater. 465, 709–715 (2018). 72K. Litzius, I. Lemesh, B. Krüger, P. Bassirian, L. Caretta, K. Richter, F. Büttner, K. Sato, O. A. Tretiakov, J. Förster, R. M. Reeve, M. Weigand, I. Bykova, H. Stoll, G. Schütz, G. S. D. Beach, and M. Kläui, “Skyrmion Hall effect revealed by direct time-resolved X-ray microscopy,” Nat. Phys. 13, 170–175 (2017). 73W. Jiang, X. Zhang, G. Yu, W. Zhang, X. Wang, M. B. Jungfleisch, J. E. Pearson, X. Cheng, O. Heinonen, K. L. Wang, Y. Zhou, A. Hoffmann, and S. G. E. te Velthuis, “Direct observation of the skyrmion Hall effect,” Nat. Phys. 13, 162–169 (2017). 137 B. Bibliography 74K. Litzius, J. Leliaert, P. Bassirian, D. Rodrigues, S. Kromin, I. Lemesh, J. Zazvorka, K.-J. Lee, J. Mulkers, N. Kerber, D. Heinze, N. Keil, R. M. Reeve, M. Weigand, B. Van Waeyenberge, G. Schütz, K. Everschor-Sitte, G. S. D. Beach, and M. Kläui, “The role of temperature and drive current in skyrmion dynamics,” Nat. Electron. 3, 30–36 (2020). 75K. Zeissler, S. Finizio, C. Barton, A. J. Huxtable, J. Massey, J. Raabe, A. V. Sadovnikov, S. A. Nikitov, R. Brearton, T. Hesjedal, G. van der Laan, M. C. Rosamond, E. H. Linfield, G. Burnell, and C. H. Marrows, “Diameter-independent skyrmion Hall angle observed in chiral magnetic multilayers,” Nat. Commun. 11, 428 (2020). 76C. Wang, D. Xiao, X. Chen, Y. Zhou, and Y. Liu, “Manipulating and trapping skyrmions by magnetic field gradients,” New J. Phys. 19, 083008 (2017). 77C. Reichhardt, D. Ray, and C. J. O. Reichhardt, “Collective transport properties of driven skyrmions with random disorder,” Phys. Rev. Lett. 114, 217202 (2015). 78C. Reichhardt and C. J. O. Reichhardt, “Noise fluctuations and drive dependence of the skyrmion Hall effect in disordered systems,” New J. Phys. 18, 095005 (2016). 79S. A. Dı́az, C. J. O. Reichhardt, D. P. Arovas, A. Saxena, and C. Reichhardt, “Fluctuations and noise signatures of driven magnetic skyrmions,” Phys. Rev. B 96, 085106 (2017). 80C. Reichhardt, D. Ray, and C. J. O. Reichhardt, “Nonequilibrium phases and segregation for skyrmions on periodic pinning arrays,” Phys. Rev. B 98, 134418 (2018). 81C. Reichhardt and C. J. O. Reichhardt, “Thermal creep and the skyrmion Hall angle in driven skyrmion crystals,” J. Phys. Condens. Matter 31, 07LT01 (2018). 82C. Reichhardt and C. J. O. Reichhardt, “Nonlinear transport, dynamic ordering, and clustering for driven skyrmions on random pinning,” Phys. Rev. B 99, 104418 (2019). 83G. Chen, “Skyrmion Hall effect,” Nat. Phys. 13, 112–113 (2017). 84E. Kalz, H. D. Vuijk, I. Abdoli, J.-U. Sommer, H. Löwen, and A. Sharma, “Collisions enhance self-diffusion in odd-diffusive systems,” Phys. Rev. Lett. 129, 090601 (2022). 85C. J. O. Reichhardt and C. Reichhardt, “Active rheology in odd-viscosity sys- tems,” EPL 137, 66004 (2022). 86P. Rieger, M. Weißenhofer, and U. Nowak, “Defect-enhanced diffusion of magnetic skyrmions,” 10.48550/arXiv.2304.12917 (2023). 138 B. Bibliography 87A. Fert, V. Cros, and J. Sampaio, “Skyrmions on the track,” Nat. Nanotechnol. 8, 152–156 (2013). 88R. Tomasello, E. Martinez, R. Zivieri, L. Torres, M. Carpentieri, and G. Finocchio, “A strategy for the design of skyrmion racetrack memories,” Sci. Rep. 4, 6784 (2014). 89X. Zhang, G. P. Zhao, H. Fangohr, J. P. Liu, W. X. Xia, J. Xia, and F. J. Morvan, “Skyrmion-skyrmion and skyrmion-edge repulsions in skyrmion-based racetrack memory,” Sci. Rep. 5, 7643 (2015). 90J. Müller, “Magnetic skyrmions on a two-lane racetrack,” New J. Phys. 19, 025002 (2017). 91O. Lee, R. Msiska, M. A. Brems, M. Kläui, H. Kurebayashi, and K. Everschor- Sitte, “Perspective on unconventional computing using magnetic skyrmions,” Appl. Phys. Lett. 122, 260501 (2023). 92D. Pinna, F. Abreu Araujo, J.-V. Kim, V. Cros, D. Querlioz, P. Bessiere, J. Droulez, and J. Grollier, “Skyrmion gas manipulation for probabilistic comput- ing,” Phys. Rev. Appl. 9, 064018 (2018). 93K. Wang, Y. Zhang, V. Bheemarasetty, S. Zhou, S.-C. Ying, and G. Xiao, “Single skyrmion true random number generator using local dynamics and interaction between skyrmions,” Nat. Commun. 13, 722 (2022). 94R. Ishikawa, M. Goto, H. Nomura, and Y. Suzuki, “Implementation of skyrmion cellular automaton using brownian motion and magnetic dipole interaction,” Appl. Phys. Lett. 119, 072402 (2021). 95T. Nozaki, Y. Jibiki, M. Goto, E. Tamura, T. Nozaki, H. Kubota, A. Fukushima, S. Yuasa, and Y. Suzuki, “Brownian motion of skyrmion bubbles and its control by voltage applications,” Appl. Phys. Lett. 114, 012402 (2019). 96Y. Jibiki, M. Goto, E. Tamura, J. Cho, S. Miki, R. Ishikawa, H. Nomura, T. Srivastava, W. Lim, S. Auffret, C. Baraduc, H. Bea, and Y. Suzuki, “Skyrmion Brownian circuit implemented in continuous ferromagnetic thin film,” Appl. Phys. Lett. 117, 082402 (2020). 97J. Lee, F. Peper, S. Cotofana, M. Naruse, M. Ohtsu, T. Kawazoe, Y. Takahashi, T. Shimokawa, L. Kish, and T. Kubota, “Brownian circuits: designs,” Int. J. Unconv. 12, 341–362 (2016). 98M. A. Brems, M. Kläui, and P. Virnau, “Information processing apparatus,” European patent disclosure EP21164676.5 (2021). 99K. Raab, M. A. Brems, G. Beneke, T. Dohi, J. Rothörl, F. Kammerbauer, J. H. Mentink, and M. Kläui, “Brownian reservoir computing realized using geometrically confined skyrmion dynamics,” Nat. Commun. 13, 6982 (2022). 139 B. Bibliography 100D. Prychynenko, M. Sitte, K. Litzius, B. Krüger, G. Bourianoff, M. Kläui, J. Sinova, and K. Everschor-Sitte, “Magnetic skyrmion as a nonlinear resistive element: a potential building block for reservoir computing,” Phys. Rev. Appl. 9, 014034 (2018). 101G. Bourianoff, D. Pinna, M. Sitte, and K. Everschor-Sitte, “Potential implemen- tation of reservoir computing models based on magnetic skyrmions,” AIP Adv. 8, 055602 (2018). 102W. Jiang, L. Chen, K. Zhou, L. Li, Q. Fu, Y. Du, and R. H. Liu, “Physical reservoir computing using magnetic skyrmion memristor and spin torque nano- oscillator,” Appl. Phys. Lett. 115, 192403 (2019). 103D. Pinna, G. Bourianoff, and K. Everschor-Sitte, “Reservoir computing with random skyrmion textures,” Phys. Rev. Appl. 14, 054020 (2020). 104R. Msiska, J. Love, J. Mulkers, J. Leliaert, and K. Everschor-Sitte, “Audio classification with skyrmion reservoirs,” Adv. Intell. Syst., 2200388 (2022). 105T. Yokouchi, S. Sugimoto, B. Rana, S. Seki, N. Ogawa, Y. Shiomi, S. Kasai, and Y. Otani, “Pattern recognition with neuromorphic computing using magnetic field-induced dynamics of skyrmions,” Sci. Adv. 8, eabq5652 (2022). 106O. Lee, T. Wei, K. D. Stenning, J. C. Gartside, D. Prestwood, S. Seki, A. Aqeel, K. Karube, N. Kanazawa, Y. Taguchi, C. Back, Y. Tokura, W. R. Branford, and H. Kurebayashi, “Task-adaptive physical reservoir computing,” 10.48550/ arXiv.2209.06962 (2023). 107M. A. Brems, K. Raab, P. Virnau, and M. Kläui, “Brownscher Reservoir- Computer mit Skyrmionen,” Phys. Unserer Zeit 54, 60–61 (2023). 108Y. Ge, J. Rothörl, M. A. Brems, N. Kerber, R. Gruber, T. Dohi, M. Kläui, and P. Virnau, “Constructing coarse-grained skyrmion potentials from experimental data with iterative Boltzmann inversion,” Commun. Phys. 6, 30 (2023). 109D. Schick, M. Weißenhofer, L. Rózsa, J. Rothörl, P. Virnau, and U. Nowak, “Two levels of topology in skyrmion lattice dynamics,” in press, 2023. 110C. Song, N. Kerber, J. Rothörl, Y. Ge, K. Raab, B. Seng, M. A. Brems, F. Dittrich, R. M. Reeve, J. Wang, Q. Liu, P. Virnau, and M. Kläui, “Commensurability between element symmetry and the number of skyrmions governing skyrmion diffusion in confined geometries,” Adv. Funct. Mater. 31, 2010739 (2021). 111T. B. Winkler, J. Rothörl, M. A. Brems, H. Fangohr, and M. Kläui, “Coarse- graining collective skyrmion dynamics in confined geometries,” 10.48550/arXiv. 2303.16472 (2023). 112J. Rothörl, S. Wettermann, P. Virnau, and A. Bhattacharya, “Knot formation of dsDNA pushed inside a nanochannel,” Sci. Rep. 12, 5342 (2022). 140 B. Bibliography 113Y. Zhao, J. Rothörl, P. Besenius, P. Virnau, and K. C. Daoulas, “Can polymer helicity affect topological chirality of polymer knots?” ACS Macro Lett. 12, 234–240 (2023). 114W. F. Brown, Micromagnetics (Krieger Publishing Company, Melbourne, FL, Dec. 1978). 115A. Aharoni, “Magnetostatic energy calculations,” IEEE Trans. Magn. 27, 3539– 3547 (1991). 116J. M. D. Coey, Magnetism and magnetic materials (Cambridge University Press, Cambridge, England, Mar. 2010). 117D. J. Griffiths and D. F. Schroeter, Introduction to quantum mechanics, 3rd ed. (Cambridge University Press, Cambridge, England, Aug. 2018). 118R. Gross and A. Marx, Festkörperphysik, 2nd ed., de Gruyter Studium (De Gruyter Oldenbourg, München, Germany, Oct. 2014). 119H. Kramers, “L’interaction entre les atomes magnétogènes dans un cristal paramagnétique,” Physica 1, 182–192 (1934). 120P. W. Anderson, “Antiferromagnetism. Theory of superexchange interaction,” Phys. Rev. 79, 350–356 (1950). 121A. Manchon, Lecture series: introduction to spintronics, https://physiqueman- chon.wixsite.com/research/spintronics-lectures, 2021. 122T. H. R. Skyrme and B. F. J. Schonland, “A non-linear field theory,” Proc. Math. Phys. Eng. Sci. 260, 127–138 (1961). 123T. Skyrme, “A unified field theory of mesons and baryons,” Nucl. Phys. 31, 556–569 (1962). 124S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, “Skyrmion lattice in a chiral magnet,” Science 323, 915–919 (2009). 125A. O. Leonov and M. Mostovoy, “Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet,” Nat. Commun. 6, 8275 (2015). 126A. O. Leonov and M. Mostovoy, “Edge states and skyrmion dynamics in nano- stripes of frustrated magnets,” Nat. Commun. 8, 14394 (2017). 127T. Gilbert, “A phenomenological theory of damping in ferromagnetic materials,” IEEE Trans. Magn. 40, 3443–3449 (2004). 128J. D. Jackson, Classical electrodynamics, 3rd ed. (John Wiley & Sons, Nashville, TN, May 2021). 141 B. Bibliography 129L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii, Statistical physics: part 2, Third Edition, Course of theoretical physics (Butterworth-Heinemann, Oxford, England, Aug. 1980). 130L. D. Landau and E. M. Lifshitz, “On the theory of the dispersion of magnetic permeability in ferromagnetic bodies,” Phys. Z. Sowjetunion 8, 153–164 (1935). 131A. Aharoni, Introduction to the theory of ferromagnetism, 2nd ed., International Series of Monographs on Physics (Clarendon Press, Oxford, England, Nov. 2002). 132M. Dyakonov and V. Perel, “Current-induced spin orientation of electrons in semiconductors,” Phys. Lett. A 35, 459–460 (1971). 133Y. K. Kato, R. C. Myers, A. C. Gossard, and D. D. Awschalom, “Observation of the spin Hall effect in semiconductors,” Science 306, 1910–1913 (2004). 134J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, “Spin Hall effects,” Rev. Mod. Phys. 87, 1213–1260 (2015). 135J. E. Hirsch, “Spin Hall effect,” Phys. Rev. Lett. 83, 1834–1837 (1999). 136J. Slonczewski, “Current-driven excitation of magnetic multilayers,” Journal of Magnetism and Magnetic Materials 159, L1–L7 (1996). 137L. Berger, “Emission of spin waves by a magnetic multilayer traversed by a current,” Phys. Rev. B 54, 9353–9358 (1996). 138O. Boulle, G. Malinowski, and M. Kläui, “Current-induced domain wall motion in nanoscale ferromagnetic elements,” Mater. Sci. Eng. R Rep. 72, 159–187 (2011). 139J. Iwasaki, W. Koshibae, and N. Nagaosa, “Colossal spin transfer torque effect on skyrmion along the edge,” Nano Lett. 14, 4432–4437 (2014). 140K. Garello, I. M. Miron, C. O. Avci, F. Freimuth, Y. Mokrousov, S. Blügel, S. Auffret, O. Boulle, G. Gaudin, and P. Gambardella, “Symmetry and magnitude of spin–orbit torques in ferromagnetic heterostructures,” Nat. Nanotechnol. 8, 587–593 (2013). 141S. Emori, U. Bauer, S.-M. Ahn, E. Martinez, and G. S. D. Beach, “Current-driven dynamics of chiral ferromagnetic domain walls,” Nat. Mater. 12, 611–616 (2013). 142Q. Shao, P. Li, L. Liu, H. Yang, S. Fukami, A. Razavi, H. Wu, K. Wang, F. Freimuth, Y. Mokrousov, M. D. Stiles, S. Emori, A. Hoffmann, J. Åkerman, K. Roy, J.-P. Wang, S.-H. Yang, K. Garello, and W. Zhang, “Roadmap of spin–orbit torques,” IEEE Trans. Magn. 57, 1–39 (2021). 143S. Woo, K. M. Song, H.-S. Han, M.-S. Jung, M.-Y. Im, K.-S. Lee, K. S. Song, P. Fischer, J.-I. Hong, J. W. Choi, B.-C. Min, H. C. Koo, and J. Chang, “Spin-orbit torque-driven skyrmion dynamics revealed by time-resolved X-ray microscopy,” Nat. Commun. 8, 15573 (2017). 142 B. Bibliography 144P. Langevin, “Sur la theorie du mouvement brownien,” Compt. Rendus 146, 530–533 (1908). 145A. Einstein, “Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen,” Ann. Phys. 322, 549–560 (1905). 146D. L. Ermak and J. A. McCammon, “Brownian dynamics with hydrodynamic interactions,” J. Chem. Phys. 69, 1352–1360 (1978). 147T. Schlick,Molecular modeling and simulation: an interdisciplinary guide, 2nd ed., Interdisciplinary Applied Mathematics (Springer New York, New York, NY, Aug. 2010). 148J. A. Anderson, J. Glaser, and S. C. Glotzer, “HOOMD-blue: a python pack- age for high-performance molecular dynamics and hard particle monte carlo simulations,” Comput. Mater. Sci. 173, 109363 (2020). 149U. Nowak, “Classical spin models,” in Handbook of magnetism and advanced magnetic materials (Wiley-Blackwell, Hoboken, NJ, July 2007). 150A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comput. Phys. Commun. 271, 108171 (2022). 151H. B. Callen and T. A. Welton, “Irreversibility and generalized noise,” Phys. Rev. 83, 34–40 (1951). 152B. Dünweg and W. Paul, “Brownian dynamics simulations without Gaussian random numbers,” Int. J. Mod. Phys. C 02, 817–827 (1991). 153D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, 2nd ed., Computational science series (Academic Press, San Diego, CA, Oct. 2001). 154C. Schütte, J. Iwasaki, A. Rosch, and N. Nagaosa, “Inertia, diffusion, and dynamics of a driven skyrmion,” Phys. Rev. B 90, 174434 (2014). 155M. von Smoluchowski, “Zur kinetischen Theorie der Brownschen Molekularbe- wegung und der Suspensionen,” Ann. Phys. 326, 756–780 (1906). 156S. J. Blundell and K. M. Blundell, Concepts in thermal physics, 2nd ed. (Oxford University Press, London, England, Sept. 2009). 157Y.-W. Li and M. P. Ciamarra, “Phase behavior of Lennard-Jones particles in two dimensions,” Phys. Rev. E 102, 062101 (2020). 143 B. Bibliography 158M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P. Bernard, and W. Krauth, “Hard-disk equation of state: first-order liquid-hexatic transition in two dimensions with three simulation methods,” Phys. Rev. E 87, 042134 (2013). 159K. Zahn and G. Maret, “Dynamic criteria for melting in two dimensions,” Phys. Rev. Lett. 85, 3656–3659 (2000). 160P. Keim, G. Maret, and H. H. von Grünberg, “Frank’s constant in the hexatic phase,” Phys. Rev. E 75, 031402 (2007). 161D. R. Nelson, “Reentrant melting in solid films with quenched random impuri- ties,” Phys. Rev. B 27, 2902–2914 (1983). 162S. Deutschländer, T. Horn, H. Löwen, G. Maret, and P. Keim, “Two-dimensional melting under quenched disorder,” Phys. Rev. Lett. 111, 098301 (2013). 163E. N. Tsiok, D. E. Dudalov, Y. D. Fomin, and V. N. Ryzhov, “Random pin- ning changes the melting scenario of a two-dimensional core-softened potential system,” Phys. Rev. E 92, 032110 (2015). 164H. Wagner, “Long-wavelength excitations and the Goldstone theorem in many- particle systems with ”broken symmetries”,” Z. Phys 195, 273–299 (1966). 165N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models,” Phys. Rev. Lett. 17, 1133–1136 (1966). 166N. D. Mermin, “Crystalline order in two dimensions,” Phys. Rev. 176, 250–254 (1968). 167P. C. Hohenberg, “Existence of long-range order in one and two dimensions,” Phys. Rev. 158, 383–386 (1967). 168W.-K. Qi, S.-M. Qin, X.-Y. Zhao, and Y. Chen, “Coexistence of hexatic and isotropic phases in two-dimensional Yukawa systems,” J. Phys. Condens. Matter 20, 245102 (2008). 169J. Hansen and I. McDonald, Theory of simple liquids: with applications to soft matter, Fourth Edition (Academic Press, Oxford, England, Feb. 2013). 170V. Ramasubramani, B. D. Dice, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, “freud: a software suite for high throughput analysis of particle simulation data,” Comput. Phys. Commun. 254, 107275 (2020). 171B. D. Dice, V. Ramasubramani, E. S. Harper, M. P. Spellings, J. A. Anderson, and S. C. Glotzer, “Analyzing particle systems for machine learning and data visualization with freud,” in Proceedings of the 18th python in science conference (2019), pp. 27–33. 144 B. Bibliography 172P. Dillmann, G. Maret, and P. Keim, “Comparison of 2D melting criteria in a colloidal system,” J. Phys. Condens. Matter 24, 464118 (2012). 173K. J. Strandburg, J. A. Zollweg, and G. V. Chester, “Bond-angular order in two- dimensional Lennard-Jones and hard-disk systems,” Phys. Rev. B 30, 2755–2759 (1984). 174K. Bagchi, H. C. Andersen, and W. Swope, “Computer simulation study of the melting transition in two dimensions,” Phys. Rev. Lett. 76, 255–258 (1996). 175W. Qi, A. P. Gantapara, and M. Dijkstra, “Two-stage melting induced by dislocations and grain boundaries in monolayers of hard spheres,” Soft Matter 10, 5449–5457 (2014). 176J. Kerr, “On rotation of the plane of polarization by reflection from the pole of a magnet,” Lond. Edinb. Dublin Philos. Mag. J. Sci. 3, 321–343 (1877). 177D. B. Allan, T. Caswell, N. C. Keim, C. M. van der Wel, and R. W. Verweij, Soft-matter/trackpy: v0.6.1, https://doi.org/10.5281/zenodo.7670439, Feb. 2023. 178J. C. Crocker and D. G. Grier, “Methods of digital video microscopy for colloidal studies,” J. Colloid Interf. Sci. 179, 298–310 (1996). 179C.-H. Sow, K. Harada, A. Tonomura, G. Crabtree, and D. G. Grier, “Measure- ment of the vortex pair interaction potential in a type-II superconductor,” Phys. Rev. Lett. 80, 2693–2696 (1998). 180A. Soper, “Empirical potential Monte Carlo simulation of fluid structure,” Chem. Phys. 202, 295–306 (1996). 181Z. Li, X. Bian, X. Yang, and G. E. Karniadakis, “A comparative study of coarse- graining methods for polymeric fluids: Mori-Zwanzig vs. iterative Boltzmann inversion vs. stochastic parametric optimization,” J. Chem. Phys. 145, 044102 (2016). 182R. Henderson, “A uniqueness theorem for fluid pair correlation functions,” Phys. Lett. A 49, 197–198 (1974). 183A. O. Leonov, T. L. Monchesky, J. C. Loudon, and A. N. Bogdanov, “Three- dimensional chiral skyrmions with attractive interparticle interactions,” J. Phys. Condens. Matter 28, 35LT01 (2016). 184A. O. Leonov, J. C. Loudon, and A. N. Bogdanov, “Spintronics via non- axisymmetric chiral skyrmions,” Appl. Phys. Lett. 109, 172404 (2016). 185L. Rózsa, A. Deák, E. Simon, R. Yanes, L. Udvardi, L. Szunyogh, and U. Nowak, “Skyrmions with attractive interactions in an ultrathin magnetic film,” Phys. Rev. Lett. 117, 157205 (2016). 145 B. Bibliography 186J. C. Loudon, A. O. Leonov, A. N. Bogdanov, M. C. Hatnean, and G. Balakr- ishnan, “Direct observation of attractive skyrmions and skyrmion clusters in the cubic helimagnet Cu2OSeO3,” Phys. Rev. B 97, 134403 (2018). 187H. Du, X. Zhao, F. N. Rybakov, A. B. Borisov, S. Wang, J. Tang, C. Jin, C. Wang, W. Wei, N. S. Kiselev, Y. Zhang, R. Che, S. Blügel, and M. Tian, “Interaction of individual skyrmions in a nanostructured cubic chiral magnet,” Phys. Rev. Lett. 120, 197203 (2018). 188S. Zhang, F. Kronast, G. van der Laan, and T. Hesjedal, “Real-space observation of skyrmionium in a ferromagnet-magnetic topological insulator heterostructure,” Nano Lett. 18, 1057–1063 (2018). 189N. Kerber, M. Weißenhofer, K. Raab, K. Litzius, J. Zázvorka, U. Nowak, and M. Kläui, “Anisotropic skyrmion diffusion controlled by magnetic-field-induced symmetry breaking,” Phys. Rev. Applied 15, 044029 (2021). 190D. Chandler, Introduction to modern statistical mechanics (Oxford University Press, New York, NY, Dec. 1987). 191COMSOL multiphysics, version 5.5, https://www.comsol.com. 192M. Mochizuki, “Spin-wave modes and their intense excitation effects in skyrmion crystals,” Phys. Rev. Lett. 108, 017601 (2012). 193I. Makhfudz, B. Krüger, and O. Tchernyshyov, “Inertia and chiral edge modes of a skyrmion magnetic bubble,” Phys. Rev. Lett. 109, 217201 (2012). 194Y. Onose, Y. Okamura, S. Seki, S. Ishiwata, and Y. Tokura, “Observation of magnetic excitations of skyrmion crystal in a helimagnetic insulator Cu2OSeO3,” Phys. Rev. Lett. 109, 037603 (2012). 195C. Reichhardt and C. J. O. Reichhardt, “Dynamics of magnus-dominated particle clusters, collisions, pinning, and ratchets,” Phys. Rev. E 101, 062602 (2020). 196I. L. Fernandes, J. Chico, and S. Lounis, “Impurity-dependent gyrotropic motion, deflection and pinning of current-driven ultrasmall skyrmions in PdFe/Ir(111) surface,” J. Phys. Condens. Matter. 32, 425802 (2020). 197J. Rothörl, M. A. Brems, T. J. Stevens, and P. Virnau, “Reconstructing diploid 3D chromatin structures from single cell Hi-C data with a polymer-based approach,” in press, 2023. 146