The Mathematics of Ethylene Oligomerisation and Polymerisation

Linear α-olefins or LAOs are produced by the catalytic oligomerisation of ethylene on a multimillion ton scale annually. A range of LAOs is typically obtained with varying chain lengths which follow a distribution. Depending on the catalyst, various types of distributions have been identified, such as Schulz–Flory, Poisson, alternating and selective oligomerisations such as ethylene trimerisation to 1-hexene and tetramerisation to 1-octene. A comprehensive mathematical analysis for all oligomer distributions is presented, showing the relations between the various distributions and with ethylene polymerisation, as well as providing mechanistic insight into the underlying chemical processes.


Introduction
Linear alpha olefins (LAOs, or 1-alkenes) find extensive use as co-monomers in olefin polymerisation and as intermediates to detergents, lubricants, and plasticisers [1][2][3]. Global consumption of LAOs grew at an average annual rate of 5.6% from 2012 to 2016, and global capacity for LAO production has been estimated at 6.2 million tonnes per year (mtpy) for 2018 [4,5]. LAOs are generally prepared via the oligomerisation of ethylene, as is shown below in Scheme 1 [6].
LAO processes may be grouped into two types: widerange and on-purpose [7]. Wide-range processes give a distribution of LAO products from C 4 to C 20+ , whilst onpurpose processes predominantly give C 4 , C 6 or C 8 LAOs [7]. Examples of catalytic systems used in industry are given in Figs. 1 and 2 [6][7][8][9].
Three of the most common types of distributions from ethylene oligomerisation and polymerisation processesnamely Schulz-Flory, selective-LAO, and Poisson-are observed with the systems in Figs. 1 and 2. In addition, a number of alternating LAO distributions with chromiumbased systems have been reported in the literature [11][12][13][14][15][16][17]. These four types of ethylene product distribution are depicted in Fig. 3.
Following on from the ground-breaking work 80 years ago by Schulz and Flory [18,19], novel mathematical approaches for analysing product distributions from ethylene oligomerisation and polymerisation processes are presented in this paper. These approaches provide straightforward, applicable methods for characterising key mechanistic features of underlying catalytic processes, and demonstrate how Schulz-Flory, alternating, selective-LAO, and Poisson distributions may be considered as parts of a connected product distribution landscape.

First-Order Oligomerisation Processes
In 1935, Schulz characterised the mathematics behind firstorder oligomerisation processes [18]. He noted the direct proportionality between mol(n)% (mol percentage for an oligomer of n ethylene units) and the product of a sequence featuring a constant propagation probability α: Shell (SHOP) 7 C 4 to C 20+ Capacity: ∼1.7 mtpy Chevron Phillips (Gulfene) 8 C 4 to C 20+ Capacity: ∼1.0 mtpy INEOS (Ethyl) 8 C 6 to C 14 Capacity: ∼0.8 mtpy The relationship between wt(n)%, mol(n)% and n is: where (RFM C 2 H 4 ) is the Relative Formula Mass of ethylene.
Using the expression for mol(n)% from Eq. 1a and condensing: Again making use of a continuous function approximation (see SI.1.2 for derivation), this implies: The following year, in 1936, Flory presented precise derivations without making use of continuous function approximations [19]. Flory noted the dependence of mol(n), the amount of moles for each oligomer of n ethylene units, on the product of (n − 1) propagation steps and 1 elimination step along with R, the total mol of oligomers produced during the reaction time: Consequently: It follows that: which gives (see SI. 1

.3 for derivation):
Comparing Eq. 1b with Eq. 3b, it is evident that both Schulz's equation and Flory's equation express: The two Eqs. 1b and 3b are subtly different, as shown in Fig. 4, because −(ln ) is the proportionality constant for Schulz and (1 − )∕ is the proportionality constant for Flory. However, in practice, this makes little difference because is typically determined by taking natural logarithms to each side of Eq. 5a, plotting ln(mol(n) %) vs. n, performing a linear regression analysis to find the gradient lnα and then exponentiating the result as a power of e [20].
The value for determined in this way will be the same in both cases.
By comparing Eq. 2c with Eq. 4b, it may be seen that: with (ln ) 2 as the proportionality constant for Schulz and (1 − ) 2 ∕ as the proportionality constant for Flory (Table 1).
In the next section, a new derivation for Eq. 3a will be presented which utilises recurrence relations, followed by the presentation of a new back-solving method to determine α and R. The analysis lays the groundwork for an analogous consideration of second-order oligomerisation processes.

