Coexistence of critical phenomena: the concept of manifold multi-spectral criticality

Prompted by the ubiquity of empirical observations of critical phenomena, often in non-equilibrium macrostates, we developed a modelling approach in which several critical phenomena coexist. Instead of a single critical point, many coexisting critical points in the system are identified, forming a one-dimensional critical manifold. Identified within our game-of-life-like heterogeneous agent-based simulation model, where agents can be created and annihilated in the presence of a catalyst, each critical point belonging to the critical manifold is associated with a multi-spectrum of critical exponents. We find this situation in non-equilibrium mixed percolation-like macrostates obeying continuous phase transitions. These macrostates are quasi-stationary, where some system characteristics are time-independent while others are not. This novel look at universality signals the existance of complexity of critical phenomena richer than described to date.

The primary motivation of this work is the everyday observation that a broad range of different critical phenomena can coexist.A necessary condition for this is the existence in the system of many critical points available for the system concerned.The system can reach these critical points following multiple independent trajectories in the multivariable system's phase space.This paradigm is a reminiscence of natural systems' ease of achieving critical macrostates.More often than not, this goes together with the observation that the system can simultaneously exhibit stationary but non-equilibrium in the statistical physics sense 1 ) and non-stationary characteristics.We will refer to such a macrostate of the system as quasi-stationary.Indeed, this work aims to model both the situation just described.
The phenomenon of criticality has been intensively studied in agent-based modelling (ABM), physics of phase transitions, and, more generally, in the science of complexity (SC) [2][3][4][5][6][7][8][9][10][11][12][13] .In physics, a multicritical behaviour has been explored for a considerable time [14][15][16][17][18][19][20][21][22][23][24][25][26][27] .Three different canonical examples are indicated below.(i) A straightforward example could be material possessing weak ferromagnetism in a critical region subjected to hydrostatic pressure 28 .Notably, the origin of the weak ferromagnetic moment is related to two opposing sublattices slightly canted concerning each other.Then the Curie/critical temperature is a critical curve, not an isolated point, depending on the pressure and terminating in the multicritical point.(ii) Another compelling example is a uniaxially anisotropic ferromagnet with longitudinal random field applied 20 .In this case-where is worth noting that both in our situation and there involved an exogenous random field-a new multicritical point was found as the intersection of three multicritical lines bicritical, tricritical, and critical end-points.(iii) Further extensive work on the topic 29 found the global phase diagram of the uniaxially anisotropic ferromagnet, which seems rather exotic.Thanks to this, the authors obtained a deep insight into the nature of the multicritical point, occurring in the five-dimensional space of the uniform longitudinal field, longitudinal random field component, transverse random field component, temperature, and degree of anisotropy.Indeed, applying an ordering field along the easy axis of magnetisation splits it into two symmetrical lines of unique bicritical points.Simultaneously the pure transverse phase of the three-dimensional space (local uniaxial random magnetic field, temperature, degree of anisotropy) spreads into two symmetrical horns.The new multicritical point is of fifth order at this stage.A particular feature of it is the decoupling of the longitudinal and transverse degrees of freedom, which holds for the associated curves of unique bicritical points.Due to the presence of random fields, the results are expected to be correct only in the vicinity of six dimensions.(iv) It is also worth noting an example characteristic for disordered systems, i.e., the existence of a multicritical point in the spin glass 30 .Generally, in the multicritical situation, the curves leading to the multicritical point are the borders of phase coexistence, and they can arise only for more than two phases 23 .
www.nature.com/scientificreports/Here we present two related issues-the first showing how our results relate to the well-known phenomenon of multicriticality.Secondly, our scenario plays out not in statistical/thermodynamic equilibrium but in quasistationary macrostates.To roughly illustrate our approach, we consider a game-of-life model of heterogeneous agents competing for survival in an exogenous catalytic random field.Agent models are increasingly explored in the context of complexity and criticality.Recent work such as by Ferreira et al 31 identified during our manuscript's revision addresses the context of human interactions under global leadership.While interesting and relevant, the model does not account for the evolution of the leader-which is the situation our work addresses.
While such an extended game-of-life model can have multiple applications, we propose a company market interpretation for the model.The content of our work is presented in the 'Manifold multi-spectral critical phenomena: results and discussion' section -concentrating on the results and followed by the 'Concluding Remarks' section.In the section 'Methodology, model, and method' the general agent model is described.There, the random state intervention acts as a catalyst, and the technological level of a leading company represents the fitness of an agent.This is discussed in more detail in the section 'Inspiration from company market' along with the pivotal market-based interpretation of our work.

Manifold multi-spectral critical phenomena: results and discussion
This section presents first-and second-order statistics, which prove the existence of a to date unobserved kind of continuous phase transition characterised by a phenomenon of multi-spectral criticality.Characteristic representation of the first-order statistics can be the average population of agents on the lattice.The second order can be the variance of this population-both in the quasi-stationary macrostate of the system.
By the 'average time-dependent population of agents, we mean the agents' time-dependent population, c(t) = N(t)/L 2 , averaged over a statistical ensemble of the macrostate's replicas obtained by Monte Carlo (MC) simulations, �c(t)� = �N(t)�/L 2 , to reduce the statistical fluctuations; it is still a time-independent quantity.Below, we present the method of determining c st (and hence c st ).This quantity can reach its pla- teau, �c(t) st � = �c st � = const .The phase transition we consider only for systems whose average number of agents/ companies, N(t) , has already relaxed to the appropriate plateau �N(t) st � = �N st � = const .In the ' Methodology, model, and method' Section below, we discuss the relaxation of the relative time-dependent average level of technology, A(t) /F(t) , showing that it also reaches its plateau, �A(t)�/F(t) st = const , for a long time, where quantity A(t) is the average technological level(fitness) of the market(system).

