\section{Introduction}

Recent years have seen an increase in efforts to understand the role of social interaction in social cognition from an embodied perspective \cite{DeJaegher-DiPaolo-2007, DeJaegher-DiPaolo-2008, DeJaegher-DiPaolo-2010, DiPaolo-2000, Froese-Gershenson-Rosenblueth-2013, Froese-2018, Candadai-et-al-2019}. It has been argued that rather than being merely an outcome of the dynamics of individual cognitive agents, social interaction can itself constitute cognition and have an effect on the individuals that partake in it\footnote{This is admittedly a contentious claim, cf. \cite{Herschbach-2012}.}. According to one specific embodied cognition account – enactivism – social interaction is defined as an active co-regulated coupling between two or more autonomous agents, where their role of interactors co-emerges with the interaction itself and their individual cognitive capacities can be reduced or augmented \cite{DeJaegher-DiPaolo-2007, DeJaegher-DiPaolo-2008, DeJaegher-DiPaolo-2010}. In this paper we focus on how different levels of coupling can influence the agent's neural and behavioral complexity.

The types of couplings distinguished in this and previous work are inspired by a well-known experiment from developmental psychology: the “double TV monitor” paradigm \cite{Murray-Trevarthen-1985}. In this experiment, 2-month-old infants interact with their mothers through a live video link. When the live video is replaced with a recorded replay of the previous actions of the mother, the infants become distressed, distracted and upset, suggesting that the reciprocity of the interaction makes a difference. That is, passive social input that is not sensitive to one’s own response is not sufficient for a positive social experience. Simulation studies described here show that it might also not be sufficient for individual cognition.

Based on the previous approaches, we propose to distinguish the following levels of coupling in dyadic interaction:
\begin{itemize}
\item 2-way or bidirectional coupling: Both agents are mutually interacting, e.g. normal interaction as a mother playing with her infant. Other examples are conversations, dancing, collaborative work, etc. \cite{DeJaegher-DiPaolo-2010}.
\item 1-way or unidirectional coupling: Active agent is in the presence of a non-interactive agent, which is showing pre-recorded behavior, e.g. “double TV monitor” experiment.
\item 0-way or no-way coupling: Active agent is not in the presence of a social partner, i.e. is alone.
\end{itemize}

This distinction is not meant to be exhaustive and future adjustments might be required. For instance, a case in which one agent is fully interactive while the other is present but staying still might fall somewhere between a 0-way and a 1-way coupling. Additionally, even within the isolated condition, one could distinguish different ways in which the agent can couple to its physical environment – bidirectionally in a full sensorimotor loop \cite{ORegan-2001} or one-directionally, whereby the agent is passively receiving stimulation from the environment. However, this work focuses on the social interaction scenario and we think it is still useful to think of levels of coupling within this domain and how they might have distinct effects on individual cognitive capacities.

\subsection{Previous work}
Evolutionary robotics (ER) has been used as a scientific tool to study minimal models of cognition. It is a methodology especially suited for embodied cognition because it allows for developing integrated sensorimotor systems which act in close coupling with their environments, i.e. ER takes into account both embodiment and situatedness in how they contribute to the solution of particular cognitive tasks \cite{Harvey-et-al-2005}. The Candadai et al. (2019) model \cite{Candadai-et-al-2019} is a minimal model of cognition in interaction based on ER. This model demonstrates as a proof of concept that social interaction transforms the neural and behavioral dynamics of embodied agents more than what is achievable by agents operating in isolation. In this model, pairs of agents were evolved in a 2-dimensional environment to maximize their individual neural entropy\footnote{It is important to emphasize here that maximizing this measure was not meant to achieve any particular adaptive outcomes. Rather, the point is to establish whether 2D brains can be evolved to become more complex and whether interaction is an enabling condition for this.}. When the agents were evolved together (interaction condition) they exhibited mutually coordinated behavior and higher individual neural entropy compared to when they were evolved alone (isolation condition). Furthermore, when agents that evolved in the interaction condition were tested in the presence of a "ghost" non-interactive partner (ghost condition), they exhibited a loss in neural and behavioral complexity. This condition is analogous to a “double TV monitor experiment” in which infants used to normal 2-way interaction mode are suddenly placed in a 1-way interaction mode.

In our previous work \cite{Resendiz-Benhumea-Froese-2020}, we have replicated the Candadai et al. (2019) model \cite{Candadai-et-al-2019} with a less constrained parameter space of the agents, leading to a broader diversity in the agent genotypes but achieving similar results in terms of neural complexity. We also provided a state-space analysis of the evolved agents’ neural network activity, which showed a single fixed-point attractor. This is consistent with previous ER work that shows that CTRNN controllers are often able to display rich dynamics despite having a single attractor because the attractor landscape is constantly shifting in the agent’s interaction with the environment \cite{Beer-2000, Buckley-2008}.

\subsection{Current work}
The Candadai et al. (2019) model \cite{Candadai-et-al-2019} demonstrated that 2-way coupling leads to higher neural complexity than 1-way coupling and 0-way coupling. However, the specific implementation adopted is open to two central criticisms.

First, it is known that continuous dynamical systems can only exhibit chaotic behavior when they have more than 2 dimensions \cite{Coddington-Levinson-1955}. This is indeed an assumption of the Candadai et al. (2019) model \cite{Candadai-et-al-2019}: the fact that the agents’ 2D brains exhibit chaotic activity is remarkable because it means that interaction can overcome the inherent limitations of their brains. However, a corollary of this is that if agents had brains with more dimensions, they could generate neural and behavioral complexity without interaction. If so, this would limit the result to very constrained systems and tell us nothing about cognition in multi-dimensional actual brains.

Second, it could be argued that a comparison between agents evolved in isolation and agents evolved in interaction is not a strong comparison that would justify far-reaching claims about the role of interaction in cognitive complexity. The isolated agents might be less complex not only because they don’t interact but because they evolve in an overall poorer environment from which they receive no input. It could be that if the environment contained some stimulation, even if not social, the agents could leverage it to create more neural complexity.

In this study, we address these two limitations and investigate whether the Candadai et al. (2019) model \cite{Candadai-et-al-2019} results still hold when we 1) include one more neuron in the internal layer of the neural architecture of the agents, i.e. we construct a 3-neuron model instead of 2-neuron model and 2) we create a new condition in which we evolve agents with a ghost agent, thereby providing them with a source of non-interactive stimulation. As a result, we implement 4 conditions based on 3 levels of coupling distinguished above, as illustrated in Fig.~\ref{fig_1},~\ref{fig_2},~\ref{fig_3}:
\begin{itemize}
\item \textbf{Interactive condition (2-way coupling):} pairs of agents are evolved in bidirectional interaction.
\item \textbf{Ghost-test condition (1-way coupling):} pairs of agents are evolved in bidirectional interaction but tested with a ghost at the end.
\item \textbf{Ghost-evolution condition (1-way coupling):} an active agent is evolved in the presence of a sufficiently complex ghost.
\item \textbf{Isolated condition (0-way coupling):} agents are evolved in isolation.
\end{itemize}

\begin{figure}[htbp]
\centering
\includegraphics[width=1.9in]{2-way_coupling.pdf}
\caption{Two-way coupling: Fully interactive condition, agents interact with each other by sending and receiving acoustic signals.}
\label{fig_1}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=1.9in]{1-way_coupling.pdf}
\caption{One-way coupling: “Live” agent (red) tested or evolved with a “ghost”
agent (blue). Live agent is sending and receiving acoustic signals. Ghost agent
is playing back pre-recorded behavior of previous trials and so sending pre-recorded signals but not receiving anything.}
\label{fig_2}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=0.59in]{0-way_coupling.pdf}
\caption{Zero-way coupling. The “isolated” agent (red agent) is evolved on its
own, therefore, the agent neither receives any input from the environment nor
interacts with other agents.}
\label{fig_3}
\end{figure}

In all conditions agents are evolved to maximize their neural complexity operationalized, for simplicity, as Shannon entropy of neural outputs. However, the maximization of predictive information (PI, \cite{Bialek-1999}), which has been successfully applied in previous self-organizing robotic systems where they exhibited a wider behavioral spectrum \cite{Martius-Der-Ay-2013}, is considered as an alternative approach for future work.

We then compare the level of neural and behavioral complexity in the last generation across all conditions. We find that despite having more powerful brains, agents that evolve in isolation exhibit lower neural and behavioral complexity. However, agents that evolve in the presence of the ghost exhibit an interesting divergence between the complexity of their neural output and behavior: their neural complexity is comparable to that of agents that evolve in interaction while their behavioral complexity is lower. We also investigate the quality of interaction between social conditions and conclude by placing our results in a broader theoretical context of embodied cognition.