Derivation of the Schulz-Flory Equation via Recurrence Relations
Catalytic oligomer isation systems which give Schulz-Flory LAO distributions propagate either via a Cossee chain growth mechanism or via a metallacyclic mechanism (Schemes 2 and 3). Both mechanisms may be generically represented by Scheme 4. Deuterium labelling studies have shown that many chromium systems operate via the metallacyclic pathway [21][22][23]. In Scheme 4: • R is the total amount in mol of oligomers produced (from T 1 onwards) during the reaction time i.e. the rate at which intermediate species arrive at M(C 2 ) n , in units of mol t −1 .
• α is the propagation probability from M(C 2 ) n to M(C 2 ) n+1 etc.
• So (1 − α) is the elimination probability from M(C 2 ) n etc. • T n is the total mol of a particular LAO with chain length (C 2 ) n produced during the reaction time (i.e. the rate at which that particular LAO is eliminated, in units of mol t −1 ), with T 1 as the starting point in the sequence.
Expressions for each LAO fraction may be found by inspection of Scheme 4, or by using a recurrence relation methodology. These two techniques are demonstrated in the centre and right columns of Table 2 respectively.
The aim is to express T n in terms of α and R. From inspection of the central column of Table 2, it may be inferred that the expression is:

Scheme 4
Generic scheme for Schulz-Flory LAO distributions. Active species shown in black, rates and eliminated products shown in blue, probabilities of propagation and elimination pathways shown in green Table 2 Expressions for T n via inspection (centre) or via recurrence relation (right)

Rate at which species arrive
Probability of eliminaƟon With constant values for α and R, it is implicitly assumed that: • The sequence begins at a suitable T 1 (e.g. a longer 1-alkene). • The oligomerisation proceeds under steady state conditions (i.e. negligible change in ethylene pressure, reaction temperature and quantity of catalytic species). • Under steady state conditions, the standing quantities of each catalytic species are constant. The sum of rates to any particular active species "box" in Scheme 4 are equal to the sum of rates from that "box". This is equivalent to saying that the sum of probabilities for any particular active species "box" is 1 e.g. for M(C 2 ) n R = R + (1 − )R ↔

Back-Solving to Determine ˛ and R, and Checking for the Best Fit
In practice, the T n terms are known (i.e. the amount in mol of each oligomer produced during the reaction time), and the aim is to find α and R. See SI.1.5 for a discussion on the method used to determine α and R, as well as a checking method to determine that the best fit has been found.

Application of the Fitting Procedure to Experiment
Bis(benzimidazole)pyridine chromium(III) chloride, in combination with MAO co-catalyst, gives a Schulz-Flory distribution of LAOs via a metallacyclic mechanism (Fig. 5) [24]. Fitting α either via a ln(mol(n)%) vs. n plot or via the recurrence relation method gives a resulting value of α = 0.84.
It is worth noting that the determination of α via the recurrence relation method offers a subtle advantage over ln(mol(n)%) vs. n plots, in that accurate fits can only be obtained for distributions with reasonably constant α values as per Scheme 4. Linear regression analysis via the ln(mol(n)%) vs. n method may permit adequate fits for data sets displaying non-constant α behaviour, leading to Schulz-Flory equations which deviate from experimental plots.

Derivations of Further Experimentally-Relevant Metrics
For polymer distributions, M n , M w , and PDI are defined as follows: These definitions can be applied to Schulz-Flory distributions. Using Eq. 3a, 4a and 4b the following expressions are obtained (see SI.1.6 for derivations):

