Meditation-Induced Coherence and Crucial Events

In this paper we emphasize that 1/f noise has two different origins, one compatible with Laplace determinism and one determined by unpredictable crucial events. The dynamics of heartbeats, manifest as heart rate variability (HRV) time series, are determined by the joint action of these different memory sources with meditation turning the Laplace memory into a strongly coherent process while exerting an action on the crucial events favoring the transition from the condition of ideal 1/f noise to the Gaussian basin of attraction. This theoretical development affords a method of statistical analysis that establishes a quantitative approach to the evaluation of the stress reduction realized by the practice of Chi meditation and Kundalini Yoga.


INTRODUCTION
The phrase normal sinus rhythm is known to be an overly simplified characterization of the true rhythmic pattern made by the series of time intervals between successive human heartbeats (Peng et al., 1993). Instead of a steady regular time interval between beats, healthy hearts manifest a heart rate variability (HRV), whose time series has a rather broad frequency spectrum. Historically, the frequency spectrum was explained to be a consequence of the firing rate of the sinus node (the heart's natural pacemaker) being controlled by signals from the autonomic (involuntary) portion of the nervous system, which has two major branches: the parasympathetic (decreases firing rate) and the sympathetic (increases firing rate). It is this continuing tug-of-war on the sinus node, with one decreasing and the other increasing the firing rate of the sinus node pacemaker cells that produces the HRV spectrum in healthy individuals. A number of techniques have been used to establish that HRV time series are fractal, for example, see Ivanov et al. (1999); Stanley et al. (1999), and consequently the frequency spectrum scales, for a review see West (2013). Herein we present a picture of the source of the HRV spectrum based on the statistics of the HRV time series.
Stress is known to disrupt the normal HRV spectrum and learning how to circumvent the debilitating influence of stress is one of the purposes of meditation. The origin of meditation is unknown, but scholars agree that it has been part of civilized culture for thousands of years, with the earliest records, circa 1,500 BCE, involving Vedantism, a Hindu tradition in India. As mentioned by Chow (2015) the Yoga Sutras of Patanjali, outlining the eight limbs of yoga, was compiled between 400 and 100 BCE. Cotemporaniously, the Bhagavad Gita was written, which discusses the philosophy of yoga and meditation.
From a physiological perspective meditation constitutes a coupling of the functionality of the heart with that of the brain and has been explored and developed for millennia. However, we are only now beginning to apply the methods of science and data analysis to this brain-heart coupling. The main purpose of the present paper is to provide a measure of reduction in the level of stress provided by meditation, which is to quantify, through the statistical analysis of HRV time series before and during meditation, precisely how much stress is alleviated by control through meditation of the heart-brain coupling.

Statistics and Memory
One of the most difficult issues confronting complexity science is the origin of 1/f -noise (Watkins, 2016(Watkins, , 2017. Is it an ergodic or a non-ergodic process? Watkins (2017) has recently pointed out that Mandelbrot, very well known for his generalization of ordinary diffusion, called Fractional Brownian Motion (FBM) (Mandelbrot and Van Ness, 1968), which has a physical origin compatible with stationarity and ergodicity, is also the author of papers (Mandelbrot, 1965(Mandelbrot, , 1967 opening a bridge between the stationary and non-stationary condition. Actually, the revisitation of the work done by Mandelbrot in the 1963Mandelbrot in the -1967 period of time leads to the creation of a connection with the Continuous Time Random Walk (CTRW) (Montroll and Weiss, 1965;Shlesinger, 2017). This is equivalent to establishing a connection with the subordination perspective (Sokolov, 2000;Thiel and Sokolov, 2017) and with the assumption that crucial events determine the observed dynamics. The crucial events are characterized by the main property that the time interval between the occurrence of consecutive crucial events is described by the waiting-time probability density function (PDF): According to Watkins (2016Watkins ( , 2017) the non-ergodic fractional renewal model of Mandelbrot (1965Mandelbrot ( , 1967 yields the inverse power law (IPL) spectrum as a function of the frequency f : The dynamics generating crucial events are renewal, which means that the sequence of waiting times {τ i } between succsive events are completely independent of one another. This important physical property is expressed by the mathematical property of the two-point correlation function: where δ i,j denotes the Kronecker function; = 1 for i = j and = 0 otherwise, in the ideal limit τ 2 << τ 2 . In this paper we focus on C(1) defined as C(1) = C(t = |i − j| = 1). The physical interpretation of this mathematical property is that the occurrence of a renewal crucial event generates a rejuvenation of the system. Consequently, the system evolves toward the occurrence of the next renewal event as if the occurrence of the earlier crucial event made the system brand new. The renewal property facilitates resolving the 1/f − paradox (Watkins, 2016(Watkins, , 2017. This paradox has been brilliantly discussed by Niemann et al. (2013) who noticed that for γ > 1(µ < 2) the spectrum S(f ) is not integrable, thereby violating the integrability condition on which the Wiener-Khinchine theorem is based. This theorem is the mathematical foundation necessary for the evaluation of a spectrum. This theorem has been sidestepped by working with a time series of length L to obtain (Lukovic and Grigolini, 2008): This result is based on the observation (Feller, 1971) that the rate R(t) of production of crucial events, all prepared at t = 0, is given by the IPL: Actually, by idenifying L in Equation (5) with 1/f , we obtain S(f ) = 1/f , which is interpreted as an indication that working with a time series of length L makes f = 1/L the smallest observed frequency, thereby settling the paradox with no need of truncating the frequency distribution, insofar as the empirical truncation is a dynamical effect of the experimental observation.
It is interesting to notice that the restriction of the IPL index to the interval 2 < µ < 3 makes S(f ) integrable, thereby generating a sort of equivalence between this process and that of FBM. To make this equivalence explicit let us focus our attention on the Langevin equation:ẋ where ξ (t) is a correlated fluctuation adopted to make the diffusion process x identical to FBM. The rate equation thereby provides a dynamical origin for FBM (Cakir et al., 2006;Tuladhar et al., 2017). This approach is, in fact, based on a rigorous Hamiltonian model for a Generalized Langevin Equation (GLE), which can be made equivalent to the generator of Fractional Gaussian Noise (FGN) used by Kou and Xie (2004) to derive a diffusive process equivalent to FBM. FBM is by construction stationary throughout the whole time region. On the other hand, the crucial events responsible for the IPLs denoted by Equations (2) and (3) generate a process that is not stationary, but one that becomes stationary in the long-time limit. Let us imagine that the time interval between consecutive crucial events is filled with either ξ = 1 or ξ = −1, as determined by a coin-tossing prescription. When µ > 2 the waiting-time PDF is used to calculate the non-stationary correlation function <ξ (t 1 )ξ (t 2 )>, which becomes stationary in the long-time limit. In fact, the coin-tossing prescription makes the correlation function of age t a , < ξ (t a )ξ (t a + τ ) >, identical to the probability that no new event occurs up to the time τ after the occurrence of an event at time t a . The brand new survival probability is proportional to 1/τ µ−1 , but for t a → ∞ the survival probability becomes proportional to 1/τ µ−2 . Although the spectrum of these fluctuations is identical to that of FBM, with scaling: in the FBM case the regression of the correlated noise ξ (t) to the origin is characterized by an exponential waiting-time distribution. In addition the corresponding function C(i, j) violates the renewal condition of Equation (4), that is, it is not a Kronecker delta function. The correlation function of FGN, as well as, the correlation function of fluctuations generated by crucial events may be extremely slow, even if S(f ) with γ < 1 does not require the settlement of any paradox, insofar as 1/f γ is integrable. There is evidence that the FBM long-time memory may be realized in many complex systems, for instance, climate, hydrology and finance, see the excellent review paper (Graves et al., 2017) to oversee the importance of this form of long memory.
Note that the FBM memory can be derived from Mori's GLE (Mori, 1965), which, in turn, is derived from Hamiltonian dynamics by projecting from the evolution of the entire Universe down to the dynamic variables of interest. This projection determines the system's memory in the sense that the intensity of the correlation function ξ (t) of a variable of interest tends to vanish for t → ∞, but the decay is sufficiently slow that the correlation function ξ (t) is not integrable. This form of memory, therefore, is a generalization of the well known Laplace determinism (Laplace, 1951), according to which given the initial conditions for all the atoms in the Universe one would be able to predict with no error the state of the Universe at any time in the future.
Herein we refer to this form of memory, including Laplace determinism, as Hamiltonian Memory. The deterministic motion given by the rate equation: is referred to as Hamiltonian memory for real . By way of contrast we refer to the non-integrable correlation function generated by the crucial event fluctuations as Crucial Memory.
How can one distinguish between these two forms of memory? The cross-correlation function of Equation (4) plays an important role in answering this question. In fact, in the case of crucial event infinite memory C(1) = 0, whereas in the case of FBM memory C(1) = 0.
Is it possible that in real systems both forms of memory coexist?
Research Bohara et al., 2017) reveals that in the case of HRV time series both forms of memory do, in fact, coexist. Although this coexistence is to some extent as paradoxical as that of wave-particle dynamics in quantum mechanics, it is not quite as surprising in biology and in the field of physiological processes. In the case of brain dynamics the concept of nested frequencies (He et al., 2010) was adopted to establish a connection between the physiological wave-like nature and the action of crucial events that are revealed by the proper method of analysis (Allegrini et al., 2015). It has to be stressed that Peng et al. (1999) revealed that meditation has the effect of generating a surprising coherent behavior in HRV time series. The connection between coherence and criticality is a subject of great interest (Termsaithong et al., 2012) that we address herein using the subordination perspective (Sokolov, 2000;Thiel and Sokolov, 2017). The connection between meditation and coherence is a subject studied by many authors moving from the adoption of positive emotions (Bradley et al., 2010) to different forms of meditation (Peng et al., , 2004. Herein we establish, with a statistical analysis of data using an approach developed in earlier work Bohara et al., 2017) that the dynamics of heartbeats result in HRV time series hosting both Hamiltonian Memory and Crucial Memory. We also prove that meditation affects both kinds of memory, turning the Hamiltonian Memory into a coherent process . Meditation also has the surprising effect of affecting Crucial Memory, as well. In fact, our analysis shows that meditation shifts the IPL index µ of Equation (1) moving it from values very close to the perfect 1/f condition, µ = 2, to values near the border with the Gaussian basin of attraction, µ = 3. The meditation-induced coherence decreases the accuracy of the method adopted in earlier work Bohara et al., 2017) to evaluate the percent of stress indicators, the Poisson events thought to be responsible for heart failure, but refine the methodology and establish that both Chi and Kundalini Yoga meditation have the effect of reducing stress, although the advanced practice of Kundalini Yoga meditation yields higher levels of reduction.

SUBORDINATION OF HARMONIC OSCILLATIONS
The subordination of a dynamic process to harmonic motion was studied by Ascolani et al. (2009). We assume that the clock's hour-hand completes a rotation from noon to noon through M clicks. The time interval between two consecutive clicks t is normalized to 1. The frequency of the regular oscillator is given by: The regular oscillation is a deterministic Hamiltonian process equivalent to the process described by Equation (9). Being cognizant of the choice t = 1, we set ≪ 1 in order for the dynamics to be as close as possible to the coherent behavior of a regular clock. Each time interval between consecutive clicks is converted into an operational time τ , whose statistics are determined by the hyperbolic PDF: which is equivalent to realizing the time interval between crucial events being given by Equation (1). Numerically these statistics are produced (Lukovic and Grigolini, 2008) by drawing a random number y with a uniform PDF on the interval [0, 1] (0 < y < 1) and converting it into the operational time using: Each realization of this procedure yields a stochastic trajectory that for the clarity of discussion is here illustrated in Figure 1. Note that we focus our attention on µ > 2, a condition making the mean value τ finite, where the brackets denote an average taken using Equation (11): The number of crucial events activated moving from noon to the long-time limit t is equal to t/ < τ >. Consequently, the number of cycles of the subordinated signal is equal to that of the coherent oscillator.
It is important to stress that, according the theoretical arguments using diffusion entropy analysis (DEA) Bohara et al., 2017), in the long-time limit the crucial events generate a scaling process with index: which is numerically close to, but different from, the FBM scaling with an index H given by Equation (8).
The evaluation of the spectrum requires averaging the solution of Equation (9) over an ensemble of realizations. For this reason we define: where ξ (t) is a single realization and the brackets denote an ensemble average. The subordinated process Y(t) can be derived from the following prescription (Ascolani et al., 2009): where ψ n (t) denotes the waiting-time PDF of a sequence of n events, the last of which occurs at time t, and (t) is the probability that no event occurs until time t far from the occurrence of the last earlier event. The evaluation of Y(t) is described in the theoretical work (Lambert et al., unpublished), which also affords rigorous mathematical arguments to evaluate the spectrum of this process. Herein we limit ourselves to adopting heuristic arguments supported by numerical results. We evaluate the Laplace transform of the subordinated harmonic process Y(t), using the renewal character of the waiting-time PDF: where the Laplace transform of f (t) is denoted f (u) . Inserting Equation (17) into the Laplace transform of Equation (16) and summing the series yields: Note that in the long-time limit, namely for u → 0, this expression reduces to: Note that in the case where the waiting-time PDF is Poisson: as shown by Ascolani et al. (2009), the subordinated harmonic solution obtained from the inverse Laplace transform of Equation (18), is: The observation of the spectrum using real data suggests introducing the concept of an effective frequency ω eff that, in the case of Poisson subordination, on the basis of Equation (21) reads: In Figure 2 we fit the numerical approach to subordination with the function: using A, λ and ω eff as fitting parameters, where we have used << 1. Note that the analytical structure of Equation (23) is the same as Equation (21). The agreement is very good except in the long-time limit. The heuristic interpretation of this numerical result is that adopting µ > 2 generates an extended time region where the process is essentially Poisson. But in the long-time limit, as suggested by Equation (19), the survival probability (t) must appear. This explains the disagreement between the fitting and the numerical result in the long-time limit. The IPL index of the survival probaility (t) is µ − 1. This is a consequence of the fact that here we are considering the brand new survival probability. The spectrum evaluation done with very extended time sequences requires the aged version of Y(t) of Equation (23) which, in turn, in the long-time limit involves the aged survival probability with µ−2. This yields, as we shall see directly, Equations (2) and (3).
For the evaluation of the spectrum S(ω) for a time series of length L we adopt the prescription: which is equivalent to: where It is convenient to point out that the correlation function, under the integral in Equation (26), is not stationary, but does become stationary in the long-time limit. We expect that if the length L of the time series is sufficiently large, the spectrum defined by Equation (24 ) reduces to the ordinary Wiener-Khinchin prescription: where ξ (t) is the aged correlation function. The aged correlation function is a damped oscillatory function, with a tail characterized by the aged survival probability, namely an IPL with index µ − 2. We expect, therefore that for f → 0 the results of Equations (2) and (3) are recovered. The numerical results of Figure 3 confirm our heuristic arguments. In fact, we see that the in the region of very small frequency the theoretical prediction γ = 3 − µ is satisfactorily confirmed. In the less precise case depicted on the bottom of the figure, we run µ = 2.8 and the numerical results correspond to an effective µ = 2.74. In the region of very large values of ω, the exponential decay of the damped oscillations generate γ = 2. Between the two IPL regions a bump appears, this being a clear signature of the effective frequency that becomes smaller by decreasing µ. This is qualitatively expected on the basis of the Poisson case of Equation (21). In fact, decreasing µ decreases the rate of event production, which corresponds to the effective frequency ω eff = r of Equation (21). Decreasing the rate of event production reduces the effective frequency.

POWER SPECTRA FROM REAL DATA
In this section we illustrate the power spectrum of Equation (24) of two individuals, C1 and C2, practicing Chi meditation and two individuals, Y1 and Y2 practicing Kundalini Yoga mediation. To make our procedure clear in Figure 4 we depict samples of the HRV time series analyzed herein, referring to individuals Y1 and C2. The data is obtained from Physionet Goldberger et al. (2000), previously analyzed using very different techniques by Peng et al. (1999).
The Chi meditation is realized by instructing the subject to breath spontaneously while visualizing the opening and closing a perfect lotus in the stomach.
The Kundalini Yoga consists of a sequence of breathing and chanting exercises while seated in a cross-legged position. In the experiment of Peng et al. (1999) the Chi meditators are novices and the Kundalini Yoga meditators are advanced practitioners. We see that the main qualitative difference between Kundalini Yoga and Chi meditation is that at the onset of medation the rate of heartbeats in the former case significantly increases, whereas in the latter it does not. However, during both types of mediation coherent-like oscillations appear.
In Figure 5 the power spectrum S(ω) for the case of Chi mediation is depicted. We see that for both subjects, C1 and C2, the transition from the pre-meditation to the meditation condition induces significant periodicity bumps. Comparing the right to the left panels, a bump appears at ω ≈ 0.34 radians per second for C1 and at ω ≈ 0.38 radians per second for C2, each induced by meditation.
Note that with the method of scaling evaluation presented in section 4 the IPL index µ of C1 is seen to change from 2.4, before meditation, to 2.55, during meditation. Similarly, the IPL index µ of C2 is observed to change from 2.28, before meditation, to 2.54, during meditation. Similar properties are found in the case of Kundalini Yoga meditation, Figure 6. In the latter case the periodicity bump during meditation for Y1 is ω ≈ 0.4 radians per second and for Y2 is ω ≈ 0.65 radians per second. Notice that there exists periodicity bumps for both Y1 and Y2 before meditation. These two subject are experienced practitioners and the fact that periodicity bumps exist, even before meditation, suggests that the practice of meditation yields permanent physiological effects. As far as the IPL index µ is concerned, the analysis done in section 4 shows that the individual Y1, before meditation, has the IPL index µ = 2.64, which moves to µ = 2.68. The individual Y2 has the IPL index µ = 2.66 before meditation and the IPL index µ = 2.72, during meditation. We shall discuss this shift more fully in section 5.
The similarity of the spectrum in Figure 6 with the idealized surrogate spectrum of Figure 3 is impressive. Actually, the top FIGURE 5 | The power spectra S(ω) of Chi meditators (C1, C2) before, on the right, and during meditation, on the left. In the left panels, the spectra for C1 and C2 meditators consist of significant bumps at ω = 0.34(rad/sec) and ω = 0.38(rad/sec) respectively, each induced by meditation.
FIGURE 6 | The power spectra S(ω) of Kundalini Yoga meditators (Y1, Y2) before, on the right, and during meditation, on the left. Notice the spectra for Y1 and Y2 meditators consist of significant bumps both before and during meditation. The bump shifts for Y1 meditator from ω = 1.89(rad/s) before meditation to ω = 0.4(rad/s) during meditation. Similarly the bump shifts for Y2 meditator from ω = 1.91(rad/s) before meditation to ω = 0.65(rad/s) during meditation. panel of Figure 6 has been realized using for the subordination process with the IPL index µ = 2.55, which is the mean value of the IPL index generated by Chi meditation and the bottom panel with the value µ = 2.8, which is the mean value generated by the Kundalini Yoga meditation. See section 4 for details.
We again note that the spectral analysis had been previously done by Peng et al. (1999), who were the first to find the meditation-induced bump. They did not, however, discuss the frequency region to the left of the bump, which according to our analysis is theoretically determined by crucial events. It is important to notice that Figure 3 is obtained by making averages using many numerical realizations. The information concerning the corresponding values of µ generating the spectrum to the left of the periodicity bump cannot be derived from the experimental spectrum. The determination of this information requires adopting the technique Bohara et al., 2017) used by us in section 4.

SEARCH OF CRUCIAL EVENTS
As mentioned earlier, the search for crucial events is carried out using the method of stripes, developed earlier Bohara et al., 2017). We observe the sequence of heartbeats recording the beat numbers in the HRV time series along the abscissa axis and the time interval between beats, T B (i), along the ordinate axis. The ordinate axis is divided into stripes of size s = T B ≈ 30 ms  and the crossing from one stripe to the nearest neighbor stripe is recorded as an event, but one that is not necessarily a crucial event. We expect, in fact, that the time interval between consecutive crucial events is filled with pseudoevents, as discussed for the surrogate data. Any event, either a pseudo-event or a crucial event, is assumed to make a random walker jump ahead by a fixed quantity equal to 1, thus generating a diffusion-like process. In spite of the fact that crucial events are rare, the long-time limit of this diffusion process is dominated by their scaling with index δ (Grigolini et al., 2001). Note that according to the scaling prediction of Equation (14), the IPL index is given by: The scaling index δ is evaluated numerically using Diffusion Entropy Analysis (DEA) (Grigolini et al., 2001). The pseudoevents imbedded between crucial events are characterized by Hamiltonian Memory, a form of memory additional to that of crucial events , as discussed earlier. The Hamiltonian Memory has the effect of violating the renewal condition C(1) = 0.  assumed that the ideal condition of healthy heartbeats, with crucial events dressed with the additional memory of pseudo-events, yields C(1) ≈ 1. HRV time series also host totally random Poisson events, with probability 1− ∈. Consequently, the evaluation of ∈ 2 is done by evaluating C(1) of Equation (4) in the case of real sequences and setting the condition:  ; Bohara et al. (2017) to the Physionet data of Peng et al. (1999). It is important to stress the different effects of different forms of meditation. Let us compare the results illustrated in Figure 7, Chi meditation, to those illustrated in Figure 8, Kundulini Yoga meditation. In both case there is a tendency for meditation to increase the IPL index µ. In the case of Chi meditation the only subject not changing µ is C7. For all the other subjects, δ moves from above the dashed line to below it. The change in µ goes from the minimal change, for C1 and C5, changing from µ = 2.47 to µ = 2.54 to the maximal change, for C2, changing from µ = 2.28 to µ = 2.54. For Kundalini Yoga meditation the minimal change is for Y1, from µ = 2.64 to µ = 2.68. The maximal change is for Y4, moving from µ = 2.66 to µ = 2.85.
Yoga meditation, in addition to increasing µ, has the surprising effect of significantly increasing ∈.  conjectured that 1− ∈ is an indicator of stress that will eventually result in heart failure. If this interpretation is valid, we would infer that yoga is an efficient method to reduce stress. Apparently Chi meditation is not as efficient in reducing stress as that of Kundulini Yoga, since for some individuals, C3, C4 and C8, ∈ 2 is decreasing rather than increasing.
Actually this latter inference might be misleading and Chi meditation also has the effect of decreasing stress, although with less intensity than Kundulini Yoga meditation. As explained in earlier work (Bohara et al., 2017;, the criterion adopted to evaluate ∈ is not an exact prescription and is based on the assumption that in the absence of Poisson events we have C(1) ≈ 1. As an effect of meditationinduced coherence this condition is violated and the value C(1) in the absence of Poisson events is much smaller than 1 and decreases as the IPL index µ comes closer to 2. This important property is illustrated in Figure 9. The adoption of a more precise way to evaluate ∈ may lead to the conclusion that Chi mediation also makes ∈ increase for all practioners.

CONCLUSIONS
In addition to the well known fact that meditation generates "exaggerated heart rate oscillations"  we find that meditation does not suppress the occurrence of crucial FIGURE 8 | DEA scaling δ, IPL index µ and ∈ 2 of the HRV time series of four different participants before and during Kundalini Yoga meditation.
FIGURE 9 | Correlation rate C(t = |i − j|) of Equation (4) for the subordinated oscillator with = 0.8 . From the top, the blue, red and green curves are for µ = 2.8, 2.5 and 2.2 respectively. As µ decreases C(1) decreases as labeled.
events. As pointed out in an earlier publication (Bohara et al., 2017) these crucial events are usually found embedded in a cloud of uncorrelated and irrelevant events, prompting us to introduce the notion of dressing. A dressed crucial events is one embedded in a cloud of Poisson events, whereas the time interval between bare crucial events is empty. The exaggerated heart rate oscillations are the form of dressing triggered by meditation. The adoption of the stripe technique of analysis (Bohara et al., 2017) made it possible to quantify the level of dressing by means of the measure ∈ 2 . This enabled them to interpret 1− ∈ as the percent of Poisson events affecting the HRV time series.  conjectured that the Poisson events are generated by stress and that an excessive amount of such events as measured by ∈ 2 , i.e., stress, may be the cause of heart failure. The results depicted in Figure 8 show that the Kundalini Yoga meditation yields, using the above interpretation, a significant reduction in the effect of stress. Consequently, this paper affords a technique of analysis of HRV that may be used to directly quantify the benefits of meditation, which are currently evaluated indirectly, through the observation of psychological health (Gard et al., 2014). Gard et al. (2014) presented a very interesting discussion of the bottom-up and top-down mechanisms activated by yoga. We find that the meditation-enhanced cognition they discuss has the important effect of making the IPL index µ move from a region close to the border µ = 2, which corresponds to ideal 1/fnoise, see Equation (5), to a region closer to the border with the Gaussian basin of attraction, µ = 3 (Annunziato and Grigolini, 2000).
Does meditation contribute to making physiological processes more complex, or less complex? We posit the plausible conjecture that meditation makes physiology more complex (Sarkar and Barat, 2008;Bhaduri and Ghosh, 2016). The observation made herein that the heartbeat dynamics, as manifest in HRV time series, is a mixture of Laplace determinism and crucial events leads us to the conclusion that meditation does, in fact, decrease the complexity of HRV. The more pronounced meditation-induced coherence is an enhancement of Laplace determinism and the significant increase of the IPL index µ signifies a tendency to move toward the Gaussian basin of attraction. The increase of ∈ 2 , with a significant reduction of the percentage of Poisson events, facilitates the creation of regular oscillations.
The transition from a region close to µ = 2 to one closer to µ = 3 is a departure from the ideal 1/f -noise and fits the observation made by psychologists (Dotov et al., 2017) that, in accordance with the phenomenological philosopher Heidegger (1962), cognition is activated by addressing a difficult task. This is in line with the observation made by Correll (2008) that the rapid response to a difficult social identification task often reveals unconscious bias; one measure of the task's difficulty being the deeper cognitive transition from 1/f -noise to white noise, corresponding to µ ≫ 1. Meditation, producing a transition from the ideal 1/f -noise condition to the Gaussian basin of attraction, which in the Poisson limit µ → ∞ yields white noise, is a sign of a difficult task that yields relaxation and stress reduction, while the difficult task of bypassing the unconscious bias may be a source of stress.
Note that phenomena such as dreamless deep sleep (Allegrini et al., 2015) and brain dynamics under anesthesia (Stramaglia et al., 2017) yield a departure from the criticality condition, and this is interpreted as favoring a lack of sensitivity to perturbation (Solovey et al., 2015) in line with other recent results (Krzemiński et al., 2017). According to Varela and his co-workers Buddhist philosophy, a source of meditation techniques similar to those analyzed in this paper, is of fundamental importance in addressing the ambitious issue of cognition (Varela et al., 2016). For this reason we believe that the present research approach may contribute to the progress of cognition science.
We found also that the long-term practice of meditation has the effect of making permanent the meditation-induced physiological changes, making the periodic bump in the spectrum become permanent.

HRV Dynamics
The nature of HRV dynamics is the object of lively debates (Glass, 2009;Sassi et al., 2009;Valenza et al., 2017). The question of whether or not the normal heart rate is chaotic (Glass, 2009) is leading investigators to refine the concept of complexity variability using the tools of non-linear dynamics (Valenza et al., 2017) and to restate the multifractal nature of heart rate variability (Sassi et al., 2009). We focus on the important role of crucial events, advocated in 2002 by . More recently the authors of Bohara et al. (2017) and Mahmoodi et al. (2017a) proved that the occurrence of crucial events activates multifractal fluctuations, thereby raising the question of what is the origin of crucial events. Herein we rely on the validity of the hypothesis that crucial events are generated by the phenomenon of Self-Organized Temporal Criticality (SOTC) (Mahmoodi et al., 2017b). The subordination theory of section 2 makes it possible to simplify the much heavier computational approach to the selforganization of interacting units, each of which is characterized by periodicity. For this reason much of our future research work will be devoted to establishing the connection between subordination and SOTC.

Heart-Brain Communication
Rapid progress is currently being made in experimentally establishing the phase synchronization between different areas of the brain (Kitzbichler et al., 2009). Phase coupling has also been established in brain-heart communication (Pfurtscheller et al., 2017), suggesting that a central command may exist. The research of Ivanov and coworkers (Bartsch et al., 2014;Liu et al., 2015;Lin et al., 2016) are of importance to properly plan the future research direction of research based on the results presented herein.
To explain why it is so, it is convenient to notice that according to the recent review paper of Mather and Thayer (Mather and Thayer, 2018) individuals with high HRV have better feelings of well-being than individuals with low HRV. This suggests a communication process with information moving from heart to several emotion regulating regions of the brain (Etkin et al., 2015), thereby implying breath control, an important ingredient of biofeedback techniques (Lehrer and Gevirtz, 2014), to be the source of therapeutic efficiency. Doubts have been raised regarding this interpretation, based on the conjectures that the enhanced oscillatory behavior may reflect "systemic fluctuations, rather than neuronal connectivity" (Tong et al., 2015). There is a connection between these doubts and the results of Liu et al. (2015) proving that the wake state corresponds to the largest number of neural links, which we conjecture to be a condition equivalent to criticality and to the generation of 1/f noise. Our results on the meditation-induced transition from values of µ close to the ideal condition µ = 2, generating 1/f noise, to µ = 3, the border with the Gaussian basin of attraction, are compatible with a significant intensity increase of coherent oscillations. Figures 2, 3 show some signs of this effect that is made much more evident by Figure 10.
Unfortunately no reliable theory yet exists to establish the true origin of this phase coupling. We are convinced that the subordination perspective rests on a process of selforganized temporal criticality (SOTC) (Mahmoodi et al., 2017b) and preliminary work has been completed that support this conjecture (Mahmoodi et al., unpublished).
It is interesting to notice that the adoption of subordination-SOTC perspective is expected to lead to results in line with other interesting findings (Lin et al., 2016) showing an alternate behavior from the brain leading heart to that of the heart leading the brain, in conflict with the impression that the emotion regulating process is based on information moving from heart to the brain.
The main result of the present paper that oscillations of increasing intensity are connected to moving from criticality to subcritical condition seems to be in conflict with the intuitive belief that meditation increases cognition, but it is not. Keep in mind that a distinction must be made between observations made during meditation and those made after mediation. Teper and Inzlicht (Teper and Inzlicht, 2013) find that meditators show a stronger executive control. This observation requires further theoretical scrutiny (Sedlmeier et al., 2018) and affords an incentive to move along the lines of our work with properly designed psychological experiments.

Breathing, Meditation, and Cognition
On the basis of results presented herein, we plan to address the problem of brain-heart communication experimentally to establish what effect meditation may have, if any, on any presumed action of the central command hypothesized by the authors of Pfurtscheller et al. (2017). The recent review (Mather and Thayer, 2018) that suggests "heart rate oscillations can FIGURE 10 | Spectra corresponding to cosine with = 0.77 subordinated to an IPL PDF with a changing µ. enhance emotion by entraining brain rhythms in ways that enhance regulatory brain networks, " implies that the central command is exerted by the heart, whose regular oscillations are stimulated by the breathing process. This entails focusing our attention on Kundalini Yoga, since Kundalini Yoga rests on breathing control . We plan to analyze EEG with the same technique as that used in this paper to assess whether or not brain dynamics and heartbeats share the same restricted range of µ values. If they do, this would suggest modeling the transmission of information from heart to the brain and back using the theoretical tool of complexity matching (Mahmoodi et al., unpublished) based on the subordination-SOTC perspective used in this paper. We expect that the adoption of this theoretical approach may make it possible for us to explain the meditation-induced enhancement of alpha waves in accordance with the observation (Chandra et al., 2017) through the use of Kriya Yoga, which is different from Kundalini Yoga.
Kundalini Yoga is known to be an efficient way to treat certain psychiatric disorders (Shannahoff- Khalsa, 2004), but the therapists using this technique are looking for further enhancement of beneficial results (Khalsa et al., 2015). We anticipate that the adoption of the analysis techniques of this paper may facilitate the progress being made in this important field of research.