The eld of complex systems is a broad scientic area which spans across dierent disciplines, such as mathematics, physics, chemistry, biology, medicine, engineering and many others. One of their most typical traits is the presence of dierent interacting components which may exhibit individually regular or chaotic (unpredictable but not random) behaviour. Such type of interactions may result in phenomena whose evolution is not trivial (i.e. complex) to predict, even when the underlying mechanisms happen to be deterministic. This complexity arises from the type of interactions, one of, which is the nonlinear terms in the dynamical equations of the models which are able to better approximate observable data. Nonlinear dynamical systems provide with theoretical concepts and tools for the analysis of a large variety of problems. For example, the investigation of systems where the motion can be chaotic, pattern formation, synchronization phenomena. Nonlinear dynamical systems are typically divided into conservative and dissipative according to whether or not the total energy of the system or some other relevant quantity is conserved. One example of a powerful tool, from this eld with many applications in dierent systems, is the so-called bifurcation analysis. Such an analysis provides an underlying theoretical explanation and model-based prediction of the qualitative variables' time evolution (periodic, quasi-periodic or chaotic) according to given model parameters (xed or time-dependent). However, when aiming to model systems with dierent components and study their long term evolution, things become more dicult and complex. Their interplay often gives rise to phenomena that can not be predicted simply by the dynamics of each component separately, its initial conditions and parameters, even when that is well understood. Hence, the emerging dynamics of the whole ensemble can be very dierent. Examples and applications of such systems can be found in dynamical astronomy, quantum mechanics, brain dynamics and elsewhere. This habilitation thesis aims to summarize a collection of recent results which span along these three research eld areas. Namely, results in (i) quantum chaotic systems, anomalous transport and localization of quantal eigenfunctions, (ii) classical chaotic systems with application in galactic dynamics and (iii) brain dynamics and synchronization phenomena with applications in medical treatments of certain diseases. In each respective subject, dierent questions and problems are posed as well as adequate tools for tackling them, which originate from the eld of nonlinear dynamical systems.

vii Résumé Le domaine des systèmes complexes est un vaste domaine scientique qui couvre diérentes disciplines, telles que les mathématiques, la physique, la chimie, la biologie, la médecine, l'ingénierie et bien d'autres. L'un de leurs traits les plus typiques est la présence de diérents composants en interaction qui peuvent présenter un comportement individuellement régulier ou chaotique (imprévisible mais pas aléatoire). Un tel type d'interactions peut aboutir à des phénomènes dont l'évolution n'est pas triviale (donc complexe") à prédire, même lorsque les mécanismes sous-jacents se trouvent être déterministes. Cette complexité émerge du type d'interactions, comme c'est le cas des termes non linéaires dans les équations dynamiques des modèles qui sont capables de mieux approcher les données observables. Les systèmes dynamiques non linéaires fournissent des concepts et des outils théoriques pour l'analyse d'une grande variété de problèmes tels que l'étude de systèmes où le mouvement peut être chaotique (imprévisible mais non aléatoire), la formation de motifs et les phénomènes de synchronisation. Les systèmes dynamiques non linéaires sont soit conservatifs soit dissipatifs, selon que l'énergie totale du système ou une autre quantité pertinente est conservée ou non. Un exemple d'outil puissant, issu de ce domaine avec de nombreuses applications dans diérents systèmes, est l'analyse dite de bifurcation. Une telle analyse fournit une explication théorique sous-jacente et une prédiction basée sur un modèle de l'évolution temporelle des variables qualitatives (périodique, quasipériodique ou chaotique) selon les paramètres des modèles donnés (xes ou dépendants du temps).

Cependant, lorsqu'on cherche à modéliser des systèmes avec diérents composants et à étudier leur évolution à long terme, les choses deviennent plus diciles et complexes. Leurs interactions donnent souvent lieu à des phénomènes qui ne peuvent être prédits simplement par la dynamique de chaque composant séparément, ni par ses conditions initiales et ses paramètres, même si cela est bien compris.

Par conséquent, la dynamique émergente de l'ensemble peut être très diérente. Des exemples et des applications de tels systèmes peuvent être trouvés dans l'astronomie dynamique, la mécanique quantique, la dynamique cérébrale et ailleurs. Cette thèse d'habilitation vise à résumer un ensemble de résultats récents qui couvrent ces trois domaines de recherche. À savoir, il en résulte (i) des systèmes chaotiques quantiques, un transport anormal et une localisation de fonctions propres quantiques, (ii) des systèmes chaotiques classiques avec application dans la dynamique galactique et (iii) la dynamique cérébrale et des phénomènes de synchronisation avec des applications dans les traitements médicaux de certaines maladies. Pour chacun de ces sujets, diérentes questions et problèmes sont posés. Pour les aborder, des outils adéquats issus du domaine des systèmes dynamiques non linéaires sont utilisés. viii ix

General introduction Nonlinear Dynamical Systems and complexity

The theoretical background and framework of dynamical systems have been and still are broadly used on a large spectrum of research areas. Applications of dynamical systems, employing dierential equations, are found in mathematics, physics, chemistry, engineering, economics, sociology, medicine etc. Dynamical systems are divided into two main categories, linear and nonlinear according to the underlying process one aims to describe. Linear dynamical systems are to a large extent simpler, well understood and have limited applications. A plethora of phenomena in nature are better and more accurately described by employing nonlinear terms in our models. Such deterministic systems are able to generate complex and often counterintuitive variables' time evolution, like for example chaotic motion (unpredictable but not random). Dynamical systems together with other research elds such as network dynamics, pattern formation, game theory, collective behavior and others, can help us to model and better understand complex interactions in nature and self-organization phenomena, i.e. the so-called complex systems [START_REF] Vemuri | Modeling of Complex Systems: An Introduction[END_REF]].

Nonlinear dynamical systems can be further divided into conservative (Hamiltonian) and dissipative. For the former, the type of evolution (qualitatively) is associated with the number of conserved quantities (e.g. total energy, angular momentum etc.) and the degrees of freedom, while for the latter one can use dissipative quantities (e.g. Lyapunov functions). Furthermore, the bifurcation analysis provides an underlying theoretical explanation and model-based prediction of the qualitative variables' time evolution (periodic, quasi-periodic or chaotic) according to given model (xed or time-dependent) parameters [START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF]].

This habilitation thesis aims to summarize a collection of recent results which span along three quite (at a quick rst glance) dierent and fascinating research eld areas. Nevertheless, there are many links among them when one models their evolution and collective behavior using concepts and properties from nonlinear dynamical systems theory. Hence, we present here results in (i) quantum chaotic systems, anomalous transport and localization of quantal eigenfunctions, (ii) classical chaotic systems with applications in galactic dynamics and (iii) brain dynamics and synchronization phenomena with applications in medical treatments of certain diseases. Evidently, in each respective subject, dierent questions and problems are posed, as well as dierent tools for tackling them.

Setting as a starting point the classical dynamical systems, there is a well-developed framework for the description of dierent natural phenomena, i.e. their evolutionary mechanisms, their potential sensitivity to initial conditions, the transition from one state to another (often associated with some model parameter, e.g. bifurcations), potential model limitations etc. Stellar dynamics is one example where physical models are used for the description of stars' motion (e.g. in galaxies), tested for supporting suciently well observed complex structural patterns as well their stability and much more.

When switching to the study of quantum systems there are evidently major dierences compared to classical systems in many aspects. Yet, they are related to each other. For example, several properties observed in quantum (chaotic) systems, such as localization of eigenfunctions and kinetic energy diusion, are directly associated with their classical corresponding model, at least up to certain time scales. In this framework, understanding the classical models can help in the better description and prediction of the quantum ones. In brain dynamics, and when setting out to model collective (e.g. electrical) activity and dynamics of networks, composed either by point neurons (micro-scale) x or mean activity of larger areas (meso/macro scale), systems of nonlinear ordinary (and partial) dierential equations also provide an adequate tool. There is a large spectrum of mathematical models which are constructed and actively (re-)validated by (new) empirical data for the description for dierent (from a physiological point of view) neuron populations whose distinct spiking patterns aect their collective (healthy or abnormal) activity. The latter is often quantied by the concept of the so-called synchronization which, roughly speaking, measures the degree of the relative rhythmic or repetitive patterns of neural activity. Switching from one rhythmic state to another can be modeled via parameter changes of the model, associated to attractors in the phase space, bifurcations and network analysis.

In this thesis, several concepts of chaos in classical systems are also mentioned and used in the quantum ones, hence, we opt to list the chapter order as follows: in Chapter 1, we start with the discussion classical chaos in astrophysics, then, in Chapter 2, we discuss a few recent ndings on quantum chaos, transport and localization and nally in Chapter 3 we focus on neural (de-)synchronization phenomena in the human brain related to healthy vs pathological activity.

Classical Chaos in Astrophysics

The main concepts of classical chaos lie on the so-called sensitivity to initial conditions in low dimensional systems and their origins go back to Henri Poincaré, the 3body problem, the fact that it is not integrable and the coexistence of quasi-periodic and chaotic dynamics in such Hamiltonian systems. In addition the basic dening property of deterministic chaos, i.e. the divergence of initially nearby phase space trajectories, is quantied by the positivity of the largest Lyapunov exponent [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF]. In astrophysics, one employs nonlinear dynamical systems to model astronomical systems, e.g. celestial bodies, planetary satellites, pulsating stars, and their evolution. To this end, one may consider the gravitational Nbody problem (GNBP), in order to model planetary systems, star clusters of various richness, galaxies and galaxy clusters and super clusters (clusters of clusters). In this framework all types of interactions, except the gravitational one, are not taken into account, as the body distances are much larger than their typical sizes and thereby are not aecting their dynamics. A special case of the GNBP, is the N = 2 case which is equivalent to the onebody (Kepler) problem. In such models, the body evolution is determined only by their gravitational self-interactions and the relative motion and patterns are self-organized. An alternative approach to the GNBP, is the use of Hamiltonian systems, namely dynamical systems that follow Hamilton's equations with conservation properties linked to symmetries of the Hamiltonian function which describes the system and its phase space (periodic, quasi-periodic and chaotic) dynamics. The Hénon & Heiles model [START_REF] Hénon | The Applicability of the Third Integral of Motion: Some Numerical Experiments[END_REF] is considered one of the rst such Hamiltonian systems, exhibiting deterministic chaos which considers point mass stars in a mean eld model approximation of an axially symmetric galaxy.

Orbits of stars are the fundamental building blocks of any galactic structure and their properties give important insight for understanding the formation and evolution of such structures [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF]. Our understanding relies signicantly on the adequacy and eciency of the models used either in the time-dependent (TD) self-consistent models or in the rather `simpler' analytical time-independent (TI) ones. It is by now well accepted that the chaotic or regular nature of orbits inuences the general stability of the N -body simulations, which is straightforwardly related to the underlying dynamics. Therefore, studying the general stability and the detailed structure of the phase (but also of the conguration) space of analytical models can be proven to be very useful, xi provided that the gravitational potentials are realistic in terms of representing density distribution proles close to those derived by simulations. In this Chapter, we summarize some recent results on TD analytical models, their, in general, good agreement with N -body simulations, the presence and evolution of periodic orbits associated with the main galaxy features like the barred shape, the distribution of regular and chaotic motion in the phase as well as conguration space.

Quantum Chaos, Diusion and Localization

Quantum chaos or wave chaos is a more recent branch of chaotic dynamical systems and its main goal is to study the properties of quantum systems which, in their classical limit, exhibit classical (deterministic) chaotic dynamics [START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF][START_REF] Haake | Quantum Signatures of Chaos[END_REF]. Applications of quantum chaos can be found in many areas, such as nuclear, atomic and molecular physics, quantum transport, mesoscopic solid-state systems, wave propagation, acoustics, quantum computers etc.

With the term quantum motion, one refers to the time evolution of the wavefunctions. Even though quantum motion of bound systems with purely discrete energy spectrum is ultimately (after suciently long time) stable and regular (almost periodic), it exhibits many features of the classical motion such as for example diusion in a chaotic domain, for time up to the Heisenberg time. The Heisenberg time, also called break time, is an important time scale in any quantum system, and is given by t H = 2π /∆E, where h = 2π is the Planck constant and ∆E is the mean energy level spacing, such that the mean energy level density is ρ(E) = 1/∆E. For time shorter than approximately t H the quantum diusion follows the classical chaotic diusion, but is stopped at larger times, just due to the interference phenomena (similar to Anderson localization occurring for particles in disordered solids). Pictorially speaking, for time up to t H the quantum system behaves as if its evolution operator has a continuous spectrum, like the classical one has in the chaotic regime, but at later times it senses the discreteness of the spectrum. If the quantum diusion stops, while the classical chaotic diusion continues, we speak about the dynamical localization, or quantum localization or Chirikov localization, rst observed in time-dependent systems [START_REF] Casati | Stochastic behavior of a quantum pendulum under a periodic perturbation[END_REF]]. The problem of quantum or dynamical localization is related to the Anderson localization model, within the framework of the tight-binding approximation, with hopping transitions between the nearest neighbors [START_REF] Fishman | Chaos, Quantum Recurrences, and Anderson Localization[END_REF]. We here summarize some recent results on dynamical localization in time-dependent periodic (Floquet) systems and diusion properties, exemplied by the quantum kicked rotator.

Modeling Neurons and the Dynamics of the Brain

Biological systems is another eld where nonlinear dynamical systems have also been used for the description and understanding of complex collective activity. One of the most important applications is in the dynamics of nerve membranes via the so-called action potential. Such a biological activity resembles the one used for the dynamics of electric circuits and ideas from the eld of nonlinear dynamical systems are also used here, associating the type of activity of individual cells/neurons or the interaction as ensembles (networks) with bifurcations of certain parameters resulting in periodic, quasi-periodic or chaotic activity.

Alan Hodgkin and Andrew Huxley pioneered the eld of dynamical neuroscience when in 1952 introduced a mathematical model of nonlinear dierential equations to describe and explain the mechanisms of the electrophysiological activity of squid giant axons (i.e. nerve membranes), based xii on their experimental data. When a squid giant axon is not sending any signal is said to be at rest. One can modulate the dynamical activity of such nerve by modifying external ionic (such as the calcium) concentrations around the axon and induce self-sustained oscillations [START_REF] Huxley | Ion Movements During Nerve Activity[END_REF][START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF], a state at which the nerve membrane acts as a neural oscillator. Linking these behaviors to dynamical systems' terminology, the resting state corresponds to a stable equilibrium point while the oscillatory to a stable limit cycle. Non-periodic (chaotic) oscillations can also be observed [START_REF] Izhikevich | Dynamical systems in neuroscience: The geometry of excitability and bursting[END_REF]]. Transitions from one state to another take place via bifurcations associated to model parameters and/or external electrical input.

When aiming to understand the collective behavior of neural systems (ensembles of neurons), besides the individual's neuron activity pattern, the concept of synchronization becomes rather relevant.

Namely, the adjustment of rhythmic activity of self-sustained periodic (neural) oscillators due to the type of (network) connectivity and interaction they have which, in a more physiologically realistic description, can be time dependent (plasticity). Such a collective activity is usually associated to the dynamics of the neurons' phases (neurons can be modeled as oscillators), their relative phase dierence which may result to phase locking, frequency entrainment or asynchronous activity [A. [START_REF] Pikovsky | Synchronization: A universal concept in nonlinear sciences[END_REF][START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF]].

Neuronal synchronization processes are relevant under normal as well as abnormal conditions. A number of brain disorders are associated with abnormal neuronal synchrony, for example Parkinson's disease (see e.g. [START_REF] Hammond | Pathological synchronization in Parkinson's disease: networks, models and treatments[END_REF]), tinnitus (see e.g. [START_REF] Elgoyhen | Tinnitus: perspectives from human neuroimaging[END_REF]) and epilepsy (see e.g. [START_REF] Wong | Cellular basis of neuronal synchrony in epilepsy[END_REF]). To specically counteract abnormal neuronal synchrony and, hence, related symptoms, a number of invasive and non-invasive techniques have been developed over the last years. Coordinated Reset (CR) stimulation is a rather recent protocol which employs basic plasticity and dynamic self-organization principles of the nervous system. Its fundamental goal is to induce long-lasting desynchronizing eects that persist cessation of stimulation. The latter are key to reducing side eects of invasive therapies such as the standard high-frequency (> 100 Hz) deep brain stimulation [START_REF] Deuschl | A Randomized Trial of Deep-Brain Stimulation for Parkinson's Disease[END_REF]] which has only acute clinical and electrophysiological eects, present only during stimulation as well as a number of adverse events, such gait disturbances and speech problems. We here summarize some recent results on ecient and optimized CR stimulation signal patterns which induce long-lasting desynchronization of abnormal activity that may occur in Parkinson's disease and tinnitus.

The thesis consists of three Chapters based largely on published papers. Each Chapter is dedicated to one of the above described topics with a separate and extensive introduction of the topic's background and general concepts, the models' presentation, a few representative results from the related published papers and a summary accompanied with perspectives for future investigation. Each Chapter is selfcontained in order to give the reader a good overview of the published work in the respective eleven +D=FJAH Classical chaos in astrophysics 1.1 Background and introduction

Orbits of starts are generally regarded as the backbone of structure in galaxies. Exploring orbital properties in general and in particular the evolution of their dynamical stability is a fundamental aspect in improving our understanding of galactic structures as a whole [START_REF] Binney | Galactic Dynamics: Second Edition[END_REF].

Our ability to explore the details of orbital stability in galaxies depends considerably on the adequacy of analytical models, which can be time-independent (TI) or time-dependent (TD). Studying the stability and the phase space structure via analytical models (see e.g. [START_REF] Manos | Regular and chaotic orbits in barred galaxies -I. Applying the SALI/GALI method to explore their distribution in several models[END_REF]) has proven to be quite useful (for a review, see [START_REF] Contopoulos | Order and Chaos in Dynamical Astronomy[END_REF]), as long as those potentials are realistic, in the sense of adequately representing the density distributions of real galaxies. It is well accepted that the chaotic or regular nature of orbits inuences the general stability of the N -body simulations, which is straightforwardly related to the underlying dynamics (self-consistent, i.e. the motions are caused from the forces generated by the masses themselves). Hence, studying the general stability and the detailed structure of the phase and conguration space of analytical models can be proven to be very useful, provided that these potentials are realistic in terms of representing density distribution proles close to those derived by simulations or extracted from observational data.

The general nature of an orbit, in conservative (Hamiltonian) TI-systems, can only be one of the following: periodic (stable or unstable ), quasi-periodic or chaotic [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF].

Nevertheless, there are cases where chaos can be characterized as weak, suggesting that orbits spend a signicant fraction of their time in conned regimes and do not ll up phase space as `homogeneously'

as the strongly chaotic ones. In these cases, the dierent rate of diusion in the phase space plays an important role, associating for example the weak chaotic motions with barred or spiral galaxy features, giving rise to a number of interesting results, see e.g. [START_REF] Athanassoula | Rings and spirals in barred galaxies -II. Ring and spiral morphology[END_REF]Athanassoula, Romero-Gómez, and Masdemont 2009;[START_REF] Athanassoula | Rings and spirals in barred galaxies -II. Ring and spiral morphology[END_REF][START_REF] Harsoula | Orbital structure in N-body models of barred-spiral galaxies[END_REF]Harsoula, Kalapotharakos, and Contopoulos 2011a;Harsoula, Kalapotharakos, and Contopoulos 2011b;[START_REF] Contopoulos | Stickiness in Chaos[END_REF][START_REF] Kaufmann | Self-consistent models of barred spiral galaxies[END_REF][START_REF] Patsis | Orbits in the Bar of NGC 4314[END_REF][START_REF] Patsis | The stellar dynamics of spiral arms in barred spiral galaxies[END_REF][START_REF] Romero-Gómez | The origin of rR1 ring structures in barred galaxies[END_REF][START_REF] Romero-Gómez | The formation of spiral arms and rings in barred galaxies[END_REF]; T. [START_REF] Bountis | Complex statistics in Hamiltonian barred galaxy models[END_REF]. There are also several results in the recent literature showing that strong local instability does not necessarily imply widespread diusion in phase space [START_REF] Cachucho | Chirikov diusion in the asteroidal threebody resonance[END_REF][START_REF] Giordano | Chaotic diusion of orbits in systems with divided phase space[END_REF]. In [START_REF] Contopoulos | Stickiness in Chaos[END_REF][START_REF] Contopoulos | Stickiness in Chaos[END_REF]] `stickiness' was studied thoroughly in 2-degrees of freedom (d.o.f.) while in [START_REF] Katsanikas | Chains of Rotational Tori and Filamentary Structures Close to High Multiplicity Periodic Orbits in a 3d Galactic Potential[END_REF]Katsanikas and Patsis 2011;[START_REF] Katsanikas | The Structure and Evolution of Conned Tori Near a Hamiltonian Hopf Bifurcation[END_REF] and in [START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF] the role of `sticky' chaotic orbits and the diusive behavior, in the neighborhood of invariant tori surrounding periodic solutions of the Hamiltonian in the vicinity of periodic orbits in conservative systems, was also studied.

Lyapunov exponents (see e.g. [C. [START_REF] Skokos | The Lyapunov Characteristic Exponents and Their Computation[END_REF] for a review and references therein) have been extensively used for the detection of chaotic motion in several dierent models. However, there are often disadvantages which hinder their use, for example their slow convergence. Several approaches exist to detect and quantify chaos, whose dierences and ecacy have been thoroughly compared and discussed in the recent literature (see [START_REF] Contopoulos | Order and Chaos in Dynamical Astronomy[END_REF][START_REF] Maone | A comparison of dierent indicators of chaos based on the deviation vectors: application to symplectic mappings[END_REF][START_REF] Maone | A comparison of dierent indicators of chaos based on the deviation vectors: application to symplectic mappings[END_REF] and references therein). In [Ch. Skokos, Gottwald, and Laskar 2016] the reader may nd a special volume with a more complete and recent review of the several broadly used chaos detection methods and their predictability as well as all the relevant references regarding their theoretical background, numerical implementation and applications in various models. We use the Generalized Alignment and (b) a case where the bar gets weaker by losing mass (see e.g. [Combes 2008a;Combes 2008b]).

Index (GALI) method [C.
There, a reliable way of using the GALI chaos detection method was used for estimating the relative fraction of chaotic vs. regular orbits in such TD potentials. We stress here that in the TD models, individual trajectories may display sudden transitions from regular to chaotic behavior and vice versa during their time evolution and in general the `sticky' behaviour, as discussed in the literature, is less pronounced. This is also the typical case in the N -body simulations where, generally speaking, the motion may also be either: (i) regular throughout the whole evolution, (ii) chaotic throughout the whole evolution, (iii) alternate between chaotic and regular motion with simultaneously orbital shape change (but not necessarily), e.g. from disc to bar like, etc.

Regarding the galaxies' evolution and formation of their several features, it is generally accepted that the most appropriate way to study them is by analyzing N -body simulations. The self-consistency of the models in this approach captures much better several details of the general dynamics. The direct application of chaos detection methods to individual orbits is still a rather dicult task while, for a large ensemble of particles, it is even harder if not unfeasible. To overcome this obstacle, mean eld potentials have been used in the literature in order to study in more detail the dynamical properties of a specic N -body simulation. These potentials are referred to as `frozen' are TI and are derived at specic snapshots of the simulations. Hence, one can apply chaos detection tools to the mean eld potential instead of the N -body simulation. For example, [START_REF] Muzzio | Spatial Structure of Regular and Chaotic Orbits in A Self-Consistent Triaxial Stellar System[END_REF] used an elliptical galaxy simulation (no bar or halo) without dissipation which collapses and eventually reaches an equilibrium state. Then, by taking a quadrupolar expansion of the frozen snapshot, they derive a stationary smooth potential. In [START_REF] Voglis | Chaotic motion and spiral structure in self-consistent models of rotating galaxies[END_REF] and references therein, the authors deal with disc galaxies, focusing mainly on the spiral structures rather than bars (no halo) while the extraction of the mean eld potential is again performed in a similar manner. Following this approach, the role of chaotic motion and diusion rate in barred spiral galaxies has also been studied [START_REF] Harsoula | Orbital structure in N-body models of barred-spiral galaxies[END_REF]Harsoula, Kalapotharakos, and Contopoulos 2011a;Harsoula, Kalapotharakos, and Contopoulos 2011b;[START_REF] Maone | A comparison of dierent indicators of chaos based on the deviation vectors: application to symplectic mappings[END_REF][START_REF] Contopoulos | Stickiness in Chaos[END_REF] while some applications to the Milky Way bar can be found in [Y. [START_REF] Wang | A new model for the Milky Way bar[END_REF]] and recently a new code for orbit analysis and Schwarzschild modelling of triaxial stellar systems was given in [START_REF] Vasiliev | A new code for orbit analysis and Schwarzschild modelling of triaxial stellar systems[END_REF]]. Nevertheless, following this approach one only derives a stationary mean eld model for an equilibrium state of the simulation under study. Furthermore, it does not incorporate an appropriate type set of parameters that would be able to describe and reproduce the time-dependencies in axis ratios, masses, pattern speed, etc. of the several components of a model, like for example the growth of the bar component or the evolution of the disc in time. Let us point out here, the fact that in all these approaches the orbits under study can be only either regular or chaotic. The latter ones may be further distinguished to strongly or weakly chaotic, depending on their diusion properties, sticky eects etc... during the whole evolution.

We here consider an N -body simulation of a disc galaxy embedded in a live halo (i.e. both the stellar disc and the dark matter halo are represented by responsive particles). Disc-halo interaction leads to the formation of a strong bar. We then measure how the galaxy components vary in time during the simulation. The time evolution of the structural parameters is provided as input to the analytical TD model we build. This `candidate' analytical mean eld potential is meant to mimic the N -body simulation evolution and more importantly to generate orbits with more similar (and in some sense `richer') morphological behaviour to those of the N -body simulation, i.e., permitting for individual orbits the interplay between regular and chaotic epochs as time evolves and providing a stable structure at the same time. Note that, in TI frozen models an orbit cannot convert from chaotic to regular. Our TD model is composed of three components (bar, disc and halo) whose parameters were tted with the N -body measurements, via the rotation curves. Note that many simplifying assumptions are made. For example, our TD model considers an (ellipsoidal) analytical bar component which is not always an excellent approximation of the shape of the actual N -body bar. Likewise, the analytical description of the halo and the disc cannot be expected to behave identically, either. However, our goal is to study the general dynamical impact in stability caused by the bar's growth in time (as it happens in the N -body simulation). Thus, by using a realistic TD model, without aiming to describe of the exact detailed dynamics yielding from the simulation, we can use chaos detection tools and quantify general trends of the fraction of regular and chaotic orbits in the phase and conguration space. Keeping this in mind, we draw (disc) initial conditions directly from the simulation and we evolve them in time with the mean eld TD potential.

The N -body and time-dependent analytical models

The N -body simulation

To serve as the base reference for the analytical model, we use one of the simulations described in [START_REF] Machado | Loss of halo triaxiality due to bar formation[END_REF]. For simplicity, we select initial conditions with a spherical halo.

The mass of the stellar disc is M d = 5 × 10 10 M ⊙ , with an exponential density prole of radial scale length R d = 3.5 kpc, and vertical scale height z 0 = 0.7 kpc:

ρ d (R, z) = M d 4πz 0 R 2 d exp ( - R R d ) sech 2 ( z z 0 ) , (1.1)
The spherical dark matter halo has a [START_REF] Hernquist | N-body realizations of compound galaxies[END_REF]] density prole and it is ve times more massive than the disc:

ρ h (r) = M h 2π 3/2 α r c exp (-r 2 /r 2 c ) r 2 + γ ′2 , (1.2)
where M h = 2.5 × 10 11 M ⊙ is the mass of the halo, γ ′ = 1.7 kpc is a core radius and r c = 35 kpc is a cuto radius. The normalisation constant α is dened by

α = {1 - √ πq exp (q 2 )[1 -erf(q)]} -1 (1.3)
where q = γ ′ /r c . For additional details on the initial conditions, see [START_REF] Machado | Loss of halo triaxiality due to bar formation[END_REF].

This is a fairly representative collisionless simulation of a strongly barred galaxy. Four snapshots of the disc particles are displayed in the upper row of Fig. 1.1. It was performed with the N -body code gyrfalcon [START_REF] Dehnen | A Very Fast and Momentum-conserving Tree Code[END_REF][START_REF] Dehnen | A Very Fast and Momentum-conserving Tree Code[END_REF]] using a total of 1.2 million equal-mass particles, with a gravitational softening length of 0.175 kpc, resulting in 0.1 per cent energy conservation. The simulation was carried out for approximately one Hubble time.

N -body simulations have been employed to study chaotic motion in simplied models of disc galaxies. For example, [START_REF] Voglis | Chaotic motion and spiral structure in self-consistent models of rotating galaxies[END_REF]] study chaos and spiral structure in rotating disc galaxies, but those galaxies are not embedded in dark matter haloes.

The connection between chaos and bars was also analysed by [START_REF] El-Zant | Dark Halo Shapes and the Fate of Stellar Bars[END_REF], with models were set up by the addition of disc, bar and halo components. They found that in centrally concentrated models, even a mildly triaxial halo lead to the onset of chaos and the dissolution of the bar in a timescale shorter than the Hubble time.

The time-dependent analytical model

We construct an analytical model that is described by its total gravitational potential

V = V B (t) + V D (t) + V H (t)
, where the three components correspond to the potentials of the bar, disc, and halo, respectively. These components will evolve in time, in accordance with the behaviour we measure from the simulation (see [Manos and Machado 2014] for a detailed motivation and justication).

(i) A triaxial Ferrers bar [START_REF] Ferrers | On the Potentials, Ellipsoids, Ellipsoidal Shells, Elliptic Laminae and Elliptic Rings, of Variable Densities[END_REF]], whose density is given by:

ρ(x, y, z) = { ρ c (1 -m 2 ) 2 if m < 1, 0 if m ≥ 1, (1.4) where ρ c = 105 32π GM B (t) abc is the central density, M B (t)
is the mass of the bar, which changes in time, and

m 2 = x 2 a 2 + y 2 b 2 + z 2 c 2 , a > b > c > 0,
V B (t) = -πGabc ρ c 3 ∫ ∞ λ du ∆(u) (1 -m 2 (u)) 3 , (1.5)
where G is the gravitational constant (set to unity),

m 2 (u) = x 2 a 2 +u + y 2 b 2 +u + z 2 c 2 +u , ∆ 2 (u) = (a 2 + u)(b 2 + u)(c 2 + u)
, and λ is the unique positive solution of m 2 (λ) = 1, outside of the bar (m ≥ 1), while λ = 0 inside the bar. The analytical expression of the corresponding forces are given in [START_REF] Pfenniger | The 3D dynamics of barred galaxies[END_REF]]. In our model, the shape parameters (i.e. the lengths of the ellipsoid axes a, b and c are) are also functions of time.

(ii) A disc, represented by the MiyamotoNagai potential [START_REF] Miyamoto | Three-dimensional models for the distribution of mass in galaxies[END_REF]: .6) where A and B are its horizontal and vertical scale-lengths, and M D (t) is the mass of the disc.

V D (t) = - GM D (t) √ x 2 + y 2 + (A + √ z 2 + B 2 ) 2 , ( 1
Here, `disc mass', refers to the stellar mass excluding the bar. As the bar grows, its mass increases at the expense of the remainder of the disc mass, such that the total stellar mass is constant: M B (t) + M D (t) = 5 × 10 10 M ⊙ . The parameters A and B are also functions of time.

(iii) A spherical dark matter halo, represented by a Dehnen potential [START_REF] Dehnen | A Family of Potential-Density Pairs for Spherical Galaxies and Bulges[END_REF]]:

V H (t) = GM H a H ×    -1 2-γ [ 1 - ( r r+a H ) 2-γ ] , γ = 2, ln r r+a H , γ = 2. (1.7)
M H is the halo mass, a H is a scale radius and the dimensionless parameter γ (within 0 ≤ γ < 3) governs the inner slope. The halo mass is constant throughout, but the parameters a H and γ are functions of time. For γ < 2 its nite central value is equal to

(2 -γ) -1 GM H /a H .
Instead of attempting to use the (disc and halo) proles from the N -body simulations, we opted to represent the bar, disc and halo using respectively the Ferrers, MiyamotoNagai and Dehnen proles. There are two reasons for such a choice. First, our approach requires analytical simplicity that could not be aorded by the proles used in the initial conditions of the numerical simulation.

Secondly, due to bar formation and evolution, the initial disc prole in the simulation soon becomes a poor representation in the inner part of the galaxy, where the bar resides. In this sense, it is not advantageous to continue using the initial proles to model later times. A MiyamotoNagai disc provides a sucient approximation for our purposes. Likewise, even though the [START_REF] Hernquist | N-body realizations of compound galaxies[END_REF] halo prole is well suited for numerical purposes, it is inconvenient from the analytical point of view. We experimented with simple logarithmic halo proles (because their rotation curves are also appropriate), but the [START_REF] Dehnen | A Family of Potential-Density Pairs for Spherical Galaxies and Bulges[END_REF]] prole was preferable, as it provided equally acceptable rotation curves, and a more satisfactory global approximation of the mass distribution. Similarly, tting the bar by a Ferrers ellipsoid is a justiable approximation. Surely, it fails to capture the N -body bar in all its complexity, particularly after the buckling instability, when the bar is substantially strong and develops the peanut-shaped feature. In general, the N -body bar will be more boxy than the Ferrers shape would allow. Nevertheless, tting an ellipsoid of the same extent allows us to obtain plausible shapes, to determine the bar orientation and to estimate its mass adequately. Ultimately, regardless of small deviations in the density proles, our goal is to obtain an analytical total potential that is approximately comparable to the overall potential of the simulation. From the simulation, we are able to measure several quantities as a function of time, which are then used to inform the analytical model.

For the bar, and following [Manos and Machado 2014], the required parameters are the bar mass, the bar shape and the bar pattern speed Ω b . First, we estimate the bar length as a function of time. This is done by measuring the relative contribution of the m = 2 Fourier component of the mass distribution as a function of radius, for each time step, and nding the radius at which the m = 2 has its most intense drop after the peak. This radius a is associated with the bar length. Then we The disc mass M D (t) is known once the bar mass has been measured, and the halo mass M H is constant. One still requires the time evolution of two disc parameters (A, B) and two halo parameters (a H , γ). This is achieved by measuring the rotation curves directly from the simulation (at each time step), and then tting the analytical v c (R) to these data. Since the disc and halo potentials are known from equations (1.6) and (1.7), we obtain their respective analytical circular velocities from Fitting Eq. (1.9) to the measured halo rotation curve, we obtain a H and γ. In the case of the disc, it is not enough to t Eq. (1.8). One must simultaneously t the MiyamotoNagai density prole to disambiguate the A + B (ignoring the inner part of the disc). When tting the disc rotation curve, we assume the total stellar mass (i.e. we take both disc and bar mass into account). Since the circular velocities rely on azimuthally averaged quantities, the presence of the bar does not greatly interfere with the quality of the ts, while its removal would lead to spurious results. The measured rotation curves (disc, halo and total), as well as the resulting tted circular velocities, are displayed in the fourth column of Fig. 1.3 (at four illustrative instants in time). Errors in the tted parameters of rotation curves were typically of about 5 per cent or less. Also shown in the rst, second and third panels of Fig. 1.3, are the disc (radial and vertical) and halo density proles. The points correspond to simulation measurements and the lines give the resulting ts.

v 2 c = R dV dR : v 2 c,D (R) = R 2 GM D [R 2 + (A + B) 2 ] 3/2 (1.8) v 2 c,H (R) = GM H r 2-γ (r + a H ) 3-γ (1.9)
One of the main arguments in favor of the adequacy of our analytical model is evidenced by the fact that its total rotation curves are in good agreement with those measured from the simulation. This indicates that the choices of proles were not unreasonable, as they result in a globally similar gravitational potential. Even if individually the densities of the components are idealized simplications, the similarity of the total potential ensures that the overall dynamical evolution should be suciently well approximated.

Finally, the resulting time evolution of the halo and disc structural parameters, measured in the manner described above, are displayed in the rst to fourth panels of Fig. 1.2. With these, the time-dependence of the analytical model is fully specied. In Table 1 in [Manos and Machado 2014], we summarize the analytical model by showing a sample of parameters for the Ferrers bar, MiyamotoNagai disc and Dehnen halo potential as tted by the N -body simulation, at four times.

Bar strength

In order to measure the bar strengths in analytical models, [START_REF] Manos | Regular and chaotic orbits in barred galaxies -I. Applying the SALI/GALI method to explore their distribution in several models[END_REF] [START_REF] Buta | A Technique for Separating the Gravitational Torques of Bars and Spirals in Disk Galaxies[END_REF][START_REF] Buta | The Distribution of Maximum Relative Gravitational Torques in Disk Galaxies[END_REF], which is a measure of the relative strength of the non-axisymmetric forces. Here, we opt instead to use a method more familiar to N -body simulations, namely measurements of the m = 2 Fourier component of the mass distribution. For the N -body simulation, we measure this component straightforwardly as a function of radius and then take the maximum amplitude to be the A 2 (see e.g. [Athanassoula, Machado, and Rodionov 2013a]. We refer to this quantity as the bar strength.

had employed the Q b parameter [R.
For the analytical model, we proceed in a way that allows us to treat it as if it could be represented by particles. We extract from the simulation a random sample of 100 000 initial conditions (i.e. positions and velocities of disc particles) at a time t 0 = 1.4 Gyr, where we use this ensemble of orbits to further study dynamical trends). The orbits of each of these `test particles' are then evolved forward in time in the presence of our time-dependent analytical potential. Their successive positions can be treated as if they were simulation particles. By stacking them at each time step, we produce the snapshots in the bottom row of Fig. 1.1. These mock snapshots display a striking resemblance to the N -body snapshots, specially bearing in mind that they were obtained by very indirect means.

While this comparison cannot be expected to yield a perfect morphological equivalence, one notices that the bar lengths are in quite good agreement, and that in both cases rings are present (although not of the same extent). The point is that the dynamics that arises from the analytical model will give rise to very similar disc and bar morphologies. In fact, the relative importance of the bar is also quite comparable, as indicated by the A 2 parameter. Analogously to the N -body case, we compute the A 2 of these mock snapshots and compare them in Fig. 1.4. We must stress here that this comparison is an a posteriori verication, i.e. the bar strength of the N -body simulation was in fact not used as an input to the analytical model. The fact the A 2 do agree well counts as a further sign of the consistency of the constructed analytical model.

It is clear, of course, that the variation of the bar strength modies the values of several parameters and yields richer information about the dynamics of a self-consistent model. N -body simulations show that in general, variations of the bar mass also change the mass ratios of the model's components, the bar shape and the pattern speed of the galaxy. Hence, if one wishes to use a mean eld potential to `mimic' a self-consistent model as accurately as possible, one should allow for all the parameters that describe the bar (together with all other axisymmetric components) to depend on time, assuming that the laws of such dependence were explicitly known. In our case, however, we adopt a simpler approach and vary only the masses of the bar and the disc, as a rst step towards investigating such models when time-dependent parameters are taken into account. Thus, we do not pretend to be able to reproduce the exact dynamical evolution of a realistic galactic simulation. Rather, we wish to understand the eects of time dependence on the general features of barred galaxy models and compare the eciency of chaos indicators like the GALI method and the Maximal Lyapunov Exponent (MLE) in helping us unravel the secrets of the dynamics in such problems. We stress that the method we introduced to construct the analytical model does not rely at all on frozen potentials. Instead, it is grounded on the detailed features of a fully time-dependent, self-consistent N -body simulation.

Unless otherwise stated, the units of the analytical model are given as: 1 kpc (length), 1000 km• sec -1 (velocity), 1 Myr (time), 2 × 10 11 M ⊙ (mass) and km• sec -1 • kpc -1 (Ω b ) while the parameter G = 1. The total mass M tot = M B (t) + M D (t) + M H is set equal to 3 × 10 11 M ⊙ and since the halo's mass M H is kept constant, the disc's mass M D (t) is varied as M D (t) = M tot -(M H + M B (t)).

Chaos detection techniques

Let us briey recall how the two main chaos detection methods used throughout the manuscript, namely the GALI and the MLE, are dened and calculated. Considering the following TD 3-d.o.f. Hamiltonian function which determines the motion of a star in a 3 dimensional rotating barred galaxy:

H = 1 2 (p 2 x + p 2 y + p 2 z ) + V (x, y, z, t) -Ω b (t)(xp y -yp x ).
(1.10)

The bar rotates around its zaxis (short axis), while the x direction is along the major axis and the y along the intermediate axis of the bar. The p x , p y and p z are the canonically conjugate momenta, V is the potential, Ω b (t) represents the pattern speed of the bar and H is the total energy of the orbit in the rotating frame of reference (equal to the Jacobi constant in the TI case).

The corresponding equations of motion are: .11) while the equations governing the evolution of a deviation vector w = (δx, δy, δz, δp x , δp y , δp z ) needed for the calculation of the MLE and the GALI, are given by the variational equations:

ẋ = p x + Ω b (t)y, ẏ = p y -Ω b (t)x, ż = p z , ṗx = - ∂V ∂x + Ω b (t)p y , ṗy = - ∂V ∂y -Ω b (t)p x , ṗz = - ∂V ∂z , ( 1
δ x = δp x + Ω b (t)δy, δy = δp y + Ω b (t)δx, δz = δp z , δp x = - ∂ 2 V ∂x∂x δx - ∂ 2 V ∂x∂y δy - ∂ 2 V ∂x∂z δz + Ω b (t)δp y , δp y = - ∂ 2 V ∂y∂x δx - ∂ 2 V ∂y∂y δy - ∂ 2 V ∂y∂z δz -Ω b (t)δp x , δ p z = - ∂ 2 V ∂z∂x δx - ∂ 2 V ∂z∂y δy - ∂ 2 V ∂z∂z δz.
(1.12)

Regarding the estimation of the value of the MLE, λ 1 , of an orbit under study we follow numerically its evolution in time together with its deviation vectors w, by solving the set of Eqs. (1.11) and (1.12) respectively. For this task we use a Runge-Kutta method of order 4 with a suciently small time step, which guarantees the accuracy of our computations, ensuring the relative errors of the Hamiltonian function (in the TI case) are typically smaller than 10 -6 . Furthermore, we need to have a xed time step in order to ensure that in the TD case the orbits vary simultaneously with the potential.

In general, the derivatives of the potential V depend explicitly on time and the ordinary dierential equations (ODEs) [Eqs. (1.11)] are non-autonomous. Hence, one has to solve together the equations for the deviation vectors [Eqs. (1.12)] with the equations of motion [Eqs. (1.11)]. Transforming the Eqs. (1.11), and consequently [Eqs. (1.12)], to an equivalent autonomous system of ODEs by considering time t as an additional coordinate (see e.g. section 1.2b [START_REF] Lichtenberg | Regular and Chaotic Dynamics[END_REF]), is not particularly helpful, and is better to be avoided as shown in [START_REF] Grygiel | Lyapunov Exponents Analysis of Autonomous and Nonautonomous Sets of Ordinary Dierential Equations[END_REF].

So, in order to compute the MLE and the GALI we numerically solve the time-dependent set of ODEs [Eq,. (1.11) and (1.12)]. Then, according to [START_REF] Benettin | Kolmogorov entropy and numerical experiments[END_REF][START_REF] Contopoulos | On the number of isolating integrals in Hamiltonian systems[END_REF][START_REF] Benettin | Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems -A method for computing all of them. I -Theory. II -Numerical application[END_REF] the MLE λ 1 is dened as: 1.13) where: 0) , (1.14) is the so-called `nite time MLE', with w(0) and w(t) being the Euclidean norm of the deviation vector at times t = 0 and t > 0 respectively. A detailed description of the numerical algorithm used for the evaluation of the MLE can be found in [C. [START_REF] Skokos | The Lyapunov Characteristic Exponents and Their Computation[END_REF].

λ 1 = lim t→∞ σ 1 (t), (
σ 1 (t) = 1 t ln w(t) w(
This computation can be used to distinguish between regular and chaotic orbits, since σ 1 (t) tends to zero (following a power law ∝ t -1 ) in the former case, and converges to a positive value in the latter. But the Hamiltonian [Eq. 1.10)] is TD in this case, which means that its orbits could change their dynamical behavior from regular to chaotic and vice versa, over dierent time intervals of their evolution. In such cases, the MLE (1.13) does not behave exactly as in TI model (presenting in general stronger uctuations) and its computation might not be able to identify the various dynamical phases of the orbits, since by denition it characterizes the asymptotic behavior of an orbit (see e.g. [START_REF] Manos | Interplay between chaotic and regular motion in a time-dependent barred galaxy model[END_REF] and also [T. [START_REF] Bountis | Application of the GALI method to localization dynamics in nonlinear systems[END_REF][START_REF] Manos | Scaling with System Size of the Lyapunov Exponents for the Hamiltonian Mean Field Model[END_REF][START_REF] Moges | On the behavior of the Generalized Alignment Index (GALI) method for regular motion in multidimensional Hamiltonian systems[END_REF] for some relevant applications). Nevertheless, we will show the MLE for a number of orbits throughout the paper, for a more global discussion of the several dynamical properties observed.

Thus, in order to avoid such problems in our study, we use the GALI method of chaos detection [C. Skokos, T. C. Bountis, and C. Antonopoulos 2007]. The GALI index of order k (GALI k ) is determined through the evolution of 2 ≤ k ≤ N initially linearly independent deviation vectors w i (0), i = 1, 2, . . . , k, with N denoting the dimensionality of the phase space of our system. Thus, apart from solving Eqs. (1.11), which determines the evolution of an orbit, we have to simultaneously solve Eqs. (1.12) for each one of the k deviation vectors. Then, according to [C. [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method[END_REF], GALI k is dened as the volume of the k-parallelogram having as edges the k unit deviation vectors ŵi (t) = w i (t)/ w i (t) , i = 1, 2, ..., k. It can be shown, that this volume is equal to the norm of the wedge product (denoted by ∧) of these vectors: GALI k (t) = ŵ1 (t) ∧ ŵ2 (t) ∧ . . . ∧ ŵp (t) . (1.15) We note that in the above equation the k deviation vectors are normalized but their directions are kept intact.

In principal and for TI systems, the GALI k (t) for regular orbits remains practically constant and positive if k is smaller or equal to the dimensionality of the torus on which the motion occurs, otherwise, it decreases to zero following a power law decay. For the chaotic ones, all GALI k (t) tend exponentially to zero with exponents that depend on the rst k LEs of the orbit [C. [START_REF] Skokos | Geometrical properties of local dynamics in Hamiltonian systems: The Generalized Alignment Index (GALI) method[END_REF]C. Skokos, T. Bountis, and C. Antonopoulos 2008]. Nevertheless, in the TD case studied in [START_REF] Manos | Interplay between chaotic and regular motion in a time-dependent barred galaxy model[END_REF] and also here, the way the theoretical estimation of the GALI's exponential rates are strongly related to the LEs, being more complicated and still open to further inquiry.

The procedure used for the detection of the several dierent dynamical epochs of the TD system is the following: We evolve the GALI k with k = 2 or k = 3 (i.e., using 2 or 3 deviation vectors) for the 2d.o.f or 3-d.o.f. cases respectively and whenever GALI k reaches very small values (i.e. GALI k ≤ 10 -8 ) we re-initialize its computation by taking again k new random orthonormal deviation vectors, which resets the GALI k = 1. We allow then these vectors to evolve under the current dynamics. For time intervals where the index decays exponentially corresponds to chaotic epochs while in the other cases to non-chaotic. The reason in doing this is that we need to follow the current dynamical stability of an orbit under-study which in principle can interplay between chaotic and regular for dierent epochs during the total time evolution. Thus, let us assume that a trajectory experiences chaotic dynamics and later on drifts to a regular regime. The volume formed by the deviation vectors will rst shrink exponentially to very small values and remain small throughout the whole evolution unless one reinitializes the deviation vectors and their volume, in order to allow them to `feel' the new current dynamics. However, when we are interested in more general dynamical trends in time (less details), we will x time intervals and we will re-initialize the deviation vectors in the beginning of each one.

Results

The 2 D.O.F. Time-Independent Case Shedding some light on the underlying TI dynamics is an important step for understanding the more complicated case of the fully 3-d.o.f. TD model where all parameters vary simultaneously in time.

By setting z, p z equal to zero at t = 0 (remaining zero at all times) in the Hamiltonian Eq. (1.10), the orbits' motion is then restricted in the 2 dimensional (x, y) space. Note that here t = 0 refers to the t 0 = 1.4 Gyr of the N -body simulation. We can then study, in a frozen potential, individual orbits and the stability of the phase space, in terms of detecting and locating chaotic and regular motion, for the several sets of the potential parameters at deferent times, as derived from the N -body simulation. We shall begin by choosing xed parameters from four time snapshots, i.e., at t = 1.4 Gyr, t = 4.2 Gyr, t = 7.0 Gyr and t = 11.2 Gyr and we integrate orbits for 10 Gyr.

In Fig. 1.5, we present the Poincaré Surface of Section (PSS) dened by x = 0, p x ≥ 0 with H = -0.19, for three typical orbits being integrated for 10 Gyr. The set of parameters for the bar, disc and halo components are chosen from the ts with the 3-d.o.f. TD Hamiltonian at t = 7.0 Gyr of the N -body simulation (see see Table 1 in [Manos and Machado 2014] for more details). The blue ( * ) points on the PSS correspond to a disc regular orbit, forming a curve by the successive intersections with the plane x = 0, having initial condition (y, p y ) = (-1.5, 0.0) with x = 0 and p x = H(x, y, p y ) (called `DR' from now on). Its projection on the (x, y)-plane is shown in the top left inset panels of Fig. 1.5 and coloured in blue. The GALI 2 for this orbit conrms that its motion is regular by oscillating to a positive value during its evolution in time as well as the MLE σ 1 following a power law decay [see second and third top inset panels of Fig. 1.5 respectively (blue), note that the axis here are in lin-log scale]. The three small black curves in the central part of Fig. 1.5, marked with (×), are formed by the successive intersections of an initial condition with (y, p y ) = (-0.4, 0.0) (we will call it `BR'). From its projection on the (x, y)-plane [rst bottom inset panel of Fig. 1.5 (black)], it is evident that it is a bar-like orbit elongated along the long x-axis. It surrounds a stable periodic orbit of period 3 in the center of these islands and its regular dynamics is clearly revealed from the evolution of the GALI 2 and the MLE σ 1 evolution [second and third top (black) inset panels of Fig. 1.5 respectively, note that the axis are again in lin-log scale]. The central scattered red points on the PSS, marked with (+), correspond to a chaotic orbit with initial condition (y, p y ) = (2.5, 0.0) (called `DC' from now on). In the three inset panels positioned vertically in the right part of the Fig. 1.5 (red), we depict its projection on the (x, y)-plane (top panel). Its GALI 2 successive and exponential decrease to zero in time (middle panel) indicates its chaotic nature. Notice that we re-initialize the deviation vectors each time the index becomes small (≤ 10 -8 ). Its MLE σ 1 , as expected, converges to a positive value (bottom inset).

Exploiting this information, we rst choose a sample of 2-d.o.f. Hamiltonian function values (for the four times mentioned above and the respective sets of parameters), from the interval of energies where the majority of the N -body simulation's particles is more probable to be found. From this In the insets we depict their projection on the (x, y)-plane together with the GALI 2 and MLE σ 1 evolution in time (see Table 1 in [Manos and Machado 2014] for the exact parameters).

sample, we select a subset of representative energies to illustrate typical phase space structures, focusing at this point on the underlying dynamics. Then, we chart the regular and chaotic regimes of the phase space with GALI 2 . In Fig. 1.6, we have used a grid of 100 000 initial conditions on the (y, p y )-plane of the corresponding PSS and we have constructed a chart of the chaotic and regular regions similar to the PSS, but with more accuracy and higher resolution, in a similar manner just like in [START_REF] Manos | Regular and chaotic orbits in barred galaxies -I. Applying the SALI/GALI method to explore their distribution in several models[END_REF], using the GALI 2 method. The dierent colour corresponds to the dierent nal value of the GALI 2 after 10 Gyr (10 000 time units) for orbits representing each cell of the grid. The yellow colour corresponds to regular orbits (and areas) where the GALI Using the above approach, we can measure and quantify the variation of the percentage of regular orbits in the phase space as the total energy increases for a specic choice of potential parameters at same xed times. The chosen values of the Hamiltonian functions cover the range of the available energy interval up to the value of the escape energy which in general is dierent.

Although the main general trend is that this percentage decreases as the energy grows, its behavior changes at high energy values where it is no more monotonic. Note that this happens for energy values H > -0.19 out of the range of N -body simulation orbits. In Fig. 1.7 we show the variation of percentages of regular motion as a function of the energy H for the dierent sets of parameters at t = 1.4 Gyr, t = 4.2 Gyr, t = 7.0 Gyr and t = 11.2 Gyr. The threshold GALI 2 ≥ 10 -8 was used to characterize an orbit as regular and GALI 2 < 10 -8 as chaotic which will also be the chaos criterion/threshold throughout this paper. We should emphasize that these percentages refer to a set of initial conditions that cover uniformly the whole 2-d.o.f. phase space. On the other hand, an ensemble of trajectories extracted from the N -body simulation does not necessarily populate `democratically' the phase space. By simply inspecting these percentages one can claim that the fraction of regular motion is systematically larger for later times, and for all energies. When looking and comparing the phase space for early times, i.e., t = 1.4, 4.2 Gyr (rst and second row) and late times, i.e., t = 7.0, 11.2 Gyr (third and fourth row) in Fig. 1.6, we may see that the central island of stability, originating bar-like orbits, is becoming larger as the time grows and this is even more evident for the relatively larger energies (see the third and fourth row from top to bottom). This indicates that the bar component becomes gradually more important and dominant. (x, y, z, p x , p y , p z ) ≈ (-4.543100, 0.499639, -0.162627, 0.048798, -0.218718, 0.002898) (we will refer to this orbit from now on as B1), which is iterated for 10 Gyr (10 000 time units) for 4, 4.2, 7.0, 11.2 Gyr are given in [Manos and Machado 2014] (Table 1).

the TD potentials mentioned above. The starting and complete set of parameters for the TD model is taken at t 0 = 1.4 Gyr by the ts with the N -body simulation. Note that in all gures' panels, we set everywhere the t 0 equal to zero instead of t 0 = 1.4 Gyr.

In the rst row of the rst block in Fig. 1.8(a,b,c,d), we show its projection on the (x, y)-plane for four successive time intervals of ∆t = 2.5 Gyr, i.e., all the three potential components V B , V D , V H are time-dependent and the total energy is not in general conserved. The colour bar next to each panel corresponds to the time (in Gyr), hence the most recent epochs of the orbit are coloured with yellow while those in the earlier ones with dark blue or black. In Fig. 1.8(e) we show its GALI 3 , capturing accurately the chaotic nature of the orbit during the rst [Fig. 1.8(a)], third [Fig. 1.8(c)] and fourth [Fig. 1.8(d)] time windows by decaying exponentially to zero. On the other hand, in the second [Fig. 1.8(b)] time window its regular (even by just looking its projected morphologically on the (x, y)-plane) behaviour is successfully revealed by the uctuates to a non-zero value of the -8 ). It turns out that the B1 begins as a regular disc-like orbit during the rst 2.5 Gyr and, as the bar starts forming and growing, it gradually evolves to a chaotic bar-like orbit until the end of the integration. We may notice how hard it is for the nite time MLE σ 1 [Fig. 1.8(f )]

to capture these dierent dynamical dierent transitions and epochs due to its time-averaged denition (see also [START_REF] Manos | Interplay between chaotic and regular motion in a time-dependent barred galaxy model[END_REF]). Furthermore, its power law decay for regular time intervals and its tendency to positive values for chaotic ones are of the same order of magnitude making it rather hard to use the temporary value σ 1 as a safe criterion of regular and weak or strong chaotic motion.

In Fig. 1.9, and in a similar manner as in Fig. 1.8, we show another characteristic disc-like orbit for most of the total of integration with initial condition:

(x, y, z, p x , p y , p z ) ≈ (-5.14416, -1.345540, 0.277956, 0.140120, -0.219648, 0.000338) (we will refer to this orbit from now on as `D1'). The evolutionary scenarios are again the same as before, i.e., in the rst row of the rst block in Fig. 1.9(a,b,c,d), we present its projection on the (x, y)plane for dierent time windows. The D1 orbit experiences a regular epoch during its rst 2.5 Gyr, then gradually becomes chaotic switching to a bar-like shape and nally becomes a chaotic but disclike now orbit. Its regular and chaotic epochs are accurately captured by the GALI 3 [Fig. 1.9(e)],

uctuating to constant value for the rst 2.5 Gyr and then successively decaying exponentially to zero for the rest of the integration. The MLE σ 1 [Fig. 1.9(f )] also reveals this dynamical evolution, by decaying with a power law for the regular part and converging to non-zero value for the three last time windows. However, here the motion does not present any further transition and/or interplay between regular and chaotic motion and again (as for the B1 orbit) the order of magnitude for the σ 1 is not varying suciently enough to lead to a safe conclusion at certain times without seeing its whole time evolution.

In [Manos and Machado 2014], it was found that the overall fraction of chaotic motion in the disc decreases as the bar grows stronger. To explore the spatially resolved evolution of chaos throughout the stellar disc, we resort to the analysis of an ensemble of orbits. From the N -body simulation, we select a sample of 1 × 10 5 disc particles at the time t 0 = 1.4 Gyr where the bar has already started to be formed and starts growing from that point on. Then, their coordinates are used as an ensemble of initial conditions to be evolved in the presence of the time-dependent analytical potential. We evolve these orbits for 10 Gyr and study their dynamical behaviour. In order to avoid confusion, from now on we reset the t 0 to be zero (starting point of our simulations). We divide the total integration time in four intervals of ∆t = 2.5

Gyr, re-initializing the GALI 3 index at the beginning of each window [Manos and Machado 2014].

The orbit is considered regular (non-chaotic) if its GALI 3 remains greater than 10 -8 during a given time window; and it is considered chaotic if it reaches GALI 3 ≤ 10 -8 . In this manner, we are able to compute fractions of chaotic motion within each time window. Additionally, at a given instant in time, we can also compute spatially resolved chaos fractions in dierent regions of the disc.

A global picture of the spatial distribution of regular and chaotic motion in the disc can be seen in the two bottom rows of Fig. 1.11, which displays the face-on and edge-on views of the ensemble of disc particles at the end of each time window, coloured by the GALI 3 index (being chaotic towards the In Fig. 1.11 and in all other such projections, the particles and orbits are displayed in the reference frame that rotates with the bar. Thus, the bar major axis always lies along the direction of the x-axis. blue, and regular towards the yellow). Some major results are already noticeable even by eye. First, the striking decrease of chaos within the bar region can be clearly seen. Secondly, even though the gap is a very low density region, it seems to hold a good portion of the chaotic orbits. Third, the outer disc as well as the ring, to a degree seem quite dominated by regular motion. Finally, another outstanding feature is the peanut-shaped view of the bar seen in the edge-on projection (sometimes called X-shaped bulge). Remarkably, particles that depart considerably from the z = 0 plane are mostly chaotic.

