Collective behavior of coupled multiplicative processes with stochastic resetting

A dynamical variable driven by the combination of a deterministic multiplicative process with stochastic reset events develops, at long times, a stationary power-law distribution. Here, we analyze how such distribution changes when several variables of the same kind interact with each other through diffusion-like coupling. While for weak coupling the variables are still distributed following power-law functions, their distributions are severely distorted as interactions become stronger, with sudden appearance of cutoffs and divergent singularities. We explore these effects both analytically and numerically, for coupled ensembles of identical and non-identical variables. The most relevant consequences of ensemble heterogeneity are assessed, and preliminary results for spatially distributed ensembles are presented.


Introduction
Besides underlying a variety of physical and chemical phenomena, which range from nuclear fission, fragmentation and coalescence to combustion and autocatalytic reactions [1][2][3], multiplicative processes play a key role in the dynamics of many systems of biological and socioeconomic nature [4][5][6][7]. The growth of organisms and populations, the expansion of urban areas, the production of goods, and the creation of wealth, among others, are driven by mechanisms with the capacity of promoting and reinforcing themselves [8]. If it were by them alone, such mechanisms would give origin to an exponential increase of the relevant quantities-number of cells or individuals, commodities, assets, money-leading to their divergence in the long term. As we all know well, however, this growth cannot last forever. In the real world, a vast host of factors, usually related to the finiteness of environmental resources, puts limits on the increase of abundances of any kind.
Mathematical modeling implements limitations to exponential growth in various ways, depending on the kind of system under study. In stochastic resetting [9][10][11][12], exponential growth is abruptly interrupted at randomly recurring times and the growing variable is taken back close to its initial value. It has been proposed as a limiting mechanism which emulates catastrophic events such as natural disasters, mass extinction, social collapses, or market-bubble bursting-the 'black swans' of economic theory [13,14]. A remarkable outcome of the combination of multiplicative processes and stochastic resetting is that, under very general conditions, it gives rise to a power-law profile in the long-time asymptotic distribution of the relevant variable. This is a desirable property, in view of the ubiquity of such kind of distributions in taxonomic abundances [15], family sizes [16], city populations [17], and wealth sharing [18], among many other examples. The occurrence of power-law distributions, in fact, is characteristic of a wide class of stochastic multiplicative models [19,20].
In recent years, the effects of stochastic resetting have been studied in connection with a variety of physical processes [12], not restricted to multiplicative growth. Its combination with both isotropic and directed diffusion [21][22][23], active-particle dynamics [24,25], and transport on networks [26] has received particular attention. Applications have covered from hydrologic phenomena [27] to RNA dynamics [28] and antiviral therapy design [29]. A thermodynamical approach to stochastic resetting has also been proposed [30]. It has been recognized, however, that coupling between dynamical elements subject to this kind of processes has been dealt with only preliminarily, in connection with interacting Brownian particles [12,31].
In this paper, we study the effects of coupling in the collective behavior of an ensemble of interacting elements whose individual dynamics is characterized by exponential growth punctuated by stochastic resetting. To put our model in a stylized socioeconomic context, we consider a population of agents, each of them possessing a variable amount of resources driven by the combination of a deterministic multiplicative process and random reset events. Stochastic resetting affects each agent independently, entailing a sudden loss of a large part of the individual resources. Thus, it does not represent global catastrophic circumstances, but rather localized events which severely disturb the state of one agent at each occurrence-such as, for instance, the failure of a personal business venture. In turn, coupling between agents acts as a form of tax levy and revenue distribution, taking a part of all resources and redistributing them over the whole population. We show that, while the stationary resource distribution of uncoupled agents follows a pure power-law function, this rather simple interaction mechanism can induce severe changes in the distribution profiles, with cutoffs and divergences at finite values of the resources. We disclose these effects in populations of identical agents and then analyze how they manifest themselves in heterogeneous ensembles. Finally, we present preliminary results on the collective behavior of spatially distributed populations.

