Erratum: ‘Representation of complex seismic sources by orthogonal moment–tensor fields’

In a recent GJI paper (Jordan & Juarez 2019; abbreviated JJ19), we developed a theoretical framework for representing complex seismic sources in terms of orthogonal moment–tensor fields. The purpose of this note is to correct a serious error in that theory discovered by the authors after publication. Specifically, the stress-glut representation theorem stated in section 2 of JJ19 is false because the tensor factorization in JJ19 eq. (18) is not exact for arbitrary stress gluts.


D E S C R I P T I O N O F T H E E R RO R
We adhere to definitions and notation in JJ19 with slight variations to handle the more general representations required by the corrected theory. In the interest of brevity, we do not repeat JJ19 references to previous work unless necessary for technical clarity.
The stress-glut density˙ describing an arbitrary seismic source is a symmetric 2nd-order tensor field defined over a finite space-time volume, x ∼ (r, t) ∈ V. Each element of this field˙ (x) is a member of S (2) 3 , the 6-D vector space of symmetric, real-valued, 2nd-order tensors of dimension 3 (Comon et al. 2008). The local source mechanism is the normalized stress-glut density,ˆ (x) ≡˙ (x)/ ˙ (x) , a unit tensor in S (2) 3 . The rank of˙ is the dimension D ≤ 6 of the subspace of S (2) 3 spanned by the measurable stress-glut elements. The field can always be represented as a sum over an orthonormal basis comprising D constant source-mechanism tensors: The pth polynomial moment of the stress-glut density˙ (x) is an ( p + 2)-order tensor linearly related to pth-order moments the of the scalar fields g α (x) by Here (x − x * ) p is the symmetric outer product tensor with Cartesian components (x k − x * k )(x l − x * l )(x m − x * m ) · · · (p factors). We assume the expansion point x * is located in V hull , the smallest convex set that contains V. For a source of maximum spatial dimension L and total duration T, time units are scaled to space units using the velocity V ≡ L/T , so that x − x * ≤ L for all x, x * ∈ V hull .
In JJ19, we considered expansions in which the leading term in (1a) is the seismic moment tensor, defined to be the zeroth moment of the stress glut, M 0 is the Aki seismic moment, andM 0 ∈ S 3 is the average source mechanism. The zeroth-degree approximation to the stress glut is given by its projection onto this average mechanism, (x) ≈M 0 g 0 (x) . (4) The scalar field g 0 (x) ≡M 0 :˙ (x) is the zeroth-degree space-time function.
Eq. (4) is widely used in seismology as a finite-source approximation. In JJ19, we sought better approximations for sources with local mechanismsˆ (x) that are variable in space and time, which we called 'mechanism complexity'. One measure of mechanism complexity is the difference between the Aki moment M 0 and the total seismic moment M T , defined in JJ19 to be the volumetric integral of the scalar moment density, The 'extra' seismic moment M T − M 0 > 0 can be a substantial fraction of the total moment for sources with high mechanism complexity; for example about 20 per cent in the case of the 2016 Kaikōura earthquake (M 7.8).
The goal of JJ19 was to improve the approximation (4) by representing complex seismic sources as sums over orthogonal momenttensor fields, whereM 0 is given by eq. (3), and the higher-degree source mechanisms are derived sequentially from the lowest-order moment of the residual fields, To describe such a 'moment-oriented representation' (MOR), we used the concept of shift-invariance. A polynomial moment of a tensor field is shift-invariant if its value does not depend on the expansion point x * . For any p ≥ 1, μ ( p) (g) is shift-invariant if and only if μ (q) (g) = 0 for all q < p. An example is μ (1) ( ˙ 1 ), the first polynomial moment of the first-degree residual field, which is shift-invariant because μ (0) ( ˙ 1 ) is, by definition, zero. In JJ19 (eq. 18), we asserted that μ (α) ( ˙ α ) can be constructed to be shift-invariant and to be uniquely factorable into the outer product of a 2nd-order mechanism tensorM α and an α-order, shift-invariant tensor J (α) ≡ μ (α) (g α ). These assertions are incorrect.
To see why, we express the α th moment of the α th residual field in its Cartesian components, This tensor, of order α + 2, is twofold symmetric in the indices i, j, which run over the 3 space dimensions, and α-fold symmetric in the indices klm · · · , which run over the 4 space-time dimensions. Therefore, μ (α) 4 admits a canonical polyadic decomposition (CPD) of the form where is an ordered set of non-negative singular values (Kolda & Bader 2009). As the terminology indicates, CPD is the generalization of matrix singular value decomposition (SVD) to tensors (De Lathauwer et al. 2000). The orthonormality conditions are Here, as in JJ19, the operator • signifies contraction over all α indices. The dimension of S (α) N , the vector space of symmetric tensors of dimension N and order α is N + α − 1 α (Comon et al. 2008).
In the special case D = 2, there is only one non-zero singular value, λ 1 and the JJ19 factorization works by settingM 1 to be U (1) 1 . Hence, section 6 of JJ19, which discusses planar fault models with a single degree of rotational freedom (D = 2), stands without correction.

