{'fraction_non_alphanumeric': 0.031214797049311384, 'fraction_numerical': 0.007812783894763618, 'mean_word_length': 4.36556833690778, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-organized criticality in boson clouds around black holes

Boson clouds around black holes exhibit interesting physical phenomena through the Penrose process of superradiance, leading to black hole spin-down. Axionic clouds are of particular interest, since the axion Compton wavelength could be comparable to the Schwarzschild radius, leading to the formation of"gravitational atoms"with a black hole nucleus. These clouds collapse under certain conditions, leading to a"Bosenova". We model the dynamics of such unstable boson clouds by a simple cellular automaton and show that it exhibits self-organized criticality. Our results suggest that the evolution through the black hole Regge plane is due to self-organized criticality.


I. INTRODUCTION
Modern particle detectors are typically associated with very high energies, like the current design of the LHC with 8 TeV center of mass energy. The need for having high energies is obvious if the particles to be discovered -like the Higgs -have a substantial rest mass or, conversely, a tiny Compton wavelength.
The present work is based on the working assumption that suitable axions exist in our Universe (we shall be more precise what "suitable" means in the body of the paper). In that case an axion cloud forms around a rotating black hole, and the combined system black hole/axion cloud can be thought of as a gigantic "gravitational atom", with the black hole playing the role of the nucleus and the axion cloud playing the role of the electrons.
As explained in [5], the dynamics of this "gravitational atom" is governed by the Penrose process of superradiance [18]. This process leads to exponential growth of the occupation numbers of the "atomic" levels, i.e., the formation of an axion Bose-Einstein condensate. The energy required for the formation of this condensate is extracted from the black hole, which thereby spins down and also loses some of its mass M . In fact, the black hole parameters (M, a) move on a Regge trajectory, where a ∈ [0, 1] is the Kerr parameter. Eventually, when the at-tractive axion self-interactions become stronger than the gravitational binding energy, the axion cloud collapsesa "Bosenova" occurs.
The aim of our present work is to map relevant aspects of the dynamics of the black hole/axion cloud system to a simple cellular automaton, and to show that it exhibits self-organized criticality. If this conclusion holds also for more realistic cellular automata it would mean that self-organized criticality governs an essential part of the dynamics of gravitational atoms. This paper is organized as follows. In section II we review some preliminaries, namely crucial properties of self-organized criticality and the black hole/axion cloud system. In section III we present the working assumptions of relevance for our cellular automaton and argue to what extent they are plausible. In section IV we present results based upon computer simulations implementing this specific cellular automaton. In section V we summarize and point to open issues.

A. Self-organized criticality
The concept of Self Organized Criticality (SOC) was proposed as an organizing principle of Nature, as it was realized that numerous spatial extended systems, in physics, biology, social media and economics, exhibit a number of properties which may be shortly characterized as flicker noise for the temporal evolution and selfsimilar (fractal) behavior for the spatial evolution [19]. Ever since this seminal paper, the analytical and numerical methods developed to be able to tackle SOC in various contexts have led to progress in understanding the nonlinear dynamics of complex, interacting systems. It is expected that SOC develops in slowly driven interaction dominated threshold systems and it has been shown that this state is an attractor of the dynamics of such systems (see for instance [20] and references therein).
Due to the complex interactions between its constituents, after a transitory phase in its evolution, the system exhibits what may be called a phase transition to a state characterized by the absence of a critical length scale, i.e., in this state all lengthscales are relevant. The archetypal SOC system is driven by an external force with a very large characteristic timescale and the system evolution is determined by local rules, i.e., only near-neighbor interactions are considered [20,21]. One important point to make is that this external driver is not artificial/conscious but is intrinsic to the physics in the system. This property of SOC is unlike the phase transitions in an Ising type system, where the phase transition is acquired by a conscious tuning of parameters to a critical value.
Because this line of research is still in its infancy and due to the associated analytical complexity, analysis of such systems relies heavily on numerical methods and so they have become intertwined in an inseparable way. It is difficult to enumerate properties of SOC without relating them to the way they are simulated in computer experiments.
Systems which exhibit SOC are discrete (in general modelled by lattice models). The interaction between constituents are usually nearest neighbor-like and governed by local rules and are modelled as interaction between sites in a lattice, such that if one site receives an information, only a few other selected sites have access to this information. The entire system is subjected to an infinitely slow external drive and due to this drive and the complex interaction, the onset of SOC usually occurs when all the sites in the lattice host a near-threshold value of some parameter. If the threshold is reached, the system becomes subjected to different rules of evolution (critical flow) which allow information to be transmitted not only locally, but across all lengthscales. If one site has reached threshold and sends information to the next site, this next site might also reach threshold and so on. The avalanches of information from one seed site to a final site (i.e. where the threshold condition is no longer satisfied) form an event. The size of an event is the total number of avalanches occurring from one seed site. Natural systems that reach SOC maintain it for an indefinite time as long as the external driver keeps acting. This equilibrium state is maintained by a balance between driving and dissipation, which occurs at the boundaries or in the bulk.
The way that SOC was historically "detected" in Nature was due to the properties exhibited by observables of the system. These observables show a powerlaw distribution of values. It was shown [20,21] that, for a system in SOC, the distribution function D(s) of some parameter s is given by where G(x) is the universal scaling function for a given universality class of systems, a is a system specific con-stant, τ is the critical exponent, b is a system dependent amplitude, L is the system spatial size and σ is determined by the physical dimension of the parameter s. A very useful feature of SOC is that if the distribution of another observable, say y, is considered, then D(y) will be equivalent to Eq. (1), with s replaced by y and different parameters a, τ , b and σ, but the same scaling function G. This is conceptually a powerful tool for observational purposes, in cases where it proves easy to record D(s), but difficult to record a perhaps more interesting distribution D(y).
An example of a system in SOC and a discussion of the properties emphasized above is given in the Appendix.

B. Massive boson instabilities near a rotating black hole
It has been shown [5,12,22] that the presence of bosons near a rotating black hole (BH) will lead to the formation of a boson cloud (BC) around the BH if a series of conditions are met. This BC-BH system is analogous to an atom. The cloud is being continuously fed with bosons due to the superradiance phenomenon. At some point, when the number of bosons in the cloud is too large, the nonlinearities in the BC lead to an unstable behavior and the cloud goes through a Bosenova-type event.
The conditions for this chain of events to occur refer to a set of parameters characterizing both the boson and the BH, {M, ω h , µ, ω R , m}, where M and ω h are the mass of the BH and the angular velocity at the horizon, respectively, and µ, ω R and m are the mass of the boson, angular velocity of the associated wavepacket and "magnetic" quantum number of the boson. The first condition requires that a boson exists such that its associated Compton wavelength, λ b , is approximately equal to the characteristic size of the BH, R g (the Schwarzschild radius). It was shown that the QCD axion can satisfy this condition [in e.g. [5]]. Superradiance, i.e., continuous feeding of the cloud with "free" bosons, the phenomenon analogous to the Penrose process for fields, occurs if the superradiance condition is met, namely if ω ∈ [0, ω h m]. Once superradiance sets in, the topology of the growing cloud is given by a frequency ω R , which satisfies the superradiance condition and is the frequency of hydrogenic-like bound levels. Growth of the BC, instability and subsequent Bosenova occur if the frequency associated to the boson as it interacts with the BH has a very small imaginary part, Γ g , such that the frequency describing the wavepacket in the potential well of the BH is ω = ω R + ıΓ g . The parameter Γ g , giving the growth of the BC, is a function of the set of parameters mentioned above and it exhibits a maximum for specific values of these parameters [5].
We are now in the position to define the "suitable axion" mentioned in the Introduction. This axion is a boson which possesses physical characteristics such that λ b ≃ R g , ω R ∈ [0, ω h m] and, in the potential well of a BH, exhibits a growth rate Γ g which is very small, ω R ≫ Γ g . The growth rate also depends on details of the bound system, such as the hydrogenic-like quantum numbers, taken here such that Γ g has its maximum possible value. We assume that this boson exists and give a general outline of the dynamics of the BC-BH system. Superradiance sets in and the number N of such bosons in the vicinity of the BH is amplified, in the linear limit, as However, during this amplification more complex processes occur. Most importantly, the nonlinearities start increasing and the bosons begin to interact with each other. The end result is that the number of axions present in the cloud is reduced, with a rate Γ d (d from dissipation), where Γ d ≪ Γ g . For our purposes, Eq. (2) can be naively modified as Even in this very simple description it is possible to recover some of the qualitative results derived rigorously in [5]. Since dissipation does not equilibrate growth, the number N of bosons grows until the back-reaction on the frequency ω R is such that the shape of the cloud departs from hydrogenic (nonlinear limit). This has been shown [5,12] to occur when where M BC is the mass of the boson cloud. The cloud grows further, but in a nonlinear regime, until the nonlinearities in the boson self interaction potential are too large and the cloud is no longer stable. A large percentage of the cloud is supposed to collapse [5,12] (critical limit). This collapse occurs on a timescale of r g of the BH, meaning that from the moment the critical limit has been reached, it takes r g units of time until a significant part of the cloud collapses into the BH.
The subtlety here is that this process occurs following an avalanche of information spreading on all lengthscales accessible to the system (i.e. with an upper bound given by the system size).
If the superradiance condition is still met, the cloud regrows within the same parameter space {M, ω h , µ, ω, m}. In a plot of the BH spin vs. M µ (a Regge plot, [5], their Fig. 3), the subsequent collapses and regrowths, with the same values of the parameters that initially satisfied the superradiance condition, resemble an almost straight line on which the systems wanders (Regge trajectories). It stays on this line for tens to hundreds of e-fold times of boson cloud growth, i.e. for tens to hundreds of Bosenovas. When the superradiance condition is no longer met with the initial parameters, the system makes a transition to another set of parameters (i.e., a new line on the Regge plot) which are again superradiant and analogous dynamics occurs.

III. CELLULAR AUTOMATON FOR BOSON CLOUD AROUND BLACK HOLES
We want to map the dynamics of BC-BH system into SOC. As was mentioned previously, the concept of SOC and the simulation of SOC are intimately connected and the only contained way to map the BC-BH dynamics to SOC is to simultaneously discuss the computer simulation associated with this analysis. SOC simulations are generally done using Monte Carlo methods on discrete grids, and we will focus on a Cellular Automata (CA) approach. The CA is a "machinery" which stores values of a parameter on a discrete grid and evolves these values in time and space according to simplistic rules designed to mimic a real phenomenon [23].
The CA used in the remainder of this paper is 2dimensional cartesian grid, of size L x × L y . While this undoubtedly simplifies the numerical effort, there is also strong physical argumentation in favor of considering a 2D grid (instead of a 3D grid). The BH has a significant spin, the quantum numbers of the atom are fixed (by assumption of the "suitable axion") and the orbits are Keplerian, which means that angular momentum is conserved thus leading to planar motion.
Three questions then arise: 1) what is the quantity which will be placed on the grid, 2) whether it is physically acceptable to discretize that quantity in this manner and 3) what are the evolution rules? We proceed by answering the first two questions here, while the third one will be the subject of Section IV.
Our simulation procedure is based on the following reasoning: if the number of axions in the cloud can be computed by some means [as in [5], their Eq. (28)], a local quantity can be defined by dividing the number of axions by the area of the cloud. Because the area of the cloud is considered to be constant, simulating the evolution of the density means simulating the evolution of the number of axions, and this is why we will refer to the quantity in the grid as "the number of bosons". This basically means that the N bosons in the L × L cloud will be accommodated into aL ×L grid, where we will from now on denote L x = L y =L. So it follows that one cell in the grid hasN = N/L 2 bosons andN will be the parameter stored in the CA.
To answer the more important question if this discretization is physically acceptable, consider that the number of bosons, N , or more precisely the N very poorly localized wavefunctions allowed by the superradiance, may be replaced by wavefunctions which peak sharply, i.e. are more localized for simulation purposes (Fig. 1). This "effective" system has the same information content from the point of view of an outside viewer looking at gravitational radiation following collapse. One will see in the following algorithm/simulations that the information released following collapse does not depend on the actual position of where the information was released. Thus, from the point of view of an external observer equipped only to record information on gravitational radiation, it is not important if an axion left the cloud by withdrawing its wavefunction from the entire cloud or only from a more localized spatial extent as in our simulation. We address now how the temporal evolution of the BC-BH system looks like in this framework. The linear dynamics of the BC-BH are implemented and are the evolution rules of the CA.N is increased at randomly chosen positions at each timestep in order to simulate superradiance andN i,j (k) is the value of the variable at the grid-cell labelled by spatial indices i and j, at time step k, as explained in detail in Section IV. We record the value of a given parameter s characterizing the system and check if the distribution function of this parameter, D(s), behaves as for a SOC system.
We summarize our assumptions • There is only one populated energy level. This is equivalent to saying that Eq. (3) expresses, in a linear approximation, the dynamics of the cloud itself. Such an assumption is justified, because the growth rate of the next available energy level is very small compared to Γ g .
• The dissipation rate Γ d is constant as the cloud density grows. By assumption, dissipation is negligible compared to growth, so even if Γ d does change, it will not change as much as to affect the dynamics [5].
• The nonlinear and critical limit are the same and occur when the parameter ǫ reaches its critical . This is supported by simulations of collapse [12].
• the size of the cloud, L, is constant and is calculated based on the linear approximation of hydrogenic states.
The mapping between BC-BH and SOC/CA proposed in this paper is summarized in Table I.

