Phase Transformation Under High Pressure Radiates as a Double Couple Deep Earthquake

Deep‐focus earthquakes (DFEs) originating at the Mantle‐Transition‐Zone (MTZ) (400–700 km) have a Double Couple (DC) radiation pattern similar to crustal earthquakes; however, their mechanism is different and governed by high pressures (15–25 GPa) at nucleation depths. We present a model of nucleation and growth of regions of phase transformation, undergoing a sudden reduction in volume (5%–10%), “volume collapse.” Successive symmetry‐breaking instabilities minimize the energy spent to move the boundary of phase discontinuity and a collapsing volume expands as a flattened pancake‐like self‐similarly expanding Eshelby ellipsoidal inclusion. At the vanishing of the M integral, expressing the balance of flows of energy across the inclusion boundary, at a critical value of the pressure, an arbitrarily small inclusion nucleates and grows at constant potential energy driven by the pressure acting on the change in volume. The inclusion develops shear eigenstrains that decompose into two DC, placing one on the basal plane to radiate without energy losses. The symmetric volume collapse radiates out as an anti‐symmetric DC, and the radiated energy is obtained as the “excess energy,” of the ambient pressure acting on the “volume collapse,” reduced by the energy consumed for the growth of the pancake surfaces, with a “pressure drop” (p0 − pcr) driving the expansion emitting a DC, even under full isotropy. The solution explains some features of the DFEs, (a) the DC radiation, (b) their large energies (the Mw 8.3 Okhotsk earthquake [2013]), (c) the absence of volumetric radiation, and (d) why they can originate in the MTZ, a long‐standing open problem.