Second-Order Oligomerisation Processes
In 2006, Tomov et al. discovered a highly active chromiumbased oligomerisation system, with activities in excess of 100,000 g mmol −1 h −1 bar −1 , which gave an alternating distribution of 1-alkenes via a metallacyclic mechanism (Fig. 6) [15]. The shape of this distribution cannot be described by first-order linear homogeneous recurrence relations with constant coefficients e.g. the Schulz-Flory equation. Variations on the first-order theme with orderly coefficients are likewise unable to accurately describe the experimentallyobserved results. For instance, in a hypothetical alternating coefficients case, the quantities of the series of oligomers 1-hexene, 1-decene, 1-tetradecene etc. always decrease with n. Yet the distribution in Fig. 6 evidently displays some sort of periodicity.
All oligomerisations may be mathematically considered in terms of recurrence relations, so from an analytic point of view it is logical to proceed from a first-order treatment to a second-order treatment. Chemically speaking, this implies a mechanism involving both single and double coordination and insertion of ethylene. The double coordination and insertion of olefins has been proposed previously. Ystenes proposed a "trigger" mechanism whereby a coordinated olefin only inserts once a second olefin has coordinated to the metal [25]. Lemos and co-workers proposed a related mechanism based on kinetic studies of Ziegler-Natta catalysts [26], and a similar mechanism has been implicated in other oligomerisation processes [27,28].
Building on the analysis from Sect. 2.1, a new "alternating equation" is derived here, which may be viewed as a direct analogue to the Schulz-Flory equation for secondorder oligomerisation processes and which accurately describes alternating distributions as per Fig. 6. This is followed by the presentation of a new back-solving method for the determination of key mechanistic parameters.

Derivation of the Alternating Equation via Recurrence Relations
In Scheme 5: • R 1 is the rate at which intermediate species arrive at M(C 2 ) n , in units of mol t −1 .
• R 2 is the rate at which intermediate species originating from before M(C 2 ) n arrive at M(C 2 ) n+1 , in units of mol t −1 .
• So the total amount in mol of oligomers, R, produced (from T 1 onwards) during the reaction time is expressed by R = R 1 + R 2 • α is the propagation probability from M(C 2 ) n to M(C 2 ) n+1 etc. • β is the propagation probability from M(C 2 ) n to M(C 2 ) n+2 etc.
• T n is the total mol of a particular 1-alkene with chain length (C 2 ) n produced during the reaction time (i.e. the (1-) etc.
rate at which that particular 1-alkene is eliminated, in units of mol t −1 ), with T 1 as the starting point in the sequence.
With reference to Scheme 5, Table 3 expresses each LAO via inspection (central column) or via recurrence relation (right column). The aim is to express T n in terms of α, β, R 1 and R 2 . In contrast to Sect. 2.1, obtaining a general expression via inspection is no longer feasible. The method of solving via a recurrence relation is much more straightforward, as will be demonstrated below.
The recurrence relation is: Following a derivation (see SI.2.1), the general solution to this expression is: [16] with: i.e. the unnormalised alternating equation.
With constant values for α, β, R 1 and R 2 , it is implicitly assumed that: • The sequence begins at a suitable T 1 (e.g. a longer 1-alkene) • The oligomerisation proceeds under steady state conditions (i.e. negligible change in ethylene pressure, reaction temperature and quantity of catalytic species) • Under steady state conditions, the standing quantities of each catalytic species are constant. The sum of rates to any particular active species "box" are equal the sum of rates from that "box". This is equivalent to saying that the sum of probabilities from any particular active species "box" is 1 e.g. for M(C 2 ) n To convert the unnormalised alternating equation, Eq. 11, into the normalised form:

Back-Solving to Determine α, β, R 1 , and R 2 , and Checking for the Best Fit
In practice, the T n terms are known (i.e. the amount in mol of each oligomer produced during the reaction time), and the aim is to find α, β, R 1 , and R 2 . See SI.2.4 for a discussion on the method used to determine α and R, as well as a checking method to determine that the best fit has been found.

Application of the Fitting Procedure to Experiment
See Sect. 2.4 for examples of fittings using the described recurrence relation method with

Rate at which species arrive
Probability of eliminaƟon

Derivations of Further Experimentally-Relevant Metrics
As per Schulz-Flory distributions, the definitions of M n , M w and PDI in Eqs. 9a, 9b, 9c can be applied to alternating distributions. Using Eqs. 11, 14 and 15, and Eqs. 12ii and 14iii in SI, the following expressions are obtained (see SI.2.5 for derivations): with: Note that with q = 0, y = 1, and r 1 = , PDI matches the expression previously obtained with Schulz-Flory distributions (see Eq. 9f): As with Schulz-Flory distributions, PDI values for alternating distributions lie between 1 and 2. This may be demonstrated by examining four limiting cases (see SI.2.6).

Connecting First-Order and Second-Order Oligomerisation Processes
By inspection of Eq. 11, two cases may be identified where the second-order general solution collapses to a first-order solution as per Eq. 10. In the first case, if β = 0, then from Eq. 10iv r 1 = and from Eq. 10v r 2 = 0 . With r 2 = 0 , Eq. 10xiii gives R 2 = 0 . Then, using Eqs. 10x and 10xi (in SI.2.1 and SI.2.4): In this sub-section, the mathematics of the second case is illustrated, followed by a method to describe the degree of similarity between second-order recurrence relations and first-order ones. This analysis has important implications for the fitting of more "gently" alternating oligomer distributions.

