{'fraction_non_alphanumeric': 0.037465362212294295, 'fraction_numerical': 0.009161341401345925, 'mean_word_length': 4.619161105815062, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
From synaptic interactions to collective dynamics in random neuronal networks models: critical role of eigenvectors and transient behavior

The study of neuronal interactions is currently at the center of several neuroscience big collaborative projects (including the Human Connectome, the Blue Brain, the Brainome, etc.) which attempt to obtain a detailed map of the entire brain matrix. Under certain constraints, mathematical theory can advance predictions of the expected neural dynamics based solely on the statistical properties of such synaptic interaction matrix. This work explores the application of free random variables (FRV) to the study of large synaptic interaction matrices. Besides recovering in a straightforward way known results on eigenspectra of neural networks, we extend them to heavy-tailed distributions of interactions. More importantly, we derive analytically the behavior of eigenvector overlaps, which determine stability of the spectra. We observe that upon imposing the neuronal excitation/inhibition balance, although the eigenvalues remain unchanged, their stability dramatically decreases due to strong non-orthogonality of associated eigenvectors. It leads us to the conclusion that the understanding of the temporal evolution of asymmetric neural networks requires considering the entangled dynamics of both eigenvectors and eigenvalues, which might bear consequences for learning and memory processes in these models. Considering the success of FRV analysis in a wide variety of branches disciplines, we hope that the results presented here foster additional application of these ideas in the area of brain sciences.


I. INTRODUCTION
Contemporary neuroscience is obsessed nowadays with detailed studies of the neuronal connections across the entire human brain. Large scale collaborative efforts [1] including the BRAIN Initiative in the USA, Brainome in China, and the BlueBrain in the European Union were launched with the objective of mapping at different resolution the connectivity of the entire brain. At a certain point theory will be desperately needed to analyze these very large maps, describing the adjacency matrix of the brain. The present work is an attempt to enter in this uncharted and challenging territory.
Under certain constraints, mathematical theory can advance predictions of the expected neural dynamics based solely on the statistical properties of their synaptic interaction matrix.
In that sense randomly connected networks of neurons are one of the classical tools of theoretical neuroscience. Only recently it was observed that the fact that the matrix of the synaptic connections is not normal (i.e., the matrix does not commute with its transpose), has dramatic consequences for the temporal dynamics of stochastic equations, which can mimic the dynamics of the network [2][3][4]. In particular, the work of Marti et. al [5] shows that increasing the symmetry of the connectivity leads to a systematic slowing-down of the dynamics and vice versa, decreasing the symmetry of the matrix leads to the speeding of the dynamics. This non-normality of the matrix not only forces matrices to have complex spectrum (which challenges several traditional tools of random matrix theory), but more importantly, its study sheds new light on the role of Bell-Steinberger [6] matrix of overlaps between the left and right eigenvectors of the connectivity matrix.
Contemporarily, the pivotal role of overlaps is understood in the simplest case of the evolution of complex Ginibre matrix -either in Smoluchowski-Fokker-Planck formalism [7,8] or in Langevin formalism [9], following the pioneering paper [10,11]. The effects of the overlaps of the Ginibre matrix for the temporal autocorrelation function of randomly connected networks was addressed analytically in the latest paper [5], confirming the numerical simulations in the weakly coupled regime of synaptic models.
In this paper we study the non-normality aspects of the popular model with excitatoryinhibitory structure [12][13][14], proposed by Rajan and Abbott in [15]. An important ingredient of this model is the introduction of the balanced condition, which stabilize the wildly fluctuating spectra of the network. In a later paper [16], the numerical study of the full non-linear dynamics in the Rajan-Abbott model has shown the emergence of a transition leading to synchronized (stationary or periodic) states. This phenomenon cannot be explained solely by the spectral features of the connectivity matrix, which motivates our study of missing non-spectral properties of non-normal networks, such as sensitivity to perturbations and transient dynamics induced by non-orthogonality of eigenvectors. Recently, it was hypothesized that the non-normality is universal in real complex networks [17].
Free random variables (hereafter FRV) theory is a relatively young mathematical theory, originating from the works of Voiculescu [18]. Partly due to the connection with large random matrices, it made in last decade a huge impact on physics [19], statistical inference [20], engineering of ICT technologies [21] and finances [22]. In brief, FRV can be viewed as a non-commutative probability theory for Big Data problems, where the information is hidden in statistical properties of eigenvalues and eigenvectors. As such, it is ideally suited for disentangling signals from noise in various kinds of complex systems. Another advantage comes form the fact, that at the operational level the formalism is simple and powerful, allowing very often to get results on the basis of "back-of-envelope" calculations.
From this perspective, it is rather bewildering, that FRV so far has not been broadly applied to the most challenging complex problem, i.e., the understanding of the brain.
Thus, in this paper we consider FRV applications to understand the neuronal networks as represented by the synaptic strength matrix. Direct application of FRV not only allows us to recover in a straightforward way well known results in the literature [15], but also to address quantitatively such issues as stability of the network with respect to perturbation and extension the existing formalisms for the heavy-tailed distributions.
The paper is organized as follows. In Sec. II we discuss two important effects caused by the non-orthogonality of eigenvectors of non-normal matrices, namely high sensitivity of the spectra and the transient behavior of the linearized dynamics. We briefly describe free probability theory in Section III, showing how it allows one to calculate the spectral density and gives an access to the eigenvector non-orthogonality. In Section IV we reframe the model introduced by Rajan and Abbott in this language. Applying the theoretical toolbox explained in Appendices we recover and generalize their main results for the unbalanced network. In doing so, we uncover the analytic formulas for the one-point eigenvector correlation function for this model, crucial for the determining its stability. Since free random variables work also in the case of heavy, spectral tails [23], we also present results for the spectra and eigenvectors of the Rajan-Abbott model adapted for the case of the Cauchy noise. We successfully confirm our analytic predictions with numerical simulations.
Further, in Section V, we show explicitly that the excitation/inhibition balance condition, not only tames the spectral outliers, but also exerts dramatic effects on non-orthogonality of eigenvectors, increasing the networks' eigenvalue condition number by several orders of magnitude. Section VI closes the paper with a summary of the main results and its implications. It also outlines main promising directions for further studies using the presented formalism. This influences not only their spectra, as the eigenvalues can be complex, but also has a strong effect on the eigenvectors. A diagonalizable non-normal matrix possesses two eigenvectors: left and right for each eigenvalue. They satisfy eigenproblems We use here physicists' "bra-ket notation", where |R i is a column and L i | is a row vector.
The scalar product is denoted as L i |R j and we define conjugated left vector |L i = ( L i |) † .
Eigenvectors are normalized to L i |R j = δ ij , but they are not orthogonal among themselves R i |R j = δ ij = L i |L j . Chalker and Mehlig introduced a matrix of scalar products of eigenvectors [10,11] Below we describe two phenomena important in neural networks, in which the nonorthogonality of eigenvectors captured in the matrix of overlaps plays an essential role.
A. Synaptic plasticity seen as perturbations of a network The synaptic strengths of real neuronal networks are not static [24]. Neural activity itself, in the course of time, allows neurons to form new connections, strengthening or weakening the existing synapses. This synaptic plasticity, on which biological learning is based, is not captured in many models. Nonetheless, the change of the synaptic strengths in a short time interval can be treated as a small additive perturbation of the initial matrix. This results in reorganization of the spectrum on a complex plane.
Considering the perturbation of the matrix X by some P , the change of the spectrum in the first order in reads The inequality follows from the Cauchy inequality and ||P || F denotes the Frobenius norm ||P || 2 F = TrP P † . This inequality is saturated (equality holds) by the rank one Wilkinson matrix P = |L i R i |. The inequality above shows that spectra of networks represented by non-normal matrices are more sensitive to changes in their connectivity. This enhanced sensitivity is driven by the non-orthogonality of eigenvectors. The quantity κ(λ i ) = √ O ii is known in the numerical analysis community as the eigenvalue condition number [25,26].

B. Eigenvector non-orthogonality in transient dynamics
Stability analysis and linear response of the dynamic systems with respect to external perturbations are among most popular methods of describing complex systems [36]. Let us consider dynamics obtained from linearization of the system in the vicinity of the fixed point Here ξ represents the external driving. Choosing it as a "spike" |ξ(t) = δ(t) |ψ(0) , we formally solve the system for t > 0 The long-time dynamics is governed by the eigenvalue with the largest real part. However, if X is non-normal, this analysis is incomplete. The behavior of the linearized dynamics can be drastically different at its early stage. In particular, the system may initially move away from the fixed point. This sometimes invalidates the linear approximation and renders the fixed point unstable, even though the linearized dynamics predicts stability.
To describe this transient dynamics, we consider the squared Euclidean distance from the fixed point, which is the squared norm of the solution (5) If we consider the "spike" |ψ(0) as a particular versor on the N -hypersphere (real or complex), averaging over all directions uniformly distributed on the hypersphere leads tō We see that all elements of the overlaps of left and right eigenvectors drive the behavior of the squared distance. First, they enhance the contributions of the eigenmodes, which is responsible for amplification of the response to the external driving. Second, since the matrix is not diagonal, they couple different eigenmodes during the evolution. This results in an interference between eigenmodes, which is reflected as an oscillatory behavior of the squared norm of the solution (see also Fig. 7). Note that for normal matrices, such effects do not exist, since left and right vectors are orthogonal and the "coupling matrix" is an identity. Recently, the transient growth was proposed as an amplification mechanism of neural signals [2,27,28].
Usually the matrix X is modeled as random. We remark that the averaging over all initial conditions is equivalent to fixing an initial vector |ψ(0) and averaging over the vectors U |ψ(0) , where U is uniformly distributed (according to the Haar measure) on the orthogonal (unitary) group. This implies that if a randomness in X is to be compatible with the average over initial conditions, the probability density function has to be invariant under unitary (orthogonal) transformations, P (X) = P (U XU † ).

A. Spectral density and eigenvector correlations
Unitarily (and orthogonally) invariant random matrices in the large size limit are described by Free Probability. Its power relies on the easiness of obtaining analytical formulas, which are very good approximations even for a relatively small size of a matrix.
An important class of matrices, the so-called bi-unitarily invariant, which generalizes the Gaussian distribution (described in Sec. III B) is important in the models of neural networks.
Despite the fact that they are genuinely non-Hermitian, due to enhanced symmetry the spectral problem is effectively one-dimensional because the spectrum is rotationally invariant on a complex plane. In this case, a powerful result holds in FRV, known as the Haagerup-Larsen theorem [29]. It states that the radial cumulative distribution function F (r) = r 0 2πρ(r )r dr , of the ensemble X can be inferred from the simple functional equation where S X (z) is so called S-transform for the ensemble X. In Appendix A we explain the probabilistic interpretation of S and we provide a simple example. Spectra of bi-unitarily invariant ensembles in large N limit are supported on either a disc or an annulus, a phenomenon dubbed "the single ring theorem" [30,31]. The inner radius of the spectrum is deduced from the condition F (r in ) = 0, while the outer one is given by F (r out ) = 1.
The applicability of free probability to non-Hermitian matrices is not limited to spectra only. It gives also access to the local averages of the overlap matrix. The one-point function associated with the diagonal elements of the overlap matrix can be calculated for any type of unitarily invariant probability [33]. For bi-unitarily invariant ensembles it takes remarkably simple form [34] O(r) = 1 The ratio of the one-point correlation function and the spectral density gives the conditional expectation of the squared eigenvalue condition number [34] Recently, the two-point function associated with off-diagonal elements of the overlap matrix has become accessible within free probability [35].

B. Example: Ginibre-Girko ensemble
We conclude this section with an example of the above construction by considering the so called Ginibre-Girko matrix X, the entries of which are independently taken from the real/complex Gaussian distribution with zero mean and 1/N variance. According to (8), we need the S-transform for X † X. This matrix belongs to the Wishart ensemble [37]. Its . This completes the calculation, since now replacing z → F (r) − 1 and using (8) we get The spectrum is therefore uniform, ρ(r) = 1 2πr dF (r) dr = 1 π , on the unit disc (F (r in ) = 0, F (r out ) = 1), reproducing Ginibre-Girko result. The eigenvector correlator comes from (10), , in agreement with [10], where it was calculated using much more laborious techniques. In the next section, we show that the same computational simplicity is preserved when considering the ensembles taking into account physiological restrictions imposed on the neural networks models.

IV. REFRAMING RAJAN-ABBOTT MODEL
The strength of synapses between all pairs of N neurons in a network is represented by the weighted adjacency (synaptic) matrix. Contrary to the Ginibre matrices, the structure of its elements is more complicated. In the minimal model [15], there are two kinds of neurons with a fraction f E N representing excitatory (E), and f I N = (1 − f E )N the remaining inhibitory (I) neurons. Their strengths are sampled from Gaussian ensembles, with means µ i and The matricial representation of the synaptic strength matrix reads therefore M + W . Here the deterministic matrix M represents the average synaptic activity. In this model it is a rank one matrix with identical rows, each containing f E N consecutive means µ E and followed by f I N consecutive means µ I . The random part W models variability across the population. It is assumed to be of the form W = XΛ, where X is the Girko-Ginibre matrix and Λ is diagonal with its first f E N elements equal to σ E and last f I N ones equal to σ I .
Several studies [14,38] show that the amount of excitation and inhibition of a neuron is the same (the so-called E/I balance) even on the scale of few milliseconds [12,13]. To incorporate this fact in the model, the balance condition is imposed on two levels. The global condition f E µ E + f I µ I = 0 means that neurons are balanced on average. This forces the last non-zero eigenvalue of M to vanish. Even in the case of a null spectrum of M , its non-normal character causes the eigenvalues of M + XΛ differ much from that of XΛ. As a result a few eigenvalues lie far beyond the spectrum of XΛ [15], see Fig. 1.
The local E/I balance is imposed on this model by demanding that the sum of strengths coupled independently to each neuron vanishes. Mathematically, the elements within each row sum to zero. This condition brings the outliers back to the disc of radius R = Explicit solutions exist for m = 2, 3, 4 types of neurons, corresponding to the quadratic, cubic or quartic algebraic equation for F (r), but other cases are easily tractable numerically. The case solved by Rajan and Abbott corresponds to the quadratic equation. Solution (13) is also equivalent to the diagrammatic construction of [39], but more explicit. The spectrum is always confined within the disc of radius r 2 out = m i=1 f i σ 2 i , as visible from the condition F (r out ) = 1.
We will argue in Sec. V that the presence of the deterministic matrix M and the balance condition exert a dramatic effect on the eigenvectors of the synaptic strength matrix.
Knowing F (r), free probability allows us to calculate via (10) also the eigenvector correlation function O(r) for its random part W . In the case of the minimal model considered by

Rajan-Abbott, it reads explicitly
where This result is inaccessible within the framework of [39].

B. Heavy-tailed noise
Cauchy noise, belonging to the regime of Lévy stable distributions, is used here as the simplest mechanism to mimic the non-Gaussianity of the realistic synaptic matrices. Since learning rules could change the initial random network structure into a small world network [40][41][42] by dynamic modification of synaptic weights, possibility of obtaining analytic benchmarks for heavy-tailed distributions is appealing. Spatial and temporal Lévy processes are omnipresent in biological time series, but the fact that they do not possess finite moments invalidates several standard tools of statistical analysis. In the case of matrices exhibiting heavy-tailed distributions of elements, the underlying mathematical structure is quite involved [43,44]. Surprisingly, also in this case FRV offers a powerful shortcut.
Application of FRV techniques for the heavy-tailed noise with α = 1 (spectral Cauchy distribution) leads (see Appendix C) to the simple result In this case, the spectrum spreads over the whole complex plane, reflecting the large fluctuation of Lévy type noises. In the case of more realistic Lévy noises, one looses the simple analytic structure presented above, but the formalism stays -the resulting equations are usually of transcendental type, but can be easily solved numerically.

V. NON-NORMALITY IN RAJAN-ABBOTT MODEL
Below we argue that imposing E/I balance not only confines the eigenvalues to a disc, but -more importantly -induces very strong non-orthogonality of eigenvectors. This in turn causes the spectra to be highly sensitive to perturbations and strengthen the transient effects.
Let us assume that the matrix W is diagonalizable. If we denote |u = (1, 1, . . . , 1) T , the E/I balance is equivalent to the fact that |u is the right eigenvector of W to the eigenvalue Since m|u = 0, m| has a decomposition into the left eigenvectors of W , except for L 1 |, Hence, the total synaptic strength matrix is decomposed as We constructed explicitly the eigenvectors of the synaptic strength matrix. The left eigenvectors are not altered when M is taken into consideration due to the E/I balance. The bi-orthogonality condition L i |R j = δ ij leaves freedom of rescaling each pair of eigenvectors by a non-zero complex number |R j → c j |R j and L j | → L j | c −1 j . These transformations allow us to set the length of left eigenvectors L j |L j = 1. The diagonal elements of the new where we have used u|u = N . In the large N limit the last term will dominate. This shows that the condition numbers grow with the size of a matrix and the effect of the matrix of averages is stronger for eigenvalues close to the origin.
Analogous reasoning for the full overlap matrix leads to the conclusion that all its elements O ij for i, j ≥ 2 are affected by the E/I balance and the deterministic matrix. The dominant term in large N is given by To study the statistics of the eigenvalue condition numbers, we performed numerical simulations by diagonalizing matrices, the random part was generated from either real or complex Ginibre ensemble. The eigenvector correlation function is juxtaposed with (14) from There is a visible mismatch between numerics for real matrices and the results from free probability, particularly evident for eigenvalues with small modulus. This fact is explained in the light of the recent result by Fyodorov [45], where he showed that the distribution of the overlap for Gaussian matrices is heavy-tailed. This distribution conditioned on real eigenvalues of the real Ginibre is so fat-tailed that even the mean does not exist, thus O(z) can be considered only outside the real axis. Being aware of this fact, we have performed further simulations only for complex matrices which do not suffer from this problem.
We studied the effect of the deterministic matrix M by juxtaposing the eigenvector correlation function in Fig. 4 and noticed the significant increase in its magnitude. This enhancement of non-normality is visible not only on the level of mean value, but also on the full distribution of the overlap (see Fig. 4 (right)).
Above conclusions are strengthened by the similar study based on Cauchy synaptic matrices. Figure 2 shows perfect agreement of our predictions with the numerics. By construction of the Rajan-Abbott local condition, spectra are unchanged. This does not hold, however, for the squared eigenvalue condition numbers -they dramatically increase (several orders of magnitude, note the scales in Figures 5 and 6). Finally, the unperturbed eigenvector correlator approaches the predicted slope (compare the predicted slope 4 to the measured 3.84). This dynamics is obtained by the linearization of the model considered in [16]. Results presented in Fig. 7 show that the deterministic connections in the network followed by the E/I balance significantly enhance the norm of the solution and the transient trajectories are present for almost all initial conditions. This is not the case if the connections were fully random. Moreover, the wild oscillations of the squared norm indicate the strong interference between the eigenmodes.
One expects these dramatic effects to be visible in the activity of individual neurons. We therefore studied the temporal dynamics of the components of the vector of neural activities The matrix is of size N = 100 with the same parameters as in Fig. 1. We chose µ = r out + 0.02 to ensure stability. Each curve corresponds to a single initial condition generated randomly from the set of vectors of unit norm.
(5) for randomly chosen initial conditions. The results, presented in Fig. 8, show that in the presence of M , the neuronal activity is not only transiently enhanced, but also more synchronized, as observed numerically in the full dynamics in [16]. This effect persistent in the non-linear model is observed as transient in the linearized dynamics. Parameters are the same as in Fig. 7.

VI. DISCUSSION AND CONCLUSIONS
In this letter we explored the use of FRV in the study of large synaptic interaction matrices. Beside recovering in a straightforward way known results on the application of random matrices to neural networks, we have addressed the issue of large fluctuations, most probably very relevant to the dynamics of learning and memory in biological neural networks [46]. Using recent results on the properties of eigenvectors in non-normal matrices, we have quantitatively linked the strength of the fluctuation of the outliers to a certain eigenvector correlator. We presented our analysis for the simplest Gaussian case, nevertheless we also pointed the way how one can consider other distributions, e.g., heavy tails. The formalism stays the same, perhaps the only difference (modulo easy Cauchy case which we solved here analytically) is that in the case of more general pdf's one has to rely on numerical solutions.
From the present results it is clear that to understand the temporal evolution of nonnormal matrix models requires considering the entangled dynamics of both eigenvectors and eigenvalues, contrary to the simple evolution of the spectra of normal matrices, for which the eigenvectors decouple in the presence of the spectral evolution [7,8,47] Our results indicate, that for the balanced networks the sensitivity of eigenvalues to additive perturbation is dramatic and increases several orders of magnitude in the networks with heavy-tailed spectrum of adjacency matrices (small worlds). Since it is commonly accepted, that spike-timing-dependent plasticity in small-world networks is a hypothetical learning mechanism (for a recent experimental study see [48]), one may worry, how synchronization of the network is possible at all. We emphasize here that the E/I balance is put into this model by hand. In the real brain the E/I balance is maintained on the scale of hundreds of milliseconds [14], and periods during which the balance is violated are not longer than few milliseconds [12,13]. More complete models of neural networks must incorporate the E/I balance as a dynamical process.
Networks adapting to the changing external conditions may change their structure in a controlled way. The high sensitivity of eigenvalues to these changes in this case might be desired, because it can facilitate the adaptation. We hypothesize that such high sensitivity in the models with dynamical E/I balance can emerge through a process a kind of selfregulated criticality [49]. Although the specifics of such process is not certain yet, there is evidence both empirical [50][51][52][53] as well as theoretical [54][55][56] of its plausibility. In addition, the connection of E/I balance with criticality was already observed at the level of neuronal avalanches analysis in EEG or MEG data [57].
Since the balance condition leads to a dramatic increase of eigenvector overlaps -conditioning the spectra -which further take crucial part in driving mechanisms of temporal evolution of the networks, one needs a powerful, stabilizing mechanism preventing the transition to the chaotic behavior. Such chaotic behavior would imply sudden and drastic reorganization of the eigenvalues leading to unwanted dynamics of the neural network.
We envision one a priori mechanism, which can tame such a behavior -it is the transient behavior. This conclusion is consistent with the model of Molino et al. [16] for non-normal Let us consider the following problem. We have two random variables x 1 and x 2 drawn from independent probability distributions p 1 (x 1 ) and p 2 (x 2 ). The distribution of the random variable s being the sum of x 1 and x 2 reads therefore One can easily unravel the convolution using the Fourier transform (characteristic function).
Thenp(k) ≡ p(s)e iks ds =p 1 (k)p 2 (k), wherep i (k) are Fourier transforms corresponding to the original densities p i (x). Note that a characteristic function generates moments of the respective distribution. We can further simplify the problem if instead of characteristic functions we consider their natural logarithms φ(k) ≡ lnp(k). Then we get the addition law, which linearizes the convolution Since φ is another generating function -this time for cumulants of the distribution -above relation means simple addition of the corresponding cumulants. The algorithm of convolution is therefore simple. First, knowing p i (x), we construct the φ i (k). Then we perform the addition law (23). Finally, we reconstruct p 1+2 (s) from the φ 1+2 (k), performing the first step in reversed order. A pedagogical and straightforward example is represented by the convolution of two independent Gaussian distributions N 1 (0, σ 2 1 ) and N 2 (0, σ 2 2 ). First step shows, that in both cases only one cumulant is non-vanishing, i.e., the second one, κ (1) The addition law and the last step of the logarithm immediately lead to the result, that the resulting distribution is also Gaussian, N 1+2 (0, σ 2 = σ 2 1 + σ 2 2 ). In free probability, the notion of independence is replaced by the notion of freeness. Two large (infinite) matrices are mutually free if their eigenvectors are maximally decorrelated, e.g., matrices X and U Y U † , where U is the Haar measure, are free.
The role of the characteristic function is played by the complex valued Green's function where ρ X (λ) is the average spectral density of the matrix X, playing here the role of probability density function in CTP. Indeed, expanding G X (z) around z = ∞ we get spectral moments M (X) k = λ k ρ X (λ)dλ. Note, that knowing G X (z) we can easily reconstruct ρ X (λ). Indeed, The role of the generating function for free cumulants is played by the so-called R-transform, The crucial relation between R(z) and G(z) reads R(G(z)) + 1 G(z) = z or G(R(z) + 1/z)) = z, i.e., R(z) is -modulo the shift 1/z -the functional inverse of the Green's function. Let us come back to the problem of spectral addition. Imagine we have now the spectral measures ρ X i (λ), corresponding to two matricial ensembles with the measures P (X i )dX i , where i = 1, 2. We are now asking, what is the spectral density of the ensemble X 1+2 = X 1 + X 2 . This is a highly non-trivial and non-linear problem, since X 1 and X 2 do not commute, but free calculus allows to solve this case in full analogy to CTP.
The algorithm is as follows: First, from ρ i corresponding to X i we construct matching G i (z) and R i (z). Then R X 1 +X 2 = R 1 (z) + R 2 (z) (26) which supersedes (23). Finally, we proceed in reverse order, reconstructing from R X 1 +X 2 (z) the Green's function G X 1 +X 2 (z), and finally the spectral density ρ X 1 +X 2 (λ). As an example, we consider the "Gaussian" in free theory, i.e., the spectral distribution the only non-vanishing cumulant of which is variance σ 2 . Thus R(z) = σ 2 z. Reconstructing Green's function gives . Taking imaginary part we reconstruct the celebrated Wigner semicircle ρ(λ) = 1 We see that the addition algorithm for two free Wigner semicircles mimics precisely the addition algorithm of two Gaussians.
Similarly to addition, one can consider multiplication laws for random variables x 1 · x 2 .
In CTP, such problem is unravelled with the help of Mellin transform, see e.g., [59]. In free calculus, the role of the Mellin transform is played by S-transform, related to R-transform as S X (z)R X (zS X (z)) = 1. The multiplication law reads S X 1 X 2 (z) = S X 1 (z) · S X 2 (z) (27) and the algorithm for multiplication follows the one for addition. However, one should be aware that the product of two symmetric (hermitian) matrices may be non-symmetric (non-hermitian). In such a case, the eigenvalues can appear on the whole complex plane, and the methods of R(z) and S(z) transforms, based on analyticity, require substantial modifications. Luckily, there exist one powerful case, governed by the Haagerup-Larsen theorem (known also as a "single ring" theorem), when analytic methods hold for complex spectra. If the complex matrix X can be decomposed as X = P U , where P is positive, U is Haar-measured and P and U are mutually free, the spectrum on the complex plane has a polar symmetry and the radial distribution can be easily read out from the singular values of X, i.e., the real eigenvalues of X † X. In mathematics, such ensembles are known as R-diagonal. To infer the information about the spectra and some correlations between left and right eigenvectors one needs only the explicit form of S X † X (z). In the case of the Ginibre ensemble X (i.e., where X ij are drawn either from real or complex Gaussian distributions), this is particularly easy, since matrix X † X is known as a Wishart ensemble.
To avoid obscure mathematics, let us recall that the Wishart ensemble is a free analogue of the Poisson distribution from classical probability [18]. This implies that all cumulants are the same, and if for convenience normalized to 1, its R transform is just, by definition, Using the above-mentioned functional relation between R and S transform we arrive at S X † X = 1 1+z . Similar techniques can be applied for generic randomness in Rajan-Abbott type models, as we show below.

B Rajan-Abbott model with Gaussian noise
We use the theorem from free probability, which states that the product of R-diagonal operator with any operator is R-diagonal [60], therefore W is subject to the Haagerup-Larsen theorem. Then W † W = ΛX † XΛ tr = X † XΛ 2 , where the last equation expresses the fact that the spectral properties are invariant under the cyclic permutations of matrices under the trace. The Green's function (resolvent) for Λ 2 reads therefore Substituting z → R Λ 2 (z)+ 1 z in (28) and using the fundamental FRV relation G(R(z)+ 1 z ) = z we arrive at Now we replace in (29) z → tS Λ 2 (t) and using the relation between S and R transforms we arrive at We note that 1 1+t is the S-transform for the Wishart ensemble (calculated above), and the multiplication law gives us the final S-transform for W † W , i.e., 1 1+t S Λ 2 (t) = S X † X (t)S Λ 2 (t) = S W † W (t), so we arrive at In the last step we substitute t → F (r) − 1 and use the Haagerup-Larsen theorem, arriving at Subtracting 1 = i f i from both sides finally gives the solution .

(33)
C Rajan-Abbott model with Cauchy noise FRV calculus is a powerful technique and the range of its applications is not confined to the basin of attraction of the Gaussian type. In particular, for random matrices X belonging to the free Lévy class (spectrum behaving like 1/λ α−1 ), the S-transform for the Wishart-Lévy matrix X † X reads S X † X (t) = 1 t(1+t) t b t/α , with b = exp[iπ(α/2 − 1)] [61]. Stability index α = 2 reproduces the Gaussian case, but a simple form can be obtained also for the Cauchy disorder α = 1. In this case S X † X (t) = − t 1+t , and when applied to (30), yields Contrary to the previous case, the spectrum is unbounded and stretches up to infinity.
Explicitly, the spectral density and the eigenvector correlator read ρ(r) = 1 2πr In the case of arbitrary α, resulting transcendental equations can be easily solved numerically.
Other type of randomness of the neural networks can also be modeled, e.g., by considering Student-Fisher spectral distributions.