Quantifying the Tetrad Effect, Shape Components, and Ce–Eu–Gd Anomalies in Rare Earth Element Patterns

Abstract

Plots of chondrite-normalised rare earth element (REE) patterns often appear as smooth curves. These curves can be decomposed into orthogonal polynomial functions (shape components), each of which captures a feature of the total pattern. The coefficients of these components (known as the lambda coefficients—$\lambda$) can be derived using least-squares fitting, allowing quantitative description of REE patterns and dimension reduction of parameters required for this. The tetrad effect is similarly quantified using least-squares fitting of shape components to data, resulting in the tetrad coefficients ($\tau$). Our method allows fitting of all four tetrad coefficients together with tetrad-independent $\lambda$ curvature. We describe the mathematical derivation of the method and two tools to apply the method: the online interactive application BLambdaR, and the Python package pyrolite. We show several case studies that explore aspects of the method, its treatment of redox-anomalous REE, and possible pitfalls and considerations in its use.

1 Introduction

Rare earth element (REE) patterns are plots of measured REE abundances in a sample, in which the concentration of each element is divided by its corresponding concentration in a reference material, a process termed normalisation (Coryell et al. 1963). The most common reference material is an average of ordinary chondrite meteorites (CI), often taken to represent an unmodified primordial composition of the solar system and the Earth. Once chondrite-normalised (CN), any deviations from a horizontal array record some geochemical fractionation process. As the geochemical properties of the REE generally vary in a smooth manner from La to Lu (Anenburg 2020 and references therein), the pattern itself is a smooth line whose shape is informative of the processes that led to the sample formation.

Quantification of various shape aspects of patterns is usually achieved using normalised element ratios. For example, $\hbox {La}_\mathrm{{N}}/\hbox {Lu}_\mathrm{{N}}$ (where subscript N indicates normalised values) measures the overall slope of a pattern, and $\hbox {La}_\mathrm{{N}}/\hbox {Sm}_\mathrm{{N}}$ might provide information on the degree of light REE (LREE) fractionation from the middle REE (MREE). However, the choice of elements is often arbitrary and chosen to support a certain hypothesis, or designed for use in simpler environments (e.g., mid-ocean ridge basalts), and can never provide information on the shape on a shorter wavelength than the one used for the ratio (Anenburg 2020). For example, the three different REE patterns shown in Fig. 1 (A, B, C) are clearly different and represent different petrogenetic processes, yet their $\hbox {La}_\mathrm{{N}}/\hbox {Lu}_\mathrm{{N}}$ are nearly identical. Additional measures, such as $\hbox {Dy}/\hbox {Dy}$* (the deviation of $\hbox {Dy}$ from a straight line connecting La and Yb) that attempts to quantify pattern concavity, are only reliable for simple patterns (Davidson et al. 2013). For example, $\hbox {Dy}/\hbox {Dy}$*$=1$ for pattern D in Fig. 1, suggesting a simple straight pattern. Yet, a quick glance at the pattern reveals it to be an artefact of its curvature, in which Dy is located precisely at the La–Yb interpolating line, simply by chance.

In cases of such complexity, descriptions of REE pattern are often supplemented by qualitative terms such as “sinusoidal”, “U-shaped”, or “spoon-shaped”. As these are loosely defined qualitative terms, they are somewhat ambiguous (i.e., one person might understand “spoon-shaped” differently from another) and unquantifiable (to what degree is a pattern spoon-shaped?). The problem is compounded when two or more qualitative shape descriptions are combined. Using normalised element ratios to describe an overall sloping REE pattern which is simultaneously spoon-shaped and sinusoidal becomes a futile endeavour. In such cases a simple visual inspection of the pattern is more informative, as “a picture is worth a thousand words”. However, visual inspections become challenging when many analyses are present and the problem of overplotting becomes serious as many individual lines are entangled into an incomprehensible jumble, from which it is not straightforward to discern trends. For instance, trends in a data set representing a rock suite in which pattern sinusoidalities increase whereas spoon-shapednesses decrease are likely to be obscured by having “a thousand pictures”.

Fig. 1

Four synthetic REE patterns generated using the following parameters: $\Lambda _\mathrm {A}{=}(2.2, 35, 300, 0, -17000)$, $\Lambda _\mathrm {B}=(4.9, 37.5, -290, 0, 0)$, $\Lambda _\mathrm {C}=(3.6, 15, 0, 8800, 0)$, and $\Lambda _\mathrm {D}=(4.6, -15, 225, 7800, 24000)$. The linear line connecting La and Yb in pattern D is plotted in light green

The above-described difficulties are greatly reduced if REE patterns are taken as mathematical functions constructed from simpler functions, each describing a shape, which we term “shape components”. The method was originally developed by O’Neill (2016) and employs least-squares fitting of coefficients to orthogonal polynomials of increasing powers, termed the “lambda shape coefficients”. Here, we explicate the method of using $\lambda$ shape components and coefficients to describe REE patterns, and expand it to include $\tau$ coefficients that describe the tetrad effect—subtle variation in REE behaviours that affects individual groups of four consecutive elements.

2 Mathematical Description of REE Patterns

2.1 Lambda Shape Coefficients

The function used by O’Neill (2016) to describe REE pattern shapes is

$$\begin{aligned} \ln {[\mathrm {REE}_i]_\mathrm {N}}=\lambda _0+\lambda _1 f^\mathrm {orth}_1 +\lambda _2 f^\mathrm {orth}_2 + \lambda _3 f^\mathrm {orth}_3+\lambda _4 f^\mathrm {orth}_4, \end{aligned}$$

(1)

where $\ln [\mathrm {REE}_i]_{\mathrm {N}}$ is the natural logarithm of the normalised element concentration. REE plots are plotted using the common logarithm (i.e., $\log _{10}[\mathrm {REE}_i]_{\mathrm {N}}$), which makes the choice of natural logarithm seem odd; the reasoning for the use of the natural logarithm is presented in O’Neill (2016). In any case, transformations between the two logarithm bases are trivial (see appendix in Anenburg 2020).

In Eq. 1, $f^\mathrm {orth}_n$ are orthogonal polynomial functions of eightfold ionic REE radii $r_i$ (Shannon 1976). The functions of increasing order are as follows

$$\begin{aligned} \begin{array}{ll} f^{\mathrm {orth}}_{1} &{}= (r_{i}-\beta ) \\ f^{\mathrm {orth}}_{2} &{}= (r_{i}-\gamma _{1})(r_{i}-\gamma _2) \\ f^{\mathrm {orth}}_{3} &{}= (r_{i}-\delta _{1})(r_{i}-\delta _2)(r_{i}-\delta _{3}) \\ f^{\mathrm {orth}}_{4} &{}= (r_{i}-\epsilon _{1})(r_{i}-\epsilon _{2})(r_{i}-\epsilon _{3})(r_{i}-\epsilon _{4}).\\ \end{array} \end{aligned}$$

(2)

Fig. 2