Quasi-stationary surface
Figure 1 shows a quasi-stationary phase diagram (multicoloured surface) in variables , q , and c st (i.e.agent activity, catalysis level and mean population of agents, respectively), at a fixed value of η (i.e.catalysis efficiency).
Figure 1.Quasi-stationary 2-dimensional surface (coloured) in variables ( , q, c st ) obtained from the simulation (at the fixed/typical value of η = 0.5 ); quantity c st denotes the mean concentration over the statistical ensemble of MC-simulation replicas.This is the snapshot picture taken from anima tion.The surface represents the quasi-stationary phase diagram embedded in the three-dimensional space.The bold black curves on the surface are related to the common critical point ( c = 0.349, q c = 0.879, �c st � = 0) , which is somewhat masked by finite-size effects.The black curves result from the intersection of mutually perpendicular planes ( , q = q c , �c st �) and ( = c , q, �c st �) with the surface-see Fig. 2. The area between these black curves marked in pink-compare projections in Figs.2-4.
By the term "quasi-stationary phase diagram," here we understand the dependence of the plateau of the relaxed average lattice population of agents, c st , on the variables and q.The population reaches a plateau in these variables (with the η variable value fixed).Below we show that the phase transition scenario plays out in the , q , and c st variables.

Manifold criticality
To reveal the phase transition, we intersect the surface in Fig. 1 with two mutually perpendicular planes: ( , q = q c = 0.879, �c st �) and ( = c = 0.349, q, �c st �) , i.e., parallel to the and q axes, respectively.The former plane contains the blue curve, while the latter includes the red one; they are shown in Fig. 2.These coloured curves are the graphical representation of formulas, Eqs.(1) and (2), respectively.We emphasise that the point ( c = 0.349, q c = 0.879, �c st � = 0) , marked with ' + , ' is an arbitrary (typical) point on the critical curve(1- dimensional critical manifold).We show this curve and '+' by top view in Fig. 3.All the other points of this curve were obtained analogously.Thus, these points correspond to the situation after subtracting the 'finite size effect' because the said power laws cut it off.The important result here is that both curves (blue and red) in Fig. 3 coincide, which allows them to be treated as one curve (with statistical dispersion).We can now systematically analyse the properties of the phase diagram surface shown in Fig. 1.
The diagram in Fig. 2 shows how the critical point, ( = c , q = q c , �c st � = 0) , (denoted by ' + ') forms.Both monotonically increasing black curves lying on (grey) mutually perpendicular planes ( , q = q c , �c st �) and ( = c , q, �c st �) , were formed as a result of their intersection with the (coloured) phase diagram surface shown in Fig. 1.Below, we show that in our model the mutual location of both planes can be selected so that the above-mentioned black curves have, in principle, a common critical point '+.' However, due to the finite size effect, the location of the '+' point may be too inaccurate.For better localisation (i.e., reducing the finite size effect), we use the following scaling laws: Eqs.(1) and (2).Indeed, points obtained in this way form coincident curves blue and red shown in Figs. 3 and 4. The smooth black dashed curve in Fig. 3 is located between the blue and red curves, thus defining an optimal 1-dimensional critical manifold.Indeed, this curve we also presented in Fig. 4.
The blue curve, fitted in Fig. 2 to the corresponding black dots lying on one ( , q = q c , �c st �) of the two grey planes, is captured by the following power/scaling law 4,32 , (1) The diagram showing the origin of a critical point marked with sign '+ ′ and located at ( = c = 0.349, q = q c = 0.879, �c st � = 0) .It is the snapshot from anima tion.Both black curves lying in mutually perpendicular (grey) planes ( , q c , c st ) and ( c , q, c st ) we took from Fig. 1 as the intersection of these planes with the (coloured) surface shown there.The blue and red curves share (to a good approximation) the '+ ′ point.The essential element here is that both curves were suitably fitted to the black curves (in the range just not below this point) using power laws obtained from Eqs. (1) and (2), respectively.Thus, this common point is a multi-spectral critical point.We have shown the course of the power curve on an analogous, exemplary plane lying between the two mentioned above in Fig. 5 below.This curve lies (accurate to the finite-size effect) in the area marked in pink on the surface in Fig. 1.The projection of this area is here (and in Figs. 3 and 4) onto a plane ( , q, �c st � = 0) marked with the same colour.Figure 3. Two close-lying curves on the ( , q, �c st � = 0) plane-blue and red small circles from the simulation with statistical errors marked.The black dashed curve defines the resultant/extrapolated curve between the above-mentioned two.For better visibility, the plot is limited to the 0.69 ≤ q ≤ 1.0 range.The common points of both curves constitute the manifold of critical points.Figure 2 shows how to create one of such points, marked by '+ ′ .We marked this (example) critical point at the intersection of mutually perpendicular thin black straight lines (horizontal and vertical).These lines form the basis of the grey planes shown in Fig. 2. The rectangular area marked in pink is the area marked with the same colour in Fig. 2. We can move this critical point along the black dashed curve and the structure shown in Fig. 2. Two more example critical points-one identified with a ' ⋄ ' symbol and one with an ' × ' symbol-are shown for illustration of this.This dashed line is described (to a good approximation) by an equation of the form = a − b q β , where a = 1.0046 , b = 0.2890 and β = 0.7933.Expanded projection of the phase diagram, represented in the three-dimensional coordinate system ( , q, c st ) in Figs. 1 and 2, to the plane ( , q, �c st � = 0) .It is supported by the projection shown in Fig. 3.The black, thin, mutually perpendicular straight lines intersecting at the point marked by'+' are the projection of mutually perpendicular planes ( c = 0.349, q, c st ) and ( , q c = 0.879, c st ) , which are oriented inwards (to keep the coordinate system clockwise).These planes are grey in Fig. 2. Between the planes we have placed here, for example, three other planes, the projections of which we have marked with diagonal straight lines in green, black, and brown.The parameter 0 ≤ k < q c marks the beginning of the green line, the parameter 0 ≤ l < c marks the beginning of the brown line, and k = l = 0 is the beginning of the black demarcation line between the above areas.For example, the green line has been drawn in the figure for k = 0.530 , and the brown line for l = 0.1 .The sign '+ ′ is the intersection of all six curves.Alternatively, the sign is the common intersection of the corresponding five planes projected onto the plane ( , q, �c st � = 0) .The z variable specifies the distance of the selected point, for example, on the green line from the beginning of this line.www.nature.com/scientificreports/ The above-given expression describes the critical behaviour (blue curve) solely in variable .It was fitted to part of the visibly rising curve of black dots, giving the critical value c = 0.3491 ∓ 0.0001 , the critical exponent ν = 0.219 ∓ 0.004 , and prefactor = 0.868 ∓ 0.004 .The mutual fit of both curves (black and blue) is suitable for almost the entire range of ≥ c .A slight deviation occurs near the critical threshold due to the finite-size effects.Also, deviations occur for large values of , which is irrelevant from the point of view of a critical phenomenon.The complementary power law, describes the critical behaviour in the q ≥ q c variable (same as above for variable).We fitted the course of the red curve in the ( = c , q, �c st �) plane to the corresponding black dots.The values of fitted quantities are as follows: q c = 0.8787 ∓ 0.0002 , prefactor Q = 0.979 ∓ 0.004 , and critical exponent β = 0.345 ∓ 0.008 , where we took the c = 0.3491 ∓ 0.0001 value obtained above.The blue and red curves share the '+ ′ point thanks to the appropriate mutual location of both perpendicular planes.The essence here is that both curves were suitably fitted to the corresponding black curves using power laws given by Eqs. ( 1) and ( 2), respectively.Thus, this common point is critical of the multi-spectral type; additional arguments are presented below.Moreover, the critical exponents ν and β characterising this critical point, controlling the boundary power laws (1) and ( 2), respectively, have offsets from the spectrum of critical exponents described later in the main text and shown in summary Fig. 10.For now, the cause of both offsets is an open question which applies to the entire 1-dimensional critical manifold.
In Fig. 3, we present two close-lying(coincident) critical manifolds (blue and red) created correspondingly by blue and red curves described in Fig. 2. That is, scaling regions of both critical manifolds overlap the entire length of the blue (shorter) curve, extending from point c (q = 0.736) = 0.892 to point c (q = 0.995) = 0.010 .Let us add that the red curve extends from point q c ( = 1) = 0.719 to point q c ( = 0.001) = 0.998 (errors in determining these quantities are up to the third significant digit after the decimal point).The space between both critical curves defines a 1-dimensional manifold of critical points (the black dashed curve).The slight split of blue and red curves results mainly from statistical errors and finite-size effects.