\section{Methods}
\subsection{Model Design}
Our model is a replication of the Candadai et al. (2019) model \cite{Candadai-et-al-2019}, except for the increased number of neurons in the internal layer of the neural architecture: three neurons instead of two.
\subsubsection{Agents and environment}
Agent bodies are circular, with a radius of 4 units. Each of them has two acoustic sensors symmetrically placed in the front of the agent at a 45$^{\circ}$ angle to the central axis; an acoustic emitter positioned in the center of the body; and two motors located in the right and left sides of the agent to allow movement in a 2-dimensional empty environment, which is an unlimited arena. Collisions are modeled as point elastic, which means, no change in the agents’ angular velocity (i.e. no friction between them) and zero net effect on their velocity vectors (i.e. energy of the complete system is conserved). This is achieved by exchanging the agents’ velocity vectors, which causes them to bounce off each other without loss of energy \cite{DiPaolo-2000}. Each agent can emit and sense acoustic signals. The strength of the acoustic signal experiences two kinds of attenuation:
\paragraph{Attenuation due to distance} The maximum strength of the emitted signal is exhibited from the source to a distance equal to 2R between the center of the bodies of the agents. The intensity of the signal obtained by each sensor is then calculated by applying the inverse square law using the distance between the sensor and the source.
\paragraph{Attenuation due to “self-shadowing” mechanism}  This attenuation occurs when the emitted signal passes within the body of the sensing agent (as the intensity of the acoustic signal is weakened due to the agent's own embodiment). It is modeled as a scaling factor over the sensory inputs in a range from 0.1 (when the sensors of the sensing agent are diametrically opposite from the source) to 1 (when the sensing agent is facing the source). The equations to calculate the distance that the emitted signal travels within the body of the sensing agent, i.e. the shielded distance (\(D_{sh}\)), are available in the Supplementary Material of the Candadai et al. (2019) model \cite{Candadai-et-al-2019}.

According to the previous points, the sensory input for each sensor of an agent is calculated using the distance between the sensor and the source (applying the inverse square law), and it is then multiplied by the “self-shadowing” attenuating factor which goes linearly from 1 (when \(D_{sh}=0\)) to 0.1 (when \(D_{sh}=2R\)).

\subsubsection{Neural architecture}
The neural architecture of each agent consists of three layers, in our previous work we called them: sensor layer, neuron layer and actuator layer. The main difference from our current model and the Candadai et al. (2019) model \cite{Candadai-et-al-2019} is presented in the neuron layer, where instead of two neurons, we use three neurons, as shown in Fig.~\ref{fig_4}

\paragraph{Sensor layer} In this layer, there are two sensor nodes with a sigmoidal activation function, whose output is given by:

\begin{equation}
o_{s}=g_{s}\sigma(I_{s}+\theta_{s})\label{eq_1}
\end{equation}

where \(gs\) is the sensory gain, \(\sigma(x)=1/(1+e^{-x})\) is the sigmoidal activation function, \(I_{s}\) is the sensory input and \(\theta_{s}\) is the bias. In this layer, both sensor nodes share common gain and bias.

\paragraph{Neuron layer} This layer is modeled as a continuous-time recurrent neural network (CTRNN) \cite{Beer-1995}. In contrast to the Candadai et al. (2019) model \cite{Candadai-et-al-2019}, this layer now consists of three fully recurrently connected neurons, which corresponds to a 3-dimensional dynamical system. Each neuron’s activity is governed by the following state equation:

\begin{equation}
    \tau_{i} \frac{dy_{i}}{dt}=-y_{i}+\sum_{j=1}^{N} w_{ij}\sigma(y_{j}+\theta_{j})+\sum_{s=1}^{2} w_{is}o_{s}
\end{equation}

where \(dy_{i}/d_{t}\) refers to the rate of change of internal state \(y_{i}\) of neuron \(i\) based on a time constant \(\tau_{i}\). The rate of change \(dy_{i}/d_{t}\) depends on the current state of the neuron, the weighted sum of outputs from other internal neurons and the total external input from the sensors. The output of each neuron based on its internal state is given by a sigmoid activation function \(\sigma(y_{i}+\theta_{j})\) where \(\theta_{j}\) refers to the neuron's bias. In this implementation, the three neurons share same time-constant and bias.

\begin{figure}[tbp]
\centering
\includegraphics[width=2.7in]{arch.pdf}
\caption{Neural architecture of the 3-neuron model based on the Candadai et al.
(2019) \cite{Candadai-et-al-2019} model, where the number of neurons in increased to three. In this
approach, the two sensor nodes share common gain and bias; the three neurons
share common time-constant and bias; and, the three actuator nodes share
common gain and bias.}
\label{fig_4}

\end{figure}
\paragraph{Actuator layer} The three internal neurons feed into the actuator layer, where the input to each actuator node is a weighted sum of the outputs of the neuron. The actuator layer contains three actuator nodes, two corresponding to the motors and one corresponding to the acoustic signal emitter. All of them are sigmoidal units with a gain and bias (but no internal state) such that the output of the actuator node \(i\), \(m_{i}\), is given by:

\begin{equation}
    m_{i}=g_{m}\sigma \left( \sum_{n=1}^{N} w_{ni}*o_{n} + \theta_{i} \right)
\end{equation}

where \(o_{n}\) is the output of the neuron, that are weighted by \(w_{ni}\), \(\theta_{i}\) is the bias and \(g_{m}\) is the gain. In this layer, the three actuator nodes share common gain and bias.

Locomotion is managed by the effective control of the two motors. Net linear velocity is given by the average of their corresponding outputs and angular velocity which rotates the agent is given by their difference divided by the radius of the agent.

\subsubsection{Evolutionary optimization}
\paragraph{Fitness function: Neural entropy}
The fitness function for the evolutionary algorithm is the Shannon entropy of neural outputs, i.e. neural entropy, which has been used as a proxy of cognitive complexity. This function does not optimize any particular task. The agents are initialized in a random position in the environment (see below) and allowed to move around over 4 trials, during which their neural activity is recorded. The neural complexity is measured as the Shannon entropy in the three-dimensional time series from the outputs of the three neurons, which are bounded in the range from 0 to 1. The output space is binned with 100 bins along each of the three dimensions, i.e. one million bins in total. Then, a 3-dimensional histogram is created using all the binning data points acquired during the 4 trials. Thus, the Shannon entropy \(H\) of the neural time series is given by:

\begin{equation}
    H=\sum_{i=1}^{100} \sum_{j=1}^{100} \sum_{k=1}^{100} -p_{ijk} log(p_{ijk})
\end{equation}

where the probability of the neural activity in a specific bin \([i,j]\), \(p_{ij}\), is given by the number of data points in that bin divided by the total number of data points. The neural entropy is then normalized to be in the range from 0 to 1 by dividing by the maximum neural entropy, i.e. \(log(100*100*100)\), when all bins are uniformly populated. Therefore, the normalized neural entropy is given by:

\begin{equation}
    \hat{H} = H/log(100*100*100)
\end{equation}

\paragraph{Genetic algorithm}
We used a real-valued genetic algorithm to optimize the parameters of the neural controllers, such as weights, gains, biases and time-constants in order to maximize the neural entropy of the agents. Each agent had 30 parameters, i.e. for N agents, the genotype consisted of 30N parameters that were initially encoded in the range \([-1, 1]\). When building the agents from each genotype to perform the 4 trials, these parameters were scaled appropriately, following the same parameter ranges as in the Candadai et al. (2019) model \cite{Candadai-et-al-2019} such that: for sensor and actuator nodes, the gains were scaled in the range \([1, 5]\) and the biases were scaled in the range \([-3, 3]\); for neuron nodes, the time-constants were scaled in the range \([1, 2]\) and the biases were scaled in the range \([-3, 3]\); all weights were scaled in the range \([-8, 8]\).

The performance evaluation of the agents was obtained according to each of the 4 conditions: interactive condition, ghost-test condition, ghost-evolution condition and isolated condition. The experimental setup for each of the 4 conditions is described in detail in the Experiments section.

For all conditions, after the performance evaluation, we generated a new population by, first, keeping an elite population of the top 4\% solutions as it is and, second, by mutating and crossing over this elite fraction to get the remainder of the solutions. Mutation was obtained by adding a zero-mean Gaussian mutation noise with variance 0.1 to the solutions and, then, crossover was obtained by swapping each parameter between a pair of solutions with a probability of 0.1.

\subsection{Experiments}
Here, we describe the implementation of the 4 conditions introduced in the Introduction:

\subsubsection{Condition 1: Interactive condition (2-way coupling)}
In this condition, we performed 100 independent runs using 96 pairs of agents that were able to interact with each other. Each pair’s agent parameters were encoded in a single genotype subjected to evolutionary search. Initially, the agents were placed at 20 units from each other. For each trial, their relative angle was modified as \([0, \pi/2, \pi, 3\pi/2]\), respectively, where both agents' heading direction was set to the right. The population was evolved for up to 2000 generations to maximize the neural entropy of both agents.

\subsubsection{Condition 2: Ghost-test condition (1-way coupling)}
In this condition, we selected the best pair of agents of the best 10 runs previously obtained in the fully interactive scenario and, then, tested them under a “ghost” condition. Red agent was the “live” agent and blue agent was the “ghost” agent. The “live” agent was able to interact with the “ghost” partner, while, the “ghost” agent was just playing back pre-recorded behavior from the previous trials in a fully interactive scenario. The “live” agent was initially positioned at a different angle (randomly selected from \([0, \pi/2, \pi, 3\pi/2]\) but different to the one it was chosen when the "live" agent was in 2-way interaction) from the “ghost” agent to avoid repeating the behavior of those previous trials. The initial distance between both agents was 20 units. We conducted 4 trials for each pair of agents, as in the fully interactive condition, and measured the normalized neural entropy of the “live” agent.

\subsubsection{Condition 3: Ghost-evolution condition (1-way coupling)}
In this condition, we selected the best blue agent of the best run in fully interactive scenario and used it as the non-interactive or “ghost” agent. We performed 10 independent runs, where only the interactive or “live” agent was evolved in the presence of the “ghost” partner, which was just playing back the pre-recorded behavior in fully interactive scenario (different for the 4 trials). The same playback was used for all runs, all agents and all generations, in order to evolve the agents to respond to a specific set of conditions and see the effects on neural and behavioral complexity through each generation. The population in each run consisted of 96 individuals, where each individual was encoding the parameters of only the “live” agent. For each trial, the "live" agent was placed 20 units from the "ghost" agent and their relative angle was modified as \([0, \pi/2, \pi, 3\pi/2]\), respectively. The population was evolved for up to 2000 generations to maximize the interaction entropy of “live” agents. Evaluation was performed the same way as for Condition 1 with fitness derived only from the “live” agent’s neural entropy.

\subsubsection{Condition 4: Isolated condition (0-way coupling)}
In this condition, we performed 10 independent runs using isolated agents that were not receiving any input and were evolved on their own to maximize their isolation entropy. Red agent was referred as the “isolated” agent. The population in each run consisted of 96 individuals, where each individual was encoding the parameters for only one agent (“isolated” agent). The population was evolved up to 2000 generations.

In all conditions, the agents' initial heading direction was set to the right.

\section{Results}
 Fig.~\ref{fig_5} shows example trajectories and neural activation of the best agents in 4 conditions from the best runs, respectively.
Complex behavior and neural activity can be clearly seen on the fully interactive condition. This complexity seems to be lost in the ghost-test condition, in line with Candadai et al. (2019) model \cite{Candadai-et-al-2019} results. The agent that exhibits complex movement trajectory when interacting with a live partner, starts to literally run in circles when the partner is non-responsive. This is the case even though the agents in the current experiment have a more complex 3-neuron brain that could in principle exhibit chaotic activity. In the isolated condition, the agent shows highly regular behavior and an oscillatory pattern of neural activity, again, in line with the results of the original model and despite a more complex brain. The most interesting for our purposes ghost-evolution condition, in which the agents are evolved in the presence of a ghost partner displays something that seems as an intermediate level of behavioral and neural complexity between the fully interactive and the isolated case.

\begin{figure*}[tbp]
\centering
\includegraphics[width=7.0in]{graph_all.pdf}
\caption{Example plots of movement trajectories and neural activity across 4 conditions. Movement trajectories in red are of a live agent, movement
trajectories in blue are of a live agent in Figure A and of ghost replay in figures B and C. Neural activity plots show output of 3 internal neurons of the live agents.}
\label{fig_5}
\end{figure*}

In order to go beyond intuitions and understand the general pattern of differences between conditions, we run 4 statistical tests comparing the best agents' neural activity and movement trajectories.

Specifically, in the first test we compared the average of the means of neural entropy of 100 agent pairs in interactive condition against the means of the best live agents in other conditions. Given unequal sample size and significantly different variances between conditions, \(F(3, 127)= 21.24\), \(p<.001\), we performed a non-parametric Kruskal-Wallis test, which showed a significant difference in neural entropy between conditions, \(H(3)= 41.7\), \(p<.001\). Focused comparisons of the mean ranks between groups showed that agents in isolated condition had significantly lower neural entropy compared to interactive and ghost-evolution conditions, as expected. However, somewhat surprisingly, the entropy in ghost-evolution condition was significantly higher than in the interactive or ghost-test conditions. This trend can also be seen in Fig.~\ref{fig_6}A.

Next, we obtained a measure of behavioral complexity of the live agents in all conditions. We recorded heading direction angles at each time point of agent trajectories, which resulted in 1D time series. We then computed sample entropy for each such time series. Since there does not seem to be a universally agreed upon measure of behavioral trajectory complexity, we used sample entropy as a measure that has been shown to be an appropriate index of complexity for biological time series more broadly \cite{Richman-2000}. We have informally validated this measure by checking that it reliably distinguishes between fully regular circular and more complex movement patterns\footnote{We originally tried to estimate behavioral complexity using image entropy calculated on generated plots of trajectories. However, this measure delivered inconsistent results, most likely because it does not consider the temporally ordered nature of movement coordinates and instead takes into account only their spatial dispersion.}. The data subjected to statistical analysis were time series sample entropy for all trials of the 10 best runs in each condition. In the interactive condition, only one of live agents was considered. Levene’s test for equality of variances was not significant and there was a significant effect of condition on the level of behavioral entropy, \(F(3,153)= 20.98\), \(p<.001\). Post-hoc Bonferroni-corrected tests showed that isolated agents’ entropy was significantly lower than entropy in all other conditions \((p<.001)\) but also ghost-evolution condition entropy was significantly lower than entropy in the interactive condition (see Fig.~\ref{fig_6}B).

As the difference in neural and behavioral complexity in the ghost-evolution condition compared to interactive condition was in opposite directions (higher and lower respectively), we run two further tests to investigate the type of interaction between the agents in 3 conditions with coupling (thus, excluding the isolated condition). In the first test we have computed the entropy of the distance between the agents while in the second their synchrony. Distance entropy was measured with a binned Shannon entropy approach adopted in the original Candadai et al. (2019) paper \cite{Candadai-et-al-2019} . Synchrony was estimated with a Dynamic Time Warping metric applied to 2-dimensional time series representing the xy-coordinates of the two agents in all trials. Only the first measure showed a significant overall effect of condition, \(F(2, 27)=4.27\), \(p<.05\) and a significant post-hoc pairwise comparison between Ghost-evolution and Ghost-test conditions. Fig.~\ref{fig_6}C and Fig.~\ref{fig_6}D show the overall trends for these measures.

\begin{figure}[htbp]
\centering
\includegraphics[width=3.2in]{all_plots_grid.pdf}
\caption{Individual and interaction measures across experimental conditions.
Figures A and B show neural and behavioral entropy of live agents. C:
distance entropy between agents. D: DTW-estimated distance between
agents (lower distance means more synchrony).}
\label{fig_6}
\end{figure}

\section{Discussion and Conclusions}
In this work, we have extended the Candadai et al. (2019) model \cite{Candadai-et-al-2019} in order to address some of its limitations. By including one more neuron in the neuron layer, i.e. three neurons instead of two, we have found that the results in terms of neural and behavioral complexity are similar to those in the original configuration. In particular, the agents’ neural entropy is still higher in the interactive condition than in the ghost-test or isolated conditions. Additionally, we have quantitatively assessed behavioral complexity in all conditions and this measure was found to be also lower in the isolated condition. This means, perhaps unsurprisingly, that a powerful brain operating in isolation, without any input is not able to achieve high levels of neural or behavioral complexity.

Our second generalization test delivered mixed insights. On the one hand, by observing example plots of neural activation (Fig. 5), it can be suggested that in the interactive condition, the neural activity of the three neurons exhibits more chaotic activity than in the other conditions, including the more stringent ghost-evolution condition. On the other hand, this suggestion is not borne out by the statistical test that shows that neural complexity in this condition is higher than in the interactive condition. This would mean that richer environment that provides constant but complex input to the agent can compensate the relative poverty of the 1-way coupling and lack of contingent response of the interaction partner. At the same time, despite a higher neural complexity, behavioral complexity in this setting is lower. Thus, agents that are evolved to maximize their neural complexity in a rich environment end up behaving in more predictable and regular manner than agents that receive specifically social stimulation.

Briefly, our findings suggests that the richness of the environment may compensate 1-way interaction in terms of neural complexity, however, this does not apply in terms of behavioral complexity. This suggests that the social world (i.e. real-time interaction between agents becoming the whole brain-body-environment-body-brain system \cite{Froese-2013}) allows for a greater repertoire of behaviors transforming our individual capacities \cite{DeJaegher-DiPaolo-2010}.

Some possible limitations may be encountered in our model based on the specific metrics that we used as a measure of neural and behavioral complexity. Further work is needed to compare the results by implementing alternative methods, e.g. evolving agents to maximize predictive information (PI), using permutation entropy and applied it to raw xy-coordinates, etc.

\section{Future Work}
In this paper, we have explored the neural and behavioral complexity of embodied agents using different levels of coupling in dyadic interaction. Future work will investigate how different modes of coupling can affect individual and interactive capacities of evolved agents. Specifically, it could be argued that an opportunity to interact with multiple partners could further enhance individual complexity. Alternatively, allowing agents to use different interactive modalities, such as distal and proximal coupling analogous to pheromone and saliva-based interactions in ants, could enrich their cognition. This will allow us to further understand how individual complexity can be generated by interaction.

\section*{Acknowledgment}
We are grateful to Alexey Yudin for sharing with us useful suggestions for the analysis of the model.

The authors thank the Scientific Computing \& Data Analysis Section (SCDA) of Research Support Division at Okinawa Institute of Science and Technology (OIST) for using their High-Performance Computing resources.

We acknowledge the help from Randall D. Beer by providing the Evolutionary Agents C++ software package v1.2.


\section{Introduction}\label{SectIntro}

Coagulation processes (also known as aggregation, agglomeration or coalescence) are amazingly common in nature. They determine a large number of phenomena, such as cloud formation, blood coagulation, milk curdling, traffic jam formation, and even planetary accretion. For decades, they were extensively investigated in several fields \cite{paper1, paper2, paper3, paper4, paper5, paper6, paper7,paper8,paper9,paper10}. A broad range of interdisciplinary applications can be found, including percolation phenomena in random graphs and complex networks \cite{2005_JPhysALushnikov, 2009_ScienceAchlioptas, 2010_PRLCosta, 2010_PRECho, 2016_PRLCho, 2016_Conv}, population genetics \cite{2005book_Hein}, pattern formation in social \cite{2014_PREMatsoukas, 2014_SciRepMatsoukas}, biological \cite{2007_EcolModelSaadi}, and man--made systems \cite{2002_PREDubovik}. Several technological applications are based on coagulation, including formation of polymers \cite{Stockmayer_1943}, food processing \cite{1981_Schmidt}, material processing \cite{Wattis_2004, Harris_2001}, and water treatment \cite{1993_Edzwald, paper11}. The aggregation formalism has recently started to play a significant role in modeling physiological processes \cite{2020_PRE_Nelson, 2020_Miangolarra}. Another interesting example of a broadly understood aggregation can be neurons accumulating into agglomerations observed in the neural network under specific conditions \cite{2018_ND_Lepek}.

In a coagulation process, particles merge irreversibly in a result of binary collisions (coagulation acts) which can be represented by the general scheme,
\begin{equation} \label{general_scheme}
   \left(i\right)+\left(j\right){{\stackrel{K\left(i,j\right)}{\longrightarrow}}}\left(i+j\right),
\end{equation}

\noindent where $\left(i\right)$ stands for a cluster of mass $i$ and $K\left(i,j\right)$ is the coagulation kernel representing the rate of the process. As the system is closed and clusters grow in time, eventually, all of the particles join to form a single giant cluster.

Deterministic (kinetic) methods for analysing coagulation processes base on the Smoluchowski aggregation equation \cite{paper12, paper13,paper14,paper15,paper16}. Requiring continuous cluster concentrations and infinite size of the system, the  Smoluchowski equation has been used to derive analytical solutions to few particular kernels, e.g., constant, $K(i,j)=const$, additive, $K(i,j) \propto i+j$, and multiplicative (product), $K(i,j) \propto ij$, and for selected initial conditions. This classic method uses a system of coupled nonlinear differential equations and provides mean--field evolution of the cluster size  distribution over time, not giving any information on the deviations from the mean solution. The equation provides an infinite--volume solution for the coagulation process and, as they are normalized with respect to the initial condition, they expire when the system moves away from the initial state. This results are failing for small systems and mature stages of the coagulation process. Therefore, explicit solutions to this equation are not ``exact solutions'' of the coagulation processes. Nevertheless, considerable literature exists on the existence and uniqueness of solutions to some general classes of kernels in Smoluchowski approach, including self--similar dynamical scaling solutions \cite{paper17,paper18,paper19,2004_Menon}. This area is still exhaustively researched \cite{Laurencot_2018, Bonacini_2019, Niethammer_2019, Banasiak_2019_Book, Canizo_2020}. On the basis of the Smoluchowski equation and the works of Ziff \cite{paper14}, a kinetic model of the gelling system consisting of two types of monomers was recently investigated \cite{2020_PRE_Nelson}. The Smoluchowski equation was also used to describe the process of protein aggregation in living cells \cite{2020_Miangolarra}.

The constant and the additive kernels are examples of the non--gelling kernels, in which the sol--gel transition is not observed when monodisperse initial conditions applied (in contrast to the product kernel). These three basic kernels are important as they became reference models of coagulation. The exact solutions to these kernels were obtained for the simplest case of monodisperse initial conditions as solutions to the master equation governing the time evolution of the probability distribution over system states \cite{paper20, 2004_PRLLushnikov, 2005_PRELushnikov, 2006_PhysD_Lushnikov, 2011_JPhysALushnikov} (stochastic approach), or through direct counting of these states \cite{grassberger1, grassberger2, 2018_PREFronczak, 2019_ROMP_Lepek} (combinatorial approach). An exact solution to the product kernel has been found for the case of arbitrary initial conditions \cite{2019_PRE_Fronczak}. The basic kernels were also revisited in a thermodynamic--focused description using the so--called linear ensembles \cite{2018_Matsoukas}.

Recently, the combinatorial approach proved its usefulness for solving more complex form of a kernel under monodisperse initial conditions, i.e., for the linear-chain kernel which included the electrorheological fluid coagulation process \cite{2021_PhysD_Lepek}. The solution involved arbitrary parameter $\alpha$ substantially increasing generality of the kernel form. The accuracy of theoretical solution tested versus numerical simulation and experimental data varied from approximate to ``exact'' depending on the parameter and coagulation time. It has been shown that for $\alpha=0$ the linear--chain kernel can be reduced to the constant kernel and reproduce its previously--known results.

\begin{center}
$\diamondsuit$
\end{center}

The idea behind the combinatorial approach bases on the fact that successive steps of the coagulation process define the space of available states. The probability distribution over the state space can be determined by studying possible growth histories of clusters using combinatorial expressions. Then, the expressions for cluster size distribution and its standard deviation are derived. For that purpose, the model uses incomplete Bell polynomials and some involved combinatorics (please see \cite{2018_PREFronczak, 2019_ROMP_Lepek}). However, such an approach gave significant advantages. The first is an explicit expression for the number of clusters of a given size at a given aggregation step, valid for any kernel. The second advantage is an expression for the standard deviation of that average. It must be emphasized that these expressions does not change when adopted to cover any new type of kernel. The only problem that is needed to be addressed involves transforming a recurrent expression for the number of possible internal states of a cluster of a particular size to a non--recursive form.

The combinatorial framework requires several assumptions: (i) monodisperse initial conditions, (ii) discrete time (i.e., time steps counted by the subsequent states of the system), (iii) one coagulation act occurring at each time step, and (iv) relatively low particle concentrations since we do not consider the possibility of having encounters among 3 particles, or several encounters at a given time.

The condition (ii) needs a thorough comment. In a real coagulating system, particles merge (through binary collisions) more or less frequently. This frequency is determined by the kernel and by the evolution stage of the system. However, an exact time between these acts remains unknown (being stochastic in nature) and is modeled in, e.g., Marcus--Lushnikov framework of aggregation. In contrast, in the combinatorial approach we use, the time is counted in coagulation acts (steps). The first coagulation act occurs at time step $t=1$, the second coagulation act occurs at time step $t=2$, and so on. Therefore, that discrete time is no more than a counter of subsequent states of the system. It must be emphasized that this timescale of binary aggregations acts is not the same as the timescale of a physical time (measured, e.g., in seconds). Binary aggregation time steps are a (nonlinear) rescaled version of the physical time. Due to this reason, the combinatorial approach does not deal with the physical time and does not give information on ``how much time such a particular process would take in real''. The discrete combinatorial time must be regarded as time steps counting subsequent states of the system.

However, it is crucial that the lack of direct relevance to the physical time is not a serious limitation. If we have sufficient information on the real system: an initial number of clusters and a ``snapshot'' of the cluster size distribution in a particular point of time we can translate this real moment of the aggregation process to the language of discrete time steps. It is because the number of clusters in the system clearly defines the time step in the combinatorial formalism. We demonstrated such a procedure comparing combinatorial predictions to the experimental data in \cite{2021_PhysD_Lepek}.

Discrete time counted as subsequent aggregation acts has another important consequence. In general, the kernel, $K$, defines the relative probabilities of particular coagulation acts which are available to occur in the regarded system state (probability distribution of merging acts). In the combinatorial framework, strictly speaking, $K \propto f(i,j)$, so all kernel formulations can be modified by multiplying by any constant number which does not change the reaction rate of the process. For simplicity and consistency with mathematical tools, however, we will use the equal sign notation throughout the work.

\begin{center}
$\diamondsuit$
\end{center}

For this work, using the combinatorial framework with recursive equations \cite{2019_ROMP_Lepek} we have studied two types of irreversible aggregation kernels with arbitrary parameters.

The first one is
\begin{equation} \label{condensation_kernel}
K\left(i,j\right) = (i+A)(j+A),
\end{equation}

\noindent where $A \ge 0$ stands for a parameter which can vary due to the choice of the particular process. According to D. Aldous \cite{paper10}, the kernel of this form represents condensation and branched--chain polymerization processes.

Condensation process is a particular case of coagulation which stands for the change of the physical state of matter from the gas phase into the liquid phase. Condensation has been researched in a number of fields, including (but not limited to) aerosol growth \cite{Drake_1972, Walter_1973, Pruppacher_1978, Wagner_1982}, water vapor condensation in cloud microphysics \cite{Long_1974, Gillespie_1975, Dziekan_2017}, and supersonic flows \cite{Wegener_1954, Yang_2017}.

In turn, branched--chain polymerization is a process of forming branched polymers which are defined as having secondary polymer chains linked to a primary backbone \cite{Stockmayer_1943}. There are several studies on the physical properties of polymers of this kind \cite{Gray_2009, Vilaplana_2010, Marshall_2013}. They are of special interest not only in industrial processing but also in biology \cite{Spouge_1983, Baskakov_2007, Fang_2011}.

The second kernel analysed in this work was the kernel combining constant and additive components in an arbitrary proportion,
\begin{equation} \label{combination_const_additve_kernel_with_B}
K\left(i,j\right) = A+B(i+j),
\end{equation}

\noindent with $A \ge 0$ and $B \ge 0$ being arbitrary parameters defining that proportion.

The kernel~(\ref{combination_const_additve_kernel_with_B}) is a part of the general kernel $ K = A + B (i + j) + Cij $, considered previosuly by Spouge \cite{Spouge_1983b}. Some stochastic solutions for the kernel $ K = A + B (i + j) $ were obtained and it was shown that for infinite systems, they are equivalent to kinetic solutions resulting from the Smoluchowski equation \cite{paper23}. In contrast to the present work, those solutions did not provide any information on the deviation from the mean behaviour, nor were validated against numerical results.

In the combinatorial framework, as mentioned above, an exact formulation of the kernel~(\ref{combination_const_additve_kernel_with_B}) is $K \propto A+B(i+j)$. As such, its right--hand side can be easily normalized to obtain $B=1$, simplifying the formulation. Then, the resulting form of the kernel to be used for further calculations is
\begin{equation} \label{combination_const_additve_kernel}
K\left(i,j\right) = A+i+j.
\end{equation}

For $A=0$, Eq.~(\ref{combination_const_additve_kernel}) simply gives the additive kernel form. In turn, if $A \gg i,j$ then the reaction rate becomes relatively constant and ceases to depend on $i$ and $j$. Such a case correspond to the constant kernel. As the combinatorial approach assumes a finite volume system and monomeric initial conditions, the condition $A \gg i,j$ can also be regarded as a condition for $A$ to be much greater than the initial number of particles (i.e., than the number of monomer units in the system).

The paper is organized as follows. Section \ref{SecBull} presents the essentials of the combinatorial approach. Section \ref{SecDerivations} provides a detailed description of our method for calculating the number of possible internal states of a cluster for the condensation kernel and for the combination of the constant and additive kernels. It includes Lagrange inversion shown in detail. Section \ref{SecComp} compares theoretical predictions to the numerical simulations. Section \ref{SecSum} discusses relation to the system size, gives concluding remarks and describes possible extensions to this work.

\section{Bullet points of combinatorial formalism} \label{SecBull}

It may be beneficial for the Reader if we briefly describe essentials of the combinatorial approach to finite coagulating systems used in this work.

\subsection{Partial Bell polynomials}

A mathematical tool being a significant part of this formalism are partial Bell polynomials (also known as incomplete or the second kind of). They provide detailed information about the partition of an arbitrary set. They are defined as
\begin{multline} \label{Bell_polynom_def}
 B_{N,k}\left(a_1,a_2,\dots ,a_{N-k+1}\right)=B_{N,k}\left(\left\{a_s\right\}\right) \\
 = N!\sum_{\left\{n_s\right\}}{   \frac{a_1^{n_1} a_2^{n_2} \dots}{ n_1!n_2! \dots (1!)^{n_1} (2!)^{n_2} \dots }   } \\ = N!\sum_{\left\{n_s\right\}}{\prod^{N-k+1}_{s=1}{\frac{1}{n_s!}{\left(\frac{a_s}{s!}\right)}^{n_s}}},
\end{multline}

\noindent where the summation is taken over all non--negative integers $\left\{n_s\right\}$ that satisfy
\begin{equation} \label{constraints}
 \sum^N_{s=1}{n_s=k} \;\;\;\; \textrm{and} \;\;\;\; \sum^N_{s=1}{sn_s=N}.
\end{equation}

Eq.~(\ref{Bell_polynom_def}) is a Diophantine equation, hence, in practice, other methods are needed to calculate values of the partial Bell polynomials (please see Appendix for an efficient recurrence equation).

\subsection{System state}

As mentioned in the Introduction, the combinatorial approach requires discrete time steps and monodisperse initial conditions (all of the clusters are monomers of size one). As a single coagulation act occurs at one time step, the total number of clusters, $k$, at time step $t$ is
\begin{equation} \label{EQ2}
 k=N-t,
\end{equation}

\noindent where $N$ is the total number of monomeric units in the system. This number does not change during the coagulation process (preservation of mass) so the state of the system at time $t$ is described by Eq.~(\ref{constraints}). In this context, $n_s\ge 0$ stands for the number of clusters of mass $s$ ($s$ is the number of monomeric units which build the cluster) and $n_1$ corresponds to monomers, $n_2$ to dimers, $n_3$ to trimers, and so on. Therefore, a state of the system at time step $t$ is described by the sequence $\left\{ n_s \right\}$,
\begin{equation} \label{EQ3}
 \mathrm{\Omega }\left(t\right)=\left\{n_1,n_2,\dots ,n_s,\dots ,n_N\right\}.
\end{equation}

\subsection{Average number of clusters of given size and its standard deviation}

Being familiar with partial Bell polynomials, we can proceed to a brief outline of the main expressions arising from the combinatorial approach used here. They are applicable for any form of the kernel. They take their origins in three observations (for details and derivations, please refer to \cite{2018_PREFronczak, 2019_ROMP_Lepek}): (i) the set of monomers can be divided into subsets in a particular number of ways, (ii) coagulation acts for a particular cluster can be distributed in different time steps, and (iii) a given cluster could be created in a particular number of ways (i.e., the number of possible histories of a cluster) and this number must be derived for the kernel of interest. By combining these expressions together, one can calculate the average number of clusters, $\left\langle n_s\right\rangle$, of a given size, $s$, as
\begin{equation} \label{ns_general}
 \left\langle n_s\right\rangle =\binom{N}{s}{\omega }_s\frac{B_{N-s,k-1}\left(\left\{{\omega }_g\right\}\right)}{B_{N,k}\left(\left\{{\omega }_g\right\}\right)},
\end{equation}

\noindent where, for simplicity,
\begin{equation} \label{omega_s}
\omega_s = \frac{x_s}{(s-1)!} \;\;\;\; \textrm{and} \;\;\;\; \left\{\omega_g\right\} =  \left\{ \frac{x_g}{(g-1)!} \right\}.
\end{equation}

Eq.~(\ref{ns_general}) describes the average number of clusters of size $s$ after $t$ steps of the aggregation process (one coagulation acts occurs in one time step). Although $t$ is not explicitly present in the equation, $k$ plays $t$'s role, as $k$ is the total number of clusters in the system and decreases linearly with time.

In Eq.~(\ref{omega_s}), please note that $\omega_s$ is a single value and depends on cluster size $s$, while $\{\omega_g\}$ is a sequence not dependent on $s$, where $g$ varies from $1$ to $N-k+1$ (i.e., to $t+1$). In other words, $g$ is a maximal size of the cluster that \textit{could} arise in the system. The number $x_s$ (or $x_g$) used in Eq.~(\ref{omega_s}) can be interpreted as the number of possible internal states of a cluster of a given size. We will describe $x_s$ in detail further in this section.

As it was mentioned before, not only was the average number of clusters obtained in the combinatorial framework, but so was the corresponding standard deviation of this average,
\begin{equation} \label{std_dev_general}
 {\sigma }_s=\sqrt{\left\langle n_s\left(n_s-1\right)\right\rangle +\left\langle n_s\right\rangle -{\left\langle n_s\right\rangle }^2}
\end{equation}

\noindent where
\begin{equation} \label{std_dev_general_addition}
 \left\langle n_s\left(n_s-1\right)\right\rangle =\binom{N}{s,s}{{\omega }_s}^2\frac{B_{N-2s,k-2}\left(\left\{{\omega }_g\right\}\right)}{B_{N,k}\left(\left\{{\omega }_g\right\}\right)}\
\end{equation}

\noindent for $2s \leqslant N$, and $\left\langle n_s\left(n_s-1\right)\right\rangle = 0$ for other cases. For short, $\binom{N}{s,s}=\binom{N}{s}\binom{N-s}{s}$.

\subsection{Number of internal states of cluster of given size}

The last expression to be presented in this section is the expression for the number of internal states of a cluster of a given size, $x_s$ (or $x_g$). It can also be interpreted as the number of possible ways to create such a cluster (the number of its ''available histories''). Here, we will use subscript $g$ to emphasize that we study a single cluster of size $g$ that \textit{could} grow in the system. In previous works \cite{2019_ROMP_Lepek,2021_PhysD_Lepek}, it has been shown that this number can be calculated as
\begin{equation} \label{xg_recursive_definition_general}
 x_g = \frac{1}{2}\sum^{g-1}_{h=1}\binom{g}{h}\binom{g-2}{h-1}x_hx_{g-h}K(g,h),
\end{equation}

\noindent where $x_h$ and $x_{g-h}$ represent the numbers of ways to create two clusters of size $h$ and $(g-h)$ that joined and became a cluster of size $g$. $K(g,h)$ is an arbitrary kernel of our interest (translated into variables $g$ and $h$). The above recurrence expression is build as follows.

When the two clusters merge the resulting cluster of size $g$ appears. We can divide the cluster of size $g$ into subclusters of size $h$ and size $(g-h)$ in exactly $\binom{g}{h}$ ways. In other words, the first binomial factor denotes the number of ways of choosing a cluster of size $h$ out of $g$ monomers.

The second binomial factor, $\binom{g-2}{h-1}$, covers the fact that coagulation acts resulting in clusters $h$ and $(g-h)$ could occupy different possible time steps. We will use an example. There are two clusters: cluster of size $(g-h)=3$ and cluster of size $h=4$. To create these clusters we need five time steps. The cluster $(g-h)$ needs two time steps to arise, while the cluster $h$ needs three time steps to arise. Coagulation acts related to the creation of the cluster of size $h$ occupy $(h-1)$ time steps out of the total number of $(g-2)$ time steps needed to create clusters of sizes $h$ and $(g-h)$, hence, the binomial factor. Therefore, coagulation acts for the $(g-h)$ cluster and coagulation acts for the $h$ cluster could happen in ten different orders (sequences) in time, as $\binom{g-2}{h-1}=10$.

The sum in Eq.~(\ref{xg_recursive_definition_general}) is taken over all possible pairs of clusters that can join to become cluster of size $g$. The factor of $\frac{1}{2}$ is used to prevent double counting of the coagulation acts.

In the next sections, we will modify Eq.~(\ref{xg_recursive_definition_general}) to describe kernels to be solved and transform it to a non--recursive form that can be used in Eqs.~(\ref{ns_general})--(\ref{omega_s}) to calculate $\left\langle n_s\right\rangle$.

\section{Deriving solutions} \label{SecDerivations}

In this Section, we will derive the number $x_g$ for the kernels defined by Eqs.~(\ref{condensation_kernel}) and~(\ref{combination_const_additve_kernel}).

\subsection{Condensation kernel} \label{subsection_condensation}

At this point, we must translate the condensation kernel, $K\left(i,j\right) = \left( i+A \right) \left( j+A \right)$, into variables $g$ and $h$ used in Eq.~(\ref{xg_recursive_definition_general}). Parameters $i$ and $j$ are the masses of two merging clusters, $g$ is the mass of the resulting cluster, and $h$ is the mass of one of the subclusters (e.g., $i=h$). The condensation kernel can be rewritten as
\begin{equation} \label{K_gh_cond}
 K\left(g,h\right) = \left( h+A \right) \left( g-h+A \right).
\end{equation}

Then, the recursive expression for $x_g$ for the condensation kernel is
\begin{equation} \label{xg_recursive_condensation}
 x_g=\frac{1}{2}\sum^{g-1}_{h=1}\binom{g}{h}\binom{g-2}{h-1}x_hx_{g-h} \left( h+A \right) \left( g-h+A \right).
\end{equation}

The strategy of solving Eq.~(\ref{xg_recursive_condensation}) to obtain a non--recurrent expression for $x_g$ is based on the generating function method and is similar as in previous case of the simple kernels \cite{2019_ROMP_Lepek}.

Expanding binomial coefficients and substituting
\begin{equation} \label{condensation_subst}
y_g=\frac{x_g(g+A)}{g!(g-1)!},
\end{equation}

we obtain
\begin{equation} \label{condensation_2}
 \frac{\left(g-1\right)}{(g+A)}y_g=\frac{1}{2}\sum^{g-1}_{h=1}{y_hy_{g-h}}.
\end{equation}

Now, we multiply both sides of Eq.~(\ref{condensation_2}) by $\sum^{\infty }_{g=1}{z^g}$ to have
\begin{equation} \label{condensation_3}
 \sum^{\infty }_{g=1}{\frac{\left(g-1\right)}{(A+g)}y_gz^g}=\frac{1}{2}\sum^{\infty }_{g=1}{\sum^{g-1}_{h=1}{\left(y_hz^h\right)\left(y_{g-h}z^{g-h}\right)}}.
\end{equation}

The generating function is defined as $G\left(z\right)\equiv \sum^{\infty }_{g=1}{y_gz^g}$. Therefore, the right--hand side of Eq.~(\ref{condensation_3}) is equal to  $\frac{1}{2}G\left(z\right)G\left(z\right)$.

At this point, we must transform the left--hand side of Eq.~(\ref{condensation_3}) to the form containing $G\left(z\right)$. These transformations give
\begin{multline} \label{condensation_4}
 \sum^{\infty }_{g=1}{\frac{\left(g-1\right)}{(g+A)}y_gz^g} =
 \sum^{\infty }_{g=1}{{\frac{g}{g+A}}y_gz^g}-\sum^{\infty }_{g=1}{\frac{1}{g+A}y_gz^g} \\
 = \frac{1}{z^A}\sum^{\infty }_{g=1}{\frac{g}{g+A}y_gz^{g+A}}-\frac{1}{z^A}\sum^{\infty }_{g=1}{\frac{1}{g+A}y_gz^{g+A}} \\ = \frac{1}{z^A}\sum^{\infty }_{g=1}{gy_g \int z^{g+A-1} dz}-\frac{1}{z^A}\sum^{\infty }_{g=1}{y_g \int z^{g+A-1} dz} \\
 = \frac{1}{z^A}\int\sum^{\infty }_{g=1}{gy_g z^{g+A-1} dz}-\frac{1}{z^A}\int\sum^{\infty }_{g=1}{y_g z^{g+A-1} dz} \\
  = \frac{1}{z^A}\int z^A \sum^{\infty }_{g=1}{gy_g z^{g-1} dz}-\frac{1}{z^A} \int z^{A-1} \sum^{\infty }_{g=1}{y_g z^g}  dz \\
  = \frac{1}{z^A}\int z^A \sum^{\infty }_{g=1}{y_g \frac{\partial}{\partial z} \left( z^{g} \right) dz}-\frac{1}{z^A} \int z^{A-1} G(z) dz \\
  = \frac{1}{z^A}\int z^A  \frac{\partial}{\partial z} G(z) dz -\frac{1}{z^A} \int z^{A-1} G(z) dz.
 \end{multline}

As we transformed both sides of Eq.~(\ref{condensation_3}), it can be rewritten as
\begin{equation} \label{condensation_5}
\int z^A \frac{\partial }{\partial z} G(z) dz - \int z^{A-1} G(z) dz = \frac{1}{2}z^A G^2(z).
\end{equation}

In the next two equations, we will use $G$ instead of $G(z)$ for short. Differentiating both sides with respect to $z$ and after elementary transformations, we have
\begin{equation} \label{condensation_6}
\frac{ \partial G}{\partial z} = \frac{1}{z} \frac{\frac{A}{2}G^2+G}{1-G}
\end{equation}

\noindent which is an ordinary differential equation with separated variables. It can be solved for an arbitrary $A$ and its solution has the form of
\begin{equation} \label{condensation_7}
\frac{G}{\left(AG+2\right)^{\frac{A+2}{A}}} = Cz.
\end{equation}

Eq.~(\ref{condensation_7}) defines implicitly the generating function $G$. However, for further considerations, we need to have an explicit series representation of $G$. It can be derived using the Lagrange inversion method (see p.~148 in \cite{paper1}). Eq.~(\ref{condensation_7}) has a form of $f\left(G\right)=F$ where $f\left(G\right)$ stands for its left--hand side and $F$ stands for its right--hand side. Applying the Lagrange inversion to the equation of that form, the series representation of the inverse function $G\left(F\right)$ can be derived. Now, we will apply such a procedure to our problem.

In the Lagrange inversion, a general form of the series representation of the inverse function is given by
\begin{equation} \label{lagrange_inversion_general}
G(F) = a + \sum^{\infty}_{n=1}{G_n \frac{ \left( F-f(a) \right)^n }{ n! } },
\end{equation}

\noindent where
\begin{equation} \label{lagrange_inversion_Gn}
G_n = \lim_{w \to a} \left[ \frac{d^{n-1}}{dw^{n-1}}  \left( \frac{w-a}{f(w)-f(a)} \right)^n  \right].
\end{equation}

Taking arbitrary $a=0$, we have $f(a=0)=0$, thus, considering Eq.~(\ref{condensation_7}), we can write
\begin{equation} \label{lagrange_inversion_a0}
G(F) = \sum^{\infty}_{n=1}{G_n \frac{F^n }{ n! } }
\end{equation}

\noindent and
\begin{equation} \label{lagrange_inversion_Gn_ours}
G_n = \lim_{w \to 0} \left[ \frac{d^{n-1}}{dw^{n-1}}  \left( Aw+2 \right)^{\frac{2+A}{A}n}  \right].
\end{equation}

The technique for finding the sequence $G_n$ is to calculate values of Eq.~(\ref{lagrange_inversion_Gn_ours}) for the first few numbers $n$, e.g., $n=1,2,3,4$. Then, basing on these results, one shall guess the expression for $G_n$. Following this scheme, in our case, we can describe the sequence as
\begin{equation} \label{Gn_ours_guessed}
G_n = 2^{\frac{2n+A}{A}} [ 2n+nA ]^{*},
\end{equation}

\noindent where asterisk denotes a special product,
\begin{equation} \label{asterisk_function_cond}
[2n+nA]^{*}=\left\{
\begin{array}{ccc}
1 & \mbox{for} & n=1, \\
\prod^{n}_{m=2}{(2n+mA)} & \mbox{for} & n \ge 2.
\end{array}
\right.
\end{equation}

Finally, from the Lagrange inversion applied to Eq.~(\ref{condensation_7}), we obtained
\begin{equation} \label{lagrange_inversion_final}
G(z) = \sum^{\infty}_{g=1}{G_g \frac{ \left( Cz \right)^g }{ g! } },
\end{equation}

\noindent where $n$ was replaced by $g$ to stay consistent with our notation and $G_g$ is given by Eq.~(\ref{Gn_ours_guessed}).

Until now, the constant $C$ remained unknown. At this point, we will calculate its value. Bearing in mind the definition of the generating function, $G(z) = \sum^{\infty }_{g=1}{y_gz^g}$, and considering Eq.~(\ref{lagrange_inversion_final}), we can write
\begin{equation} \label{comparison_to_definition}
\sum^{\infty}_{g=1}{G_g \frac{ \left( Cz \right)^g }{ g! } } = \sum^{\infty }_{g=1}{y_gz^g}.
\end{equation}

Comparing the first elements of these sums we have $G_1 \frac{ Cz }{ 1!} = {y_1 z}$. Therefore,
\begin{equation} \label{comparison_first_coefficients}
G_1 C = y_1.
\end{equation}

From the initial condition, $x_1 \equiv 1$, and considering the substitution, Eq.~(\ref{condensation_subst}), we obtain $y_1=A+1$. On the other hand, $G_1$ is given by Eq.~(\ref{Gn_ours_guessed}). In consequence,
\begin{equation} \label{constant_C}
C = 2^{-\frac{A+2}{A}}(A+1).
\end{equation}

Now, considering the respective coefficients in Eq.~(\ref{comparison_to_definition}), $G_g \frac{ C^g }{ g! } = {y_g}$, and after elementary transformations, we have
\begin{equation} \label{yg_final_cond}
y_g = \frac{(A+1)^g}{2^{g-1}g!} [2g+gA]^{*} .
\end{equation}

The last step to derive $x_g$ is going back from the substitution, Eq.~(\ref{condensation_subst}). Finally, we obtain the number of ways to create a cluster of size $g$ for the condensation kernel as
\begin{equation} \label{xg_final_cond}
x_g = \frac{(A+1)^g (g-1)!}{2^{g-1}(A+g)} [2g+gA]^{*},
\end{equation}

\noindent where the special function, $[ . ]^*$, is given by Eq.~(\ref{asterisk_function_cond}).

Eq.~(\ref{xg_final_cond}) together with Eqs.~(\ref{ns_general}) and~(\ref{omega_s}) fully defines the average number of clusters of a given size at a given time step, $\langle n_s \rangle$, for the condensation kernel.

\subsection{Combination of constant and additive kernels}

In this Section, we will derive the number of possible internal states of a cluster of given size $g$, $x_g$, for the kernel (\ref{combination_const_additve_kernel}).

Similarly as in the case of the condensation kernel, we must translate the kernel of interest, this time $K\left(i,j\right) =  A+i+j$, to the language of $g$ and $h$ used in the general expression for $x_g$, Eq.~(\ref{xg_recursive_definition_general}). Bearing in mind that $i$ and $j$ are the masses of two merging clusters, $g$ is the mass of the resulting cluster and $h$ is the mass of one of the sub-clusters (e.g., $i=h$), we can rewrite the kernel as
\begin{equation} \label{K_gh_comb}
 K\left(g,h\right) = A + h + (g-h) = A + g.
\end{equation}

Hence, using Eq.~(\ref{xg_recursive_definition_general}) the recurrent expression for $x_g$ for this kernel is
\begin{equation} \label{xg_recurrent}
 x_g=\frac{1}{2}\sum^{g-1}_{h=1}\binom{g}{h}\binom{g-2}{h-1}x_hx_{g-h} \left( A + g \right).
\end{equation}

Expanding binomial factors and substituting
\begin{equation} \label{substitution} y_g = \frac{x_g}{g!(g-1)!},
\end{equation}

\noindent we obtain
\begin{equation} \label{combination_2}
 \frac{\left(g-1\right)}{(g+A)}y_g=\frac{1}{2}\sum^{g-1}_{h=1}{y_hy_{g-h}}.
\end{equation}

Eq.~(\ref{combination_2}) has the same form as the analogous equation for the condensation kernel with the only difference in the substitution, Eq.~(\ref{substitution}). Thus, we will use the solution from \ref{subsection_condensation}. However, because of the substitution we need to calculate $C$ once again. Now, from the initial condition ($x_1 \equiv 1$), we have $y_1=1$, so the constant is
\begin{equation} \label{constant_C_again}
C = 2^{-\frac{A+2}{A}}.
\end{equation}

Therefore,
\begin{equation} \label{yg_final_comb}
y_g = \frac{1}{2^{g-1}g!} [2g+gA]^{*},
\end{equation}

\noindent where the asterisk denotes a special product given by Eq.~(\ref{asterisk_function_cond}).

Finally, going back from the substitution, Eq.~(\ref{substitution}), we obtain the number of ways to create a cluster of size $g$ for the kernel (\ref{combination_const_additve_kernel}) as
\begin{equation} \label{xg_final_comb}
x_g = \frac{(g-1)!}{2^{g-1}} [2g+gA]^{*},
\end{equation}

\noindent where, again, the product function, $[ . ]^*$, is given by Eq.~(\ref{asterisk_function_cond}).

Eq.~(\ref{xg_final_comb}) together with Eqs.~(\ref{ns_general}) and~(\ref{omega_s}) fully defines the average number of clusters of a given size at a given time step, $\langle n_s \rangle$, for the combination of the constant and additive kernels.

\section{Theoretical results compared to numerical simulations} \label{SecComp}

\begin{figure*}[ht!]
\includegraphics[width=1\textwidth]{condensation_Fig.png}
\caption{Theoretical calculations of the average number of clusters of size $s$, $\langle n_s \rangle$, for the condensation kernel, $K(i,j)=(A+i)(A+j)$, with different values of $A=\{0, 3, 10, 30, 100, 10^6 \}$ compared to the numerical simulation results. Solid lines represent theoretical predictions based on the combinatorial equations, Eqs.~(\ref{ns_general_final})--(\ref{omega_g_condensation_final}).  Although theoretical results are defined only for integer $s$, solid lines are used as guidelines for eyes. Circles, squares and triangles correspond to the results obtained by numerical simulation. The number of monomers in the system was $N=200$. Three stages of the coagulation process are presented: $t=20$ (beginning of the process, $\nicefrac{t}{N}=0.1$, squares, black), $t=100$ (the half-time, $\nicefrac{t}{N}=0.5$, circles, red) and $t=180$ (late stage of the process $\nicefrac{t}{N}=0.9$, triangles, blue). For $A=0$, the results (both theoretical and numerical) fully correspond to the results obtained for the product kernel. For $A \gg N$ (here, $A=10^6$), the results correspond accurately to the known results for the constant kernel. For each case, ${10}^5$ independent simulations were performed.}
\label{Figure_condensation}
\end{figure*}

\begin{figure*}[ht!]
\includegraphics[width=1\textwidth]{combination_Fig.png}
\caption{Theoretical calculations of the average number of clusters of a given size, $\langle n_s \rangle$, for the combination of constant and additive kernels, $K(i,j)=A+i+j$, with several values of $A=\{0,3,10,30,100,10^6\}$, compared to the numerical simulation results. Solid lines represent theoretical predictions based on the combinatorial expressions, Eq.~(\ref{ns_general_final_comb})--(\ref{omega_g_combination_final}). Although theoretical results are defined only for integer $s$, solid lines are used as guidelines for eyes. Circles, squares and triangles correspond to the results obtained by the numerical simulation. The number of monomers in the system was $N=100$. Three stages of the coagulation process are presented: $t=30$ (beginning stage of the process, $\nicefrac{t}{N}=0.3$, squares, black), $t=70$ (moderate time, $\nicefrac{t}{N}=0.7$, circles, red) and $t=95$ (a very late stage of the process, $\nicefrac{t}{N}=0.95$, triangles, blue). Subsequent values of $A$ represent different proportions of the constant and additive parts involved in the kernel. For $A=0$, the results (both theoretical and numerical) fully correspond to the results obtained previously for the additive kernel. For $A \gg N$ (here, $A=10^6$), the results correspond accurately to the known results for the constant kernel. For each case, ${10}^5$ independent simulations were performed.}
\label{Figure_combination}
\end{figure*}

\begin{figure*}[ht!]
\includegraphics[width=1\textwidth]{stddev_Fig.png}
\caption{Standard deviation predicted by the theoretical model, Eqs.~(\ref{std_dev_general}) and~(\ref{std_dev_general_addition}), vs. standard deviation obtained by the simulation for the condensation kernel and for the combination of the constant and additive kernels. For both, $A=10$. Solid and dashed lines represent combinatorial results for the average number of clusters of a given size, $\langle n_s\rangle$, and for its standard deviation, respectively; squares and circles represent, respectively, the average number of clusters and standard deviation obtained by simulation. Left column: the condensation kernel, $N=200$, (a) $t=20$, $\nicefrac{t}{N}=0.1$, (b) $t=100$, $\nicefrac{t}{N}=0.5$, and (c) $t=180$, $\nicefrac{t}{N}=0.9$. Right column: the combination of the constant and additive kernels (for this case, $N=100$), (d) $t=30$, $\nicefrac{t}{N}=0.3$, (e) $t=70$, $\nicefrac{t}{N}=0.7$, and (f) $t=95$, $\nicefrac{t}{N}=0.95$. Upper inset figures: plots of variance divided by mean, $ \sigma^2/ \langle n_g\rangle$. Again, solid lines for theoretical predictions and circles for numerical data. Lower inset figures: plots of $\langle n_s\rangle$ for the highest $s$ (logarithmic scale). For each case, ${10}^5$ independent simulations were performed.}
\label{Figure_stddev}
\end{figure*}

In this Section, we will compare theoretical solutions arising from the combinatorial framework which were derived in Section~\ref{SecDerivations} to numerical results obtained by the simulation.

The numerical simulation term refers to a simulation of the aggregating system of $N$ particles where the information on the number of clusters and on the masses of these clusters are stored in the program memory and updated step by step in order with subsequent aggregation events. For each aggregation event, two clusters are selected randomly from the distribution suitable for the kernel which is simulated. The crucial part of the simulation is implementing the reaction kernel, $K$. In the previous work on the linear--chain kernel, an exhaustive description on how to implement an arbitrary kernel in such a simulation has been given (Appendix in \cite{2021_PhysD_Lepek}). The code used for the simulations, as well as for theoretical calculations, is available in the Internet \cite{cpp_libraries}.

\subsection{Final expressions for condensation kernel}

To have a final expression for the average number of clusters of a given size at a given time step, $\left< n_s \right>$, we need to use Eqs.~(\ref{ns_general}) and~(\ref{omega_s}) with $x_g$ specified by Eqs.~(\ref{xg_final_cond}) and~(\ref{asterisk_function_cond}). In effect, we obtain a set of self--contained expressions to calculate $\left< n_s \right>$ for the condensation kernel, $K(i,j)=(A+i)(A+j)$. For clarity, we will write these expressions explicitly,
\begin{equation} \label{ns_general_final}
 \left\langle n_s\right\rangle =\binom{N}{s}{\omega }_s\frac{B_{N-s,k-1}\left(\left\{{\omega }_g\right\}\right)}{B_{N,k}\left(\left\{{\omega }_g\right\}\right)},
\end{equation}

\noindent where the time is counted by the total number of clusters remaining in the system, $k=N-t$, and
\begin{equation} \label{omega_s_condensation_final}
\omega_s = \frac{(A+1)^s}{2^{s-1}(A+s)} [2s+sA]^{*},
\end{equation}
\begin{equation} \label{omega_g_condensation_final}
\left\{\omega_g\right\} =  \left\{ \frac{(A+1)^g}{2^{g-1}(A+g)} [2g+gA]^{*} \right\}
\end{equation}

\noindent and $[ . ]^*$ denotes the special product function given by
\begin{equation} \label{asterisk_function_cond_final}
[2g+gA]^{*}=\left\{
\begin{array}{ccc}
1 & \mbox{for} & g=1, \\
\prod^{g}_{h=2}{(2g+hA)} & \mbox{for} & g \ge 2.
\end{array}
\right.
\end{equation}

Using the above set of equations we can calculate cluster size distribution for any stage of the system's evolution. The only parameter here is an arbitrary real value $A \ge 0$.

Taking advantage of a known relation for Bell polynomials \cite{Comtet_1974},
\begin{multline} \label{Bell_relation_1}
B_{N,k}\left(abz_1,ab^2z_2,\dots,ab^{N-k+1}z_{N-k+1}\right) \\ =
 a^kb^NB_{N,k}\left(z_1,z_2,\dots,z_{N-k+1}\right),
\end{multline}

\noindent we can simplify the set of Eqs.~(\ref{ns_general_final})--(\ref{omega_g_condensation_final}) to the form of
\begin{equation} \label{ns_general_final_simplified}
 \left\langle n_s\right\rangle =\binom{N}{s} \frac{[2s+sA]^{*}}{(A+s)} \frac{B_{N-s,k-1}\left(\left\{ \frac{[2g+gA]^{*}}{A+g} \right\}\right)}{B_{N,k}\left(\left\{\frac{[2g+gA]^{*}}{A+g} \right\}\right)}.
\end{equation}

In Figure~\ref{Figure_condensation}, the theoretical solutions arising from our combinatorial equations were plotted against the results of the numerical simulation for several values of $A=\{0,3,10,30,100,10^6\}$ to show the performance of the solutions in respect to $A$. The total number of monomers in the system was $N=200$. For each $A$, we present three stages of the aggregation process. These are: $t=20$ (when the aggregation begins, $\nicefrac{t}{N}=0.1$), $t=100$ (half--time, $\nicefrac{t}{N}=0.5$), and $t=180$ (a late stage of the process, $\nicefrac{t}{N}=0.9$). The case of $A=0$ corresponds to the product kernel which is the widely--known example of a gelling kernel. It can be observed that, for $t=20$ and $t=100$, the theoretical predictions model the results of the simulation with a high precision. For the latest stage of the process, this accuracy significantly falls with the theoretical curve only preserving the giant cluster appearing in the system. The results for $A=0$, both theoretical and numerical, fully correspond to the results obtained previously for the product kernel in \cite{2019_ROMP_Lepek}. In the case of $A=3$ and the latest time, $t=180$, the theoretical solution behaves similarly to the numerical data, thus, the combinatorial solution can be regarded as approximate. For $A=3$ and earlier times, the theoretical results are perfectly accurate. As well, for higher values of $A$ the coagulating system is modeled with an excellent precision for any time of the process. For $A \gg i,j$, the solution (so the numerical simulation) converges to the case of the constant kernel.

Additionally, using Eqs.~(\ref{std_dev_general}) and~(\ref{std_dev_general_addition}) we can easily obtain standard deviation of the average numbers of clusters. For the condensation kernel, these theoretical predictions of the standard deviation were compared to the numerically--obtained results in Figure~\ref{Figure_stddev}.

\subsection{Final expressions for combination of constant and additive kernels}

Again, to obtain a final expression for the average number of clusters of a given size, $\left< n_s \right>$, we need to use Eqs.~(\ref{ns_general}) and~(\ref{omega_s}) with $x_g$ specified by Eqs.~(\ref{xg_final_comb}) and~(\ref{asterisk_function_cond}). Eventually, we have a set of self--contained expressions for the kernel of the form $K(i,j)=A+i+j$. Once again, let us write these expressions explicitly,
\begin{equation} \label{ns_general_final_comb}
 \left\langle n_s\right\rangle =\binom{N}{s}{\omega }_s\frac{B_{N-s,k-1}\left(\left\{{\omega }_g\right\}\right)}{B_{N,k}\left(\left\{{\omega }_g\right\}\right)},
\end{equation}

\noindent where
\begin{equation} \label{omega_s_combination_final}
\omega_s = \frac{1}{2^{s-1}} [2s+sA]^{*},
\end{equation}
\begin{equation} \label{omega_g_combination_final}
\left\{\omega_g\right\} =  \left\{ \frac{1}{2^{g-1}} [2g+gA]^{*} \right\},
\end{equation}

\noindent and $[ . ]^*$ denotes, again, the special function given by Eq.~(\ref{asterisk_function_cond_final}). Analogously as for the condensation kernel, we can calculate cluster size distribution for an arbitrary real parameter $A \ge 0$ at arbitrary time of the system evolution.

Again, using the relation (\ref{Bell_relation_1}), the set of Eqs.~(\ref{ns_general_final_comb})--(\ref{omega_g_combination_final}) can be simplified to the form of
\begin{equation} \label{ns_general_final_comb_simplified}
 \left\langle n_s\right\rangle =\binom{N}{s} [2s+sA]^{*} \frac{B_{N-s,k-1}\left(\left\{ [2g+gA]^{*} \right\}\right)} {B_{N,k}\left(\left\{[2g+gA]^{*} \right\}\right)}.
\end{equation}

In Figure~\ref{Figure_combination}, average numbers of clusters of a given size, $\left< n_s \right>$, are presented for several values of $A=\{0,3,10,30,100,10^6\}$. There are three data series contained in each plot, corresponding to three different stages of the system evolution. The initial number of monomers, thence, the total number of monomeric units in the system was $N=100$. The first series ($t=30$, $\nicefrac{t}{N}=0.3$) corresponds to the early stage of the coagulation process. The second series ($t=70$, $\nicefrac{t}{N}=0.7$) corresponds to the later stage of the process. The third series presents the cluster size distribution at the end of the process ($t=95$, $\nicefrac{t}{N}=0.95$), only five steps before the moment when all of the particles join into one single cluster. For $A=0$, the results (both theoretical and numerical) fully correspond to the results obtained previously for the additive kernel \cite{2019_ROMP_Lepek}. For $A \gg N$, the process can be regarded as usual constant kernel coagulation. In Figure~\ref{Figure_combination}, we have used $A=10^6$ and the results correspond accurately to the known result for the constant kernel. For all values of $A$, as well as for any time, the combinatorial predictions proved to be ``exact'' solutions of the process.

In Figure~\ref{Figure_stddev}, we present standard deviation of the cluster size distribution given by the combinatorial approach, Eqs.~(\ref{std_dev_general}) and~(\ref{std_dev_general_addition}), for $A=10$ and several stages of the system evolution, $t=30,70,95$. We compared them to the standard deviation calculated for the data obtained by the simulation. Combinatorial estimates followed numerical deviation with perfect precision for any time of the process.

\begin{figure}[ht]
\includegraphics[scale=0.6]{sizedep_Fig.png}
\caption{Absolute difference between theoretical predictions and numerical data normalized by system size, $\Delta_N = \frac{1}{N} \sum_s{ | \langle n_s \rangle^{theory} - \langle n_s \rangle^{simul} | }$, for $N = \{  2^5, 2^6, 2^7, 2^8, 2^9 \}$. Blue squares represent condensation kernel, red triangles represent combination of the constant and additive kernels. Dashed lines are only guidelines for eyes. For both kernels, $A=10$ and $\nicefrac{t}{N}=0.5$ were used. For each case, $10^5$ independent simulations were performed, excepting $N=2^9$ when it was $10^4$.}
\label{Figure_sizedep}
\end{figure}

\begin{figure*}[ht]
\includegraphics[width=0.85\textwidth]{summary_Fig.pdf}
\caption{Some of the kernel forms solved in the combinatorial approach can be reduced to a form corresponding to the constant kernel. These are following kernels: previously solved linear--chain kernel, $ K = \left(i^{-1} + j^{-1} \right)^{\alpha} $, condensation kernel, $ K=(A+i)(A+j) $, and a linear combination of the constant and additive kernels, $ K = A + i + j $. In the case of the latter two, for the constant value $ A =0 $, one obtains previously known kernels, respectively, product and additive. In this way, theoretical solutions were obtained for all "intermediate" processes between the constant kernel and the product, additive, and linear--chain kernels.}
\label{Figure_Summary}
\end{figure*}

\section{Discussion and conclusion} \label{SecSum}

For this work, using combinatorial framework with recursive equations we have studied the two types of irreversible aggregation kernels with arbitrary parameters, the condensation kernel, $K(i,j) = (A+i)(A+j)$, and the combination of the constant and additive kernels, $K(i,j)=A+i+j$. Expressions for the average number of clusters of a given size and for the corresponding standard deviation were obtained and plotted against the numerical simulation results.

What are the advantages of the above analysis over a conventional solutions of the rate equations? As far as we know there were no “exact” solution of the rate equation with the condensation kernel. On the other hand, some stochastic solutions for the kernel $ K = A + B (i + j) $ were obtained \cite{Spouge_1983b} and it was shown that, for infinite systems, they are equivalent to kinetic solutions of the Smoluchowski equation \cite{paper23}. In contrast to the present work, they did not provide any information on the deviation from the mean behaviour, nor were validated by the numerical results. Additionally, it is worth to note that the discreteness of the combinatorial approach does not diminish the
generality of the solutions as it allows to obtain continuous--time results by the asymptotic analysis (see, e.g., Eqs.~(60)--(62) in \cite{2018_PREFronczak}).

While presenting results in the Figures we used rather small sizes of the systems. However, the combinatorial framework allows to analyze systems of any size, including large ones. The problem is that obtaining relevant numerical simulation results for the size larger than $N = 400$ particles for an arbitrary kernel is computationally expensive due to the fact that we need to repeat the simulation millions of times to obtain reliable statistics. This is the main reason why we did not compare combinatorial and numerical results for larger sizes. To investigate possible dependency on the size and to be sure that quality of the theoretical predictions does not decrease with increasing system size we performed additional simulations for several values of $N$. In Figure~(\ref{Figure_sizedep}), we present an absolute difference between numerical and theoretical results as a function of the system size. No alarming observations were made.

Computational complexity of the expressions including Bell polynomials also may be an issue for calculations. However, in some cases (e.g., electrorheological coagulation \cite{2021_PhysD_Lepek}), it was possible to reduce Bell polynomials to much simpler form. Then, one can easily analyze much larger system sizes. Additionally, an efficient formula to calculate partial Bell polynomials is provided in the Appendix. Nonetheless, we suppose computational issues to be diminishing in future due to the increasing computational power of machines.

Some of the kernel forms solved in the combinatorial approach so far, i.e., linear--chain, $ K = \left(i^{-1} + j^{-1} \right)^{\alpha} $, condensation, $ K = (A + i) ( A + j) $, and a linear combination of the constant and additive kernels, $ K = A + i + j $, can be reduced to a form corresponding to the constant kernel. For the latter two and $ A=0 $, we obtain the previously known kernels, respectively, product and additive. In this way, theoretical solutions were obtained for all "intermediate" processes between the constant kernel and the product, additive, and linear--chain kernels (as shown in the Figure~\ref{Figure_Summary}). However, the quality of the combinatorial solutions decreases with increasing $ \alpha $. In the case of a condensation kernel, accuracy of the combinatorial solutions decreases once the gel point is exceeded. For the combination of the constant and additive kernels, the combinatorial solutions remain ''exact'' for any $ A $ and for any stage of the aggregation. Obtaining explicit theoretical solutions for such a wide continuous space of kernel forms can be perceived as a novelty in a research on the coagulation processes.

Obviously, several questions still remain unanswered. Why do combinatorial expressions give exact results for some kernels while approximate for other? Is it possible to modify the equations to cover arbitrary initial conditions or the post--gel phase as it was done in the for the product kernel in the Marcus--Lushnikov approach \cite{2019_PRE_Fronczak}? Are there other kernels which can be solved in the combinatorial framework with arbitrary parameters and extend the view in Figure~\ref{Figure_Summary}? Answering any of these questions would be a serious step forward in the coagulation theory.

\acknowledgments
This work has been supported by the National Science Centre of Poland (Narodowe Centrum Nauki) under grant no.~2015/18/E/ST2/00560 (A.F. and M.\L.).

\section*{Appendix}

As incomplete Bell polynomial definition is given by a Diophantine equation (which is NP--hard problem) other expression must be used to calculate its value in practice. An efficient way to calculate these polynomials is recurrence relation,
\begin{multline} \label{Bell_recursion}
 B_{n,k} \left( a_1, a_2, \dots, a_m, \dots, a_{n-k+1} \right) \\ = \sum_{m=1}^{n-k+1} \binom{n-1}{m-1} a_{m} B_{n-m,k-1},
\end{multline}

\noindent where $B_{0,0} = 1$, $B_{n,0} = 0$ for $n\ge1$, and $B_{0,k} = 0$ for $k\ge1$.

Some computational environments implement Bell polynomials (e.g., Wolfram Mathematica). For some others, a relevant code can be easily found in the Internet (e.g., Matlab). A problem in applying theoretical predictions is the explosion of digits of Bell polynomials for $N>70$ as this number immediately exceeds the precision available in standard programming environments. Therefore, calculations need to be performed with the help of arbitrary precision computation packages. An implementation for C++ environment using an arbitrary precision package can be found in \cite{cpp_libraries}.