Third order interactions shift the critical coupling in multidimensional Kuramoto models
Ricardo Fariello Marcus A. M. de Aguiar

	






Departamento de Ciências da Computação,
	Universidade Estadual de Montes Claros, 39401-089, Montes Claros, MG, Brazil 
	


Instituto de Física Gleb Wataghin, Universidade Estadual de Campinas, Unicamp 13083-970, Campinas, SP, Brazil


	
The study of higher order interactions in the dynamics of Kuramoto oscillators has been a topic of intense recent
research.  Arguments based on dimensional reduction using the Ott-Antonsen ansatz show that such interactions usually
facilitate synchronization, giving rise to bi-stability and  hysteresis. Here we show that three body interactions shift the critical coupling for 
synchronization towards higher values in all dimensions, except D=2, where a cancellation occurs. 
After the transition, three and four body interactions combine  to facilitate synchronization. 
We show simulations in D=3 and 4 to illustrate the dynamics.
	


Third order interactions shift the critical coupling in multidimensional Kuramoto models
    Marcus A. M. de Aguiar 
    April 29, 2024
========================================================================================





§ INTRODUCTION


Pairwise interactions are not enough to describe the dynamics of many complex systems <cit.>. 
In several cases actions take place at the level of groups of agents, as found in neuroscience <cit.>, 
ecology <cit.>, biology <cit.> and social sciences <cit.>.
Efforts to account for interactions beyond the usual two body terms have generated a large body of literature in several areas, 
particularly in propagation of epidemics <cit.> and synchronization 
<cit.>.

In the specific case of the Kuramoto model, it has been shown that higher order terms create bi-stable regions in parameter space, 
and therefore hysteresis, a feature that does not exist in the original model with pair interactions <cit.>. 
For fully connected networks and Lorentzian distribution of natural frequencies, where the Ott-Antonsen ansatz <cit.> can be applied, 
theory predicts that the critical coupling for synchronization does not change, but a line of saddle-node bifurcation marks the appearance 
of a new synchronized equilibrium that co-exist with the asynchronous state. Moreover,
above the critical coupling, higher order terms facilitate synchronization.

Here we consider the multidimensional version of the Kuramoto model, as proposed in <cit.>. 
We first derive the higher order corrections in D dimensions that generalize the corresponding interactions in D=2. Then, we show that,
generally, the critical coupling for synchronization changes when three-body interactions are taken into account, making the transition 
from disordered to ordered states more difficult. Only in two dimensions, corresponding to the original Kuramoto model, this dependence 
disappears. Above the modified critical coupling third and fourth order terms combine to facilitate synchronization. These
features are demonstrated using mean field arguments and are verified with numerical simulations, since the Ott-Antonsen ansatz 
does not extend to dimensions larger than 2.






§ HIGHER ORDER INTERACTIONS IN THE MULTIDIMENSIONAL KURAMOTO MODEL


In the Kuramoto model, a set of N oscillators, represented only by their phases θ_i are coupled according to the equations

    θ̇_i    =   ω_i + k_1/N∑_j=1^N sin(θ_j-θ_i)


where  ω_i are their natural frequencies, selected from a symmetric distribution g(ω). Three and four body interactions 
are introduced as


    θ̇_i    =   ω_i + k_1/N∑_j=1^N sin(θ_j-θ_i) + 
    		k_2/N^2∑_j,k=1^N sin(2θ_j-θ_k -θ_i)
       +   k_3/N^3∑_j,k,m=1^N sin(θ_j+θ_k -θ_m-θ_i)


where i = 1, ..., N and k_1, k_2 and k_3 are the coupling constants for pairwise (1-simplex), triplets (2-simplex) and quadruplets (3-simplex) interactions.  
We note that the specific form of the three and four body interactions in Eq.(<ref>) follows the choice in <cit.>, although other, more
symmetric forms, have also been considered  <cit.> .

Using the definition of the order parameter


    z = 1/N∑_j e^i θ_j = r e^i ψ


we obtain


    1/N∑_j e^i (θ_j - θ_i) = r e^i (ψ - θ_i),



    1/N∑_k e^i (θ_k + θ_i - 2θ_j) = r e^i (ψ + θ_i - 2θ_j)


and


    z^2 z^* e^-iθ_i= 1/N^3∑_j,k,m e^i (θ_j +θ_k -θ_m- θ_i) = r^3 e^i (ψ - θ_i).




