### Derivation of the size of the giant component

Let us consider a system of nodes bearing 2*N* different types of half-edges. The probability that a node is characterised by the half-edge count vector ({varvec{k}}=(k_1,dots ,k_{2N})) is given by the multivariate degree distribution (u({varvec{k}})). The permutation matrix *P* defines a pairing on the 2*N* half-edge types, with (P_{i,j}=1) if half-edge types *i* and *j* form a bond and (P_{i,j}=0) otherwise.

Let the *i*-excess degree distribution (u_i({varvec{k}})) be the probability that a node connected to a randomly sampled half-edge of type *i* has excess configuration ({varvec{k}}) (thus not counting the already sampled half-edge). It is given by

$$begin{aligned} u_i({varvec{k}})=(k_i+1)frac{u({varvec{k}}+varvec{e_i})}{{mathbb {E}}[k_i]}, end{aligned}$$

(12)

where (varvec{e_i}) are the standard basis vectors.

Let ({varvec{z}}=(z_1,dots ,z_{2N}),;|z_i|<1,; z_iin mathbb C), we define the generating functions of the degree distribution

$$begin{aligned} U({varvec{z}}): {mathbb {C}}^{2N} rightarrow {mathbb {C}} end{aligned}$$

and the *i*-excess degree distribution

$$begin{aligned} U_i({varvec{z}}): {mathbb {C}}^{2N} rightarrow {mathbb {C}} end{aligned}$$

as

$$begin{aligned} U({varvec{z}})&=sum _{{varvec{k}}>0} {varvec{z}}^{{varvec{k}}} u({varvec{k}}),\ U_i({varvec{z}})&=sum _{{varvec{k}}>0} {varvec{z}}^{{varvec{k}}} u_i({varvec{k}})=frac{1}{{mathbb {E}}[k_i]}frac{partial U({varvec{z}})}{partial z_i}, ;i=1,dots ,2N end{aligned}$$

(13)

Note that the vector power in these definitions and in what follows is evaluated as

$$begin{aligned} {varvec{z}}^{{varvec{k}}}:=prod _{i=0}^{2N} z_i^{k_i}. end{aligned}$$

Let (W(x): {mathbb {C}} rightarrow {mathbb {C}}) and (W_i(x): {mathbb {C}} rightarrow {mathbb {C}}) for (i=1,dots 2N) be a family of complex functions satisfying the system of (2N+1) equations:

$$begin{aligned}{}& W(x)=xU({varvec{P}}varvec{omega }(x)),\&W_i(x)=xU_i({varvec{P}}varvec{omega }(x)),,i=1,dots ,N end{aligned}$$

(14)

with (varvec{omega }(x):=left( W_1(x),dots ,W_{2N}(x)right) ^top ). Here (W_i(x)) play an auxiliary role as we are mainly concerned with the properties of *W*(*x*). The theory for coloured random graphs^{34} interprets *W*(*x*) as the generating function for the probability *w*(*n*) that a uniformly at random chosen node is part of a finite component of size *n*, that is (W(x)=sum _{n>0} x^{n} w(n)) with (xin {mathbb {C}},,|x|le 1.) Therefore, the probability that a randomly sampled node is not in any finite-sized component but in the giant component is given by

$$begin{aligned} g_text {node}=1-W(1), end{aligned}$$

(15)

where *W*(1) can be calculated by substituting (x=1) in Eq. (14), yielding

$$begin{aligned} W(1)&=sum _{{varvec{k}}>0} u({varvec{k}}) (varvec{Ps})^{{varvec{k}}}={mathbb {E}}[(varvec{Ps})^{{varvec{k}}}],; end{aligned}$$

(16)

with ({varvec{s}}=varvec{omega }(1)), satisfying system of equations