Globally coupled ensembles and the dynamics of an uncoupled agent
We consider a population of N agents, each of them described by a one-dimensional variable x i (t) > 0, which represents the agent's resources. As advanced in the introduction, resources change due to three mechanisms: • Exponential growth at rate λ i > 0; • Stochastic resetting to x i = u i > 0, with frequency q i ; • Resource redistribution. This form of diffusion-like coupling between agents is implemented by taking a fraction of resources from each agent at rate a i , so that a total amount T(t) = i a i x i (t) is collected per time unit. This amount is then redistributed, assigning a fraction b i to each agent, with i b i = 1. The collection rates a i give a direct measure of the strength of coupling. The variable x i (t) satisfies the stochastic equatioṅ for i = 1, . . . , N, where the three terms on the right-hand side respectively represent the mechanisms listed above. In the second term, is a Poisson process, with δ(t) being Dirac's delta. It describes reset events at times t i,k (i = 1, . . . , N, k = 1, 2, 3, . . .), occurring with frequency q i [27]. The prefactor (u i − x i ) insures that resetting occurs from the current value of x i to u i . The third term in (1) globally couples the evolution of different agents, through the mean-field-like quantity T(t).
Considering, as in any mean-field approach, that T(t) is a prefixed quantity whose value has to be self-consistently determined once the solution for each x i (t) has been found, we can write the Chapman-Kolmogorov equation for the probability distribution where is the velocity of probability flow along x i . If, as we assume from now on, v i f i → 0 for both x i → 0 and x i → ∞, and for all t 0, (3) preserves the total probability, and the integral in the last term on its right-hand side equals unity. The general procedure for solving (3) is to find a piecewise solution in the intervals 0 < x i < u i and u i < x i , where the term proportional to δ(x i − u i ) vanishes, and then match the solutions at each side of u i using the delta function to define a boundary condition. The delta function, in fact, determines a finite discontinuity in the whole solution, at all times. If agents are decoupled from each other, with a i = 0 for all i, (3) can be fully solved for any initial condition [11]. In this case, the values of λ i and u i could be arbitrarily fixed without generality loss, by an appropriate rescaling of time and resources. In our analysis, however, we explicitly keep the two parameters to encompass the cases where they vary from agent to agent. Here, we are interested in the long-time asymptotic solution, which reads for x i u i , and 0 otherwise. Indeed, the probability vanishes for 0 < x i < u i due to the combined effect of reset events, which sustainedly transfer probability to u i , and a velocity of probability flow, v i = λ i x i , which is positive everywhere. Above u i , on the other hand, f st i (x i ) is a power-law stationary distribution with an exponent controlled by the parameters q i and λ i . As expected, its tail becomes fatter when the multiplication factor λ i grows and/or when the reset frequency q i decreases. The moments of the distribution f st i (x i ), are well defined up to a critical exponent γ c = q i /λ i . For γ > γ c , they diverge. Note, in particular, that the stationary mean value x i st (γ = 1) is finite and positive if λ i < q i . In the following, as needed to obtain our analytical results, we assume that this condition holds. The nonexistence of low-order stationary moments is directly related to the absence of self-averaging properties, which jeopardizes the statistical predictability of the system [33,34].
The left panel of figure 1 shows a typical numerical realization of the solution to (1) for an agent in a decoupled ensemble, with q i = 2.5, λ i = 1, and u i = x i (0) = 1. The integration step is δt = 10 −3 . The evolution consists of periods of variable length during which x i grows exponentially, punctuated by instantaneous events at which x i is reset to u i . Dots in the right panel show a normalized histogram of some 10 7 successive values of x i along the same realization, and the straight line represents the corresponding Chapman-Kolmogorov prediction (6).
In the following sections, we first study the effect of coupling on the stationary distribution of resources in a homogeneous ensemble, where all agents have the same individual parameters. Next, we generalize these results to the non-homogeneous case.