Taking the imaginary part of these equations we can rewrite the Kuramoto model as


    θ̇_i = ω_i + k_1 r sin(ψ - θ_i) +  k_2 r 1/N∑_j=1^N sin(2θ_j-θ_i-ψ) + k_3 r^3 sin(ψ-θ_i).



Following <cit.> we also define a second order parameter


    z_2 = 1/N∑_j e^2 i θ_j = p_2 e^i ξ_2


and the related vector


    p⃗_2 = p_2 (cosξ_2,sinξ_2)


so that (<ref>) becomes


    θ̇_i = ω_i + k_1 r sin(ψ - θ_i) +  k_2 r p_2 sin(ξ_2-θ_i-ψ) + k_3 r^3 sin(ψ-θ_i).


We use p_2 and ξ_2 for the module and phase of z_2 as the vector r⃗_2 will be reserved for a later definition.

Eq.(<ref>) can be generalized to more dimensions if we can write it in terms of 
the unit vectors σ⃗_⃗i⃗ = (cosθ_i,sinθ_i) ≡ (σ_ix,σ_iy) 
<cit.>.  Computing σ̇_ix = - θ̇_i sinθ_i and using Eq.(<ref>) we find


    σ̇_ix   =    - sinθ_i[ ω_i + k_1 r sin(ψ - θ_i) +  k_2 r p_2 sin(ξ_2-θ_i-ψ) + k_3 r^3 sin(ψ-θ_i)]  
       ≡    -ω_i σ_iy + k_1 T_1x + k_2 T_2x + k_3 T_3x.


The contributions of the interaction terms are computed as follows:


    T_1x    =     - r sinθ_i sin(ψ - θ_i) 
       =    - r  sinψsinθ_i cosθ_i + r sinθ_i^2 cosψ
       =     r cosψ - r [ sinψsinθ_i  +  cosθ_i cosψ] cosθ_i 
       =    r_x -(r⃗·σ⃗_i) σ_ix


where we defined the vector


    r⃗ = (rcosψ,rsinψ)


representing the center of mass of the system.

The third order term is T_2x =  - r p_2 sinθ_i sin(ξ_2-ψ - θ_i), 
which becomes identical to T_1x if we replace r → r p_2 and ψ→ξ_2 - ψ. Therefore we find


    T_2x   =    r_2x -(r⃗_2 ·σ⃗_i) σ_ix


where

    r⃗_2 = r p_2(cos(ξ_2-ψ), sin(ξ_2-ψ)) = (r_x p_2x+r_y p_2y, r_x p_2y-r_y p_2x).


Similarly we obtain


    T_3x   =    r_3x -(r⃗_3 ·σ⃗_i) σ_ix


where

    r⃗_3 = r^3 (cosψ, sinψ) = r^2 r⃗


with r_3 = r^3 and ψ_3 = ψ. 

After doing a similar calculation for the σ_iy we can write the dynamical equations in vector form as 


    d σ⃗_⃗i⃗/d t = 𝐖_i σ⃗_⃗i⃗ + ∑_j=1^3 k_j [r⃗_j -(σ⃗_i ·r⃗_j) σ⃗_i]


where we identify r⃗_1 ≡r⃗ and  𝐖_i is the anti-symmetric matrix 


    𝐖_i = ( 
    	[    0 -ω_i;  ω_i    0 ]).



In order to extend the model to higher dimensions it is essential that all vectors r⃗_j can be written in terms of the unit vectors σ⃗_j.
Clearly


    r⃗ = 1/N∑_k=1^N σ⃗_k


and


    r⃗_3 = r_2^2 r⃗.


To write r⃗_2 in terms of σ⃗'s we first realize that


    p_2x = 1/N∑_j cos 2θ_j = 1/N∑_j [σ_jx^2 - σ_jy^2] = 2/N∑_j σ_jx^2 - 1 
    	    ≡  2⟨σ_x^2 ⟩ - 1  = 1 - 2 ⟨σ_y^2 ⟩


and


    p_2y = 1/N∑_j sin 2θ_j = 2/N∑_j  σ_jxσ_jy≡  2⟨σ_x σ_y ⟩ .