$$begin{aligned} s_i&=W_i(1)=U_i(varvec{Ps})\&=sum _{{varvec{k}}>0} u_i({varvec{k}}) (varvec{Ps})^{{varvec{k}}}\&=frac{1}{{mathbb {E}}[k_i]} sum _{{varvec{k}}>0} (k_i+1) u(k_i+varvec{e_i}) (varvec{Ps})^{{varvec{k}}}\&=frac{{mathbb {E}}[k_i(varvec{Ps})^{{varvec{k}}-varvec{e_i}}]}{{mathbb {E}}[k_i]}, end{aligned}$$

(17)

for (i=1,dots ,2N). These results were stated in Eqs. (7) and (8).

### Derivation of the weight-average size of finite components

We derive the expression of the weight-average size of finite components, defined as

$$begin{aligned} langle s rangle _w=frac{frac{d}{dx}W(x)|_{x=1}}{1-g_text {node}}, end{aligned}$$

(18)

and normalised to the system size without the giant component. The evaluation of (frac{d}{dx}W(x)) results in

$$begin{aligned} frac{dW(x)}{dx} bigg |_{x=1}&=U(varvec{Pomega }(x))|_{x=1}+\&quad+bigg (frac{partial U({varvec{z}})}{partial z_1} frac{d(varvec{Pomega }(x))_1}{dx} + dots \&quad+frac{partial U({varvec{z}})}{partial z_N} frac{d(varvec{Pomega }(x))_N}{dx} bigg )bigg |_{x=1,{varvec{z}}={varvec{s}}}\&=W(1)+left( frac{partial U({varvec{z}})}{partial z_1},dots ,frac{partial U({varvec{z}})}{partial z_N}right) bigg |_{{varvec{z}}={varvec{s}}}{varvec{P}}{varvec{y}}\&=1-g_text {node}+{varvec{s}}^top {varvec{D}}{varvec{P}}{varvec{y}}. end{aligned}$$

(19)

with ({varvec{y}}=frac{d}{dx}{varvec{w}}(x)|_{x=1}) and ({varvec{D}}=text {diag}{{mathbb {E}}[k_1],dots ,{mathbb {E}}[k_{2N}]}). The elements of vector ({varvec{y}}=(y_1,dots ,y_N)^top ), where (y_i=frac{partial }{partial x}W_i(x)|_{x=1}), are given by

$$begin{aligned} y_i&=U_i(varvec{Pomega }(x))|_{x=1}+\&quad+bigg (frac{partial U_i({varvec{z}})}{partial z_1} frac{d(varvec{Pomega }(x))_1}{dx} + dots \&quad+frac{partial U_i({varvec{z}})}{partial z_N} frac{d(varvec{Pomega }(x))_N}{dx} bigg )bigg |_{x=1,{varvec{z}}={varvec{s}}}\&=s_i+left( frac{partial U_i({varvec{z}})}{partial z_1},dots ,frac{partial U_i({varvec{z}})}{partial z_N}right) bigg |_{{varvec{z}}={varvec{s}}}{varvec{P}}{varvec{y}}. end{aligned}$$

(20)

Utilising the matrix function ({varvec{X}}({varvec{s}})) with its elements

$$begin{aligned} X_{i,j}({varvec{s}})&=frac{partial U_i({varvec{z}})}{partial {z_j}}bigg |_{{varvec{z}}={varvec{s}}}\&=frac{{mathbb {E}}[(k_ik_j-delta _{i,j}k_i){varvec{s}}^{{varvec{k}}-varvec{e_i}-varvec{e_j}}]}{{mathbb {E}}[k_i]}, end{aligned}$$

(21)

with (i,j=1,dots ,2N), Eq. (20) is written in its matrix form ({varvec{y}}={varvec{s}}+{varvec{X}}({varvec{s}})varvec{Py}). Thus, ({varvec{y}}) is expressed explicitly as

$$begin{aligned} {varvec{y}}&=[{varvec{I}}-varvec{XP}({varvec{s}})]^{-1}{varvec{s}}. end{aligned}$$

(22)

Having Eqs. (19) and (22) in hand, Eq. (18) is written as