Homogeneous ensembles
In homogeneous ensembles, all agents have the same individual parameters, so that we can drop the index i of all quantities involved in (3) and (4). Moreover, the redistribution coefficient is b i ≡ b = N −1 for all agents, and the last term on the right-hand side of (4) becomes withx(t) the arithmetic mean of resources over the ensemble. In order to find a stationary solution to the Chapman-Kolmogorov equation (3) for homogeneous ensembles, which is expected to stand for the long-time asymptotic distribution of resources of any single agent, two assumptions are necessary. First, we need thatx(t) becomes independent of time as t → ∞-discarding, for instance, sustained drift or oscillations. Second, the asymptotic value ofx, which is an average taken over the ensemble, must coincide with the average of resources over the stationary distribution for a single agent, x st as defined in (7) for γ = 1, in order to pose the self-consistency relation from which the value ofx is derived. This ergodic hypothesis implies in turn that we neglect stochastic fluctuations inx, which amounts to assume that we are working with an infinitely large ensemble. Note that the validity of these assumptions has to be assessed by studying the full temporal evolution of the distribution f(x, t). Below, we provide numerical evidence supporting it.
In the homogeneous case, the general stationary solution to (3) for any form of v(x) and x = u reads where the arbitrary constant K can vary between the various domains of x where the solution is being computed. In any specific solution, K has to be tuned for each interval in order to comply with the boundary conditions imposed by the delta-like term in (3), by the form of v(x), and by the overall normalization of Inspection of the functional form of v(x) makes readily clear that it is necessary to consider two separate cases. If coupling is sufficiently weak, a < λ, whatever the value ofx > 0, v(x) is positive. The situation is similar to the case of a single uncoupled agent (section 2) and the problem has to be first separately solved in the regions 0 < x < u and x > u. The two solutions must then be matched at x = u. On the other hand, for strong coupling, a > λ-and assuming that, as expected, x > u-the velocity changes sign at where, however, f st (x) is expected to be continuous. In this case, the solution has three pieces to be matched, in the domains 0 < x < u, u < x < x 0 , and x > x 0 . The solutions for the two cases are explicitly given in the following.

Weak coupling: a < λ
Much as in the case of an uncoupled agent, the combination of reset events and an everywhere positive velocity This distribution is correctly normalized and complies with the boundary condition (5). Note that (12) is a straightforward generalization of (6) for a = 0. The self-consistency relation (9) holds if cf (7). The stationary averagex, which is well defined if λ < q (see section 2), is independent of the collection rate a, as expected from the fact that the total amount of resources is not affected by the redistribution process.
For large x, f st (x) exhibits a power-law tail whose exponent strongly depends on a. As a increases and resource redistribution becomes more effective, the steepness of f st (x) grows, indicating a less dispersed distribution of individual resources around the (a-independent) mean valuex. Accordingly, the critical order for the convergence of the moments defined in (7) increases as is not a pure power-law function, but develops a shoulder just above x = u. The width of this shoulder, W, can be defined as the distance between u and the point of intersection of the extrapolation of the large-x power-law tail and f st (u). This gives As an example of this case, the right-hand panel in the uppermost row of figure 2 shows the distribution f st (x) for q = 2.5, λ = 1, u = 1, and a = 0.5. The curve corresponds to the analytical prediction (12) and dots stand for a normalized histogram built from a numerical realization for 1000 coupled agents. Excellent agreement between both data sets validate the assumptions made in the analytical calculation. The dotted line In all cases, q = 2.5, λ = 1, and u = 1. Right: analytical (curves) and numerical (dots) results for the stationary resource distribution f st (x). Numerical integration of (1) was performed as explained in the main text for figure 1 (section 2). In the panel to the right of the uppermost row, the dotted straight line illustrates the construction used to define the width W of the shoulder in the curve, equation (14).
illustrates the geometrical construction that defines the width W of equation (14). In the panel to the left, we show a portion of the time evolution of resources for a single agent in the same numerical realization.
It is worthwhile noticing that the form of f st (x) in (12) has a well-defined limit for a → λ: In the limit, thus, there is a change in the functional form of the resource distribution, with the power-law tail replaced by a purely exponential decay in x. This implies, in particular, that the moments x γ are defined for any order γ. The change from power-law to exponential decay anticipates a significant modification in f st (x) as we switch from weak to strong coupling.

