Repositorio Institucional de Documentos

Tesis Doctoral Complex Systems: Nonlinearity and Structural Complexity in spatially extended and discrete systems (Sistemas Complejos: Nolinealidad y Complejidad Estructural en sistemas espacialmente extendidos y discretos) Memoria presentada por Jes´ us G´ omez Garde˜ nes en la Facultad de Ciencias de la Universidad de Zaragoza para optar al grado de Doctor. LUIS MARIO FLOR´ IA PERALTA, Profesor Titular del Departamento de F´ ısica de la Materia Condensada de la Universidad de Zaragoza CERTIFICA que la presente memoria, “Complex Systems: Nonlinearity and Structural Complexity in spatially extended and discrete systems”, ha sido realizada en el Departamento de F´ ısica de la Materia Condensada de la Universidad de Zaragoza bajo su direcci´ on, y autoriza su presentaci´ on para que sea calificada como Tesis Doctoral. Zaragoza, 17 de Octubre de 2006 Fdo: Luis Mario Flor´ ıa Peralta YAMIR MORENO VEGA, Investigador contratado “Ram´ on y Cajal” del Instituto de Biocomputaci´ on y F´ ısica de los Sistemas Complejos de la Universidad de Zaragoza CERTIFICA que la presente memoria, “Complex Systems: Nonlinearity and Structural Complexity in spatially extended and discrete systems”, ha sido realizada en el Departamento de F´ ısica de la Materia Condensada de la Universidad de Zaragoza bajo su direcci´ on, y autoriza su presentaci´ on para que sea calificada como Tesis Doctoral. Zaragoza, 17 de Octubre de 2006 Fdo: Yamir Moreno Vega Contents Agradecimientos v Resumen vii Introducci´ on . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Parte I: Localizaci´ on en redes no lineales de Schr¨ odinger . . . . . . . . . . viii Parte II: Estructura y propagaci´ on sobre redes complejas . . . . . . . . . . . ix Parte III: Din´ amica no lineal en redes complejas . . . . . . . . . . . . . . . . x Publicaciones relacionadas con la Tesis doctoral . . . . . . . . . . . . . . . . xii 1 Introduction: Complex Systems? 1 I Intrinsic localization in nonlinear Schr¨ odinger lattices 9 2 Discrete Breathers and Nonlinear Schr¨ odinger lattices 13 2.1 The Salerno Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Discrete space-time symmetries: (p,q) resonant states . . . . . . . . . . 19 2.3 Discrete Breathers numerics . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.1 Newton continuation . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.2 Floquet stability analysis . . . . . . . . . . . . . . . . . . . . . . 22 3 Discrete Breathers in 1D Nonlinear Schr¨ odinger lattices 25 3.1 DB’s in the standard Salerno Model . . . . . . . . . . . . . . . . . . . . 26 3.1.1 The structure of the solution . . . . . . . . . . . . . . . . . . . . 29 3.1.2 The background amplitude . . . . . . . . . . . . . . . . . . . . . 32 3.1.3 Floquet analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.2 Particle perspective on DB’s . . . . . . . . . . . . . . . . . . . . . . . . 39 3.2.1 Collective variables theory. . . . . . . . . . . . . . . . . . . . . . 39 3.2.2 Energy balance governs mobility. . . . . . . . . . . . . . . . . . 41 3.2.3 Oscillating breathers . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2.4 Validity and limitations of particle perspective . . . . . . . . . . . 48 3.3 DB’s in the Salerno Model with competing nonlinearities . . . . . . . . . 49 3.3.1 Continuum limit . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.2 Pinned discrete breathers . . . . . . . . . . . . . . . . . . . . . . 51 3.3.3 Moving discrete breathers . . . . . . . . . . . . . . . . . . . . . 57 3.4 Conclusions and Prospective Remarks . . . . . . . . . . . . . . . . . . . 61 i ii CONTENTS 4 Discrete Breathers in 2D Nonlinear Schr¨ odinger lattices 65 4.1 The Salerno model in two dimensions . . . . . . . . . . . . . . . . . . . 66 4.2 Pinned discrete breathers . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.2.1 Pinned DB’s in the standard Salerno model . . . . . . . . . . . . 70 4.2.2 Pinned DB’s in the Salerno model with competing nonlinearities . 74 4.3 Discrete vortex breathers . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.3.1 Vortex crosses . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.3.2 Vortex squares . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.3 Bound states of discrete vortex crosses . . . . . . . . . . . . . . . 86 4.4 Mobile discrete breathers . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.4.1 Structure and stability of (1,0,1) fixed points. . . . . . . . . . . . 88 4.4.2 Unstable manifold behaviour and ubiquity of pulson states. . . . . 91 4.5 Conclusions and Prospective Remarks . . . . . . . . . . . . . . . . . . . 94 II Structure and Dynamics of Complex Networks 97 5 Network Structure and Generation 101 5.1 Describing Complex Networks . . . . . . . . . . . . . . . . . . . . . . . 101 5.1.1 Basic definitions . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.1.2 Single nodes properties . . . . . . . . . . . . . . . . . . . . . . . 102 5.1.3 Network properties . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.1.4 Looking at Networks Mesoscale . . . . . . . . . . . . . . . . . . 108 5.2 Overview of network generation models . . . . . . . . . . . . . . . . . . 110 5.2.1 Random graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.2.2 Small-world networks . . . . . . . . . . . . . . . . . . . . . . . 112 5.2.3 Scale-Free networks . . . . . . . . . . . . . . . . . . . . . . . . 113 5.3 Global versus local knowledge . . . . . . . . . . . . . . . . . . . . . . . 116 5.3.1 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.3.2 Network properties . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.4 Interpolation between Random and Scale Free Networks . . . . . . . . . 121 5.4.1 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.4.2 Network growth and degree evolution . . . . . . . . . . . . . . . 123 5.4.3 Network properties . . . . . . . . . . . . . . . . . . . . . . . . . 127 6 Propagation through Complex Networks 131 6.1 Epidemic spreading and Immunization . . . . . . . . . . . . . . . . . . . 131 6.1.1 Modeling epidemic spreading . . . . . . . . . . . . . . . . . . . 132 6.1.2 Epidemic spreading in general complex networks . . . . . . . . . 135 6.1.3 Immunization strategies . . . . . . . . . . . . . . . . . . . . . . 138 6.1.4 Covering based Immunization . . . . . . . . . . . . . . . . . . . 141 6.2 Information transmission and Jamming . . . . . . . . . . . . . . . . . . . 150 6.2.1 Information dynamics on Networks . . . . . . . . . . . . . . . . 151 6.2.2 Shortest path routing . . . . . . . . . . . . . . . . . . . . . . . . 156 6.2.3 Congestion-aware routing . . . . . . . . . . . . . . . . . . . . . 160 6.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 CONTENTS iii III Nonlinear Dynamics of Complex Networks 173 7 Activatory-Inhibitory interactions in Networks 177 7.1 Modeling biological networks . . . . . . . . . . . . . . . . . . . . . . . 177 7.1.1 Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178 7.1.2 Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 7.2 Regulatory dynamics in scale-free networks . . . . . . . . . . . . . . . . 187 7.2.1 The model: basic dynamical features. . . . . . . . . . . . . . . . 187 7.2.2 Statistical characterization of island’s dynamics and structure. . . 199 7.2.3 Structure inside dynamical islands . . . . . . . . . . . . . . . . . 204 7.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209 8 Synchronization on Complex Networks 211 8.1 The Kuramoto model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212 8.1.1 Solution to the Kuramoto model . . . . . . . . . . . . . . . . . . 214 8.1.2 Synchronization in complex networks . . . . . . . . . . . . . . . 216 8.2 Synchronization in local scale-free networks . . . . . . . . . . . . . . . . 219 8.3 Homogeneous versus heterogeneous topologies . . . . . . . . . . . . . . 222 8.4 Synchronization in structured networks . . . . . . . . . . . . . . . . . . 233 8.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237 9 Conclusions 239 Appendices 243 A Computation of Periodic Orbits . . . . . . . . . . . . . . . . . . . . . . . 243 B Linear Stability of Periodic Orbits . . . . . . . . . . . . . . . . . . . . . . 246 Bibliography 251 Agradecimientos En primer lugar quiero dar las gracias a mis dos directores de Tesis, Luis Mario Flor´ ıa y Yamir Moreno, por su constante motivaci´ on, apoyo y cordialidad durante estos ´ ultimos a˜ nos. Su ayuda y consejo han sido y ser´ an siempre ´ utiles tanto en lo cient´ ıfico como en lo personal. Tambi´ en quisiera agradecer a los miembros del grupo de F´ ısica estad´ ıstica y no lineal, Fernando Falo, Juanjo Mazo, Jos´ e Luis Garc´ ıa Palacios, Pedro Mart´ ınez, Santiago Cuesta, Diego Prada y David Zueco, por haber logrado mantener un ambiente familiar dentro del trabajo. Gracias tambi´ en a aquellos, Alex Arenas, Sa´ ul Ares, Alan R. Bishop, Albert D´ ıazGuilera, Vito Latora y Franz Mertens, que me han acogido durante las estancias y visitas que se han producido a lo largo de esta Tesis. Tambi´ en quiero agradecer a los que se desplazaron hasta nuestra Universidad, tanto los anteriormente citados como Stefano Boccaletti, Michel Peyrard, ´ Angel S´ anchez y Giorgios Tsironis. Con todos ellos he disfrutado tanto discutiendo ampliamente mis trabajos de Tesis como de su amabilidad. Gracias a todos los miembros del departamento de F´ ısica de la Materia Condensada y del Instituto de Biocomputaci´ on y F´ ısica de Sistemas Complejos (BIFI), en particular, a Jos´ e Luis Alonso y Luis Mart´ ın por tener siempre las puertas abiertas. Gracias a mis compa˜ neros, y ya amigos, durante toda esta larga etapa, David, Diego, Javi, Iv´ an, Pablo y Rom´ an, con los que, adem´ as de juntarnos innumerables fines de semana, hemos discutido algo de F´ ısica de vez en cuando. No me olvido de aquellos que empezaron conmigo hace ya 9 a˜ nos, C´ esar, Diego, Guillermo y Sanvi. Tambi´ en me gustar´ ıa recordar a todos los miembros del Club Atl´ etico Chirikov, con ellos las tardes de los viernes y noches de los jueves fueron siempre apasionantes, tanto en la victoria como en la derrota. Gracias tambi´ en a mis amigos de “Marianistas”, Diego, Fernando, Javi, Laura, Luis, Mar, Oscar y To˜ no, a los que, a pesar de estar durante a˜ nos f´ ısicamente separados la mayor parte de las ocasiones, agradecer´ e siempre su imperturbable amistad y su siempre c´ alido saludo. Gracias por estar ah´ ı. Es de ley dar gracias a quien ha prestado el apoyo econ´ omico necesario para que esta Tesis fuera llevada a cabo. Gracias tanto al Consejo Superior de Investigaciones Cient´ ıficas como al Ministerio de Educaci´ on Cultura y Deporte por las becas y contratos que he disfrutado desde que comenc´ e mi labor investigadora. Por ´ ultimo, quisiera dedicar esta Tesis en primer lugar a mi familia y muy especialmente a mis padres, Jes´ us y Lola. Ellos, con todo el cari˜ no del mundo, me han ense˜ nado el camino para ser persona. Esta Tesis, como tantas otras cosas, se la debo a ellos. Y por ´ ultimo (y no menos importante), a Ana, por haber querido recorrer este camino conmigo y haberme sufrido y apoyado durante los peores momentos de esta Tesis. v xii RESUMEN •La identificaci´ on de subestructuras con din´ amica independiente correspondiente a autoregulaci´ on y “subclusters” de sincron´ ıa. •Analizar qu´ e efectos tiene la perturbaci´ on sobre distintos elementos de la red para la determizaci´ on de nodos (o grupos de nodos) relevantes. •Establecer una diferenciaci´ on entre las distintas clases de topolog´ ıa seg´ un los resultados de los tres puntos anteriores y, de esta manera, intentar determinar qu´ e par´ ametros topol´ ogicos son los causantes del diferente comportamiento. Publicaciones relacionadas con la Tesis doctoral Parte de los resultados relatados en cada parte de esta Tesis se han publicado en los siguientes art´ ıculos: •Parte I: [1] “Mobile Localization in nonlinear Schr¨odinger lattices”. J. G´ omez-Garde˜ nes, F. Falo and L.M. Flor´ ıa, Physics Letters A 332, 213 (2004). [2] “Nonintegrable Schr¨odinger discrete breathers”. J. G´ omez-Garde˜ nes, L.M. Flor´ ıa, M. Peyrard and A.R. Bishop, Chaos 14, 1130 (2004). [3] “Discrete Breathers in Two-Dimensional Anisotropic Nonlinear Schr¨odinger lattices”. J. G´ omez-Garde˜ nes, L.M. Flor´ ıa and A.R. Bishop, Physica D 216, 31 (2006). [4] “Solitons in the Salerno Model with competing nonlinearities”. J. G´ omez-Garde˜ nes, B.A. Malomed, L.M. Flor´ ıa and A.R. Bishop, Phys. Rev. E 73, 036608 (2006). [5] “Discrete solitons and vortices in the two-dimensional Salerno model with competing nonlinearities”. J. G´ omez-Garde˜ nes, B.A. Malomed, L.M. Flor´ ıa and A.R. Bishop, Phys. Rev. E 74, 036607 (2006). •Parte II: [6] “Local versus Global knowledge in the Barab´asi- ´ Albert scale-free model”. J. G´ omez-Garde˜ nes and Y. Moreno, Phys. Rev. E 69, 037103 (2004). [7] “Improved routing strategies for Internet traffic delivery”. P. Echenique, J. G´ omez-Garde˜ nes and Y. Moreno, Phys. Rev. E 70, 056105 (2004). RESUMEN xiii [8] “Distance-d covering problems in scale-free networks with degree correlations”. P. Echenique, J. G´ omez-Garde˜ nes, Y. Moreno and A. V´ azquez, Phys. Rev. E 71, 035102(R) (2005). [9] “Dynamics of jamming transitions in complex networks”. P. Echenique, J. G´ omez-Garde˜ nes and Y. Moreno, Europhysics Letters 71, 325 (2005). [10] “Immunization of Real Complex Communication Networks”. J. G´ omez-Garde˜ nes, P. Echenique and Y. Moreno, Eur. Phys. J. B 49, 259 (2006). [11] “From Scale-free to Erd¨os-R´enyi Networks”. J. G´ omez-Garde˜ nes and Y. Moreno, Phys. Rev. E 73, 056124 (2006). •Parte III: [12] “On the robustness of complex heterogeneous gene expression networks”. J. G´ omez-Garde˜ nes, Y. Moreno and L.M. Flor´ ıa, Biophysical Chemistry 115, 225 (2005). [13] “Michaelis-Menten dynamics in complex heterogeneous networks” J. G´ omez-Garde˜ nes, Y. Moreno and L.M. Flor´ ıa, Physica A 352, 265 (2005). [14] “Scale-Free topologies and Activatory-Inhibitory interactions” J. G´ omez-Garde˜ nes, Y. Moreno and L.M. Flor´ ıa, Chaos 16, 15114 (2006). [15] “Current trends in the modeling of biological networks”. Y. Moreno, L.M. Flor´ ıa and J. G´ omez-Garde˜ nes, AIP Conference Proceedings 851 ”From Physics to Biology: The interface between experiment and computation: Bifi 2006 II International Congress” 150 (2006). [16] “Synchronization of networks with variable local properties”. J. G´ omez-Garde˜ nes and Y. Moreno, cond-mat/0608309 [to appear in Int. J. Bif. Chaos (2007)]. [17] “Paths to Synchronization on Complex Networks”. J. G´ omez-Garde˜ nes, Y. Moreno and A. Arenas, cond-mat/0608314 [submitted]. Chapter 1 Introduction: Complex Systems? The greatest challenge today, not just in cell biology and ecology but in all of science, is the accurate and complete description of complex systems. Scientist have broken down many kinds of systems. They think they know most of the elements and forces. The next task is to reassemble them, at least in mathematical models that capture the key properties of the entire ensembles. Edward O. Wilson [1]. This thesis covers the analysis of two fundamental ingredients for the correct modeling of real macroscopic systems: Nonlinearity and Structural Complexity. The study of systems where these ingredients are present is systematically related to the field of “Physics of Complex Systems”. It is not easy however to find a formal definition of what a complex system is and most books on the matter submit the reader from the very beginning to some illustrative examples of complex phenomena rather than establishing the general principles and characteristics of a complex system. The standard classification of the wide range of studied natural systems into physical fields is mainly based on the energy range involved for their description, namely, (ranging from higher to lower energy) particle physics; nuclear physics; molecular atomic and optical physics; soft and condensed matter. With this classification in mind the question about what field the Physics of Complex System belongs to naturally arises. However, complex systems are common to a number of physical disciplines belonging to different “energy ranges” so it is difficult to place them into a single physical compartment. In addition, one can find lots of examples of systems called “complex” outside these traditional branches of physics in chemistry, biology, ecology, and social and economical sciences. Then, rather than becoming a particular physical field, the physics of complex systems has emerged as an interdisciplinary subject. What are the unifying characteristics of complexity in the phenomena studied by such a diverse branch of scientific disciplines? One can state that the fingerprint of a complex system is reflected by the display of organization without a central organizing principle. This collective organizational behaviour is not usually explained by decomposing a complex system into its parts and analyzing their isolated properties. In this sense, the Physics of Complex Systems is a new way for analyzing systems where complex phenomena are displayed rather than a new physical field. 1 2 CHAPTER 1. INTRODUCTION: COMPLEX SYSTEMS? One of the first attempts to show the need of a Physics of Complex Systems is due to Philip W. Anderson in his celebrated article “More is different” [2]. In this article, Anderson tells about the concept of broken symmetry reflected when one moves from a small system to a macroscopic one. In this transition it may happen that some symmetries of the single systems, that determine their physical behaviour, are lost when embedded in a bulk of many systems and unexpected phenomena occur (Emergence). In this latter situation the knowledge of the physical laws governing each single building block is in many cases not enough to explain the collective behaviour of the big system. Emergentism versus Reductionism This new way of thinking in physics is strongly related with the Emergentism philosophy. The advent of Emergentism philosophy constituted a punch at the old-fashioned Reductionist movement that led the theory of science for decades and, in particular, the way of physics during the XIX century and the first half of the last century. Emergentism states that the observed phenomena are classified into different levels of description and that each one of these levels is independent in the sense that each has its own laws. The emergence of such levels is the result of the increase of the problem’s complexity. On the other hand, reductionists assume the unity of science so that a hierarchical organization is established among disciplines: •Chemistry is based on Physics. •Fundamental Biology is based on Chemistry. •Psychology is based on Biology. •Sociology is based on Psychology. •Political science and Anthropology are both based on Sociology. Whereas the first two of these reductions were commonly accepted, it was not the case for the last ones yielding a big controversy. For example, aspects of evolutionary psychology and socio-biology are rejected by those who claim that complex systems are inherently irreducible. On the contrary, strong reductionists claim that the behavioural sciences should become truly scientific disciplines by being based on genetic biology arguments. Examples of this long-standing controversy can be found in several scientific forums. Perhaps, the most fruitful ones, in the sense of the number of rationales on the subject, are found in the context of life sciences: the mind-body problem, the explanation of Darwinian evolution theory, etc... As the theory of science evolved, the Reductionism-Emergentism controversy progressively affected each stage of the above chain of reductions, arriving up to the very first level. In fact, it was Karl Popper, one of the godfathers of the modern theory of science, who stated that Chemistry was not reducible to Physics [3]. The big controversy arrives to Physics under the term “Complex Systems” questioning the possibility of explaining every physical system only in terms of the properties of its constituents (Constructionist hypothesis). As mentioned above, a great deal of physical systems is said to display complex behaviour, i.e. they display some new phenomena 3 that cannot be predicted by only looking at their parts. This impossibility goes beyond the limitations related to the large amount of elements that are involved; there are cases where the properties of the isolated elements seem to be violated by the physical description of the macroscopic behaviour. The concept of symmetry breaking is thus central for a proper (in physical terms) definition of what is called emergent phenomena in Physics. The aforementioned article by P. W. Anderson [2] introduced this new concept in “one of the early manifestos for this infinitely quiet revolution” leading the new way of looking at physical phenomena. Let us review the most salient examples of the so-called complex behaviour. Symmetry breaking arguments have been recursively invoked to explain unexpected experimental discoveries in the field of condensed matter physics. The most representative examples are superconductivity, superfluidity, liquid crystals and antiferromagnetism. The physical explanation of each of these collective phenomena has contributed to the growth of conceptually new frameworks in many-body physics. It was Landau [4, 5] the first to formulate phase transitions as processes where symmetry reductions occurred. This point of view allowed him, after the theoretical explanation of ferromagnetism, to predict antiferromagnetism [6, 7] by generalizing the idea of spin rotation symmetry breaking. Subsequently, ideas on gauge symmetry breaking led to the explanation of superfluidity [8, 9] and superconductivity [10]. Anderson [2] thus claimed that symmetry breaking in many-body systems is a generalized phenomenon yielding different emergent collective behaviours depending on the particular type of broken symmetry. Therefore, it seems reasonable that, since the principles that govern the behaviour of a macroscopic system of particles appear defined by the whole system, these phenomena should be studied separately to those found for more elementary levels of description. The success of the holistic description of simple (in terms of the knowledge of the governing physical laws) systems where collective behaviour shows up, paved the way for the search of new conceptual schemes to explain other emergent phenomena at higher complexity levels. This search was carried out by the study of mathematical models that account for different emergent behaviours such as the appearance of dynamically coherent states, spontaneous localization or pattern formation in extended systems. A second set of motivating phenomena is revealed by the observation of self-similar (fractal) patterns in natural systems, that is viewed as the emergence of self-organization behaviour. These examples are different forms of the dynamical or spatio-temporal order that appears in seemingly different systems of a large number of interacting elements in nature. The breakthroughs in the description of the first class of phenomena is closely related to the advances in the studies of nonlinear dynamical systems. From the very first soliton theory [11] in continuous nonlinear systems and the discovery of localized states in nonlinear chains [12], we have seen that coherent structures emerge from large scale nonlinear models, possess their own entity (particle-like behaviour, well-defined life times, characteristic interaction patterns, etc...). In fact, no matter the complexity of the underlying equations, spatial or temporal coherent structures are many times described with the help of a low dimensional phase space. The concept of self-organized criticality, introduced by Per Bak, Kurt Wiessenfeld and Chao Tang [13, 14], constitutes one of the best explanation of nature complexity and, perhaps, it represents one of the major conceptual achievements of the physics of complexity. Self-organized criticality tries to capture the essential ingredients to explain the 4 CHAPTER 1. INTRODUCTION: COMPLEX SYSTEMS? critical-like behaviour (manifested by observations of fractal structure and power-laws) of many natural systems without a central controller unit. The original idea was to describe the dynamics of sandpiles, accounting for the avalanches that may happen when grains are progressively incorporated, by means of a simple cellular automaton model. The success of this simple model was seized to relate the model to a variety of phenomena where criticality was already observed, like earthquakes, forest fires, epidemics and, indeed, evolution theory (relating it to the the theory of punctuated-equilibrium [15]). Reading again Anderson in [16] (written nearly 20 years after his treatise on complexity) and reviewing the most relevant examples of complex behaviour up to now, it is evident that emergence is progressively being accepted as a necessary ingredient to face all the new phenomena that has appeared in the last decades under the name of Complexity. A number of scientific institutes and groups are nowadays contributing to the growth of complexity physics which is still in its infancy. Centers with a long history tackling complexity such as the “Sante Fe Institute”, the “Max Plank Institute for Physics of Complex Systems” at Dresden, the “Complex Systems group and the Center for Nonlinear Studies” at Los Alamos, the “New England Complex Systems Institute” at Boston, etc... and those of much recent creation like our “Institute of Biocomputation and Physics of Complex Systems” at Zaragoza are actively contributing to this growth. In summary, the physics of complex systems tries to explain emergent phenomena without losing the sight of the whole system (unlike the fully reductionist way of doing, which destroys the systemic level). It is then important to keep in mind that, although physical systems have a clear hierarchical ordering (nobody doubts that a system is composed by its parts) “each level can require a whole new conceptual structure” [2] (at least for our limited way of thinking), making thus impossible bridging the gaps by the systematic use of a bottom-up approaches. What are the essential ingredients of a Complex System? All the examples listed above are labeled as complex phenomena due to the impossibility to explain them by studying in isolation the parts of the systems where they occur. The behaviour of these systems is thus intrinsically new with respect to the properties of the single elements of the system. It is clear that a complex system is composed of a large number of elements. However, a big ensemble of building blocks is not enough by itself to guarantee the emergence of unexpected phenomena such as those described above (long range correlations from short-range interactions, localization in extended systems, self-organization and adaptability, etc...). Then, it is important to distinguish the attribute complex from complicated. Airplanes, computers or swiss clocks are examples of complicated systems made up of a large number of pieces. They all conform a directed cause-effect chain so that the malfunctioning of a single piece stops the whole system. They are also designed by an external agent (humans) for a unique function. What are then the key ingredients of a system for the observation of complex phenomena? Before being tempted to answer this question, it comes the doubt about whether it is reasonable to expect complex systems to have a well defined number of characteristic properties. Unavoidably, our mind tends naturally to overuse classification, which is the 5 natural (hardwired) way of thinking we manage for everyday’s life. However, the holistic roots of emergent phenomena makes it intrinsically difficult to have well defined boundaries between what is a real complex system and what is not. (In fact, we have up to now defined a complex system by their behavioural rather than structural properties.) On the other hand, the number of complex phenomena observed numerically and experimentally gives us some hints for unveiling some recurrent structural ingredients. Let us briefly summarize the most relevant ones: •Nonlinearity: It is clear that only a few natural systems can be described by means of linear relations. The need for a nonlinear modeling of the interactions is clearly seen by the nonlinear response of many real systems to perturbations. In fact, nonlinearity also appears when first principles equations are obtained for very simple systems. •Non regular structure of the interactions: The network of interactions is of utmost importance. It is revealed a high diversity in the amount of connections that a single element of the system has with combination of short-range and long-range links. In addition, loops are very usual in real systems. •Surroundings do matter: Many complex systems are open. They may share a balance between dissipated energy and incoming energy flux with the surroundings in order to achieve dynamical stability. Although neither complete nor precise (they can be presented alone or all together), the properties listed above are shared by a number of systems where complex phenomena is observed. Evidently, natural complex systems, like e.g. a protein, incorporate all these complexity levels, but, on the other hand, synthetic models incorporating only one of them are able to reproduce the relevant phenomena. It is seen that one of these ingredients in synthetic models, specially in the case of nonlinear systems with many degrees of freedom, can lead the system towards complexity. It is also worth stressing that, as well as the difficulties for explaining complex phenomena, all the above properties incorporate additional mathematical and computational difficulties. Despite the efforts for unveiling the attributes of a complex system, the question about what are the essential ingredients that generate complexity remains open, specially when it is clear that a variety of dynamical mechanisms can produce self-similar structures. Our “Complex Systems” The association of several ingredients of complexity in synthetic models is even more interesting than the study of systems where only one source of complexity exists. One would expect the observation of new emergent phenomena, different from those related to any of the sources of complexity. This expectation motivates our studies in this Thesis. In particular, we will focus on extended systems of interacting elements where two sources of complexity are present, namely, nonlinearity and/or structural complexity. As explained above, the emergence of coherent structures in extended nonlinear systems has been studied since decades. Our concern in this Thesis is to dedicate a first part to the study of localization in nonlinear homogeneous lattices. In particular, we will 6 CHAPTER 1. INTRODUCTION: COMPLEX SYSTEMS? address the study of localized states in one and two-dimensional nonlinear Sch¨ odinger lattices. These states, usually termed intrinsic localized modes or discrete breathers are time periodic, spatially localized and are seen as ubiquitous solutions to a number of homogeneous nonlinear lattices (for it the attribute “intrinsic” in their denomination). Besides, nonlinear Sch¨ odinger lattices are seen as paradigmatic equations of importance for several branches of physics like Bose-Einstein condensates or nonlinear optics. In this first part we will be specially interested in the mobility of such coherent structures. The main difference with classical solitons in continuum equations relies in the absence of continuum space translational symmetry, that makes the finding of such solutions non trivial. Besides, we will study other types of coherent states, like discrete vortices or oscillating discrete breathers, in order to have a complete description of the behaviour of localized solutions in this important class of lattices. The second part of this Thesis will concern the study of complex networks,i.e. extended systems of interacting elements where the patterns of connections between them is random. The study of this class of systems has been traditionally ascribed to graph theory. However, the recent discoveries on the self-similar character of the structure of connections in many real (social, biological, technological, etc...) systems have led to a burst in the activity of the so-called physics of complex networks. There is a lot of important consequences of the scale-free behaviour of real networks, like robustness under random perturbations, the “small-world” effect, absence of threshold for epidemics spreading and a complete new behaviour for most dynamical processes that take place on top of them. The self-similar patterns of connections and its ubiquity in nature lead to the conclusion that a large amount of systems share the same self-organizational principles. However, the question about the mechanism that drives the evolution of networks to these common structural patterns is still unsolved. We will focus on both the modeling of network growth and the study of several simple dynamics of interest in human-made (technological) scale-free networks. The study of complex network structure and the analysis of simple dynamics on top of scale free graphs try to unveil the improvements that a heterogeneous pattern of connections provides to the deployment of the network’s function. However, these two elements, function (dynamics) and structure, are many times presented to us entangled. That is, the growth and time evolution of the network of interactions (that determines its scalefree feature) is performed at the same time the system develops its function. Therefore, the structure is the result of a kind of selective process that drives to the most efficient architecture. In this case, the study of how network grows is similar to the problem of finding those network architectures that are the most efficient for its functioning. Besides, most of the dynamics of real systems are seen to be nonlinear and, therefore, the analysis of systems of elements with both nonlinear and random interactions comes as necessary. Our main purpose in the third part of the Thesis is to analyze two systems of this kind and, therefore, approach to the problem on the Structure-Function relation. We will analyze this relation in two biologically relevant systems, namely, a scale-free network with activatory-inhibitory (Michaelis-Menten type) interactions and the Kuramoto model of phase oscillators on top of different network architectures. In these two studies we do not pretend to find the definite answer to the Structure-Function problem, but to discuss new tools and discover new phenomena that could lead to a better understanding of this relation. 7 The studies on this Thesis are thus separated in three parts depending on the sources of complexity involved in their description: nonlinearity (Part I), structural complexity (Part II) and (finally Part III) both. Along this thesis we will face problems associated to the emergence of coherent structures, self-similar structural patterns and finally selforganization of dynamical patterns. Therefore, the concept of emergence will be the recurrent idea behind our studies. 14 CHAPTER 2. DISCRETE BREATHERS AND NONLINEAR SCHR ¨ ODINGER LATTICES observation of these localized states was found to be generic of homogeneous nonlinear lattices and therefore different from that due to presence of Anderson modes as a consequence of the existence of any inhomogeneity (defect or impurity) of the harmonic lattice. Observations of genuine discrete breather solutions were mainly based on the numerical simulations of the nonlinear dynamics. The question on their existence as true solutions of the nonlinear lattice remained unsolved until R.S. Mackay and S. Aubry [22] established the theorem for the existence of discrete breather solutions. This theorem is based on the concept of anti-integrability (developed by S. Aubry for studying the Frenkel-Kontorova model [49, 50]) or, applied to general lattices, the anticontinuum limit. This concept refers to the limiting case when there is no coupling between adjacent sites of the lattice so that the system is composed of a set of independent oscillators whose dynamics is governed by their corresponding on-site potentials. Then, considering the state where a single oscillator evolves following an orbit of frequency ωbwhile the remaining sites are in the rest state one can ask whether this state of energy confinement would remain when the coupling between sites is adiabatically incorporated. The continuability of discrete breathers from the uncoupled limit implies two conditions •Non-resonance condition: The oscillation frequency and its harmonics must rely outside the phonon band of the lattice at the rest state nωb,ω(q)∀q∈[−π/2, π/2] (n=1,2, ...) (2.2) •Anharmonicity condition: The on-site potentials governing the dynamics of the isolated sites must be nonlinear so that the frequencies, ωb, of their orbits fulfills ∂ωb/∂I,0, where Iis the action. The proof of the existence theorem is based on the implicit function theorem and provides a practical way for constructing localized solutions of the type (2.1). After the rigorous formulation of the existence conditions of discrete breathers several questions arised. From one hand, its stability and robustness in noisy environments has been studied in detail [51] since their experimental observation and potential applications to real systems implies relative large life times. Another hot topic is the issue of their mobility. Taking into account the general form (2.1) of a pinned discrete breather one would expect their mobile counterparts to have the form Φn(t)=f(n−vb−x0)exp(iωbt),(2.3) with f(n−vb−x0)∼exp[Γ|n−vbt−x0|] when n→ ±∞. The possibility of transferring energy packets across lattices opens the door to a wide range of applications in nonlinear optics, solid state and soft matter physics. However, since the continuous translational invariance is broken due to discreteness, the computation of pinned discrete breathers of the form (2.1) does not guarantee the success in constructing mobile localized states like (2.3) by means of a change of the reference system. Different approaches have been used for studying this problem ranging from the “kicking” method [52–55] (where a static solution is perturbed with the so-called pinning or marginal mode in order to make it move) to analytical approximations were continuous variables (collective coordinates) accounting for the localization center are introduced [56–60]. Our approach to this problem tries 2.1. THE SALERNO MODEL 15 to generalize the method of continuation for pinned breathers to obtain mobile solutions. For this purpose we start with make use of the concept of (p,q) resonant states that will allow us to unify the problem of finding both mobile and pinned discrete breathers. In this chapter we introduce the set of nonlinear Schr¨ odinger equations that we study along the two forthcoming chapters as well as to summarize the basic definitions and techniques used for characterizing breather solutions to these equations. We will start in section 2.1 describing the Salerno model [61] which provides a two-parametric family of nonlinear Schr¨ odinger lattices. In section 2.2 we address the definition of the concept of (p,q) resonant solutions to which general discrete breather solutions belong. Finally, in section2.3 the basic technique to obtain and characterize discrete breathers are summarized. 2.1 The Salerno Model The continuous nonlinear Schr¨odidinger equation (NLS) constitutes a key tool for a number of fields as diverse as the study of Bose-Einstein condensates (where the mean field approximation is of the NLS-type, the Gross-Pittaevskii equation), the study of nonlinear (Kerr type) optical fibers, molecular chains (where Davydov solitons are studied), etc... Besides, the NLS equation is specially interesting for nonlinear physics since it appears when considering the lowest order of nonlinearity for any dynamical equation on a dispersive medium where energy is conserved. The most general form of this equation is i˙ Φ(x,t)=−▽2Φ(x,t)−γ|Φ(x,t)|2Φ(x,t),(2.4) where Φ(x,t) is a complex field which, in the context of Bose-Einstein condensates, accounts for the macroscopic wave-function of the condensate. The parameter γaccounts of the competence between the dispersive (Laplacian term) and the nonlinear parts. This cubic nonlinear equation posses the singular property of being integrable. The integrability was probed by means of the Inverse Scattering Method (ISM) technique [11, 62, 63] in [64] providing a family of nonlinear waves. The Discrete Nonlinear Schr¨ odinger equation The physical relevance of the NLS equation along with its integrable character make it one of the most studied models by nonlinear physicists during the last decades. Besides, discretizations of this equation are also of great interest. The natural the discretization of eq. (2.4) yields the so-called standard discrete nonlinear Schr¨odinger equation (DNLS) [23], i˙ Φn=−C(Φn+1+ Φn−1)−γ|Φn|2Φn,(2.5) where Φnis now a complex variable, the parameter Camounts the nearest neighbor coupling, and γis the strength of the nonlinearity. The above discretization does not conserve the integrability of the continuous model (2.4) although the wide applicability to physical fields is preserved. In particular, the DNLS equation is particularly relevant for •Dynamical description of Bose-Einstein condensates trapped in a periodic potential well (optical trap) [33–35, 65–68]. 16 CHAPTER 2. DISCRETE BREATHERS AND NONLINEAR SCHR ¨ ODINGER LATTICES •Pulse dynamics in nonlinear waveguide arrays [28, 29, 69–72]. •Adiabatic approximation of the Holstein polaron [23]. •Excitation dynamics in biopolymers lattices [26]. The dynamics governed by the DNLS equation (2.5) is derived from the Hamiltonian H=−CX n (ΦnΦn+1+ ΦnΦn+1)−γ 2X n|Φn|4,(2.6) where Φndenotes the complex conjugate of Φn. Both variables, Φnand Φn, are canonically conjugated with the usual Poisson structure {U,V}=X n"U Φn V Φn−V Φn U Φn#.(2.7) Beside, the DNLS equation has a second integral of the dynamics, namely the norm, N=X n|Φn|2,(2.8) that in the context of Bose-Einstein condensates accounts for the total number of bosons whereas for waveguides arrays it is the total power of the beam. The Ablowitz-Ladik equation AnotherimportantdiscretizationofthecontinuousNLSequationistheso-calledAblowitzLadik equation (AL). Although this lattice is not so physically relevant as the usual discretization, DNLS equation (2.5), the AL model preserves the integrability of its continuous counterpart (2.4). In fact, the AL model is an extremely exceptional example of an integrable nonlinear lattice that was discovered by M.J. Ablowitz and J.F. Ladik in 1976 [73, 74] by means of the ISM in its discrete version [75, 76]. The AL model reads as follows i˙ Φn=−(Φn+1+ Φn−1)C+γ 2|Φn|2,(2.9) where again Φn(t) is a complex probability amplitude, the parameter Camounts the nearest neighbor coupling, and γis the strength of the nonlinearity. The nonlinear term in the AL equation is of the intersite type and hence differs with its counterpart in the DNLS model which is an onsite nonlinearity. The AL model (2.9) has a deformed Poisson structure defined by {U,V}=X n"U Φn V Φn−V Φn U Φn#1+γ 2|Φn|2,(2.10) and the conserved Hamiltonian is H=−CX n (ΦnΦn+1+ ΦnΦn+1).(2.11) 2.1. THE SALERNO MODEL 17 The integrability of the AL equation results in an infinite number of conserved quantities. Along with the Hamiltonian the two conserved magnitudes of lowest order in {Φn}are N=2 γX n ln(1 +γ 2|Φn|2),(2.12) P=iX n (ΦnΦn+1−ΦnΦn+1),(2.13) which are the norm and the momentum respectively. The integrable AL equation possesses a two-parameter family of exact breather solutions of the form Φn(t)=s2 γsinhβsech[β(n−x0(t))] × exp [i(α(n−x0(t)) + Ω(t))] .(2.14) As can be observed the solutions possess the continuous spatial symmetry x0→x0+ǫ and hence analytic mobile breather solutions with a similar form to eq. (2.3) are available for this exceptional lattice. The two parameters of this breather family can be chosen to be the breather frequency ωband velocity vb, vb=˙x0=2sinhβsinα β(2.15) ωb=˙ Ω = 2coshβcosα+αvb,(2.16) where −π≤α≤πand 0 < β < ∞. The AL moving breather (instantaneous) profile interpolates between the rest state Φn=0 of the lattice (at n→ ±∞) in an exponentially localized region around x0(t), while traveling with velocity vb. The Salerno model In the above two equations, DNLS (2.5) and AL (2.9), the self-focussing effect of local nonlinearity balanced by the opposite effect of the dispersive coupling makes possible the existence of localized periodic solutions (breathers) of the discrete field, where the profile of |Φn|decays exponentially away from the localization center: Φn(t)=|Φn|exp[iωb(t))] .(2.17) In the uncoupled limitC→0 of the DNLS equation, also known as the anti-integrable or anti-continuous limit, discrete breathers can be easily constructed by selecting a periodic oscillation Φn0(t) of frequency ωb=γ|Φn0|2at site n0and Φn=0 for n,n0. These solutions can be uniquely continued (we will see the procedure below) to nonzero values of the coupling C, and constitute the one-parameter family of immobile on-site breathers of the DNLS equation. Unfortunately, the continuation from the uncoupled limit does not provide solutions where the localization center moves along the lattice with velocity vb(as for the AL case), i.e, mobile discrete breathers. On the other hand, the connection between the integrable 18 CHAPTER 2. DISCRETE BREATHERS AND NONLINEAR SCHR ¨ ODINGER LATTICES (though physically limited) AL equation and the physically relevant (though nonintegrable) DNLS equation is provided by the model originally introduced by M. Salerno in [77], i˙ Φn=−(Φn+1+ Φn−1)hC+µ|Φn|2i−2νΦn|Φn|2.(2.18) The above lattice provides a Hamiltonian interpolation between the standard DNLS equation (2.5), for µ=0 and ν=γ/2, and the integrable AL lattice (2.9) when µ=γ/2 and ν=0. In the following we will set the value of γ=2. The Hamiltonian of the Salerno equation is given by H=−CX n (ΦnΦn+1+ ΦnΦn+1)−2ν µX n|Φn|2 +2ν µ2X n ln(1 +µ|Φn|2),(2.19) which contains the AL and DNLS Hamiltonian for the above limits. The Poisson structure of the Salerno model takes the form {U,V}=X n"U Φn V Φn−V Φn U Φn#1+µ|Φn|2,(2.20) which, for µ,0, takes the same functional form as that of the Ablowitz-Ladik equation, eq. (2.10), and in the limit µ=0 it becomes the standard Poisson structure according to that of the DNLS limit, eq. (2.7). In addition to the Hamiltonian, this equation possesses, for any value of the parameters µand ν, the following conserved norm N=1 µX n ln(1 +µ|Φn|2).(2.21) While the SM was originally introduced in a rather abstract context, it has recently found direct physical realization, as an asymptotic form of the Gross-Pitaevskii equation describing a Bose-Einstein condensate of bosonic atoms with magnetic momentum trapped in a deep optical lattice [78]. In that case, the onsite nonlinearity is generated, as usual, by collisions between atoms, while the intersite nonlinear terms account for the long-range dipole-dipole interactions. This latter interaction may be attractive (µ > 0) or repulsive (µ < 0), if the external magnetic field polarizes the atomic momentum along the lattice or perpendicular to it, respectively. The continuation of the family (both pinned and mobile) discrete breathers from the AL integrable limit allows numerical observations of the interplay between the integrable term, weighted by the parameter µ, and the nonintegrability, weighted by ν. We will inspect the effects that the combination of these two nonlinearities, with both similar and opposite signs, has on discrete breathers dynamics. 2.2. DISCRETE SPACE-TIME SYMMETRIES: (P,Q) RESONANT STATES 19 2.2 Discrete space-timesymmetries: (p,q)resonant states In order to unify the problem of finding pinned and mobile discrete breathers by means of continuation methods we start defining the concept of (p,q) resonant states. Suppose that a frequency ωb=2π/Tbis given, we will say that a solution Φ = {Φn(t)}is (p,q)resonant with respect to the reference frequency ωb, if the following condition holds, for all nand t: Φn(t)= Φn+p(t+qTb).(2.22) After q Tb-periods, these solutions repeat the same profile but displaced by plattice sites. In more technical terms, these (p,q) resonant solutions are fixed points Φof the operator LpTq=M(2.23) (M−I)Φ = 0,(2.24) where Land Tare, respectively, the lattice translation and the Tb-time evolution operator L{Φn(t)}={Φn+1(t)}(2.25) T{Φn(t)}={Φn(t+Tb)}.(2.26) We now consider some examples of (p,q) resonant solutions with respect to the frequency ωb; the first example is simply provided by the family of plane wave solutions of eq. (2.18): Φn(t)=Aexp[i(kn −ωt)] .(2.27) It is easily seen, by inserting (2.27) in eq. (2.18), that the values of ω,kand |A|define a surface in the three-dimensional space, the nonlinear dispersion relation surface ω(k,A) (see figure 2.1): ω=−2[1 +µ|A|2]cosk−2ν|A|2.(2.28) Note that due to the nonlinear character of the eq. (2.18), the frequency ωdepends on both wave number kand amplitude |A|of the plane wave. 0 0.5 1 1.5 2 |A|2-3 -2 -1 0 1 2 3 k -4 -2 0 2 ωFigure 2.1: Plot of the nonlinear dispersion relation surface of nonlinear plane waves, eq. (2.28), as a function of the amplitude Aand the wave number kof the plane wave. The values of µand νare fixed to 0.5. 20 CHAPTER 2. DISCRETE BREATHERS AND NONLINEAR SCHR ¨ ODINGER LATTICES One can easily determine those plane waves that are (p,q) resonant with respect to ωb: the eq. (2.22) imposes the following condition on ωand k ω ωb =1 qp 2πk−m,(2.29) where mis any arbitrary integer. These planes in the 3-d space (ω,|A|,k) intersect the dispersion relation surface at (in general) several one-parameter families (branches) kj(|A|), in the first Brillouin zone (−π≤k≤π). If we are not interested in unreasonably large (and not interesting) amplitude values |A|of the plane waves, the number of branches is finite: one can see that for fixed values of all the parameters (p,q,ωb,ν,µ), there is a finite number of branches in the limit |A| → 0; there is also a well defined (parameter dependent) threshold value of the amplitude at which a pair of new branches (tangent bifurcation) appear (i.e. these plane waves can only resonate with ωbfor amplitudes above some threshold value). Thus, by a suitable bounding of the amplitude, for each couple (p,q) one finds a finite number, s, of branches of (p,q) resonant plane waves. (Note also that this number diverges when p/qtends to an irrational). A different, and highly nontrivial, example of (p,q) resonant solutions is provided by the solitary waves (2.14) of the AL lattice. From eq. (2.16) it is clear that the choice 2πvb/ωb=p/qselects a (p,q) resonant solitary wave with respect to the frequency ωb, i.e. a breather solution where the two time scales involved, given by its frequency ωband velocity vb, are commensurate. The set of velocity values of resonant AL breathers is dense and any AL moving breather is a limit of some sequence of resonant ones. Note also that immobile breathers are (0, 1) resonant with respect to the frequency ωb. In the integrable limit, the plane waves and the AL breathers are both exact independent solutions. Integrability makes possible that the initial localization of energy is maintained with time evolution, without decaying away by exciting radiation. It is a well established result that (even far away from this integrable limit) immobile discrete breathers remain exact solutions of the lattice dynamics. Our concern in the next sections is the question of moving discrete breathers away from integrability in eq. (2.18). In order to study them, we will focus on (p,q) resonant solutions. The motivating of this restriction comes from its accessibility to numerics. First we will motivate the numerical (Newton) method that allow us to study these solutions with an adequately high precision. 2.3 Discrete Breathers numerics We introduce here the numerical techniques that we have used. As a whole, one could refer to them as the (SVD-) regularized Newton method. They do not constitute a novel method in ”discrete Breather numerics”, as they have been already used, e.g. in [79] to refine moving breathers of Klein-Gordon lattices obtained by other numerical means (see, by contrast, [80]). From the methodological side, what is novel here is the systematic use of them in the investigation of the family of moving Schr¨ odinger breathers reported below in 3.1. To some extent, the presentation here is self-contained but for further details on these techniques we refer to the Appendices and the proposed bibliography. First, in 2.2 we 2.3. DISCRETE BREATHERS NUMERICS 21 introduce the notion of (p,q) resonant solution, providing some illustrative examples. The (SVD) regularized Newton algorithm is presented in 2.3.1, and finally in 2.3.2 we briefly explain the basics of Floquet stability analysis. 2.3.1 Newton continuation A well-known numerical procedure to obtain exact periodic solutions of nonlinear lattices is the Newton continuation [22, 53, 79, 81]. The different practical implementations of this procedure work very successfully when, for example, one obtains numerically exact immobile discrete breathers of eq. (2.18), from the uncoupled limit µ=0 and C=0, where exact periodic discrete breathers are trivially constructed. The iteration of the Newton operator Tconverges rapidly to its fixed point (i.e. the solution to be computed) provided the starting point, ˆ Φ0, is close enough, and the solution of the following system of linear equations is a well-posed problem: (DT −1)(Φn−Φn+1)=[T −I]Φn,(2.30) where DTis the Jacobian matrix of the Newton operator, and Φn(the n-th iteration solution of (2.30)) converges quadratically to the fixed point solution. By adiabatic change of a model parameter, one constructs a uniquely continued exact fixed point solution for each parameter value, using each time, as starting point of the Newton iteration, the solution previously computed. The matrix (DT − 1) must be invertible, in order to uniquely compute Φn+1. Degeneracies associated with the +1 eigenvalues of DT, if any, have to be removed in order to obtain a unique fixed point solution. When continuing immobile (time periodic) discrete breathers of eq. (2.18), a convenient prescription is commonly used, namely to restrict the operator action to the subspace of time-reversible solutions (see Appendix A and [53, 81]). This provides a practical way of removing degeneracies, allowing unique continuation of immobile discrete breathers. However, for the continuation of general (p,q) resonant solutions (of which periodic solutions are only the particular case p=0 and q=1), one has to use M=LpTqas the Newton operator. One has also to deal with the degeneracies of M, and imposing time-reversibility could, in this case, be too restrictive, since in general (p,q) resonant solutions are not time-reversible. A well-known solution to the problem of removing degeneracies when no clear restrictions are available, is provided by the so-called singular value decomposition (SVD) [53, 79, 82, 83] of the matrix (DLpTq−1) : (DLpTq−1) =J=PVQ ,(2.31) where P,Vand Qare 2N×2Nsquare matrices. Pand Qare orthogonal matrices and Vis diagonal (vjδij) with possibly null (zero) elements, called singular values, associated with the null space of J(the subspace that is mapped to zero Jx =0). The columns of Pwhose same-numbered elements vjare nonzero are an orthonormal set of basis vectors that span the range of J(the subspace reached by this matrix). The rows of Qwhose samenumbered elements vjare zero are an orthonormal basis for the null space of J. One can 22 CHAPTER 2. DISCRETE BREATHERS AND NONLINEAR SCHR ¨ ODINGER LATTICES numerically use this SVD decomposition, checking the (numerical) vectors spanning the null space to identify degeneracies, and using at iteration steps the pseudoinverse matrix Q∗ˆ V−1P∗,(2.32) where ˆ V−1is diagonal with elements 1/vjfor vj,0 and 0 for vj=0. The convergence criterion for the fixed point solution is that X j[T −I]Φn+1j<N·10−16 ,(2.33) where N is the size of the lattice, i.e. the solutions obtained along the two forthcoming chapters can be regarded as exact up to machine precision. As a judicious test of our numerical codes, we have used both procedures (reduction to time-reversible subspace and SVD decomposition) to obtain immobile discrete breathers (for which both methods are valid) of the Salerno model. Both agree, up to the highest possible accuracy, from the uncoupled limit up to the A-L limit (and viceversa). 2.3.2 Floquet stability analysis A very useful outcome of the numerical Newton method of computing solutions of eq. (2.18) is the Jacobian matrix of the Newton operator, usually called the Floquet matrix F. This matrix is the linear operator associated with the linear stability problem (see Appendix B and [84]) of the fixed point solution. Indeed, the Jacobian Fof the Newton operator M F=DM(2.34) maps vectors in the tangent space of the solution (small initial perturbations ~ ǫ(0) of the fixed point solution) into their TM-evolved vectors, i.e. ~ ǫ(TM), after a period of M. That is: ~ ǫ(TM)=F~ ǫ(0) ,(2.35) The Floquet matrix of a Hamiltonian system is real and symplectic, so the Floquet eigenvalues λcome in quadruplets, λ, 1/λ, ¯ λ, 1/¯ λ. The necessary condition for the stability of the solution is that all the eigenvalues lie on the unit circle of the complex plane, |λ|=1. To illustrate the Floquet analysis of (p,q) resonant solutions of the NLS lattice (2.18), we now obtain the Floquet spectrum of modulational instabilities of a (p,q) resonant plane wave, Φn(t)=Aexpi(kn −ωt).(2.36) One has to investigate the evolution of small perturbations, in both amplitude and phase, of the plane wave Φn(t)=(A+In)expi(kn −ωt+ϕn),(2.37) where we assume that the perturbation parameters are small compared with those of the plane wave solution. Introducing expression (2.37) in (2.18) and considering the following form for the perturbations {In, ϕn}: In(t)=Iexpi(Qn −Ωt) ϕn(t)=ϕexpi(Qn −Ωt) (2.38) 2.3. DISCRETE BREATHERS NUMERICS 23 1 1.02 1.04 1.06 1.08 1.1 -0.5 0 0.5 1 1.5 |λ|exp(-TM) θFloq / TM A2=0.01 A2=0.02 A2=0.03 A2=0.04 A2=0.05 Figure 2.2: Plot of the modulus of the unstable Floquet eigenvalues |λ|(corresponding to the positive values of ℑ(Ω) in eqs. (2.44) and (2.45)), versus the Floquet angle, θFloq. Both quantities are conveniently normalized to the period of the map TM. The amplitude of the excursion of |λ| and the range of values of θFloq for which |λ|>1 grow as the amplitude Aof the plane wave is increased. The parameters in eq. (2.18) are µ=ν=0.5 and the wave number of the plane wave is k=0.5. we obtain the dispersion relation for the perturbation parameter Ω: [Ω−2(1 +µA2)sinksin Q]2=16(1 +µA2)× sin2Q/2cosk[(1 +µA2)sin2Q/2cosk −µA2cosk−νA2],(2.39) as obtainedin [85, 86]. Fromthe above expression onederivesthe values ofΩ(A,Q,k;ν, µ) for the modulational perturbations. When the parameter Ωhas a nonzero imaginary part, i.e. the right-hand side of (2.39) is negative, the plane wave (A,k) becomes unstable under the corresponding modulational (Q) perturbation, whose amplitude will grow exponentially fast in the linear regime (tangent space). Modulational perturbations (2.38) correspond to eigenvectors {In, ϕn}of the Floquet matrix: In(t+TM)=exp(−iΩTM)In(t) (2.40) ϕn(t+TM)=exp(−iΩTM)ϕn(t) (2.41) with associated Floquet eigenvalues exp(−iΩTM). The real part of Ωgives the angle in the complex plane, θFloq =−ℜ(Ω)TM,(2.42) while the imaginary part ℑ(Ω) gives the modulus of the Floquet eigenvalue, |λ|=exp(ℑ(Ω)TM),(2.43) 30 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES -0.75 -0.5 -0.25 0 0.1 0 20 40 60 80 100 ℜ(Φn) -0.75 -0.5 -0.25 0 0.1 0 20 40 60 80 100 ℑ(Φn) -150 -300 -450 0 20 40 60 80 100 site ϕn 10-4 10-3 10-2 10-1 100 0 20 40 60 80 100 |Φn|2 Site Site (d) (a) (b) (c) Figure 3.3: Instantaneous profile of a (1,1) resonant breather with ωb=2.678 and vb=0.426; the nonintegrable parameter is ν=1.0 (standard DNLS equation). (a) Real part, (b) imaginary part, (c) modulus and (d) phase. The resonant condition for the harmonic composition of the background gives the contribution of three plane waves. The existence of these plane waves is revealed by the modulation of the extended tail in the modulus profile (c). |Φn|2 |A|2 120 160 200 time 60 80 100 site 10-3 10-2 10-1 100 101 Figure 3.4: Time evolution of |Φn|2profile of a mobile discrete Schr¨ odinger breather. The frequency of the solution is ωb=5.050 and the velocity is vb=0.804. Note that the background is composed by a single plane wave with amplitude A. The nonintegrable parameter of eq. (2.18) is ν=0.2. 3.1. DB’S IN THE STANDARD SALERNO MODEL 31 10-3 10-1 101 103 105 0 0.5 1 1.5 2 2.5 3 3.5 4 S(ω) ω 0 1 2 3,4 S( ) ω (b) -1 -0.5 0 0.5 1 1.5 -3 -2 -1 0 1 2 3 ωph/ωb kph 4 0 5 2 3 1 6 (a) Figure 3.5: (a) Plot of the graphical solving of the resonant condition (in the Aj→0 limit) for a (1,2) resonant breather with ωb=2.384 and vb=0.189. (b) Power Spectrum S(ω) of the background of this solution at ν=1.0. From (a) eq. (2.29) gives the contribution of seven plane waves (j=0, ..., 6) but only five (j=0, ..., 4) of them are visible due to the difference of orders of magnitude between the amplitudes |Aj|. The agreement between the resonant condition equation (for the fitted value of Aj) and the frequencies observed in S(ω) is up to machine accuracy. localized core needs for its motion to ”surf over” a specific extended state of radiation (see figure 3.4): (ˆ Φbackg)n(t)= s−1 X j=0 Ajexpi(kn −ωjt).(3.6) We note that among the members of the (s-parameter) continuous family of (p,q) resonant plane waves (see Section I), the fixed point solution contains only a particular member (Aj,ωj) from each branch (see figure 3.5.a). This selection varies smoothly with the (adiabatic) continuation parameter ν. In particular, the amplitude modulus |Aj|selected increases smoothly from its zero value at the integrable limit (ν=0), for both signs of ν. If the bare core of a fixed point solution (i.e. after subtraction of the background) is taken as initial condition for a direct numerical integration of the equations of motion, one observes radiative losses, along with the corresponding changes in shape, velocity, etc. of the localized moving core. The motion of the bare localized core (not anymore a solution) excites extended states of the lattice. Thus, regarding the exact fixed point solution, one could say that radiative losses of the running core are exactly canceled out when the localized core runs, with specific velocity, on top of the specific linear combination of (Aj, ωj) resonant plane waves (3.6). A complementary numerical observation is the following: Taking as initial condition for a direct integration of the equations of motion (2.18), a superposition of an immobile discrete breather and the background of a (p,q) resonant mobile breather, it evolves into a moving discrete breather, with approximate velocity vb=(pωb)/(2πq). One thus would say that the background promotes breather translational motion with adequate velocity. In the next section 3.1.2, a connection between background characteristics and the particle perspective (i.e. the Peierls-Nabarro barrier of collective variable theories), will be established in order to further illuminate the physical description of discrete breather mobility. 32 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES Whatever physical perspective one may prefer, the numerical fact is that the generic structure of the fixed point solution is given by the superposition (3.3). Not too far from ν≃0, where the amplitudes Ajof the fixed point background have small values, one can carefully check that if the bare core is given as a starting guess for Newton iteration, this converges well to the exact complete solution (core +background), by developing the specific selection of Ajamplitudes. This confirms the robustness of the numerics. Though previous observations of nondecaying tails of numerically accurate mobile discrete breathers in Klein-Gordon lattices [53] and/or (solitary) traveling waves [95] in self-focusing equations had been reported (see also the interesting discussions on this issue in [80] and [96]), no systematic study on those tails and their role is available. However we clearly see that they are an essential part of the exact solution. As argued in the introductory section, the translational motion of a discrete breather introduces a new time scale. In a nonintegrable context, this fact unavoidably implies resonances with plane wave band spectra, and an exact self-sustained moving DB solution could only exist on top of a developed resonant background. This seems to have been (with a few exemptions) not fully appreciated in most of current literature on mobile breathers, where the background is most often either ignored or deliberately suppressed. A notable feature of the plane wave content of the background ˆ Φbackg is that the amplitude modulus |Aj|in (3.6) differ by orders of magnitude, i.e. |A1| ≫ |A2| ≫ |A3|..., so that only a few frequencies are dominant for most practical purposes (see figure 3.5.b). In other words, the extended background associated to a spatially localized moving core is, in turn, strongly localized in the reciprocal (k-space) lattice. The possible relevance of this observation is further discussed below in the concluding section. 3.1.2 The background amplitude Inordertocharacterizethespecificfeaturesofthenonintegrablemotionofdiscretebreathers, we focus here on the (perhaps) most remarkable among those features: the background amplitude of the uniquely continued fixed point. How does it evolve along the continuation path in parameter space? For positive values of νwe have followed the line in parameter space (figure 3.1) µ+ ν=1 (see equation (2.18)), while for negative values, we took the path µ−ν=1. Note that taking this latter path is similar to studying staggered breathers in the former one due to the staggering transformation reported above. We do not expect other paths to make important differences. As stated earlier, near ν≃0, the amplitude grows from its zero value (at the integrable limit) for both signs of this parameter, for it is a nonintegrable effect. However, for larger values of nonintegrability |ν|the background amplitude evolution shows some important differences for the two signs of ν. In figure 3.6 we plot the background amplitude (modulus) of the (1, 1) resonant fixed point, versus the continuation parameter ν, for three different values of the breather frequency ωb. For ν > 0, one observes that the amplitude steadily increases with νbefore continuation stops (i.e. Newton iteration ceases to converge beyond a certain maximum ν value). Note that the amplitude grows faster for higher values of the frequency, and that the continuation stops (correspondingly) at a smaller value of ν. This may suggest that the failure of fixed point continuation is related to a somewhat excessive growth of the background amplitude, an issue that will be discussed later. 3.1. DB’S IN THE STANDARD SALERNO MODEL 33 0 0.005 0.01 0.015 0.02 0.025 -0.4 -0.2 0 0.2 0.4 |Φbackg|2 ν abc Figure 3.6: Background amplitude versus νfor three different (1/1) resonant breathers with frequencies: (a) ωb=5.65, (b) ωb=4.91, (c) ωb=4.34. Note the two different behaviours: for positive values of ν|Φbackg|2is a monotonous increasing function of νwhile for the negative part it shows smooth rises and falls. For ν < 0, after an initial growth the background amplitude decreases down to almost negligible values around ν≃ −0.3, then grows and again decreases close to zero at ν≃ −0.39, and so on, in progressively narrower intervals with larger peak amplitude, until continuation stops. Most noticeable is the fact that the intervals neither depend on the breather frequency ωbnor on the breather velocity vb. Why do background amplitudes decay so dramatically at those regions in parameter space? An important hint is presented in the next section, where the Floquet stability analysis of immobile discrete breathers will show a coincident situation of mirror-symmetry breaking (and its absence for positive ν values). For other values of pand qthat we have numerically investigated, the same features of the background amplitude variation as shown in figure 3.6 are qualitatively reproduced. 3.1.3 Floquet analysis On the basis of the general arguments given in [84, 97], the Floquet spectra of immobile DB in the thermodynamic limit, N→ ∞, consists of two components: the (continuous) Floquet spectrum of the asymptotic state of the solution (rest state), and a discrete part associated with spatially localized eigenvectors. The continuous part is composed by small amplitude (linear) plane waves, the so-called phonons. However, for mobile DB the asymptotic state of a (p,q) resonant fixed point solution is a superposition of plane waves, the background ˆ Φbackg. From this, one should expect the Floquet spectrum of a (p,q) resonant DB being composed of two components: the discrete (spatially localized eigenvectors) and a continuous part associated with the linear stability of the background plane waves. The continuous part of the Floquet spectrum should reflect the same results of the modulational instability analysis of section 2.3.2. In particular, this means that any 34 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES 10 0 0.2 0.4 0.6 0.8 1 ωb ν 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 00000000000000000000000 11111111111111111111111 11111111111111111111111 11111111111111111111111 11111111111111111111111 11111111111111111111111 11111111111111111111111 Figure 3.7: Continuation diagram of (1,1) resonant breathers as a function of the frequency ωb. The end of the numerical continuation, νmax(ωb), is represented by the line with dots. The region where mobile breathers suffer from core instabilities is limited by the shaded area. modulational instability a plane wave may suffer will be also an instability of a fixed point solution whose background contains this plane wave. In the future we will refer to any instability of the continuous part of the Floquet spectrum as background instability. Any instability from the discrete part is a core instability. First we focus on core instabilities. For this we turn attention to the continuation of mobile (p,q) resonant breathers. Figure 3.7 shows in the ν−ωbplane (dotted line), the values νmax(ωb) where the numerical continuations stop due to non convergence of Newton iteration for p=1, q=1 and ν > 0. As it was remarked above, the continuation stop is associated with the rapid increase of the background amplitude shown in figure 3.5. Only low frequency breathers, for which the background amplitude increases more slowly, can be numerically continued all the way to the standard DNLS equation. The linear stability analysis of (p,q) resonant breathers yields a well defined region in the ν−ωbdiagram where core instabilities appear. There is an island inside the continuation region of figure 3.7, where the Floquet spectra contain a real eigenvalue λ > 1. We observe the evolution of this Floquet eigenvalue (and its complex conjugate) as the parameter νis increased in figure 3.8.a, for a (1, 1) breather of frequency ωb=2.678. Here the angle (θFloq) in the complex plane is plotted versus ν. The interval of constant zero angle corresponds to the section (constant ωb) of the instability island in figure 3.7. Along the whole continuation path, the profile of the corresponding unstable eigenvector is localized. An example of this profile inside the instability island is shown in figures 3.8.b and 3.8.c, where one observes that the localized instability shows a decaying background along the direction opposite to the motion. The decay rate increases as the modulus of the eigenvalue grows and decreases again when λreturns to the unit circle. On the other hand, the stable Floquet eigenvector associated with 1/λ shows a wing decaying along the mirror symmetric direction. The direct integration of the equation of motion reveals that the unstable solution experiences a pinning after a transient of regular motion 3.1. DB’S IN THE STANDARD SALERNO MODEL 35 (a) -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0 0.1 0.2 0.3 0.4 0.5 θFloq ν (b) -0.5 -0.25 0 0.25 0 20 40 60 80 100 120 ℜ(εunst) (c) 0 0.25 0.5 0 20 40 60 80 100 120 site ℑ(εunst) Figure 3.8: (a) Floquet angle evolution of the spectra of a (1,1) resonant breather with ωb= 2.678. The thick trajectory corresponds to the localized eigenvector that becomes unstable (θFloq = 0 interval). Instantaneous profile of the real (b) and imaginary (c) part of the Floquet unstable eigenvector of a (1,1) breather with ωb=3.207 and ν=0.26. The decaying tails along the direction opposite to the motion reveals the energy loss that the unstable eigenvector causes to the solution. 36 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES 0 -0.3 -0.6 -0.9 -1.2 -1.5 -1 -0.5 0 0.5 1 k0 ν Stable Unstable Figure 3.9: Modulational instability existence diagram for a plane wave with wave number k0∈[−π/2,0]. This diagram fixes the region where mobile discrete breathers with a background composed of only one plane wave do not suffer from background instability. with velocity vb=p/(qTb). After the solution pins at site n, its core center oscillates around this site. The trapping of the unstable MB could be interpreted as a result of the energy losses that the growth of the linearly unstable perturbation induces on the solution. Returning to the instability island shown in the diagram of figure 3.7, some final observations are worth summarizing: (i) there is a range of frequencies where mobile breathers of the standard DNLS equation (ν=1) suffer from this instability; (ii) very high frequency breathers do not experience this instability (in the short range where they can be continued); (iii) very low frequency breathers are stable all the way up to ν=1. We turn now to background instabilities. Once we know the plane wave content (k0, k1,..) of a (p/q)-resonant fixed point, we can know whether the solution is subject to MI or not and, if it is unstable, what are the harmful perturbations (Q). This problem is not so simple because we cannot know a priori the plane wave content if we do not have the amplitudes of each one (2.29). However, we can derive a necessary condition for not having MI if we consider that, from (2.29), the background is always composed of at least one plane wave (m=0) with k0between [−π/2,0]. From this we can simplify the analysis of the background stability to the k0plane wave stability as a necessary condition for the MB stability. For this we calculate, for each νand k, the value of the right-hand side of (2.39) for all the range of Q( [−π, π] ) and A. If this value is always positive the plane wave with this k0is free from modulational instabilities at this point of the model (2.18) with parameter ν. From this extensive exploration we obtain, see figure 3.9, the region in the k−νplane where MI is present. In the range of νbetween [−1,−0.5] there is no modulational instability for single plane waves of any value of kbetween [−π/2,0], and in particular for k0. However, this does not guarantee that moving breathers are free from these instabilities in this region, unless the background has only one plane wave (as is sometimes the case). On the contrary, in the region ν > 0 any moving breather suffers such instabilities. The transition area in the region ν∈[−0.5,0] presents MI depending on which k0we have. For the 3.1. DB’S IN THE STANDARD SALERNO MODEL 37 (a) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 ℑ(λ) ℜ(λ) (b) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 ℑ(λ) ℜ(λ) (c) -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 ℑ(λ) ℜ(λ) Figure 3.10: Floquet spectra of (1,1) resonant breathers: (a) for ωb=4.348, vb=0.692 and ν=0.08 the spectra shows the core (localized) instability; (b) for ωb=6.610, vb=1.052 and ν=0.07 the spectra shows the background (modulational) instability (also present but not visible in (a)); (c) for ωb=4.348, vb=0.692 and ν=−0.39 the solution is linearly stable. 38 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES -1 -0.5 0 0.5 1 -0.305 -0.303 -0.301 -0.299 ξ ν one-site DB two-site DB two-site DB 18 12 6 046 48 50 52 54 |Φn|2 18 12 6 046 48 50 52 54 |Φn|2 18 12 6 046 48 50 52 54 |Φn|2 -1 -0.5 0 0.5 1 -0.392 -0.391 -0.39 -0.389 ν one-site DB two-site DB two-site DB two−site DB DB one−site DB asymmetric (b) (a) Figure 3.11: Graphical representation of the two first symmetry breaking bifurcations for ν < 0. The quantity ξin the vertical axes of both figures is defined, referred to the one-site breather, as the difference between the modulus |Φ|of the two sites adjacent to the maximum (|Φmax|), i.e. ξ∼ |Φmax−1| − |Φmax+1|. For the one-site DB ξ=0 and for the two-site DB ξ=1, for this ξis conveniently normalized with the the difference between Φmax and Φmax±1. The continuous lines represent the regions where the static solutions are linearly stable while the discontinuous ones represent the unstable regions. The modulus profile of the three immobile coexisting solutions are plotted in the central insets for ωb=6.215 and ν=−0.3012. range where no plane-wave with kbetween [−π/2,0] is subject to MI we can assure that if there is only one contribution, k0, to the background the corresponding MB solution is stable. For example, this is the case for (1/1) resonant breathers if ωb>4 and for (1/2) resonant breathers if ωb>8.46. The Floquet spectra of a moving breather satisfying these requirements is plotted in figure 3.10.c. After the analysis of both types of instabilities eventually experienced by moving Schr¨ odinger breathers, we finally report on a most relevant numerical fact revealed by the Floquet analysis of the family of standard immobile discrete breathers for ν < 0 (or, similarly, the family of staggered immobile discrete breathers for ν > 0). Near ν≃ −0.3 an immobile two-site DB experiences a mirror symmetry-breaking (pitchfork) bifurcation becoming linearly unstable. When approaching the bifurcation point, two conjugate Floquet eigenvalues quickly approach +1, where they meet, and then separate along the real axis. The eigenvector associated to the unstable λ > 1 Floquet eigenvalue is localized and odd-symmetric, and is termed the symmetry-breaking or depinning mode φdep. We recall here that the background of an immobile breather is the rest state ˆ Φ = 0, whose continuous spectrum consists of small amplitude (linear) plane waves. The depinning mode, on the other hand, is a localized core instability of the immobile breather, favoring a translation of the core center. For a smaller value of ν≃ −0.39 there is another symmetry-breaking bifurcation where the two-site breather becomes stable, again interchanging the stable character with the one-site. The corresponding bifurcation diagram for these two symmetry breaking transitions is plotted in figure 3.11. In the first symmetry breaking bifurcation, two unstable mirror-asymmetric immobile breathers emerge from the bifurcation point, progressively evolve toward the (stable) two- 3.2. PARTICLE PERSPECTIVE ON DB’S 39 site breather, and finally collide in a new pitchfork bifurcation from where a unstable two-site breather emerges. The net result is an inversion of stability between one-site and two-site immobile breathers. Around the narrow interval of νvalues where these two bifurcations occur, the energies of the three types of breathers involved (one-site, two-site, and asymmetric) have very small differences. From a particle perspective, this should make the breather motion easier. It is precisely in this same narrow interval where (see 3.1.2) we observe that the background amplitude of moving breathers becomes negligible. This is not a coincidence as we will argue in 3.2. 3.2 Particle perspective on discrete breathers The appealing framework and success of collective variable approaches (see e.g. [56–60, 98]) to the problem of nonintegrable motion of discrete breathers relies on the fidelity of a particle-like description of these field excitations that they provide. In these approaches, the effective dynamics of only a few degrees of freedom (e.g. the localization center, and the spatial width of the state, etc...in some instances [65, 68]) replaces the whole description of the moving localized state. Though unable to account for all the nonintegrable features, perturbative collective variable theories of NLS lattices provide a sensible physical characterization of important features of the nonintegrable mobility of localized solutions, like the emergence [99] of a Peierls-Nabarro barrier to motion. Here we summarize the main results of this particlelike description and compare them with the behaviour of numerically exact (p,q) resonant moving breathers. Our goal is twofold: to acquire a correct physical understanding of the numerical facts, and then to make an assessment of validity and intrinsic limitations of collective variable approaches. 3.2.1 Collective variables theory. A presentation of the particle perspective on moving Schr¨ odinger breathers near the AL integrable limit can be found in [58, 59] (see also [56, 57, 60, 98]), where the interested reader will find the relevant formal aspects of the theory. Using the integrable solitary wave (2.14) as an ansatz for the moving breather solution in the perturbed AL lattice, ν,0 and small in (2.18), one considers the parameters α, β,x0and Ωas dynamical variables (variation of constants). The time evolution of these parameters in the perturbed lattice is governed by: ˙x0=2sinαsinhβ β(3.7) ˙ Ω = 2cosαcoshβ+α˙x0+g(β) (3.8) ˙ β=0 (3.9) ˙α=−ν∞ X s=0 8π3sinh2β β3sinh(π2s/β)sin(2πsx0) (3.10) 46 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES 3.2.3 Oscillating breathers The emergence of the Peierls barrier and the behaviour of the background amplitude illustrate the physical interpretation of this background as a (p/q)-resonant energy support to overcome the barrier to motion. We now confirm this statement searching for another kind of solution: oscillating breathers. These solutions are predicted by collective coordinates approaches and are a consequence of the loss of translational invariance out of the integrable limit. Following the above interpretation of the background role one can imagine certain solutions with a background amplitude not high enough for surpassing the Peierls barrier and allowing travel along the lattice. In terms of a well defined potential, considering the particle perspective, the center of these localized solutions would be oscillating between (n−1/2) and (n+1/2) for the unstaggered ones or between nand (n±1) for the staggered ones. From our perspective, the oscillating breathers are solutions with two frequencies: the internal one of the breather (ωb) and the one corresponding to the oscillatory motion (ωosc). Once again, we have a problem dealing with two time scales and consequently we have to impose that the two frequencies are commensurate pωb=qωosc. The fixed point problem is now associated with the map: TqTbΦn(t)= Φn(t) (3.17) We cannot, however, develop the Newton iteration scheme in a similar way as for mobile breathers. There is no longer any family of oscillating breathers providing a good start point for the continuation (they are intrinsic solutions of the nonintegrable regime because they appear as the Peierls barrier emerges). The way to obtain a good ansantz (as Cretegny and Aubry already used to find mobile breathers in Klein-Gordon lattices [53]) is to perform a small perturbation of the static solution (pinned at a site n) with the depinning internal mode: Φansantz n= Φstatic n(ωb)+ǫδφdep n(3.18) The dynamics of the perturbed solution for small enough values of ǫshows the oscillating behaviourexpectedandforlargeenough values ofǫthebreather starts to move. Obviously in both cases the motion finishes after a transient due to radiation (they are not exact solutions). Tuning the parameter ǫwe search for those oscillatory transients whose ωosc is resonant with the breather frequency ωb. The transient is much more stable when the nonintegrable parameter νis very small, close to the AL limit. We first search here for a good initial guess for the method and then obtain the exact solution of the map (3.17). Once the exact solution is obtained for a small ν, we can perform the continuation to higher values in the same way as we did for mobile solutions. In figure 3.16.a we show the evolution of the amplitude of oscillation as νis increased from 0.05 to 0.18. The amplitude of the oscillation is represented by the phase portrait of the localization center of the breather defined as in (3.16). The continuation reflects that the amplitude of the oscillation, for a fixed value of ωosc, grows with ν. In figure 3.16.b the density plot of |ˆ Φn|2is shown as a function of time, revealing the oscillating pattern of the solution. The existence of exact oscillating breathers is a consequence of the existence of a Peierls barrier. The structure of these solutions reveals the existence of a background (resonant with the map) whose amplitude grows as ν(and consequently the amplitude of 3.2. PARTICLE PERSPECTIVE ON DB’S 47 (a) -0.12 -0.08 -0.04 0 0.04 0.08 0.12 40.5 40.75 41 41.25 41.5 V0 X0 (b) Figure 3.16: (a) Evolution of the localization center x0of an exact (1/18)-oscillating breather for different values of ν: 0.05, 0.06,.., 0.18. The internal frequency is ωb=3.086. The amplitude of the oscillation of x0increases with νrevealing the nonlinear character of the motion for the highest values of ν.(b) Density plot of the time evolution of |ˆ Φn|2for the above oscillating breather. 48 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES oscillation) is increased. This is the picture we expected from the role played by the interaction background-core in the energy balance during motion. The monotonous growing behaviour of the background versus the oscillation amplitude, strongly suggests that if the amplitude of the former is increased the solution will be able to translate steadily. This has been checked by direct numerical integration, because no exact solutions connecting the oscillating with the mobile ones can be obtained due to the different maps employed to obtain both types of solutions. However, the existence of a background in the exact oscillating breather solutions and its behaviour with the amplitude of the breather oscillations are fully consistent with the interpretation of the results obtained for the mobile solutions. 3.2.4 Validity and limitations of particle perspective The most basic result of the perturbative collective variable theories away from the integrable regime is the existence of a Peierls-Nabarro potential function of the core (collective variable) center. It expresses (in particle-like terms) that the breather position is no longer indifferent because the continuous translational invariance has been broken. From this also naturally comes the existence of oscillating breathers. We have seen how our numerics fully confirm the qualitative validity of these predictions. A further prediction concerns the phase portrait’s transition studied in [57]. Despite the fact that our end of continuation is correlated with the equipotential lines profile of the numerical PN barriers, and the phase portrait transition is also related to their sudden growth, no clear connection (between transition and end of continuation) can be established. The end of continuation is itself sensibly interpreted as a numerical consequence of the sudden increases of the amplitude background, and does not imply neccesarily the existence of any global phase portrait transition. However, in some respects the perturbative collective variable theory is clearly incomplete: For example, it is unable to predict the observed localized (core) instability bifurcation and the observed symmetry breaking transitions for ν < 0. These bifurcations could easily appear in a theory with (at least) two variables (a dimer) experiencing the Peierls-Nabarro potential, which would demand an improved perturbative ansatz. This improved ansatz must coincide in the integrable limit with the AL solution. One can use the numerical results to guide the construction of such an improved ansatz. In this respect, the following observation may be relevant. The parameter βof the AL solution determines both the amplitude and the width of the localized pulse. However, our numerical estimates of these breather characteristics for immobile breathers show clearly that, for fixed value of ωb, the breather width is independent of ν, while the amplitude varies with it. In other words, away from integrability, width and amplitude of the (immobile) breather are no longer a single collective variable. Beyond any other limitation of the perturbative collective variable theory, the background (an indispensable part of the exact solution) is absent in the perturbative ansatz, and it cannot appear later in that context. A complete theory of (nonlinear Schr¨ odinger) breather motion should somehow incorporate the background in the ansatz itself. If correct, it should then predict that the background amplitude grows from zero with the nonintegrability parameter ν, and (ideally) so on with all the numerically observed behaviours. One possible way to develop the analytical approach could be to use the method presented in [100]. In this scheme, eq. (3.15) may play an important role, for it provides the energy 3.3. DB’S IN THE SALERNO MODEL WITH COMPETING NONLINEARITIES 49 balance governing the translational motion of the breather core. In other words, our results show that the core energy is not an invariant of motion and this requires the existence of a finely tuned background whose nonlinear interaction with the core compensate the core energy variations. 3.3 Discrete breathers in the Salerno model with competing nonlinearities In the above sections we have mainly focused on the study of mobile discrete breathers. In fact, the characterization of usual (non-staggered) pinned discrete breathers along the standard (µ > 0 Salerno path was already considered in previous works [58, 59, 89, 90] concluding that eq. (2.18) gives rise to pinned discrete breathers at all values of the DNLS parameter ν, and all positive values of the AL coefficient, µ. As already mentioned above if νis negative one can make it positive by means of the staggering transformation, and hence study those staggered pinned discrete breather along the SM with ν > 0 (and hence finding the symmetry breaking bifurcation reported in section 3.1.3). However, the sign of µcannot be altered. In particular, the proper AL model (ν=0) with µ < 0 does not give rise to localized solutions. The latter circumstance suggests considering soliton dynamics in the SM with µ < 0, i.e., with competing nonlinearities, which is the subject of the present section 1. In order to study the SM with competing nonlinearities, it is necessary to redefine the conserved norm (2.21) and Hamiltonian (2.19) by N=1 µX n ln1+µ|Φn|2,(3.19) H=X n"−ΦnΦ∗ n+1+ Φn+1Φ∗ n−2ν µ|Φn|2+2ν µ2ln1+µ|Φn|2#.(3.20) Whereas the Poisson structure of the standard Salerno model (eq. 2.20) remains valid. The above redefinitions of the norm (3.19) and Hamiltonian (3.20) are introduced in order to remain valid when h1+µ|Φn|2itakes negative values at some sites, due to the use of µ < 0. In this section we will study the existence and characterization of both pinned and mobile discrete breathers when these two competing (on-site self-focusing and inter-site self-defocusing) nonlinearities coexist in the Salerno model. In 3.3.1 a continuum approximation (CA) of the Salerno model is used in order to investigate the behaviour of the discrete breathers when µ < 0 in an analytical form. It is found that, although they might exist in a semi-infinite band of frequencies (as occurs for the above studied case µ > 0), they actually occupy a finite band, with an solution (peakon) at the edge of the band. After this calculations a family of discrete breathers is constructed for µ < 0 in section 3.3.2 by means of a continuation of these pinned solutions from the standard DNLS limit (µ=0) where they are easily obtained. The continuation results show that they form a family of 1Remind that the SM with µ < 0 is also physically relevant for it describes the repulsive case for the longrange dipole-dipole interactions in a Bose-Einstein condensate of bosonic atoms with magnetic momentum trapped in a deep optical lattice as introduced in section 2.1. 50 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES regular pinned discrete breathers, including a peakon-like one, similar to what was found in the CA, but discrete breathers extend beyond the peakon in the form of a novel solution termed cuspon that we will characterize in this part. In section 3.3.2, the pinned breather stability is explored by means of both standard Floquet analysis and direct simulations, with the conclusion that only a small part of the family is unstable. Two-breathers bound states are reported in 3.3.2, where it is demonstrated that stability exchange between inphase and out-of-phase states occurs at a point where the bound breathers are peakons. For what concerns to mobile breathers we show in section 3.3.3 that they can be continued up to a critical strength of the inter-site self-defocusing nonlinearity. 3.3.1 Continuum limit To introduce the continuum approximation (CA) in eq. (2.18), we define Φ(x,t)≡ e2itΨ(x,t), and expand Ψn±1≈Ψ±Ψx+(1/2)Ψxx, where Ψis now treated as a function of the continuous coordinate x, which coincides with nwhen it takes integer values. After that, the continuum counterpart of eq. (2.18) is derived, iΨt=−2(1−|µ|)|Ψ|2Ψ−1−|µ||Ψ|2Ψxx ,(3.21) where we have set ν= +1 and µ < 0, in order to inspect the interesting region. Equation (3.21 ) conserves the norm and Hamiltonian, which are straightforward counterparts of expressions (3.19) and (3.20), Ncont =1 µZ+∞ −∞ dx ln1−|µ||Ψ|2,(3.22) Hcont =Z+∞ −∞ "|Ψx|2+2 1 |µ|−1!|Ψ|2+2 µ2ln1−|µ||Ψ|2#.(3.23) Localized solutions to eq. (3.21) are sought as Ψ(x,t)=U(x)eiωt, with a real function U(x), this solutions are usually referred to as envelope solitons in the continuum context [101]. The localized envelope U(x) obeys the equation d2U dx2=ω−2(1−|µ|)U2 1−|µ|U2U,(3.24) which may give rise to solitons, provided that ω > 0 and |µ|<1. The absence of soliton solutions for |µ|>1 implies that if the intersite self-defocusing, accounted for by µ < 0, is stronger than the onsite self-focusing, the self-trapping of solitons is impossible in the CA. Equation (3.24) can be cast in the form U′′ xx |µ|−1−1=−W′(U),(3.25) where the effective potential W(U) is W=−1 2U2−1−Ω 2|µ|ln1−|µ|U2,with Ω≡|µ|ω 1−|µ|; (3.26) the expansion of the potential (3.26) for U2→0 yields W≈h−ΩU2+|µ|(1−Ω)U4i 2.(3.27) 3.3. DB’S IN THE SALERNO MODEL WITH COMPETING NONLINEARITIES 51 This form of the equation shows that solitons exist in a finite band of frequencies, 0 < Ω<1, rather than in the entire semi-infinite band, Ω>0, where the linearization of equation (3.24) produces exponentially decaying solutions that could serve as the solitons’ tails. The reduction of the semi-infinite band to a finite one is typical for soliton families in models with competing nonlinearities, such as the cubic-quintic NLS equation [102]. Further, it follows from the divergence of potential (3.26) at U2=1/|µ|that the solitons’s amplitude A, which is a monotonously increasing function of Ω, is smaller than 1/p|µ| for 0 <Ω<1, and A=1/p|µ|at Ω = 1. Solitons can be found in an explicit form near the edges of the existence band: at small ω(i.e., small Ω), U(x)≈pω(1−|µ|)sech√2ωx,(3.28) while precisely at the opposite edge of the band, Ω = 1, the exact solution is a peakon, Upeakon(x)=1/p|µ|exp−p(1/|µ|)−1|x|.(3.29) In other words, at a given frequency ω, the peakon solution is found at |µ|=µp≡1/(1+ω).(3.30) Note that norm (2.21) of the peakon is Npeakon =π2/[6p|µ|(1 −|µ|)] ,(3.31) and its energy is also finite. Close to this point, i.e., for 0 <1−Ω≪1, the solution is different from the limiting form (3.29) in a narrow interval |x|.p|µ|/(1−|µ|)(1 −Ω), where the peak is smoothed. Finally let us remark that the CA based on eq. (3.21) is valid if the intrinsic scale of all continuum solutions, that may be estimated through the curvature of the soliton’s profile at x=0 as l∼1/qU′′ xx/U, is large, l≫1 (recall the lattice spacing is 1 in the present notation). According to eq. (3.29), the latter condition implies (1/|µ|)−1≪1 (i.e., strictly speaking, the CA applies in the case when the competing nonlinearities in the SM nearly cancel each other). It is relevant to note that, in the standard version of the SM (previously studied in sections 3.1 and 3.2), with µ > 0, the CA presented here give rise to pinned envelope solitons in the entire semi-infinite band, ω > 0 and then consistent with the exact solutions obtained for the discrete model in the above sections and earlier works [58, 59, 89, 90]. 3.3.2 Pinned discrete breathers In order to find exact pinned ((0, 1) resonant) discrete breather solutions in a numerical form, we look for solutions to eq. (2.18) which are localized and time periodic with frequency ωb=2π/Tb(that is related to ωin the continuum equation (3.21) by ωb≡ω−2). Pinned solutions are widely known for the DNLS limit (µ=0) since they can be obtained both continuing the analytical AL pinned beathers along the standard SM (as previously done in sections 3.1 and 3.2) and by the continuation from the anticontinuum limit, C=0, of the DNLS equation (2.5). It is then possible to make a numerical continuation of such 52 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES solutions for µ < 0 by adiabatic changes of the model parameter µand successive applications of the shooting methods in order to obtain the numerically exact pinned discrete breather for a given frequency ωband µ. In general all the pinned solutions were computed starting from the DNLS limit, µ=0, and increasing |µ|at a fixed value of ωb. The continuations were performed using an increment δ(|µ|)=10−2at each step, or smaller if higher accuracy was needed. As shown in the previous section, the breather family in the continuum equation (3.21) ends with the peakon solution (3.29). To compare the numerically determined shape of the discrete breathers with the feasible peakon limit, we fitted the breathers’ tails to the asymptotic form |Φn|=Aexp[−Γ(|n−n0|)],(3.32) with constant A,Γ, and n0, which follows from the linearized equation (2.18) for large |n|. This procedure yielded the decay rate, Γ = Γ(µ, ωb), amplitude, A=A(µ, ωb) (and the center’s position n0), as functions of parameters µand ωbof the pinned breather family. Once A(µ, ωb) and n0were found, we defined γ(µ, ωb)≡A−|Φn0|to measure a deviation of the true discrete soliton from a conjectured peakon shape obtained by formal extension of the tail inward. In figure 3.17.a we show the evolution of γproduced by several continuations of the discrete breather solutions (at different frequencies ωb). We define µp(ωb) as a value of µat which an exact discrete peakon of internal frequency ωbis found, that we realize as vanishing of γ(µ, ωb)at µ=µp. In figure 3.17.b we plot the evolution of the breather’s amplitude as the continuation is performed. It is observed that the amplitude increases with |µ|, reaching the predicted value, 1/p|µ|, at the exact peakon solution. A noteworthy result, evident from figure 3.17, is the persistence of discrete breathers beyond the peakon limit (which means continuability of the solutions to γ < 0). The apparent intersection of different curves at one point in figure 3.17.a is a spurious feature (see the inset in the figure): an accurate consideration shows that the curves actually intersect at close but different points. In contrast, the intersection of the curves in figure 3.17.b indeed happens at a single point, which corresponds to discrete breathers taking the peakon shape. Figure 3.18 displays typical examples of the numerically found discrete breathers. It demonstrates that the solutions corresponding to γ < 0 are cuspons, with a superexponential shape, that do not exist in the continuum equation (3.21). The discrete character of the SM with the competing nonlinearities allows this new type of solution, as happens with the quasi-collapsing states in the standard DNLS equation in two dimensions (see next chapter). Cuspon solutions continue into the region of |µ|>1, where the CA yields no breathers, but, due to the sharp change of the solution with the increase of |µ|, finding numerical solutions at larger values of |µ|becomes increasingly more difficult. In figure 3.19.a we compare the line of the existence of the peakons in the continuum limit, and the actual location of discrete peakons. It is seen that the agreement between the CA and numerical findings is good for smaller |ωb|(in this case, the discrete breathers are broad), while at larger |ωb|the discrete breathers are narrow, hence the agreement with the CA deteriorates. 3.3. DB’S IN THE SALERNO MODEL WITH COMPETING NONLINEARITIES 53 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0 0.5 1 1.5 2 γ |µ| ωb=-2.01 ωb=-2.09 ωb=-2.25 ωb=-2.51 ωb=-2.86 ωb=-3.33 -0.2 -0.175 -0.15 -0.125 -0.1 1.06 1.02 0.98 0.94 0.9 |µ| 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 |Φn0|2|µ| |µ|/|µp| ωb=-2.01 ωb=-2.04 ωb=-2.09 ωb=-2.16 ωb=-2.25 ωb=-3.08 ωb=-4.70 1 0.8 0.6 0.4 0.2 00.2 00.4 0.6 0.8 1 γ µ| | µ| | µ| p| | Φ n 0 | 2 | µ | (b) (a) Figure 3.17: (a)The mismatch with the peakon shape, γ, as a function of |µ|, for discrete breathers found at different frequencies ωb. (Note in the inset that there is no common intersection of all the curves). (b)The breather’s amplitude vs. |µ|. The axes are rescaled to show that the amplitude of the peakon solutions (attained at |µ|=|µp|) are equal to 1/p|µ|, as predicted by the continuum approximation. Floquet analysis Performing the linear stability analyses of pinned breathers it is found that these solutions are linearly stable along the whole µ-continuation, except for a relatively small region, as shown in figure 3.20.a. The entire instability island in the (ωb,|µ|) plane is displayed in figure 3.19.b. The instability displayed is revealed by a Floquet multiplier leaving the unit circle at +1 (harmonic bifurcation). The eigenvector associated to this multiplier show a localized profile around the pinned solution. Note, in particular, that the peakon and cuspon solutions are stable. The stability of the discrete breathers was also checked by direct simulations of perturbed (along the unstable direction given by the Floquet eigenvector whose Floquet multiplier is λ > 1) solutions, using the full equation (2.18). The 54 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES 10-10 10-8 10-6 10-4 10-2 100 80 100 120 140 160 180 200 220 |Φn|2 n µ=-0.956 µ=-2.640 µ=-0.300 10-2 10-1 100 165 150 135 |Φ | n2 Figure 3.18: Generic examples of three different types of discrete breathers of frequency ωb= 2.091: a quasi-continuous sech-like solution at |µ|=0.3, an exact peakon at |µ|=0.956, and a cuspon at |µ|=2.64. |µ|=0.884 − 0.4 0.2 0 7 6 5 4 3 2 (a) (b) (c) µpp µ CL p µ ωb | | | | | | Figure 3.19: (a)The value of µp, at which the soliton assumes the peakon shape: the prediction of the continuum approximation, equation (3.30) (solid curve), and numerical results for discrete breathers (dots). The small region where the pinned discrete breathers are found to be unstable (for that purpose, the vertical axis shows |µ|, rather than |µp|) is also shown. The inset displays the relative difference between the numerically found values of |µp|and the predicted ones, |µCL p|, provided by the continuum approximation. (b)Zoom of the area in the ωb,µpplane where the instability island is located. (c)Norm of the discrete breathers vs. their frequency for |µ|=0.884. 3.3. DB’S IN THE SALERNO MODEL WITH COMPETING NONLINEARITIES 55 1.2 1 0.8 0.6 0 100 200 300 400 500 time n0 n0+1 n0+2 |Φ ( )| |Φ ( )| t 0 n n λF (0) Φn | | Φn| time 3002001000 400 500 1.2 1 0.8 0.6 | | (b) (a) |µ| (t)| Figure 3.20: (a)The absolute value of the Floquet multipliers, {|λ|}, for the linearization of perturbations around discrete solitons, is shown vs. |µ|, at several fixed values of the frequency (which are chosen so as to make the instability intervals well separated). (b) A robust pulson generated from an unstable breather at |µ|=0.922. results of these simulations corroborated the predictions of the linear stability analysis. Direct simulations of the evolution of perturbed unstable breathers, a typical example of which is displayed in figure 3.20.b, show that, after a transient stage, a localized pulson (showing simultaneous width and amplitude oscillations) is formed. The pulsons are (quasi-)periodic in time, and persist indefinitely. This behaviour resembles that found in the ordinary two-dimensional DNLS equation in quasi-collapsing states that will be described in next chapter. A necessary stability condition for soliton families in models of the NLS type may be provided by the Vakhitov-Kolokolov (VK) criterion [103]: if the norm Nof the breather is known as a function of its frequency ωb, the breathers can be stable against small perturbations with real eigenvalues, provided that dN/dωb<0. Although the applicability of the VK criterion to the present model has not been proven (and counter-examples are known, when solutions predicted by the criterion to be unstable are actually stable [104]), it is relevant to test the criterion here, numerically computing N(ωb)according to eq. (3.19). The result is that the VK criterion precisely explains the stability and instability of the discrete breathers, except for the cuspons (see below), as shown in figure 3.19.c. A noteworthy feature of the N(ωb)dependence is a divergence of the total norm due to the infinite contribution of the central site to expression (3.19) in the case of the exact peakon solution, with Φn0 2=1/|µ|. An example of the N(ωb)dependence showing the divergence is plotted in figure 3.21. As concerns the cuspons, whose amplitude exceeds the critical value, 1/p|µ|, the norm (3.19) converges for them, and features a positive 62 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES (i) The set of “s” plane waves which take part in the background of a (p,q)−resonant discrete breather with internal frequency ωbis derived by the simple selection rule for the wave-numbers kj ω(kj,Aj) ωb =1 qp 2πkj−m,(3.35) i.e. only the plane waves which are (p,q)-resonant with the internal period of the breather can be components of {Φbckg n(t)}. The number of solutions of (3.35) fixes “s”. (ii) The amplitudes {Aj}of the nonlinear plane waves differ by orders of magnitude yielding a localization in the k-space. (iii) There exist a strong positive correlation between the amplitude of the background and the strength of the Peierls-Nabarro barrier arising from the periodic lattice. This correlation is particularly clear when symmetry breaking transitions occur for the also studied case of ν < 0 and µ > 0, and reflects the link between non-integrability and the existence of the background dressing of the mobile core. Another interesting effect is obtained for the SM with competing nonlinearities. In this case the continuation of mobile breathers of a given frequency stops near the divergence of the Peierls-Nabarro barrier for pinned breathers with the same frequency. (iv) Finally, the interpretation of the correlation described in (iii) is reinforced from a study of the energy evolution of the mobile core: There is an energy balance brought by the background when the core moves along the lattice. In particular, it can be observed how the core energy oscillates periodically so that it takes the maximum energy value when the core visits the inter-site configuration. This extra energy periodically obtained by the core is provided by the interaction backgroundcore, with the energy maximum clearly related to the background amplitude. It is worth stressing that the most relevant predictions of perturbative collective variable theory are confirmed by our numerical results, which show the existence of PeierlsNabarro barriers to breather translational motion. Furthermore, the existence of exact oscillating breather solutions for the standard SM is numerically confirmed. They are found to contain an extended background whose amplitude is typically much smaller than for mobile breathers. The correlation between the Peierls-Nabarro barrier EPN (computed from immobile breathers) and the amplitude background of moving breathers correctly suggests that the background has a role in the energy balance required to overcome the barriers to translational motion. The interpretation is also fully consistent with the observations on the background amplitude behaviour of spatially oscillating anchored breathers in the standard SM. Currently used effective particle (collective variable) theories are thus seen as intrinsically incomplete, because core energy is not an invariant of motion. Any sensible improved approach must adopt equation (3.33) as starting point for improved perturbative ansatzes, and we hope that our work will stimulate further studies along these lines. Numerically exact moving discrete breathers with an infinitely extended tail of small amplitude were already observed in some cases for Klein-Gordon lattices with Morse 3.4. CONCLUSIONS AND PROSPECTIVE REMARKS 63 µ ν AL DNLS Standard SM ✗ Staggered behaviour SM CompetingNonlin. Mobile Breathers CUSPONS ✗ SBB ? Figure 3.27: Schematic plot of discrete breather’s existence diagram in the (ν,µ)-plane. For the standard Salerno model µ > 0 we have continued pinned breathers along the path µ+ν=1 for (ν > 0) from the Ablowitz-Ladik integrable lattice to the DNLS limit. Mobile breathers can be continued all the way (from AL to DNLS eqs.) along this path provided their frequencies are small enough. For the standard Salerno model with ν < 0 we have performed the continuation along the path µ−ν=1 finding a symmetry breaking bifurcation for pinned breathers at ν≃ −0.3 that prevents continuing both pinned and mobile solutions far beyond this point and then excluding the possibility of reaching the DNLS equation with ν < 0. The staggering transformation between the two regions (ν > 0 and ν < 0) of the standard Salerno model implies that staggered breathers cannot be continued to the DNLS limit with ν > 0 because of a symmetry breaking bifurcation. For the Salerno model with competing nonlinearities (µ < 0) we find a transition from smooth peaked pinned breathers to cuspon states where the energy is hyperlocalized around the breather center. Cuspon breathers appear as stable solutions of the dynamics. Mobile breathers in the competing Salerno mobile cannot be continued beyond this transition point since the Peierls-Nabarro energy diverges at this point. The staggering transformation implies that these latter results applies for staggered breathers when µ < 0 and ν < 0. However, this region could not be explored for typical breather states since, on one hand, the purely AL lattices with µ < 0 does not admit localized states as solutions (preventing the continuation from ν > 0 and µ < 0 ) and, on the other hand, the continuation from ν < 0 and µ < 0 stops near the symmetry breaking bifurcation as reported above. Then, this region of the Salerno model with competing nonlinearities remains unexplored and apparently forbidden for our continuation methods. potential by Cretegny and Aubry [53], however no investigation of the background of these exact solutions is reported, so they were able to ”..suggest that generally a strictly localized breather cannot propagate without radiating energy”. Our systematic study of the NLS lattices allows us to go further by showing that the extended background (here fully characterized) plays an important and subtle role in the translational motion of the localized core. Indeed, it is an indispensable part of the exact solution in the nonintegrable regime. Exact mobile localization only exists over finely tuned extended states of the 64 CHAPTER 3. DISCRETE BREATHERS IN 1D NONLINEAR SCHR ¨ ODINGER LATTICES nonlinear lattice. Mobile ”pure” (i.e. rest state background) localization must be regarded as very exceptional [96]. Before concluding this chapter, it is worth commenting on some of the differences between the Newton continuation of fixed points that we use in this chapter, and other important recent approaches to breather numerics. The work by Ablowitz et al [107] uses discrete Fourier analysis to obtain a nonlinear nonlocal integral equation, from where the ” ... soliton is thus viewed as a fixed point of a nonlinear functional” (sic) in the Fourier transformed space of functions. Following these authors, their results seem to differ from those of early pioneering work [108] (nowadays textbook material [23]) ”in which a continuous traveling solitary waves were reported using Fourier series expansions with finite period L while assuming convergence as L → ∞” (sic). Ablowitz et al term continuous a solution that can be defined offthe lattice points, which they see as ”necessary when discussing traveling waves in lattices” (sic), and disagree with some conclusions reported in the earlier works. The (”orthodoxy matters”) discussion above helps us to clarify how differently our numerical approaches ”sees” the discrete Schr¨ odinger breather problem: The very concept of a variable defined offthe lattice points is intrinsically alien to our discrete approach, which neither needs of it nor excludes its eventual consideration. In contrast to those views (but not at all in logical opposition), we consistently view the thermodynamical limit (N→ ∞) in lattice space, much in the sense used e.g. by Serge Aubry in his celebrated work on the Frenkel-Kontorova ground state problem [49]: The infinite size limit is built up from a subsequence of PBC (finite) lattices for which the limit is well defined. This will make the Fourier-transformed k-space continuum. Closer to our approach in some respects, though technically different in many others, is the formal approach purposed recently by James and collaborators [109, 110]. It is also worth mentioning that these results have been reproduced recently for other kind of solutions (dark breathers) [111] and have constituted [112] a “(negative) result” about the impossibility of constructing “exponentially localized fundamental (single-humped) moving discrete solitons” in the nonintegrable part of the Salerno model. There are, at very different levels, several open questions to further research. From a technical point of view, it is important to analyze carefully the irrational limit p/q→σ, of the solutions. In particular, in this limit the number of resonant plane wave branches tends to a continuum and one could (or not) expect that exponential localization in the reciprocal lattice persists in that limit. This can be addressed numerically, though systematic investigations may require some efforts in optimizing the time efficiency of current numerical schemes. An important issue regarding applications is the phenomenology of multibreather states. In particular, studies on collisions of a pair of breathers may find in this study of exact mobility a useful reference in order to deal with the complexities that emerge from the many time-length scales involved in these physically relevant phenomena. Much simpler multibreather states, e.g. train-like chains of (moderately) separated moving breathers could also be investigated. Not least, the perspective and results presented here may be of some interest to studies of the effects of coupling to (nonthermal and/or thermal) radiation baths in the breather and multibreather states of nonlinear lattices [51] and the practical manipulation and patterning of localized ”hot spots” by external fields [113]. Chapter 4 Discrete Breathers in two-dimensional Nonlinear Schr¨ odinger lattices Given the ubiquity of such breathers in discrete nonlinear physical systems (which exist on essentially all length scales), these nonlinear excitations are likely to be important in many physical phenomena, including melting, fracture, and the buckling and folding of biopolymers. They may also prove useful in technologies ranging from ’smart’ materials with tunable collective responses to light-induced, alloptical switches and networks. With the acquisition of this new animal, the nonlinear ’zoo’ has become an altogether more interesting place. David K. Campbell in [114]. The study of two-dimensional nonlinear Schr¨ odinger lattices has attracted much attention [115, 116] in recent years due to the new phenomena emerging when the dimensionality of the lattice is increased. Some examples of these new features are the existence of vortex-breathers [117] which supports energy flux, the appearance of an energy threshold for the creation of discrete breathers [118–122] and the ubiquity of an instability (the quasi-collapse) of some discrete breather solutions leading to a highly localized pulson state [123–128]. These theoretical efforts have their counterpart in recent advances in the field of nonlinear optics. The studies of two-dimensional arrays of coupled nonlinear waveguides allow the experimental observation of those effects studied theoretically. Specially relevant is the recent experimental breakthrough (theoretically designed in [129]) by Fleisher et al [72, 130], where a two-dimensional array of nonlinear waveguides is induced in a photosensitive material. This technique provides a clear experimental verification of the two-dimensional discrete breather existence in this system. In particular, besides the observation of standard discrete breathers, these works reported the first observations of staggered discrete breathers. Our study in this chapter focus on the computation of numerically exact discrete breathers in two-dimensional anisotropic nonlinear Schr¨ odinger lattices, i.e. where the couplings in the two spatial directions are not necessarily equal. The use of the shooting methods introduced in section 2.3.1, and redefined here for the two-dimensional case in section 4.1, allow us to find these solutions and analyze their structural and stability properties. Both pinned and mobile discrete breathers are studied. In the latter case we will study only the ones whose motion is along one axis of the lattice. The analysis of 65 66 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES the numerically exact solutions help to shed light on some features of the properties and stability of localized solutions reported in previous works. After introducing in section 4.1 the two-dimensional anisotropic Salerno lattice and provide explanations on the implementation of the numerical procedures used to study the dynamics of 2D discrete breathers, we will focus on pinned ones. The analysis of the results on pinned discrete breathers for anisotropic nonlinear Schr¨ odinger lattices is reported in section 4.2 for both the standard version of the SM (section 4.2.1) and that with competing nonlinearities (section 4.2.2). In both studies we present the numerical computations of the fixed point norm, as a function of three parameters: breather frequency, transversal coupling, and nonlinearity (see below). They show, as anticipated, the socalled quasi-collapse transition. In these studies we present numerically computed sectors of the bifurcation surface and take a brief look at the nonlinear dynamics on the unstable manifold, whose typical trajectories have been called pulson states. Early numerical work on the 2D quasi-collapse phenomena in isotropic lattices was reported in [127, 128] and [125]. A three-year-old account of the ”state of knowledge” on 2D Schr¨ odinger lattices can be found in Section six of [131]. Interestingly, for the case of competing nonlinearities a transition to 2D cuspon states is also found. In this region of the Salerno model we have also studied the existence and stability of in-phase and out-of-phase bound states of pinned breathers motivated by the results obtained in the previous chapter for the 1D case (section 3.3.2). As introduced above, a new class of breathing solutions are possible in the 2D model: discrete vortices [117]. We investigate vortex breathers of two types, vortex crosses and vortex squares, in section 4.3 (in the framework of the isotropic model). The analysis of their linear stability reveals parametric stability regions (which turn out to be rather narrow) for the vortices, and helps to identify various bifurcations (including a generic Hamiltonian Hopf bifurcation) responsible for their destabilization. Direct simulations demonstrate that the instability transforms the vortices into ordinary breathers in the case of the standard Salerno model, and into vortex pulsons, that keep the vortical topology, in the most interesting case of competing nonlinearities. Finally, we have also introduced bound states of vortex crosses and analyze their stability. Mobile solutions are finally reported in section 4.4. For this type of solutions we have focused on a single type of mobile breather, namely those moving along the direction of stronger lattice coupling constant. The structure of each of these mobile exact discrete breathers is that of a localized moving core superimposed on a specific extended state of resonant small amplitude radiation, the background. An extensive Floquet stability analysis of this type of solutions is performed in two sectors of the three-dimensional parameter space, revealing the existence of two different transitions. The tangent space eigenvectors associated to each of the transitions are presented, and the relation of the unstable manifold trajectories to pulson states is analyzed afterwards. 4.1 The Salerno model in two dimensions Motivated by the results reported in the last chapter our aim here focus on extending the continuation scheme for calculating exact discrete breathers in higher dimensional systems. In particular we focus on the following two-dimensional nonlinear Schr¨ odinger 4.1. THE SALERNO MODEL IN TWO DIMENSIONS 67 lattice i˙ Φnm =−C1(Φn+1,m+ Φn−1,m)+C2(Φn,m+1+ Φn,m−1)(1+µ|Φn,m|2)−2νΦn,m|Φn,m|2(4.1) This lattice can be viewed as the two-dimensional Salerno model. The two coupling parametersC1and C2provide a technical advantage for numerics (see below), but they are also introduced for theoretical and experimental interest. The possibility of controlling the ratio between the two linear couplings of the two transversal directions has been studied in various works as a way of analyzing how the intrinsic 2D phenomena (such as the quasi-collapse) emerge. In fact, for C1<< C2,µ=0 and ν,0 equation (4.1) describes a set of weakly coupled nonlinear waveguide arrays and can be considered as a case of “intermediate dimensionality”. This extreme has been studied experimentally in [132] and using perturbative methods in [133]. On the other hand, this equation incorporates, as two particular limits, the physically relevant standard two-dimensional DNLS equation (µ=0, ν,0) and the two-dimensional counterpart of the AL lattice (µ > 0, ν=0) which is not integrable. The continuation between these two limits provides a useful tool for studying the interplay between the on-site and inter-site nonlinearities in the 2D case. Moreover, the anisotropy (or freedom in the values of the coupling parameters C1and C2) allows to include an integrable model among the members of the family of nonlinear lattices described by eq. (4.1). That is, for ν=0, Ci=0 and Cj,0 one obtains a set of integrable AL 1D chains. In this sense, every 2D model included by eq. (4.1) is connected with this integrable model where analytic discrete breathers are available. The SM (4.1) may find a straightforward physical application as a discrete model for the BEC of dipolar atoms trapped in a deep two-dimensional optical lattice [78]; in that case, as stated for the 1D Salerno model, assuming that a strong magnetic field aligns the momenta parallel (perpendicular) to the lattice plane, and the condensate is strongly confined in the vertical direction, one will again deal with the dipole-dipole attraction (repulsion), i.e. µ > 0 (µ < 0) in eq. (4.1). Similarly to the 1D version of the Salerno model eq. (4.1) has two dynamical invariants, the Hamiltonian H=−C1X n,mΦn,mΦn+1,m+ Φn+1,mΦn,m −C2X n,mΦn,mΦn,m+1+ Φn,m+1Φn,m −2ν µX n,m|Φn,m|2+2ν µ2X n,m ln1+µ|Φn,m|2,(4.2) and, due to the phase invariance of the equations of motion, the following norm (4.1) N=1 µX n,m ln1+µ|Φn,m|2.(4.3) Note that we have included here the needed redefinition in the logarithmic terms of both quantities in order to manage with a correct description of the dynamical invariants within the Salerno model with competing nonlinearities. In the same manner as in the 1D case we will focus on a special set of 2D discrete breathers. For this, we have to generalize the definition (2.22) introduced in section 2.2 for 68 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES a (p,q) resonant solution in the 1D model to the 2D case. In this context, discrete breathers solutions are characterized by three time scales. Namely, one associated with the internal oscillation ωband the other two derived from the translation of the localization center, i.e. its velocity~ vb=(vx,vy). The subset of 3-tuples (ωb,~ vb) that fulfill the (px,py,q)-resonance condition vx2π ωb =px q(4.4) vy2π ωb =py q,(4.5) (where px,pyand qare integers) denote the breather solutions that can be obtained with our continuation method. These solutions are those that after qperiods of the internal frequency, ˆ Φ(t0+qTb), translates pxand pylattice sites in the xand ydirection of the square lattice, respectively, i.e. ˆ Φn,m(t0)=ˆ Φn+px,m+py(t0+qTb),(4.6) where, again, PBC are applied ΦNx+1,m= Φ1,m,Φ0,m= ΦNx,m,Φn,Ny+1= Φn,1and Φn,0= Φn,Ny(with Nxand Nybeing the lattice size in the xand ydirection respectively). Consequently, a (px,py,q)-resonant state will be a solution of the following set of equations F(px,py,q,ωb,ν,C1,C2)h{ˆ Φn,m(t0)}i=Lpy yLpx xTq (ωb,ν,C1,C2)h{ˆ Φn,m(t0)}i={ˆ Φn,m(t)},(4.7) where the operators Liare the lattice translation in the i-direction, Lx[{Φn,m(t0)}]={Φn+1,m(t0)},(4.8) Ly[{Φn,m(t0)}]={Φn,m+1(t0)}.(4.9) Besides, T(ωb,ν,C1,C2)is the time evolution operator given by equation (4.1) over one period Tb=2π/ωb,T(ωb,ν,C1,C2)[{Φn,m(t0)}]={Φn,m(t0+Tb)}.(4.10) In order to illustrate the 2D time scales resonance let us to consider the plane wave solutions of equation (4.1): Φn,m(t)=Aexp[i(kxn+kym−ωt)]. These solutions possess the following nonlinear dispersion relation ω(~ k,A)=2(C1coskx+C2cosky)(1 +µA2)−2νA2.(4.11) Hence, we can obtain the subset of plane waves which are (px,py,q)-resonant with some time scale τ(i.e. after a time qτthey have translated pxand pysites in the xand ydirection, respectively). Each member of these subsets will be labeled by the pair ~ k=(kx,ky) and from the condition (4.7) it follows that the corresponding set of values of~ kfor each family will satisfy the relation ω(~ k,A)=1 qτ~ p·~ k−m 2π,(4.12) where mis an integer and ~ p=(px,py). In figure 4.1 the corresponding values of ~ kare represented for two resonances of type (px=1,py=0,q=1) and (px=1,py=1,q=1). 4.1. THE SALERNO MODEL IN TWO DIMENSIONS 69 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 C2=0.1 C2=0.5 C2=1.0 C2=1.5 -3 -2 -1 0 1 2 3 -3 -2 -1 0 1 2 3 C2=0.1 C2=0.5 C2=1.0 C2=1.5 ky kxkx (a) (b) Figure 4.1: Wave numbers, ~ k=(kx,ky), of the (1,0,1) (a) and (1,1,1) (b) resonant plane waves for m=0 (see equation (4.12)). Different values of C2, while C1is fixed (C1=1), are shown. The reference time scale for the resonance is set to τ=2.4315 (ω=2.584). The method used for solving equation (4.7) for each resonant 3-tuple (ωb,~ vb) is the same as in the 1D case, already described in section 2.3.1. Then, the implicit function theorem assures that a fixed point solution of a map (4.7) given by ~ ξ=(px,py,q,ωb,ν, C1,C2) can be obtained provided that (i) the Jacobian of the operator F~ ξ[{Φn,m(t0)}]−I is invertible, and (ii) we know a fixed point of a map corresponding to an infinitesimally close set of parameters, ~ ξ−δ~ ξ=(px,py,q,ωb−δωb,ν−δν,C1−δC1,C2−δC2). As explained in section 2.3.1 the first demand can be satisfied using a singular value decomposition (SVD) of the Jacobian in order to obtain the pseudo-inverse operator. On the other hand, when the second condition is fulfilled convergence of the Newton-Raphson iterative scheme is guaranteed. For this, we start with a sufficiently good trial solution, {Φ0 n,m(t0)}and solve the equation {δΦ0 n,m(t0)}=−DF~ ξh{Φ0 n,m(t0)}i−1·F~ ξh{Φ0 n,m(t0)}i,(4.13) in order to obtain {Φ1 n,m(t0)}={Φ0 n,m(t0)}+{δΦ0 n,m(t0)}. We iterate these calculations to the desired convergence, and then the solution, {ˆ Φn,m(t0)}, is obtained. In our numerics this is the case when F~ ξh{Φi n,m(t0)}i<N·10−16 ,(4.14) (where Nis the total number of sites in the square lattice) is fulfilled. Once the solution is found we use it as the following trial solution, {Φ0 n,m(t0)}, for solving the map (4.7) corresponding to the next set of parameters ~ ξ′=~ ξ+δ~ ξ. There are two possible paths for developing the continuation method depending on the choice of the starting point of the continuation. One possibility is to start from the full anti-continuum limit, C1=C2=0, where a pinned breather solution of frequency ωbis written as ˆ Φn,m(t)=δn,n0δm,m0rωb 2νexp(iωbt).(4.15) Starting from the above solution, we can perform the continuation increasing the parameters C1and C2as usual, and so obtain the whole family of (px=0,py=0,q=1) resonant 70 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES discrete breathers. An alternative path starts from the one-dimensional limit, C2=0. The choice of this second limit (which implies taking as the very initial trial solution of the continuation the whole set of 1D solutions obtained in the previous chapter) is justified when seeking mobile solutions. As stated above, this limit offers the possibility of studying strongly anisotropic lattices as a controlled interpolating situation between one and two dimensions. On the other hand, employing this strategy we can only obtain those solutions which are (px=p,py=0,q)-resonant, i.e. the two-dimensional continuation of those one-dimensional (p=px,q)-resonant discrete breathers. Hence, the solution from which we start is ˆ Φn,m(t)=δm,m0ˆ Φ1D n(t),(4.16) where ˆ Φ1D n(t) is the corresponding (p=px,q)-resonant one-dimensional solution. In what follows we will employ both continuation paths when we study the case of pinned breathers (section 4.2), and we will show that the results obtained are the same when approaching the same limit (the standard two-dimensional DNLS). 4.2 Pinned discrete breathers We first focus on the characterization of pinned ((0,0,1)-resonant) discrete breathers for the standard Salerno model (with special attention to the DNLS equation) in section 4.2.1 and for the SM with self-defocusing inter-site nonlinearity in section 4.2.2. 4.2.1 Pinned discrete breathers in the standard Salerno model As we have discussed, we can choose two different starting points for the continuation of (0,0,1)-resonant fixed points (pinned breathers) of equation (4.7): (i) the full anticontinuum (AC) limit (C1=C2=0), or (ii) the (one-dimensional, 1D) limit of uncoupled chains (C1,0,C2=0), where they were obtained in the previous chapter from continuation along the standard 1D Salerno model by increasing values of the parameter νfrom the one-dimensional AL lattice (2.9). As a test for our codes, we have checked that both paths arrive to the same solution. In fact, unique continuations can proceed along any path on the plane of parameters (C2, ν) that we have explored. Early works [123–125] on the isotropic two-dimensional standard DNLS equation analyzed the so-called quasi-collapse instability of pinned discrete breathers, i.e. the condensation of all the energy into a few modes in discrete nonlinear systems, which corresponds to the onset of a singularity (wave collapse) [126] in multidimensional continuum models. Subsequent numerical works [127] extended these studies to the isotropic 2D Salerno lattice and addressed the question of how the instability is affected by the presence of impurity lattice sites. As expected, our results further corroborate the existence of quasi-collapse instabilities in the anisotropic case: The phase diagram in parameter space (ωb,C2,ν) consists of two regions (stable and unstable) separated by the surface of transition. As we perform the continuation of breather solutions across the parameter space we scan the Floquet stability of the computed solution. In figure 4.2 we present the two stability transition curves in the plane (ωb,C2,ν=1), i.e. the function Cth 2(ωb), corresponding to the two different continuation starts. The continuation from the AC limit is made through the path C1=C2 4.2. PINNED DISCRETE BREATHERS 71 0 0.2 0.4 0.6 0.8 1 5 4 3 2 1 0 C2th ωb C1=C2 C1=1 C ω 2 th b Figure 4.2: Evolution of the threshold value of the coupling parameter, Cth 2, as a function of the frequency, ωb, for two different continuations starts. The values of Cth 2limit the region where pinned discrete breathers are linearly stable (unstable for C2>Cth 2). The instability yields a hyperlocalized state (quasi-collapse). The continuation from the fully uncoupled limit (C1=C2=0) (filled circles) is performed using the path C1=C2. For the continuation (bold circles) from the 1-dimensional limit (C1=1, C2=0) the coupling in the new direction C2is progressively increased. and the one from the 1D limit is made at C1=1. The convergence of the two paths at C2=1 is clearly seen. The Vakhitov-Kolokolov criterion [103] for stability of the pinned discrete breather solution derived and used for the 2D DNLS in [124, 125], ∂N ∂ωb!C2,ν >0,(4.17) is of a very general character and our numerics illustrate it clearly. On the other hand, the Floquet stability analysis detects the dimensionality (and a basis in tangent space) of the unstable linear manifold associated with the quasi-collapse instability that these exact discrete breathers experience for some parameter values. We have computed numerically, for a fine grid of ωbvalues and a coarser grid of C2and ν, the function N(ωb,C2, ν), from which we show some sectors in figures 4.3 and 4.4. In figure 4.3 we show the numerically computed norm (4.3) as a function of the breather frequency N(ωb), for three different values of the transversal coupling C2, and a fixed value of ν=1 (anisotropic DNLS limit). We observe the existence of a minimum value, min N(ωb)=Nth ,0, which is thus seen as an excitation threshold for the creation of these solutions. The position of the minimum ωth b(C2), which naturally increases with C2, separates the stable and unstable branches of pinned breathers. Breathers corresponding to values of ωbwhere N(ωb) has a negative slope are unstable: This is shown in the insets, where the Floquet spectra of two representative examples of pinned discrete Schr¨ odinger breathers are plotted in the complex plane. Note that the high accuracy of the numerical solution allows an unprecedented detailed Floquet analysis of the instability, paving the way to rigorous analytical characterizations of the quasi-collapse unstable 78 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES Bound states of discrete breathers In addition to isolated pinned breathers, we have studied their bound states when µ < 0. Two types, in-phase and π-out-of-phase, of pairs of identical breathers, with the same frequency ωband different distances between them, has been analyzed. For this purpose, we first continued these solutions, at µ=0, from the anti-continuum limit up to the 2D DNLS equation (C=ν=α=1), and then decreased the value of µinto the region of competing nonlinearities (µ < 0). At the same time, the linear stability analysis of these periodic solutions was performed by the numerical computation of their Floquet spectra. We have computed two different patterns of bound states of breathers. The first type consists of two discrete breathers with their centers, (n(j) 0,m(j) 0), with j=1,2, lying on the same lattice axis (so that n(1) 0=n(2) 0or m(1) 0=m(2) 0), whereas for the second type of bound states the centers are related by n(1) 0=n(2) 0±dand m(1) 0=m(2) 0∓d,i.e. they are aligned along a diagonal of the lattice. In figure 4.9 we show the absolute value of the Floquet multipliers as a function of µfor in-phase and out-of-phase bound states, aligned along a lattice axis for the case of ωb=7.0, with three different values of the distance between breather centers in the pair. Results of similar computations for the diagonal-aligned bound states with ωb=8.0 are shown in figure 4.10. As in the 1D version of the model (see section3.3.2), for µ=0 in-phase bound states are linearly unstable (the more unstable the closer breathers are in the pair), while out-of-phase pairs are stable. As observed in figures 4.9 and 4.10, at µ=−0.3 for the pattern of the first type (ωb=7.0), and at µ=−0.25 for the second one (ωb=8.0), the in-phase bound states become stable regardless of the distance between breathers. Simultaneously, out-of-phase states become unstable, also regardless of the 1.8 1.4 1 0.6 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ d=3 1.8 1.4 1 0.6 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ d=3 d=5 1.8 1.4 1 0.6 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ d=3 d=5 d=7 1.4 1.2 1 0.8 0.6 |λ| d=3 1.4 1.2 1 0.8 0.6 |λ| d=3 d=5 1.4 1.2 1 0.8 0.6 |λ| d=3 d=5 d=7 Figure 4.9: The absolute value of the Floquet multipliers as a function of µfor in-phase (top) and out-of-phase (bottom) axis-aligned bound states of breathers with ωb=7 (C=1). The figure shows cases when the two breather centers are separated by d=3, 5 and 7. It can be observed that, irrespective of the value of d, the stability interchange occurs at µ=−0.3. 4.3. DISCRETE VORTEX BREATHERS 79 4 3 2 1 0 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ d=1 4 3 2 1 0 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ d=1 d=2 4 3 2 1 0 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ d=1 d=2 d=3 5 4 3 2 1 0 |λ| d=1 5 4 3 2 1 0 |λ| d=1 d=2 5 4 3 2 1 0 |λ| d=1 d=2 d=3 Figure 4.10: The same as in the previous figure for in-phase (top) and out-of-phase (bottom) diagonal-aligned bound states of breathers with ωb=8 (C=1). Shown are the results for the states with the separation between the breather centers d=1, 2 and 3. It can be observed that for all these states undergo the stability exchange at µ=−0.25. separation between breather centers. The same stability exchange between inand out-of-phase states was observed in the 1D case, where it occurs at the value of µat which the discrete breather solution is a peakon. However, here in the 2D case the discrete breathers in the pair are cuspons on both sides of the stability-exchange point. Nevertheless, we find, for both types of the bound states, that the values of µat this point is exactly the same at which the cuspon’s norm, N(ωb, µ), diverges (see the previous subsection). In other words, the stability interchange between inand out-of-phase bound states is associated with the divergence of the breather norm N(ωb, µ), rather than to the appearance of a peakon (contrary to the 1D case, where the emergence of a peakon and norm divergence occur simultaneously). As a conclusion, although the divergence of the norm does not switch the stability of single pinned discrete breathers, it marks the stability border of bound states of breathers, regardless of their size and orientation relative to the lattice. 4.3 Discrete vortex breathers A natural generalization of the fundamental discrete breathers are discrete vortices, which are well-known solutions of the ordinary 2D DNLS model [117]. A vortex is characterized by the phase circulation around its center, ∆θ, that must be a multiple of 2π. Hence they may be labeled by an integer number (vorticity, or topological charge), S≡∆θ/(2π). In this section, we consider vortices only in the isotropic model (C=1), with the purpose of analyzing their behaviour at both the standard SM (µ > 0) and at the SM with competing nonlinearities (µ < 0). In the framework of the 2D DNLS model, influence of the lattice anisotropy on fundamental and vortical discrete breathers was studied in [134]). 80 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2 1 0 0464 6 2 1 2 0 02466 0 2 4 8 −1 −2 1 0 2 −2 −1 0 1 2 20 nm nm Re( Re(Φn,m) Φn,m) −2 −1 −2 −1 Figure 4.11: Two examples of fundamental (|S|=1) discrete vortices. Profiles of the real part of the square vortex with M=1 and vortex cross are shown in the top and bottom panels, respectively. Both solutions are found for µ=−0.4 and ωb=7.0 (as noted in the text, we fix C=1 for the vortex solutions). We will construct two types of vortices, on-siteand off-site-centered ones (alias vortex crosses and vortex squares), both with |S|=1. Vortex squares are characterized by the number of lattice bonds, M, that each side of the square comprises; in this section, we only deal with M=1. Two examples of these two species of the solutions are plotted in figure 4.11. 4.3.1 Vortex crosses In order to construct fundamental (|S|=1) vortex crosses centered around the lattice site (n0,m0), we start with the anticontinuum (C1=C2=0) DNLS (µ=0) limit. The corresponding seed pattern includes nonzero fields Φn0,m0+1=−iΦn0+1,m0=−Φn0,m0−1=iΦn0−1,m0=pωb/2.(4.18) Then, by adiabatically increasing the inter-site coupling (Newton continuation in Cusing C1=C2=C), we reach the isotropic DNLS model, and start the continuation to positive values of the inter-site nonlinearity, µ. Performing the continuation in Cat µ=0, we have found that, for low-frequency vortex solutions, there is a critical value, Cc, that depends 4.3. DISCRETE VORTEX BREATHERS 81 8 7 6 5 2 1 0-1-2 µ 8 7 6 5 2 1 0-1-2 µ 10 7.5 5 2.5 0 |λ| ωb Unstable Unstable Stable =8.0 b ω µ 2 1 0 -1 -2 2 1 0-1-2 Im(λ) Re(λ) 2 1 0 -1 -2 Im(λ) (a) (c) (b) (d) µ=0.46 µ=−0.3 Figure 4.12: (a) The absolute value of the Floquet multipliers as a function of µfor a vortex cross with ωb=8. Two bifurcations can be inferred from the Floquet distributions in panels (c) and (d): a Hamiltonian Hopf bifurcation at µ=0.46, and a harmonic bifurcation at µ=−0.3. A similar set of two bifurcations is found at other frequencies. The entire stability diagram is displayed in panel (b), showing a narrow stability region. (As noted in the text, we fix C=1 for the vortex solutions). on frequency ωb, at which a Hamiltonian Hopf bifurcation (HHB) [135] occurs and the vortex solution turns unstable for C>Cc(ωb). This phenomenon was already reported in previous works [117, 134]. Higher-frequency vortex solutions, which are stable in the DNLS equation in the considered range of parameters, undergo destabilization through a bifurcation of the same type as a result of the continuation in µ, at C=1. The Hamiltonian-Hopf character of the bifurcation can be seen in figure 4.12.c, which shows the Floquet spectrum after the bifurcation: it is seen that a quadruplet of complex eigenvalues λjexit the unit circle. After this (first) bifurcation, further bifurcations of the same type occur at increasing values of µ, as observed in the right part of figure 4.12.a. Similar to what was reported in Ref. [117] for the DNLS model, in direct simulations unstable vortex crosses evolve into onsite-centered fundamental discrete breathers (with S=0) by transferring almost all the energy to one of the sites which originally formed the cross. The corresponding instability border (for C=1) in the (µ,ωb) plane is depicted by the right curve of figure 4.12.b. More interesting is the case of µ < 0. In this regime, we have found that fundamental vortex crosses experience another bifurcation, with a quadruplet of Floquet eigenvalues leaving the unit circle at λ= +1 (the so-called harmonic bifurcation). With the decrease of µ, the corresponding two pairs of the eigenvalues move along the real axis in the opposite direction, until each pair breaks up, as shown in figure 4.12.d. The unstable eigenvectors, 82 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES (a) 10-6 10-5 10-4 10-3 10-2 10-1 100 0 2 4 6 8 10 n 0 2 4 6 8 10 m 0.5 0.4 0.3 0.2 0.1 0 |δΦ*n,m| (b) 10-6 10-5 10-4 10-3 10-2 10-1 100 0 2 4 6 8 10 n 0 2 4 6 8 10 m 0.5 0.4 0.3 0.2 0.1 0 |δΦ**n,m| (c) 7 6 5 4 3 80 60 40 20 0 |Φn,m|2 time (n0,m0+1) (n0,m0-1) (n0+1,m0) (n0-1,m0) Figure 4.13: (a) and (b) Intensity profiles of the unstable Floquet eigenvectors, δΦ∗δΦ∗∗, corresponding to the bifurcation at µ=−0.3 (for C=1) of the vortex cross with ωb=8, see figure 4.11.d. (c) Time evolution of the lattice field at sites around the center of the same unstable vortex solution. Pulsonic dynamics of the amplitudes is observed, without decay of the vortex pattern. 4.3. DISCRETE VORTEX BREATHERS 83 δΦ∗and δΦ∗∗, associated with this bifurcation are plotted in figure 4.13.a and 4.13.b (in this notation, ∗does not stand for complex conjugation). The shape of each eigenvector reveals strong localization at two opposite sites of the vortex cross, each one separately breaking the spatial symmetry (2D isotropy) of the original solution. Adding a small perturbation to the solution along one unstable direction causes oscillations of the amplitudes around the vortex center, as shown in figure 4.13.c. Such behaviour persists at longer times; in fact, the vortex pattern does not disappear but rather suffers irregular modulations of its local amplitudes. This picture of the instability development supplements the stability diagram for the fundamental vortex crosses, which is displayed in figure 4.12.b in the (µ,ωb) plane (as noted above, for the isotropic model, with C=1). Note that the border of the instability which transforms the vortex cross into its oscillatory counterpart (the left curve in the figure) stays in the µ < 0 region, even for large frequencies. Therefore, unlike the HHB described above, this instability is dominated by the competition between the selfdefocusing inter-site and self-focusing on-site nonlinearities. A further insight into the nature of this bifurcation is provided by the observation that it coincides exactly with the divergence of norm N(ωb, µ) of the discrete breather (and of the vortex cross solution), and thus it coincides with the stability interchange between in-phase and out-of-phase bound states analyzed above in section 4.2.2. Regarding the vortex cross as made up of two (perpendicular) out-of-phase bound states of breathers (say, left-right and top-bottom), one would be tempted to interpret the quadruplet of eigenvalues leaving the unit circle at +1 as the two pairs of eigenvalues that signal the simultaneous instability of both out-of-phase bound states. At least, this interpretation would explain the fact that a quadruplet of eigenvalues simultaneously leave the unit circle at +1, and it is fully consistent with the shape of the Floquet eigenvectors in figure 4.13. This interpretation suggests that the bifurcation of vortex crosses occurring in the left part of figure 4.12.b is the same one experienced by out-of-phase bound pairs of breathers in figure 4.9 (for separation d=1). In any case, a noteworthy numerical finding is that these bifurcations (of bound states and vortex crosses) not only coincide but are also characterized by the divergence of the breather norm. 4.3.2 Vortex squares We have also studied the smallest (M=1) vortex squares carrying S=1vorticity. Forthis purpose, we have performed the continuation of the corresponding solution family, starting from a configuration with nonzero components Φn0,m0=−iΦn0,m0+1=−Φn0+1,m0+1= iΦn0+1,m0=√ωb/2 in the anticontinuum limit, eq. (4.18). As in the case of the vortex cross, we have first performed the continuation in the coupling constant Cto obtain the corresponding solutions for the DNLS model (C=1, µ=0). Again, for low-frequency vortex squares, we have observed an HHB at some critical value of C. For high-frequency solutions, a bifurcation of the same type is observed when the continuation is performed from the DNLS model to values µ > 0. In figure 4.14.a, one can observe this bifurcation for the vortex square with ωb=8. The corresponding HHB (see figure 4.14.c) occurs with a quadruplet of the Floquet eigenvalues leaving the unit circle. The behavior of the unstable solution is the same as for the vortex cross, and, after a transient, a regular breather with S=0 emerges at one of corner sites of the former vortex square, while the field at 84 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES 8 7 6 5 2 1 0-1-2 µ 8 7 6 5 2 1 0-1-2 µ 10 7.5 5 2.5 0 Unstable b µ Unstable ω |λ| Stable ωb=8.0 2 1 0 -1 -2 2 1 0-1-2 Im(λ) Re(λ) 2 1 0 -1 -2 Im(λ) (a) (b) (c) (d) µ=0.14 µ=−0.04 Figure 4.14: (a) The absolute value of the Floquet multipliers as a function of µfor a vortex square of minimum size (M=1) with ωb=8 (C=1). Two bifurcations are revealed by Floquet distributions in panels (c) and (d). At µ=0.14, we find a Hamiltonian Hopf bifurcation, whereas at µ=−0.04 a quadruplet of eigenvalues leave the unit circle and start a short trip to +1, from where they leave the unit circle again. The entire stability diagram is represented in panel (b), showing a narrow stability region. three other corners nearly vanishes (i.e. the energy mainly concentrates at a single site of the initial vortex structure). With the continuation of the vortex square to µ < 0, we have again (as in the case of vortex crosses) found that the solutions suffer a destabilizing bifurcation different from that at µ > 0. However, the bifurcation for µ < 0 (see figure 4.14.d) is also different from its counterpart for the vortex cross (which was displayed above in figure 4.12.d). At some value µ < 0, a quadruplet of Floquet multipliers leave the unit circle, to return to it at +1. After this brief excursion, they immediately leave the unit circle again, and instability grows with |µ|. Unlike its counterpart for the vortex cross, this bifurcation does not correspond to the interchange of stability for the bound state of breathers analyzed in 4.2.2, which actually occurs at a lower value of µ, where the vortex square is already unstable. However, it is remarkable that precisely at this value of µthe quadruplet of eigenvalues outside the unit circle meet instantaneously at +1, so that the vortex square is marginally stable at that point. Profiles of unstable eigenvectors, δΦ∗and δΦ∗∗, are shown in figure 4.15.a and 4.15.b. Each one is localized at two non-adjacent corners of the plaquette where the vortex square is located. The dynamics triggered by the original solution being perturbed by this δΦ∗(or equivalently δΦ∗∗) is displayed in figure 4.15.c. Again (as in the case of the vortex cross), the vortex pattern is not destroyed (in contrast with the unstable behavior at µ > 0). 4.3. DISCRETE VORTEX BREATHERS 85 (a) 10-6 10-5 10-4 10-3 10-2 10-1 100 0 1 2 3 4 5 6 7 n 0 1 2 3 4 5 6 7 m 0.5 0.4 0.3 0.2 0.1 0 |δΦ*n,m| (b) 10-6 10-5 10-4 10-3 10-2 10-1 100 0 1 2 3 4 5 6 7 n 0 1 2 3 4 5 6 7 m 0.5 0.4 0.3 0.2 0.1 0 |δΦ**n,m| (c) 3.9 3.8975 3.895 3.8925 3.89 60 80 100 120 140 160 180 200 220 |Φn,m|2 time (n0,m0) (n0-1,m0+1) (n0,m0+1) (n0-1,m0) 0 0+1 (n ,m ) (n ,m ) 0+1 (n ,m ) (n ,m ) 0+1 0 00 0+1 Figure 4.15: (a) and (b) Intensity profiles of the unstable Floquet eigenvectors, δΦ∗and δΦ∗∗, corresponding to the bifurcation, at µ=−0.04 (for C=1), of the vortex square with ωb=8, shown in figure 4.13.d. (c) The time evolution of the lattice field at the vortex-square’s corners for the same unstable solution. The simulations reveal periodic evolution of the amplitudes with a clear sequence of energy transfer between the adjacent sites following the same pattern as the current flux in the original vortex solution. 86 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES Instead, the lattice field at the vortex-square sites develops a periodic pulsonic behavior, in which at least two frequencies can be identified. One of the frequencies accounts for periodic transfer of energy between four corners of the square vortex, following the same path as the flux current: (n0,m0)→(n0,m0+1) →(n0+1,m0+1) →(n0+1,m0)→(n0,m0)→... (4.19) Another noteworthy feature of the dynamics in this case is that the total amount of energy that is periodically transferred between neighboring sites varies, also in a regular periodic fashion, thus giving rise to the second frequency. Again (as happened for the vortex cross), the instability observed at µ < 0 induces a pulsonic dynamics of the lattice amplitudes but, in the present case, the dynamics is much more regular. An intriguing numerical observation is that the value of µat which the quadruplet of eigenvalues meet at +1 (so that the vortex-square solution momentarily becomes marginally stable) occurs exactly when the breather norm diverges. The entire stability diagram for the fundamental vortex squares is presented in figure 4.14.b. Again, we find a narrow stability region for low-frequency vortex squares that expands as the frequency increases. 4.3.3 Bound states of discrete vortex crosses As a first step towards the characterization of the stability of more complex 2D arrangements of vortices, we have studied two types of bound states of vortex crosses, with the vortex centers aligned along a lattice axis (say, the x-direction). In the two types of the bound state, the vortices have equal or opposite vorticities, see figures 4.16.a and 4.16.b. Both types of solutions were studied on the isotropic Salerno lattice with competing nonlinearities (C=ν=1 and µ < 0), and were numerically obtained by the continuation at µ=0 from the anticontinuum limit (C=0), followed by the a second continuation in the direction of negative inter-site nonlinearity µ. The Floquet spectrum of the solution was also numerically computed along the continuation path. At µ=0, bound states of vortices with equal vorticities are stable, while those with opposite vorticities are unstable. To explain this numerical observation, one has to realize that the right-most member of the breather set forming the left vortex, and its left-most counterpart in the right vortex are out-of-phase (in-phase) in the former (latter) case, see figures 4.16.a and figures 4.16.b. Then, the stability analysis of bound states of breathers reported above in section 4.2.2 suggests that the stability of the bound states of vortices is actually dominated by the stability of the local bound state of the two constituent breathers (one from each vortex) that are in the closest proximity. This analysis is further validated by comparison of unstable Floquet eigenvalues for the bound state of vortices with opposite vorticities and those for the bound state of in-phase breathers (for the corresponding values of the frequency and separation between the centers). Whenµdecreases, a destabilizingbifurcation occurs, as expected, in theequal-vorticity bound state, precisely at the same value of µwhere the simultaneous instability of the vortex cross (in section 4.3.1) and the out-of-phase bound state of ordinary breathers occurs. By inspection of the Floquet spectrum for the bound state of vortices, one can clearly identify pairs of eigenvalues associated with each of these instabilities that take place simultaneously at this bifurcation point. It clear that the stability of bound states of discrete 4.3. DISCRETE VORTEX BREATHERS 87 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 -1 -0.8 -0.6 -0.4 -0.2 0 |λ| µ out-of-phase in-phase (c) (b) (a) Equal Vorticities Opposite Vorticities Figure 4.16: A schematic representation of the in-phase (a) and out-of-phase (b) bound states of vortices with S=1, in the limit of C=0, ν=1, µ=0. Vectors stand for the instantaneous values of Φn,min the complex plane, with Φn,m=√ωb/2. These solutions are continued in C up to C=1, and then continued in µ. Panel (c) shows the evolution of the Floquet multipliers as a function of µwhen µ < 0. The results correspond to ωb=8 and the distance between the two vortex centers is set to be d=5 (as seen in (a) and (b)). vortex and that of single vortices in the SM with µ < 0 is related to the behaviour found for bound states of two pinned breather solutions. The decomposition of any complex solution in terms of this latter building blocks is clearly of importance. 94 CHAPTER 4. DISCRETE BREATHERS IN 2D NONLINEAR SCHR ¨ ODINGER LATTICES 4.5 Conclusions and Prospective Remarks We have studied here the dynamics of exact numerical discrete breathers (pinned, vortical and mobile ones) in a two-dimensional anisotropic nonlinear Schr¨ odinger lattices. These solutions are computed from a set of uncoupled 1D chains into increasing non-zero values of the coupling in the transversal direction in order to reach the 2D limit. It is convenient review the most salient results in order to have a compact picture of the 2D behavior of discrete breathers. Pinned breathers.- We have performed an extensive exploration in the parameter space (ωb,C2,ν) of breather frequency, transversal coupling and Salerno parameter, by computing the Floquet spectra of the numerical solutions. Both the 1D solutions of the standard and the competing SM have been continued into the 2D regime in order to see the effects of the dimensional increment. In particular we have found the link between the unstable behaviour found for certain breather frequencies in the 1D competing Salerno model (whose pulsonic character resembled those of the well known 2D unstable solutions) and the quasi-collapse instability that appears for low frequency breathers when the coupling in the transverse direction is incorporated. Furthermore, we have analyzed the dynamics on the quasi-collapse unstable manifold, where the unstable breather experiences a shift in frequency towards the (higher) value of the stable breather with the same norm. The excess of energy is coherently transferred to oscillations of the breather width, so that the resulting pulson state is characterized by two frequencies. We have also recovered the 2D counterpart of the 1D cuspons and peakons for the 2D SM with competing nonlinearities. Again these hyperlocalized states are stable Finally, the stability analysis of in-phase and out-of-phase bound states of breathers in the isotropic lattice reveals that there is a stability interchange between both types of bound states, precisely at the same value of the intersite-nonlinearity parameter (µ) where the breather norm diverges as happened for the 1D model. Vortex breathers.- In addition to fundamental breathers, discrete vortices of two types, crossand square-shaped ones, have also been constructed, and their stability regions identified. In direct simulations, unstable vortices in the standard 2D Salerno model of the ordinary type transform into regular breathers, while in the model with the competing nonlinearities the instability turns vortices into localized vortical pulsons, without destroying their topological character. It is then worth mentioning the ubiquity of this pulsonic attractors of the dynamics in the model. Regarding the stability of bound states of vortex crosses, we have shown that it is determined by the stability of the local bound state of two constituent breathers (forming the two vortices) which are in the closest proximity. Mobile breathers.- We have studied discrete breathers moving along the strong coupling direction for the standard 2D SM. These solutions are composed of an exponentially localized core on top of an extended background which is itself the finite sum of a finite set of nonlinear 2D plane waves. The time scales associated with these plane waves are resonant with the core internal frequency as happens in the 1D case. In particular, the background chooses a finite set of plane waves from a continuous family of resonant solutions. The Floquet analysis of these mobile discrete breathers reveals the existence of two distinct types of instability. One is the counterpart, for mobile breathers, of the quasi-collapse experienced by pinned breathers. The other instability occurs in a region of parameter space where pinned breathers are linearly stable. The analysis of the dynam- 4.5. CONCLUSIONS AND PROSPECTIVE REMARKS 95 ics on the unstable manifold show that the excess of energy is partly transferred to a small moving pulse, ejected from the center of localization, which justifies the designation of a fission instability. However, part of the energy excess is also transferred to width oscillations. The appearance of pulson states far from the quasi-collapse regime indicates that the tendency to allocate energy in the form of width oscillations is a general 2D feature, not exclusively associated to quasi-collapse instabilities. We leave the question on mobility of 2D discrete breathers in an arbitrary lattice direction. The results obtained here shed light about how this mobility can be obtained. In fact, our experiences show that mobility of pinned breathers can be induced based on the existence of the extended background in the numerically exact mobile solution. On the other hand, the results obtained here and the aforementioned future work may help to design and better understand recent numerical experiments reported in [136], concerning the interaction between high amplitude pinned breathers and mobile ones. These experiments provides a possible way for routing and blocking mobile discrete breathers via the interaction with the high amplitude pinned ones, resulting in a plausible implementation of logical functions. Part II Structure and Dynamics of Complex Networks 97 p Presentation of Part II The second part of the Thesis is devoted to the study of the structure of complex networks. Traditionally, physics has focused on systems where the underlying topology of elements’ interactions is described by regular lattices such as those studied in the preceeding part. However, in the recent years, physicists have started to look to those systems where the interactions among constituents reflect the abstract relations between pairs of elements rather than being determined by the proximity in a physical space. These relations can be determined by the existence of monetary transactions between banks in economic networks, or cooperative and friendship relations between individuals in social networks, or assemblies of different molecules working together to develop cellular tasks in biological networks, etc... From the highest to the lowest level of description we find complex networks of interactive elements that cannot be described by regular patterns of connections. The growing interest in the characterization of the above systems has led to the emergence of the so-called network science [137]. Let us review the development of this new interdisciplinary field. One can settle the first steps of network science with the works on graph theory [138, 139] in the middle of last century. The most remarkable result is the theoretical analysis of a random network by the mathematicians Paul Erd¨ os and Alfr´ ed R´ enyi [140, 141]. However, networks where interactions among elements are completely random are a coarse-grained approach to real networked systems, assuming a homogeneous disorder in what concerns the patterns of connections. The burst in the study of complex networks came with the advent of the XXI century along with the development of the Internet and the World Wide Web. This development has provided a large amount of data-sets for unveiling the relations established among industrial companies, institutions, scientists, etc... Besides, the explosion of human mobility (provided by the increase of accessible infrastructures and transportation companies) and the boom of new telecommunications tools (mobile phones, instant messaging services, etc...) has hugely facilitated the stablishment of new agent networks with a high global character. These two ingredients, the emergence of new networked systems and the high accessibility to data-sets describing them, constituted an unprecedented opportunity for scientist to analyse their topological features. The analysis of real complex networks revealed that seemingly different systems share a common property when looking to the distribution of the number of connections that the elements of the networks have. It is found [142–154] that most of the networks present a power law functional shape for this statistical quantity and, therefore, they differ from that accounted by Paul Erd¨ os and Alfr´ ed R´ enyi, where all the elements present a similar number of neighbours. Besides the surprising fact that real networks (accounting for many 99 100 different types of interactions) share the scale-free character, understanding the (common or not) origin of this internal organization has become a challenging question for many researchers. The above astonishing findings lead physicists to construct simple models of network growth in order to reproduce the “universal” properties found for real networks. In this sense, the models developed by Duncan J. Watts and Steven H. Strogatz [155], and AlbertL´ aszl´ o Barab´ asi and R´ eka Albert [143] deserve special mention. Another direction of research has been focused by the search for statistical measures of network topology in order to handle efficiently the large amount of available data-sets and characterize those networks they represent with a few meaningful indicators. The purpose of these two types of studies differ strongly from those of concern of traditional graph theory (where rigorous theorems of difficult real application are proved) and, at the same time, are methodologically far from the meticulous system characterization performed by biologists (who tend to overpay attention to single element details to catalogue systems so that a unitary analysis of different systems becomes difficult). The statistical point of view and the unitary approach to the problem of network characterization performed by physicists have clearly taken advantage over other disciplines under the name of statistical physics of complex networks2. In chapter 5 we will briefly review the most important tools for characterizing network structure and present two model of synthetic network generation. The next step of network science has been to look to network dynamics. Most networks are not composed of mere static objects but, on the contrary, their elements develop a function. This function can be as simple as being routers for the transfer of entities among their elements, or as complicated as being regulatory agents of some internal dynamical processes performed at each network node. It is important to difference two kind of studies to the problem of network dynamics. On one hand, given that real networks can be described by a set of statistical measures and that one can stablish subsets of networks which are qualitatively similar in terms of these quantities, it is therefore interesting to find how to implement dynamical processes on top of the network in order to take advantage of these topological characteristics. This search for efficient algorithms is not only motivated for practical purposes but it is also interesting for studying the interplay between dynamics and the underlying networked substrate. In this sense, the models developed for constructing synthetic networks are a useful benchmark for studying this interplay before applying the results to real networks. This first kind of studies are therefore interesting for networks whose dynamics can be modeled or modified. This is the case of technological, logistic and certain social networks, e.g. information and transportation networks are susceptible targets of this kind of studies. In chapter 6 we will deal with two related problems, namely the interplay between network structure and, first, the performance of immunization strategies aimed at stopping the spread of epidemics and, second, the routing policies for information dynamics. As we mentioned above there is a second kind of studies on network dynamics. Instead of varying the dynamical properties assuming a fixed substrate in this second class of studies both structure and dynamics are asummed to co-evolve towards a stable state where system’s performance is optimal. We will leave the discussions on this interesting problem for part III and now let us focus on the study of the structure of networks and its influence over the dynamical performance. 2Interesting reviews and tutorials on the subject are found in [156–163]. Chapter 5 Network Structure and Generation There may well be no useful parallel to be drawn between the way in which complexity appears in the simplest cases of many-body theory and chemistry and the way it appears in the truly complex cultural and biological ones, except perhaps to say that, in general, the relationship between the system and its parts is intelectually a one-way street. — Philip W. Anderson in More is Different [2]. This chapter is devoted to the description of the structure of networks and the modelling of their growth and evolution. We have seen in the preceding pages that a great amount of empirical data about the patterns of connections among the constituents of social, technological, logistic and biological networks is nowadays available. At this point several questions arise such as How similar real networks are? or Is there any common feature (regularity) between networks with the same function? In order to answer (if possible) these kind of questions we have to define some properties that would allow us to give a quantitative and qualitative description of the architecture of networks. After defining these magnitudes we will briefly describe some important models of synthetic networks that try to capture some of the ingredients observed when analysing the “native” ones. We will round offthe chapter with a deep analysis of two models of network design that will be employed along the forthcoming parts of this Thesis. 5.1 Describing Complex Networks Before defining the magnitudes employed for describing the networks we give a brief account of some formal definitions and notations inherited from classical graph theory. Afterwards, we list and explain the most used quantities for characterizing networks’ structure at local, global and mesoscopic levels. 102 CHAPTER 5. NETWORK STRUCTURE AND GENERATION 5.1.1 Basic definitions We start with the formal definition of a network (or, in mathematicals terms, a graph) G(V,E) as an ordered pair of set of sets: a non null set Vof elements called nodes (or vertices) and another set Eof pairs (i,j), with i,j∈V, called links (or edges, arcs) that denote that nodes iand jare connected. Normally one imposes that i,jso that selfconnections are avoided. We will denote by Nand Lthe cardinalities of the sets Vand Erespectively. Along with this definition we can also consider that the elements in Eare ordered pairs (i,j),(j,i). In this case we will talk of a directed network (or digraph). It is also very common to assign weigths (numbers) to the edges so that we have a weighted (or valued) network. The cardinality of Vand Ecan tell us about the nature of the graph. Taking into account that the maximal cardinality of Eis N 2we will talk about a sparse network when L ≪ N2and a dense one when L ∼ N2. Asubgraph G′(V′,E′) of G(V,E) is a graph such that V′⊆Vand E′⊆E. A subgraph is said to be maximal with respect to a given property if it cannot be extended by adding elements either to V′or E′without loosing its property. Aditionally, we say that a subgraph G′is induced by Gif E′contains all the pairs (i,j)∈Ewith i,j∈V′. The complement of a graph G(V,E) is the graph G′(V′,E′) (sometimes denoted G) so that V=V′but whose edge set E′consists of the edges not present in E. Then a graph G(V,E∪E′) will be a complete graph, i.e. every node will be connected with the rest of the N−1 nodes. In order to manage with a graph one usually label with natural numbers the elements of Vso that i=1, ..., N. Different assignation of labels to the elements of Vyield isomorphic networks and the topological properties are not affected. Formally, two graphs G1(V1,E1) and G2(V2,E2) are said isomorphic when one can stablish a bijective relation φ:V1→V2 that preserve the connections, i.e. if (i,j)∈E1then (φ(i), φ(j)) ∈E2. One can represent the graph with the so-called adjacency matrix Awhose elements are aij =1 if (i,j)∈Eand aij =0 otherwise. This matrix will be symmetric for undirected graphs but, in general, this is not the case for digraphs. In the case of weighted networks one can replace the non zero elements of Aby the weights of the corresponding links in order to obtain a complete representation of the graph. The analysis of the adjacency matrix will give the topological characterization of the networks. At present, there is a large amount of different magnitudes used for characterizing the networks architecture. However, the more specific field we study the more properties we will find for describing these complex topologies. Then, we will only emphasize on those quantities which are of general use and therefore will be used along the works described in this Thesis. We will divide the magnitudes depending on the scale involved for their definition. From our point of view this is a useful definition since local (microscopic) or global (macroscopic) properties play a key role depending on the kind of dynamics placed on top of the network. 5.1.2 Single nodes properties Local Magnitudes We will refer to local quantities when one takes into account a node iand its neighbours, Γi. Obviously, the first local property is the degree of a node i,ki, which is the cardinal 5.1. DESCRIBING COMPLEX NETWORKS 103 of the set Γi,i.e. the number of nodes which iis linked to or, in terms of the adjacency matrix ki= N X j=1 aij .(5.1) If one is considering a directed network one will talk about the “in-degre” of a node i, kin i, and its “out degree”, kout i, which are the number of incoming and outgoing links that a node shares with its neighbours. Again, we can obtain such quantities from the adjacency matrix by kin i= N X j=1 aij and kout i= N X i=1 aij .(5.2) Another interesting local measure is the so-called clustering coefficient of a node which measures the number of connections among the neighbours of a node, ei. This quantity is usually normalized to one by dividing by its maximum value ki 2so that it measures the probability that two neighbours jand mof a node i(aij =aim =1) are also linked to each other (ajm =1). A formal expression of the clustering of a node iwith ki neighbours is ci=2PN j,m=1aijaimajm ki(ki−1) .(5.3) A third local property arise when looking at the degree of the neighbours of a given node (here we assumme that this information is on the local horizon of the nodes) so that we can define the average nearest neighbours degree. We can write this quantity as knni=PN j=1aij PN m=1ajm PN j=1aij .(5.4) Global Magnitudes Now we will define two properties that are defined taken into account pairs of nodes that are not necessarily linked and are then influenced by the topology of the whole graph. These measures are the closeness centrality and the betweeness centrality. As we will see in section 6.2 these two magnitudes will play a key role when dealing with problems of propagation through networks. First of all we define the distance between two elements of the network, dij, as the length of the geodesic that goes from node ito j. In principle, one can observe more 00 00 00 11 11 11 00 00 00 11 11 11 00 00 00 11 11 11 00 00 00 11 11 11 00 00 00 11 11 11 000 000 000 111 111 111 00 00 00 11 11 11 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 000 000 000 111 111 111 00000 00000 00000 00000 11111 11111 11111 11111 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0000000 0000000 0000000 0000000 0000000 1111111 1111111 1111111 1111111 1111111 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 00 00 00 11 11 11 0 0 0 1 1 1 0 0 0 1 1 1 00 00 00 00 11 11 11 11 0011 0011 00 00 00 11 11 11 0 0 0 0 1 1 1 1 00 00 11 11 0 0 0 1 1 1 00 00 00 00 00 11 11 11 11 11 0 0 0 1 1 1 00 00 11 11 00 00 00 00 00 00 11 11 11 11 11 11 0 0 1 1 0 0 0 1 1 1 0 0 1 1 0 0 0 1 1 1 0011 0 0 0 1 1 1 00 00 00 11 11 11 00 00 00 11 11 11 00 00 00 11 11 11 00 00 00 11 11 11 00 00 00 11 11 11 000 000 000 111 111 111 00 00 00 11 11 11 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 000 000 000 111 111 111 00000 00000 00000 00000 11111 11111 11111 11111 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0000000 0000000 0000000 0000000 0000000 1111111 1111111 1111111 1111111 1111111 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 00 00 00 11 11 11 00 00 00 00 11 11 11 11 0 0 0 0 1 1 1 1 00 00 11 11 0 0 0 1 1 1 00 00 00 00 00 11 11 11 11 11 0 0 0 1 1 1 00 00 00 00 00 00 11 11 11 11 11 11 0 0 0 1 1 1 0 0 0 1 1 1 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 00000 00000 00000 00000 00000 00000 11111 11111 11111 11111 11111 11111 000000 000000 000000 000000 000000 000000 000000 111111 111111 111111 111111 111111 111111 111111 00 00 00 00 00 00 00 00 11 11 11 11 11 11 11 11 00000 00000 11111 11111 Figure 5.1: Two different situation in terms of the clustering of the striped node. In the first case (left) clustering is 0 while for the second example (right) the probability of finding two connected neighbours raises to 2/7. 110 CHAPTER 5. NETWORK STRUCTURE AND GENERATION (a) (b) Figure 5.5: (a) Schematic picture of a set of 4 communities (surrounded by dashed circles). The density of inner-links between nodes of the same community is much larger than that of the links with members of the rest of the network. (b) The 7 possible 4-nodes motifs. graph which is a subgraph of G. An example of all the possible 4-node connected graphs is illustrated in figure 5.5.b. The concept of motifs (originally introduced by Uri Alon and coworkers [173–177]) was employed for studying the finding of recurrent patterns of interconnections between a small number of nodes in biological and other networks. In order to obtain a quantitative description for the appearance of the significant motifs in a graph G, one makes use of matching algorithms for counting the total number of occurrences of each n-node subgraph Min the original graph and in the randomized ones. Then, one can define the statistical significance of a given motif Mby some score function, like the so-called Z-score [174, 176] ZM=nM−hnrand Mi σrand nM ,(5.15) wherenMand nrand Marethe number of times the subgraph Mappears in Gand inits randomized counterpart repectively. σrand nMis the standard deviation of the number of appearances in the randomized network ensemble. 5.2 Overview of network generation models In this section we briefly account for several important network models. Nowadays there is a huge number of ways for generating complex networks that try to capture the properties of real graphs. Many of them are variations of the models we present below since they represent seminal works on the matter. For a complete review on currents trends in network modeling we refer the reader to [156, 157, 161, 162]. 5.2. OVERVIEW OF NETWORK GENERATION MODELS 111 5.2.1 Random graphs We call random graphs to those network where the links between nodes are randomly distributed3. In their seminal work in the subject, Erd¨ os and R´ enyi [140] (ER for short) proposed a method for the construction of random graphs with Nnodes and Llinks: Starting from Nisolated nodes, pairs of randomly chosen nodes are linked avoiding self and multiple connections. This process is stoped when Llinks have been stablished. A single graph obtained using the above recipe is one of the  N 2! Lpossible equiprobable realizations. The set of all these possible realizations is called the set of uniform random graphs with Nnodes and Llinks, GER N,L. In a random graph the probability that two given nodes are linked is L/N 2. Another possible strategy for constructing random graphs is to sample every pair of nodes and with probability 0 <p<1 link them. This procedure defines a different set called random binomial graphs,GER N,p, that contains graphs with different number of total links Lbeing pL(1 −p)N 2−L (5.16) the probability that a graph belonging to GER N,phas Llinks. Then the average number of links of a graphs in this set is pN 2. The two sets (uniform and binomial) of random graphs are tightly related to the canonical and grand canonical ensembles of the equilibrium statistical mechanics when one looks at the number of edges as the number of particles in the system. Both ensembles converge to the same set in the thermodynamic limit N→ ∞ when approached keeping hkifixed (wich is equivalent to fix 2L/Nand p(N−1) in the uniform and the binomial sets repectively). The structural properties of the ER graphs vary as a function of p, showing a dramatic change at the critical probability pc=1/Nthat corresponds to hki=1. In particular: •If p<pc, the size of the giant component of the graph is of the order O(ln N) graph and there is no graph component with more than one loop. •If p=pc, the size of the giant component goes with O(N2/3). •If p>pc, the graph has a giant component with a number of loops that scales as O(N) and there is no other graph component with more than O(lnN) elements neither with more than one loop. This transition (characterized by Erd¨ os and R´ enyi in [141]) is strongly related with the percolation transition studied in the theory of critical phenomena [178]. In ER graphs the probability that a node has kneighbours follows the binomial distribution P(k)= N−1 k!pk(1 −p)N−1−k,(5.17) 3In fact, all the network models analysed in this chapter are strictly random in the sense of the mechanism adopted for their construction. However, the term random graph is overused in the literature for calling Erd¨ os-R´ enyi networks. 112 CHAPTER 5. NETWORK STRUCTURE AND GENERATION 0 1 0.05 p Figure 5.6: Three kind of networks obtained used the Watts-Strogatz method starting from a regular one-diemnsional network where every node is linked to its 4 nearest neighbours. For p=0 we obtain the regular network. For small values of p≪1 we have “Small World” networks from the small amount of reassigned links. Finally, for p=1 random networks are obtained. that for hkifixed and N→ ∞ tends to the Poisson distribution P(k)=hkik k!exp(−hki).(5.18) Erd¨ os and R´ enyi graphs are uncorrelated since the links are launched at random independently of the degree of the nodes. As a consequence, P(k′|k) and knn(k) are independent of k. Concerning to the connectivity properties of an ER graph, when p>lnN/Nnearly all the generated graphs are composed of one single component and the average path length takes values around hki=lnN/ln(pN)=ln N/hkibecause locally the ER topology is viewed as a tree like structure where a single node has hkineighbours, hki2nodes at distance 2,... Finally, since pis the probability of two nodes sharing a link there will be pk(k−1)/2 links among the neighbours of a node of degree kso that the clustering coefficient goes as c=p=hki/N, and then it vanishes in the thermodynamic limit. 5.2.2 Small-world networks In 1998 Watts and Strogats (WS) proposed a method of graph constrution that allows to obtain networks with a high clustering coefficient andsmall average path leghth [155]. We have seen that ER graphs present a small value of Land that cvanishes when N→ ∞. On the other hand, a regular network with connections to first, second, third,... next nearest neighbours presents a high value of the clustering coefficient joined with large values of L. In some sense, the WS model interpolates smoothly between these two topologies. The WS procedure starts from a ring (see figure 5.6) where every node is symmetrically linked to its 2mnext nearest nodes so that there is L=mN links. Then, every link is considered and with probability pit is substituted by another link that connects one of the original nodes with a new one chosen at random. Note that for p=0 we maintain the original regular topology whereas for p=1 an ER random graph is generated. In figure 5.7 it is shown how for a range of pvalues the WS model generates networks with 5.2. OVERVIEW OF NETWORK GENERATION MODELS 113 Figure 5.7: Evolucion of the clustering coefficient and the average shortest path length as a function of p. Note that near p=4·10−3the clustering remains comparable to the values of the regular network whereas Lhas decreased significatively. both the small world property (due to shortcuts added when p,0) and high clustering coefficient (inherited from the regular topology), two characteristics shared by a number of real networks. This result reveals that the clustering coefficient is very robust under link reasignation whereas Lrapidly decreases when a few shortcuts are incorporated. Analytical calculations on the transition observed in the WS model are found in [179– 182]. It has been shown that the appearance of the small world character as pincreases is not a phase transition but a crossover phenomenon. The characteristic length satisfy the scaling relation L(N,p)=N f(Np) where f(x)∼(cif x≪1 ln x xif x≫1(5.19) Besides, in [181] the authors found the analytical expressions for the clustering and the degree distribution as a function of the control parameter p c(p)=3(m−1) 2(2m−1)(1 −p)3(5.20) P(k;p)= min(k−m,m) X i=0 m i!(1 −p)ipm−i(pm)k−m−i (k−m−i)! exp(−pm),(5.21) the last equation (5.21) is valid provided k≥motherwise P(k<m;p)=0. The WS model was later modified by Newman and Watts in order to solve the possible formation of disconnected graphs of the network as shortcuts were incorporated. Then, they proposed to add new links between randomly chosen nodes instead of making the rewiring process [183]. They considered every node and with probaibility pa link was stablished with any other node of the network so that the average number of shortcuts added is pN. 5.2.3 Scale-Free networks There are a large number of models that reproduce the power law functional form for the degree distribution. However, we will focus here on those models that incorporates the growing character present in real networks, where the amount of nodes grows with time, to the formulation of the model. These models usually consist of an initial small subset of nodes to which new nodes are sequentially incorporated by launching new links over 114 CHAPTER 5. NETWORK STRUCTURE AND GENERATION 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 000 000 000 000 000 111 111 111 111 111 000 000 000 000 000 111 111 111 111 111 000 000 000 000 111 111 111 111 000 000 000 000 111 111 111 111 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 11 11 11 11 00 00 00 00 00 00 00 11 11 11 11 11 11 11 0000 0000 0000 1111 1111 1111 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 00 00 00 00 00 00 00 00 11 11 11 11 11 11 11 11 0000 0000 0000 1111 1111 1111 00000 00000 00000 00000 00000 00000 00000 00000 00000 11111 11111 11111 11111 11111 11111 11111 11111 11111 000111 00 00 11 11 00 00 00 00 00 00 00 00 00 00 11 11 11 11 11 11 11 11 11 11 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 00000 00000 11111 11111 0000000 0000000 0000000 0000000 0000000 0000000 1111111 1111111 1111111 1111111 1111111 1111111 000 000 000 000 000 000 000 000 111 111 111 111 111 111 111 111 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 000 000 000 000 000 000 000 000 000 000 000 000 000 111 111 111 111 111 111 111 111 111 111 111 111 111 000 000 000 000 000 000 000 000 000 000 000 000 000 111 111 111 111 111 111 111 111 111 111 111 111 111 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 00000000 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 11 000000 000000 000000 000000 000000 111111 111111 111111 111111 111111 000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000 000000000 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 111111111 00000000 00000000 00000000 00000000 00000000 11111111 11111111 11111111 11111111 11111111 00000000000 00000000000 11111111111 11111111111 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 000 000 000 000 000 000 000 000 000 000 000 000 000 000 111 111 111 111 111 111 111 111 111 111 111 111 111 111 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 00000 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 11111 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 000000 000000 000000 000000 000000 000000 000000 000000 000000 000000 000000 000000 111111 111111 111111 111111 111111 111111 111111 111111 111111 111111 111111 111111 000000000 000000000 000000000 111111111 111111111 111111111 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 1111111 00 00 00 00 00 00 00 11 11 11 11 11 11 11 0 0 0 0 0 0 0 1 1 1 1 1 1 1 000000 000000 000000 111111 111111 111111 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 0000 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 1111 000 000 000 000 000 000 000 000 000 000 111 111 111 111 111 111 111 111 111 111 00 00 00 00 00 00 00 00 11 11 11 11 11 11 11 11 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 0000000000 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 1111111111 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 Figure 5.8: Schematic representation of the grwoth process. At each time steps a new node is incorporated to the network core linking to m=3 nodes that already belong to the core. those nodes that already take part of the network (see figure 5.8). In particular, the work by Barab´ asi and Albert (BA) in 1999 [143] suposed an important breakthrough to the problem of finding the roots of the SF behaviour of real networks and had the growing process as a key ingredient of their formulation. The BA models works starting from an initial core of m0isolated nodes. At each time step a new node is incorporated to the network by launching m≤m0links over the already existing ones so that the network core grows linearly in time. The probability that one of the core nodes, i, receives a link from the new node is proportional to its degree, ki, Πpa i=ki PN(t) j=1kj ,(5.22) where N(t) is the number of nodes that form the network core at time t,N(t)=m0+t−1. Besides, the total number of links at time tevolves as L(t)=mt. The above rule for node selection was termed preferential attachment and favours that a node with more links than others will increase its connectivity at a higher rate (this is usually referred to as richer gets richer). Obviously, the soonest a node is incorporated to the network core the most connected it will be at larger times. The solution of the BA model was found by the same authors by means of a mean field approximation4[143, 187]. In this formulation the connectivity of a node i,ki, is considered as a real continuous and time-derivable variable. Considering that new nodes are uniformly incorporated in time and that they attach mnew links, we can write the evolution equation of kias ∂ki ∂t=mΠpa i(ki)=mki PN(t) j=1kj =ki 2t,(5.23) with the initial condition ki(ti)=m, and tibeing the time when node iwas added to the network core. The solution to eq. (5.23) is ki(t)=mt ti1/2 .(5.24) 4Other solutions to this important model have been found solving the rate equation for the connectivity distribution [184–186]. 5.2. OVERVIEW OF NETWORK GENERATION MODELS 115 10-5 10-4 10-3 10-2 10-1 100 100101102103 P(k) k k−3 Figure 5.9: Degree distribution, P(k), for a BA network with N=5·104. The linear fit of the data in the log-log plot yield an exponent for the power law of γ=3. In order to obtain the degree distribution we first take the cumulative distribution, P(ki< k). From (5.24) one obtains P(ki<k)=Pti≥m2 k2t.(5.25) Finally, since the adition of new nodes is performed uniformly, the probability of finding at time ta node that was incorporated to the network at time tiis P(ti;t)=1/(m0+t). Then, the above probability (5.25) can be written as P(ki<k)=1−m2t k2(m0+t),(5.26) so that the degree distribution yiels P(k)=∂P(ki<k) ∂k=2m2t m0+tk−3.(5.27) Taking the limit when t→ ∞ we obtain the power law P(k)=2m2k−3with the exponent γ=3. In figure 5.9 we show the numerical results for the degree distribution when the BA is implemented. Analytical calculations accounting for other magnitudes have been performed. For example in [188] the authors showed that the shortest path length is smaller than that observed for ER graphs. In particular L∼log N/log(log N). Besides, the clustering coefficient in BA networks vanishes in the thermodynamic limit as happened for ER graphs. However, although the clustering decay is seen to be slower, c∼N−3/4, than in the case of ER networks, it represents a major weakness of the BA model. Variations of the preferential attachment rule (5.22) has been broadly studied after the BA model appeared. These variations try to obtain more flexible models in order to grow networks with other characteristics similar to those found in real network for which the BA fails to reproduce (existence of high clustering, the presence of degree corelations, the variety of exponents found for the ower law distribution, etc ...) while keeping the SF character induced by the preferential attachment. Some examples of these variations can be found in [186, 189–199]. 116 CHAPTER 5. NETWORK STRUCTURE AND GENERATION 5.3 Global versus local knowledge In this section, we revisit one of the main assumptions of the Barab´ asi-Albert model: the preferential attachment rule. We study a model in which the PA rule is applied to a neighborhood of newly created nodes and thus no global knowledge of the network is assumed. We numerically show that global properties of the BA model such as the connectivity distribution and the average shortest path length are quite robust when there is some degree of local knowledge. In contrast, other properties such as the clustering coefficient and degree-degree correlations differ and approach the values measured for real-world networks. As explained in Sec. 5.2.3 the first scale-free network model, introduced by Barab´ asi and Albert, postulated that there are two fundamental ingredients of many real networks [143, 187]: their growing character and the preferential attachment (PA) rule. The preferential attachment rule considers that the probability that an old node links to newly added nodes is proportional to its degree k(see eq. 5.22). However, the BA model assumes that one knows the connectivity of all nodes when a new node links to the network. This is clearly an unrealistic assumption. This drawback of the model construction has not passed unnoticed and many models have been introduced to produce scale-free networks and to test whether or not the basic assumptions of the BA recipe are necessary conditions to build up these networks [156, 161]. There are some models in which the PA rule is limited to a neighborhood due to geographic constraints [200], or where its linear character is investigated [201]. In the model decribed here, we adopt a different perspective. Our aim is to test to what extend the global character of the PA rule in the original BA model is important. We introduce a model in which the PA is applied only to a neighborhood of the newly added node depending on the value of a variable which measures the affinity between different nodes. By going down from the BA limit of the model to the the limit where all nodes are distinct, we test to what extend the global knowledge of each node’s connectivity is fundamental to get a scale-free graph. Through numerical simulations we find that in a wide range of the model parameters, average quantities such as the connectivity distribution and the shortest path length are not affected by the use of local knowledge of the network whereas other properties like the clustering coefficient are more sensitive to local details. 5.3.1 The model The model is defined in two layers. The first discriminates among all the nodes by assigning to each node at the moment of its creation a parameter aiwhich measures how close or distinct a given node is from the rest of the elements that compose the network. Then, we apply the preferential attachment rule in the neighborhood defined by nodes with common affinities. Specifically, the network is constructed by repeated iteration of the following rules: (i) Start from a small core of nodes, mo, linked together. Assign to each of these monodes a random affinity aitaken from a probability distribution, P(a). In what follows, we will use for simplicity a form for P(a) uniformly distributed between (0,1). 5.3. GLOBAL VERSUS LOCAL KNOWLEDGE 117 (ii) At each time step, a new node jwith a random affinity ajis introduced and linked to mnodes already present in the network according to the rules specified below. (iii) Search through all nodes of the network verifying whether or not the condition ai−µ≤aj≤ai+µis fulfilled, where µis a parameter that controls the affinity tolerance of the nodes. The nodes that satisfy the affinity condition are grouped in a set Aas potential candidates to gain new links. (iv) Apply the preferential attachment rule to the set A5,i.e., when choosing the nodes to which the new vertex links, we impose that the probability that vertex iconnects to the new node depends on its connectivity such that Π(ki)=ki Ps∈Aks.(5.28) (v) Finally, repeat steps (ii)-(iv) t times such that the final size of the network is N= mo+t. It is worth mentioning that the inclusion of the affinity parameter ais not a mere artifact. Indeed, most real systems are formed by non-identical elements and thus it is natural to assume that although a given node could have a large connectivity a newly created element will not link to that node because they have very little in common. This feature is clearly manifested in social networks like the WWW −where individuals bookmark different web pages accordingly to their “affinity”−or the scientist citation network [144]. In this way, it is very unlikely to find a citation in a condensed matter paper referring to a paper wrote by a psychologist. Additionally, the same argument can be translated to biological networks such as predator-prey webs or protein-protein interaction networks. Obviously, when µis large enough as to dilute the first layer of the model, we recover the BA model. The problem then consists of determining to what extend the local preferential attachment will give the same results, or in other words, does the knowledge of the entire network substantially contribute to the properties observed in the BA networks? 5.3.2 Network properties We have performed extensive numerical simulations of the model described in the preceding section. In all cases, the numerical results have been obtained after averaging over at least 500 iterations varying the system size from 103up to 1.2×104nodes. We first generate the BA network by setting the parameter µto its maximum value such that the preferential attachment applies to the entire set of nodes and then tune µin order to systematically reduce its value and therefore the size of the set Ato which the second choice eq. (5.28) is applied. Figure 5.10 shows the number of nodes with connectivity kfor several values of µ. It turns out that irrespective of the range to which the preferential attachment is applied the stationary probability of having a node with connectivity kis the same as for the BA model, namely, Pk∼k−γwith γ≈3. This result could be intuitively understood by noting 5In case that the number of elements in the set Ais smaller than mwe just add a link to all nodes in A without applying the PA rule. 118 CHAPTER 5. NETWORK STRUCTURE AND GENERATION 100 101 102 103 104 101102103 NP(k) k µ=1 µ=0.7 µ=0.2 µ=0.1 µ=0.04 Figure 5.10: Number of nodes with connectivity kfor different values of µ. The size of the network is N=104nodes and mo=m=3. The power-law distribution has an exponent equal to 3. Note that the BA limit corresponds to µ=1. that although the rules for the network generation has been changed at a local level, from a global perspective the average properties should not change radically. To realize this point, think of the network as being made up of different small components, as given by the affinity constraint, each of which is constructed following the BA algorithm. It is then clear that for large system sizes, each graph will follow the power law distribution Pk∼k−3and so will be for the entire network. The above argument applies only to average global properties, but there is nothing that guarantees a priori that the components of the network will link together in such a way that other properties will not be affected. This is the case of the average shortest path length L. As already introduced, complex networks show the noticeable property, known as small-world property, that the average path length increases at most with the logarithm of its size. We expect that for high values of µthe network is composed by a unique giant component and no fragmentation arises. When the range to which the affinity criterion is applied decreases, the network will gradually loose its compactness and will stretch approaching a one-dimensional structure with some small components. Further reduction of µprovokes the break down of the network in many isolated clusters. Figures 7.4 and 5.12 substantiate this picture. Figure 7.4 represents the ratio between the average path length obtained for different values of µand that of the BA network, for several system sizes. As µrestricts the PA range, the network undergoes a transition characterized by a growth of L(µ) an eventually becomes fragmented giving rise to an infinite shortest path length. We note here that although the results shown in the figure have been obtained for a uniform distribution of affinity values ai, the qualitative behavior does not change for other probability distributions and only the value at which the transition is 5.3. GLOBAL VERSUS LOCAL KNOWLEDGE 119 1 1.4 1.8 2.2 2.6 0 0.2 0.4 0.6 L(µ) / L(BA) µ N=2⋅103 N=4⋅103 N=6⋅103 N=8⋅103 N=1⋅104 N=1.2⋅104 Figure 5.11: Ratio between the average shortest path length for different µvalues, L(µ), and that of the BA network (L(1)) for several system sizes. The horizontal line marks the BA limit. A transition from graphs fulfilling the small-world property to a regime in which networks break down in many small pieces rising the value of L(µ) is observed. See the text for further details. (a) Pajek (b) Pajek (c) Pajek (d) Pajek Figure 5.12: Graph representations of four networks produced with different values of µ. The values of µcorrespond to (a) µ=1, (b) µ=0.2, (c) µ=0.1, (d) µ=0.04. Each network is made up of N=500 nodes. 126 CHAPTER 5. NETWORK STRUCTURE AND GENERATION Figure 5.17: Model A. Monte Carlo simulation (points) versus mean field (lines) results for ˆ kpa(t=N) as a function of the birth time t0for different values of α. The parameters of the model were N=105and A=m=m0=1. The statistics of the Monte Carlo simulations were performed using 104networks for each value of α. the incoming PA degree of a node i,ˆ kpa i dˆ kpa i dt=(1 −α)mˆ kpa i+A (1 −α)mt +AΩ(t),(5.35) (with the initial condition ˆ kpa i(ti 0)=0). Obviously, in the limit α=0 we recover the mean field equation for the Generalized Dorogovtsev model [186] (which, when A=m, describes the Barab´ asi-Albert model). For α,0 the influence of the uniform random linking is evident from the presence of Ω(t). The number of nodes that start to have the above dynamics at some time t0is dΩ(t)/dtevaluated at time t=t0which for α,0 is not constant as we have seen in the previous calculation of Ω(t). The solution of (5.35) is then given by ˆ kpa i(t=N) A=−1+exp(1 −α)mZN ti 0 dt (1 −α)mt +AΩ(t).(5.36) We have solved numerically eq. (5.36) in order to obtain ˆ kpa i(t=N) (or kpa i(t=N)= ˆ kpa i(t=N)+αm) as a function of ti 0. This function, along with the number of nodes which are incorporated to the connected set at time ti 0=t0, gives the degree distribution of the PA links. We have compared the results given by eq. (5.36) for different values of αwith the corresponding ones obtained by performing Monte Carlo simulations of the model (averaging over 104networks for each value of α). The results, plotted in figure 5.17, show a very good agreement for the mean field model and the numerical network construction. As expected, the sooner a node is incorporated to the connected set the higher its final PA degree. However, as discussed above, one can observe that this gain of the oldest nodes becomes less important when the value of αgrows due to the combination of two effects: 5.4. INTERPOLATION BETWEEN RANDOM AND SCALE FREE NETWORKS 127 (i) the application of the PA rule becomes less frequent and (ii) the fast growth of the connected set tends to make more homogeneous the PA probability of the nodes. MODEL B In the second proposal the two different linking processes are completely independent. For this, we consider that ΠPA iis a linear function of the (incoming and outgoing) links that appear as a product of the application of the PA rule. Then, kpa iwill be zero until it launches its αmPA links over the rest of the nodes, i.e. regardless of ku i. Then, the mean field equation for the evolution of kpa iis given by dkpa i dt=(1 −α)mkpa i 2(1 −α)mt +m0,(5.37) with the initial condition kpa i(ti 0)=(1−α)mand ti 0being the time when node ilaunches its mlinks. Solving the above equation yields kpa i(t)=(1 −α)m"t ti 0#1/2 .(5.38) Because the nodes launch their links at a constant rate (one node per time step), it is easy to obtain the degree distribution P(kpa) P(kpa)=2(1 −α)2m2(kpa)−3,(5.39) which is simply a power law distribution with a Barab´ asi-Albert exponent regardless of the value of α. On the other hand, the relative weight of the power law with respect to the Poisson distribution in the total degree distribution P(k) will be obviously affected by α (as the prefactor in the above equation suggests). 5.4.3 Network properties In this section we discuss the transition from SF to ER networks in terms of the global topological features of the networks. We have performed Monte Carlo simulations of the two models and compared how the relevant topological measures evolve as a function of α. We are interested in obtaining how the different correlations between the uniform and PA linking rules affect several structural measures. To do this, we have studied the behavior of three magnitudes that behave very different in the two known limiting cases (SF and ER networks), namely: the degree distribution P(k), the average shortest path length hLiand the second moment of the degree distribution hk2i. Degree distribution - The degree distribution evolution is clearly different for the two models. In figure 5.18 we have plotted the degree distribution and the rank-degree relation for both models. The rank-degree relation provides a useful tool for analysing the degree heterogeneity of the networks [204] and thus it is helpful when looking at the transition between ER and SF networks. As can be observed from figures 5.18(a) and 5.18(c) the correlated model A shows a smooth transition from the power law (α=0) to the Poisson distribution (α=1). The decay of the tails (k>> 1) of the degree distribution and the rank-degree relation becomes progressively faster as αgrows revealing the decrease of 128 CHAPTER 5. NETWORK STRUCTURE AND GENERATION Figure 5.18: Monte Carlo results for the degree distribution P(k) and rank-degree relation for several values of α. (a) and (c) show the results for model A revealing a progressive increase of the tails decaying rate when α→1. The results for model B ((b) and (d)) show how the decaying rate is not affected by α. The networks were generated with the following parameters N=105and m=m0=3 (A=3 for model A). the exponent of Ppa(kpa) as expected from the results obtained by the analytical insights developed for model A. For model B the transition is completely different as it is shown in figures 5.18(b) and 5.18(d). In both representations the decaying rate of the tails is independent of αand the transition to the Poisson distribution is much more apparent for low values of k. In this sense one can conclude that highly connected nodes persist along the transition of model B while for model A the heterogeneity is progressively lost. Average shortest path length - The different evolution of the degree distributions observed above suggests to look at how the average shortest path length behaves along the two paths of interpolation. It is well known that the existence of high degree nodes makes the network more compact due to the possibility of finding shortcuts between nodes going through the hubs. Hence, the persistence of highly connected nodes determines the small diameter of the scale-free network. The results obtained are shown in figure 5.19(a). As expected, the average shortest path length as a function of αgrows slower for model B because the probability of finding hubs is higher than for the networks generated using model A for the same value of α. Second moment of P(k) - In order to obtain a quantitative measure of the evolution of the degree heterogeneity for the two models it is convenient to measure the second moment of the degree distribution, hk2i. This magnitude diverges (in the thermodynamic limit N→ ∞) for scale-free networks with exponents between 2 and 3. So, we expect a decrease of the heterogeneity on the path to ER graphs. As can be observed from 5.4. INTERPOLATION BETWEEN RANDOM AND SCALE FREE NETWORKS 129 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 α Model A 1 1.5 2 2.5 3 3.5 0 0.2 0.4 0.6 0.8 1 α Model A Model B 0.8 0.85 0.9 0.95 1 Model A 0.8 0.85 0.9 0.95 1 Model A Model B <L> <k > 2 2 <k > ER <L>ER (a) (b) Figure 5.19: Average path length (a) and second moment of the degree distribution (b) as a function of α. Both quantities are represented normalized by their respective values in the ER limit. The results clearly manifest the two different transitions of the models regarding the heterogeneity evolution along the interpolating path. The averaged networks had the following parameters N=104and m=m0=3 (A=3 for model A). figure 5.19(b), model A shows a faster decrease of hk2ias expected from the study of the degree distribution while for model B the transition is much smoother revealing again the persistence of highly connected nodes along the path to the ER limit. As for other properties like the clustering coefficient and degree-degree correlations we have checked that they remain unchanged irrespective of the value αand wheter model A or B is implemented. The present model provides a useful tool to study the influence of the degree of heterogeneity in dynamical processes of different kinds just as the Watts-Strogatz model have proved to do so in the transition from regular to random structures. In particular, there exist open questions in phenomena such as the synchronization of coupled oscillators [205] where this kind of model could be particularly relevant to explore the system’s behavior in the region where homogeneous and heterogeneous architectures coexist. This question will be deeply analysed in section 8.3. Chapter 6 Propagation through Complex Networks The better a simulation is for its own purposes, by the inclusion of all relevant details, the more difficult it is to generalize its conclusions for other species. For the discovery of general ideas in ecology, therefore, different kinds of mathematical descriptions, which may be called models, are called for. Whereas a good simulation should include as much detail as possible, a good model should includes as little as possible. — J. Maynard Smith in Models in Ecology [206]. In this chapter we will focus on two of the main dynamical processes studied on top of complex networks, namely, the analysis of Epidemic spreading and Information dynamics. The interest of studying these problems is twofold. First, the simplicity of the description of the two procceses allow for analytical results, heuristic insights and extensive numerical simulations in order to explain the role that the underlying topology has on the dynamics. Then, one of the advantages of studying these dynamics is that the simple formulation of the models (usually expressed by means of linear rules) used for their description does not mask the effects of the topological complexity. Besides, one can realize by looking at the literature that a great number of the networks whose characterization is available (mostly due to the simplicity for unveiling the links between their components) can be regarded as either technological or logistic networks. Then, the study of epidemics and information propagation is justified for practical purposes. 6.1 Epidemic spreading and Immunization The history of the studies on epidemic spreading starts with the first works by epidemiologist at the beginning of the 20th century [207]. However, the burst in the mathematical modeling of disease transmission took place in the middle of the 20th century by the formulation of a large variety of models (interesting books on the matter are [208–211]) 132 CHAPTER 6. PROPAGATION THROUGH COMPLEX NETWORKS aimed at reproducing the evolution patterns of the number of casualties and infected people during epidemic periods (see figure 6.1). Recently, the attention has been redirected to the spreading of informatic viruses. The interest in this field has been coupled to the availability of data about potential transmission networks (like the internet or peer-to-peer networks). The development and deployment of a digital immune system to prevent technological networks from the spreading of viruses and to minimize the damage produced by intentional attacks are in the root of recent research efforts [158, 212–222]. In this section we will first introduce two general models (SIR and SIS) that describe the spread of epidemics on homogeneous systems. Then, we will turn our attention to the disease transmission in heterogenous substrates and the performance of different immunization strategies will be compared. Finally we will report on a new immunization strategy based on the covering problem of complex networks. The performance of this new algorithm depends on the local structure of the network. We will implement this strategy along with the afore mentioned in order to compare their results when deployed on top of real networks. 6.1.1 Modeling epidemic spreading There are many different models to describe the epidemic transmission problem. However, nearly all of them are variations of some general and coarse grained models like the SIR (Susceptible-Infected-Removed) and SIS (Susceptible-Infected-Susceptible). The different variations have to do with an increased complication of the models to study particular diseases. In order to focus on the importance that the topology of the network has on the spreading of a disease we will deal with the most simplified descriptions of the epidemics dynamics. The SIR model The SIR model was introduced by Kermack and Mc Kendrick in 1927 [207] to explain the rapid rise and fall in the number of infected patients observed in epidemics such as the plague (London 1665-1666, Bombay 1906) and cholera (London 1865). This model was 200 175 150 125 100 75 50 25 0 13 12 11 10 9 8 7 6 5 4 3 2 1 0 Deaths Week Figure 6.1: Histogram of the number of deaths due to the influenza-pneumonia epidemic during the 1968-1969 winter in New York. This extremely damaging influenza was named the “Hong Kong flu” due to the place where it started. 6.1. EPIDEMIC SPREADING AND IMMUNIZATION 133 recovered by the work of Anderson and May [223] after being oversought for decades. The SIR model is a typical example of the so-called compartmental models. In this class of models the elements are viewed as parts of several groups (or compartments) so that the evolution equations are referred to the number of elements of each group. The SIR model describes the spreading of infectious diseases in which each individual can be either immunized or dead after the contagion. Following this assumption we can classify the population into three different groups: •Susceptible: Those healthy people who have not been infected and thus are likely to contract the disease in the future. •Infected: People who has been contagied and are currently suffering the effects of the disease. They can infect Susceptible people in the course of their disease. •Recovered and Removed: Composed by people that finally died due to the disease or recovered and got immunized. Then, individuals can change their state by means of the jumps between the three compartments, S→I→R. The dynamical rules accounting for the flux among the three states determine a set of differential equations for the densities of the population groups s(t)=S(t)/N,i(t)=I(t)/Nand r(t)=R(t)/N. The change rate for susceptible elements is always negative and proportional to the number of contacts among infected and susceptible elements. We call λthe probability that one susceptible individual gets infected in one contact, then we can write ds dt=−λhkis(t)i(t).(6.1) The evolution of the proportion of infected individuals, i(t), has two contributions, one positive −˙s(t) and one negative accounting for the recovering (or death) rate of the infected individuals di dt=λhkis(t)i(t)−µi(t),(6.2) where µis the recovering (or death) rate that corresponds to the inverse of the average disease time for an individual. Taking into account that r(t)=1−s(t)−i(t) the last evolution equation for the recovered density is dr dt=µi(t).(6.3) The above formulation of the model equations assumes the homogeneous mixing hypothesis that considers that the set of susceptible people with whom an infected individual establishes contacts is taken at random within the whole population. This is manifested in the constant value for the number of contacts hkiso that the approach is only valid for homogeneous networked systems. Along with this assumption we have considered homogeneity in the agent characteristics so that λand µ(although seen as averages) are meaningful. This model is seen as a mean field approximation to the epidemic spreading problem. One can modify the SIR model by adding more compartments (like e.g. in models for VIH propagation where a set of people suffering an incubation period, or more 134 CHAPTER 6. PROPAGATION THROUGH COMPLEX NETWORKS 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Epidemic percolation *r <k> λµ Healthy state 0.2 0.4 0 0.6 0.8 1 λ/µ =0.5/<k> λ/µ =1/<k> λ/µ =2/<k> r* 1/ λ λ λ Figure 6.2: Graphical solution of r⋆given by last expression in eq. (6.4). It can be observed the emergence of the second solution (corresponding to the epidemic percolation) for λ/µ > hki−1. The inset shows a qualitative plot of the phase transition for the SIR model. technically a latency period, should be distinguished from those who have the disease already diagnosed) or by considering that the time scale involved is slow enough so that additional terms accounting for natural birth and death rates should be incorporated. We can rescale conveniently the above equations (µ→1, t→µt,λ→λ/µ) in order to have a single control parameter λto study the behaviour of the model. The question to answer in the SIR model has to do with the conditions under which a small infectious seed leads to a significant fraction of individuals at the recovered state when the steady state (i(t)=0) is reached for the whole system. One can translate this question in terms of a general bond percolation problem since an epidemic is set when there is a large enough fraction of “occupied” bonds (those that the epidemy used to spread) to lead to the formation of a network component whose size scales with the size of the graph (signaling an epidemic percolation). In fact, there exist an exact mapping between both problems [224, 225] that can be used as a powerful tool for solving the epidemic spreading across general networks [226, 227]. In order to answer this question we consider the initial conditions s(0) =1/N,i(0) =r(0) =0 and we look for the value r⋆=limt→∞ r(t) It is easy to notice that r⋆=1−s⋆because i⋆is necessarily null. Then, dividing eq. (6.1) by eq. (6.3) to get rid of i(t) we have ds dr=−λhkis(t)−→ s(t)=exp(−λhkir(t))−→ r⋆=1−exp−λhkir⋆.(6.4) Last equation has always r⋆=0 as a solution (no epidemic percolation) and if R0= λhki>1 there is a second one with r⋆,0 corresponding to a significative spread of the disease. R0is usually termed as the effective reproductive rate and its physical meaning is clear: it corresponds to the average secondary infections produced when a single infected individual is introduced in a healthy population. If this number is greater than one the disease reaches a non null fraction of the population. From the phase transitions point of view one would speak about r⋆as the appropriate order parameter and λas the control parameter. It can be obtained that at the critical point, λ=hki−1,r⋆behaves as r⋆∼ (λ−hki−1) so that the critical exponent is 1 as expected from a mean field treatment. 6.1. EPIDEMIC SPREADING AND IMMUNIZATION 135 The SIS model The SIS model was originally introduced by Hethcote and Yorke [228] in 1984 for describing the propagation of gonorrhea and has been largely used for studying the transmission of tuberculosis. These diseases have several common attributes that make it different from other infections. The most important difference is that the infection does not confer immunity to recovered subjects so that the SIR model is no longer valid. The SIS model considers only two compartments composed of susceptible and infected nodes so that a continuous flux between both compartments is allowed, S↔I. Then, the relevant equation for the SIS model is di dt=λhki[1−i(t)]i(t)−µi(t),(6.5) where we have assumed s(t)=1−i(t) (i.e. no deaths associated to the disease are considered). Again we can rescale the equation in order to obtain µ=1. The SIS model is analogous to the SIR model in what refers to the existence of a epidemic transition. However in the SIS model the two regimes are differentiated by whether the disease persist indefinitely in the population (due to the fact that subjects can be reinfected many times) or not. Then, imposing di/dt=0 in (6.5) we obtain two different steady states, one with i(t)=0 and the second one, for the case λ > hki−1, with i(t)=(λ−1/hki)/λ corresponding to the an endemic state. Due to the manifest analogy between SIR and SIS models we will focus on the SIR formalism in the forthcoming discussions about the behaviour of epidemic spreading on heterogenous networks. It has been shown that the same qualitative results are also obtained for both models when studying more complex topologies. 6.1.2 Epidemic spreading in general complex networks The homogeneous mixing hypothesis assumed in the previous discussion cannot be applied to many real systems like technological networks where their components do not interact with a similar number of network elements. Then, it is necessary to incorporate the ingredient of heterogeneity to the problem of epidemic spreading. In the SIR model (and analogously for the SIS formulation) the different compartments S(t), I(t) and R(t) are now characterized by the subsets {Sk(t)},{Ik(t)}and {Rk(t)}labeled by the connectivity kof their components [212–214, 216, 227]. These variables are normalized so that Sk(t)+Ik(t)+Rk(t)=1, then S(t)=NX k P(k)Sk(t),I(t)=NX k P(k)Ik(t) and R(t)=NX k P(k)Rk(t).(6.6) This second classification of the compartment elements into degree classes allows to take into account the heterogeneous character of real networks. At a mean field description the evolution of these new magnitudes satisfy the following Chapter 9 Conclusions We want to conclude with a brief summary of the most relevant results obtained along the three parts of this Thesis. We want to stress here their relevance as well as some prospective research that arises in the light of these results. In the first part we have studied the issue of intrinsic localization (discrete breathers) in nonlinear Schr¨ odinger anharmonic lattices (described by the Salerno model). The major achievement of these studies is the computation of exact mobile localized modes. For these computations, it was important to develop a generalized continuation method, that can be thought of, as the natural extension of those employed for computing standard (pinned) discrete breathers. The generalized continuation method allows to obtain, in a highly systematic way, families of mobile, oscillating and vortical discrete breathers (as well as pinned ones). The problem on the existence of mobile discrete breathers has been extensively discussed after the theory for pinned localized modes was successfully developed. The use of collective variable methods and numerical simulations of perturbed pinned solutions lacks the precision required to obtain general arguments about the possibility of having mobile localized states in nonlinear lattices. However, the computation of mobile discrete breathers in this Thesis is neither unbiased (i.e. based on a priori ansatzes) neither suffers from low numerical accuracy. On the contrary, our continuation procedure computes mobile discrete breather solutions as exact fixed points solutions of a map and, therefore, the unique requirement is that the Jacobian of the map is invertible so that the iterative method converges to the desired solution. Concerning mobile breathers, our results point out that, except for integrable and other exceptional (e.g. vanishing Peierls-Nabarro barrier) situations, mobile localized states in nonlinear Schr¨ odinger lattices are described by a localized part (the core) and an asymptotically extended background composed of plane waves. We have obtained numerical evidences of the importance of this extended background in the core mobility. In particular, we have shown how the Peierls-Nabarro barriers that a mobile breather experiences periodically in its motion across the lattice are surpassed with the help of the energy balance core-background. In this sense, we have observed that the higher is the PeierlsNabarro barrier, the higher is the energy flux between core and background and the higher is the amplitude of the extended background. These observations reveal the essential role of the background in core mobility pointing out that collective variables approaches are incomplete when considering only those degrees of freedom relative to the core. 239 240 CHAPTER 9. CONCLUSIONS The study of the whole Salerno lattice, both in its oneand two-dimensional versions, has provided a detailed account of the existence and properties of discrete breather solutions and, at the same time, it has proved the versatility of the continuation method. Several questions remain open after these studies. In particular, it would be desirable to perform a deep analysis on the mobility of two-dimensional discrete breathers and more exotic solutions like trains of discrete breathers or vortical states. It would be also interesting to apply the continuation methods to other important nonlinear lattices such as Klein-Gordon or Peyrard-Bishop-Dauxois models. Finally, it is also interesting to develop a collective variable theory accounting for the relevant ingredients of the mobile solutions found. In the second part we have studied the structure of complex networks and the performance of propagation dynamics on top of them. Several results have been obtained for each of these two issues. We have first presented two models of network construction that provide two network families where only a few topological characteristic vary significantly among the members of these families. The purpose of these models is to provide a useful tool for analyzing the role that these changing structural properties have on the performance of different network dynamics. In fact, these models have been used for this purpose along the Thesis. Whereas one model preserves the scale-free character of the degree distribution while the clustering coefficient and average characteristic path length are varied, the second one allows the degree distribution to interpolate between the scale-free and Poisson distributions while other magnitudes remain comparable. This latter model would be very useful to shed new light on the roots of the different phenomena found when dealing with homogeneous and heterogeneous topologies. The studies on the propagation dynamics on networks have been focused on two important dynamics, namely, the SIR model for epidemic spreading and the analysis of coarse-grained information routing models. The main objective of these two studies is to analyze the efficiency of different routing and immunization algorithms depending on the substrate topology. For the studies performed on epidemic spreading the main results concern the development of a new immunization method based on the d-covering problem. We have implemented an heuristic method for finding the nearly smallest subset of network’s elements that should be covered so that every node in the network has at least one covered node at a distance less than or equal to d. The results show that, depending on the degree-degree correlations of the network, the obtained solution is very different. Besides we have shown that the obtained solution of the d-covering problem, when thought of as immune nodes to a SIR epidemics, yields a very efficient algorithm compared to those already existing in the literature. The efficiency of the d-covering subset of immune nodes also depends strongly on the correlations of the network when a SIR-like epidemics is studied. The study of information routing dynamics also yields interesting results. In particular, the main result concerns the study of a congestion-aware strategy for the routing of information packets across the network. The use of shortest-path strategies in scale-free networks lead to a fast congestion of highly connected nodes and hence to the development of jamming for low levels of injected information. We have obtained a more robust routing policy making use of an effective distance that takes into account the congestion level of the network at the local scale. However, the shift of the onset of jamming 241 is achieved at the expenses of a sudden growth of the congestion levels at the jammed phase. We have explained the microscopic origins of this first order like phase transition as a consequence of the effective fragmentation of the network. This fragmentation is due to the formation of dynamic walls composed of those nodes that do not allow to receive data packets from their neighbors, due to their high level of local congestion. We have seen in this second part several examples on the relation of topological characteristics like the clustering coefficient, the average path length and the degree-degree correlations on the development of two simple model dynamics of importance for scalefree networks. Besides, the modeled algorithms for immunization and data routing have been designed for a nearly optimal deployment when applied on top of highly heterogeneous networks. The nature of the simple dynamics studied here and their application to human-made real systems allows our models to be reliable in these kind of networks. This is not anymore the case of real biological networks where both topology and (nonlinear) dynamics are imposed to the system. This extreme has been analyzed in the third part of the Thesis. The third part of this thesis is devoted to the study of nonlinear dynamics on top of complex networks. It is thought of as the confluence of the above two parts because it applies the tools obtained from the studies on nonlinear localization in homogeneous lattices and the analysis of complex networks structure. Along with the results obtained in this part, the mixed use of these tools constitutes a, somewhat, novel feature since the study of nonlinear complex networks is still in its infancy. We have studied two different nonlinear systems: a Michaelis-Menten regulatory dynamics (where activatory and inhibitory terms compete) and the paradigmatic Kuramoto model of coupled phase oscillators. In these two problems, related to diverse natural systems, the main purpose is to unveil the relation between the networked structure of the systems and the function they fulfill. The search for this “Structure-Function” connection is based on the assumption that the evolution of the real networks is the result of a kind of optimization for the performance of their function. Then, a first step is to analyze coarse-grained synthetic systems modeled by relevant nonlinear dynamics. The studyof activatory-inhibitory regulation, modeledbymeansofageneralMichaelisMenten dynamics between the network nodes, allows to approach the problem of genegene regulation. In this sense, some important results are related to the experimental observations of this kind of systems. The first important result concerns the fragmentation of the network into independent dynamical clusters while the rest of the network remains at the rest (zero activity) state. The dynamics of these dynamical islands show a very rich dynamical behavior: steady, periodic and chaotic states. When these emergent dynamical clusters of self-sustained non-zero activity are considered as networks defined by its nodes and the links among them, new topological features, different from those of the underlying network, are obtained. In this regard, the most important finding is a clustering coefficient for the dynamical islands much higher than that of the substrate network (a Barab´ asi-Albert scale-free network). A second important result is obtained when looking at the observed bifurcations. Periodic clusters display either period doubling or tripling bifurcations on their route to chaos. Analyzing the shape of the Floquet eigenvectors associated to these bifurcations it is possible to determine which nodes are responsible for the transition from the old (period 1) attractor to the new (period 2 or 3) one. This method allows us to observe that, differently from other processes on networks, nodes’ substruc- 242 CHAPTER 9. CONCLUSIONS tures and not single nodes are responsible for the evolution of the dynamical clusters. The second objective of this third part is the study of the synchronization paths in networks of Kuramoto phase oscillators. This study is performed on a variety of networked substrates, namely, the two networks families introduced in the first part of the Thesis and structured networks. For this purpose, we have introduced a new order parameter that allows to unveil the local patterns of the synchronized clusters that emerge as the coupling strength is increased. In this sense, the main result is obtained when comparing the synchronization paths in homogeneous and heterogeneous networks. The results point out that the emergence of the giant synchronized cluster for Erd¨ os-R´ enyi networks is the result of the coalescence of multiple small synchronized clusters. This simultaneous clusters’ collapse is thus produced in a narrow region for the coupling parameter so that the degree of global synchronization is rapidly increased from zero near the synchronization onset. On the other hand, for scale-free networks the process is described by a gradual growth of the giant synchronization cluster. This synchronization cluster, organized around the central hubs of the networks, grows by incorporating more and more synchronized nodes as the coupling is increased. As a result, the onset of global synchronization occurs earlier (in terms of the coupling strength) than in the Erd¨ os-R´ enyi case. However, the global coherence in scale-free networks grows, far from the synchronization onset, at a much slower rate than in the case of Erd¨ os-R´ enyi graphs due to the (one-by-one) additive growth for the synchronized cluster. These two works constitute interesting examples on the interplay between Function and Structure. In the first study it is clear that nonlinear dynamics shows up an emergent structure with new topological characteristics. For the second work it is shown that, depending on the underlying structure, radically different patterns of synchronization are obtained. Therefore, the importance of tackling the combined study of both structural complexity and nonlinear dynamics is clear, since a separate analysis would be incomplete. The mutual influence observed thus prevent from going from one to the other or vice versa. The continuation of the presented work would be carried following different directions. Perhaps, the most ambitious direction is to make one step further into the analysis of real networks. For example, the availability of gene expression data (although one must be careful and selective with the large amount of experimental data sets) motivates the study of real gene regulatory networks in order to apply the tools and results found in this part. At the same time, other kind of relevant nonlinear dynamics, like models of neural activity (Hopfield, integrate and fire, etc...), could be also studied by means of similar techniques in order to obtain more examples on the interplay of structure and dynamics. The results presented in this Thesis are intended to analyze and understand several phenomena displayed when two essential ingredients of complexity are present (both separated and combined). It would take still a long time before the understanding of simple dynamics and models allows to go one step farther and attack the unification of these two ingredients in order to have a framework within which one can solve the “StructureFunction” problem. In this sense, the present work provides some examples and tools about how the problem could be first tackled. As all research work that does not suffice to fully unravel the features of a given problem, our work also motivates further studies on this interesting question that we intend to pursue in the years to come. Appendices We want to add two appendices about the computation and stability characterization of periodic orbits since these solutions have extensively appeared throughout this Thesis. Although for particular types of periodic solutions and dynamical systems, the results reported here can be further extended we have tried to briefly summarize the essential features about these two issues. A Computation of Periodic Orbits The computation of periodic solutions to a set of Ncoupled nonlinear differential equations ∂~ x ∂z=~ F~ ξ[~ x],(A.1) where ~ xare the variables of the system and ~ ξdenote the parameters of the particular equations, can be formulated as a problem of finding the solution of a system of Nnonlinear equations, with Nvariables xi(i=1,...,N) and a set of mparameters ξi(i=1,...,m), of the form ~ G(~ x;~ ξ)=0.(A.2) As we introduced in section 2.2, let us consider a periodic solution ~ x~ ξas a fixed point solution of a N-dimensional map M M~ yn=~ yn+1,(A.3) where the map Mcan be constructed using the z-evolution operator (zis usually time or space) given by equations (A.1) over a (time or space) period Twhen we are looking for z-periodic solutions M=T~ ξ,T,(A.4) or a combination of an index (lattice) displacement and a z-evolution operators when looking for combined periodicities as in (p,q)-resonant states for time and lattice displacement, eq. (2.24), M=LpTq ~ ξ,T.(A.5) Given the particular definition of the map M, the desired periodic solution will satisfy eq. (A.2) in the form ~ G(~ x;~ ξ)=M~ x~ ξ−~ x~ ξ=0.(A.6) One is typically interested in a particular solution corresponding to a special choice of the parameters ~ ξbut, on the other hand, the only available solution corresponds to a 243 244 APPENDICES simplified version of the system corresponding to ~ ξ0. In these cases the solution can be found by means of a homotopy procedure [296]: given a known solution, ~ x~ ξ0, to some special choice of parameters, ~ ξ0, the solution, ~ x~ ξ, is computed via the computation of intermediate solutions to a chain of equations with parameters, ~ ξ0→~ ξ1→... → ~ ξn−1→~ ξn=~ ξ. The latter path in parameter space is conveniently used so that every intermediate solution can be found. There are several methods used for solving each step in the chain of equations and nearly all of them make use of the solution found for the latter system as the ansatz for the analytical or numerical methods used at each step. The homotopy strategy is based on the implicit function theorem that assures the existence of an unique solution ~ x~ ξn, so that ~ G(~ x~ ξn;~ ξn)=0, when there exist a solution ~ x~ ξn−1 (~ G(~ x~ ξn−1;~ ξn−1)=0) and ~ ξnbelongs to an open set centered at ~ ξn−1. The conditions that must be fulfilled are: (i) ~ G(~ x;~ ξ) is continuous on an open set centered at (~ x~ ξn−1,~ ξn−1). (ii) The Jacobian determinant of ~ G(~ x;~ ξ) evaluated at (~ xn−1 ~ ξ,~ ξn−1) is non-null, DethD~ xG(~ x~ ξn−1;~ ξn−1)iij=Det ∂Gi(~ x;~ ξ) ∂xj(~ x~ ξn−1,~ ξn−1) ,0.(A.7) In order to satisfy the local convergence conditions of the theorem, the homotopic computation must be carried by dividing the path toward the desired solution into as much intermediate steps as necessary. In this way, the solution for ~ x~ ξnwould not differ very much to that for ~ x~ ξn−1so that expressing the new solution as ~ x~ ξn=~ x~ ξn−1+~ ∆one would write ~ G(~ x~ ξn−1+~ ∆;~ ξn)=0=~ G(~ x~ ξn−1;~ ξn)+D~ xG(~ x~ ξn−1;~ ξn)~ ∆ + . . . , (A.8) neglecting those higher order terms in ~ ∆. From the above expression one can obtain the difference ~ ∆between the old and the desired solution by just computing the inverse of the Jacobian matrix D~ xG ~ ∆ = −[D~ xG(~ x~ ξn−1;~ ξn)]−1·~ G(~ x~ ξn−1;~ ξn).(A.9) Due tothe errorat the truncationin Taylor expansion (A.8), one mustiterate thisprocedure until the desired convergence (bounded by machine precision) is reached,  ~ G(~ x~ ξn;~ ξn)< ǫ. For this purpose, one uses as the new trial solution the one obtained by the last computation of ~ ∆. Calling ~ yi ξnthe trial function used at the ith stage of the iterative computation of solution ~ xξnand ~ ∆ithe obtained solution of eq. (A.9) at this stage, a schematic picture of the whole process for computing ~ x~ ξnfrom the initial ansatz ~ x~ ξn−1(the solution of the A COMPUTATION OF PERIODIC ORBITS 245 before equation in the homotopy chain) is Stage 1 ~ x~ ξn−1=~ y1 ~ ξn A.9 −→ ~ ∆1−→ ~ y2 ~ ξn=~ y1 ~ ξn+~ ∆1 Stage 2 ~ y2 ~ ξn A.9 −→ ~ ∆2−→ ~ y3 ~ ξn=~ y2 ~ ξn+~ ∆2 . . .. . .. . . Stage m~ ym ~ ξn A.9 −→ ~ ∆m−→ ~ ym+1 ~ ξn=~ ym ~ ξn+~ ∆m . . .. . .. . . Stop when  ~ G(~ yi ~ ξn;~ ξn)< ǫ −→ ~ x~ ξn=~ yi ~ ξn. The convenience of using this iterative process relies on its quadratic convergence but, on the other hand, one must posses a good ansatz for the initial trial function (since no global convergence is assured) and hence a homotopic continuation to the desired solution is required. For the particular situation when one is interested in the computation of purely zperiodic solutions (such as discrete breathers for the case of time periodic solutions), the equation to solve would be written as G~ x~ ξ(z0);~ ξ=T~ ξ,T~ x~ ξ(z0)−~ x~ ξ(z0)=0,(A.10) where z0stands for the z-origin of integration. Therefore, eq. A.9 will take the form ~ ∆ = −nhD~ xT~ ξn,T~ x~ ξn−1(z0)i−Io−1·G~ x~ ξn−1(z0);~ ξn,(A.11) where the matrix D~ xT~ ξ,h~ x(z0)is computed integrating from z0to z0+hthe equations that are obtained by deriving eq. (A.1) respect to the initial conditions, ~ x(z0),  ∂D~ xT~ ξ,z~ x(z0) ∂zij =∂2xi(z) ∂z∂xj(z0)= N X k=1 ∂hF~ ξ~ x(z)ii ∂xk(z) ∂xk(z) ∂xj(z0)= = N X k=1Aik ∂xk ∂xj(z0)=hA·D~ xT~ ξ,z~ x(z0)iij ,(A.12) with the initial condition D~ xT~ ξ,z0~ x(z0)=I,i.e. integrating the linearized equations around the solution, ~ x~ ξ(z). The matrix D~ xT~ ξ,h~ x~ ξ(z0)provides a map between an initial perturbation of the solution, ~ δx~ ξ(z0), and its evolution up to z0+h, ~ δx~ ξ(z0+h)=hD~ xT~ ξ,h~ x~ ξ(z0)i·~ δx~ ξ(z0).(A.13) For a z-periodic solution the elements of matrix Ain eq. (A.12) are z-periodic functions with the same period Tand hence D~ xT~ ξ,qT ~ x~ ξ(z0)=hD~ xT~ ξ,T~ x~ ξ(z0)iq,(A.14) 246 APPENDICES with qinteger. This implies that it is enough to integrate the evolution of the linearized equations over a period Tand obtain D~ xT~ ξ,T~ x~ ξ(z0)for characterizing the time evolution of the perturbations after an integer number of periods. For a time periodic solution the socalled Floquet matrix, F=D~ xT~ ξ,T~ x~ ξ(t0), contains all the information about the linear stability of a periodic solution (we will discuss it later in Appendix B). As the implicit function theorem states, the key point for being able to compute the solutions to eq. (A.10), is that the Jacobian is invertible so that eq. (A.11) can be solved. In the particular case of time periodic solutions this implies that the eigenvalues of the Floquet matrix, F, must be different from +1. In other words, if the spectra of the Floquet matrix contains perturbations that are also time periodic with the same period of the solution, T, their associated eigenvalue will be +1 and therefore we could not solve eq. (A.11). The existence of such perturbations causes the degeneracy of the linear problem since given any solution of period Tone can construct another one by adding any combination of such perturbations and therefore the solution is not unique. Although the existence of degeneracies depends on the particular dynamical system, when we are dealing with autonomous dynamical systems, such as eq. (A.1), there always exist one Floquet eigenvalue +1. The associated eigenvector is related with the time translation invariance of the solutions, ~ δx~ ξ(t)=˙ ~ x~ ξ(t), and acts translating the instant solution across its path in the phase space. However, the non invertibility of the Jacobian does not prevent us from implementing the continuation scheme. The common way used to solve this problem is to restrict the rank of the Jacobian matrix to a subspace orthogonal to its kernel. This restriction does not influence the efficency of the method since the kernel vectors correspond to directions in tangent space that convert the solution into itself. There exist several ways for restricting the Jacobian rank that depend on the particular properties of the periodic solution. For example, in the case of time-reversible orbits, i.e. those which are invariant under the transformation R(~ q,~ p,t)→(~ q,−~ p,−t) (where ~ qand ~ pdenote the two sets of canonically conjugated variables of the system), one can fix the time origin without loss of generality setting ~ p=0 [81]. With this restriction we prevent from perturbations inside the same periodic manifold and, as a plus, we have reduced the degrees of freedom to the half, N/2. However, this tricky method does not allow to compute other kind of periodical orbits which are not time-reversible (as mobile discrete breathers) and the use of other methods such as the Singular Value Decomposition [82, 83] is required. B Linear Stability of Periodic Orbits Before analysing the stability analysis for periodic solutions let us briefly focus on the stability characterization of general orbits of dynamical systems. The linear system of equations defined in (A.12) defines the most general tool for characterizing the stability of dynamical systems solutions: the Lyapunov exponents which are the eigenvalues, {µj} (j=1, ..., N), of matrix Athat define the system of linear differential equations for the evolution of linear perturbations. In particular, the general definition of Lyapunov exponents can be expressed in terms of the eigenvalues of matrix D~ xT~ ξ,t~ x~ ξ(t0)(which for B LINEAR STABILITY OF PERIODIC ORBITS 247 x(t) x(t+h) x(t+2h) |d| |d | |d | |d| |d| 1 2 Figure B.1: Schematic representation of the rescaling procedure for the computation of the largest Lyapunov exponent µ1. period Tsolutions and t=Tis the Floquet operator), {λj(t)}, as µj=lim t→∞ 1 tln|λj(t)|.(B.1) Although there are several techniques for computing Lyapunov exponents, it is somehow a hard task since it implies large integration times in order to get accurate values independent of the time origin choice. However, one is most of the times interested in the value of the largest Lyapunov exponent (say µ1) which is the easiest to calculate due to the tendency of any perturbation to grow towards the direction associated to the largest Lyapunov exponent (see [296, 333]). The largest Lyapunov exponent, µ1, tell us whether the solution repeals nearby orbits (perturbations), µ1>0, and thus the solution is regarded as chaotic. In the case of stable periodic orbits of autonomous dynamical systems the maximum Lyapunov exponent is always 0, this corresponds to the Floquet eigenvalue +1 associated to the time translational invariance. The largest Lyapunov exponent can be expressed as µ1=lim t→∞ 1 tlnD~ xT~ ξ, t~ x~ ξ(t0)~ δx~ ξ(t0)  ~ δx~ ξ(t0)  .(B.2) This expression, which makes uses of the ratio of separation of an initial perturbation ~ δx~ ξ(t0) after large times, turns out to be helpful for computing µ1. The computational method consists in making a perturbation of the solution with a tangent vector of arbitrary direction and modulus |d|, and follow the evolution of the perturbed orbit for a time interval, h. Then, the distance, |d1|between the original and the perturbed orbits is measured. At the same time h, the perturbed solution is varied by preserving its direction but being rescaled to |d|. This process is iterated for a number of times (see figure B.1) so that a set of distances {|di|}with (i=1, ..., n) is collected. Finally, averaging these measures one obtains the largest Lyapunov exponent µ1=1 nh n X i=1 ln|di| |d|.(B.3) 254 BIBLIOGRAPHY [55] M. O  M. J  A. E, Enhanced mobility of strongly localized modes in waveguide arrays by inversion of stability, Phys. Rev. E 67, 056606 (2003). [56] A. A. V and Y. B. G, On the motion of solitons in discrete molecular chains, Theor. Math. Phys. 68, 873 (1987). [57] C. C, Y. S. K, O. K, and K. H. S, Moving localized modes in nonlinear lattices, Phys. Rev. B 47, 14228 (1993). [58] D. C  A. R. B  N. G-J, Localized states in discrete nonlinear Schr¨ odinger equations, Phys. Rev. Lett. 72, 591 (1994). [59] D. C  A. R. B  N. G-J, Perturbation theories of a discrete, integrable nonlinear Schr¨ odinger equation, Phys. Rev. E 53, 4131 (1996). [60] R. S. M and J. A. S, Effective hamiltonian for travelling discrete breathers, J. Phys. A 35, 3985 (2002). [61] M. S, Quantum deformations of the discrete nonlinear Schr¨ odinger equation, Phys. Rev. A 46, 6856 (1992). [62] R. K. B  P. J. C (), Solitons, Springer, (Berlin, Ge), 1980. [63] G. E,Solitons, Springer, (Berlin, Ge), 1981. [64] V. E. Z and A. B. S, Exact theory of two-dimensional self-focusing and onedimensional self-modulation of waves in nonlinear media, Sov. Phys. JETP 34, 62 (1972). [65] A. S, S. F, S. G, and S. R. S, Quantum coherent atomic tunneling between two trapped Bose-Einstein condensates, Phys. Rev. Lett. 79, 4950 (1997). [66] F. D, S. G, L. P. P, and S. S, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999). [67] J. D, J. E. S, D. L. F, C W. C, L. A. C, J. C, L. D, E. W. H, K. H, W. P. R, S. L. R, B. I. S, and W. D. P, Generating solitons by phase engineering of a Bose-Einstein condensate, Science 287, 97 (2000). [68] A. T and A. S, Self-focusing and defocusing in waveguide arrays, Phys. Rev. Lett. 86, 2353 (2001). [69] H. S. E, Y. S, R. M, A. R. B, and J. S. A, Discrete spatial optical solitons in waveguide arrays, Phys. Rev. Lett. 81, 3383 (1998). [70] R. M, U. P, J. S. A, H. S. E, and Y. S, Dynamics of discrete solitons in optical waveguide arrays, Phys. Rev. Lett. 83, 2726 (1999). [71] R. M, H. S. E, Y. S, M. S, and J. S. A, Selffocusing and defocusing in waveguide arrays, Phys. Rev. Lett. 86, 3296 (2001). [72] J. W. F, M. S, N. K. E, and D. N. C, Observation of two-dimensional discrete solitons in optically induced nonlinear photonic lattices, Nature 422, 147 (2002). BIBLIOGRAPHY 255 [73] M. A and J. L, A nonlinear difference scheme and inverse scattering, Stud. Appl. Math. 55, 213 (1976). [74] M. A and J. L, Nonlinear differential-difference equations and Fourier analysis, J. Math. Phys. 17, 1011 (1976). [75] K. M. C and M. K, A discrete version of the inverse scattering problem, J. Math. Phys. 14, 594 (1973). [76] K. M. C, On discrete inverse scattering problems. II, J. Math. Phys. 14, 916 (1973). [77] M. S, Quantum deformations of the discrete nonlinear Schr¨ odinger equation, Phys. Rev. A 46, 6856 (1992). [78] Z. D. L, P. B. H, L. L, J. Q. L, and W. M. L, Magnetic soliton and soliton collisions of spinor Bose-Einstein condensates in an optical lattice, Phys. Rev. A 71, 053611 (2005). [79] T. C and S. A, Spatially inhomogeneous time-periodic propagating waves in anharmonic systems, Phys. Rev. B 55, R11929 (1997). [80] S. F and K. K, Moving discrete breathers?, Physica D 127, 61 (1999). [81] J. L. M´ ı  S. A, Breathers in nonlinear lattices: numerical calculation from the anticontinuous limit, Nonlinearity 9, 1501 (1996). [82] G. S,Linear Algebra, Academic Press, New York, 1980. [83] H. P, S. A. T, W. T. V, and B. P. F,Numerical Recipes, Cambridge University Press, New York, 1992. [84] J. L. M´ ı, S. A, and L. M. F´ ı, Intrinsic localized modes: Discrete breathers. existence and linear stability, Physica D 113, 283 (1998). [85] Y. S. K and M. P, Modulational instabilities in discrete lattices, Phys. Rev. A 46, 3198 (1992). [86] Y. S. K and M. S, Modulational instabilities in the discrete deformable nonlinear Schr¨ odinger equation, Phys. Rev. E 49, 3543 (1994). [87] M. J, S. A, Y. B. G, P. L. C, and K. O. . R, Dynamics of breathers in discrete nonlinear Schr¨ odinger models, Physica D 119, 115 (1998). [88] P. G. K, K. O. R, and A. R. B, The discrete nonlinear Schr¨ odinger equation: A survey of recent results, Int. Journal of Modern Physics B 15, 2833 (2001). [89] D. H, N. G. S, H. G, and G. P. T, Spatial properties of integrable and nonintegrable discrete nonlinear Schr¨ odinger equations, Phys. Rev. E 52, 255 (1995). [90] D. H and G. P. T, Wave transmission in nonlinear lattices, Phys. Rep. 307, 333 (1999). [91] J. L. M´ ı  F. F  P. J. M´ ı  L. M. F´ ı, Discrete breathers in dissipative lattices, Phys. Rev. E 63, 66603 (2001). 256 BIBLIOGRAPHY [92] P. J. M´ ı  M. M  L. M. F´ ı  F. F, Dissipative discrete breathers: Periodic, quasiperiodic, chaotic, and mobile, Chaos 13, 610 (2003). [93] D. Z, P. J. M´ ı, L. M. F´ ı, and F. F, Mode-locking of mobile discrete breathers, Phys. Rev. E 71, 036613 (2005). [94] S. V. D, P. G. K, B. A. M, and D. J. F, Two-soliton collisions in a near-integrable lattice system, Phys. Rev. E 68, 56603 (2003). [95] D. B. D, J. C. E, H. F, and J. A. D. W, Solitons on lattices, Physica D 68, 1 (1993). [96] S. F, Y. Z, and K. K, Moving lattice kinks and pulses: An inverse method, Phys. Rev. E 59, 6105 (1999). [97] J. L. M´ ı and S. A, Finite size effects on instabilities of discrete breathers, Physica D119, 163 (1998). [98] K. K, Perturbative study of classical Ablowitz-Ladik type soliton dynamics in relation to energy transport in α-helical proteins, Phys. Rev. E 61, 5839 (2000). [99] Y. S. K and D. K. C, Peierls-nabarro potential barrier for highly localized nonlinear modes, Phys. Rev. E 48, 3077 (1993). [100] R. B and M. P, Discreteness effects on a sine-Gordon breather, Phys. Rev. B 43, 8491 (1991). [101] T. D and M. P,Physics of Solitons, Cambridge University Press, (Cambridge, UK), 2006. [102] K. I. P, D. I. P, and I. V. T, Self-action of light beams in nonlinear media: soliton solutions, Opt. Quant. Electr. 11, 471 (1979). [103] M. G. V and A. A. K,Sov. J. Radiophys. Quant. Electr. 16, 783 (1973). [104] I. M. M, B. V. G, R. D, and B. A. M, Finite-band solitons in the Kronig-Penney model with the cubic-quintic nonlinearity, Phys. Rev. E 71, 016613 (2005). [105] T. K, P. G. K, and B. A. M, Stability of multiple pulses in discrete systems, Phys. Rev. E 63, 036604 (2001). [106] D. E. P, P. G. K, and D. J. F, Stability of discrete solitons in nonlinear Schr¨ odinger lattices, Physica D 212, 1 (2005). [107] M. J. A, Z. H. M, and G. B, Methods for discrete solitons in nonlinear lattices, Phys. Rev. E 65, 026602 (2002). [108] H. F  (M. R  M. P .), Nonlinear coherent structures in physics and biology, Springer, (Berlin, Ge), 1991. [109] Y. S and G. J, Numerical computation of travelling breathers in Klein-Gordon chains, Physica D 204, 15 (2005). [110] G. J and Y. S, Travelling breathers with exponentially small tails in a chain of nonlinear oscillators, Commun. Math. Phys. 257, 51 (2005). BIBLIOGRAPHY 257 [111] B. S´ -B and M. J, Exact numerical solutions for dark waves on the discrete nonlinear Schr¨ odinger equation, Phys. Rev. E 71, 036627 (2005). [112] T. R. O. M  A. R. C  P. G. K  J. C, Radiationless Traveling Waves in Saturable Nonlinear Schr¨ odinger Lattices, Phys. Rev. Lett. 97, 124101 (2006). [113] D. C  A. R. B  N. G-J  M. S, Electric-field-induced nonlinear Bloch oscillations and dynamical localization, Phys. Rev. Lett. 74, 1186 (1995). [114] D. K. C, Nonlinear physics: Fresh breather, Nature 432, 455 (2004). [115] P. G. K, K. O. R, and A. R. B, Two-dimensional discrete breathers: Construction, stability, and bifurcations, Phys. Rev. E 61, 2006 (2000). [116] P. G. K, K. O. R, and A. R. B, Localized excitations and their thresholds, Phys. Rev. E 61, 4652 (2000). [117] B. A. M and P. G. K, Discrete vortex solitons, Phys. Rev. E 64, 026601 (2001). [118] S. F, K. K, and R. S. MK, Energy thresholds for discrete breathers in one-, two-, and three-dimensional lattices, Phys. Rev. Lett. 78, 1207 (1997). [119] M. I. W, Excitation thresholds for nonlinear localized modes on lattices, Nonlinearity 12, 673 (1999). [120] M. K, Energy thresholds for discrete breathers, Phys. Rev. Lett. 92, 104301 (2004). [121] M. K, Dimension dependent energy thresholds for discrete breathers, Nonlinearity 17, 1923 (2004). [122] G. K, K. O. R, and A. R. B, Delocalizing transition of Bose-Einstein condensates in optical lattices, Phys. Rev. Lett. 89, 030402 (2002). [123] V. K. M, S. L. M, I. V. R, and S. K. T, Two-dimensional solitons in discrete systems, JETP Lett. 60, 829 (1994). [124] E. W. L, K. H. S, and S. K. T, Stability of discrete solitons and quasicollapse to intrinsically localized modes, Phys. Rev. Lett. 73, 1055 (1994). [125] E. W. L, K. H. S, V. K. M, S. L. M, I. V. R, and S. K. T, Instability of two-dimensional solitons in discrete systems, JETP Lett. 62, 677 (1995). [126] J. J. R and K. R, Blow-up in nonlinear Schr¨ odinger equations I: A general review, Phys. Scr. 33, 481 (1986). [127] P. L. C, Y. B. G, and K. O. R, Dynamics in discrete twodimensional nonlinear Schr¨ odinger equations in the presence of point defects, Phys. Rev. B 54, 900 (1996). [128] P. L. C, Y. B. G, V. K. M, S. L. M, K. O. R, J. J. R, I. V. R, and S. K. T, Discrete localized states and localization dynamics in discrete nonlinear Schr¨ odinger equations, Phys. Scr. 67, 160 (1996). 258 BIBLIOGRAPHY [129] N. K. E, S. S, D. N. C, J. W. F, and M. S, Discrete solitons in photorefractive optically induced photonic lattices, Phys. Rev. E 66, 046602 (2002). [130] J. F, T. C, M. S, N. K. E, and D. N. C, Observation of discrete solitons in optically induced real time waveguide arrays, PRL 90, 023902 (2003). [131] J. C. E  M. J  (L. V´ , R. S. MK, M. P. Z .), Proceedings of Localization and Energy Transfer in Nonlinear Systems, World Scientific, Singapore, 2003. [132] D. C, S. B-A, R. M, J. S. A, H. S. E, Y. S, and D. R, Stable soliton vs collapse dynamics of short laser pulses in nonlinear structures with intermediate dimensionality, HAIT J. Sci. Eng. 1, 363 (2004). [133] T. P  F. L  (D. B. D  J. C. E .), Proceedings of the Confererence on Nonlinear Coherent Structures in Physics and Biology, WWW: http://www.ma.hw.ac.uk/solitons/procs/. [134] P. G. K, D. J. F, R. M. C-G, B. A. M, and A. R. B, Discrete solitons and vortices on anisotropic lattices, Phys. Rev. E 72, 046613 (2005). [135] J. C.   M, Hamiltonian Hopf bifurcation with symmetry, Nonlinearity 3, 1041 (1990). [136] D. N. C and E. D. E, Blocking and routing discrete solitons in twodimensional networks of nonlinear waveguide arrays, Phys. Rev. Lett. 87, 233901 (2001). [137] A. L. B´ , Taming complexity, Nature Physics 1, 68 (2005). [138] B. B,Random Graphs, Academic Press, London, 1985. [139] C. B,The Theory of Graphs, Dover, New York, 2001. [140] P. E ¨  and A. R´ , On random grahs, Pub. Math. Debrecen 6, 290 (1959). [141] P. E ¨  and A. R´ , On the evolution of random graphs, Pub. Math. Inst. Hung. Acad. Sci. 5, 17 (1960). [142] M. E. J. N, The structure of scientific collaboration networks, Proc. Nat. Acad. Sci. 98, 404 (2001). [143] A. L. B´  and R. A, Emergence of scaling in random networks, Science 286, 509 (1999). [144] S. R, How popular is your paper? An empirical study of the citation distribution, Eur. Phys. J. B 4, 131 (1998). [145] H. J, B. T, R. A, Z. N. O, and A. L. B´ , The large-scale organization of metabolic networks., Nature 407, 651 (2000). [146] M. F, P. F, and C. F, On power-law relationships of the Internet topology, Comput. Commun. Rev. 29, 251 (1999). BIBLIOGRAPHY 259 [147] L. A. N. A, A. S, M. B, and H. S,Proc. Nat. Acad. Sci. 97, 11149 (2000). [148] A. V  R. P-S  A. V, Large-scale topological and dynamical properties of the internet, Phys. Rev. E 65, 066130 (2002). [149] M. E. J. N, S. F, and J. B,Phys. Rev. E 66, 035101 (2002). [150] P. S, S. D, A. C, P. A. S, and G. M,Phys. Rev. E (2003). [151] R. G´ , L. D, A. D-G, F. G, and A. A,Phys. Rev. E 68, 065103 (2003). [152] R. G´   S. M  A. T  L. A. N A,Proc. Nat. Acad. Sci. 102, 7794 (2005). [153] R. A  H. J  A. L. B´ , Internet: Diameter of the World-Wide Web, Nature 401, 130 (1999). [154] H. E  L. I. M  S. B, Scale-free topology of e-mail networks, Phys. Rev. E 66, 035103 (2002). [155] D. J. W and S. H. S, Collective dynamics of small-world networks, Nature 393, 440 (1998). [156] S. N. D and J. F. F. M, Evolution of networks, Adv. Phys. 51, 1079 (2002). [157] M. E. J. N, The structure and function of complex networks, SIAM Review 45, 167 (2003). [158] R. P-S and A. V,Evolution and Structure of the Internet, Cambridge University Press, Cambridge, U.K., 2004. [159] S. B  H. G. S (.), Handbook of Graphs and Networks: From the Genome to the Internet, Wiley-VCH, Berlin, 2002. [160] S. N. D and J. F. F. M,Evolution of Networks, Oxford University Press, Oxford, 2003. [161] R. A and A. L. B ´ , Statistical mechanics of complex networks, Rev. Mod. Phys. 74, 47 (2002). [162] S. B, V. L, Y. M, M. C, and D.-U. H, Complex networks: Structure and dynamics, Physics Reports 424, 175 (2006). [163] E. B-N  H. F  Z. T (.), Complex Networks, Springer, (Berlin, Ge), 2004. [164] Provided by Clip2 Distributed Search Solutions. [165] Autonomous System level representation of the Internet as of April 16th 2001. Provided by the National laboratory for applied networ research (NLANR), National Science Foundation, http://moat.nlanlr.net. 260 BIBLIOGRAPHY [166] Router level graph representation of the Internet. Mapping the internet within the scan project, Information Sciences Institute, http://www.isi.edu/div7/scan/. [167] S. M, The small world problem, Psycol. Today 2, 60 (1967). [168] P. S. D, R. M, and D. J. W, An experimental study of search in global social networks, Science 301, 827 (2003). [169] M. E. J. N, Scientific collaboration networks I. Network construction and fundamental results, Phys. Rev. E 64, 016131 (2001). [170] M. R, I. F, and A. I, Mapping the Gnutella network, IEEE Internet Computing 6, 50 (2002). [171] S. W and K. F,Social Networks Analysis, Cambridge University Press, Cambridge, 1994. [172] F. R, C. C, F. C, V. L, and D. P, Defining and identifying communities in networks, Proc. Nat. Acad. Sci. 101, 2658 (2004). [173] S. S-O, R. M, S. M, and U. A, Network motifs in the transcriptional regulation network of Escherichia coli, Nature Genetics 31, 64 (2002). [174] R. M, S. S-O, S. I, N. K, D. C, and U. A, Network motifs: Simple building blocks of complex networks, Science 298, 824 (2002). [175] S. M and U. A, Structure and function of the feed-forward loop network motif, Proc. Nat. Acad. Sci. 100, 11980 (2003). [176] R. M, S. I, N. K, R. L, S. S-O, I. A, M. S, and U. A, Superfamilies of designed and evolved networks, Science 303, 1538 (2004). [177] N. K, S. I, R. M, and U. A, Efficient sampling algorithm for estimating subgraph concentrations and detecting network motifs, Bioinformatics 20, 1746 (2004). [178] J. J. B, N. J. D, A. J. F, and M. E. J. N,The Theory of Critical Phenomena. An Introduction to the Renormalization Group, Clarendon Press, Oxford, UK, 1992. [179] M. E. J. N and D. J. W, Renormalization group analysis of the small-world network model, Phys. Lett. A 263, 341 (1999). [180] M. B and L. A. N. A, Small-world networks: Evidence for a crossover picture, Phys. Rev. Lett. 82, 3180 (1999). [181] A. B and M. W, On the properties of small-world network models, Eur. Phys. J. B 13, 547 (2000). [182] M. E. J. N, C. M, and D. J. W, Mean-field solution of the small-world network model, Phys. Rev. Lett. 84, 3201 (2000). [183] M. E. J. N and D. J. W, Scaling and percolation in the small-world network model, Phys. Rev. E 60, 7332 (1999). BIBLIOGRAPHY 261 [184] P. L. K, S. R, and F. L, Connectivity of growing random networks, Phys. Rev. Lett. 85, 4629 (2000). [185] P. L. K, G. J. R, and S. R, Degree distributions of growing networks, Phys. Rev. Lett. 86, 5401 (2001). [186] S. N. D, J. F. F. M, and A. N. S, Structure of growing networks with preferential linking, Phys. Rev. Lett. 85, 4633 (2000). [187] A. L. B´  and R. A, Mean-field theory for scale-free random networks, Physica A272, 173 (1999). [188] B. B´  and O. R,Preprint. Department of Mathematical Sciences, University of Memphis (2002). [189] R. A and A. L. B ´ , Topology of evolving networks: Local events and universality, Phys. Rev. Lett. 85, 5234 (2000). [190] G. B and A. L. B´ , Competition and multiscaling in evolving networks, Europhys. Lett. 54, 436 (2001). [191] G. B and A. L. B´ , Bose-Einstein condensation in complex networks, Phys. Rev. Lett. 86, 5632 (2001). [192] E. M. J, M. G, and M. E. J. N, Structure of growing social networks, Phys. Rev. E 64, 046132 (2001). [193] P. H and B. J. K, Growing scale-free networks with tunable clustering, Phys. Rev. E 65, 026107 (2002). [194] K. K and V. M. E´ ı, Highly clustered scale-free networks, Phys. Rev. E 65, 036123 (2002). [195] S. N. D, A. V. G, and J. F. F. M, Pseudofractal scale-free web, Phys. Rev. E 65, 06122 (2002). [196] G. S´ , M. A, and J. K´ , Structural transitions in scale-free networks, Phys. Rev. E 67, 056102 (2003). [197] E. R and A. L. B´ , Hierarchical organization in complex networks, Phys. Rev. E67, 026112 (2003). [198] A. V´ , M. B ˜ , Y. M, R. P-S, and A. V, Topology and correlations in structured scale-free networks, Phys. Rev. E 67, 046111 (2003). [199] M. A. S, M. B˜ , and A. D´ ı-G, Competition and adaptation in an internet evolution model, Phys. Rev. Lett. 94, 038701 (2005). [200] S. M, M. B, H. E. S, and L. A. N. A, Truncation of power law behavior in scale-free network models due to information filtering, Phys. Rev. Lett. 88, 138701 (2002). [201] P. L. K and S. R, Organization of growing random networks, Phys. Rev. E 63, 066123 (2001). 262 BIBLIOGRAPHY [202] G. C, A. C, P. D. L. R´ ı, and M. A. M˜ , Scale-free networks from varying vertex intrinsic fitness, Phys. Rev. Lett. 89, 258702 (2002). [203] Z. L, Y. C. L, N. Y, and P. D, Connectivity distribution and attack tolerance of general networks with both preferential and random attachments, Phys. Lett. A 303, 337 (2002). [204] E. F. K, Revisiting ”scale-free” networks, BioEssays 27, 1060 (2005). [205] G. G and R. L, Synchronous neural activity in scale-free network models versus random network models, Proc. Nat. Acad. Sci. 102, 9948 (2005). [206] J. M S,Models in Ecology, Cambridge University Press, Cambridge (UK), 1974. [207] W. O. K and A. G. MK, A contribution to the mathematical theory of epidemics, Proc. Roy. Soc. Lond. A 115, 700 (1927). [208] N. T. J. B,The Mathematical Theory of Infectious Diseases and its Applications, Hafner Press, New York, 1975. [209] R. M. A and R. M. M,Infectious Diseases in Humans, Oxford University Press, Oxford, 1992. [210] J. D. M,Mathematical Biology, Springer, (Berlin, Ge), 1993. [211] O. D and J. H,Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation, Wiley, New York, 2000. [212] R. P-S and A. V, Epidemic spreading in scale-free networks, Phys. Rev. Lett. 86, 3200 (2001). [213] R. M. M and A. L. L, Infection dynamics on scale-free networks, Phys. Rev. E 64, 066112 (2001). [214] R. P-S and A. V, Epidemic dynamics and endemic states in complex networks, Phys. Rev. E 63, 066117 (2001). [215] M. B˜  and R. P-S, Epidemic spreading in correlated complex networks, Phys. Rev. E 66, 047104 (2002). [216] Y. M, R. P-S, and A. V, Epidemic outbreaks in complex heterogeneous networks, Eur. Phys. J. B 26, 521 (2002). [217] R. P-S and A. V, Immunization of complex networks, Phys. Rev. E 65, 036104 (2002). [218] Y. M, J. B. G´ , and A. F. P, Epidemic incidence in correlated complex networks, Phys. Rev. E 68, 035103(R) (2003). [219] R. C, S. H, and D. -A, Efficient immunization strategies for computer networks and populations, Phys. Rev. Lett. 91, 247901 (2003). [220] P. H, Efficient local strategies for vaccination and network attack, Europhys. Lett. 68, 908 (2004). BIBLIOGRAPHY 263 [221] P. E, J. G´ -G˜ , and Y. M, Improved strategies for internet traffic delivery, Phys. Rev. E 71, 035102(R) (2005). [222] J. G´ -G˜ , P. E, and Y. M, Immunization of real complex communication networks, Eur. Phys. J. B 49, 259 (2006). [223] R. M. A and R. M. M, Population biology of infectious diseases: Part i, Nature 280, 361 (1979). [224] P. G, On the critical behavior of the general epidemic process and dynamical percolation, Math. Biosci. 63, 157 (1983). [225] L. M. S, C. P. W, I. M. S, C. S, and J. K, Percolation on heterogeneous networks as a model for epidemics, Math. Biosci. 180, 293 (2002). [226] C.M and M. E. J. N, Epidemics and percolation in small-world networks, Phys. Rev. E 61, 5678 (2000). [227] M. E. J. N, Spread of epidemic diseases on networks, Phys. Rev. E 64, 016128 (2002). [228] H. W. H and J. A. Y, Gonorrhea transmission dynamics and control, Lecture Notes in Biomathematics 56 (1984). [229] M. E. J. N, Ego-centered networks and the ripple effect, Soc. Netw. 25, 83 (2003). [230] M. G and D. J,Computers and Intractability: A Guide to the Theory of NPcompleteness, Freeman, San Francisco, 1979. [231] R. P-S, A. V´ , and A. V, Dynamical and correlation properties of the internet, Phys. Rev. Lett. 87, 258701 (2001). [232] Y. M, M. N, and A. F. P, Dynamics of rumor spreading in complex networks, Phys. Rev. E 69, 066130 (2004). [233] M. T, H. T, and T. S, Critical behaviors and 1/Φnoise in information traffic, Physica A 233, 824 (1996). [234] M. T, H. T, and K. F, Dynamic phase transition observed in the internet traffic flow, Physica A 277, 248 (2000). [235] M. A  M and A. L. B´ , Fluctuations in network dynamics, Phys. Rev. Lett. 92, 028701 (2004). [236] R. V.S´ and S. V, Information transfer and phase transitions in a model of internet traffic, Physica A 289, 595 (2001). [237] S. V and R. V. S´ , Self-organized critical traffic in parallel computer networks, Physica A 312, 636 (2002). [238] S. V and R. V. S´ , Internet’s critical path horizon, Eur. Phys. J. B 38, 245 (2004). [239] T. O and R. S, Phase transition in a computer network traffic model, Phys. Rev. E58, 193 (1998).