IV. RESULTS
We address now the third question posed in the previous Section by formulating the evolution rules for our CA. Let L be the cloud radius, i.e., the r c quantity from Eq. (11) of [5] L ≡ r c = (n + l + 1) 2 (µR g ) 2 R g , with n = 0, l = 1, µR g = 0.42 and µ = 3 · 10 −11 eV . n and l are quantum numbers associated with the hydrogenic state of the cloud and these particular values assure maximum growth. In this case, the constant [24] most effective instability rate occurs when [5] Γ g = 1.5 · 10 −7 rg −1 .
The solution of the differential equation (2) is given by where N i is the initial value of the variable N . For simulation purposes, temporal evolution is recorded in discrete timesteps, ∆t. The solution at timestep k may be rewritten as If ∆t is taken to be of the order of one gravitational radius and the notationΓ g = Γ g ∆t is used, the equation for the evolution of the system in natural time units is formally written as The number of bosons added to the cloud due to superradiance in one natural timestep is given by Since for the QCD axion e kΓg ∼ = 1 and eΓ g − 1 ∼ = 1.5 · 10 −7 , the dynamics equation used as a basis for the simulation is It is useful to restate this equality in terms of the critical number of bosons, N c . We will do this by assuming that the initial number of bosons is some fraction x ∈ (0, 1) of the critical number, N c is easily computable as N c = M BH ǫ/µ = 2.4 × 10 16 (all quantities are calculated for mass expressed in eV). So the equation for the temporal evolution of the total number of bosons in the cloud can be written as where ∆N = 3.6x · 10 9 . This is to be simulated on a 100 × 100 grid with the help of the parameterN . One may immediately writē for all sites labelled i, j in the grid and at each timestep k. This is not general enough because so far we have no reason to believe that the spatial growth of the number of bosons is homogeneous. To circumvent this problem, the same total amount ∆N will be added to the cloud, but to a limited number of randomly chosen sites, i.e. the update rule for the CA will read for n randomly chosen sites, with N 0 = ∆N/n. If n = 100 thenN 0 = 3.6x · 10 7 . For the discretized case the critical parameter isN c = 2.4 · 10 12 . Dissipation must also be taken into account, but at a much smaller rate. This is modelled asN for m randomly chosen sites, with m ≪ n. For n = 100 we consider m = 1. For future reference we note that in numerical simulations the important quantity is actually the ratioN c /N 0N To summarize, for a QCD axion with µR g = 0.42, initial cloud population with x = 95% and after rescaling the values ofN c andN 0 so thatN 0 is of order unity (for practical reasons) we get the following parameters used in the simulationN c = 280000,N 0 = 4 n = 100, m = 1.
In words, at each time step k, a number 100 units are added to the cloud andN grows locally. The local enhancement of axion number is done randomly, i.e., 100 randomly chosen grids will receive one unit. At the same timestep, 1 unit dissipates from a randomly chosen site and it is assumed that this unit escapes to infinity (it is not fed to the BH). Adding new axions at randomly chosen positions is a nontrivial assumption, but it seems well justified based on two considerations: the observations of laboratory Bosenovas which emphasize the stochastic nature of this process [24] and the previous statement that, until the critical limit is reached, the system evolves in a linear regime, Eq. (3).
At each timestep k all the grids are verified to see whether the critical densityN c is reached. Once the criticality condition is met, the gravitational potential is no longer important and transfer of information is due to the nonlinear interaction of the bosons. The next question is "how many neighbors does the site interact with?" so that information of collapse could be transported to that number of neighbors. We will assume it interacts with 3 sites, namely those three which are nearest neighbors and closer to the BH than the current site. The critical condition and the critical flow are implemented as The result of this simulation is the distribution D(s) of event sizes (Fig. 2). This is exactly how an archetypical SOC system behaves [19] and it is clear that after a growth period Self Organized Criticality is reached and the system stays in this state. This is a first important result.