The First-Order Surface in R 2 /R 1 , α and β Space
Equation 20 features three variables, R 2 /R 1 , and β, which are physically interpretable as per Scheme 5. Conveniently, then, it is possible to plot in 3D space the combinations of R 2 ∕R 1 , and β values which give rise to first-order behaviour. Furthermore, since r 2 does not feature in Eq. 22 and since β is not necessarily 0, Eq. 10iv (Supplementary Information SI.2.1) now operates as a parametric equation where different combinations of α and β may give the same result for r 1 . A visualisation of the first-order surface in R 2 /R 1 , and space, constructed by lines along which r 1 is constant, is shown in Fig. 7 and 8. Note that first-order distributions following Eq. 18 are restricted to the alpha axis with β = 0 and R 2 /R 1 = 0.

Degree of Alternation
By sight, it is evident that second-order distributions can possess varying degrees of alternating character. To arrive at a definition for the degree of alternation in a second-order distribution, we need to consider the pictorial representation of Eq. 11. Figure 9 illustrates how an alternating distribution is mathematically constructed by adding a constantly decreasing first-order recurrence relation to an alternating first-order recurrence relation (recall that with 0 ≤ α < 1, 0 ≤ β < 1, and 0 < (1 − α − β) < 1, then 0 < r 1 < 1, The degree of alternation is defined as: Absolute values are taken in the numerator to avoid cancellation of terms on summation (refer to Fig. 9) Fig. 7 Visualisation of the first-order surface in R 2 ∕R 1 , and space along the axis (left) and along the axis (right). Colours correspond to particular values of r 1 (key centre) Deg. of Alt. may be considered as a percentage).
As r 1 and r 2 may be expressed in terms of α and β, Eq. 23b gives the degree of alternation as a function of R 2 /R 1 , and β. In other words, every point in R 2 /R 1 , and β space has an associated degree of alternation. A visualisation of the degree of alternation map in R 2 /R 1 , and β Fig. 10. Note that although 0 ≤ R 2 ∕R 1 , for convenience Fig. 10 runs only from 0 ≤ The red areas in Fig. 10 represent those areas with the largest degree of alternation, i.e. where R 2 ∕R 1 or are highest. The light blue layer corresponds to the first-order surface shown in Figs. 7 and 8, where distributions display no degree of alternation. This layer acts as a crossing-point, above which distributions display "normal" alternating behaviour, and below which distributions display "inverse" alternation. Figure 11 illustrates these results. Note that in Fig. 11, and remain constant while R 2 ∕R 1 varies. While experimentally so far only "normal" alternating distributions  have been found, the discovery of "inverse" alternating behaviour is theoretically possible.

Implications for Second-Order Fitting with "Gently" Alternating Distributions
"Gently" alternating distributions display significantly less c 2 r n 2 character than c 1 r n 1 character. Such distributions lie close to the first-order surface in R 2 ∕R 1 , and space as illustrated in Figs. 7 and 8. In back-solving to determine , , R 1 , and R 2 , "gently" alternating distributions present a unique difficulty in the fitting process when experimental error is taken into account.
The first-order surface in Figs. 7 and 8 is constructed by lines along which r 1 is constant and along which distributions appear identical. Second-order distributions which are similar in appearance, but which lie close to the first-order surface, may accurately be described using quite different R 2 ∕R 1 , and values.
By way of illustration, consider the hypothetical secondorder distribution shown in Fig. 12, as may be detected via GC from 1-hexene onward. If, due to experimental error, a quantity of 1-hexene evaporates prior to GC measurement, then the back-solving method produces a different fit for the data set slightly perturbed at C 6 .
Comparing the fitted parameters in Fig. 12 with those in Fig. 13, there are significant differences between values for , , R 2 ∕R 1 , c 2 and r 2 . The only reliable metric for "gently" alternating distributions is the first-order c 1 r n 1 term. In borderline cases where a second-order distribution lies too close to the first-order surface to obtain reliable values for , , R 1 , and R 2 , pressure studies (see Sect. 2.4) and/or repeat experiments will be necessary.