Multi-spectral critical behaviour of c st
We can now extend the data scaling with Eqs.(1) and ( 2) concerning boundary planes ( , q = q c , c st ) and ( = c , q, c st ) , respectively, to any oriented plane ( , q, c st ) between them perpendicular to ( , q, c st = 0) plane (see Fig. 4 for illustration).The basis of the ( (q), q, c st ) plane, e.g., the green straight line placed on ( , q, c st = 0) plane, we can present in variables ( , q) in the form of the linear relation between them, where 0 ≤ k < q c and k ≤ q ≤ 1 .For convenience, we only use non-negative values for the k directional parameter-the consequences of which we discuss below.Equation ( 3) is obtained assuming that the green straight line should pass through two points: ( = 0, q = k) and ( c , q c )-we marked this latter point in Fig. 4 by sign '+ ′ .It is the same critical point we marked in Figs. 2 and 3.It should be emphasised that the point marked by '+ ′ is only an example of a critical point on the dashed critical curve in Fig. 3.We might choose any other point on this curve for our consideration such as the ones denoted with ⋄ and ×.
The basis of the other ( , q( ), c st ) plane is the brown curve in Fig. 4. Its equation in the variables ( , q) takes the following linear form, where 0 ≤ l < c and l ≤ ≤ f = l + c −l q c , having f = (q = 1) .For convenience, we have limited the l directional parameter range to non-negative values, similar to the k parameter.This type of double parameterisation is the price we pay for using non-negative parameters k and l.We obtained Eq. ( 4), assuming that the straight brown line should pass through two points: ( = l, q = 0) and ( c , q c ) .The diagonal black straight line meeting Eqs. ( 3) and ( 4) for k = l = 0 is the boundary line separating both methods of parameterisation discussed above.Eqs. ( 3) and ( 4) are intended to offer useful parameterisation for straight lines passing through ( c , q c ) to define (below) the compact variable z conveniently.The plane ( (q), q, c st ) based, for example, on the green straight line shown in Fig. 4 intersects the (multicoloured) surface presented in Fig. 1, giving the black dots shown in Fig. 5.