FIG. 2:
The distribution of event size, D(s) for a QCD-axion cloud in superradiance around a BH. The flow at criticality is given by the rules (19). The cloud was initialized at 95% of its critical boson numberNc.
To test whether or not the fact that the system reaches SOC is rule-dependent, the system evolution is studied for different rules at criticality (different critical flows). The difference is that apart from the three units allowed to move to other sites, an additional number of boson units completely leave the cloud. This was done for 3 units, 6 units [Eq. (20), Fig. 3] and 12 units. For the case of six units, this is written and implemented as In all cases SOC is reached in a time period that depends on how close N i is to the critical value. The number of events needed to relax the system changes.
A more realistic simulation has to account for when the superradiance condition is no longer fulfilled. It is then we expect that SOC to cease, followed by a transitory period, and afterwards the system will again settle into a SOC state characterized by the new parameters. Due to the limited computational capabilities and the very slow growth rate, we use the fact that the system will eventually be in SOC and initialize the automaton very close to its critical parameter, x = 0.99.
To compute how much mass the BH has lost because of superradiance, we assume as an example that superradiance stops when 10% of the BH mass is lost, i.e. superradiance stops for M BH = M but immediately starts again for M BH = 0.9M . All the parameters characteristic to the system, and thus to the simulation, change accordingly. Simulations including the change in superradiance parameters were performed for all 4 sets of rules at criticality. The same qualitative results were obtained in all simulations: there is SOC on both sides of the transition, the transition is sharp and the ratio of the number of events prior to transition and following transition is approximately 2.
When exiting the superradiance condition for the parameter set governed by M , superradiance for the parameter set governed by = 0.9M instantly begins according to our rules. The results of this is thatN , which has a value aroundN M c is compared to the lower value of N 0.9M c and avalanches occur until the new value for the critical parameter is reached and SOC is installed and maintained.

