Chaos in hydrodynamic BL Herculis models.
Abstract
We present nonlinear, convective, BL Hertype hydrodynamic models that show complex variability characteristic for deterministic chaos. The bifurcation diagram reveals a rich structure, with many phenomena detected for the first time in hydrodynamic models of pulsating stars. The phenomena include not only period doubling cascades en route to chaos (detected in earlier studies) but also periodic windows within chaotic band, typeI and typeIII intermittent behaviour, interior crisis bifurcation and others. Such phenomena are known in many textbook chaotic systems, from the simplest discrete logistic map, to more complex systems like Lorenz equations.
We discuss the physical relevance of our models. Although except of period doubling such phenomena were not detected in any BL Her star, chaotic variability was claimed in several higher luminosity siblings of BL Her stars – RV Tau variables, and also in longerperiod, luminous irregular pulsators. Our models may help to understand these poorly studied stars. Particularly interesting are periodic windows which are intrinsic property of chaotic systems and are not necessarily caused by resonances between pulsation modes, as sometimes claimed in the literature.
keywords:
hydrodynamics – methods: numerical – chaos – stars: oscillations – stars: variables: BL Herculis1 Introduction
Chaotic dynamics is present in many astrophysical systems and stellar variability is not an exception, although in this case, chaos was studied mostly in the context of hydrodynamic models of large amplitude pulsators. Buchler & Kovács (1987) and Kovács & Buchler (1988) found a chaotic behaviour in their radiative typeII Cepheid models (W Vir and RV Tau). In depth analysis of chaos in these models was conducted by Serre, Kolláth & Buchler (1996) and Letellier et al. (1996). Also Buchler & Moskalik (1992) found chaotic behaviour in two sequences of radiative BL Hertype models, however did not analyse the phenomenon. Recently chaotic behaviour was reported in convective hydrodynamic models of BL Her stars (Smolec & Moskalik, 2012) and RR Lyrae stars (Plachy, Kolláth & Molnár, 2013).
On observational side chaos was detected in type II Cepheids of RV Tau type (R Scuti and AC Her; Buchler, Kolláth & Serre, 1996; Kolláth et al., 1998) and in several semiregular variables (Buchler, Kolláth & Cadmus, 2004), and in Miratype variable (Kiss & Szatmáry, 2003).
In this paper we report on the chaotic behaviour we have found in a sequence of nonlinear convective BL Her models. For the first time in stellar pulsation modelling we clearly demonstrate the appearance of dynamical phenomena well known and common to classical chaotic systems, both discrete (e.g. logistic or Hénon maps) and continuous (e.g. Rössler or Lorenz equations). In all these systems the basic route to chaos is through a cascade of period doubling bifurcations also present in our models and in earlier studies of radiative typeII Cepheid models (Kovács & Buchler, 1988). In addition, our models display a full wealth of dynamic behaviour characteristic for deterministic chaos. Within chaotic band we find several windows of nonchaotic variation (windows of order), with stable period limit cycles. These windows are either extremely narrow or relatively large. In the latter case the periodic window is preceded by the typeI intermittent behaviour, till the periodic cycle is born through the tangent bifurcation. This periodic cycle again undergoes a series of period doubling bifurcations en route to chaos. The interior crisis bifurcations, in which separate chaotic bands merge, leading to the abrupt increase of the attractors volume are detected in our models, as well as crisis induced intermittency and typeIII intermittency.
Chaos was not detected in any BL Her star so far. In our opinion however, these models are important for several reasons. (i) Chaos does occur in largerluminosity typeII Cepheids – RV Tau stars, as well as in semiregular variables. Our models may shed more light on variability of these poorly studied stars. (ii) We initiated the survey of nonlinear convective pulsation models of typeII Cepheids extending to the highest luminosities (RV Tau domain) in which chaotic variability is expected, as previous radiative models and observations indicate. In the present paper we introduce and test the methods to study chaos in such models. (iii) The striking similarity between our hydrodynamic models of pulsating stars and even the simplest chaotic systems, like logistic map, is noteworthy, indicating that many very different systems may share the same dynamical properties. (iv) Finally, although the chaotic behaviour was not detected in any BL Her star so far, we cannot exclude such possibility in the future. We note that the period doubling effect in these stars was predicted by Buchler & Moskalik (1992), based on radiative hydrodynamic models, but it took 20 years to discover the effect in the first star of this type (Soszyński et al., 2011; Smolec et al., 2012).
In Section 2 we summarize the properties of the logistic map, which will help us better understand phenomena occurring in our hydrodynamic models. The reader familiar with the chaos theory may safely skip this Section. In Section 3 we briefly describe the computation and basic properties of the models. In the next Sections we present detailed analysis of the models showing both chaotic and periodic variation, including discussion of the largest Lyapunov exponents (Section 5). We discuss our results in Section 6 and comment on observability of the chaotic phenomena in Section 7. Summary in Section 8 close the paper.
Initial results of this study were reported in the conference proceedings of IAU Symposium No 301 (Precision Asteroseismology), Smolec & Moskalik (2014).
2 Logistic map
In this Section we briefly discuss the properties of the logistic map, which is the simplest 1D discrete system showing deterministic chaos. The map is defined as:
(1) 
where is a parameter. We are interested in a range as then the iterates of eq. (1) are bounded, . For iterations converge to 0 and for they diverge. The great advantage of the logistic map is its simplicity – most of the properties, e.g. bifurcation points, fixed points and their stability, may be computed analytically. Full understanding of the mechanisms behind the observed behaviours is thus possible. At the same time logistic map displays nearly all types of behaviours characteristic for deterministic chaos in more complex, continuous and higher dimension systems, that are also present in our hydrodynamic models. In the case of our models analytical approach is not possible and we are left with the complex output of pulsation code. Comparison with results discussed in this Section allows for a better understanding of our hydrodynamical models. The properties of the logistic map described below may be found in numerous textbooks (e.g. Petigen, Jürgens & Saupe, 2004) and original articles (e.g. May, 1976), and are given here without derivation.
Depending on the value of the iterations of any initial (trajectory) either converge to a periodic cycle, or not, and then the chaotic attractor is present. In Fig. 1 we show the bifurcation diagram – a possible longterm values of (initial iterations omitted) as a function of . The plot is a stack of greyscaled histograms: for each value of , we computed several thousand iterations of eq. (1) and calculated the probability with which the iterations fall to one of the 200 bins into which the interval was divided. Grey bands of chaos are clearly visible, as well as windows of order in which stable periodic cycles are present. The right part of the figure shows the zoom of the largest, period3 cycle window.
To visualise the attractor it is useful to construct the first return maps, i.e. plots of vs. (omitting the initial transient). For the discussion below it is instructive to consider the evolution of return map as is increased. In a range it is shown in the animation that may be found in the online version of this article as additional supporting information.
The consecutive iterates of eq. (1) may be constructed geometrically, using the plot of vs. as illustrated in the left panel of Fig. 2 for (thin red line). For initial value one plots the vertical line towards the and from that point the horizontal line towards the diagonal, and then repeats the procedure to get a trajectory. In the discussed case it converges to the fixed point, (period1 cycle), which is located at the crossing of and the diagonal. Its stability depends on the slope of at the intersection. If the modulus of the slope is less than 1 the iterations converge towards and the period1 cycle is stable. It is the case for : iterations of any converge towards . For the slope is steeper than 1, the period1 cycle becomes unstable and stable period2 cycle () is born through the period doubling (pitchfork) bifurcation. The iterates of eq. (1) alternate between two values. The appearance of the period doubling bifurcation is illustrated in the middle panel of Fig. 2 with the help of vs. plot. For intersects the diagonal at a single point corresponding to a stable period1 cycle. At the bifurcation point () is tangent to the diagonal, and for larger , intersects the diagonal at three points. The middle point corresponds to the unstable period1 cycle, which is a degenerate case of period2 solution. The two other points have the same slopes and correspond to the stable period2 cycle.
Following the same scenario, at another period doubling bifurcation occurs giving rise to stable period4 cycle (and period2 cycle loses its stability). The cascade of period doubling bifurcations will finally lead to chaos. The extent of the domain of stable period cycle () decreases as increases and the ratio approaches the Feigenbaum constant () which is a universal constant for many other systems too (Feigenbaum’s universality, see Feigenbaum, 1983). Beyond (the accumulation point) chaos appears for the first time. The described period doubling route to chaos is schematically illustrated in the upper part of Fig. 3 with stable and unstable cycles marked with black solid and red dashed lines, respectively.
Beyond the accumulation point the chaotic domain extends, which however, is densely packed with windows of order – stable period cycles. To understand how these are born we focus our attention on the most prominent period3 window (zoomed in the right part of Fig. 1) and the iteration of (rightmost panel of Fig. 2). Before the period3 cycle is born the system is chaotic and the three vertices of approach the diagonal. At/beyond the bifurcation () they touch/intersect the diagonal and give rise to a pair of period3 cycles, of which one is stable and one is unstable, as is apparent from the analysis of the slopes of at the intersection with the diagonal. The stable period 3cycle will soon undergo a series o period doubling bifurcations en route to chaos, just as described in the previous paragraphs, and as is clearly visible in the right panel of Fig. 1. This scenario is schematically plotted in the bottom part of Fig. 3.
Two interesting phenomena occur at the edges of the just discussed period3 window (and other periodic windows as well). Before the tangent bifurcation occurs we observe the intermittency, illustrated in Fig. 4. The bottom panels show the consecutive iterations of eq. (1) for three different values of and the top panels show the third return maps, i.e. plots of vs. . Obviously the points fall along the curve. Just before the bifurcation (left and middle panel) evolution of the system is characterized by long intervals of almost periodic behaviour interrupted by shorter bursts of chaos. This is the intermittent behaviour first analysed for the Lorenz equations by Mannevile & Pomeau (1979), followed by an indepth analysis by Pomeau & Mannevile (1980). As control parameter increases the almost periodic intervals become longer (cf. left and middle panels in Fig. 4) up to a critical value of at which tangent bifurcation occurs and stable period3 cycle is born (and unstable period3 cycle as well; right panel in Fig. 4).
To discuss the appearance of intermittency we focus attention on the left panel of Fig. 4 (inset). Constructing the trajectory for each third iterate of eq. (1) geometrically, one must fall into the intermittent channel – a narrow region between the vertex of and the diagonal – in which iterations must be trapped for a while (thin zigzag in the inset). Similar channels are also present at the two other vertices approaching the diagonal and iterates of eq. (1) fall consecutively into the three channels. Closer to the bifurcation, narrower the channels and longer the iterations are trapped within, with apparently more periodic variation (middle panel of Fig. 4, inset). As the iterations leave the intermittent channels a chaotic burst is observed.
In a broader context, intermittency is one of the routes to chaos (for a review see Eckmann, 1981), characterized by sporadic switching between qualitatively different behaviours, the laminar (periodic) behaviour and chaotic bursts. Intermittency is associated with a bifurcation in which stable periodic cycle becomes unstable. Depending on the bifurcation in which the stability is lost, three types of intermittency were distinguished by Pomeau & Mannevile (1980). The typeI intermittency is related to a tangent bifurcation in which stable and unstable limit cycles collide and both vanish. It is the case for the just described intermittency in the logistic map. TypeII intermittency is related to the subcritical Hopf bifurcation (in the Hopf bifurcation the stationary solution bifurcates into periodic orbit)^{1}^{1}1We follow the convention adopted e.g. in Seydel (2010) and call the bifurcation subcritical if the stable solution exists on one side of the bifurcation point only. If the stable solution exists on either side, the bifurcation is supercritical.. TypeIII intermittency appears together with subcritical period doubling bifurcation, in which the stable limit cycle collides with the unstable perioddoubled cycle. A summary of these scenarios may be found in the original paper by Pomeau & Mannevile (1980) and e.g. in Becker et al. (1999). Hydrodynamic models to be discussed in the coming sections display both typeI and typeIII intermittency.
The stable period3 cycle born in the tangent bifurcation will soon undergo a series of period doubling bifurcations which will create three separate chaotic bands as is well visible in Fig. 1 and in Fig. 5 (left and middle panels), which is the same as Fig. 4 except we focus on the right edge of the period3 window. These three bands will merge in the interior crisis bifurcation, first described by Grebogi, Ott & Yorke (1982), in which the volume of the attractor changes suddenly. The crisis bifurcation occurs as the three chaotic bands hit the unstable period3 cycle (born in the tangent bifurcation) which scatters the trajectories into previously unvisited regions – Fig. 5 the rightmost panel. This is well visible in the right panel of Fig. 1 in which unstable cycles of period 3 (and period 6) are plotted with the dashed lines and in Fig. 5 where location of unstable fixed points of period3 is indicated with ’+’ signs.
Directly after the occurrence of crisis bifurcation, crisis induced intermittency is observed (well visible in the bottom right panel of Fig. 5). The trajectory is confined in the region of the former, precrisis attractor, with sporadic excursions out of it.
3 Hydrodynamic models
All hydrodynamic models analysed in this paper were computed with the Warsaw nonlinear convective pulsation codes (Smolec & Moskalik, 2008). Numerical parameters of the models (zoning) are the same as in our previous papers (section 3 of Smolec et al. (2012) and Smolec & Moskalik (2012)). The physical parameters of the models and parameters of the turbulent convection model (Kuhfuß, 1986) are the same as in Smolec & Moskalik (2012). In particular , and . We focus on a single sequence of models with the same luminosity, and varying effective temperature, , which is our control parameter through this paper. The models cover a strip extending over K, from K to K (corresponding periods of the fundamental mode are d and d). The maximum temperature difference between the consecutive models is only K and in the most interesting domains the difference is as small as K. In Fig. 6 we show the location of our models in the HR diagram (thick horizontal line), together with the location of models that we have studied in Smolec & Moskalik (2012) that show periodic and quasiperiodic modulation of pulsation (thin horizontal lines).
In nonlinear computations, the initial static model was perturbed with the velocity profile, and was integrated for at least pulsation cycles with time steps per pulsation cycle. By default, all the models were initialized in the same manner, with velocity eigenfunction of the fundamental mode, scaled to match the 4.5 km s surface velocity. To check for the possible dependence of results on the initialization, in particular to check whether e.g. two attractors are possible for the same model (hysteresis), several models were initialized in a different manner (with larger surface velocity or with a mixture of the fundamental mode and first overtone eigenfunctions). In all considered cases, the computations converged to the same attractor, only the length of the initial transient phase was different. To check the longterm stability of the computed attractors (either periodic or chaotic), we have computed up to pulsation cycles for a few models. In all cases the attractor emerging from the calculation of first pulsation cycles remained stable.
We note that in these models the eddyviscous dissipation is strongly decreased, (for equations see Smolec & Moskalik, 2008). It results in significant pulsation amplitudes and, as our treatment of radiation is very simple (diffusion approximation) and model mesh is fixed, in erratic light curves with spurious spikes (see section 4 in Smolec & Moskalik (2012) for detailed discussion of this point). Therefore in this paper we analyse the radius variation only, which is smooth. In Fig. 7 we plot a section of typical time series for model showing chaotic variability ( K).
4 Analysis of the BL Her hydrodynamic models
In this Section we analyse the radius variation of our hydrodynamic models.
4.1 Bifurcation diagram
Deterministic chaos is present in our models beyond doubt. The conclusion is unavoidable once the bifurcation diagram is plotted – Fig. 8. The diagram is constructed in a similar way as in the case of logistic map (Fig. 1). For each value of our control parameter – the effective temperature – we computed the probability that the maximum radii, , falls into one of the bins into which the range of radius variation in our models () was divided. The stack of greyscaled histograms is displayed in Fig. 8 and shows a striking similarity to classical chaotic systems. One can pick other parameter than maximum radius to construct the bifurcation diagram, but results are qualitatively the same. As compared to Fig. 1 the bifurcation diagram looks rough which results from smaller resolution in the control parameter and decreased number of bins along vertical axis, necessary to have a reasonable statistics in each bin.
On both sides of the computation domain singleperiodic (cycle1) pulsation in a fundamental mode is present. The chaotic bands are reached through a series of period doubling bifurcations both from the cool and the hot side. Qualitatively the same scenario is observed for the Gauss (or mouse map), see e.g. Hilborn (2000). The chaotic bands are separated with periodic windows of order. The largest period3 window is nearly a onetoone copy of the just discussed counterpart seen in case of the logistic map. Within chaotic bands pronounced structures are clearly visible as well – the probability of hitting particular bins by maximum radii is not equal. Some values of maximum radii are clearly preferred as indicated with darker bands, migrating across the bifurcation diagram as control parameter changes.
4.2 Period doubling cascade en route to chaos
Our hydrodynamic models display two pronounced period doubling cascades that lead from a singleperiodic fundamental mode pulsation (period1 cycle) to chaos. The first cascade extends on the cool side of the computation domain. We clearly observe the appearance of period2 cycle (first in the K model and present up to K model) and period4 cycle (first in the K model and present also in the K model). Further periodic cycles are not resolved in our model grid with K step in effective temperature. The maximum radii of the K model are bounded within eight well separated chaotic bands, which, as effective temperature is increased, merge smoothly into one chaotic band (at K).
The other period doubling cascade, on the hot side of the computation domain, extends over much larger temperature range, but otherwise it is a mirror image of the just discussed cascade. The consecutive period doubling bifurcations occur with decreasing control parameter. For such situation, the terms periodhalving or inverse period doubling cascade are in use. Here we describe the route from order (period1 cycle) to chaos as effective temperature is decreased. Period2, period4, period8 and period16 cycles are all clearly detected. They appear for the first time in models with temperatures K, K, K and K, respectively. Period8 and period16 variation is detected in only one model each. In K model maximum radii are bounded within four separate chaotic bands.
It is interesting to check whether the lengths of the domains of the consecutive period cycles follow the Feigenbaum scaling. Our model grid is too coarse to exactly pinpoint the bifurcation points and hence only a rough estimates are possible. Assuming that the bifurcation occurs halfway between the neighbouring period and period models we get (for the cascade on the hot side): K, K and K for the extents of the period2, period4 and period8 domains, respectively (see Section 2). The ratios (assuming K error in the estimate of the bifurcation point) are K and K not significantly different from the Feigenbaum constant () toward which the ratio converges as for the logistic map and other iterated maps, too (Feigenbaum, 1983; Collet, Eckmann & Lanford, 1980; Collet, Eckmann & Koch, 1981). Feigenbaum scaling is also observed with periodic solutions of the ordinary differential equations (Seydel, 2010). For the cascade on the cool side we may only estimate which is K.
4.3 Case studies 1: chaotic models
Before we describe the phenomena that shape the bifurcation diagram, we present a more detailed analysis for three chaotic models ( K, K and ). Their location in the bifurcation diagram (Fig. 8) is shown with black arrows. Nearly all of the computed chaotic models were analysed in the same manner as presented below.
There is not much to learn from the timeseries alone. As an example in Fig. 7 we show radius variation over pulsation cycles for K model, and indeed, no obvious regularity can be noticed. Much more useful are return maps for maximum radii, ie. plots of vs. . These are Poincaré maps with surface of section defined by at maximum expansion phase. The maps are shown in Fig. 9 for the discussed models. The points do not populate the plot in a random manner but fall along a characteristic, albeit rather complex shape, which evolves as the effective temperature changes. This evolution is illustrated with animation that may be found in the online version of this article as additional supporting information (see also other maps for chaotic models of different effective temperatures – grey dots in Figs. 11 and 13). The complex shape, as compared to analogous maps for some classical systems (e.g. tent map for the maximum values in the Lorenz system, see Lorenz (1963) or in Petigen, Jürgens & Saupe (2004)) is not surprising. Our system is much more complex and return map is only a 2D projection of complex dynamics occurring on a higher dimension manifold. The insets in Fig. 9 provide insight into the finestructure of the chaotic attractor. Although the numerical noise does not allow to show the cascade of such zoomins into smaller and smaller regions we conclude that the attractor’s structure is most likely fractal and the attractor is strange.
Chaos clearly manifests in the Fourier spectra, which are plotted in Fig. 10. In each case a time series for pulsation cycles was analysed with Period04 software (Lenz & Breger, 2005). The spectra were prewhitened with the frequency of the fundamental mode and its harmonics – location of these frequencies is shown with dashed lines. For two models also additional highestsignal frequencies in the resulting spectra were prewhitened and the result is also shown in Fig. 10.
In each case wide bands of signal are present in the spectra, which are characteristic signature of chaos. One large band is present for the K model (Fig. 10, top panel), without any obvious structure. To the contrary, in the case of the K model (Fig. 10, middle panel), wide bands are concentrated around and its harmonics. The cause of this difference becomes clear once the bifurcation diagram is analysed (Fig. 8). In the case of K model there are no obvious preferred values, or ranges of values, for the maximum radii. For the K model, the maximum radii fall preferentially into three ranges which manifest in Fig. 8 as three darkgrey bands within one large chaotic domain extending between K and K. The signal at and its harmonics is obviously not coherent which successive prewhitening (see Fig. 10) shows.
The situation is slightly different for the K model (Fig. 10, bottom panel). After prewhitening with the fundamental mode and its harmonics strong signal is present at subharmonic frequencies. It results from the twoband structure of the attractor clearly visible in the return map (Fig. 9, right panel). The maximum radii alternate between the two bands and vary chaotically within each of them. Therefore, signal at and its harmonics is strong and highly coherent but, after prewhitening with these frequencies, only very wide bands of signal without any obvious structure are present in the frequency spectra.
4.4 Periodic windows of order
Within the chaotic band (Fig. 8) seven windows, in which much more ordered behaviour is observed, including strictly periodic variation, may be identified. Below they are briefly described in the order they appear as effective temperature is increased. The return maps for some of the models are plotted in Figs. 11, 12 and 13. In each case the neighbouring chaotic model (of lower effective temperature) is plotted with grey points for reference. The most interesting period3 window is discussed in more detail in Section 4.5 and two interesting models from period6 and period9 windows are discussed in Section 4.6