Multi-spectral critical relaxations
In Fig. 6 we present the time-dependence of the population of agents on the lattice, c(t) = N(t)/L 2 obtained by MC simulations (for details see 'Methodology, model, and method' Section below and Suppl ement ary Infor mation.The curves c(t) obtained were determined for the points lying on the green diagonal straight line shown in Fig. 4.They result from averaging the corresponding curves over a statistical ensemble of 400 replicas i.e., single independent repetitions of a given system simulation.As shown in Fig. 6, the plateaus are reached at different rates.Plateaus for small and large companies'(agents') populations on the lattice have the highest rates.On the other hand, plateaus for intermediate populations (i.e., in the vicinity of z c ) are reached slower. (2) The relaxation time (as the reciprocal of the rate) has a local maximum in the variable z.Here we only discuss the properties of the system after it has reached a plateau for different values of the plateau.A striking feature of the relaxation spectrum presented in Fig. 6 is the apparent dependence of its plateaus on z, that is, both on the level of state intervention(catalyst), q, and on the level of company(agent) activity, .The main goal of this work is investigating this relationship.8) by fitting it to data (black dots) from computer simulations.These data were additionally averaged over the statistical ensemble of replicas of the initial simulation (numerical experiment).Therefore, we denote them by c st ; hence, the plot's vertical axis has such a notation.From this fit we obtained = 0.990 ∓ 0.001 and ζ = 0.309 ∓ 0.002 with k = 0.530 and z c = 0.494 , i.e., q c = 0.879 and c = 0.349 .These parameters locate the multi-spectral critical point as shown (by the ′ + ′ sign) in Figs. 2-4.The green curve is drawn along the green straight line in Fig. 4, starting from '+ ′ .Fits for points marked with ⋄ and × in www.nature.com/scientificreports/In Fig. 7, the integrated non-linear relaxation time dependence on z variable is shown at fixed η .The relaxa- tion time we can define as follows, where t 0 ≥ t ip , while t ip locates in time an inflection point (see red dots on curves in Fig. 6), while the non-linear relaxation function φ(t) takes the form, where brackets . . .mean the averaging over a time-dependent statistical ensemble of replicas.Thus, we have defined the global relaxation mode for the company market (where Eq. ( 24)) defines an example of the local mode).

Spectral multicritical growth and slowing down
The result shown in Fig. 7 is the pillar of this work and shows that there is z = z c for which the integrated non- linear relaxation time, τ , diverges in the z variable.Thus, it reveals this variable's critical growth and slowing down effects.This is supported by the observation that this divergence is of the power-law nature.Namely, solid lines on the plot are described by the following (dichotomous) formula, where the sign '− ′ concerns the left edge of the peak (i.e., the situation for z < z c described by the red curve) and the sign '+ ′ concerns the right edge (i.e., the situation when z > z c described by the blue curve), exponents ξ ∓ (> 0) are the critical dynamic ones, A ∓ denote prefactors/amplitudes, and B ∓ denote the background level.In Table 1 we have compiled the obtained values of quantities (exponents, prefactors, and background coefficients).Figure 6 presents the critical threshold we took from fit.
The plots in Fig. 8 show the spectra of both critical exponents ξ ∓ .The existence of these spectra confirms the spectral multicritical nature of the continuous phase transformations in the system.
( gives the solid lines.We can talk about the multicritical feature of this effect because we have the divergence of τ for both (red and blue curves) critical 1-manifolds presented in Fig. 3.This is what justifies the name "multicritical dynamics/kinetics/relaxation" we observe.Fits for points marked with ⋄ and × in Fig. 3  www.nature.com/scientificreports/

Phase classification
The dependence of the integrated non-linear relaxation time τ (see Fig. 7) enables us quantitatively to differentiate populations of agents/companies on the lattice.These correspond to distinct areas in the quasi-stationary phase diagram shown in Fig. 1 and in its projection shown in Fig. 3.These areas can be intuitively described as follows:

Low
The low population/concentration area (LP, defined by the almost horizontal part of the phase diagram surface in Fig. 1), where the probability of meeting two agents is residual-agents behave there as independent of each other, e.g., like free particles in dilute lattice gas.Although the agents almost do not block each other's displacements between lattice sites, the intensity of their mutual interactions is low.This means that their effective activity roughly expressed by c st is also tiny.We refer to this area of the phase diagram low-concentric or low-effective activity, whereby the term ' effective activity' means the intensity/ frequency of interactions between agents leading to mergers of companies or the creation of new companies (in the form of spin-offs).

High
The high population region (abbreviated by HP) is separated from the LP by a narrow transition region due to finite size effects-this region gets narrower the greater the lattice size.In this area, there is a rapid increase in the population of agents, resulting in a significant increase in their effective activity.This area can be divided into two parts: 1.An area of multiscaling, i.e., multi-spectral critical region, where the power-law increase of agents lattice population is in force-there is a multi-channel mutual interaction of agents, including an autocatalytic interaction/dependence.This area is the scaling region.2.An area of high population where the probability of the agents meeting is high, i.e., outside the scaling area.There, the number of vacancies on the lattice is relatively small.We use the above classification in the subsections below, i.e., the behaviour of the order parameter (or agent population on the lattice, see Fig. 1).Since we are dealing with a continuous phase transition, we discuss singularities of variance at multi-spectral critical points.We confirm the presence of singularities in this phase transition, i.e., transition between the inactive and active phase.

Divergence of variance
The green curve in Fig. 5 represents the fit to the black dots of the c st vs. z ≥ z c given by one-dimensional scaling (power law),   7)).The k and l control parameters are in the ranges 0 ≤ k ≤ q c and 0 ≤ l ≤ c , respectively.There is a large variability of the exponent ξ − compared to the exponent ξ + .Spectra are shown for all the three critical points marked with ⋄ , + and × in Fig. 3 and the z range on both sides of the criticality threshold z c -the left plot on the left (i.e., for z < z c ) and the right one on the right side (i.e., for z ≥ z c ).
Vol where the compact scaling variable, z, is given for 0 ≤ k < q c by the following formula, with its critical value z c = (q c − k) 2 + 2 c , and obtained from Eq. ( 9).Thus, the compact variable z measures the distance of a given point located on the green straight line from its beginning (the point ( = 0, q = k, �c st � = 0) ) along this line.To simplify the notation, we do not introduce additional indexing of the z variable and the z c and z f parameters using the directional parameter k and l.
For the complementary case of 0 ≤ l < c we obtain with its critical value z c = ( c − l) 2 + q 2 c .As before, the z variable measures the distance of a given point placed on the straight brown line defined by the direction parameter l from the beginning of this line (i.e., from the point ( = l, q = 0, �c st � = 0) ) along its direction.
In this case, we can formulate a one-dimensional power law complementary to Eq. ( 8), where both the prefactor and the multi-spectral critical exponent ζ depend explicitly on l instead of k.Next, (from Eqs. ( 11) and ( 4)) we obtain, The following continuity boundary condition applies, We limit ourselves to two opposite areas of the ( , q) variables, namely to (0 ≤ q < q c ) × (0 ≤ < c ) and (q c ≤ q ≤ 1.0) × ( c ≤ ≤ 1.0) .It is worth emphasising that the proposed multiscaling laws, Eqs. ( 8) and ( 12), describe universality classes indexed with the k and l parameters.The plot in Fig. 10 shows the non-uniform sub-spectrum of the ζ exponent.The borderline scaling laws described by Eqs.(1) and ( 2) are complementary to Eqs. ( 8) and (12) .
To summarise this subsection, we note that in our case, the order of the phase transition remains unchangedwe are dealing with a continuous phase transition.The essence of our phase transition, different from multicriticality, is shown in Fig. 4. The dashed line is the critical manifold, the same as in Fig. 3.