10.1029/2023JB026797
2 of 16 The "transformational faulting" model for example, from metastable olivine to denser wadsleyite and to ringwoodite, was first introduced by Green and Burnley (1989), as densified "anticracks," also Frohlich (1994).The Okhotsk earthquake (8.3 Mw at 609 km depth) was attributed by Schubnel et al. (2013), to an instantaneous phase change from olivine to spinel.The transformation faulting model is assumed in the literature to be driven by the deviatoric stresses.However, the estimated deviatoric stresses at these depths are in the range 25-100 MPa, -for example, Dziewonski and Anderson (1981) and Dziewonski et al. (1981) and not large enough to drive earthquakes of such a magnitude.For the Okhotsk (2013), earthquake the deviatoric stresses were assumed at 150 MPa, and did not provide enough energy to explain their large radiated energy, giving an efficiency coefficient >1 (Ye et al., 2013).Seismological observations, and laboratory analog experiments, show that a reduction of ∼5% to 10% in volume (density increase) associated with the high-pressure transformation (as the atoms find equilibrium positions closer to each other in the crystal lattice undergoing a phase transformation) can produce substantial acoustic emissions, for example, in analog experiments by Meade andJeanloz (1989, 1991), and Schubnel et al. (2013), see also Frohlich (1994Frohlich ( , 2006)).Moreover, the large energy released by the pressure acting on the change in volume has to be accounted for.This raises the question of what role the pressure plays in the initiation and propagation of these acoustic emissions, and the DFEs, and why it does not hinder their initiation as intuitively expected.These questions are addressed here in a dynamic Eshelby inclusion model of a sudden change in volume under high pressure in full isotropy, as the Randall model ( 1964), but, which, due to instabilities, yields results consistent with the observed radiation of a predominantly DC seismic source, and shows that sufficiently high pressure indeed drives the DFEs (acting on the change in volume), providing their large energies.
The hypothesis that pressure-induced phase transformations with sudden volume collapse cause deep-focus seismicity has been advocated over the years (Benioff, 1963;Randall, 1964).It was already stated by Randall (1964) in his abstract that " such phase transition (sudden volume change), even for small change in density, offers a much more concentrated source of seismic energy than does concentrated faulting," and this is at the core of the fundamental physics of the phenomenon that produces these large earthquakes.However, there has been a puzzle, since, a sudden volume change should produce radial displacements, which were not observed, while, on the contrary, were observed shear motions predominantly of a Double-Couple seismic source, with little or no volumetric component (e.g., Frohlich, 1994Frohlich, , 2006)).These observations have thrown doubt on the model of sudden volume change, which received no further consideration and neither did the effect of pressure.This paper helps resolve the puzzle, as the unstable mode of the Randall problem (Figure 1) of sudden volume collapse under high pressure, by considering the energy spent to expand the boundary of phase discontinuity.
The phenomenon of an expanding inclusion of phase change differs from fracture.In the Meade and Jeanloz (1989) experiments, where the specimen is subjected to 70 GPa and acoustic emissions are recorded, it is described as "displacive phase transformation, not fracturing or cracking of the samples… due to the very high pressure above the brittle to ductile transition."This is a new defect, the flattened densified dynamically expanding Eshelby inclusion with transformation strain   *  under high pressure; it models the mechanics of the anticrack of Burnley and Green (1989) as an expanding region of phase change.The phase transformation is modeled as a dynamically expanding Eshelby (1957) ellipsoidal inclusion with transformation strain, or eigenstrain, which was first introduced in geophysics as a source for seismicity by Backus and Mulcahy (1976), called "stress glut," also Segall (2010) for wider applications to geophysics, Meng et al. (2012) for numerical applications.The transformation strain was introduced by Eshelby in the famous thought experiment, where an element is taken outside the material, it undergoes a (transformation) strain   *  , and, in order to be inserted back into the hole (undeformed configuration), corresponding stresses are applied induced by imaginary body forces (see also, Eshelby, 1975).One may then think of the eigenstrains as distributed body forces inside the volume of the flattened pancake.
The Eshelby problem, that has the important constant stress property in the interior of the ellipsoid (-which allows the inclusion to move as a whole-), assumes that the material moduli (at the macroscopic scale) are uniform and that the transformation strain is also uniform inside the ellipsoid.Although there is variation in the moduli in geophysical materials it is assumed that the variation in the moduli within a thin layer at a given depth in the subducting zone, where the phase transformation occurs, is small.In the deep earthquakes analog experiment the geological material (olivine) at the micrometer scale consists of grains of different shapes of size <10 μm (e.g., Schubnel et al., 2013, Figure 4), but this does not affect the homogenized modulus at the macroscale, which can vary slowly, so that the assumption of uniform material moduli in the Eshelby inclusion is approximately satisfied.The uniform eigenstrain assumption, corresponding to uniform densification in the phase transformation, is reasonable if the material is approximately uniform undergoing the same phase transformation.The ellipsoid can expand in MARKENSCOFF 10.1029/2023JB026797 3 of 16 regions where the Eshelby property holds: Brown (2012), states that, "as it expands and encounters obstacles, the ellipsoid can elongate or change both its axial ratios and the orientation of its axes in any way, so long as it remains ellipsoidal."The model of Green and Burnley (1989) assumes that transformation faults self-organize forming a bigger fault.Also, the subducting lithosphere is assumed cold (e.g., Ye et al., 2013), so no temperature effects would enter the constitutive relation.Whether the elasticity Eshelby solution applies at these depths may possibly be questioned.As the eigenstrains in the Eshelby problem are equivalent body forces and the solution is obtained based on the dynamic Green's function (for a point load), the dynamic Eshelby problem should be applicable as long as the dynamic Green's function is; it was assumed to hold by Knopoff andRandall (1970), andRandall (1964), in their models of deep earthquakes due to phase transformation.Experiments on single-crystal olivine under high pressures by Mao et al. (2015) confirm elastic wave propagation with the moduli dependent on the pressure (at a given depth).
The deep earthquakes, modeled here as regions of dynamically expanding phase transformations, are governed by the energetics of the moving phase boundaries, that dictate the shape and the deformation as to minimize the energy spent to move the boundary of phase discontinuity in a given stress field, with the shape, in turn, determining the radiation.The assumption that the deformation fields are those of a self-similarly expanding Eshelby inclusion provides properties that are intricately related to the energetics.Modeled as regions of self-similarly expanding Eshelby inclusions (Ni & Markenscoff, 2016a), allow the maximum radiation of energy, since, in the interior of the self-similarly expanding inclusion the particle velocity is zero (a lacuna), and no kinetic energy is lost inside; this also allows the phase transformation to take place under equilibrium conditions.Thus, the input  Randall (1964): Stable solution of sudden volume change preserving symmetry and emits volumetric radiation; "it is a more concentrated source of energy than transformation faulting" (Randall).(b) Under high pressure, and high densification   *  , in full isotropy, successive symmetry breaking instabilities minimize the energy spent to move the boundary of phase discontinuity which expands planarly as (c) a densified concentrated "blade-like" inclusion containing large shear transformation strains.The eigenstrains decompose (b) into two DC shear sources: the DC1 with eigenstrains on planes parallel to the basal plane radiating through the expanding perimeter without energy loss, while the DC2, on perpendicular planes, assumes the expansion loss.(d) The direction of the displacement around the pancake is indicated corresponding to pressure P (positive and negative) and shear S waves, emanating the radiation of a DC.10.1029/2023JB026797 4 of 16 energy goes to nucleate/grow the boundary of the expanding pancake, and the rest radiates out as seismic energy.To obtain the radiated energy is an objective of this study.Markenscoff (2019aMarkenscoff ( , 2020Markenscoff ( , 2021a) ) assumed the deformation fields to be those of a self-similarly expanding Eshelby ellipsoidal inclusion with uniform volume collapse eigenstrains and considered the energetics of moving phase boundaries according to Noether's theorem (1918) of the calculus of variations in a variable domain; they determine the shape, that, in turn determines the radiation.She showed that the problem of Randall (1964), exhibits an instability (Figure 1), under high pressure, and acquires a flattened "pancake-like" shape, as to minimize the energy spent as the inclusion grows large, rather than the symmetry-preserving spherical shape (Markenscoff, 2020).To note, the pancake shape is a shape that avoids the boundary becoming a sink or source of energy, as an alternative possibility to the symmetry-preserving spherical one (Markenscoff, 2019a).The absence of volumetric radiation was explained (as canceled by the pressure), and due to symmetry in material properties, geometry, and loading a symmetric compensated linear vector dipole (CLVD) shear seismic source was obtained (Markenscoff, 2021a).However, the outward radiation was not studied in Markenscoff (2021a) and the question remained why the vast majority of the observations for the DFEs are those of a DC seismic source type.
The explanation for this fundamental feature of DFEs was obtained here by obtaining a successive instability, which minimizes the energy spent to expand the perimeter as well, so that only one DC radiates through the perimeter without additional energy losses.There will be some CLVD radiation if the boundary of the ellipsoid is not totally flat, but the higher the pressure the flatter is the ellipsoid, consistent with the fact that the deeper earthquakes are DC (e.g., Okhotsk, 2013).The balance of flows of energies through a moving phase boundary separating a region with eigenstrain is developed here (Figure 3): when at a critical pressure the outward flow cancels the inward one, and there is no interaction energy lost to create a phase boundary, nucleation can occur, and an arbitrarily small inclusion can grow at constant potential energy, while the excess energy is radiated into the bulk (Figure 3).This illustrates the physics of the M integral of theoretical physics (Budiansky & Rice, 1973;Jackiw, 1972).Markenscoff (2020Markenscoff ( , 2021a) ) obtained (respectively for the sphere and the pancake) the critical pressure at which an arbitrarily small inclusion will grow at constant potential energy at constant rate of densification.For an earthquake to occur, the critical nucleation pressure p cr must be smaller than the corresponding ambient pressure p 0 at the depth at which the earthquake is nucleated; the difference (p 0 − p cr ) provides the "stress drop," here a "pressure drop" acting on the change in volume eigenstrain   *  .The critical pressure to nucleate a pancake-shape inclusion is (Markenscoff, 2021a) where ν is Poisson's ratio (equal to minus the ratio of the lateral strain over the longitudinal) and μ the shear modulus (of the untransformed material in the Eshelby solution).The values for Table 1 were obtained based on Abramson et al. (1997), Dziewonski andAnderson (1981), andMao et al. (2015).
Note.The depth of 609 km corresponds to the Okhotsk earthquake (2013) and the phase transformation is assumed to be the one indicated with the orange arrow.The "pressure drop" (p 0 − p cr ) acting on the change in volume   *  drives the expansion.For these phase transformations in the MTZ the ambient pressure exceeds the critical and the DFE can occur.1 the values of the critical nucleation pressure are lower than the ambient pressure p 0 at the corresponding depth of the earthquake and the radiated energy is obtained in the Methods/Results section.Moreover, it is seen that the critical pressure exceeds (and cancels) the (tensile) mean stress due to the volume collapse in the three transformations indicated in Table 1, so that there is no volumetric radiation.Change in shear modulus and bulk modulus in the phase transformations are not accounted for here, but the "Eshelby property" (constant interior stress) and the energetics of the M integral are also valid in full anisotropy (Willis, 1971) and in dynamic expansion.For the values in Table 1, the "pressure drop" is p 0 − p cr = 2.9 and 7.03 GPa, respectively, and acting on the collapsed volume drives the earthquake.For the Okhotsk earthquake at 609 km (2013) the pressure drop is 8.5 GPa.Whether, for larger densifications and larger moduli of the geological materials below 700 km, the critical pressure for nucleation and growth becomes higher than the ambient pressure, not allowing an earthquake to nucleate, has to be examined depending on the phase transformation, and, this might be one explanation of why the earthquakes stop below a certain depth in the Mantle.The important result obtained in this paper is the radiated energy, which is computed for the fault data of the Okhotsk (2013) earthquake and shows the ample energy provided by this model.

Under High Pressure a Self-Similarly Expanding Flattened Eshelby Ellipsoidal Inclusion With "Volume Collapse" Emits the Radiation of a Double Couple Even Under Full Isotropy
A solution of mechanics is obtained for a dynamic phase transformation with "volume collapse" in full isotropy under high pressure exhibiting a successive instability to the one obtained in Markenscoff (2021a).
The analysis is based on the self-similarly expanding Eshelby ellipsoidal inclusion with uniform eigenstrain (Ni & Markenscoff, 2016a, 2016b) with constant axes speeds 1/s i (with s i denoting the slowness) with 1/s 3 ≪ 1/s 1 , s 1 t/s 3 t = a 3 /a 1 → 0, satisfying the conservation of linear momentum (see Text S1 in Supporting Information S1).The displacement vector u i is continuous, and with the summation convention implied for repeated indices i, j = 1, 2, 3, we have the following relations for the strains, stress and strain energy density for the Eshelby (1957) inclusion problem (Table 2).
The displacement is continuous, and for the static ellipsoid the traction is continuous across the interface, while the hoop stress components are discontinuous.In the presence of inertia, for the expanding ellipsoid, the traction is also discontinuous (jump conditions in Equations 1 and 2).To the transformation due stress σ ij is added the pre-stress p 0 δ ij to obtain the total stress.
The fields in the pancake-like flattened ellipsoid were obtained in Markenscoff (2021a) following an asymptotic analysis where the eigenstrains are assumed large, the strains small, the pressure large, the deviatoric pre-stresses small, the strain energy density large and the total energy of the very thin inclusion is finite.The large collapsing volume produces to the leading order a large tensile mean stress   = − *  , where K is the bulk modulus in the untransformed state, which is canceled by the ambient pressure if larger (Table 1).For the material not to break in tension the superposed compressive pre-stress must be large enough to cancel the tension, which implies that the pancake solution is valid for this order of ambient pressures, consistent with the range of the deep earthquakes, as exhibited in Table 1.Assuming elastic strains of the order of 10 −3 or smaller, for a change in volume of the order   *  ∼ 5%-10%, the ratio of the axes speeds/lengths (small to large axis) of the ellipsoid should be of the order s 1 /s 3 ∼10 −4 or smaller.To note that in analog experiments (Wang et al., 2017) the thickness of the faults was 100 nm while the length could be millimeters (making for a ratio of 10 −4 ); however the pressures are 2-5 GPa in the Schubnel, 2013, experiments.
The total strain The elastic strain eij is the total minus the transformation strain The strain energy density   , which constitute a symmetric center of shear, expected from the overall symmetry of the problem.As the deviatoric pre-stress is of smaller order than the pressure, the shear stresses that are of the order of the eigenstrains are very large, producing a large distortional strain energy density.This explains the generation of very large shear stresses, a question posed in Ye et al. (2013).The Eshelby property of constant interior strain (valid only for the ellipsoid, Markenscoff, 1998a) is obtained also in self-similar expansion, and, consequently, in the interior domain the particle velocity must be zero.Indeed, the Eshelby ellipsoid expanding self-similarly (with constant subsonic axes speeds) emits P, S, and M waves, first presented in Burridge and Willis (1969), for the self-similarly expanding elliptical crack, and in Ni and Markenscoff (2016b), for the expanding inclusion, see Text S1 in Supporting Information S1.The M waves cancel the P and S in the interior domain where the particle velocity is zero (shown explicitly Ni & Markenscoff, 2016a) and is based on analytic properties alone (Markenscoff, 1998b).More generally, the interior domain can be predicted mathematically to be a lacuna, a region in which the solution vanishes, which is a topological property of hyperbolic systems of partial differential equations (Atiyah et al., 1970), qualitatively described by Burridge (1967Burridge ( , 1971)), as a region "of absolute quiet in a sea of waves."The M waves are emitted by the boundary of the expanding ellipsoid and satisfy the Hadamard jump conditions across the moving interface (Markenscoff, 2021b).As shown in Burridge and Willis (1969) the M waves have the boundary of the expanding ellipsoid as their degenerate wave-front, and they yield Rayleigh waves when the crack expands at the Rayleigh wave-speed.For the spherical inclusion the M waves were obtained in Markenscoff (2021b), where it was shown that the dynamic Eshelby problem properties also hold for a Newtonian fluid in a phase transformation with different viscosity, as the static Eshelby property also holds for a Newtonian fluid (Eshelby, 1957).One may note here that the Eshelby ellipsoidal inclusion is a very rich solution with all the components of eigenstrains, with three of them that yield cracks in a particular limit (Burridge & Willis, 1969;Markenscoff, 2019aMarkenscoff, , 2019b;;Mura, 1982), and the others giving new defects of inclusions of phase transformations that manifest themselves in phenomena such as this one.
Across a moving surface of discontinuity the jump conditions in the traction and particle velocity express the conservation of linear momentum and compatibility, respectively, and are given in terms of the normal boundary velocity  ̇ by the Hadamard jump conditions (1) where the double brackets denote jumps across the surface of discontinuity.From the jump conditions we can determine the external to the inclusion limits.We extend to dynamics a procedure described for statics in Mura (1982, Equation 6.4), and write the unknown jump of the displacement gradient, at a point with normal   on the moving interface, in the form where the vector λ k is to be determined making use of the Hadamard jump conditions 1 and 2; it is obtained in Text S2 in Supporting Information S1 as where, for isotropic material, This is a point-wise solution for the fields approaching in the limit the moving boundary from the exterior.
The normal n at a point x on the surface of the ellipsoid is expressed in terms of the three axis lengths as ).The expressions 5a and 5b applied to the flattened surfaces of the ellipsoid under consideration here, with the normal boundary velocity tending to zero as  ̇ → 0 , to the leading order, give the static solution, which is consistent with experiments in failure waves (e.g., Clifton, 1993).However, by the above procedure we cannot obtain the fields in the limit of the sharp tip perimeter (Figure 1), which would require the evaluation of the limit of the full solution for the ellipsoid of Ni and Markenscoff (2016a), being a topic of future investigation.
We present now the non-unique decomposition of a CLVD into two DC (e.g., Julian et al., 1998;Knopoff & Randall, 1970) and show why the CLVD may choose to decompose in a way that minimizes the losses due to the interaction energy to expand the perimeter.A CLVD can be decomposed asymmetrically into the two DC, shown as an example of a non-unique possibility in Knopoff and Randall (1970), depicted in Figure 2, as with