Period9 window is present for models in a range K (Fig. 11). In fact this window is not strictly periodic but display a complicated internal structure. Within this window models were computed with a smaller 0.5 Kstep in effective temperature. For all models nine bands are clearly visible, which are either very wide ( K, K) or very narrow ( K) indicating a possible strict period9 cycle. In the case of K model each of the nine bands is split and thus 18 bands are apparent. A detailed analysis of the K model (that shows a rare case of typeIII intermittency) is presented in Section 4.6.

Period6 window is present for one model of K (Fig. 13a). The neighbouring K models display one chaotic band. Certainly our model grid lacks resolution to provide more insight into the dynamic scenarios within such a narrow window. In return map six very small bands rather than points are present, and thus model is not strictly periodic. The periodic model might be located just a tiny fraction of Kelvin away. This remark also applies to other very narrow windows discussed below.

Period5 window is present for one model of K (Fig. 13b). The neighbouring K models display one chaotic band. Detailed analysis shows that in fact for this model the maximum radii form five narrow chaotic bands which may join during the chaotic bursts for several pulsation cycles. Two such bursts happened within cycle integration of the model and are illustrated in Fig. 14. In this model we also detect a signature of typeIII intermittency: switching between period5 cycle and period10 cycle. Since the effect is barely visible for this model, we postpone its discussion to Section 4.6, in which we present a more clear example of typeIII intermittency in one model in period9 window.