The Effect of Pressure in Second-Order Steady-State Oligomerisation
In Sect. 2.2, a catalytic scheme (Scheme 5) for the production of alternating 1-alkene oligomer distributions was presented. The terms M C 2 n , M C 2 n+1 etc. group together all intermediates involving n ethylene units. In order to investigate the effect of changing ethylene pressure on and , it is necessary to propose a more detailed scheme with the individual elementary reaction steps that occur within each active species "box" (i.e. between each type of M C 2 n species). This so-called candidate scheme is kept as general as possible, for application to any processes proceeding via a metallacyclic mechanism (Scheme 6): Following a derivation (see SI.4.1), it can be shown that the ratio between the probabilities for double ethylene insertion β and single insertion α is defined as: So long as k 2 , k 7 ≠ 0 , the relationship between the probabilities and and ethylene pressure P et can be expressed as:    Table 4 and Fig. 14.
The fitting of the data in Table 4 with Eq. 24b gives values for g 1 = 0.177 and g 2 = 0.176 . Note that g 2 is equivalent to the ratio of rate constants for single insertion vs. double insertion from a double-ethylene associated intermediate. There is a good degree of correlation between experimentally determined and the fitted values of ∕ in Fig. 14 (goodness of fit R 2 = 0.974).
A final note should be added, that for ethylene oligomerisations at lower pressures with these catalysts, it is possible that the most abundant oligomeric product, 1-butene in this case, is incorporated during the oligomerisation reaction. We have observed this with some of the most active catalysts, resulting in a minor distribution of 2-ethyl-1-alkenes with the opposite alternating pattern, alongside the major distribution of 1-alkenes with the normal alternating behaviour [17].

Steady-State Oligomerisation with Nonconstant Propagation Probabilities
In assuming constant values for the propagation probabilities and in Schemes 4 and 5 (Sects. 2.1 and 2.2 respectively), recurrence relations with constant coefficients are generated.
The solution of these may be expressed in equation format (see Eqs. 8 and 11). If, however, the propagation probabilities adopt nonconstant values in the catalytic cycle, then recurrence relations with variable coefficients are generated. The aim in such cases is to determine a unique value for each unknown parameter, and/or to elucidate the key mechanistic features of the system.
In this sub-section, a method for handling first-order recurrence relations with variable coefficients is presented. The method finds application in oligomerisation systems where changes at some point in the catalytic cycle. Second-order systems with nonconstant propagation probabilities will be discussed later, with a particular focus on selective-tetramerisation.