Next, for any vector σ⃗_j  we define the diadic matrix Λ_j = σ⃗_j r⃗^T and the
vector


    1/N∑_jΛ_jσ⃗_j = 
       1/N∑_j(
    	[ σ_jx r_x σ_jx r_y; σ_jy r_x σ_jy r_y;          ])
    	(
    	[ σ_jx; σ_jy ]) =
    	1/N∑_j(
    	[ r_xσ_jx^2 + r_yσ_jxσ_jy; r_xσ_jxσ_jy + r_yσ_jy^2 ])


which can also be written as


    (
    	[ r_x⟨σ_x^2  ⟩ + r_y⟨σ_xσ_y⟩; r_x⟨σ_xσ_y⟩ + r_y⟨σ_y^2  ⟩ ])	= 1/2(
    	[ r_x (1+p_2x) + r_y  p_2y;  r_x p_2y + r_y (1-p_2x) ])	 =1/2 (r⃗ + r⃗_2)


where we used Eqs.(<ref>) and (<ref>). Therefore


    1/N∑_jΛ_jσ⃗_j  
    	= 1/N∑_j (r⃗·σ⃗_j) σ⃗_j = 1/2(r⃗+r⃗_2)
	

or


    r⃗_2 =  2/N∑_j (r⃗·σ⃗_j) σ⃗_j  - r⃗.
	



Eq.(<ref>) can now be extended to higher dimensions by simply considering unit vectors σ⃗_i in D-dimensions, rotating on the surface of the corresponding (D-1) unit sphere 
<cit.>. Particles are now represented by D-1 spherical angles, generalizing the single phase θ_i of the original model. The matrices 𝐖_i become D × D anti-symmetric matrices containing the D(D-1)/2 natural frequencies of each oscillator. 
The equations up to fourth order can be summarized as follows:


    d σ⃗_⃗i⃗/d t = 𝐖_i σ⃗_⃗i⃗ + ∑_j=1^3 k_j [r⃗_j -(σ⃗_i ·r⃗_j) σ⃗_i]


with


    r⃗_1 = r⃗ = 1/N∑_k=1^N σ⃗_k,



    r⃗_2 =  2/N∑_j (r⃗·σ⃗_j) σ⃗_j  - r⃗
	

and


    r⃗_3 = r_2^2 r⃗.







§ PHASE TRANSITION ON THE SPHERE AND HIGHER DIMENSIONS


The transition to synchronization in D=3  is discontinuous in the case of all-to-all pair interactions, 
with r=0 if k_1 < 0 and r(0^+) =1/2 <cit.>. The effects of k_2 and k_3 can be estimated with a 
mean-field calculation of the vectors r⃗_2 and r⃗_3.
Using (<ref>)  we write r⃗_2 = ρ⃗_2 - r⃗ with


    ρ⃗_2 =  2/N∑_j (r⃗·σ⃗_j) σ⃗_j 
    	                        =  (1/N∑_j σ⃗_j σ⃗_j^T ) 2 r⃗
	

where the dyadic matrices are


    σ⃗_j σ⃗_j^T = (
    	[   σ_jx^2 σ_jxσ_jy σ_jxσ_jz; σ_jyσ_jx   σ_jy^2 σ_jyσ_jz; σ_jzσ_jx σ_jzσ_jy   σ_jz^2 ]).


For asynchronous states we can assume that the components of the σ⃗_j are uncorrelated and that the average of  
σ_jx^2, σ_jy^2 and 
σ_jz^2  is  1/3, since they sum to one. Thus, the average in Eq.(<ref>) is 1/3 times de identity matrix and
ρ⃗_2 = 2 r⃗ /3 and r⃗_2 = ρ⃗_2 - r⃗ = - r⃗/3.  Therefore, the net effect of the third order correction 
is to shift k_1 → k_12≡  k_1 - k_2/3. The first order phase transition that would occur at 
k_1=0 is now shifted to the line k_12=0, or to


    k_1 = k_2/3.


To the right of this critical line the σ⃗_j's start to correlate until they all point in the direction of r⃗. In this region
we can approximate σ⃗_j ≈r⃗ and get


    ρ⃗_2 =   (1/N∑_j σ⃗_j σ⃗_j^T ) 2 r⃗≈ 2 r^2 r⃗
	

so that r⃗_2 = r^2 r⃗.