Divergence of susceptibilities
The scaling, Eqs.(1) & (2), and multiscaling Eqs. ( 8) & ( 12), laws discussed above lead to divergence of susceptibility at the criticality threshold ( c , q c , �c st � = 0) .We define susceptibility, χ , as a partial derivative, where y = and q, respectively, for expressions Eqs.(1) and ( 2), and y = z for expressions Eqs. ( 8) and ( 12).This definition results in susceptibility divergence because the values of the critical exponents ν( c , q c ), β( c , q c ) (8) c st (z; l, q c , c ) =�(l, q c , c ) www.nature.com/scientificreports/ ,ζ(k, c , q c ) , and ζ(l, c , q c ) are less than 1 in the full range of c , q c , k, and l.Exemplary (typical) values of the ν and β exponents are given in the subsection above, and the ζ exponents in Fig. 10 above (while the single value of ζ is given in the caption under Fig. 5).The course of susceptibility χ vs. z is similar to the course of variance in Fig. 9.The critical exponent controlling the descent of the right slope of the susceptibility peak has the form of exponent 1 − ζ > 0 .The figure also applies to each point of the critical manifold.

Diverging fluctuations
In the above subsections, we indicated the conditions necessary for continuous phase transitions to appear.This was the behaviour of the population-or the order parameter.We also noticed the divergence of the susceptibilities directly related to the above mentioned behaviour.In this subsection, we indicate the singular behaviour of the variance of the population, Var[c st ] , (see Fig. 9 for details) as a sufficient condition for the Figure 9.The diverging variance Var[c st (z; k, c , q c )] vs. z.The singularity of variance at the critical point z c = 0.494 , as marked in Fig. 5, is well seen.Parameters characterising the singular component in Eq. ( 16) are as follows: µ = 1.286 ± 0.065 and V = 0.00010 ∓ 0.00002 .While the parameters characterising the linear background in the mentioned formula are as follows: a = 0.0050 ∓ 0.0007 and b = 0.0032 0.0004 .This result and the one presented in Fig. 5 indicate that we are dealing with a continuous phase transition at z c .Fits for points marked with ⋄ and × in Fig. 3 look analogical and are not shown.8) and ( 12)), which depend on the parameters k and l, respectively.Spectra are shown for all the three critical points marked with ⋄ , + and × in Fig. 3.
existence of a continuous phase transition and multi-spectral critical phenomena.A variance sub-spectrum consisting of critical exponents co-creates the multi-spectra, which characterises the multi-spectral criticality (see Fig. 11 for details).
The multi-scaling hypothesis indexed, for example, by the k parameter takes the following form for z ≥ z c : where z c = z c (k, c , q c ) and z f = z f (k, c , q c ) appeared in the description of Eq. ( 9) and in Eq. ( 10), respectively.Moreover, the background linear function k, c , q c ) = −a(k, c , q c ) • z + b(k, c , q c ) , where coefficient a > 0 .
To get the analogous scaling hypothesis indexed by the l parameter, it is formally enough to replace the k parameter with the l parameter everywhere in the above equation.In Fig. 9, we present the diverging dependence of the variance Var[c st (z; k, c , q c )] vs. z for the multi-spectral critical point, which in Figs.2-4 has been tagged with '+ ′ .Each critical point ( c , q c , �c st � = 0 ) belonging to the critical manifold is associated with a family of vertical planes.These contain the power laws described by equations Eqs. ( 1) & ( 2) and ( 8) & (12).Hence, each family of these planes is associated with a given level/order of the statistic to form a spectrum of critical exponents.One can move from one plane to another by continuously changing the values of the k or l parameters.Therefore, each critical point along the critical manifold possesses a spectrum of critical exponents and for this reason the critical manifold becomes multi-spectral.