In order to quantify in more detail these results, we measure the fraction of chaotic orbits as a function of time in each region (i.e. at the end of each time window, we obtain the number of particles having GALI 3 ≤ 10 -8 in a region divided by the total number of particles within that region). The result is shown in Fig. 1.12. The fraction of chaotic motion within the bar drops from nearly 40 per cent to less than 10 per cent. The outer disc remains essentially regular, with a non-zero but negligible appearance of chaos throughout the evolution. The fact that the ring region undergoes an initial increase in chaos can be ascribed in good measure to the edges of the bar. The gap region displays some interesting behaviour. Between the rst and second time windows, the gap becomes depleted in terms of total number of particles, but at the same time its fraction of chaotic motion increases. From then on, it decreases, but the gap continues to be the region holding the highest local fraction of chaotic motion in the disc. The large amount of chaotic motion seen in the gap region is not unexpected. In fact it is well known that orbits that oscillate between the Lagrangian points L 1 and L 2 are unstable and therefore the transition zone between the bar and the disk is expected to be chaotic, see e.g. [START_REF] Harsoula | Orbital structure in N-body models of barred-spiral galaxies[END_REF]Athanassoula, Romero-Gómez, and Masdemont 2009]. Instead of considering the state (chaotic or regular) of each particle at specic instants the ends of the four time windows, we here account for the changes of their state. For example, one given orbit that was found to be regular at the end of the evolution might have been chaotic at the beginning, or it might have been continuously regular. In either case, where did such particles originate? Do the particles that change dynamical state (and those that don't) share a common locus at the beginning of the evolution? To explore these issues, we will examine separately the orbits that change dynamical behaviour and those that do not. This will allow us, in a sense, to map the origins and the destiny of face-on and edge-on vies of the orbits that are permanently regular (cyan), and those that are permanently chaotic (red). Bottom rows: face-on and edge-on views of the orbits that start regular and end chaotic (green), and those that start chaotic and end regular (purple).

Orbits that are permanently regular. As regards morphology, the permanently regular ones are qualitatively unremarkable, in the sense that they occupy almost any region of the galactic disc. There is thus little qualitative distinction between them and the entire ensemble and they merely map the normal evolution of the galactic disc as a whole. The only noticeable structures that are not quite covered by these orbits are the gap region, and, vertically, the peanut. (Due to the method employed to create the initial conditions in the N -body simulation of [START_REF] Machado | Loss of halo triaxiality due to bar formation[END_REF], there is a residual transient seen as a vague spiral pattern at t = 0 and it subsides on a short time scale.)

Orbits that are permanently chaotic. The permanently chaotic orbits, on the the other hand, display peculiar features. They are tightly restricted to the region of the bar, and partially to the gap. Indeed, they spend nearly the entire evolution conned within this region. There is not one single permanently chaotic orbit to be found in the outer regions of the disc. At the instant t = 0, these particles whose future destiny is to be permanently chaotic are initially located within a reasonably well-dened ring, i.e. they are mostly found within 2 kpc < r < 4 kpc.

For the remainder of the orbits, we will focus on two regimes: those that start regular and end chaotic (6.1 per cent), and those that start chaotic and end regular (20.8 per cent), regardless of the intermediate states (i.e. the transitions in the second and third time windows). Finally, there is a subset of orbits (7.9 percent) that do undergo two changes of regime, but nevertheless nish as they started; these will be disregarded. Let us consider now the two cases where the nal state diers from the initial state (lower rows of Fig. 1.13):

Orbits that begin regular and end chaotic. These start at t = 0 from a similar locus as the permanently chaotic, but here the ring is slightly larger and more diuse. This subset also includes some orbits very close to the origin (r < 0.5 kpc) in the beginning, which are not present in the permanently chaotic case. The volume occupied by these orbits contracts gradually, but they are more extended than the permanently chaotic ones, encompassing the region of the gap at later times as well. They are also vertically extended, being the major contributors to the structure of the peanut.

In fact, this is the only subset of particles which signicantly populates the peanut, in the regions of about 2.5 kpc < |z| < 5 kpc of height. The gap region, and mainly the peanut, are regions where chaos is important. However, it is only the initially regular and nally chaotic orbits that depart considerably from the plane. The orbits that were already chaotic from the beginning do not visit such heights.

Orbits that begin chaotic and end regular.

Interestingly, the initial locus of this subset is approximately the complement of the previous case. Here, the orbits at t = 0 occupy the region internal to the ring dened by the previous case, while avoiding the very centre. In the third row of Fig. 1.13, the purple points overlap green points in the t = 2.5 -10 Gyr frames. But in the t = 0 frame, the purple points ll precisely an empty region. Subsequently, the initially chaotic orbits evolve to be essentially part of the bar and end regular.

Origin and morphology of boxy orbits

Regular and chaotic two-dimensional and three-dimensional orbits of stars in such models of a galactic potentials provide with boxy orbital structures, which are the largest part of the bar galaxy component.

Boxiness in the equatorial plane is associated either with quasi-periodic orbits in the outer parts of stability islands, or with sticky orbits around them, which can be found in a large range of energies.

In the 3D barred models in [Patsis and Katsanikas 2014a;[START_REF] Patsis | The phase space of boxy-peanut and X-shaped bulges in galaxies -I. Properties of non-periodic orbits[END_REF]] it has been suggested that inner boxy features can be built by means of quasi-periodic orbits at the edges of the stability islands of the x1 family, as well as with sticky orbits just beyond the last invariant torus around the stable x1 periodic orbit (see e.g. [START_REF] Manos | Probing the local dynamics of periodic orbits by the generalized alignment index (GALI) method[END_REF] for applications of the GALI index for periodic orbits of dierent stability and their vicinity). It has been also proposed that such orbits support boxiness both in face-on, as well as in edge-on projections at the central region of the bar (about within half the way to the end of the bar). A similar dynamical phenomenon was leading to the boxy features on the galactic plane in the bars of 2D barred-spiral models in [START_REF] Tsigaridi | Morphologies introduced by bistability in barred-spiral galactic potentials[END_REF]. By considering snapshots of the N -body simulation and for the orbital study we treat each snapshot independently, as an autonomous (TI) Hamiltonian system and calculate regular, sticky and chaotic orbits and hence examine the degree of chaoticity of the bar-supporting orbits.

In [Chaves-Velasquez et al. 2017], the location of boxy orbits in diagrams was indicated for certain TI analytical models for snapshots of the TD one. They are always found about the transition region from order to chaos. By perturbing such orbits in the vertical direction a class of 3D non-periodic orbits, which have boxy projections both in their face-on and side-on views, was found. Furthermore, the kind of orbits supporting double boxy morphologies in the successive snapshots was investigated, and the manner they evolve in time, i.e. from the model of the earlier snapshot to the model for the nal one.

In a rotating Ferrers bar, the elliptical periodic orbits of the main families are characterized by a single non-zero initial condition along the minor axis of the bar, namely their position along the y-axis in our models. The curve of zero velocity (ZVC) in a (E J , y 0 ) diagram separates the region where orbital motion is allowed from the region where it is not. Since the main family consists of direct periodic orbits, only the y 0 > 0 part of such a diagram is of interest for us. An (E J , y 0 ) diagram is the projection of a complete (E J , y 0 , p y 0 ) gure with all possible initial conditions. However, it is sucient for describing the properties of the orbits we present below. The line that gives the y 0 initial condition of the main family of periodic orbits is the characteristic curve of the model. Since we want to study chaoticity in a large range of energies, we have created such (E J , y 0 ) diagrams for the potentials of the three snapshots we study. In order to calculate the degree of chaoticity of the planar orbits around the main family of periodic orbit as we move from the center of the system towards corotation, we use the GALI 2 index. We have used the GALI 2 index to color-code each point in the allowed region in the (E J , y 0 ) areas. The shade of the color indicates the GALI 2 index that a given orbit, i.e. a point in the (E J , y 0 ) diagram, has at the end of the integration. In other words, the color of an (E J , y 0 ) point indicates if the orbit with y 0 initial condition at E J will lead to regular (large log 10 (GALI 2 ) values) or chaotic (very small log 10 (GALI 2 ) values) motion. At the borders between these regions we nd points with intermediate log 10 (GALI 2 ) values, which correspond to sticky chaotic orbits.

In Fig. 1.14, we show the chaoticity of the planar orbits on the equatorial plane of a TI model based on the 2nd snapshot (at t=4.2 Gyr from Fig. 1.1). We sample the GALI 2 index at two time windows. First after time t 1 , corresponding to 1 Gyr and then after time t 2 , corresponding to 10 Gyr.

In this way we investigate both the relatively short-term as well as the long-term behaviour of the orbits. Darker shades indicate more chaotic orbits. The color for each orbit is determined according to its log 10 (GALI 2 ) value and is taken from the colour bars given to the right of the gures.

Summary and outlook

The adequacy of the TD model and its verication is shown by (i) the similarity of the rotation curve with the one from the N -body simulation, ensuring the global dynamics are well approximated, (ii) the morphological similarity of the disc and bar features with those of the N -body simulation. and, (iii) the very good quantitative agreement of the length and strength of the bar in the resulting mock snapshots with the N -body bar. Such comparisons indicate that the TD model is able to adequately capture both the dynamics and the morphology of the barred galaxy model in question. Regarding the TD analytical model, we established stability trends in terms of estimating the amount of regular and chaotic motion in dierent time-windows. In this case, we used a more realistic set of initial conditions coming directly from the simulation itself and iterated them under the constructed TD potential. Even though, we did not manage to span the whole orbital richness of our ensemble of initial conditions, we have rather achieved to give a avor of the possible evolutions for individual trajectories. It turns out that the complete set of orbits tends to become relatively more `regular' in time.

By directly following the dynamical evolution of ensembles of orbits within the analytical potential, we were able to calculate the fractions of chaotic and regular motion resolved in both time and space.

With this information we could evaluate not only the global trends in time, but also across several regions of the galactic disc and of the halo, associating them with distinct morphological features.

We scrutinized the dierent changes of regime during the evolution, tracing the types of orbits back to their common origins.

The time-dependence of the analytical model ensures rather realistic dynamical transitions similar to an N -body simulation, i.e. bar formation and growth, development of a ring, a dynamical halo, etc. At the same time, this setup serves ideally to apply the GALI chaos detection method and in this way to determine the current chaotic (or otherwise) dynamical state of any given orbit at a xed time interval, something that it is extremely hard to do in an N -body simulation. Moreover, the TD model has the advantage over a derived frozen potential from an N -body simulation, from the point of view that it incorporates smooth dynamical evolution via its time-dependent parameters and allows us to follow the both dynamical and morphological transitions.

A recipe for building two-and three-dimensional boxy structures in rotating bars was presented in [Chaves-Velasquez et al. 2017]. The basic idea is based on rst starting with the planar backbone of periodic orbits for building a bar, namely with the well known x1 family. Then, instead of populating the model with regular quasi-periodic orbits encountered in the immediate neighborhood of the periodic orbit, one considers either periodic orbits close to the last Kolmogorov-Arnold-Moser torus (see e.g. [START_REF] Contopoulos | Order and Chaos in Dynamical Astronomy[END_REF]) or, more eciently, the sticky chaotic orbits that surround the islands of stability, as they appear in the surfaces of section. In this way, such a selection of orbits secures a boxy morphology on the plane.

In a 3D model, when we eject out of the plane particles that follow the 2D boxy orbits by adding a p z = 0 perturbation, we nd that there is always a ∆p z range of perturbations for which all three projections of the 3D orbits are boxy. A remarkable property of these sticky boxy orbits is the formation of an X feature embedded in the bar in the face-on projections. For as long as we have the usual ellipses of the x1 family (or the x1-tree in 3D models according to [Ch. [START_REF] Skokos | Orbital dynamics of three-dimensional bars I. The backbone of three-dimensional bars. A ducial case[END_REF]) in a rotating bar, we can nd a class of boxy 2D and 3D orbits. They are sticky chaotic orbits as their GALI 2 index indicates and they can support the bar, or a part of the bar, for many Gyr.

Observational features that can be reproduced by using such orbits as building blocks, can rstly be the boxy-or peanut-shaped bulges in the central parts of the bars. In these cases in the face-on views of the galaxies, we will observe boxy isophotes in their central parts, inside the bar, as in the sample of galaxies presented by [START_REF] Erwin | Peanuts at an angle: detecting and measuring the three-dimensional structure of bars in moderately inclined galaxies[END_REF]. On the other hand, the work in This work focused on one particular galaxy simulation, but the method may be extended to dierent galaxy types, taking as input the results of other N -body simulations. Here we focused on a strongly barred galaxy to maximise the eects we wished to explore. Clearly a natural extension would be to compare the present results with alternative galaxy models of varying bar strengths, disc masses, halo proles, etc. Such a systematic exploration would reveal to which parameters the dynamical stability is most sensitive. For example, given the high fraction of chaotic motion found in the inner halo, the question arises as to what is the role of the dark matter prole in determining that behaviour. One might explore whether a more cuspy inner halo would help or rather hinder the rise of regular motion. A further development would be the inclusion of models containing gas [START_REF] Patsis | Stellar and gas dynamics of late-type barred-spiral galaxies: NGC 3359, a test case[END_REF]] and star formation. To this end, the hydrodynamical simulations of [Athanassoula, Machado, and Rodionov 2013b] would be ideally suited, since they already oer a systematic grid of models for galaxies with dierent halo triaxialities and dierent initial gas fractions, thus resulting in a variety of bar strength evolutions. More broadly, models derived from a fully cosmological hydrodynamical simulation of galaxy formation would oer an even more realistic scenario (e.g. [START_REF] Valluri | Halo Orbits in Cosmological Disk Galaxies: Tracers of Formation History[END_REF]) than the usual models of isolated galaxies. Finally, a specic issue that merits further analysis is the behaviour of the X-shaped (or boxy/peanut) bulge [Patsis and Katsanikas 2014b;[START_REF] Patsis | The phase space of boxy-peanut and X-shaped bulges in galaxies -I. Properties of non-periodic orbits[END_REF], particularly in light of the recent interest in the kinematics and structure of the Milky Way's own bulge [START_REF] Saito | Mapping the X-shaped Milky Way Bulge[END_REF][START_REF] Zoccali | The GIRAFFE Inner Bulge Survey (GIBS). I. Survey description and a kinematical map of the Milky Way bulge[END_REF].

In

+D=FJAH

Complex phenomena in quantum systems [START_REF] Casati | Stochastic behavior of a quantum pendulum under a periodic perturbation[END_REF][START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF][START_REF] Haake | Quantum Signatures of Chaos[END_REF]]. Through the Fourier transform connection between the time and energy, the dynamical localization reects itself also in the time-independent eigenfunctions, both in the eigenstates of the Floquet operator in time-periodic systems and in the eigenfunctions of the time-independent, classically chaotic, systems. Namely, if all classical transport times like the diusion time (time necessary to occupy the entire classically available chaotic part of the phase space) are all shorter than the Heisenberg time t H , we nd extended eigenstates. In the contrary case, one observes localized eigenstates.

Another important aspect of quantum chaos is the statistics of the energy spectra of classically chaotic quantal systems. One of the main cornerstones in the development of quantum chaos [START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF][START_REF] Haake | Quantum Signatures of Chaos[END_REF] is the nding that in classically fully chaotic, ergodic, autonomous Hamiltonian systems with purely discrete spectrum the uctuations of the energy spectrum around its mean behavior obey the statistical laws described by the Gaussian Random Matrix Theory (RMT) [START_REF] Mehta | Random Matrices[END_REF][START_REF] Guhr | Random-matrix theories in quantum physics: common concepts[END_REF], provided that we are in the suciently deep semiclassical limit. The latter condition means that all relevant classical transport times are smaller than Heisenberg time t H . This statement is known as the Bohigas -Giannoni -Schmit (BGS) conjecture [START_REF] Bohigas | Characterization of Chaotic Quantum Spectra and Universality of Level Fluctuation Laws[END_REF][START_REF] Casati | On the connection between quantization of nonintegrable systems and statistical theory of spectra[END_REF]. Since ∆E ∝ f , where f is the number of degrees of freedom (the dimension of the conguration space), we see that for suciently small -value the stated condition will always be satised. Alternatively, xing the , we can go to high energies such that the classical transport times become smaller than t H . The role of the antiunitary symmetries that classify the statistics in terms of dierent ensembles of RMT, namely Gaussian Orthogonal Ensemble (GOE), Gaussian Unitary Ensemble (GUE) or Gaussian Symplectic Ensemble (GSE), has been elucidated in [START_REF] Robnik | False time-reversal violation and energy level statistics: the role of anti-unitary symmetry[END_REF][START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF][START_REF] Haake | Quantum Signatures of Chaos[END_REF][START_REF] Mehta | Random Matrices[END_REF]]. The theoretical foundation for the BGS conjecture has been initiated rst by Berry [START_REF] Berry | Regular and irregular semiclassical wavefunctions[END_REF]], and later further developed by Richter and Sieber [START_REF] Sieber | Correlations between Periodic Orbits and their R?le in Spectral Statistics[END_REF], arriving nally in the almost-nal proof proposed by the group of F. Haake [START_REF] Müller | Semiclassical Foundation of Universality in Quantum Chaos[END_REF]Heusler et al. 2004;[START_REF] Müller | Semiclassical Foundation of Universality in Quantum Chaos[END_REF][START_REF] Müller | Periodic-orbit theory of universal level correlations in quantum chaos[END_REF]].

On the other hand, if the system is classically integrable, then Poisson statistics applies (see e.g. [START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF][START_REF] Haake | Quantum Signatures of Chaos[END_REF]] and references therein, and for recent advances [START_REF] Robnik | On spectral statistics of classically integrable systems[END_REF]). In the mixed type regime, where classical regular regions coexist in the classical phase space with chaotic regions, being a typical KAM-scenario which is the generic situation, the so-called Principle of Uniform Semiclassical Condensation (of the Wigner functions of the eigenstates; PUSC) applies, see e.g. [START_REF] Berry | Regular and irregular semiclassical wavefunctions[END_REF]]. Consequently the Berry-Robnik statistics [START_REF] Berry | Semiclassical level spacings when regular and chaotic orbits coexist[END_REF]Prosen and Robnik 1999a] is observed, again under the same semiclassical condition stated above, requiring that t H is larger than all classical transport times. For more details about the mixed type regime, see [START_REF] Berry | Semiclassical level spacings when regular and chaotic orbits coexist[END_REF][START_REF] Prosen | Energy level statistics in the transition region between integrability and chaos[END_REF]Prosen and Robnik 1994a;Prosen and Robnik 1994b;[START_REF] Prosen | Semiclassical energy level statistics in the transition region between integrability and chaos: transition from Brody-like to Berry-Robnik behaviour[END_REF]Prosen and Robnik 1999b;[START_REF] Grossmann | On level spacing distributions for 2D non-normal Gaussian random matrices[END_REF] , while some more recent developments in the eld were published in [START_REF] Batisti¢ | Dynamical localization of chaotic eigenstates in the mixed-type systems: spectral statistics in a billiard system after separation of regular and chaotic eigenstates[END_REF][START_REF] Batisti¢ | The intermediate level statistics in dynamically localized chaotic eigenstates[END_REF]Batisti¢ and Robnik 2013a;Batisti¢ and Robnik 2013b].

If the couplings between the regular eigenstates and chaotic eigenstates become important, due to the dynamical tunneling, we can use the ensembles of random matrices that capture these eects [START_REF] Vidmar | [END_REF]]. As the tunneling strengths typically decrease exponentially with the inverse eective Planck constant, they rapidly disappear with increasing energy, or by decreasing the value of the Planck constant. In such cases the regular and chaotic eigenstates can be separated and the dynamical localization in the chaotic eigenstates can be studied (see [START_REF] Prosen | Proceedings of the International School of Physics Enrico Fermi, Course CXLIII[END_REF]] for a review and references therein). In such a situation it turns out that the Wigner functions of the chaotic eigenstates no longer uniformly occupy the entire classically accessible chaotic region in the classical phase space, but are localized on a proper subset of it. Indeed, this has been analyzed with unprecedented precision and statistical signicance in [START_REF] Batisti¢ | Dynamical localization of chaotic eigenstates in the mixed-type systems: spectral statistics in a billiard system after separation of regular and chaotic eigenstates[END_REF]Batisti¢ and Robnik 2013a;Batisti¢ and Robnik 2013b] for the case of mixed type systems. The important discovery is that the level spacing distribution of the dynamically localized chaotic eigenstates in periodic, as well as time-independent systems, is exceedingly well described by the Brody distribution (see e.g. [START_REF] Brody | Randommatrix physics: spectrum and strength uctuations[END_REF]] and references therein) with the Brody parameter values β BR within the interval [0, 1],

where β BR = 0 yields the Poisson distribution in case of the strongest localization, and β BR = 1 gives the Wigner surmise (2D GOE, as an excellent approximation of the innite dimensional GOE). The Brody distribution was found to t the empirical data much better than the distribution function proposed in [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] and the references therein, characterized by the parameter β IZ . It is well known that the Brody distribution so far has no theoretical foundation, but there are many empirical results showing that we have to consider it seriously thereby being motivated for seeking its physical foundation.

In this chapter, we review the work published mainly in [Manos and Robnik 2013;[START_REF] Batisti¢ | The intermediate level statistics in dynamically localized chaotic eigenstates[END_REF]Manos and Robnik 2014;Manos and Robnik 2015], where the quantum kicked rotator (QKR) [START_REF] Casati | Stochastic behavior of a quantum pendulum under a periodic perturbation[END_REF]] was explored from the classical point of view (the standard map, SM), its quantum analog was analyzed using the N -dimensional model of Izrailev, and the semiclassical connection between the two pictures was considered. We shall treat the cases with the classical dimensionless kick parameter K in the range K ∈ [START_REF] Manos | Survey on the role of accelerator modes for anomalous diusion: the case of the standard map[END_REF]35], and for some purposes even up to K = 70, and in the end shall focus on the case K = 10, which is the most chaotic one in the sense that it is fully chaotic with minimal (in fact practically undetected) regular regions among all cases K ∈ [START_REF] Manos | Survey on the role of accelerator modes for anomalous diusion: the case of the standard map[END_REF]70] and among them best exhibits normal diusion. Izrailev's N -dimensional model introduced and discussed in [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] is treated for various N ≤ 3000 values, which in the limit N → ∞ tends to the QKR. Due to the niteness of N the observed (dimensionless) localization length of the eigenfunctions in the space of the angular momentum quantum number does not possess a sharply dened value, but has a certain distribution instead. Its reciprocal value is almost Gaussian distributed. This might be expected on the analogy with the nite time Lyapunov exponents in the Hamiltonian dynamical systems. In order to corroborate the theoretical ndings on this topics we perform the numerical analysis of the nite time Lyapunov exponents in the SM (classical kicked rotator), especially the decay of the variance.

Indeed, in the Shepelyansky picture [START_REF] Shepelyansky | Localization of quasienergy eigenfunctions in action space[END_REF]] the localization length can be obtained as the inverse of the smallest positive Lyapunov exponent of a nite 2k-dimensional Hamiltonian system associated with the band matrix representation of the QKR, where k is the quantum kick parameter (to be precisely dened below) and where N plays the role of time.

However, unlike the chaotic classical maps or products of transfer matrices in the Anderson tightbinding approximation, where the mean value of the nite time Lyapunov exponents is usually equal to their asymptotic value at innite time and the variance decreases inversely with time, as we also carefully checked, here the distribution is found to be independent of N : It has a nonzero variance even in the limit N → ∞. This happens because the quantum kicked rotator at N = ∞ cannot be exactly modeled with nite bandwidth (equal to 2k) band matrices, but only approximately. Hence, the underlying Hamiltonian system of the Shepelyansky picture has a growing dimension with N , implying asymptotically an innite set of Lyapunov exponents and behavior dierent from the nite dimensional Hamiltonian systems. The observation of the distribution of the localization length around its mean value with nite variance also explains the strong uctuations in the scaling laws of the kicked rotator, such as e.g. the entropy localization measure as a function of the theoretical scaling parameter Λ, to be discussed below. On the other hand, the two dierent empirical localization measures, namely the mean localization length as extracted directly from the exponentially localized eigenfunctions and the measure based on the information entropy of the eigenstates, are perfectly well linearly connected and thus equivalent. Therefore these results call for a rened theory of the localization length in the quantum kicked rotator and similar systems, where we must predict not only the mean value of the inverse localization length but also its (Gaussian) distribution, in particular the variance [Manos and Robnik 2015].

The classical and quantum kicked rotator models

The classical kicked rotator model: the standard map

The kicked rotator was introduced in [START_REF] Casati | Stochastic behavior of a quantum pendulum under a periodic perturbation[END_REF]] and is one of the key model systems in classical and quantum chaos, especially for time-periodic (Floquet) systems. The Hamiltonian function is:

H = p 2 2I + V 0 δ T (t) cos θ. (2.1)
Here p is the (angular) momentum, I the moment of inertia, V 0 is the strength of the periodic kicking, θ is the (canonically conjugate, rotation) angle, and δ T (t) is the periodic Dirac delta function with period T . Between the kicks the rotation is free, therefore the Hamilton equations of motion can be immediately integrated, and thus the dynamics can be reduced to the SM, or so-called Chirikov-Taylor map, given by:

   p n+1 = p n + V 0 sin θ n+1 , θ n+1 = θ n + T I p n , (2.2)
where the quantities (θ n , p n ) refer to their values just immediately after the n-th kick. Then, by introducing new dimensionless momentum P n = p n T /I, we get:

{ P n+1 = P n + K sin θ n+1 , θ n+1 = θ n + P n , (2.3)
where the system is now governed by a single classical dimensionless kick parameter K = V 0 T /I, and the mapping is area preserving.

The generalized diusion process of the SM [Eq. ( 2.

3)] is dened by:

(∆P ) 2 = D µ (K)n µ , (2.4)
where n is the number of iterations (kicks), and the exponent µ is in the interval [0, 2], and all variables P , θ and K are dimensionless. Here D µ (K) is the generalized classical diusion constant. In the case µ = 1 we have the normal diusion, and D 1 (K) is then the normal diusion constant, whilst in the case of anomalous diusion we observe subdiusion when 0 < µ < 1 or superdiusion if 1 < µ ≤ 2.

In the case µ = 2 we have the ballistic transport which is associated strictly with the presence of accelerator modes.

In the case of the normal diusion µ = 1 the theoretical value of D 1 (K) is given in the literature, e.g. in [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]]:

D 1 (K) = { 1 2 K 2 [1 -2J 2 (K) (1 -J 2 (K))] , if K ≥ 4.5 0.15(K -K cr ) 3 , if K cr < K ≤ 4.5 , (2.5)
where K cr 0.9716 and J 2 (K) is the Bessel function. Here we neglect higher terms of order K -2 .

However, there are many important subtle details in the classical diusion further discussed below.

The dependence of the diusion constant for the growth of the variance of the momentum on K is very sensitive, and described in the theoretical result [Eq. (2.5)], and fails around the period 1 accelerator mode intervals

(2πn) ≤ K ≤ √ (2πn) 2 + 16, (2.6)
n any positive integer. In these intervals for the accelerator modes n = 1 we have two stable xed points located at p = 0, θ = π -θ 0 and p = 0, θ = π + θ 0 , where θ 0 = arcsin(2π/K). There are two unstable xed points at p = 0, θ = θ 0 and p = 0, θ = 2π -θ 0 . For example, in the case K = 6.5 we have θ 0 ≈ 1.31179. Moreover, as the diusion might even be anomalous, we have recalculated the eective diusion constant D eff = (∆P ) 2 /n numerically, which in general is not equal to the D µ dened in Eq. (2.4). In Fig. In Fig. 2.2 we show the variance of the momentum P in the SM [Eq. ( 2.

3)] with K = 6.5 (red crosses) where small islands and accelerator mode of period 1 are present and K = 10.0 (blue stars)

where the phase space is fully chaotic for the same initial conditions as in Fig. 2.1 as a function of the discrete time n (number of iterations), in log-log representation. The two slopes associated with dierent types of diusion are µ(K = 6.5) = 1.61252 (dotted), µ(K = 10.0) = 0.991334 (solid) with standard deviation errors ±0.01271 (0.7881%) and ±0.0009537 (0.0962%) respectively.

In this manner we have calculated the diusion exponent µ for all K on the interval K ∈ [K cr , 70] and the result is shown in Fig. 2.3. We show the diusion exponent µ as a function of K after n = 5000 iterations, using a ne grid of 314 × 314 (≈ 100000) initial conditions on the plane (θ, P ) = (0, 2π).

The µ is calculated by the slopes, of the lines of the variance of the momentum P as a function of iterations, as it is described in [Eq. (2.4)] and for a grid of cells on the entire phase space. The intervals on the black horizontal line µ = 0.9 indicate the intervals of stable accelerator modes of period 1 [Eq. (2.6)]. All intervals of K with exponent µ ≈ 1 are associated with normal diusion processes. The large peaks (appearing mainly for K > 2π marked with full black circles) reect the anomalous diusion due to accelerator modes [of period 1, being located inside the intervals predicted by the Eq. (2.6)]. However, there is a number of relatively smaller peaks for K < 2π (more clearly presented in the inset panel of Fig. 2.3), whose origin is accelerator modes of higher period as we will see later, and also for 2π < K < 4π, both these sets are marked with empty circle. With the symbol (×) we mark few typical examples, close to those peaks, for which the diusion is normal and are studied in detail in [Manos and Robnik 2014].

All the large peaks for K > 2π, marked with full black circles in Fig. 2.3, correspond to regimes with accelerator modes of period 1 and they decrease monotonically as a power law:

f (x) = ax b , (2.7)
where a = 2.41645 and b = -0.195896 [blue dotted line in Fig. 2.3, with asymptotic standard error ±0.04294 (1.777%) and ±0.00537 (2.741%) respectively] indicating that for K > 70 their eect decreases signicantly. On the other hand, the size of the successive accelerator modes of period 1 intervals decays with a power law dened simply and analytically by the Eq. (2.6).

In order to understand the eect of the presence of accelerator modes in the diusion and transport properties of the phase space in the SM, we rst picked an, as much as possible, representative sample of K-values. In more detail, we included in our test-cases all the K-values which correspond to all the main peaks appearing in Fig. 2.3 with K > 2π together with a few cases from the `plateaus' of this curve. Furthermore, we took into account the peaks occurring for 1 K 2π (see the empty black circles in the inset zoom in Fig. 2.3) which are associated with accelerator modes of higher periodicity, as it will be seen thereupon. The case with K = 3.8, whose µ value is ≈ 1, is chosen for comparison reasons from the plateau and as it turns out has no accelerator modes in its phase space causing anomalous diusion. Here we should stress that we repeated the same procedure for larger number of iterations n = 10000 and it turns out that the exponent µ has well converged to the values shown in Fig. 2.3.

The distribution of the momenta in the case of normal diusion is found to be perfect Gaussian, In the case of normal diusion (Gaussian statistics) for the above quantities, one expects to nd α = 2 for the Lévy stable distribution and diusion exponent µ = 1, while in the general case we nd other values. We have calculated the µ exponent emerging from a small box/ensemble of initial conditions in the phase space and thus produced the µ-landscape in the phase space of the SM. Furthermore, we have also employed the GALI-method [Ch. Skokos, T. Bountis, and Ch.

Antonopoulos 2007; C. [START_REF] Skokos | The Smaller (SALI) and the Generalized (GALI) Alignment Indices: Ecient Methods of Chaos detection[END_REF] to calculate the GALI-index in order to identify the regular and chaotic regions in the phase space of the SM, which also quanties the degree of chaos (indirectly). The GALI-landscape and the µ-landscape are found to correspond very well to each other, where it seems that the µ-plot contains more information than GALI-plot [Manos and Robnik 2014]. Hence, the SM exhibits normal diusion for most of the K values on the interval [K cr , 70],

except for the accelerator mode intervals where anomalous diusion is observed, with the exponent µ typically being larger than 1. Using these plots and the described methodology we found that the case K = 10 is the closest to full chaos (no regular islands present) and exhibits the normal diusion for all initial conditions in the phase space of the SM.

The quantum kicked rotator model

The quantum kicked rotator (QKR) is the quantized version of Eq. ( 2.1), namely:

Ĥ = - 2 2I ∂ 2 ∂θ 2 + V 0 δ T (t) cos θ.
(2.8)

The Floquet operator F acting on the wavefunctions (probability amplitudes) ψ(θ), θ ∈ [0, 2π), upon each period (of length T ) can be written as (see [START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF]):

F = exp ( - iV 0 cos θ ) exp ( - i T 2I ∂ 2 ∂θ 2
) ,

(

where now we have two dimensionless quantum control parameters: .10) which satisfy the relationship K = kτ = V 0 T /I, K being the classical dimensionless control parameter of Eq. (2.3). By using the angular momentum eigenfunctions:

k = V 0 , τ = T I , ( 2
|n = a n (θ) = 1 √ 2π exp(i n θ), (2.11)
where n is any integer, we nd the matrix elements of F , namely:

F m n = m| F |n = exp ( - iτ 2 n 2 ) i n-m J n-m (k), (2.12)
where J ν (k) is the ν-th order Bessel function. For a wavefunction ψ(θ) we shall denote its angular momentum component (Fourier component) by:

u n = n|ψ = ∫ 2π 0 a * n (θ)ψ(θ) dθ = 1 √ 2π ∫ 2π 0 ψ(θ) exp(-inθ) dθ.
(2.13)

The QKR has very complex dynamics and spectral properties. As the phase space is innite (cylinder), p ∈ (-∞, +∞), θ ∈ [0, 2π), the spectrum of the eigenphases of F , denoted by ϕ n , or the associated quasienergies ω n = ϕ n /T , can be continuous, or discrete. It is quite well understood that for the resonant values of τ : τ = 4πr q , (2.14) q and r being positive integers without common factor, the spectrum is continuous, as rigorously proven by Izrailev and Shepelyansky [START_REF] Izrailev | Quantum resonance for a rotator in a nonlinear periodic eld[END_REF], and the dynamics is (asymptotically) ballistic, meaning that starting from an arbitrary initial state the mean value of the momentum p increases linearly in time, and the energy of the system E = p2 /(2I) grows quadratically without limits. For the special case q = r = 1 this can be shown elementary. Such behavior is a purely quantum eect, called the quantum resonance. Also, the regime of quadratic energy growth manifests itself only after very large time, which grows very fast with the value of the integer q from Eq. (2.14), such that for larger q this regime practically cannot be observed.

For generic values of τ /(4π), being irrational number, the spectrum is expected to be discrete but innite. But the picture is very complicated. Casati and Guarneri [START_REF] Casati | Non-recurrent behaviour in quantum dynamics[END_REF] have proven that for τ /(4π) suciently close to a rational number, there exists a continuous component in the quasienergy spectrum. So, the absence of dynamical localization for such cases is expected as well. Without a rigorous proof, we nally believe that for all other (good) irrational values of τ /(4π)

we indeed have discrete spectrum and quantum dynamical localization. In such case the quantum dynamics is almost periodic, and because of the eective niteness of the relevant set of components u n and of the basis functions involved, just due to the exponential localization (see below), it is even eectively quasiperiodic (eectively there is a nite number of frequencies), and any initial state returns after some recurrence time arbitrarily close to the initial state. Thus the energy cannot grow indenitely.

Localization and Diusion Properties

The asymptotic localized eigenstates are exponentially localized. The (dimensionless) theoretical localization length in the space of the angular momentum quantum numbers is given below, and is equal (after introducing some numerical correction factor α µ ) to the dimensionless localization time t loc [Eq. (2.16)]. We here denote it by L unlike in reference [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] and [Manos and Robnik 2013]. Therefore, an exponentially localized eigenfunction centered at m in the angular momentum space [Eq. (2.11)] has the following form:

|u n | 2 ≈ 1 L exp ( - 2|m -n| L ) , (2.15)
where u n is the probability amplitude [Eq. ( 2.13)] of the localized wavefunction ψ(θ). The argument leading to t loc in Eq. ( 2.16) originates from the observation of the dynamical localization (see e.g. [START_REF] Casati | Stochastic behavior of a quantum pendulum under a periodic perturbation[END_REF]), and is well explained in [START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF]], in case of normal diusion µ = 1, whilst for general µ we gave a theoretical argument in [Manos and Robnik 2013]. We shall denote σ = 2/L, and will later on determine the σ's directly from the individual numerically calculated eigenstate.

But where one is able to observe such phenomena (spectral statistics, namely Brody-like level spacing distribution) analogous in the quantum chaos of time-independent bound systems with discrete spectrum? To see these eects the system must have eectively nite dimension, because in the innite dimensional case we simply observe Poissonian statistics. Truncation of the innite matrix F mn in Eq. ( 2.12) in tour de force is not acceptable, even in the technical case of numerical computations, since after truncation the Floquet operator is no longer unitary.

The only way to obtain a quantum system which shall in this sense correspond to the classical dynamical system [Eqs. (2.1), (2.2) and ( 2.

3)] is to introduce a nite N -dimensional matrix, which is symmetric unitary, and which in the limit N → ∞ becomes the innite dimensional system with the Floquet operator [Eq. (2.9)]. The semiclassical limit is k → ∞ and τ → 0, such that K = kτ = constant. As it is well known [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]], for the reasons discussed above, the system behaves very similarly for rational and irrational values of τ /(4π). Such a N -dimensional model [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] will be introduced below.

Following [Manos and Robnik 2013] and references therein, we nd that the dimensionless Heisenberg time (also called break time or localization time and denoted by t loc in units of kicking period T ) is equal to the dimensionless localization length L (see Eq. ( 2.20 below). Here we repeat briey the main steps to derive this result:

L ≈ t loc = ( α µ D µ (K) τ 2 ) 1 2-µ . (2.16)
where α µ is a numerical constant to be determined empirically, and in case of normal diusion µ = 1 is close to 1/2. Since this semiclassical approach and derivation is quite important, we repeat the arguments given in [Manos and Robnik 2013].

The generalized diusion process of the SM [Eq. ( 2.

3)] is dened by Eq. (2.4). As the real physical angular momentum p and P are connected by P = pT /I we have for the variance of p the following equation:

(∆p) 2 = I 2 T 2 D µ n µ .
(2.17)

Now we argue as follows: The general rule in quantum chaos is that the quantum diusion follows the classical diusion up to the Heisenberg time (or break time, or localization time), dened as:

t H = 2π ∆E , (2.18)
where ∆E is the mean energy level spacing. In our case we have the quasienergies and ∆E = ∆ω, where ∆ω = ∆ϕ/T , and ∆ϕ is the mean spacing of the eigenphases. This might be estimated at the rst sight as ∆ϕ = 2π/N , but this is an underestimate, as eectively we shall have due to the localization only L levels on the interval [0, 2π). Therefore ∆ϕ = 2π/L and we nd:

t H = 2πT ∆ϕ = T L. (2.19)
Since T is the period of kicking, and t H is the real physical continuous time, we get the result that the discrete time [number of iterations of Eq. ( 2.3) at which the quantum diusion stops], the localization time t loc is indeed equal to the localization length in momentum space, i.e.:

t loc ≈ L.

(2.20)

Since our derivation is not rigorous, we use the approximation symbol rather than equality, in particular as the denition depends linearly on the denition of the Heisenberg time. Now the nal step: By inspection of the dynamics of the Floquet quantal system [Eqs. (2.8),(2.9)] one can see (the derivation can be found in [START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF]) that the value of the variance of the angular momentum at the point of stopping the diusion (∆p) 2 is proportional to 2 L 2 , and to achieve equality we introduce a dimensionless numerical (empirical,correction) factor α µ by writing (∆p) 2 = 2 L 2 /α µ , which on the other hand must be equal just to the classical value at stopping time t loc , namely equal to (I/T ) 2 D µ L µ . From this it follows immediately Eq. (2.16). The numerical constant α µ is found empirically by numerical calculations, for instance in the literature the case K = 5 with µ = 1 is found to be α 1 = 0.5 (however, we nd numerically α 1 = 0.45, taking into account Eq. ( 2.16) when studying the model's localization properties). Thus, we have the theoretical formula for the localization length in the case of generalized classical diusion [Eqs. (2.4),(2.17)], which we use in dening the scaling parameter Λ below.

The Floquet (Izrailev) model

The motion of the QKR [Eq. (2.8)] after one period T of the ψ wavefunction can be described also by the following symmetrized Floquet mapping, describing the evolution of the kicked rotator from the middle of a free rotation over a kick to the middle of the next free rotation, as follows:

ψ(θ, t + T ) = Û ψ(θ, t), (2.21) Û = exp ( i T 4I ∂ 2 ∂θ 2
) exp

( -i V 0 cos θ ) exp ( i T 4I ∂ 2 ∂θ 2
) .

Thus, the ψ(θ, t) function is determined in the middle of the rotation, between two successive kicks.

The evolution operator Û of the system corresponds to one period. Due to the instant action of the kick, this evolution can be written as the product of three non-commuting unitary operators, the rst and third of which correspond to the free rotation during half a period Ĝ(τ /2) = exp

( i τ 4 ∂ 2 ∂θ 2
)

, τ ≡ T /I, while the second B(k) = exp(-ik cos θ), k ≡ V 0 / , describes the kick. Like before, we have only two dimensionless parameters, namely τ and k, and K = kτ = V 0 T /I. In the case K ≡ kτ 1 the motion is well known to be strongly chaotic, for K = 10 certainly without any regular islands of stability, as mentioned, and also there are no accelerator modes, so that the diusion is normal (µ = 1). We have carefully checked that the case K = 10 is the closest to the normal diusion µ = 1 for all K ∈ [K cr , 70]. The transition to classical mechanics is described by the limit k → ∞, τ → 0 while K = const. We shall consider the regimes on the interval 3 ≤ k ≤ 20, but will concentrate mostly on the semiclassical regime k ≥ K, where τ ≤ 1.

In order to study how the localization aects the statistical properties of the quasienergy spectra, we use the model's representation in the momentum space with a nite number N of levels [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]], which we refer to as Izrailev model:

u n (t + T ) = N ∑ m=1 U nm u m (t), n, m = 1, 2, ..., N . (2.22)
The nite symmetric unitary matrix U nm determines the evolution of an N -dimensional vector, namely the Fourier transform u n (t) of ψ(θ, t), and is composed in the following way:

U nm = G nn ′ B n ′ m ′ G m ′ m (2.23) where G ll ′ = exp ( iτ l 2 /4
) δ ll ′ is a diagonal matrix corresponding to free rotation during a half period T /2, and the matrix B n ′ m ′ describing the one kick has the following form:

B n ′ m ′ = 1 2N + 1 × (2.24) 2N +1 ∑ l=1 { cos [ ( n ′ -m ′ ) 2πl 2N + 1 ] -cos [ (n ′ + m ′ ) 2πl 2N + 1 ]} × exp [ -ik cos ( 2πl 2N + 1 )]
.

The Izrailev model [Eqs. (2.22-2.24)] with a nite number of states is considered as the quantum analogue of the classical SM on the torus with closed momentum p and phase θ, where U nm describes only the odd states of the systems, i.e. ψ(θ) = -ψ(-θ), provided we have the case of the quantum resonance, namely τ = 4πr/(2N +1), where r is a positive integer, as in Eq. (2.14). The matrix (2.24) is obtained by starting the derivation from the odd-parity basis of sin(nθ) rather than the general angular momentum basis exp(inθ).

Nevertheless, we shall use this model for any value of τ and k, as a model which in the resonant and in the generic case (irrational τ /(4π)) corresponds to the classical kicked rotator, and in the limit N → ∞ approaches the innite dimensional model [Eq. ( 2 The dierence of behavior between the generic case and the quantum resonance shows up only at very large times, which grow fast with (2N + 1), as explained above. It turns out that also the eigenfunctions and the spectra of the eigenphases at nite dimension N of the matrices that we consider do not show any signicant dierences in structural behavior for the rational or irrational τ /(4π), which we have carefully checked. Indeed, although the eigenfunctions and the spectrum of the eigenphases exhibit sensitive dependence on the parameters τ and k, their statistical properties are stable against the small changes of τ and k. This is an advantage, as instead of using very large single matrices for the statistical analysis, we can take a large ensemble of smaller matrices for values of τ and k around some central value of τ = τ 0 and k = k 0 , which greatly facilitates the numerical calculations and improves the statistical signicance of our empirical results. Therefore our approach is physically meaningful. Similar approach was undertaken by Izrailev (see [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] and references therein). In Figs. we can get the rst quantitative empirical localization measure to be discussed and used later on. The new nding in [Manos and Robnik 2015] is that σ has a distribution, which is close to the Gaussian (but cannot be exactly that, because σ is a positive denite quantity). It does not depend on N and survives the limit N → ∞. Therefore also L has a distribution whose variance does not vanish in the limit N → ∞.

Dynamical (Chirikov) Localization

Following [Manos and Robnik 2013] and [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] we introduce another measure of localization.

For each N -dimensional eigenvector of the matrix U nm the information entropy is:

H N (u 1 , ..., u N ) = - N ∑ n=1
w n ln w n , (2.25) where w n = |u n | 2 , and

∑ n |u n | 2 = 1.
In case of the random matrix theory being applicable to our system [Eqs. (2.21) and (2.22-2.24)],

namely the Circular orthogonal ensemble (COE) (or GOE), due to the isotropic distribution of the eigenvectors of a COE of random matrices, we have the probability density function of |u n | on the interval [0, 1]:

w N (|u n |) = 2Γ(N/2) √ πΓ((N -1)/2) (1 -|u n | 2 ) (N -3)/2 .
(2.26)

It is easy to show that in the limit N → ∞ this becomes a Gaussian distribution:

w N (|u n |) = √ 2N π exp ( - N |u n | 2 2
) , (2.27) and the corresponding information entropy [Eq. (2.25)] is equal to:

H GOE N = ψ ( 1 2 N + 1 ) -ψ ( 3 2
) ln

( 1 2 N a ) + O(1/N ), (2.28)
where a = 4 exp(2-γ) ≈ 0.96, while ψ is the digamma function and γ the Euler constant ( 0.57721...). For a uniform distribution over M states w n = 1/M we get H N ≈ log M , and thus M ≈ exp(H N ). Thus, we get the insight that the correct measure of localization must be proportional to exp(H N ), but properly normalized, such that in case of extendedness (GOE/COE) it is equal to N .

Therefore the entropy localization length l H is dened as:

l H = N exp ( H N -H GOE N
) .

(2.29) Indeed, for entirely extended eigenstates l H = N . Thus, l H can be calculated for every eigenstate individually. However, all eigenstates, while being quite dierent in detail, are exponentially localized, and thus statistically very similar. Therefore, in order to minimize the uctuations one uses the mean localization length d ≡ l H , which is computed by averaging the entropy over all eigenvectors of the same matrix (or even over an ensemble of similar matrices):

d ≡ l H = N exp ( H N -H GOE N
) .

(2.30)

The localization parameter β loc is then dened as:

β loc = d N ≡ l H N .
(2.31)

The parameter that determines the transition from weak to strong quantum chaos is neither the strength parameter k nor the localization length L, but the ratio of the localization length L to the size N of the system in momentum p: .32) where L ≈ t loc , the theoretical localization length Eq. (2.16), was derived in [Manos and Robnik 2013].

Λ = L N = 1 N ( α µ D µ (K) τ 2 ) 1 2-µ , ( 2
Λ is the scaling parameter of the system. The relationship of Λ to β loc is discussed e.g. in [START_REF] Casati | Scaling behavior of localization in quantum chaos[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] and further developed in [Manos and Robnik 2013] (Sect. VII).

Here we just summarize by showing the empirical scaling property of Λ versus β loc in Fig. 2.6, where the approximate analytical description is given by the function:

β loc = γΛ 1 + γΛ , γ = 4.04, (2.33)
which is similar to the scaling law in in [START_REF] Casati | Scaling behavior of localization in quantum chaos[END_REF]], but not the same. Namely, the value γ = 4.04 diers somewhat from γ ≈ 3.2 in [START_REF] Casati | Scaling behavior of localization in quantum chaos[END_REF], where β loc is plotted versus x ≈ 4Λ.

In producing this plot we have used Eqs. (2.16),(2.32) for Λ, which is just a rough estimate. Indeed, as we shall show below, following [Manos and Robnik 2015], Λ is not a number in a given system, at xed K, k and N , but has a distribution, whose reciprocal is approximately Gaussian distributed, and Eq. (2.16) is just a rough estimate of the mean value of Λ. Therefore we should not be surprised any more to see large uctuations in the scaling law of Fig. 2.6, an observation entirely unexplained so far, but claried in [Manos and Robnik 2015]. The statistical properties of the localization measure will be discussed below.

β loc Λ K=5 K=7 K=10 K=12 K=14 K=17 K=20 K=25 K=30 K=35 f(Λ)=γΛ/(1+γΛ), γ=4.04

Level spacing distribution: P (S), W (S) and U (S)

To study the eigenvalue statistics of quantum Floquet systems and quantum maps one considers the eigenphases ϕ n ∈ [0, 2π) dened by λ n = e iϕn . In such case the spectral unfolding procedure is very easy, as the mean level density is N/(2π), i.e. the mean level spacing is 2π/N . The histogram of the level spacing distribution P (S) is the distribution of the spacings S n := N 2π (ϕ n+1 -ϕ n ) with n = 1, ..., N and ϕ N +1 := ϕ 1 , in the bins of certain suitable size ∆S. The factor N/2π ensures that the average of all spacings S n is 1, and thus P (S) is supported on the interval [0, N ], and its upper limit goes to ∞ when N → ∞.

The cumulative distribution W (S), or integrated level spacing distribution, preserving the full accuracy of all numerical eigenvalues/spacings, useful especially when the number of levels N is small, is dened as:

W (S) = ∫ S 0 P (x)dx ≡ #{n|S n ≤ S} N .
(2.34)

Finally, we shall use also the so-called U -function (see Appendix A in [Manos and Robnik 2015]):

U (S) = 2 π arccos √ 1 -W (S).
(2.35)

The U -function has the advantage that its expected statistical error δU is independent of S, being constant for each S and equal to δU = 1/(π √ N s ), where N s is the total number of objects in the W (S) distribution. The numerical pre-factor 2/π in Eq. (2.35) is determined in such a way that U (S) ∈ [0, 1] when W (S) ∈ [0, 1]. The U -function is an excellent and rened criterion used to assess the goodness of the t of the theoretical level spacing distribution.

The important special level spacing distributions that we are using are the following: X The Poisson distribution:

P Poisson (S) = e -S , W Poisson (S) = 1 -e -S .
(2.36) X The COE or GOE distribution:

P COE (S) ≈ P Wigner (S) = π 2 exp ( - π 4 S 2 ) , W COE (S) = 1 -exp ( -πS 2 4
) .

(2.37) X The Brody distribution, [START_REF] Brody | Randommatrix physics: spectrum and strength uctuations[END_REF]:

P BR (S) = C 1 S β exp ( -C 2 S β+1
) , (2.38) where the two parameters C 1 and C 2 are determined by the two generic normalization conditions that must be obeyed by any P (S):

∫ ∞ 0 P (S)dS = 1, S = ∫ ∞ 0 SP (S)dS = 1, (2.39)
thus with S = 1 being the mean distance between neighboring levels (after unfolding). Hence:

C 1 = (β + 1)C 2 , C 2 = [ Γ ( β + 2 β + 1
)] β+1 , (2.40) where Γ(x) denotes the Gamma function. In the strongly localized regime β = 0 we observe Poissonian statistics while in the fully chaotic one β = 1 and the RMT applies. The Brody cumulative level spacing distribution is:

W BR (S) = 1 -exp(-C 2 S β+1
).

(2.41) X Izrailev distribution: In [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF][START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]], Izrailev suggested the following distribution in order to describe the intermediate statistics, i.e. the non-integer β in the following PDF could be associated with the statistics of the quasienergy states with chaotic localized eigenfunctions, also approximating the level spacing distribution arising from the Dyson COE joint probability distribution [START_REF] Porter | Statistical theory of spectra: uctuations[END_REF]]:

P IZ (S) = A ( 1 2 πS ) β exp [ - 1 16 βπ 2 S 2 - ( B - 1 4 πβ ) S
] ,

(2.42)

where again the two parameters A and B are determined by the two normalization conditions 1 = S = 1 given above.

Of course, we must be fully aware of the fact that both, Brody and Izrailev distributions, are approximations. It is clear that at β BR = 1 we get precisely Wigner surmise [Eq. (2.37)], which is the exact GOE only for two-dimensional Gaussian random matrices, and thus only an (excellent) approximation for the innite dimensional GOE case. Indeed, if we try to t the exact innite dimensional GOE level spacing distribution with the Brody distribution, we do not get β BR = 1, but instead β BR = 0.953, see [START_REF] Brody | Randommatrix physics: spectrum and strength uctuations[END_REF]]. Also, we should mention that Izrailev et al. have published an improved distribution function [Casati, Izrailev, and Molinari 1991], which we have also tested, and is dened by:

P new IZ (S) = AS β (1 + BβS) f (β) exp[- π 2 16 βS 2 - π 2 (1 - β 2 S)] (2.43)
where f (β) =

2 β (1-β 2 ) β
-0.16874 and A, B are the normalization parameters. We found (see below) that in our applications it is even worse than the original version [Eq. (2.42)].

Lévy stable distribution

The physical origin and relevance of the Lévy stable distribution to this kind of problems, like the SM, is well summarized in e.g. [START_REF] Zaslavsky | The Physics of Chaos in Hamiltonian Systems[END_REF][START_REF] Zaslavsky | Self-similarity, renormalization, and phase space nonuniformity of Hamiltonian chaotic dynamics[END_REF][START_REF] Klafter | Lévy statistics in a Hamiltonian system[END_REF][START_REF] Geisel | Accelerated diusion in Josephson junctions and related chaotic systems[END_REF][START_REF] Zaslavsky | Hierarchical structures in the phase space and fractional kinetics: I. Classical systems[END_REF]. In probability theory, an α-Lévy skew stable distribution is a four parameter family of continuous probability distributions, characterizing the location, scale, skewness and kurtosis. Following [START_REF] Nolan | Stable Distributions -Models for Heavy Tailed Data[END_REF]],