Period3 window with subsequent period doubling cascade extends between K (period3 cycle) and K (three chaotic bands). This largest window is discussed in detail in the next Section.

Period6 window with subsequent (inverse) period doubling cascade. The window extends between K, at which three chaotic bands are present, till K at which period6 cycle, reached through the inverse period doubling cascade, ceases to exist. This window is in fact a mirror image of the just discussed period3 window, except the (inverse) period doubling cascade seems to be truncated at the cease of period6 cycle, which is followed by chaos instead of stable period3 cycle. The scenario around K is uncommon and will be discussed in more detail in Section 4.6 and in Section 6.

The latter two windows represent a clear example of period3 bubble (or remerging Feigenbaum tree, Bier & Bountis, 1984). The outlook at the bifurcation diagram (Fig. 8) reveals that the two just discussed windows are in fact tightly connected. The scenarios at the cool and hot side of the chaotic band separating the windows and extending between K are not only mutual mirror images. In fact the three chaotic bands formed at the two sides of the chaotic domain do not disappear as they merge into one chaotic band (in the interior crisis bifurcation, see next Section), but sustain their identity and smoothly merge within the chaotic domain as the darkgrey bands in Fig. 8 indicate.

Period6 window is present for one model of K (Fig. 13c). The neighbouring K models display two chaotic bands. The ones at the cool side merge into one band at K.
4.5 Intermittency and crisis at period3 window
In this Section we focus attention on the largest period3 window extending between K and K, and its direct vicinity. For the most interesting temperature ranges, at the edges of the window, the models were computed with a very small, Kstep in effective temperature. The corresponding part of the bifurcation diagram displays a striking similarity to the bifurcation diagram of the logistic map, Fig. 1, and bifurcation diagrams of many other dynamical systems (e.g. Rössler system, see in Petigen, Jürgens & Saupe, 2004). This similarity, and the analysis presented in this Section, lead to conclusion that in these systems the same dynamical phenomena lead to the appearance of the period3 window and its subsequent evolution to chaos.
In Fig. 15 (top two panels) we show the values of maximum radii during the pulsation cycles in two models directly preceding the appearance of the period3 window. Intermittency is clearly visible. For the K model the intervals with apparently periodic, period3 variation are short, last by up to 50 pulsation cycles, and are interrupted by much longer intervals of chaos. As the period3 window is approached the almost periodic behaviour dominates, as is the case for the slightly hotter model of K. The long intervals of period3 behaviour are only sporadically interrupted with much shorter bursts of chaos. This behaviour is characteristic for typeI intermittency. The appearance of intermittency and the birth of the period3 window are further illustrated in Fig. 16. It shows the third return maps (top) and small sections of radius variation for three different models at the onset of period3 window. The Figure is a counterpart of Fig. 4 for the logistic map. The intermittent behaviour is well visible in the models directly preceding the period3 window. The return maps show the formation of the intermittent channels (one is zoomed in the insets for K and K models), which get narrower as effective temperature is increased and bifurcation point is approached. As system evolves through the intermittent channels the apparently periodic variation is observed. For the hotter model, the density of points at the very narrow intermittent channel is high indicating long intervals of periodic variation. This is also illustrated with the sections of radius variation (bottom panels of Fig. 16) which were chosen to show the almost periodic variation interrupted with short chaotic burst (at d for K model). In between K and K the intermittent channel touches the diagonal and period3 cycle is born in the tangent bifurcation (rightmost panel in Fig. 16).
The ensuing scenario follows the one depicted schematically in the lower part of Fig. 3 and discussed in more detail for the logistic map (Section 2). As effective temperature is increased, the period3 cycle undergoes a series of period doubling bifurcations en route to chaos. The first bifurcation takes place at K and leads to period6 cycle. The model at K is a period12 cycle. Further period doublings are not resolved in our computations and the following models display a 6band chaos and finally 3band chaos.
The three chaotic bands merge in an interior crisis bifurcation which is illustrated with the help of Fig. 17, which should be analysed together with its counterpart for the logistic map – Fig. 5. For the first two models 3banded chaos is present – the values of maximum radii are bounded in three, well separated ranges, which is clearly visible both in the return maps and in the sections of radius variation. At the tangent bifurcation leading to the appearance of period3 window the unstable period3 cycle was also created. Its exact location in the return map cannot be computed as was the case for logistic map, but the corresponding three points must be located at the diagonal. As effective temperature is increased (cf. left and middle panels of Fig. 17) the three chaotic bands expand, approaching the diagonal and the anticipated unstable period3 cycle, location of which may now be easily guessed. At the interior crisis bifurcation the threeband chaotic attractor collides with the unstable period3 cycle and expands into oneband chaotic attractor.
Directly after the crisis bifurcation, the maximum radii still fall preferentially into the ranges defined by the three chaotic bands. It is well visible in the return map in the rightmost panel of Fig. 17 (higher density of points along the three extended arches) and in the corresponding section of radius variation. In addition in Fig. 15 (bottom panel) we plot the values of maximum radii during the consecutive pulsation cycles for the same model. For many pulsation cycles the system evolves in a phasespace defined by the former threeband chaotic attractor and only sporadically gets scattered over a larger space. This behaviour is called crisis induced intermittency and ceases as effective temperature is increased, albeit, still the probability of maximum radii falling into one of the three bands is larger thorough the full chaotic domain separating the period3 and period6 windows ( K; see Fig. 8).
4.6 Case studies 2: periodic models
In this Section we discuss in more detail two models displaying periodic variation: K and K models. In the bifurcation diagram (Fig. 8) their location is indicated with green arrows. The first model is located within period6 window while the second is located within period9 window. Figs. 18 and 19 show the frequency spectra for the two models. For clarity, these are limited to range.
As arrow in Fig. 8 indicates K model is located at the border of chaotic band and period6 window. Whether it is in fact a period3 cycle variation, just a moment before the period doubling bifurcation, or period6 cycle, just after the bifurcation, is not clear from the bifurcation diagram. Analysis of return map, which consists of three slightly elongated clumps extending over less than does not provide the clue, either. Before we discuss the Fourier spectrum for the considered model, we briefly describe other method that may be used to resolve the issue, which we find particularly useful for similar models at the direct vicinity of bifurcation points. We analyse growth of the maximum kinetic energy of the model, from one pulsation cycle () to other (). The corresponding kinetic energy growth rate, which we define as , summed over pulsation cycles, , should be close to zero for the period limit cycle. The method has thus additional advantage of pointing whether variations converge to limit cycle, or not. In the bottom panel of Fig. 20 we plot both and for the considered model. For a reference, in the top panel we plot for a singlyperiodic K model located well beyond the chaotic domain. The mean value of for this model is zero, as expected. The small () fluctuations around the mean result partly from the interpolation (the time step is not an integer part of the period) and partly from the numerical noise which is relatively high in our models with low eddy viscosity (see Section 3). Analysis of the bottom panel of Fig. 20 clearly points that the K model display period6 cycle behavior rather than period3. This is further supported with the Fourier spectrum (Fig. 18). After prewhitening with the fundamental mode and its harmonics (top panel of Fig. 18) a very strong signal is present at and its harmonics. After prewhitening we clearly detect signal at (and harmonics), which is however, three order of magnitude weaker (middle panel of Fig. 18). After the next prewhitening, signal is still present, which indicates that period6 cycle is not strictly periodic, but displays irregular variation with very small amplitude. We note that only a fraction of Kelvin away chaotic band extends.
The other model ( K) is not strictly periodic either, which is clear already from the bifurcation diagram and first return map (Fig. 11, brown points), which consists of nine separate and extended bands. Accordingly, strong signal at (and harmonics) is present (top panel in Fig. 19). After prewhitening, strong signal at and its harmonics is present, which is not coherent however, as successive step of prewhitening shows (bottom panel of Fig. 19). Fig. 21, showing the maximum radii over the consecutive pulsation cycles, explains the origin of the signal. The model switches intermittently between period9 and period18 cycle: starting from period9 cycle the subharmonic cycle grows up to some amplitude and then model switches back to cycle9 behaviour. This is typeIII intermittency (see Pomeau & Mannevile, 1980, and Section 2), a result of collision between stable period9 cycle and unstable perioddoubled (period18) cycle arising from subcritical period doubling bifurcation. TypeIII intermittency is rarely observed in dynamical systems and in most cases concerns switching between period1 cycle and period2 cycle (e.g. Dubois, Rubio & Berge, 1983). The very similar behaviour as in the case of K model was observed in electric circuits (see Thamilmaran & Lakshmanan, 2002, their fig. 21).
5 Largest Lyapunov exponents
One of the key features of chaotic system is its sensitivity to initial conditions, which can be quantified using Lyapunov exponents. For the chaotic attractor, the two initially nearby trajectories diverge at an exponential rate given by the largest Lyapunov exponent, . To determine for our models we used the code and algorithm developed by Rosenstein, Collins & De Luca (1993) based on the original algorithm of Grassberger & Procaccia (1983). Below we briefly describe the underlying idea and refer the reader to original papers for details.
We first reconstruct the attractor dynamics using the method of delays. We use the radius values equally spaced in time, , to construct delay vectors, , in the embedding space:
is the embedding dimension and is time delay (or lag). The algorithm finds the close neighbours in the phase space, two points and with sufficiently small distance between, . For chaotic system the distance should grow exponentially in time, so that , corresponding to the largest Lyapunov exponent. Analysis of results for many pairs of close neighbours leads to a robust estimate of .
For each model we have analysed only a section of the computed time series, of length , expressed below as a multiple of the fundamental mode pulsation period. The default length of the analysed data is pulsation cycles. The embedding dimension cannot be smaller than the physical dimension of the attractor. Provided is high enough, its exact value should not affect the determined values of Lyapunov exponents. We find no systematic differences in the largest Lyapunov exponents computed assuming and , and pick the latter value as the default (see also next paragraph). For the time delay we follow Rosenstein, Collins & De Luca (1993) and use for which autocorrelation function drops to of its initial value.
Before running computations for all models we checked the sensitivity of to chosen values of , and , for one chaotic model of K. For the chosen dense time sampling (d, more than points per pulsation cycle) the lag resulting from autocorrelation function is . In the top row of Tab. 1 we show computed assuming default values of , and . Next, we have varied the values of these parameters and report the resulting largest Lyapunov exponents in the following rows of Tab. 1. We find that determination of is robust and does not depend on exact values of the discussed parameters, provided that embedding dimension and lag are sufficiently high. Also, the chosen default length of the analysed timeseries is sufficient, which we further verified computing for all chaotic models, assuming the default values for and , and two different lengths of timeseries, the default one, of pulsation cycles, and the shorter, of pulsation cycles. The results are plotted in the bottom part of the bifurcation diagram displayed in Fig. 8, with thick black and thin blue lines respectively. Only for a few models a noticeable difference is present. We conclude that with pulsation cycles the attractor dynamics is well probed.
[d]  
The positive values of the largest Lyapunov exponents clearly establish the chaotic nature of our models. The typical values of (Fig. 8) vary between d and d. The largest values are slightly above d and are calculated for models between K and K. Within chaotic bands variation of with effective temperature is not smooth. Large drops of towards zero are seen, as expected, at the edges of periodic windows (for strictly periodic model ). Smaller drops within chaotic domain may indicate a nearby periodic window, which was not detected because of too coarse resolution in effective temperature.
It is interesting to compare the computed values of to those determined for two RV Tau stars and one Miratype variable. Results obtained for these stars are collected in Tab. 2 together with determination for one radiative chaotic model of Kovács & Buchler (1988). Ranges of rather than single values are given as determination of from real stellar data is much more sensitive to the values of embedding dimension, lag, etc. The reader is referred to tables in original papers (references are given below Table) for details.
The values of largest Lyapunov exponent determined for our models are typically two orders of magnitude larger than values determined for stars or a radiative model listed in Tab. 2. On the other hand, the pulsation periods of our models (always between d) are order or two orders of magnitude shorter (see third column in Tab. 2). In the context of pulsating stars however, it is more appropriate to use the length of a single pulsation cycle (pulsation period) as a unit of time. When expressed in units of inverse pulsation cycles rather than d, the values of largest Lyapunov exponents of our models are comparable to the values determined for RV Tau/Miratype stars. They are larger than for RV Tau type stars, but a factor smaller than for Mira type star.
star/type  [ d]  period [d]  Ref. 