Strong coupling: a > λ
When the effect of redistribution overcomes the exponential growth of individual resources (a > λ), as advanced above, the velocity of probability flow is negative for x > x 0 , with x 0 given by (11). Since there is no mechanism that transports probability from small values of x toward x > x 0 , the long-time asymptotic probability distribution vanishes in that interval. Therefore, we have f st (x) = 0 for both x < u and x > x 0 .
In the intervening interval, u x x 0 , f st (x) is still given by (10) or, equivalently, by (12), which we now conveniently rewrite as For a just above λ, the exponent in (16) is a large positive number. This implies that the probability distribution monotonically decreases for u < x < x 0 , and vanishes at x = x 0 . In particular, f st (x) is continuous at that point. However, the derivative of f st (x) at x 0 switches abruptly from 0 to −∞ when a crosses from below the critical value a 1 = λ + q/2, and is finite and negative for a = a 1 . This implies a sudden change in the profile of the distribution just below x 0 : for a < a 1 , it matches smoothly with the vanishing solution at x > x 0 , while for a a 1 there is a sharp discontinuity in its slope. The second and third rows in figure 2 respectively illustrate the two cases for a = 1.8 and 3.2, with q = 2.5, λ = 1, and u = 1. Note, in the respective panels to the left, how the fast growth of x i (t) after each reset event is soon moderated by the action of resource redistribution, becoming increasingly slower as time goes on. This effect, which turns out to be consistent with the cutoff at x 0 in the stationary distribution f st (x), is stronger for larger values of a.
Abrupt transitions are not yet exhausted, though. At a second critical value of the collection rate, a 2 = λ + q, the exponent in (16) crosses 0 and becomes negative. For a = a 2 , f st (x) is constant in the interval u < x < x 0 . For a > a 2 , meanwhile, the distribution is a monotonically increasing function which diverges at x = x 0 . The divergence, however, is integrable for all a > a 2 , as the exponent is always greater than −1. An example of this case is illustrated for a = 4 in the lowermost row of figure 2, with all the other parameters equal as above.
As a keeps growing and the exponent in (16) approaches −1, an increasingly large fraction of the total probability concentrates just below the divergence at x = x 0 . At the same time, x 0 becomes closer and closer tō x; cf (11). In fact, when the redistribution rate increases and coupling neatly dominates the dynamics, we expect that individual resources accumulate around their ensemble average. In the limit of large a, their distribution f st (x) becomes a (one-sided) approximant to Dirac's delta.

Non-homogeneous ensembles
We now turn back the attention to the Chapman-Kolmogorov equation (3) in the case where all parameters are generally different from agent to agent. Following mutatis mutandis the same procedure as for homogeneous ensembles, the stationary solution for the resource distribution of agent i reads cf (12). Here, we have assumed that T(t) = i a i x i (t) is a constant T for asymptotically long times. Along the same lines as for the homogeneous case, we estimate this constant by using the average of each x i over the withx i computed as in (9). In the product b i T which appears in (17) Replacing into (18) yields a self-consistency equation for T, whose solution is This expression makes it possible to compute T in terms of the individual parameters of agents, thus completing the solution to our problem.
For a sufficiently large ensemble, the sums in (20) can be estimated by means of integral approximations, assuming that the individual parameters are drawn from a distribution g (q, λ, u, a, b). Namely, g(q, λ, u, a, b)F(q, λ, u, a, b), (21) where F (q, λ, u, a, b) stands for the summands. Whether the integral in (21) can be explicitly performed, of course, depends crucially on the form of g (q, λ, u, a, b), which encompasses all the information on the distribution of parameters and their mutual correlations. Since the simultaneous consideration of all these statistical properties is impractical, in the following we limit ourselves to the analysis of two representative cases in which T can be given an explicit expression.