for a random variable X with distribution function F (x), the characteristic function is dened by ϕ

(u) = E exp(iuX) = ∫ ∞ -∞ exp(iux)dF (x)
. Then, a random variable X is stable if and only if X δ = aZ + b, with a > 0, b ∈ R and Z is a random variable with characteristic function:

E exp(iuZ) = { e -|u| α [1-iβ tan πα 2 (signu)] , α = 1 e -|u|[1+iβ tan 2 π (signu)] log |u| , α = 1 , (2.44)
where 0 < α ≤ 2 and -1 ≤ β ≤ 1 (the symbol δ = indicates that both expressions have the same probability law). We then adopt the parametrization k = 0, S(α, β, γ, δ; 0) for which the random variable X given by:

X δ = { γ(Z -β tan( πα 2 )) + δ, α = 1 γZ + δ, α = 1 , ( 2.45)
has characteristic function: .46) where Z = Z(α, β) is dened as described in Eq. (2.44), α ∈ (0, 2] is the index of stability or characteristic exponent, β ∈ [-1, 1] the skewness parameter, γ > 0 the scale parameter and δ ∈ R location parameter. For the ts with data we used the Stable Distribution package of Mathematica [START_REF] Rimmer | Stable Distributions in Mathematica[END_REF]. Two important special cases are the Gaussian distribution with α = 2 and the Cauchy-Lorentz with α = 1 which are the only ones with an explicit closed formula.

S(α, β, γ, δ; 0) ≡ E exp(iuX) = { e iuδ-γ α |u| α (1+iβ(-1+|uγ| 1-α )sign(u) tan( πα 2 )) , α = 1 e iuδ-γ|u|(π+2iβ log(|uγ|))sign(u) π , α = 1 , ( 2

Results

In Fig. 2.7, we present the probability distribution function of the momentum P after n = 5000 iterations (black lled circles) and the ts with the α-Lévy stable distribution (solid blue line) for a sample of K-values associated with the principal peaks (presence of accelerator modes of period 1) of the Fig. 2.3, i.e., K = 6.50,11.90,13.20,18.95. In Fig. 2.7(a),(d),(e),(f ) α is generally not equal to 2 and two K-values without accelerator modes, i.e., K = 7.0, 10.0, are shown in Fig. 2.7(b),(c) where α is equal to 2. Here, and for comparison reasons with the best-t function depicted in the gures, the Gaussian distribution (red dashed line) is derived by the S(α, β, γ, δ; 0) probability distribution by setting α = 2 and keeping all the remaining parameters the same as given by the ts generally for α = 2. The total number of initial conditions is ≈100000 (314 × 314) on a uniform grid in the entire phase plane (θ, P ) ∈ [0, 2π] × [0, 2π]. For larger K-values with accelerator modes of period 1, their eect in the diusion process is becoming gradually weaker, as can also be seen by the decay of the µ value in Fig. 2.3, and the probability distributions tend to Gaussian-like ones. We have also calculated the χ 2 -test for all the density distributions of Fig. 2.7. However, the derived quantitative information turns out, not to dier much to what the panels show. Let us also stress that the goodness of the several ts is not the main point here but rather the deviation from the Gaussian prole which is evident when the eect of the accelerator mode is relatively strong as shown in panel (a) and expressed by α parameter value far from 2. Moreover, due to our main motivation, i.e., to focus on the relevant quantum time scales, we do not expect to capture very accurately the tails, which are expected to follow the theoretical curve at much larger times. More examples and the details regarding the whole set of parameters can be found in [Manos and Robnik 2014]. Note that, for K < 2π the t was done for an ensemble in a cell around the origin (θ, P ) = (0, 0) instead of a grid uniform in the entire phase space, in order to exclude the data coming from islands of stability whose momenta do not diuse at all and mix up the distribution. The probability distribution function of the momentum P after n = 5000 iterations (black lled circles) and the ts with the stable α-Lévy distribution for a sample of K-values aswith the principal peaks (presence of accelerator modes of period 1) of the Fig. 2.3, i.e., K = 6.5, 11.9, 13.2, 18.95 (a),(d),(e),(f ) where α is generally not equal to 2 and two K-values without accelerator modes, i.e., K = 7.0, 10.0 (b),(c) where α is equal to 2. Here, the Gaussian distribution (red dashed line) for comparison reasons with the best-t function, depicted in the gures, is derived by the S(α, β, γ, δ; 0) probability distribution by setting α = 2 and keeping all the rest parameters same as the given by the ts in general α = 2. The total number of initial conditions is ≈100000 (314 × 314) on a uniform grid in the entire phase plane (θ, P ) ∈ [0, 2π] × [0, 2π]. For larger K-values with accelerator modes of period 1, their eect in the diusion process is becoming gradually weaker, as can be seen also by the decay of the µ value in Fig. 2.3, and the probability distributions tend to Gaussian-like ones.

In Fig. 2.8(a) we present the outcome of the calculation of the GALI 2 on the whole phase plane (θ, P ) ∈ [0, 2π]×[0, 2π] for K = 3.1, for 500×500 initial conditions uniformly distributed. Each initial condition is colored according to the color scale seen at the right side of the panel. For chaotic orbits, having small GALI 2 value (≈ 10 -8 ) are colored black, while the yellow color corresponds to regular motion, found here to be ≈ 13.52% of the whole plane, with high -close to zero -values (the color bar is in a logarithmic scale). Thus, we can clearly identify even tiny regions of regular motion which are not easily seen in the phase space portraits given often by simple Poincaré surface of sections.

Having located the stable region of the phase space, the next point of interest is to distinguish among them, those that are due to accelerator modes from those due to islands of stability.

The distinction can be eciently achieved with the use of the diusion exponent µ color-plot in Fig. 2.8(b), where we rst consider again a grid of 500 × 500 cells (on the entire phase space (θ, P ) ∈ [0, 2π] × [0, 2π]) with 50 × 50 initial conditions in each and evolve all of them together for n = 5000 iterations. Then, for each ensemble of each cell separately, we calculate numerically the diusion coecient D µ (K) as a function of the iterations n and perform a t procedure just like in Fig. 2.2 to calculate the diusion exponent µ [represented in the color bar of Fig. 2.8(b)] that characterizes the dierent kind of diusion of this small area. Depending on the relative location of each ensemble one may expect to nd: (i) normal diusion (µ = 1) inside chaotic regimes without the presence of accelerator modes, (ii) subdiusion (0 < µ < 1) inside islands of stability, (iii) superdiusion (1 < µ < 2) inside chaotic regimes with the presence of accelerator modes in the phase space and (iv) ballistic transport (µ ≈ 2) inside and in the very close vicinity of accelerator modes.

It turns out that the stable regions around (π/2, 0), (3π/2, 0), (π/2, π), (3π/2, π) and (π/2, 2π), (3π/2, 2π) are indeed islands of stability since their diusion exponent µ is smaller than 1. The remaining tiny stable areas correspond to stable higher period accelerator modes with µ ≈ 2.

Here we manage to locate the accelerator modes of higher periods 2,3,4,... which in general are not so easy to be calculated analytically. In

Analysis of the level spacing distribution of numerical spectra

For the numerical calculations and results regarding the spacing distributions P (S) (and W (S)) for the eigenphases ϕ j , we have considered a range of 41-values of the quantum parameter k (= 2, 3, ..., 42) keeping xed the classical parameter K = 7, where the phase space is fully chaotic (see [Manos and Robnik 2013]). In order to ameliorate the statistics we considered a sample of 161 matrices U nm of size N = 398 (≈ 64, 000 elements), in a similar manner as in [Izrailev 1988]) with slightly dierent values of k (with the step size ∆k = ±0.00125 k). For some samples we reached up to 641 matrices U nm of size N = 398 acquiring qualitatively the same results.

For the ensemble of M = 641 matrices of size N = 398, in case K = 7 and k = 11, using the χ 2 best tting procedure (described in more detail below) we found β BR = 0.421 and all three representations clearly show excellent agreement with the best tting Brody distribution. In Fig. 2 always the outer one) PDFs by using the U -function and W -distribution. Thus, in case of the ideal tting the data would lie on the abscissa. In this case, based on the P (S) we show the histogram. It is seen that Brody distribution is better tting the data than the Izrailev distribution. Since the deviations are really small, the statistical signicance very high, we plot in Fig. 2.9b the dierences U (data) -U (Brody/Izrailev) versus U BR . Thus if data are on the abscissa the agreement is perfect. As can be seen, the deviations from that are really small, and clearly smaller for Brody. In Fig. 2.9c we show the ne dierences of W (data) -W (Brody/Izrailev) versus W BR , and again we clearly see that Brody is better. Finally, in Fig. 2.9d we show the same thing as in Fig. 2.9c, except not versus W but versus S instead.

Residues and χ 2 test:

In the best tting procedure we have calculated both the residues and the χ 2 as follows:

X PDFs residues: R PDFs =

∑ N i=1 (P BR,IZ (i) -data(i)) 2 . X χ 2 : χ 2 PDFs = ∑ N i=1 (P BR,IZ (i)-data(i)) 2 P BR,IZ (i)
.

In Fig. localization measures vary in the semiclassical limit of the increasing value of the quantum parameter k, at xed classical parameter K = 10. Indeed, the theoretical estimate of L in Eq. ( 2.16), at xed K, and remembering k = K/τ , shows that approximately the mean value of the localization length should increase quadratically with k, or equivalently, the slope σ should decrease inversely quadratically with k. This prediction is observed, and is demonstrated in Fig. 2.16. This is also in agreement with the prediction based on the tight-binding approximations in reference [Kottos et al. 1996] [Eq. (6)]. We give, in [Manos and Robnik 2015] (Table I), the mean slope σ and the standard deviation of σ, as well as the mean value of the related quantity 2/l H and its standard deviation for various k = k 0 = 3, 4, 5, . . . , 19, for each of them taking seven nearby values of k, namely k = k 0 ± jδk, where j = 0, 1, 2, 3 and δk = 0.00125, for matrices of dimension N = 3000.

Each histogram for all k 0 was tted with the Gaussian distribution and then the mean values and the standard deviations were extracted. All four quantities decrease to zero with increasing k, meaning In the former case the behavior is roughly as 1/k 2 , in agreement with the theoretical estimate 1/k 2 of Eq. ( 2.16), and in the latter case also like 1/k 2 , surely not as the theoretical estimate 1/k based on the Lyapunov exponents method in [Kottos et al. 1996] [Eq. (9)].

that in the semiclassical limit the localization lengths monotonically increase to innity, so that in this limit we have asymptotically extended states (no localization), and their standard deviation also goes to zero as 1/k 2 , which is dierent from the tight-binding approximations in reference [START_REF] Kottos | Scaling properties of Lyapunov spectra for the band random matrix model[END_REF] [Eq. ( 9)].

Next we studied how does the distribution of the localization measure σ behaves as a function of the dimension N of the Izrailev model Eqs. (2.22-2.24). Since in the limit N → ∞ the model converges to the innitely dimensional quantum kicked rotator, we would at rst sight expect that following the Shepelyansky picture [START_REF] Shepelyansky | Localization of quasienergy eigenfunctions in action space[END_REF]] σ should converge to its asymptotic value, which is sharply dened in the sense that the variance of the distribution of σ goes to zero inversely with N . Namely, at xed K and k Shepelyansky reduces the problem of calculating the localization length to the problem of the nite time Lyapunov exponents of the approximate underlying nite dimensional Hamilton system with dimension 2k. The localization length is then found to be equal to the inverse value of the smallest positive Lyapunov exponent. In our case, the dimension of the matrices N of the Izrailev model plays the role of time. As it is known, and analyzed in detail in [Manos and Robnik 2015], the nite time Lyapunov exponents have a distribution, which is almost Gaussian, and its variance decays to zero inversely with time. Thus on the basis of this we would expect that the variance of σ decays inversely with N . However, this is not what we observe. In [Manos and Robnik 2015] (Table II) we clearly see that at constant K = 10 and k = 10 the mean value of σ is constant and obviously equal to its asymptotic value of N = ∞, while the variance of σ does not decrease with N , as 1/N , but is constant instead, independent of N . This is in disagreement with the banded-matrix models of the tight-binding approximations and thus disagrees with the Eq. ( 9) of reference [START_REF] Kottos | Scaling properties of Lyapunov spectra for the band random matrix model[END_REF],