$$begin{aligned} langle s rangle _w=1+ frac{{varvec{s}}^top {varvec{D}}{varvec{P}}[{varvec{I}}-varvec{XP}({varvec{s}})]^{-1}{varvec{s}}}{1-g_text {node}}, end{aligned}$$

(23)

which is stated in Eq. (9).

### Mathematical description: living polymerisation of polymer chains

In this section, we introduce a master equation for monomer species in living polymerisation. First we will discuss the formulation without colours, followed by the monomer PBE with *N* colours. Living polymerisation is a special type of chain-growth polymerisation characterised by a very fast initiation and suppressed termination of radicals. Hence, only propagation is modelled explicitly. We consider a monomer with one functional group (vinyl group), which allows the formation of at maximum two bonds per monomer and therefore leads to linear strands of polymers.

#### Monomer equation for living chain polymerisation

The monomer equation describes the concentration of the *monomer* species over time by a system of ODEs, rather than the concentrations of polymers of different size as is the case with the conventional PBE. We consider a system with the initial fraction of active monomers units to be set to (R_0=0.01). During the propagation reaction, the active monomer unit connects with an unreacted monomer and a covalent bond is formed. Due to the asymmetry in the reactants (vinyl+radical), the formed connection is considered asymmetric as well. This asymmetry is represented as the orientation of a directed edge. Monomer species (mathcal M_{k_1,k_2}) are characterised by the number and direction of edges they bear—the number of in-edges (k_1) and out-edges (k_2).

In order to describe living polymerisation, we need to distinguish between five monomer species: (1) free monomers ({mathcal {M}}_{0,0}), (2) active monomers ({mathcal {M}}^*_{0,0}) with zero (half-) edges, (3) active monomers ({mathcal {M}}_{1,0}) with one in-edge and zero out-edges, (4) dead ends ({mathcal {M}}_{0,1}), which initially were active sites and have already connected to one monomer with zero in-edges and one out-edge, and (5) consumed monomers (mathcal M_{1,1}) with one in-edge and one out-edge. The polymerisation process is described by the following reaction mechanism:

$$begin{aligned}&{mathcal {M}}^*_{0,0}+{mathcal {M}}_{0,0}xrightarrow {k_text {p}}{mathcal {M}}_{0,1}+{mathcal {M}}_{1,0},\&{mathcal {M}}_{1,0}+{mathcal {M}}_{0,0}xrightarrow {k_text {p}}mathcal M_{1,1}+{mathcal {M}}_{1,0}. end{aligned}$$

(24)

Active monomer units ({mathcal {M}}^*_{0,0}) or ({mathcal {M}}_{1,0}) connect to free monomers ({mathcal {M}}_{0,0}) at rate (k_text {p}). The product species consist of one consumed monomer (mathcal M_{0,1}) or ({mathcal {M}}_{1,1}) and a new active monomer (mathcal M_{1,0}). This reaction mechanism translates into the following system of ODEs for monomer concentrations, which we denote as ({mathcal {M}}_{i,j}(t)):

$$begin{aligned}&dot{{mathcal {M}}}_{0,0}(t)= -k_text {p} {mathcal {M}}_{0,0}(t) [{mathcal {M}}_{1,0}(t) + {mathcal {M}}^*_{0,0}(t) ], \&dot{{mathcal {M}}}^*_{0,0}(t)= -k_text {p} {mathcal {M}}_{0,0}(t) {mathcal {M}}^*_{0,0}(t), \&dot{{mathcal {M}}}_{0,1}(t)= k_text {p} {mathcal {M}}_{0,0}(t) {mathcal {M}}^*_{0,0}(t),\&dot{{mathcal {M}}}_{1,0}(t)= k_text {p} {mathcal {M}}_{0,0}(t) {mathcal {M}}^*_{0,0}(t),\&dot{{mathcal {M}}}_{1,1}(t)=k_text {p} {mathcal {M}}_{0,0}(t) mathcal M_{1,0}(t), end{aligned}$$

(25)