Distribution in a i
We have seen in section 3 that, for a homogeneous ensemble, different values of the collection rate a give origin to qualitatively diverse asymptotic resource distributions. It thus makes sense, in a non-homogeneous ensemble, to consider first the effect of a distribution in a i , keeping the other parameters equal for all agents. Within a continuous approximation as in (21), and supposing that a i is drawn from a distribution g(a), we get From (19), the average of individual resources is This result shows that heterogeneity in the parameter a i has a disparate impact on the resources of different agents. In fact, the average of individual resources is no more independent of a i as it was in homogeneous ensembles [cf (13)], but monotonically decreases as a i grows. This can be interpreted in terms of the fact that the larger a i the greater the amount of resources taken from agent i. Meanwhile, upon redistribution, all agents get the same amount. It is interesting to note that, as expected from the combined fact that redistribution does not affect the total amount of resources and that agents are identical in all parameters except a i , the mean value ofx i over the whole ensemble coincides with that of the homogeneous case: In other words, although the average of individual resourcesx i depends on the specific value of a i , total resources are independent of the form in which a i is distributed.
In the case where g(a) is a narrow symmetric distribution centered at a certain valueā, the integral I a can be approximated using Sommerfeld's expansion [35], where σ a is the standard deviation of a and μ (2n) a are higher-order even moments of g(a). To the first significant correction in the expansion, (23) yields We see that at the center of the distribution, i.e. for a i =ā,x i is smaller than the average value expected for the same parameters. Only for a i <ā − σ 2 a /(q − λ +ā) is the individual average of resources larger than the result (13) for homogeneous ensembles.
When the assumption that g(a) is a narrow symmetric distribution does not hold, (23) shows that the overall inverse dependence ofx i on a i described in the preceding paragraph is still valid. However, details on howx i compares with (13) for each a i depend on the specific value of the integral I a . Here, we do not pursue this analysis further but just point out that, for a large class of distributions g(a), the computation of the integral can be performed using calculus techniques in the complex plane, or exploiting the fact that it coincides with the double Laplace transform of g(a), namely, I a = L 2 [g].
The left panel of figure 3 shows, as symbols, the values ofx i computed along 10 4 time units in a numerical simulation of an ensemble of 1000 agents. The individual parameters a i are uniformly distributed over the interval (0.5, 1.5), while q = 2.5, λ = 1, and u = 1. Note that, for this choice, some of the agents in the ensemble have power-law resource distributions (a i < λ), while others have their resources distributed over a finite interval (a i > λ). In the same plot, the curve shows the first-order Sommerfeld approximation (26). Although the distribution of a i is relatively wide, the approximation turns out to be excellent.

Distribution in λ i
Mirroring the analysis of the previous section, we now consider the case where the ensemble is nonhomogeneous with respect to the multiplicative parameter λ i , while all the other parameters are identical between agents. Assuming that, in the continuous limit, λ i is distributed according to a function h(λ) over the interval (0, λ max ), with λ max < q (see section 2), we find The average of individual resources turns out to bē According to this result,x i increases as λ i grows, representing the fact that the resources of agents with stronger multiplication must be relatively larger on the average. When h(λ) is a narrow symmetric distribution centered ath, proceeding as in section 4.1 yields to the first significant order in σ λ , the standard deviation of λ. At the center of the distribution, λ i =λ,x i is now larger than the value expected for a homogeneous ensemble with the same parameters. A distribution in λ i implies that resources are being created at different speed for each agent. A consequence of this fact is that the mean value ofx i over the ensemble is no more the same as for the homogeneous case or for the case with a distribution in a i only, and now depends on h(λ). To the leading order in the Sommerfeld approximation, we have Thus, for this particular form of h(λ), total resources over the ensemble are larger than when the multiplicative factor equals the average valueλ for all agents. Moreover, in contrast with all the cases analyzed previously, total resources now depend on the value of the coupling strength a. Specifically, when (30) is valid, they decrease as a grows, meaning that redistribution entails an overall resource loss. Actually, it turns out that the decrease of total resources for increasing a occurs independently of the form of h(λ). In fact, from (30) we find ∂ ∂a with Σ 2 the variance of (q − λ + a) −1 over h(λ). In other words, whatever the distribution of the multiplicative factor λ, total resources shrink as coupling between agents becomes stronger. The equality in the last side of (31) only holds if Σ 2 = 0, i.e. when h(λ) is a delta function.
In figure 3, the right panel shows numerical results forx i in a simulation of an ensemble where λ i is uniformly distributed in (0.5, 1), with q = 2.5, a = 0.75, and u = 1. This ensemble, again, comprises both agents whose resources are distributed following power-law functions (λ i > a) and agents with resources bounded to finite intervals (λ i < a). As in the case of a non-homogeneous ensemble with distributed a i , the first-order Sommerfeld expansion, shown as a curve, provides a remarkably good approximation.