First-Order Recurrence Relations with Nonconstant Propagation Probabilities
Experimental LAO distributions which deviate from Schulz-Flory have been observed, for example chromium(III) bis(carbene)pyridine systems reported by McGuinness and co-workers ( Fig. 15) [29]. The observation has been rationalised by less favourable elimination at the early stages of metallacyclic growth [29]. In Fig. 15, is variable before C 8 , then adopts a constant value from C 8 onwards.
Data from C 8 (n = 4) onward may be treated with the back-solving method described in Sect. 2.1 to find the constant value and R, the total mol of C 8+ oligomers produced during the reaction time (i.e. the rate at which intermediate To find ′ (i.e. the value at C 6 ) and R ′ , two distinct equations involving the two unknown variables may be identified (note that T ′ 1 is the experimental data point for 1-hexene): These equations sum to give: Given that there are two distinct equations for two unknowns ( ′ and R ′ ), a unique solution does exist. R ′ may be found with Eq. 25c, then substituted into Eqs. 25a or 25b to find ′ . The second non-constant term T ′ 2 (i.e. C 4 ) may now be appended to Scheme 7, and the value at C 4 found using the same methodology. The values calculated for Fig. 15 via this methodology are 0.86, 0.65 and 0.54 from C 4 , C 6 and C 8+ respectively. In other words, for this particular chromium catalyst, the termination probability is less favourable at the early stages and increases as the metallacycle grows from a metallacyclopentane to-heptane andnonane. Once the metallacycle has reached a certain size, the termination probability becomes constant.
Computational studies employing the M06L functional on simplified model compounds have indicated that a singleand double-coordination mechanism operates (see Fig. 17) [30][31][32][33]. Elimination from a chromacyclopentane intermediate to give 1-butene is calculated to be energetically unfavourable [30]. Competing processes to 1-hexene and 1-octene formation diverge at a singly-coordinated chromacyclopentane intermediate 7, with 1-hexene formed predominately via single insertion  [30]. Given the varying differences in energy between the early intermediates in Fig. 17, the selective-tetramerisation process in Fig. 16 may be mathematically considered as a second-order recurrence relation with variable coefficients. Data from C 16 onwards corresponds to the tail end of an alternating distribution, where and values are constant by virtue of negligible geometry and energy differences between larger chromacyclic species. For the purpose of the present discussion, setting aside the difficulties in fitting gently alternating distributions as outlined in Sect. 2.3, the region from C 16 onwards may conceivably be treated with the back-solving method described in Sect. 2.2 to find the constant values for and as well as In an effort to find ′ and ′ (i.e. the and values at C 14 ) and R ′ 1 and R ′ 2 , three distinct equations involving the four unknown variables may be identified (note that T ′ 1 is the experimental data point for 1-tetradecene): These equations sum to give: Scheme 7 Catalytic scheme for Fig. 15 at the boundary of the variable and constant regions. Active species shown in black, rates and eliminated products shown in blue Given that there are three distinct equations for four unknowns ( ′ , ′ , R ′ 1 and R ′ 2 ), a unique solution does not exist in this case. If no further information is available/applied from e.g. theoretical studies, an assumption (such as � � = ) is required as a fourth equation to determine unique solutions for the four unknown variables.

Connecting Schulz-Flory Distributions and Poisson Distributions
The mathematics behind Schulz-Flory ethylene oligomer distributions was first accurately described in 1936 [19].
In 1940, Flory went on to give a kinetic derivation of Poisson's equation for polymer distributions obtained from a monofunctional monomer (ethylene oxide) [34]. Yet the mathematical connection between Schulz-Flory distributions and Poisson distributions has not been fully elucidated. In this section, the relationship between the two is made explicit from first principles. It should be noted that the treatment presented below is described with ethylene as the monomer, but is generally applicable to oligomerisation and polymerisation processes involving the sequential incorporation of monomers in forming linear product chains.

First Principles
The mechanics of a first-order ethylene oligomerisation or polymerisation process whereby propagation proceeds etc.
Scheme 8 Catalytic scheme for Fig. 16 at the boundary of the variable , and constant , regions via metal alkyl or metallacyclic species is summarised in Scheme 9. In this model and the following mathematical treatment, we assume that the rate constants k p and k e are constant over time and independent of chain length. Furthermore, all active species must start the polymerisation instantaneously and there is no catalyst deactivation. M(n) represents catalyst intermediates involving n ethylene units (e.g. M(0) = metal hydride species, M(1) = metal ethyl or metallacyclopropane species, M(2) = metal butyl or metallacyclopentane species etc.), k p is the rate constant of propagation, k e is the rate constant of elimination, and L is the total precatalyst loading.
A general expression for M(n) , as a function of time, may be derived (see SI.5.1): The label "non-steady-state" refers to time dependent terms, whilst the label "steady-state" refers to time independent terms. It is worth noting again that Scheme 9 assumes k p and k e are constant over time and independent of M(n) , i.e. no distinction in propagation or elimination kinetics is made between the various intermediate groups in the catalytic scheme. In addition, by combining M(n) species in forming the kinetics equations, it is implicitly assumed that a single species within each M(n) group dominates. For instance, the M(0) group includes both metal hydrides with a vacant site and ethylene-associated metal hydrides. This assumption is no longer required for the Schulz-Flory case, given that steady-state kinetics operate.

Derivation of the Schulz-Flory Equation for First-Order Steady-State Oligomerisation
Referring to Scheme 9, the scenario where steady state is achieved from t 0 (i.e. instantaneously from the beginning of reaction) with k e ≠ 0 corresponds to the limiting case of firstorder steady-state oligomerisation. In such a case, by ignoring non-steady-state terms, Eq. 27simplifies to: Non-steady-state term Steady-state term (27) Given that elimination occurs at a constant rate, the interest lies only in the sum of linear alpha olefins of chain length n ethylene units continuously eliminated from t 0 to t end : Equation 29a is the unnormalised Schulz-Flory equation, and Eq. 29b is the normalised Schulz-Flory equation. Note in Eq. 29a that the total mol of oligomers produced during the reaction time R = L k e t end . Given that both Lk 2 e t end k p (where L and t end are known) and may be fitted, k p and k e may be inferred.
Lk 2 e t end k p

Derivation of Poisson's Equation for First-Order Living Polymerisation
Referring again to Scheme 9, setting k e = 0 , i.e. there is no chain termination process, corresponds to the limiting case of first-order non-steady-state living polymerisation.
In such a case, Eq. 27 simplifies to: Given that no termination occurs throughout, the interest lies only in the distribution of intermediate species at t end (i.e. the end of reaction), which would give n-alkene products on thermal termination (as per the INEOS (Ethyl) process [8]) or n-alkane products upon hydrolysis. At this point: is Poisson's equation scaled by the precatalyst loading L , and Eq. 31b is Poisson's equation scaled by the factor 1 − e − −1 (which quickly approaches 1 as increases). Note that is the mean number of ethylene units per chain. Given that may be fitted (where t end is known), k p may be inferred.

Non-Steady-State Living Polymerisation with Catalysed Chain Growth
In 1952, Ziegler et al. reported the Aufbaureaktion involving the formation of long chain trialkylaluminium species via stepwise insertion of ethylene monomers into Al-C bonds (starting with triethylaluminium, for instance) [35]. The effect of transition metals, such as colloidal Ni, in suppressing chain growth in the Aufbaureaktion to selectively form 1-butene was subsequently discovered [36]. This prompted a screening of the d-block transitional metals for other metal effects and, in the case of titanium, this led to the first reported transition metal catalysed ethylene polymerisation [37,38]. In their early work, Ziegler et al. described the process as a transition metal catalysed Aufbaureaktion, with the alkylaluminium functioning as the polymerisation catalyst and the transition metal compound functioning as the co-catalyst [37]. Later it was found that the polymer chain grows on the transition metal centre, rather than on aluminium. More recently however, several processes have been reported where the polymer chain is transferred during the course of the reaction from a transition metal to a main group metal such as aluminium, magnesium or zinc [39][40][41][42][43][44]. If this chain transfer is very fast and reversible, then it appears as if the polymer chains grow on the main group metal centres instead. Such a reaction is generally described as a catalysed chain growth (CCG) process or a coordinative chain-transfer polymerisation (CCT), or, to employ the original terminology, a transition metal catalysed Aufbaureaktion [44][45][46].
Tremendous advances have been made in CCG during the past 2 decades, but the mathematics of generating product distributions via CCG has not yet been reported. Such a treatment from first principles is presented here and then applied next.

First Principles
Scheme 10 summarises the mechanics of a first-order ethylene oligomerisation or polymerisation process whereby propagation proceeds via metal alkyl or metallacyclic species, and chain transfer to and amongst a population of non-propagating agents is possible. A similar caveat as in Sect. 3.1 applies in that in this model and the following mathematical treatment, we assume that the rate constants k p , k e and k ct are constant over time and independent of chain length. Furthermore, all active species must start the polymerisation instantaneously and there is no catalyst deactivation.
M(n) represents transition metal intermediates involving n ethylene units, CT(n) represents the population of chain transfer species bearing an alkyl chain of length n ethylene units, k(M(n)) p is the rate constant of propagation for M(n) , k(M(n)) e is the rate constant of elimination for M(n) , and k ct 's are rate constants of chain transfer between species. Scheme 10 is not labelled with individual k ct 's given the general nature of the present treatment.
It is possible to show (see SI.5.2) that: This expression for change does not show dependence on chain transfer terms. In other words, the underlying M(n) system drives towards the distribution that it would adopt in the absence of chain transfer agent. However, the nature of instantaneous steady state for CT species slows down the overall rate at which the M(n) system evolves vs. the case with no chain transfer.

Numerical Simulations
Two illustrative examples are provided to demonstrate the preceding mathematics. Numerical simulations of evolving oligomerisation or polymerisation systems which involve chain transfer, as per Scheme 10, are possible with Microsoft Excel (see SI.5.3 for details).

The Fitting Process
The preceding discussion demonstrated how Poisson distributions or Schulz-Flory distributions are still generated for first-order systems involving CCG. Section 2.1 provided the details for modelling first-order steady state oligomerisation systems, which by extension may be applied to first-order oligomerisation systems involving CCG. Section 3.1 gave a derivation to the Poisson distribution for living polymerisation systems, ending at Eqs. 31a and 31b. Furthermore, k p may now be interpreted as the effective propagation rate constant for the entire system. The process for modelling living polymerisation systems with CCG is described next, along with derivations of the experimentallyrelevant metrics M n , M w and PDI . For illustration, a comparison is also made between an experimentally determined result vs. the fitted data.

Back-Solving to Determine L and λ, and Checking for the Best Fit
In practice, the mol(n) t end values in Eq. 31a are known (i.e. mol of each oligomer produced during the reaction time), and the aim is to find L and . See SI.5.4 for a discussion on the method used to determine and R , as well as a checking method to determine that the best fit has been found.

Derivations of Further Experimentally-Relevant Metrics
For a distribution of n-alkenes, the following expressions can be derived (see SI.5.5 and SI.5.6): As increases, the factor 1 − e − approaches 1. Under such circumstances, Eqs. 35 and 37 simplify.  With high values of , PDI tends towards unity, as per the literature [47].

Comparison of an Experimentally Determined Result vs. Fitted Data
CCG polymerisations have been performed, using ZnEt 2 as the CT agent with 2,6-diacetylpyridinebis(2,6-diisopropylanil)FeCl 2 /MAO (experimental details, see SI.5.7). Termination of the polymerisation reaction is carried out by hydrolysis, resulting in a mixture of long-chain alkanes. Given a highly symmetric Poisson distribution of the resulting n-alkane oligomers, quantitative 13 C{ 1 H} NMR analysis may be used to estimate the mean chain length by comparing integrals for CH 3 with those for CH 2 . For the given experiment, quantitative 13 C{ 1 H} NMR spectroscopy gives a CH 3 :CH 2 ratio of 2:45.38, implying a mean chain length between 46 and 48 carbons (see Fig. 24).
Experimental characterisation of the alkane distribution in such mixtures is non-trivial. The average molecular weight is typically too low for accurate GPC analysis, and the longer alkanes are often insufficiently volatile for conventional quantitative GC-FID analysis. Only the first part of the n-alkane oligomer distribution is usually  Table 5) vs. experimental GC data 1 3 quantifiable via GC. These signals have been used here as inputs for the fitting model. Following the procedure previously described in this sub-section, the values in Table 5 are found.