with the initial conditions being set as follows: (mathcal M_{0,0}(0)={mathcal {M}}_0(1-R_0)), ({mathcal {M}}^*_{0,0}(0)=mathcal M_0R_0), ({mathcal {M}}_{0,1}(0)=0), ({mathcal {M}}_{1,0}(0)=0) and ({mathcal {M}}_{1,1}(0)=0).

We quantify the progress of the polymerisation by the conversion of functional groups, (chi (t)), which, in the case of linear polymers, is equivalent to the conversion of monomers:

$$begin{aligned} chi (t)=1-frac{{mathcal {M}}(t)}{{mathcal {M}}(0)}. end{aligned}$$

(26)

The two-variate degree distribution (u(k_1,k_2)) at any time point *t* is extracted from the solution of Eq. (25) by setting:

$$begin{aligned} u(k_1,k_2)=frac{{mathcal {M}}_{k_1,k_2}(t)+mathcal M^*_{k_1,k_2}(t)}{sum _{k_1,k_2}[{mathcal {M}}_{k_1,k_2}(t) +mathcal M^*_{k_1,k_2}(t)]}. end{aligned}$$

(27)

#### ODEs with colours for living chain polymerisation

We extend the system of ODEs from the previous section to incorporate the formation time of edges in the species description, which allows to compute the coloured degree distribution. Consider *N* time intervals (Delta t_c=t_{c+1}-t_c) with (c=1,dots ,N). To each newly formed edge we assign colour *c* if its formation time (tin [t_{c-1},t_c)). The monomer species (mathcal M_{{varvec{k}}}) are now indexed by the count vector

$$begin{aligned} {varvec{k}}=(k_1,dots ,k_{2N}), end{aligned}$$

(28)

where the index (k_{2(i-1)+1}) denotes the count of in-edges of colour *i* and (k_{2i}) the count of out-edges of colour *i*, for (i=1,dots ,N). In the case of linear polymer chains, the count vector elements are restricted to (k_{2(i-1)+1},,k_{2i}in {0,1}), which does not hold true for a polymer network. We group the count vectors by the total number of in-edges and out-edges irrespectively of their colour:

$$begin{aligned} {varvec{k}}_{m,n}={ {varvec{k}} : sum _{i=1}^N k_{2(i-1)+1} =m, ; sum _{i=1}^N k_{2i}=n }. end{aligned}$$

Let the elements of vectors (Delta {varvec{k}}_{1,0}(c)) and (Delta {varvec{k}}_{0,1}(c)) be given by

$$begin{aligned}&[Delta {varvec{k}}_{1,0}(c)]_i=delta _{i,2(c-1)+1},\&[Delta {varvec{k}}_{0,1}(c)]_i=delta _{i,2c},\ end{aligned}$$

(29)

with (delta _{i,j}) being Kronecker delta, and let

$$begin{aligned} Theta (t,c)= [theta (t-t_c)+theta (t_{c+1}-t)], end{aligned}$$

with (theta (t)) the Heaviside step function. For all (c=1,dots ,N) representing time intervals, the reaction mechanism of Eq. (24) is extended using the multi-index notation:

$$begin{aligned}&{mathcal {M}}^*_mathbf{0}+{mathcal {M}}_{ mathbf{0}}xrightarrow {k_text {p}} {mathcal {M}}_{ Delta {varvec{k}}_{0,1}(c)} + {mathcal {M}}_{Delta {varvec{k}}_{1,0}(c)},\&{mathcal {M}}_{{varvec{k}}}+{mathcal {M}}_{ mathbf{0}}xrightarrow {k_text {p}} {mathcal {M}}_{{varvec{k}+Delta {varvec{k}}_{0,1}(c)}} + {mathcal {M}}_{ Delta {varvec{k}}_{1,0}(c)},;{{varvec{k}}}in {varvec{k}}_{1,0} end{aligned}$$

(30)

The latter reaction mechanism translates into a system of ODEs as follows. For monomer species that have no half-edges, we have:

$$begin{aligned} dot{{mathcal {M}}}_mathbf{0}(t)&= -k_text {p} {mathcal {M}}_mathbf{0}(t)left[ {mathcal {M}}^*_{ mathbf{0}}(t) + sum _{ varvec{kappa }in {varvec{k}}_{1,0}} {mathcal {M}}_{ varvec{kappa }}(t)right] ,\ dot{{mathcal {M}}}^*_{ mathbf{0}}(t)&=-k_text {p} {mathcal {M}}_mathbf{0}(t) {mathcal {M}}^*_{ mathbf{0}}(t),\ end{aligned}$$

(31)

For monomers with one out-edge of arbitrary colour, we have

$$begin{aligned} dot{{mathcal {M}}}_{{varvec{k}}}(t)=k_text {p} {mathcal {M}}_mathbf{0}(t) {mathcal {M}}^*_mathbf{0}(t)Theta (t,c), end{aligned}$$

(32)

for all ({{varvec{k}}} in {varvec{k}}_{0,1}). For monomer species with one in-edge of arbitrary colour, we have:

$$begin{aligned} dot{{mathcal {M}}}_{{varvec{k}}}(t)=&-k_text {p} {mathcal {M}}_mathbf{0}(t) {mathcal {M}}_{{varvec{k}}}(t) \&+ k_text {p}{mathcal {M}}_mathbf{0}(t) left[ {mathcal {M}}^*_mathbf{0}(t) + sum _{varvec{kappa }in {varvec{k}}_{1,0}} {mathcal {M}}_{ varvec{kappa }}(t) right] Theta (t,c), end{aligned}$$

(33)

for all ({varvec{k}} in {varvec{k}}_{1,0}). For monomer species with one in-edge of arbitrary colour and one out-edge of colour *c* we have:

$$begin{aligned} dot{{mathcal {M}}}_{{varvec{k}}+Delta {varvec{k}}_{0,1}(c)}(t)= k_text {p} {mathcal {M}}_mathbf{0}(t) {mathcal {M}}_{{varvec{k}}}(t) Theta (t,c), end{aligned}$$

(34)

for all ({varvec{k}}in {varvec{k}}_{1,0}). The initial conditions are given by ({mathcal {M}}_mathbf{0}(0)={mathcal {M}}_0(1-R_0)), ({mathcal {M}}^*_mathbf{0}(0)={mathcal {M}}_0R_0) and ({mathcal {M}}_{ {varvec{k}}}(0)=0) for ({varvec{k}} ne mathbf{0}).

From the solution of the ODE system (29)–(34), the 2*N*-variate degree distribution at time *t* is computed:

$$begin{aligned} u({varvec{k}})=frac{{mathcal {M}}_{{varvec{k}}}(t)+mathcal M^*_{{varvec{k}}}(t)}{sum _{{varvec{k}}} [mathcal M_{{varvec{k}}}(t)+{mathcal {M}}^*_{{varvec{k}}}(t)]}, end{aligned}$$

(35)

where conversion (chi (t)) is defined as

$$begin{aligned} chi (t)=1-frac{{mathcal {M}}_mathbf{0}(t)}{{mathcal {M}}_mathbf{0}(0)}, end{aligned}$$

(36)

See Table 1 for the numerical values used in the examples.

### Mathematical description: free-radical photopolymerisation

We consider the photopolymerisation of mono- and diacrylates and derive the differential equations for the concentration of monomer species at *N* discrete time intervals in a similar way as in the preceding section. The concentrations are used to construct the coloured degree distribution.

#### ODEs with colours for free-radical photopolymerisation

The reaction mechanism for polymer species can be found in Ref.^{13}. We characterise the state of a monomer unit ({mathcal {M}}_{v,r,{varvec{k}}}) by its number of vinyl groups *v*, radicals *r*, and the edge count vector ({varvec{k}}=(k_1,dots ,k_{2N})) where *N* is the number of time intervals. The covalent bond formed during polymerisation is represented by a directed edge oriented from the initial radical site to the initial vinyl site. In the case of photopolymerisation of mono- and diacrylates, the elements of the count vector are restricted to ( k_iin {0,1,2}, , i=1,dots ,2N), as every monomer unit can form at maximum two in-edges and two out-edges.