and also disagrees with the Shepelyansky picture. The reason is that the associated Shepelyansky's

Hamilton system is only approximate construction, because with increasing N the matrix elements of the Floquet propagator (matrix) outside the diagonal band of width 2k become important, and thus the dimension of the Hamilton system cannot be considered nite, constant and equal to 2k, but increases with N . As a consequence we have the constant value of the variance of σ, and thus constant variance of the localization length L = 2/σ, and therefore the localization length has a distribution with nonvanishing variance even in the limit N = ∞. This is precisely the reason why the semiclassical prediction of the localization length in Eq. ( 2.16) fails in detail and we nd strong uctuations in the plot of L against the 2/σ of Figs. 2.6 and 2.13. The proper theory of the localization length must predict its distribution rather than just its approximate mean value. The problem of quantum or dynamical localization is related to the Anderson localization model, within the framework of the tight-binding approximation, with hopping transitions between the nearest neighbors only [START_REF] Fishman | Chaos, Quantum Recurrences, and Anderson Localization[END_REF].

In [Manos and Robnik 2015], we have numerically analyzed the behavior of the nite time Lyapunov exponents for a classical Hamilton system exemplied by the SM, following Fujisaka [START_REF] Fujisaka | Statistical Dynamics Generated by Fluctuations of Local Lyapunov Exponents[END_REF]] and Ott [START_REF] Ott | Chaos in Dynamical Systems[END_REF]] and also for the 2 × 2 random transfer matrices of the tight-binding approximation to describe the Anderson localization. In both cases we have shown that the distribution of the positive Lyapunov exponent is excellently described by a Gaussian distribution, whose mean value converges with time to the asymptotic value of the innite time, while the variance decays inversely with time t (the number of iterations in the case of the SM), and n, the number of random matrices in the product. The latter are random unimodular transfer matrices of the tight-binding approximation, of the form:

T n = ( W -1 1 0 ) (2.47)
where W = E -E 0 n is drawn from a distribution, dened by a given model. E is the eigenenergy of the system, and E 0 n is the uctuating on-site potential. We have tested three quite dierent distributions for W , namely Gaussian, box distribution and the Cauchy-Lorentz distribution [Manos and Robnik 2015], and found that the shape of the distribution of the positive Lyapunov exponent for any n (=number of matrices in the product) depends very weakly on the overall shape of the W -distribution, while the mean value and the variance depend only on the variance of W . Indeed, the evidence for the predicted decay of the variance of the nite time Lyapunov exponents is overwhelming, as shown in Fig. 2.17, where we plot the standard deviation as a function of time in log-log plot, showing that it decays inversely with the square root of time.

In the context of our Izrailev model the dimension N of the matrix plays the role of time. The width of the diagonal band is equal to 2k. to the asymptotic value, while the variance does not decay to zero, but rather remains constant, independent of N as clearly demonstrated in [Manos and Robnik 2015]. From this we conclude that even in the limit N → ∞ the localization length has a certain distribution with nonvanishing variance, or more precisely, its inverse (the slope σ) has an almost Gaussian distribution with nonvanishing variance. We believe that this is the cause of the strong uctuations observed for example in Figs. 2.6 and 2.13(a).

Summary and outlook

We have reviewed some recent results [Manos and Robnik 2013;[START_REF] Batisti¢ | The intermediate level statistics in dynamically localized chaotic eigenstates[END_REF]; found, in agreement with previous results, e.g. see [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]] and references therein. However, even after a great numerical eort in maximally improving the statistical relevance of this scaling law, large uctuations around the mean value have been observed. Also, we have shown, that the Brody distribution describes the energy level spacing distribution very well, and the spectral Brody parameter β BR , determining the level spacing distribution, was found to be linearly related to β loc .

The main conclusion of our analysis is the empirical fact based on numerical computations of the eigenfunctions of the N -dimensional Izrailev model, that the localization length has a distribution with nonvanishing variance not only for nite N , but even in the limit N → ∞. This is the reason, we believe, for the strong uctuations in the scaling laws which involve the empirical localization measures and the theoretical semiclassical value of the localization length. In the Shepelyansky picture [START_REF] Shepelyansky | Localization of quasienergy eigenfunctions in action space[END_REF]] this might seem to be a contradiction, but the resolution of the puzzle is that in the limit of large N the nite dimensional Hamiltonian system extracted from the Floquet propagator of the quantum kicked rotator is not good enough, and therefore the matrix elements outside the main diagonal band of width 2k play a role, even if they are small, but nevertheless plentiful, making the Hamiltonian system eectively innite dimensional, with innitely many Lyapunov exponents. This nding is a challenge for the improved semiclassical theory of the localization length, to derive and explain the discovered distribution function.

On the other hand, the simple model of the Anderson localization based on the tight-binding approximation, with only nearest neighbor interactions, described by the product of 2 × 2 unimodular matrices, has a nite dimension, as the transfer matrices are exactly two-dimensional, and therefore the variance vanishes in the limit of large times n (number of matrices in the product) as 1/n. The same conclusion applies to such a model with a nite number of interacting neighbors. Indeed, according to e.g. [Kottos et al. 1996] the variance of σ should vanish as V ar(σ) ∝ 1/(N k 2 ), but our work shows that in the quantum kicked rotator this is not observed: the variance does not depend on N , and decays with k faster than 1/k 2 , namely as 1/k 4 . Thus, we found some important dierences between the dynamical localization in the quantum kicked rotator and the Anderson tight-binding model of localization, and the Shepelyansky picture, which rest upon the banded matrix models with rigorously nite bandwidth.

Perspectives

The problem of calculating the distribution of the localization length (or its inverse) in the semiclassical framework is open for the future work. Also, the theoretical derivation of the Brody distribution to explain the level spacing distribution of the energies [START_REF] Batisti¢ | Dynamical localization of chaotic eigenstates in the mixed-type systems: spectral statistics in a billiard system after separation of regular and chaotic eigenstates[END_REF]Batisti¢ and Robnik 2013a;Batisti¢ and Robnik 2013b] in time-independent systems, and of the quasienergies [START_REF] Izrailev | Quantum localization and statistics of quasienergy spectrum in a classically chaotic[END_REF]Manos and Robnik 2013;[START_REF] Batisti¢ | The intermediate level statistics in dynamically localized chaotic eigenstates[END_REF]

+D=FJAH !

Complex neuronal activity patterns in the human brain

Background and introduction

Synchronization of oscillations is a generic mechanism in animate and inanimate systems (see e.g.

[A. [START_REF] Pikovsky | Synchronization: A universal concept in nonlinear sciences[END_REF]). In fact, oscillators of qualitatively dierent type may share fundamental synchronization mechanisms. Synchronization processes may occur within as well as between dierent systems of the human body, e.g. between heartbeat intervals and respiratory cycles [START_REF] Bartsch | Phase transitions in physiologic coupling[END_REF]. Neuronal synchronization processes are relevant under normal as well as abnormal conditions. A number of brain disorders are associated with abnormal neuronal synchrony, for example Parkinson's disease (see e.g. [START_REF] Lenz | Single unit analysis of the human ventral thalamic nuclear group -Tremor-related activity in functionally identied cells[END_REF][START_REF] Nini | Neurons in the globus pallidus do not show correlated activity in the normal monkey, but phase-locked oscillations appear in the MPTP model of parkinsonism[END_REF][START_REF] Hammond | Pathological synchronization in Parkinson's disease: networks, models and treatments[END_REF]), tinnitus (see e.g. [START_REF] Ochi | Eects of quinine on neural activity in cat primary auditory cortex[END_REF][START_REF] Eggermont | Maladaptive Neural Synchrony in Tinnitus: Origin and Restoration[END_REF][START_REF] Elgoyhen | Tinnitus: perspectives from human neuroimaging[END_REF]), migraine (see e.g. [START_REF] Angelini | Steady-State Visual Evoked Potentials and Phase Synchronization in Migraine Patients[END_REF][START_REF] Bjørk | Quantitative EEG Power and Asymmetry Increase 36 h Before a Migraine Attack[END_REF]) and epilepsy (see e.g. [START_REF] Wong | Cellular basis of neuronal synchrony in epilepsy[END_REF]). Neuronal dynamics and, in particular synchronization processes, crucially depend on the patterns and types of neuronal connections [START_REF] Sporns | Networks of the Brain[END_REF]]. For instance, according to computational studies it makes a signicant dierence whether neurons interact through gap-junctions or synapses. This is relevant for the emergence of dierent kinds of synchronization patterns [START_REF] Esfahani | Stimulus-dependent synchronization in delayedcoupled neuronal networks[END_REF] and epileptic seizures [START_REF] Volman | Gap junctions and epileptic seizurestwo sides of the same coin?[END_REF].

Connectivity and function are strongly connected and may undergo plastic changes throughout the life course [START_REF] Hübener | Neuronal plasticity: beyond the critical period[END_REF]. The timing pattern of neuronal activity may strongly determine the strength of neuronal connections [START_REF] Bliss | Long-lasting potentiation of synaptic transmission in the dentate area of the anaesthetized rabbit following stimulation of the perforant path[END_REF]. Spike-timing-dependent plasticity (STDP) is a pivotal mechanism by which neurons adapt the strength of their synapses to the relative timing of their action potentials (see e.g. [H. [START_REF] Markram | Spike-timing-dependent plasticity: A comprehensive overview[END_REF]).

Based on seminal experimental studies [START_REF] Markram | Regulation of synaptic ecacy by coincidence of postsynaptic APs and EPSPs[END_REF]] a series of computational studies focused on how adaptive coupling and activity dependent synaptic strength inuence the collective neuronal dynamics was performed. In the presence of STDP a plethora of qualitatively dierent stable dynamical regimes emerge [START_REF] Zhang | Synaptic plasticity induced transition of spike propagation in neuronal networks[END_REF]]. In fact, multistability (the coexistence of multiple attractors in the phase space) is a typical feature of neuronal networks and oscillator networks equipped with STDP. Multistability was found in dierent neural network models comprising different STDP models, e.g. in phase oscillator networks with both symmetric and asymmetric phase dierence-dependent plasticity, a time continuous approximation of STDP [START_REF] Maistrenko | Multistability in the Kuramoto model with synaptic plasticity[END_REF] as well as in phase oscillator networks with STDP and in dierent types of neuronal networks with STDP [Pster and P. Tass 2010; O.V. [START_REF] Popovych | Depending on the T s value, more (or less) ON-cycles may be administered within a xed time interval[END_REF]Ebert, Hauptmann, and P. Tass 2014] and other types of neural network models (e.g. [START_REF] Proix | Individual brain structure and modelling predict seizure propagation[END_REF]] and references therein).

A number of computational studies were dedicated on desynchronizing synchronized ensembles or networks of oscillators or neurons [START_REF] Winfree | The Geometry of Biological Time[END_REF][START_REF] Tass | A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations[END_REF] [START_REF] Silchenko | Impact of acoustic coordinated reset neuromodulation on eective connectivity in a neural network of phantom sound[END_REF] as well as abnormal cross-frequency coupling [START_REF] Adamchic | Psychometric Evaluation of Visual Analog Scale for the Assessment of Chronic Tinnitus[END_REF] within a tinnitus-related network of brain areas (see [Manos, Zeitler, and P. A. Tass 2018a] and references therein).

So far, the pre-clinical [P. [START_REF] Tass | Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation[END_REF]P. A. Tass, Qin, et al. 2012a] and clinical [P. [START_REF] Tass | Counteracting tinnitus by acoustic coordinated reset neuromodulation[END_REF]Adamchic, Hauptmann, et al. 2014] [START_REF] Tass | Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation[END_REF]P. Tass, Adamchic, et al. 2009;P. A. Tass, Qin, et al. 2012a;Adamchic, Hauptmann, et al. 2014].

Investigation goals, hypotheses and protocols

In this chapter, we set out to investigate the impact of the CR stimulation frequency and intensity on the eects during stimulus delivery (so-called acute eects ), on transient eects emerging directly after cessation of stimulation (so-called acute after-eects ), and on eects outlasting cessation of stimulation (so-called sustained after-eects ). The ultimate goal of this study is to improve the calibration of CR stimulation, in particular, by providing computationally generated predictions that can be tested in subsequent pre-clinical and clinical studies. The computational study presented here is organized around three hypotheses: A particular closed loop embodiment of CR stimulation, periodic stimulation with demandcontrolled length of high-frequency pulse train, is basically a time-shifted entrainment of the dierent neuronal subpopulations with stimulus intensities adapted to the amount of undesired synchrony. Accordingly, the duration of a stimulation cycle was selected to be reasonably close to the mean period of the synchronized neuronal oscillation. In STDP-free networks of Kuramoto [START_REF] Kuramoto | Chemical oscillations, waves, and turbulence[END_REF]] and FitzHugh-Nagumo [START_REF] Fitzhugh | Impulses and Physiological States in Theoretical Models of Nerve Membrane[END_REF] X Hypothesis #2: Dierent embodiments of CR stimulation may dier with respect to eect size and robustness. In a series of computational studies (see e.g. [P. A. Tass 2003b; P. A.

Tass 2003a] and references therein) and in all pre-clinical [P. [START_REF] Tass | Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation[END_REF]P. Tass, Adamchic, et al. 2009] and clinical studies [P. A. Tass, Qin, et al. 2012b;[START_REF] Adamchic | Psychometric Evaluation of Visual Analog Scale for the Assessment of Chronic Tinnitus[END_REF][START_REF] Damchic | Linking the Tinnitus Questionnaire and the subjective Clinical Global Impression: which dierences are clinically important? In[END_REF]Adamchic, Hauptmann, et al. 2014;[START_REF] Adamchic | Reversing pathologically increased EEG power by acoustic coordinated reset neuromodulation[END_REF][START_REF] Silchenko | Impact of acoustic coordinated reset neuromodulation on eective connectivity in a neural network of phantom sound[END_REF]] performed so far, CR was applied either with xed sequences or rapidly varying sequences (RVS), where the sequence was randomly varied from cycle to cycle. In a recent computational study, it was shown that at intermediate stimulation intensities the CR-induced anti-kindling eect may signicantly be improved by CR with slowly varying sequences (SVS), i.e. by appropriate repetition of the sequence with occasional random switching to the next sequence [START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF]].

However, this study was not performed for a larger range of CR stimulation frequencies. By denition, SVS CR stimulation features signicantly more periodicity of the stimulus pattern.

Accordingly, the dependence of resonance and/or anti-resonance eects on the CR stimulation frequency might be more pronounced for SVS CR as opposed to RVS CR.

X Hypothesis #3: Pronounced acute eects might provide a necessary, but not sucient condition for pronounced sustained after-eects. In a pre-clinical study in Parkinsonian monkeys with CR-DBS delivered at an optimal and a less favorable intensity, it was shown that long and pronounced acute therapeutic after-eects coincide with long-lasting, sustained after-eects [P. A. Tass, Qin, et al. 2012a]. However, according to computational studies the relationship between acute after-eects and sustained long-lasting eects might be more involved, at least for particular parameter combinations [START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF].

Related to these hypotheses, to assess the robustness of CR stimulation against initial network conditions we performed our numerical simulations for dierent network initializations. Here, we did not systematically vary the stimulation duration. Rather, based on a pre-series of numerical simulations, we here used a xed stimulation duration that is reasonably short, but nevertheless enabled to robustly achieve an anti-kindling for properly selected values of stimulation frequency and intensity. In fact, our goal was to nd stimulation parameters enabling short, but notwithstanding eective CR stimulation. Keeping the stimulation duration at moderate levels may be benecial for applying the CR approach to dierent invasive as well as non-invasive stimulation modalities. For instance, standard high frequency (HF) DBS, i.e. permanent electrical high-frequency pulse train stimulation delivered to dedicated target areas through implanted depth electrodes, used for the treatment of, e.g. Parkinson's disease [START_REF] Benabid | Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus[END_REF][START_REF] Krack | Five-Year Follow-up of Bilat-eral Stimulation of the Subthalamic Nucleus in Advanced Parkinson's Disease[END_REF][START_REF] Deuschl | A Randomized Trial of Deep-Brain Stimulation for Parkinson's Disease[END_REF] may cause side eects. If side eects are caused by stimulation of non-target tissue, they may be reduced by adapting the spatial extent of the current spread to the target's anatomical borders by appropriate electrode designs as introduced, e.g. by [START_REF] Martens | Spatial steering of deep brain stimulation volumes using a novel lead design[END_REF][START_REF] Dijk | A novel lead design enables selective deep brain stimulation of neural populations in the subthalamic region[END_REF][START_REF] Bour | Directional Recording of Subthalamic Spectral Power Densities in Parkinson's Disease and the Eect of Steering Deep Brain Stimulation[END_REF], in particular, to spatially tailor stimuli by means of directional DBS (see e.g. [START_REF] Contarino | Directional steering[END_REF]). However, some side eects may, at least partly, be caused by stimulating the target region itself (see e.g. [START_REF] Jahanshahi | Parkinson's Disease, the Subthalamic Nucleus, Inhibition, and Impulsivity[END_REF]]. Accordingly, no matter how precisely stimuli are delivered to DBS targets, the amount of stimulation should be decreased as much as possible.

HF DBS may not only cause side eects by electrical current spreading outside of the target region, but also by chronic stimulation of the target itself or by functional disconnection of the stimulated structure (see e.g. [START_REF] Ferraye | Eects of subthalamic nucleus stimulation and levodopa on freezing of gait in Parkinson disease[END_REF]). Accordingly, it is key to reduce the integral stimulation current. Electrical CR stimulation of the STN employs signicantly less current compared to HF DBS [P. A. Tass, Qin, et al. 2012b;Adamchic, Hauptmann, et al. 2014;[START_REF] Wang | Coordinated Reset Deep Brain Stimulation of Subthalamic Nucleus Produces Long-Lasting, Dose-Dependent Motor Improvements in the 1-Methyl-4-phenyl-1,2,3,6-tetrahydropyridine Non-Human Primate Model of Parkinsonism[END_REF]]. However, to further improve the CR approach, in a previous computational study the spacing principle [START_REF] Ebbinghaus | Memory; A Contribution to Experimental Psychology[END_REF] was used to achieve an anti-kindling at subcritical intensities, i.e. particularly weak intensities rendering permanently delivered CR stimulation ineective [O. V. Popovych, Xenakis, and P. A. Tass 2015]. According to the spacing principle [START_REF] Ebbinghaus | Memory; A Contribution to Experimental Psychology[END_REF], learning eects can be improved by repeated stimuli spaced by pauses as opposed to delivering a massed stimulus in a single long stimulation session. Spaced CR stimulation at subcritical intensities might possibly be applied to CR DBS. However, for clinical applications, in particular, for non-invasive applications of CR stimulation, such as acoustic CR stimulation for tinnitus or vibrotactile stimulation for PD, it is crucial to achieve therapeutic eects within a reasonable amount of time. Applications of non-invasive medtech devices typically rely on the patients' compliance and should favorably require short stimulation durations. Accordingly, we here set out to apply the spacing principle to CR stimulation at supercritical intensities, i.e. intensities that enable an anti-kindling for moderate stimulation duration and properly selected stimulation frequencies. The overall goal is to design short-term dosage regimen that improve CR stimulation ecacy, while keeping the integral amount of stimulation as well as the overall duration of the protocols at comparably low levels.

To come up with favorable combinations of stimulation parameters, in our numerical study we used dierent data analysis methods, e.g. macroscopic measures assessing the average amount of the population's synchrony and synaptic connectivity. These measures are appropriate to demonstrate relevant stimulation eects, such as stimulation-induced transitions from pronounced neuronal synchrony to desynchronized states.

To establish a pharmacological treatment for clinical use, in humans typically a 4-phase sequence of clinical trials is performed [START_REF] Friedman | Fundamentals of Clinical Trials[END_REF]. In pre-clinical studies pharmacokinetic, toxicity and ecacy are studied in non-human subjects. In human-studies (phase I ) safety and tolerability of a drug are studied in healthy volunteers. Proof of concept studies (phase IIA) determine whether a drug can have any ecacy, whereas dose-nding studies (phase IIB ) are performed to reveal optimum dose at which a drug has biological activity with minimal side-eects. Eectiveness and the clinical value of a new intervention are studied in a randomized controlled trial (phase III ), compared with state of the art treatment, if available. Finally, post-marketing surveillance trials (phase IV ) are performed to detect rare or long-term adverse eects within a much larger patient population and over longer time periods. There might also be combinations of dierent phases.

In principle, this 4-phase pattern is also valid for medical technology, e.g. neuromodulation technologies. However, if neuromodulation technologies aim at the control of complex dynamics of e.g. neural networks, dierent parameters and dosage regimens may have complex, non-linear and even counterintuitive eects, see e.g. [O. V. Popovych, Xenakis, and P. A. Tass 2015]. The latter work illustrates how computational modelling can be used to generate hypotheses for dose-nding studies.

In general, performing dose-nding studies simply by trial and error may be impossible because of the substantial parameter space to be tested, with trial durations and related costs getting out of hands.

The development of proper dosage strategies and regimens enables favorable compromises between therapeutic ecacy and detrimental factors such as side-eects or treatment duration. This is relevant, e.g. for the development of pharmaceutical [START_REF] Dash | Pharmaceutics: Basic Principles and Application to Pharmacy Practice[END_REF] or radiation therapy [START_REF] Symonds | Walter and Miller's Textbook of Radiotherapy: Radiation physics[END_REF]. DBS is the standard treatment of medically refractory movement disorders [START_REF] Benabid | Long-term suppression of tremor by chronic stimulation of the ventral intermediate thalamic nucleus[END_REF][START_REF] Krack | Five-Year Follow-up of Bilat-eral Stimulation of the Subthalamic Nucleus in Advanced Parkinson's Disease[END_REF][START_REF] Deuschl | A Randomized Trial of Deep-Brain Stimulation for Parkinson's Disease[END_REF]]. The clinical [START_REF] Temperli | How do parkinsonian signs return after discontinuation of subthalamic DBS?[END_REF] and electrophysiological [START_REF] Kühn | Highfrequency stimulation of the subthalamic nucleus suppresses oscillatory beta-activity in patients with Parkinson's disease in parallel with improvement in motor performance[END_REF]] eects of standard HF DBS occur only during stimulation and cease after stimulation oset.