Since r⃗_3 = r^2 r⃗, it does not contribute when r=0 and, therefore, does not affect the critical line k_2=3k_1.
In the synchronous region it contributes positively, facilitating synchronization. 
Together with the third order term, that acts similarly, the effect of the higher order terms is to shift the curves 
r=r(k_1) to the left by making k_1 → k_1eff = k_1 - r^2(k_2 + k_3) if k_1eff > k_2/3, making synchronization 
more effective (larger values of r) than the original model with pair interactions. If we call r(k_1,k_2,k_3) the 
equilibrium value of the order parameter for fixed values of the coupling constants, then


    r(k_1,k_2,k_3) = {[                       0       f     k_1 < k_2/3; r(k_1+r^2(k_2+k_3),0,0)     f     k_1 > k_2/3 .;                         ].


Conversely, all curves collapse onto r(k_1,0,0) if shifted to the left according to:


    r(k_1-r^2(k_2+k_3),k_2,k_3) = r(k_1,0,0)


in the region where r ≠ 0.

Using similar arguments we see that the analogue of Eq.(<ref>) in D dimensions results in (1/D) 𝐈 in the disordered region,
where 𝐈 is the D × D identity matrix. This implies that ρ⃗_2 = (2/D) r⃗ and 
r⃗_2 = -[(D-2)/D] r⃗.  In D=2 there is no displacement in the critical point <cit.>, 
but it shifts the transition in all other dimensions. The critical coupling goes from k_1c to k_1c + k_2 (D-2)/D.

In odd dimensions we expect the transition to be discontinuous, with the critical line shifted. In even dimensions
the transition remains continuous, but ever steeper as k_2 and k_3 increases.









§ NUMERICAL RESULTS


In this section we show numerical simulations of Eqs.(<ref>) for D=3 and D=4. We consider N=2000
particles and choose the matrices of natural frequency 𝐖_i such that all entries are Gaussian distributed 
around ω=0 with unit mean square deviation Δ=1. See <cit.> for phase diagrams with pair 
interactions and different values of Δ and ω. Initial conditions are uniformly distributed over the corresponding 3D
or 4D sphere and the equations are integrated up to t=2000. The order parameter r is computed as the average
over the last 500 time units. 

Figures <ref> and <ref> show simulations for D=3. 
Fig. <ref>(a) shows heatmaps of the order parameter r as a function of k_1 and k_2 for k_3=0 whereas
Fig. <ref>(b) shows results for k_3=2. In both panels the critical line from disorder to synchronization follows k_2 = 3 k_1 and is
independent of k_3. The phase transition is always discontinuous, as in the case with k_2=k_3=0, but the
size of the discontinuity from r=0 to r ≠ 0 increases from 1/2 at k_2=k_3=0
to about 0.9 when k_2=4 and k_3=0. For k_3 > 0 the jump increases even faster.

Fig. <ref>(a) shows r as a function of k_1 for some values of k_2 and k_3. As k_2 increases the critical
point moves towards larger values and the jump at the critical point increases. Fig. <ref>(b) shows the map 
r(k_1,k_2,k_3)  → r(k_1-r^2(k_2+k_3),k_2,k_3), which shifts the curves to the left (when r ≠ 0). All
curves fall approximately on the same curve r(k_1,0,0) as predicted by Eq.(<ref>).

Similar plots for D=4 are shown in Figs. <ref> and <ref>. In this case the phase transition is continuous when
k_2=k_3=0, with k_1c≈ 1.5. Again, the value of k_3 does not affect the shift in the critical coupling, which is
moved from k_1c to approximately k_1c + k_2/2, as shown in Fig. <ref>(a) for k_3=0 and  Fig. <ref>(b) for k_3=2. 
The approximation is still accurate, even is the neighborhood of the critical line. To the right of the critical line the system achieves 
higher values of synchronization as compared to the case k_2=k_3=0, leading to a discontinuity that increases as the intensity of 
higher order interactions increases, as shown in Fig. <ref>(a). Fig. <ref>(b) shows that shifting the curves according to Eq.(<ref>) 
also makes them all collapse into r(k_1,0,0).
 










§ CONCLUSIONS


The effect of higher order interactions on the dynamics of complex systems has been a hot topic of research in 
the past years. Although this type of many body action depends in a great measure on the topology of 
connections between agents, mean field estimates can serve as a guide to understand the changes expected
when these interactions are included. Here we considered the multidimensional version of the Kuramoto model
as introduced in <cit.> and derived the functional form of third and fourth order corrections that
are natural extensions of their counterparts in the usual, D=2, Kuramoto model. 

We have shown that the third order correction has a special role in displacing the critical coupling constant towards higher 
values in all dimensions, except D=2. After the critical points, higher order terms tend to facilitate synchronization in
a simple way, shifting the order parameter r(k_1,k_2,k_3) to r(k_1+r^2(k_2+k_3),0,0).
The shift causes the phase transition to become discontinuous also in even dimensions, where the original model
predicts continuous phase transitions.

We have not investigated the possibility of bi-stable regions, which do occur in D=2 <cit.>. That
would require scanning the space of initial conditions for each value of the coupling parameters, since no version of
the Ott-Antonsen ansatz is available for D>2. This is an interesting next step in this research.




	This work was partly supported by FAPESP, grant 2021/14335-0 (ICTP‐SAIFR) and CNPq, grant 301082/2019‐7. 








22
	
	
	urlstyle[Battiston et al.(2020)Battiston, Cencetti, Iacopini, Latora, Lucas,
	Patania, Young, and Petri]battiston2020networks
	Federico Battiston, Giulia Cencetti, Iacopo Iacopini, Vito Latora, Maxime
	Lucas, Alice Patania, Jean-Gabriel Young, and Giovanni Petri.
	Networks beyond pairwise interactions: Structure and dynamics.
	Physics Reports, 874:0 1–92, 2020.
	
	[Petri et al.(2014)Petri, Expert, Turkheimer, Carhart-Harris, Nutt,
	Hellyer, and Vaccarino]petri2014homological
	Giovanni Petri, Paul Expert, Federico Turkheimer, Robin Carhart-Harris, David
	Nutt, Peter J Hellyer, and Francesco Vaccarino.
	Homological scaffolds of brain functional networks.
	Journal of The Royal Society Interface, 110
	(101):0 20140873, 2014.
	
	[Sizemore et al.(2018)Sizemore, Giusti, Kahn, Vettel, Betzel, and
	Bassett]sizemore2018cliques
	Ann E Sizemore, Chad Giusti, Ari Kahn, Jean M Vettel, Richard F Betzel, and
	Danielle S Bassett.
	Cliques and cavities in the human connectome.
	Journal of computational neuroscience, 44:0 115–145,
	2018.
	
	[Ganmor et al.(2011)Ganmor, Segev, and Schneidman]ganmor2011sparse
	Elad Ganmor, Ronen Segev, and Elad Schneidman.
	Sparse low-order interaction network underlies a highly correlated
	and learnable neural population code.
	Proceedings of the National Academy of sciences, 1080
	(23):0 9679–9684, 2011.
	
	[Grilli et al.(2017)Grilli, Barabás, Michalska-Smith, and
	Allesina]grilli2017higher
	Jacopo Grilli, György Barabás, Matthew J Michalska-Smith, and Stefano
	Allesina.
	Higher-order interactions stabilize dynamics in competitive network
	models.
	Nature, 5480 (7666):0 210–213, 2017.
	
	[Sanchez-Gorostiaga et al.(2019)Sanchez-Gorostiaga, Bajić, Osborne,
	Poyatos, and Sanchez]sanchez2019high
	Alicia Sanchez-Gorostiaga, Djordje Bajić, Melisa L Osborne, Juan F Poyatos,
	and Alvaro Sanchez.
	High-order interactions distort the functional landscape of microbial
	consortia.
	PLoS Biology, 170 (12):0 e3000550, 2019.
	
	[Benson et al.(2016)Benson, Gleich, and Leskovec]benson2016higher
	Austin R Benson, David F Gleich, and Jure Leskovec.
	Higher-order organization of complex networks.
	Science, 3530 (6295):0 163–166, 2016.
	
	[de Arruda et al.(2020)de Arruda, Petri, and Moreno]de2020social
	Guilherme Ferraz de Arruda, Giovanni Petri, and Yamir Moreno.
	Social contagion models on hypergraphs.
	Physical Review Research, 20 (2):0 023032,
	2020.
	
	[Iacopini et al.(2019)Iacopini, Petri, Barrat, and
	Latora]iacopini2019simplicial
	Iacopo Iacopini, Giovanni Petri, Alain Barrat, and Vito Latora.
	Simplicial models of social contagion.
	Nature communications, 100 (1):0 2485, 2019.
	
	[Jhun et al.(2019)Jhun, Jo, and Kahng]jhun2019simplicial
	Bukyoung Jhun, Minjae Jo, and B Kahng.
	Simplicial sis model in scale-free uniform hypergraph.
	Journal of Statistical Mechanics: Theory and Experiment,
	20190 (12):0 123207, 2019.
	
	[Vega et al.(2004)Vega, Vázquez-Prada, and
	Pacheco]vega2004fitness
	Yamir Moreno Vega, Miguel Vázquez-Prada, and Amalio F Pacheco.
	Fitness for synchronization of network motifs.
	Physica A: Statistical Mechanics and its Applications,
	343:0 279–287, 2004.
	
	[Berec(2016)]berec2016chimera
	Vesna Berec.
	Chimera state and route to explosive synchronization.
	Chaos, Solitons & Fractals, 86:0 75–81, 2016.
	
	[Skardal and Arenas(2019)]skardal2019abrupt
	Per Sebastian Skardal and Alex Arenas.
	Abrupt desynchronization and extensive multistability in globally
	coupled oscillator simplexes.
	Physical review letters, 1220 (24):0 248301,
	2019.
	
	[Skardal and Arenas(2020)]skardal2020higher
	Per Sebastian Skardal and Alex Arenas.
	Higher order interactions in complex networks of phase oscillators
	promote abrupt synchronization switching.
	Communications Physics, 30 (1):0 218, 2020.
	
	[Moyal et al.(2024)Moyal, Rajwani, Dutta, and Jalan]sarika2024
	Bhuwan Moyal, Priyanka Rajwani, Subhasanket Dutta, and Sarika Jalan.
	Rotating clusters in phase-lagged kuramoto oscillators with
	higher-order interactions.
	Phys. Rev. E, 109:0 034211, Mar 2024.
	10.1103/PhysRevE.109.034211.
	URL <https://link.aps.org/doi/10.1103/PhysRevE.109.034211>.
	
	[Biswas and Gupta(2024)]biswas2024symmetry
	Dhrubajyoti Biswas and Sayan Gupta.
	Symmetry-breaking higher-order interactions in coupled phase
	oscillators.
	Chaos, Solitons & Fractals, 181:0 114721, 2024.
	
	[Ott and Antonsen(2008)]Ott2008
	Edward Ott and Thomas M. Antonsen.
	Low dimensional behavior of large systems of globally coupled
		oscillators.
	Chaos, 180 (3):0 1–6, 2008.
	ISSN 10541500.
	10.1063/1.2930766.
	
	[Chandra et al.(2019)Chandra, Girvan, and Ott]2019continuous
	Sarthak Chandra, Michelle Girvan, and Edward Ott.
	Continuous versus discontinuous transitions in the d-dimensional
	generalized kuramoto model: Odd d is different.
	Physical Review X, 90 (1):0 011002, 2019.
	
	[Lipton et al.(2021)Lipton, Mirollo, and Strogatz]lipton2021kuramoto
	Max Lipton, Renato Mirollo, and Steven H Strogatz.
	The kuramoto model on a sphere: Explaining its low-dimensional
	dynamics with group theory and hyperbolic geometry.
	Chaos: An Interdisciplinary Journal of Nonlinear Science,
	310 (9):0 093113, 2021.
	
	[Ashwin and Rodrigues(2016)]ashwin2016hopf
	Peter Ashwin and Ana Rodrigues.
	Hopf normal form with sn symmetry and reduction to systems of
	nonlinearly coupled phase oscillators.
	Physica D: Nonlinear Phenomena, 325:0 14–24, 2016.
	
	[Barioni and de Aguiar(2021)]barioni2021complexity
	Ana Elisa D Barioni and Marcus AM de Aguiar.
	Complexity reduction in the 3d kuramoto model.
	Chaos, Solitons & Fractals, 149:0 111090, 2021.
	
	[Fariello and de Aguiar(2024)]Fariello2024a
	Ricardo Fariello and Marcus AM de Aguiar.
	Exploring the phase diagrams of multidimensional kuramoto models.
	Chaos, Solitons & Fractals, 179:0 114431, 2024.
	



	
