{'fraction_non_alphanumeric': 0.028126680025414203, 'fraction_numerical': 0.011460827916524119, 'mean_word_length': 4.942064759692173, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transient Chaos Generates Small Chimeras

While the chimera states themselves are usually believed to be chaotic transients, the involvement of chaos behind their self-organization is not properly distinguished or studied. In this work, we demonstrate that small chimeras in the local flux dynamics of an array of magnetically coupled superconducting quantum interference devices (SQUIDs) driven by an external field are born through transiently chaotic dynamics. We deduce analytic expressions for small chimeras and synchronous states which correspond to nonchaotic attractors in the model. We also numerically study the bifurcations underlying the multistability responsible for their generation. Transient chaos manifests itself in the short term flux oscillations with erratically fluctuating amplitudes, exponential escape time distribution and irregular dependence of the escape time to initial conditions. We classify the small chimera states in terms of the position of the non-synchronized member and numerically construct their basin of attraction. The basin is shown to possess an interesting structure consisting of both ordered and fractal parts, which again can be attributed to transient chaos.


I. INTRODUCTION
The intriguing coexistence of synchronous and asynchronous groups in a collection of non-locally interacting identical oscillators, known as chimera or multichimera states [1,2], are as enigmatic as the Greek mythological creature with incongruous body parts that they are named after. Following the 2002 discovery of these states in the paradigmatic Kuramoto-Sakaguchi model [3], they have been speculated to be related to a diverse range of contexts. Examples include neural dynamics during unihemispheric sleeps of some birds and mammals [1,4,5,6] and the onset of epileptic seizures in human beings [7], fatal loss of coordination between human cardiac cells [1] to the laminarturbulent coexistence in the Couette flow [8]. Apart from these, chimeras are also numerically and experimentally found in diverse dynamical systems [1] ranging from mechanical [9,10,11,12,13], spin-torque [14] and chemical oscillators [1,15,16,17,18], certain neuron models [19,20,21,22,23,24] to coupled chaotic discrete maps [25,26,27]. A relatively recent addition, which will be the subject of our study in the later sections, is the dynamics of magnetic fluxes threading superconducting quantum interference devices (SQUIDs) in a quantum metamaterial [28,29,30], where chimera and multichimera states, predicted from numerical simulations, are presently awaiting experimental verification. Such ubiquity of chimeras, however, does not overshadow their mystery. Far from that, as it turns out [1,31], the morphology of the chimera states are as thought-provoking as are the conditions for their emergence, destruction and their dynamical properties, bifurcations and stability. In fact, the apparently irregular oscillations of the asynchronous oscillators in a chimera state render the system chaotic. This chaotic behavior of chimera states has been characterized by studies on Lyapunov spectra for varied interaction ranges of oscillators and system sizes [32,33,34], sensitivity to initial conditions [35] as well as typical bifurcation scenarios and Poincare sections [36]. By nature, this chaos can be spatial [32], temporal [35] or a subtle coexistence of both [36]. But it has also been shown that in finite size systems, chimera states have a finite lifetime [15,32,34,37], after which they collapse to a globally synchronous state, which implies that these states are examples of transiently chaotic states.
Transient chaos is a feature of certain dissipative systems, where the system evolves erratically for a finite time (called the 'escape time') due to the presence of non-attracting chaotic sets and finally settles down to one of the available (chaotic or regular) attractors which then determines the asymptotic dynamics [38,39]. Such transiently chaotic systems typically exhibit positive finite-time Lyapunov exponents [38,39], exponential distributions of escape times [38,39], irregular dependence of the escape time on initial conditions [38,39] and fractal basins of attraction [38,39]. There are ample evidences of real systems exhibiting such behavior [38,39], e.g., chemical reactions [38], NMR voltage signals [38,40] and driven pendulum [38] and the motion of particles of a fluid advecting in the wake of a cylinder [38,41] to mention few. In this regard, while the chimera states themselves are established as transiently chaotic dynamical states, there are no detailed studies relating the generation of chimera states to transient chaos. In this work, we have drawn that connection explicitly for the case of small chimeras. While it is a common experience from numerical simulations that chimeras are formed after some initial transient dynamics of the system which is typically initiated with a random configuration, there is no existing rigorous study of such dynamics. This is the problem we address in the current work and we explicitly show that such dynamics is a chaotic one over a finite time scale, which ultimately gives way to chimeras or synchronous states.
The system we shall focus on an artificially fabricated medium made of a lattice of magnetically coupled identical SQUIDs, each of which is basically a superconducting ring interrupted by a single Josephson junction. The material as a whole interacts with an incident electromagnetic field in some well-defined manner which, owing to the high nonlinearity and multistability [42,43,44,45] of the individual SQUIDs, can be tuned by varying experimentally accessible parameters of the system. The chimeras predicted in the flux oscillations of this material also show some interesting unique features, which make them distinct from most of those found in other dynamical models. Firstly, owing to the high degree of multistability of the individual SQUIDs, the chimera states are very sensitive to initial conditions and model parameter values [28,29,30]. Very slight changes can result in a state with a different amount of synchronization, which may not be a chimera at all. The chimera states here also exhibit the rare characteristic that they are seen to appear in even nearest-neighbor-coupled SQUIDs [29], where in most [1] (though not all [19,20,21,46]) other cases, they emerge in systems with long-range but non-global interaction. In the morphology also, these chimeras are peculiar; because they maintain one order of magnitude difference between the amplitudes of oscillations in synchronous and asynchronous groups, with the former usually having lower amplitudes of oscillation of its members [28]. Apart from these reasons, the study of chimeras in this system is further encouraged by the fact that the fabrication of such materials has been recently possible with Hypres Nb/AlO x /Nb junction process on silicon substrates (with the SQUID superconducting transition temperature at T c = 9.2K) [45,47,48], allowing for experimental probes to their peculiar dynamics. Thus, it is worthwhile to numerically investigate the nature of chimera states in more detail for this system.
In this work, we first systematically constructed the numerical phase diagrams for the various collective states of the system to identify the parameter regime for chimera states and understand the mechanism behind their emergence. For studies on transient chaos, we mainly focus on small system sizes. There we show that while the chimeras and synchronized states correspond to regular attractors whose analytic expressions we have derived, the route of the system towards them involves chaotic dynamics on finite time scales. A numerically efficient distinction between the initial chaotic dynamics and the asymptotic dynamics near the attractors in this system is facilitated by the fact that the latter dynamics is characterized by a steady amplitude of oscillation. Thus we have used the temporal fluctuations of oscillation amplitudes to determine the escape time for a specific initial condition. For varied initial conditions, the escape times vary erratically, with their distribution shown to be an exponential one. The final chimera states are distinguished by the position of the non-synchronized oscillator in the lattice. We plot the basins of attraction for each of them, which are shown to be fractal with additional special properties. Thus we obtain all the major signatures of transient chaos in the dynamics preceding the emergence of chimeras, namely, the short term flux oscillations with erratically fluctuating amplitudes, exponential escape time distribution, irregular dependence of the escape time to initial conditions and fractal basins of attraction of the attractors. Thus we establish that chimeras in a small-sized metamaterial ('small chimeras') are born through transient chaos. In future, this study can be generalized to other models and larger lattice sizes and the universality of these results be tested.
This study may turn out to be the key to understanding the mechanism of self-organization of chimeras in larger systems where the chimera state may itself be a chaotic transient [32,35]. In those systems, if we can quantitatively differentiate between the dynamics preceding a fully developed chimera state and the chaotic dynamics during the evolution of a fully developed chimera state, then we may try to find signatures of transient chaos in the former part as well. If we succeed, we can compare the natures of transient chaos in those two parts of dynamics and establish the chimera-sustaining systems as unique chaotic systems exhibiting sequentially transient dynamics and search for the underlying attracting and nonattracting invariant sets. In fact, in that respect, the typical observation for a chimera state that the position of the synchronous and asynchronous parts in the lattice is strongly dependent on the initial condition [3,35] is encouraging enough to probe the involvement of transient chaos. This observation seems to be in tune with our findings that the basins of attraction of small chimeras are fractal ones, which also lead to a loss of predictability of the system's final state in terms of the position of the non-synchronized oscillator. From more practical points of view, understanding how the chimera states emerge may lead to their efficient control. More importantly, we may be able to predict, delay or even inhibit their emergence altogether in cases where chimera states and partial synchronization are unwanted, for example, in epileptic seizures [7], and power grid desynchronization events [1]. In summary, this study, directly relating chaos and birth of chimeras, provides a useful understanding of the mechanisms at play behind the genesis of one of the most important and ubiquitous collective dynamical states.
The organization of this work is as follows: section II provides an introduction to the model, gives examples of a number of possible dynamical states and explains the phase diagrams of the model, so that we are able to characterize the system completely and choose the pertinent parameter regimes for small chimeras. The next section is devoted to analytic and numerical studies on smallest chimera states: their formation, classification, dynamics, basins of attraction and multiple manifestations of transient chaotic dynamics in them. Finally, section IV summarizes our findings and mentions few directions of future study.

A THE MODEL EQUATIONS AND COLLECTIVE STATES
In this section, we introduce the mathematical model of the metamaterial which we shall use below. The material is composed of a one-dimensional linear array of equally spaced N identical SQUIDs driven by an external magnetic field. The dynamical variables are the supercurrents I n (n = 1, ..., N ) in the individual SQUID rings and magnetic fluxes φ n (n = 1, ..., N ) threading them (Fig. 1). A current flowing through a SQUID ring generates a magnetic field around it, which contributes to the fluxes enclosed by the other rings and thus the SQUIDs engage in an effective long-range dipole-dipole interaction. The flux and current in each SQUID are further related and can be described by the the resistively and capacitively shunted junction (RCSJ) SQUID model [49]. These two facts enable one to write a set of coupled dimensionless equations for the fluxes only as [28] where n = 1, 2.., N , and φ e (t) = φ ac cosΩt is the flux enclosed by each ring due to an external time-dependent driving magnetic field and the coupling matrix for the inverse-cubic interaction between the SQUIDs is defined aŝ Here λ 0 is the coupling coefficient between nearest neighboring SQUIDs which depends on the lattice constant and other geometrical factors of the metamaterial, β is a parameter determining whether a SQUID is hysteretic or not and γ is the loss coefficient. One may consult Refs. [28] for the derivation of Eq. (1) and the normalization units of various quantities used in those equations. In numerical simulations of Eqs. (1) and (2), we initialize the fluxes with values randomly and uniformly chosen from the specified symmetric interval [−φ R /2, φ R /2] and also choose the initial flux velocity to be zero all the fluxes, i.e.,φ n (0) = 0 ∀n. In actual experiments, as suggested in [28], the random flux configuration can be achieved by thermal quenching, where the whole set-up is cooled through the superconducting transition temperature T c of the SQUIDs. We numerically solve Eqs. (1) and (2) using the fourth order Runge-Kutta method with a constant step size of dt = 0.02 unit as in Ref. [28]. The solutions show that, starting from a random flux configuration on the lattice, one usually obtains a number of different ordered states, e.g., (1) synchronous, where all the fluxes coherently oscillate as a mass, (2) cluster, where also the oscillations are mostly synchronous except for a few fluxes at scattered positions throughout the lattice, and (3) chimera and multichimera, where the fluxes spontaneously separate into one or more synchronous groups separated by asynchronous flux groups of similar sizes. Fig. 2 illustrates the diversity of these collective spatiotemporal pattterns.

B ORDER PARAMETERS AND PHASE DIAGRAMS
In order to have a complete idea about the behavior of the system and possible dynamical states in various parameter ranges, we need to identify the collective states of the system. In this section, we identify the parameter regimes corresponding to chimera, multichimera and synchronized states and elucidate the transitions between them, so that we gain an insight into how chimera states are formed and broken. For the necessary systematic identification and statistical characterization of the various collective modes of the system mentioned above, we will use three order parameters. They have been successfully used in earlier works to detect and characterize chimera states in varied range of systems, e.g., from coupled Kuramoto [50], and Lorentz oscillators [51] to Mackey-Glass delayed system [51] and Hidmarsh-Rose neurons [52]. The first of these measures is the traditionally used complex Kuramoto order parameter [53] R(t). At a particular instant of time, this is calculated from a particular flux configuration using the following equation, The magnitude of this quantity, averaged over a sufficiently large time interval during equilibrium (occurring typically after 5000 time units for N = 64) encodes the amount of global coherence in the system. In fact, |R| is 0 for asynchronous state, 1 for globally synchronous state and increases monotonically in between. Thus, for chimera, multichimera and cluster states, the system is only partially synchronized, so |R| takes some intermediate value between 0 and 1. However, since Eq. (3) shows that |R| is unchanged under arbitrary permutations of the fluxes φ i s, most of which obviously destroy their local correlations and chimera order, hence an intermediate value of |R| does not necessarily imply a chimera. To overcome this difficulty, two new order parameters were devised in earlier works [51,52], which are capable of detecting and distinguishing chimera and multichimera states by calculating the discontinuity of local order in a given flux configuration. These are called strength of incoherence (SI) S and discontinuity measure (DM) η. To calculate them, as in the earlier works [51,52], we start by defining Next, we bin the total N fluxes into N g bins of equal sizes N b = N/N g and calculate the time-averaged variance of z i for each bin by the usual equation: with k = 1, 2, ..., N g being the bin index. This quantity measures the fluctuation of fluxes in a single bin over the To identify the bins with a relatively lower fluctuation of z (corresponding to an enhanced local coherence), one defines a variable λ for each bin by the equation [51,52] which is 1 for a bin having variance σ lower than a sufficiently small threshold δ and 0 for other bins. Summing up this binary measure of local order over the bins, we can then define the strength of incoherence S as [51,52] so that, for complete flux synchrony S = 0, for complete asynchrony S = 1 and 0 < S < 1 implies chimera or multichimera states. To further differentiate between the last two states, we count the number of discontinuities separating coherent and incoherent regions in the flux profile as characterized by different values of λ. Thus the discontinuity measure is defined by summing up the differences between λ values of neighboring bins [51,52] using the following equation along with the periodic boundary condition λ Ng +1 = λ 1 . For both completely synchronous and completely random flux configurations, where the flux profile shows no discontinuity in the coherence profile and thus λ is same for all the bins, η takes the value 0. However, it is 1 for chimera state and a larger integer for multichimera states since both of them involve coexistence of domains with different amounts of coherence. Apart from these states, we also encountered cluster states. They are, however, harder to distinguish by these measures, and one often has to resort to more direct (often manual) ways [51] for their detection. In our case, we identified them by the fact that unlike chimeras, they do not show significant lowering of |R| from 1, but the SI and DM values do not correspond to synchronous states. Using this observation alongwith the three order parameters, we now have a computationally efficient way of systematically identifying various collective states from the dynamics of the flux profiles in the system.
To deduce the nature of collective states at different parameter values we plot the phase diagrams in β − γ plane for the three different order parameters |R|, SI and DM in Fig. 3. These show, for sufficiently low values of the damping parameter γ, the system is in an asynchronous state for any value of β. This can be explained with our observation in the next section that the time required for the system to reach equilibrium increases when the damping parameter is lowered and thus for smaller γ, steady states such as chimeras and synchronous states are not formed within a reasonable time. For higher γ values, as is seen from Fig. 3, the system follows a specific path towards its journey to synchronization through the growth of discontinuity in the system as revealed by the measure η (Fig. 3c). Fig. 3 shows, for sufficiently low values of β, the system is in chimera states. If β is increased further, discontinuities are introduced within the chimera states. The chimera gets fragmented and synchronous parts appear within its asynchronous portions and the system moves towards multichimera states as indicated by the increasing DM values. On increasing β more, the number of discontinuities in the flux profile increases further and system goes into cluster states with a large number of synchronous portions. Finally those synchronous parts merge together, reducing the number of discontinuities and leading the system to a completely synchronous state. This route shows how increasing the SQUID hysteresis parameter β affects the chimera state and ultimately destroys it. This transition, inferred from the amount of discontinuities in the flux profile as measured by η, is also confirmed and visualized from the results of direct simulation which we plotted in Fig. 2b-2f. In the next section, we shall reconcile this phenomena with the bifurcation diagram of the system with respect to β where we shall relate the formation of these discontinuities to the distribution of the fluxes among multiple possible solutions when the system becomes multistable for certain ranges of β. Thus, the phase diagrams studied here will guide us to find the suitable parameter regimes to look for chimera states with a smaller system size. Furthermore, as we showed, the phase diagrams clarify the dependence of the formation and annihilation of chimeras on the hysteresis parameter. This will complement our study on similar aspects in smaller systems with the bifurcation diagram discussed in the next section.
Since our system of SQUIDs is known to exhibit pronounced multistability [28,29,30], for consistency we used the same random initial flux configuration for all the points in the phase diagrams. To complement this, we have studied the multistability of the system explicitly to check the robustness of the phase diagrams. For that purpose, histograms of equilibrium values of |R| are plotted for different initial flux configurations (Fig. 4). The histograms show that, away from the various phase boundaries in Fig. 3c, the system is monostable and thus the prediction of various collective states distinguished by η in Fig. 3c and hence the mentioned route towards synchronization is statistically robust. only near the phase boundaries the system can exhibit multistability manifested by a wide distribution of the possible equilibrium values of |R| in Fig 4. In this respect, we note that the histograms of the statistical order parameter |R| will fail to detect another kind of multistable behavior if two microscopically distinct states yield identical values of the order parameter. For example, two chimera states with different lattice portions synchronized can have the same |R| value but are distinguishable. For small system size, this issue is studied in detail in the next section.

A CLASSIFICATION AND DYNAMICS OF POSSIBLE SMALL CHIMERA STATES
After the discovery and thorough analytic and experimental study of chimera states in a wide variety of systems, there is now a recent interest in realizing these (or at least similar) states in systems of minimum possible size. In this respect one usually adopts the notion of the so-called weak chimera states, defined by Ashwin and Burylko [54,55] as the states where two or more individuals of a system are frequency synchronized and the rest (one or more) oscillate with a different frequency than the synchronized group. Weak chimeras of this type have been numerically demonstrated in systems of three, four or sometimes a few more identical units for the cases of inertial Kuramoto model [56], Hansel-Mato-Meunier model [54] and Lang-Kobayashi-type model of delay-coupled lasers [57]; and experimentally shown in coupled pendula [58] and laser systems [59]. In this section, we shall show that the system of SQUIDS modeled by Eq. (1) with N = 3 adds to this list. In fact, evidently the system of 3 SQUIDs can be completely synchronized, or they can have chimera states, which may then be distinguished and classified by the position of the non-synchronized SQUID in the lattice. In that respect, Fig. 5 shows that we can obtain all of these possibilities from simulations of Eq. (1) with various initial conditions and parameter values. The plots also show that while the chimeras, cluster states and synchronous states arising in larger systems (as those shown in Fig. 2 and Ref. [28]) typically exhibit larger oscillation amplitudes for asynchronous fluxes and smaller oscillation amplitudes for synchronous fluxes, this is not necessarily the case for small chimeras. Indeed, as we found through direct simulations, large amplitude (magnitude around 1) and globally synchronous oscillation for φ ac = 0.01 can occur upto a lattice size of N = 10. This shows that smaller metamaterial lattices do sustain coherent oscillations with amplitudes much larger than that of the driving flux. This inspires us to seriously consider the system size while designing real materials.
Having seen the plethora of dynamical states in the numerical simulations of small-sized SQUID systems, we now aim for their analytical investigation. In that respect, the particular set of equations governing the dynamics of the model allows one to explicitly construct some exact relations for small chimera and synchronous states. As an example, let us first consider the case when the fluxes in SQUIDs 1 and 2 are phase-locked but the flux through third one does not synchronize with them. This is a genuine symmetry-broken solution since Eq. (1) is symmetric under exchange of the first and the third SQUID. Writing Eq. (1) explicitly with φ 1 = φ 2 = φ 0 (t) = φ 3 , Now, for φ 1 (t) = φ 2 (t), the right side of Eqs. (8) and (9) must be identical. This leads to the following algebraic relation of the flux of the synchronized pair with that of the third one where we have the symmetry of the matrix Λ −1 to simplify the expression. Using the explicit form of the coupling matrix Λ for N = 3, we rewrite Eq. (11) in a more compact form as As noted earlier, the particular form of the coupling matrix λ makes Eqns. (8) and (10) symmetric under the exchange of fluxes in SQUIDs 1 and 3. So, the corresponding relation for the chimera where SQUIDs 2 and 3 are synchronized is identical to Eq. (12) and reads A similar straightforward analysis fails for the third type of chimera state, but that type does exist, as seen from direct simulations (Fig. 5). However, unlike the first two types of chimeras just described, this solution preserves the symmetry of Eq. (8)- (10). Nonetheless, Eqs. (12)- (13) show that these exact results for the simplest chimera state contain all the elements of the Ashwin-Burylko criterion described above and thus theoretically confirms the possibility of small chimeras in this model. A similar analysis can be continued to systems of four and more SQUIDs in future, where there will be more diverse types of chimeras as found in Ref [57]. The above analysis for chimeras can be further continued to deduce the condition for complete synchronization. For this, we set φ 1 = φ 2 = φ 3 = φ 0 in Eqs.
where we have expanded the sine in Eq. (14) upto linear order for small φ ac and used the fact that γ << Ω so that tan −1 Ω γ ≈ π 2 . It is easy to check that the two sides of Eq. (15) are of the same order for most typical parameter values that experimentalists work with. Hence ideal or approximate synchronization is quite common in the system and they coexist with chimera states in almost all of the parameter regime studied here. This inspires us to explore the basins of attraction for different chimera and synchronized states and this is done below.
The above analytic results fail to predict the stability of the various states. To check their stability, we take resort to plotting the numerical bifurcation diagram (Fig. 6) corresponding to flux oscillation amplitudes. This is done at the equilibrium regime when the amplitudes of all the flux oscillations become steady over time. As we see from the figure, there is only one stable amplitude possible for the flux oscillations in the SQUIDs when the hysteresis parameter β has sufficiently large values. This situation corresponds to completely synchronized states with all the three fluxes oscillating with very small amplitudes. As β is decreased from this value, the system undergoes a bifurcation which generates an additional branch of larger possible amplitudes for the fluxes. With lowering of β both the branches further bifurcate and thus the system becomes multistable with respect to the possible oscillation amplitudes. This continues until as β approaches zero, when the lower branch disappears and only one solution branch remains. This multistability opens up the opportunity for the amplitudes of the three fluxes to stay on different branches and thus yield a chimera state. For larger systems, such multistability for intermediate values of the hysteresis parameter can result in the formation of discontinuities in the flux profiles as indicated in the phase diagrams of Fig. 3. Thus the bifurcation diagram showing the mechanism for generation of small chimeras may actually be valid for larger system sizes as well. On the other hand, this inspires one to look for the special features of small chimeras, which we discuss below, in larger systems as well.

B SIGNATURES OF TRANSIENT CHAOS
The synchronized and chimera states described in the last subsection act as attractors of the system, since they represent the asymptotic steady dynamics. As we have found, they are nonchaotic for a wide regime of parameter values, like the ones plotted in Fig. 5. However, an investigation to their full time series (Fig. 7) reveals that initially the dynamics is apparently irregular with randomly and widely fluctuating amplitudes of oscillation of the fluxes. This leads us to suspect the involvement of transient chaos in the emergence of chimera and synchronous states for the system. To confirm this, we carry out a number of studies. The simplest signature of transient chaos is the statistics of escape times [38,39] and dependence of the same on initial conditions. However, to obtain the escape time, one needs to divide the total time evolution of the system into two distinct parts corresponding to initial transient dynamics and asymptotic dynamics near the attractors. In principle, these two parts can be distinguished by finite time Lyapunov spectra which typically shows different characteristics for the two different regimes [38]. However, the current system allows for a more straightforward and numerically efficient way for this distinction. This is based on the observation that (Fig. 7) the amplitudes of flux oscillations are stable for sufficiently long time only when the system is close to the attractor and show pronounced fluctuations during the transient period. As a quantitative measure, we have detected the standard deviation of a sufficiently large (here 20) number of consecutive amplitudes and detected escape only when this is lower than a small threshold value (here 0.005). The escape time thus found is plotted against varied initial conditions. The resulting plot (Fig. 8a) indicates the irregular dependence of the escape time on initial conditions, as is the usual case [38] for transiently chaotic systems. In particular, the discrete points on the plot with long escape times are the set of points which belong to the immediate neighborhood of the chaotic set responsible for the transient chaos. We also find the distribution of escape times for initial conditions over a chosen region in phase space (see Appendix A for details of the region). As another demonstration of transient chaos, for sufficiently long transients, the asymptotic part of the distribution is shown to be well-fitted by an exponential curve (Fig. 8b) yielding the time scale for the escape to the attractors. The time scale is dependent on the system parameters, in particular, on the damping parameter γ which directly controls the rate of relaxation of the system. This relation, as is evident from Fig. 8c, can be fitted to a power law, with an exponent of approximately 0.92. This increase of the time scale on lowering the damping parameter can also be used to explain our earlier observation that for sufficiently small values of γ, the system fails to reach any collective state in reasonable time. This is yet another example of how the insights from the dynamics of small metamaterial systems can be successfully applied to larger system sizes as well.

C BASIN OF ATTRACTION FOR A SYSTEM OF THREE SQUIDS
In the previous sections, we mentioned that the small chimera states and synchronized states are possible attractors at certain parameter values. However, the system governed by Eq. (1) is well-known to show pronounced multistability [28,29,30] which we explained with the bifurcation diagram of Fig. 6. For larger systems this can correspond to different initial conditions leading to a huge possibility of different chimera states, which are distinguished by the positions of the asynchronous and synchronous oscillations in the lattice and possibly nonchimera (e.g. synchronous or cluster) states. But at the level of only three SQUIDs, as explained in the earlier subsection, the only five possibilities are complete synchronization, three types of chimeras and the asynchronous states. In this section, we show that in this small system, multistability plays a dramatic role in determining the final collective state. To demonstrate this, in Fig. 9a and 9b, we plot the basins of attraction of the system projected in the planes of initial conditions spanned by φ 1 (0) − φ 2 (0) (with φ 3 (0) fixed at 0) and φ 1 (0) − φ 3 (0) (with φ 2 (0) fixed at 0) both withφ i (0) = 0 for i = 1, 2, 3 (the actual phase space of the second-order system Eq. (8)- (10) with N = 3 is six-dimensional). To plot the basins, we detected a chimera if |R| calculated for a state and averaged between times 4800 − 5000 was between 0.4 − 0.8. A greater time-averaged |R| was marked as synchronous, while a lesser was marked as some other state. In this respect, we note that, unlike the case for bigger systems, |R| can now successfully identify a small chimera since it yields an intermediate value only if any two of the fluxes are synchronous and that always corresponds to a chimera state.
The basins of attraction are seen to be composed of both regular and apparently random components and have interesting structures. The regular parts would correspond to predictability of final states, while for the other parts, a minute change in the initial conditions may lead to a very different final state. The relatively regular parts in the plots occur in the vicinity of the origin and along the lines on which two initial fluxes are equal, except for the line φ 1 (0) = φ 3 (0) (Fig. 9a and 9b). In the former case, all three fluxes start with magnitudes very close to each other, so that they ultimately synchronize. In the latter cases, any two of the fluxes start with values close to each other, such that the initial configurations themselves resemble a chimera state. However, as is seen from the figures, this may or may not ultimately lead to a chimera state since they can be totally synchronized also. Interestingly, the initial conditions responsible for synchronized and chimera states in such cases are seen to be arranged in a regular manner with well-defined periodicity which also slowly varies with the distance from the origin ( Fig. 9a and 9b). This regularity vanishes sufficiently away from the origin, which indicates that, for initial fluxes of significantly larger values, the system is not predictable at all, even if two of the initial fluxes are equal. This extreme sensitivity to initial conditions is a well-known property of chaotic systems [60,61], where these are responsible for fractal basin boundaries. In those cases, the trajectories of the system behave chaotically near the vicinity of a nonattracting chaotic set in the phase space and, finally escape it via different routes and consequently terminate on different attractors. So, this behavior can be yet another manifestation of the involvement of initial transient chaos in the system's final relaxation to an attractor. This idea can be carried over to systems of larger sizes with a richer diversity of attractors and the involvement of transient chaos to the generation of chimera states in other systems can be studied. In fact, it is well-known [3,35] that, for spatially extended systems, the position of the synchronous portion in the lattice is strongly dependent on the initial condition. This might be an indication that transiently chaotic dynamics is involved in shaping a chimera's spatial profile in any general system.
The study of fractal basins in dissipative systems is a relatively less explored topic with generic predictions being done only recently [62,63,39]. In this respect, a particularly interesting type of fractal basin is the so called Wada basin where basins corresponding to three or more attractors meet at the same boundary [61,62]. This interesting topological structure has been discovered in certain chaotic systems [61,62,64] and conditions for its appearance are being developed [61,62]. In our system, we check the fractality and Wada property of our basins through direct numerical simulation. To do this, we have adapted the method in Ref [62], while a different method is discussed in Ref. [65]. To confirm the fractal and Wada properties, we have selected a portion of the basin boundary and zoomed in parts of it in successive steps. As is seen from the typical example in Fig. 9c, as we keep on increasing the resolution, an apparent boundary between two basins actually shows portions of other basins inside it at higher resolution and it continues until very small scales. To ensure the numerical precision, we choose a step size of dt = 0.002 in doing this. However, continuing this to arbitrary fine length scales would really be an interesting theoretical study, but not of interest from experimental considerations where one has limited control over initial flux adjustments. Nonetheless, the studies in this section establish the system as an experimentally realizable example of dissipative system showing sensitive dependence on initial conditions and fractal Wada basins, thereby encouraging the experimental study of dissipative chaos.

IV. CONCLUSIONS
In the preceding sections, we described how the small chimera states emerge out of transient chaos in a model of nonlinear metamaterial. These studies, in one hand, encourage to study the mechanism behind emergence of chimeras in other models, and in the other, establish the superconducting metamaterial as a unique system to experiment with nonlinear dynamics. We now conclude by mentioning some limitations of the study and offering possible directions of future works.
The effective model used to simulate the material has many idealizations and accommodations of realistic effects may have significant effect on collective dynamics. Firstly, our considerations here neglected any imperfections in the system coming from, e.g, critical current fluctuations of the SQUIDs [66], intrinsic low-frequency flux noises due to spin-diffusion [67] or local magnetic impurities [68], and other 1/f noise sources [69]. Furthermore, the characterizations of all the individual SQUIDs in the material may not be identical and so ideally there should be a distribution of the β,γ and λ 0 parameters. The external magnetic field may not be temporally coherent and its spatial variations may not be negligible. Furthermore, as was shown in Ref. [29], considerations of nonzero initial flux velocities allow more diversity and multistability in the dynamics with prominent multichimera states containing varied number of synchronous groups. We neglected this possibility throughout our studies.
Our present work complements the identification of chimera states as transiently chaotic states themselves and shows that even in the self-organization of these states, transient chaos is involved. To check if this relation is robust, one needs to study other chimera-sustaining systems and discern the characteristics of the transients leading to chimera states there. Further sophisticated studies may involve investigations of the finite-time Lyapunov spectra at various times. This would be particularly needed if the chimera states themselves are chaotic and then one may like to observe the differences in the Lyapunov spectra corresponding to the transient part, the fully developed chimera states and (if present) the state obtained after the chimera state collapses. For extended systems, the transient chaos may be spatial, temporal or spatiotemporal and thus offer much richer dynamical possibilities. Since currently both chimeras and chaos enjoy a sound establishment as prominent contemporary research topics and both exhibit a certain level of ubiquity in natural and model systems, we hope that the newly found relation between them will open up novel directions of future work in the field of nonlinear dynamics.