# Examining inverse generative social science to study targets of interest

Thomas Chesney  Asif Jaffer  Robert Pasley

###### Abstract

We assess an emerging simulation research method—Inverse Generative Social
Science (IGSS) (Epstein et al., 2023)—that harnesses the power of evolution by
natural selection to model and explain complex targets.

Drawing on a review of recent papers that use IGSS, and by applying it in two
different studies of conflict, we here assess its potential both as a
modelling approach and as formal theory.

We find that IGSS has potential for research in studies of organistions. IGSS
offers two huge advantages over most other approaches to modelling. 1) IGSS
has the potential to fit complex non-linear models to a target and 2) the
models have the potential to be interpreted as social theory.

The paper presents IGSS to a new audience, illustrates how it can contribute,
and provides software that can be used as a basis of an IGSS study.

###### keywords:

agent-based modeling , genetic programming , simulation , conflict model

††journal: arXiv

[inst1]organization=Nottingahm University Business School,
addressline=Wollaton Road, city=Nottingham, postcode=NG8 1BB, country=UK

[inst2]organization=Karachi Institute of Business Administration,
addressline=University Road, city=Karachi City, postcode=Sindh 75270,
country=Pakistan
, state=/*authors listed alphabetically

##  1 Introduction

In response to (Melnyk et al., 2023) calling for a refresh of simulation
studies organisation studies, we assess an emerging simulation research
method—Inverse Generative Social Science (IGSS) (Epstein et al., 2023)—that
harnesses the power of evolution by natural selection to model and explain
complex targets. IGSS has the ability to explore a vast search space of
variables and their relationships seeking non-linear models to fit to data.
More importantly the best models found have the potential to act as formal
theory and as such bring with them the many advantages thereof, while also
answering calls for more formal theory in studies of organisations (Chesney,
2021; Vancouver et al., 2018). Drawing on a review of recent papers that use
IGSS, and by applying it in two different studies of conflict, we here assess
its potential both as a modelling approach and as formal theory.

IGSS marries the power of evolutionary computing (Holland, 1992)—specifically
here we use genetic programming (Koza et al., 1992)—and agent-based modelling
(Gilbert and Troitzsch, 2005). Genetic programming is used to evolve micro-
specifications—hereafter referred to as _rules_ —that when inserted into an
agent model’s code, dictate individual agents’ behaviour. When such a model is
run, if it generates macro behaviours that have been observed in a target of
interest then under the paradigm of generative or ‘bottom up’ social science
(Epstein and Axtell, 1996), that target is considered to have been explained.
We evaluate this view in Section 3. Before that, in Section 2 we outline IGSS.
Then in Section 4 we present two IGSS studies of conflict. A discussion of the
work is offered in Section 5.

##  2 Overview

A typical IGSS study will progress as follows. There exists a target for which
an explanation is sought. There also exist data on that target—or such data is
collected—called the reference dataset. The reference dataset is input to
software that evolves rules using genetic programming. The goal is to find
rules that when implemented and run in an agent-based model will simulate the
target such that it behaves as described in the reference dataset. This
progression is unpacked over the proceeding paragraphs.

The target is a phenomenon or situation of interest that is being studied. The
reference dataset stores observations on it capturing some interesting aspect
of how it behaves. It is this behaviour that is to be explained by evolved
micro rules. The reference dataset is quantitative and contains variables that
have been identified as being theoretically relevant, or at least are
suspected of such. This is the same as any quantitative research dataset.
Indeed the reference dataset may have been collected previously for another
study.

Agent-based modelling and genetic programming are well known in the
organisations and decision support literature (see for example: (Hajmohammad
and Shevchenko, 2020; Shevchenko et al., 2020; Braune et al., 2022; Shady et
al., 2022; Chesney et al., 2017; Zaffar et al., 2019)). Agent-based modelling
simulates a target allowing it to be experimented on, explored and observed
(Chesney, 2021). Genetic programming is the automatic evolution of computer
code (Banzhaf et al., 1998). To combine them requires genetic programming to
evolve rules that will then be used to dictate agent behaviour in an agent
model. (In addition the agent model will likely be required during evolution
to test evolved rules for fitness.)

Genetic programming starts with a population of random (and therefore almost
certainly meaningless) rules–the first generation of rules. For example:

  1. 1.

IF agentAttribute1 $>$ agentAttribute2 THEN behaviour1

  2. 2.

IF agentAttribute3 == globalVariable1 THEN behaviour2

  3. 3.

IF agentAttribute1 $/$ agentAttribute2 $<$ agentAttribute3 THEN behaviour1

  4. 4.

IF agentAttribute4 != globalVariable2 THEN behaviour3

In this example Behaviours 1, 2 and 3 are each one possible action that an
agent might take. Agent attributes are variables such as an agent’s age or a
summary of the situation they are currently in, perhaps the number of other
agents that are close by. Global variables hold the same value for all agents.
Examples could include a count of the number of agents, the current simulated
weather, or the legal regulations under which the agents are operating.

Next each of the rules is tested for fitness. Fitness is the answer to the
question: if agents follow this rule, how close will the output of the agent
model be to the reference dataset? A common fitness metric is the squared
difference between agent model data and reference data. It could be (and often
will be) that to calculate fitness the agent model must be run which can make
the process of completing an IGSS study take many weeks of computer
processing. Three genetic operations are then used repeatedly to create the
next generation of rules:

  1. 1.

One of the rules (chosen with probability proportional to its fitness, so most
likely a rule with high fitness) is allowed to reproduce and simply moves into
the next generation.

  2. 2.

One of the rules (again one most likely that has high fitness) mutates
randomly and the mutated rule moves into the next generation. A mutation means
that part of the rule is changed at random.

  3. 3.

Two of the rules (as before two that most likely have high fitness) breed and
the two child rules move into the next generation. This breeding is known as
crossover which means that part of the first rule and part of the second rule
are swapped with each other.

These genetic operations allow for the essential ingredients of evolution to
take place: reproduction, mutation and natural selection, which are known to
be powerful tools for finding solutions (Nowak and Rauh, 2006). The second
generation of rules is then tested for fitness. A third generation is created
as the second was and so on until a stopping condition is reached.

So far few IGSS studies have been conducted but all have used the approach
outlined above with small variation (the most notable of which is probably
using genetic algorithms rather than genetic programming but the overall goal
is the same). Greig et al. (2023) investigate flocking behaviour in birds.
Their reference dataset comes from a single run of an ‘ideal’ agent-based
model which was written based on observations of birds in flight (Reynolds,
1987). Miranda et al. (2023) studying field irrigation decisions by farmers
use reference data collected from an experiment using human participants. Vu
et al. (2023) use genetic programming to match survey data in a study of
alcohol consumption.

Creating a suitable test for fitness is often a challenging part of an AI
effort but here IGSS has a huge advantage. Given that the goal is matching—not
exceeding or extrapolating from or any other goal that might be assigned to an
AI—the reference dataset, the necessary existence of a reference dataset means
the fitness function can often simply be some measure of closeness, possibly
Euclidean distance or as was suggested earlier, mean squared error.

This overview glosses over a lot of complexity–conducting an IGSS study is not
a trivial matter. For example Vu et al. (2023) test fitness against multiple
objectives, not just one. The IGSS modeller must also develop the agent model
that will use the evolved rules and decide on parameters such as the agent
model’s scale, what type of agents to include and the numbers of each, other
initial parameter settings, and environment settings such as size and shape.
The modeller may also need to make decisions about the discrete time nature of
agent-modelling, specifically the point at which fitness is tested. This
decision is very much context dependent. The bird flocking model of Greig et
al. (2023) for example is neverending: the agent birds are constantly flying
with no end destination and therefore fitness can be tested at any point. This
is not the case in the study by (Vu et al., 2023) of alcohol use and they test
for fitness at a set point in time. Then after all of this, the final rule or
rules that have been produced will need to be pruned and interpreted.

##  3 Code as theory

We will use two phrases later that will be useful to define upfront. _Language
of the target_ refers to descriptions of behaviour using terms that have
meaning in the context of the target situation. If the target is salary
negotiation, the language of the target might include ‘wages’, ‘hardball’,
‘play hard to get’ and ‘hold your nerve’. _Language of the target theory_
refers to descriptions of behaviour using terms that have meaning in the
theory being used to explain the target. This time for salary negotiations
these might include ‘shared pot’, ‘utility’, ‘bounded choice’, ‘rational’ and
so forth.

The word theory in social science is very imprecise, something that has long
been criticised (Sutton and Staw, 1995). Table 1 shows some of the ways it is
used and the list is certainly incomplete. It is unlikely that this paper—or
any other—will convince readers to agree on a meaning for theory, but we can
make clear the definition that we use: theory is an explanation of data. Those
data come from observing a target. Theory abstracts out elements of interest,
as defined by the researcher, from a target and explains them. The theory will
not attempt to explain more than these elements. (This is one reason why
theory can be so poor when used for prediction.) Theory can appear in many
formats. One format that will be relevant to our discussion later is computer
code. (Again, perhaps not all readers will agree with this but a justification
will be given at that time.)

The purpose of a social scientist is to contribute to theory–to explain. In
most scientific endeavours, theory appears twice. When using an agent model as
a research method, the first is that theory is used to develop a model, to
dictate what micro behaviour is implemented. It is however rare that theory
will be available to guide all micro behaviour. To make a simulation work,
often additional behaviour—possibly even from additional theories—will have to
be melded together. Discussing this in relation to using simulations for
social science, Van Hemel et al. (2008, p.98) say: _The use of such terms as
‘theory’, ‘framework’, ‘model’, and ‘paradigm’ in psychology and the social
sciences is as informal as the models themselves. One person‘s conceptual
model is another person‘s theory or framework…In psychology and the social
sciences, theorizing about a problem typically begins with verbal conceptual
models, which then may be elaborated and adjusted over time as relevant
empirical data accumulate. Formal mathematical models, computational models,
statistical models, etc. rely on verbal conceptual models to specify variables
and relations among them, although a host of extra assumptions and plausible
estimates are typically needed to translate a verbal theory into a workable
implementation._

The exception to this is exploratory research where theory is not used, or not
relied on as much, in building a model. Instead, when used in agent modelling,
researchers are keen to explore different ideas they have for what might be
going on in a target. These ideas do not have to be grounded firmly in theory
and are implemented as micro behaviour to allow a researcher to observe
whatever macro behaviour emerges from them.

The second time theory appears in scientific endeavours is as an output. The
contribution a scientist makes is to theory–developing a new theory; adapting,
confirming, or validating an existing theory. Theory is both an input and an
output; theory goes into science, theory comes out of science.

Although rare in social science, it is possible that theory will be formal
(Adner et al., 2009). It is recognised that formalising theory brings many
advantages (Edwards, 2010; Davis et al., 2007; Vancouver et al., 2016, 2010;
Farrell and Lewandowsky, 2010; Busemeyer and Diederich, 2010; Lomi and Larsen,
2001; Grand et al., 2016). Formats for this include formal logic, mathematics,
and computational models (which because they are used to simulate data
ultimately means they are computer code). Examples of formal logic are
difficult to find, see (Kamps, 1999). Mathematics is of course common if we
include results from linear regression and structured equation modelling,
although (Sutton and Staw, 1995) for example would not consider this alone to
be theory. Mathematical modelling work such as that used by for example (Gupta
and Chutani, 2020) is common although this method does blur with computational
models and there is a question over whether writers who use this approach
actually mean for their mathematical models to be interpreted as theory
(though perhaps this is not an essential characteristic of theory). An entire
field is devoted to computational models of social behaviour, for examples see
the Journal of Artificial Societies and Social Simulation.

Vancouver et al. (2018) demonstrate the process of turning imprecise theories
into computational models and discuss why this is valuable. By implementing a
theory as a computational model such as an agent-based model it can be tested
and explored. Running a model generates data that can be compared with
observations. This will demonstrate—prove in fact—that a model explains—or
does not explain—those observations (it doesn’t prove that it is the _correct_
explanation, only that it is _an_ explanation). Implementing a theory as a
computational model makes an imprecise social science theory formal and
precise, and allows for rigorous testing (Edwards, 2010; Edwards and Berry,
2010; Adner et al., 2009).

IGSS sits in the generative social science paradigm as described by Epstein
and Axtell (1996). They say, “ _we consider a given macrostructure to be
‘explained’ by a given microspecification when the latter’s generative
sufficiency has been established_ ” (Epstein and Axtell, 1996, p.177). Under
IGSS, the final chosen rule that has been selected from the evolved rules
(which will have high fitness but might be the most fit–this is explained
later in this section) should explain the target. If theory is an explanation,
then this rules should at least have the potential to be considered theory–we
would argue that it _is_ theory.

Vancouver et al. (2018) would likely agree. Perhaps Sutton and Staw (1995)
would too. Although they neither approve of or discount computer code as being
theory, citing a lack of agreement about whether a model can constitute theory
(p.371), agent model code does seem to fit well with their description of what
strong theory is. Agent model code answers the question of why behaviours
emerge. The code makes explicit the connections among elements in the target.
The models present causal relationships and identify the timing of events.
Code presents underlying processes that allow us to understand “systematic
reasons for a particular occurrence or nonoccurrence” and present the
microprocesses involved. A counter argument to this last point might be that
agent-based models do not concern inter-relationships among macro level
variables (although we think this is incorrect–the output from agent models do
allow us to examined macro level variables). Sutton and Staw (1995, p.378)
finish their section with a quotation from Karl Weick: a good theory
“explains, predicts, and delights”. There is certainly something delightful
about observing a phenomena emerge ‘in front of your eyes’ on a computer
screen while running an agent model.

At least in the sense of explaining social processes of interest, agent-based
models can be considered as theory (Chesney, 2021) although as we say this we
remember that not all readers will agree on what theory is. An agent model’s
code can be pruned to test and observe how robust the output is, thus
potentially simplifying a theory. In addition, a model/theory can be explored
to generate hypotheses that have not been thought of previously, and
computational theories can be tested in ways that would be expensive,
unethical, or impossible otherwise (Chesney, 2021; Vancouver et al., 2018,
2010). As such agent-based models offer real possibilities for theory
advancement that deserve more attention in this regard, see for example
(Oliveira and Secchi, 2021).

We now consider how this idea of formal theory fits with IGSS. IGSS is best
seen as exploratory research. Exploratory research is used to study a
phenomenon that has not been studied previously. It is used to identify
research questions for further examination and hypotheses that will later be
tested. An exploratory researcher will likely get to prioritise these
questions and lay out the future research agenda. As for statistics,
exploratory cannot use null hypothesis significance testing. To do so would be
a form of data dredging. (To explain why briefly: if 20 random ideas are
tested using a significance threshold of 5%, it would be expected that 1 of
them would be observed because 5% equates to a 1 in 20 chance. Therefore if 20
ideas that might be real effects are explored in an agent model, but in fact
all of them are just random, one of them will probably wrongly be identified
as being real.) Instead, exploratory research relies on graphical methods,
descriptive statistics such as measures of central tendency, spread and
correlation, and qualitative approaches. Many of these, such as interviews and
focus groups, are inappropriate for agent models. Instead, when traditional
agent modelling is used for exploratory research, it involves exploring ranges
and combinations of parameter settings, and less frequently exploring a small
number of competing rules.

The jurisdiction of IGSS is this later space–exploring competing sets of
rules, many more than could possibly be explored manually. Indeed, IGSS can be
considered ‘turbo charged’ exploratory research. The consequence of this is
that IGSS should firmly be considered exploratory research. Indeed if our
interest was in implementing existing theory as an agent-based model then we
would not have need for IGSS. Even if theory is guiding the design of the
agents’ environment, or guiding the selection of which agent breeds to include
or the numbers of each to have, eventually the genetic programming algorithm
must take over and learn micro rules on its own, at which point it will not be
being guided by existing theory–it will be _exploring_ its search space. The
consequence is that theory as an input is less relevant to IGSS than to other
applications of agent models.

There are a infinite number of micro rules that will produce any macro
behaviour perfectly and it would be impossible to decide between them simply
by qualitatively or quantitatively observing the output of a model. When
discussing this, Chesney (2021, p.45) uses the analogy of rolling a die:
_Imagine you are holding a loaded die. If it is thrown on the ground, six
comes up most frequently, followed by five, then four, and so on down to one.
If however it is thrown against a certain slope, the frequency of getting a
six is reduced. In fact if you throw the die now, each side comes up about 16
times in 100, exactly as expected from a normal die. Is the die now uniform?
Or more simply, is it right? The system of the die’s use has now produced a
fair die even though we have not actually changed what most people would see
as the ‘correct’ part of the system that should have been changed._

This issue is dealt with explicitly in extant IGSS studies. Greig et al.
(2023) recognise the problem thus: “it is important to clarify that different
solutions in program space can potentially map to the same solution in logical
space” and state that the solution will take a “new set of non-trivial tools”.
Miranda et al. (2023) list their 10 best rules which are quite different from
each other and there is no obvious way to choose between them other than
closeness to the reference dataset, which as with the die metaphor does not
necessarily give the ‘correct’ solution. Interestingly they examine 10,000
rules using factor analysis to try to identify common elements (they refer to
them as ‘strategies’) found within them. This is a very exciting approach and
is discussed further later.

Vu et al. (2023) acknowledge the problem and deal with it in part by supplying
their genetic algorithm with a “library of theory building blocks” which are
“entities, attributes and mechanisms”. These come from an existing
computational model (Buckley et al., 2022) which is itself based on the Theory
of Planned Behaviour (Ajzen, 1991). As part of their discussion they state
that the fact that the genetic algorithm uses all of the building blocks
demonstrates that they “are all of theoretical importance for determining
individuals’ intentions” to behave in a certain way (in their case drinking
alcohol). They then (Paragraph 6.4) go on to list additional empirical
evidence as support of the building blocks’ importance. This of course gives
one of the ways in which the infinity of models is reduced, by giving the IGSS
algorithm access to theoretically valid constructs to work with. The pity in
this is that there may exist constructs that hadn’t yet been considered by
human researchers that the AI will not be given opportunity to find.

Epstein et al. (2023) and Epstein (2014) discuss the concept of cognitively
plausible agents. This means that fitness is not enough to judge evolved
rules. A rule has to be plausible. It must be possible to interpret it
meaningfully in either the language of the target or the language of the
target theory and in Section 5 we discuss this in relation to our models from
Section 4. Epstein et al. (2023) points out that a purely rationale agent for
example would not be cognitively plausible and suggest that agent emotion,
deliberation and social influence must contribute to interpreting and even
guiding evolved behaviour.

Genetic programming is favoured for IGSS rather than many other forms of AI as
it produces output that has strong potential to be meaningfully interpreted
(although it is not the only possibility–a decision tree is another
candidate). Contrast this with the output of an algorithm which trains a
neural network which would be an interconnected series of weights something
that could not possibly be interpreted meaningfully.

Meaning |  Example |  Reference
---|---|---
academic discipline |  gender theory |  Jule (2014)
research method or paradigm of study |  grounded theory |  Martin and Turner (1986)
topic of study and research method |  game theory |  Camerer (2011)
a belief |  marxism or marxist theory |  Singer (2018)
an umbrella term for a group of theories |  theory of evolution |  Nowak and Rauh (2006)
an explanation of data |  this is the definition favoured here |  Chesney (2021)
answers to normative questions |  theory of how society should be |  Ougaard (2013)
Table 1: Uses of the word ‘theory’ in social science.

###  3.1 Metatheory

We end this section with a comment on a possibility that may be of interest to
researchers. Rather than generating a theory, IGSS generates _candidate_
theories, competing alternative explanations. This suggests a role for IGSS in
the study of metatheory. Metatheory is theory on theory, or theorising about
the development and refinement of theories. A theory explains why certain
behaviours are observed in a target. There exists an underlying process that
creates those behaviours. It is this process that is being explained as it
reveals itself in the target. With a traditional agent-based model, agents are
told what this process is. With IGSS this behaviour is evolved in what was
earlier called ‘turbo charged’ exploratory research. This is valuable in and
of itself but it will also allow for the study of how rules evolve; such a
study would contribute to metatheory, a theory about the development of a
theory, an explanation of where rules come from. We can picture IGSS as
building a model of models to examine under what conditions a certain
observation that is explained by a certain theory emerges.

For example, the factor analysis carried out by Miranda et al. (2023) might be
considered to be an early instance of this. Factor analysis is the authors’
term and might be misleading. They use ‘factor’ to refer to a set of commands
that are applied to a number of parameters to produce a return value
(Gunaratne and Garibay, 2020). In social science terms examples of factors
would a construct such as personality, age, height, or a behaviour strategy
such as playing tit-for-tat. From the 10,000 rules that were evolved, they
assess how each factor and how joint factors (factor interaction) impact on
model fitness by measuring joint contribution. This involves training a random
forest regressor on the factor presence to fitness data and assessing this
using tree analysis software (Saabas, 2019). The most important factors then
are the theoretical elements that when included lead to a theory being
successful. We suggest that this can be considered the beginnings of a study
of metatheory.

##  4 Studies of conflict

Understanding and dealing with conflict is an important part of research into
organisations. To illustrate IGSS and guide our discussion, in this section we
1) evolve behaviour in a hawk-dove game and 2) re-create the well known civil
disobedience model.

###  4.1 Hawk dove model

For the hawk-dove agent model, in each time period agents move to a location
upon which is a certain amount of a desirable resource. Locations can
accommodate up to two agents at a time. Each agent must assess their situation
and decide how much of the resource to try to take. If there is conflict
between two agents over a resource then both get zero, otherwise they take
what they wanted. One run of the simulation lasts 100 time periods. The
decision of how much resource to try to take is the rule that is evolved using
IGSS with the format:

IF condition THEN take x ELSE take y

Agents initially choose a location at random and will return to it unless they
receive zero, in which case they choose another location at random. An agent’s
situation is recorded as a number of variables that are then used to construct
the rule. These are:

  1. 1.

amount of resource on my current location, $resource$

  2. 2.

number of agents on my current location, $agents$

  3. 3.

amount of resource I took in the previous round, $previousTook$

  4. 4.

amount of resource on my location in the previous round, $previousResource$

  5. 5.

number of agents on my location in the previous round, $previousAgents$

  6. 6.

number of agents in the simulation, $totalAgents$

  7. 7.

total amount of resource I have taken, $totalResource$

In this study we created reference datasets manually to capture various
scenarios of interest. We looked at: equality (where all agents end up with
the same total amount of resource, roughly the maximum possible for each) and
several forms of inequality (using the wealth distribution in the UK described
by The Equality Trust as a guide). The study was used to answer the question:
how would all members of a society have to act in order to achieve this
distribution of wealth? The code is written in R and is available at
github.com/ThomasCNotts/RIGSS.

####  4.1.1 Results

To achieve near perfect equality, IGSS found the following rule:

  * IF (previousResource - previousResource) * (previousResource - previousResource) $>=$ (previousTook - resource) - (totalResource - agents) THEN take = 1

Evolved rules almost always have to be pruned to remove redundancy. In this
case, the rule prunes down to ‘always take 1’ which given the amount of
resource and the number of agents in the simulation is a simple way to achieve
equality. (The rule actually reduces to ‘if zero is greater than a negative
number then take 1’.)

For inequality we started with a reference dataset that produces the
distribution of wealth shown in the right panel of Figure 1. To try to match
this IGSS found the rule:

  * IF (previousTook - totalAgents) != (agents AND previousTook) THEN take = 1 ELSE take = 9

which with pruning reduces to:

  * IF previousTook $>=1$ THEN take = 1 ELSE take = 9

and produces the distribution shown in the left panel of Figure 1.

Figure 1: The left panel shows wealth inequality with two extremes, split
between two halves of society. The right panel shows the inequality produced
by IGSS. Figure 2: A reference dataset that is closer to the real wealth
inequality in the UK.

###  4.2 Civil disobedience model

The second illustrative model was written in NetLogo which is perhaps less
well known than R. A specialist agent development environment, NetLogo is free
to download. The civil disobedience agent model (Epstein, 2002) models
violence flaring up in a densely populated environment. A version called
Rebellion.nlogo is included in the NetLogo models library. Citizen agents are
either happy or unhappy with their government, and are worried about arrest if
they protest with violence. They calculate their political grievance based on
an individual perceived hardship parameter and a government legitimacy global
parameter, and calculate their chance of arrest based on the number of other
people rebelling and the presence of police. If an agent’s grievance exceeds
the risk of arrest by a small threshold, the agent decides to rebel.
Occasionally, given certain initial parameters, there is a city-wide riot; but
what other mirco rules would lead to the same behaviour? An IGSS study is a
suitable way to explore this question.

There are three rules governing behaviour:

  1. 1.

Rule M: if there are any new sites free, citizens move to a new site if they
are not in jail; police move to a new site.

  2. 2.

Rule A: citizens who are not in jail decide if they should riot. The rule for
this is: riot if (grievance - risk-aversion * estimated-arrest-probability)
$>$ threshold.

  3. 3.

Rule C: if police sees a rioter or rioters, they arrest one of them randomly.

The civil disobedience model has an advantage for researchers who want to run
their first IGSS study in that the behaviour executed in each tick (the
NetLogo term for a discrete unit of time) is independent of events that have
happened in the previous tick–agents do not have to remember or factor in
their past. This means that we can generate a static reference dataset and use
this throughout learning. The alternative would be having to run the model
every time we wanted to test the fitness of a generated rule as was the case
with the hawk dove model which slows running considerably. To handle the A–M–C
order of rules, we evolve each rule separately, one at a time.

We ran the model for 40 time periods using three different combinations of
global variables, giving us data covering 120 time periods. For each time
period we collected data on every agent, recording the agent variables. Since
some internal agent states change more than once during a time period we
captured certain agent variables at three different steps in each period.
Essentially all data—including any hidden states—that determines behaviour was
recorded in a csv file which became the reference dataset. A total of 122,160
data points were generated with 38 variables making up each point. These data
were given as input to a second NetLogo program that implemented the genetic
programming algorithm.

As before genetic programming was used to generate a candidate rule of the
form IF condition THEN action. Due to bias in the reference dataset, in that
in most time periods many more agents choose not to riot than choose to riot,
we did not use a simple accuracy measure to test fitness. Instead we used a
balanced accuracy measure as a fitness measure. This takes into account both
sensitivity (true positive rate) and specificity (true negative rate) to
provide a more comprehensive assessment of the model’s ability to classify
instances from different classes correctly.

Rule M, the move rule, is relatively simple and the algorithm discovered a
rule with perfect fit quickly. It is not qualitatively identical to the actual
Rule M in the original Rebellion.nlogo but is quantitatively the same:

  * IF 0 $>=$ jail-term0 AND threshold $<$ free-neighborhood0 [set movement-calculated 1]

Rule A is by far the most complicated of the three and it took much longer to
achieve 98.07% fitness:

Rule A:

  * IF ((movement-tracker0 = 0 AND movement-tracker1 = 1 OR active-binary $>$ 0) OR government-legitimacy * government-legitimacy $<$ perceived-hardship) AND (count cops-on neighborhood) * estimated-arrest-probability = jail-term [set active? True]

In the language of the target:

  * IF (((agent did not move last time AND agent moved this time) OR agent is active since last time) OR government-legitimacy squared $<$ hardship) AND count police * arrest probability = jail term of the agent THEN riot

Like Rule M, a version of Rule C was discovered with perfect fitness:

Rule C:

  * if breed = cops [enforce]

We created a new model with the new evolved rules and compared the results
with the original model. The results are shown in Figure 3. Qualitatively the
two would be difficult to distinguish.

Figure 3: The left panel shows the output from the original Rebellion.nlogo;
the right panel shows the same output from Rebellion featuring our evolved
rules.

##  5 Discussion

In 2009, Frigg and Reiss (2009) published an examination of claims made about
the philosophy of computer simulation arguing that far from demanding a new
‘metaphysics, epistemology, semantics and methodology’, simulation raises few
if any new philosophical problems. As IGSS is exploratory research it shares a
philosophy of science with exploratory research. The fact that IGSS generates
candidate theories quantitatively at a larger scale than ‘manual’ exploratory
research does not in itself bring fresh philosophical challenges. However
there is much to get excited about. IGSS offers two huge advantages over most
other approaches to modelling. 1) IGSS has the potential to fit complex non-
linear models to a target and 2) the models have the potential to be
interpreted as social theory.

One of the most interesting prospects is that the candidate theories IGSS
generates can be _non-linear_. Not only that but IGSS allows for the
possibility of a candidate theory to be discovered made up of multiple non-
linear behaviours belonging to different types of agents thus introducing
heterogeneity between different types of agents. This can be compared with
model discovery in the physical sciences where AI has long been used to search
for linear and non-linear ordinary differential equations that yield an
observable relationship between inputs and outputs (see for example Schmidt
and Lipson (2009)). The downside of this is that the search space will be
enormous. While the review in this paper has highlighted several approaches to
reducing it, this remains a avenue for future research.

Our illustrative examples raise many points to note about IGSS. Starting with
specific results before moving on to make general comments about IGSS, for the
hawk-dove model, while our system did not find a good fit for either of our
inequality models, it did find a rule that leads to a smoother (and indeed
therefore more realistic) unequal distribution of wealth (right side panel,
Figure 1). The rule that gave allowed this to emerge exhibited extremes; in
the language of the target it is that agents take either very little or a lot:

  * IF previousTook $>=1$ AND agents $>=1$ THEN take = 1 ELSE take = 9

According to the preceding argument, this rule can explain wealth inequality.
This should be explored further in a fresh agent model to probe the
relationships between extremes and wealth distribution. (Strictly as this
would still be exploratory research, it would not have to be a fresh agent
model although a fresh model would be needed if hypothesis testing were going
to be used–see Chesney (2021, p.104). Indeed a brief exploration of this rule
in our agent model reveals that any values where the second is greater than
the first will lead to inequality, but the values above maximise it.)

Looking at the civil disobedience model, Rule A as presented in the previous
section may seem like it it is not explainable and this is a good case where
the rule must be interpreted in the language of the target or the target
theory. To recap the evolved version of the rule is:

  * IF ((movement-tracker0 = 0 AND movement-tracker1 = 1 OR active-binary $>$ 0) OR government-legitimacy * government-legitimacy $<$ perceived-hardship) AND (count cops-on neighborhood) * estimated-arrest-probability = jail-term [set active? True]

This can be interpreted in the language of the target as:

  * IF (((agent did not move last time AND agent moved this time) OR agent is active since last time) OR government-legitimacy squared $<$ hardship) AND count police * arrest probability = jail term of the agent THEN riot

These parameters can be further interpreted. For example, arrest probability =
jail term is effectively arrest probability = 0. This rule then starts to make
more sense. Human interpretation like this will be essential to make sense of
complicated rules and as we have said before, this sort of interpretation is
what a social scientists does so we should not fear it.

More generally, from our studies it is apparent that when running an IGSS
study, there are a number of decisions that must be made almost arbitrarily.
Each of these decisions could be a subject for further research. For example,
the choice of primitives to give to the genetic programming algorithm to work
with. Primitives are the mathematical operators from which the rules are built
and there is very little theoretical guidance on what these should be. Other
examples are the size of the initial population of rules and the probabilities
of each genetic operation (that is, the chance of reproduction, mutation or
crossover happening to a rule). In our hawk-dove model, the number of
generations was chosen arbitrarily. Each of these decisions could potentially
impact on results and yet there currently exists no systematic way to decide
on all of them other than ‘rules of thumb’. We ran our hawk dove model for 100
time periods, the civil disobedience for 40–again these are arbitrary
decisions that we suspect will have little bearing on results but this should
be studied further.

Lastly we give some suggestions for possible IGSS studies in the area of
decision support and organisations. IGSS could be used to model individual-
level decision-making in supply chains (e.g., purchasing behavior, supplier
relationship management) and simulate how these micro-behaviors aggregate to
produce macro-level outcomes like market fluctuations, supply chain
disruptions, or resilience. This could help predict and understand complex
phenomena that emerge from individual interactions within the supply chain.
IGSS could explore novel and potentially more efficient or ethical supply
chain configurations that wouldn’t be easily discovered through traditional
optimization approaches. Finally IGSS could model how individual agents adapt
their behavior when faced with changes like new regulations, economic crises,
or technological advancements providing insight into potential unintended
consequences of policies or predict how supply chains react to disruptions,
aiding in better preparedness and policy design.

We finish with a final observation. An agent model is not a necessary part of
an IGSS study. A reference dataset can be generated in any number of ways,
only one of which is an agent model. A genetic programming algorithm could be
run on data produced in an experiment, a survey, or an existing dataset to
produce rules. It may be natural though not essential to then implement those
rules in an agent model to demonstrate that they lead to the behaviour
observed in the target; this is of course the idea behind generative social
science. The main point is, that when IGSS software becomes easier to use,
there will be many existing datasets that have be used to create many accepted
theories that could be explored with IGSS to produce alternative theories that
can then be assessed by researchers.

## References

  * Adner et al. (2009) Adner, R., Polos, L., Ryall, M., Sorenson, O., 2009.  The case for formal theory.  Academy of Management Review 34, 201–208.
  * Ajzen (1991) Ajzen, I., 1991.  The theory of planned behaviour.  Organizational Behavior and Human Decision 50, 179–211.
  * Banzhaf et al. (1998) Banzhaf, W., Nordin, P., Keller, R.E., Francone, F.D., 1998.  Genetic programming: an introduction: on the automatic evolution of computer programs and its applications.  Morgan Kaufmann Publishers Inc.
  * Braune et al. (2022) Braune, R., Benda, F., Doerner, K.F., Hartl, R.F., 2022.  A genetic programming learning approach to generate dispatching rules for flexible shop scheduling problems.  International Journal of Production Economics 243, 108342.
  * Buckley et al. (2022) Buckley, C., Field, M., Vu, T.M., Brennan, A., Greenfield, T.K., Meier, P.S., Nielsen, A., Probst, C., Shuper, P.A., Purshouse, R.C., 2022.  An integrated dual process simulation model of alcohol use behaviours in individuals, with application to us population-level consumption, 1984–2012.  Addictive behaviors 124, 107094.
  * Busemeyer and Diederich (2010) Busemeyer, J.R., Diederich, A., 2010.  Cognitive modeling.  Sage.
  * Camerer (2011) Camerer, C.F., 2011.  Behavioral game theory: Experiments in strategic interaction.  Princeton university press.
  * Chesney (2021) Chesney, T., 2021.  Agent-Based Modelling of Worker Exploitation.  Springer.
  * Chesney et al. (2017) Chesney, T., Gold, S., Trautrims, A., 2017.  Agent based modelling as a decision support system for shadow accounting.  Decision Support Systems 95, 110–116.
  * Davis et al. (2007) Davis, J., Eisenhardt, K., Bingham, C., 2007.  Developing theory through simulation methods.  Academy of Management Journal 32, 480–499.
  * Edwards (2010) Edwards, J.R., 2010.  Reconsidering theoretical progress in organizational and management research.  Organizational Research Methods 13, 615–619.
  * Edwards and Berry (2010) Edwards, J.R., Berry, J.W., 2010.  The presence of something or the absence of nothing: Increasing theoretical precision in management research.  Organizational Research Methods 13, 668–689.
  * Epstein and Axtell (1996) Epstein, J., Axtell, R., 1996.  Growing artificial societies.  Brookings Institution Press.
  * Epstein (2002) Epstein, J.M., 2002.  Modeling civil violence: An agent-based computational approach.  Proceedings of the National Academy of Sciences 99, 7243–7250.
  * Epstein (2014) Epstein, J.M., 2014.  Agent_Zero: toward neurocognitive foundations for generative social science.  Princeton University Press.
  * Epstein et al. (2023) Epstein, J.M., Garibay, I., Hatna, E., Koehler, M., Rand, W., 2023.  Special section on” inverse generative social science”: Guest editors’ statement.  Journal of Artificial Societies & Social Simulation 26.
  * Farrell and Lewandowsky (2010) Farrell, S., Lewandowsky, S., 2010.  Computational models as aids to better reasoning in psychology.  Current Directions in Psychological Science 19, 329–335.
  * Frigg and Reiss (2009) Frigg, R., Reiss, J., 2009.  The philosophy of simulation: hot new issues or same old stew?  Synthese 169, 593–613.
  * Gilbert and Troitzsch (2005) Gilbert, N., Troitzsch, K., 2005.  Simulation for the social scientist second edition.  Open University Press.
  * Grand et al. (2016) Grand, J.A., Braun, M.T., Kuljanin, G., Kozlowski, S.W., Chao, G.T., 2016.  The dynamics of team cognition: A process-oriented theory of knowledge emergence in teams.  Journal of Applied Psychology 101, 1353.
  * Greig et al. (2023) Greig, R., Major, C., Pacholska, M., Bending, S., Arranz, J., 2023.  Learning interpretable logic for agent-based models from domain independent primitives.  Journal of Artificial Societies and Social Simulation 26.
  * Gunaratne and Garibay (2020) Gunaratne, C., Garibay, I., 2020.  Evolutionary model discovery of causal factors behind the socio-agricultural behavior of the ancestral pueblo.  Plos one 15, e0239922.
  * Gupta and Chutani (2020) Gupta, V., Chutani, A., 2020.  Supply chain financing with advance selling under disruption.  International Transactions in Operational Research 27, 2449–2468.
  * Hajmohammad and Shevchenko (2020) Hajmohammad, S., Shevchenko, A., 2020.  Mitigating sustainability risk in supplier populations: an agent-based simulation study.  International Journal of Operations & Production Management 40, 897–920.
  * Holland (1992) Holland, J.H., 1992.  Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence.  MIT press.
  * Jule (2014) Jule, A., 2014.  Gender theory, in: Michalos, A. (Ed.), Encyclopedia of Quality of Life and Well-Being Research, Springer.
  * Kamps (1999) Kamps, J., 1999.  On criteria for formal theory building: Applying logic and automated reasoning tools to the social sciences, in: AAAI/IAAI, pp. 285–290.
  * Koza et al. (1992) Koza, J.R., et al., 1992.  Genetic programming.  MIT Press Cambridge.
  * Lomi and Larsen (2001) Lomi, A., Larsen, E.R., 2001.  Dynamics of organizations: Computational modeling and organizational theories.  The MIT Press.
  * Martin and Turner (1986) Martin, P.Y., Turner, B.A., 1986.  Grounded theory and organizational research.  The journal of applied behavioral science 22, 141–157.
  * Melnyk et al. (2023) Melnyk, S.A., Thürer, M., Blome, C., Schoenherr, T., Gold, S., 2023.  (re)-discovering simulation as a critical element of om/scm research: call for research.  International Journal of Operations & Production Management .
  * Miranda et al. (2023) Miranda, L., Baggio, J., Ozmen Garibay, O., 2023.  Evolutionary model discovery of human behavioral factors driving decision-making in irrigation experiments.  JASSS: Journal of Artificial Societies and Social Simulation 26.
  * Nowak and Rauh (2006) Nowak, K.L., Rauh, C., 2006.  The influence of the avatar on online perceptions of anthropomorphism, androgyny, credibility, homophily, and attraction.  Journal of Computer-Mediated Communication 11, 153–178.
  * Oliveira and Secchi (2021) Oliveira, N., Secchi, D., 2021.  Theory building, case dependence, and researchers’ bounded rationality: An illustration from studies of innovation diffusion.  Sociological Methods & Research , 0049124120986201.
  * Ougaard (2013) Ougaard, M., 2013.  What is Theory. CBS Press, Copenhagen. chapter What is Theory in Political Science.
  * Reynolds (1987) Reynolds, C.W., 1987.  Flocks, herds and schools: A distributed behavioral model, in: Proceedings of the 14th annual conference on Computer graphics and interactive techniques, pp. 25–34.
  * Saabas (2019) Saabas, A., 2019.  Treeinterpreter software.  available at GitHub URL: https://github.com/andosa/treeinterpreter.
  * Schmidt and Lipson (2009) Schmidt, M., Lipson, H., 2009.  Distilling free-form natural laws from experimental data.  science 324, 81–85.
  * Shady et al. (2022) Shady, S., Kaihara, T., Fujii, N., Kokuryo, D., 2022.  A novel feature selection for evolving compact dispatching rules using genetic programming for dynamic job shop scheduling.  International Journal of Production Research 60, 4025–4048.
  * Shevchenko et al. (2020) Shevchenko, A., Pagell, M., Lévesque, M., Johnston, D., 2020.  Preventing supplier non-conformance: Extending the agency theory perspective.  International Journal of Operations & Production Management 40, 315–340.
  * Singer (2018) Singer, P., 2018.  Marx: A very short introduction.  Oxford University Press.
  * Sutton and Staw (1995) Sutton, R.I., Staw, B.M., 1995.  What theory is not.  Administrative science quarterly , 371–384.
  * Van Hemel et al. (2008) Van Hemel, S.B., MacMillan, J., Zacharias, G.L., et al., 2008.  Behavioral Modeling and Simulation: From Individuals to Societies.  National Academies Press.
  * Vancouver et al. (2016) Vancouver, J.B., Li, X., Weinhardt, J.M., Steel, P., Purl, J.D., 2016.  Using a computational model to understand possible sources of skews in distributions of job performance.  Personnel Psychology 69, 931–974.
  * Vancouver et al. (2010) Vancouver, J.B., Tamanini, K.B., Yoder, R.J., 2010.  Using dynamic computational models to reconnect theory and research: Socialization by the proactive newcomer as example.  Journal of management 36, 764–793.
  * Vancouver et al. (2018) Vancouver, J.B., Wang, M., Li, X., 2018.  Translating informal theories into formal theories: The case of the dynamic computational model of the integrated model of work motivation.  Organizational Research Methods 23, 238–274.
  * Vu et al. (2023) Vu, T.M., Buckley, C., Duro, J.A., Brennan, A., Epstein, J.M., Purshouse, R.C., 2023.  Can social norms explain long-term trends in alcohol use? insights from inverse generative social science.  Journal of artificial societies and social simulation: JASSS 26.
  * Zaffar et al. (2019) Zaffar, M.A., Kumar, R.L., Zhao, K., 2019.  Using agent-based modelling to investigate diffusion of mobile-based branchless banking services in a developing country.  Decision Support Systems 117, 62–74.


[a]G. Morlino

# Contribution to Galactic cosmic rays from young stellar clusters

S. Menchiari  E. Amato  N. Bucciantini

###### Abstract

The origin of Galactic cosmic rays (CR) is still a matter of debate. Diffusive
shock acceleration (DSA) applied to supernova remnant (SNR) shocks provides
the most reliable explanation. However, within the current understanding of
DSA several issues remain unsolved, like the CR maximum energy, the chemical
composition and the transition region between Galactic and extra-Galactic CRs.
These issues motivate the search for other possible Galactic sources.
Recently, several young stellar clusters (YSC) have been detected in gamma
rays, suggesting that such objects could be powerful sources of Galactic CRs.
The energy input could come from winds of massive stars hosted in the clusters
which is a function of the cluster total mass and initial mass function (IMF)
of stars. In this work we evaluate the total CR flux produced by a synthetic
population of YSCs assuming that the CR acceleration occurs at the termination
shock of the collective wind resulting from the sum of cluster’s stellar
winds. We show that the spectrum produced by YSC can significantly contribute
to energies $\gtrsim 100$ TeV if the diffusion inside the wind-blown bubble is
Bohm-like and the spectral slope is harder than the one produced by SNRs.

##  1 Introduction

The origin of Galactic Cosmic Rays (GCRs) in the knee region is still debated
in the community. Recently, young and massive stellar clusters (YMSCs) have
been suggested as alternative candidate sources to supernova remnants (SNRs).
The energy input could come from the cluster stellar winds which provide,
during a lifetime of several million years, an energy comparable with respect
to SNR outputs [1]. In addition, YMSCs could offer favorable conditions for
particle confinement. The stellar winds, in fact, generate a wind-blown bubble
around the cluster with a typical size of tens of pc, where the magnetic
turbulence may be enhanced with respect to the Galactic interstellar medium
(ISM), increasing the particle diffusion time and, consequently, the
possibility to achieve very high energies. Besides energetic considerations,
acceleration of particles from the wind of massive stars appears to be a
necessary component to explain the 22Ne/20Ne anomaly in CR composition [2]:
the enhancement observed with respect to Solar abundances requires
acceleration from a carbon-enriched medium rather than from the standard ISM
at a few percent level. The relative contribution of these different source
populations (SNRs and YMSC) to the observed CRs is yet to be clarified, in
particular in the region across the _knee_ around few PeVs.

In recent years, diffuse gamma-ray emission has been detected in coincidence
with many young massive stellar clusters (YMSC) by several gamma-ray
facilities, like Fermi-LAT, H.E.S.S. [3] and LHAASO. These findings strongly
supporting the idea that some acceleration mechanism is taking place there.
Detailed morphological and spectral analysis suggest a continuous injection of
possible hadronic origin [4, 5]. The Cygnus cocoon has even been observed at
the highest energies ever probed by gamma-ray astronomy: a 1.4 PeV photon was
detected from LHAASO [7], strengthen the idea that stellar clusters may act as
PeVatrons.

The exact location where particle acceleration takes place in YMSCs is still
unclear. For compact and young systems, the so-called wind termination shock
(WTS) [8] is expected to be strong enough to enable particle acceleration at
such high energies [9]. Alternatively, stochastic acceleration might be driven
by the highly turbulent environment of the cluster, particularly in its core
[10], further amplified once SN explosions start to occur [11]. In this work
we want to estimate the contribution of YMSCs to the Galactic CR flux before
supernova start to explode. Hence, we will consider only acceleration at the
cluster WTS. For each SC, we build a stellar population, which properties in
terms of wind speed and mass loss rate are used to determine the size of the
wind-blown bubbles, as described in §2. Then, we will integrate the WTS
contribution over the entire population of Galactic SCs. However, the SC
population is reasonably well know only within $\sim 2$ kpc from the Sun and,
to overcame this lack of information, we will build a synthetic population of
SC based on properties of local clusters, as explained in §3. Finally, in §4
we summarize the acceleration model of [8] and we discuss the results in §5.
The reader is also referred to a companion work presented in the same
conference [12] where the same approach is used to estimate the SC
contribution to the Galactic diffuse gamma-ray emission.

##  2 Properties of stellar winds and wind-blown bubbles

To describe the bubble structure around a SC, for each star we need two
quantities, the wind velocity and the mass loss rate and the star age. In the
following we will only deal with main sequence stars, neglecting the
contribution from the stellar final stages, like Wolf-Rayet, which, however,
may contribute up to $\sim 30\%$ of the wind power in a SC [1]. In the stellar
wind theory, the wind velocity is generally written as [15]

$v_{\rm w,\star}=C(T_{\rm eff})v_{\rm esc}=C(T_{\rm
eff})\left[2gR_{\star}\,(1-\Gamma)\right]^{1/2}$ (1)

where $v_{\rm esc}$ is the escape speed form the star, $g$ is the surface
gravity and $R_{\star}$ the stellar radius. The factor
$\Gamma=L_{\star}/L_{\rm Edd}$ takes into account the reducing effect of
Thomson scattering on the gravitational potential. The wind velocity is in
general larger than $v_{\rm esc}$ due to the radiation pressure from the star.
Such an effect is accounted for by the function $C$ which depends on the
effective surface temperature of the star which is, in turn, estimated using
the Boltzmann law: $T_{\rm eff}=\left[\frac{L_{\star}}{4\pi
R_{\star}^{2}\sigma_{b}}\right]^{1/4}$, $\sigma_{b}$ being the Boltzmann
constant. $C$ ranges from 1 for $T_{\rm eff}<10^{4}$ K up to 2.65 for $T_{\rm
eff}>2\times 10^{4}$ K. In this range we assume a linear increase with $T_{\rm
eff}$.

The stellar mass loss rate is a rather difficult quantity to constraint from
observations. Stellar evolution codes provide results which, however, depends
on several quantities like the metallicity and the the stellar rotation. For
the sake of simplicity, here we use the approximate model by [16] where the
mass loss rate depends only on the stellar luminosity, mass and radius, and
reads [for a comprehensive discussion see also 4]:

$\dot{M}_{\star}=9.55\times
10^{-15}(L_{\star}/L_{\odot})^{1.24}(M_{\star}/M_{\odot})^{0.16}(R_{\star}/R_{\odot})^{0.81}\,{\rm
M_{\odot}\,yr^{-1}}\,.$ (2)

The stellar luminosity is only a function of its mass and it is taken from [4]
and consists of a smoothed broken power law mixing two different empirical
mass-luminosity relations. The relation between stellar radius and mass is
given to first approximation by
$R_{\star}/R_{\odot}=0.85(M_{\star}/M_{\odot})^{0.67}$ [17]. We notice that
the uncertainty in the mass-ratio relation translates into an uncertainty of
only $\sim 15\%$ in the final cluster’s luminosity. For the purpose of this
work, the properties of stellar winds can be considered almost stationary
during the main sequence lifetime, which last

$\log{\left(T_{\rm age}/\rm
yr\right)}=6.43+0.825\,\left[\log{(M_{\star}/120\rm M_{\odot})}\right]^{2}\,.$
(3)

After such a time, we assume that the stellar wind do not contribute anymore
to the SC wind. We stress again that the subsequent explosion of SNe from
massive stars is neglected.

For the initial mass function (IMF) of star inside a cluster, we adopt the
distribution from [18] which is a broken power-law in several mass range,
$\xi_{\star}=A_{i}(M_{\rm cl})M_{\star}^{k_{i}}$, which reduces to a Salpeter
IMF for $M>M_{\odot}$ with $k_{i}=2.35$. The IMF for each cluster is
normalized to give the SC total mass, i.e.
$\int_{M_{\star,\min}}^{M_{\star,\max}}M\xi_{\star}(M_{\star},M_{\rm
cl})dM=M_{\rm c}$. The minimum and maximum stellar mass are assumed to be
$M_{\star,\min}=0.08M_{\odot}$ which is related to the minimum theoretical
mass to support significant nuclear burning and $M_{\star,\max}=150$ M⊙ that
is the maximum observed stellar mass.

Once the stellar distribution is fixed, the properties of the SC collective
wind can be estimated from the mass and momentum flux conservation,
integrating over all stellar winds. Hence, the final SC mass loss rate and
wind speed are:

$\displaystyle\dot{M}_{\rm c}(M_{\rm
c})=\int_{M_{\star,\min}}^{M_{\star,\max}}\dot{M}_{\star}\,\xi_{\star}(M_{\star},M_{\rm
c})\,dM_{\star}$ (4) $\displaystyle v_{\rm w,c}=\frac{1}{\dot{M}_{\rm
c}}\int_{M_{\star,\min}}^{M_{\star,\max}}\dot{M}_{\star}v_{\rm
w,\star}\,\xi_{\star}(M_{\star},M_{\rm c})\,dM_{\star}$ (5)

The bubble structure is described using the classical solution for adiabatic
expansion [19]. Defining the cluster age $t$ and the cluster luminosity
$L_{\rm w,c}=\frac{1}{2}\dot{M}_{\rm c}v_{\rm w,c}^{2}$, as well as the ISM
mass density $\rho_{0}$, the position of the TS is located at

$R_{\rm
s}(t)=48.6\,\left(\frac{\rho_{0}}{\mathrm{cm}^{-3}}\right)^{-0.3}\left(\frac{\dot{M}_{\rm
c}}{10^{-4}\mathrm{M}_{\odot}\mathrm{yr}^{-1}}\right)^{0.3}\left(\frac{v_{\rm
w,c}}{1000\,\mathrm{km\,s}^{-1}}\right)^{0.1}\left(\frac{t}{10\,\mathrm{Myr}}\right)^{0.4}\,\mathrm{pc}$
(6)

while the bubble radius is

$R_{\rm
b}(t)=174\,\left(\frac{\rho_{0}}{\mathrm{cm}^{-3}}\right)^{-0.2}\left(\frac{L_{\rm
w,c}}{10^{37}\,\mathrm{erg\,s}^{-1}}\right)^{0.2}\left(\frac{t}{10\,\mathrm{Myr}}\right)^{0.6}\mathrm{pc}\,.$
(7)

Because stellar clusters born from giant molecular clouds, the local ISM
density where their winds expand is usually denser than the average Galactic
ISM. Here we assume a reference value of $\rho_{0}=10$ protons/cm3 identical
for all SCs.

##  3 Stellar cluster distribution

The SC distribution in mass, time and position in the Galaxy is defined as

$\xi_{\rm c}(M_{\rm c},T,\vec{r})=\frac{dN}{dTdM_{\rm c}d\Sigma}\,,$ (8)

such that $\xi_{\rm c}dTdM_{\rm c}$ is the number of SC with initial masses
$[M_{\rm c},M_{\rm c}+dM_{\rm c}]$ formed per unit surface $d\Sigma$ of the
Galactic disk, at the position $\vec{r}$, in the time interval $[T,T+dT]$.
Following [20] we assume that the distribution is factorized in time and mass.
Moreover, we also assume that the distribution can be factorized in space such
that the cluster initial mass function (CIMF) depends only by the distance $r$
from the Galactic Centre through a normalization factor. Hence, we can write

$\xi_{\rm c}(M_{\rm c},T,\vec{r})=\psi(T)f_{\rm c}(M_{\rm c})\rho_{\rm c}(r)$
(9)

where $\psi$ is the SC formation rate (CFR), $f_{\rm c}(M_{\rm c})$ is the
CIMF and $\rho_{\rm c}(r)$ is the cluster radial distribution, normalized to
be unity at the Sun location, $r_{\odot}=8.5$ kpc. To get the present
distribution of SC we should integrate in time Eq.(9), however, here we are
interested in describing only the young population of SC, with an age
$\lesssim 10$ Myr, because for larger ages the wind power drops to negligible
values and the production of CRs becomes irrelevant. [20] showed that the
present SC distribution in the solar neighborhood is compatible with a
formation rate roughly constant during the last $\sim 50$ Myr. Hence the CFR
can be assumed constant. Its value can be derived from the work by [21] who
obtained a surface star formation rate in the solar neighbourhood of $350\,\rm
M_{\odot}\,Myr^{-1}\,kpc^{2}$ for clusters’ mass between $100\,\rm M_{\odot}$
and $3\times 10^{4}\,\rm M_{\odot}$. This corresponds to an average CFR of
$\bar{\psi}=0.63\,\rm kpc^{-2}\,Myr^{-1}$. For the CIMF, we follow [20] which
derived for the SC population in the solar neighborhood the following broken
power-law:

$f_{\rm c}(M_{\rm c})=\,\begin{cases}k_{1}\,M_{\rm c}^{-1.63}\hskip
14.22636pt{\rm for}\hskip 14.22636ptM_{\rm c,\min}<M_{\rm c}<M_{\rm c}^{*}\\\
k_{2}\,M_{\rm c}^{-1.24}\hskip 14.22636pt{\rm for}\hskip 14.22636ptM_{\rm
c}^{*}<M_{\rm c}<M_{\rm c,\max}\end{cases}$ (10)

where $M_{\rm c,\min}=2.5\,\rm M_{\odot}$, $M_{\rm c,\min}=6.3\times
10^{4}\,\rm M_{\odot}$, and $M_{\rm c}^{*}=100\,\rm M_{\odot}$. The constants
$k_{1}$ and $k_{2}$ are obtained from the continuity at $M_{\rm c}^{*}$ and
from the normalization condition $\int_{M_{\rm c,\min}}^{M_{\rm c,\max}}f_{\rm
c}(M_{\rm c})dM_{\rm c}=1$. Due to strong stellar light extinction in the
Galactic plane, the spatial distribution of SC is known with sufficient
accuracy only in the solar neighborhood (for a distance $\lesssim 2$ kpc from
the Sun [20]). As a consequence $\rho_{\rm c}$ should be derived from some
other proxy. Here we use the distribution of pulsars as derived by [22], which
reads

$\rho_{c}(r)=\left(\frac{r+r_{\odot}}{2r_{\odot}}\right)^{1.64}\exp{\left[-4.01\,\frac{r-r_{\odot}}{r_{\odot}+0.55{\rm
kpc}}\right]}$ (11)

where $r_{\odot}=8.5$ kpc is the Sun position. On top of the radial
distribution, we also account for the distribution inside the Galactic spiral
arms, using the same procedure that [23] adopted to evaluate the distribution
of SNRs. The spiral structure is realized by choosing a galactocentric
distance $r$ from Eq.(11) and than by choosing randomly an arm. The polar
angle is then determined so that the cluster lies in the centroid of the arm.
The actual position of the SC is computed by applying a correction to the
galactocentric distance drawn from a normal distribution centered at zero with
a standard deviation $0.07\,r$. Figure 1 shows the result of a single
realization of SC population with an age younger than 3 Myr, in terms of
radial distribution from the GC and position in the Galactic disk. The total
number of clusters results to be $\simeq 300$.

(a) (b)

Figure 1: Distribution of a single realization of stellar cluster population
with age $<3\,\rm Myr$, as a function of Galactocentric radius (a) and in the
Galactic plane (b). The solid line in (a) show the pulsar distribution from
Eq.(11) for comparison.

##  4 Cosmic ray acceleration

Following [8], we assume that the acceleration of particles only occurs at the
WTS developed by the SC wind. The (relativistic) energy distribution function
of accelerated particles located at the shock position can be written as

$f_{s}(E)=\xi_{\rm cr}\,\frac{L_{\rm
w,c}}{m_{p}c^{2}}\frac{1}{\Lambda_{p}}\left(\frac{E}{m_{p}c^{2}}\right)^{-s}\,e^{-\Gamma(E)}\,,$
(12)

where the normalization constant $\Lambda_{p}$ is defined such that the CR
luminosity of the system is equal to $\xi_{\rm cr}$ times the wind luminosity,
i.e. $L_{\rm cr}\equiv\int Ef_{s}(E)u_{2}\,dE=\xi_{\rm cr}L_{c,w}$, $u_{2}$
being the downstream wind speed. The solution in Eq. (12) has the typical
power-law term $\propto E^{-s}$, found for the case of plane parallel shocks,
plus an additional exponential function which accounts for the effects due to
the spherical geometry and to the escaping of particles from the bubble
boundary, which determine the maximum energy of the system. The exponential
function has a complicated expression which, however, can be approximated by
the following formula:

$e^{-\Gamma(E)}\simeq\left[1+A(E/E_{\max})^{b}\right]\,e^{-k(E/E_{\max})^{c}}\,.$
(13)

where $A$, $b$, $k$ an $c$ are fitting parameters, while $E_{\max}$ is the
nominal maximum energy defined by the condition that the upstream diffusion
length is equal to the shock radius, i.e. $D_{1}(E_{\max})/u_{1}=R_{s}$,
$D_{1}$ being the upstream diffusion coefficient. The diffusion properties
inside the bubble represents the most uncertain parameter of the system. Again
following [8] we parameterise the diffusion coefficient as
$D=v/3\,r_{L}^{\delta}\,\left(r_{L}/L_{c}\right)^{1-\delta}$, where $r_{L}$ is
the Larmor radius, while $L_{c}\simeq 1$ pc, is the coherence length-scale of
the magnetic field, assumed to be of the order of the SC core radius. The
exponent $\delta$ is equal to 1/3, 1/2 and 1 for Kolmogorov, Kraichnan and
Bohm-like diffusion, respectively. Finally, the magnetic field, $\delta B$,
used to evaluate the Larmor radius, is estimated assuming that a fraction
$\eta_{B}$ of the wind luminosity at the shock, is converted into magnetic
pressure, namely $(\delta B^{2}/4\pi)\,4\pi
R_{s}^{2}v_{w}=\eta_{B}\,\dot{M}\,v_{w}^{2}/2$. $\eta_{B}$ is expected to be
of the order of few percent.

On a very general ground, the spectral slope of accelerated particles is
determined by the effective compression ratio at the shock, $\sigma$, which
includes the velocity of the scattering turbulence. Using hybrid simulations,
it has been recently shown [24] that the downstream turbulence is, in general,
more effective than the upstream one in determining the slope, hence here we
will include only such an effect. In a parametric form, we can write the mean
velocity of the waves downstream as
$\bar{v}_{A,2}=\chi_{A}\,\sqrt{11/8}\,\eta_{B}^{1/2}\,v_{w}$, where $v_{A}$ is
the Alfvén speed and $\chi_{A}=0$ for waves that are symmetrically moving in
all directions. The comporession ratio is than written as

$\sigma=\frac{u_{1}}{u_{2}+v_{A,2}}=\frac{\sigma}{u_{2}+\sqrt{11/8\,\xi_{B}}\,\chi_{A}\,u_{1}}\,.$
(14)

The value of the parameter $\chi_{A}$ obtained from numerical simulations is
of the order of few tens of percent. However, in such simulations the magnetic
field amplification is only due to the CR streaming, while in the case of
stellar winds the magnetic field is more probably determined by MHD
instabilities. Hence it is not clear whether the results by [24] can be
straightforwardly applied to our case. As a consequence, we will take
$\chi_{A}$ as a free parameters.

Now we do have all the ingredients to evaluate the CR flux produced by SCs.
However, one last caveat need to be addressed. The instantaneous CR luminosity
is formally obtained from the particle flux escaping from the wind bubble,
which reads $\phi_{\rm esc}=4\pi R_{b}^{2}\,\left[D\partial f/\partial
R\right]_{R=R_{b}}$, and the total amount of CR injected by a single SC should
be given by the integral during its lifetime, $\int\phi_{\rm esc}dT$. However,
one can easy realize that such a contribution is always negligible with
respect to the amount of particles stored inside the bubble, which is
approximately given by $4\pi/3f_{s}R_{b}^{3}$. This apparent inconsistency is
given by the fact that the solution provided by [8] is stationary and does not
account for the time evolution of the wind-bubble after the adiabatic phase,
when the bubble will fade out, releasing all CR stored during the acceleration
phase. This issue can be solved replacing the escaping flux with the flux
injected inside the bubble, i.e. $\phi_{s}=4\pi R_{s}^{2}\,u_{2}f_{s}$. [8]
have shown that $\phi_{\rm esc}$ and $\phi_{s}$ differ only slightly at
energies $\gtrsim E_{\max}$.

##  5 Results and discussion

Using the approximation provided in Eq.(9) we can write the CR flux at the
present time $T$ injected by the entire population of SC as follows:

$Q_{\rm
SC}(E,T)\simeq\bar{\psi}\,\int_{0}^{T}\int_{M_{c,\min}}^{M_{c,\max}}f_{c}(M_{c},t)\,\phi_{s}(M_{c},t)\,dt\,dM\times\int_{0}^{R_{\rm
disk}}\rho_{c}(r)\,dr\,.$ (15)

In the present work, rather then performing the analytical integral, we
proceed generating a synthetic population of SC and then we sum up over the
contribution of all individual clusters. In Figure 2 we show the CR spectrum
injected by the synthetic population shown in Figure 1. Two different cases
are considered: $\chi_{A}=0$ and $\chi_{A}=0.1$. The CR acceleration
efficiency and magnetic amplification efficiency are fixed to $\xi_{\rm
cr}=0.05$ and $\eta_{A}=0.05$. The CR spectrum is shown for three different
assumption of the diffusion coefficient (Bohm, Kraichnan and Kolmogorov) and
is compared with the spectrum of CR injected by SNRs. Notice that all Figures
only show the proton CR component. Heavier elements are not discussed here.
One can see that for Bohm diffusion, the $E_{\max}$ reaches PeV energies, and
the SC contribution dominates the CR spectrum above $\sim 100$ TeV. The case
with $\chi_{A}=0$ show a harder spectrum ($s=2.03$) than the one having
$\chi_{A}=0.1$ ($s=2.18$).

Notice that the contribution from SNRs is estimated using some
simplifications. We assume that all SN explode into an uniform medium with
density 0.1 cm-3. The kinetic explosion energy is fixed to $10^{51}$ erg and
the ejecta mass to $5\,\rm M_{\odot}$. The particle acceleration and escape is
than calculated using the model by [26] which include the magnetic field
amplification due to streaming instability (even though in a simplified
approach). As one can see from Figure 2, the effective maximum energy produced
by SNRs is only $\simeq 50$ TeV (similar results are obtained using more
refined models like [27] for the same SNR parameters).

There are two important aspects to discuss. The relative normalization between
the SCs and the SNRe contributions and the different slopes. The normalization
of the CR flux produced by SNRs is obtained taking the Salpeter IMF for the
entire Galaxy stellar population and normalizing it to a SFR of $2\,\rm
M_{\odot}\,yr^{-1}$ [28]. Than we assume that all stars with $M>8\,\rm
M_{\odot}$ explode as SNe. This approach gives a rate of $\simeq 1$ SN every
century, a factor roughly 2 smaller (but still compatible within the errors)
than the one estimated from the combined evidence from external galaxies and
from the observation of historical SNe in our Galaxy [29]. Then, we assume
that the SNR acceleration efficiency is the same as the WTS, namely $5\%$.
These assumptions gives the power ratio between SCs and SNRs equal to $P_{\rm
SC}/P_{\rm SNR}\simeq 8\%$. If one account for the uncertainties of several
parameters, such a value ranges between $1\%$ and $30\%$. Given the smaller
power injected by SCs with respect to SNRs, the possibility that they are
responsible for the CR spectrum at PeV energies strongly depends on two
quantities: the diffusion coefficient in the bubble and the value of the
spectral slope. Concerning the former, only a diffusion close to Bohm-like
seems able to produce a maximum energy substantially larger than the one
produced by SNRs. In addition, the spectral slope needs to be harder than the
one produced by SNR, like the case shown in Figure 2(a). If these two
conditions are realized, than SCs could be responsible for the observed CRs at
PeV energies. On the other hand, steeper diffusion and/or slope similar (or
steeper) than the one produced by SNRs would question the role of SCs, at
least in the simplified model presented here.

(a) (b)

Figure 2: CR spectra injected by stellar clusters, compared with the one
injected by SNRs (red-dashed lines). Panels 2(a) and 2(b) show the cases
corresponding to $\chi_{A}=0$ and $\chi_{A}=0.1$, respectively. Different
lines are calculated with different assumptions for the diffusion coefficient
in the wind-bubble, as labelled.

## References

  * [1] C. J. Cesarsky & T. Montmerle, Space Sci. Rev., 36 (1983) 173;
  * [2] M. Casse & J. A. Paul, _ApJ_ , 258 (1982) 860C;
  * [3] F. Aharonian, R. Yang & E. de Ona Wilhelmi, _Nature Astronomy_ 3 (2019) 561;
  * [4] S. Menchiari, PhD thesis (2023) eprint arXiv:2307.03477;
  * [5] P. Blasi & G. Morlino, emphMNRAS 523 (2023) 4015;
  * [6] Z. Cao et al. [LHAASO Coll.], Nature 594 (2021) 33;
  * [7] X. Y. Wang et al. [LHAASO Collaboration], _Chinese Physics C_ 46 (2022) 3;
  * [8] G. Morlino, P. Blasi, E. Peretti & P. Cristofari, _MNRAS_ 504 (2021) 6069;
  * [9] S. Gupta, B. B. Nath, P. Sharma & D. Eichler, _MNRAS_ 493 (2020) 3159;
  * [10] D. V. Badmaev, A. M. Bykov & M. E. Kalyashova, _MNRAS_ 517 (2022) 2818;
  * [11] T. Vieu & B. Reville, _MNRAS_ 519 (2023) 136;
  * [12] S. Menchiari, G. Morlino, E. Amato, N. Bucciantini, PoS(ICRC2023) 649 (2023)
  * [13] T. Cantat-Gaudin et al. [Gaia Coll.], _A &A_ 640 (2020) A1;
  * [14] S. Gabici, F. Aharonian & S. Casanova, _MNRAS_ 396 (2009) 1629;
  * [15] R.-P. Kudritzki & J. Puls , _ARA &A_ 38 (2000), 613;
  * [16] H. Nieuwenhuijzen & C. de Jager, _A &A_ 231 (1990) 134;
  * [17] O. Demircan & G. Kahraman, _Ap &SS_ 181 (1991) 313;
  * [18] P. Kroupa _MNRAS_ 322 (2001)231;
  * [19] R. Weaver, R. McCray, J. Castor, P. Shapiro, R. Moore, _ApJ_ 218 (1977) 377;
  * [20] A. E. Piskunov, A. Just, N. V. Kharchenko, et al. _A &A_, 614 (2018) A22;
  * [21] H. J. G. L. M. Lamers & M. Gieles, _A &A_ 455 (2006) L17;
  * [22] D. R. Lorimer _IAU Symposium 2018_ (2004) preprint arXiv:0308501;
  * [23] P. Cristofari, S. Gabici, T. B. Humensky et al. _MNRAS_ 471 (2017) 201;
  * [24] D. Caprioli, C. Haggerty, & P. Blasi, _ApJ_ 905 (2020) 2;
  * [25] M. A. Miville-Deschenes, N. Murray & E. J. Lee, _ApJ_ 834 (2017) 57;
  * [26] G. Morlino & S. Celli, _MNRAS_ 508 (2021) 6142;
  * [27] P. Cristofai, P. Blasi & E. Amato, _APh_ 123 (2020) 102492;
  * [28] D. Elia, S. Molinari, E. Schisano et al., _ApJ_ 941 (2022) 162;
  * [29] G. A. Tammann, W. Loeffler & A. Schroeder, _ApJS_ 92 (1994) 487;