𝑘𝑘𝑘𝑘
The total potential energy of a mechanical system is the Gibbs free energy in isothermal processes, or the Helmholtz free energy in adiabatic ones, and is defined as the strain energy minus the potential energy of the loading mechanism on the boundary.In the case of an inclusion with volume collapse under loading mechanism   (0)  , we write (Eshelby, 1961;Mura, 1982, Equations 13.9-13.16)the total potential energy of the system as the sum, E tot = W 0 + W 1 + ΔW, where W 0 is the potential energy in the absence of eigenstrain, W 1 the potential energy due to the inclusion alone, which is the strain energy in the whole body due to the inclusion, which is  − ∫ Ω  *    (a positive quantity), and ΔW is the remaining part of the potential energy, due to the interaction of the loading mechanism (the pre-stress) with the inclusion, "appropriately called interaction energy" (Eshelby, 1961), and it is given by (Eshelby, 1961;Mura, 1982, Equation 13.16) We note that, for negative volume collapse eigenstrains and compressive applied stress, the first term in Equation 7 is negative and can overcome the second positive term at a critical pre-stress   (0)  .The presence of the second term indicates that upon removal of the applied loading an inclusion with phase transformation of volume collapse will not go back to the untransformed configuration.Figure 3 below illustrates physically the interaction energy for a moving phase boundary.
We then examine the interaction energy in view of the decomposition of the CLVD into the two DCs in Equation 11,