The reaction mechanism is summarised as follows:

Photoinitiation of the initiator (I_2) and formation of two initiator radicals

*I*per initiator:$$begin{aligned} I_2 xrightarrow {k_text {d}} 2I. end{aligned}$$

Initiation of a vinyl group by reaction with the initiator radical:

$$begin{aligned} I+ {mathcal {M}}_{v,r,{varvec{k}}} xrightarrow {vk_text {ini}} {mathcal {M}}_{v-1, r+1,{varvec{k}}}. end{aligned}$$

Propagation (including crosslinking) by the reaction of the radical with a vinyl group and formation of an edge of colour

*c*:$$begin{aligned}&{mathcal {M}}_{v,r,{varvec{k}}}+{mathcal {M}}_{v’,r’,{varvec{k}}’}xrightarrow {k_text {p}} \&,,,,,,,,,,,mathcal M_{v,r-1,{varvec{k}}+Delta {varvec{k}}_{0,1}(c)} + mathcal M_{v’,r’,{varvec{k}}’+Delta {varvec{k}}_{1,0}(c)}. end{aligned}$$

Termination by disproportionation by the reaction of two radicals:

$$begin{aligned}&{mathcal {M}}_{v, r,{varvec{k}}} + mathcal M_{v’,r’,{varvec{k}}’} xrightarrow {rr’k_text {td}} mathcal M_{v, r-1,{varvec{k}}} + {mathcal {M}}_{v’, r’ -1,{varvec{k}}’}. end{aligned}$$

The corresponding system of ODEs reads as:

$$begin{aligned} begin{pmatrix} dot{{mathcal {M}}}_{v,r,{varvec{k}}}(t)\ {dot{I}}_2(t)\ {dot{I}}(t)\ end{pmatrix} = begin{pmatrix} G_{v,r,{varvec{k}}}(t)\ -k_text {d} I_2(t), \ 2k_text {d} I_2(t)-k_text {i} I(t) c_v(t)\ end{pmatrix}, end{aligned}$$

(37)

where

$$begin{aligned} G_{v,r,{varvec{k}}}(t)&:= k_text {ini}I(t)big [ (v+1) {mathcal {M}}_{v+1,r-1,{varvec{k}}}(t) – v {mathcal {M}}_{v,r,{varvec{k}}}(t) big ]\&quad +k_text {p}c_r(t)big [(v+1) {mathcal {M}}_{v+1,r-1,{varvec{k}}-Delta {varvec{k}}_{1,0}(c)}(t) \&quad – v {mathcal {M}}_{v,r,{varvec{k}}}(t) big ] \&quad + k_text {p}c_v(t)big [(r+1) {mathcal {M}}_{v,r+1,{varvec{k}}-Delta {varvec{k}}_{0,1}(c)}(t) \&quad – r {mathcal {M}}_{v,r,{varvec{k}}}(t) big ] \&quad +2k_text {td}c_r(t)big [(r+1) {mathcal {M}}_{v,r+1,{varvec{k}}}(t) – r {mathcal {M}}_{v,r,{varvec{k}}}(t) big ]. end{aligned}$$

(38)

The radical concentration (c_r(t)) and vinyl concentration (c_v(t)) in the system are given by

$$begin{aligned}&c_r(t)=sum _{v,r,{varvec{k}}} r{mathcal {M}}_{v,r,{varvec{k}}}(t),\&c_v(t)=sum _{v,r,{varvec{k}}} vmathcal M_{v,r,{varvec{k}}}(t). end{aligned}$$

(39)

The polymerisation progress is quantified by the vinyl conversion

$$begin{aligned} chi (t)=1-frac{c_v(t)}{c_v(0)}, end{aligned}$$

(40)