R Cyg / Miratype  430  1  
R Sct / RV Tau  70  2  
AC Her / RV Tau  37.5  3  
D5200 / model  11.5  4 
References: (1) Kiss & Szatmáry (2003) (tab. 2), (2) Buchler, Kolláth & Serre (1996) (tab. 1), (3) Kolláth et al. (1998) (tab. 1), (4) Serre, Kolláth & Buchler (1996) (tab. 1).
6 Discussion
The nonlinear stellar pulsation equations we solve, form a much more complex system than classical chaotic systems discussed in the literature. Yet the resulting dynamical scenario is qualitatively similar to that arising from the iteration of the simplest logistic map (compare the bifurcation diagrams in Figs. 1 and 8 and animations attached with the online version of this paper as supporting information). Most of conclusions about the origin of dynamical phenomena found in our models are drawn based on the analogies between our models and simpler systems for which strict analytical reasoning is possible. However, for some of the phenomena we detect, we do not find a satisfactory analogy. The appearance of period6 window, extending between K and K, is one of them.
The scenario that is expected and that is encountered in many other systems, and is also present in the period3 window (between K and K) is the following. First, a stable period3 cycle is born together with unstable period3 cycle. The stable branch undergoes a series of period doubling bifurcations to form three chaotic bands, which finally collide with the unstable period3 cycle to form one chaotic band (Section 4.5). To the contrary, in the window at K, a stable period6 cycle, which looks like two stable period3 cycles born very close to each other, emerges from the chaos. We are not aware of any bifurcation that may lead to such scenario and of any other system showing such behaviour.
One of the possibilities is that in fact a pair of stable and unstable period3 cycles is born, as expected, and stable cycle immediately undergoes a period doubling bifurcation. To check this, we have computed additional models in the interesting temperature range (with Kstep in effective temperature), but they display either chaos or a period6 behaviour (in Section 4.6 we analysed one of these models). If the proposed scenario indeed takes place it must occur in a temperature range narrower than K.
The situation at K looks even more complex. It seems that at this temperature we deal with discontinuity – the bifurcation diagram (Fig. 8) divides into two parts apparently decoupled from each other. On the cool side we clearly see a gradual evolution of the chaotic bands, divided, from time to time, by periodic windows. This gradual evolution seems to continue till K, but not for the band extending at higher effective temperatures. This behaviour may result from the coexistence of two attractors in the system. Note that by default we initialized all the models along a sequence in the same manner (Section 3). It is possible that such initialization leads to different attractors for models cooler/hotter than K. To check the possible existence of other stable attractor(s) (with different basins of attraction), we have repeated the computation for many models, but with several different initializations. In all cases however, we finally arrived at the same attractor as in the case of our default initialization. At the moment the cause and nature of bifurcation we observe at K remains unclear.
The other phenomenon we have not discussed yet, is the appearance of chaos itself. The chaotic bands appear through a well understood period doubling route; the question is about the trigger. In the case of radiative models of Buchler & Kovács (1987) and Kovács & Buchler (1988), analysis of the Floquet coefficients clearly shows that the first period doubling bifurcation is caused by the halfinteger resonance between pulsation modes ( between fundamental and second overtone modes; Moskalik & Buchler, 1990). The following cascade en route to chaos was not analysed. Moskalik & Buchler (1990) analysed a toy model of parametrically driven oscillator and showed that the first period doubling in such system results from the resonance and the following cascade is a result of increasing nonlinearity. We note that the nonlinearity may be the only cause of period doublings in classic chaotic systems void of internal resonances.
We do not have appropriate tools (Floquet coefficients) to proof that halfinteger resonance is responsible for the period doubling of singleperiodic pulsation we observe at the cool and the hot sides of our computation domain. The closest resonances are (Fig. 6) on the hot side of the computation domain and, and on the cool side. The loci of the resonance, causing the period doubled pulsation detected in a single BL Her star (Smolec et al., 2012) is located more than 300 Koff the hot side of the computation domain considered here and likely plays no role. We cannot exclude the possibility that nonlinearity is the only cause of the observed behaviours.
The appearance of periodic windows is also very interesting and important for stellar pulsation studies. As in the case of perioddoubled pulsation, resonances were also invoked as a possible explanation. In a recent study Plachy, Kolláth & Molnár (2013) proposed that a resonance between the fundamental mode and the first overtone is responsible for a period20 cycle behaviour they found in one of their RR Lyrae models (their model H; for other model they propose a resonance). In this case however, a caution is needed, as inferences about the role of resonance are not based on firm theoretical grounds. The presence of first overtone cannot be deduced from the frequency spectrum. The role of resonance is most likely^{2}^{2}2Plachy, Kolláth & Molnár (2013) do not discuss in detail how the connection between the periodic pulsation and resonances is made. deduced based on approximate coincidence of the model’s location in the HR diagram with the loci of the resonance determined with linear pulsation periods. Since pulsation periods change in the nonlinear regime, and finetuning of such highorder resonance is difficult, claims on the possible role of resonances must be supported with other (dynamical) arguments. We are not aware of any studies showing that such highorder resonances may indeed have a noticeable effect on stellar pulsations.
In a simpler explanation, period behaviour detected in the models, is an intrinsic property of nonlinear, chaotic system. The K neighbours of the discussed model of Plachy, Kolláth & Molnár (2013) show chaos and thus situation corresponds to periodic window within chaotic band, just as we report in this paper (Section 4.4), and as is found in many chaotic systems void of resonances. In chaotic systems the spectrum of periodic windows is dense (it is one of the key properties of chaotic systems; e.g. Petigen, Jürgens & Saupe, 2004), but most of the windows are extremely narrow. Based on extreme similarity of our bifurcation diagram (Fig. 8), to bifurcation diagrams for other systems, we conclude that also in the case of our computations the spectrum of periodic windows is most likely dense, but most of the windows are extremely narrow (in effective temperature). With the default K resolution of our model computation only few of the windows were detected (and most of them are narrower than K). Studying linear periods we find no tight connection between location of the periodic windows and the loci of highorder resonances between low order pulsation modes. We conclude that existence of periodic windows is an intrinsic property of nonlinear system studied in this paper. There is no need to invoke resonances to explain them.
7 Observability of the chaotic phenomena
Phenomena discussed in this paper are certainly exotic in the context of BL Her stars. Except of period doubling phenomenon detected in one star (Soszyński et al., 2011; Smolec et al., 2012), BL Her stars are singleperiodic pulsators. We recall however, that existence of perioddoubled BL Her star was predicted years before its discovery (Buchler & Moskalik, 1992). The future detection of other, more complex dynamical phenomena in these stars cannot be excluded.
Period doubling and irregular brightness variations are common in higher luminosity siblings of BL Her stars, namely in more luminous type II Cepheids of RV Tau type and in semiregular variables. Period doubling is a characteristic feature of RV Tau stars which, in addition, show irregular variation. Chaos was reported in several stars of these types (Buchler, Kolláth & Serre, 1996; Kolláth et al., 1998; Buchler, Kolláth & Cadmus, 2004; Kiss & Szatmáry, 2003). Our models suggest that other phenomena, intrinsic to chaotic dynamics, may also occur in these groups of stars. Particularly interesting would be the discovery of period behaviour (with other than 2) and of intermittency.
Period behaviour may arise either within a period doubling cascade or within chaotic band. A possible period4 behaviour in RV Tau type star, supporting the perioddoubling transition to chaos scenario, was reported by Pollard et al. (2000). Unfortunately, our unpublished analysis of the OGLEIII data for this star does not confirm the detection. Models indicate that in the case of period doubling cascade the domains of period behaviour get narrower as increases, making the detection less probable for larger . Also most of the periodic windows within chaotic domain are very narrow. Nevertheless, in our opinion, the firm detection of period4 behaviour in RV Tau stars is only a matter of time.
From observers point of view, the strongest evidence for period behaviour would be the presence of additional signal in the frequency spectrum, at and its harmonics. The difficulty arises for longperiod variables, as very long timeseries is necessary to get sufficient resolution in the subharmonic part of the frequency spectrum (i.e. in a range ). In addition, the effect may be obscured by irregular variability on top of period behaviour, as is commonly observed in RV Tau stars (see eg. lightcurves in Soszyński et al., 2011). In these stars, a sporadic switching of the deep and shallow minima also occurs (Wallerstein, 2002) which is yet another difficulty in the analysis. Therefore, inspection of the timeseries, folding the light curve with multiple of the basic period, are invaluable tools to search for the effect. This however requires not only long timeseries, but also well sampled timeseries. In this respect the projects aimed at longterm and regular monitoring of the sky or individual stars, such as OGLE (Udalski et al., 2008), EROS (e.g. Beaulieu et al., 1995), ASAS (Pojmanski, 2002), or programs led by amateur associations of variable star observers (e.g. AAVSO) are extremely important and should be continued as long as possible. It is worth noting that the observations of the only stars rigorously analysed for the presence of chaotic dynamic (i.e. those from Tab. 2) were all collected by amateur astronomers and in all cases covered more than 10 years (up to a century for R Cyg).
Based on our rather restricted model survey we cannot predict the expected amplitude of brightness alternations within period cycle. In perioddoubled RV Tau stars the amplitude of alternations vary from a hundredth of magnitude to few tenths of magnitude. More precise the photometry, larger the chance to detect the smaller effect.
Longlasting and frequent observations of the stars are also crucial for the possible detection of intermittency. Models indicate that the intervals of apparently periodic/chaotic behaviour may last many tenths of pulsation cycles which, as pulsation periods of the luminous stars are long, requires a very long monitoring. Again, intermittency is expected in the narrow domains close to the edges of periodic windows and thus probability of its detection is certainly very small.
8 Summary
The BL Her models discussed in this paper fall along a single stripe of constant luminosity in the HR diagram and cover a range of only K. Yet they display a wealth of dynamical behaviours characteristic for deterministic chaos. Many of the discussed phenomena are detected for the first time in the context of stellar pulsation models. It was possible because our model survey was dedicated to study such phenomena – a tiny step in effective temperature, sometimes as small as K (and K max), allowed to follow the dynamical evolution of the system from singleperiodic pulsation, through period doubling cascade to well developed chaotic regime, and back to singleperiodic pulsation. The chaotic regime turned out to be a goldmine of interesting dynamical phenomena. We found several periodic windows (with cycle3, 5, 6, 7 and 9 behaviours). We stress that the existence of periodic windows is not related to resonances among pulsation modes, but is an intrinsic property of a chaotic system. At the edges of the largest period3 and period6 windows we have found intermittent behaviour and crises bifurcations. Particularly interesting is intermittency – a sporadic switching between two qualitatively different behaviours. In typeI intermittency intervals of apparently periodic (period, in general) behaviour are interrupted with bursts of chaos. In typeIII intermittency the oscillations switch between two periodic cycles. In our models we detected switching between period9 and period18 cycles.
Detection of the discussed phenomena in the stars would be extremely interesting, however it requires a long and regularly sampled time series. The already available data from projects such as OGLE offer the best opportunity for a successful search.
Acknowledgments
Model computations presented in this paper were conducted on the psk computer cluster in the Copernicus Centre, Warsaw, Poland. This research is supported by the Polish Ministry of Science and Higher Education through Iuventus+ grant (IP2012 036572) awarded to RS. PM is supported by the Polish National Science Centre through grant DEC2012/05/B/ST9/03932.
References
 Beaulieu et al. (1995) Beaulieu J.P., et al., 1995, A&A, 303, 137
 Becker et al. (1999) Becker J., Rödelsperger F., Weyrauch Th., Benner H., Just W., Cenys A., 1999, Phys. Rev. E, 59, 1622
 Bier & Bountis (1984) Bier M., Bountis T.C., Phys. Lett. A, 104, 239
 Buchler & Kovács (1987) Buchler J.R., Kovács G., 1987, ApJL, 320, L57
 Buchler & Moskalik (1992) Buchler J.R., Moskalik P., 1992, ApJ, 391, 736
 Buchler & Kolláth (2011) Buchler J.R., Kolláth Z., 2011, ApJ Lett., 731, 24
 Buchler, Kolláth & Serre (1996) Buchler J.R., Kolláth Z., Serre T., 1996, ApJ, 462, 489
 Buchler, Kolláth & Cadmus (2004) Buchler J.R., Kolláth Z., Cadmus R.R., 2004, ApJ, 613, 532
 Collet, Eckmann & Lanford (1980) Collet P., Eckmann J.P., Lanfold III O.E., 1980, Comm. Math. Phys., 76, 211
 Collet, Eckmann & Koch (1981) Collet P., Eckmann J.P., Koch H., 1981, J. Stat. Phys., 25, 1
 Dubois, Rubio & Berge (1983) Dubois M., Rubio M.A., Berge P., 1983, Phys. Rev. Lett., 51, 1446
 Eckmann (1981) Eckmann J.P., 1981, Rev. Modern Phys., 53, 643
 Feigenbaum (1983) Feigenbaum M.J. 1983, Physica D, 7, 16
 Grassberger & Procaccia (1983) Grassberger P., Procaccia I. 1983, Physica D, 9, 189
 Grebogi, Ott & Yorke (1982) Grebogi C., Ott E., Yorke J.A. 1982, Phys. Rev. Lett., 48, 1507
 Hilborn (2000) Hilborn R.C., 2000, “Chaos and Nonlinear Dynamics. An Introduction for Scientists and Engineers”, second Edition, Oxford University Press
 Kiss & Szatmáry (2003) Kiss L.L., Szatmáry K., 2003, A&A, 390, 585
 Kolláth et al. (1998) Kolláth Z., Buchler J.R., Serre T., Mattei J., 1998, A&A, 329, 147
 Kovács & Buchler (1988) Kovács G., Buchler J.R., 1988, ApJ, 334, 971
 Kuhfuß (1986) Kuhfuß R., 1986, A&A, 160, 116
 Lenz & Breger (2005) Lenz P., Breger M., 2005, Comm. in Asteroseismology, 146, 53
 Letellier et al. (1996) Letellier C., Gouesbet G., Soufi F., Buchler J.R., Kolláth Z., 1996, Chaos, 6, 466
 Lorenz (1963) Lorenz E.N, 1963, J. of the Atm. Sci., 20, 130
 Mannevile & Pomeau (1979) Mannevile P., Pomeau Y., 1979, Phys. Lett. A, 75, 1
 May (1976) May R.M., 1976, Nature, 261, 459
 Moskalik & Buchler (1990) Moskalik P., Bucher J.R., 1990, ApJ, 355, 590
 Plachy, Kolláth & Molnár (2013) Plachy E., Kolláth Z., Molnár L., 2013, MNRAS, 433, 3590
 Pojmanski (2002) Pojmanski, G. 2002, Acta Astron., 52, 397
 Pollard et al. (2000) Pollard K.R., et al., 2000, ASP Conf. Ser., 203, 89
 Pomeau & Mannevile (1980) Pomeau Y., Mannevile P., 1980, Comm. Math. Phys., 74, 189
 Petigen, Jürgens & Saupe (2004) Petigen H.O., Jürgens H., Saupe D, 2004, “Chaos and Fractals”, secon edition, SpringerVerlag, New York
 Rosenstein, Collins & De Luca (1993) Rosenstein M.T., Collins J.J., DeLuca C.J. 1993, Physica D, 65, 117
 Serre, Kolláth & Buchler (1996) Serre T., Kolláth Z., Buchler J.R., 1996, A&A, 311, 845
 Seydel (2010) Seydel R., 2010, “Practical Bifurcation and Stability Analysis”, third edition, Springer, New York
 Smolec & Moskalik (2008) Smolec R., Moskalik P., 2008, Acta Astron., 58, 193
 Smolec & Moskalik (2012) Smolec R., Moskalik P., 2012, MNRAS, 426, 108
 Smolec & Moskalik (2014) Smolec R., Moskalik P., 2014, IAU Symp. No 301, p. 489
 Smolec et al. (2012) Smolec R., et al., 2012, MNRAS, 419, 2407
 Soszyński et al. (2011) Soszyński I., et al., 2011, Acta Astron., 61, 285
 Thamilmaran & Lakshmanan (2002) Thamilmaran K., Lakshmanan M., 2002, Int. J. of Bifurcation and Chaos, 12, 783
 Udalski et al. (2008) Udalski A., Szymanski M., Soszynski I., Poleski R., 2008, Acta Astron., 58, 69
 Wallerstein (2002) Wallerstein G., 2002, PASP, 114, 689
SUPPORTING INFORMATION
Additional Supporting Information may be found in the online version of this article:
Animation 1: The animation shows the evolution of first return map for logistic equation as control parameter is increased from to . Note that smaller step is used within chaotic regime.
Animation 2: The animation shows the evolution of first return map for BL Her models discussed in the paper. Note the different, smaller effective temperature step at the edges of period3 window.