Illustration of progressive addition of shape components, using a $\lambda _1$ and subsequent b $\lambda _2$, c $\lambda _3$, and d $\lambda _4$. Dark blue lines are the full patterns described by the $\Lambda$ vector in each panel header, light blue lines are the $\lambda _nf_n^\mathrm {orth}$ components added in each step, and light green lines show the previous full pattern to which the component has been added. For example, $\Lambda _\mathrm {b} = \Lambda _\mathrm {a}+\lambda _2$

The parameters $\beta$, $\gamma _n$, $\delta _n$, and $\epsilon _n$ are pre-calculated constants, such that each $f^\mathrm {orth}_n$ describes a fixed, predetermined shape (see Appendix). The use of ionic radii is preferred to atomic number, as it widens the distance between each LREE relative to HREE. First, each LREE often varies in its properties from its neighbour more than the HREE: in most natural samples, the fractionation degree of Ce differs from Pr more than Ho differs from Tm. This behaviour follows that of ionic radii, and by using them, the greater fractionation degree is accommodated by lower-order shape components. Second, O’Neill (2016) empirically demonstrated that using ionic radii results in a superior fit compared to using atomic numbers. The drawback is that the mathematical treatment becomes slightly more complex, but this is not an issue in nearly all applications if the available software described below is made use of. We follow the convention used by O’Neill (2016) of using the eightfold coordinated “ionic radius” from Shannon (1976), itself based on an earlier structure determination of simple lanthanide fluorides (Greis and Petzel 1974). Whether these values are sufficiently accurate and admissible as representative values for the entire field of geochemistry is still debatable, while being mindful that the “ionic radius” of an element varies based on the local environment and that ions are not rigid balls (Liu et al. 2013; Gibbs et al. 2013). There is ongoing work attempting to refine REE ionic radius values (Gagné 2018), and it is possible that slightly different radii will be used in future versions of the tools presented here.