𝑖𝑖𝑖𝑖
. We compute the interaction energy of the basal plane DC1 (1/2) (−1,1,0) with the pre-stress   (0)  =  0  stresses: the stresses σ ij are those in the "pancake" inclusion due to densification as given by Equation 6in Markenscoff (2021a) Thus, the DC1 (1/2) (−1,1,0) has zero interaction energy with the pancake's self stresses, as well with the symmetric pressure so that the interaction energy in Equation 7 vanishes.By contrast, for the CLVD eigenstrains and the pancake stresses, the second term in Equation 7 would be a positive quantity (spent energy).Moreover, it can be seen that, in the limit, as the perimeter expands with normal in the x 1 direction, the DC2 (eigenstrains   *
In the in-plane DC1 (1/2) (−1,1,0) is expressing the physics that the interaction energy of the inclusion with the pre-stress field is minimized, which is developed in the sequel and shown in Figure 3.In Equation 6a or (1, −1,0) in Equation 6b, the eigenstrains   * 11 = − * 22 give an equivalent shear   * 12 eigenstrain at 45°, and the question is along which radial direction of the circular pancake it will occur, since the axes are taken in arbitrary directions in full radial symmetry (Figure 2).The direction of the maximum shear pre-stress   (0)  12 on the Ox 1 x 2 plane, no matter how small, gives the direction of the perturbation; it is not driving (to the leading order) the transformation of volume collapse, which is driven by the pressure acting on the change in volume.The DC2 on the vertical planes are secondary deformation systems with the obtained instability showing how they are produced in a flattened ellipsoid with densification.These secondary non-radiative shear localization systems appear to be consistent with the "tangled dislocations" observed in the first experiments in olivine amorphization at 56 GPa (Jeanloz and Ahrens, 1977) and also with the experimental findings of both shear nanobands (NSB) and densified lenses in olivine phase transformations in Wang et al. (2017); they were brought up in the discussion in Zhan (2020, Figure 5), stating that "the NSBs are filled with fine-grained spinel accommodating the seismic slips."