Concluding remarks
To the best of our knowledge, this type of rich scaling phenomena along a nonlinear manifold has to date not been observed in the science of complexity, including agent-based modelling, entailing econophysics and sociophysics.More specifically, we built a numerical, heterogeneous agent-based model inspired by the market of competing companies in the presence of random state interventionism.This interventionism serves as a catalyst operating in an egalitarian manner, i.e., regardless of the technological level(fitness) of companies and market share(rank).We examined the behaviour of heterogeneous agents in the semi-stationary macrostates of the system.We observed in these macrostates a new type of non-equilibrium continuous phase transitions, generating the set of critical thresholds in the form of a multi-spectral 1-dimensional critical manifold.This manifold's essential feature is that every critical point is associated with at least four different types of sub-spectra of critical exponents responsible for (i) the power-law behaviour of the order parameter, (ii) diverging susceptibility, (iii) variance divergence, and (iv) growth & slowing-down effects.
Every critical exponent characterises a different curve in the macrostates' space on which the system can reach a given multi-spectral critical point.For example, it can be characteristic of the growth & slowing-down effects concerning the system's relaxation to semi-stationary macrostates in the scaling area.As the system approaches ( 16) Figure 11.Spectra of the µ exponent (appearing in Eq. ( 16)), which show significant oscillations.This is related to the multi-spectral critical point ( c , q c , �c st � = 0) marked by '+ ′ in Figs. 2-4.Spectra are shown for all the three critical points marked with ⋄ , + and × in Fig. 3.
the critical threshold, the observed divergence of the relaxation time within the scaling region is one characteristic and the telling system behaviour for a continuous phase transition.
We chose to study an egalitarian variant of our model, i.e., the one in which no agent is preferred initially (we treat this as the rule of no preference).As it transpired, this variant has the broadest spectrum of populations' plateaus and, what is the most significant, it leads to the multi-spectral critical phenomena depending on the compact variable z (containing the exogenous random level of catalyst (intervention), q, and endogenous single agent's(company's) conditional activity level ).The fitness of the agent (e.g., technological level of the company) controls the local dynamics/kinetics of the agent in our model [33][34][35][36] -the agent/company tries to achieve the leader's fitness (technological level), competing for survival.We treat this fitness of the leader as a deterministic, time-dependent exogenous field acting on the system.
From the point of view of the thermodynamics of the lattice gas or Ising magnetic system, the variable c st can play the role of the order parameter-concentration/population of lattice gas particles or magnetisation.There is well-known relation between Ising spin, s, and the (lattice gas) particle occupation number, c (equals 0 or 1), in the form s = c − 1/2., could be an analogon of the thermodynamic activity, q the analogon of the hydrostatic pressure applied to the system, albeit this analogy is incomplete because q is random.In contrast, the hydrostatic pressure is not.The analogy which fully holds here is that both quantities are exogenous and η is the analogon of a permeability.We have shown that the phase transition scenario plays out in the , q , and c st variables (at a fixed η).
We would like briefly to address the differences and similarities of our situation relative to those containing multicriticality.At the outset, the critical difference stems from multi-spectral phenomena occurring in quasistationary macrostates.At the same time, multicriticality occurs in macrostates of statistical equilibrium of systems.However, in our case, we are dealing only with a single 1-dimensional critical manifold.In physics, even when dealing with multicriticality, it either appears for the boundary points of a single critical manifold or as a result of connecting or crossing various critical manifolds.These points have different properties from those identified in our case.
We focused on systems with curves of critical points and not isolated critical points.Several critical phenomena relate to the n-dimensional critical manifolds, where n ≥ 1 .The critical curves sometimes intersect with other transition curves or branches of different phases or phase borders, generating multicritical intersection points.Multicritical points can also locate/terminate on the (n − 1)-dimensional border of the n-dimensional critical manifold.However, systems are still critical at this border.Multicritical points are unique points in thermodynamic or other systems' parameter/variable space with a continuous phase transition.A single critical exponent no longer describes the critical behaviour at a multicritical point.Moreover, the critical exponents can change abruptly at these points as the control parameter/variable is varied.Hence, the points on the borders usually belong to a different class of universality than those realised within a bulk part of the critical manifold.
Furthermore, in our situation, the semi-critical curves (intersecting the 1-dimensional critical manifold) lie entirely in one phase on the phase diagram.Each semi-critical curve has a power law associated with it, and many such curves end at the same critical point.This means many critical exponents (i.e., a continuous local spectrum of critical exponents) characterise a given critical point-we call it a multi-spectral critical point.Indeed, the current definition of the multicritical point does not cover our situation.However, the multicriticality known so far bears similarities to our situation because multicritical and multi-spectral critical points are characterised by many critical exponents, and not by a single one, as in the case of the usual critical point.
Criticality as a generic phenomenon is a signature of complex systems that are the scene of the struggle of opposites.Balancing these opposites is necessary (at least) for a quasi-stationary macrostate to define.However, more is needed to obtain the critical curve or-in general-the n-dimensional critical manifold.For example, for a market of competing companies involving random state interventionism, we deal (in the scaling region) with a single balance condition generating the law of companies' conservation.This law describes balancing the companies' number emerging in the market with the number of companies disappearing.In our model, two stochastic processes (expressed in terms of probability currents) are responsible for the disappearance of companies from the market: (i) a merger of companies, as a result of which two companies become one, and (ii) bankruptcy of companies, i.e., the disappearance of companies from the market.In the quasi-stationary macrostate, both these destructive processes are balanced by the company creation process.Two companies can create (with a given probability) a spin-off, i.e., a new company related to them.However, on the curves of critical points (which we mainly consider), the single balance condition mentioned above breaks down into two special balance conditions corresponding to different balancing paths: (i)The first path corresponds to (a) balancing the central part of the probability current of company creation with (b) the central part of the probability current of its annihilation (associated with the merger of companies and the full probability current of bankruptcy of companies).(ii)The second path is related to balancing the residual parts of the currents of the probabilities of creating and annihilating companies (already without the participation of bankruptcies).The transition from the single balance condition to two specific balance conditions can be seen as a symmetry-breaking phenomenon characterising critical phenomena considered in this work.
Our approach is a modern variant of the game of life 37,38 combined with growth theory 38 and references therein, together with innovation spread/diffusion 39 .Thus, our work directly relates to the collective dynamics of opinion formation for acceptances of innovations 40 ) developed in the agent-based modelling and catalyst-based competitive survival model context.We dealt with two types of the non-equilibrium (irreversible) competitive game of life: (i) agent-based catalytic game of life, where the catalyst may be the effectiveness of state intervention, and (ii) agent-based autocatalytic game of life.The catalyst enters this game indirectly through the company's technological level(fitness), while the descendant (created by a pair of agents) is the autocatalyst.This can be related by analogy to critical phenomena in catalytic and autocatalytic irreversible chemical reactions 41,42 .
Our results show critical divergence of variance along the critical manifold associated with dysregulation of the agent's system at critical levels of the catalyst.In the particular case of the market scenario, this warns against From the point of view of econophysics, the various multi-spectral critical points represent the macrostates of the system with unlimited risk because then we deal with divergent fluctuations and susceptibilities in variables q, , and z.In effect, egalitarian state intervention either does little if it is below the critical thresholds even for a relatively large state intervention level, q, or it introduces a rapid increase in fluctuations and hence uncertainty and risk in the market when it passes the multi-spectral critical threshold.Our work has universal hallmarks because it bridges and extends different branches of the science of complexity.Notably, these include the rapidly developing innovation diffusion/spreading discipline, the underresearched area of dissemination of innovative technologies in the market of competing companies and the potentially meaningful applications to chemical and biological reaction systems.