In [Manos, Zeitler, and P. A. Tass 2018a], we studied the inuence of the CR stimulation frequency and the intensity on the outcome of CR stimulation with RVS and SVS [START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF]. CR stimulation consists of sequences of stimuli delivered to each sub-population. For RVS CR stimulation, the CR sequence is randomly varied from one CR stimulation period to another [Peter A. [START_REF] Tass | Long-term anti-kindling eects of desynchronizing brain stimulation: a theoretical study[END_REF]. Conversely, SVS CR stimulation is characterized by repeating a sequence for a number of times before randomly switching to the next sequence. In [Manos, Zeitler, and P. A. Tass 2018a], we demonstrated that the ecacy of singleshot CR stimulation with moderate stimulation duration depends on the stimulation parameters, in particular, on the intensity as well as the relationship between CR stimulation frequency and intrinsic ring rates.

RVS CR stimulation turned out to induce pronounced long-lasting desynchronization, e.g. at weak

intensities and CR stimulation frequencies in a certain range around the neurons' intrinsic ring frequencies. In contrast, SVS CR stimulation enabled even more pronounced anti-kindling, however, at the cost of a signicantly stronger dependence of the stimulation outcome on the CR stimulation frequency. Dosage regimen design is an integral part of pharmacokinetic methodology, aiming at an optimization of drug delivery and eects [START_REF] Williams | Dosage regimen design: pharmacodynamic considerations[END_REF]]. By a similar token, we hypothesize that appropriate dosage regimens might further enhance the ecacy of RVS and SVS CR stimulation.

To probe dierent dosage regimens, we considered dierent stimulation singleshot and multishot CR stimulation protocols. Protocols A and B have identical integral stimulation duration, whereas Protocols C and D may require less stimulation.

Protocol A: Spaced multishot CR stimulation with xed stimulation parameters. Instead of one singleshot CR stimulation we deliver the identical CR shot ve times, where the duration of each single pause equals the duration of each identical singleshot. Intersecting singleshot stimuli by pauses to increase stimulation ecacy, resembles the so-called spacing principle, a learning-related mechanism that is well-established in psychology [START_REF] Ebbinghaus | Memory; A Contribution to Experimental Psychology[END_REF], education [START_REF] Kelley | Making long-term memories in minutes: a spaced learning pattern from memory research in education[END_REF], and neuroscience [START_REF] Naqib | Molecular determinants of the spacing eect[END_REF]. According to the spacing principle, learning eects can be enhanced by delivering a stimulus in a spaced manner, as opposed CR stimulation and, especially, SVS CR stimulation has pronounced periodic characteristics.

Accordingly, the CR stimulation frequency turned out to be a sensitive parameter, in particular, for SVS CR stimulation (see [Manos, Zeitler, and P. A. Tass 2018a]). For this reason, for stage (i)

of Protocol C and D we perform a demand-controlled variation of the CR stimulation frequency to prevent from, e.g. unfavorable resonances or phase locking dynamics. Note these demand-controlled changes of the CR stimulation frequency are mild and hardly change the networks' ring rates.

In [START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF], we tested the performance of the dierent Protocols A-D by selecting unfavorable stimulation parameters, which render CR stimulation ineective according to [Manos, Zeitler, and P. A. Tass 2018a]. By design, Protocols C and D work well for all parameter pairs (K, T s ) related to eective singleshot CR stimulation. In that case, CR stimulation actually ceases due to lack of demand. Note, in all four stimulation protocols we keep the stimulation intensity xed. Only Protocols C and D require feedback of the stimulation outcome.

The Hodgkin-Huxley spiking neuron model

In this study we use the conductance-based Hodgkin-Huxley neuron model [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF] for the description of an ensemble of spiking neurons. The set of equations and parameters read (see [START_REF] Hansel | Phase Dynamics for Weakly Coupled Hodgkin-Huxley Neurons[END_REF][START_REF] Popovych | Depending on the T s value, more (or less) ON-cycles may be administered within a xed time interval[END_REF]):

C dV i dt = I i -g N a m 3 i h i (V i -V N a ) -g K n 4 i (V i -V K ) -g l (V i -V l ) + S i + F i (3.1) dm i dt = α m (V i )(1 -m i ) -β m (V i )m i (3.2) dh i dt = α h (V i )(1 -h i ) -β h (V i )h i (3.3) dn i dt = α n (V i )(1 -n i ) -β n (V i )n i (3.4)
The variable V i with i = 1, . . . , N , describes the membrane potential of neuron i, and:

α m (V ) = (0.1V + 4)/[1 -exp(-0.1V -4)], (3.5
)

β m (V ) = 4 exp[ -V -65 18 ], (3.6
)

α h (V ) = 0.07 exp[ -V -65 20 ], (3.7
)

β h (V ) = 1 [1 + exp(-0.1V -3.5)]
, (3.8) α n (V ) = 0.01V + 0.55 [1 -exp(-0.1V -5.5)] ,

(3.9)

β n (V ) = 0.125 exp[ -V -65 80
].

(3.10)

The total number of neurons is N = 200, while g N a = 12 mS/cm 2 , g K = 36 mS/cm 2 , g l = 0.3 mS/cm 2 are the maximum conductance per unit area for the sodium, potassium and leak currents respectively. The constants V N a = 50 mV, V K = -77 mV and V l = -54.4 mV refer to the sodium, potassium and leak reversal potentials respectively. C is the constant membrane capacitance (C = 1 µF/cm 2 ), I i the constant depolarizing current injected into neuron i, determining the intrinsic ring rate of the uncoupled neurons. For the realization of dierent initial networks, we used random initial conditions drawn from uniform distributions, i.e. I i ∈ [I 0 -σ I , I 0 + σ l ] ( I 0 = 11.0 µA/cm and σ l = 0.45 µA/cm ), h i , m i , n i ∈ [0,1] and V i ∈ [-65, 5] mV. The initial values of the neural synaptic weights c ij are picked from a normal distribution N (µ c = 0.5 µA/cm 2 , σ c = 0.01 µA/cm 2 ). Hence, in this setup the neurons are not identical. The S i term refers to the internal synaptic input of the neurons within the network to neuron i, while F i represents the current induced in neuron i by the external CR stimulation. In addition, in order to model variations of the model parameters, we add a sinusoidal external current input of the form I var = A • sin(2π • f • t) to the right-hand side of Eq. (3.2), where f and A are the frequency and the amplitude of the signal respectively.

Network and neuron coupling description

The N = 200 spiking Hodgkin-Huxley neurons are placed on a ring and the N s = 4 stimulations sites are equidistantly placed in space at the positions of neurons i = 25, 75, 125, 175. The neurons interact via excitatory and inhibitory chemical synapses by means of the postsynaptic potential (PSP) s i which is triggered by a spike of neuron i [START_REF] Gerstner | A neuronal learning rule for sub-millisecond temporal coding[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience: The geometry of excitability and bursting[END_REF]] and modelled using an additional equation:

ds j dt = 0.5 (1 -s j ) 1 + exp [-(V j + 5) /12]
-2s j . (3.12) where V r,j is the reversal potential of the synaptic coupling (20 mV for excitatory and 40 mV for inhibitory coupling), and c ij is the synaptic coupling strength from neuron j to neuron i. There are no neuronal self-connections within the network (c ii = 0 mS/cm 2 ). The variable:

S i = N -1 N ∑ j=1 (V r,j -V i ) c ij |M ij |s j ,
M ij = ( 1 - d 2 ij σ 2 1 ) exp ( - d 2 ij ( 2σ 2 2
) ) (3.13) describes the spatial prole of coupling between neurons i and j and is of a Mexican hat-type [START_REF] Cowan | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF][START_REF] Dominguez | A Spiking Neuron Model of Cortical Correlates of Sensorineural Hearing Loss: Spontaneous Firing, Synchrony, and Tinnitus[END_REF][START_REF] Rocha | Linking the Response Properties of Cells in Auditory Cortex with Network Architecture: Cotuning versus Lateral Inhibition[END_REF] with strong short-range excitatory (M ij > 0) and weak long-range inhibitory interactions (M ij < 0). Here d ij = d|i -j| is the distance between neurons i and j, while:

d = d 0 N -1 (3.14)
determines the distance on the lattice between two neighboring neurons within the ensemble. d 0 is the length of the neuronal chain d 0 = 10. σ 1 = 3.5 , and σ 2 = 2.0. In order to limit boundary eects, we consider that the neurons are distributed in such a way that the distance d ij is taken as:

d • min (|i -j|, N -|i -j|) for i, j > N/2.

Spike-Timing-Dependent Plasticity

The synaptic weights c ij are dynamical variables that depend on the time dierence, ∆t ij = t i -t j , between the onset of the spikes of the post-synaptic neuron i and the pre-synaptic neuron j, denoted by t i and t j , according to [START_REF] Bi | Synaptic modications in cultured hippocampal neurons: Dependence on spike timing, synaptic strength, and postsynaptic cell type[END_REF][START_REF] Popovych | Depending on the T s value, more (or less) ON-cycles may be administered within a xed time interval[END_REF]:

∆c ij =    β 1 e -∆t ij γ 1 τ , ∆t ij ≥ 0 β 2 ∆t ij τ e ∆t ij γ 2 τ , ∆t ij < 0 (3.15)
with parameters β 1 = 1 , β 2 = 16 , γ 1 = 0.12, γ 2 = 0.15, τ = 14 ms and δ = 0.002, According to the value of ∆t ij , the synaptic weight c ij is updated in an event-like manner, i.e. we add or subtract an increment δ ċij for excitatory or inhibitory connections respectively, with learning rate δ > 0 every time a neuron spikes. Furthermore, we restrict the values of c ij on the interval [0,1] mS/cm 2 for both excitatory and inhibitory synapses, ensuring in this way that their strengthening or weakening remains bounded. The maximal inhibitory synaptic weight cmax was set to be 1 in all our stimulations.

The time window of the plasticity is adjusted with respect to the intrinsic ring rate of the neuron population in order to exhibit multistability. There, dierent time-windows (via dierent choices of parameters) were selected for the STDP for two dierent neuron models, i.e. one with bursting neurons (FitzHugh-Rinzel) and one for spiking neurons (Hodgkin-Huxley). In our simulations, the STDP tends to simply stabilize the ongoing ensemble evolution and does not, by itself, (de-)synchronize the network. The parameters were, in general, chosen such that the ratio

∆t ij γ 1,2 τ
is normalized, and the plasticity takes place within a time interval associated with the spiking period of the individual neurons. We analyzed two additional cases for small variation of the plasticity time-window (τ = 12 and τ = 16) and obtained very similar general eects. The selected xed value τ = 14, used throughout the entire study, also allows us to compare our results with previously published studies.

Coordinated Reset Stimulation

Coordinated Reset (CR) stimulation was applied to the neuronal ensemble of N spiking Hodgkin-Huxley neurons. This was done sequentially via N s (= 4 in this study) equidistantly spaced stimulation sites [P. A. Tass 2003a]: one stimulation site was active during T s /N s , while the other stimulation sites were inactive during that period. After that another stimulation site was active during the next T s /N s period. All N s stimulation sites were stimulated exactly once within one stimulation ON-cycle. Therefore, the duration of each ON-cycle is T s (in ms). The spatiotemporal activation pattern of stimulation sites is represented by the indicator functions ρ k (t)(k ∈ [1, . . . , N ]), taking the value 1 when the k th stimulation site is active at t and 0 else.

G stim (t) = t -t k τ stim e -(t-t k )/τ stim , t k ≤ t ≤ t k+1 .
(3.16)

Here τ stim = Ts 6Ns denotes the time-to-peak of G stim , and t k is the onset of the k th activation of the stimulation site. Note that the period (or frequency) through the τ stim parameter of the CR stimulation determines not only the random onset timing of each single signal but also its temporal duration. The spatial spread of the induced excitatory postsynaptic currents in the network is dened by a quadratic spatial decay prole given as a function of the dierence in index of neuron i and the index x k of the neuron at stimulation site k: (3.17) with d the lattice distance between two neighboring neurons as dened in Eq. (3.14) and σ d = 0.8 the spatial decay rate of the stimulation current. Thus the total stimulation current is expressed by the following equations: (3.18) where V r = 20 mV denotes the excitatory reversal potential, V i the membrane potential of neuron i, K the stimulation intensity, and ρ, G, D as described above. For the RVS CR stimulation, sequences are randomly chosen from a set of N s !(= 24) dierent possible sequences during each ON-cycle (an example is shown in Fig. 3.1A for CR stimulation period T s = 10 ms for the rst 90 ms of an activated CR period). Each newly drawn sequence is indicated by a dierent color and lasts exactly one ON cycle. On the other hand, for the SVS-l CR stimulation, one rst randomly picks a sequence and repeats it l times before switching to another

D (i, x k ) = 1 1 + d 2 (i -x k ) 2 /σ 2 d ,
F i = [V r -V i (t)] • K Ns ∑ k=1 D (i, x k ) ρ k (t) G stim (t) ,

Macroscopic measurements and statistical tools

The synaptic weights, being aected by the STDP and the dierent intrinsic periods of the neurons, change dynamically in time. One ecient way to measure the strength of the coupling within the neuronal population at time t is given by the following synaptic weight (averaged over the neuron population):

C av (t) = N -2 ∑ i,j
sgn(M ij )c ij (t), (3.19) where M ij is dened in Eq. (3.13) and sgn is the sign-function. Furthermore, one may additionally measure the degree of the neuronal synchronization within the ensemble, using the order parameter [START_REF] Haken | Synergetics[END_REF][START_REF] Kuramoto | Chemical oscillations, waves, and turbulence[END_REF]:

R(t) = |N -1 ∑ j e iφ j (t) | (3.20)
where φ j (t) = 2π(t-t j,m ) t j,m+1 -t j,m for t j,m ≤ t < t j,m+1 is a linear approximation of the phase of neuron j between its m th and (m + 1) th spikes at spiking times t j,m and t j,m+1 . R (t) is inuenced by the synaptic weights, as the latter are time dependent due to the STDP. The order parameter R measures the extent of phase synchronization in the neuronal ensemble and takes values between 0 (absence of in-phase synchronization) and 1 (perfect in-phase synchronization).

In our numerical calculations, we estimate C av [see Eq. (3.19)] and R av . The latter quantity is averaged over the last 100 • T s . Whenever we plot the order parameter versus time, we determine the moving average < R > over a time window of 400 • T s , because of the presence of strong uctuations. For the statistical description and analysis of the non-Gaussian distributed C av and R av data (n = 11 samples), we use the median as well as the Inter-Quartile Range (IQR).

Results

Initial simulation setup