The Energetics of Moving Phase Boundaries and the Radiated Energy Due To "Volume Collapse" Under High Pressure
We will obtain in this section more generally the energetics of a moving surface of phase discontinuity.The surface of the inclusion containing transformation strain   *  is a surface of discontinuity of the displacement gradient, and consequently of the strain energy density, and, also, in dynamics, discontinuity of the traction and particle velocity (Equations 1 and 2), which determine the net flow of energy across the moving interface.We consider the equation of motion   −  2 ∕ 2 = 0 and, as, for example, in Freund (1990, p. 224), we take the inner product with the particle velocity  ∕ , and obtain the energy-momentum equation where U is the stress work is the kinetic energy density at any material point.
We note that (Eshelby, 1970, Equation 54) the energy flow vector is   = − ̇ , and that the energy flow in the direction of the positive normal is   = − ̇ , so that for increase of energy in the volume the energy flow is inwards.
We consider (Figure 3) an arbitrary volume of material possessing a quantity F, where F = (U + T), the sum of the strain and kinetic energy densities, and containing an internal surface Σ separating two phases, V + and V − , and moving with an arbitrary velocity   ̇⃗  .We apply the Reynolds Transport Theorem (e.g., Aris, 1962, pp. 85-86) as developed for fluids (with  ⃖ ⃗  being the particle velocity) and write the total time derivative, which equates the rate of change of the total F possessed by the volume the change of F produced the volume minus the flux of F through the surface.We apply the R.T.T. to the volumes V − and V + separately (Figure 3), as The first integral on the right-hand-side (RHS) of Equation 10expresses the rate of increase of (U + T) energy produced in the volume; this differs in the two volumes as in V − there is transformation strain affecting all quantities.The application of the divergence theorem to Equation 9 writes this term as a surface integral producing a jump contribution across the surface of discontinuity Σ. Equation 9 can be integrated for any arbitrary time interval (t 1 , t 2 ), and hence the time integration can be omitted so that the expressions are rates of energy flow.
In the limit as Σ − → Σ + → Σ we write in the sum of the terms Equations 10a and 10b the contribution of the jumps of the fluxes of energy across the moving interface, and their directions, obtained as The variation of the total energy of the body with the expansion of the surface of discontinuity is obtained as a surface integral over the moving boundary Σ in Equation 10c, demonstrating the effect of the energy-momentum tensor (e.g., Eshelby, 1970Eshelby, , 1975) ) in Equation 9on the energy flux through the boundary.With the application of the Hadamard Equations 1 and 2, Equation 10c acquires the expression obtained in Text S1 in Supporting Information S1 as Equation S3.5 where  ̇ is the normal outward velocity at a point of the moving surface of phase change.
The expression in Equation 11 is the same as the one given in statics by Eshelby (1970, Equation 43), also Eshelby (1977), for the force on an interface and as obtained by Markenscoff (2019a).The integrand in Equation 11 is the product of the "driving force" =<  >  *  (  = 1 2 3) , acting along the outward normal to the phase boundary, times the boundary velocity  ̇ .With the self-stresses of the inclusion the "driving force" in Equation 11 is negative, showing that there is a net energy flux into the inclusion from the Σ + surface, hence, showing that there is energy lost inside the moving phase boundary between the surfaces Σ − and Σ + in the limit.This is expected in the presence of inertia for the same reason as in statics, where, as Eshelby writes (1961, Equations 2.30-2.32), in the Eshelby thought experiment, when the transformed material is inserted back in the untransformed shape, and welded in the body, energy is extracted from the elastic body (still with zero energy in the matrix) so as to relax to zero the layer of the body force traction   = −    = − *   on the boundary of the inclusion as it acts on the relaxation displacement.
However, with the application of loading, and specifically the pressure pre-stress p 0 δ ij , the driving force on the RHS of Equation 11 with the total stress , and it can become zero for a critical nucleation pressure, at which the outward flow due to the positive driving force   *  cancels the inward flow due to the self-stresses, so that indeed, there is no interaction energy of the expanding inclusion with the pre-stress field of the pressure, and an arbitrarily small inclusion can nucleate quasi-statically and grow at constant potential energy, as can be shown by relating the interaction energy to the potential energy.Equation 11 is the expression derived in Markenscoff (2014Markenscoff ( , 2020) ) for the M integral for an inclusion (e.g., Budiansky & Rice, 1973) for invariance of the potential energy functional under an infinitesimal transformation in scaling  λ , with     =   , (see Text S3 in Supporting Information S1), The M integral for the self-similarly expanding surface of discontinuity is equal to minus the change of the potential energy of the system with scaling, which is equal to the change in the interaction energy, and Equation 11 is obtained in Equation S3.6 in Text S3 in Supporting Information S1 as Figure 3. Energy flows across the moving phase boundary of an inclusion with eigenstrain under pressure from the bulk into the inclusion and from the inclusion into the bulk.The inclusion extracts energy from the bulk to accommodate its own self-stresses in the matrix, indicated in blue, while the pressure acting on the volume collapse   *  produces energy flow out of the inclusion boundary into the bulk indicated in red.At a critical pressure, when the two flows are equal and cancel each other, there is obviously no interaction energy between the inclusion with eigenstrain and the pressure loading mechanism, which allows an arbitrarily small inclusion to nucleate and grow at constant potential energy, as the variation in the interaction energy is related to the variation in the potential energy, This is the phenomenon of the growth of an inclusion governed by the vanishing of the M integral (Equation 12).It obtains the radiated energy.MARKENSCOFF 10.1029/2023JB026797 11 of 16 Equation 12 is an equation of motion for the expanding ellipsoid determining the axis speeds.
The radiated energy E R for surface earthquakes, modeled with the mechanics of moving cracks/faults, was originally obtained by Kostrov (1974), subsequently by Husseini and Randall (1976); in Rudnicki and Freund (1981) it was given in their Equation 46as ) , which is the input energy minus the energy lost in the expansion process of the propagating crack.Similarly here, the radiated energy rate is equal to the "excess energy," of the rate of input energy of the pressure acting on the "volume collapse" minus the losses to nucleate the expanding pancake, which is the energy to generate the faces of the pancake; integrated over time the total radiated energy is obtained as Analogously to the "stress drop," in this volume collapse under high pressure seismic phenomenon, the difference can be called "pressure drop" and drives the expansion acting on the reduction in volume.As we do not have a solution for the equation of motion giving the expansion speed for a given "pressure drop" we cannot estimate the radius of the fault, and in Table 3 we take the radius of it from the literature (Ye et al., 2013).For a plane phase boundary with dilatational eigenstrain the equation of motion relating the input energy to the boundary velocity was obtained by Markenscoff (2010aMarkenscoff ( , 2010b)).
The radiated energy according to Equation 13 is analogous to the "excess strain energy" W 0 of Kanamori (1977) for great earthquakes, which is the work done by the fault surfaces during quasi -static fault extension.The evaluation of Equation 13with the values shown in Table 1 for the phase transformation assumed for the Okhotsk earthquake gives the radiated energy in Table 3.In Ye et al. (2013) the radiated energy is estimated at 1.5 × 10 17 J.Thus, the model predicts energies that can explain the observation with a significant margin that could account also for dissipation damping.Losses due to damping are not considered here.for two fault radii a 1 , with 50 km being a conservative estimate.With axes ratio 10 it would be 3.75 × 10 18 J.