Methodology, model, and method
The agent-based model used in this work can take multiple variants and configurations 36,43 yielding itself to adaptation to the conditions imposed by a given market situation.We describe a variant of the model in which we observed multi-spectral criticality, not identified in the previous variants we investigated.Below, we only cover the elements of the model required for interpreting the results obtained.Detailed complementary information is provided in the Supplementary Information accompanying this work.

Inspiration from company market
While the agent model can have many interpretations and applications, our motivation came from agent-based modelling of the competing companies' markets in the presence of state intervention 36,43 .We are studying the market of competing companies as a primary market.In particular, the stock exchange would not exist without the market for competing companies/agents struggling to survive, which determine the dynamics/evolution of the stock exchange.State intervention in competing companies' markets can tremendously impact the economy, mainly through technological growth.Understanding of the role of state interventionism in the free market economy is still in its infancy, although this problem is as old as the free market itself [44][45][46] .

Model substrate
We set up the model as follows.First, a two-dimensional regular lattice network of size L × L is created (where L is the linear size of the lattice).Each site i of the lattice represents a discrete location where a company(agent) can, but does not have to, reside.A pair of time-dependent functions describe the local state of the occupied site i: 1.The first function defines the company's(agent's) technology level(fitness), A i (t) at site i.
2. The second function represents the company market share(agent rank in the system) ω i (t) at this site.
The company's market share corresponds to the agent's instantaneous rank in the agent community system formulated in the ABM language.At the same time, agent's fitness represents its temporary ability to survive in that community, just as the company's technological level assesses its temporary ability to survive on the market.Thus, the time-dependent microstate of the occupied lattice site is defined by the functions' pair (A i (t), ω i (t)) , i = 1, 2, . . ., N(t) , which are the endogenous micro characteristics; N(t) is the endogenous macro characteristics specifying the number of agents present on the lattice at time t.An unpopulated lattice site is in the microstate of (0, 0).Upon initialisation, the lattice sites are populated with probability c 0 subject to the condition of single occupancy per site.We chose c 0 = 0.8 , meaning that, on average, the lattice starts with 80% of sites occupied.This particular choice allows us to compare our results with previous work 35,43 (for q → 0 ).Each occupied site is also assigned the technology level(fitness) A i (0) (drawn from a uniform distribution) and assigned equal shares ω i (0) = 1 c 0 L 2 to ensure egalitarian initial condition.The share(rank) ω i (t) describes what percentage of the market(system) a given company(agent) at lattice site i holds at any time t.However, the leader defines an exogenous field; therefore, it is not included in the occupancy count.The normalisation condition, holds for any given time t.Only N(t) lattice sites are populated at time t; the remaining L 2 − N(t) sites are empty.The model assumes that the progression of technology depends on the (global) technological growth rate of the leader.Companies in different countries can copy the leader's solutions to increase technology.It constitutes the technology diffusion stimulated by the outer fields.The leader's technology is assumed to grow according to Eq. ( 19) below, following Moore's macroeconomic law of exponential growth 48 .The parameter σ is the growth rate at which the leader's technology, F(t), increases.So, we have, where we assume F(0) = 1 for further calculations for standardisation reason.

Mutual interaction and fitness (technological level)
An important assumption of our model is that, apart from interventionism, the technological level is a decisive factor in companies' survival probability and success.This follows a similar class of models [33][34][35][36] .
Let us first assume there is no other company in the nearest vicinity of a given company located in i th lattice site.It may result from a fluctuation or when the concentration of agents on the lattice c st ≤ 0.5 .By c st , we denote the concentration of agents on the lattice in the quasi-stationary macrostate of the system.Then, we can propose the linear first-order recursive equation describing the local (discrete) stochastic dynamics of the level of technology with multiplicative noise in the form: where Here, t ≥ 1 (we measure it in MCS/site units), and A i represents the time-dependent company's technology level in the i th lattice site, while F is the leader's time-dependent technology level.The situation of the system corresponding to the left branch of the binary tree of life (shown in Fig. 12), i.e., the successful intervention consisting of drawing a random number r < q , we describe by the upper branch in Eq. ( 21).The lower branch in Eq. ( 21) refers to the right branch of the binary tree of life of a subbranch of ' Activity' (not 'Bankruptcy').Thus, the increase in the level of technology in Eq. ( 20) depends on the distance between the company's level of   η, r < q 1, q ≤ r ≤ 1.

Figure 13.
Example relationships of the mean population c st vs. η for five typical q values at fixed typical = 0.9 .Such a choice of these values , q, η covers both phases-the bottom three curves belong to LP and the top two to HP. Technical note: filled circles indicate the simulation results and thin solid lines have been drawn between them to guide eye for interpolation.There is a clear dependence of c st on q which was to be expected in the light of Fig. 1 and only a negligible, residual dependence on η-the catalyst level.www.nature.com/scientificreports/Summarising, we have shown the way from the control variables, ( , q, η) , of the model through the time- dependent technology levels(fitnesses), A i (t) , to the basic survival probabilities of companies(agents), p i (t) , on the lattice site i.It is the technological level(fitness) which controls the local dynamics/kinetics of the company(agent)-the company(agent) tries to achieve the technological level(fitness) of a leader, competing with other companies(agents) for survival in the market-therefore, except for the leader, the company's(agent's) technological levels(fitnesses) are the central endogenous characteristics of our model.