Model Outputs
There is a high degree of correlation between calculated and experimental values in Table 5 (goodness of fit R 2 = 0.996), and the calculated mean chain length of 46.6 carbons is within the range previously specified by quantitative 13 C{ 1 H} NMR spectroscopy. This method therefore allows the entire product distribution to be determined from a subset of experimental data obtained by GC analysis.

Conclusions and Summary
A variety of different distributions-including Schulz-Flory, alternating, selective-LAO, and Poisson distributions-of linear hydrocarbon products can be obtained from ethylene oligomerisation and polymerisation processes. Oligomerisation processes are steady-state in nature, and can be mathematically described using recurrence relations. First order oligomerisation processes propagate via sequential single-insertions of ethylene, whereas second order oligomerisation processes propagate via sequential single-and double-insertions of ethylene. From an analytical point of view, an important consideration with regards to second order processes is the degree of alternation exhibited by a given distribution, which is generally more pronounced at higher ethylene pressure. The degree of alternation for a particular distribution may be defined as the relative ratio of alternating vs. Schulz-Flory character. An overview of the key analytical results from first order and second order processes with constant probabilities of propagation is given in Table 6.
Examples of first order and second order oligomerisation processes which involve varying probabilities of propagation have been identified from the literature. Such systems include selective oligomerisation catalysts such as the wellknown chromium-based tetramerisation catalysts. An analytical method is available for handling first order processes with varying probabilities of propagation. However, there is a theoretical limitation to a similar approach for analogous second order processes.
The connection between distributions from first order steady-state oligomerisations (Schulz-Flory) and first order non-steady-state polymerisations (Poisson) is made through a consideration of kinetic differential equations. The addition of a chain transfer agent (e.g. ZnEt 2 ) to such systems results in a slowing down of the overall propagation rate. Accurate modelling of a CCG polymerisation system featuring ZnEt 2 as the CT agent has been achieved, allowing for the derivation of several key mechanistic parameters.
The mathematics of ethylene oligomerisation presented here can be easily applied to other systems where chain growth proceeds through an oligomerisation process. Examples of such systems would be the oligomerisation of higher α-olefins to poly(α-olefin) oils, the dehydropolymerisation of aminoboranes, chain growth reactions in the Fischer-Tropsch process or the formation of oligonucleotides. By fitting the experimental results to the formulas described herein for the various oligomerisation scenarios, valuable information may be extracted regarding the underlying reaction mechanisms of the oligomerisation process.