Discussion
The deep earthquakes are modeled as instabilities of nucleation and growth in a dynamic phase transformation with densification ("volume collapse") under high pressure.Successive instabilities occur as to minimize the energy spent to move the boundary of phase discontinuity.An unstable mechanics solution under high pressure is obtained for the problem of Randall (1964) of a sudden volume change exhibiting here a new instability that produces a DC seismic source (driven by the high pressure) (Figure 1).
The shape (axis speeds) is dictated by the energetics of moving surfaces of phase discontinuity, and instabilities minimize the energy spent in the nucleation and growth process of the phase boundary in a stress field.There have not been such energy considerations in geophysics with regard to phase transformation models (e.g., Knopoff & Randall, 1970).Knopoff and Randall (1970) investigated models of seismic radiation from suddenly occurring phase transition and proposed model B for change in density and/or bulk modulus, and model C for a sudden change in shear modulus in the presence of uniaxial strain.They assumed intuitively spherical expansion of the seismic source concluding that model B would give isotropic radiation of P waves, and that model C would give a CLVD.However, by considerations of the "equilibrium shape" (Eshelby, 1970;Markenscoff, 2019a) of the boundary of the region of phase change, the problem of Randall (1964), and model B of Knopoff and Randall (1970) exhibit an instability (Figure 1), under high pressure, and acquire a flattened "pancake-like" shape, as to minimize the energy spent as the inclusion grows large, rather than the symmetry-preserving spherical shape (Markenscoff, 2020).Due to the symmetry in material properties, geometry, and loading, a symmetric CLVD shear seismic source was obtained (Markenscoff, 2021a).However, the outward radiation was not studied in Markenscoff (2021a) and the question remained that the vast majority of the observations for the DFEs are of a DC seismic source type.This fundamental issue was resolved here by obtaining a successive instability, which minimizes the energy spent to expand the perimeter as well.Markenscoff (2020Markenscoff ( , 2021a) obtained a critical pressure at which an arbitrarily small inclusion with transformation strain due to sudden densification can nucleate and grow at constant potential energy, driven by the pressure acting on the change in volume   *  .At this critical pressure the flow of energy from the inclusion into the matrix and from the matrix into the inclusion is balanced, indicating that there is no interaction between the inclusion and the matrix (Figure 3), which allows an arbitrarily small inclusion to nucleate and grow at no expenditure of energy.In an instability, the ellipsoidal inclusion with volume collapse eigenstrains   *  assumes a flattened pancake-like shape allowing it to grow large at less energy expense.The flattened shape creates deviatoric eigenstrains; the mean eigenstrain of volume collapse induces a tensile mean stress in the inclusion that is canceled by the pressure that is high enough at the depth of the DFEs, as shown in Table 1 for the corresponding phase transformations; the strain energy density in the inclusion is predominantly distortional.The interior of the inclusion is a lacuna with zero particle velocity and no kinetic energy dissipated inside, and where the phase transformation can take place under equilibrium conditions.Due to symmetry in geometry and loading the eigenstrains are of a symmetric shear center, CLVD (CLVD in geophysics) as obtained in Markenscoff (2021a).The emitted radiation was not analyzed in Markenscoff (2021a) and it was assumed to be of the CLVD type as the eigenstrain in the inclusion.However, a CLVD would consume more energy to grow and expand rather than the DC obtained here.Some specific comments on the solution are: • Regarding the non-unique decomposition of the CLVD to two DCs, it was suggested by Knopoff and Randall (1970), Equation 12, that the non-uniqueness be removed by assuming a symmetric decomposition as (−1/2,−1/2,1) = (1/2) (0,−1,1) + (1/2) (−1,0,1).It should be noted that in spherical expansion, the CLVD of Model C in Knopoff and Randall (1970), produces zero interaction energy between the fields of the CLVD and the DC in the direction in which it is proposed that their axes be taken parallel.Thus, for a CLVD due to their Model C (change in shear modulus due to uniaxial strain), the proposed parallelism to remove the non-uniqueness in the CLVD to DC decomposition is justified.However, if applied here, this decomposition would resist the radiation through the perimeter of a flattened ellipsoid with volume collapse requiring non-zero interaction energy.• The nucleation pressure for the inclusion with transformation strain is obtained at the value for which the interaction energy of the inclusion with the loading mechanism of the applied loading (pressure) vanishes.This is shown by deriving the flow of energies across a moving surface of phase discontinuity from the matrix into the inclusion and from the inclusion into the matrix in Equations 11 and 12.The inclusion extracts energy from the matrix to be accommodated into the bulk, while the pressure acting on the volume collapse 10.1029/2023JB026797 13 of 16 eigenstrains produces energy flow into the bulk.The mismatch of the two is the interaction energy, and, if it is zero, then the inclusion can grow at constant potential energy according to the M integral.It also exemplifies the importance of the concept of eigenstrains in accounting for the energetics of defects in a pre-stress field governing their growth and shape, and valid also in fluids and a variety of phenomena in continua, including pressurized magma in volcanic eruptions.The radiated energy is obtained according to Equation 13, with the input energy due to the ambient pressure p 0 times the change in volume   *  , reduced by the energy spent for the nucleation and growth of the surface of phase discontinuity, which are the surfaces of the pancake; again, there is no additional energy loss in expanding the perimeter with a DC placed on the basal plane.Due to the lacuna property of the self-similarly expanding ellipsoid, there is no kinetic energy dissipation in the interior domain, which allows the phase transformation to take place under equilibrium conditions, and allows more energy to radiate out.To note that, unlike fracture, which introduces a length scale spoiling the self-similarity (see discussion in Burridge & Willis, 1969), in the expanding inclusion self-similarity is preserved and no characteristic length is introduced (any damping losses would be proportional to the boundary velocity).The driving force for the expansion of the inclusion under pressure is provided by the "pressure drop" analogous to a "stress drop" for cracks.The radiated energy is calculated on the basis of Equation 13 when the radius of the fault is taken as estimated in the literature (in Ye et al., for the Okhotsk, 2013 quake) with values shown in Table 3.Thus, the model predicts energies that can explain the observation with a large enough margin that could account also for losses due to dissipation damping.
• In this symmetric model, the pancake is assumed circular; however, if the basal plane eigenstrains were unequal   * 11 ≠  * 22 (elliptical shape, two different axis speeds), then an additional DC would be produced on the basal plane as shown in Markenscoff (2021a, Equation 4b).
• Changes in the shear and bulk modulus can be treated by the same Eshelby Tensor (also, see Markenscoff, 2019aMarkenscoff, , 2020)); here only the change in density was considered as to demonstrate the shear generation even in the absence of change in shear modulus.As deep earthquakes involve the phase transformation of anisotropic materials (e.g., Li et al., 2018) the treatment can be extended to the anisotropic dynamic Eshelby problem by Fourier transforms (Willis, 1971).All steps in the current analysis, including the M integral, are valid in anisotropy, but here only material isotropy is considered as it captures the main physical features of the deep earthquakes, as due to the very high pressure acting on the reduction in volume (due to increase in density).• In the presence of heat conduction (but no other non-mechanical transfer modes) U in Equation 9 can be replaced by ∂U/∂t = ρ∂e/∂t + ∂q j /∂x j where e is the internal energy per unit mass and q j are the components of the heat flux vector.Then Equation 9is regrouped as  ( ̇ − ) − ( +  )∕ = 0 , which, with the divergence theorem, yields the energy flux expression of Freund (1990) and Kostrov and Nikitin (1970).The expression for the energy flux can separate into a recoverable and a dissipative part (Freund, 1990, p. 230).
• The deformation fields are those of the dynamic Eshelby ellipsoidal inclusion problem, self-similarly expanding, which is a lacuna with zero interior particle velocity.As the lacuna phenomenon and the dynamic Eshelby inclusion properties are also valid for a Newtonian fluid transforming into another Newtonian fluid (Bilby et al., 1975, including inertia effects shown in Markenscoff, 2021b), the current analysis can be extended to hydrous self-similarly expanding ellipsoids with transformation strain, and analogously analyze instabilities in dehydrating poroelastic solids coupled with volume collapse under high pressure, that can model intermediate earthquakes (70-350 km), for example, Ferrand et al. (2017), Jung et al. (2004), and Meade and Jeanloz (1991).It is conjectured that with the porosity depending on the eigenstrain of volume collapse, as with plastic porosity in Segall and Rice (1995), the critical instability pressure will be decreased inducing seismicity at lower critical pressures.