V. SUMMARY AND OPEN ISSUES
Starting from rigorous analytical results describing the dynamics of a boson-cloud-black-hole (BC-BH) system, we derived a simple map from this dynamics to a cellular automaton. Within our working framework (see assumptions at the end of Section III and Table I) we have shown that the BC-BH system reaches self-organized criticality (SOC) and stays there as long as the superradiance condition is fulfilled. Since laboratory observations of the growth and collapse of Bose-Einstein condensates lead to a similar picture -albeit only for short timescales, see for instance Fig. 3 in [24] -we have some confidence that this is a physical property of our model and not an artifact of inadequate modelling. When the superradiance condition is no longer valid, the evolution rules change quantitatively and the simulated system exhibits a transition after which it again settles into SOC. By varying the rules we checked that this is a rule-independent result (in the sense that it is recovered for all the sets of rules that we have employed) which gives some confidence that this result holds for actual axion cloud near black holes.
The quantitative details of the transition are rule dependent (Figs. 2-4). For the first set of transition rules, Eq. (19), the temporal width of the transition period is quite large, 25 × 10 4 r g . For the second set of transition rules, Eq. (20), the temporal width of the transition is of the order 10 4 r g . The ratio of the medium numerical value of D(s) before and after transition is approximately 2 for four different evolutions of the system at criticality. The open issue is whether or not these aspects are related to underlying physics or are numerical artifacts.
We assumed that a new superradiance domain is instantly valid, which led to very large number of events associated with the transition when compared to the SOC state (Fig. 4). An open question is whether or not this assumption is physically justified. This should be investigated in future work.
We expect that observational data of gravitational waves from a large range of masses for BH will be able to shed light on at least two very important issues. One is the physics at transition and observational data will help devise more realistic rules for suitable cellular automata. Also, an appropriate observational database will help in determining the scaling function associated with this phenomena. Starting from a set of initial conditions, the pile grows and in the mean time the slope of the pile increases. It is the value of this slope, as seen in the Eqs. (A1), (A2), that will at some point reach a critical value. The meaning of the critical value is that if this value is reached and one more grain is added, material will slide. In the critical state the system is just barely stable with respect to further perturbations. Information about the perturbation will be known across all lengthscales in the system.
An event consists of the seed (i.e. placing of the grain) and the subsequent information (grain) transfer or avalanche until complete relaxation. Each event has an associated time T (i.e. the time it takes the perturbation to die out), a cluster size s (i.e. the total number of slidings induced by the perturbation) and a total energy (i.e. release of energy, as transfers occur based on first principles that the final state is an energetically lower state). For this simple model, the distribution of event sizes D(s) is given by [20] where L is the system size and θ(x) is the Heaviside stepfunction. To put this in the form of Eq. (1), we rewrite it as and by identification a = 1, τ = 1, b = 1, σ = 1 and the scaling function is G(x) = xθ (1 − x).