Spatially extended ensembles
As a final case of study, we relax the condition of globally coupled ensembles introduced in section 2 and consider a system of spatially distributed agents where resource redistribution occurs between neighbors. The preliminary numerical results presented below correspond to a linear array of N identical agents with periodic boundary conditions. Resource redistribution is homogeneous over a symmetric neighborhood of size 2K + 1 around each agent, including K neighbors to each side. The diffusive nature of this form of coupling between agents is apparent. The evolution equation for the resources of agent i, analogous to (1), readṡ where P i (t) is again a Poisson process representing reset events at rate q. Now,x i is the average of resources over the neighborhood N i of agent i, given bŷ We have performed numerical simulations of (32) for arrays of N = 1000 agents with K = 1, . . . , 10. Attention has been focused on the long-time stationary resource distribution of individual agents and on spatial correlations between individual resources at each time. As in section 3, we have assumed that individual resource distributions can be equivalently estimated from both temporal series and ensemble averages. Moreover, for each parameter set, spatial correlations in the stationary state were supposed to depend on the distance between agents only. Symbols in figure 4 represent the numerical estimation of the individual stationary resource distributions for the values of a already considered in section 3 (see figure 2) and four values of K, as indicated in the legend. Solid curves, in turn, correspond to the analytical result for globally coupled homogeneous ensembles. Firstly, we see that in all cases the distributions exhibit fatter tails than for the globally coupled situation, and that the larger the value of K the fatter the tail. This is due to the fact that, as the number of agents in each interacting group decreases, the statistical fluctuations in the reset times within the group become stronger. On the average, this allows for larger values of the individual resources.
For weak coupling (a = 0.5), the distributions seem to preserve their power-law decay for large x. For small x, differences between distributions for various K are comparable to those observed for other values of a, but they are less visible in the plot due to the logarithmic scale. For strong coupling (a > λ = 1), while individual resources of globally coupled agents are limited to a finite interval, distributions for the spatially extended ensemble protrude beyond the upper limit, with a tail that reaches much higher values of x. However, our numerical results have not made possible to unambiguously discern the functional form of the tail. For both weak and strong coupling, as expected, resource distributions for spatially extended ensembles approach their global-coupling counterparts as K grows. Meanwhile, for any given value of K the overall difference with the corresponding distribution for the globally coupled ensemble grows with the coupling strength a. This points to the fact that the more efficient the process of resource redistribution, the larger the effect of localized coupling between agents with respect to all-to-all interactions.
Instantaneous correlations between individual resources in spatially extended ensembles were quantified by the standard Pearson coefficient, wherex(t) and σ 2 (t) are the average and variance of resources over the whole ensemble. In our numerical simulations, C(t) was measured at regular intervals during a total of 10 8 time steps, and then averaged over time. The resulting correlation measure was recorded as a function of the distance d between pairs of agents, C(d). Figure 5 shows the numerical estimation of C(d) for the same parameter sets as in figure 4. Overall, correlations grows with a, as coupling between agents becomes stronger. In all cases, moreover, C(d) exhibits a slowly decaying, high plateau of width K, revealing the relatively similar values of resources within the neighborhood of each agent. Beyond d = K, the correlation drops abruptly and keeps decaying, now faster than on the plateau. For a = 0.5, due to the rapid decay of C(d), we were not able to determine a well-defined functional form for large d. Such low correlations can be ascribed to the fact that, for weak coupling, individual  figure 4) and resource redistribution is not able to establish significant correlations between interacting agents.
For strong coupling, on the other hand, our results suggest that C(d) decays exponentially for large d (see the straight lines in figure 5). The slope of the exponential tails decreases rapidly as K grows, roughly as K −1 . In contrast, it is relatively insensible to variations in a. In physical systems, exponentially decaying time and space correlations have been attributed to a variety of features, from finite-range interactions and spectral gaps [36,37] to mixing and chaotic dynamics [38]. Among stochastic processes, the classical Ornstein-Uhlenbeck process is the paradigm for such kind of correlations (in time [32]). Their origin in the present system of coupled multiplicative agents with reset events remains an open problem deserving future consideration.