Conclusions
The DFEs (400-700 km) are modeled as events involving phase transformations under high pressure for which successive instabilities occur as to minimize the energy spent to move the boundary of phase discontinuity.The instabilities can resolve the "mystery" of how the earthquake can initiate and radiate out unhindered by the high pressures, and explain their large energies.It was shown that a suddenly collapsing volume (densification of 5%-10%) under high pressure, even in full isotropy can radiate out the energy as an anti-symmetric DC shear

Figure 1 .
Figure 1.In a series of instabilities a sudden change in volume under high pressure and under full isotropy produces a self-similarly expanding flattened Eshelby ellipsoidal inclusion that radiates as a Double Couple (DC) deep earthquake.(a) Problem ofRandall (1964): Stable solution of sudden volume change preserving symmetry and emits volumetric radiation; "it is a more concentrated source of energy than transformation faulting"(Randall).(b) Under high pressure, and high densification   *  , in full isotropy, successive symmetry breaking instabilities minimize the energy spent to move the boundary of phase discontinuity which expands planarly as (c) a densified concentrated "blade-like" inclusion containing large shear transformation strains.The eigenstrains decompose (b) into two DC shear sources: the DC1 with eigenstrains on planes parallel to the basal plane radiating through the expanding perimeter without energy loss, while the DC2, on perpendicular planes, assumes the expansion loss.(d) The direction of the displacement around the pancake is indicated corresponding to pressure P (positive and negative) and shear S waves, emanating the radiation of a DC.