A complete REE pattern is thus described by the linear combination of the shape coefficients ($\lambda _n$) acting as scaling factors for their respective shape components ($f_n^\mathrm {orth}$). Specifically, $\lambda _0$ represents the average REE abundance, $\lambda _1$ represents the linear slope, $\lambda _2$ represents quadratic curvature, $\lambda _3$ represents cubic curvature (i.e., sinusoidality), and $\lambda _4$ represents the quartic curvature (i.e., W-shape), essentially providing finer control that, when combined with the previous coefficients, allows for various spoon shapes and the like. Figure 2 shows an example of a REE pattern constructed by the gradual addition of shape coefficients. An interactive online app, ALambdaR (https://lambdar.rses.anu.edu.au), allows easy visualisation of the effect of shape coefficients on REE patterns, and readers are encouraged to experiment with it in order to gain an intuitive understanding of the method.

2.2 Tetrad Effects

2.2.1 Background and Previous Studies

In addition to the smooth variation in REE geochemical properties across the lanthanide row, there are cases where the REE are divided into four segments (“tetrads”), each consisting of four elements: (i) La–Ce–Pr–Nd, (ii) Pm–Sm–Eu–Gd, (iii) Gd–Tb–Dy–Ho, and (iv) Er–Tm–Yb–Lu (Kawabe 1992). Within each tetrad, the middle two elements are fractionated relative to the elements at the edges (e.g., Tb and Dy are fractionated relative to Gd and Ho). The fractionation can be positive or negative, which is manifested as a series of ridges or troughs on a REE pattern, occasionally called M-type and W-type, respectively. Various combinations of these types can also be observed (Azizi et al. 2020). The tetrad effect stems from quantum mechanical effects on the electronic structure and energies of the 4f subshell (Kawabe 1992). Note that these M-type and W-type tetrads are distinct from the quartic $f_4^\mathrm {orth}$ of a smooth pattern described by $\lambda _4$, often referred to as M-shaped or W-shaped.

Previous attempts to quantify the tetrad effect included the ratio of the middle elements from their values on an interpolated line connecting the edge elements (Irber 1999), or calculation of the standard deviation of this ratio (Monecke et al. 2002; Abedini et al. 2020). These methods suffer from three major problems. Firstly, the magnitude of the tetrad effect can only be calculated if all four elements are measured and do not present anomalies. This eliminates tetrad 2, as it includes Pm, which does not occur naturally, and Eu, which is commonly anomalous because of its two oxidation states (Eu²⁺ and Eu³⁺). This can also lead to elimination of tetrad 1, as Ce is also occasionally anomalous because of its two oxidation states, Ce³⁺ and Ce⁴⁺, with the latter being common in low-temperature environments; this is especially problematic because the tetrad effect is most likely to occur in surface environments and highly evolved magmatic systems (Bau 1997) up to the low hundreds of degrees Celsius. Secondly, these methods are not suitable for REE patterns which exhibit curvature, namely the vast majority. These interpolation methods will report non-zero tetrads even if the pattern had not been influenced by tetrad effects, simply because curved lines deviate from straight lines (Bau 1997), as these methods essentially calculate the degree of precisely that deviation. Thirdly, tetrads are not symmetric. They arise from the electronic structure of the elements and the stability of empty (La), half-full (Gd), and full (Lu) f shells (Kawabe 1992; Peppard et al. 1969). Additionally, quarter-full f shells should also contribute to the tetrad effect, but the number of f orbitals (14) is not divisible by four, and true quarter-filled shells will only be realised in non-existing elements that occur between Nd and Pm, and Ho and Er (Nugent 1970). This mismatch leads to asymmetry which is not captured by simple linear interpolation on a normalised REE pattern plot. Finally, the two methods based on standard deviations (Monecke et al. 2002; Abedini et al. 2020) may erroneously quantify analytical noise as a tetrad effect if the middle elements are shifted to the same direction.

A superior method for tetrad effect quantification was developed by Minami and Masuda (1997). It involves least-squares fitting of four quadratic curves to the four tetrads, mathematically constraining the positions of the curves such that they intersect correctly at the quarter-filled non-existent element positions using the method of Lagrange undetermined multipliers. The method also solves the problem of missing elements—they are simply excluded from the fitting procedure, using the remaining elements instead. Finally, the method partially solves the issue of curvature by dividing the entire REE pattern into four linear slopes unto which parabolas are superimposed. Unfortunately, the method has largely remained unadopted by the geochemical community, likely due to its mathematical complexity and the lack of accessible computational tools to apply it.

2.2.2 Tetrad Functions Compatible with $\lambda$ Shape Coefficients

A major limitation for tetrad fitting is that the method of Minami and Masuda (1997) cannot be easily applied in conjunction with the REE shape coefficients described above, as the least-squares matrix construction is different. Minami and Masuda (1997) only use a subset of the REE for the fitting of each tetrad, whereas $\lambda$ fitting considers all elements simultaneously. Furthermore, as the tetrad fitting already includes some degree of curvature (albeit in low resolution), there is a conflict between that and some of the $\lambda$ coefficients. Here we adapt their method so that it is compatible with the O’Neill (2016) $\lambda$ shape coefficients, and allows compatible fitting of both shape and tetrad coefficients.

Each tetrad is initially represented by a negative parabola that rises above the x-axis at its respective location and is zero elsewhere, essentially the function

$$\begin{aligned} f^\mathrm {tet}(x)=\max (0,1-x^2)=\frac{1-x^2 + \sqrt{\left( 1-x^2\right) ^2}}{2}. \end{aligned}$$

(3)

Four such parabolas have to be generated, each shifted laterally so that their x-axis intersections occur at the following positions on an atomic number axis: (i) La and Nd–Pm midpoint: 57 and 60.5, (ii) Nd–Pm midpoint and Gd: 60.5 and 64, (iii) Gd and Ho–Er midpoint: 64 and 67.5, and (iv) Ho–Er midpoint and Lu: 67.5 and 71. Their maxima are located between the intersections at 58.75, 62.25, 65.75, and 69.25, respectively. Consequently, $1-x^2$ is transformed to form all four tetrads via $f^\mathrm {tet}\left( g(z)\right)$ where $g(z) = \frac{\left( z-z_0\right) }{1.75}$, $z_0$ are the maxima positions, and 1.75 is half of the parabola width (3.5/2) at $f^\mathrm {tet}=0$. Finally, as the $\lambda$ method expects ionic radii rather than atomic number, a third function that transforms r to z is required. Previous attempts to parameterise r as a function of z have found agreement with Slater’s model for atomic shielding (Raymond et al. 2010; Seitz et al. 2007), but here we use a simple polynomial function $h(r)=z$ that maps the Shannon (1976) radii to atomic number. The polynomial order can be arbitrarily large, and we find that a seventh-order polynomial is sufficient to accurately generate parabolas at the correct locations when the parabolas are calculated using $f^\mathrm {tet}\left( g\left( h\left( r\right) \right) \right)$.

Each tetrad is then described by $\tau _n f^\mathrm {tet}_n$, where $\tau _n$ are each of the four tetrad coefficients, with positive values indicating a negative parabola (i.e., a ridge, or M-type tetrad), and negative values indicating a positive parabola (i.e., a trough, or W-type tetrad).

2.3 Combining $\lambda$ and $\tau$ Shape Coefficients

The sum of all functions described so far leads to the complete function that describes REE patterns

$$\begin{aligned} \ln {[\mathrm {REE}_i]_\mathrm {N}}= & {} \lambda _0+\lambda _1 f^\mathrm {orth}_1 +\lambda _2 f^\mathrm {orth}_2 + \lambda _3 f^\mathrm {orth}_3+\lambda _4 f^\mathrm {orth}_4 \nonumber \\&+\tau _1 f^\mathrm {tet}_1+\tau _2 f^\mathrm {tet}_2 +\tau _3 f^\mathrm {tet}_3+\tau _4 f^\mathrm {tet}_4.\end{aligned}$$

(4)

Two examples of complete REE patterns that include both $\lambda$ and $\tau$ terms are given in Fig. 3.

Fig. 3

Complete REE patterns showing a an M-type tetrad effect and b a W-type tetrad effect (bottom). Dark blue lines are the smooth patterns derived from $\lambda _n$ terms, light blue lines are the tetrads derived from $\tau _n$ terms, light green lines are the quadratic curves that underlie the tetrad model, and purple lines are the complete REE patterns (i.e., $\Lambda +\mathrm {T}$)

Fitting data to this form is a straightforward procedure in most modern data analysis software packages. There are currently two tools available for fitting a function of this form (Eq. 4) to REE patterns. The first is BLambdaR, an interactive online tool written in the Shiny package for R, using R’s native linear algebra routines (Chang et al. 2020; R Core Team 2020; Anderson et al. 1999; Wilkinson and Rogers 1973; Chambers 1992). It is accessible at https://lambdar.rses.anu.edu.au/blambdar/. The second is pyrolite, a Python package which includes functions for fitting lambdas and tetrads amongst a suite of tools tailored to transformation and visualisation of geochemical data (Williams et al. 2020).

2.4 Treatment of Uncertainties

Naturally, the fitting of often noisy real-world data using $\lambda$ and $\tau$ coefficients is never perfect. Several numerical indicators can be used to assess the fit quality. The first indicator is the coefficient of determination ($r^2$), a number ranging from 0 to 1 which is commonly interpreted as the proportion of variability in the data explained by the model, with higher $r^2$ values indicating lower unexplained variance. Whilst useful for linear fits, $r^2$ has the problem that it increases for every fitted parameter, even if said parameter does not actually improve the fit. As our model may include up to nine parameters (five $\lambda$ and four $\tau$), this may be an issue of concern. Therefore, we use the adjusted $r^2$

$$\begin{aligned} r^2_\mathrm {adj}=1-\left( 1-r^2\right) \times \frac{n-1}{n-p}, \end{aligned}$$

(5)

where n is the number of fitted REE and p is the number of fitted coefficients. Note that the denominator often includes a $-1$ term in the statistical literature to account for an intercept parameter, but in our case this parameter ($\lambda _0$) is included in p, rendering decrementation by one unnecessary. $r^2_\mathrm {adj}$ attempts to account for the potentially spurious increase of $r^2$ when many parameters are fitted, and provide a less biased estimate of unexplained variance (Freund et al. 2010). 

The second indicator, which measures the goodness of fit, is the reduced chi-squared

$$\begin{aligned} \chi ^2_\nu =\frac{\chi ^2}{\nu }=\frac{\sum \left( \ln {[\mathrm {REE}_i]_\mathrm {N}}-\ln \widehat{[\mathrm {REE}_i]_\mathrm {N}}\right) ^2}{\nu \times s(\ln [\mathrm {REE}])*^2},\end{aligned}$$

(6)

where $\nu$ indicates degrees of freedom (i.e., fitted variables minus fitted parameters), the numerator is the sum of squared residuals, the hat symbol indicates expected value from fit, and s indicates the assumed uncertainty on REE measurements as a fraction (i.e., 5% uncertainty should be taken as $s=0.05$). See Bevington and Robinson (2003) and the Appendix in O’Neill (2016) for more details. In general, lower $\chi ^2_\nu$ is better. However, given that we know that there are errors in the data because analytical methods are not foolproof and unavoidable uncertainties exist—$s(\ln [\mathrm {REE}])*^2>0$—it is very unlikely for $\chi ^2_\nu$ to be too low, as that would imply much lower uncertainties in the REE data than is possible according to instrumental and other sources of error. $\chi ^2_\nu$ is then an indicator of this issue: $\chi ^2_\nu =1$ indicates that the fitted residuals of the model are in accordance with the assumed data uncertainty, $\chi ^2_\nu >1$ indicates that the model does not do a good job of fitting or that we have underestimated data uncertainties, and $\chi ^2_\nu <1$ indicates that we are overfitting the model (i.e., it fits a curve with no actual significance) or that we have overestimated data uncertainties. Geochronologists and other geoscientists might be more familiar with an alternative term for $\chi ^2_\nu$, the mean squared weighted deviation, or MSWD (Powell et al. 2002). Individual $\chi ^2_\nu$ values may not be particularly informative themselves, as a single measurement may be exceptionally accurate or inaccurate just by chance. However, when integrated along a large database, all individual $\chi ^2_\nu$ values are themselves $\chi ^2$-distributed with a theoretical maximum that corresponds to the actual s, allowing one to estimate an overall representative uncertainty on REE data across the database (O’Neill 2016).

We also include statistics on individual parameters: the standard error (SE) and p-value. The SE is calculated by taking the square root of the diagonal of the variance-covariance matrix derived from the measured data (Freund et al. 2010; Bevington and Robinson 2003). The ratio of the coefficients and their respective SE is taken, from which a double-sided probability of that value or larger under a t distribution with $\nu$ degrees of freedom is computed as the p-value. Better fits result in lower p-values. In hypothesis testing, p-values smaller than 0.05 are often taken as a rubber stamp for statistical significance, or whether there is a real influence by a certain fitted parameter. We would caution against such a simplistic interpretation in this case, as no hypothesis testing is conducted per se, and the physical interpretation of the shape coefficients is debatable (see below). There is an ongoing discussion on the use and abuse of p-values (Amrhein et al. 2017; McShane et al. 2019; Greenland 2019), and we recommend that they are taken as a semi-quantitative indicator of whether a certain pattern exhibits a certain shape. Values as low as 0.0001 can safely be interpreted as “the shape definitely exists”, and values as high as 0.3 may require more scrutiny. This is particularly important for tetrads: if they cannot be consistently and significantly fitted, it is best to leave them out of the model, as their inclusion breaks $\lambda$ orthogonality, as discussed in the next section.

2.5 Orthogonality and Lack of Thereof During Tetrad Fitting

One principal strength of the $\lambda$ shape parameters and their respective shape components is predicated upon their orthogonality. All components are independent of each other. For instance, bending of a linearly sloping pattern in a fashion consistent with $f_2^\mathrm {orth}$ would increase the magnitude of $\lambda _2$ with no influence on $\lambda _1$. Another striking example, which we encourage readers to attempt on their own using the online tools, is to generate REE data with $\Lambda =(0,0,0,0,-20000)$ in ALambdaR so that the pattern consists essentially of an M-shaped $f_4^\mathrm {orth}$. Next, feed the REE data to BLambdaR, which correctly extracts the initial $\Lambda$ vector (permitting minor deviations due to rounding errors, etc.). Unticking the option to fit $\lambda _4$ leads BLambaR to counterintuitively report $\Lambda =(0,0,0,0)$, as if the fitted pattern is a simple horizontal line. This behaviour provides an excellent demonstration of orthogonality: any $f_n^\mathrm {orth}$ polynomial is essentially invisible to all other orthogonal polynomials, and allows the method to be used consistently across the range of all naturally occurring REE patterns. A certain shape (e.g., quadratic curvature: $\lambda _2$) is intercomparable across all other patterns, regardless of the degree of their slope and sinusoidality, and other shapes described by $\lambda _{\ne 2}$.

The above-described property of orthogonality breaks down when tetrad components ($\tau _nf_n^\mathrm {tet}$) are added to the least-squares fitting. The quadratic $f^\mathrm {tet}_n$ functions are only defined in their respective four-element tetrads and are consequently independent of each other. That is, excluding or including a certain tetrad will not directly affect the other tetrads. However, they are not orthogonal with respect to the $f^\mathrm {orth}_n$ functions, and cannot be made so. The consequence is that when both $\lambda$ and $\tau$ coefficients are sought in a single least-squares fit, the fitted coefficients do not follow the expected mutual orthogonal independence. This can be easily seen for the above pattern ($\lambda _4=-20000$) by ticking all $\tau _n$ while $\lambda _4$ remains unticked. The pattern contains no tetrads and no other $\lambda$-shapes—because that is how we prepared it, yet the fitting leads to the following spurious parameters: $\Lambda =(-0.115, -1.59, -12.1, 666)$ and $\mathrm {T}=(0.57, 0.043, -0.01, 0.15)$. The problem leads to a concern: how can the fitted parameters be trusted and interpreted when attempting to quantify the tetrad effect? An attempt to address these issues is given in the case studies below.

3 Physical Significance of $\lambda$ Shape Coefficients

Accurately measured REE patterns represent actual fractionation processes experienced by the analysed sample. An underlying assumption for least-squares fitting is that real data are well represented by the fitted model. This leads to the question: are the predetermined $f^\mathrm {orth}_n$ shape functions meaningful, and do they directly represent actual natural processes? The physical significance of $\lambda _0$ is trivial, as it represents the overall average chondrite-normalised abundance of REE in a pattern and does not affect the shape.

Fig. 4

Three REE patterns, shifted horizontally. $\Lambda _\mathrm {A}=(1, 0, 300, 0, 0)$, $\Lambda _\mathrm {B}=(1.1, 14.5, 300, 0, 0)$, $\Lambda _\mathrm {C}=(1.5, -18, 300, 0, 0)$

REE patterns essentially reflect REE partitioning and segregation between two or more phases. This partitioning is invariably related to crystal-chemical constraints, and is often modelled using the lattice strain model (Blundy and Wood 2003). Quadratic fits to measured REE partition coefficients in $D-r_i$ space are usually accurate (where D is the concentration ratio between two phases), particularly for minerals that contain a single crystallographic site that can host REE. High-quality partition data often reveal minor asymmetry, which is successfully modelled by an additional cubic $r_i^3$ term. The physical significance of $\lambda _1$ and $\lambda _2$ is therefore a proxy for simple lattice strain partitioning, as they similarly describe quadratic curves in ratio $-r_i$ space. The parabola defined by $\lambda _2f^\mathrm {orth}_2$, on its own, has an extremum close to Eu (pattern A in Fig. 4), and it can be shifted horizontally by varying $\lambda _0$ and $\lambda _1$ (as any parabola defined by $y(x)=ax^2+bx+c$ can be transformed in $x-y$ space by varying b and c), demonstrated by patterns B and C in Fig. 4. However, this is a simplification because the hypothetical patterns in Fig. 4 would have formed by fractionation of a single-REE-site mineral from an unfractionated CI chondrite pattern. Real REE patterns are unlikely to be that simple, because they are formed by fractionation of numerous minerals, each having more than one REE-hosting crystallographical site, from an already-fractionated source. In essence, a currently measured REE pattern reflects the entire lattice strain partitioning pedigree since the formation of the solar system. Modelling this history is theoretically possible, but arduous and most likely non-unique, leading to accurate reconstruction being unachievable. Nevertheless, fractionated natural patterns will inevitably be smoothly curved, leading to the use of $\lambda _3$ and $\lambda _4$. There is arguably no direct physical significance to the cubic $f^\mathrm {orth}_3$ and quartic $f^\mathrm {orth}_4$, but these higher-order polynomials represent fine-tuning and approximate small variations that stem from the cumulative parabolas, of varying heights, widths, and positions, encapsulated in a single REE pattern. The lack of physical significance does not, however, detract from the usefulness of the method, as it does an excellent job of describing various shape features of REE patterns. The method also has practical applications in petrogenetic modelling by parameterising REE crystal/melt partition coefficients as $D=\delta _0+\delta _1f_1^\mathrm {orth}+\cdots$ and using these coefficients in petrogenetic process vectors ($\Psi$) and coefficients ($\psi _n$). See O’Neill (2016) for more details and examples.

4 Quantification of Cerium and Europium Anomalies

A common feature of REE patterns is the presence of anomalies: the deviation of the redox-sensitive elements Ce and Eu from the smooth pattern exhibited by all other elements. The presence, sign, and magnitude of Ce and Eu anomalies are useful indicators of various petrogenetic processes. For example, Ce and Eu anomalies in zircon are indicative of the magma oxidation state from which the zircon crystallised (Loucks et al. 2020 and references therein). Anomalies are typically quantified by taking the ratio of the observed normalised concentration of Ce and Eu to a calculated normalised concentration derived by the geometric mean of the two neighbouring elements, often termed Ce* and Eu* (geometrically being the linear interpolation between the two elements on a logarithmic plot): $\mathrm {Ce/Ce^*_{LI}=\sqrt{La_N\times Pr_N}}$ and $\mathrm {Eu/Eu^*_{LI}=\sqrt{Sm_N\times Gd_N}}$, where subscript LI indicates linear interpolation. A major disadvantage of this method is the assumption of linearity, namely the question of whether the geometric mean indeed represents the true concentration of Ce and Eu had there been no redox effects. Whilst this method may give acceptable results for relatively straight patterns, it rapidly becomes inaccurate when patterns are strongly curved. For example, the REE pattern shown in Fig. 5 is a typical zircon-like pattern, with positive Ce and negative Eu anomalies. Ce anomalies in zircon are a useful petrogenetic indicator (Burnham and Berry 2014; Trail et al. 2015), and many approaches to zircon Ce anomaly quantification have been published in recent decades (Zhong et al. 2019). Visually, the difference between the interpolated line and the non-anomalous pattern looks negligible—the difference is smaller than the symbol size and the line plots on top of the circles. However, this is misleading, as the y-axis covers six orders of magnitude. We indicate the expected value of Ce from the smooth curve constructed from shape components as $\mathrm {Ce^*_{\lambda \tau }}$. The pattern has been generated with a $\mathrm {Ce/Ce^*_{\lambda \tau }}=40$, but calculating the anomaly using the interpolated line results in $\mathrm {Ce/Ce^*_{LI}}=47.5$, an overestimation of about 19%. Likewise, $\mathrm {Eu/Eu^*_{\lambda \tau }}=0.6$ (the true Eu anomaly), whereas the calculated anomaly via interpolation is $\mathrm {Eu/Eu^*_{LI}}=0.667$ , an underestimation of about 11%.

Fig. 5

A zircon-like REE pattern generated using the parameters $\Lambda =(4.6, -63, -545, -1400, \hbox {10500})$ in light blue, and the same pattern with Ce and Eu anomalies in dark blue. Lightest blue are linear lines connecting La and Pr, and Sm and Gd. The pattern is truncated at Tb for clarity, and the full pattern can be obtained by inputting the $\Lambda$ vector into ALambdaR

Following the above example, more accurate anomaly quantifications are conducted by taking the ratio of measured Ce and Eu to that value expected from the least-squares fit, instead of linear interpolation. Clearly, this requires that anomalous Ce and Eu are not used for the actual fitting. It should be noted that for this method to succeed, accurate values for all REE are crucial, particularly for elements at the extreme edges (i.e., La, Yb, and Lu) when calculating a Ce anomaly. In some materials, this proves difficult, with Lu having the lowest chondritic abundance of all the REE. Additionally, La is notoriously difficult to measure in minerals that reject LREE, most notably zircon (Zou et al. 2019; Burnham 2020). In such cases, we recommend a visual inspection of the measured pattern, fitted line, and statistical indicators as given by pyrolite and interactively by BLambdaR, to assess whether an anomaly determination is reliable.

5 Case Studies

5.1 Gadolinium Anomaly and Tetrad Effect in Seawater

REE compositions and patterns of seawater are a useful proxy for studying various ocean processes such as mixing of reservoirs, sediment load, hydrothermal vent input, etc. Our understanding of these processes might be improved if shape coefficients are used to describe seawater. As an example, we took the representative seawater data “BATS” of van de Flierdt et al. (2012) and fitted $\lambda$ coefficients to it, treating Ce and Eu as anomalous. The match between the fitted line and the data is close, but not perfect (adjusted $r^2$ is 0.95 and 0.90 for shallow and deep waters, respectively, Fig 6a, Online Resource A1). Adding $\tau$ coefficients markedly improved the fit, with adjusted $r^2$ increasing to 0.997 and 0.985, respectively (Fig 6b, Online Resource A2).

Fig. 6

Seawater data (circles) and least-squares fit (lines). a Lines fit using $\lambda _n$ coefficients. b Lines fit using $\lambda _n$ and $\tau _n$ coefficients. Data from van de Flierdt et al. (2012)

It may be argued that the improved fit is least-squares trickery, because increasing the amount of parameters invariably leads to better fits and greater $r^2$. However, there are several lines of evidence to suggest the effect is real. First, we are using the adjusted $r^2$ which accounts for spurious increases in determination caused by additional parameters. Second, all $\tau$ coefficients are consistently negative (i.e., reflecting a W-type tetrad effect) and showing roughly the same magnitude, except $\tau _1$ which is known to be stronger in geological environments (Fig. 7). Had the additional $\tau$ coefficients described noise and spuriously improved the least-squares fit, this consistency would have been unlikely. Third, the uncertainty range of the $\tau$ coefficients is well away from zero (Fig. 7), and the $\tau$ significance p-values for deep water—consisting of more abundant REE, hence more likely to be accurately and precisely measured—are all under 0.012. For less accurate shallow water, the $\tau _1$ p-value is 0.068, whereas the others are under 0.04.

Fig. 7

Tetrad coefficients for seawater (BATS) and PAAS. Uncertainty shown is 1 standard error

Whether seawater has a tetrad effect or not has been previously debated (McLennan 1994; Kawabe et al. 1998), but our modelling shows that, when reliable REE data are normalised using reliable CI values, a significant tetrad effect is indeed observed, in agreement with recent observations (Ernst and Bau 2021). Seawater is commonly normalised to PAAS (Post-Archean Australian shale) rather than CI. Fitting $\tau$ coefficients to accurately measured PAAS data of Pourmand et al. (2012) reveals that PAAS does not contain any statistically significant tetrad effect, with p-values ranging from 0.24 to 0.96 (Fig. 7, Online Resource A3). Therefore, the choice of either CI or PAAS for normalisation is unlikely to generate, obscure, or erase tetrad effects. The tetrad fitting also shows that the minor positive Gd anomaly observed in Fig. 6 is the result of the tetrad effect and represents a natural phenomenon. A positive Gd anomaly in seawater, lakes, or the rivers that feed them is interpreted as anthropogenic contamination, and the anomaly is usually quantified using deviations of Gd from linear interpolation of other REE on a PAAS-normalised plot (Kulaksız and Bau 2011a), or by using non-orthogonal linear or polynomial least squares (Kulaksız and Bau 2011b; Müller et al. 2002). As these methods do not take overall pattern curvature or tetrad effects into consideration, subtle anomalies which derive from the tetrad effect might be erroneously ascribed to anthropogenic contamination. We recommend that least-squares fits to REE patterns using $\lambda$ and $\tau$  parameters be examined to check for consistency with the linear interpolation or fitting methods.

5.2 Tetrad Effects in the Toongi Rare Metal Deposit

The Toongi Deposit in NSW, Australia, hosts substantial resources of Zr, Hf, Nb, Ta, REE, and Y. The mineralisation is hosted in eudialyte and various later alteration products such as bastnäsite, vlasovite, milarite, catapleiite, and gaidonnayite (Spandler and Morris 2016). These minerals often exhibit unusual and complex REE patterns, and Spandler and Morris (2016) qualitatively noted the presence of the tetrad effect in the patterns they obtained from the Toongi minerals.

Fig. 8

Tetrad $\tau$ coefficients for Toongi minerals. a to d show the first to fourth tetrad coefficients, respectively. Values fitted with $\lambda _4$ are shown in light green, and values fitted with $\lambda _4$ omitted are shown in dark blue. ca, catapleiite; ga, gaidonnayite

Figure 8 shows $\tau$ coefficients derived from the Spandler and Morris (2016) data (see Online Resource A4 and A5). As discussed above, there is a clash between $\tau$ and $\lambda$ shape components, particularly when the patterns are highly curved at the edges, in which case the same features are fitted by both $\tau _1$, $\tau _4$, and $\lambda _4$. To investigate the issue, we fitted the data with $\lambda _4$ (light green) and without $\lambda _4$ (dark blue). It is immediately evident that eudialyte has small positive tetrads, which are negligible compared to the later minerals. This is consistent with the inference of Spandler and Morris (2016) that it is a primary mineral, unaffected by lower-temperature hydrothermal processes. For minerals other than eudialyte, there is a clear difference between the $\lambda _4$ inclusive and exclusive values of $\tau _1$and $\tau _4$. Which values are correct? The two middle tetrads—$\tau _2$ and $\tau _3$—are informative as they remain mostly unchanged irrespective of whether $\lambda _4$ is fitted or not. Therefore, we expect the external tetrads to follow the same behaviour, which is only achieved when $\lambda _4$ is not fitted. In contrast, $\tau _1$ and $\tau _4$ are diminished in magnitude and reversed in sign when $\lambda _4$ is included, a result inconsistent with the middle tetrads. Therefore, in this case of strongly curved LREE and HREE, we proffer that $\lambda _4$ is spuriously preventing the correct identification of the tetrads.

The LREE tetrads, $\tau _1$ and $\tau _2$, are strongly positive in all secondary minerals. One possibility is that LREE in the Toongi secondary assemblage are dominated by bastnäsite and other REE carbonates. Although Spandler and Morris (2016) do not provide full REE data for these minerals, partial patterns based on La, Ce, Pr, Nd, and Sm suggest the presence of a negative $\tau _1$ for REE carbonates. The secondary minerals shown in Fig. 8 are LREE-poor relative to the MREE and HREE, and their positive $\tau _1$ and $\tau _2$ may complement the LREE-rich carbonates. The HREE tetrad behaviour is less clear. The HREE are likely dominated by vlasovite, and its positive $\tau _3$ and $\tau _4$ complement the slightly negative $\tau _3$ and $\tau _4$ observed in milarite.

The above discussion is not meant as a comprehensive study of REE partitioning and the tetrad effect in Toongi, but serves as an example of a range of possibilities made available by the use of quantitative $\lambda$ and $\tau$ determinations. It also shows an example in which it is best to exclude $\lambda _4$ from the least-squares fitting.

5.3 $\lambda _4$ and Ce Anomalies in the Troodos Ophiolite

In our final case study, we examine REE patterns from oceanic plagiogranites collected from the Troodos ophiolite in Cyprus by Anenburg et al. (2015). A representative pattern is shown in Fig. 9. All $\tau$ coefficients are highly insignificant, regardless of whether $\lambda _4$ is included in the fit, and therefore there is patently no tetrad effect exhibited by the pattern (Online Resources A6 and A7), although by simple visual inspection of the pattern, it might appear as if there is a tetrad effect. This pattern is remarkable by how well it is described by $\lambda _4$ (Online Resource A8). Its magnitude is rather large (>10,000) compared to some of the lower-order coefficients ($\lambda _2\approx 4$ and $\lambda _3\approx 200$). $\lambda _4$ is often less precisely determined than the lower-order coefficients (O’Neill 2016), but in this case $\mathrm {SE}(\lambda _4)<20\%$, whereas $\mathrm {SE}(\lambda _2)>40\%$ and $\mathrm {SE}(\lambda _3)>20\%$.

Fig. 9

REE data for a Troodos ophiolite plagiogranite (dark blue dots), $\lambda$ polynomial fit to the data, excluding Eu and Ce (dark blue line), and the $\lambda _4$ component of the fit (light blue line)

The peculiar dominance of $\lambda _4$ can be explained as a combination of the three most REE-rich minerals that occur in the plagiogranites: hornblende, zircon, and epidote–allanite. REE patterns for the Troodos hornblende are characterised by two positively straight segments: a steep segment from La to Sm, and an almost horizontal segment from Gd to Lu (Gillis 2002), such that it would be well described by $\lambda _1<0$ and $\lambda _2<0$. Epidotes and allanites are LREE-rich (Anenburg et al. 2015), and zircon is HREE-rich (Chen et al. 2020), such that small amounts of these accessory minerals allow the lightest LREE and heaviest HREE to poke out of the hornblende-dominated pattern. This serves as an example of the lack of distinct physical significance for $\lambda _4$. The central wide ridge records the maximum of hornblende curvature, and the two steep peripheral segments record allanite and zircon (Fig. 9). However, there is no reason to assume that $\lambda _4$ is particularly efficient in describing hornblende, allanite, or zircon patterns. It is just coincidental that the combinations of three lower-order $f^\mathrm {orth}_n$ and their lateral translations were accurately described by $\lambda _4$.

The presumably coincidental yet excellent $\lambda _4$ fit raises an interesting question. What is the nature of the Ce anomaly in the pattern shown in Fig. 9? Is it a negative anomaly, as Ce plots below the La–Pr interpolation line ($\mathrm {Ce/Ce^*_{LI}}<1$), or is it positive, as Ce plots above the fitted line ($\mathrm {Ce/Ce^*_{\lambda \tau }}>1$)? The REE patterns of a rock sample reflect the REE patterns of its constituent minerals. The plagiogranites described by Anenburg et al. (2015) are LREE-poor, with Ce concentrations of around 5 ppm. This low Ce content is very sensitive to Ce-rich materials, such as zircon. Plagiogranite-hosted zircon from Troodos is exceptionally Ce-rich, with Ce/Ce* often reaching values greater than 100 (Chen et al. 2020; Morag et al. 2020), and as such is expected to contribute substantial Ce to the whole rock budget (along with trace contributions from epidote and allanite). This anomalous contribution of Ce would not be part of a smooth curve, and propagate the anomaly to the whole rock, under the assumption that Ce-rich zircon did not fractionate earlier and impart a measurable anomaly to the melt. This case study shows—in addition to the examples shown earlier—that Ce anomalies obtained using the deviation from the polynomial fit are more accurate where REE patterns exhibit curvature and better reflect geological and mineralogical observations relative to the deviation from the interpolation line. In cases of high curvature such as the one presented here, the inferred anomaly even has the opposite sign.

6 Concluding Remarks

We present a method for parameterising REE patterns, combining the orthogonal polynomials of O’Neill (2016) with adapted tetrad polynomials based on Minami and Masuda (1997). Our method provides uncertainty estimates and is implemented in two forms: an interactive online app written in R and Shiny (BLambdaR), and a Python package that includes our $\lambda$ and $\tau$ fitting algorithm with other geochemical data analysis tools (pyrolite). The case studies show that there is no one single fitting recipe that applies universally, and we recommend that different fitting options be explored. The case studies also show some of the considerations that apply when treating real-world data and how they might be approached.

Appendix

1.1 Orthogonal Polynomial Root Derivation

The process for deriving the orthogonal polynomials (Eq. 2) is explained in Bevington and Robinson (2003). $\beta$ is calculated by solving ${\scriptstyle \sum }(r_i-\beta )=0$, leading to $\beta =\frac{1}{N}{\scriptstyle \sum }{r_i}$, which is the mean of all ionic radii, where $N=14$ (i.e., the number of REE).

$\gamma _n$, $\delta _n$ and $\epsilon _n$ require solving of equation systems of increasing complexities (Bevington and Robinson 2003)

$$\begin{aligned}&\begin{array}{l} {\scriptstyle \sum }{(r_i-\gamma _1)(r_i-\gamma _2)}=0\\ {\scriptstyle \sum }{r_i(r_i-\gamma _1)(r_i-\gamma _2)}=0,\\[6mm] \end{array}\\&\begin{array}{l} {\scriptstyle \sum }{(r_i-\delta _1)(r_i-\delta _2)(r_i-\delta _3)}=0\\ {\scriptstyle \sum }{r_i(r_i-\delta _1)(r_i-\delta _2)(r_i-\delta _3)}=0\\ {\scriptstyle \sum }{r_i^2(r_i-\delta _1)(r_i-\delta _2)(r_i-\delta _3)}=0,\\[6mm] \end{array}\\&\begin{array}{l} {\scriptstyle \sum }{(r_i-\epsilon _1)(r_i-\epsilon _2)(r_i-\epsilon _3)(r_i-\epsilon _4)}=0\\ {\scriptstyle \sum }{r_i(r_i-\epsilon _1)(r_i-\epsilon _2)(r_i-\epsilon _3)(r_i-\epsilon _4)}=0\\ {\scriptstyle \sum }{r_i^2(r_i-\epsilon _1)(r_i-\epsilon _2)(r_i-\epsilon _3)(r_i-\epsilon _4)}=0\\ {\scriptstyle \sum }{r_i^3(r_i-\epsilon _1)(r_i-\epsilon _2)(r_i-\epsilon _3)(r_i-\epsilon _4)}=0. \end{array} \end{aligned}$$

The solutions in Bevington and Robinson (2003) only work for equispaced $r_i$, but that is not the case here: the difference in ionic radii is not constant (e.g., the difference between La and Ce is greater than Yb and Lu). O’Neill (2016) solved these non-linear systems for $\gamma _n$, $\delta _n$, and $\epsilon _n$ using Excel’s “Solver” function that numerically finds the values that satisfy the equations. However, an exact analytical solution is achieved by application of Vieta’s formula, which allows the conversion of the above systems into simple polynomial equations whose roots are the sought-after parameters. Starting with the example of $\gamma _n$, the systems are constructed by rewriting the above equations in a different form

$$\begin{aligned}\begin{aligned}&{\scriptstyle \sum }{r_i^2}-{\scriptstyle \sum }{r_i(\gamma _1+\gamma _2)}+{\scriptstyle \sum }{(\gamma _1\gamma _2)}=0\\&{\scriptstyle \sum }{r_i^3}-{\scriptstyle \sum }{r_i^2(\gamma _1+\gamma _2)}+{\scriptstyle \sum }{r_i(\gamma _1\gamma _2)}=0. \end{aligned} \end{aligned}$$

By defining $\gamma _1+\gamma _2=a$ and $\gamma _1\gamma _2=b$ and noting that ${\scriptstyle \sum }{r_i^n}$ are constants, the systems simplify to the following problem, written in matrix form

$$\begin{aligned} \begin{pmatrix} -{\scriptstyle \sum }{r_i} &{} N \\ -{\scriptstyle \sum }{r_i^2} &{} {\scriptstyle \sum }{r_i} \end{pmatrix} \begin{pmatrix} a\\ b \end{pmatrix} = \begin{pmatrix} -{\scriptstyle \sum }{r_i^2}\\ -{\scriptstyle \sum }{r_i^3} \end{pmatrix}. \end{aligned}$$

Once solved, a and b are then used as coefficients in the following quadratic equation

$$\begin{aligned} \gamma ^{2}-a\gamma +b=0, \end{aligned}$$

the roots of which are the two parameters $\gamma _1$ and $\gamma _2$.

Similarly, the system for $\delta _n$ is rewritten as

$$\begin{aligned}\begin{aligned}&{\scriptstyle \sum }{r_i^3}-{\scriptstyle \sum }{r_i^2(\delta _1+\delta _2+\delta _3)}+{\scriptstyle \sum }{r_i(\delta _1\delta _2+\delta _1\delta _3+\delta _2\delta _3)}-{\scriptstyle \sum }{(\delta _1\delta _2\delta _3)}=0\\&{\scriptstyle \sum }{r_i^4}-{\scriptstyle \sum }{r_i^3(\delta _1+\delta _2+\delta _3)}+{\scriptstyle \sum }{r_i^2(\delta _1\delta _2+\delta _1\delta _3+\delta _2\delta _3)}-{\scriptstyle \sum }{r_i(\delta _1\delta _2\delta _3)}=0\\&{\scriptstyle \sum }{r_i^5}-{\scriptstyle \sum }{r_i^4(\delta _1+\delta _2+\delta _3)}+{\scriptstyle \sum }{r_i^3(\delta _1\delta _2+\delta _1\delta _3+\delta _2\delta _3)}-{\scriptstyle \sum }{r_i^2(\delta _1\delta _2\delta _3)}=0. \end{aligned} \end{aligned}$$

Define $\delta _1+\delta _2+\delta _3=a$, $\delta _1\delta _2+\delta _1\delta _3+\delta _2\delta _3=b$, and $\delta _1\delta _2\delta _3=c$, and solve the linear system

$$\begin{aligned} \begin{pmatrix} -{\scriptstyle \sum }{r_i^2} &{}\quad {\scriptstyle \sum }{r_i} &{}\quad -N \\ -{\scriptstyle \sum }{r_i^3} &{}\quad {\scriptstyle \sum }{r_i^2} &{}\quad -{\scriptstyle \sum }{r_i}\\ -{\scriptstyle \sum }{r_i^4} &{}\quad {\scriptstyle \sum }{r_i^3} &{}\quad -{\scriptstyle \sum }{r_i^2} \end{pmatrix} \begin{pmatrix} a\\ b\\ c \end{pmatrix} = \begin{pmatrix} -{\scriptstyle \sum }{r_i^3}\\ -{\scriptstyle \sum }{r_i^4}\\ -{\scriptstyle \sum }{r_i^5} \end{pmatrix}. \end{aligned}$$

a, b, and c are then used in the cubic equation

$$\begin{aligned} \delta ^3-a\delta ^2+b\delta -c=0 \end{aligned}$$

to solve for the roots $\delta _1$, $\delta _2$, and $\delta _3$.

Finally, the system for $\epsilon _n$ is rewritten as

$$\begin{aligned} \begin{aligned}&{\scriptstyle \sum }{r_i^4}-{\scriptstyle \sum }{r_i^3(\epsilon _1+\epsilon _2+\epsilon _3+\epsilon _4)}+{\scriptstyle \sum }{r_i^2(\epsilon _1\epsilon _2+\epsilon _1\epsilon _3+\epsilon _1\epsilon _4+\epsilon _2\epsilon _3+\epsilon _2\epsilon _4+\epsilon _3\epsilon _4)}\\ {}&\quad - {\scriptstyle \sum }{r_i(\epsilon _1\epsilon _2\epsilon _3+\epsilon _1\epsilon _2\epsilon _4+\epsilon _1\epsilon _3\epsilon _4+\epsilon _2\epsilon _3\epsilon _4)}+{\scriptstyle \sum }{(\epsilon _1\epsilon _2\epsilon _3\epsilon _4)}=0\\&{\scriptstyle \sum }{r_i^5}-{\scriptstyle \sum }{r_i^4(\epsilon _1+\epsilon _2+\epsilon _3+\epsilon _4)}+{\scriptstyle \sum }{r_i^3(\epsilon _1\epsilon _2+\epsilon _1\epsilon _3+\epsilon _1\epsilon _4+\epsilon _2\epsilon _3+\epsilon _2\epsilon _4+\epsilon _3\epsilon _4)}\\ {}&\quad -{\scriptstyle \sum }{r_i^2(\epsilon _1\epsilon _2\epsilon _3+\epsilon _1\epsilon _2\epsilon _4+\epsilon _1\epsilon _3\epsilon _4+\epsilon _2\epsilon _3\epsilon _4)}+{\scriptstyle \sum }{r_i(\epsilon _1\epsilon _2\epsilon _3\epsilon _4)}=0\\&{\scriptstyle \sum }{r_i^6}-{\scriptstyle \sum }{r_i^5(\epsilon _1+\epsilon _2+\epsilon _3+\epsilon _4)}+{\scriptstyle \sum }{r_i^4(\epsilon _1\epsilon _2+\epsilon _1\epsilon _3+\epsilon _1\epsilon _4+\epsilon _2\epsilon _3+\epsilon _2\epsilon _4+\epsilon _3\epsilon _4)}\\ {}&\quad -{\scriptstyle \sum }{r_i^3(\epsilon _1\epsilon _2\epsilon _3+\epsilon _1\epsilon _2\epsilon _4+\epsilon _1\epsilon _3\epsilon _4+\epsilon _2\epsilon _3\epsilon _4)}+{\scriptstyle \sum }{r_i^2(\epsilon _1\epsilon _2\epsilon _3\epsilon _4)}=0\\&{\scriptstyle \sum }{r_i^7}-{\scriptstyle \sum }{r_i^6(\epsilon _1+\epsilon _2+\epsilon _3+\epsilon _4)}+{\scriptstyle \sum }{r_i^5(\epsilon _1\epsilon _2+\epsilon _1\epsilon _3+\epsilon _1\epsilon _4+\epsilon _2\epsilon _3+\epsilon _2\epsilon _4+\epsilon _3\epsilon _4)}\\ {}&\quad -{\scriptstyle \sum }{r_i^4(\epsilon _1\epsilon _2\epsilon _3+\epsilon _1\epsilon _2\epsilon _4+\epsilon _1\epsilon _3\epsilon _4+\epsilon _2\epsilon _3\epsilon _4)}+{\scriptstyle \sum }{r_i^3(\epsilon _1\epsilon _2\epsilon _3\epsilon _4)}=0. \end{aligned}\end{aligned}$$

 Define $\epsilon _1+\epsilon _2+\epsilon _3+\epsilon _4=a$, $\epsilon _1\epsilon _2+\epsilon _1\epsilon _3+\epsilon _1\epsilon _4+\epsilon _2\epsilon _3+\epsilon _2\epsilon _4+\epsilon _3\epsilon _4=b$, $\epsilon _1\epsilon _2\epsilon _3+\epsilon _1\epsilon _2\epsilon _4+\epsilon _1\epsilon _3\epsilon _4+\epsilon _2\epsilon _3\epsilon _4=c$, and $\epsilon _1\epsilon _2\epsilon _3\epsilon _4=d$, and solve the linear system

$$\begin{aligned} \begin{pmatrix} -{\scriptstyle \sum }{r_i^3} &{}\quad {\scriptstyle \sum }{r_i^2} &{}\quad -{\scriptstyle \sum }{r_i} &{}\quad N\\ -{\scriptstyle \sum }{r_i^4} &{}\quad {\scriptstyle \sum }{r_i^3} &{}\quad -{\scriptstyle \sum }{r_i^2} &{}\quad {\scriptstyle \sum }{r_i}\\ -{\scriptstyle \sum }{r_i^5} &{}\quad {\scriptstyle \sum }{r_i^4} &{}\quad -{\scriptstyle \sum }{r_i^3} &{}\quad {\scriptstyle \sum }{r_i^2}\\ -{\scriptstyle \sum }{r_i^6} &{}\quad {\scriptstyle \sum }{r_i^5} &{}\quad -{\scriptstyle \sum }{r_i^4} &{}\quad {\scriptstyle \sum }{r_i^3} \end{pmatrix} \begin{pmatrix} a\\ b\\ c\\ d \end{pmatrix} = \begin{pmatrix} -{\scriptstyle \sum }{r_i^4}\\ -{\scriptstyle \sum }{r_i^5}\\ -{\scriptstyle \sum }{r_i^6}\\ -{\scriptstyle \sum }{r_i^7} \end{pmatrix}. \end{aligned}$$

a, b, c, and d are used in the quartic equation

$$\begin{aligned} \epsilon ^4-a\epsilon ^3+b\epsilon ^2-c\epsilon +d=0, \end{aligned}$$

which is solved for the roots $\epsilon _1$, $\epsilon _2$, $\epsilon _3$, and $\epsilon _4$.