with (c_v(0)) the initial vinyl concentration. The remaining initial concentrations are given as (I_2(0)), *I*(0), and (mathcal M_{v,r,{varvec{k}}}(0)=0) for all (v,,r,,{varvec{k}}) except for ({mathcal {M}}_{2,0,{varvec{k}}_{0,0}}(0)=c_text {di}mathcal M_0) and (mathcal M_{1,0,{varvec{k}}_{0,0}}(0)=c_text {mono}{mathcal {M}}_0). The fractions (c_text {di}) and (c_text {mono}) refer to the fractions of di- and monovinyls in the systems. All numerical values are presented in Table 1.

After solving the system of ODEs, the degree distribution is computed:

$$begin{aligned} u({varvec{k}})=frac{sum _{v,r}mathcal M_{v,r,{varvec{k}}}(t)}{sum _{v,r,{varvec{k}}}mathcal M_{v,r,{varvec{k}}}(t)}. end{aligned}$$

(41)

#### Special case: cross-linked living polymerisation

Living polymerisation is a simplified case of the free-radical polymerisation: As living polymerisation does not include initiation and termination reaction, the respective rate constants are set to zero, (k_text {ini}=0, ,k_text {td}=0). Also, initiator and initiator radicals are not present in the system, hence (I_2(0)=0), (I(0)=0). After applying these changes, the system of ODEs (37) simplifies to

$$begin{aligned} dot{{mathcal {M}}}_{v,r,{varvec{k}}}(t)&= k_text {p} c_r(t) big ((v+1) {mathcal {M}}_{v+1,r-1,{varvec{k}}-Delta {varvec{k}}_{1,0}(c)}(t) \&quad – v {mathcal {M}}_{v,r,{varvec{k}}}(t) big ) \&quad + k_text {p} c_v(t) big ((r+1) {mathcal {M}}_{v,r+1,{varvec{k}}-Delta {varvec{k}}_{0,1}(c)}(t) \&quad – r {mathcal {M}}_{v,r,{varvec{k}}}(t) big ). end{aligned}$$

(42)

The initial conditions are given by

({mathcal {M}}_{2,0,{varvec{k}}_{0,0}}(0)=c_text {di}(1-R_0)^2{mathcal {M}}_0),

({mathcal {M}}_{1,1,{varvec{k}}_{0,0}}(0)=c_text {di}2R_0(1-R_0)){mathcal {M}}_0),

({mathcal {M}}_{0,2,{varvec{k}}_{0,0}}(0)=c_text {di}R_0^2{mathcal {M}}_0),

({mathcal {M}}_{1,0,{varvec{k}}_{0,0}}(0)=c_text {mono}(1-R_0){mathcal {M}}_0),

({mathcal {M}}_{0,1,{varvec{k}}_{0,0}}(0)=c_text {mono}R_0{mathcal {M}}_0),

({mathcal {M}}_{v,r,{varvec{k}}}(0)=0) elsewhere,

with (c_text {di}) and (c_text {mono}) the fraction of di- and monofunctional monomers and (c_text {di}+c_text {mono}=1).

### Randomised matrix (tilde{{varvec{M}}})

Matrix (tilde{{varvec{M}}}) is defined in a similar manner as ({varvec{M}}) in Eq. (4),

$$begin{aligned} {tilde{M}}_{i,j}=frac{{tilde{mu }}_{i,j}}{{tilde{mu }}_{j}}-delta _{i,j},;i,j=1,dots ,2N end{aligned}$$

(43)

with ({tilde{mu }}_{i,j}) defining the second mixed moments and ({tilde{mu }}_{j}) the first moments of a 2*N*-variate degree distribution with randomised colours. The randomised degree distribution is based on the uncoloured directed degree distribution (u(k_1,k_2)) with every half-edge having equal probability of having any colour *i* with (i=1,dots ,N), which results in a multinomial distribution.

In order to mathematically describe the randomised system, we define the uncoloured directed degree distribution as the distribution that is obtained from an originally coloured degree distribution by neglecting the colour, but retaining the orientation:

$$begin{aligned} u(x,y)=sum _{Omega _{x,y}} u({varvec{k}}), end{aligned}$$

(44)

where the summation is performed over

$$begin{aligned} Omega _{x,y}={mathbf{k} |sum _{i>0} k_{2(i-1)}=x, sum _{j>0} k_{2j}=y } end{aligned}$$

containing all configurations of vector (mathbf{k}) with *x* in-edges and *y* out-edges regardless their colour.

The first moments of the randomised degree distribution are given by

$$begin{aligned} {tilde{mu }}_{2(i-1)}&=sum _{k_1,k_2}sum _{k_1’=0}^{k_1} k_1′ u(k_1,k_2) \× left( {begin{array}{c}k_1\ k_1’end{array}}right) left( frac{1}{N}right) ^{k_1′} left( frac{N-1}{N}right) ^{k_1-k_1′}, end{aligned}$$

(45)

$$begin{aligned} {tilde{mu }}_{2i}&=sum _{k_1,k_2}sum _{k_2’=0}^{k_2} k_2′ u(k_1,k_2) \× left( {begin{array}{c}k_2\ k_2’end{array}}right) left( frac{1}{N}right) ^{k_2′} left( frac{N-1}{N}right) ^{k_2-k_2′}. end{aligned}$$

(46)

In order to characterise the second mixed moments, we distinguish between the second moments with same half-edge types,

$$begin{aligned} {tilde{mu }}_{2(i-1),2(i-1)}=&sum _{k_1,k_2}sum _{k_1’=0}^{k_1} k_1’^2 u(k_1,k_2) \× left( {begin{array}{c}k_1\ k_1’end{array}}right) left( frac{1}{N}right) ^{k_1′} left( frac{N-1}{N}right) ^{k_1-k_1′}, end{aligned}$$

(47)

$$begin{aligned} {tilde{mu }}_{2i,2i}=&sum _{k_1,k_2}sum _{k_2’=0}^{k_2} k_2’^2 u(k_1,k_2) \× left( {begin{array}{c}k_2\ k_2’end{array}}right) left( frac{1}{N}right) ^{k_2′} left( frac{N-1}{N}right) ^{k_2-k_2′}, end{aligned}$$

(48)

the second moments with same half-edge orientation,

$$begin{aligned} {tilde{mu }}_{2(i-1),2(j-1)}=&sum _{k_1,k_2}sum _{k_1’=0}^{k_1} sum _{k_1”=0}^{k_1′} k_1”(k_1′-k_1”) u(k_1,k_2) \× left( {begin{array}{c}k_1\ k_1’end{array}}right) left( frac{1}{N}right) ^{k_1′} left( frac{N-1}{N}right) ^{k_1-k_1′}, end{aligned}$$

(49)

$$begin{aligned} {tilde{mu }}_{2i,2j}=&sum _{k_1,k_2}sum _{k_2’=0}^{k_2} sum _{k_2”=0}^{k_2′} k_2”(k_2′-k_2”) u(k_1,k_2) \× left( {begin{array}{c}k_2\ k_2’end{array}}right) left( frac{1}{N}right) ^{k_2′} left( frac{N-1}{N}right) ^{k_2-k_2′}, end{aligned}$$

(50)

and the second moments with different types and orientation,

$$begin{aligned} {tilde{mu }}_{2(i-1),2j}=&sum _{k_1,k_2}sum _{k_1’=0}^{k_1}sum _{k_2’=0}^{k_2} k_1’k_2′ u(k_1,k_2) \× left( {begin{array}{c}k_1\ k_1’end{array}}right) left( frac{1}{N}right) ^{k_1′} left( frac{N-1}{N}right) ^{k_1-k_1′}\× left( {begin{array}{c}k_2\ k_2’end{array}}right) left( frac{1}{N}right) ^{k_2′} left( frac{N-1}{N}right) ^{k_2-k_2′}, end{aligned}$$

(51)

with ({tilde{mu }}_{2(i-1),2j}={tilde{mu }}_{2j,2(i-1)}).