Summary and conclusion
We have presented a statistical characterization of the stationary collective behavior of a set of coupled stochastic variables governed by the combination of a multiplicative process and random resetting. While the former induces exponential growth of the variables, the latter inhibits their long-time divergence. In the framework of a stylized socioeconomic model, each variable is interpreted as the amount of resources of an agent. Interaction between agents is introduced by means of a mean field, diffusion-like coupling scheme, in which a part of the resources of each agent is collected and the total is redistributed all over the ensemble. Our focus was put on the effects that this redistribution process has on the individual 'wealth' of agents, as well as on the total amount of resources in the system. Our analytical results, derived from a mean field-like approach for an infinitely large ensemble, were validated by means of numerical simulations.
If coupling is absent, the interplay between multiplication and resetting determines that, in the long run, individual resources are distributed following a pure power-law profile, with an exponent controlled by the multiplicative factor and the frequency of reset events. Our main result is to have disclosed a series of abrupt changes in the resource distribution as coupling is turned on and its strength increases. If the frequency at which resources are collected is less than the multiplication rate, the power-law distribution subsists but it is increasingly steeper, transforming into an exponential when the two rates become equal. Beyond this point, the tail of the distribution disappears and resources are confined to a finite interval. For stronger coupling, across a second critical point, the derivative of the distribution at the upper end of the interval changes from zero to infinity. Finally, after a third transition, the distribution develops a singularity at that end, leading to a sharp accumulation of resources around their mean value as coupling becomes definitely dominant.
When the ensemble is homogeneous-namely, when all agents have the same parameters-the mean value of individual resources is independent of the coupling strength, in spite of the disparate profiles that the resource distribution can adopt. This result is a direct consequence of the conservation of total resources during redistribution among identical agents. The situation changes, however, when non-homogeneous ensembles are considered. We have shown that, if the rate of resource collection is different between agents, the larger this rate the smaller the average of individual resources, as expected. Nevertheless, total resources over the ensemble are still independent of how collection rates are distributed. On the other hand, when multiplication rates are different from agent to agent, those agents with larger rates naturally have higher average resources, but total resources consistently diminish as the strength of coupling grows. This remarkable collective effect of resource redistribution seems to indicate that collecting resources from 'wealthy' agents, which are those with higher multiplication rates, may benefit less-favored agents upon redistribution, but produces an overall depletion of resources. In view of the minimalistic nature of our model, however, any implication on the functioning of real economic systems must obviously be drawn with caution.
Preliminary numerical results on the spatially extended version of the same model, where resource redistribution occurs in a limited neighborhood of each agent, have been presented for a linear array of agents with periodic boundary conditions and different neighbor sizes. We have shown that, in this version, individual resource distributions do not seem to be bounded to a finite interval when coupling is strong, as in the case of global coupling. This may be an indication that the cutoff found when agents are globally coupled is a specific feature of that case, in the limit of an infinitely large ensemble. The differences between distributions with local and global coupling grow with the coupling strength and, as expected, diminish for larger neighborhoods. We remark that the Chapman-Kolmogorov approach used in the case of global coupling cannot be straightly extended to spatially extended systems, due to the strong local fluctuations that preclude application of the hypotheses needed to apply the same analytical procedure. Further exploration of one-dimensional ensembles as well as analysis of other topologies (multi-dimensional arrays, networks) are worth considering in subsequent work. Both for globally coupled and spatially extended systems, statistical properties such as ergodicity, self-averaging, and predictability-which we have not addressed here-also deserve future attention, not only because of their intrinsic interest but also in view of their relevance to many applications, in particular, in economic theory [14,39,40]. Straightforward extensions of the model may incorporate agentto-agent cooperative or competing mechanisms, resource-dependent reset probabilities and collection rates (such as in progressive taxing systems), and coupling between reset events affecting different agents, among many other variants.
In summary, we have shown that a rather simple form of coupling between stochastic variables driven by multiplicative processes and random reset events gives rise to a rich variety of collective statistical behavior. This stimulates study of other kinds of interaction in similar models for complex socioeconomic systems.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.