For each initial network of N = 200 non-identical-neurons and parameter conditions (or simply network), we apply RVS and SVS CR signals (dierent per network). For each network, the initial conditions for each neuron were randomly drawn from random distributions as given in the Hodgkin-Huxley Spiking Neuron Model subsection. We start the simulation with an equilibration phase, which lasts 2 s. Later on, we evolve the network under the inuence of STDP (which will be present until the end of the simulation). We then integrate the network for 60 s with STDP without any external stimulation yet, where a rewiring of the connections takes place, resulting in a strongly synchronized state. Next, we apply CR stimulation for 128 s (resetting the starting time to t = 0 s). During this CR-on period three stimulation ON-cycles (the stimulation is on) alternated with two OFF-cycles (the In order to probe and chart the CR stimulation intensity and frequency parameter space, we restrict the CR stimulation intensity to values in the interval (K ∈ [0.20, ..., 0.50]). This particular choice is based on our previous experience and numerical studies (see e.g. [O.V. [START_REF] Popovych | Depending on the T s value, more (or less) ON-cycles may be administered within a xed time interval[END_REF][START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF]) where it was found that weaker intensities were not able to suciently desynchronize the neuron ensemble while larger intensities did not signicantly improve (or sometimes even worsen) the outcome of RVS and SVS CR stimulation signals. We then set an initial-central value for the CR stimulation period (that denes the initial/starting frequency) which in principle is selected close to the intrinsic ring rate of the strongly synchronized network. In this case, and before applying the CR stimulation, the intrinsic ring rate of the network is ≈71 Hz which corresponds to T s ≈ 14 ms. Hence, we begin with the CR stimulation period T 0 = 16 ms which gives an initial stimulation frequency f 0 = 1/T 0 (in a similar manner just like in [O.V. [START_REF] Popovych | Depending on the T s value, more (or less) ON-cycles may be administered within a xed time interval[END_REF][START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF] and adjusted to a value close to the intrinsic one).

Then we dene such a period interval [T smin , T smax ] in ms (T s : integer) that allows us to create an approximately equidistant grid on the frequency space: f stim ∈ [25%f 0 , ..., 175%f 0 ]. This initial T 0value is also well studied for dierent types of signal patterns aiming to optimize the CR eect with the use of dierent type of CR stimulation sequences (see e.g. [START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF]). Then, we dene the ratio (%) of CR sequence frequency per ON-cycle (f stim ) over the frequency of the reference stimulation frequency (f 0 = 62.5 Hz, T 0 = 16 ms) as r 0 = (f stim /f 0 ) • 100 and we end up in studying the intensity and frequency-ratio (K, r 0 ) parameter space. In Table 1 in [Manos, Zeitler, and P. A. Tass 2018a], we show the conversion between the stimulation frequency-ratio and period.

For comparison reasons, we also give the corresponding ratios r int (%) of CR stimulation frequency per ON-cycle (f stim ) over the frequency of the intrinsic ring rate of the network frequency (f int = 71.4 Hz, T int = 14 ms) without any external stimulation.

Impact of CR Stimulation Duration and Signals on Dierent Initial Networks

Before presenting the core of our ndings, let us rst start by discussing how the RVS CR stimulation duration aects the long-lasting anti-kindling of dierent initial randomly chosen networks. In Fig. 3.2, we show the evolution of the mean synaptic weight C av as a function of time for dierent total CR-on time durations: t = 64 s (Fig. 3.2A), t = 128 s (Fig. 3.2B), and t = 256 s (Fig. 3.2C). 128 s is the standard CR-on period used throughout the paper. The CR stimulation intensity is K = 0.20, and the period T s = 10 ms. A general trend appears in this sequence of panels, i.e. the longer the CR stimulation lasts, less spread of the C av regarding the long-lasting anti-kindling eect is observed after stimulation oset. This is shown in Fig. 3.2D with boxplots. The last box (corresponding to t = 256 s of total CR-on period) has no outliers and shows a more uniform long-lasting eect (as shown in Fig. 3.2C) for all 11 network initializations, not only during the CR-on period but also afterwards during the CR-o period. However, there is no statistically signicant decrease of the median of the C av from t = 64 s to t = 128 s (right-sided Wilcoxon rank sum test, p = 0.0209, 5% signicance level). Moreover, the median value of the C av does not change signicantly between t = 128 s (Fig. 3.2B) and t = 256 s (Fig. 3.2C, both-sided Wilcoxon rank sum test, p = 0.8955). Hence, the intermediate stimulation duration t = 128 s provides fairly good results. Furthermore, for considerably larger stimulation durations the anti-kindling is typically, but not always more pronounced. From a clinical standpoint, it is desirable to achieve reasonably pronounced stimulation eects without excessive stimulation durations. Accordingly, in this computational study we choose t = 128 s as total CR-on time, and t = 256 s as total CR-on/o time.

For the dierent simulations, we use dierent random initial networks and CR signals. For the sake of generality, we do not pick any optimal combination of random initial network and RVS CR times the length of the box (above/below). There is no statistically signicant dierence between the data sets at t = 128 s and t = 256 s (p = 0.8955 two-sided Wilcoxon rank sum test). The total CR-on/o time is twice as long as the CR-on period. (E) An identical RVS CR stimulation signal (the one of network 1) was used for all 11 initial networks for t = 128 s [comparison with (B)]. In all cases, the CR stimulation intensity is K = 0.20 with period T s = 10 ms. stimulation signal that would induce a favorable or biased behavior. This is to assess whether CR eects are robust with respect to dierent initial conditions. Fig. 3.2B shows a typical example where 11 dierent random stimulation signals where applied to 11 dierent initial networks during the CR-on period, with CR stimulation intensity K = 0.20 and stimulation period T s = 10 ms.

The CR-on/o period lasts 128 ms respectively.

During the CR-on period the mean synaptic weights C av evolve in a similar manner for all networks, with little deviations between the dierent curves. They reach approximately the same small value at the end of the CR-on period. The latter corresponds to weak excitatory synaptic connectivity and, in most cases in this paper, to globally well-desynchronized states. However, the post-stimulation dynamics of C av may be quite diverse.

Some networks retain their weak average connectivity while others, like network 2 and 9 (Fig. showing that the dispersion around the median of the C av results is very small in large parts of the parameter plane. In contrast, small IQRs are found only for small R av , in regions with strong desynchronization. Figs 3.3A and 3.3B display two main bands in the (K, r 0 )-parameter space associated with small dispersion: The rst band is aligned along the horizontal axis, for weak stimulation intensities (i.e. K = 0.20 and K = 0.25) and stimulation frequencies greater than 40% of the standard f 0 corresponding to a stimulation period of T 0 = 16 ms. The second band runs along the vertical stimulation intensity K axis, and for relatively high frequencies, i.e. for f stim = 160%f 0 (T s = 10 ms) and f stim = 175%f 0 (T s = 9 ms) which correspond to ≈ 155% and ≈ 140% of the ring rate of the synchronized neurons, respectively. For these (bottom-horizontal and right-hand-side-vertical bands) the dispersion around the median values is quite small for both C av and R av (Figs. 3.3C and 3.3D). In addition, the vertical stripe at the reference frequency value f 0 (100%, T s = 16 ms), studied in [Zeitler and P. A. For the stage (i) control, the variation of the CR stimulation frequency is not adapted to frequency characteristics of the neuronal network. Rather a minor variation of the CR stimulation frequency is performed to make a fresh start with the subsequent single CR shot. These minor changes of the CR stimulation frequency do not lead to changes of the neurons' intrinsic ring rates of more than ±3%. Due to the stage (iii) control, the demand-controlled shutdown of CR stimulation, the maximum integral stimulation duration of Protocol C and D can reach the level of Protocols A and B, but may well fall below. We use the order parameter to assess the amount of synchronization. We here present a few representative results regarding Protocols A and C only. The complete analysis for all Protocols can be found in [START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF].

Protocol A: Spaced multishot CR stimulation with xed stimulation parameters.

For this stimulation protocol all stimulation parameters are kept constant (Fig. The spacing protocol does not cause an improvement of the stimulation outcome in this case, too [START_REF]Chaos Detection and Predictability[END_REF][START_REF] Popovych | What Can Simulations Contribute to Neuroimaging Data Analytics[END_REF]. Case II:

(K, T s ) = (0.20, 28).

Protocol C: Spaced multishot CR stimulation with demand-controlled variation of stimulation period T s and intensity.

We study the stimulation outcome of only ve symmetrically spaced consecutive single CR shots with stimulation period T s and intensity varied according to a three-stage control scheme. To this end, for both RVS CR and SVS CR (not shown here) stimulation we consider two unfavorable parameter pairs of xed CR stimulation period and intensity, respectively. One example refers to cases where CR stimulation induces acute eects, but no long-lasting desynchronizing eects. The other example concerns the case where CR stimulation causes neither acute nor long-lasting desynchronizing eects in a reliable manner. We consider a regular and a random type of demand-controlled variation of the CR stimulation period T s . Note, in both cases the CR stimulation period is not adapted to frequency characteristics of the network. We consider the time courses of the time-averaged order parameter < R > and R av , the order parameter averaged over a window of length 100 • T s at the end of pause.

Demand-controlled regular variation of the CR stimulation period and demand-controlled variation of the intensity: At the end of each pause we calculate the order parameter R av averaged over a window of length 100 • T s . We vary the CR stimulation period and intensity according to the amount of synchrony, based on a three-stage control scheme:

1. Insucient desynchronization: If R av > 0.4, we decrease the CR stimulation period of the subsequent RVS shot by T s (j + 1) = T s (j) -1 ms, where the index j stands for the j-th CR shot. As lower bound we set T s = 9 ms (corresponding to ≈ 156% of the intrinsic ring rate), in order to avoid undesirably high CR stimulation frequencies. In a previous computational study the latter turned out to be unfavorable for desynchronization (see [Manos, Zeitler, and P. A. Tass 2018a]). As soon as T s reaches its lower bound of 9 ms, it is reset to T s (1) + 1 ms.

2. Moderate desynchronization: If 0.2 ≤ R av ≤ 0.4, we preserve the CR stimulation period for the subsequent CR shot: T s (j + 1) = T s (j), where the index j denotes the j-th CR shot. 0.2 ≤ R av ≤ 0.4 is considered to be indicative of a desynchronization eect.

3. Sucient desynchronization: If R av < 0.2, the CR stimulation is suspended for the subsequent shot by setting K = 0 for the next shot and until 0.2 ≤ R av . R av < 0.2 is considered a sucient desynchronization.

Spaced multishot RVS CR stimulation with demand-controlled regular variation of the stimulation period T s and demand-controlled variation of the intensity: In both Cases (I and II) this protocol reliably induces a desynchronization for all networks tested (Figs. In realistic biological systems intrinsic (model) parameters typically vary over time. These variations may be of complex dynamical nature [see e.g. [START_REF] Timmer | Pathological tremors: Deterministic chaos or nonlinear stochastic oscillators[END_REF][START_REF] Yulmetyev | Regular and stochastic behavior of Parkinsonian pathological tremor signals[END_REF]]. To obtain some indication as to whether Protocol C is robust against low-amplitude intrinsic variations of the neuronal ring rates, we added a low-amplitude term I var = A • sin(2π • f • t) to the right-hand side of Eq. (3.1). In the stimulation-free case, I var causes variations of the neurons' ring rates in the order of ±3% and no qualitative changes of the network dynamics. For dierent frequencies f this type of variation does not signicantly aect the long-term desynchronization outcome of Protocol C (f = 0.004 Hz, 4 Hz and 20 Hz, see Figs. 3.8 and 3.9). By the same token, the neuronal ring rates are not signicantly altered by the additional periodic force.

Note, this is not intended to be a comprehensive study of the impact of periodic forcing of arbitrary frequency on the spontaneous or stimulation-induced dynamics of the model network under consideration. Rather, the slow oscillatory forcing is meant to model slow physiological modulatory processes in an illustrative manner. In the extreme case of f = 0.004 Hz the slow oscillatory modula-

Summary and outlook

By systematically varying the CR stimulation frequency and intensity and comparing the stimulation outcome of the two dierent CR protocols, RVS and SVS CR stimulation [Manos, Zeitler, and P. A. Tass 2018a], RVS CR proved to be more robust with respect to variations of the stimulation To the best of our knowledge, in our study in a plastic network the CR stimulation frequency and intensity were systematically varied for the rst time to investigate the impact on acute and long-lasting stimulation outcome. Remarkably, pronounced acute desynchronization (as measured by means of the standard order parameter) does not necessarily lead to long-lasting desynchronization.

On the one hand, this nding might inspire future computational and pre-clinical studies aiming at specically designing stimulation protocols for long-lasting (as opposed to acute) desynchronization.

On the other hand, this nding is signicant for the development of clinical calibration procedures [START_REF]Chaos Detection and Predictability[END_REF][START_REF] Popovych | What Can Simulations Contribute to Neuroimaging Data Analytics[END_REF].

for CR stimulation, see [START_REF] Adamchic | Acute eects and after-eects of acoustic coordinated reset neuromodulation in patients with chronic subjective tinnitus[END_REF].

We cannot expect a stimulation technique to be generically eective, irrespective of the neural network model used. Nevertheless, stepwise adding further physiologically and anatomically relevant features to the neural network models employed may help to generate specic predictions and, ultimately, to further improve stimulation protocols and dosage regimes. In that sense, the nding that RVS CR stimulation at weak to moderate intensities and stimulation frequencies adapted to the neurons' intrinsic ring rates causes a desynchronization in neural network models without STDP [Lysyansky, O. V. Popovych, and P. A Tass 2011] and with STDP as shown in this study, is relevant and, in fact, in accordance with pre-clinical ndings [J. [START_REF] Wang | Coordinated Reset Deep Brain Stimulation of Subthalamic Nucleus Produces Long-Lasting, Dose-Dependent Motor Improvements in the 1-Methyl-4-phenyl-1,2,3,6-tetrahydropyridine Non-Human Primate Model of Parkinsonism[END_REF]P. A. Tass, Qin, et al. 2012a]. Furthermore, the fact that SVS CR stimulation might even be more eective, but requires more careful parameter adaptation may guide future development of calibration techniques as put forward in [START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF].

We demonstrated that over a wide range of stimulation parameters favorable acute eects do not automatically lead to favorable long-lasting, sustained after-eects. This is in agreement with a computational study in the same model, but performed in only a restricted parameter range [START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF], as well as with an electroencephalograms (EEG) experiment performed in tinnitus patients [START_REF] Adamchic | Acute eects and after-eects of acoustic coordinated reset neuromodulation in patients with chronic subjective tinnitus[END_REF]]. To characterize stimulation induced ef-fects, we here used the average synaptic weight and the average amount of neuronal synchrony. These macroscopic quantities enabled us to eectively investigate the impact of variations of stimulation parameters on the stimulation outcome.

In neural networks without STDP tested so far, CR stimulation works at higher intensities as well, see e.g. [START_REF] Lysyansky | Desynchronizing anti-resonance eect of m:n ON-OFF coordinated reset stimulation[END_REF]. In that case, pronounced cluster states are induced, but coherent synchrony is reliably suppressed. This is not the case in the neural network model with STDP studied here. [Manos, Zeitler, and P. A. Tass 2018a]). This is relevant for two reasons:

(i) Since eective CR stimulation does not require phase resets time-locked to the individual stimuli, further computational studies should elucidate whether it makes sense to calibrate CR stimuli for preclinical and clinical applications by selecting stimulus parameters that favorably achieve phase resets.

Corresponding results might be relevant for the design of calibration procedures and, in addition, challenge existing patents that are based on selecting parameters that optimally achieve phase resets of the stimuli delivered to the individual sub-populations. (ii) By the same token, our results do not only challenge current hypotheses on the mechanism of CR stimulation, but also fundamental patents in the eld of invasive as well as non-invasive CR stimulation. Accordingly, future computational studies should focus on the mechanism of action of CR stimulation in networks with STDP in order to actually understand and possibly improve anti-kindling protocols.

Our goal was to accomplish an anti-kindling in a way as robust as possible, complying with clinically motivated constraints. For instance, striving for anti-kindling induced at minimal stimulation We computationally showed that a spacing with rigid ve-shot timing structure, but exible, demand-controlled variation of stimulation frequency and intensity (Protocol C) provides a short-term dosage regimen that signicantly improves the long-term desynchronization outcome of RVS and SVS CR stimulation. At the end of each pause between CR shots, the stimulus after-eect is assessed. If the desynchronization is considered to be insucient, a mild variation of the CR stimulation frequency is performed to possibly provide a better t between network and CR stimulation frequency, without actually adapting the stimulation frequency to frequency characteristics of the network stimulated. If desynchronization is considered to be moderate, the subsequent CR shot is delivered with parameters unchanged. If desynchronization is sucient, CR stimulation is suspended during the subsequent shot.

Intriguingly, in the vast majority of parameters and networks tested, this short-term dosage regimen induces a robust and reliable long-lasting desynchronization. This protocol might be a candidate especially for non-invasive, e.g. acoustic [P. [START_REF] Tass | Counteracting tinnitus by acoustic coordinated reset neuromodulation[END_REF] or vibrotactile [START_REF] Syrkin-Nikolau | Coordinated reset vibrotactile stimulation shows prolonged improvement in Parkinson's disease[END_REF][START_REF] Tass | A model of desynchronizing deep brain stimulation with a demand-controlled coordinated reset of neural subpopulations[END_REF], applications of CR stimulation to increase desynchronization ecacy, while keeping the stimulation duration at moderate levels.

Demand-controlled variation of CR stimulation frequency and intensity (Protocol D) alone (i.e.

without inserting pauses) is not sucient to signicantly improve the outcome of RVS and SVS stimulation [START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF]. Hence, introducing pauses signicantly improves the eect of the demand-controlled variation of CR stimulation frequency and intensity. In principle, stimulation parameters other than the CR stimulation frequency might be varied depending on the stimulation outcome. However, in this study we have chosen to vary the CR stimulation frequency, since the latter turned out to be a sensitive parameter, especially for SVS CR stimulation (see [Manos, Zeitler, and P. A. Tass 2018a]). In fact, the short-term dosage regimen with demand-controlled variation of stimulation parameters (Protocol C) might help to turn SVS CR stimulation in a method that causes a particularly strong anti-kindling in a robust and reliable manner.

Protocol C does not require a direct adaption of the CR stimulation frequency to measured quantities reecting frequency characteristics of the network. We have chosen this design, since it might be an advantage not to rely on specic biomarker-type of information. For instance, in the case of PD a number of relevant studies were devoted to closed-loop DBS (see e.g. [START_REF] Graupe | Adaptively controlling deep brain stimulation in essential tremor patient via surface electromyography[END_REF][START_REF] Rosin | Closed-Loop Deep Brain Stimulation Is Superior in Ameliorating Parkinsonism[END_REF][START_REF] Little | Adaptive deep brain stimulation in advanced Parkinson disease[END_REF][START_REF] Rosa | Adaptive deep brain stimulation in a freely moving Parkinsonian patient[END_REF]). A relevant issue in this context is the availability of a biomarker adequately reecting the individual patient's extent of symptoms [START_REF] Beudel | Adaptive deep brain stimulation in Parkinson's disease[END_REF][START_REF] Kühn | Innovations in deep brain stimulation methodology[END_REF]. In fact, it is not clear whether low or high frequency beta band oscillations might be more appropriate as biomarker-type of feedback signal [START_REF] Beudel | Adaptive deep brain stimulation in Parkinson's disease[END_REF]. For several reasons, beta band oscillations might possibly not be an optimal feedback signal (see e.g. [START_REF] Kühn | Innovations in deep brain stimulation methodology[END_REF]). Enhanced beta band oscillations are not consistently found in all PD patients [START_REF] Beudel | Adaptive deep brain stimulation in Parkinson's disease[END_REF][START_REF] Kühn | Innovations in deep brain stimulation methodology[END_REF]. The clinical score of PD patients might more appropriately be reected by the power ratio of two distinct bands of high frequency oscillations around 250 Hz and 350 Hz [Özkurt et al. 2011]. Appropriate biomarkers might depend on the patient phenotype [START_REF] Quinn | Beta oscillations in freely moving Parkinson's subjects are attenuated during deep brain stimulation[END_REF]]: In tremor dominant (compared to akinetic rigid) PD patients resting state beta power may decrease during tremor epochs [START_REF] Bronte-Stewart | The STN beta-band prole in Parkinson's disease is stationary and shows prolonged attenuation after deep brain stimulation[END_REF][START_REF] Quinn | Beta oscillations in freely moving Parkinson's subjects are attenuated during deep brain stimulation[END_REF]. By a similar token, theta and beta oscillations interact with highfrequency oscillations under physiological [START_REF] Yanagisawa | Regulation of Motor Representation by Phase-Amplitude Coupling in the Sensorimotor Cortex[END_REF] as well as pathological [START_REF] Yang | Beta-Coupled High-Frequency Activity and Beta-Locked Neuronal Spiking in the Subthalamic Nucleus of Parkinson's Disease[END_REF] conditions. Also, quantities assessing the interaction of brain oscillation, e.g. phase amplitude coupling might be used as biomarker to represent the amount of symptoms [START_REF] Beudel | Adaptive deep brain stimulation in Parkinson's disease[END_REF].

Also, activity in the beta band might be relevant for compensatory purposes, as recently shown in a parkinsonian monkey study with sensorimotor rhythm neurofeedback [START_REF] Philippens | Sensorimotor rhythm neurofeedback as adjunct therapy for Parkinson's disease[END_REF].

It might be another potential advantage for clinical applications that the three-stage control of the proposed short-term dosage regimen (Protocol C) could possibly be approximated by scores reecting the patient's state or the amount symptoms. A simple three-stage rating of the patient's state (bad, medium, and good) might replace the feedback signal-based stages (i), (ii), and (iii). Assessments of the patient's state might be performed in a pause after a CR shot. Depending on the rating, the CR stimulation frequency or intensity of the subsequent CR shot may be varied. In particular, for non-invasive application of CR stimulation a non-invasive assessment of the stimulation eect might straightforwardly be realized.

Intrinsic variations of sucient size might naturally mimic variations of the relationship between CR stimulation frequency and intrinsic neuronal ring rates as introduced on purpose in Protocol C.

Accordingly, already the purely spaced stimulation without demand-controlled variability (Protocol A) might display some variability of the relationships between intrinsic ring rates and CR stimulation frequency simply due to the intrinsic variability. However, at least with the frequencies 0.004 Hz, 4

Hz and 20 Hz in the low-amplitude term I var = A • sin(2π • f • t) added to the right-hand side of Eq. (3.1), we were not able to observe any substantial improvement of the desynchronizing outcome of Protocol A (Fig. 3.9 from Supplementary Fig. 2 in [START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF]).

However, more physiological patterns of ring rate modulations might have a more signicant impact on the stimulation outcome of Protocol A. In future studies typical variations of the signals relevant to a particular pre-clinical or clinical application might be taken into account to further improve desynchronizing short-term dosage regimen. The additional periodic forcing considered here was meant to illustrate the stability of the suggested control approach. However, future studies could also provide a detailed analysis of the interplay of one or more periodic inputs and noise, thereby the focus will now be on a feedback-adjusted modulation of synaptic patterns to induce long-lasting therapeutic eects. Currently, clinical proof of concept (phase IIa ) is available for deep brain CR stimulation for the therapy of Parkinson's disease [Adamchic, Hauptmann, et al. 2014] and acoustic CR stimulation for the treatment of chronic subjective tinnitus [P. [START_REF] Tass | Counteracting tinnitus by acoustic coordinated reset neuromodulation[END_REF]]. In addition, promising rst in human (phase I ) data are available for vibrotactile CR stimulation for the treatment of Parkinson's disease showing pronounced and highly signicant sustained therapeutic eects [START_REF] Syrkin-Nikolau | Coordinated reset vibrotactile stimulation shows prolonged improvement in Parkinson's disease[END_REF]. For the clinical development of these treatments it is mandatory to perform dose-nding studies (phase IIb) to reveal optimal stimulation parameters and dosage regimens, for comparison see [START_REF] Friedman | Fundamentals of Clinical Trials[END_REF]. The latter are required to get properly prepared for large ecacy (phase III ) trials [START_REF] Friedman | Fundamentals of Clinical Trials[END_REF].

Since CR stimulation modulates complex neuronal dynamics, dose-nding studies are sophisticated, since stimulation parameters as well as dosage patterns have to be chosen appropriately. Selecting appropriate stimulation parameters and dosage regimens by trial and error may neither be eective nor aordable, since it would require a huge number of patients. In contrast, this work illustrates the important role of computational medicine in generating hypotheses for dose-nding studies.

Specically, we show that spacing (i.e. adding pauses in between stimulation epochs) as well as moderate and unspecic parameter variations adapted in the course of the therapy are not sucient to overcome limitations of CR stimulation. Intriguingly, the combination of both, spacing plus adaptive moderate parameter variation increases the robustness of the stimulation outcome in a signicant manner. This computational prediction can immediately be tested in dose-nding studies and, hence, help to optimize the CR therapy, shorten the development time and reduce related costs.

Perspectives

Besides the spike timing-dependent plasticity (modulating relative synaptic strength connectivity between neurons on a time scale of seconds to minutes) there is also the structural plasticity (deletion pre-existing or generation new synapses in order to homeostatically adapt the ring rate of the neurons on a time scale of days to months), see e.g. [START_REF] Diaz-Pier | Automatic Generation of Connectivity for Large-Scale Neuronal Network Models through Structural Plasticity[END_REF]]. The latter one is not extensively studied while combining both types of plasticity is a rather complicated task as one should design models that can operate at dierent time scales. An interesting future task could be to investigate the role of structural plasticity in explaining long-term deep brain stimulation eects. In this framework, we aim to set up a rather detailed network and incorporate physiologically relevant mechanisms, such as dynamical evolution of synaptic connections between neurons, in order to optimize the stimulation This type of investigation is not able to take place at a clinical setting and hence one needs to employ theoretical and computational studies in order to ultimately improve the quality of life of patients and minimize side eects of the deep brain stimulation after surgery.

Nonlinear dynamical systems have been recently employed in the study of the human brain also at larger spatio-temporal scales. The human brain is a complex neural network able to self-organise into dierent emergent states crucial for its healthy and pathological functioning. At a meso-and macro-scale level, a brain connectivity network consists of nodes (grey matter regions) and edges.

Edges can represent white matter tracts in structural networks connectome, the so-called structural connectivity, or correlations between for example, two blood-oxygen-level-dependent (BOLD) or EEG time series in functional networks, the so-called functional connectivity (FC), see e.g. [START_REF] Sporns | The Human Connectome: A Structural Description of the Human Brain[END_REF]. Interdisciplinary approaches using concepts from nonlinear dynamics, physics, bi-ology and medicine, have allowed us to understand in more depth how the human brain functions. A dynamical-model approach, can provide links between attractors, bifurcations, synchronisation patterns and empirical neuroimaging data, such as EEG, BOLD functional magnetic resonance imaging (fMRI), etc. By choosing adequate model parameters (e.g. via a parameter sweep exploration), it is feasible to build customised virtual brain activity for individual subjects. These parameters can serve as dynamical biomarkers and predictors of dierent brain states (healthy vs diseased) and behavioural modes [O. V. [START_REF] Popovych | What can computational models contribute to neuroimaging data analytics?[END_REF]]. Following such concepts, the virtual epileptic patient has been recently proposed [START_REF] Jirsa | The Virtual Epileptic Patient: Individualized whole-brain models of epilepsy spread[END_REF], where such medical-treatment approaches using personalised mathematical models for epileptic patients has been illustrated. By better understanding the pathological activity as opposed to health, one aims to optimize appropriately the treatment.

Furthermore, such brain networks nodes and edges are not uniquely dened. Typically, the voxel parcellation is based on certain criteria from the analysis of the neuroimaging data. These criteria may be based for example on structural neural information of the areas or their functionality or coactivation during a task, see e.g. [START_REF] Stanley | Dening nodes in complex brain networks[END_REF][START_REF] Thirion | Which fMRI clustering gives good brain parcellations?[END_REF][START_REF] Eickho | Imaging-based parcellations of the human brain[END_REF]. Ultimately, one is able to reduce the dimensionality of the brain data by merging hundred thousands of voxels from the high-resolution neuroimaging data into a few brain regions. So far, there is no unied brain parcellation. One important benet for such a result would be the interpretability and comparability of the results for dierent subjects and studies.

There are many ways to design such a parcellation map (brain atlas). Recently, we investigate two paradigmatically distinct cases [Popovych O. V., Kyesam J., Manos T., Diaz-Pier S., Hostaedter F., Schreiber J., Yeo, B. T. T. and Eickho S. B., Impact of brain parcellation and empirical data on modeling of the resting-state brain dynamics, (submitted) 2020]: (i) a parcellation of the cortex regions according to its folding properties, e.g., into gyral-based parcels encircled by tracing from the depth of one sulcus to another [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]] and (ii) a parcellation approach based on the brain function, where the patterns of the resting-state FC can be used to group the voxels (or vertices) into parcels of similar connectivity [START_REF] Schaefer | Local-Global Parcellation of the Human Cerebral Cortex from Intrinsic Functional Connectivity MRI[END_REF]]. Thus, we can model and investigate the eects of brain parcellation in great detail with a dynamical system, e.g. coupled phase oscillators suggested for modeling cortical oscillations and resting-state BOLD dynamics [START_REF] Breakspear | Generative models of cortical oscillations: neurobiological implications of the Kuramoto model[END_REF][START_REF] Cabral | Role of local network oscillations in restingstate functional connectivity[END_REF][START_REF] Ponce-Alvarez | Resting-State Temporal Synchronization Networks Emerge from Connectivity Topology and Heterogeneity[END_REF][START_REF] Fukushima | Comparison of uctuations in global network topology of modeled and empirical brain functional connectivity[END_REF].

In this approach, each node is modeled by a phase oscillator (e.g. [START_REF] Kuramoto | Chemical oscillations, waves, and turbulence[END_REF]), the total number of brain node is assigned by the given brain atlas, the relative connectivity strength and delay signal propagation can be estimated by the empirical neuroimaging data (for each subject) as well as estimates about the nodes' frequency. All the above can be used as inputs into the dynamical system and allow for a relatively small free set of parameters (physiologically relevant, such as global coupling, delay scaling etc.) to be explored and tuned in a way that can describe and ultimately explain adequately the dynamics of certain resting-states or brain activity during tasks as well as distinguish between healthy and pathologically states.

  with a, b and c being the semi-axes of the ellipsoidal bar. The corresponding bar potential is:

Figure 1 . 1 :

 11 Figure 1.1: Snapshots of the N -body simulation (upper panels) at four dierent times, displaying stellar density, on the same range, projected on the xy plane. Each frame is 40 by 40 kpc. To illustrate bar lengths and shapes, we overlay ellipses (which are not isophotal ts). The lower panels display the result of evolving an ensemble of initial conditions in the presence of the constructed TD analytical potential.

)Figure 1 . 2 :

 12 Figure 1.2: Time evolution of the halo, disc and bar parameters, measured from the N -body simulation (points) and supplied to the analytical model (tted polynomials). First and second panel: parameters of the Dehnen halo prole. Third: parameters of the MiyamotoNagai disc. Fourth: bar mass and disc mass. Fifth: semi-major axes of the bar. Sixth: bar pattern speed.

Figure 1 . 3 :

 13 Figure 1.3: (Colour online) Parameters of the Dehnen halo prole and of the MiyamotoNagai disc are obtained by tting the circular velocity curves at each time. Here we display four dierent times. The fourth column exhibits total, halo and disc circular velocity curves. The rst and second columns show the disc density proles (radial and vertical, respectively). The third column has the halo density prole. Points come from measurements of the N -body simulation, while lines are tted proles.

Figure 1 . 4 :

 14 Figure 1.4: Bar strength, measured from the simulation (line) and from the analytical model (points).

A 2

 2 is the maximum relative contribution of the m = 2 Fourier component of the mass distribution in the disc.

Figure 1 . 5 :

 15 Figure 1.5: The Poincaré Surface of Section dened by x = 0, p x ≥ 0 with H = -0.19, for three typical orbits (two regular and one chaotic) being integrated for 10 Gyr. The set of parameters for the bar, disc and halo components are chosen from the ts with the 3-d.o.f. TD Hamiltonian at t = 7.0 Gyr of the N -body simulation. In the insets we depict their projection on the (x, y)-plane together with the GALI 2 and MLE σ 1 evolution in time (see Table1in[Manos and Machado 2014]

The 3 D

 3 .O.F. Time-Dependent Case In Fig. 1.8 we show the evolution of an orbit from the ensemble of the N -body simulation with initial condition:

Figure 1 . 6 :

 16 Figure 1.6: The (in)stability map for the 2-d.o.f. frozen potential case using a grid of ≈ 100 000 initial conditions on the PSS and integrating them for 10 Gyr. The colour-bar represents the nal GALI 2 values of each initial condition in the end of the iteration. The rst row refers to a set of potential parameter given at t = 1.4 Gyr for Hamiltonian values H = -0.31, -0.28, -0.25, -0.22, the second row at t = 4.2 Gyr for H = -0.28, -0.25, -0.22, -0.19, the third row at t = 7.0 Gyr for H = -0.28, -0.25, -0.22, -0.19 and the fourth row at t = 11.2 Gyr for H = -0.28, -0.25, -0.22, -0.19.The yellow (light-grey in b/w) colour corresponds to regular orbits where the GALI 2 oscillates around to relatively large positive values, the black color represents the chaotic orbits where GALI 2 tends exponentially small values, while the intermediate colors in the colour-bars between the two represent `weakly chaotic or sticky' orbits. The exact set of parameters used at each t = 1.4, 4.2, 7.0, 11.2 Gyr

Figure 1 .Figure 1 . 8 :

 118 Figure 1.8: The 3-d.o.f. orbit B1 dierent projection on the (x, y)-plane in dierent time windows is depicted in the rst (four-panel) row (from top block-part to the bottom respectively) and the colour bar corresponds to the time (in Gyr). Its GALI 3 and MLE σ 1 evolution in time is shown just below them. Note that the orbit starts as a disc-like and gradually its shape turns to barred, displaying the bar's growth through the parameters of the Ferrers' potential.

Figure 1 . 9 :Figure 1 . 10 :

 19110 Figure1.9: Same as in Fig.1.8 but for the 3-d.o.f. disc-like orbit D1 . Note, that here the disclike pattern slightly varies from cases to case. The dierent degree of chaoticity can be accurately captured by the frequency and fast decay to zero of the GALI 3 [Fig.1.9(e)], indicating that the orbits is relatively `strongly chaotic' . This information can not be revealed in such a way by the MLE σ 1 shown in Figs.1.9(f ).

Figure 1 . 11 :

 111 Figure 1.11: Upper rows: face-on and edge-on views of the ensemble of orbits at the end of each time window. Lower rows: face-on and edge-on views of projected GALI 3 indices for the disc ensemble.Chaotic orbits are those with GALI 3 ≤ 10 -8 . The face-on views also display the ellipses used to dene the regions referred to as: bar, ring, gap and outer disc. Each frame is 20 kpc wide, and the particles are displayed in the reference frame that rotates with the bar.

Figure 1 .

 1 Figure 1.12: Fraction of chaotic motion as a function of time, measured within dierent regions of the disc. These regions are schematically indicated in the right-hand panels, and more clearly detailed in Fig. 1.11.

  regular and chaotic motion. By selecting those orbits which are permanently regular (64.3 per cent) and those which are permanently chaotic (only 0.9 per cent). The remainder (34.8 per cent) change their nature at least once during the evolution. Let us consider rst those orbits that do not undergo any change of regime throughout the evolution (upper rows of Fig. 1.13).

Figure 1 .

 1 Figure 1.13: Evolution of regular/chaotic regimes separated in dierent congurations. Upper rows:

Figure 1 . 14 :

 114 Figure 1.14: The chaoticity of the planar orbits on the equatorial plane of model based on the 2nd snapshot (see Fig. 1.1) is given color coded in (E J , y 0 ) diagrams. The colour of each orbit (each point in the gures) corresponds to the value of the log 10 (GALI 2 ) quantity calculated for it and is taken from the color-bar on the right hand side of the gures. In (a) we calculate log 10 (GALI 2 ) for t 1 = 1 Gyr, while in (b) for t 2 = 10 Gyr. In both gures the zero velocity curve is indicated with ZVC. The continuous black line in the region where motion is allowed is the characteristic of the main family. Capital letters (A,B,...,F) and arrows pointing to the points C and D are used for facilitating the description of the evolution of the curve (see[Chaves-Velasquez et al. 2017] for more details). We observe that in general the orbits with the smaller GALI 2 index in (a), which reach values log 10 (GALI 2 ) -5, become strongly chaotic in (b). However, in (b) appear also additional features indicating chaotic behavior, that are absent in (a). Such features are the dark blue tails above the characteristic for -0.27 E J -0.17. The six heavy dots at E J = -0.2 indicate the initial conditions of the orbits we use to demonstrate the relation between GALI 2 and their morphology (see [Chaves-Velasquez et al. 2017] for more details). Arrows point to the 1st and 6th of them. The ve heavy, yellow dots at E J =-0.206, -0.195, -0.18798, -0.17 and -0.162 indicate the initial conditions of the boxy orbits we presented in [Chaves-Velasquez et al. 2017].

[

  Chaves-Velasquez et al. 2017] indicates that in cases of slow rotating bars as in the TD models, the 3D boxy structure may constitute a major part of the bar. The presence of the X feature in the face-on views of the orbits, as well as the presence of a ring surrounding the bar, raises the question whether a dynamical mechanism as the one proposed by[START_REF] Tsigaridi | Morphologies introduced by bistability in barred-spiral galactic potentials[END_REF], acts in galaxies like IC 5240 presented in [R. J.[START_REF] Buta | The de Vaucouleurs Atlas of Galaxies[END_REF].PerspectivesSuch an analytical TD model has the advantage of oering a fully time-dependent and astrophysically realistic galaxy model, as indicated by the fact that it was successful in recovering several dynamical and morphological features of a barred galaxy. Similar adequate TD models, exhibiting similar dynamical evolution like N -body simulations, are potentially suited to a broad class of applications in galactic dynamics (e.g. double bars, central mass concentrations and even bar dissolution). Once the output of an N -body simulation has been modeled into a time-dependent analytical potential, a variety of analyses could then be undertaken, particularly regarding orbital studies. Work that generally relied on highly simplied (and usually frozen) analytical potentials could take advantage of more astrophysically realistic galaxy models. This bridging would aord an approximation to the richness of detail of an N -body simulation, at a lower computational cost and with the versatility of a simple analytical formulation.

  [Marostica D. A., Machado R. E. G., Manos T. and Athanassoula E., From stellar to halo bars: quantifying the dark matter response, (submitted), 2020.], we aim to characterise the structure and kinematics of the halo bar, with the goal of correlating them with the properties of the stellar bar. Finally, and regarding families of orbits that can contribute to boxy morphologies, one can extend the work in [Chaves-Velasquez et al. 2017] (in frozen time snapshots) and study their robustness in fully TD analytical models [Manos T., Skokos Ch., Patsis P. A., Orbits supporting bars in 3D, rotating, non-autonomous Hamiltonian systems (in preparation), 2020].

Figure 2 . 1 :

 21 Figure 2.1: The classical diusion constant D eff = (∆P ) 2 /n for the SM as a function of K (δK = 0.05) for three discrete times n, i.e., the number of the iterations of the SM, n = 1000 (lower red dashed line), n = 5000 (intermediate blue solid) and n = 10000 (upper black dot-dashed). In the background we have plotted the classical diusion constant D 1 (gray dotted line) [Eq. (2.5)]. The presence of accelerator modes at certain intervals of K (and the sticky objects around) generate anomalous diusion which is rendered by peaks. Here we used ≈100000 (314 × 314) initial conditions uniformly distributed in a grid on the entire phase space [0, 2π] × [0, 2π].

Figure 2 . 2 :

 22 Figure 2.2: The variance of the momentum P in the SM [Eq.(2.3)] with K = 6.5 (red crosses) where small islands and accelerator mode of period 1 are present and K = 10.0 (blue stars) where the phase space is fully chaotic for the same initial conditions as in Fig.2.1 as a function of the discrete time n (number of iterations), in log-log representation. The two slopes associated with dierent types of diusion are µ(K = 6.5) = 1.61252 (dotted), µ(K = 10.0) = 0.991334 (solid) with standard deviation errors ±0.01271 (0.7881%) and ±0.0009537 (0.0962%) respectively.

Figure 2 . 3 :

 23 Figure 2.3: The diusion exponent µ as a function of K after n = 5000 iterations and for ≈100000 (314 × 314) initial conditions on the plane (θ, P ) = (0, 2π). The intervals on the black horizontal line µ = 0.9 indicate the intervals of stable accelerator modes of period 1 [Eq.(2.6)]. All intervals of K with exponent µ ≈ 1 are associated with normal diusion processes. The large peaks (appearing mainly for K > 2π marked with full black circles) reect the anomalous diusion accelerator modes (mainly of period 1). The smaller peaks for K < 2π (more clearly presented in the inset panel) originate by accelerator modes of higher period together with those for 2π < K < 4π marked with empty circle and a few typical examples close to those peaks [marked with the symbol (×)], for which the diusion is normal, are studied thoroughly in[Manos and Robnik 2014]. The blue dotted line corresponds to the power law which describes the decay of the exponent µ of the main peaks' amplitude due to accelerator modes of period 1.

  whilst for anomalous diusion a strong departure from the Gaussian distribution is observed, being well tted by a stable Lévy distribution, characterized by the parameter α ∈ [0, 2]. The details are given in reference[Manos and Robnik 2014]. For each one of the K-values of the nonlinearity kick parameter of Fig.2.3, we have performed a thorough study by calculating and comparing the following quantities:(a) The index of stability α-parameter of the Lévy stable distribution. (b) The diusion exponent µ as described in Eq.(2.4).

  .21)], restricted to the symmetry class of the odd eigenfunctions. It is of course just one of the possible discrete approximations to the continuous innite dimensional model.

  2.4 and 2.5 (from [Manos and Robnik 2013]), we show the examples of strongly exponentially localized eigenstates by plotting the natural logarithm of the probabilities w n = |u n | 2 versus the momentum quantum number n, for two dierent matrix dimensions N . By calculating the localization length L from the slopes σ = 2/L of these eigenfunctions using Eq.(2.15)

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: (a) A sample of strong localized eigenstates for K = 7, r = 222, k ≈ 2.00 and N = 398 (b) Same for K = 7, r = 444, k ≈ 2.00 and N = 796.

Figure 2 . 6 :

 26 Figure 2.6: The parameter β loc vs. Λ for 161 × 398 elements for various values of K and for a wide range of k values, where the scaling law [Eq. (2.33)] is shown with the black line.

Figure 2 . 7 :

 27 Figure 2.7: The probability distribution function of the momentum P after n = 5000 iterations (black lled circles) and the ts with the stable α-Lévy distribution for a sample of K-values as-

Fig. 2 .Figure 2 . 8 :

 228 Figure 2.8: (a) The GALI 2 for K = 3.1 with 50 × 50 initial conditions on a grid 500 × 500 on the entire phase space (θ, P ) ∈ [0, 2π] × [0, 2π]. (b) the diusion exponent µ for the same kick parameter value for 50 × 50 initial conditions on a 500 × 500 cell grid of the entire phase space calculated after n = 5000 iterations (see text for more details). (c) Two examples (marked with an empty square and empty circle in panel b) following dierent diusion processes: a trajectory transported ballistically (with P < 0) by the eect of an accelerator mode of period 4 (•) and one oscillating between islands of stability ( ) of period 4 too.

Figure 2 . 9 :

 29 Figure 2.9: Intermediate statistics (panel (a)) for distribution P (S) (histogram -black line) of the model tted with distribution P BR (blue-solid line), P IZ (red-dashed line pointed by the arrow in the inset gure) and P new IZ (black-dotted line) for M × N = 641 × 398, K = 7 and k = 11 (see text for discussion). In panels (b),(c),(d) we show the dierence of the numerical data and the best tting Brody (blue-thick line), old Izrailev (red-medium line) and Izrailev new/improved (black-thin line;

Figure 2 .Figure 2 . 11 :

 2211 Figure 2.10: χ 2 for P BR and P IZ for N = 398 (a), N = 796 (b) and N = 4000 (c) for various values of the best ts for β BR (black line) and β IZ (gray line) using 150 subintervals in S.

Figure 2 .Figure 2 . 14 :

 2214 Figure 2.13: (a) L versus 2/ σ for matrices of dimension N = 1000 (crosses and solid t line) and for matrices of dimension N = 3000 (stars and dashed t line), for 7 nearby values of k, namely k = k 0 ± jδk, where j = 0, 1, 2, 3 and δk = 0.00125, for k 0= 3, 4, 5, ..., 19. (b) We plot the mean value of 2/(N σ ) versus β loc for k 0 = 3, 4, 5, . . . , 19 and 7 matrices of dimension N = 3000 with k = k 0 ± jδk, where j = 0, 1, 2, 3 and the step size δk = 0.00125.

Figure 2 .Figure 2 .

 22 Figure 2.15: Histograms of l H in (a) and 2/l H in (b) for the system k = 10 described by the matrices of dimension N = 3000. In both cases we show the Gaussian best t.

Figure 2 .

 2 Figure 2.17: The standard deviation of the positive nite time Lyapunov exponents for the SM (stars)and for the product of random transfer matrices with the box distribution of W (empty boxes), as a function of time in log -log presentation, and their best ts. The slope is exactly -1/2.

  Manos and Robnik 2014; Manos and Robnik 2015] on the dynamical localization in the N -dimensional Izrailev model. The analysis of the classical system (SM) and of the quantum kicked rotator has been performed for many dierent values of the classical kick parameter K on the interval [K cr , 70] ≈ [0.97, 70], of the quantum parameter k and matrix dimensions N in the interval [400, 3000]. The aspects of classical generally anomalous diusion have been studied and the important relevance of the accelerator modes were elucidated, and the semiclassical approximation for the average localization length L has been derived. The entropy localization measure l H [Eq. 2.29] has been calculated, and the corresponding parameter β loc = l H /N was dened. The scaling law of β loc versus Λ = L/N was

  in time-periodic systems, of chaotic eigenstates, is still open for the future.Furthermore, in the classical kicked rotor model, one can extend these ndings on diusion properties and potential scaling laws for 1D and 2D lattices and investigate how the fraction of the chaotic, regular motion along with the presence of accelerator modes may impact systems with more complex congurations, e.g. when varying the system size, the kick parameter per rotor and coupling type and strength within the lattice[Moges H., Manos T. and Skokos Ch., Anomalous diusion and chaotic motion in standard map lattices (in preparation), 2020].

  ]. The clinical need for stimulation techniques that cause desynchronization irrespective of the network's initial state [P. A. Tass 2001], thereby being reasonably robust against variations of system parameters and, hence, not requiring time-consuming calibration, motivated the design of the so-called Coordinated Reset (CR) stimulation [P. A. Tass 2003b; P. A. Tass 2003a]. CR stimuli aim at disrupting in-phase synchronized neuronal populations by delivering phase resetting stimuli typically equidistantly in time, separated by time dierences T s /N s , where T s is the duration of a stimulation cycle, and N s is the number of active stimulation sites [P. A. Tass 2003b; P. A. Tass 2003a]. The spatiotemporal sequence by which all stimulation sites are activated exactly once in a CR stimulation cycle is called the stimulation site sequence, or briey sequence. Taking into account STDP in oscillatory neural networks qualitatively changed the scope of the desynchronization approach: Computationally it was shown that CR stimulation reduces the rate of coincident ring and, mediated by STDP, also decreases the average synaptic weight, ultimately preventing the network from generating abnormally increased synchrony. This anti-kindling, i.e., unlearning of abnormally strong synaptic connectivity and of excessive neuronal synchrony, causes long-lasting sustained eects that persist cessation of stimulation. As shown computationally, anti-kindling can robustly be achieved in networks with plastic excitatory and inhibitory synapses, no matter whether CR stimulation is administered directly to the soma or through synapses [O.V. Popovych and P. Tass 2012; P. A. Tass and O. Popovych 2012]. In line with these computational ndings, long-lasting CR-induced desynchronization and/or therapeutic eects were accomplished with invasive as well as non-invasive stimulation modalities.Long-lasting desynchronization was induced by electrical CR stimulation in rat hippocampal slices rendered epileptic by magnesium withdrawal [P.[START_REF] Tass | Long-lasting desynchronization in rat hippocampal slice induced by coordinated reset stimulation[END_REF]]. Electrical CR deep brain stimulation (DBS) caused long-lasting therapeutic after-eects in parkinsonian non-human primates. Bilateral therapeutic after-eects for at least 30 days were caused by unilateral CR stimulation delivered to the subthalamic nucleus (STN) of parkinsonian MPTP monkeys for only 2 h per day during 5 consecutive days. In contrast, standard permanent high-frequency deep brain stimulation did not induce any sustained after-eects (see e.g. [P.[START_REF] Tass | Counteracting tinnitus by acoustic coordinated reset neuromodulation[END_REF]). In patients with Parkinson's disease electrical CR-DBS delivered to the STN caused a signicant and cumulative reduction of abnormal beta band oscillations along with a signicant improvement of motor function[Adamchic, Hauptmann, et al. 2014]. Non-invasive, acoustic CR stimulation was developed for the treatment of patients suering from chronic subjective tinnitus [P. A.[START_REF] Tass | Unlearning tinnitus-related cerebral synchrony with acoustic coordinated reset stimulation: theoretical concept and modelling[END_REF] P. A. Tass, Qin, et al. 2012a]. In a proof of concept-study acoustic CR stimulation caused a statistically and clinically signicant sustained reduction of tinnitus symptoms [P. A.Tass, Qin, et al. 2012a;[START_REF] Adamchic | Psychometric Evaluation of Visual Analog Scale for the Assessment of Chronic Tinnitus[END_REF]] together with a concomitant decrease of abnormal neuronal synchrony[START_REF] Adamchic | Reversing pathologically increased EEG power by acoustic coordinated reset neuromodulation[END_REF], abnormal eective connectivity [A. N.

  proof of concept studies for invasive and non-invasive CR stimulation were driven by computationally derived hypotheses and predictions. Theoretically predicted phenomena and mechanisms, such as long-lasting stimulation eects [Peter A. Tass and Majtanik 2006; O.V. Popovych and P. Tass 2012; Hauptmann and P. A. Tass 2007], cumulative stimulation eects [Hauptmann and P. A. Tass 2009], and improvement by weak stimulus intensity [Lysyansky, O. V. Popovych, and P. A Tass 2011], were veried based on dedicated theory-driven study protocols for pre-clinical and clinical proof of concepts [P. Tass, A.

X Hypothesis # 1 :

 1 Due to the inherently periodic structure of CR stimulation the relation between CR stimulation frequency and the spontaneous neuronal ring rates (prior to stimulation) matters. Periodic delivery of CR stimuli with xed sequence basically constitutes a time-shifted entrainment of the dierent neuronal subpopulations [P. A. Tass 2003b; P. A. Tass 2003a].

  model neurons the impact of CR stimulation intensity and frequency on the desynchronizing outcome of CR was studied in detail [Lysyansky, O. V. Popovych, and P. A Tass 2011].

  to administering a massed stimulus in a single long stimulation session. Computationally, it was shown that subcritical CR stimulation at subcritical (ineective) intensities may become eective if intersected by rather long pauses and delivered suciently often, e.g. eight times [O. V.[START_REF] Popovych | The Spacing Principle for Unlearning Abnormal Neuronal Synchrony[END_REF]. However, shorter pauses were not sucient. As yet, spaced CR stimulation at supercritical intensities was not studied. Here, we focus on comparably short stimulation protocols. Accordingly, we use CR stimulation of sucient intensity and deliver ve single CR shots intersected by pauses. We consider a symmetric dosage regimen, with identical duration of single shots and pauses.Protocol B: Long singleshot CR stimulation with xed stimulation parameters. To assess the impact of the spacing principle, as a control condition we simply stimulate ve times longer, without any pause and with stimulation parameters kept constant. Protocol B is shorter, but employs the same integral stimulation duration as Protocol A.Protocol C: Spaced multishot CR stimulation with demand-controlled variation of the CR stimulation frequency and intensity. As in Protocol A, we deliver spaced CR stimulation comprising ve identical CR shots, intersected by pauses, where all shots and pauses are of equal duration. However, at the end of each CR shot we monitor the stimulation eect and perform a three-stage control: (i)If no pronounced desynchronization is achieved, the CR stimulation frequency of the subsequent CR shot is mildly varied by no more than ±3%. (ii) If an intermediate desynchronization is observed, the CR stimulation frequency remains unchanged and CR stimulation is continued during the subsequent shot. (iii) If a pronounced desynchronization is achieved, no CR stimulation is delivered during the subsequent shot. Note, for stage (i) we do not adapt the CR stimulation frequency to a measured quantity. We consider two dierent variation types employed for stage (i): with regular and with random variation of the CR stimulation frequency. Regular variation means to increase or decrease the CR stimulation frequency in little unit steps. In contrast, random variation stands for randomly picking the CR stimulation frequency from a restricted interval.Protocol D: Long singleshot CR stimulation with demand-controlled variation of the stimulation frequency. To assess the specic pausing-related impact of the evolutionary spacing principle, as a direct control condition we perform Protocol C without pauses. To this end, we string ve CR shots together, without pauses, and evaluate the stimulation eect at the end of each CR shot. If no pronounced desynchronization is achieved, the CR stimulation frequency is slightly varied by no more than ±3% for the subsequent CR shot. During each single CR shot stimulation parameters are kept constant. Only from one CR shot to the next the CR stimulation frequency can be varied. Overall, Protocol D is shorter than Protocol C, but uses the same integral stimulation duration as in Protocols A-C.

( 3 . 11 )

 311 Initially we draw s i ∈ [0 , 1 ] (randomly from a uniform distribution). The coupling term S i from Eq. (3.1) (see [O.V. Popovych and P. Tass 2012]) contains a weighted ensemble average of all postsynaptic currents received by neuron i from the other neurons in the network:

  The stimulation signals induced single brief excitatory post-synaptic currents. The evoked timedependent normalized conductances of the postsynaptic membranes are represented by α-functions as follows [O.V. Popovych and P. Tass 2012]:

Figure 3 . 1 :

 31 Figure 3.1: Time evolution of CR stimulation signals. (A) RVS CR stimulation signal with period T s = 10 ms for the rst 90 ms of an activated CR period. The vertical lines indicate the successive ON-and OFF cycles and the temporal distance between two successive vertical lines correspond to the period T s of each cycle (every stimulation site is activated exactly once during the ON cycles). A change of color indicates a change of sequence. (B) SVS-4 CR stimulation signal with the same period but here the total time spans up to two completed ON-and OFF cycles (∼ 125 ms) while the color changes as a new sequence is drawn.

  stimulation is o ) as in the example stimulation signals shown in Fig 2. Each ON-and OFF-cycle lasts Ts. After 128 s the CR stimulation ceases permanently and we continue the evolution of the CR-o period for extra 128 s.

Figure 3 . 2 :

 32 Figure 3.2: Impact of the total CR-on time on the mean synaptic weight C av for dierent initial random networks and RVS CR. (A) Time evolution of the C av for dierent total CR-on time durations, t = 64 s, (B) t = 128 s (this is the standard CR-on period used throughout the paper) and (C) t = 256 s. In all these cases, 11 dierent initial networks were stimulated with dierent RVS CR stimulation signals during the CR-on period. The thick red horizontal lines indicate the CR-on/o stimulation periods (the end is marked with a vertical gray line) while the horizontal gray dashed lines are visual cues for mutual comparison. (D) Boxplots of the mean synaptic weights presented in (A)-(C), showing the median values (black lines within the boxes). The box frames depict the middle 50%, the upper and lower whiskers the rst and last 25% respectively while the outliers (black dots) are set as 1.5

Fig. 3 .

 3 Fig.3.2B: The outcome at the end of the CR-on period is fairly uniform, while the post-stimulation dynamics of C av is diverse. Replacing one random external stimulation signal by another one may improve the long-term outcome in some cases (e.g. network 8 green dotted line), but worsen the outcome in others (e.g. network 3 blue solid line). These plots indicate that both the random initialization of the network and the dierent stimulation signals during the CR-on period impact on the nal outcome at the end of the CR-o period in a complex manner.

Figure 3 . 3 :

 33 Figure 3.3: Global overview of the mean synaptic weight and synchronization at the end of the CRo period using RVS CR stimulation. (A) Median of the mean synaptic weight C av , (B) median of the order parameter R av (11 dierent random initial network congurations and 11 dierent RVS CR random signals). Long-lasting anti-kindling is achieved in all dark regions as indicated by the corresponding color-bars. Panels (C) and (D) show the dispersion around these median values by plotting their IQR respectively. All IQR values being close to zero indicate that the middle 50% of the distribution are very close to the median value.

Figure 3 . 4 :

 34 Figure 3.4: Fine-grained T s -period grid analysis for RVS CR stimulation at intensity K = 0.20. (A) Boxplots of C av (mean synaptic weight) and (B) R av (order parameter) for xed and weak stimulation intensity K = 0.20 for a ner sample on the T s integer value interval at the end of the CR-o period. The red and green dots indicate the reference stimulation period T 0 = 16 ms and intrinsic ring rate period T int = 14 ms, respectively.

  Case I: RVS CR stimulation induces a desynchronization during the CR shots (Fig.3.6A), but no reliable, long-lasting desynchronization in the subsequent pauses (Fig. 3.6B). The spacing protocol with ve identical RVS CR shots does not signicantly improve the desynchronizing outcome of a single RVS CR shot. In fact, in the boxplots the large dispersion around the median value remains almost unchanged in the course of this protocol (Fig. 3.6B). Case II: Neither during the RVS CR shots nor during the subsequent pauses a sucient desynchronization is observed (Figs. 3.6C,D).

Figure 3 . 5 :

 35 Figure 3.5: Schematic summary of the CR stimulation protocols. (A) Protocol A: Spaced multishot CR stimulation with xed stimulation parameter. (B) Protocol B: Long singleshot CR stimulation with xed stimulation parameters. (C) Protocol C: Spaced multishot CR stimulation with demandcontrolled variation of the CR stimulation frequency and intensity. (D) Protocol D: Long singleshot CR stimulation with demand-controlled variation of the stimulation frequency and intensity (see text).

Figure 3 . 6 :

 36 Figure 3.6: Protocol A: Spaced multishot RVS CR stimulation with xed stimulation period T s . (A,C) Time evolution of the order parameter < R > averaged over a sliding window during 5 consecutive RVS CR shots with xed CR stimulation period. Dierent colors correspond to dierent networks. Stimulation parameters are unfavorable for anti-kindling in Case I (A,B) and Case II (C,D) (see text). (A,C) The horizontal solid red lines indicate the CR shots, while the horizontal dashed grey lines serve as visual cues. Spacing is symmetrical, i.e. CR shots and consecutive pauses are of the same duration. (B,D) Boxplots for R av , averaged over a window of length 100 • T s at the end of each pause, illustrate the overall outcome for all tested 11 networks. Case I: (K, T s ) = (0.[START_REF]Chaos Detection and Predictability[END_REF][START_REF] Popovych | What Can Simulations Contribute to Neuroimaging Data Analytics[END_REF]. Case II: (K, T s ) = (0.20, 28).

Figure 3 . 7 :

 37 Figure 3.7: Protocol C: Spaced multishot RVS CR stimulation with demand-controlled regular variation of the stimulation period T s and with demand-controlled variation of the intensity. (A,C) Time evolution of the order parameter < R > averaged over a sliding window during 5 consecutive RVS CR shots. If R av , the order parameter averaged over a window of length 100 • T s at the end of a pause, exceeds 0.4, the CR stimulation period of the subsequent RVS shot is decreased by T s → T s -1 ms (see text). Stimulation parameters are unfavorable for anti-kindling in Case I (A,B) and Case II (C,D) (see text). (A,C) The horizontal solid red lines indicate the CR shots, while the horizontal dashed grey lines serve as visual cues. Spacing is symmetrical, i.e. CR shots and consecutive pauses are of the same duration. (B,D) Boxplots for the time-averaged order parameter R av at the end of each pause, illustrate the overall outcome for all tested 11 networks. Case I: (K, T s ) = (0.30, 11). Case II: (K, T s ) = (0.20, 28).

Figure 3 . 8 :

 38 Figure 3.8: Protocol C in the presence of intrinsic variations of the ring rates caused by a modulatory low-amplitude current inputI var = A • sin(2π • f • t), with A = 1, f = 0.004 Hz (A, B), f = 4 Hz (C,D) and f = 20 Hz (E, F). Spaced multishot RVS CR stimulation with demand-controlled random variation of the stimulation period T s and with demand-controlled variation of the intensity. The low-amplitude variation I var is active during the entire simulations, respectively. Its low amplitude A = 1 ensures that the dynamics of the network is not drastically aected. (A, C, E) Time evolution of the order parameter < R > averaged over a sliding window during 5 consecutive RVS CR shots respectively. If R av at the end of a pause exceeds 0.4, the CR stimulation period of the subsequent SVS shot is decreased by T s → T s -1 ms (see text). (B, D, F) Boxplots for the time-averaged order parameterR av at the end of each pause, illustrate the overall outcome for all tested 11 networks respectively. The horizontal solid red lines indicate the CR shots, while the horizontal dashed grey lines highlight the two control thresholds (see text). Case I stimulation parameters are unfavourable for anti-kindling: (K, T s ) = (0.[START_REF]Chaos Detection and Predictability[END_REF][START_REF] Popovych | What Can Simulations Contribute to Neuroimaging Data Analytics[END_REF].

Figure 3 . 9 :

 39 Figure 3.9: Protocol A in the presence of intrinsic variations of the ring rates caused by a modulatory low-amplitude current input I var = A • sin(2π • f • t), with A = 1, f = 0.004 Hz (A, B), f = 4 Hz (C, D) and f = 20 Hz (E, F). Spaced multishot RVS CR stimulation with xed stimulation period T s . Same series of simulations and analysis as in Fig.3.8. (A, C, E) Time evolution of the order parameter < R > averaged over a sliding window during 5 consecutive RVS CR shots respectively. (B, D, F) Boxplots for the time-averaged order parameter R av at the end of each pause, illustrate the overall outcome for all tested 11 networks respectively. Spacing is symmetrical, i.e. CR shots and consecutive pauses are of the same duration. Case I stimulation parameters are unfavourable for anti-kindling: (K, T s ) = (0.30, 11).

  For both RVS CR and SVS CR, for several parameters tested the long-term outcome deteriorates with increasing stimulation intensity. Accordingly, based on our results, in pre-clinical and clinical applications stimulation at higher intensities should be avoided. Another important aspect refers to the more pronounced periodicity of SVS CR pattern. In previous papers (lacking a wider scan of the parameter space), SVS CR stimulation appeared to be superior to RVS CR stimulation [Zeitler and P. A. Tass 2015; Zeitler and P. A. Tass 2016]. Here, we showed that SVS CR stimulation decisively depends on the appropriate choice of the stimulation frequency. This sensitivity may signicantly reduce the performance in the presence of biologically realistic variations of the neuronal ring rates and might, hence, be the very reason, why the outcome of SVS CR stimulation is signicantly better for smaller numbers of sequence repetitions. In contradiction to the results obtained in networks without STDP [P. A. Tass 2003b; P. A. Tass 2003a; Lysyansky, O. V. Popovych, and P. A Tass 2011], CR stimulation may cause a fullblown anti-kindling without any phase resets of the subpopulations time locked to the corresponding stimuli (see S4 Fig. in

  intensities led to the computational development of spaced CR stimulation [O. V. Popovych, Xenakis, and P. A.Tass 2015] and two-stage CR stimulation with weak onset intensity[START_REF] Zeitler | Anti-kindling Induced by Two-Stage Coordinated Reset Stimulation with Weak Onset Intensity[END_REF]. The motivation behind these developments was to avoid side eects by substantially reducing stimulation intensities. In this context, it might turn out to be benecial that RVS CR stimulation causes sustained after-eects over a wide range of stimulation frequencies even at weak intensity (Fig.3.3). Accordingly, RVS CR stimulation might provide an appropriate stimulation protocol, in particular, if applied in an open-loop manner, without the ability to calibrate the stimulation parameters, especially the stimulation frequency by adapting it to the dominant peaks in the frequency spectrum of electrophysiological signals such as local eld potentials or EEG signals.By comparing spaced CR stimulation with xed stimulation parameters (Protocol A) and massed, continuous CR stimulation with equal integral duration (Protocol B) with a exible spaced CR stimulation with demand-controlled variation of CR stimulation frequency and intensity (Protocol C), and with a exible non-spaced CR stimulation with demand-controlled variation of CR stimulation frequency and intensity (Protocol D), we demonstrated[START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF]] that Protocol C enables to signicantly improve the long-term desynchronization outcome of both RVS and SVS CR stimulation, even at comparatively short integral stimulation duration. Remarkably, spacing alone (Protocol A) is not sucient to provide an ecient short-term dosage regimen. In fact, in particular cases vefold longer stimulation duration might even be more ecient than ve consecutive single CR shots with identical integral stimulation duration, at least for RVS CR stimulation[START_REF] Manos | How stimulation frequency and intensity impact on the long-lasting eects of coordinated reset stimulation[END_REF]. The low performance of pure spacing (Protocol A) might be due to the low number of single CR shots, here ve, as opposed to slightly larger numbers of CR shots, say eight, tested for the case of subcritical CR stimulation before [O. V. Popovych, Xenakis, and P. A.Tass 2015]. However, more important might be the approx. fty-fold longer stimulation and pause duration used for the spaced subcritical CR stimulation protocol. The long spaced subcritical CR stimulation protocol might be benecial for invasive application, such as DBS, and help reduce side-eects by substantially reducing stimulation current intake of the issue.

  focusing on stochastic resonance and related phenomena (e.g. [A. S.[START_REF] Pikovsky | Coherence Resonance in a Noise-Driven Excitable System[END_REF]).The number of stimulation sites and CR stimulation spatial decay was based on[START_REF] Lysyansky | Desynchronizing anti-resonance eect of m:n ON-OFF coordinated reset stimulation[END_REF]. In accordance to that study, adding more stimulation sites does neither lead to qualitatively dierent results nor does it improve the stimulation outcome.The short-term dosage regimen proposed here provides a closed-loop CR stimulation concept that enables to signicantly increase the robustness and reliability of the stimulation outcome. Our results motivate to further improve the CR approach by closed loop or feedback-based dosage regimen. Compared to the computationally developed initial concept of demand-controlled CR-induced desynchronization of networks with xed coupling constants [P. A. Tass 2003b; P. A. Tass 2003a],

  protocols[Manos T., Diaz-Pier S. and Tass P.A. Long-term desynchronization by coordinated reset stimulation in a neural network model with synaptic and structural plasticity, (in revision) 2020].

  Skokos, T. C. Bountis, and C. Antonopoulos 2007; Manos, T. Bountis, and C. Skokos 2013; C. Skokos and Manos 2016]. In [Manos, T. Bountis, and C. Skokos 2013], a study was carried out focusing on the dynamics of a barred galaxy model containing a disc and

a bulge component. Considering a TD analytical model extending a TI one

[START_REF] Manos | Studying the global dynamics of conservative dynamical systems using the SALI chaos detection method[END_REF][START_REF] Manos | Regular and chaotic orbits in barred galaxies -I. Applying the SALI/GALI method to explore their distribution in several models[END_REF][START_REF] Bountis | Complex statistics in Hamiltonian barred galaxy models[END_REF]

whose mass parameters of the bar and disc potential vary linearly as functions of time (the one at expense of the other). Two very general conceivable cases in barred galaxies were analyzed: (a) a model where the mass of bar grows, considering a common trend found in N -body simulations due to the exchange of angular momentum (see e.g.

[START_REF] Athanassoula | Morphology, photometry and kinematics of N -body bars -I. Three models with dierent halo central concentrations[END_REF][START_REF] Athanassoula | What determines the strength and the slowdown rate of bars? In[END_REF]

  Gyr, t = 4.2 Gyr, t = 7.0 Gyr and t = 11.2 Gyr. index. Note that the plot is in lin-log scale and the deviation vectors are re-initialized, by taking

							t = 1.4 (Gyr)
	100					t = 4.2 (Gyr)
							t = 7.0 (Gyr)
							t =11.2 (Gyr)
	of regular motion	40 60 80
	%
		20
		0	-0.35	-0.3	-0.25	-0.2	-0.15	-0.1
					energy
	Figure 1.7: Percentages of regular motion for the 2-d.o.f. frozen case as a function of the energy H
	for the dierent sets of parameters at t = 1.4

again k new random orthonormal deviation vectors, each time the GALI 3 becomes very small (i.e. GALI 3 ≤ 10

  2.1 Background and introductionQuantum chaos, or wave chaos, is the study of the phenomena in the quantum domain which corre-

spond to the classical chaos in the Hamiltonian systems

[START_REF] Stöckmann | Quantum Chaos -An Introduction[END_REF][START_REF] Haake | Quantum Signatures of Chaos[END_REF]

]. Although quantum motion (time evolution of the wavefunctions) of bound systems with purely discrete energy spectrum is ultimately (after a suciently long time, asymptotically) stable and regular, in fact almost periodic, it exhibits many features of the classical motion such as e.g. diusion in a chaotic domain, for times up to the Heisenberg time. The Heisenberg time, also called break time, is an important time scale in any quantum system, and is given by t H = 2π /∆E, where h = 2π is the Planck constant and ∆E is the mean energy level spacing, such that the mean energy level density is ρ(E) = 1/∆E. For times shorter than approximately t H the quantum diusion follows the classical chaotic diusion, but is stopped at larger times, just due to interference phenomena, which occur due to the wave nature of the underlying system, and are typically destructive.

Pictorially speaking, for times up to t H the quantum system behaves as if its evolution operator has a continuous spectrum, like the classical one has in the chaotic regime, but at later times it senses the discreteness of the spectrum. If the quantum diusion stops, while the classical chaotic diusion continues, we speak about the dynamical localization, or quantum localization or Chirikov localization, rst observed in time-dependent systems

  t we get β BR =0.421, β IZ (old)=0.416 and β IZ (new)=0.376, and based on the W (S) t we get β BR =0.421, β IZ (old)=0.401, β IZ (new)=0.350.

  2.10 we show three examples for the χ 2 as a function of the tted parameter β BR or β IZ .It is clearly demonstrated that Brody t is signicantly better than Izrailev. In Fig.2.11 we show the

	0.25							0.3								0.35
	0.2	Izrailev					0.25								0.3
																0.25
								0.2		Izrailev
	0.15															0.2
	χ 2						χ 2	0.15								χ 2
	0.1															0.15
								0.1
																0.1
	0.05							0.05								0.05
					Brody								Brody
	0							0								0
	0.1	0.2	0.3	0.4 β BR (PDF) & β IZ (PDF) 0.5 0.6	0.7	0.8	0.9	0	0.1	0.2	0.3	0.4 β BR (PDF) & β IZ (PDF) 0.5 0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9

β BR (PDF) & β IZ (PDF) Izrailev Brody

  Shepelyansky reduces the problem of the localization

	-0.6					MLEs (standard map)
	-0.8					f(x)=0.504*x+0.563 MLEs (random matrix)
	-2.2 -2 -1.8 -1.6 -1.4 -1.2 -1 Log 10 (standard deviation)					g(x)=0.504*x+0.383
	-2.4
	-2.6	1	1.5	2	2.5 Log 10 (t)	3	3.5	4

  Tass 2015], but with a t = 64 s CR-on period, is also associated with robust long-lasting antikindling and desynchronization for all CR stimulation intensity values K. Along the lines of our model analysis, the horizontal band corresponds to pronounced desynchronizing outcome at favorably weak CR stimulation intensities within a range of stimulation frequencies. However, we have to keep in mind that the discrete grid is not very dense. Hence, in order to investigate whether this conclusion is justied, we calculated C av and R av for all the integer period T s values for K = 0.20, ranging from f stim = 175%f 0 (T s = 10 ms) to f stim = 40%f 0 (T s = 40 ms). Fig 3.4 shows this ne-grained analysis. The boxplot for C av is shown in Fig 3.4A, and for the R av in Fig 3.4B. Note, in this gure the horizontal axis shows the CR stimulation period instead of the frequency. And it is sorted from larger to smaller values for an easier comparison between the two representations. The red and green dots indicate the reference stimulation period T 0 = 16 ms and intrinsic ring rate period T int =14 ms respectively. For T s ∈ [9, • • • , 24 ms] we observe robust anti-kindling eects. In contrast, for T s ∈ [25, • • • , 28 ms] many networks tend to be in a synchronized state, while for T s ∈ [29, • • • , 38 ms] the anti-kindling is found to be robust again, before nally reaching the largest T s value where the CR stimulation signals

	Another region with similar
	characteristics lies at the center of Figs 3.3A and 3.3B for intermediate stimulation intensity and
	frequency values.
	At a rst glance, among those two bands in Figs 3.3A and 3.3B, where dark color dominates
	suggesting long-lasting anti-kindling after cessation of CR stimulation, the horizontal band seems
	especially intriguing.

are not eective at all. In summary, at weak stimulation intensities favorable stimulation outcomes are achieved within wide ranges of the stimulation frequency. For further analyses of stimulation induced eects observed in particular ranges of the stimulation intensity/frequency parameter plane,

  In the majority of networks tested, these parameters did neither lead to acute nor long-lasting desynchronization after administration of a single CR shot (Fig.3.4) where T s = 28 ms corresponds to ≈ 50% of the intrinsic ring rate (or ≈ 36 Hz)]. For both cases, we investigate the order parameter < R > averaged over a sliding window for 11 dierent networks (marked with dierent color/line types, Figs 3.6A,C). Boxplots of the order parameter R av averaged over a window of length 100 • T s at the end of each pause demonstrate the overall stimulation outcome for alltested 11 networks (Figs. 3.6B,D).

3.5)

. Accordingly, the CR stimulation period T s remains constant, too. We study the stimulation outcome of only ve symmetrically spaced consecutive single CR shots. To this end, for both RVS CR and SVS CR (not shown here) stimulation we consider two unfavorable parameter pairs of xed CR stimulation period and intensity, respectively. One example refers to cases where CR stimulation induces acute eects, but no long-lasting desynchronizing eects. The other example concerns the case where CR stimulation causes neither acute nor long-lasting desynchronizing eects in a reliable manner.

RVS CR stimulation: Case I: (K, T s ) = (0.

[START_REF]Chaos Detection and Predictability[END_REF][START_REF] Popovych | What Can Simulations Contribute to Neuroimaging Data Analytics[END_REF]

. At a stimulation duration of 128 s these parameters caused only an acute, but no long-lasting desynchronization in the majority of networks studied (Fig.

3

.6), where T s = 11 ms corresponds to ≈ 127% of the intrinsic ring rate (or ≈ 91 Hz)]. Case II: (K, T s ) = (0.20, 28).

papers listed in the List of publications" section. All journals gave permission to use these papers in their published versions and to be inserted after each respective Chapter. However, that choice would substantially increase the thesis' size, and for this reason we opted for a more dense and concise style, still enabling an independent reading. A collection of all the aforementioned papers can be found either on the respective online journals' web-pages or alternatively on my personal web-page: https://sites.google.com/site/thanosmanos/Publications. At the end of the thesis one can nd the full bibliography for all chapters.

Acknowledgments

I would like to acknowledge and thank all my collaborators over the last recent years. I have enjoyed the personal and scientic interaction with all of them. The experience I have acquired working with them has helped me to learn a lot and further broaden my scientic research interests into new to me elds.

I would also like to thank Dr. Constanza ROJAS-MOLINA and Dr. Charalampos SKOKOS for reading early drafts of this manuscript and giving me their feedback which ultimately improved its quality.

I am also very grateful to Dr. Astero PROVATA, Dr. Rainer KLAGES and Prof. Andrey SHILNIKOV for having agreed to review this habilitation thesis as well as Prof. Alessandro TORCINI for agreeing to supervise it.

Finally, I wish to thank Meriem MATOUK for her help in the amelioration of the French of the Résumé.

Cergy-Pontoise

Thanos (Athanasios) MANOS September, 2020

List of publications

The results presented in this thesis are based on material published in the following articles (in chronological order).

The distribution of the localization measures

Following [Manos and Robnik 2015], we here present the results about the distribution of the localization measures. We restricted our analysis exclusively to the case K = 10, as this case is the closest to the normal diusion regime µ = 1. First we demonstrate that the localization measures L = 2/σ and l H are very well dened, linearly related and thus equivalent. In Fig. 2.12 we showed this in the diagram of the mean σ versus 2/ l H , where both averagings are over all eigenfunctions for matrices of dimension N = 3000, for 7 nearby values of k around k 0 , namely k = k 0 ± jδk, where j = 0, 1, 2, 3 and δk = 0.00125, for k 0 = 3, 4, 5, . . . , 19. In the Fig. 2.13, we show the relationship of the theoretical L in Eq. (2.16) and the mean value of the empirical 2/ σ for k 0 = 3, 4, 5, ..., 19. It is clearly seen in Fig. 2.13(a) that there are strong uctuations which we attribute to the fact that 2/σ has a certain distribution with nonvanishing variance, to be presented and described below, and that the theory of L resulting in Eq. (2.16) is too simple, as it corresponds only roughly to the value of 2/ σ . On the other hand, in Fig. 2.13(b) we see again that the two empirical localization measures are exactly linearly related. We should mention that in the cases of larger k > 19 the slopes σ are so small, and the localization too weak, that we cannot get reliable results. Therefore we limit ourselves to the interval 3 ≤ k ≤ 19.

We have thereby demonstrated that the empirical localization measures are well dened, while the theoretical prediction for their mean values is not good enough. The reason is that the localization measures of a given xed system (with xed K = 10 and k) have a distribution with nonvanishing variance, which is out of the scope of current semiclassical theories, as they do not predict this distribution and the corresponding variance. This nding, as the central result of the paper [Manos and Robnik 2015], is demonstrated in Fig. 2.14. The distributions are clearly seen to be close to a Gaussian, but cannot be exactly that, as σ is always a positive denite quantity. Its inverse, the localization length equal to L = 2/σ, has a distribution whose empirical histograms are much further away from a Gaussian, so that in this sense σ is the fundamental quantity. Indeed, it corresponds to the nite time Lyapunov exponent known in the theory of dynamical systems.

As l H and 2/σ are equivalent localization measures, the former one is expected also to have a distribution, which we demonstrate in the histograms of Fig. 2.15. We have also analyzed how the