C O R R E C T I O N T O T H E M O M E N T -O R I E N T E D R E P R E S E N TAT I O N
The basis for an orthogonal moment-tensor representation is an orthonormal set of up to six source mechanisms. There are various ways to construct a basis set when the stress glut is known. Multilinear principal component analysis (PCA, De Lathauwer et al. 2000;Jolliffe & Cadima 2016) provides the most efficient expansion in a least-squares sense, because each mechanism added to the basis set successively minimizes the Euclidean norm of the stress-glut residual. We distinguish the results of PCA from the MOR by using a different notation for the basis elements and scalar fields (cf. eq. 1): In PCA,N α is chosen to maximize the projection of the resid- onto the basis element, which is equivalent to minimizing the squared norm of δ˙ α+1 (x), Here δ˙ α+1 The members of this basis set can be computed as the 2nd-order symmetric eigentensors of the outer-product integral, The tensor G is 4th-order and has the elastic-tensor symmetries, and its terms can be ordered by their non-zero eigenvalues, ν 2 0 ≥ ν 2 1 ≥ · · · ≥ ν 2 D−1 > 0, which are the squared norms of the PCA scalar fields, ν 2 α = V h 2 α (x) dx. In particular, the zeroth-degree term in the PCA expansion maximizes ν 2 0 over all choices ofN 0 . This sizemeasure is not necessarily equal to the square of the Aki moment, which instead maximizes V g 0 (x) dx over all choices ofM 0 . Hence, the directionN 0 in S (2) 3 is not necessarily aligned with the seismic moment tensorM 0 , and, in general, h 0 (x) = g 0 (x) and δ˙ 1 (x) = ˙ 1 (x). The PCA basis set of orthogonal mechanisms is an incomplete parametrization of the source model; to set up the inverse problem, we must also parametrize the scalar fields defined by the basis set. In a multifault rupture, these scalar fields will be non-zero only on an unknown (and probably fractal) set of rupture surfaces embedded in the source space-time volume-the 'support set' of the stress glut. In some cases, the spatial geometry of a multifault rupture is known well enough from local geological and geophysical data to admit a plausible discretization of the spatial support. Hamling et al. (2017) and Wang et al. (2018), for example, derived source models for the 2016 Kaikōura earthquake using an inferred geometry of up to 12 fault segments of various sizes and orientations. Situations where the rupture geometry is well constrained are untypical, however.
Owing to the uncertainty of the support set, we seek to parametrize the stress glut in terms of its low-order moments μ ( p) (˙ ). These moment tensors are integrals over the source volume and thus insensitive to the geometrical details of the rupture geometry. Moreover, they provide a natural parametrization for the teleseismic study of source complexity, because the low-order (small p) terms dominate the low-frequency source radiation. Projecting Downloaded from https://academic.oup.com/gji/article-abstract/222/2/1333/5856841 by guest on 14 June 2020 these tensors onto the basis set in eq. (1) generates the moment array μ ( p) (g α ). When considering this moment array as the collection of source parameters, we abbreviate its subdiagonal elements ( p < α) by ε ( p) α , its diagonal elements ( p = α) by J (α) , and those above the diagonal ( p > α) by μ ( p) α . Table 1 displays this moment-array notation for D = 4, 0 ≤ p ≤ 3, a sufficiently general case for our discussion.
The first basis element of an MOR (α = 0) is defined by the source monopole ( p = 0) or centroid moment tensor (CMT). The CMT maximizes the Aki seismic moment M 0 = J (0) / √ 2 and concentrates this moment at the source centroid, x 0 = μ (1) (˙ )/J (0) . We expand the stress-glut polynomial moments about this centroid, which we take to be the coordinate origin, so that μ (1) 0 = 0. The normalized second moment μ (2) 0 /J (0) then specifies the characteristic space-time ellipsoid of g 0 ; twice the square-root of its largest eigenvalue is the characteristic length of the zeroth-degree field, L 0 , which we presume to be comparable to the outer source dimension L. The first-order residual field, obtained by subtracting the CMT projection from the stress-glut density, has no monopole moment: μ (0) ( ˙ 1 ) = 0. All zeroth-order moments of higher degree therefore vanish; that is ε (0) α = 0 for α ≥ 1, as shown in Table 1. In JJ19, we took the first-degree mechanism (α = 1) to maximize the dipole moment ( p = 1) of the first-degree residual field, the second-degree mechanism (α = 2) to maximize the quadrupole moment ( p = 2) of the second-degree residual field, and so on. This optimization minimizes the residual moment μ (α) ( ˙ α+1 ): In JJ19, we incorrectly asserted that each step in this sequential orthogonalization procedure completely removes the target multipole from the residual field; for example the dipole moment of the second-order residual field is identically zero. For a general source, this is not true; the subdiagonal terms in the moment array (i.e. those with 1 ≤ p < α, labelled ε ( p) α in Table 1) will be non-zero and must be retained to properly represent the wavefield and pose inverse problems.
The displacement wavefield u(x) can be expressed as the integrated contraction of˙ (x ) with a third-order stress-displacement Green's tensor H(x; x ). Expanding H in a Taylor series about x * represents the wavefield as a summation over the moment array μ ( p) (g α ): The factor in brackets is the fully contracted (scalar) operator μ ( p) klm··· (g α ) ∂ k ∂ l ∂ m · · · (p derivatives). We assume the displacement wavefield is band-limited with a centre wavelength L, which we use to non-dimensionalized the derivatives,∇ p x * = L p ∇ p x * . To nondimensionalize the moments, we use the characteristic source size J (0) and length L 0 , The wavefield for each degree can then be written In the point-source limit, L 0 /L → 0, the expansion point x * becomes indistinguishable from the centroid x 0 , and the non-zero moments with the smallest p will tend to dominate the sum. The zeroth-degree wavefield can thus be approximated by the monopole term, and the first-degree field by the dipole term: However, the leading terms of the higher-degree fields u α , α ≥ 2, are not source multipoles of degree α, as claimed in JJ19 (eqs 39 and 40), but are typically dipoles corresponding to the terms ε (1) α . This correction affects the composition of inverse problems for sources with high degrees of mechanism complexity (D ≥ 3). We consider the Bayesian estimation of moment-array elements for a particular source from its observed wavefield u obs . These data are inverted by a sequential orthogonalization procedure that increases the maximum degree of the target MOR at each step; for example at the αth step, a vector of data functionals f α = F α (u obs ) is inverted for an MOR up to degree α conditional on the results for α − 1. As the notation indicates, we allow the data set to vary from step to step. For example, as α increases, the data set may be augmented with seismic observations at smaller wavelengths to better resolve the higher-order moments.
The first step (α = 0) estimates the CMT parameters (up to 10) from f 0 using the source monopole (17) as the wavefield model. The values of (M 0 , J (0) , x 0 ) that maximize the posterior model probability will approximate the true source monopole. In the second step (α = 1), the data-residual vector f 1 = F 1 (u obs − u 0 ) is inverted using the dipole form (18) of the first-degree wavefield ( p = 1). The first-degree mechanismM 1 and the dipole vector J (1) that maximize the posterior model probability conditional on f 1 , (M 0 , J (0) , x 0 ), and the orthogonality constraintM 0 :M 1 = 0 will approximate the largest term in the CPD of μ (1) ( ˙ 1 ). From eq.
If the data warrant going further (i.e. if the weight of evidence from u obs supports D ≥ 2), we proceed with the next step (α = 2), which is to invert the residual vector f 2 = F 2 (u obs − u 0 − u 1 ) for moments of order p ≤ 2. The corrected model for the second-degree wavefield includes a first-order moment, Accounting for the orthogonality conditions (M 0 :M 2 =M 1 : specify second-order tensor J (2) . The norm ε (1) 2 is maximized by the CPD mechanism U (1) 2 at the singular value λ (1) 2 , whereas J (2) is maximized by the CPD mechanism U (2) 2 at λ (2) 2 . If the data provide sufficient information, the estimateM 2 obtained by maximizing the posterior model probability conditional on f 2 and the parameter priors (posteriors from the previous step) will lie somewhere between these two mechanisms, depending on the relative sizes of the two terms in (19).
We define the direction cosines cos θ = U (1) 2 : U (2) 2 and cos ψ = M 2 : U (2) 2 , assuming without loss of generality that 0 < θ ≤ π/2; then, the estimate can be written, where 0 ≤ ψ ≤ θ is a parameter that dials between the two CPD mechanisms. The mechanism that best fits the data set will approximate the one that maximizes the norm of the second-degree wavefield. The norms of the two terms are Using the scaling in eq. (15), we find If the normalized moments are comparable in magnitude, ε (1) 2 ≈ Ĵ (2) , then the dipole term will dominate (ψ ≈ θ), owing to the extra wavelength factor in the denominator of the quadrupole term, in which caseM 2 ≈ U (1) 2 . Of course, if the true rank of˙ were D = 3, then U (1) 2 will be identical to U (2) 2 (θ = 0). This inversion procedure can be continued to higher values α and p until weight of evidence from u obs does not support an increase in the number of source parameters. Of course, the number of new parameters proliferates as we go higher degree (e.g. 20 are needed just to specify J (3) ), and their contributions to the low-frequency wavefield decrease as (L 0 /L) p . From a practical standpoint, the utility of a moment-based approach to the source inverse problem is likely to be limited to low degree, although some complex sources, such as the 2016 Kaikōura earthquake, may require expansions as high as D = 5 to fully characterize the stress glut.

C O R R E C T I O N T O N U M E R I C A L R E S U LT S
In JJ19, we illustrated the MOR using a realization of Graves & Pitarka's (2016) kinematic source model. The realization is an M w 6.85 rupture simulation on a vertical fault with a strikeslip mechanism perturbed by local variations in the fault strike, rake and dip ( fig. 3 of JJ19). The Aki moment of this event is M 0 = 2.13 × 10 19 Nm, and its characteristic outer scale derived from μ (2) 0 is L 0 = 47.8 km. We recomputed the MOR for this realization using the CPD optimization algorithm of section 2, which we implemented using the Tensorlab 3.0 software package (Vervliet et al. 2016). For the purposes of comparison with JJ19, we construct the MOR using the JJ19 optimization of eq. (13); that is by taking ψ = 0 in eq. (20). The g α fields, shown in Figs 1(a)-(e), differ slightly (by a few per cent or less) from those in fig 4 of JJ19, owing to an incorrect assumption in the JJ19 algorithm that effectively fixedM α J (α) to be the largest diagonal component of μ (α) ( ˙ α ), rather than the projection that has the largest singular value. The mechanism tensorsM α are also nearly the same, exceptM 3 , which is reversed in fig 4 of JJ19 owing to a sign error in the plotting code. Table 2 displays the non-dimensionalized moment array of the MOR, μ ( p) (g α ) , evaluated at the event centroid, and it also list the seismic moment fractions for each degree α. In JJ19, we defined the fractional momentsḾ α in terms of the direction cosines of the local source mechanism, cos θ α (x) =M α : ˙ (x). Substituting the identity α cos θ 2 α (x) = 1 into eq. (5), we obtained the totalmoment sum, For sources with mechanism complexity, the zeroth-degree fractional momentḾ 0 = V m(x)cos 2 θ 0 (x)dx differs from the Aki moment M 0 = V m(x) cos θ 0 (x)dx. Our corrected results yield moment fractionsF α ≡Ḿ α /M T close to those given in JJ19, withF 2 showing the largest difference (0.006). The seismograms and radiation patterns calculated for the higher-degree terms are also very similar to those given in figs 5 and 6 of JJ19. Fig. 1 compares the MOR and PCA representations of the Graves-Pitarka realization. The leading term in the PCA repre-sentationN 0 h 0 is nearly identical toM 0 g 0 , and μ (0) (h 0 )/ √ 2 differs from the Aki moment μ (0) (g 0 )/ √ 2 by less than three parts in 10 4 . The two optimization methods also yield similar moment-tensor fields at higher degrees, although their ordering is not the same. In particular, the second-and third-degree fields are switched-M 2 g 2 ≈N 3 h 3 andM 3 g 3 ≈N 2 h 2 -which reflects the inequalities μ (2) (g 2 ) > μ (3) (g 3 ) and h 2 > h 3 . The PCA orders the orthogonal fields in a decreasing sequence of fractional moments, whereas the MOR does not. Aside from this ordering difference, the two expansions are essentially the same.
These representations have a physical interpretation related to the composition of the source model. The Graves-Pitarka realization involves mechanism rotations due to out-of-plane strike variations and in-plane rake and dip variations. These angular variations are constructed to be statistically independent (Graves & Pitarka 2016), and they are well characterized by von Mises distributions with zero means and standard deviations of 10 • , 15 • , 10 • , respectively ( fig. 3  of JJ19). In JJ19, we showed that double-couple distributions with a single rotational degree of freedom generate moment fields with a similar set of orthogonal mechanisms. Applying the von Mises model to the Graves-Pitarka realization ( fig. 1 of JJ1) yields fractional moments with expected values of 0.130, 0.066 and 0.030, respectively. Even though the standard deviation of the strike variations is equal to that of the dip variation (10 • ), its fractional moment is about four times larger, owing to the greater effect of the out-ofplane rotations on this measure of mechanism complexity.
Out-of-plane strike variations generate a double-couple source mechanism that is rotated by 45 • about the vertical axis relative to the average strike-slip mechanism ( fig. 2 of JJ19). The first-degree term,M 1 ≈N 1 , approximates this mechanism, and the first-degree moment fractions, 0.108 for the MOR and 0.113 for PCA (Tables 2  and 3), are close to the expected value of 0.130. Likewise, the mechanismM 3 ≈N 2 corresponds to rake variations, and it yields Downloaded from https://academic.oup.com/gji/article-abstract/222/2/1333/5856841 by guest on 14 June 2020   Values in boldface were incorrectly assumed to be zero in JJ19.  Owing to the statistical fluctuations, the correspondence ofM 2 ≈N 3 to dip variations is less precise (e.g. the CLVD fraction ofM 2 is almost 40 per cent), but this field accounts for only a small fraction of the total moment (0.035 and 0.028, respectively). This fraction is again fairly well predicted by the von Mises mean (0.030). In JJ19, the higher-degree residuals ε (α) β>α defined by the MOR of eq. (13) were incorrectly assumed to be zero. The size of these residuals is measured by the ratio [Note that (24) does not imply term-by-term equality; that is λ In the Graves-Pitarka example, these residuals are not zero, but they are relatively small: r 1 = 0.146, r 2 = 0.120 and r 3 = 0.008, indicating that the J (α) terms dominate the expansion of μ ( p) ( ˙ α ). However, J (α) generally will not dominate the wavefield u α . From Table 2, for example, we see that ε (1) 2 ≈ 0.6(L/L 0 ) J (2) and ε (1) 3 ≈ 0.7(L/L 0 ) 2 J (3) , so in this case the dipole terms will be the largest if the source is small compared to a wavelength.

D I S C U S S I O N A N D C O N C L U S I O N S
The MOR derived in JJ19 is incomplete, because it includes only rank-one approximations to the residual stress-glut moment Downloaded from https://academic.oup.com/gji/article-abstract/222/2/1333/5856841 by guest on 14 June 2020 μ (α) ( ˙ α ). Here we have corrected the MOR by adding the higherdegree terms ε (α) β>α to the moment array. For α ≥ 2, these terms need not be small and should be included in the representation of sources with high degrees of mechanism complexity.
The revised MOR provides a theoretical framework for inverting low-frequency wavefield observations for stress-glut integrals that are insensitive to the geometrical details of the rupture geometry. In the point-source (low-frequency) limit, the zeroth-degree wavefield u 0 can be approximated by a CMT monopole with a mechanism M 0 and a scalar moment J (0) , and the first-degree wavefield u 1 by a mechanismM 1 orthogonal toM 0 and a dipole with vector moment J (1) . Expanded to p = 2, the second-degree wavefield u 2 includes a quadrupole with tensor moment J (2) , but also a dipole with vector moment ε (1) 2 . The latter will tend to dominate the radiation at low frequencies and will thus determineM 2 .
We have outlined the use of a moment-oriented representation in a progressive Bayesian inversion scheme that increases the maximum degree of the target MOR at each step, thus building up the joint probability distribution of the moment-array parameters conditional on the observations. At the αth step, the data residual vector f α = F α (u obs − α−1 β = 0 u β ) is inverted for (M α , J (α) , ε (1≤ p<α) α ). As α and p increase, the parameters proliferate, and their contributions to the wavefield decrease, limiting our ability to estimate the complete parameter set.
When to terminate the moment expansion will depend on the the effective rank of the stress glut and the resolving power of the data. For sources with simple mechanisms, such as ruptures on nearplanar faults, the effective rank will be close to one, and the zerothdegree (CMT monopole) representation (4) will be an adequate approximation of the source. For complex sources, higher-degree moments will be necessary to fully characterize the stress glut and account for the total moment. For example, we have computed both the MOR and PCA representation of the Wang et al. (2018) model of the 2016 Kaikōura earthquake and find that the total-moment fractions are significant (>5 per cent) for all members of a complete basis of deviatoric mechanisms, thus requiring a representation of rank D = 5 (Juarez & Jordan 2019).

A C K N O W L E D G E M E N T S
We thank Jeff McGuire for a helpful review. This research was supported by the Southern California Earthquake Center (SCEC) under National Science Foundation Cooperative Agreement EAR-1033462, U.S. Geological Survey Cooperative Agreement G12AC20038, and grants from the W. M. Keck Foundation and Pacific Gas & Electric Company. The SCEC contribution number for this paper is 9997.