Figure 4 .
Figure 4. Expanded projection of the phase diagram, represented in the three-dimensional coordinate system ( , q, c st ) in Figs.1 and 2, to the plane ( , q, �c st � = 0) .It is supported by the projection shown in Fig.3.The black, thin, mutually perpendicular straight lines intersecting at the point marked by'+' are the projection of mutually perpendicular planes ( c = 0.349, q, c st ) and ( , q c = 0.879, c st ) , which are oriented inwards (to keep the coordinate system clockwise).These planes are grey in Fig.2.Between the planes we have placed here, for example, three other planes, the projections of which we have marked with diagonal straight lines in green, black, and brown.The parameter 0 ≤ k < q c marks the beginning of the green line, the parameter 0 ≤ l < c marks the beginning of the brown line, and k = l = 0 is the beginning of the black demarcation line between the above areas.For example, the green line has been drawn in the figure for k = 0.530 , and the brown line for l = 0.1 .The sign '+ ′ is the intersection of all six curves.Alternatively, the sign is the common intersection of the corresponding five planes projected onto the plane ( , q, �c st � = 0) .The z variable specifies the distance of the selected point, for example, on the green line from the beginning of this line.Then, for instance, z c = Figure 4. Expanded projection of the phase diagram, represented in the three-dimensional coordinate system ( , q, c st ) in Figs.1 and 2, to the plane ( , q, �c st � = 0) .It is supported by the projection shown in Fig.3.The black, thin, mutually perpendicular straight lines intersecting at the point marked by'+' are the projection of mutually perpendicular planes ( c = 0.349, q, c st ) and ( , q c = 0.879, c st ) , which are oriented inwards (to keep the coordinate system clockwise).These planes are grey in Fig.2.Between the planes we have placed here, for example, three other planes, the projections of which we have marked with diagonal straight lines in green, black, and brown.The parameter 0 ≤ k < q c marks the beginning of the green line, the parameter 0 ≤ l < c marks the beginning of the brown line, and k = l = 0 is the beginning of the black demarcation line between the above areas.For example, the green line has been drawn in the figure for k = 0.530 , and the brown line for l = 0.1 .The sign '+ ′ is the intersection of all six curves.Alternatively, the sign is the common intersection of the corresponding five planes projected onto the plane ( , q, �c st � = 0) .The z variable specifies the distance of the selected point, for example, on the green line from the beginning of this line.Then, for instance, z c =

Figure 5 .
Figure 5.The dependence of the order parameter c st vs. z.The solid curve (in green) we obtained from the formula Eq. (8) by fitting it to data (black dots) from computer simulations.These data were additionally averaged over the statistical ensemble of replicas of the initial simulation (numerical experiment).Therefore, we denote them by c st ; hence, the plot's vertical axis has such a notation.From this fit we obtained = 0.990 ∓ 0.001 and ζ = 0.309 ∓ 0.002 with k = 0.530 and z c = 0.494 , i.e., q c = 0.879 and c = 0.349 .These parameters locate the multi-spectral critical point as shown (by the ′ + ′ sign) in Figs.2-4.The green curve is drawn along the green straight line in Fig.4, starting from '+ ′ .Fits for points marked with ⋄ and × in Fig 3 look analogical and are not shown.

Figure 6 .
Figure 5.The dependence of the order parameter c st vs. z.The solid curve (in green) we obtained from the formula Eq. (8) by fitting it to data (black dots) from computer simulations.These data were additionally averaged over the statistical ensemble of replicas of the initial simulation (numerical experiment).Therefore, we denote them by c st ; hence, the plot's vertical axis has such a notation.From this fit we obtained = 0.990 ∓ 0.001 and ζ = 0.309 ∓ 0.002 with k = 0.530 and z c = 0.494 , i.e., q c = 0.879 and c = 0.349 .These parameters locate the multi-spectral critical point as shown (by the ′ + ′ sign) in Figs.2-4.The green curve is drawn along the green straight line in Fig.4, starting from '+ ′ .Fits for points marked with ⋄ and × in Fig 3 look analogical and are not shown.

Figure 7 .
Figure 7.The dependence of the integrated non-linear relaxation time, τ , vs z at fixed value of η = 0.5 .A visible divergence of τ at the critical threshold z c = 0.4937516 defines our model's multicritical nonlinear growth and slowing down effects in the z variable.The black dots result from the MC simulation, while the formula Eq. (7) gives the solid lines.We can talk about the multicritical feature of this effect because we have the divergence of τ for both (red and blue curves) critical 1-manifolds presented in Fig.3.This is what justifies the name "multicritical dynamics/kinetics/relaxation" we observe.Fits for points marked with ⋄ and × in Fig.3look analogical and are not shown.

Figure 8 .
Figure 8. Spectra of both multicritical exponents ξ ∓ (appearing in Eq. (7)).The k and l control parameters are in the ranges 0 ≤ k ≤ q c and 0 ≤ l ≤ c , respectively.There is a large variability of the exponent ξ − compared to the exponent ξ + .Spectra are shown for all the three critical points marked with ⋄ , + and × in Fig.3and the z range on both sides of the criticality threshold z c -the left plot on the left (i.e., for z < z c ) and the right one on the right side (i.e., for z ≥ z c ).

Figure 10 .
Figure 10.The plots of the spectra of the ζ exponent (appearing in Eqs.(8) and (12)), which depend on the parameters k and l, respectively.Spectra are shown for all the three critical points marked with ⋄ , + and × in Fig.3.

( 17 )
N(t) i=1 ω i (t) = 1.over the system of all the individual technology levels(fitnesses) weighed by the respective shares(ranks) for each i.
t)A i (t).
∓ 12.65268 ∓ 0.55088 76.0835 ∓ 4.7216 improper state interventions, especially for the market states in the multi-spectral critical points.