Figure 2 .
Figure 2. Symmetry-breaking instability in the deformation of the collapsing volume under pressure in full isotropy inside an Eshelby flattened ellipsoidal inclusion containing uniformly distributed centers of eigenstrain: the non-unique decomposition of the compensated linear vector dipole (CLVD) deformation C (−1/2,−1/2,1)   *  into two Double Couple (DC) (either the ones in green, or in red color) allows to place one, the DC1 (1/2) (1,−1,0), on planes parallel to the basal plane, with zero interaction energy with the symmetric pancake stresses, and allowing the expansion of the perimeter without losses.The DC2 (−1,0,1) is a non-radiative shear localization that consumes the energy expense to nucleate and grow the pancake surfaces.This explains the DC radiation.
, the eigenstrains are those on the basal DC1   1∕2)(1 + )∕(1 − 2) *  , and the integration is over the volume Ω of the ellipsoid, yielding for the second first term in Equation7

in the Phase Transformation Indicated by the Arrows, and With the Ambient Pressure Being p 0 at the Depth of the Occurrence
of the Earthquake in the Mantle-Transition-Zone (MTZ); the Mean Stress σ m in the Volume Collapse Transformation Is Tensile, Given by  − *  , Where K is the Bulk Modulus in the Untransformed Configuration, and Is Canceled in the Above Transformations by the Critical Pressure so That No Volumetric Radiation Is Emitted; Values Are Indicated at the Discontinuities of 410 and 5 of 16For the transformations indicated in Table

Table 2
Definitions of Strain, Stress and Strain Energy Density for an Eshelby Inclusion